Effective self-similar expansion of a Bose-Einstein condensate: free space vs confined geometries

We compare the exact evolution of an expanding three-dimensional Bose-Einstein condensate with the evolution obtained from the effective scaling approach introduced in Ref. [1]. This approach, which consists in looking for self-similar solutions to be satisfied on average, is tested here in different geometries and configurations. We find that, in case of almost isotropic traps, the effective scaling reproduces with high accuracy the exact evolution dictated by the Gross-Pitaevskii equation for arbitrary values of the interactions, in agreement with the proof-of-concept of Ref. [2]. Conversely, it is shown that the hypothesis of universal self-similarity breaks down in case of strong anisotropies and trapped geometries.


I. INTRODUCTION
Self-similarity is a remarkable property that plays a key role for describing the dynamics of several ultracold atomic systems.It characterizes the free expansion and the collective excitations of Bose-Einstein condensates both in the noninteracting and the hydrodynamic regimes [3][4][5][6][7][8][9][10], the expansion of a one-dimensional Bose gas in the mean-field Thomas-Fermi regime and in the Tonks-Girardeau regime [11], of a superfluid Fermi gas [9,12,13], and of a thermal cloud [14].Beside its conceptual appeal -the fact that a system evolves maintaining the same exact shape throughout the whole dynamics -, self-similarity also implies a strong simplification in the numerical treatment, allowing for a dramatic reduction of the mathematical complexity of the equation governing the system: instead of having to deal with partial differential equations, one can obtain the evolution of the system by solving a set of ordinary differential equations for the time-evolution of the scaling parameters that characterize the rescaling of the coordinates.
In view of this numerical simplification, approximate self-similar behaviors have been conjectured by several authors for dealing with problems that otherwise would be impossible to tackle.Effective scaling approaches have been employed as approximate solutions for describing the collective excitations of a trapped Bose gas [1], the expansion of an interacting Fermi gas [15], and of quantum degenerate Bose-Fermi mixtures [16].This approach, originally proposed in Ref. [1], consists in using a self-similar ansatz for the evolution of the system in the hydrodynamic regime, then requiring the hydrodynamic equations to be satisfied on average, after integration over the spatial coordinates.s Recently, this effective approach was tested for the case of a freely expanding quasi-one-dimensional Bose-Einstein condensate (BEC), by comparing the approximate solution with the exact evolution of the system as obtained from the solution of the Gross-Pitaevskii (GP) equation.Remarkably, in this case it was found that the effective scaling approach is indeed very accurate in reproducing the exact evolution, for arbitrary values of the interactions [2].
In the present paper we extend this analysis to the case of a three-dimensional (3D) BEC [17], for which we consider the expansion in free space and in trapped geometries (namely, in cylindrical and planar waveguides).In the former case (free space), we find that the effective scaling approach for the expansion of a sphericallysymmetric condensate is rather accurate even for intermediate values of the interaction -in between the noninteracting and hydrodynamic limits where the scaling is exact -, similarly to what found for the quasi-1D case [2].Deviations from this behavior are observed instead for prolate and especially oblate condensates, signalling that the scaling approach becomes less accurate when the expansion along certain directions is faster than along the others, causing local variations of the density that do not conform to the hypothesis of self-similarity.Likewise, self-similarity may also break down in the presence of a residual trapping.In fact, the density profile along the trapped directions is determined by the interplay of the confining potential and the nonlinear term, and it therefore gets modified as the system expands.In this scenario, the scaling approach can be safely employed only in the noninteracting limit (for obvious reasons) and in the hydrodynamic limit.In the latter case, a necessary condition is also that the evolution time should be sufficiently short such that the initial Thomas-Fermi (TF) profile along the trapped directions is preserved.This is conceptually different from the case of the evolution in free space, where e.g.once a TF profile is fixed by the initial conditions it is then maintained during the whole evolution, regardless of the variation in the local value of the nonlinear term.
The paper is organized as follows.In Sec.II we review the general procedure for constructing the effective scaling approach within the hydrodynamic formulation of the GP equation, as well as the form of the GP equation in the rescaled coordinate system.Then, in Sec.III we apply this approach to the case of a freelyexpanding condensate in 3D, considering first the case of a spherically-symmetric BEC, and then the cases of oblate and prolate geometries.In Sec.IV, as an example of confined geometries, we examine the expansion in a waveguide, first within the GP approach, and then by using the 1D non-polynomial Schrödinger equation (NPSE) of Ref. [18].There we also comment on the case of a planar waveguide.Final considerations are drawn in the conclusions.

II. SCALING APPROACH
Let us consider a Bose-Einstein condensate described by a wave function ψ(x, t) that evolves according to the Gross-Pitaevskii (GP) equation [5] i with m being the particle mass, g = 4π 2 a/m the interaction strength (with a the s-wave scattering length), and V (x, t) a generic trapping potential, that may depend on time.In particular, we shall consider the case of a time-dependent harmonic potential of the form and we shall focus on the expansion dynamics of the condensate following the release of the confinement along some direction.Eq. ( 1) can be transformed into a system of two coupled hydrodynamic equations by means of a Madelung transformation, ψ(x, t) ≡ n(x, t) exp {iS(x, t)} [5], where v ≡ ( /m)∇S, and represents the so-called quantum pressure term.Exact Scaling.The scaling approach consists in looking for solutions characterized by a density profile that evolves self-similarly as with n 0 (x) being fixed by the initial conditions, and with all the time dependence being contained within the scaling parameters λ i (t) (here the term j λ −1 j is just a volume normalization factor, with i, j = 1, 2, 3 labelling the spatial directions).In addition, the continuity equation (3) yields the following scaling of the velocity field: By introducing the ansatz (6) into the hydrodynamic equation ( 4), one gets m λi In the noninteracting regime (g = 0) and in the so-called hydrodynamic regime (negligible P ) the above expression can be factorized in two terms, one depending on the spatial coordinates and the other on the time coordinate alone, so that the scaling ansatz (6) represents an exact solution.
Effective Scaling.Contrarily, when such a factorization is not possible, one may follow a different approach by requiring the self-similarity to be satisfied on average, as discussed in Refs.[1,2,15,16].First, it is convenient to transform the expression in Eq. ( 7) by integrating over the i-th coordinate [19], where q(x, t) ≡ P (x, t) + V (x, t) + gn 0 (x i /λ i (t))/ j λ j .Then, we multiply Eq. ( 8) by n 0 (x i /λ i (t)) and we get rid of the coordinate dependence by integrating over the volume.This sort of averaging procedure -that constitutes the essence of the effective scaling approach [1,2] -yields an equation for the scaling parameters depending only on the time variable.Specific expressions for each case considered in this paper will be discussed in the following sections.
Scaled Gross-Pitaevskii equation.The scope of this work is to test the accuracy of the effective scaling approach in reproducing the exact evolution by quantifying the deviation of the scaling dynamics from the exact solution of the GP equation in (1).To this end, it is convenient to rewrite the GP equation in terms of the rescaled wave function, defined as [3] ψ(x, t) = 1 j λ j (t) with φ(x, 0) = n 0 (x).Then, plugging the above expression for ψ(x, t) in the Gross-Pitaevskii equation ( 1) yields where we have defined ∇ j = ∂/∂(x j /λ j ).Notice that the partial derivative ∂ t only operates on the explicit dependence on time in φ(•, t), which accounts only for the deviation from the self-similar solution.The latter can be measured by defining the following fidelity [2] that can take values in the range F ∈ [0, 1], the upper bound (F = 1) corresponding to exact self-similar solutions.This quantity provides a quantitative estimate of the accuracy of the effective approach.
In the following, it is also useful to introduce dimensionless units, which can be conveniently defined in terms of natural oscillator units of a suitable frequency ω 0 , namely τ ≡ ω 0 t, with the spatial coordinates being expressed in units of a ho ≡ /(mω 0 ), and so on.As for the interaction strength, we define g ≡ 4πN a/a ho , by including also the number of atoms N in the definition.From now on, all the quantities will be considered as dimensionless, unless otherwise stated.

III. EXPANSION IN FREE SPACE A. A spherical condensate
As a first example, let us consider the case of BEC released from a spherically-symmetric harmonic trap, with ω(0) = 1 and ω(τ ) = 0 for τ > 0.Here the scaling ansatz takes the following form (see Appendix A 1) where the scaling parameter λ(τ ) evolves according to and the parameters A(g) and B(g) -which are fixed by the initial conditions -fulfill the sum rule A(g)+B(g) = 1 [2], see Fig. 1.Manifestly, the above formula in Eq. ( 14) interpolates between the non-interacting (λ −3 term) and TF (λ −4 term) limits [5].
The behavior of the scaling parameter λ(τ ) as a function of time (for different values of g), and its final value has a slight minimum at the crossover region between the noninteracting and hydrodynamic regimes.The first behavior is directly connected to the power law in the scaling equation ( 14): in the TF limit the RHS behaves like ∝ λ −4 [10], and this implies a slower growth of λ (> 1) as compared to the ∝ λ −3 behavior of the noninteracting limit.This reflects the fact that as g is increased, the density distribution gets wider and this narrows the momentum distribution, thus reducing the contribution of kinetic energy to the condensate expansion.As for the non-monotonic behavior at intermediate values of g, we notice that the minimum of λ f corresponds to the min- imum of the fidelity, shown in Fig. 3.This figure shows that the accuracy of the effective scaling can be very high, with a fidelity always above 98%, for any value of the interaction strength (τ ≤ 10).This approach seems therefore quite robust, similarly to the 1D case discussed in Ref. [2].
Its reliability can also be appreciated from the behavior of the density profiles, which are shown in Fig. 4 for different values of g.Only small deviations can be seen, especially around the minimum of fidelity (at g ≡ g * ).
It is also interesting to notice that for g < g * the scaling ansatz produces a slightly larger central value for the density, whereas the opposite happens for g > g * .FIG. 4. Comparison of the radial density profiles of a spherical condensate obtained from the numerical solution of the full GP equation (solid) and from the scaling argument (dashed) for different interaction strengths, at the final evolution time τ f = 10 (the corresponding fidelity is marked by the black circles in Fig. 3).Notice the different range of horizontal axis in the two plots.

B. Pancake-and cigar-shaped condensates
Here we consider the case of a condensate that expands from a cylindrically-symmetric harmonic trap where ω ⊥,z (τ ) ≡ 0 for τ > 0. As unit frequency here we choose the geometric average of the trapping frequencies in the three spatial directions and we define the parameter α = ω ⊥ (0)/ω z (0) characterizing the trap anisotropy.In this case, the self-similar ansatz for the density takes the following form The corresponding scaling equations are (see Appendix A 2 for the derivation and the definition of the various constants) where the constants A i and B i fulfill the following sum rules Behavior of the constants Ai and Bi (i = 1, 2, 3) in Eqs. ( 18) and ( 19) as a function of the interaction strength g.Here α = 10.Their behavior as a function of the interaction strength g, see Fig. 5, guarantees that the known results [5] corresponding to the exact solutions in the non-interacting (g ≪ 1) and Thomas Fermi limits (g ≫ 1) are recovered.
In the following, we shall compare the case of prolate and oblate geometries, with α = 10 and α = 0.1 respectively, with the previously considered spherical case (α = 1).The corresponding values of the fidelity after an expansion time τ f = 10 are shown in Fig. 6.This figure shows that, in the presence of a trap anisotropy, the assumption of a self-similar evolution is less effective for intermediate values of the interactions, especially for the case of oblate geometries (the case with α = 0.1 in the figure).A qualitative argument for explaining this behavior is the following.When a pancake-shaped condensate is released from the trap, the initial expansion is characterized by a fast dilatation along the axial direction, whereas the radial dynamics is much slower.The former causes a fast local variation of the density, which drives the radial expansion out of the self-similar regime.This happens due to the fact that the radial dynamics soon decouples from that of the mean field term, contrarily to what occurs e.g. in the expansion of a spherical condensate.This is confirmed by Fig. 7(a), where we plot the behavior as a function of time of the condensate radial and axial widths, σ ⊥ (τ ) ≡ r 2 ⊥ τ and σ z (τ ) ≡ z 2 τ , respectively.This figure shows that the value of σ z (τ )/σ z (0) remains very close to one in the rescaled coordinate system we are using, indicating that the evolution of the condensate along the axial direction is well reproduced by the hypothesis of self-similarity.However, the expansion in the transverse plane is faster than what is dictated by the self-similar approach, signaling that the latter substantially underestimates the actual value of σ ⊥ (τ )/σ ⊥ (0).
In the opposite case of a cigar-shaped condensate, the fidelity substantially improves (though not at the level of the spherically-symmetric case), see the curve for α = 10 in Fig. 6.The reason for this is that here the dynamics is driven by two of the three spatial directions, so that the breaking of the self-similarity mostly affects one direc- tion only, namely the axial one.Indeed, Fig. 7(b) shows that the scaling for the transverse width is almost exact, whereas the axial width increases faster that what is predicted by the self-similar expansion.This confirms that the self-similar approach underestimates the expansion along the "slow" direction.
The above behavior has a direct consequence on the value of the aspect ratio AR(τ ) ≡ σ ⊥ (τ )/σ z (τ ), a quantity that is typically used to characterize the condensate expansion [5,15,16].Consistently with the scaling approach, here it is convenient to consider the ratio AR(τ )/AR(0).In case of an exact self-similar expansion, the former quantity is expected to be constant and equal to one in the rescaled coordinate system we are using, AR(τ )/AR(0) ≡ 1 for any τ .The deviations from the ideal behavior are shown in Fig. 8.We see that in the spherical case the self-similar approach provides very accurate results, whereas in case of oblate and prolate condensates significant deviations are observed intermediate values of the interaction (up to 50%), as one may expect from what was discussed previously.
To conclude this section, we remark that the scaling approach can be still very useful from the numerical point of view, even when it is not able to reproduce the exact evolution.Indeed, it allows us to replace the usual GP equation with the one in Eq. (10), which evolves the system in a rescaled coordinate system, thus greatly reducing the size of the computational grid in case of a free expansion.

IV. CONFINED GEOMETRY: EXPANSION IN A WAVEGUIDE
Let us now turn to the case of a condensate initially confined by the same cylindrically-symmetric potential in Eq. (15), that is now allowed to expand in the presence of the radial confinement, following the switch-off of the axial potential [ω z (τ ) = 0 for τ > 0].This case represents a generalization of the quasi-1D case considered in Ref. [2], the difference being that here the condensate profile along the radial direction is not restricted to that of the harmonic oscillator ground state -several transverse excited states can be occupied, depending on the values of g and α.The 1D mean field limit of Ref. [2] is recovered for g ≪ 4π √ α, see Ref. [20].

A. Scaling of the GP equation
Considering that the expansion takes place along the axial direction, in this section we express frequencies in units of the axial frequency, namely ω z (0) = 1, ω ⊥ (0) ≡ α [21].In the following we shall use α = 10, as an example.In the present case, the self-similar ansatz for the density takes the same form as in Eq. ( 17), with Eq. ( 18) being replaced by where the constants A i and B i are those previously defined (modulo a straightforward rescaling, owing to the different choice of ω 0 ).The behavior of the scaling parameters λ z and λ ⊥ as a function of time (for different values of g) and as a function of g (after an expansion time τ = 10) are shown in Fig. 9, in the left and right panels respectively.Their behavior corresponds to a contraction of the radial size, and to a dilatation in the axial direction, as one would expect for a condensate expanding along the waveguide.Similarly, the monotonic behavior of λ z (τ f ) as a function of g -which again interpolates between the the two limiting behaviors (non-interacting, up to g ≈ 10 −2 , and TF for g 10 3 ) -, agrees with the naive expectations.Instead, the local minimum displayed in the behavior of λ ⊥ (τ f ) as a function of g, see Fig. 9d, signals a failure of the self-similarity assumption, since the radial size is expected to decrease monotonically by increasing g.Indeed, by looking at the density cuts along the axial and radial directions in Fig. 10, it is manifestly evident that for intermediate values of the interactions the effective scaling approach breaks down.In these figures, the dashed profile corresponds to the initial state -which is equivalent to the self-similar solution in the rescaled coordinate system -, and the solid line is the solution of the scaled GP equation at τ f = 10.It is evident that for intermediate values of g, where the radial scaling parameter becomes too small (see Fig. 9), the actual density profile is characterized by a strong contraction along the axial direction (and by a corresponding radial broadening) with respect to what a perfect scaling would predict (dashed lines in the figure).As anticipated in the introduction, this behavior is due to the FIG.10.Expansion in a waveguide: comparison of the density profiles as given by the numerical solution of the GP equation (solid) and by the scaling approach (dashed) for different interaction strengths and an evolution time of τ = 10.Density cuts along the z-axis at r ⊥ = 0 (left), and along the r ⊥ -axis at z = 0 (right).For intermediate values of g, the scaling solution overestimates greatly the central density and underestimates the radial width.The fidelity corresponding to each value of g considered here is marked by the black circles in Fig. 11.
fact that self-similarity is explicitly violated in the presence of a residual trapping.In fact, the density profile along the trapped directions is determined by the interplay between the confining potential and the nonlinear term, and it therefore gets modified as the condensate expands: it cannot be maintained during the expansion along the waveguide [22].In this respect, also in the case of g > 10 4 -for which the scaling approach appears reasonably accurate, see Figs. 9, 10 -, we expect this approach to eventually break down in the asymptotic regime of very low densities [23].
Indeed, the above behavior can be inferred by looking at the fidelity, which is shown in Fig. 11 as a function of g, at different evolution times.As a reference, here we also show the value of the fidelity at τ = τ f for the case of a shallower radial confinement, namely α = 1.Then, by comparing the two cases at τ = τ f = 10, for α = 1 and 10 (the black and yellow continuous lines, respectively), it is clear that in the former case the system is able to evolve self-similarly even for large values of g (in the time range considered here), whereas this is not the case for α = 10.The reason behind this is that for a lower initial local density, the system takes longer to exit the self-similar regime (for values of g for which the scaling is exact in the quasi-1D limit).Consequently, we expect the self-similarity to also break down eventually for α = 1, though at larger expansion times than those in the case of α = 10.
Therefore, we can conclude that the effective scaling approach to the GP equation for the waveguide expansion works only if one of the following two conditions is satisfied: i) short evolution times, when the radial profile has not yet deviated from the initial one (that is, the local density has not decreased too much), or ii) the system is in the 1D mean field limit [2,20], g ≪ 4π √ α, where the density profile is always characterized by a Gaussian transverse profile (for arbitrary times).

B. Effective 1D model (NPSE)
Here we consider an alternative approach, by means of the 1D non-polynomial Schrödinger equation (NPSE) of Ref. [18], that provides an effective 1D description of the condensate evolution in the presence of a transverse confinement.This approach consists in factorizing the condensate wave function into a radial and an axial component, ψ(r ⊥ , z) ≡ φ(r ⊥ )f (z), where φ(r ⊥ ) is a Gaussian of width σ(z, τ ) = 1 + 2ã|f (z, τ )| 2 , and f satisfies the following non-polynomial Schrödinger equation with all the quantities being dimensionless [ã ≡ aN/ /(mω z )], and with both φ and f normalized to unity.Remarkably, besides the obvious advantage of having to deal with just a one-dimensional equation, this approach is also appealing because the radial profile is effectively described by a Gaussian wave function, so that the self-similarity is not explicitly violated.However, there is a price to pay, given that the non-polynomial character complicates the derivation of the scaling equations, as it will be clear in the following.
By proceeding as in the previous cases, we first transform Eq. ( 23) into the corresponding hydrodynamic form (by writing f = √ ne iS ), and then we make use of the  usual scaling ansatz for the density The corresponding scaling for the wave function reads , t e i 2 λ(t) with ϕ(z, 0) = n 0 (z).This procedure yields where q(τ ) ≡ P (0, τ ) + α 1 + 3ãn 0 (0)/λ 1 + 2ãn 0 (0)/λ , with ξ = z/λ and ∂ z = λ −1 ∂ ξ .Then, the next step would be to multiply Eq. ( 26) by n 0 (z/λ(t)) and to get rid of the coordinate dependence by integrating over the volume, as we did for Eq. ( 8).However, the presence of the density n(z, τ ) in the denominator of the nonlinear term does not permit to obtain an equation with an explicit dependence on λ(τ ).Then, in order to proceed as in the previous section, we have to make further approximations, taking either the strongly-interacting (or high-density) limit or the weakly-interacting (or low-density) limit, as discussed in Appendix B.
Weakly-interacting limit.For 2ãn 0 ≪ 1, the equation for the scaling parameter takes the form with A + B = 1 (see Appendix B).It is important to remark that this sum rule is obtained by means of the corresponding truncated expression of the NPSE [obtained by replacing the nonlinear term with the expansion in Eq. (B2)] and not from the full NPSE in Eq. ( 23).In the following, the parameters A and B will be calculated from the full equation, so that the regime of validity of the weakly-interacting approximation can be directly controlled by looking at the value of the sum A + B. The latter is plotted as a function of g in Fig. 12a, showing that the the sum rule is satisfied up to ã ≃ 10 −1 .In principle this result can be easily improved -namely, the regime of validity extended to higher values of g -by considering higher orders in the expansion (B2).
In the present case, the approximate self-similar solution is obtained from Eq. (25), along with the solution of Eq. ( 29).The corresponding fidelity in reproducing the exact solutions of the NPSE ( 23) is reported in Fig. 12b.Different density profiles are shown in Fig. 12c.Remarkably, within the validity of the present approximation, ã 10 −1 , the fidelity stays above the 99% threshold (in the time range considered here, τ ≤ 10).Therefore, one may expect that a similar level of accuracy can be extended to stronger interactions by considering higher orders in the expansion (B2), as discussed above.
Strongly-interacting regime.In the opposite limit, 2ãn 0 ≫ 1 we have with A 1 + A 2 + A 3 = 1 (see Appendix B).Similarly to the previous case, this sum rule is obtained by means of the truncated NPSE, here with the nonlinear term being replaced by the expansion in Eq. (B3).In this case, see Fig. 13a, the sum rule is fulfilled for ã 10 1 , indicating that the present approximation is not reliable for lower values of the interaction.Again, one can improve this result by considering higher orders in the expansion (B3).As for the fidelity, shown in Fig. 13b, we find that it is surprisingly lower than what one would expect, even within the regime of validity of the present approximation (ã 10 1 ).This is possibly due to the fact that the condition 2ãn 0 ≫ 1 can only be fulfilled in the bulk of the profile, and not in the tails where the density drops.In any case, we have verified that the depletion of the fidelity comes mainly from the phase of the wave function, that cannot be reproduced with sufficient accuracy by the scaling ansatz.Conversely, from the density profiles shown in Fig. 13c, we can see that the density evolves almost self-similarly even down to ã = 10, where the fidelity is already pretty low.In this respect, the approximate self-similar approach turns out to be more accurate and effective for the NPSE (despite the complexity of having to treat the strong-and weak-interaction regimes separately), as compared to the case of the GP equation considered in Sec.IV A.

V. CONCLUSIONS
We have compared the exact evolution of an expanding 3D Bose-Einstein condensate with that obtained by means of an effective scaling approach based on the assumption of a self-similar dynamics.This approach consists in looking for self-similar solutions to be satisfied on average, by integrating the hydrodynamic equations over the spatial coordinates [1,2,15,16].Different geometries and configurations have been considered.
In case of isotropic traps, we have found that the hypothesis of self-similarity -which is exact in the noninteracting and hydrodynamic limits [5] -is indeed rather accurate even for intermediate values of the interactions.This result provides an extension to three dimensions of the 1D proof-of-concept discussed in Ref. [2].Conversely, we have found that significant deviations from the exact evolution characterize the expansion of prolate and oblate condensates.This behaviour originates from the fact that when some of the directions expand much faster than the others, they produce local variations of the density that break the hypothesis of self-similarity.Accordingly, we have found that self-similarity also breaks down for the expansion in a waveguide, due to the presence of the residual trapping.Indeed, the balance between the effect of the confining potential and that of the nonlinear term changes during the expansion owing to the local variations of the density.This in turn affects the density profile along the trapped directions, thus resulting in a behavior that is not self-similar.In the specific case of a cylindrical waveguide, we have shown that the effective scaling approach provides reliable results only for short evolution times, when the radial profile has not yet deviated from the initial profile (that is, when the local density has not decreased too much), or if the system is in the so-called 1D mean field limit [2,20], where the density profile is always characterized by a Gaussian transverse profile.
For the waveguide expansion we have also analyzed the case of the non-polynomial Schrödinger equation of Ref. [18], which has the advantage of describing the radial profile in terms of a Gaussian wave function, so that self-similarity is not explicitly violated, at least.This comes with a cost, since the non-polynomial character of the equation forces us to split the scaling equations in two different regimes, for weak and strong interactions.However, the accuracy of this approximation can be controlled by means of the sum rules that the parameters of the model have to satisfy.Here we have considered the lowest order scaling equations, finding satisfactory results in terms of the fidelity, compared to the case of the GP equation (within the regimes of validity of each approximation).We expect that a similar level of accuracy can be extended to other interaction regimes by considering higher orders in the expansion of the nonlinear term.
Finally, we remark that although the hypothesis of selfsimilarity is justified only for certain geometries and configurations, the scaling approach can still be very useful for simulating a free expansion from the numerical point of view.Indeed, even if the scaling is only approximate, solving the dynamical equations in a rescaled coordinate system can greatly reduce the size of the computational grid.

10 −1 10 0 10 1 10 2 10 3 10 4 10 5 10 6 F 10 FIG. 3 .
FIG.3.Behavior of the fidelity as a function of the interaction strength g, at different evolution times, for a spherical condensate.The black circles correspond to the density profiles shown in Fig.4.

10 FIG. 6 .
FIG.6.Behavior of the fidelity as a function of the interaction strength g, at τ f = 10, for different values of the trap anisotropy α.

10 −1 10 0 10 1 10 2 10 3 10 4 10 5 10 6 (FIG. 9 .
FIG. 9. Expansion in a waveguide: behavior of (a,b) the axial and radial scaling parameters as a function of time for different interaction strengths, and (c,d) of the corresponding values at τ f = 10 as a function of the interaction strength g.The red dots in (c,d) correspond to the values of g considered in (a,b).

10 FIG. 11 .
FIG. 11.Expansion in a waveguide: plot of the fidelity as a function of the interaction strength, at different evolution times (for α = 10).The solid black line corresponds to the case with α = 1, at τ = τ f = 10.The vertical dashed lines represent the value 4π √ α for α = 1, 10 (from left to right, respectively -see text).The black circles correspond to the density profiles shown in Fig. 10

FIG. 12 .
FIG. 12. Scaling approach for the NPSE in the weaklyinteracting limit.(a) Behavior of the constants A and B as a function of the interaction strength.(b) Behavior of the fidelity at different evolution times as a function of the interaction strength.(c) Comparison of the density profiles as given by the numerical solution of the NPSE (solid) and by the scaling argument (dashed) for different interaction strengths, at τ f = 10 [the corresponding fidelity is marked by the solid circles in panel (b)].

5 FIG. 13 .
FIG. 13.Scaling approach for the NPSE in the stronglyinteracting limit.(a) Behavior of the constants Ai (i = 1, 2, 3) as a function of the interaction strength.(b) Behavior of the fidelity at different evolution times as a function of the interaction strength.(c) Comparison of the density profiles as given by the numerical solution of the NPSE (solid) and by the scaling argument (dashed) for different interaction strengths, at τ f = 10 [the corresponding fidelity is marked by the solid circles in panel (b)].