1 Introduction↩︎

In the last 40 years, laser cooling of atoms to very low temperatures [1] has been a crucial experimental achievement in quantum and solid state physics, with several Nobel prices awarding the development of techniques in this subject. Among them, cold atoms in dissipative optical lattices [2] still play an important role, and Sisyphus cooling in particular is still being widely used to cool atoms below the Doppler temperature limit.

Cold atoms offer an invaluable setup to study directed transport, and many experimental realizations of dissipative [3]–[8] and nondissipative [9], [10] ratchets [11], [12] have been demonstrated on them. Additionally, they are also known to exhibit unsual transport behaviour beyond Boltzmann–Gibbs statistical mechanics [13], usually referred as a regime of infinite density [14]–[17] in which the probability distributions are non-normalizable due to the broken ergodicity.

This paper focused in cold atoms in standard setups of optical lattices associated with Sisyphus cooling. They are created by counterpropagating laser beams, and can be tuned experimentally with a high level of control. Here they are studied theoretically at the level of the semiclassical equations for the atomic phase space densities, which allows a detailed description of the Sisyphus mechanism.

Special attention is paid to atomic transport. A perturbing probe beam is introduced [2], [18], [19] to break the system symmetries [12], [20] and put the atoms into directed motion. The probe perturbation excites atomic density waves, which are referred to as Brillouin propagating modes [19] due to the resemblance to acoustic waves rippling through a dense fluid, even though in the present case the optical lattice is dilutely occupied and the mode is sustained in a medium of non-interacting particles due to the interaction with the laser beams.

Previous research has focused mostly on one propagating mode, the one with same frequency and wave number as the propagating perturbation, ignoring the effect of other excited modes, propagating or not. Here we address this issue by deriving a series expansion of the current on the amplitudes of the excited atomic density waves. Additionally, the expansion directly reveals a singularity that is identified with the threshold to the regime of infinity density. The threshold values are in agreement with previous analytical results [14], [15] obtained from a simplified, approximate Fokker-Planck equation [21] based on space and Zeeman sublevel averaging of the semiclassical equations—thus with uncontrolled approximations. The method presented here does not make such approximations, being based on the semiclassical equations it takes explicitly into account the microscopic origin of Sisyphus cooling, thus providing desirable support for previous results on infinity densities.

The paper is organized as follows. In Sec. 2 the system models studied are described in detail. The new method, relating the current and higher moments to the atomic density modes, is presented in Sec. 3. The application of the method to the system models, its analytical results and numerical validation, is in Sec. 4. The transition to the regime of infinity density, coming out naturally from a singularity in the analytical results, is discussed in in Sec. 5. Finally, Sec. 6 ends with the main conclusions.

2 System models↩︎

First, we consider the simplest model of a dissipative optical lattice with Sisyphus cooling, a one-dimensional (1D) system generated by atoms with a transition \(J_g= 1/2 \rightarrow J_e=3/2\), mass \(m_a\), illuminated by two counter-propagating laser fields with orthogonal linear polarizations. This setup generates a 1D optical lattice, so-called lin\(\perp\)lin configuration [22]. The atom, in the Zeeman sub-level of the atomic ground-state \(+\) or \(-\), experiences the optical potential \[U_{\pm}(x)= \frac{U_0}{2}[-2\pm \cos(2k_l x)], \label{eq:upm1D}\tag{1}\] where \(x\) is the laser beam propagation axis, \(k_l\) the laser field wave vector, and \(U_0\) the optical lattice depth.

Further details of the dissipative optical lattice are given, in the semi-classical approximation for weak laser intentities, by the following [23] coupled Fokker-Planck equations for the phase space density \(P_{\pm}(x,p,t)\) of atoms in the ground state Zeeman sublevel \(|\pm\rangle=|J_g=1/2,M_g=\pm1/2\rangle\) at the position \(x\) with momentum \(p\), \[\begin{align} \left[ \frac{\partial}{\partial t} + \frac{p}{m_a}\frac{\partial}{\partial x} -U_{\pm}^\prime(x)\frac{\partial}{\partial p}+F_{\pm}(x,t)\frac{\partial}{\partial p}\right] P_{\pm}= \nonumber\\ -\gamma_{\pm}(x) P_{\pm}+\gamma_{\mp}(x) P_{\mp} +\frac{\partial^2}{\partial p^2}\left[D_0 P_{\pm}\right], \label{eq:fp} \end{align}\tag{2}\] where \(\gamma_{\pm}\) is the transition rate between the ground state sublevels, \(D_0\) is a noise strength describing the random momentum jumps that result from the interaction with the photons, and \(F_\pm(x,t)\) is generally a non-conservative force—coming, for example, from radiation pressure, or, as a second example, being an arbitrary time-dependent driving force \(F(t)\) that can be generated by phase modulating one the lattice beams [24].

Equation (2 ) is complemented by the normalization condition, \[\int dx \int \!dp\left[ P_{-}(x,p,t) + P_{+}(x,p,t) \right] =1. \label{eq:norm:fp}\tag{3}\]

A quantity of special interest is the average atomic current, which measures the directed motion, being defined as \[\begin{align} &\langle v \rangle = \lim_{t\rightarrow\infty} \frac{\langle [x(t)-x(0)] \rangle }{t} = \nonumber\\ &\lim_{t\rightarrow\infty} \frac{ \int_0^t \!\! dt' \int \!\! dx \int \!\! dp \, \frac{p}{m_a}[P_+(x,p,t')+P_-(x,p,t')]}{t}. \label{eq:currdef} \end{align}\tag{4}\]

2.1 1D lin\(\perp\)lin setup↩︎

In the 1D lin\(\perp\)lin configuration without driving there is no radiation pressure, i.e. \(F_\pm(x,t)=0\), the optical potential is given by (1 ), with \[U_0=-2\hbar \Delta'/3>0, \label{eq:u0:def}\tag{5}\] where \(\Delta'\) is the light-shift per beam for a closed transition having a Clebsch-Gordan coefficient equal to unity. Furthermore, the transition rates between the internal states are given by \[\gamma_{\pm}(x)=g_0\pm g_1\cos(k_0 x), \label{eq:rates}\tag{6}\] where \(k_0=2k_l\), \(g_0=\Gamma^\prime/9\) and \(g_1=g_0\), with \(\Gamma^\prime\) being the photon scattering rate per lattice beam.

Note that in writing (2 ) we are not explicitly considering an extra diffusive term that comes up in the semiclassical approximation [23], \(\partial^2 D_{\mp\pm} P_{\mp}/\partial p^2\), which describes a further momentum kick when the transitions take place, and has to be corrected ad hoc to avoid artificial singuralities produced by the semi-classical approximation [15]. When corrected, the effect of that term is nevertheless small, and commonlly neglected, unless shallow optical potentials are considered [15]—more about this issue in Sec. 5.

In addition, we are not explicitly considering any state or spatial dependence of the noise coefficient \(D_0\), because their effect is observed to be small in the simulation results reported in this paper. Neglecting \(D_{\mp\pm}\) and taking the space average of the remaining diffusion coefficient in the original semi-classical equations [23] yields \(D_0=35\hbar^2k_l^2\Gamma^\prime/90\).

The application of an additional weak probe beam generically produces extra small contributions to all the functions above: the optical potential, radiation pressure forces, transition rates and noise terms, though the most relevant contribution usually is the one in the optical potential, being the others negligible. We consider here a probe beam that is polarized parallel to the 1D counter-propagating lattice beam, which yields the following optical potential \[\begin{align} U_{\pm}(x)= \frac{U_0}{2}[-2\pm \cos(k_0 x) + \varepsilon_p \cos(k_0 x-\delta_p t+\phi_p) ], \nonumber\\ \label{eq:upm1D:v1b} \end{align}\tag{7}\] where \(\varepsilon_p=2E_p/E_0\) is twice the ratio between the electric field of the probe \(E_p\) to that of the underlying optical lattice \(E_0\), \(\delta_p\) is the probe frequency detuning and \(\phi_p\) the probe phase.

2.2 3D lin\(\perp\)lin setup↩︎

In addition to the above 1D lin\(\perp\)lin configuration, we will study also the one-dimensional system model that arises in the standard 3D-lin\(\perp\)lin configuration [2], after neglecting movement in the two directions perpendicular to the direction of interest, usually taken as the \(x\)-axis [18], [19], [25]–[27]. A weak \(y\)-polarized probe is added in the \(z-\)direction [19], [27], yielding an extra term in the optical potential which is the superposition of two sinusoidal potentials travelling in the \(x\)-direction, each with opposite velocity and similar shape as in the previous 1D system model. More specifically, by formally considering \(y=z=0\) and one of the travelling probe drives, the optical potential in each sublevel of the ground state is given by [2], [27] \[\begin{align} U_\pm(x)=\frac{U_0}{2}\Big[-\frac{3}{2}-\frac{1}{2}\cos(2 k_0 x) \pm \cos(k_0 x) \nonumber \\ + \varepsilon_p\cos(k_0 x -\delta_p t+\phi_p) \Big], \label{eq:Upm:3D} \end{align}\tag{8}\] with transition probability rates between them given by \[\gamma_{\pm}(x)=g_0\pm g_1\cos(k_0 x) + g_2\cos(2 k_0 x) , \label{eq:gammapm}\tag{9}\] where \(k_0=k_l\sin\theta_x\), \(U_0=-16\hbar \Delta'_0/3\), \(\Delta'_0\)(\(<0\)) is the light-shift per lattice field, \(g_0=2\Gamma_s/3\), \(g_1=8\Gamma_s/9\), \(g_2=2\Gamma_s/9\), and \(\Gamma_s\) is the photon scattering rate per lattice beam, and now \(\varepsilon_p=E_p/(2E_0)\). Like before, for the sake of simplicity, we are neglecting the probe contribution to the transition rates, the radiation forces, and noise terms, being observed their effect to be small in the simulations. The equation to solve is also (2 ), with the noise strength \(D_0=5\hbar^2k_0^2\Gamma_s/18\).

3 Fourier mode theory↩︎

We develop in this section a theory that will allow us to visualize the atomic density modes excited by the probe and their precise contribution to the current.

First, the atomic Fourier modes are defined from the phase space density \(P_{\pm}(z,p,t)\) by means of the following Fourier transform \[\mathcal{P}^\pm_{\omega,k,q}=\frac{1}{T_p}\int_0^{T_p}\!\!\! dt \, e^{-i\omega t}\int \!dx \, e^{i k x}\int \!dp \,e^{i pq/\hbar} P_{\pm}(x,p,t), \label{eq:pfourier:def}\tag{10}\] where \(T_p=2\pi/\delta_p\) is the time period introduced by the probe. The fact that the driving probe is time periodic allows us to focus only on solutions which has the same periodicity, i.e. to \[\omega=l \delta_p, \label{eq:om:discrete}\tag{11}\] with \(l\) integer—thus neglecting transitory dependencies on a specific initial condition, which are expected to die out after a transient time interval, leaving the periodic solution. In addition, real phase space densities implies \[\mathcal{P}^\pm_{-\omega,-k,-q}=\left( \mathcal{P}^{\pm}_{\omega,k,q}\right)^*, \label{eq:pfourier:def:real}\tag{12}\] where \(*\) denotes complex conjugate.

Now, let us consider a generic 1D setup determined by (2 ) and (3 ), with a periodic optical potential \[U_\pm(x+2\pi/k_0)=U_\pm(x),\] and space periodic transition rates, \[\gamma_\pm(x+2\pi/k_0)=\gamma_\pm(x),\] and time and space periodic driving, \[F_\pm(x+2\pi/k_0,t)=F_\pm(x,t+2\pi/\delta_p)=F_\pm(x,t), \label{eq:force:period}\tag{13}\] for all \(x, t\).

Using (10 ), the Fokker-Planck Eq. (2 ) is transformed into \[\begin{align} & i\omega \mathcal{P}^\pm_{\omega,k,q} -\frac{\hbar k}{m_a}\frac{\partial}{\partial q}\mathcal{P}^\pm_{\omega,k,q} \nonumber\\ & -\frac{iq}{\hbar}\Big( \sum_n\mathcal{F}^{(0)\pm}_n \mathcal{P}^\pm_{\omega,k-nk_0,q} + \sum_{l,m} \mathcal{F}^{p\pm}_{l,m} \mathcal{P}^\pm_{\omega-l\delta_p,k-mk_0,q} \Big) = \nonumber \\ & -\sum_n \Big( \gamma_n^{\pm} \mathcal{P}^\pm_{\omega,k-nk_0,q} - \gamma_n^{\mp} \mathcal{P}^\mp_{\omega,k-nk_0,q} \Big) \nonumber\\ & -\frac{q^2}{\hbar^2} D_0 \mathcal{P}^\pm_{\omega,k,q} \,\,, \label{eq:fp:new} \end{align}\tag{14}\] where we have taken advantage of the mentioned periodicity, allowing us to write \[\begin{align} -\frac{\partial U_\pm}{\partial x} = \sum_n \mathcal{F}^{(0)\pm}_n e^{-i n k_0 z}, \tag{15}\\ F_\pm(z,t) = \sum_{l,m} \mathcal{F}^{p\pm}_{l,m} e^{-i ( m k_0 z -l\delta_p t)`}, \tag{16} \end{align}\] and similarly to the transition rates \[\gamma_\pm(z) = \sum_n \gamma^{\pm}_n e^{-i n k_0 z}. \label{eq:gamdef}\tag{17}\] In the problems considered here, there are no force bias, i.e. \[\mathcal{F}^{(0)\pm}_0=\mathcal{F}^{p\pm}_{0,0}=0, \label{eq:nobias}\tag{18}\] and the states are symmetric: \[\gamma_0^\pm=\gamma_0. \label{eq:g0}\tag{19}\]

From (10 ), it is clear the normalization condition (3 ) is translated into \[\mathcal{P}^-_{0,0,0}+\mathcal{P}^+_{0,0,0}=1. \label{eq:norm:fp:new}\tag{20}\]

Once \(\mathcal{P}^\pm_{\omega,k,q}\) is known, the phase space density is obtained back by means of the inverse transform, \[P_{\pm}(x,p,t)=\frac{1}{4\pi^2\hbar}\sum_l \int \!dk \int \!dq \, e^{i (-k x+tl\delta_p )}e^{-i pq/\hbar} \mathcal{P}^\pm_{l\delta_p,k,q}.\] We are specially interested in the marginal probability \(P_{\pm}(x,t)=\int dp P_\pm(x,p,t)\), which can be obtained from \(\mathcal{P}^\pm_{\omega,k,q}\) with \(q=0\), \[P_{\pm}(x,t)= \frac{1}{2\pi}\sum_l \int \!dk \, \mathcal{P}^\pm_{l\delta_p,k,0} e^{-i (k x- t l \delta_p)}. \label{eq:prob:mar}\tag{21}\] The normalization condition (20 ) and the equations of motion (14 ) only requires excitation of a discrete set of possible wave number values, \[k=n k_0 \label{eq:k:discrete}\tag{22}\] where \(n\) is an integer.

Note equation (21 ) expresses the atomic density as a sum of plane waves, atomic density modes, each moving with velocity \((l/n)\delta_p/k_0\). The case \((l,n)=(1,1)\) has been explicitly referred to as a Brillouin-like propagation mode [2], [18], [19], because of the analogy with the resonances produced by light scattering on propagative modes in fluids, such as sound waves.

Assuming, like in (10 ), that the time \(t=0\) is after all the transients have already died out, the atomic current (4 ) can be written in terms of the mode amplitudes as \[\begin{align} \langle v \rangle &=& \frac{\delta_p}{2\pi} \int_{0}^{2\pi/\delta_p} \!\! dt \int \!\! dx \int \!\! dp \, \frac{p}{m_a}[P_+(x,p,t)+P_-(x,p,t)] \nonumber\\ &=& \lim_{q\rightarrow 0}\frac{\hbar \delta_p }{i m_a 2\pi}\frac{\partial}{\partial q} \int_0^{2\pi/\delta_p} \!\! dt \int dx\int dp \, e^{i(k x-\omega t+pq/\hbar)} \nonumber\\ &{ }&\quad \quad\times [P_+(x,p,t)+P_-(z,p,t)] \Big|_{\omega=k=0} \nonumber\\ &=& \frac{\hbar}{i m_a}\lim_{q\rightarrow 0}\left( \frac{\partial \mathcal{P}^+_{0,0,q}}{\partial q}+\frac{\partial \mathcal{P}^-_{0,0,q}}{\partial q} \right). \label{eq:currdef:four} \end{align}\tag{23}\]

The equations (14 ) and (20 ) can be used to find \(\partial \mathcal{P}^\pm_{0,0,0}/\partial q\) in terms of the mode amplitudes \(\mathcal{P}^\pm_{\omega,k,0}\), thus providing the exact contribution of each excited mode to the atomic current, yielding one of the central analytical results of this paper, in which the current (4 ) is written as \[\langle v \rangle =\sum_{l,n} v[l,n] =\sum_{l,n} v^+_{l,n}\mathcal{P}^+_{l\delta_p,nk_0,0} + v^-_{l,n}\mathcal{P}^-_{l\delta_p,nk_0,0}, \label{eq:vmodes}\tag{24}\] with explicit analytical expressions for the coefficients \(v^\pm_{l,n}\) for the specific system at hand, as reported in Sec. 4.

The same procedure can be applied to compute any moment of the velocity distribution, \(\langle v^n \rangle\), with \(n\) a positive integer, since following similar steps as in (23 ) one easily obtains \[\begin{align} \langle v^n \rangle = \left(\frac{\hbar}{i m_a}\right)^n \lim_{q\rightarrow 0}\left( \frac{\partial^n \mathcal{P}^+_{0,0,q}}{\partial q^n}+\frac{\partial^n \mathcal{P}^-_{0,0,q}}{\partial q^n} \right), \label{eq:moment} \end{align}\tag{25}\] with \(q\)-derivatives that can be found in terms of the amplitudes of the atomic waves in a similar fashion as before.

3.1 Basic symmetry↩︎

Due to the vectorial nature of both the driving force \(F_\pm(x,t)\) and the current \(\langle v\rangle\), the inversion transformation \(x\rightarrow-x\) yields an inverted current [28], a fact which is behind the —easy to verify— statement that the solution of the problem with conservative force \[\widehat{\mathcal{F}}^{(0)\pm}_n = -\mathcal{F}^{(0)\pm}_{-n},\] and driving force \[\widehat{\mathcal{F}}^{p\pm}_{l,m}= -\mathcal{F}^{p\pm}_{l,-m}, \label{eq:oppforce}\tag{26}\] is just \[\widehat{\mathcal{P}}^\pm_{\omega,k,q}=\mathcal{P}^\pm_{\omega,-k,-q},\] which, by virtue of (23 ), yields \[\langle \widehat{v}\rangle= -\langle v\rangle. \label{eq:vinver}\tag{27}\] Note if the potential is spatially symmetric, then \[\mathcal{F}^{(0)\pm}_n = -\mathcal{F}^{(0)\pm}_{-n}, \label{eq:spasymm}\tag{28}\] and to produce a current we need a symmetry breaking probe.

3.2 Sketch of the calculation of the current in terms of mode amplitudes↩︎

From (23 ), the calculation of the current is reduced to find an expression for \(\partial \mathcal{P}^\pm_{0,0,0}/\partial q\). However, this is not a trivial task, because the term proportional to \(\partial/\partial q\) in (14 ) vanishes for \(k=0\).

To proceed, we focus in (14 ) for \(\omega=k=0\), and Taylor expand the dependency on \(q\), having for small \(q\), \[\mathcal{P}^\pm_{0,0,q} = \mathcal{P}^\pm_{0,0,0} + q \frac{\partial \mathcal{P}^\pm_{0,0,0}}{\partial q}+ \frac{q^2}{2} \frac{\partial^2 \mathcal{P}^\pm_{0,0,0}}{\partial q^2}+\ldots \label{eq:taylor}\tag{29}\] After summing the \(-\) and \(+\) forms of Eq. (14 ) for \(\omega=k=0\) and inserting the expansion (29 ), we find, in each order, with \(n\) a non-negative integer, \[\begin{align} &0& = \sum_n \left( \mathcal{F}^{(0)+}_{n'} \mathcal{P}^+_{0,-n'k_0,0} + \mathcal{F}^{(0)-}_{n'} \mathcal{P}^-_{0,-n'k_0,0} \right) \nonumber\\ & & + \sum_{l,m} \left( \mathcal{F}^{p+}_{l,m} \mathcal{P}^+_{-l\delta_p,-mk_0,0} + \mathcal{F}^{p-}_{l,m} \mathcal{P}^-_{-l\delta_p,-mk_0,q} \right), \\ & & \frac{(n+1)D_0}{i\hbar} \left( \frac{\partial^n \mathcal{P}^+_{0,0,0}}{\partial q^n}+\frac{\partial^n \mathcal{P}^-_{0,0,0}}{\partial q^n} \right) = \nonumber \\ & & \sum_{n'} \left( \mathcal{F}^{(0)+}_{n'} \frac{\partial^{n+1} \mathcal{P}^+_{0,-n'k_0,0}}{\partial q^{n+1}} + \mathcal{F}^{(0)-}_{n'} \frac{\partial^{n+1} \mathcal{P}^-_{0,-n'k_0,0}}{\partial q^{n+1}} \right) \nonumber\\ & & + \sum_{l,m} \Big( \mathcal{F}^{p+}_{l,m} \frac{\partial^{n+1} \mathcal{P}^+_{-l\delta_p,-mk_0,0}}{\partial q^{n+1}} \nonumber\\ & & \quad \quad \quad + \mathcal{F}^{p-}_{l,m} \frac{\partial^{n+1} \mathcal{P}^-_{-l\delta_p,-mk_0,0}}{\partial q^{n+1}} \Big). \label{eq:curr:sol} \nonumber\\ \end{align}\tag{30}\] Eq. (30 ) with \(n=1\) facilitates the task to calculate the current, because, by means of (23 ), it provides a useful expression in terms of derivatives with respect to \(q\) of amplitudes with either \(k\ne0\) or \(\omega\ne0\). These derivatives can be calculated by using the following equations, obtained by using the Taylor expansion \[\frac{ \partial^{n'} \mathcal{P}^\pm_{\omega,k,q}}{\partial q^{n'}} =\frac{\partial^{n'} \mathcal{P}^\pm_{\omega,k,0} }{\partial q^{n'}} + q \frac{\partial^{1+n'} \mathcal{P}^\pm_{\omega,k,0}}{\partial q^{1+n'}}+ \frac{q^2}{2} \frac{\partial^{2+n'} \mathcal{P}^\pm_{\omega,k,0}}{\partial q^{2+n'}}+\ldots \label{eq:taylor2}\tag{31}\] with \(n'\) a non-negative integer, in (14 ), \[\begin{align} &\frac{\hbar k}{m_a}\frac{\partial \mathcal{P}^\pm_{\omega,k,0} }{\partial q} &= i\omega \mathcal{P}^\pm_{\omega,k,0} \nonumber \\ & &+ \sum_n \left( \gamma^\pm_n \mathcal{P}^\pm_{\omega,k-nk_0,0} - \gamma^\mp_n \mathcal{P}^\mp_{\omega,k-nk_0,0} \right), \tag{32}\\ &\frac{\hbar k}{m_a}\frac{\partial^2 \mathcal{P}^\pm_{\omega,k,0} }{\partial q^2} &= i\omega \frac{\partial \mathcal{P}^\pm_{\omega,k,0}}{\partial q} -\frac{i}{\hbar} \sum_n \mathcal{F}^{(0)\pm}_n \mathcal{P}^\pm_{\omega,k-nk_0,0} \nonumber \\ & & -\frac{i}{\hbar} \sum_{l,m} \mathcal{F}^{p\pm}_{l,m} \mathcal{P}^\pm_{\omega-l\delta_p,k-mk_0,0} \nonumber\\ & &+ \sum_n \left( \gamma^\pm_n \frac{\partial \mathcal{P}^{\pm}_{\omega,k-nk_0,0} }{\partial q} - \gamma^\mp_n \frac{\partial \mathcal{P}^{\mp}_{\omega,k-nk_0,0} }{\partial q}\right). \nonumber\\ \tag{33} \end{align}\] The direct use of (32 )–(33 ) allows us to calculate most terms, except \(\partial \mathcal{P}^{\mp}_{\omega,0,0}/\partial q\), with \(\omega\ne0\). In that case, that first derivative can also be found from (33 ), but it requires linear inversion of a matrix of range \(2\times2\), because similar first derivative terms also appears in the right hand side. Performing that calculation we obtain \[\begin{align} &\frac{\partial \mathcal{P}^\pm_{\omega,0,0} }{\partial q} &= \frac{1}{2\gamma_0 i \omega-\omega^2}\Bigg[ \frac{\gamma_0 i}{\hbar} \Big[ \nonumber\\ & & \sum_n \left( \mathcal{F}^{(0)+}_n \mathcal{P}^+_{\omega,-nk_0,0} + \mathcal{F}^{(0)-}_n \mathcal{P}^-_{\omega,-nk_0,0} \right) \nonumber\\ & &+ \sum_{l,m} \left( \mathcal{F}^{p+}_{l,m} \mathcal{P}^+_{\omega-l\delta_p,-mk_0,0} + \mathcal{F}^{p-}_{l,m} \mathcal{P}^-_{\omega-l\delta_p,-m k_0,0} \right) \Big] \nonumber \\ & & -\frac{\omega}{\hbar} \Big[ \sum_n \mathcal{F}^{(0)\pm}_n \mathcal{P}^\pm_{\omega,k-nk_0,0} +\sum_{l,m} \mathcal{F}^{p\pm}_{l,m}\mathcal{P}^\pm_{\omega-l\delta_p,k-mk_0,0} \Big] \nonumber \\ & &- i\omega \sum_{n\ne0} \left( \gamma^\pm_n \frac{\partial \mathcal{P}^\pm_{\omega,-nk_0,0}}{\partial q} - \gamma^\mp_n \frac{\partial \mathcal{P}^\mp_{\omega,-nk_0,0}}{\partial q} \right) \Bigg]. \label{eq:fp95sol:dp} \end{align}\tag{34}\]

Applying the above steps, it is then straightforward to arrive to an explicit expression for the current \(\langle v\rangle\) —or an arbitrary moment \(\langle v^n\rangle\) using similar steps— in the form (24 ) for a particular system. The application to a system with Sisyphus cooling, such as those presented in Sec. 2 reveals the presence of the regime of infinity density, as discussed in Sec. 5.

3.3 Calculation of the mode amplitudes in the computer simulation↩︎

Numerical solutions of the equation (2 ) are obtained by generating a large number of individual atomic trajectories [23] using a stochastic algorithm [29].

For convenience, let us denote the atomic mode amplitudes for \(q=0\) as \[P_\pm[l,n]=\mathcal{P}^\pm_{l\delta_p,n k_0,0}.\] They are the Fourier transform, Eq. (10 ), of the atomic density \(P_\pm(x,t)\), \[\begin{align} P_\pm[l,n]=\frac{\delta_p}{2\pi}\int_0^{2\pi/\delta_p}\!\!\! dt \, e^{-i l \delta_p t}\int \!dx \, e^{i n k_0 x} P_{\pm}(x,t). \nonumber \\ = \lim_{l'\rightarrow\infty}\frac{\delta_p}{2\pi l'}\int_0^{2\pi l'/\delta_p}\!\!\! dt \, e^{-i l \delta_p t}\int \!dx \, e^{i n k_0 x} P_{\pm}(x,t). \nonumber \\ \label{eq:pfourier:def2} \end{align}\tag{35}\] The atomic density can be expressed in terms of atomic trajectories \(x^{\pm}_j(t)\) using the Dirac delta function \[P_{\pm}(x,t)= \frac{1}{N}\sum_{j=1}^N \delta[x-x^{\pm}_j(t)] . \label{eq:dens:dirac}\tag{36}\] Thus, inserting (36 ) into (35 ) yields \[P_\pm[l,n]= \lim_{l'\rightarrow\infty} \frac{\delta_p}{2\pi l' N}\sum_{j=1}^{N}\int_0^{2\pi l'/\delta_p}\!\!\! dt \, e^{i [n k_0 x^{\pm}_j(t) - l \delta_p t] }. \label{eq:pfourier:def3}\tag{37}\] Since in the simulation individual atomic trajectories are generated, Eq. (37 ) can be readily implemented to compute the atomic mode coefficients, involving an average over many atomic trajectories and many probe periods. At a particular time, the atom is in a definite state \(+\) or \(-\), and thus, it only contributes to the corresponding density via (37 ).

4 Application to the system models↩︎

Our goal is to express the current as (24 ), i.e. an expansion of the atomic mode amplitudes \(P_\pm[l,n]\).

In the system 1D lin\(\perp\)lin, the force and rate constants (15 ) and (17 ), defining the optical potential (1 ) and transition rates (6 ), respectively, are given by \[\begin{align} \mathcal{F}^{(0)+}_{\pm1}=\mp\frac{F_0}{2i}, \quad \mathcal{F}^{(0)-}_{\pm1}=-\mathcal{F}^{(0)+}_{\pm 1} \tag{38}\\ \gamma_{\pm 0}^\pm=g_0, \quad \gamma_{\pm1}^+= \frac{g_1}{2}, \quad \gamma_{\pm1}^-=-\gamma_{\pm1}^+. \tag{39} \end{align}\] with \(F_0=k_0 U_0/2\).

In the system 3D-lin\(\perp\)lin, Eqs. (8 )–(9 ), in addition to (38 )–(39 ), there are the following coefficients \[\begin{align} \mathcal{F}^{(0)+}_{\pm2}=\pm\frac{F_0}{2i}, \quad \mathcal{F}^{(0)-}_{\pm2}=\mathcal{F}^{(0)+}_{\pm 2} \tag{40}\\ \gamma_{\pm2}^+= \frac{g_2}{2}, \quad \gamma_{\pm2}^-=\gamma_{\pm2}^+. \tag{41} \end{align}\]

In both cases, the probe addition to the potential has the same form, and the driving force coefficients (16 ) are given by \[\mathcal{F}^{p\pm}_{1, 1} = -\frac{F_p e^{-i\phi_p}}{2i}, \quad \mathcal{F}^{p\pm}_{-1, -1} = \frac{F_p e^{+i\phi_p}}{2i}, \label{eq:fp:v1b}\tag{42}\] where \(F_p=U_0\varepsilon_p k_0/2\).

It can be checked that the solution, of both systems, possesses the following atomic state symmetry \[\mathcal{P}^-_{l\delta_p,n k_0,q}= (-1)^{n-l} \mathcal{P}^+_{l\delta_p,n k_0,q}. \label{eq:pminus:b}\tag{43}\] The proof involves verifying that (43 ) provides indeed a valid solution of (14 )—because of the specific way the sign of each term is changed in each atomic state—and invoking the uniqueness of the solution. This general, simple relation for the mode amplitudes \(\mathcal{P}^-\) in terms of the coefficients in the other state, \(\mathcal{P}^+\) allows us to focus only on the densities of one atomic sub-level.

4.1 System 1D lin\(\perp\)lin↩︎

By following the steps presented in Sec. 3.2, we obtain the following expansion of the current in terms of mode amplitudes \[\begin{align} \langle v \rangle_\mathrm{1D} = \frac{m_a}{m_a F_0 g_1 - 2 D_0 k_0} \Bigg[ -\frac{\text{Im}\left[P_+[0,1]\right] 8 F_0 g_0^2}{k_0} \nonumber \\ +\frac{\text{Im}\left[P_+[0,2]\right] F_0 \left(-4 g_0 g_1+F_0 k_0/m_a\right)}{k_0} \nonumber \\ +\frac{\text{Im}\left[e^{i\phi_p}P_+[1,2]\right] 2 F_0 F_p }{m_a} \nonumber \\ +\frac{ \text{Im}\left[e^{i\phi_p}P_+[1,1]\right] 2 F_p \delta_p^2}{k_0} \nonumber \\ +\frac{\text{Im}\left[e^{i2\phi_p}P_+[2,2]\right] F_p^2 }{m_a} \Bigg]. \label{eq:opcionB95kp95ne95k0} \end{align}\tag{44}\] Note the only contribution that depends explicitly on the probe frequency \(\delta_p\) is the Brillouin mode \((l,n)=(1,1)\), the same propagating mode as the probe potential, all other modes depend explicitly only on the optical potential force \(F_0\), probe amplitude \(F_p\), or transition rates \(g_0\) and \(g_1\).

All terms in (44 ) share a common denominator, \(m_a F_0 g_1 - 2 D_0 k_0\), which produce a singularity when is zero, and it will be discussed in Sec. 5.

The expression (44 ) is validated numerically in Fig. 1, where we depict the contribution of each mode, \(v[l,n]\), after numerically calculating the mode amplitudes \(P_+[l,n]\) in the simulations using (37 ). As a test, we also show the sum of all contributions and compare it with the current calculated directly in the simulation from its definition (4 ), confirming the analytical calculation. Reduced units are defined, in all simulations reported in thus paper, such that \(m_a=\hbar=k_l=1\).

Previous research [18], [19], [25], [26] has focused on one propagating mode, the mode \((l,n)=(1,1)\), which travels at speed \(\delta_p/k_0\) and is obviously excited by the probe potential, having the same frequency and wave number. Here this mode is confirmed to provide the most important contribution to the current. However, other modes are also excited, and they even produce a noticeable contribution, as illustrated in Fig. 1. The atomic mode \((1,2)\), which is also propagating to the right, though at the smaller speed \(\delta_p/(2k_0)\), as well as the non-propagating mode \((0,2)\), produce contributions to the current in the opposite direction. The total current is, however, dominated by the mode \((1,1)\) and goes in the same direction as the propagating probe potential.

4.2 System 3D-lin\(\perp\)lin↩︎

In the system 3D–lin\(\perp\)lin, the additional terms (40 )–(41 ) produce more mode contributions to the current, appearing both more contributing terms for the previously discussed modes, and new contributions from other modes. \[\begin{align} \langle v \rangle_\mathrm{3D} = \frac{m_a}{m_a F_0 g_1 - 2 D_0 k_0} \Bigg[ \nonumber \\ -\frac{\text{Im}\left[P_+[0,1,0]\right] F_0\left( 8 g_0^2 {\color{blue}-4g_2^2/3+F_0 k_0/(2m_a)} \right) }{k_0} \nonumber \\ +\frac{\text{Im}\left[P_+[0,2,0]\right] F_0 \left(-4 g_0 g_1 {\color{blue}-8g_1 g_2/3} +F_0 k_0/m_a\right)}{k_0} \nonumber \\ \color{blue}+\frac{\text{Im}\left[P_+[0,3,0]\right] F_0 \left(-16 g_0 g_1/3 - 2 g_2^2 -3F_0 k_0/(2m_a)\right)}{k_0} \nonumber \\ \color{blue}+\frac{\text{Im}\left[P_+[0,4,0]\right] F_0 \left(-2 g_1 g_2/3 +F_0 k_0/(2m_a)\right)}{k_0} \nonumber \\ \color{blue}-\frac{\text{Im}\left[P_+[0,5,0]\right] F_0 2 g_2^2 }{3 k_0} +\frac{\text{Im}\left[e^{i\phi_p}P_+[1,1,1]\right] 2 F_0 F_p }{m_a} \nonumber \\ \color{blue} -\frac{\text{Re}\left[e^{i\phi_p}P_+[1,-2,1]\right] F_0 F_p }{2 m_a} \color{blue} +\frac{\text{Re}\left[e^{i\phi_p}P_+[1,2,1]\right] 3 F_0 F_p }{2 m_a} \nonumber \\ +\frac{ \text{Im}\left[e^{i\phi_p}P_+[1,0,1]\right] 2 F_p \delta_p^2}{k_0} +\frac{\text{Im}\left[e^{i2\phi_p}P_+[2,0,2]\right] F_p^2 }{m_a} \Bigg]. \label{eq:3Dlin} \nonumber \\ \end{align}\tag{45}\]

Like in the previous system, the biggest contributor to the current is the propagating mode \((1,1)\), producing movement in the same direction as the propagating probe potential, here with the help of non-propagating mode \((0,1)\), which is also present in the system 1D lin\(\perp\)lin —also producing there a current in the same direction, but with a much less appreciable contribution. New contributions from the non-propagating modes \((0,3)\), \((0,4)\) and \((0,5)\)—due to the extra mode-connecting terms in (8 ) and (9 )—are observed, all pointing to the opposite direction of propagation. A similar feature is observed in the new modes \((1,-1)\) and \((1,3)\), which together with the propagating mode \((1,2)\) produce current in the opposite direction. The overall current is much smaller—about four times smaller—than the one produced alone by the Brillouin mode \((1,1)\) due to those counter-propagating contributions.

5 Into the regime of infinite density↩︎

The fact that the optical potential force \(\mathcal{F}_{\pm1}^{(0)\pm}\), Eq. (38 ), inverts its sign when the atomic state is changed is crucially coupled to the similar state-alternating sign in the transition rate \(\gamma_{\pm1}^\pm\), Eq. (39 ), being responsible for the celebrated Sisyphus cooling—a kind of dissipative mechanism due to the transitions taking place with higher probability when the atoms are at the top of the optical potential barriers, thus the atom spends most of the time climbing hills, like the punishment inflicted in the king Sisyphus of ancient Greek mythology. However, it is known the potential should not be too shallow [21], or the dynamics may exploit in the form of superdiffusion [8] and non-normalizable probability densities [13], [15]. Indeed, the calculations presented in Sec. 4 indicate that the current diverges to infinity, in both studied systems, when \[F_0 g_1-2D_0 k_0/m_a = 0, \label{eq:curr:inf}\tag{46}\] because the left hand side of (46 ) appears as a common factor in the denominator in the formulas for the current, (44 ) and (45 ). It is a consequence of the appearance of a term proportional to \(\partial \mathcal{P}^+_{0,0,0}/\partial q +\partial \mathcal{P}^-_{0,0,0}/\partial q\) in the right hand side of (30 ) (with \(n=1\)), which has to be taken to the left hand side in order to complete the calculation.

In fact, the same thing happens in any moment \(\langle v^n\rangle\). Using (25 ) and (30 ), and looking for the terms in the right hand side that yield a contribution containing \(\partial^n \mathcal{P}^+_{0,0,0}/\partial q^n +\partial^n \mathcal{P}^-_{0,0,0}/\partial q^n\), i.e. the terms proportional to \(\mathcal{F}^{(0)\pm}_{+1}\gamma^\pm_{-1}\) and \(\mathcal{F}^{(0)\pm}_{-1}\gamma^\pm_{+1}\), yields the common factor, \([1-F_0 g_1 m_a/(D_0 k_0(n+1)]^{-1}\), and thus the following singularity condition, \[F_0 g_1-(1+n)D_0 k_0/m_a = 0, \label{eq:moment:inf}\tag{47}\] which in terms of the potential depth is written as \[U_0 = \frac{2(1+n) D_0}{g_1 m_a}. \label{eq:U0:inf}\tag{48}\] Though this result has been obtained from the semiclassical equations (2 ), which do not explicitly contain all the diffusive terms predicted by the semiclassical derivation, that is, the terms \(\partial^2 D_{\mp\pm} P_{\mp}/\partial p^2\) and the spatial dependence of the noise coefficient, the same analysis in the more complex system model shows that the exact same result (47 ) is obtained if one considers \(D_0\) as the spatial average of all diffussive terms.

This singularity can be identified with the threshold for the transition into an anomalous regime, in which the \(n\) moment of the velocity distribution, \(\langle v(t)^n\rangle\), diverges for long times. The result (48 ) is in perfect agreement with the prediction [14], [15] for the system 1D lin\(\perp\)lin obtained from a simplified, approximate Fokker-Planck equation [21] for shallow lattices, when the noise stregth \(D_0\) contains the spatial average of both diffusive terms \(D_{\mp\pm}\) and \(D_{\pm\pm}\) [23], i.e. \(D_0=(35+6)\hbar^2 k_l^2\Gamma'/90\). This is remarkable, as our method uses the full semiclassical equations (2 ), and hence, unlike the approximate Fokker-Planck equation, it does explicitly takes into account the microscopic origin of Sisyphus cooling, including the periodicity of the lattice potential and transitions between the atomic sublevels.

The only assumption made by the current method is the existence of the Taylor expansion (31 ), which is numerically confirmed by the simulations presented in the paper in the regime above the singularity. In contrast, in shallow lattices—i.e. with potential wells \(U_0\) below the threshold given by (48 )—this assumption does not hold, giving rise to the regime of infinite densities, with diverging moments and non-normalizable densities.

In the system 3D-lin\(\perp\)lin, to accurately cover the shallow potential regime, \(D_0\) should also include the spatial average the diffusive term \(D_{\mp\pm}\), thus \(D_0=(1+1/5)5\hbar^2k_0^2\Gamma_s/18\).

6 Conclusions↩︎

Starting from the semiclassical equations for the atomic phase densities, we have presented an exact method to calculate the precise contribution to the current of the excited atomic density waves in cold atoms under dissipative optical lattices. Explicit analyitcal expressions are provided for the popular system 1D lin\(\perp\)lin under a simple—but realizable—symmetry breaking perturbation, as well as the one-dimensional model that results from focusing in a specific direction in the 3D-lin\(\perp\)lin system.

The analytical results are validated with numerical simulations of the stochastic atomic trajectories associated with the semiclassical equations. They show that several atomic modes, not just the one associated with the perturbation, contribute relevantly to the directed motion.

Additionally, the analytical solution predicts a singularity for each moment of the velocity distribution, which is identified with the threshold for the transition to the regime of infinite density. The threshold values of the present work are identical to the ones obtained from an approximated Fokker-Planck equation of the semiclassical equations, derived by spatial and Zeeman level averaging the latter, in the system 1D lin\(\perp\)lin, thus providing a solid ground to previous analytical results.

Finally, it is worth mentioning that the proposed method can easily be applied to more complex systems and driving perturbations.

References↩︎