Retroreflection and diffraction of a Bose–Einstein condensate by evanescent standing wave potential

The characteristic of the angular distributions of accelerated Bose–Einstein condensate (BEC) atoms incidence on the surface is designed using the mathematical modeling method. Here, we proposed the idea to study retroreflection and diffraction of a BEC from an evanescent standing wave potential (ESWP). The ESWP is formed by multiple reflections of the laser beam from the surface of the prism under the influence of gravity. After BEC’s reflection and diffraction, the so-called BEC’s density rainbow patterns develop due to the interference which depends on the surface structure which we model with the periodic decaying evanescent field. The interaction of accelerated bosonic atoms with a surface can help to demonstrate surface structures or to determine surface roughness, or to build future high spatial resolution and high sensitivity magnetic-field sensors in two-dimensional systems.


Scientific Reports
| (2020) 10:20674 | https://doi.org/10.1038/s41598-020-77597-8 www.nature.com/scientificreports/ produced by totally internally reflecting a laser beam at the surface of a refractive medium and then retroreflecting the light back along its original path as shown in Fig. 1. The evanescent field decreases exponentially in the direction normal to the surface and is modulated sinusoidal along the surface as depicted in Fig. 1 23 .
The underlying two-dimensional model for a BEC in the ESWP trap geometry is explained in "Theoretical model", where we also provide estimates for experimentally realistic parameters. A detailed study of the reflection and diffraction of the BEC wavepacket is provided in "Dynamics of the BEC". And, in "Conclusion" we give a brief summary and draw our conclusions. Finally, in "Method and numerical technique" we describe a method and numerical technique in detail.

Theoretical model
We start working out a two-dimensional model for the dynamics of the BEC in the ESWP. Assuming that both thermal and quantum fluctuations can be neglected, we arrive at a mean-field theory with the two-dimensional time-dependent Gross-Pitaevskii equation (GPE) 31 On the right-hand side the first term represents the kinetic energy of the atoms with mass m, the last term stands for their two-particle interaction, where its strength G = G 3D /(a y √ 2π) , here G 3D = N4πa 2 /m is related to the s-wave scattering length a, m denotes a mass of the bosons, and a y = ( /mω y ) 1/2 defines the width of the oscillator along the y-axis. For this pancake configuration of the trap potential, we have in mind ω y >> ω ⊥ , here ω y defines the frequency of the trap along the y−axis and the ω ⊥ defines the frequency of the radial component, which has to be small enough in order to ensure a quasi two-dimensional setup 31 . In comparison to a conventional transmission grating, the atomic reflection grating has to achieve two functions: reflection and diffraction. For this purpose, we proposed here an evanescent standing wave potential. A reflection optical potential with a decaying field along the z-direction and a periodic modulation along the x-axis in the presence of gravity. This kind of periodicity can be realized by overlapping two evanescent waves with opposite wave-vectors as shown in Fig. 1. The BEC atoms falling off on this optical potential will experience a reflection of the z-component of their momentum, this reflection is only possible when their kinetic energy does not increase the potential height, this condition helps BEC to avoid to penetrate into the potential. Furthermore, V(x, z) in Eq. (1) represents the underlying anharmonic potential energy of our system: Here, g represents the gravitational acceleration, and V 0 = Ŵ 3 I 0 /(8π 2 cδ 3 ) describes the strength of the evanescent field, where Ŵ is the natural line-width of the Cs atoms, defines the wavelength of the EW, I 0 denotes the peak intensity of the EW, and δ 3 corresponds to the detuning frequency of the hyperfine sub-level F = 3 of the Cs atom band. While, the dimensionless quantity η is related to the magnitude of the wave vector of the standing wave 32 , and 1/κ = /4π n 2 sin 2 β − 1 defines the decay length, where n represents the refractive index of the medium and β denotes the angle of incidence of the laser beam. When η → 0 , the leftover anharmonic potential (2) can be approximated around its minimum z min 0 = 1 κ ln( V 0 κ mg ) with the axial frequency ω z = √ gκ .
In order to have a concrete set-up in mind, we will use for our analysis the following parameter values which stem from the GOST experiment 20,30 . We denote the number N of Cs atoms, where the s-wave scattering length amounts to a = 440 a 0 with the Bohr radius a 0 . The inverse decay length amounts to κ = 6.67 × 10 6 m −1 as www.nature.com/scientificreports/ the EW is produced by a far-detuned laser with wavelength = 852 nm in a nearly round spot on the surface of the prism with reflective index n = 1.45 , and the angle of incidence of the laser light is β = 49.2 o . The natural line width of the Cs atom is Ŵ = 2π × 5.3 MHz , the detuning between the hyperfine sub-levels is δ 3 = 2π × 1 GHz , and the peak intensity of the EW amounts to I 0 = 9.6 × 10 7 W/m 2 , so the strength of the EW is of the order V 0 = 0.96 × k B K , where k B is Boltzmann's constant and the axial harmonic frequency amounts to ω ⊥ = 2π × 1.2 kHz.
In the following considerations we assume ω y = 2π × 12.0 kHz , which fulfills the quasi two-dimensional condition ω y >> ω ⊥ for the ESWP experiment. As the BEC does not penetrate very far into the repulsive ESWP, it is hardly influenced by the Van-der-Waals potential, which follows from the above values of the GOST experiment 30 . In order to do the numerical simulation, we use the dimensionless equation as here, we consider dimensional spatial coordinates z = κz and x = κx , the dimensionless time t = t(mg)/ κ and the two-particle dimensionless interaction strength G = N2 √ 2πãk/ã y , here, ã = aκ being a dimensionless s-wave scattering length, ã y = a y κ is the dimensionless oscillator length along the y-axis, and k = 2 κ 3 /(g m 2 ) defines the dimensionless kinetic energy constant. Additionally, we measure energies in gravitational energy mg/κ and dimensionless frequency describe as ω z = ω z κ/(mg) and dimensionless strength of the evanescent field represents as Ṽ 0 = κV 0 /(mg). By keeping in mind the last experimental values, we have the following dimensionless quantities: the dimensionless s-wave scattering length amounts to ã = 0.033 , the dimensionless optical decaying strength is given by Ṽ 0 = 906 , the dimensionless kinetic energy strength is k = 0.066 , the dimensionless radial frequency amounts to be ω r = 1.303 , and dimensionless two-particle interaction yields G = 0.086 N. We have normalized the wavefunction as |�(x,z)| 2 dxdz = 1 . To avoid complexity, we drop all parameter tildes for rest of the manuscript.

Dynamics of the BEC
Initially, we trapped a BEC in a two-dimensional trap at the top of the mirror as shown in Fig. 1a. Numerically, we achieved this process by letting the initial trapped potential as V (x, z) = x 2 + (z − z 0 ) 2 , here z 0 = 10 . At the time t = 0 , the initial trapping potential is switched off and condensate is released into the previously described standing evanescent potential Eq. (2) and the BEC expands ballistically (interacting gas) in potential. Here, numerically we avoid penetration of atoms into the prism by consider a very high potential for the negative z-axis i.e., V (z < 0) → ∞ . During the bouncing of the BEC, the atoms experience a weakly confining potential along the x-direction and relatively linear along the z-direction due to the gravitational field.
We demonstrate the diffraction patterns by using simple geometry. In our model, the potential plane is well approximated by a sinusoidal exponentially decaying field, which can be modeled here as a reflective grating for the BEC atoms. For quantum scattering the periodicity interval d leads to Brag scattering for constructive interference with the azimuthal angle ψ and the de Broglie wavelength db = ( /(2mE 0 )) 1/2 . Above equation describes the angular positions of diffraction spots ψ of order n. For small angle approximation sin(ψ) ≈ ψ , the Bragg condition is given by △ψ ≈ db /d , which has been verified experimental for the Fast Atom diffraction 33 . The azimuthal exit angle ψ is associated to the deflection angle θ in the detection plane via tan(ψ) = tan(φ in ) sin(θ) as shown in Fig. 1c, by using the small angle approximation sin(ψ) = sin(θ) sin(φ in ) the above Eq. (4) takes the form Here, db⊥ defines as a de Broglie wavelength with respect to the dynamics of projectiles in a plane normal to axial strings as shown in Fig. 1c. We notice from Eq. (5) that the azimuthal splitting of diffraction spots △ψ for a given de Broglie wavelength db is independent of the incident angle φ in . The diffraction spots △ψ helps to obtain the periodicity d of the surface potential. The external periodic potential frequency leads to a blazing effect similar to an optical grating, which enhances the intensity of finite diffraction order.
As we already discussed in "Introduction" that the quantum treatment for the retroreflection and diffraction of a Bose-Einstein condensate by standing evanescent wave does not permit analytical solutions, due to the anharmonicity and nonlinear interaction of BEC atoms. Therefore, we tackle this problem by using a numerical technique called time-splitting spectral method [34][35][36][37][38] , which perfectly depicts the description of the reflection and diffraction phenomenon. Therefore, although in this scientific report, we give up analytical approach, but we obtain a detailed insight into the dynamic process of the BEC. The evolution of the BEC probability density as it approaches the evanescent wave and then gets reflected as shown in Fig. 2a-l. In this special case, the amplitude of a standing wave is zero, therefore, the BEC start bouncing at the exponentially decaying potential. In Fig. 2, at the BEC density graph, additionally we plot a white line that mimics the periodic reflecting potential. The mean position of the BEC exhibits a periodic bouncing as shown in Fig. 2, which is natural as this is the clean case and ideally there are not any impurities present at the prism surface. As we can see that after letting η = 0 , our proposed model becomes a one-dimensional case as can be seen from Eq. (2), therefore to avoid ballistic expansion of the BEC cloud along the x-axis, we explicitly let a harmonic potential x 2 /2 along the x-axis. If we do not www.nature.com/scientificreports/ consider the harmonic potential along the x-axis, we notice that we lose all the BEC atoms in the x-axis, therefore to see the retroreflection of the BEC in this special case, we need to let harmonic potential along the x-axis. This way, we record bouncing of the BEC on a hard surface. We notice that the width ( σ x ) of the BEC along the x-axis is constant as predicted in Figs. 2 and 6b, this is quite obvious as we know that along the x-axis the potential is harmonic. We would like to mention that the width of the BEC along the x-axis is model by standard deviation. As we make a small increase in the dimensionless amplitude η = 0.01 of the standing wave as shown in Fig. 3a-l, the BEC temporal density starts to spread during the flight back as can be seen in Fig. 3d-f. This shows that the prism surface is not smooth anymore, therefore the BEC wavepacket spreads, by studying the spread of the BEC we can predict the shape of the surface structure similar to Surface Physics, where many experimentalists use the Fast Atom Diffraction technique to study surface patterns 33,39 . After the second reflection from the prism, atoms indicate interference maxima and minima of the temporal density of the BEC as shown in Fig. 3l at dimensionless time t = 60 . This interesting phenomenon can be explained, the interference patterns are  www.nature.com/scientificreports/ produced by overlapping the part of the BEC wavefunction that is still approaching to the maximum height with the part that is being diffracted previously. These maxima and minima of density wave-packets are quite visible when we further increase the dimensionless amplitude of the standing wave potential, let say η = 0.1 as shown in Fig. 4. In the panel Fig. 4c, the BEC wavefunction is interacting with the standing wave antinode at x = 0 . In this case, there is a circular wave-front formed around the standing wave antinodes, so some of the portion of the wavefunction is diffracted back earlier as compared to the other portion which still touching the minima of the standing wave. The superposition of the circular wave-fronts with each other generates maxima and minima of the temporal density, the temporal density gets weaker and blurred from the middle as we move far from the standing wave as predicted in Fig. 4f at dimensionless time t = 30 . As the BEC wave function further evolves in time, we can see a superposition of four contributions of the BEC wavepacket, as can be seen in the panel (k) of Fig. 4 and a rainbow kind of density structure appears. The rainbow pattern has all the information of the surface structure, thus we see four maxima of the density patterns, therefore we can claim that there are four periodic structures are present at the prism surface.  www.nature.com/scientificreports/ Interestingly, the BEC wavepacket starts exhibiting the typical Fresnel or near-field features of such functions, these features are very similar to the light when it passes from a diffraction grating. As we increase further the standing wave amplitude to η = 0.9 as shown in Fig. 5a-l, the circular wave-fronts of BEC wavepacket reaped around the standing wave antinodes as shown in Fig. 5c. We would like to mention here, our mimicking periodic one-dimensional white potential is not describing the exact Physics of the ESWP as gravity starts playing a very important role, therefore, it is getting hard to mimic a two-dimensional potential (Fig. 1b) as a onedimensional case. It can be seen from the density rainbow structure as presented in Fig. 5e that end-part of the BEC wavepacket take more time to reflect back as compared to the middle part. The wavefront of the BEC could not find the perfect plan surface instead the surface is quite rough. So by comparing the clean surface case with the dirty periodic surface case, someone can provide a very good estimate of the periodic surface structure. We note that the mean position of the BEC wavepacket along the z-axis is periodic for all the different dimensionless periodic amplitudes as shown in Fig. 6a, however the mean-position of the BEC decrease with every bounce which means the loss of energy as depicted in Fig. 6a. So a natural question arises where does this energy go, in this closed system, we notice that the total energy of the system is conserved however, this decrease in height is happened due to the spread of the wavepacket in the x-axis. We observe that the standard-deviation (width) of the BEC wavepacket along the x-axis increases linearly with time and the growth of the linear curve increase with the increase of the dimensional periodic strength η as shown in Fig. 6b. For a special case, η = 0 the width of the BEC does not increase as the BEC is bounded in the x-axis due to the harmonic confinement, therefore we do not observe any decrease in mean-position of the BEC wavepacket along the z-axis. For η = 0.9 , the BEC wavepacket losses its coherence due to the loss of a definite phase as it bounces on a quite high periodic potential surface and on every bounce a relative phase is added. Therefore, as the periodic potential strength increases the Quantum decoherence increases and we see quite dull interference patterns as shown in Fig. 5l, thus for high values of η , we can not track any information from these interference patterns.

Conclusion
In conclusion, we have demonstrated a simple retroreflection and diffraction technique for ultra-cold atoms, where the diffraction patterns can be controlled with the experimental tuneable parameters. We note that the spreading of BEC in the vertical direction is quite arduous. Atoms that originated at the top have a longer bouncing time-period than those started at the center. This effect is quite small during the first few bounces, however, it increases as the number of bounces started to increase. The bright and dark interference fringes can be seen, however, the coherence of the BEC started to lose after few bounces as in very bounce BEC atoms get phase due to path traveling along the z-axis. The prospect of exploitation of the effect of rainbow scattering of the BEC wavepacket to examine the structural features of the surface is examined. Fast Atom Diffraction from surfaces is one of the standard methods for solid surface structures near to the room temperature. However, our technique can be used to find out the surface structures at ultra-cold temperature. Our technique will open new prospects in the field of BEC diffraction and reflection from the prism surface. Nowadays magnetic-field sensors 40 are efficient of making measurements with both good field sensitivity and high spatial resolution. For instance, with low sensitivity magnetic force microscopy 41 permits investigation of magnetic structures with a spatial resolution in the nanometre range, although atomic magnetometers and SQUIDs facilitate to investigate extremely sensitive magnetic-field measurements, but at low resolution. In this respect, Wildermuth et al. investigated one-dimensional Bose-Einstein condensates in a microscopic field-imaging technique that combines high spatial resolution (within 3 micrometres) with high field sensitivity (300 picotesla) 42 . Our proposed model findings can help experimentalists to build future high spatial resolution and high sensitivity magnetic-field sensors in two-dimensional systems.

Method and numerical technique
We derive the Gross-Pitaevskii equation by using the Hamilton principle of least action with the action functional The dynamics of two-dimensional condensate can be observed when the longitudinal trap frequency ω y is much higher than the radial trap frequency ω r . If this condition holds then the condensate wave function can be factorized where φ(y) = 1 π 1/4 a 1/2 y e −y 2 /(2a 2 y ) . By substituting Eq. (8) into Eq. (6) and integrate the over y-axis, the effective two-dimensional Lagrangian density To determine the two-dimensional time-dependent Gross-Pitaevskii equation, we use the Euler-Lagrangian equation where, r = (x, z) . By using the two-dimensional Lagrangian density Eq. (9) and after some algebra, the twodimensional time-dependent Gross-Pitaevskii equation reads (1) 43 .
In order to determine the dynamics of the condensate wave function, we solve the two-dimensional timedependent Gross-Pitaevskii equation (1) in imaginary time numerically by using the split-operator method [34][35][36][37][38] . It is worthwhile to mention here that we use the imaginary time simulation technique that helps us to find a suitable ground state wave function, for any random initial guess, for a given interaction strength. The ground state wave function then served as an initial condition for the rest of the dynamic cases. In our case, we use initially Gaussian as a guess which by simulating in imaginary time results in a ground state for a specific interaction strength. For numerical simulations, we perform discretization of the dimensionless GP Eq. (3), with space step △ x =△ z = 0.083 and the time step dt = 0.001.