On a nonlinear Schr{\"o}dinger equation for nucleons in one space dimension

We study a 1D nonlinear Schr{\"o}dinger equation appearing in the description of a particle inside an atomic nucleus. For various nonlinearities, the ground states are discussed and given in explicit form. Their stability is studied numerically via the time evolution of perturbed ground states. In the time evolution of general localized initial data, they are shown to appear in the long time behaviour of certain cases.


Introduction
This paper is concerned with the study of solutions to a nonlinear Schrödinger (NLS) type equation which, in a specific non-relativistic limit proper to nuclear physics, describes the behavior of a particle inside the atomic nucleus.This equation is, at least formally (see Appendix A), deduced from a relativistic model involving a Dirac operator and, in space dimension d = 1, is given by (1) i∂ where φ ∈ L 2 (R, C) is a function that describes the quantum state of a nucleon (a proton or a neutron), α ∈ N * is a strictly positive integer and a > 0 is a parameter of the model.Note that equation (1) is Hamiltonian and has a conserved energy (2) Solitary wave solutions for this equation can be constructed by taking φ(t, x) = e ibt ϕ(x) with ϕ a real positive square integrable solution to the stationary equation The reasoning that the solution can be chosen to be real is the same as for the standard NLS equation.
Positive square integrable solutions of (3) can be seen as ground states solutions of (1) since they are minimizers of (2) among all the functions belonging to (4) where f + denotes the positive part of any function f .As shown in [4], this is the appropriate way to define ground states for this energy.Indeed, on the one hand, by adapting the arguments of [4], it can been shown that the energy E is not bounded from below in the set ϕ ∈ H 1 (R), R |ϕ| 2 = 1 .On the other hand, if ϕ ∈ X, one can show (see [4]) that |ϕ| 2 ≤ 1 a.e. in R. As a consequence, for any ϕ ∈ X, E[ϕ] ≥ − a α+1 .
In this paper we will prove the existence of solitary waves solutions to (1) for any value of α ∈ N * , give them in explicit form, and study numerically their stability as well as the time evolution of more general initial data with |φ| < 1.
In the physical literature, the most relevant case is given by α = 1 which leads to the cubic nonlinear Schrödinger type equation (5) i∂ Nevertheless, it could be mathematically interesting to investigate also the behavior of solutions for other power nonlinearities as for example the quintic nonlinearity (α = 2) which corresponds to the L 2 critical case for the usual NLS equation.For the latter equation, it is known that initial data with a mass larger than the ground state can blow up in finite time, see for instance [7] and references therein.
To our knowledge, the above model was mathematically studied for the first time in [3], where M.J. Esteban and S. Rota Nodari consider the equation which is the generalization of (3) for any spatial dimension d ≥ 1 and for α = 1.In particular, the existence of real positive radial square integrable solutions has been shown whenever a > 2b.Note that solutions to (6) do not have a simple scaling property in the parameter b as ground states for the standard NLS equation.This makes it necessary to study several values of b in this context.This result has then been generalized in [8], where the existence of infinitely many squareintegrable excited states (solutions with an arbitrary but finite number of sign changes) of ( 6) was shown in dimension d ≥ 2.
In [4] (see also [6]), using a variational approach the existence of solutions to ( 6) is proved without considering any particular ansatz for the wave function of the nucleon and for a large range of values for the parameter a. Finally, in [6], M. Lewin and S. Rota Nodari proved the uniqueness, modulo translations and multiplication by a phase factor, and the non-degeneracy of the positive solution to (6).The proof of this result is based on the remark that equation ( 6) can be written in terms of u = arcsin(ϕ) as simpler nonlinear Schrödinger equation.
The same can be done for (3).Indeed, by taking u := arcsin(ϕ α ), one obtains In Appendix B, we generalize the results of [6] for any α ∈ N * in spatial dimension 1 by proving the following theorem.
The paper is organized as follows: In Section 2 we derive the explicit form of solutions to (3) for any α ∈ N * whenever a > (α + 1)b > 0, and we show their behavior for various values of the parameters.The computation presented in Section 2 is justified in the Appendix B where the proof of Theorem 1 is done.In Section 3, we outline the numerical approach for the time evolution of initial data according to (1).This code is applied to perturbations of the ground states for various values of the nonlinearity parameter α and for initial data from the Schwartz class of rapidly decreasing functions.In Section 4, we discuss the generalization of the model in higher space dimension.Finally, a formal derivation of the equation ( 1) is presented in Appendix A).

Ground states
In this section we construct ground state solutions to the equation ( 1) and show some examples for different values of the parameters.
First of all, equation ( 3) can be integrated once to give where we have used the asymptotic behavior of ϕ for x → ∞.Putting ψ := ϕ −2α , we get from (9), which has for a = (α + 1)b the solution Here x 0 is an integration constant reflecting the translation invariance in x of the ground state and ψ(x 0 ) is chosen in order to have a C 1 solution to (10) defined for any x ∈ R. Using the translation invariance, we will assume in the following that the maximum of the solution is at x = 0, and then we put x 0 = 0.The solution to equation (10) for a = (α + 1)b leading to the wanted asymptotic behavior of ϕ will not be globally differentiable.
Summing up, with (11) we get the ground states for 0 < (α + 1)b < a in the form Let us point out that this construction will be further justified in Appendix B where Theorem 1 is proven.
As a concrete example we show the solutions (12) for a = 9 and various values of b < a/(α + 1).The solutions for α = 1 can be seen in Fig. 1.With b → a/2, the solutions become broader and broader and have a larger maximum.The peak near 1 becomes also flatter.For b = 4.499, the maximum is roughly at 0.9999 and almost touches on some interval the line 1.
For α = 2, 3 we get in the same way the figures in Fig. 2. It can be seen that the higher nonlinearity has a tendency to lead to more compressed peaks as in [1].But due to the missing scaling invariance of the ground states here, it is difficult to compare them.

Numerical study of the time evolution
In this section we study the time evolution of initial data for the equation (5).We study the stability of the ground state and the time evolution of general initial data in the Schwartz class of smooth rapidly decreasing functions for various parameters.
The results of this section can be summarized in the following Conjecture 2. The ground states of equation ( 5) are asymptotically stable if the perturbed initial data satisfy |φ(x, 0)| < 1.  3.1.Time evolution approach.The exponential decay of the stationary solutions makes the use of Fourier spectral methods attractive.Thus we define the standard Fourier transform of a function u, and consider the x-dependence in equation ( 1) in Fourier space.
The numerical solution is constructed on the interval x ∈ L[−π, π] where L > 0 is chosen such that the solution and relevant derivatives vanish with numerical precision (we work here with double precision which corresponds to an accuracy of the order of 10 −16 ).The solution φ is approximated via a truncated Fourier series where the coefficients φ are computed efficiently via a fast Fourier transform (FFT).This means we treat the equation, ( 14) and approximate the Fourier transform in ( 14) by a discrete Fourier transform.
The study of the solutions to ( 14) is challenging for several reasons: first it is an NLS equation which leads to a stiff system of ODEs if FFT techniques are used.Since a possible definition of stiffness is that explicit time integration schemes are not efficient, the use of special integrators is recommended in this case.But most of the explicit stiff integrators for NLS equations, see for instance [5] and references therein, assume a stiffness in the linear part of the equations.However, here the second derivatives with respect to x appear in nonlinear terms.Since we are interested in a fourth order method (for the accuracy needed for |φ| ∼ 1 as discussed below), we will nonetheless use an explicit approach (the standard explicit fourth order Runge-Kutta method) with very small steps for stability reasons, but also since the accuracy is needed in some of the studied examples.This appears to be more efficient than an implicit scheme where a nonlinear equation has to be solved iteratively at each time step.
An additional problem of equation ( 1) is the singular term for |φ| → 1.Since the equation is focusing, it is to be expected that for initial data with modulus close to 1 it will be numerically challenging since the focusing nature of the equation might lead for some time to even higher values of |φ|.The accuracy of the solution is controlled as in [5]: the decrease of the Fourier coefficients indicates the spatial resolution since the numerical error of a truncated Fourier series is of the order of the first neglected Fourier coefficients.The error in the time integration is controlled via conserved quantities.We use the energy (2) which is a conserved quantity of ( 1), but which will numerically depend on time due to unavoidable numerical errors.In the examples below, the relative energy is always conserved to better than 10 −6 since the stability conditions as discussed above enforce the use of very small time steps.
We test the numerical approach at the example of the ground state.Concretely we consider the ground state solution for α = 1, a = 9, b = 4.4 as initial data.We use N = 2 10 Fourier modes in x for x ∈ 5[−π, π] and N t = 10 5 time steps for t ∈ [0, 1].Note that though the ground state solution is stationary, it is not time independent.We compare the numerical and the exact solution, i.e., the solution to (12) times e ibt at the final time t = 1.This difference is of the order of 10 −14 as shown in Fig. 3.The relative conservation of the energy (2) is during the whole computation of the order of 10 −14 .This shows that the ground state can be numerically evolved with an accuracy of the order 10 −14 , and that the conservation of the numerically computed energy indicates at accuracy of the time integration.3.2.Perturbations of ground states.We first consider the stability of the ground states constructed in the previous section.To this end we perturb it first in the form φ(x, 0) = λϕ(x), where λ ∼ 1.We use N = 2 11 Fourier modes and N t = 5 * 10 5 time steps for t ∈ [0, 0.25], i.e., more than a whole period of the perturbed ground state.In Fig. 4 we show the solution for the perturbed ground state with λ = 0.99.It can be be seen that after a short phase of focusing a ground state with slightly larger maximum than the initial data is reached.In addition there is some radiation towards infinity.The reaching of a ground state is even more obvious from the L ∞ norm of the solution shown on the left of Fig. 5.Note that a comparison with the ground state corresponding to the presumed final state of the time evolution in Fig. 4 is difficult since in contrast to ground states of the standard NLS equation, the ground states here do not have a simple scaling property in b.The Fourier coefficients of the solution at the final time on the right of Fig. 5 indicate that the solution is fully resolved in x.If we perturb the same ground state as in Fig. 4 with a factor λ > 1 (such that ||λϕ|| ∞ < 1), we observe a similar behavior as can be seen in Fig. 6.As is more obvious from the L ∞ norm in Fig. 7 on the left, a ground state with slightly lower maximum than the initial data is quickly reached.The decrease of the modulus of the Fourier coefficients on the right of Fig. 7 indicates that the numerical error in the spatial resolution is of the order of 10 −8 .This shows that there are stronger gradients to resolve in this case than in Fig.The same initial data as above are perturbed with a localized perturbation of the form φ(x, 0) = ϕ(x) ± 0.001 exp(−x 2 ).The resulting L ∞ norms of the solutions to (5) for these initial data are shown in Fig. 8.In both cases the L ∞ norm appears to approach a slightly smaller ground state (for the − sign in the initial data) respectively slightly larger ground state (for the + sign).In any case the ground states appear to be stable also in this case.
We repeat the experiments of Fig. 4 and Fig. 6 for α = 2, i.e., a higher nonlinearity.The L ∞ norms for the perturbed ground states can be seen in Fig. 9. Again a ground state with slightly smaller maximum state is reached for λ = 0.99, and with a slightly larger maximum for λ = 1.001.The oscillations which appear for larger times are due to us studying the perturbations in a periodic setting and not on R.  of the computational domain and interacts after some time with the ground state which leads to periodic excitations of the latter.Note that the quintic nonlinearity is L 2 critical in the standard NLS equation, i.e., solutions to initial data of sufficient mass blow up in finite time.Here no blow-up is observed, ||φ|| ∞ < 1 for all times.But in addition the ground state appears to be stable also for the Gaussian perturbations of Fig. 8   The same experiments as above are shown for an even higher nonlinearity α = 3 in Fig. 10.Once more the ground states appear stable (also for Gaussian perturbations not shown here).Note that for the standard NLS a septic nonlinearity would be supercritical which again would lead to a blow-up of initial data of sufficiently large mass.

3.3.
Schwartz class initial data.An interesting question in this context is whether these stable ground states appear in the long time behavior of solutions to generic localized initial data.To address this question we consider initial data of the form φ(x, 0) = µ exp(−x 2 ) with 0 < µ < 1, again for a = 9.We use N = 2 12 Fourier modes for x ∈ 5[−π, π] and N t = 5 * 10 5 time steps for the indicated time intervals.In Fig. 11 it can be seen that the L ∞ norm of the solution appears to oscillate around some asymptotic values, and that some radiation is emitted towards infinity.
The former effect is more visible on the left of Fig. 12 where the L ∞ norm of the solution is shown.Since there is no dissipation in the system, the final ground state will be only reached asymptotically.The small oscillations in the L ∞ norm for larger times are again due to the fact that the problem is treated as being periodic.Radiation emitted towards infinity reenters because of this periodicity on the other side of the computational interval if it leaves on one side.This     5) with a = 9 for the initial data φ(x, 0) = 0.9 exp(−x 2 ).
is more visible on the right of Fig. 12, where the solution at the final time of the computation is shown.The situation is similar for higher nonlinearity.In Fig. 13 we show the case α = 2. On the left for µ = 0.9, the L ∞ norm of the solution appears to be decreasing.The initial data seem to be simply radiated towards infinity.On the right we show the case µ = 0.96 where the L ∞ norm seems to oscillate around some final state with non-vanishing L ∞ norm, presumably a ground state.This would indicate that the ground state can appear in the long time behavior of generic solutions.
For even higher nonlinearity α = 3, we consider in Fig. 14 the case µ = 0.9 on the left and µ = 0.99 on the right.In both cases the initial data appear to be simply radiated away to infinity.Note that this does not imply that ground states cannot be observed in the long time behavior of solutions in this case, it just shows that this is not the case in the studied examples.For instance it is possible that Gaussian initial data do not have a sufficiently broad maximum in this case.1) for α = 3 and the initial data φ(x, 0) = 0.9 exp(−x 2 ) on the left, and φ(x, 0) = 0.99 exp(−x 2 ) on the right.

Outlook: analysis of the model in higher space dimension
It seems interesting to investigate the model described in this paper in higher space dimension, d = 3 being the most relevant case from a physical point of view.
On the one hand, to generalize the model for any dimension d > 1, one can simply replace ∂ x in equation ( 1) by the operator ∇.This leads to a quasilinear Schrödinger equation of the form ).On the other hand, at least in dimension d = 2 and d = 3, another possibility is to formally derive the equation of the model by following the arguments presented in Appendix A and by taking as starting point the Dirac equation in dimension 1 < d ≤ 3.This will lead to a slightly more complicated quasilinear Schrödinger equation.
In both cases, solitary wave solutions for any α ∈ N * are expected to exist and one can investigate their behavior analytically and numerically.However, the study of the time-dependent equation from an analytical point view seems much more involved.
From a numerical point of view, a stiff time integrator would be recommended for higher dimensions in order to overcome stability constraints.A straight forward, but computationally expensive approach would be an implicit scheme.More interesting would be Rosenbrock-type integrators based on a linearization of the equation after the spatial discretisation, and an eponential integrator for the Jacobian of the resulting system.This can be done efficiently by using so-called Leja points, see for instance [2] and references therein.
We leave this generalization to future work.

Appendix A. Formal derivation of the non-relativistic model
Consider the following nonlinear Dirac equation in one space dimension with κ 1 and κ 2 positive constants and α ∈ N * .Here Ψ = (ψ, ζ) is a 2-spinor that describes the quantum state of a nucleon of mass m, and σ 1 and σ 3 are the Pauli matrices given by In nuclear physics, the interesting regime is when the parameters κ 1 and κ 2 behave like m, whereas κ 1 − κ 2 stays bounded.More precisely, let κ 1 = θm and κ 1 − κ 2 = λ with θ and λ positive constants.As a consequence, the nonlinear Dirac equation ( 15) can be written as ( 16) Hence, by writing φ(t, x) = e imt ψ(t, x) and χ(t, x) = e imt ζ(t, x), we obtain As usual, in the non-relativistic regime, the lower spinor χ is of order 1/ √ m.Hence, we have to perform the following change of scale with a = 2λ θ , and F, G defined by Finally, denoting ε = 1 m the perturbative parameter, we obtain In particular, when ε = 0, we have which leads at least formally to the time-dependent quasilinear Schrödinger equation Appendix B. Existence of positive solutions to the stationary equation In this appendix, we prove Theorem 1 and we make rigorous the construction of ground states presented in Section 2.
As in [3], we write the stationary equation (3) as a system of first order ODEs for any strictly positive integer α.
The existence of solutions to (19) is an immediate consequence of the Cauchy-Lipschitz theorem.More precisely, we have the following lemma.
A key ingredient of the proof is to remark that system (19) is the Hamiltonian system associated with the energy As a consequence, to have a complete description of the dynamical system, it is enough to analyze the energy levels of (20), i.e. the curves in the (χ, ϕ)-plane defined by Γ c = {(χ, ϕ)|H(χ, ϕ) = c}.As a consequence of Remark 6, we have the following lemma.
Hence, let a > (α + 1)b.In this case we are able to prove the existence and uniqueness (modulo translations) of a positive solution ϕ to of (3) such that lim x→±∞ ϕ(x) = 0.
Remark 11.As shown in Section 2, a straightforward computation leads to the explicit formula for ϕ.In particular, for any fixed x 0 ∈ R, Finally, we conclude by proving the non-degeneracy of the solution ϕ.The linearized operator at our solution ϕ is defined by Our goal is to prove that in L 2 , ker L = span{ϕ }.

Figure 3 .
Figure 3. Difference of the numerical solution to the equation (5) for initial data being the ground state solution for a = 9 and b = 4.4, and the exact solution for t = 1.

Figure 5 .
Figure 5. On the left the L ∞ norm of the solution of Fig. 4, on the right the modulus of the Fourier coefficients of the solution at the final time.

Figure 7 .
Figure 7. On the left the L ∞ norm of the solution of Fig. 4, on the right the modulus of the Fourier coefficients of the solution at the final time.

Figure 8 .
Figure 8. L ∞ norms of the solution of (5) for the initial data φ(x, 0) = ϕ(x) ± exp(−x 2 ), on the left for the minus sign, on the right for the plus sign in the initial data.

Figure 12 .
Figure 12.On the left the L ∞ norm of the solution of Fig. 11, on the right the solution at the final time.