Diffractive focusing of a uniform Bose-Einstein condensate

We propose a straightforward implementation of the phenomenon of diffractive focusing with uniform atomic Bose-Einstein condensates. Both, analytical as well as numerical methods not only illustrate the influence of the atom-atom interaction on the focusing factor and the focus time, but also allow us to derive the optimal conditions for observing focusing of this type in the case of interacting matter waves.


Introduction
When focusing an electromagnetic wave, the position and the strength of the focus are typically controlled by a lens which imprints a position-dependent phase on the incoming wave. However, focusing is possible even without a lens, namely by employing the concept of diffractive focusing, where the focus is a consequence of the non-trivial shape of the initial wave function. In this article we extend this concept to non-linear matter waves and show how it can be experimentally realized with Bose-Einstein condensates (BECs).
Already in 1816 A. J. Fresnel realized [1,2] that light passing through a circular aperture creates a bright spot on the symmetry axis before it starts to spread. This type of focusing is nowadays known as diffractive focusing and was originally described for light by the two-dimensional (2D) paraxial Helmholtz equation [3]. The same effect occurs for matter waves in one or multiple dimensions. Indeed, since the equations for electromagnetic fields within the paraxial approximation have a form similar to the Schrödinger equation of a free particle, diffractive focusing was studied [4,5] and successfully observed for water waves and plasmons [6], as well as for atoms [7], electrons [8], neutrons [9], and molecules [10]. In all these cases, the effect of diffractive focusing manifests itself [11] provided the initial wave function is (i) a real-valued one and (ii) has a non-Gaussian shape. Moreover, this kind of focusing can be very useful for the collimation of waves, such as water waves and x-rays, for which no ordinary lenses exist.
In contrast to the studies mentioned above, in the present article we explore this phenomenon for an atomic BEC [12,13] in a regime where the atom-atom interaction plays a key role. Our paper has a twofold objective: (i) We generalize the diffractive focusing effect to the case of interacting waves, and (ii) we demonstrate that, in this regime, a rather straightforward implementation is possible.
For this purpose, we consider an atomic BEC with a rectangular initial wave function emulating the case of previous studies that analyzed matter wave diffraction out of a rectangular slit [4,5,6]. In the laboratory, such shapes can be realized with BECs confined to so-called box potentials leading to uniform ground-state-density distributions due to the repulsive atom-atom interactions. The required potentials can for instance be generated by blue-detuned light sheets, in some cases combined with higher-order Laguerre-Gaussian (LG) laser beams [14,15,16].
When the trap is switched-off, the freely evolving rectangular matter wave undergoes diffractive focusing in complete analogy to a matter wave being diffracted by a rectangular aperture in the near-field or paraxial regime. The particularity of the realization we are discussing in this study is that there is no need for a dedicated aperture, but the box potential itself acts as the aperture and forms the required non-Gaussian initial state. Hence, the size of the aperture is given by the characteristic length of the box potential or by the size of the BEC itself.
This straightforward implementation of diffractive focusing occurs in all spatial dimensions where the initial wave function is close to having a rectangular shape. For the sake of simplicity, we study this effect in a quasi-1D BEC configuration that can easily be generalized to 3D. Our analytical and numerical methods, based on the Gross-Pitaevskii equation (GPE) [17,18] as well as the dynamics of Wigner functions in phase space [19], help to identify the optimal conditions for observing this type of focusing and for defining several benchmarks such as the focusing factor or the focus time.
We emphasize that our results apply not only to the physics of Bose-Einstein condensates [20], but also to other physical systems, whose dynamics is governed by the so-called cubic Schrödinger equation, for example, to nonlinear fiber optics [21] and deep water surface water waves of moderate steepness [22]. In order to study the effect of the interaction, that is the magnitude of the nonlinear term in the cubic Schrödinger equation, we propose here to observe diffractive focusing of interacting waves with a uniform BEC. Indeed, for this system all required experimental techniques already exist.
Our article is organized as follows. In Section 2, we derive the effective 1D GPE starting from the 3D GPE and introduce an optical potential by using a laser in a higher-order LG mode. We then focus in Section 3 on finding the optimal parameters for this potential to obtain a nearly rectangular ground state. Moreover, we study in detail the effect of the atom-atom interaction on the main features of diffractive focusing. Furthermore, we compare the results of our quasi-1D model to a full 3D treatment to prove their validity. In Section 4, we explain the effect of diffractive focusing in the interacting case by studying the dynamics of the Wigner function in phase space. We conclude in Section 5 by summarizing our results and discussing further interesting avenues.
Detailed calculations are presented in three Appendices. In Appendix A we derive the effective 1D GPE for the longitudinal wave function and the effective 1D interaction strength in the two limits of almost non-interacting and weakly interacting BECs. We then evaluate in Appendix B the chemical potential and the energy of a BEC in the relevant external potentials. Furthermore, Appendix C presents the Thomas-Fermi wave function for the optical trapping potential.

Theoretical foundations
In order to demonstrate the effect of diffractive focusing with a BEC we consider a quasi-1D setup that contains all relevant aspects and allows for an elementary presentation of the core features. To this end, we effectively freeze the dynamics in two dimensions, and analyze the focusing of an appropriately shaped wave function in the third dimension.
In this section we first introduce the effective 1D GPE that governs the dynamics of the BEC. We then present a special form of the box-shaped trapping potential based on a higher order LG mode.

From a 3D to a quasi-1D Bose-Einstein condensate
To arrive at a quasi-1D BEC consisting of N atoms of mass m, we start from the 3D GPE [20] i for the macroscopic wave function ψ = ψ (r, t) which is normalized according to the condition where r ≡ (x, y, z) is the position vector with the Cartesian coordinates x, y and z.
Here we assume that the atoms are interacting via a contact potential whose strength is determined by the s-wave scattering length a s , resulting in the interaction constant In addition, the external potential consists of the harmonic trap in the transverse directions determined by the trap frequencies ω x and ω y , as well as the box potential V Box = V Box (z, t) along the z-axis enforcing a rectangular ground state. Throughout this paper we consider the case where the longitudinal characteristic length L z of the external potential is much larger than the transverse one L ⊥ ≡ /mω ⊥ with ω ⊥ ≡ √ ω x ω y . In Appendix A and Appendix B we show that in this limit there is no dynamics in the transverse direction, that is in the x-y plane, as long as N a s L z . Based on these assumption we derive in Appendix A the effective 1D GPE [23] i for the wave function ϕ = ϕ(z, t) along the z-direction. Here the effective interaction strengthg is determined by the interaction constant g, Eq.
(3), the number of particles N , and the coupling parameter c ⊥ , which in general is not a constant and originates from the non-linear coupling between the transverse (x-y plane) and the longitudinal (z-axis) dynamics.
In Appendix A we derive also analytical expressions for c ⊥ in two limiting cases, namely for (i) almost non-interacting and (ii) weakly interacting atoms. If 0 ≤ N a s L ⊥ L z , the atom-atom interaction is so small that there is no coupling between the dynamics in the transverse and the longitudinal degrees of freedom, resulting in the expression (8) In the opposite limit, if L ⊥ N a s L z , we can apply the Thomas-Fermi approximation [20] for the ground state of the 3D GPE, Eq. (1), and arrive at the formula for the parameter c ⊥ .
We conclude by noting that in section 3.3 we perform full 3D numerical simulations based on Eq. (1), to test the validity of the effective 1D description. We find that for the parameters considered in this article the 1D GPE, Eq. (6), with the coupling parameter c ⊥ given by Eq. (9) describes correctly the dynamics of the BEC.

Optical box potential
We realize the box potential V Box by a LG laser beam, more precisely by the radially symmetric LG l 0 mode with the intensity profile [15] where l = 0, 1, 2, ... is the order of the mode, w 0 and P are the waist and the power of the beam, respectively. Here ρ ≡ y 2 + z 2 measures the distance from the beam axis, as depicted in Fig. 1. By choosing the waist w 0 of the LG beam to be much larger than the characteristic lengths of the harmonic trap along the transverse directions, that is when w 0 √ l L ⊥ , as shown in Fig. 1, we can neglect the cylindrical symmetry of the beam profile and use the effective 1D intensity profile I l (ρ = z) along the z-axis instead.
In addition, if the laser frequency is far blue-detuned from the atomic resonance, the associated optical dipole potential reads [15,24] V l (z) = Γ 2 8∆ Here Γ , ∆ and I s denote the decay rate, the detuning, and the saturation intensity, respectively.
With the help of the explicit expression, Eq. (10), of the intensity, we find with for the trapping potential caused by the LG mode. If the chemical potential of the ground state is much lower than the maximum V l (z l ) of the trapping potential located at z l ≡ w 0 l/2, we can approximate Eq. (12) around the potential minimum at z = 0 and obtain the power-law (13) for the trapping potential V l .
In the following we choose our parameters such that this power-law approximation is valid and Eq. (13) can be used to describe the box potential. Moreover, we note that similar box-like trapping potentials can also be realized by combining appropriate Hermite-Gauss beams with Gauss beam endcaps [14], or by employing blue-detuned painted potentials [16].

Diffractive focusing
In this section we first identify the parameters of the LG potential V l , Eq. (13), which allow us to create a quasi-1D BEC ground state wave function with the desired rectangular shape. By taking this state as the initial one, we then study in detail the effect of the atom-atom interaction on the diffractive focusing. Finally, we compare the results of our quasi-1D model to a full 3D treatment.

Preparation of the initial state of rectangular shape
We start with the discussion of the ground state of a quasi-1D BEC in the LG potential given by Eq. (13). To be specific, we throughout our article consider 87 Rb atoms [25,26] and emphasize that the values of all relevant parameters, listed in Table 1, are accessible in a state-of-the art experiment. According to Table 1 we find for the ratio N a s L ⊥ ∼ = 245, which implies that the atoms are indeed weakly interacting. Since in this case the parameter c ⊥ is given by Eq. (9), we obtain from Table 1 the value With the help of the imaginary-time propagation method [27] we have solved Eq. (6) numerically, and have obtained the wave function ϕ l = ϕ l (z, t) for the ground state of a BEC in the LG box potential V l given by Eq. (13) for different values of l, namely l = 2, 6, 10, 12, as shown in Fig. 2.
For increasing values of l the potential V l becomes more rectangular, that is flatter at z = 0 and steeper at the edges as displayed in Fig. 2 (a). Since for the ground state we consider the natural scattering length of 87 Rb we fulfill the condition of the Thomas-Fermi approximation [20] and the ground-state wave function has the form of the inverse potential, as shown in Appendix C. Thus, when the potential gets more rectangular, the ground state is more rectangular, too, as apparent from Fig. 2 (b), because the interaction between the atoms enforces a more homogeneous density distribution within the box potential.
Indeed, the fidelity between the ground state wave function ϕ l and the normalized wave function of a rectangular shape with the same amplitude h l ≡ max(ϕ l ), serves as our measure of the rectangularity of the ground state. As displayed in Fig. 3 for the fidelity F as a function of l, for l ≥ 10 the fidelity reaches 99% and we therefore consider the case l = 10 in the remainder of our article.

Influence of atom-atom interaction
Now we are in the position to analyze the effect of the atom-atom interaction on diffractive focusing. We start with the ground state of the BEC in the trapping potential V ⊥ + V 10 and then switch off only the LG potential V 10 while simultaneously changing the scattering length a s from its initial value, a s = 100 a 0 , to its final one a (f) s via a Feshbach resonance [28], where a 0 corresponds to the Bohr radius.
The resulting time evolution in the z-direction, calculated from Eq. (6) with a splitstep algorithm [29], is shown in Fig. 4 for two different values of the final scattering length, namely for a  Since the initial wave function ϕ(z, 0) has a nearly rectangular shape, we characterize the focusing effect by the focusing factor which describes the increase of the amplitude of the probability density during the dynamics in comparison with its initial value.  for four different values of l, with the initial state being the ground state of the complete potential V = V ⊥ (x, y) + V l (z) (solid lines), or the corresponding rectangular state (dashed lines), defined by Eq. (17), respectively. For a growing interaction strength the focusing factor and the focal time decrease rapidly. Hence, diffractive focusing is only visible if the atom-atom interaction is very weak during the dynamics. Figure 5 (b) shows that for a given l the focal time of the system depends strictly on the size of the initial profile, and the results for the ground-state wave function resemble the ones for the corresponding rectangular state.
One the other hand, Fig. 5 (a) shows that a appropriate fidelity, as defined by Eq. (16), is necessary to reduce the deviations in the predictions of the focus factors corresponding to ϕ l and ϕ  . Influence of the atom-atom interaction on the phenomenon of diffractive focusing of a uniform BEC apparent in the time evolution of the normalized 1D density distribution |ϕ(z, t)| 2 /|ϕ(0, 0)| 2 . The initial state ϕ(z, 0) is the ground state of the trapping potential V = V ⊥ + V 10 and the dynamics (a, b) takes place after switching off only the longitudinal potential V 10 and changing the scattering length a s from its initial value a s = 100 a 0 to its final one a ϕ l (z, 0) of l = 10 (red curve), which provides a fidelity of more than 99%, as shown in Fig. 3.
Furthermore, Fig. 6 presents a contour plot of the focusing factor f for different values of the initial and final effective interaction strength, Eq. (7),g (i) andg (f) with l = 10. For small values ofg (i) → 0 (a s → 0) the focusing effect reduces drastically, as the ground state wave function approaches the one of a particle in the box potential, that is the non-Gaussian shape vanishes.
On the other hand, we recognise that values ofg (i) = 3000 ω z already result in similar behaviour as our reference case ofg (i) = 17129.71 ω z , which corresponds to a s = 100 a 0 . Here the interactiong is used since we undergo a transition of the two limiting cases (i) almost non-interacting (smallg) and (ii) weakly interacting atoms (largeg), and in general the relation betweeng and a s is unknown. Moreover, for any valuesg (i) > 0, the focusing factor f decreases for growing values ofg (f) , as displayed in Figs. 5 and 6.
To summarize, in order to observe the phenomenon of diffractive focusing for a quasi-1D BEC, we require (i) a large initial interactiong (i) for preparing a nearly rectangular state, and (ii) a small final interactiong (f) to have a measurable effect during the dynamics. Hence, experimentally the use of Feshbach resonances [28] is mandatory to tune the atom-atom interaction in the desired way. As mentioned at the beginning of this discussion we use 87 Rb, but this effect takes place for any BEC with low temperature, for example, 39 K with a wide Feshbach resonance [30].

Justification of the quasi-1D approximation
We conclude this discussion of the phenomenon of diffractive focusing in a BEC by briefly examining to what degree the effective 1D GPE given by Eq. (6) describes the free propagation of the quasi-1D BEC. For this purpose, we again start from the ground state of the BEC in the trapping potential V = V ⊥ (x, y) + V 10 (z) and then switch off The focusing factor f depends only weakly on the initial interactiong (i) . However, a vertical cut (red line) through the contour plot at a fixed value ofg (i) = 17129.71 ω z corresponding to a s = 100 a 0 (red line) reveals that f depends strongly on the final interaction lengthg (f) , in complete agreement with the dependence shown in Fig. 5 (a). The reference case from Fig. 4 (b,d) is marked by a red cross. The corresponding scattering length giving rise to the interaction strength g can be calculated by inverting Eq. (7) and choosing an appropriate value of c ⊥ . the LG potential V 10 , while simultaneously changing the scattering length to its final value a (f) s . First, we solve Eq. (6) for the wave function ϕ = ϕ(z, t) numerically and obtain the time dependence of the normalized 1D density |ϕ(0, t)| 2 / max z |ϕ(z, 0)| 2 at z = 0.
Here we consider the two cases of almost non-interacting and weakly interacting atoms depicted in Fig. 7 by the orange and blue curve, respectively. For these limits the parameter c ⊥ , Eq. (A.12), is determined by the transverse wave function Φ 0 = Φ 0 (x, y), Eq. (A.7), and given by Eqs. (8) and (9), respectively.
Then, we perform the full 3D numerical simulation of the GPE given by Eq. (1) for the wave function ψ = ψ(r, t). The time evolution of the normalized integrated density  1). For all cases, the initial state is given by the ground state of the trapping potential V ⊥ + V 10 and the dynamics occurs after switching off only the longitudinal potential V 10 and instantaneously changing the scattering length a s from its initial value a s = 100 a 0 to its final one a is displayed by the green curve in Fig. 7. As a result, for our trap configuration, with the relevant parameters listed in Table  1, the quasi-1D approximation is very reliable and in excellent agreement with the results of the full 3D simulation.
A comparison between the curves corresponding to (i) almost non-interacting atoms (orange line) and (ii) weakly interacting atoms (blue line) with the full 3D curve reveals that the ground state obtained within the Thomas-Fermi approximation and leading to the interaction parameter c ⊥ given by Eq. (9), is more accurate in describing the dynamics of the system. This statement holds true even for small values of a (f) s , as long as the ground state was created for large values of the initial scattering length a s . We note that due to the change of the scattering length at t = 0 the transverse wave function does not describe the ground state anymore and the system undergoes collective excitations. In our 3D simulations we observed such breathing oscillations in the transverse direction. However, as a consequence of the large anisotropy of the trapping potential, the time scale of the transverse dynamics is much shorter compared to the longitudinal motion and, thus, the influence of these fast oscillations on the slower longitudinal dynamics mostly averages out for the parameters considered in this article.

Diffractive focusing viewed from Wigner phase space
This section illuminates the phenomenon of diffractive focusing of a BEC from quantum phase space. For this purpose we first recall the essential ingredients of the Wigner formulation [19] of quantum mechanics and then study classical trajectories in the absence and the presence of an atom-atom interaction. This elementary approach provides us with a deeper insight into the dynamics of the Wigner function for an interacting matter wave.

Wigner function essentials
The Wigner function [19] corresponding to the wave function ϕ = ϕ(z, t) is defined as where p is the momentum. Integration of W over p, or over z yields the relations connecting the marginals of W to the probability density distributions |ϕ(z, t)| 2 and |φ(p, t)| 2 in position and momentum space, respectively [19]. Herẽ is the momentum representation of the wave function ϕ = ϕ(z, t).
Although the Wigner function W is normalized, that is these properties do not imply that W is always positive. Indeed, the Wigner function is a quasi-probability distribution [19] and its negative parts reflect the quantum features of the system under consideration.

Classical trajectories
Instead of deriving and solving the dynamical equation for the Wigner function corresponding to the 1D GPE, Eq. (6), we obtain the time-dependent Wigner function directly from the definition, Eq. (20), of W in terms of the time-dependent wave function ϕ = ϕ(z, t) determined by solving Eq. (6) numerically. In order to visualize the dynamics in phase space, we take a point {z, p} in phase space and find the "classical trajectories" {Z(t), P (t)} governed by the Hamilton equations d dt subjected to the initial conditions Z(0) ≡ z and P (0) ≡ p. Here the classical Hamiltonian corresponds to the 1D GPE, Eq. (6), without the trapping potential. We emphasize that the use of Eqs. (25), (26) and (27) implies the knowledge of the wave function ϕ = ϕ(z, t) at all times obtained by numerically solving Eq. (6).

Time evolution without atom-atom interaction
Before we consider the case of interacting particles, we first recall [5,31] the interactionfree dynamics (a  (22), the integration over the momentum or the position variable provides us with the position distribution |ϕ(z, t)| 2 (lower sub-figure), or the momentum distribution |φ(p, t)| 2 (left sub- figure). Figure 8 brings out most clearly the origin of the phenomenon of diffractive focusing. Indeed, at t = 0, the Wigner function W = W (z, p; 0) exhibits both positive and negative values. During the free expansion, t > 0, the parts of the Wigner function corresponding to p > 0 (p < 0) move to the right (left) along the straight lines {z + pt, p}, displayed in Fig. 8 (d) by different colors for four different initial points in phase space. These lines are parallel to the z-axis and describe the free classical motion, resulting from Eqs. (25), (26) and (27) At the time of focusing, t = t f , the position distribution |ϕ(z, t f )| 2 features a narrow maximum at z = 0. Indeed, integration over p in Eq. (21) for fixed position z, yields a maximum only at the values of z, which correspond to the maximal values of W = W (z, p, t f ), displayed by dark red color. According to Fig. 8 (a)-(d), this is the line z = 0 in phase space. In other words, focusing takes place at z = 0, because at t = t f all negative parts of the initial Wigner function W (z, p; 0) have moved away from the p-axis [4,5]. However, they now subtract from the positive parts of the wings and make the distribution in space even narrower.
For t > t f , the negative parts of W (z, p, 0) have moved further away from the line z = 0. Since the positive parts of the original Wigner function are at lower momenta than the negative ones, they move slower and are therefore left at the center of the phase space. They are the origin of the spreading of the position distribution |ϕ(z, t)| 2 .

Time evolution with atom-atom interaction
Next, we discuss the nonlinear time evolution of the Wigner function in the case of a non-vanishing atom-atom interaction, namely for the final value of the scattering length a (f) s = 0.58 a 0 . Figure 8 (e) presents the same initial Wigner function W (z, p; 0) as Fig.  8 (a).
In contrast to the case of no atom-atom interaction, Figs. 8 (f)-(h) indicate that the parts of the Wigner function corresponding to p > 0 (p < 0) do not only move to the right (left) but also move up (down). This effect can be explained as follows.
For short times, t t f , we can neglect in Eq. (6) the kinetic energy term (− 2 /2m)∂ 2 /∂z 2 compared to the interaction termg|ϕ(z, t)| 2 and arrive at the nonlinear equation which is not easy to solve. However, we note that Eq. (28) conserves the quantity |ϕ(z, t)| 2 , that is ∂|ϕ(z, t)| 2 /∂t = 0, resulting in the simplified equation with the solution Thus, for t t f , the wave function ϕ = ϕ(z, t) picks up only a positiondependent phase determined by the initial distribution |ϕ(z, 0)| 2 and the effective interaction strengthg. The gradient −gt∂|ϕ(z, 0)| 2 /∂z of this phase defines the increase in momentum, which is a function of z.
This result also follows from Eqs. (25) and (27) and is displayed in Fig. 8 (h) by the classical trajectories corresponding to different initial points in phase space. Hence, due to this increase in momentum, the negative parts of the Wigner function get deformed and move faster away from the center compared to the caseg = 0, resulting in the focus appearing at earlier times. This behavior is also confirmed by Fig. 4 (b).
For t > t f , the position distribution |ϕ(z, t)| 2 spreads further, as shown in Fig. 8  (h), and reduces its amplitude. Thus, the interaction termg|ϕ(z, t)| 2 in Eq. (6) gets smaller compared to the kinetic energy term, and the evolution of the wave function ϕ = ϕ(z, t) can be described solely by the Schrödinger equation. This effect is illustrated in Fig. 8 (h) by the classical trajectories, which in the long-time limit are again parallel to the z-axis.
A closer look at the momentum distribution of Fig. 8 (h) reveals that two maxima are forming symmetrically around the origin of phase space due to the non-vanishing interaction. The momenta at which they occur directly depend on the final interaction strengthg (f) of the system and increases for largerg (f) . For very long times the peaks smear out into a large momentum distribution at the center, however, if the interaction is set to zero beforehand, then the double-peak structure could be preserved.
In summary, we emphasize that diffractive focusing originates from the negative parts of the Wigner function [4,5]. According to the Hudson theorem [32], a pure state with an initial Gaussian profile has a positive Wigner function at any point of the phase space, and therefore the state does not show diffractive focusing at all. A similar behavior occurs for any classical state. Hence, diffractive focusing for a given BEC gives us an opportunity to check whether this BEC is prepared in a non-Gaussian or non-classical state.

Conclusions and outlook
In this article we have studied the phenomenon of diffractive focusing of interacting matter waves employing analytical as well as numerical methods. We have proposed a straightforward implementation of this effect with an atomic BEC confined by a box-like trap as realized for instance in Ref. [16]. The interaction of the atoms forming a BEC leads to the non-linearity in the GPE and is an essential ingredient in the preparation of a rectangular wave function, previously studied [4,5,6] in the context of Schrödinger waves, or the paraxial approximation in optics, and obtained by a rectangular slit.
As benchmarks, we have identified the focusing factor and the focus time which both are functions of the strength of the atom-atom interaction. These measures allow us to derive the optimal conditions for observing this type of self-focusing of a BEC. Having identified the origin of diffractive focusing for interacting matter waves, illuminated by the time evolution of the Wigner function in phase space, we conclude that the cleanest realization occurs when the atom-atom interaction is switched off during the dynamics by a magnetic Feshbach resonance.
For the sake of simplicity we have restricted our treatment to a quasi-1D case. However, the effect of diffractive focusing takes also place for higher dimensions [33] and could be realized with a 3D box potential [16] generated by blue-detuned laser light. Indeed, the focus factor achieves the value 4 for a cylindrically symmetric rectangular shape in two dimensions [1], whereas it is 1.8 for the rectangular initial profile in one dimension [4,5]. Moreover, diffractive focusing crucially depends on the initial profile, that is a more non-Gaussian or non-classical initial wavefunction results in stronger focusing. Thus, it would be useful to find the optimal initial profile giving rise to the best focusing. The problem of finding such an optimal state determined solely by the atom-atom interaction is highly non-trivial due to the non-linearity of the dynamical equation.
We conclude by emphasizing that the results presented here can be immediately applied to other physical systems, whose dynamics is governed by the Gross-Pitaevskiitype equation, that is the cubic Schrödinger equation, for instance, to nonlinear optics [21] and deep water surface water waves of moderate steepness [22]. Moreover, the diffractive focusing can be used to generate bright sources of matter waves for dedicated applications in precision measurements [34,35]. A more detailed discussion of these points goes beyond the scope of this article and has to be postponed to a future publication.

Appendix A.1. Decoupling of transverse and longitudinal dynamics
To derive an equation governing the dynamics of a quasi-1D BEC which consists of N atoms of mass m, we start from the 3D GPE [20] i for the BEC wave function ψ = ψ (r, t) with the external potential being the sum of a harmonic trap in the transverse directions determined by the trap frequencies ω x and ω y , and a box potential V Box = V Box (z, t) yielding trapping in the z-direction.
In this article we consider a highly anisotropic cigar-shaped trapping geometry defined by the relation where L z is the longitudinal characteristic length of the external potential and L ⊥ ≡ L x L y is the transverse one with L x ≡ /mω x and L y ≡ /mω y . In this case, as shown in Appendix B, the total energy per particle of the BEC is approximately given by the relation with ω ⊥ ≡ √ ω x ω y .
Hence, for 0 ≤ N a s L z , (A.6) the total energy per particle E/N is much smaller than the characteristic energy scale ω ⊥ of the transverse direction, making it impossible to drive collective excitations in that direction as long as the energy of the system is conserved. Indeed, we effectively freeze out the transverse dynamics [36,37,38] Consequently, the total wave function can be approximated by the product of the real-valued wave function Φ 0 = Φ 0 (x, y) describing the ground state in the transverse direction, the wave function ϕ = ϕ (z, t) along the z-direction, and a time-dependent phase factor, where the constant ε 0 shall be determined later as to simplify the equations. Moreover, the function Φ 0 is chosen to be normalized, that is When we insert our ansatz, Eq. (A.7), into the 3D GPE, Eq. (A.1), we obtain the identity Finally, we multiply both sides of Eq. (A.9) from the left by Φ 0 , integrate over x and y, and arrive at the non-linear equation  According to Eqs. (A.4) and (A.6), there exist two distinct cases where both Φ 0 andg can be calculated analytically, namely the limit of almost non-interacting atoms, 0 ≤ N a s L ⊥ L z , and the case of weakly-interacting atoms, L ⊥ N a s L z . In the next sections we consider these two situations.

Appendix A.2. Almost non-interacting atoms
For almost non-interacting atoms with 0 ≤ N a s L ⊥ L z , we neglect the interaction term gN |ψ| 2 in the 3D GPE (A.1), and the equation becomes approximately separable in the coordinates x, y, and z. As a result, the transverse wave function Φ 0 in the ansatz, Eq. (A.7), for ψ coincides with the wave function of the ground state of a 2D harmonic oscillator.
By inserting Eq. (A.14) into the definitions for c ⊥ , Eq. (A.12), and ε 0 , Eq. (A.13), and performing the integration over x and y, we obtain the explicit expressions 15) for the parameter c ⊥ and for the constant ε 0 . This approach to quasi-1D BECs has already been discussed in similar ways by other groups [36,37,38]. Our results exactly coincide with their findings for the same order of approximation.

Appendix A.3. Weakly interacting atoms
In the case of weakly interacting atoms, that is for L ⊥ N a s L z , the interaction term gN |ψ (r, t) | 2 in the 3D GPE (A.1) is the leading one and we can apply the Thomas-Fermi approximation [20] by neglecting the kinetic term when determining the ground-state wave function. Starting from the stationary solution we thus obtain the ground-state wave function Here Θ denotes the Heaviside function and is the chemical potential derived in Appendix B when the box potential is approximated by infinitely high potential walls separated by 2L z . As a result, within the Thomas-Fermi approximation, the total wave function φ = φ(x, y, z) given by Eq. (A.18) is again the product of the normalized transverse wave function .21) and the longitudinal wave function  .13), and neglecting the second-order derivatives over x and y, we arrive at the formula for the constant ε 0 , where µ is given by Eq. (A.19). We emphasize that the expression, Eq. (A.23), for c ⊥ is still obtained in the regime where the motion along the transverse direction is effectively frozen out (N a s L z ), but in contrast to the previous case the interaction between the particles is taken into account when determining the shape of the transverse ground-state wave function Φ 0 . The comparison between the solutions of the effective 1D GPE (A.10) with c ⊥ given by Eq. (A.23), and the 3D GPE (A.1), presented in Fig. 7, shows that, for the parameters considered in this article, our expression, Eq. (A.23), for c ⊥ describes the dynamics more accurately than the standard formula, Eq. (A.15), corresponding to weakly interacting atoms.
We conclude this discussion by noting that the case of even stronger atom-atom interaction, when L z N a s , can be treated in a similar way. Indeed, the dynamics along the transverse direction is then much faster compared to the longitudinal direction due to the relation L ⊥ L z . Here one can perform the adiabatic approximation [37,38] to factorize the total wave function and to describe the longitudinal dynamics of the quasi-1D BEC.
In order to derive an analytical expression for µ, we approximate V Box by the potential of infinitely high walls where b 2 x ≡ 2µ/(mω 2 x ) and b 2 y ≡ 2µ/(mω 2 y ). By introducing the polar coordinates x ≡ b x r cos θ and y ≡ b y r sin θ with 0 ≤ r ≤ 1 and 0 ≤ θ ≤ 2π, we arrive at for the chemical potential of a BEC being confined by an infinitely high box potential along the z-axis and two harmonic potentials along the x-and y-direction. We emphasize that Eq. (B.7) is valid for arbitrary length scales L ⊥ and L z of the external potentials.

Appendix B.2. Energy
A similar calculation can be performed to find the total energy [20] E = N dxdydz V (x, y, z) |φ (x, y, z)| 2 + g 2 |φ (x, y, z)| 4 (B.8) of a BEC within the Thomas-Fermi approximation. However, the more convenient approach consists of inserting the result, Eq. (B.7), for the chemical potential into the definition [20]  for the total energy per particle within the Thomas-Fermi approximation.
Appendix C. Thomas-Fermi wave function for optical trapping potential In this appendix we derive the wave function for ground state of the potential V l , Eq. reads [20] ϕ T F (z) = µ T F − V l (z) g (C.6)