Optomechanical transport of cold atoms induced by structured light

Optomechanical pattern forming instabilities in a cloud of cold atoms lead to self-organized spatial structures of light and atoms. Here, we consider the optomechanical self-structuring of a cold atomic cloud in the presence of a phase structured input field, carrying orbital angular momentum. For a planar ring cavity setup, a model of coupled cavity field and atomic density equations describes a wide range of drifting modulation instabilities in the transverse plane. This leads to the formation of rotating self-organized rings of light-atom lattices. Using linear stability analysis and numerical simulations of the coupled atomic and optical dynamics, we demonstrate the presence of macroscopic atomic transport corresponding to the pattern rotation, induced by the structured pump phase profile


I. INTRODUCTION
The spontaneous emergence of spatiotemporal order is a prominent feature of physical systems driven far from equilibrium [1]. Optical systems provide valuable platforms for studying spatial self-organization, where the self-sustained patterns are typically encoded in the internal excitations of a nonlinear medium [2,3]. Besides their fundamental interest, optical spatial structures offer potential applications in information processing such as optical memories and registers [4].
Recent experiments have provided paradigmatic examples of spatial self-organization in transverse nonlinear optics with cold atoms such as density, electronic and magnetic spatial ordering [5][6][7]. In the first case, the bunching of atoms due to optomechanical forces can provide positive feedback and lead to spatial instabilities [8]. Other mechanisms for pattern formation in cold atoms rely instead on optical pumping and internal state nonlinearities. In the regime of quantum degeneracy, the non-equilibrium phase transition corresponding to the onset of spatial self-organization has been interpreted as a quantum phase transition [9,10]. Multimodal configurations have been shown to enrich significantly the physical scenario, allowing the realization of frustrated interactions and supersolid states with atomic gases [11,12].
The introduction of orbital angular momentum of light (OAM) and its relative ease of generation and control paved the way to a plethora of applications in optical manipulation and information technologies [13,14]. Focusing on cold atoms, the use of OAM modes enables transfer of angular momentum to both motional and internal atomic degrees of freedom [15]. Macroscopic rotation of cold atomic gases is also achieved in the dispersive regime by means of rotating ring lattice trapping potentials [16,17]. Such geometries and related configurations, when applied to cold atoms, offer the possibility to engineer classical or quantum transport [18] and to simulate topological condensed matter effects such as fractional quantum Hall physics [19].
In this work we extend previous studies of optomechanical instabilities in cold atoms, allowing the input field to have * giuseppe.baio@strath.ac.uk a spatially dependent phase structure and carrying non-zero OAM. Externally controlled phase gradients are known to induce drifting pattern dynamics, as demonstrated with hot atomic vapors [20,21], and can be applied to control the motion of dissipative solitons [22][23][24]. Moreover, drifting dynamics induced by OAM was observed for patterned states in a single mirror system with photorefractive nonlinearity [25]. For a cold atomic gas with optomechanical nonlinearity, we show that such dynamics, corresponding to the motion of coupled structures of light and atoms, is capable of inducing controllable atomic mass transport. Our atomic ensemble is assumed to be kept at constant temperature with optical molasses, providing strong damping to the center of mass dynamics. Moreover, dealing with density redistribution-induced interactions, atoms behave as linear scatterers. This opens the possibility of extending the present results to soft matter systems such as suspensions of colloidal particles [26][27][28].
The paper is organised as follows. In Sec. II we review the known features of optomechanical self-structuring in a ring cavity for a spatially homogeneous input pump. In Sec. III we discuss drifting pattern dynamics induced by OAM in the pump profile, which leads to the formation of self-organized rotating ring lattices. Finally, in Sec. IV, we compare our predictions with numerical results from particle dynamics simulations, showing the presence of macroscopic optomechanical transport corresponding to the pattern rotation.

II. THE MODEL
In this section we discuss a theoretical model, adapted from Refs. [29,30], describing the transverse dynamics of the cavity field and the density of the atomic ensemble. A linear stability analysis provides the threshold condition for transverse optomechanical modulation instabilities (MI) of the homogeneous (flat) stationary states.

A. Model Equations
We consider a cold thermal cloud of two-level atoms within a planar ring cavity geometry of effective length L, sketched in Fig. (1). Theoretical models describing the coupled dynam- plitude A I (r) and phase exp (ilϕ) drives a single longitudinal cavity mode. The ring cavity has effective length L and mirror transmittivity τ. The intra-cavity field E(r,t) interacts with a cloud of two-level atoms of thickness l, average atomic density N 0 and temperature T . ics of cavity field and density modulations of a cold atomic gas involve the paraxial wave equation for a single longitudinal mode of the cavity field E(r,t) coupled to the atomic density distribution N(r,t) = N 0 n(r,t), where N 0 is the average density of the sample and r denotes the transverse spatial coordinate [30]. The field equation is derived in the slowly varying envelope, rotating wave and mean field approximations as follows [31]: where θ is the cavity-pump detuning, describing the linear shift on the cavity field, A I (r) represents a spatially dependent pump and diffraction is described by the transverse Laplacian ∇ 2 ⊥ . Here we focus on a pump profile of the form A I (r) = A I (r) exp (ilϕ), where A I (r) is a generic radial profile and r, ϕ are the radial and azimuthal coordinates respectively. In addition, κ = cτ/L represents the cavity linewidth and √ a = λ L/4πτ the diffraction length, where λ is the light wavelength, τ the mirror trasmittivity and c the speed of light in vacuum. The response of the two level saturable medium 1 is parametrized in terms of the complex susceptibilityγ = 2C(1 + i∆), with ∆ = 2δ /Γ, where δ is the light detuning with respect to the atomic resonance. Finally, the parameter C = b 0 /2τ(1 + ∆ 2 ) is known as cooperativity parameter and it depends on the optical thickness of the sample at resonance b 0 , τ and ∆. For our description of atomic motion, we focus on the limit when the population of the excited atomic state is negligible. This is captured by the saturation parameter 2 s = |E(r,t)| 2 . In general, the resulting dipole force 1 The characteristic dependence is obtained by adiabatic elimination of the atomic variables. 2 The cavity field E is rescaled by the saturation intensity at detuning ∆, i.e.
[I sat (1 + ∆ 2 )] 1 2 , where I sat represents the saturation intensity at resonance. f dip acting on the atom center of mass can be derived from the AC Stark potential: . (2) Note the dependence of Eq. (2) on the detuning ∆, such that the atom is trapped in high (low) intensity regions when ∆ < 0 (∆ > 0). Assuming the presence of optical molasses, we adopt a classical treatment of positions r j and momenta p j (i, j = 1, . . . , N), leading to the canonical set of stochastic differential equations: where m is the atomic mass, γ describes friction (momentum damping) and ξ j (t) is a stochastic variable such that ξ j (t) = 0 and ξ j (t)ξ j (t ) = 2D p δ j j δ (t − t ), with D p being the diffusion constant in the space of momenta. In the thermodynamic limit (N → ∞), the above equations map to the Chandrasekhar equation for the phase space distribution function f (r, p,t) in the 1-particle picture [27,32]: where β = 1/k B T . Note that, in the presence of damping, the collisional contribution in Eq. (4) is replaced by the term ∇ p ( f p) + m β ∇ 2 p f , describing relaxation in momentum space. In the limit of strong viscous damping (more generally, when t 1/γ), the momentum p is eliminated from the dynamics and one derives the Smoluchowski equation for the density distribution n(r,t), by direct integration of Eq. (4) in momentum space, namely: where diffusive dynamics, as an effect of the optical molasses, is characterized by the constant D = 1/γmβ [33]. Clearly, in the case of hot atoms (β → 0), the diffusive term dominates and one has n(r,t) ≈ 1, i.e. the homogeneous state. Generally speaking, MI in our system is due to a competition between two-level and density-driven nonlinearities. In what follows, we perform a linear stability analysis, showing that density inhomogeneities alone are sufficient to achieve (optomechanical) MI. The features of the resulting self-organized states in the presence of OAM are discussed in Sec. (III).

B. Linear Stability Analysis
The family of homogeneous stationary states of our system is simply found imposing that all derivatives in both Eqs. (1,5) vanish. For the cavity field, i. e. ∂ E 0 /∂t = 0 and ∇ 2 ⊥ E 0 = 0 in Eq. (1), this yields: with n 0 = 1. A consequence of the two-level nonlinear term in Eq. (6) is that the homogeneous field intensity |E 0 | 2 , as a function of the pump intensity A I , is not necessarily single valued [31]. It is easy to see that, in the low saturation limit, our model Eqs. (1) and (5) read as follows: where time and space coordinates are rescaled as t = κt, r = r/ √ a, D = D/κa, while the constant σ =hΓ∆/4k B T denotes an optomechanical coupling strength. Moreover, Eq. (6) now simply reads A I = [1 + i(θ + 2C∆)] E 0 , so that there is always a one-to-one correspondence between the values of |A I | 2 and |E 0 | 2 in the low saturation limit. Perturbations around steady states values E 0 and n 0 are introduced in Fourier space as δ E(q,t ) = a e i(q·r +νt ) + b * e −i(q·r +ν * t ) , δ n(q,t ) = c e i(q·r +νt ) + c * e −i(q·r +ν * t ) . Thus, defining the variable x = (a, b, c) T , the linearized system obtained from Eqs. (7,8) can be cast in the usual matrix form (M−νI)x = 0, where M is the following matrix: By imposing the marginal stability condition, i.e. ν = 0 or det(M) = 0, one finds the following analytical threshold: leading to the critical wavevector q 2 c = 1 − (θ + 2C∆). Here we study the linear stability properties by spanning two different parameters regions. The two instability diagrams in Fig.  (2) show the dependence of the linear growth rate (GR), given by -Im[ν], as a function of |E 0 | 2 vs. detuning ∆ (Fig.2 -Top) 3 and |E 0 | 2 vs. temperature T (Fig.2 -Bottom). The value of Im[ν] is computed by estimating the minimum imaginary part among the eigenvalues of M in Eq. (9). Finally, the homogeneously pumped system spontaneously selects an hexagonal patterned state for values of |E 0 | close to threshold [34]. We remark here that, assuming low saturation, purely optomechanical structures occurs in the proximity of the threshold only since, for stronger pumping, electronic nonlinearity becomes relevant. 3 Only the case ∆ > 0 is shown as we have symmetric behaviour for ∆ < 0.

III. PATTERN DYNAMICS WITH STRUCTURED PHASE
In this section we analyze the effect of the spatially dependent pump A I (r ) on the transverse geometry and dynamics of the patterned states of our system. Before presenting numerical results, we show how drifting pattern dynamics arises from Eqs. (7,8), based on the argument from Refs. [22,35].

A. Drifting pattern dynamics
In order to include OAM in our system, the pump rate A I (r ) must have an azimuthally dependent phase exp(ilϕ), where l ∈ Z. Indeed, the factor exp (ilϕ) generates l inter-twined phase fronts and a non-trivial topological structure due to the phase singularity at r = 0 [14]. In general, any phase profile including the factor exp (ilϕ) is responsible for a nonvanishing contribution of OAM through the transverse plane. Let us further simplify the analysis recalling that the stationary solution of Eq. (8) is given by the Gibbs distribution [29]: The steady states of the coupled system above MI can in principle be obtained by integrating numerically Eq. (7), simply eliminating the density distribution. Thus, Eq. (11) acts back as the nonlinear term in the field equation, namely [30]: where we emphasize here the dependence of n eq (r ) on the field intensity. Therefore, the coupled system of Eqs. (7,8) reduces now to Eq. (12) only. Let us consider the local transformation E(r ,t ) =Ẽ(r ,t ) exp(ilϕ). The corresponding equation forẼ(r ,t ) reads now as follows: One can recognize that the LHS of Eq. (13) has the form of a covariant derivative D t = ∂ ∂t + v dr · ∇ ⊥ , i.e. the expression of a time derivative in a locally co-moving reference frame, according to the velocity v dr (r ) = 2l∇ ⊥ ϕ = 2l r φ. Thus, in the presence of OAM, the family of patterned solutions are seen in the laboratory frame drifting with a velocity v dr (r ) [22]. Furthermore, such an argument carries implications also for the geometry of the patterned state in our case. Indeed, the RHS of Eq. (13) shows a spatially dependent contribution to the decay rate and cavity detuning. Therefore, the spatial rigidity of the patterned phases is broken according to the spatial symmetry of the pump A I (r ) and the stationary solution is retrieved in a rotating reference frame defined (in polar coordinates) byφ(r ,t ) = ϕ − ω(r )t , with the angular frequency ω(r ) = 2l/r 2 . This corresponds to a differential rotation of the concentric patterned solutions at a fixed radius, as found in Ref. [35]. In what follows, we support the above observation by integrating Eq. (12) numerically.

B. Numerical Results
We shall restrict here our analysis to a 2D radial "top hat" profile with rapidly vanishing tails. A commonly used example is given by the following: FIG. 3: Formation of self-organized ring lattices with a red detuned pump intensity 5% above threshold and l = 1. The 2D intensity I(r) (left) and density n eq (r) are shown at different integration times. Parameters are chosen as follows: θ = 4, C∆ = −3.5 (∆ < 0), σ = −25. Note that atoms bunch towards the maxima of intensity.
where the constant η > 0 determines the side steepness of the flat part with radius r 0 . The advantage of using Eq. (14) is twofold: the used boundary conditions do not affect numerical results and allow us to investigate the transverse 2D dynamics induced by the structured phase factor exp (ilϕ) on a wide area inside the integration domain.
We integrate Eq. (12) using a Fourier split-step method. A domain of 10 critical wavelengths is discretized in 256×256 points and time step dt = 5 × 10 −3 . Fig. (3) shows an example of optomechanical self-structuring with OAM index l = 1. The presence of a helical phase exp (ilϕ) in the pump profile (14) causes the field to vanish at r = 0 and induces the formation of diffractive rings whose intensity exceeds the MI threshold. Such rings are spaced by a critical wavelength, visible from Fig. (3) when t = κt = 20. Since the intensity between contiguous rings is below the MI threshold, each ring undergoes a 1D angular instability independently, leading to a set of concentric ring lattices structures [35]. Moreover, by varying η and r 0 , one can also influence the features of the patterned states. For example, increasing the steepness with the parameter η may lead to outer diffraction rings faster achieving azimuthal MI. For any balanced choice of such parameters, selfstructuring dynamics always leads to a set of ring light-density lattices, each one rotating independently from the others. Similar behaviour is found with higher order OAM modes (l > 1).
The self-organized light-atom ring lattices display differential rotation, in accordance to the angular frequency ω(r ) = 2l/r 2 . To show this, we consider the field in Eq. (12) at a fixed radius r = r 0 , mapping the variation of ϕ ∈ [0, 2π] to an effective 1D model with periodic boundary conditions and an input pump of the form A I (ϕ) = A I e ilϕ . The same argument is repeated in Sec. IV for studying atomic transport. Fig. (4) shows a comparison, for l = 1, 2, 3, of the rotation frequency estimated at several radial distances and showing good agreement with the predicted ω(r ). Analogous rotational behaviour can be expected for the case of bistable optomechanical dissipative solitons, as those found in [36]. OAM or other structured phase pump profiles can thus be used in this case to control the motion and the effective interactions of such localized dissipative structures [23,24].
In the rest of this section, we discuss the effect of varying the diffusion coefficient D on the ring lattice rotational dynamics. In order to take into account this dependence and obtain realistic values of the rotation frequency comparable with experiments, we numerically integrate the coupled cavity Eq. (7)  case at fixed radius 4 . Note that, since we have D = 1/γmβ , any variation of D at constant temperature is inversely proportional to a variation in the momentum damping coefficient γ, in the particle Eqs. (3). Results are shown in Fig. (5), where values of the estimated angular frequency are plotted against the rescaled diffusion coefficient D = D/κa. For higher values of D (decreasing the friction γ), one observes a faster ring lattice rotation speed which saturates to a value roughly proportional to the index l. On the other hand, by increasing the friction γ, the rotation of the ring lattices created by the dipole interaction between atoms and light is affected by a drag, due to slower atomic diffusion. In the limit of very small diffusion (D 1), this drag can compensate the effect of the phase gradient introduced by the pump OAM and the concentric ring lattices would be expected to remain stationary. The above observations have profound consequences from the point of view of atomic center-of-mass dynamics since drifting self-organized light-density modulations induce transport of the self-trapped atoms, as shown in the next section.

IV. ATOMIC TRANSPORT
In this section we demonstrate OAM induced atomic transport corresponding to the rotation of the self-organized ring lattices. This is done by simply addressing the average momentum of the atomic distribution which quantifies the macroscopic rotational motion of the atoms, occurring at the onset of self-organization. We restrict ourselves to the 1D case i.e. atoms confined in a ring geometry at a fixed radius within the cavity 2D transverse domain. Such a case can be realized experimentally by appropriately using Laguerre-Gauss beams of a given OAM in order to excite a single patterned ring [15].

A. Average momentum
We start our observations by recalling that, at equilibrium, the Chandrasekhar Eq. (4) for a dilute 1D atomic gas at constant temperature T is solved by the Maxwell-Boltzmann distributionf (p), namely [37]: where p is the average (bulk) momentum around which the distributionf (p) is centered. Since we are interested in the strong friction limit, we assume that the total phase space distribution has the following factorized form [29,37]: where r varies over a finite interval, depending on the chosen radius. Macroscopic transport in this regime can be addressed by looking at the average momentum p , whose time evolution is obtained multiplying Eq. (4) by the momentum p and integrating over momentum space [38]: where f dip is the dipole force in Eq. (2), expressed in the low saturation limit. We now rearrange Eq. (17), according to the following simple manipulations. Firstly, we have that p 2 ∂ f ∂ r d p = ∂ ∂ r n p 2 , where p 2 = (mβ ) −1 is the second moment off (p) in Eq. (15).
Atomic currents can thus be obtained integrating Eq. (18) along the 1D azimuthal domain. However, the first term in its RHS vanishes in our case 5 so that, after integration, the average momentum finally reads: This quantifies the steady state mass current along a rotating ring in the overdamped limit, from the knowledge of n(r,t) only. In what follows, we compare p(t) density with the ensemble averages obtained from simulating particle dynamics.

B. Particle model
Since a critical wavenumber q c is selected at the MI threshold of Eq. (7), it is convenient to re-define here atomic positionr j = q c r j and momentap = p j /hq c , so that the set of particle Eqs. (3) in 1D now reads [39]: where ω p =hq c /2mκ defines characteristic energy scales of an atom oscillating in a standing wave potential of modulation q c . In the overdamped case, one imposes momentap j at their steady state values, leading to the following epression of the average momentum or atomic current: i.e. simply obtained as an ensemble average. Most importantly, this provides an expression of the relative constants in Eq. (5) in terms of microscopic parameters: where the diffusion D and σ have been expressed in terms of the (scaled) momentum spread ζ 2 p . Note that D p =γζ 2 p . We now integrate numerically the particle Eqs. (20) in the overdamped limit. Those are coupled to the field equation in 1D through a density profile n(r,t ), numerically reconstructed at all times from the particle positionsr j (t ).
A typical outcome of our simulations is presented in Fig.  (6), for the case of 5 × 10 4 atoms. Consistently with the 2D numerical integration in Sec. III, we observe drifting pattern dynamics in the 1D azimuthal geometry, i.e. ring lattice rotation at a fixed radius, induced by the OAM carried by the input pump. The results displayed in Fig. (7) show that, at the onset of self-organization 6 , the atoms spontaneously develop net mass current along 1D rings, in the presence of OAM. This FIG. 6: Numerical results from the 1D particle simulation for a set of 5 × 10 4 atoms. (a) κt vs.r plot showing the time evolution (κt max = 60) of the 1D cavity field intensity I(r). (b) Selfordering of the density n(r,t), numerically reconstructed from the particle trajectories. For blue detuned light (∆ > 0), atoms bunch in the minima of I(r). The OAM dependent slope of the patterned phase in this graph corresponds to the rotation of the pattern in a ring geometry. Parameters are chosen as follows: OAM index l = 3, Γ κ =ω p =γ = 1, ζ p = 0.707, ∆ = 100, such that D ≈ 2, σ ≈ 25.
corresponds to the slope (drift) of the patterned field intensity I(r,t ) and atomic density n(r,t ), displayed in the t vs.r plots in Fig. (6). In particular, in Fig. (7 -Top), we compare the values of the current p j (t ) (in units ofhq c ) obtained by means of both the above approaches, i.e. from Eq. (19) (density) and Eq. (21) (particles). In both cases (l = 3, 4), one observes a exponential-like increase of p j (t ) , which eventually evolves to fluctuations around a steady state value. Moreover, the two series show a highly correlated behaviour. Finally, as shown in Fig. (7 -Bottom), in order to obtain more precise values of the steady state currents p j (t ) for different OAM indices l = 1, . . . , 4 , we averaged over a set of 10 simulations with the same set of parameters for each value of l. As a final remark, we observe from Fig. (7 -Bottom), that the steady state current values increases non linearly with the OAM index l. This might suggest that the optomechanical transfer of OAM considered here, i.e. from the cavity field to the atoms, occurs more efficiently with higher order OAM modes input. However, a detailed analysis of this mechanism in such a configuration is out of the aim of this work and will be investigated in future studies.

V. CONCLUDING REMARKS
Clouds of cold atoms in optical cavities, under the action of a coherent beam of light carrying OAM, can spontaneously form concentric rings of rotating light-atom lattices in the transverse plane, by means of optomechanical self-structuring (See Fig. (3)). The rotation of these spatio-temporal structures is due to the gradient of the transverse phase distribution externally imposed on the pump beam. The observed angular velocity is directly proportional to the OAM and inversely proportional to the square of the ring radius (See Fig. (4)). Moreover, for increasing friction and decreasing diffusion, the atoms introduce a drag that slows down the rotation of the patterned rings (See Fig. (5)). By means of numerical evidences, we showed that such ring lattice rotation can sustain macroscopic atomic motion along azimuthal 1D rings in the 2D transverse domain (See Figs. (6, 7)). Therefore, optomechanical transport and azimuthal mass currents induced by the OAM in the structured pump have been predicted for a cold thermal gas in the assumption of overdamped motion.
A possible generalization of the case considered here is that of rotating optomechanical instabilities without momentum damping [40]. Similarly, extensions to the dynamics of optomechanical dissipative solitons and their interaction in the presence of an OAM carrying pump are also of interest for future theoretical studies and experimental realizations. Finally, our results can be applied to the case of a quantum degenerate gas, in order to study vortex formation and turbulent regimes [41]. The OAM transfer mechanism discussed here represents a potential platform to induce quantized persistent currents of ultra-cold atoms, where such states are indeed an ideal candidate for the realization of atomtronic devices in circular 1D atomic traps [14,15].

ACKNOWLEDGMENTS
We are grateful to T. Ackemann and S. Franke-Arnold for helpful discussions. All authors acknowledge financial support from the European Training Network ColOpt, which is funded by the European Union (EU) Horizon 2020 programme under the Marie Skłodowska-Curie Action, grant agreement 721465.