The Influence of s-wave interactions on focussing of atoms

The focusing of a rubidium Bose-Einstein condensate via an optical lattice potential is numerically investigated. The results are compared with a classical trajectory model which under-estimates the full width half maximum of the focused beam. Via the inclusion of the effects of interactions, in the Bose-Einstein condensate, into the classical trajectory model we show that it is possible to obtain reliable estimates for the full width half maximum of the focused beam when compared to numerical integration of the Gross-Pitaevskii equation. Finally, we investigate the optimal regimes for focusing and find that for a strongly interacting Bose-Einstein condensate focusing of order 20 nm may be possible.


I. INTRODUCTION
Atom lithography is a technique where the gradient forces applied by laser fields on a beam of atoms are used to direct the atoms into nanostructures deposited on a plane surface 1,2 . It has been the topic of considerable study in the field of atom optics over the last two decades as it provides a scheme of writing nanometer structures directly onto a substrate in parallel process 3 . The main application in the development of nano-lithography could be the race for an increased density of transistors in the computer chips 4 . The possibility of using light to deposit feature sizes of a few nanometers was first suggested in 1987 by Balykvin and Letokhov 5,6 . Later in 1991, Mc-Clelland and Scheinfein 7 proposed a particle optics approach in which the atom lens created by the laser light can focus an atomic beam down to a surface. The principle was experimentally demonstrated by Timp 1 in 1992 using Na deposited on Si via a standing light wave, and was similarly followed by McClelland et al. 2 to focus Cr in the presence of a 1D lattice. Atom lithography continued to develop to two-and three-dimensional nanostructures in a single process utilizing the selectivity of atom-light interaction (using 2D or 3D lattices) 8 . There have also been demonstrations depositing Yb 9 and Fe 10,11 . Almost all the experiments accomplished so far, have used an oven source of atoms in which the beam is collimated with an aperture followed by a transverse laser cooling process 2 before traveling through a focusing potential. The traditional study of focusing in optical lattices uses the classical trajectories model of atomic motion when travelling through the light 7 . The principle is based on atom-light interactions, resulting from the dipole force 12 which causes neutral atoms to become manipulated with a near resonant laser light [12][13][14] . Ideally, this model predicts almost perfect focusing in the absence of interactions.
Nevertheless, there are some advantages in using a Bose-Einstein condensate (BEC) of neutral atoms as the source of atom deposition. Since the de-Broglie wavelength of a gas of atoms is of the order of the mean field distance between particles, an ultra cold source such as a BEC would bring atoms to wavelengths of ∼ nm or pm for nano-Kelvin temperatures 15 resulting in an excellent collimation of the beam of atoms as well as a high flux density 16 . Using a BEC source can also reduce effects such as chromatic aberration and angular divergence 2 , with the longitudinal and transverse velocity distributions typically being much lower when incident on the surface compared to those resulted from thermal sources.
However, for BEC's, the effect of interactions must be accounted for by introducing the s-wave interaction between atoms. The investigation of the matter wave focusing dynamics of a trapped BEC for both an interacting and non-interacting case was conducted theoretically by D. Muuray and P. Ohberg in 2004 17 . In their work, they derived the time to focus as a function of the focusing strength for a 3D Bose-condensed cloud. The significance of this research is to take atomic interactions into account when focusing a free propagating BEC and to scale its effect on the broadening of the focal spot sizes and peak densities achievable in realistic nano-lithography experiments.
The principle of atom lithography using a BEC is schematically depicted in Fig 1. The cloud of 87 Rb is initially confined by a harmonic trap. Turning off the trap, the released condensate expands whilst propagating along the vertical axis. It then encounters the optical lattice being focused by its nodes along the horizontal axis resulting in a large periodic array of nano-structures that are separated by λ/2 where λ is the wavelength of lattice.
In this paper, the focal properties of an optical lattice are derived by the classical equations. Using this model, we estimate the Full Width at Half Maximum (FWHM) as well as peak density of 87 Rb focused structures, which can be further used to determine the possible contributions of structure broadening such as the spherical and diffraction aberrations. Following this, the Gross-Pitaevskii Equation (GPE) is introduced and its application for the purpose of BEC lithography is studied. Conducting various numerical simulations via the GPE for different 87 Rb BEC and focusing potential parameters, we estimate resultant focused structure resolutions and peak densities. Not only is the key role of the FIG. 1. Schematic representation of atom deposition using a 87 Rb BEC focused by an optical lattice. The condensate is initially confined by a harmonic trap. Once the BEC is released from the trap, it starts propagating whilst expanding along the falling (longitudinal) direction. An optical potential with a sinusoidal configuration along the transverse axis focuses the BEC resulting in the periodic atomic distribution.
s-wave interaction within the BEC investigated through the GPE, but also by taking the transverse and longitudinal velocity distributions of the condensate, the classical trajectories model is developed to consider the influence of atomic interactions on focused profiles.

II. THE CLASSICAL TRAJECTORIES APPROACH
When atoms are exposed to a potential field, an atomic dipole moment is induced by the electric field. The induced dipole moment interacts with the gradient of the field resulting in a dipole force gradient 18 applied towards the nodes or anti-nodes of the periodic field. The resultant dipole potential on stationary atoms is well established 18 . However, when atoms carry a momentum, the created optical dipole potential behaves as a focusing thick lens on neutral moving atoms 19 : where ∆ denotes the detuning of the laser frequency from the atomic resonance, is the Plank constant, γ is the natural linewidth of the atomic transition (spontaneous decay rate from the excited level) and I s is the saturation intensity related to the atomic D 2 line transition, 5 2 S 1/2 −→ 5 2 P 3/2 , for 87 Rb. While a red detuned laser light from resonance, ∆ < 0, in 87 Rb D 2 line, directs the falling atoms towards the nodes of the standing potential, a blue detuned light, ∆ > 0 brings the atoms to the anti-nodes of the focusing potential (∆ = ω L − ω 0 where ω L and ω 0 are respectively the frequency of incident laser light and resonance between 5 2 S 1/2 and 5 2 P 3/2 ). We consider a geometry for the optical potential such that the laser intensity profile, I(x, z), shapes a mask with a periodic scheme of multiple focusing lenses along the direction of focusing atoms. Hence, a Gaussian standing potential is considered which contains of a sinusoidal behavior of the intensity along the x-axis (the focusing direction) and a Gaussian envelope function along the zaxis, the direction of falling atoms. This optical lattice potential is represented by where I 0 is the maximum intensity of the spatially varying Gaussian profile, σ z is the radius of the beam at 1/e 2 value of the maximum intensity and k = 2π/λ is the wavenumber. Since almost no force is applied to atoms along the y direction compared to the other two directions, neglecting this causes the focusing potential, Eq. (1), to become independent of y.
In order to scale the parameters involved in an optimal focusing potential, we exploit the classical trajectories approach 2 . Neglecting the y axis due to the symmetry of problem, the classical equations of motion for atomic trajectories are given by Using conservation of energy, one can combine Eqs. (3) and (4) solving for x as function of z where E 0 represents the total energy of each atom and x = dx dz . To obtain the focal properties of the lattice, we consider the paraxial approximation 20 in which the sinusoidal term of the potential consisting of a multi-node wave along the x axis is converted to a single node harmonic part. In addition, the approximation neglects aberrations, i.e. the trajectories are assumed to be perfectly parallel to the z axis when falling towards the lattice. These can be implemented by U (x, z) E 0 , dx dz 1 and kx 1. Applying these, Eq. (5) converts to Now considering sin(kx) ∼ kx, sin 2 (kx) ∼ 0 and cos(kx) ∼ 1, one obtains the following second order differential equation According to the relation between the maximum intensity and the corresponding value of power in a standing wave Gaussian beam 20 , I 0 = 8P 0 /πσ 2 z , the required power value of the harmonic potential to focus atoms at any desired spots along the focal axis (z-axis) is achieved as function of the potential factors and atoms' initial kinetic energy, where ξ = q 2 σ 2 z is a dimensionless parameter. We now proceed with a specific example for focusing ultra-cold 87 Rb atoms using the paraxial approximation displayed in Fig. 2. To meet the requirements of approximation, the potential wavelength is chosen to be 400 times greater than the actual 87 Rb D 2 line wavelength, λ = 400λ D2 = 312 µm where λ D2 = 780.027 nm. This forms a potential with a harmonic distribution along the x axis. While the potential radius is adjusted to σ z = 100 µm, ξ = 5.37, 3.37, 2.37, and 1.37 corresponding to the potential power of P 0 = 43.018 , 26.996, 18.986, and 10.975 µW respectively provided that the initial velocity of atoms and the detuning from resonance are selected as v i = v z = 1 cm/s and ∆ = 200 GHz. For each value of ξ or P 0 , two 87 Rb atomic trajectories falling symmetrically from x i = ±λ/4 = ±78 µm landing at x f = 0 and a particular z f are plotted. As a case in point, for ξ = 5.37, solving Eq.(7) yields the focal point optimally placed at z f = 0 and x f = 0 (the center of potential). However, reducing power values results in the focal point being located at z < 0 below the potential center.
We note that this is an ideal scenario in which for a particular value of P 0 , all atoms are focused exactly at one spot and no aberration is considered. However, in practice there always exists a spherical aberration in optics while atoms are entering a thick lens 21 . This arises from the fact that atoms traveling farther to the focal axis will experience a weaker force than those that are closer to this axis. This effect can be considered through the exact numerical solution of Eq.(5). As an illustration, the atomic trajectories for ultra-cold 87 Rb starting at z i = 20 µm whilst moving at a longitudinal velocity of v z = 1 cm/s through the potential of a radius size of σ z = 10 µm focused at z f = 0 (ξ = 5.37) have been numerically simulated and are depicted in Fig 3. Here, we choose a lattice whose wavelength is 16 Cross-section view (top view) for the 87 Rb classical atomic trajectories for different values of ξ and consequently different harmonic potential powers, P0, when applying the paraxial approximation. The solid and dashed curves respectively represent the area above (z > 0) and below (z < 0) the center of potential (z = 0). The two gray-shaded oval areas illustrate the harmonic potential intensity profile, I(x, z), which are darker for higher intensity regions. The maximum intensity value for ξ = 5.37, 3.37, 2.37, and 1.37 is respectively I0 = 1.095 × 10 4 , 6.874 × 10 3 ,4.834 × 10 3 , and 2.794 × 10 3 W/m 2 . A choice of λ = 400λD 2 = 312 µm for the potential wavelength necessitates the two adjacent peaks (located at x = ±λ/4 = ±78 µm, z = 0) to be apart by λ/2 = 156 µm along the x axis. The trajectories corresponding to ξ = 5.37 (brown curves) are focused at x f = z f = 0 while setting ξ = 3.37, 2.37, and 1.37 shifts the focus points to z < 0. Parameters used in the simulation are: σz = 100µm, ∆ = 200 GHz, γ = 37 MHz, Is = 16.5 W/m 2 , m = 1.44 × 10 −25 Kg and vi = vz = 1 cm/s. times larger than the actual wavelength of 87 Rb D 2 line, λ = 16λ D2 = 12.48 µm. This is practical in realistic experiments through the use of a Spatial Light Modulator (SLM) 22,23 . Aligning the laser detuning to ∆ = 200 GHz, the required optimal value of lattice power and maximum peak intensity to focus the atoms to the center of potential (z = 0) are calculated as P 0 = 0.068 µW and I 0 = 1.752 × 10 3 W/m 2 . Unlike the paraxial solution, for each lattice node here, the atoms do not fully land at the same focal point, and this process is identical for all lattice nodes implying the spherical aberration. Hence, one could evaluate the linewidth of the created structure and estimate the broadening contribution arising from this effect inside every lattice slit, D = λ/2 = 6.24 µm. For instance, we have selected the lattice central node, from x = −λ/4 = −3.12 µm to x = λ/4 = 3.12 µm, and plotted a histogram distribution for 62400 trajecto- However, the spherical aberration is not the only reason for profile broadening. In a lens like lattice, there are finite adjacent apertures positioned in a distance of λ/2 apart along the entire lens. Since atoms exhibit wave-like behavior, they interfere together with their de-Broglie wavelengths, λ dB , while crossing through the apertures. This causes a limit on the linewidth of structures known as the diffraction effect. The angular resolution produced by the diffraction is estimated by Rayleigh Criterion 24,25 where the factor β = 1.22 accounts for a circular aperture 24 while β = 0.88 is dedicated to a rectangular (or cylindrical) aperture 20,26,27 . The de-Broglie wavelength of atoms is defined as is the most probable longitudinal velocity. The angular resolution in Eq. (10) can be converted to the diffractionlimited FWHM by where θ is considered to be a very small angle. Since each aperture has a rectangular structure to the incident atomic beam, (∆x) diff would be For the previous example, for the atoms with a mass of m = 1.44 × 10 25 Kg falling at v z = 1 cm/s through a lattice of wavelength of λ = 312 µm, the diffraction-limited FWHM is estimated as (∆x) diff = 0.3589 µm. In this calculation, the focal length (the difference between the coordinates of the focal point and the principal plane location when focusing at z = 0), f = 5.531 µm, σ z = 10 µm and λ = 12.48 µm. According to Eq. (12), this broadening is proportionally related to the focal length and is inversely associated with the velocity of atoms as well as the size of apertures. Hence, to minimize (∆x) diff , one needs to either use a relatively shorter focal length, a higher longitudinal velocity or a larger potential wavelength. Eventually, taking into account the total structure broadening predicted by the classical trajectories model in absence of atomic interactions, and considering the results for (∆x) sph of the stated example, we can estimate the total value of FWHM via the classical trajectories model for the cases, z f = 0 (with σ z = 10 µm), which is evaluated as (∆x) Classical = 0.495 µm.

III. THE GROSS-PITAEVSKII EQUATION METHODOLOGY
Using a BEC 28 as a source of ultra-cold atoms brings several advantages to atom deposition as it can significantly reduce the linewidth of longitudinal and transverse velocity distributions providing excellent coherence and collimation for the atomic beam as well as offering relatively small de Broglie wavelengths, high peak densities and quality spatial modes 15,16,29 . In this section, we take into account the atomic interaction when focusing a free propagating 87 Rb BEC. The impact of interactions on the broadening of the nano-focal spot sizes and peak densities are estimated, which is accomplished through the use of the GPE 30,31 .
The three-dimensional time dependent GPE modeling the dynamics of a BEC is represented by 30,32-34 where ψ(r, t) indicates the BEC wavefunction at different times of propagation, V ext is the time-dependent external potential applied on the BEC, m is the atomic mass and is the Planck constant. The interactions between atoms within the cloud are considered using a non-linear mean field potential, V mean , which estimates the averaged exerted potential on any particular atom by all other atoms given by 29,[35][36][37] where quantifies the atomic interactions, |ψ(r, t)| 2 describes the atomic density, and a s is the s-wave scattering length. The value a s can be practically tuned from a s > 0 (repulsive interactions) to a s < 0 (attractive interactions) utilizing Feshbach resonance 38 .

A. The BEC Ground State
In order to generate the BEC ground state wavefunction at t = 0, we assume that the condensate is initially confined via a harmonic trapping potential defined by where ω x , ω y , ω z represent the harmonic trap frequencies along the x, y, z axes respectively. For our purpose, we consider a cylindrical (cigar-shaped) condensate with two radial axes associated with the two tight trap frequencies, ω y and ω z , and one axial axis corresponding to the weak trap frequency, ω x , where ω y , ω z > ω x . This configuration allows the BEC to be elongated along the x axis compared to the y and z axes. The Thomas-Fermi solution 39,40 to Eq. (14) can be used as an initial function to Eq. (14) to acquire the exact solution for the ground state wavefunction of system, ψ g (r, t = 0), which itself is exploited as an initial condition to Eq. (14) to achieve the propagation state, ψ(r, t > 0). The process of calculating ψ g (r, t = 0) and ψ(r, t > 0) is numerically conducted using an imaginary and real time step respectively via the Embedded Runge-Kutta scheme along with adaptive Fourier split-step size 41 . We consider 10 5 atoms trapped by ω x = 2π × 10 Hz and ω y,z = 2π × 70 Hz producing a cylindrical BEC elongated along the x axis. Once the harmonic trap, V trap , is turned off (ω i = 0, i = x, y, z), the condensate starts expanding due to the s-wave interaction between atoms.

B. The Time Dependent Focusing Potential
For simplicity of numerical calculation for the evolving BEC through a focusing potential, we assume that the BEC is located in a stationary frame while the optical lattice is situated in a moving frame approaching the BEC along the z axis. Hence, V lattice would be dependent on time by z(t) = 1 2 gt 2 + v 0 t, which is the varying distance as a function of time following the free falling method where g and v 0 are the gravity and initial velocity kick respectively. Thus, once the confining potential is switched off at t > 0, the optical lattice potential [see Eq. (1)] is switched on so that where z 0 denotes the initial distance between the center of lattice and the center-of-mass of condensate along the z axis. As a result, combining Eqs. (14), (15), (16), and (18), the dynamics of a focused BEC at t > 0 would be given by the following equation To make a direct comparison between the output of the GPE and classical trajectories models, we exploit the same example in section II in which the 87 Rb BEC  . For this case, provided that the laser detuning is set to ∆ = 200 GHz, the optimal lattice power and maximum peak intensity values to focus the BEC at z = 0 are required as P 0 = 0.068 µW and I 0 = 1.752×10 3 W/m 2 respectively. Applying the stated factors along with the relevant data for the mass, saturation intensity and spontaneous emission rate associated with 87 Rb D 2 line in Eq. (19), the BEC focusing dynamics is numerically achieved, and the full process as well as results are shown in Figs 5(a-f).
Here we have considered two different cases to examine the focal spot sizes and peak densities. Firstly, the BEC is strongly interacting and the s-wave scattering length is chosen as a s = 100a 0 leading to a repulsive BEC [see

C. The Variation of BEC and Potential Factors
In this section, we consider more examples to investigate the influence of altering the BEC longitudinal velocity and lattice radius size on the deposited profiles. In addition to the previous instance discussed in section III B for the focused profile with σ z = 10 µm, v z = 1 cm/s, two more rounds of simulations are considered accounting for both interacting and non-interacting BEC's. While parameters such as ξ = 5.37, λ = 12.48 µm, ∆ = 200 GHz remain unchanged, we consider σ z = 20 µm, v z = 1 cm/s and σ z = 10 µm, v z = 2 cm/s. The results for FWHM and central peak density values are summarized in Table I. We notice that doubling the lattice radius size for the same BEC velocity leads to a reduction almost by half in the structure peak, and an increase by a factor of two in the structure linewidth for both a s = 100a 0 and a s = 0. Furthermore, for the two cases, increasing the condensate velocity for a given potential radius results in a decline by half in the FWHM and a growth by a factor of two in the peak density. This arises from the fact that the potential peak intensity is directly proportional to the square of BEC velocity whereas it is inversely proportional to the square of lattice radius size, I 0 ∝ v 2 z /σ 2 z .

IV. VELOCITY DISTRIBUTION OF A BEC
Due to the s-wave interactions between the atoms within a condensate, the velocity of atoms does not remain constant over time. Hence, a distribution is formed in both the longitudinal (along the z axis) and transverse (along the x axis) velocity profiles so that each one encompasses an associated peak representing the most probable velocity. Taking a Fast Fourier Transform (FFT) directly from the BEC density profile in position space, one is able to gain the BEC density profile in momentum space at different times of propagation. Since inter atomic interactions cause the condensate to expand (for a s > 0) over time the density profile in momentum space also spreads.
As an illustration, we have studied separately the BEC transverse and longitudinal velocity distributions at t = 2 ms in Figs 6(a-c). Here, the BEC is kicked by v z = 1 cm/s when being released from the trap. Neglecting the gravity acceleration in BEC's falling, the linewidth for each distribution is derived by applying a Gaussian fit [see Figs 6(b, c)]. The values of FWHM for the transverse and longitudinal profiles at t = 2 ms are estimated as ∆v x = 0.0264 cm/s and ∆v z = 0.1791 cm/s respectively. It is clear that the tight trap frequency along the z axis (ω z = 2π × 70 Hz) causes a significantly enhanced widening velocity profile along the falling axis compared to the horizontal axis (∆v z ∆v x ). In sections IV A and IV B, we explore the possibility of implementing the information from ∆v z and ∆v x to the classical trajectories model to predict the resultant structure broadenings.

A. The Chromatic Aberration in Classical Trajectories Model
In optics, the chromatic aberration occurs because lenses have different refractive indices for different wave- the resultant broadening along the x axis is given by where is the angle between the incident atomic beam and the focal axis, and D is the lens slit size (which is the distance between two adjacent peaks in an optical lattice) given by D = λ/2. Since a change in velocity, v z , would vary the focal length, f , one can write where the variation of focal length with respect to the longitudinal velocity of atoms, df dvz (or the kinetic energy of atoms) can be broken down as The dimensionless parameter, ξ, in Eq. (23) is a function of E 0 and consequently a function of v z [see Eq. (9)], represented by where C is a constant coefficient. Now, taking into account that and using Eq. (23), Eq. (22) converts to Finally, substituting Eqs. (26) and (21) into Eq. (20), the broadening in the focal spot sizes resulting from longitudinal velocity spread is described as where the data for df /dξ for a range of desired focal lengths is derived from Fig 8 estimated by the paraxial approximation in section II. We now calculate the contribution of the chromatic aberration in broadening the focused structure for the example mentioned in sections II and III B. According to Fig 8,to focus the atoms to the center of potential, ξ = 5.37, one would obtain dF/dξ = −0.0249 (F is the dimensionless focal length). Then, for a lattice with a size of σ z = 10 µm, one would calculate df /dξ = σ z (dF/dξ) = −0.249 µm. Hence, given a lattice wavelength of λ = 12.48 µm, a focal length of f = 5.531 µm (for focusing to the center of the potential with σ z = 10 µm, and ξ = 5.37), the longitudinal velocity and linewidth, v z = 1 cm/s and ∆v z = 0.1791 cm/s at t = 2 ms (for the BEC falling from z i = 20 µm to z f = 0), the resultant broadening magnitude is (∆x) c = 0.405 µm.
It is clear that the most significant factor in the broadening is due to the longitudinal velocity linewidth. As a result, to reduce (∆x) c , one would better choose a relatively short propagation time (or a small z 0 ) for the condensate to prevent the longitudinal velocity profile from a high spread. Moreover, from Eq. (27), one can understand that larger v z values can also result in less broadening in the focused structures, which is in an agreement with what is predicted by the diffraction effect [see Eq. (12)].

Atom (1)-Collimated
Atom (2)-Not collimated (2) (1) FIG. 9. A schematic illustration of two atoms with an angle of θ with respect to each other crossing through a thick lens with two principal planes being focused at different points.
The focal length, f , is considered for atom (1) collimated to the z axis while f1 is associated with atom (2) moving under an angle of θ with respect to the z axis. The virtual object of a size of ∆x0 is located at z = −z0 on the left side of the lens, and the created image of a size of (∆x)a appears on the right side of the lens.

B. The Angular Divergence in Classical Trajectories Model
Another broadening which appears in the deposited structures results from the a divergence in the falling beam of atoms. The divergence is sourced from the transverse velocities of the atoms within the cloud, v x = 0. As shown in Fig 9, we consider two atoms approaching a thick lens. Atom (1) is collimated to the focal axis (z axis), and atom (2) is moving towards the lens with an angle of θ with respect to the focal axis. The size of the virtual object caused by the two atoms located at z = −z 0 before the lens is chosen as ∆x 0 . Considering f and f 1 , respectively, as the associated focal lengths to atom (1) and (2), the created image is demagnified by (∆x) a . This is represented by the following equation Given that the object size is obtained by the linewidth arising from the beam divergence would be estimated as where the beam divergence, is directly taken from the transverse velocity linewidth, ∆v x , and the most probable longitudinal velocity, v z . Inserting Eq. (31) into Eq. (30), one would obtain As a case in point, we refer to the example in sections II and III B. A transverse velocity linewidth of ∆v x = 0.0264 cm/s [see Fig. 6(c)] and a probable longitudinal velocity of v z = 1 cm/s result in θ = 0.0264 rad at t = 2 ms for the BEC falling from z i = 20 µm focusing to the center of lattice (ξ = 5.37) of a radius size of σ z = 10 µm. In such an instance, considering the associated focal length, f = 5.531 µm, the angular divergence broadening is estimated as (∆x) a = 0.146 µm.
Clearly, a longer time of BEC propagation gives the atoms more chance to interact within the cloud. Thus, to reduce ∆v x and consequently (∆x) a , one can choose to utilize a fall time for the focusing process.

V. GPE AND CLASSICAL TRAJECTORIES AGREEMENT
As stated thus far, the influence of the s-wave interaction between atoms appears as the longitudinal and transverse velocity distributions in the classical calculations. That being said, one can define the following equivalency where which can be rewritten as following using Eq.

VI. NUMERICAL INVESTIGATION OF FOCUSING
We now study in more detail the impact of the focusing potential aperture size (or the potential wavelength) as well as its radius on the deposited BEC structures. We consider a 87 Rb BEC trapped by ω x = 2π × 10 Hz, and ω y = ω z = 2π × 70 Hz, whose center-of-mass is at z 0 = 500 µm from the lattice center at z = 0. In the results presented, the freely propagating and interacting (a s =  100a 0 ) condensate is aimed to be focused optimally at z = 0 (ξ = 5.37). We initially consider a lattice radius size of σ z = 10 µm while three various potential wavelengths are employed, λ = 16λ D2 , 8λ D2 , and 4λ D2 . The results are extracted for a range of initial momentum kicks differing from p = 16 k to p = 128 k (where k is the wavenumber associated with the 87 Rb D 2 line defined by k = 2π/λ D2 ) is applied to the BEC. As a result, for an optimal focus, the lattice power and peak intensity are adjusted accordingly based on the BEC kinetic energy for any certain momentum kick. Figs. 10(a, b) display the associated outputs for the focused BEC linewidths and peak densities respectively. The FWHM in each numerical calculation is achieved by taking an average over the linewidth of the structure peaks whose height are greater than 1/e of the highest central peak.
Analyzing the results, firstly, one can conclude that imparting higher momentum kicks to the BEC reduces the focal spot sizes and enhances the peak densities for any potential wavelength. However, larger λ values necessitate higher potential powers and peak intensities resulting in the superior structure resolutions and heights (see the blue stars in Figs. 10(a, b) compared to the red triangles and green crosses). Furthermore, one can understand that the FWHM and peak density data points for λ = 16λ D2 , 8λ D2 , and 4λ D2 tend to a steady state at relatively large magnitudes of initial momentum kick [i.e. p ≥ 96 k, see Figs. 10(a, b)]. This implies the characteristic factors of a focused profile (including the resolution and peak density) become independent of the BEC momentum kick at relatively high longitudinal velocities. At this point, the profile resolution also becomes independent of the lattice wavelength while the focused peak density remains strongly dependent on the size of the momentum kick.
In Fig. 11(a, b), a fixed value of λ = 16λ D2 is considered while three different lattice radii, σ z = 10, 20 and 40 µm are selected. Since an expansion in the potential radius size leads to a decline in the power and peak intensity values, it is expected larger FWHMs will be obtained as well as lower peak densities for σ z = 40 µm (indicated by the green crosses) compared to those of σ z = 10 µm (indicated by the blue stars). As an instance, choosing σ z = 10, 20 and 40 µm would respectively result in (∆x) int GPE = 19.68, 23.27 and 48.83 nm for the resolutions, and 4.068 × 10 4 , 2.692 × 10 4 and 1.501 × 10 4 atoms/µm 2 for the peak densities at p = 128 k. Moreover, the profile characteristic factors tend to a steady state for every potential size at larger momentum kicks [i.e. p ≥ 96 k, see Figs. 11(a, b)].
Finally, the focused 87 Rb BEC profile is studied through a focusing lattice of λ = λ D2 . The output FWHM and peak density values are plotted as a function of potential power in Figs. 12(a, b)   (from P = 0.3 to 1.83 mW) regardless of an optimal power value for any momentum kick, P = P λD 2 (n k). In other words, the focused profiles are examined at z = 0 even when they are not at their optimal focus stage. According to Figs. 12(a, b), for each value of momentum kick, a rise in the focusing potential power results in an enhanced resolution as well as a higher peak density. For instance, setting P = 1.83 mW for p = 256 k produces a focused profile with (∆x) int GPE = 51.25 nm and a highest peak of 2975 atoms/µm 2 . However, increasing the value of beam radius (from σ z = 100 to 200 µm) broadens the corresponding profile resolutions creating shorter peaks, which is indicated in Figs. 13(a, b). In such a case, for P = 1.83 mW and p = 256 k, the resolution and highest peak density respectively reduce to (∆x) int GPE = 85.87 nm and 1726 atoms/µm 2 .
It is also noticeable that for any magnitude of momentum kick applied to the BEC, the focused profile characteristic factors become independent of the potential power at larger P values. As a result, in this case, there exists a threshold in power value for improving the profile resolution and peak density meaning that one may not expect to obtain a considerable improvement in focal spot sizes when using P > 2 mW [see Figs. 12(a, b) and Figs. 13(a, b)]. We note that the threshold range, P > 2 mW, is only the case for a lattice of λ = λ D2 as it could be reduced by a factor of ten to P 16λD 2 (128 k) = 0.39 mW for a lattice of λ = 16λ D2 , which could lead to a focused profile of (∆x) int GPE = 19.68 nm [see Fig. 11(a)].

VII. CONCLUSIONS
In this paper, we summarized the classical trajectories approach and its applications in focusing ultra-cold 87 Rb atoms. The structure of an appropriate focusing potential as well as its focal properties for an optimal focus were analyzed using the paraxial solution to the equations of atomic motion. Moreover, we derived a relation between the required value of potential power and the desired focus point on the focal plane along with employing a relevant example. The spherical aberration on broadening of focused structures was explored through numerical solutions to the classical equations. We then calculated the structure linewidth resulting from the diffraction contribution. We inferred that an increase in lattice aperture size and atomic longitudinal velocity or a reduction in the focal length results in a lower FWHM of diffraction as in this case atoms would have less chance to interfere together.
We then moved on to consider the influence of s-wave interactions between the atoms on focusing of 87 Rb condensate using the GPE. It is found that the resultant structure linewidths are significantly impacted by s-wave interactions. Moreover, we concluded that either reducing the lattice radius size or increasing the BEC initial kinetic energy could improve the characteristic factors (resolutions and peak densities) of a focused 87 Rb profile. Further to this, the contribution for both the BEC longitudinal and transverse velocity distributions in broadening the structure size was explored via the classical trajectories model. We found a reliable agreement between the GPE approach and a classical trajectories model which includes the interaction contribution.
Finally, conducting a variety of numerical simulations with different lattice radii and wavelengths, we noticed that narrower and higher density structures arise from relatively lower potential radius sizes and greater wavelengths. It was shown that nano-meter structures are achievable in principle. The example of this is a resultant structure of (∆x) int GPE = 19.68 nm via a focusing lattice of σ z = 10 µm and λ = 16λ D2 . Furthermore, it was inferred that applying higher momentum kicks necessitating greater optimal potential powers and peak intensities cause superior profile resolutions and peak fluxes. However, for a certain lattice wavelength, there is a threshold point in potential power where for larger values than this point, the focused profile characteristic factors are no longer dependent on power magnitudes tending to a steady state. It was observed that exploiting a relatively smaller lattice wavelength increases the threshold power whilst causing a destructive impact on the structure resolution and peak density. For instance, the approximate threshold power for a lattice of λ = λ D2 and λ = 16λ D2 is about P λD 2 = 1.83 mW and P 16λD 2 = 0.39 mW resulting in (∆x) GPE