Atomic soliton transmission and induced collapse in scattering from a narrow barrier

We report systematic numerical simulations of the collision of a bright matter-wave soliton made of Bose-condensed alkali-metal atoms through a narrow potential barrier by using the three-dimensional Gross–Pitaevskii equation. In this way, we determine how the transmission coefficient depends on the soliton impact velocity and the barrier height. Quite remarkably, we also obtain the regions of parameters where there is the collapse of the bright soliton induced by the collision. We compare these three-dimensional results with the ones obtained by three different one-dimensional nonlinear Schrödinger equations. We find that a specifically modified nonpolynomial Schrödinger equation is able to accurately assess the transmission coefficient even in a region in which the usual nonpolynomial Schrödinger equation collapses. In particular, this simplified but very effective one-dimensional model takes into account the transverse width dynamics of the soliton with an ordinary differential equation coupled to the partial differential equation of the axial wave function of the Bose–Einstein condensate.

Localized soliton-like structures, known as bright matter-wave solitons, can be generated in Bose-Einstein Condensates (BEC) with attractive interatomic interactions.Since the first experimental realization of such structure about two decades ago 1 , dynamics of matter-wave bright solitons in an attractive BEC have been intensely studied, both at the quantum level and using the Gross-Pitaevskii equation (GPE).Current experimental capabilities offer an unprecedented opportunity to test many-body theories in an ultracold Bose gas, and matter-wave solitons are an excellent target for the predictions 2,3 .Moreover, several technological applications are based on the possibility of generating and manipulating such kind of coherent structures: some of them are interferometry 4 even beyond the quantum limit, and quantum-enhanced metrology 5,6 .The remarkable analogy of models based on the 3D-GPE with the equations of motion occurring in optics with Kerr media allowed to study common aspects on the same ground, such as the Hong-Ou-Mandel experiment 7 .
In a typical setup of an atomic interferometer using solitons, a matter-wave soliton is prepared in a quasi-1D trap, i.e. a confining potential made of a strong radial component and a weak or absent axial component.This setup experimentally allows the creation of cigar-shaped condensates in the case of repulsive interparticle interaction and in the noninteracting case.By using attractive interparticle interactions that are obtainable for example by using Feshbach resonances, one can generate a matter-wave soliton.In the latter case, solitons loaded into quasi-1D traps have typical axial widths that are comparable to the radial potential characteristic length.To achieve interference, the soliton is set into motion by phase imprinting, and it collides with a narrow potential barrier set by a narrow laser beam, acting analogously to an optical beam splitter, designed to be able to split the number of atoms into two even solitonic packets.The resulting two solitons are then recombined in a later stage, through the same barrier.After the first splitting, split solitons may achieve a differential phase shift, thus allowing the observation of interference in the recombined packet.
From the theoretical point of view, the study of quantum matter-wave solitons was carried out mostly in 1D, where the many-body wave function is well known to have an exact solution by Bethe ansatz 8 .This procedure relies on having very strong radial confinement, and it is not sensible to the 3D dynamics, thus losing details about the transverse degrees of freedom that are especially important near the point of GPE collapse in which the reduction of the size of the condensate brings it to a regime in which other interaction effects start to be non-negligible, like three-body interactions causing depletion from the trap 9 .
The stability of solitons in quasi-1D harmonic traps is highly nontrivial, as the cubic GPE has a critical dimension equal to 2. Instability is in the form of a collapse, also known in the mathematical literature as nonlinear blow-up 10 .In one dimension, the collapse is prevented by the Vlasov-Petrishev-Talanov theorem 10 .GPE collapse due to an arbitrary attractive interaction potential can be triggered in this context by loading into the trap a suitably high number of particles 11 .Moreover, the radial anisotropy of the trap can play a role in the critical number of particles for collapse 12,13 .It is fundamental to remark that, even with a purely 1D model, by adding a barrier-like external potential to the 1D-GPE the problem becomes non-integrable, and requires approximate methods to be tackled.
The 3D-GPE dynamics represents a useful tool not only as an approximation of the full quantum dynamics but also as providing signatures of soliton entanglement across the barrier 14 .In fact, the discontinuity in the reflection coefficient indicates the possibility of creating "Schrödinger cat" states, exploring quantum entanglement phenomena 14 .The 3D-GPE model predicts a peculiar behavior of the transmission coefficient with the barrier: at low values of the velocity and the barrier height it is a discontinuous function of the parameters.This aspect has been investigated 14 and is frequently referred to as the particle behavior of the impinging soliton.
Various dimensional reduction schemes for the 3D-GPE have been proposed, 2,[15][16][17] .In this work, we compare three schemes of dimensional reduction with full 3D simulations.1D effective equations have a great computational advantage in the description of the dynamics and are routinely used in studies of atomic interferometers.The simplest one is the 1D-GPE, obtained by imposing a fixed transverse wave function as the lowest energy eigenstate of the transverse harmonic potential.An improved model is called nonpolynomial Schrödinger equation (NPSE) 17 , which is based on assuming the transverse width of the trial as a variational parameter and obtaining the equation of motion as Euler-Lagrange (EL) equations.Furthermore, in the original NPSE formulation, derivatives of the transverse width parameter present in the Lagrangian are neglected, and the corresponding equation is algebraic.So we also consider the non-approximated version of the NPSE, which we call NPSE+ for brevity.Previous work 18 highlighted the behavior of the transmission coefficient and the barrier-induced collapse with the 3D-GPE and the NPSE, studying regions in the barrier height vs. number of particles plane.We investigate the differences in the collision process varying the velocity, focusing on the transmission coefficient and the onset of barrier-induced collapse.

Gross-Pitaevskii equations
The model is based on the Hartree approximation for bosons 8 , using which it is possible to derive the following Lagrangian, called Gross-Pitaevskii Lagrangian, for the field ψ(r,t), representing the wave function of the Hartree product state, with all the particles in the same single-particle quantum state, where U is the external potential, N is the number of particles, and g is the contact potential, which can be linked to the s-wave scattering length a s with the expression The associated EL equation is the 3D-GPE: Standard dimensional reduction of the 3D-GPE in a tight transverse harmonic potential relies on the assumption that the transverse degree of freedom of the wave function is frozen to the ground state of the harmonic potential, as we will briefly review now.Let the external potential be written as where ω ⊥ is the (isotropic) strength of the potential, and V is the axial part of the potential.The role of anisotropy on the transverse potential was studied in 12,13 .The strength ω ⊥ naturally sets a characteristic length scale l ⊥ = h/(mω ⊥ ).Let us assume that the wave function is composed of a constant Gaussian transverse part φ , which is the ground state of the transverse harmonic potential, and a time-varying axial component f as where  This is physically justified when the interaction energy is much smaller than the energy difference between the first excited state and the ground state of the transverse potential.Inserting this ansatz into Eq.( 3), and integrating along the transverse coordinates, one obtains the corresponding wave equation, called the 1D-GPE: where we defined g 1D = g(N − 1)/(2πl ⊥ ).The expression of g 1D can be used to define a normalized nonlinear parameter γ = (N − 1)|a s |/l ⊥ .Using Eq. ( 7), for g 1D < 0 the ground states of the axial problem are constituted by stable solitons, and no collapse is expected for any interaction strength.Instead, for the 3D-GPE case, there exists a critical nonlinear parameter for the existence of stable solitons.The value, above which static wave function collapse is expected, is about γ c ≈ 0.67.

Variational ansatz and NPSE
As shown in 17 , a better approximation is to consider the separation of the total wave function in a transverse Gaussian component with non-constant transverse width σ (x,t) and to find the equation of motion using a variational principle.The resulting equation is the NPSE.The function φ in the ansatz Eq. ( 5) is substituted by a more general where σ is a function to be determined as a variational parameter.In 17 , the calculations were done assuming that derivatives of σ are negligible.By keeping these derivatives terms, it is possible to write the following effective 1D Lagrangian (detailed calculations are given in the Methods section) The corresponding EL equations for f and σ provide the solution to the variational problem and are obtained as shown in 19 : 3/9 We will refer to the above coupled equations as NPSE+.By neglecting the derivative of σ in Eq. ( 9), one obtains another effective 1D Lagrangian 17 , whose EL equations, called NPSE, correspond to We remark that the NPSE+, as opposed to the NPSE, respects the variational principle, so the corresponding ground state energy is bound to be greater or equal to the true ground state energy of the 3D-GPE.
We will use the 3D-GPE as a reference equation, and compare the predictions of the 1D-GPE, the NPSE, and the NPSE+.The axial densities of the ground state solutions are shown in Fig. 1, where we have set the nonlinear parameter to a very high value γ = 0.65, near the 3D-GPE static collapse value γ c .We notice that the 1D-GPE fails to represent accurately the axial wave function, NPSE+ and NPSE have similar accuracy.NPSE+ has the additional cost of the computation of the solution of the transverse width differential equation coupled to the axial wave function partial differential equation.

Generalization of bound on splitting energy
The soliton splitting event can be verified only for specific ranges of the transmission coefficient 20 , depending on the initial soliton velocity.These values can be computed by analyzing energy conservation in the splitting event.The interplay of kinetic energy and internal energy of the solitons during the scattering event has been discussed in [20][21][22] , by using Lieb-Liniger 8 energies E G pertaining to the soliton internal degrees of freedom in the total Hamiltonian.Imposing energy conservation, the kinetic energy of the initial soliton must satisfy where N is the number of atoms in the initial soliton, and n is the one in the transmitted soliton.This is the condition that must be satisfied for having a splitting event of transmission coefficient T = n/N.In our case, the internal energy of the soliton can be computed either numerically or analytically.The chemical potential of the stationary solution for the NPSE can be obtained 17 from the implicit relation and selecting only the stable branch of the solutions, i.e. the one satisfying the Vakhitov-Kolokolov criterion ∂ ∂ n µ < 0. For the other equations used in this work, it is possible to obtain numerically the value of µ(n) from stationary state solutions.Using the values of the chemical potential, we are able to write the ground state energy of the nonlinear wave equation corresponding to an N-particle soliton as it is possible to obtain different ranges of transmission coefficients that are accessible for a given value of the initial kinetic energy.

Soliton solutions
We review some properties of the solitonic ground state of the equations, comparing them.We study the highly nonlinear regime, in which γ = 0.65.Soliton solutions in this case are stable for the 3D-GPE and the NPSE and NPSE+ for γ < γ c 10, 16 .Soliton solutions are shown for all the equations in Fig. 1 and their transverse width in Fig. 3.The simulations show a better agreement of the NPSE+ equation with respect to the 3D-GPE in the transverse width.The computation of the NPSE+ transverse width is obtained by iteratively solving Eq. ( 10) and then Eq. ( 11).In the solution of the ordinary differential equation Eq. ( 11) we apply Dirichlet boundary conditions corresponding to the vanishing of the axial wave function at infinity.The computation of the 3D-GPE transverse width is done by a least square fit on the radial distribution of the wave function, namely defined as   where The chemical potential of the solitons is shown in Fig. 2, and the corresponding energy is shown in Fig. 4. We remark that, being the NPSE only an approximation of the true variational solution, its chemical potential is allowed to be less than the 3D-GPE chemical potential, thus becoming less than the bound set by the variational principle.As expected, the difference becomes more pronounced at high nonlinearities, implying a more localized wave function, where the terms proportional to the derivatives of σ in Eq. ( 10) and ( 11) become more relevant.

Scattering from a narrow barrier
We are interested in computing the transmission coefficient for various velocities and barriers.We assume energy in units of hω ⊥ , time in units of ω −1 ⊥ , and length in units of l ⊥ .The barrier is Gaussian, centered in x = 0, and it is parametrized by the peak value parameter b, The width w is fixed to w = 0.5.In our simulations, the velocity ranges in v ∈ [0.1, 1.0], and the barrier in b ∈ [0.0, 1.0].In particular, by setting a sufficiently high γ, for example near to γ c , namely γ = 0.65, we analyze the onset of barrier-induced collapse, happening for high soliton velocity and high barrier height, as shown in Fig. 5. Our results show that the vanishing of the transverse width predicted by the NPSE, suggesting a barrier-induced collapse, is a very weak indicator of an actual collapse.Instead, the NPSE+ collapsing region is not due to a vanishing of the transverse width, but to a sudden concentration of the axial density in smaller and smaller regions, like in the 3D-GPE case.In fact, we have set the numerical threshold of the collapse to the detection of a single probability per site greater than 0.3.The abrupt change in the local maximum density we observe between stable solutions and collapsing ones justifies the validity of this criterion.In the region of parameters we investigated, the NPSE+ is not collapsing, as reported in the transmission functions at constant velocity in Fig. 6, so it is ineffective in predicting the barrier-induced collapse present in the 3D-GPE.Results obtained by the familiar 1D-GPE, are collapse-free.The comparison of the transmission coefficient versus the barrier height with fixed velocity reported in Fig. 6 shows that, quite remarkably, the NPSE+ can describe accurately the transmission coefficient in the non-collapsing region.

Discussion
In this article, we have presented a numerical study of the collision of a bright matter-wave soliton with a narrow potential barrier using the three-dimensional Gross-Pitaevskii and three dimensionally reduced versions of it.We investigated how the choice of dimensional reduction impacts the description of some features of the process, namely the transmission coefficient and the onset of barrier-induced collapse, also using the familiar one-dimensional Gross-Pitaevskii.We first reviewed the  ground state properties given by all the schemes, highlighting the role of the variational transverse width.Then we compared the scattering properties: our results show that by using the NPSE in a regime of high barrier height and high velocity it fails to describe the 3D dynamics due to the the vanishing of the transverse width of the solution.In such cases the collapse phenomena in the 3D solutions are absent, and the NPSE is not capturing the correct dynamics.
Our main result is that by adopting a slight modification of the NPSE, by using the true variational solution with the NPSE ansatz called NPSE+, we can predict the transmission factor and the dynamics of the transverse width more accurately, even though the collapse phenomenon is not captured by this effective equation.We believe the present work is a valuable contribution to the field of matter-wave soliton interferometry and quantum measurement, as results can be used to predict the dynamics of experimentally accessible scenarios.For example, the NPSE+ can be used for modeling interferometric experiments in highly nonlinear regimes where the determination of the transverse width is important.In the setting of a quasi 1D harmonic trap corresponding to a transverse frequency ω ⊥ /2π = 254Hz, loaded with about 28 × 10 3 atoms of 7 Li, analogously to a past experiment 23 , the nonlinear regime we have studied is achieved for an s-wave scattering length of a s ≈ −5.52 × 10 −11 m.

Numerical methods
The time-marching scheme we use for all the simulations is the split-step Fourier method (SSFM), adopting Strang splitting of the nonlinear and linear part of the evolution operator.The SSFM is well-known to be accurate to the second order in time and to every order in space, thus being highly efficient in the spatial discretization 24,25 .The drawback of the method -or the feature, depending on the point of view -is to natively implement periodic boundary conditions.The implementation of absorption boundaries in the context of this method is still possible but not straightforward 26,27 .We assume the field to be localized away from the boundary in order to neglect this problem.In our setup, we use a unit of energy of hω ⊥ , a unit of time of ω −1 ⊥ and a unit of length of l ⊥ , constituting the natural units for the (isotropic) harmonic confinement.In these units, we consider for the 1D simulations a total length of L = 40, with a grid of N = 512 points.In the 3D simulations, we use a grid of (N x , N y , N z ) = (512, 40, 40) points, with total lengths of (L x , L y , L z ) = (40, 10, 10).The time step in both setups is chosen to be h t = 0.01.These parameters have been proven to give a total truncation error in the L ∞ norm of the order of 10 −4 in 1D solitonic ground state solutions and 3D linear problems with anisotropic three-dimensional harmonic trap.
The ground state solutions are computed using an imaginary-time propagation method.We point out that some modifications of this method are available under the name of normalized gradient-flow methods 28 .

Derivation of the NPSE+
Following 16 , we write the 3D Lagrangian We are interested in integrating along the transverse coordinates without neglecting the terms proportional to ∂ ∂ x σ and ∂ 2 ∂ x 2 σ .The novel terms arise from ih ∂ ∂t ( f φ ) and h 2m ∇ 2 ( f φ ).By separating the derivatives, we have Integrating the term proportional to ∂ ∂t σ gives 0, as one may realize by looking at the symmetry of its prefactor.However, the term proportional to ∂ 2 ∂ x 2 φ gives a non-null contribution to the 1D Lagrangian.The integration gives By considering the Euler-Lagrange equations, we recover Eq. ( 10) and Eq. ( 11).

) 2/ 9 ( a )
Soliton solutions compared.(b) Zoom applied to the top part of the soliton solutions.

Figure 1 .
Figure 1.Comparison of the axial ground state wave function f .The space coordinate x is in units of l ⊥ = h/(mω ⊥ ), the characteristic length of transverse harmonic confinement of frequency ω ⊥ .The nonlinear parameter is set to γ = (N − 1)|a s |/l ⊥ = 0.65.The three-dimensional Gross-Pitaevskii equation is the red solid line, the nonpolynomial Schrödinger equation without the corrections is the green dotted line, the one with the corrections is the green dashed line.The one-dimensional Gross-Pitaevskii equation is the dash-dot grey line.

Figure 2 .
Figure 2. Chemical potential µ as a function of the nonlinear parameter γ. µ is in units of hω ⊥ .The three-dimensional Gross-Pitaevskii equation is the red solid line, the nonpolynomial Schrödinger equation without the corrections is the green dotted line, and the one with the corrections is the green dashed line.The one-dimensional Gross-Pitaevskii equation is the dash-dot grey line.

Figure 3 .
Figure 3.Comparison of the transverse width parameter σ .Both σ and x are in units of l ⊥ = h/(mω ⊥ ).The colors are as in Fig. 2.

Figure 4 . 8 TFigure 5 .
Figure 4. Energy E as a function of the nonlinear parameter γ, assuming a trap geometry such that l ⊥ /|a s | = 2 × 10 4 , that corresponds to a critical particle number N c ≈ 13400.E is in units of hω ⊥ .The three-dimensional Gross-Pitaevskii equation is the red solid line, the nonpolynomial Schrödinger equation without the corrections is the green dotted line, and the one with the corrections is the green dashed line.The one-dimensional Gross-Pitaevskii equation is the dash-dot grey line.

Figure 6 .
Figure 6.Transmission coefficient T versus barrier height b with fixed velocity v. v that is in units of hω ⊥ /m, and b is in units of hω ⊥ .The three-dimensional Gross-Pitaevskii equation is the red solid line, the nonpolynomial Schrödinger equation without the corrections is the green dotted line, and the one with the corrections is the green dashed line.The one-dimensional Gross-Pitaevskii equation is the dash-dot grey line.