Estimating scattering potentials in inverse problems with Volterra series and neural networks

Inverse problems often occur in nuclear physics, when an unknown potential has to be determined from the measured cross sections, phase shifts or other observables. In this paper, a data-driven numerical method is proposed to estimate the scattering potentials, using data, that can be measured in scattering experiments. The inversion method is based on the Volterra series representation, and is extended by a neural network structure to describe problems, which require a more robust estimation. The Volterra series method is first used to describe the one-dimensional scattering problem, where the transmission coefficients, and the phase shifts are used as inputs to determine the unknown potentials in the Fourier domain. In the second example the scattering process described by the radial Schrödinger equation is used to estimate the scattering potentials from the energy dependence of the phase shifts, where neural networks are used to describe the scattering problem. At the end, to show the capabilities of the proposed models, real-life data is used to estimate the 3S1\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${}_{}^3 S_1$$\end{document}NN potential with the neural network approach from measured phase shifts, where a few percent relative match is obtained between the measured values and the model calculations.


Introduction
Inverse problems in low-energy nuclear physics are related to the determination of (generally) local scattering potentials from observables such as scattering phase shifts, transmission coefficients, total-and differential cross sections, spin observables etc. [1][2][3][4][5]. These problems are generally ill-conditioned, and in many cases even ill-defined [6,7] and require a careful examination of the actual problem on which the inversion is done. In low-energy nuclear physics, usually the interaction potential is sought from elastic scattering data, a e-mail: balassa.gabor88@gmail.com (corresponding author) which could be measured in colliding experiments [8]. At high enough energies, where inelastic channels could open, the potentials could also contain imaginary parts, that takes into account the absorption of the incoming flux, thus able to describe non-elasticity in the scattering experiments. Nuclear interactions are generally very complicated in nature, however at low energies it can be described fairly well, with optical potential models, which could take into account the elastic-and inelastic interactions as well [9,10].
There are several methods for estimating the interaction potentials from measured data that can be used in different scenarios. One such possibility is the original Gelfand-Levitan-Marchenko method, which uses the energy dependence of the phase shifts of a specific partial wave, to determine the scattering potential [11]. Another type of inversion is the Newton-Sabatier method, which solves the inversion problem at fixed energy from measured phase shifts of the different partial waves [12]. Apart from the fixed-angular momentum and fixed-energy cases, there also exist mixed case inversion techniques, where the phase shifts of several angular momenta, at different energies are used to obtain the potentials [13].
In this work, two different methods are proposed, which are both able to estimate the interaction potential from measurable data in quantum scattering problems. The first method uses a non-causal Volterra series representation to make a connection between the observables and the potentials in the Fourier domain, while the second method gives a sequential method to estimate the potentials at discrete space-points with neural networks, using only the energy dependent phase shifts at fixed angular momenta. The methods described here are purely data-driven and the identified systems are constrained to an operating range, in which the training has been done. If the differential equations, which govern the inverse systems are known, then it could be possible to analytically express e.g. the Volterra kernels [14], thus giving a more robust estimation, in a much wider operating range. The inverse system, however is rarely known, and is usually highly nonlinear, in which case it is necessary to constrain ourselves into a well-defined range, where the nonlinearities could be locally modeled.
The paper is organized as follows. In Sect. 2 the necessary mathematical tools are summarized, which contains the Volterra series representation of nonlinear dynamical systems, and the applied neural network model as well. After the short summary, in Sect. 3 the one-dimensional scattering problem is addressed, where a non-causal Volterra representation is used to estimate the scattering potentials in the Fourier domain, using the phase shifts and the transmission coefficients as inputs. In Sect. 4 the low-energy elastic scattering is addressed, where a neural network model is used to describe the interaction potential in low-energy nuclear scattering, from the energy dependent phase shifts as inputs. The proposed neural network model is then used in Sect. 5 to describe the low-energy neutron + proton scattering data, giving an approximate 2% relative match between the measured-, and the obtained phase shifts calculated from the potential after inversion. At the end, Sect. 6 concludes the paper, giving some further possibilities on how to use and extend the proposed methods.

Mathematical models
In this section the necessary mathematical tools are summarized, including the causal and non-causal Volterra series, and the multilayer perceptron (MLP) type neural network structures, which will be both used in the later sections. Both mathematical models are able to identify nonlinear dynamical systems, and were used in many engineering applications, where the underlying mathematical models are too complex or not known exactly [15][16][17].

Volterra series
In linear system theory, where a physical system can be described by linear differential equations, the system can be characterized by a simple convolution integral, where an input u(t) is convolved with a so-called transfer function h(t), giving the y(t) output of the system as follows [18]: where y(t) is the output of the system, u(τ ) is the input, while h(τ ) is the transfer function. Here, the value t is some parameter e.g. time, position, or any other variable on which the inputs, and outputs could depend. This simple representation can be extended even if the physical system is described by a system of linear differential equations, in which case a transfer matrix can be introduced [19]. In the linear case the behavior and the stability of the system is fully characterized by the transfer function, thus the main problem is to identify the h(t) kernel. In real-life, however, many physical systems are nonlinear in nature, and the simple linear representation is simply not valid. In this case many simplifying assumptions can be made e.g. one could take the best linear approximation (BLA) of a weakly nonlinear system in the least-squares sense [20,21]. The best linear approximation could be enough in many practical applications, however if the nonlinearities are dominant, other methods are necessary to adequately model the nonlinear system. One of these methods is the Volterra series representation of nonlinear dynamical systems [22,23], which is the extension of the linear theory by introducing higher order polynomial terms in the convolutional integral, and can be written as: where h n (τ 1 , . . . , τ n ) is the nth order Volterra kernel of the system. In the most general case the integrals are going from −∞ to ∞, in which case the model is called non-causal. In this case all the future and past values are used to describe the system, however, as it was mentioned before, the parameter t does not have to be time, but can be position or any other variable on which the system depends. If the system is causal, then the integral limits can be changed to [0, ∞], which greatly simplifies the identification problem. As it will be shown later, the non-causal representation is a very useful tool in the addressed problems, which is also shown in [24], where a discrete, non-causal Volterra model is used to describe inverse quantum mechanical problems. In practical applications, the system is usually discretized and the mathematical model could be cast into a system of linear equations, where the coefficient matrix contains the inputs with different lags, while the unknown parameter vector contains the identifiable kernel functions at the discrete points [25]. The identification of the kernel functions h 1 , h 2 , . . . , h n are quite the cumbersome task, if the order, and the memory of the system is large, in which case alternative methods exist to identify the kernels [26,27]. One particularly interesting Fig. 1 The schematic view of a multiple input-multiple output, feedforward, multilayer perceptron model, with h hidden layers method is the neural network approach, where different neural network structures are used to obtain the nth order Volterra kernels by Taylor expanding the activation functions around the bias values [28].
In the applications in this paper, multiple input and single output (MISO) systems are in focus, in which case the Volterra representation changes to: where u i (t) is the ith input variable, while the H N {·} multiinput kernel functions are expressed as: As it can be seen, in the MISO case, the inputs can be mixed together giving higher order mixed terms e.g the second order h 2 (τ 1 , τ 2 )u 1 (t − τ 1 , t − τ 2 ), which could be included in the identification process.

Neural networks
The neural network model is a different approach, which can be used to describe nonlinear dynamical systems. In this case the system is described by many interconnected layers built up by so-called neurons, which could have nonlinear or linear activation functions as well [29,30]. In this paper only feed-forward multilayer perceptron structure is used, therefore in this section, only that will be described in more detail. The schematic view of a fully-connected, feed-forward MLP network can be seen in Fig. 1, where for the sake of generality the multiple input and multiple output (MIMO) case is shown. The number of inputs are set to N , while the number of outputs are set to M. The mapping from u 1 , u 2 , . . . , u N to y 1 , y 2 , . . . , y M are described by the block structure shown in Fig. 1, where each line corresponds to a weight parameter, and each rectangular block corresponds to a neuron, which is decribed by an activation function and a bias parameter. The first h columns are the so-called hidden layers, while the last layer labeled by o is the output layer. The input-output relationship can be expressed in a simple compact form if one introduces vector and matrix notations, in which case the output vector y = [y 1 y 2 . . . y M ] T can be expressed as: where N h = [n h1 n h2 . . . n hk h ] T is a vector of the activation functions at the hidden-layer h (or at the output layer o), T is the vector of bias paremeters, and W h is a weight matrix constructed by the elements w h i j . This feed-forward structure is usually trained by the back propagation technique [31], which is a very efficient way in calculating the derivatives, that is needed for solving the optimization problem for the unknown weights and biases. A very good property of the neural network approach is its generalization capability, which is the ability to predict the output to previously unseen input data. With a careful construction of the network, very good generalization properties can be achieved [32] for different network structures as well.

One-dimensional scattering with a non-causal Volterra model
The first problem, which is addressed with the Volterra formalism will be the one-dimensional quantum scattering problem, where an incoming particle, described by the wave function φ(x) = e ikx is scattered on a bounded potential V (x), giving a scattered wave Ψ (x). The scattering process can be described by the Lippmann-Schwinger equation [33], which in one-dimension has the following form: where G(x, x ) is the Green function of the Schrödinger operator in one-dimension [34], k is the wave number, while φ(x) is the incoming-, and Ψ (x) is the outgoing wave function.
Assuming that the potential goes to zero at some finite distance, after discretization the integral equation can be rewritten into a linear system of equations, then can be solved by standard numerical methods. In the inverse problem, the unknown potential is sought from observables, which can be measured in scattering experiments. In [24] the above prob- Fig. 2 The block diagram of the inverse system lem was examined with a non-causal Volterra model using the real and imaginary parts of the scattered wave as inputs of the system. In an experimental setup however, the scattered wave is not accessible and the typical observables are the asymptotic phase shifts and/or transmission/reflection coefficients [35,36]. In this section, a non-causal Volterra system is used to describe the scattering potential in the Fourier domain, using the scattering phase shifts, and transmission coefficients as inputs. The block diagram of the system in question can be summarized in Fig. 2, where the input-output relationship is described by a dynamical system, which is essentially nonlinear. The nonlinear nature of the inverse scattering problem will be more apparent in Sect. 4, where the phaseshift → potential relation can be described by a first order nonlinear differential equation. In this section, a closed form of the dynamical system will not be derived, and a purely data-driven technique is proposed, which is able to describe the inverse problem sketched in Fig. 2. The two inputs of the system are the k-dependent phase-shifts, and transmission coefficients, which is defined as the ratio of the transmitted intensity to the incident intensity: where T (k) is the transmission coefficient, Ψ T R is the transmitted wave function, Ψ I N is the incoming wave function. An example for the system inputs and outputs can be seen in Fig. 3, where the potential V (x) and the corresponding Δφ(k) and T (k) functions are also shown.
It is apparent from the block diagram in Fig. 2, that the output of the system is not V (x), but instead its Fourier transform V (k), which is a more natural choice if the inputs of the system are all k-dependent quantities. The non-causal nature of the applied Volterra system means, that after discretization, the output V (k i ) depends on the previous the current k i , and the following values k i+1 , . . . k i+M as well, which can be summarized as follows: where Δφ is the phase shift, T is the transmission coefficient, H  To identify the kernels, 12,000 training samples have been generated by solving the Lippmann-Schwinger equation, which will be the operating range of the identified system. The potentials used for training were uniformly generated random samples, in the predefined range, which is a common technique in system identification to excite the system [37]. There are more refined methods for generating excitation signals e.g. random phase multisines [38], which could have been also used in this case as well. In addition, only symmetric potentials were generated so that the Fourier transform of V (x) will be purely real.
Here, theh = c = 1 natural system of units is used, and the distance is measured in inverse energies. The Fourier transform for the potential has been calculated between k ∈ [0, 20] MeV, with Δk = 0.2 MeV, so the Volterra system is discretized into 100 discrete points in momentum space. The choice of this discretization in momentum space corresponds to δx ≈ 0.163 MeV −1 in position space, which is approximately the maximum resolution, what we can get after inverse Fourier transforming the obtained V (k) potential. 1 To each wave-number k i , one Volterra system is identified. This does not take away from the generality of the system, as it is always possible to extend the model by introducing an extra input (in this case the wave number k), so that only one Volterra model is needed to describe the inverse system. It is, however, much simpler to identify separate Volterra systems for each wave number, due to the complicated cross-terms, which would arise otherwise. Identifying the kernels, it turned out that second order kernels, with M = 35 memories, and first order cross-terms were sufficient to describe every system, giving a few percent relative error at all of the wave numbers, for which the mean squared error (which is used in the optimization process, and is defined in Eq. 9) can be seen in Fig. 4.
where N is the number of samples, V true is the true potential, and V est is the estimated potential given by the second order Volterra model, which has been identified at k j . In Fig. 5 the results can be seen for a randomly generated test potential, discretized for 50 points between x ∈ [−4, 4] MeV −1 in position space, where a remarkable agreement is achieved between the estimated and the true potentials in Fourier-, and also in position space, after inverse Fourier transforming the obtained V (k) potential. For better visibility, the true potential in Fig. 5 is linearly interpolated between the 100 inversion points in the Fourier space, and also between the 50 discrete points in position space, which are used in the inverse Fourier transform. The model works for more complicated test potentials as well, which can be seen in Fig. 6. Here, some larger discrepancies can be seen near k ≈ 6 MeV, which corresponds to the largest MSE values in Fig. 4, and therefore is expected. With more memories and/or including higher order terms and cross-terms, the model can be further improved.
In this operating range, the nonlinearities are quite negligible, as only second order systems were enough to cover them, with a good accuracy. If the magnitude of the potential is larger, or the potential covers a wider range in space, the simple second order approximation might not be enough, and one has to include higher order terms as well. It was also important to include, at least the first order cross terms, in which case the combinations of the two inputs at different k i -s e.g. Δφ(k i )T (k j ), are also appearing in the Volterra representation.
In the next section, we go a step further, and only the phase shift information will be used at a wider operating 1 The maximum resolution is the consequence of the Nyquist criteria. range, where a different approach is more suitable to solve the inverse problem in low-energy elastic scattering problems, which arises in low-energy nuclear physics.

Low-energy elastic scattering with neural networks
In this section, a new data-driven method will be introduced, which is able to determine the scattering potentials in lowenergy elastic scattering experiments, when the only observables are the energy dependent s-wave phase shifts. In this regard, the inversion is a fixed angular momenta method, where the energy dependence of the phase shifts are used at fixed angular momentum to describe the inverse scattering problem. The method is universal, and it could be easily extended to different angular momenta as well, which is a future goal of the proposed method. In the previous section a Volterra model is used to describe the inverse problem, while here the method is extended with the use of a feed-forward neural network structure, which is able to grasp higher order nonlinearities. After describing the neural network model, and showing its working principles through simple examples, it will be used to describe real-life scattering data in Sect. 5.
Here, first a short introduction is given of the theory of lowenergy quantum scattering, where we solely focus on spherically symmetric potentials, but now in three-dimensions, so that the method can be readily used and compared with real experiments as well. In low-energy nuclear scattering the two-body scattering system with masses m 1 and m 2 can be described by the three-dimensional time-independent Schrödinger equation: where μ = m 1 m 2 /(m 1 + m 2 ) is a reduced mass of the two colliding particles, Ψ (r) is the normalized wave func- is the Laplace operator, while V (r) is the potential, and E is the energy. If we only consider spherically symmetric potentials and separate the radial and the angular parts as Ψ (r) = R l (r )Y m l (θ, φ) the problem can be simplified into solving the radial Schrödinger equation, which has the form of: where u l (r ) = r R l (r ) is the radial wave function, and l is the angular momentum of the system. In scattering experiments, we are interested in the asymptotic behavior of the wave functions, which can be expressed as the sum of an incoming plane wave e i kr and an outgoing spherical wave, weighted by the scattering amplitude f (k, θ), which is related to the differential cross section as [39,40]: The scattering amplitude can be expressed with the scattering phase shifts by doing a partial wave expansion, giving the following form: where δ l is the phase shift of the lth partial wave, and P l is the lth order Legendre polynomial. In scattering experiments the differential-, and total cross sections are measured, and the phase shifts are determined by fitting the different partial waves [41,42]. The forward problem is therefore to solve the Schrödinger equation, and by fitting the expected asymptotic wavefunction, to obtain the phase shifts. The Schrödinger equation, which describes the problem is linear, however the inverse problem, which appears in measurements are nonlinear, as it can be seen from the so-called Variable Phase Approximation (VPA) [43][44][45][46], which is a first order nonlinear differential equation, which relates a compact and local scattering potential to the asymptotic phase shifts as: where δ l (r ) is the accumulated phase shift at r , while U (r ) = 2μV (r )/h 2 , and k = (2μE/h 2 ) 1/2 . The j l and n l functions are the Ricatti-Bessel and Ricatti-Neumann functions. The initial condition for the VPA equation is δ l (r = 0) = 0, which corresponds to a zero phase shift, when the potential still not disturbs the incoming wave. The measured, accumulated phase shift will be the asymptotic δ l (r → ∞) value. This equation can be used in atomic and nuclear physics to obtain the scattering phase shift for a given potential [47,48], however for higher angular momenta the numerical solution is quite difficult, and a very precise solver e.g. fourth or fifth order Runge-Kutta method is necessary.
In this paper, only the s-wave (l = 0) scattering is considered, in which case the VPA equation becomes: In the inverse problem, the asymptotic δ(r → ∞) values at different energies are used as inputs to describe the unknown potential V (r ).
Knowing that the phase shifts can be described directly by a first order nonlinear differential equation, one way for inversion would be to give a closed form to the direct problem (if its possible) at different energies, then solve the nonlinear system of equations with a suitable method. In ideal case, the closed form would be the analytical solution of Eq. (15), however due to the nonlinearities, it is not a trivial problem to solve. In the system identification side, it is possible to estimate the behavior of the nonlinear dynamical system in a well defined operating range, with either Volterra series or neural networks. This is, however, not an efficient way to do the inversion, because one still has to solve a system of nonlinear equations to obtain the potentials.
In the Volterra system theory there are methods on how to obtain the inverse kernels from the identified direct kernels [49], while in the neural network case, there are several methods exist for inversion [50][51][52]. Most of these methods require a numerical inversion scheme, where e.g. gradient methods are used to determine the inverse solution. One especially interesting method is the use of an invertible neural network structure [53,54], where the inverse of the system could be uniquely determined from the forward network. This could be a convenient way to solve these kind of problems, because the forward process is usually much simpler to train, however in this case a very careful training is needed if the inverse system is badly conditioned. Another interesting method, which could be tested in the future is the use of onedimensional convolutional neural networks (1D CNN) [65], which could have nice properties e.g. insensitivity to small perturbations in the input data, and due to its convolutional form it could be well suited to describe dynamical systems as well. In the followings only multilayer perceptron type neural networks are used to describe the dynamical systems, and the usage of other types of networks is left to future works.
The method, which will be proposed, identifies the inverse system sequentially, using the phase shifts as inputs and the potentials as outputs. The schematic view of the inversion process can be seen in Fig. 7, where Δφ refers to the energy dependent s-wave phase shifts, while V (r ) is the potential. In the method the position space (distance from the origin r ) is discretized at M different space points, denoted by r 1 , r 2 , r 2 , . . . , r M , where r 1 is the closest point to the origin, while r M is the furthest. The inversion will be given in these so-called inversion points, which ideally spans the whole space, where the V (r ) potential differs significantly Fig. 7 The schematic view of the inversion process, where each N i block corresponds to a fully connected MLP network from zero. In this sense, we assumed that the potential goes to zero (or at least close to zero) at a finite distance, which is a reasonable assumption in nuclear physics. In practice, it is usually possible to give, at least a crude estimation to the range of the potential, thus this information could be easily put into the identification process. If nothing is known from the potential, then it also should not be a problem, however the identification process could take a longer time, because one have to span a larger region, with possibly more inversion points.
The method starts from the assumption, that a closed form can be given to the potential at the last inversion point V (r M ), using only the phase shifts as inputs. To describe this relation, a multilayer perceptron model is identified, which is marked as N 1 in Fig. 7. For the next point, V (r M−1 ), the phase shifts, and also the previously estimated V (r M ) potential value is used as inputs for the second dynamical system N 2 . The V (r M−1 ) potential still only depends on the measured phase shifts, but now through N 1 as well. In this way the identification process turned out to be much simpler. The method can be generalized for the following points as well, where the last inversion point V (r 1 ) will depend on the phase shifts, and all of the previously estimated potentials as well, which also depend on the input phase shifts, so ultimately we will arrive at a closed form between the phase shifts and the potentials in the inversion points. The method can be summarized as follows: V (r 1 ) = N M (Δφ 1 . . . , Δφ k , V (r M ), . . . , V (r 2 )), (18) where each N i (·) expression represents a multilayer perceptron model. The number of free parameters can be expressed by counting the weight and bias parameters at every MLP used in the inversion points as: where M is the number of inversion points, h (i) j is the number of neurons in the j'th hidden layer of the i'th MLP, k is the number of phase shifts used as inputs, while n (i) is the number of hidden layers in the i'th MLP.
To test the method in a simple scenario, let us assume, that the sought potential is bounded between V (r ) ∈ [−10, 0] MeV in r ∈ [0, 6] fm, and dies out at r > 6 fm, while for the mass, we fix m = m 1 = m 2 = 1 GeV. These values give a good benchmark to examine the method, as the nonlinearities in this case are becoming non-negligible, thus a nonlinear dynamical model is necessary. Here, instead of the wave number k, we introduce the laboratory kinetic energy T lab as: where T lab is the laboratory kinetic energy. In the following examples T lab is going to be used as the input parameter for the phase shifts, so Δφ(k) → Δφ(T lab ) exchange is used, to be able to compare the model results directly with the measurements in Sect. 5. To do the identification, 10,000 training-, and 2000 validation samples has been generated in the predefined range, where the phase shifts has been calculated by the VPA equation shown in Eq. (15). The identification has been done for different systems containing different number of hidden layers and/or neurons to be able to determine an optimal system complexity for the actual problem. The inputs for the system were the laboratory kinetic energy dependent phase shifts, Δφ(T lab ), from 1 to 99 MeV, with ΔT lab = 5 MeV steps, so in overall 20 points at T lab = [1, 6, 11, 16, . . . , 99] MeV. The number of inversion points is set to M = 20 between 0 and 5.713 fm, with Δr = 0.3007 fm. To train the system the Levenberg-Marquardt back propagation algorithm is used [63] sequentially for every inversion point one after another, starting from the last point at the largest distance. The inputs and outputs were normalized between [−1, 1] in every neural block, and to avoid overfitting, the training is stopped after the validation error starts to increase, or the gradient reaches a minimal value set to be 10 −7 .
To make a first guess for the complexity, four different systems were trained and compared by calculating the relative errors in the inversion points for 2000 test samples, which can be seen in Fig. 8, where the relative errors are shown in Three of the four systems consists only one hidden layer h (i) , with 2, 5, 20 neurons in each of the neural blocks N i , while in the last system two hidden layers were used, with 20 neurons in each. The activation functions were tanh(·) nonlinearities in the hidden layers, and a linear relationship in the output layer. For the first inversion point in r = 5.713 fm only the input phase shifts could be used for teaching, therefore one could expect a larger error, and a slower learning rate, which can be seen in the calculated relative errors as well. It can be seen that with increasing complexity, the relative errors are decreasing in each of the inversion points. With using 10 neurons and one hidden layer, the relative error is just a few percent in almost all of the points. In all cases, the largest errors can be seen in the last few points at small distances below 0.5 fm, which is the consequence of the large sensitivity of the inversion in this range. With using two hidden layers (or more neurons), the errors could be further reduced, however to achieve a very good inversion it is not necessary to use two hidden layers in each of the inversion points. The final model, which will be used in the followings consists one hidden layer, with [10] neurons in the first 16 inversion points r 20 , r 19 , . . . , r 5 , and [20,20] neurons in two hidden layers in the last 4 inversion points r 4 , . . . , r 1 , where the errors are the largest. Following Fig. 7 this means,  that N 1 ,N 2 ,. . .,N 16 are MLP's with one hidden layer and 10 neurons, while N 17 ,N 18 ,N 19 , and N 20 are MLP's with two hidden layers, each having 20 neurons. The relative errors in the inversion points for this system can be seen in Fig. 9, where again, a few percent relative error has been obtained in all of the points. In Fig. 10 the training process can be followed, where the dependence of the mean squared error on the number of epochs determined from the normalized inputs/outputs are shown in four different inversion points r 20 , r 13 , r 5 , r 1 , where the training and validation loss is shown as well. The figures clearly show the previously mentioned behavior regarding the time to reach an optimal solution. The first inversion point (r 20 ) needs more than 10 2 epochs to reach an MSE, which corresponds to a few percent relative error in that point, while in the subsequent points a convergence can be achieved relatively fast, after 10-20 epochs. Near the last inversion point (r 1 ) the MSE of the converged solution is a little larger than in the previous points, which is also expected from the previously determined relative errors in Fig. 9.
To further examine the model, the accuracy has been calculated using different number of training samples for the previously determined model, which can be seen in Fig. 11. The accuracy shows a sharp rise at a few hundred samples, then a slowly saturating region after a few thousand samples, which is the usual shape one will get with multilayer perceptron models [64]. This shows us, that in this system, with only a few thousand training samples a good generalization can be obtained with the proposed model.
To show the model capabilities through some examples, the results for two randomly generated potential compared with the original potentials can be seen in Fig. 12. The inversion technique is capable to give satisfying results within a few percentage of relative error, even with this simple neural network structure. It is worth mentioning, that the discrepancies, especially at small r -s, could be improved, by setting up more inversion points, and/or using more phase shifts as  inputs, however in that case the time of training and the complexity of the model could also increase.

Estimating the 3 S 1 N N scattering potential with neural networks
After introducing the neural network model through a simple example, let us proceed to the main application of the inversion method, which is the determination of the scattering potential in low-energy elastic S-wave N N scattering from the experimental phase shifts. In particular, we are interested in the 3 S 1 channel, in which case one bound state, the deuteron, appears in the measured phase shifts [55,56]. In the identification process, there are no restrictions to the number of possible bound states, however, knowing apriori the number of bound states, it could be possible to restrict the training samples to only consist potentials, which can give a specific number of bound states. In this case the training and the identification of the system could be greatly simplified, however the obtained model will be more restricted. Before generating the training samples, it is necessary to gather some information about the identifiable system, as it could greatly simplifies the identification process. In this case, it is known, that the effective range of the scattering process should be a few Fermi in position space, however it is also known, that strong repulsive, and Coulomb-like terms could also appear in the potentials. As it was described in the previous section, the first inversion point (furthest from the origin) does not have to be zero, but it is desirable, to be at least close to zero, as in that case, the identification is much simpler, and a simpler model could be used to describe that point. In this case, we can assume, that at approximately r ≈ 3 fm, the potential is negligible, so for first approximation, we set the range of the inversion to r ∈ [0, 3] fm, with M = 30 inversion points, starting from r 1 = 0 fm, and ending at r 30 = 2.9986 fm, which corresponds to Δr = 0.1034 fm.
The operating range for the inversion has been set, so the model could describe a bounded attractive and repulsive potential in the whole range, and also a highly repulsive term e.g. a Coulomb-, or a nonperturbative strong repulsive core, at low r -s near r = 0. The potentials are also assumed to go to zero at r ≈ 3-4 fm. Furthermore, it would be possi-ble to include infinities at r = 0 as well, in which case the inversion point at r = 0 fm has to be excluded. In this case, we assumed that the potential at r < 0.5 fm is a finite value between −500 and 20, 000 MeV, while at r > 0.5 fm the values are constrained in the range of V (r ) ∈ [−500, 500] MeV, and go to zero at 3-4 fm. In the generation of the excitation functions, or in other words the V (r ) training potentials, it is necessary to include the previous assumptions, while keeping the randomness of the generated functions. This is realized with the help of Piecewise Hermite Interpolating Polynomials (PCHIP) [57], which is a polynomial interpolation method, that is able to give continuous functions in a predefined range, so it is possible to generate V (r ) functions, where |V (r )| < V max constraint can be easily satisfied. To do the interpolation a number of control points has to be given, which are used as fixed points, in which between the interpolation is done. To keep the randomness of the potentials, the control points are chosen randomly between r c ∈ [0, r max ] for each V i (r ) excitation, by first generating a random number N , which gives the number of control points, then generating N number of r i points, which will be the place of the control points. At the final step a random number is generated for each control point in the operating range of V (r ), then we interpolate between the control points. Two randomly generated potential, which respects our previous assumptions are shown in Fig. 13, where it can be seen, that at small distances the potential could be very large, while in the region r > 3 fm it goes to zero. It can also be seen, that if the number of control points are large, then the potential keeps its randomness, that we desire in an identification process, where one optimally wants to excite all of the frequencies until the Nyquist limit. To train the neural network, 10,000 samples have been generated with different number of control points in each sample from N = 2 up to N = 18, while for validation another 2000 samples have been generated, which are also used in the training process to avoid overfitting. Following the reasoning of the previous example, the complexity of the network has been chosen, so that it could describe the system with a few percent relative errors in all of the inversion points. The final structure for all of the inversion points are MLP's with two hidden layers, having 30 neurons in each, with tanh(·) activation functions, and a linear output layer. In this model, every inversion point is described by this structure, so every N i (·) (i = 1, . . . , 30) has the same MLP structure, which has to be trained sequentially, starting from the first inversion point r 30 .
The training process for two representative inversion points (r 30 , r 23 ) can be followed in Fig. 14, where the expected behavior is seen again. The first inversion point r 30 , which only depends on the input phase shifts, the training process is slower and need more epochs to reach convergence, than what is observed in the following inversion points, which depend on the previous inversion points and the input phase shifts as well.
The obtained relative error in the 30 inversion points for 2000 test samples after training are shown in Fig. 15. The largest errors again occur in the first few inversion points at the largest distances, and at the last few inversion points at small distances, which was the observed behavior in the last example as well. In overall a few percent error could be achieved with this model. An example of the model capabilities to a randomly generated potential can be seen in Fig. 16.
After the inverse system is identified in the desired operating range, it is possible to estimate the low-energy N N scattering potentials from the measured phase shifts. In this  section the neutron+proton scattering is investigated between T lab ∈ [1, 300] MeV laboratory kinetic energies, for which the measured phase shifts are taken from [58]. After inversion, the obtained potential and the laboratory kinetic energy dependent, re-calculated phase shifts can be seen in Fig. 17, where a very good match has been achieved between the original-, and the re-calculated phase shifts. The averaged relative error, through the measured points, between T lab ∈ [1, 300] MeV is approximately 2%. The obtained potential has a strong repulsive core at small distances, while at larger distances a moderate attractive force emerges. This picture is comparable with the usual nuclear potential models for s-wave scattering [59,60], where the strong repulsive core is the consequence of the strong force at small distances, while the medium and large distance attractive terms are related to the various meson exchange forces [61,62]. It is worth mentioning, that we did not have to make any distinction for scattering with-, and without bound states, and could include them both in the identified model. The identified system in this form can only be used to describe l = 0 scatterings, however by adding the angular momenta to the inputs, it can be easily extended to describe potentials, with different angular momenta as well.

Conclusions
In this paper, two data-driven inversion methods are proposed, which are both able to describe inverse quantum mechanical scattering systems, in a well-specified operating range, using nonlinear system identification techniques. In these problems, the scattering potential is sough from the known observables e.g. phase-shifts, and transmission coefficients, using finite dimensional, truncated Volterra systems, and feed-forward neural networks. Two main structure is proposed, which are tested in two scattering scenarios. Firstly, a non-causal Volterra system is used to describe the onedimensional quantum scattering problem, where the inputs were the energy dependent phase shifts, and the transmission coefficients, while the output was the Fourier transformed scattering potential. The system was designed, so that in every k i (in the range, constrained by the Nyquist limit), one non-causal Volterra system is identified, which used the input functions at k i−M , k i−M+1 , . . . , k i , . . . , k i+M , defined by the memory of the system. To test the method, the operating range has been set to V (x) ∈ [0, 1] MeV, with a test mass of m = 1 MeV. With these parameters the nonlinearity of the system could be modeled by second order Volterra series, with M = 35 memories, by including only the first order cross-terms of the input phase shifts, and transmission coefficients.
For the second problem, a more realistic model is proposed, where the low-energy elastic scattering on spherically symmetric potential is addressed, with neural networks. In this model, only the energy dependent phase shifts were used as inputs to the inverse system, while the potential in position space is used as the output. To solve the inversion, a sequential neural network architecture is proposed, which consists one MLP network at each inversion point, which are connected by their outputs. To test the model, the operating range has been set to T lab ∈ [1, 99] MeV, and V (r ) ∈ [−10, 0] MeV, with test masses of m 1 = m 2 = 1 GeV. The network has been trained by generating a sufficient number of training samples, with the help of the Variable Phase Equation, then tested to randomly generated continuous potentials, giving good results in each case, with a few percent relative error.
The neural network construction then used to describe the low-energy elastic neutron+proton scattering from measured phase shifts between T lab ∈ [1, 300] MeV. To train the network, the operating range has been set to be able to describe the short distance repulsive-, and the larger distance attractive terms as well, where the training samples were constructed with the help of Piecewise Cubic Hermitian Interpolating Polynomials, which were able to easily constraint the potential functions between the prespecified range. The results of the inversion shows a physically reasonable potential, with a short range repulsive, and a long range attractive part. The phase shifts were calculated back from the obtained potential, and compared to the measured phase shifts, giving good results, with around 2% relative error, averaged through the number of measurement points.
The methods proposed here, were both able to solve the inversion problem in quantum scattering, at a predefined operating range. As the inverse system is essentially nonlinear, by extending the operating range, the nonlinearities will have greater roles, and a more complex system has to be identified, which in the case of the Volterra model, means more memories, and higher order terms, with more cross-terms as well, while in the neural network structure, it means more hidden-layers, with more neurons. This method is advantageous, if the full nonlinear dynamical equations, which governs the system, are not known, or one only wants to operate the system in a well-defined range. In this case one option would be to identify a model (Volterra/Neural Network) for the forward system, than by generating a sufficient number of samples, identify the inverse system as well, with one of the described methods. Another possibility would be to directly identify the inverse system by using only the measured data, which could be more problematic if there are no sufficient number of observables to work with. In low-energy quantum scattering experiments the dynamical equations, which governs the system are well-known, therefore there is no need to identify a black box system, and one can utilize the apriori knowledge in quantum mechanics to build a sufficient model even for more complex systems as well.
In the case, that the dynamical equations are known precisely, the Volterra kernels can be calculated analytically and the higher order terms can be obtained more easily. As the Variable Phase Equation describes the forward problem, it could be possible to analytically calculate the Volterra kernels for the forward process, then invert the obtained Volterra system, to be able to estimate the potentials from the phase shifts, which could be a future application of the proposed methods.