Multimode Collective Atomic Recoil Lasing in Free Space

Abstract

Cold atomic clouds in collective atomic recoil lasing are usually confined by an optical cavity, which forces the light-scattering to befall in the mode fixed by the resonator. Here we consider the system to be in free space, which leads into a vacuum multimode collective scattering. We show that the presence of an optical cavity is not always necessary to achieve coherent collective emission by the atomic ensemble and that a preferred scattering path arises along the major axis of the atomic cloud. We derive a full vectorial model for multimode collective atomic recoil lasing in free space. Such a model consists of multi-particle equations capable of describing the motion of each atom in a 2D/3D cloud. These equations are numerically solved by means of molecular dynamic algorithms, usually employed in other scientific fields. The numerical results show that both atomic density and collective scattering patterns are applicable to the cloud’s orientation and shape and to the polarization of the incident light.

1. Introduction

In the original configuration of the Collective Atomic Recoil Lasing (CARL, [1,2]) an ensemble of N initially randomly distributed atoms set in the arm of a ring cavity scatters photons in the reverse direction of an incident laser injected along the cavity axis [3,4]. The process arises due to a collective instability, where each atom is affected by the optical field scattered by the other atoms. After the production of ultra-cold atomic samples, a similar phenomenon was first demonstrated using a cigar-shaped Bose–Einstein Condensate (BEC) [5] and later using a cold, thermal gas [6], named Superradiant Rayleigh Scattering. In [5], superradiant scattered light was observed to propagate along the major axis of the atomic cloud, simultaneous with the formation of a matter-wave momentum grating in the cloud. Some features of this phenomenon have been described by single-mode/mean-field models similar to that of the (CARL) [7,8,9,10,11,12]. These mean-field models are appropriate in certain specific cases, where there is a well-defined propagation axis and consequently, to a good approximation, a single spatial mode, e.g., in a single-mode cavity or in a highly elongated sample where its major axis defines an ‘end-fire mode’ [7] which dominates the emission direction. In general, however, for arbitrary shapes of atomic ensembles, many spatial modes are involved simultaneously in the collective scattering process.

In a recent paper, we presented a model able to describe collective scattering of light by an atomic cloud into multiple vacuum modes, when this is irradiated by an optical field [13]. For the sake of simplicity, in the model a scalar field was assumed, neglecting polarization effects and near-field terms in the scattered light. The model’s equations, which follows the particles’ trajectories with the sole dependency of their positions, were implemented in a molecular dynamics (MD) code. The response of a two-dimensional cold atomic gas was proven to be sensitive to the cloud’s shape and to its major axis’ orientation with respect to the propagation axis of the incident field. Any change in these features triggers distinct pattern formation in the cloud, due to the atomic redistribution caused by different scattering events. The atomic rearrangements vary from unidimensional gratings to organizations in two dimensions, occurring when the cloud’s major axis is aligned or not with the propagation axis of the pump field, respectively. The scattering was proven to befall across the major axis of the cloud with no need for an optical resonator. Patterns similar to those observed in [5] with a BEC in the momentum space have been shown in [13], but in the real space and with thermal atoms. In addition, the capabilities of the code were demonstrated by showing the formation of a density grating when the study of collective light-scattering was extended to a three-dimension atomic cloud.

Specifically, three main shapes and orientations were investigated in [13] with different outcomes. On the one hand, when the cloud is elongated perpendicularly to the incident field, the atomic scattering occurs mainly across the perpendicular axis (the longest path), which combined with the pump translates in a 2D density pattern formation. On the other hand, the grating formations observed in both 2D and 3D profiles, elongated along the pump propagation axis, exhibit analogies with the ones observed in either the collective atomic recoil lasing (CARL) [1], which enhances the backscattered light, or in the free-electron laser (FEL) [14]. In addition, a third case with a spherical cloud was investigated. With this configuration, any density grating was expected due to the lack of any longest scattering path. However, the peculiarity of this case is that an electrostrictive force, generated by collective scattering, slightly deforms the initially spherical shape. Then, the cloud is stretched along the pump’s propagation direction and compressed along the perpendicular direction, by outlining two preferred scattering paths along its edges. The resulting two-dimensional density grating differs from the one observed in a perpendicular elongated cloud, since the light is scattered at an angle larger than ${90}^{\circ }$ with respect to the pump’s incident axis.

In the present paper, we extend the former scalar model of [13] to a full vectorial model of light. The polarization effects are shown to be relevant in the scattering process: since the radiation electric field is transverse, propagation along the pump’s polarization direction is not allowed. The study is substantially extended by considering 3D collective scattering with different distributions and orientations of the atomic cloud. Moreover, we investigate the three-dimensional intensity spatial profiles, for a better understanding of the grating formation for each configuration. The vectorial model, likewise its scalar version, has been derived from a multimode theory, where the vacuum radiation modes are adiabatically eliminated. The resulting N coupled equations, where each atom of the cloud is subdued to the radiation force produced by the other $N-1$ atoms, are now also accounting for the polarization effects of the incident and scattered fields. These new equations, although more complex due to the presence of addition short-range terms, are also eligible to be implemented in MD codes. The algorithm used in the current letter has been adapted to work with the rather new open-source programming language developed by the MIT, Julia. The former algorithm employed in the scalar study, Verlet velocity, has been substituted with its close relative leapfrog algorithm, or in this case the variant “Verlet leapfrog”, included in one of Julia’s packages [15]. The vectorial model is introduced in Section 2, together with a brief description of the employed MD algorithm, the chosen inter-particle interaction, and the rescaling of the motion equations. These results, for the 2D and 3D geometries, are discussed in Section 3.1 and Section 3.2, respectively. Further details of the analytical derivation of the model can be found in the Appendix A.

2. Vectorial Model of Collective Scattering

2.1. Motion Equations

Following the line applied in our previous work [13], we consider a gas with N two-level atoms driven by a laser field propagating along the z axis, but in this case with linear polarization along the x or y axis. The pump field has a frequency ${\omega }_0=c{k}_0$, with wavevector ${\mathit{k}}_0=k_0\widehat{\mathit{z}}$ where $k_0=2\pi /{\lambda }_0$, and a Rabi frequency ${\mathsf{\Omega }}_0=d{E}_0/\hslash$, where $E_0$ is the electric field and d the atomic dipole. A far detuning ${\Delta }_0={\omega }_0-{\omega }_a\gg \mathsf{\Gamma }$ is chosen between ${\omega }_0$ and the atomic frequency ${\omega }_a$, where $\mathsf{\Gamma }=d^2{k}_0^3/2\pi {ϵ}_0\hslash$ is the transition linewidth. In such far-detuned limit, and when a dilute gas is deemed, multiple scattering of photons can be neglected, thus reducing the scattering process in the cloud to single-scattering. Multiple scattering is linked to large optical thickness $b(\Delta )=b_0/[1+{(4\Delta /\mathsf{\Gamma })}^2]$, where $b_0=aN/{\left(k_{0}R\right)}^2$ is the on-resonance optical thickness [16], with R as the transverse radius of the atomic cloud and where $a=2$ is for a scalar field and $a=3$ or a polarized pump. Consequently, each vacuum mode $\mathit{k}$ is populated when an incident photon with mode ${\mathit{k}}_0$, is scattered into the vacuum mode by an atom of the cloud. The atom experiences a recoil momentum $\mathit{p}=\hslash ({\mathit{k}}_0-\mathit{k})$ every time it scatters a photon. The scattered optical field with mode $\mathit{k}$ interferes with the incident field with mode ${\mathit{k}}_0$, generating a dipole force proportional to the atomic recoil momentum.

By summing over all the radiated modes, we obtain a set of N coupled equations describing how each jth atom endures the radiation force exerted by each of the other atoms. When the polarization vector of the optical field is taken into account, the equations for the position ${\mathit{r}}_j$ and the momentum ${\mathit{p}}_j$ of each jth atom is (see Appendix A):

$${\dot{\mathit{r}}}_j=\frac{{\mathit{p}}_j}{m},$$

(1)

$$\begin{array}{cc}\hfill & {\dot{\mathit{p}}}_j=\mathsf{\Gamma }\hslash k_0{\left(\frac{{\mathsf{\Omega }}_0}{2{\Delta }_0}\right)}^2\sum _{m\ne j}\{\left(\widehat{\mathit{z}}-{\widehat{\mathit{r}}}_{jm}\right){\mathsf{\Pi }}_{jm}\frac{sin\left[k_0(r_{jm}-z_{jm})\right]}{k_0{r}_{jm}}+\left[\widehat{\mathit{z}}{Q}_{jm}-{\widehat{\mathit{r}}}_{jm}({\mathsf{\Pi }}_{jm}+Q_{jm})\right]\frac{cos\left[k_0(r_{jm}-z_{jm})\right]}{{\left(k_0{r}_{jm}\right)}^2}\hfill \\ \hfill & -\left.\left(\widehat{\mathit{z}}-3{\widehat{\mathit{r}}}_{jm}\right)Q_{jm}\frac{sin\left[k_0(r_{jm}-z_{jm})\right]}{{\left(k_0{r}_{jm}\right)}^3}+3{\widehat{\mathit{r}}}_{jm}{Q}_{jm}\frac{cos\left[k_0(r_{jm}-z_{jm})\right]}{{\left(k_0{r}_{jm}\right)}^4}\right\},\hfill \end{array}$$

(2)

where

$${\mathsf{\Pi }}_{jm}=\sum _{\alpha ,\beta }[{\delta }_{\alpha \beta }-{\left({\widehat{r}}_{jm}\right)}_{\alpha }{\left({\widehat{r}}_{jm}\right)}_{\beta }]{\widehat{e}}_{0\alpha }{\widehat{e}}_{0\beta }^{*}=1-({\widehat{\mathbf{r}}}_{jm}·{\widehat{\mathbf{e}}}_0)({\widehat{\mathbf{r}}}_{jm}·{\widehat{\mathbf{e}}}_0^{*}),$$

(3)

$$Q_{jm}=\sum _{\alpha ,\beta }[{\delta }_{\alpha \beta }-3{\left({\widehat{r}}_{jm}\right)}_{\alpha }{\left({\widehat{r}}_{jm}\right)}_{\beta }]{\widehat{e}}_{0\alpha }{\widehat{e}}_{0\beta }^{*}=1-3({\widehat{\mathbf{r}}}_{jm}·{\widehat{\mathbf{e}}}_0)({\widehat{\mathbf{r}}}_{jm}·{\widehat{\mathbf{e}}}_0^{*}).$$

(4)

where ${\widehat{e}}_{0\alpha }$ and ${\widehat{e}}_{0\beta }$ are the transverse components of the polarization vector ${\widehat{\mathbf{e}}}_0$ of the pump. The atomic mass is m and ${\mathit{r}}_{jm}={\mathit{r}}_j-{\mathit{r}}_m$, with ${\widehat{\mathbf{r}}}_{jm}={\mathit{r}}_{jm}/r_{jm}$ (where $m\ne j$).

The oscillating force represented by Equation (2) couples each atom j with all the other ones along two direction components: one coinciding with the external field $\widehat{\mathit{z}}$ and one connecting atom j to any other atom m, ${\widehat{\mathbf{r}}}_{jm}$. The force has a finite range, with four terms decreasing with the atom separation as $1/{\left(k_0{r}_{jm}\right)}^n$, within the range $1\le n\le 4$. The complete derivation of these equations is presented in the Appendix A. In the derivation, the population of the excited state has been neglected, due to the large detuning or weak field (linear optic regime). Furthermore, by assuming $\mathsf{\Gamma }\gg {\omega }_{rec}$ (where ${\omega }_{rec}=\hslash k_0^2/2m$ is the recoil frequency), the internal degrees of freedom of the atoms are adiabatically eliminated. Finally, the radiated field is also adiabatically eliminated by using the Markov approximation, with the resulting Equation (2).

2.2. Bunching Equation

The “bunching” factor is a reminiscence of the free-electron laser [14] and, since CARL is labeled as its atomic version, it is an important parameter to quantitatively study the density grating formation. This parameter allows measurement of how strong the periodic modulation of the atomic density along a certain direction is. Its value ranges from zero, meaning that the atoms are uniformly distributed, to one, when the atoms are maximally organized in bunches with period $\mathsf{\Lambda }=2\pi /|{\mathit{k}}_0-\mathit{k}|$. In the backward scattering $\mathit{k}=-{\mathit{k}}_0$ that period is $\mathsf{\Lambda }=\pi /k_0={\lambda }_0/2$. The bunching factor is described by the following expression:

$$M(\mathit{k},t)=\frac{1}{N}\sum _{j=1}^N{e}^{i({\mathit{k}}_0-\mathit{k})·{\mathit{r}}_j\left(t\right)},$$

(5)

where ${\mathit{r}}_j\left(t\right)$ is the position of a particle at a certain time and ${\mathit{k}}_0$ and $\mathit{k}$ are, in that order, the wave vector of the incident field and the scattered modes. The scattered intensity in the direction $\mathit{k}$ is

$$I_s\left(\mathit{k}\right)=I_1{N}^2{sin}^2\psi {\left|M(\mathit{k},t)\right|}^2,$$

(6)

with $I_1=(\hslash {\omega }_0\mathsf{\Gamma }/16\pi r^2){({\mathsf{\Omega }}_0/2{\Delta }_0)}^2$ as the single-atom Rayleigh scattering intensity and $\psi$ as the angle between the pump polarization axis and the scattering direction. The derivation of Equation (6) can be followed in Appendix A.3 of Appendix A.

2.2.1. Leapfrog Algorithm

It is unclear who introduced the leapfrog algorithm, but some early references to this method can be found in the first volume of The Feynman Lectures on Physics [17] from the 1960s. Subsequently, it was popularized in the 1980s by the hands of Hockney and Eastwood [18]. In the leapfrog technique, the velocity and the position are not calculated in the same time-step, there is a half-step lag between them. The study of the evolution of an all-atom system starts by setting the initial position for each particle belonging to that system. Subsequently, the velocity is calculated (leaps) over the position, then the velocity leap over the position, and so on until the desired time steps are acquired. Check Figure 1 for a graphical summary.

Usually, to reach a numerical solution with this method, the equations need to be transformed into iterative expressions that allow calculation of the position and the velocity of each particle at every time-step. This is performed by substituting the position and velocity within those expressions with $\alpha =({\alpha }_{n+1}+{\alpha }_n)/2$ and $v_{\alpha }=(v_{\alpha ,n+1/2}+v_{\alpha ,n-1/2})/2$, respectively, and their derivatives with $\dot{\alpha }=({\alpha }_{n+1}-{\alpha }_n)/\Delta t$ and $\dot{v_{\alpha }}=(v_{\alpha ,n+1/2}-v_{\alpha ,n-1/2})/\Delta t$, apiece, being $\alpha =x,y,z$. The final expressions are rearranged to calculate ${\alpha }_{n+1}$ and $v_{\alpha ,n+1/2}$, thus producing two iterative equations suitable for the algorithm. However, as regarded earlier, Julia gives access to packages of predefined functions of symplectic algorithm, which includes the hybrid “leapfrog Verlet” algorithm [15] among others. This fact represents a great advantage when using Julia, because the coding process can be greatly accelerated, since there is no need for adapting the differential equations into iterative expressions.

The main advantage of employing the symplectic leapfrog algorithm is that it is less expensive, computationally speaking, than other numerical methods. It allows calculation of the drift of the particles within a system with high level of accuracy, and with very little compromise in the general picture of the energy conservation. In terms of memory usage, this numerical method requires less storage than the well-known RK4 method from the Runge–Kutta family, whose introduction dates previous to the computation times (between the 1890s and the early 1900s).

2.2.2. Minimum Inter-Particle Distance Definition

To avoid the complex interaction arising when the atoms get too close to each other, a short-range repulsive potential needs to be implemented in the system dynamics. Traditionally, the common option is to employ the Lennard-Jones pair potential [19]. This simple potential can accurately characterize the repulsive-attractive interactions caused by the van der Waals forces between two neutral atoms. When it is applied in the equations describing CARL, it allows the counteraction of the singularities arising from pairs of atoms getting too close to each other.

In the present study, instead of calculating the short-range potential between every two atoms at each time-step, we have adopted another less orthodox option. The idea is to add a constant minimum distance between a pair of atoms to soften the strong repulsion between them, making the computation calculation lighter and faster. The concept comes from astrophysics, where it is employed in the numerical resolution of the gravitational force in stellar systems, e.g., galaxies or star clusters. Such astrophysical systems are not too different from an atomic cloud computationally speaking. Therefore, we decided to implement the pioneer Plummer model [20], whose original idea was to avoid the singularity points located in black holes in systems of galaxies. This idea is simplified using a similar concept to avoid the singularity-like reaction generated by the terms referring to shorter distance ($1/r_{jm}^n$ with $n\ge 2$) in Equation (2). A limit or cut-off $\xi$ is defined, as the numerical scattering length, and added into the inter-particle distance

$$r_{jm}\to \sqrt{r_{jm}^2+{\xi }^2}$$

(7)

to force the existence of a minimum value, hence avoiding the smaller values that can fragment the system.

The parameter $\xi$ barely affects the behavior of the particles of the system when they are separated by large distances and it recreates the repulsive effect generated due to van der Waals forces. In reality, whenever two particles encounter each other, they experience an overlapped crossing, which can be explained as a virtual elastic collision.

2.2.3. Scaling Model’s Equations

The N coupled equations of motion and the bunching need to be scaled to avoid increasing the numerical error. This can be achieved by working with position, momentum, and time dimensionless variables with smaller absolute values. To have a fair comparison with the previous scalar study, we adopt the change of variables proposed in Ref. [13]. The position vector for every atom is scaled with the wave vector of the incident light, becoming ${\mathit{r}}^{\prime }=k_0\mathit{r}$. The momentum variable is redefined as ${\mathit{p}}^{\prime }=\mathit{p}/\hslash k_0$, where $\hslash k_0$ is the momentum photon. The time is redefined as $t^{\prime }={\omega }_{r}t$, where ${\omega }_r=\hslash k_0^2/2m$ is the recoil frequency. All variables with a “prime” notation will be from now on refereed as variables without prime, i.e., $r^{\prime }$ will become r, but still in representation of $k_{0}r$. With these units, the force on the jth atom is proportional to the constant

$$A=\frac{\mathsf{\Gamma }}{{\omega }_r}{\left(\frac{{\mathsf{\Omega }}_0}{2{\Delta }_0}\right)}^2.$$

(8)

The distance between two atoms $r_{jm}$, has actually been further rescaled employing the singularity-avoiding parameter to $r_{jm}\to {\left(r_{jm}^2+{\xi \prime }^2\right)}^{1/2}$; the cut-off parameter also being rescaled as ${\xi }^{\prime }=k_0\xi$ (later simply regarded as $\xi$).

Even with a large detuning, it is possible to make an extra approximation for the wave vectors of the field and the atomic transition, $\left|\mathit{k}\right|\approx |{\mathit{k}}_0|$. By adding such approximation to the one already applied to the atoms’ positions, assuming ${\mathit{k}}_0$ directed along the z axis and writing the wave vector $\mathit{k}$ in polar coordinates, $\mathit{k}=k_0(sin\theta cos\varphi ,sin\theta sin\varphi ,cos\theta )$, the bunching in Equation (5) becomes:

$$M\left(t\right)=\frac{1}{N}\sum _{j=1}^N{e}^{-i[(sin\theta cos\phi )x_j+(sin\theta sin\phi )y_j-(1-cos\theta )z_j]},$$

(9)

where the dimensionless position variable for each atoms ${\mathit{r}}_j=(x_j,y_j,z_j)$ carries the implicit dependency on time. The bunching for the two-dimensional case can be obtained by setting $\phi =0$.

2.2.4. Parameters and Settings Applied in the Results

Real atomic clouds evolve in three-dimensional free space, but in the simulations they are also considered to be systems of two dimensions. The more realistic representations contain $N={10}^4$ and the bi-dimensional examples $N=5\times {10}^3$. Making reference to a gas of ⁸⁷Rb, we choose the following parameters adopted in the simulations: the D2-line transition at ${\lambda }_0=780$ nm has a decay rate $\mathsf{\Gamma }=3.8\times {10}^7{\mathrm{s}}^{-1}$ and the recoil frequency is ${\omega }_r=\hslash k_0^2/2\mathrm{m}\approx 2.3\times {10}^4{\mathrm{s}}^{-1}$. Moreover, we set $A=1$, giving the relationship between the detuning and the pump field as ${\mathsf{\Omega }}_0/\left|{\Delta }_0\right|\approx 1/20$.

A main distinction has been executed while setting the value of the singularity-avoiding parameter: in the vectorial model, due to the presence of higher order short-range terms, we set $\xi =1$, whereas the value employed for the scalar model in [13] was $\xi ={10}^{-2}$. When the additional short-range terms in Equation (2) are manually suppressed (forcing them to be zero) the value used for the vectorial case can approximate the one used for the scalar one. When a smaller value for this parameter is chosen, the cloud experiences a slow but constant expansion, due to a large value attained by the force, caused by the small distance existing between atoms.

In the same manner, we delimit the simulation step for the time evolution, which has been fixed to $\delta t={10}^{-3}{{\omega }_r}^{-1}$. The total simulation time has been set to $t=0.3{{\omega }_r}^{-1}$, yielding a total of 300 steps. The step used in [13] is roughly 6 times smaller, but using such time interval simply increases the simulation time with no visual impact in the results.

3. Results

There are several variables that can be adjusted to obtain different results when keeping the system parameters fixed, such as the pump frequency or the detuning, among others. For instance, we can talk about the space dimensions of the gas (2D or 3D); the cloud’s orientation with respect to the incident laser beam, which can accommodate different directions for the polarization vector; the shape of the cloud, being elliptical, which includes many degrees of elongation, or circular; and the density distribution of the atoms in the cloud (e.g., Gaussian or uniform). To give some hierarchy to such a variety, the central concept that will allow the characterization of such a complex system will be the space dimensions of the systems: 2D and 3D. The same numerical random seed has been introduced when defining the initial positions for the atoms inside both the flat and the volumetric clouds. The incident field’s polarization is linear, perpendicular to its propagation axis (z axis) and parallel to x or y axis. Restricting the number of system configurations available even further, only the off-axis scattering (elongated 2D/3D clouds with their major axis perpendicular to the incident field propagation direction) is here investigated. The analysis of this configuration allows visualization of the effects derived from having the polarization field parallel or perpendicular to the major axis of the cloud.

3.1. Polarized Light-Scattering from a 2D Atomic Cloud

3.1.1. Pump Polarization Parallel to the Major Axis of an Elliptical Cloud: Backscattering and 1D Grating Formation

When the cloud stretches along the x axis, the atoms experiences a total momentum change of $\hslash k_0(\widehat{\mathit{z}}\pm \widehat{\mathit{x}})$ [13], i.e., depending whether the photon is emitted downward (positive sign) or upward (negative sign). When the fact that light cannot be scattered along the direction of its polarization is added to this calculated scattering event, the only possible scattering direction for the current case is the z axis. This suppresses the $\widehat{x}$ component of the total recoil momentum. This suppression has obvious implications both in the grating generated in the cloud and the bunching profile, due to the combined action of the external and scattered fields.

The results shown in Figure 2 exhibit remarkable changes in the grating pattern generated in the cold vapor, in comparison with the scalar case, caused by the scattering restriction imposed by the polarization vector. The polarization blocks the scattering across the transverse direction (x axis), leaving the light with the only choice but to scatter along the propagation axis of the pump field (z axis). The initial distribution, shown in panel (a), is modified into a pure 1D grating pattern depicted in plot (b) of this figure. The grating is a consequence of the two only possible atomic momentum changes: some atoms remain in place, with a net kick close to zero, while others endure a double momentum kick, with a momentum variation of $2\hslash k_0$.

The bunching profiles shown in the polar plot of Figure 2d is consistent with the grating patterns appearing in panel (b). The intensity pattern manifests more similarities with the one generated in the scalar case with an elliptical cloud parallel to the propagation axis, than with the analogous case that presents the same cloud orientation: only the lobes representing the forward and backwards scattering direction are visible. The backscattering lobe is somewhat distorted and not as clear as in the scalar model, although it preserves some symmetry. The main cause for a deformation of this kind may come from the dipole-dipole force which starts to push the atoms towards the right of the plot. The force pushes with greater strength the atoms located close to the origin at $x=0$, where the optical thickness is larger, displacing them towards positive value of z faster than the ones located at the tips of the cloud (with larger $|k_{0}x|$ value). Summarizing, the radiation force transforms the expected straight strips of bunches of atoms, into arcs. These are equidistantly separated with a periodic space between them of $\approx \lambda /2$.

3.1.2. Pump Polarization Perpendicular to the Major Axis of An Elliptical Cloud: Off-Axis Scattering

We repeat the previous simulation, but now with the pump field polarization directed along the y axis. The new system does not experience any propagation restriction, because the polarization is perpendicular to the plane containing the atomic cloud. Therefore, a bi-dimensional grating similar to that obtained with the scalar model is expected to be observed for this system.

The grating portrayed in Figure 3b is not very different than the one obtained for the scalar system (see [13]). As predicted, pump polarization along the y axis does not affect the dynamics in the $(x,z)$ plane. The slight differences observed between the vectorial and scalar models arise from the short-range interaction terms. In fact, the force equation has now some additional short-range terms ($1/r_{jm}^n$ with $n>2$), and the terms that are typified in both models ($1/r_{jm}^n$ with $n\le 2$), are not identical. Moreover, the cut-off parameter has been augmented to $\xi =1$, instead of the $\xi ={10}^{-2}$ from the scalar study. The adjustment of this parameter does not affect the validity of the results, but its alteration certainly introduces a slight modification in the way particles interact with each other.

Analyzing the bunching factor in plot (d) of the same figure, it is possible to confirm that the system eventually responds as in the scalar case, providing a similar 2D pattern: two large lobes appear along the positive and negative directions of the longest scattering path (x axis). An earlier time snapshot is presented in Figure 3a,c, showing the effects of a non-negligible thickness that allows the generation of a transient 1D grating. After the preliminary formation of a longitudinal grating, the scattering direction is defined along the main axis of the ellipsoidal cloud (i.e., the x axis), hence reshaping the 1D grating into a 2D grating.

In both the simulations presented so far, we observe some atoms slowly drifting outward from the center of the ellipse. For instance, observing Figure 3, it can be identified how the cloud in (b) has expanded its area in comparison with that in (a). This is due to the ’collisions’ between the atoms: when two atoms become to close each other, the short-range terms in the force give them a kick, which is numerically controlled by the cut-off parameter, set to $\xi =1$. To control numerically the stability of the simulation, we have assumed a cut-off parameter $\xi$ larger than the one adopted for the scalar case in ref. [13] ($\xi ={10}^{-2}$). If this value were to be applied to the present simulations with the vectorial model, the cloud would suffer a stronger and faster repulsion, making the 2D grating impossible to appear.

3.2. Scattering from 3D Cloud of a Linearly Polarized Field

3.2.1. Pump Field Perpendicular to Ellipsoid’s Major Axis with In-Plane Polarization

The first three-dimensional case, like in the 2D system, is investigated with the pump linearly polarized along the x axis, which forbids the scattering along that axis. Since in 3D space there are still two more propagation directions that are allowed, it is expected that the light forms a series of atomic rods (composed of several hundreds of atoms each) along the axis of the polarization vector. Again, the result is due to the combination of the photons transferred from the pump to the scattered modes along the y and z axes. Such scenario can be discerned in Figure 4b, where these atomic rods are not completely straight, displaying an arched shape, as in the 2D system, due to the effect of an electrostrictive force [13]. Since the optical thickness is larger on the central ($y,z$) plane (located at $x=0$), the central area endures a greater force from the pump field than the tips of the cloud, which gives to the atomic rods their characteristic arched shape. The complete suppression of the bunching along the polarization axis, can be observed in Figure 4c, where the lobes along the polarization directions (x axis) are completely missing.

3.2.2. Pump Field Perpendicular to Ellipsoid’s Major Axis with Out-of-Plane Polarization

The external optical field is here considered with the polarization vector rotated ${90}^{\circ }$ with respect to the previous case, i.e., along the y axis. As the scattering is forbidden along the pump polarization y axis, tilted planes of particles are now visible in two chosen time instants: faintly drawn in Figure 5a and well defined in panel (b) of the same figure. The snapshot in (a) shows a grating composed of these tilted planes in formation. The planes with a positive slope are more visible because the atomic recoil is first directed downwards (towards negatives value of the x axis), which is confirmed by the upwards lobe of the bunching factor in panel (c). The titled planes with a negatives slope are also represented in a second snapshot (b) at a later time. Here they are dimly traced due to the latter appearance of the downwards scattering, which is represented by the downward lobe of the bunching factor in panel (d).

In addition, we do not observe any lobe in the bunching factor along the y axis either, as expected, since it is the pump polarization direction. The asymmetry observed in the lobes of bunching factor along the x-axis—see Figure 5c—is due the random initial conditions of the numerical simulations. This asymmetry is less evident at later time instant—Figure 5d—when the collective scattering process becomes dominant.

3.2.3. Pump Field Parallel to Ellipsoid’s Major Axis

We display the pattern emerging in a characteristic “cigar-shaped” cloud of cold atoms, which elongates along the propagation axis of the pump. The cloud evolved in the presence of a linearly polarized pump field, whose initial distribution and bunching intensity profile are shown in Figure 6, is shown in Figure 7. The plots in this figure are at the same time for a cloud with same shape and initial conditions, with snapshots (a) and (c) for a pump polarized along x and snapshots (b) and (c) for a pump polarized along y.

A similar behavior than the one observed in ref. [13] with the scalar model in the 2D and 3D cases, is here seen in the bunching profile of the vectorial model. For both pump polarization cases, the bunching snapshots in Figure 7c,d show a large lobe emerging backward along the z axis (pump field axis), together with the forward lobe already present at the time $t=0$ in Figure 6b and due to the diffraction of the pump field [13]. The backward lobe is the well-known signature of the collective recoil lasing effect.

There are some interesting features that can be highlighted from the snapshots (a) and (b) of Figure 7. First, both panels show the cloud atomic density redistribution into a longitudinal 1D grating. Secondly, the polarization effects are visible in both the cases: the left panel shows horizontal grating due to the suppressed scattering in the vertical direction, whereas the right panel shows a vertical grating, consequence of the forbidden scattering along the horizontal direction. Lastly, both clouds undergo a squeezing effect on the radial traverse direction with respect to the z axis, stronger at the right-end of both ensembles. This electrostriction effect is a phenomenon typically observed in dielectrics, which can change their shape whenever they are subject to a strong electric field.

In the collective recoil lasing this electrostrictive effect was already observed by the scalar model in [13], with a 2D and a 3D atomic distribution. The current analysis confirms the presence of this effect with a more realistic model taking into account the light polarization. This force is not yet well understood and further investigation is needed for a satisfactory explanation. For instance, it appears that the effect rides along the cloud like a wave: it starts at the right tip of the cloud, as displayed in Figure 7a,c, and it travels backwards along the z axis. After this sort of backwards-propagating wave has passed a section, this section is not susceptible to the polarization vector of the field anymore, but it is still reactive to the gradient of the radiation force. There are some studies which have experimentally observed similar results, as the one presented in [21].

3.2.4. Scattering from a Spherical Atomic Distribution

As with previous cases, a density grating also emerged when a spherical cluster of atoms is probed (see Figure 8), although no tilted planes are defined for such case. This comes because of the symmetry of the spherical initial cloud, so that no initial preferred scattering directions are present, and the atoms start to scatter randomly, except that in the direction of the polarization vector. We observed that two dynamics occurring at the same time in this case: the first one moves the inner atoms towards the shell and the outer atoms inwards; the second one distributes the particles into bunches, avoiding the polarization vector axis. On the one hand, the first action seems to be due to the electrostrictive force (outer atoms moving inwards) and the scattering effect (inner atoms moving outwards, thanks to a higher optical thickness along the z axis). The combination of these two processes deforms the edges of the cloud into an ’egg-like’ shape. In addition, this effect forces the atoms to be forward-displaced faster the closer they are to the outer shell, due to the definition of a preferred scattering path.

On the other hand, the second dynamics arises because of the photon-scattering, which starts to occur across this newly defined longest path along the external shell of the cloud. The scattering helps sketching arched rods of atoms on the outer surface vertically—Figure 9a—or horizontally—Figure 9b—oriented, depending on the polarization vector of the pump field. If the pump were to be set as a scalar field as in ref. [13], a collection of concentric circles could be spotted. However, in reality, such circles would represent the edges of shallow circular cones, which have their center located on the z axis, with a value smaller than the values of their edges.

Observing the bunching formations, in Figure 9c,d, it is possible to see that the projections onto the ($x,y$) plane present an elongated shape along y and x, respectively, as direct consequence of having a linearly polarized pump field. When the scalar model is used, these projections appear as a clover with four leaves, and they represent a segmentation of the spherical cloud into smaller spherical bunches of atoms following a 3D pattern.

4. Conclusions

The main scope of this work is to show the density grating formation caused by collective atomic recoil lasing, triggered in a cold atomic cloud when a polarized external optical field is applied. The results attained are similar to the ones observed with a previous scalar model [13], with the exception that when the scattering direction coincides with the polarization vector of the pump field, the scattered field is suppressed. Minor differences are seen between the patterns achieved for the scalar and the vectorial models, due to the presence of different short-range terms in either model. If the cold atomic vapor is shaped into a spherical distribution, the simulations succeed in clearly showing the electrostrictive nature of the force generated by the collective scattering process: the cloud suffers an elongation along the pump propagation axis, hence establishing the longest scattering path along the edge or the outer shell of a 2D or 3D systems, respectively. However, the experimental observation of CARL when the cold atomic vapor is shaped into a spherical distribution remains a difficult task, compared to the more favorable ellipsoidal case. Nevertheless, the role of the electrostrictive force creating some preferred direction where collective scattering set up, it is intriguing and should deserve a deeper investigation.

The present model of equations has also some numerical advantages. In fact, since the model only depends on the particles’ positions, it allows its implementation into molecular dynamics (MD) algorithm, typically employed in other scientific fields. The new open-access programming language Julia has revealed to be of great instrument to implement such MD codes, making it swift when using the built-in Julia algorithms, while allowing the accurate following of the trajectory of each particle and respecting the energy conservation of the system.