Absorption imaging of a quasi 2D gas: a multiple scattering analysis

Absorption imaging with quasi-resonant laser light is a commonly used technique to probe ultra-cold atomic gases in various geometries. Here we investigate some non-trivial aspects of this method when it is applied to in situ diagnosis of a quasi two-dimensional gas. Using Monte Carlo simulations we study the modification of the absorption cross-section of a photon when it undergoes multiple scattering in the gas. We determine the variations of the optical density with various parameters, such as the detuning of the light from the atomic resonance and the thickness of the gas. We compare our results to the known three-dimensional result (Beer-Lambert law) and outline the specific features of the two-dimensional case.


Introduction
The study of cold atomic gases has recently shed new light on several aspects of quantum many-body physics [1,2,3,4]. Most of the measurements in this field of research are based on the determination of the spatial density of the gas [5]. For instance one can use the in situ steady-state atomic distribution in a trapping potential to infer the equation of state of the homogenous gas [6]. Another example is the time-of-flight method, in which one measures the spatial density after switching off the trapping potential and allowing for a certain time of ballistic expansion. This gives access to the momentum distribution of the gas, and to the conversion of interaction energy into kinetic energy at the moment of the potential switch-off.
To access the atomic density n(r), one usually relies on the interaction of the atoms with quasi-resonant laser light. The most common method is absorption imaging, in which the shadow imprinted by the cloud on a low intensity probe beam is imaged on a camera. The simplest modelling of absorption imaging is based on a mean-field approach, in which one assumes that the local value of the electric field driving an atomic dipole at a given location depends only on the average density of scatterers. One can then relate the attenuation of the laser beam to the column atomic density n (col) (x, y) = n(r) dz along the line-of-sight z. The optical density of the cloud D(x, y) ≡ ln[I in (x, y)/I out (x, y)] is given by the Beer-Lambert law D BL (x, y) = σ n (col) (x, y), where σ is the photon scattering cross-section, and I in (resp. I out ) are the incoming (resp. outgoing) intensity of the probe laser in the plane xy perpendicular to the propagation axis. For a closed two-level atomic transition of frequency ω 0 = ck 0 , σ depends on the wavelength λ 0 = 2π/k 0 associated to this transition and on the detuning ∆ = ω − ω 0 between the probe light frequency ω and the atomic frequency: Here Γ represents the natural line width of the transition (i.e., Γ −1 is the natural life time of the excited state of the transition). Eq. (2) assumes that the intensity of the probe beam is much lower than the saturation intensity of the atomic transition. Quasiresonant absorption imaging is widely used to measure the spatial distribution of atomic gases after a long time-of-flight, when the density has dropped sufficiently so that the mean-field approximation leading to Eq. (1) is valid. One can also use absorption imaging to probe in situ samples, at least in the case where σ n (col) is not very large so that the output intensity is not vanishingly small. This is in particular the case for low dimensional gases. Consider for example a 2D gas, such that the translational degree of freedom along z has been frozen. For a probe beam propagating along this axis, one can transpose the Beer-Lambert law of Eq. (1) by simply replacing the column density by the surface density n (2D) of the gas. This 2D Beer-Lambert law can be heuristically justified by treating each atom as a disk of area σ that blocks every photon incident on it. In an area A σ containing N = An (2D) 1 randomly placed atoms, the probability that a photon is not blocked by any of the disks is (1 − σ/A) N ≈ exp(−σn (2D) ).
In a quasi-2D gas there is however an important limitation on the optical densities to which one may apply the Beer-Lambert prediction of Eq. (1). Already for σ 0 n (2D) = 1 the mean interparticle distance is only 0.7 λ 0 and one may expect that the optical response of an atom strongly depends on the precise location of its neighbours. More precisely the exchange of photons between closely spaced atoms induces a resonant van der Waals interaction that significantly shifts the atomic resonance frequency with respect to its bare value ω 0 . The optical density of the gas at resonance may then be reduced with respect to Eq. (1), and this was indeed observed in a series of experiments performed with a degenerate 87 Rb gas [7,8].
The general subject of the propagation of a light wave in a dense atomic sample, where multiple scattering plays an essential role, has been the subject of numerous experimental and theoretical works (see e.g. [9,10] in the context of cold atoms, and [11] for a review). Here we present a quantitative treatment of the collective effects that appear when a weak probe beam interacts with a quasi-2D atomic gas. We consider an ensemble of N atoms at rest with random positions and we investigate the transmission of quasi-resonant light by the atom sheet. We model the resonance transition between the atomic ground (g) and excited (e) states by a J g = 0 ↔ J e = 1 transition. We present two equivalent approaches; the first one is based on the calculation of the field radiated by an assembly of N dipoles, where each dipole is driven by an external field plus the field radiated by the N −1 other dipoles; the second one uses the standard T matrix formalism of scattering theory. We show that in both cases the optical density of the medium can be determined by solving the same 3N × 3N linear system. A similar formalism has been previously used for the study of light propagation in small 3D atomic samples, in the presence of multiple scattering (see e.g. [12,13,14,15,16,17,18,19]). However its application to quasi-2D samples has (to our knowledge) not yet been investigated, except in the context of Anderson localisation of light [13]. Our numerical calculations are performed for N = 2048 atoms, which is sufficient to reach the 'thermodynamic limit' for the range of parameters that is relevant for experiments. We show in particular that even for moderate values of σ 0 n (col) , the optical density is notably reduced compared to what is expected from the Beer-Lambert law (e.g., more than 20 % reduction for σ 0 n (col) = 1). We investigate how the absorption line shape is modified by the resonant van der Waals interactions and we also show how the result (1) is recovered when one increases the thickness of the gas, for a given column density n (col) .
The paper is organised as follows. In section 2, we detail the modelling of the atomlight interaction with the two-level and rotating wave approximations. Then we explain the principle of the calculation for the absorption of a weak probe beam crossing the atom slab (section 3). The ensemble of our numerical results are presented in section 4. Finally in section 5 we discuss some limitations to our model and draw some concluding remarks.

The electromagnetic field
We use the standard description of the quantised electromagnetic field in the Coulomb gauge [20], and choose periodic boundary conditions in the cubic-shaped quantisation volume V = L x L y L z . We denote a q,s the destruction operator of a photon with wave vector q and polarisation s (s⊥q). The Hamiltonian of the quantised field is and the transverse electric field operator reads E(r) = E (+) (r) + E (−) (r) with and The wave vectors q are quantised in the volume V as where n i is a positive or negative integer.

The atomic medium
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 ω 0 the energy difference between e and g. The atomic Hamiltonian is thus (up to a constant) The restriction to a two-level approximation is legitimate if the detuning ∆ between the probe and the atomic frequencies is much smaller than ω 0 . The modelling of this transition by a J g = 0 ↔ 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 ↔ J e = 3 for Rb atoms in [7,8]), but are more complex to handle [21,22] and they are thus out of 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.

The atom-light coupling
We treat the atom-light interaction using the electric dipole approximation (length gauge), which is legitimate since the resonance wavelength of the atoms λ 0 is much larger than the atomic size. We write the atom-light coupling as: where D j is the dipole operator for the atom j. We will use the rotating wave approximation (RWA), which consists in keeping only the resonant terms in the coupling: where h.c. stands for Hermitian conjugate. Here D (+) j represents the raising part of the dipole operator for atom j: where d is the electric dipole associated to the g − e transition andû α 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 (Lamb shift) and in the fact that the excited state e acquires a non-zero width Γ For simplicity we will incorporate the Lamb shift in the definition of ω 0 . Note that the proper calculation for this shift requires that one goes beyond the two-level and the rotating wave approximations. 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 atomlight 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 were also taken into account. The small parameter associated to the RWA is ∆/ω 0 , which is in practice in the range 10 −6 − 10 −9 ; the RWA is thus an excellent approximation.
Formally the use of the electric dipole interaction implies to add to the Hamiltonian an additional contact term between the dipoles (see e.g. [23,12]). 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 Sec. 4) do not depend on the size of the excluded volume, and we can safely omit the additional contact term in the present work.

Interaction of a probe laser beam with a dense quasi-2D atomic sample
We present in this section the general formalism that allows one to calculate the absorption of a quasi-resonant laser beam by a slab of N atoms. We address this question using two different approaches. The first one maps the problem onto the collective behaviour of an assembly of N oscillating dipoles [12]. The equation of motion for each dipole is obtained using the Heisenberg picture for the Hamiltonian presented in section 2. It contains two driving terms, one from the incident probe field and one from the field radiated by all the other dipoles at the location of the dipole under study. The steady-state of this assembly of dipoles is obtained by solving a set of 3N linear equations. The second approach uses the standard quantum scattering theory [24], which is well suited for perturbative calculations and partial resummations of diagrams. We suppose that one photon is incident on the atomic medium and we use resummation techniques to take into account the multiple scattering events that can occur before the photon emerges from the medium. The relevant quantity in this approach is the probability amplitude T ii that the outgoing photon is detected in the same mode as the incident one [14,17], and we show that T ii is obtained from the same set of equations as the values of the dipoles in the first approach.

Wave propagation in an assembly of driven dipoles.
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 ( =û x ). Since we consider a J g = 0 ↔ 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 polarisation. 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 polarisation 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: This equation can be integrated between the initial time t 0 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: 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 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 E(r, t) =Ē(r) e −iωt + c.c., and we obtain after some algebra (see e.g. [12,25]) where we set and The function g α,α (kr) is identical to the one appearing in classical electrodynamics [26], 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] ( where the reduced detunig δ = 2∆/Γ has been defined in Eq. (2) and u j,j = k(r j − r j ). This can be written with matrix notation where the 3N vectors |X and |Y are defined by 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 [27], 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 it will be the case in this work, see Sec. 4) have been studied in [27,28,29,30]. Eq. (15) has a simple physical interpretation: in 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 [31] who named it, together with Eq. (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.

Absorption signal
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 = − /2 and z = /2. We can consider two experimental setups to address this problem. The first one, represented in Fig. 1a, consists in measuring after the atomic sample the total light intensity with the same momentum k = kû 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 = > /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 Fig. 1b, 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 favored 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. We assume that the lens in Fig. 1a 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: where the integral runs over the lens area. Since the incident probe beam is supposed to be linearly polarised along x, we calculate the x component of the field in F . Plugging the value of the field given in Eqs. (12,17) we obtain the transmission coefficient 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 1 . More precisely the domain in the lens plane contributing to the integral for the dipole j is essentially a disk of radius . When this small disk is entirely included in the lens aperture, i.e., the larger disk of radius R centered on x = y = 0, we obtain We use the result (20) for all atoms, which amounts to neglect edge effects for the dipoles located at the border of the lens, and we obtain: with n (col) = N/πR 2 and where the coefficient Π is defined by 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 T is known the optical density is obtained from As an example, consider the limit of a very sparse sample where multiple scattering does not play a significant role (σ 0 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:

Light absorption as a quantum scattering process
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 and we consider the initial state where all atoms are in their ground state and where a single photon of wave vector k = kû z and polarisation =û x is incident on the atomic medium with |G ≡ |1 : g, 2 : g, . . . , N : g . The state |Ψ i is an eigenstate of H 0 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 |G ⊗ |q, s 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 |E j,α = |1 : g, . . . , j : e α , . . . , N : g ⊗ |vac , where |vac stands for the vacuum state of the electromagnetic field. The subspace spanned by the states |E j,α 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 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 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 0 , 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 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 We define the displacement operator and use the general result [20] where the effective Hamiltonian H eff is For the following calculations, it is convenient to introduce the dimensionless matrix [M ] proportional to the denominator of the right hand side of (34): 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 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 ω 0 , and its imaginary part reads: 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 already introduced for the field radiated in r j by a dipole located in r j : where the tensor g α,α is defined in Eqs. (13)(14). Suppose now that the atoms are uniformly distributed over the transverse area L x L y of the quantisation volume. We set n (col) = N/(L x L y ) and we rewrite the expression (30) of the desired matrix element T ii as where the coefficient Π has been defined in (22). The result (40) combined with (21) leads to which constitutes the 'optical theorem' for our slab geometry, since it relates the attenuation of the probe beam 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 where 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 non-Hermitian [32,33]. For a short distance r between the atoms (kr 1), the leading term in h 1 (u) and h 2 (u) is u −3 and the energies (real parts of the eigenvalues) of the six eigenstates vary as ∼ ± Γ/(kr) 3 (resonant dipoledipole 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 [34].
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 [29,30]. More precisely, one could think of solving the linear system (16), or equivalently calculating T ii in Eq. (30), by using the expansion of the column vector |Y defined in Eq. (17) on the left (|α j ) and right ( β j |) eigenvectors of H eff . Then one could inject this expansion in 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.

Beyond the sparse sample case: 3D vs. 2D
For a sparse sample, we already calculated the optical density at first order in density (Eq. (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 Eq. (34), which gives Consider the case of a resonant probe δ = 0 for simplicity. The result T ≈ 1 − σ 0 n (col) /2 obtained for a sparse sample in Eq. (24) corresponds to the first term [P/(z − H 0 )] 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 1). We start our study by adding the first term (n = 1) in the expansion (42) to the zero-th order term already taken into account in Eq. (24). This amounts to take 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. Injecting the first two terms of the expansion (42) into (40), we obtain We now have to average this result on the positions of the atoms j and j . There are N (N − 1) ≈ N 2 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 −1 , the leading term in g xx is the e iu /u contribution of h 1 (u) in Eqs. (13)- (14). We thus arrive at 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 get (cf. Eq. (20)) For a thick gas (k 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 of 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: where we recognize the first three terms of the power series expansion of T = exp(−σ 0 n (col) /2), corresponding to the optical density D = σ 0 n (col) .
Double scattering diagrams for a 2D gas ( = 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 Eq. (43) can now be written If we keep only the terms varying as e iu /u in h 1 and h 2 (Eq. (14)), we can calculate analytically the integral in (47) and find the same result as in 3D, i.e., iσ 0 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 1 , h 2 , since the terms in e ikr /(kr) 3 , corresponding to the short range resonant van der Waals interaction, are actually dominant. Therefore the expansion of the transmission coefficient 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. For a 3D gas there are also corrections to the second term in Eq. (45) due the 1/r 3 contributions to h 1 and h 2 . 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 −3 , to be compared with the value ∼ n (col) k −2 of the second term in Eq. (45).
Therefore one can have simultaneously n (3D) k −3 1 and n (col) k −2 1, if the thickness of the gas along z is 1/k.

Absorption of light by a slab of atoms
In order to study quantitatively the optical response of a quasi-2D gas, we have performed a Monte Carlo calculation of the transmission factor T given in Eq. (21), and of the related optical density D = ln |T | −2 . We start our calculation by randomly drawing the positions of the N atoms, we then solve numerically the 3N × 3N linear system (16), and finally inject the result for the N dipoles in the expression of T . The atoms are uniformly distributed in a cylinder of axis z, with a radius R and a thickness . The largest spatial densities considered in this work correspond to a mean inter-particle distance ≈ k −1 . Around each atom we choose a small excluded volume with a linear size a = 0.01 k −1 . We varied a by a factor 10 around this value and checked that our results were essentially unchanged. Apart from this excluded volume we do not include any correlation between the positions of the atoms. This choice is justified physically by the fact that, in the case of large phase space densities which motivates our study, the density fluctuations in a 2D Bose gas are strongly reduced and the two-body correlation function g 2 (r, r ) is such that g 2 (r, r) ≈ 1 [35].
In this section we first determine the value of N that is needed to reach the 'thermodynamic limit' for our problem: for a given thickness , D should not be an independent function of the number of atoms N and the disk radius R, but should depend only of the ratio N/πR 2 = n (col) . We will see that this imposes to use relatively large number of atoms, typically N > 1000, for the largest spatial densities considered here. All subsequent calculations are performed with N = 2048. We then study the dependence of D with the various parameters of the problem: the column density n (col) , the thickness of the gas , and the detuning ∆. In particular we show that for a given n (col) we recover the 3D result (1) when the thickness is chosen sufficiently large.

Reaching the 'thermodynamic limit'
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 2 and we investigate how D depends on N either for a strictly 2D gas ( = 0) or for a gas extending significantly along the third direction ( = 20 k −1 ). 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 D is a function of n (col) only.
The results are shown in Fig. 2a and Fig. 2b, where we plot D as a function of N . We perform this study for four values of the density n (col) , corresponding to σ 0 n (col) = 0.5, 1, 2 and 4. Let us consider first the smallest value σ 0 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 −3 and we find that the calculated optical density is independent of N (within standard error) already for N 100, for both values of . Consider now our largest value σ 0 n (col) = 4; for a strictly 2D gas ( = 0), D reaches an approximately constant value independent of N for N 1000. For σ 0 n (col) = 4 and a relatively thick gas ( = 20 k −1 , blue squares in Fig. 2b), 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 remaining part of the paper we will show only results obtained with N = 2048 atoms for column densities not exceeding σ 0 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%.

Measured optical density vs. Beer-Lambert prediction
We now investigate the variation of the optical density D = ln |T | −2 as function of the column density of the sample n (col) , or equivalently of the Beer-Lambert prediction D BL = n (col) σ. We suppose in this section that the probe beam is resonant (∆ = 0), and we address the cases of a strictly 2D gas ( = 0) and a thick slab ( = 20 k −1 ).
Consider first the case of a strictly 2D case, = 0, leading to the results shown in Fig. 3a. We see that D differs significantly (∼ 25%) from D BL already for D BL around 1. A quadratic fit to the calculated variation of D for σ 0 n (2D) < 1 (continuous red line) The discrepancy between D and D BL increases when the density increases: for D BL = 4, the calculated D is only ≈ 1.4. For such a large density the average distance between nearest neighbours is ≈ k −1 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. Consider now the case of a thick sample, = 20 k −1 (Fig. 3b). The calculated optical density is then very close to the Beer-Lambert prediction over the whole range that we studied. This means that in our chosen range of optical densities, the mean-field approximation leading to D BL is satisfactory as soon as 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 D BL when the thickness of the gas increases. We show in Fig. 4 the variation of D as function of for three values of the column density corresponding to D BL = 1, 2 and 4. An exponential fit D = α + β exp(− / c ) to these data for 2 k −1 ≤ ≤ 20 k −1 gives a good account of the observed variation over this range, and it provides the characteristic thickness c needed to recover the Beer-Lambert law. We find that c ≈ 3.0 k −1 for D BL = 1, c ≈ 3.5 k −1 for D BL = 2, and c ≈ 4.4 k −1 for D BL = 4. Remark. For the largest value of the column density considered here (n (col) σ 0 = 4) we find that D increases slightly above the value D BL when is chosen larger than 20 k −1 (upper value considered in Fig. 4). We believe that this is a consequence of the edge terms that we neglected when approximating Eq. (19) by Eq. (21). These terms become significant for D BL = 4 because for our atom number N = 2048, the sample radius R ≈ 55 k −1 is then not very large compared to its thickness for 20 k −1 . In order to check this assumption, we also calculated numerically the result of Eq. (19) (instead of Eq. (21)) for practical values of the parameters (position and radius) of the lens represented in Fig. 1a. The results give again D ≈ D BL , but now with D remaining below D BL . Since our emphasis in this paper is rather put on the 2D case, we will not explore this aspect further here.

Absorption line shape
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 D with the detuning of the probe laser. We show in Fig. 5a the results for a strictly 2D gas ( = 0) for n (col) σ 0 = 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 σ 0 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) σ 0 = 4. Finally we note the emergence of large, non-symmetric wings in the absorption profile. This asymmetry is made more visible in Fig. 5b, where we show with full blue squares the same data as in Fig. 5a for n (col) σ 0 = 4, but now plotting D/D BL as function of δ. For a detuning δ = ±15, the calculated optical density exceeds the Beer-Lambert prediction by a factor 4.1 (resp. 2.8) on the red (resp. blue) side.
In order to get a better understanding of these various features, we give in Fig. 5b two additional results. On the one hand we plot with empty blue squares the variations of D/D BL for a thick gas ( = 20 k −1 ) with the same column density n (col) σ 0 = 4. There are still some differences between D and D BL in this case, as already pointed out in [17], but they are much smaller than in the = 0 case. This indicates that the strong deviations with respect to the Beer-Lambert law that we observe in Fig. 5a are specific 2D features. On the other hand we plot with black stars the variations of D/D BL for a 2D gas ( = 0) in which we artificially increased the exclusion radius around each atom up to a = k −1 instead of a = 0.01 k −1 (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 manybody features resulting from multiple scattering of photons among larger clusters of atoms. The comparison of the results obtained for a = 0.01 k −1 and a = k −1 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 −1 , 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 −1 , 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) 3 . A given detuning δ can thus be associated to 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 + |1 : e x ; 2 : g )/ √ 2. If the pair is aligned along the x axis, this state has an energy ω 0 − 3 Γ/2(kr) 3 , hence it is resonant with red detuned light such that δ = −3/(kr) 3 . If the pair axis is perpendicular to x, the state written above has an energy ω 0 + 3 Γ/4(kr) 3 , hence it is resonant with blue detuned light such that δ = 3/2(kr) 3 . 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 , is smaller than the value r red = (3/|δ|) 1/3 for the opposite value −δ. Since the probability density for the pair distance is P(r) ∝ 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 with r is even stronger (P(r) ∝ 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 references 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.

Summary
We have presented in this paper a detailed analysis of the scattering of light by a disordered distribution of atoms in a quasi-two dimensional geometry. The particles were treated as fixed scatterers and their internal structure was modeled as a two-level system, with a J = 0 ground state and a J = 1 excited state. In spite of these simplifying assumptions the general trend of our results is in good agreement with the experimental finding of [7], where a variation of the measured optical density similar to that of Fig. 3 was measured.
Several improvements in our modeling can be considered in order to reach a quantitative agreement with theory and experiment. The first one is to include the relatively complex atomic structure of the alkali-metal species used in practice, with a multiply degenerate ground state; this could be done following the lines of [21,22]. A second improvement consists in taking into account the atomic motion. This is in principle a formidable task, because it leads to a spectacular increase in the dimension of the relevant Hilbert space. This addition can however be performed in practice in some limiting cases, for example if one assumes that the particles are tightly bound in a lattice [36,37]. When the atom-light interaction is used only to probe the spatial atomic distribution of the gas, neglecting the particle motion should not be a major problem. Indeed the duration of the light pulse is quite short (∼ 10 microseconds only). Each atom scatters only a few photons in this time interval and its displacement is then smaller than the mean interatomic spacing for the spatial densities encountered in practice. The acceleration of the atoms under the effect of resonant van der Waals interaction should also have a minor effect under relevant experimental conditions. Finally another aspect that could be valuably studied is the interaction of the gas with an intense laser beam [38]. One could thus validate the intuitive idea that saturation phenomena reduce the effects of resonant van der Waals interactions [39,8], and are thus helpful to provide a faithful estimate of the atomic density from the light absorption signal.