Absorption imaging of a quasi-two-dimensional gas: a multiple scattering analysis

We use the standard description of the quantized electromagnetic field in the Coulomb gauge [22], and choose periodic boundary conditions in the cubic-shaped quantization volume ${\cal V}=L_x L_yL_z$. We denote by a_q,s the destruction operator of a photon with wave vector q and polarization s (s⊥q). The Hamiltonian of the quantized field is

$$\begin{equation} H_{{\rm F}}=\sum_{\boldsymbol q,\boldsymbol s} \hbar cq\, a_{\boldsymbol q,\boldsymbol s}^\dagger a_{\boldsymbol q,\boldsymbol s} \end{equation} \tag{ 3 }$$

and the transverse electric field operator reads E(r) = E⁽⁺⁾(r) + E⁽⁻⁾(r) with

$$\begin{equation} \boldsymbol E^{(+)}(\boldsymbol r)={\rm i}\sum_{\boldsymbol q,\boldsymbol s}\sqrt{\frac{\hbar c q}{2\varepsilon_0 {\cal V}}} \,a_{\boldsymbol q,\boldsymbol s}\,{\rm e}^{{\rm i}\boldsymbol q\cdot \boldsymbol r}\,\boldsymbol s, \end{equation} \tag{ 4 }$$

and E⁽⁻⁾(r) = (E⁽⁺⁾(r))^†. The wave vectors q are quantized in the volume ${\cal V}$ as q_i = 2πn_i/L_i, i = x,y,z, where n_i is a positive or negative integer.

We consider a collection of N identical atoms at rest in positions r_j, j = 1,...,N. We model the atomic resonance transition by a two-level system with a ground state |g〉 with angular momentum J_g = 0 and an excited level of angular momentum J_e = 1. We choose as a basis set for the excited manifold the three Zeeman sublevels |e_α〉, α = x,y,z, where |e_α〉 is the eigenstate with eigenvalue 0 of the component J_α of the atomic angular momentum operator. We denote by ℏω₀ the energy difference between e and g. The atomic Hamiltonian is thus (up to a constant)

$$\begin{equation} H_{{\rm A}}=\sum_{j=1}^N \sum_{\alpha=x,y,z} \hbar \omega_0 \, |j:e_\alpha\rangle \langle j:e_\alpha| . \end{equation} \tag{ 5 }$$

The restriction to a two-level approximation is legitimate if the detuning Δ between the probe and the atomic frequencies is much smaller than ω₀. The modelling of this transition by a $J_g=0\leftrightarrow J_e=1$ transition leads to a relatively simple algebra. The transitions that are used for absorption imaging in real experiments often involve more Zeeman states ($J_g=2 \leftrightarrow J_e=3$ for Rb atoms in [7, 8]), but are more complex to handle [23, 24], and they are thus beyond the scope of this paper. However, we believe that the most salient features of multiple scattering and resonant van der Waals interactions are captured by our simple level scheme.

We treat the atom–light interaction using the electric dipole approximation (length gauge), which is legitimate since the resonance wavelength of the atoms λ₀ is much larger than the atomic size. We write the atom–light coupling as

$$\begin{equation} V=-\sum _j \boldsymbol D_j \cdot \boldsymbol E(\boldsymbol r_j) , \end{equation} \tag{ 6 }$$

where D_j is the dipole operator for the atom j. We will use the RWA, which consists in keeping only the resonant terms in the coupling:

$$\begin{equation} V\approx -\sum _j \boldsymbol D_j^{(+)} \cdot \boldsymbol E^{(+)}(\boldsymbol r_j) + {\rm h.c.}, \end{equation} \tag{ 7 }$$

where h.c. stands for the Hermitian conjugate. Here D⁽⁺⁾_j represents the raising part of the dipole operator for atom j:

$$\begin{equation} \boldsymbol D_j^{(+)}= d \sum _{\alpha=x,y,z} |j:e_\alpha\rangle \langle j:g|\, \boldsymbol {\hat u_\alpha}, \end{equation} \tag{ 8 }$$

where d is the electric dipole associated with the g–e transition and $\boldsymbol {\hat u_\alpha }$ is a unit vector in the direction α.

When a single atom is coupled to the electromagnetic field, this coupling results in the modification of the resonance frequency (the Lamb shift) and in the fact that the excited state e acquires a non-zero width Γ

$$\begin{equation} \Gamma= \frac{d^2\omega_0^3}{3\pi\varepsilon_0\hbar c^3}. \end{equation} \tag{ 9 }$$

For simplicity we will incorporate the Lamb shift into the definition of ω₀. Note that the proper calculation for this shift requires that one goes beyond the two-level and the RWAs. The linewidth Γ, on the other hand, can be calculated from the above expressions for V using the Fermi golden rule.

The RWA provides a very significant simplification of the treatment of the atom–light coupling, in the sense that the total number of excitations is a conserved quantity. The annihilation (resp. creation) of a photon is always associated with the transition of one of the N atoms from g to e (resp. from e to g). This would not be the case if the non-resonant terms of the electric dipole coupling D⁽⁺⁾_i·E⁽⁻⁾ and D⁽⁻⁾_i·E⁽⁺⁾ were also taken into account. For a single atom, the small parameter associated with the RWA is Δ/ω₀, which is, in practice, in the range 10⁻⁶–10⁻⁹. For a collection of atoms, the RWA is still an excellent approximation if the shifts introduced by resonant van der Waals interactions are small compared to atomic frequencies. In practice, this requires that the interatomic distances remain very large compared to the atom size (the Bohr radius).

Formally, the use of the electric dipole interaction implies to add to the Hamiltonian an additional contact term between the dipoles (see, e.g., [12, 25]). This term will play no role in our numerical simulations because we will surround the position of each atom by a small excluded volume, which mimics the short-range repulsive interaction between atoms. We checked that the results of our numerical calculations (see section 4) do not depend on the size of the excluded volume, and we can safely omit the additional contact term in this work.

In this section, we assume that the incident field is prepared in a coherent state corresponding to a monochromatic plane wave E_L  e^i(kz−ωt). We choose the polarization to be linear and parallel to the x-axis ($\boldsymbol \epsilon =\boldsymbol {\hat u_x}$). Since we consider a $J_g=0 \leftrightarrow J_e=1$ transition, this choice does not play a significant role and we checked that we recover essentially the same results with a circular polarization. Note that the situation would be different for an atomic transition with larger J_g and J_e since optical pumping processes would then depend crucially on the polarization of the probe laser.

The amplitude E_L is supposed to be small enough that the steady-state populations of the excited states e_j,α are small compared to unity. This ensures that the response of each atomic dipole is linear in E_L; this approximation is valid when the Rabi frequency dE_L/ℏ is small compared to the natural width Γ or the detuning Δ.

Using the atom–light coupling (6), the equations of motion for the annihilation operators a_q,s in the Heisenberg picture read

$$\begin{equation} \dot a_{\boldsymbol q,\boldsymbol s}(t)=-{\rm i} \,cq\, a_{\boldsymbol q,\boldsymbol s}(t)+\sqrt{\frac{cq}{2\hbar \varepsilon_0 {\cal V}}} \sum_{j'} \boldsymbol s^*\cdot \boldsymbol D_{j'}(t)\,{\rm e}^{-{\rm i}\boldsymbol q\cdot \boldsymbol r_{j'}}. \end{equation} \tag{ 10 }$$

This equation can be integrated between the initial time t₀ and the time t, and the result can be injected in the expression for the transverse field to provide its value at any point r:

$$\begin{eqnarray} &&\fl E_{\alpha}(\boldsymbol r,t)=E_{{\rm free},\alpha}(\boldsymbol r,t)+\sum_{j',\alpha'}\sum_{\boldsymbol q,\boldsymbol s} \int_0^{t-t_0}{\rm d}\tau\; \frac{cq}{2\varepsilon_0 {\cal V}} \left[ {\rm i} D_{j',\alpha'}(t-\tau)\,{\rm e}^{{\rm i}\boldsymbol q\cdot (\boldsymbol r-\boldsymbol r_{j'})-{\rm i}cq \tau}\,s_\alpha s^*_{\alpha'} + {\rm h.c.}\right],\nonumber\\ \end{eqnarray} \tag{ 11 }$$

where E_free stands for the value obtained in the absence of atoms. We now take the quantum average of this set of equations. In the steady-state regime the expectation value of the dipole operator D_j(t) can be written as d_j e^−iωt + c.c., and the average of E_free(r,t) is the incident field E_L  e^i(kz−ωt) + c.c. We denote the average value of the transverse field operator in r as $\langle \boldsymbol E(\boldsymbol r,t)\rangle =\bar {\boldsymbol E} (\boldsymbol r)\,{\rm e}^{-{\rm i}\omega t}+{\rm c.c.}$, and we obtain after some algebra (see, e.g., [12, 27])

$$\begin{equation} \bar{E}_{\alpha}(\boldsymbol r)=E_{{\rm L}}\,\epsilon_{\alpha}{\rm e}^{{\rm i} kz}+\frac{k^3}{6\pi\varepsilon_0}\sum_{j',\alpha'}g_{\alpha,\alpha'}(\boldsymbol u_{j'})\,d_{j',\alpha'}, \end{equation} \noindent \tag{ 12 }$$

where we set u_j = k(r − r_j) (with k ≈ k₀),

$$\begin{equation} g_{\alpha,\alpha'}(\boldsymbol u)=\delta_{\alpha,\alpha'}h_1(u)+\frac{u_\alpha u_{\alpha'}}{u^2}h_2(u) \end{equation} \tag{ 13 }$$

and

$$\begin{equation} h_1(u)=\frac{3}{2}\,\frac{{\rm e}^{{\rm i}u}}{u^3}(u^2+{\rm i}u-1),\quad h_2(u)=\frac{3}{2}\,\frac{{\rm e}^{{\rm i}u}}{u^3}(-u^2-3{\rm i}u+3) . \end{equation} \noindent \tag{ 14 }$$

The function g_α,α'(kr) is identical to the one appearing in classical electrodynamics [28], when calculating the field radiated in r by a dipole located at the origin.

We proceed similarly for the equations of motion for the dipole operators D⁽⁻⁾_j and take their average value in steady state. The result can be put in the form [12]

$$\begin{equation} (\delta+i) d_{j,\alpha} + \sum_{j'\neq j,\;\alpha'} g_{\alpha,\alpha'}(\boldsymbol u_{jj'}) d_{j',\alpha'}= - \frac{6\pi\varepsilon_0}{k^3}\, E_{{\rm L}}\,\epsilon_{\alpha}{\rm e}^{{\rm i} kz_j}, \end{equation} \noindent \tag{ 15 }$$

where the reduced detuning δ = 2Δ/Γ has been defined in equation (2) and u_j,j' = k(r_j' − r_j). This can be written with matrix notation

$$\begin{equation} [M] |X\rangle=|Y\rangle, \end{equation} \noindent \tag{ 16 }$$

where the 3N vectors |X〉 and |Y 〉 are defined by

$$\begin{equation} X_{j,\alpha}=-\frac{k^3}{6\pi \epsilon_0 E_L}d_{j,\alpha}, \quad Y_{j,\alpha}=\epsilon_{\alpha}\,{\rm e}^{{\rm i}k z_{j}}, \end{equation} \noindent \tag{ 17 }$$

and where the complex symmetric matrix [M] has its diagonal coefficients equal to δ + i and its off-diagonal coefficients (for j ≠ j') given by g_α,α'(u_jj'). This matrix belongs to the general class of Euclidean matrices [29], for which the (i,j) element can be written as a function F(r_i,r_j) of points r_i in the Euclidean space. The spectral properties of these matrices for a random distribution of the r_i's (as will be the case in this work, see section 4) have been studied in [29–32].

Equation (15) has a simple physical interpretation: in the steady state each dipole d_j is driven by the sum of the incident field E_L and the field radiated by all the other dipoles. This set of 3N equations was first introduced by L L Foldy in [33] who named it, together with equation (12), 'the fundamental equations of multiple scattering'. Indeed for a given incident field, the solution of (16) provides the value of each dipole d_j, which can then be injected in (12) to obtain the value of the total field at any point in space.

From the expression of the average value of the dipoles, we now extract the absorption coefficient of the probe beam and the optical density of the gas. We suppose that the N atoms are uniformly spread in a cylinder of radius R along the z-axis and located between z = −ℓ_z/2 and z = ℓ_z/2. We can consider two experimental setups to address this problem. The first one, presented in figure 1(a), consists in measuring after the atomic sample the total light intensity with the same momentum $\boldsymbol k=k \boldsymbol {\hat u_z}$ as the incident probe beam. This can be achieved by placing a lens with the same size as the atomic sample, in the plane z = ℓ' > ℓ_z/2 just after the sample. The light field at the focal point of the lens F gives the desired attenuation coefficient. We refer to this method as 'global', since the field E(F) provides information over the whole atomic cloud. One can also use the setup sketched in figure 1(b), which forms an image of the atom slab on a camera and provides a 'local' measurement of the absorption coefficient. In real experiments local measurements are often favoured because trapped atomic sample are non-homogeneous and it is desirable to access the spatial distribution of the particles. However, for our geometry with a uniform density of scatterers, spatial information on the absorption of the probe beam is not relevant. Therefore we only present the formalism for global measurements, which is simpler to derive and leads to slightly more general expressions. We checked numerically that we obtained very similar results when we modelled the local procedure.

Zoom In Zoom Out Reset image size

We assume that the lens in figure 1(a) operates in the paraxial regime, i.e. its focal length f is much larger than its radius R. We relate the field at the image focal point of the lens to the field in the plane z = ℓ' just before the lens:

$$\begin{equation} \boldsymbol E(F) =-\frac{{\rm i} \,{\rm e}^{{\rm i}kf}}{\lambda_0 f} \int_{\cal L} \boldsymbol E (x,y,\ell')\;{\rm d}x\,{\rm d}y , \end{equation} \tag{ 18 }$$

where the integral runs over the lens area. Since the incident probe beam is supposed to be linearly polarized along x, we calculate the x component of the field in F. Plugging the value of the field given in equations (12) and (17), we obtain the transmission coefficient

$$\begin{equation} {\cal T}\equiv \frac{\left. E_x(F)\right|_{{\rm with\; atoms}}}{\left.E_x(F)\right|_{{\rm no\; atom}}} = 1-\frac{{\rm e}^{-{\rm i}k\ell'}}{\pi R^2}\sum_{j,\alpha} X_{j,\alpha}\int_{\cal L}g_{x,\alpha}[k(\boldsymbol r-\boldsymbol r_j)]\;{\rm d}x\,{\rm d}y. \end{equation} \tag{ 19 }$$

This result can be simplified in the limit of a large lens by using an approximated value for the integral appearing in (19). We suppose that kℓ' ≫ 1 so that the dominant part in g_x,α is the e^iu/u contribution to h₁. More precisely, the domain in the lens plane contributing to the integral for the dipole j is essentially a disc of radius $\sqrt {\lambda (\ell '-z_j)}\sim \sqrt {\lambda \ell '}$ centred on (x_j,y_j). When this small disc is entirely included in the lens aperture, i.e. the larger disc of radius R centred on x = y = 0, we obtain

$$\begin{equation} \int_{\cal L}g_{x,\alpha}[k(\boldsymbol r-\boldsymbol r_j)]\;{\rm d}x\,{\rm d}y \approx \frac{3{\rm i}\pi}{k^2}\delta_{x,\alpha}\,{\rm e}^{{\rm i}k(\ell'-z_j)}. \end{equation} \tag{ 20 }$$

We use the result (20) for all atoms, which amounts to neglecting edge effects for the dipoles located at the border of the lens, and we obtain

$$\begin{equation} {\cal T}= 1-\frac{{\rm i}}{2}\sigma_0 n^{({\rm col})} \Pi , \end{equation} \tag{ 21 }$$

with n^(col) = N/πR² and where the coefficient Π is defined by

$$\begin{equation} \Pi=\frac{1}{N}\sum_j X_{j,x}\,{\rm e}^{-{\rm i}kz_j} . \end{equation} \tag{ 22 }$$

This coefficient captures the whole physics of multiple scattering and resonant van der Waals interactions among the N atoms. Indeed one takes into account all possible couplings between the dipoles when solving the 3N × 3N system [M]|X〉 = |Y 〉. Once ${\cal T}$ is known the optical density is obtained from

$$\begin{equation} {\cal D}\equiv \ln \left|{\cal T}\right|^{-2} . \end{equation} \tag{ 23 }$$

As an example, consider the limit of a very sparse sample where multiple scattering does not play a significant role (σ₀n^(col) ≪ 1). All non-diagonal matrix elements in [M] are then negligible and [M] is simply the identity matrix, times i + δ. Each X_j,x solution of the system (16) is equal to e^ikz_j/(i + δ), and we obtain as expected

$$\begin{equation} \sigma_0 n^{({\rm col})}\ll 1: {\cal T}\approx 1-\frac{1}{2(1-{\rm i}\delta)}\sigma_0 n^{({\rm col})},\quad {\cal D}\approx \frac{\sigma_0 n^{({\rm col})}}{1+\delta^2} . \end{equation} \tag{ 24 }$$

In order to study the attenuation of a weak probe beam propagating along the z-axis when it crosses the atomic medium, we can also use quantum scattering theory. The Hamiltonian of the problem is

$$\begin{equation} H=H_0+V,\quad H_0=H_{{\rm A}}+H_{{\rm F}}, \end{equation} \tag{ 25 }$$

and we consider the initial state where all atoms are in their ground state and where a single photon of wave vector $\boldsymbol k=k\boldsymbol {\hat u_z}$ and polarization $\boldsymbol \epsilon =\boldsymbol {\hat u_x}$ is incident on the atomic medium

$$\begin{equation} |\Psi_{{\rm i}}\rangle =|{\cal G}\rangle \otimes |\boldsymbol k, \boldsymbol \epsilon\rangle , \end{equation} \tag{ 26 }$$

with $|{\cal G}\rangle \equiv |1:g,\ 2:g,\ \ldots , N:g\rangle$. The state |Ψ_i〉 is an eigenstate of H₀ with energy ℏω. The interaction of the photon with the atomic medium, described by the coupling V , can be viewed as a collision process during which an arbitrary number of elementary scattering events can take place. Each event starts from a state $|{\cal G}\rangle \otimes |\boldsymbol q, \boldsymbol s\rangle$ and corresponds to

• (i)  
  The absorption of the photon in mode q,s by atom j, which jumps from its ground state |j:g〉 to one of its excited states |j:e_α〉. The state of the system is then

  $$\begin{equation} |{\cal E}_{j,\alpha}\rangle=|1:g,\ \ldots ,j:e_\alpha,\ldots,N:g\rangle \otimes |{\rm vac}\rangle, \end{equation} \tag{ 27 }$$

  where |vac〉 stands for the vacuum state of the electromagnetic field. The subspace spanned by the states $|{\cal E}_{j,\alpha }\rangle$ has dimension 3N.
• (ii)  
  The emission of a photon in the mode (q',s') by atom j, which falls back into its ground state.

Finally, a photon emerges from the atomic sample, and we want to determine the probability amplitude to find this photon in the same mode |k,〉 as the initial one.

The T matrix defined as

$$\begin{equation} T(E)=V+V\frac{1}{E-H+{\rm i}0_+}V, \end{equation} \tag{ 28 }$$

where 0₊ is a small positive number that tends to zero at the end of the calculation, provides a convenient tool to calculate this probability amplitude. Generally,

$$\begin{equation} T_{ if}=\langle \Psi_{ f}|T(E_i)|\Psi_{ i}\rangle \end{equation} \tag{ 29 }$$

gives the probability amplitude to find the system in the final state |Ψ_f〉 after the scattering process. The states |Ψ_i〉 and |Ψ_f〉 are eigenstates of the unperturbed Hamiltonian H₀, with energy E_i. Here we are interested in the element T_ii of the T matrix, corresponding to the choice |Ψ_f〉 = |Ψ_i〉. Using the definition (28), we find

$$\begin{equation} T_{ ii}=\frac{\hbar \omega d^2}{2\varepsilon_0 {\cal V}}\sum_{j,j'} {\rm e}^{{\rm i}k (z_{j}-z_j')} \;\langle {\cal E}_{j',x}| \frac{1}{\hbar \omega-H+{\rm i}0_+} |{\cal E}_{j,x}\rangle . \end{equation} \tag{ 30 }$$

We now have to calculate the (3N) × (3N) matrix elements of the operator 1/(z − H), with z = ℏω + i0₊, entering into (30). We introduce the two orthogonal projectors P and Q, where P projects on the subspace with zero photon, and Q projects on the orthogonal subspace. We thus have

$$\begin{equation} P |{\cal E}_{j,\alpha}\rangle= |{\cal E}_{j,\alpha}\rangle, \quad P|{\cal G}\rangle \otimes |\boldsymbol k, \boldsymbol \epsilon\rangle=0, \end{equation} \tag{ 31 }$$

$$\begin{equation} Q |{\cal E}_{j,\alpha}\rangle= 0, \quad Q|{\cal G}\rangle \otimes |\boldsymbol k, \boldsymbol \epsilon\rangle=|{\cal G}\rangle \otimes |\boldsymbol k, \boldsymbol \epsilon\rangle.\qquad \end{equation} \tag{ 32 }$$

We define the displacement operator

$$\begin{equation} R(z)=V+V\frac{Q}{z-QH_0Q-QVQ}V \end{equation} \tag{ 33 }$$

and use the general result [22]

$$\begin{equation} P\frac{1}{z-H}P=\frac{P}{z-H_{{\rm eff}}}, \end{equation} \tag{ 34 }$$

where the effective Hamiltonian H_eff is

$$\begin{equation} H_{{\rm eff}}=P\left(H_0+R(z) \right)P . \end{equation} \tag{ 35 }$$

For the following calculations, it is convenient to introduce the dimensionless matrix [M] proportional to the denominator of the right-hand side of (34):

$$\begin{equation} [M]_{(j',\alpha'),(j,\alpha)}= \frac{2}{\hbar \Gamma}\langle {\cal E}_{j',\alpha'}| z-H_{{\rm eff}}|{\cal E}_{j,\alpha}\rangle . \end{equation} \tag{ 36 }$$

It is straightforward to check⁵ that for z → ℏω this matrix coincides with the symmetric matrix appearing in (16). Indeed the matrix elements of R(z) are

$$\begin{equation} \langle {\cal E}_{j',\alpha'}| R(z) |{\cal E}_{j,\alpha}\rangle=\frac{\hbar d^2}{2\varepsilon_0 {\cal V}}\sum_{\boldsymbol q,\boldsymbol s} cq\, s_{\alpha}^* s_{\alpha'}\; \frac{{\rm e}^{{\rm i}\boldsymbol q\cdot (\boldsymbol r_{j'}-\boldsymbol r_j)}}{z-\hbar \omega}, \end{equation} \tag{ 37 }$$

which can be calculated explicitly. For j = j', the real part of this expression is the Lamb shift that we reincorporate in the definition of ω₀, and its imaginary part reads

$$\begin{equation} \langle {\cal E}_{j,\alpha'}| R(z) |{\cal E}_{j,\alpha}\rangle=-{\rm i}\frac{\hbar \Gamma}{2}\,\delta_{\alpha,\alpha'}. \end{equation} \tag{ 38 }$$

For j ≠ j', the sum over (q,s) appearing in (37) is the propagator of a photon from an atom in r_j in internal state |e_α〉, to another atom in r_j' in internal state |e_α'〉. This is nothing but (up to a multiplicative coefficient) the expression that we have already introduced for the field radiated in r_j' by a dipole located in r_j:

$$\begin{equation} \langle {\cal E}_{j',\alpha'}| R(z) |{\cal E}_{j,\alpha}\rangle=-\frac{\hbar \Gamma}{2} g_{\alpha,\alpha'}(\boldsymbol u_{j,j'}), \end{equation} \tag{ 39 }$$

where the tensor g_α,α' is defined in equations (13) and (14).

Suppose now that the atoms are uniformly distributed over the transverse area L_xL_y of the quantization volume. We set n^(col) = N/(L_xL_y) and we rewrite the expression (30) of the desired matrix element T_ii as

$$\begin{equation} \frac{ T_{ ii}L_z}{\hbar c}=\frac{1}{2N}\sigma_0 n^{({\rm col})} \;\sum_{j,j'}{\rm e}^{{\rm i} k (z_{j}-z_{j'})} \;[M^{-1}]_{(j,x),(j',x)}=\frac{1}{2}\sigma_0 n^{({\rm col})} \Pi, \end{equation} \tag{ 40 }$$

where the coefficient Π has been defined in (22). The result (40) combined with (21) leads to

$$\begin{equation} {\cal T}=1-{\rm i}\frac{ T_{ ii}L_z}{\hbar c}, \end{equation} \tag{ 41 }$$

which constitutes the 'optical theorem' for our slab geometry, since it relates the attenuation of the probe beam ${\cal T}$ to the forward scattering amplitude T_ii.

The emergence of resonant van der Waals interactions is straightforward in this approach. Let us consider for simplicity the case when only N = 2 atoms are present. The effective Hamiltonian H_eff is a 6 × 6 matrix that can be easily diagonalized and its eigenvectors, with one atom in |e〉 and one in |g〉, form in this particular case an orthogonal basis, although H_eff is a non-Hermitian [34, 35]. For a short distance r between the atoms (kr ≪ 1), the leading term in h₁(u) and h₂(u) is u⁻³ and the energies (real parts of the eigenvalues) of the six eigenstates vary as ∼ ± ℏΓ/(kr)³ (resonant dipole–dipole interaction). The imaginary parts of the eigenvalues, which give the inverse of the radiative lifetime of the states, tend either to Γ or 0 when r → 0, which correspond to the superradiant and subradiant states for a pair of atoms, respectively [36].

For N > 2 the eigenvectors of the non-Hermitian Euclidean matrix H_eff are, in general, non-orthogonal, which complicates the use of standard techniques of spectral theory in this context [31, 32]. More precisely, one could think of solving the linear system (16), or equivalently calculating T_ii in equation (30), by using the expansion of the column vector |Y 〉 defined in equation (17) on the left (|α_j〉) and right (〈β_j|) eigenvectors of H_eff. Then one could insert this expansion into the general expression of the matrix element T_ii, to express it as a sum of the contributions of the various eigenvalues of H_eff. However, the physical discussion based on this approach is made difficult by the fact that since H_eff is non-Hermitian, the {|α_j〉} and the {|β_j〉} bases do not coincide. Hence the weight 〈β_j|Y 〉〈Y |α_j〉 of a given eigenvalue in the sum providing the value of T_ii is not a positive number, and this complicates the interpretation of the result.

As a final comment, we note that our approach does not make use of the quite common concept of self-energy. In the context of the propagation of waves in a random medium [11], the self-energy Σ appears by averaging the relation G = G₀ + G₀V G over disorder, where G = 1/(z − H) and G₀ = 1/(z − H₀). The average $\tilde G$ of G then satisfies $\tilde G= G_0+G_0\Sigma \tilde G$, where Σ is an infinite sum of irreducible diagrams. Here we choose a different but equivalent route, which consists in solving exactly the scattering problem for N atoms and 1 photon for a given realization of the disorder, i.e. a given draw {r_j} of the positions of the atoms, obtain the attenuation ${\cal T}$ of the sample, and subsequently average this result over several realizations of the disorder (see section 4).

For a sparse sample, we have already calculated the optical density at first order in density (equation (24)) and the result is identical for a strictly 2D gas and a thick one. The approach based on quantum scattering theory is well suited to go beyond this first-order approximation and look for differences between the 2D and 3D cases. The basis of the calculation is the series expansion of equation (34), which gives

$$\begin{equation} P\frac{1}{z-H}P = \frac{P}{z-H_0} +\sum_{n=1}^{\infty} \frac{P}{z-H_0} \left( PR(z)P \frac{1}{z-H_0}\right)^n . \end{equation} \tag{ 42 }$$

Consider the case of a resonant probe δ = 0 for simplicity. The result ${\cal T}\approx 1-\sigma _0 n^{({\rm col})}/2$ obtained for a sparse sample in equation (24) corresponds to the first term [P/(z − H₀)] of this expansion. Here we investigate the next order term and explain why one can still recover the Beer–Lambert law for a thick (3D) gas, but not for a 2D sample.

Double scattering diagrams for a thick sample (kℓ_z ≫ 1). We start our study by adding the first term (n = 1) in the expansion (42) to the zeroth-order term already taken into account in equation (24). This amounts to taking into account the diagrams where the incident photon is scattered on a single atom and those where the photon 'bounces' on two atoms before leaving the atomic sample. Inserting the first two terms of the expansion (42) into (40), we obtain

$$\begin{equation} \frac{T_{ ii}L_z}{\hbar c}=\frac{1}{2}\sigma_0 n^{({\rm col})} \left[ -i + \frac{1}{N} \sum_j \sum_{j'\neq j} {\rm e}^{{\rm i}k(z_j-z_{j'})}g_{xx}(\boldsymbol u_{jj'})\right] . \end{equation} \tag{ 43 }$$

We now have to average this result on the positions of the atoms j and j'. There are N(N − 1) ≈ N² couples (j,j'). Assuming that the gas is dilute so that the average distance between two atoms (in particular |z_j − z_j'|) is much larger than k⁻¹, the leading term in g_xx is the e^iu/u contribution of h₁(u) in equations (13) and (14). We thus arrive at

$$\begin{equation} \frac{T_{ ii}L_z}{\hbar c}=\frac{1}{2}\sigma_0 n^{({\rm col})} \left[-{\rm i}+ \frac{3N}{2k} \langle {\rm e}^{{\rm i}k(z-z')} \frac{{\rm e}^{{\rm i}k |\boldsymbol r-\boldsymbol r'|}}{|\boldsymbol r-\boldsymbol r'|} \rangle \right], \end{equation} \noindent \tag{ 44 }$$

where the average is taken over the positions r and r' of two atoms. We first calculate the average over the xy coordinates and we obtain (cf equation (20))

$$\begin{equation} \frac{T_{ ii}L_z}{\hbar c}=\frac{1}{2}\sigma_0 n^{({\rm col})} \left[-{\rm i} + \frac{{\rm i}}{2}\sigma_0 n^{({\rm col})} \langle {\rm e}^{{\rm i}k (z-z')} \,{\rm e}^{{\rm i}k |z-z'|}\rangle \right] . \end{equation} \noindent \tag{ 45 }$$

For a thick gas (kℓ_z ≫ 1) the bracket in this expression has an average value of ≈1/2. Indeed the function to be averaged is equal to 1 if z < z', which occurs in half the cases, and it oscillates and averages to zero in the other half of the cases, where z > z'. We thus obtain the approximate value of the transmission coefficient:

$$\begin{equation} k\ell_{z } \gg 1: {\cal T}=1-{\rm i}\frac{T_{ ii}L_z}{\hbar c}\approx 1-\frac{1}{2}\sigma_0 n^{({\rm col})} + \frac{1}{8} \left( \sigma_0 n^{({\rm col})} \right)^2 , \end{equation} \noindent \tag{ 46 }$$

where we recognize the first three terms of the power series expansion of ${\cal T}=\exp (-\sigma _0 n^{({\rm col})}/2)$, corresponding to the optical density ${\cal D}=\sigma _0 n^{({\rm col})}$.

Double scattering diagrams for a 2D gas (ℓ_z = 0). When all atoms are sitting in the same plane, the evaluation of the second-order term (and the subsequent ones) in the expansion of T_ii in powers of the density is modified with respect to the 3D case. The calculation starts as above and the second term in the bracket of equation (43) can now be written as

$$\begin{equation} \frac{1}{N} \sum_j \sum_{j'\neq j} g_{xx}(u_{jj'})=n^{(\rm 2D)} \int g_{xx}(\boldsymbol u)\;d^2u. \end{equation} \noindent \tag{ 47 }$$

If we keep only the terms varying as e^iu/u in h₁ and h₂ (equation (14)), we can calculate analytically the integral in (47) and find the same result as in 3D, i.e. iσ₀n^(2D)/4. If this was the only contribution to (47), it would lead to the Beer–Lambert law also in 2D, at least at second order in density. However, one can check that a significant contribution to the integral in (47) comes from the region u = kr < 1. In this region, it is not legitimate to keep only the term in e^ikr/kr in h₁, h₂, since the terms in e^ikr/(kr)³, corresponding to the short-range resonant van der Waals interaction, are actually dominant. Therefore the expansion of the transmission coefficient ${\cal T}$ in powers of the density differs from (46), and one cannot recover the Beer–Lambert law at second order in density. Calculating analytically corrections to this law could be done following the procedure of [12]. Here we will use a numerical method to determine the deviation with respect to the Beer–Lambert law (see section 4.2).

Remark 1. For a 3D gas there are also corrections to the second term in equation (45) due the 1/r³ contributions to h₁ and h₂. However these corrections have a different scaling with the density and can be made negligible. More precisely their order of magnitude is ∼n^(3D)k⁻³, to be compared with the value ∼n^(col)k⁻² of the second term in equation (45). Therefore, one can have simultaneously n^(3D)k⁻³ ≪ 1 and n^(col)k⁻² ≳ 1 if the thickness ℓ_z of the gas along z is ≫1/k.

We start our study by testing the minimal atom number that is necessary to obtain a faithful estimate of the optical density. We choose a given value of n^(col) = N/πR² and we investigate how ${\cal D}$ depends on N either for a strictly 2D gas (ℓ_z = 0) or for a gas extending significantly along the third direction (ℓ_z = 20 k⁻¹). We consider a resonant probe for this study (Δ = 0). We vary N by multiplicative steps of 2, from N = 8 up to N = 2048 and we determine how large N must be so that ${\cal D}$ is a function of n^(col) only.

The results are shown in figures 2(a) and (b), where we plot ${\cal D}$ as a function of N. We conduct this study for four values of the density n^(col), corresponding to σ₀n^(col) = 0.5,1,2 and 4. Let us consider first the smallest value σ₀n^(col) = 0.5. For each value of N we perform a number of draws that is sufficient to bring the standard error below 2 × 10⁻³ and we find that the calculated optical density is independent of N (within standard error) already for N ≳ 100, for both values of ℓ_z. Consider now our largest value σ₀n^(col) = 4; for a strictly 2D gas (ℓ_z = 0), ${\cal D}$ reaches an approximately constant value independent of N for N ≳ 1000. For σ₀n^(col) = 4 and a relatively thick gas (ℓ_z = 20 k⁻¹, blue filled circles in figure 2(b)), reaching the thermodynamic limit is more problematic since there is still a clear difference between the results obtained with 1024 and 2048 atoms. This situation thus corresponds to the limit of validity of our numerical results. In the rest of this paper, we will show only results obtained with N = 2048 atoms for column densities not exceeding σ₀n^(col) = 4. The number of independent draws of the atomic positions (at least 8) is chosen such that the standard error for each data point is below 2%.

Zoom In Zoom Out Reset image size

We now investigate the variation of the optical density ${\cal D}=\ln |{\cal T}|^{-2}$ as a function of the column density of the sample n^(col), or equivalently of the Beer–Lambert prediction ${\cal D}_{{\rm BL}}=n^{({\rm col})}\sigma$. We suppose in this section that the probe beam is resonant (Δ = 0), and we address the cases of a strictly 2D gas (ℓ_z = 0) and a thick slab (ℓ_z = 20 k⁻¹).

Consider first the case of a strictly 2D case, ℓ_z = 0, leading to the results shown in figure 3(a). We see that ${\cal D}$ differs significantly (∼25%) from ${\cal D}_{{\rm BL}}$ already for ${\cal D}_{{\rm BL}}$ around 1. A quadratic fit to the calculated variation of ${\cal D}$ for σ₀n^(2D) < 1 (continuous red line) gives

$$\begin{equation} {\cal D} \approx {\cal D}_{{\rm BL}}\left(1-0.22\, {\cal D}_{{\rm BL}} \right) . \end{equation} \tag{ 48 }$$

The discrepancy between ${\cal D}$ and ${\cal D}_{{\rm BL}}$ increases when the density increases: for ${\cal D}_{{\rm BL}}=4$, the calculated ${\cal D}$ is only ≈1.4. For such a large density the average distance between nearest neighbours is ≈k⁻¹ and the energy shifts due to the dipole–dipole interactions are comparable to or larger than the linewidth Γ. The atomic medium is then much less opaque to a resonant probe beam than in the absence of dipole–dipole coupling.

Zoom In Zoom Out Reset image size

Consider now the case of a thick sample, ℓ_z = 20 k⁻¹ (figure 3(b)). The calculated optical density is then very close to the Beer–Lambert prediction over the whole range we studied. This means that in our chosen range of optical densities, the mean-field approximation leading to ${\cal D}_{{\rm BL}}$ is satisfactory if the sample thickness exceeds a few optical wavelengths λ = 2π/k.

It is interesting to characterize how the optical density evolves from the value for a strictly 2D gas to the expected value from the Beer–Lambert law ${\cal D}_{{\rm BL}}$ as the thickness of the gas increases. We show in figure 4 the variation of ${\cal D}$ as a function of ℓ_z for three values of the column density corresponding to ${\cal D}_{{\rm BL}}=1,2$ and 4. An exponential fit ${\cal D}=\alpha + \beta \exp (-\ell _{z }/\ell _{{\rm c}})$ to these data for 2 k⁻¹ ⩽ ℓ_z ⩽ 20 k⁻¹ gives a good account of the variation observed over this range, and provides the characteristic thickness ℓ_c needed to recover the Beer–Lambert law. We find that ℓ_c ≈ 3.0 k⁻¹ for ${\cal D}_{{\rm BL}}=1$, ℓ_c ≈ 3.5 k⁻¹ for ${\cal D}_{{\rm BL}}=2$ and ℓ_c ≈ 4.4 k⁻¹ for ${\cal D}_{{\rm BL}}=4$.

Zoom In Zoom Out Reset image size

Remark 2. For the largest value of column density considered here (n^(col)σ₀ = 4), we find that ${\cal D}$ increases slightly above the value of ${\cal D}_{{\rm BL}}$ when ℓ_z is chosen larger than 20 k⁻¹ (the upper value considered in figure 4). We believe that this is a consequence of the edge terms that we neglected when approximating equation (19) by (21). These terms become significant for ${\cal D}_{\rm BL}=4$, because for our atom number N = 2048, the sample radius R ≈ 55 k⁻¹ is then not very large compared to its thickness for ℓ_z ≳ 20 k⁻¹. In order to check this assumption, we also calculated numerically the result of equation (19) (instead of equation (21)) for practical values of the parameters (position and radius) of the lens represented in figure 1(a). The results give again ${\cal D}\approx {\cal D}_{\rm BL}$, but now with ${\cal D}$ remaining below ${\cal D}_{\rm BL}$. Since our emphasis in this paper is rather placed on the 2D case, we will not explore this aspect further here.

Resonant van der Waals interactions manifest themselves not only in the reduction of the optical density at resonance, but also in the overall line shape of the absorption profile. To investigate this problem we have studied the variations of ${\cal D}$ with the detuning of the probe laser. We show in figure 5(a) the results for a strictly 2D gas (ℓ_z = 0) for n^(col)σ₀ = 1,2 and 4. Several features show up in this series of plots. First we note a blue shift of the resonance, which increases with n^(2D) and reaches Δ ≈ Γ/4 for σ₀n^(2D) = 4. We also note a slight broadening of the central part of the absorption line, since the full-width at half-maximum, which is equal to Γ for an isolated atom, is ≃1.3Γ for n^(col)σ₀ = 4. Finally, we note the emergence of large, non-symmetric wings in the absorption profile. This asymmetry is made more visible in figure 5(b), where we show with full blue squares the same data as in figure 5(a) for n^(col)σ₀ = 4, but now plotting ${\cal D}/{\cal D}_{\rm BL}$ as a function of δ. For a detuning δ = ± 15, the calculated optical density exceeds the Beer–Lambert prediction by a factor of 4.1 (resp. 2.8) on the red (resp. blue) side.

Zoom In Zoom Out Reset image size

In order to gain a better understanding of these various features, we give in figure 5(b) two additional results. On the one hand we plot with empty blue squares the variations of ${\cal D}/{\cal D}_{\rm BL}$ for a thick gas (ℓ_z = 20 k⁻¹) with the same column density n^(col)σ₀ = 4. There are still some differences between ${\cal D}$ and ${\cal D}_{\rm BL}$ in this case, as already pointed out in [18], but they are much smaller than in the ℓ_z = 0 case. This indicates that the strong deviations with respect to the Beer–Lambert law that we observe in figure 5(a) are specific 2D features. On the other hand, we plot with black stars the variations of ${\cal D}/{\cal D}_{\rm BL}$ for a 2D gas (ℓ_z = 0) in which we artificially increased the exclusion radius around each atom up to a = k⁻¹ instead of a = 0.01 k⁻¹ (blue full squares) for the other results in this paper. This procedure, which was suggested to us by Robin Kaiser, allows one to discriminate between effects due to isolated pairs of closely spaced atoms, and many-body features resulting from multiple scattering of photons among larger clusters of atoms. The comparison of the results obtained for a = 0.01 k⁻¹ and a = k⁻¹ suggests that the blue shift of the resonance line, which is present in both cases, is a many-body phenomenon, whereas the large-amplitude wings with a blue–red asymmetry, which occurs only for a = 0.01 k⁻¹, is rather an effect of close pairs.

This asymmetry in the wings of the absorption line in a 2D gas can actually be understood in a semi-quantitative manner by a simple reasoning. We recall that for two atoms at a distance r ≪ k⁻¹, the levels involving one ground and one excited atom have an energy (real part of the eigenvalues of H_eff) that is displaced by ∼ ± ℏΓ/(kr)³. A given detuning δ can thus be associated with a distance r between the two members of a pair that will resonantly absorb the light. To be more specific let us consider a pair of atoms with kr ≪ 1, and suppose for simplicity that it is aligned either along the polarization axis of the light (x) or perpendicularly to the axis (y). In both cases the excited state of the pair that is coupled to the laser is the symmetric combination $(|1:g;2:e_x\rangle +|1:e_x;2:g\rangle )/\sqrt 2$. If the pair is aligned along the x-axis, this state has an energy ℏω₀ − 3ℏΓ/2(kr)³; hence it is resonant with red detuned light such that δ = −3/(kr)³. If the pair axis is perpendicular to x, the state written above has an energy ℏω₀ + 3ℏΓ/4(kr)³; hence it is resonant with blue detuned light such that δ = 3/2(kr)³. This clearly leads to an asymmetry between red and blue detuning; indeed the pair distance r needed for ensuring resonance for a given δ > 0, r_blue = (3/2|δ|)^1/3 k⁻¹, is smaller than the value r_red = (3/|δ|)^1/3 k⁻¹ for the opposite value −δ. Since the probability density for the pair distance at small r is ${\cal P}(r)\propto r$ in 2D for randomly drawn positions, we expect the absorption signal to be stronger for −δ than for +δ. In a 3D geometry the variation of the probability density at small r is even stronger ( ${\cal P}(r)\propto r^2$), but it is compensated by the fact that the probability of occurrence of pairs that are resonant with blue detuned light is dimensionally increased. For example, in our simplified modelling where the pair axis is aligned with the reference axes, a given pair will be resonant with blue detuned light in 2/3 of the cases (axis along y or z) and resonant with red detuned light only in 1/3 of the cases (axis along x). This explains why the asymmetry of the absorption profile is much reduced for a 3D gas in comparison to the 2D case.