text:
# Nonlinear magneto-optical rotation in optically thick media

## I Introduction

Resonant nonlinear magneto-optical rotation (NMOR) Gawlik (1994); Budker et al. (1999) has been the subject of extensive theoretical and experimental studies because it provides a very sensitive way to measure magnetic fields (see, e.g., Ref. Budker et al. (2000a)). While NMOR experiments are usually carried out with vapor samples of moderate optical thickness (no more than $`2`$ absorption lengths), the intense research of recent years on electromagnetically-induced transparency (EIT; see, e.g., Ref. Harris (1997) for a review) has motivated investigation of optically thick media Scully and Fleischhauer (1992); Fleischhauer and Scully (1994).

Recently, NMOR in the vicinity of the D1 and D2 lines by an optically thick vapor of rubidium was studied Sautenkov et al. (2000); Novikova et al. (2001); Matsko et al. (2001); Novikova and Welch (2002); Matsko et al. (2002). The authors of Ref. Novikova et al. (2001) used a $`5`$cm-long vapor cell containing 87Rb and a laser beam tuned to the maximum of the NMOR spectrum near the D1 line (laser power was $`2.5`$ mW, beam diameter $`2`$ mm). They measured maximum (with respect to laser frequency and magnetic field) polarization rotation $`\phi _{max}`$ as a function of atomic density. It was found that $`\phi _{max}`$ increases essentially linearly up to $`n3.5\times 10^{12}\mathrm{cm}^3`$. At higher densities, $`d\phi _{max}/dn`$ decreases and eventually becomes negative. The maximum observed rotation was $`10`$rad with the applied magnetic field of $`0.6`$G. (This shows the effect of power broadening on the magnetic field dependence; at low light power, the maximum rotation would occur at a magnetic field about an order of magnitude smaller.) Rotation slope $`d\phi /dB|_{B=0}`$, an important parameter for magnetometry, was also measured as a function of atomic density in Ref. Novikova et al. (2001).

The goal of the present contribution is to analyze the scaling of optimal magnetometric sensitivity with respect to optical thickness of the sample. We analyze NMOR for a simple system—an isolated $`J=1J=0`$ transition—for which analytical solutions for the density matrix are readily obtained. We assume that the transverse thickness of the sample is small, so that the trapping of spontaneously-emitted radiation can be neglected. (This assumption may not be justified in actual experiments with optically thick media Matsko et al. (2001).) We also neglect mixing between different velocity groups due to velocity-changing collisions. Such mixing is important in the presence of a buffer gas, or for anti-relaxation-coated cells.

A theoretical treatment of NMOR in dense media has been carried out earlier in Refs. Sautenkov et al. (2000) and Fleischhauer et al. (2000). The present analysis partially overlaps with this earlier work and extends it in several ways: we extend the treatment to the case of large Doppler broadening and the intermediate regime where power and Doppler widths are comparable; provide a qualitative discussion of the contribution to optical rotation as a function of optical depth in the medium; discuss the role of the coherence and the alignment-to-orientation conversion effects Budker et al. (2000b); discuss scaling of magnetometric sensitivity with optical thickness of the medium; and finally, give a general argument on scaling based on the observation that an optimized nonlinear optical rotation measurement is a way to carry out a spin measurement of an atomic ensemble with sensitivity given by fundamental quantum noise limits Budker et al. (2000a); Fleischhauer et al. (2000).

## II Description of the density matrix calculation

The calculation is based on a standard density matrix approach. The time evolution of the atomic density matrix $`\rho `$ is given by the Liouville equation (see, e.g., Ref. Stenholm (1984)):

$`{\displaystyle \frac{d\rho }{dt}}={\displaystyle \frac{1}{i\mathrm{}}}[H,\rho ]{\displaystyle \frac{1}{2}}\{\mathrm{\Gamma },\rho \}+\mathrm{\Lambda },`$ (1)
where the square brackets denote the commutator and the curly brackets the anti-commutator, and the total Hamiltonian $`H`$ is the sum of light-atom interaction Hamiltonian $`H_L=\stackrel{}{d}\stackrel{}{E}`$ (where $`\stackrel{}{E}`$ is the electric field vector, and $`\stackrel{}{d}`$ is the electric dipole operator), the magnetic field-atom interaction Hamiltonian $`H_B=\stackrel{}{\mu }\stackrel{}{B}`$ (where $`\stackrel{}{B}`$ is the magnetic field and $`\stackrel{}{\mu }`$ is the magnetic moment), and the unperturbed Hamiltonian $`H_0`$. $`\mathrm{\Gamma }`$ is the relaxation matrix (diagonal in the collision-free approximation)

$`\xi Jm\left|\mathrm{\Gamma }\right|\xi Jm=\gamma +\gamma _0\delta (\xi ,\xi _e),`$ (2)
where $`\gamma `$ is the ground state depolarization rate (e.g., due to transit of atoms through the laser beam), $`\gamma _0`$ is the spontaneous decay rate from the upper state, and $`\xi `$ represents the quantum number distinguishing the ground state ($`\xi _g`$) from the excited state ($`\xi _e`$). $`\mathrm{\Lambda }=\mathrm{\Lambda }^0+\mathrm{\Lambda }^{repop}`$ is the pumping term, where the diagonal matrix

$`\xi _gJ_gm\left|\mathrm{\Lambda }^0\right|\xi _gJ_gm={\displaystyle \frac{\gamma \rho _0}{\left(2J_g+1\right)}}`$ (3)
describes incoherent ground state pumping ($`\rho _0`$ is the atomic density), and

$`\xi _gJ_gm\left|\mathrm{\Lambda }^{repop}\right|\xi _gJ_gm^{}=\gamma _0{\displaystyle \underset{m_e,m_e^{},q}{}}J_g,m,1,q|J_e,m_eJ_g,m^{},1,q|J_e,m_e^{}\rho _{\xi _eJ_em_e\xi _eJ_em_e^{}},`$ (4)
describes repopulation due to spontaneous relaxation from the upper level (see, e.g., Ref. Rautian and Shalagin (1991)). Here $`\mathrm{}|\mathrm{}`$ are the Clebsch-Gordan coefficients.

The electric field vector is written (see, e.g., Ref. Huard (1997))

$`\stackrel{}{E}`$ $`={\displaystyle \frac{1}{2}}[E_0e^{i\varphi }(\mathrm{cos}\phi \mathrm{cos}ϵi\mathrm{sin}\phi \mathrm{sin}ϵ)e^{i\left(\omega tkz\right)}+c.c.]\widehat{x}`$
$`+{\displaystyle \frac{1}{2}}[E_0e^{i\varphi }(\mathrm{sin}\phi \mathrm{cos}ϵ+i\mathrm{cos}\phi \mathrm{sin}ϵ)e^{i\left(\omega tkz\right)}+c.c.]\widehat{y},`$
where $`\omega `$ is the light frequency, $`k=\omega /c`$ is the vacuum wave number, $`E_0`$ is the electric field amplitude, $`\phi `$ is the polarization angle, $`ϵ`$ is the ellipticity (arctangent of the ratio of the major and minor axes of the polarization ellipse), and $`\varphi `$ is the overall phase. By substituting (LABEL:lightfield) into the wave equation

$`\left({\displaystyle \frac{\omega ^2}{c^2}}+{\displaystyle \frac{d^2}{dz^2}}\right)\stackrel{}{E}={\displaystyle \frac{4\pi }{c^2}}{\displaystyle \frac{d^2}{dt^2}}\stackrel{}{P},`$ (6)
where $`\stackrel{}{P}=Tr(\rho \stackrel{}{d})`$ is the polarization of the medium, the absorption, rotation, phase shift, and change of ellipticity per unit distance for an optically thin medium can be found in terms of the density matrix elements (these expressions are given in Ref. Rochester et al. (2001)). Once the solutions for the density matrix are obtained, we will perform an integration to generalize the result to media of arbitrary thickness.

## III The Doppler-free case

We consider the case of a $`J=1J=0`$ transition, and linearly-polarized incident light, with a magnetic field directed along the light propagation direction (Faraday geometry). Using the rotating wave approximation, the solution of Eq. (1) is obtained, and from this, analytic expressions for thin medium absorption and rotation are found. These expressions can be simplified by assuming that $`\gamma \gamma _0`$. We first consider the case where the power-broadened line width is much greater than the Doppler width. In this case, the absorption coefficient per unit length $`\alpha `$ is found to be

$`\alpha {\displaystyle \frac{\alpha _0}{(2\mathrm{\Delta }/\gamma _0)^2+2\kappa /3+1}},`$ (7)
where $`\mathrm{\Delta }`$ is the light-frequency detuning from resonance,

$`\kappa ={\displaystyle \frac{d^2E_0^2}{\mathrm{}^2\gamma \gamma _0}}`$ (8)
is the optical pumping saturation parameter, and

$`\alpha _0{\displaystyle \frac{1}{6\pi }}\lambda ^2n`$ (9)
is the unsaturated absorption coefficient on resonance, where $`\lambda `$ is the transition wavelength and $`n`$ is the atomic density. For $`\mathrm{\Omega }\gamma `$, where $`\mathrm{\Omega }=g\mu B`$ is the Larmor frequency ($`g`$ is the Landé factor, and $`\mu `$ is the Bohr magneton), the slope of optical rotation per unit length, $`d\phi /\left(d\mathrm{\Omega }dx\right)`$, (proportional to rotation for small magnetic fields) is found to be

$`{\displaystyle \frac{d\phi }{d\mathrm{\Omega }dx}}{\displaystyle \frac{1}{\gamma }}{\displaystyle \frac{\alpha _0}{(2\mathrm{\Delta }/\gamma _0)^2+2\kappa /3}},`$ (10)
where we have neglected linear optical rotation. In general, optical rotation can be induced by either linear dichroism or circular birefringence. Analysis of the steady-state polarization of the ground state shows that both processes contribute here: the contribution due to linear dichroism induced in the medium is given by

$`{\displaystyle \frac{d\phi }{d\mathrm{\Omega }dx}}|_{dichr.}{\displaystyle \frac{1}{\gamma }}{\displaystyle \frac{1}{1+\left(2\mathrm{\Delta }/\gamma _0\right)^2}}{\displaystyle \frac{\alpha _0}{(2\mathrm{\Delta }/\gamma _0)^2+2\kappa /3}},`$ (11)
and the contribution due to circular birefringence (arising due to alignment-to-orientation conversion in the presence of both the magnetic field and the strong electric field of the light Budker et al. (2000b)) is given by

$`{\displaystyle \frac{d\phi }{d\mathrm{\Omega }dx}}|_{biref.}{\displaystyle \frac{1}{\gamma }}{\displaystyle \frac{\left(2\mathrm{\Delta }/\gamma _0\right)^2}{1+\left(2\mathrm{\Delta }/\gamma _0\right)^2}}{\displaystyle \frac{\alpha _0}{(2\mathrm{\Delta }/\gamma _0)^2+2\kappa /3}}.`$ (12)
The sum of the two contributions produces the Lorentzian line shape of Eq. (10).

We now generalize the formulas for absorption and rotation to the case of thick media, and find the magnetometric sensitivity. In the Doppler-free case we are presently considering, we can further simplify the expressions by assuming $`\mathrm{\Delta }=0`$ (we will have to include non-zero detunings in the discussion of the Doppler-broadened case below), and $`\kappa 1`$. (We will see from the final result that this holds everywhere in the medium when the input light power is optimal.) Generally, the medium in the presence of a magnetic field produces light ellipticity as well as rotation; however, the ellipticity is an odd function of detuning and is zero on resonance. The rate of change of the saturation parameter as light travels through the medium is given by

$`{\displaystyle \frac{d\kappa (x)}{dx}}`$ $`=\alpha \kappa (x)`$
$`{\displaystyle \frac{3}{2}}\alpha _0,`$ (13)
so solving for the saturation parameter as a function of position,

$`\kappa (x)\kappa _0{\displaystyle \frac{3}{2}}\alpha _0x,`$ (14)
where $`\kappa _0`$ is the saturation parameter at $`x=0`$.

The contribution to the small-field optical rotation of the “slice” of the medium at position $`x`$ is found by substituting Eq. (14) into Eq. (10):

$`{\displaystyle \frac{d\phi }{d\mathrm{\Omega }dx}}(x)`$ $`{\displaystyle \frac{1}{\gamma }}{\displaystyle \frac{\alpha _0}{2\kappa _0/3\alpha _0x}}.`$ (15)
This contribution is plotted as a function of position in Fig. 1. The plot illustrates that the part of the medium near its end contributes significantly more to the overall rotation than the part near its beginning. This is because light power, and correspondingly, the light broadening of the resonance, is lower at the end of the medium \[the ratio of the rotation per unit length is approximately $`\kappa (\mathrm{})/\kappa _01(3/2)\alpha _0\mathrm{}/\kappa _0`$\].

Integrating over the length $`\mathrm{}`$ of the medium gives the total slope:

$`{\displaystyle \frac{d\phi }{d\mathrm{\Omega }}}`$ $`={\displaystyle _0^{\mathrm{}}}{\displaystyle \frac{d\phi }{d\mathrm{\Omega }dx}}(x)𝑑x`$
$`{\displaystyle \frac{1}{\gamma }}{\displaystyle _0^{\mathrm{}}}{\displaystyle \frac{\alpha _0}{2\kappa _0/3\alpha _0x}}𝑑x`$
$`={\displaystyle \frac{1}{\gamma }}\mathrm{ln}\left({\displaystyle \frac{\kappa _0}{\kappa _03\alpha _0\mathrm{}/2}}\right).`$ (16)
The slope (16) is plotted as a function of $`\alpha _0\mathrm{}`$ in Fig. 2.

The photon shot-noise-limited magnetometric sensitivity $`\delta B`$ 11endnote: 1In the general case, there may exist an additional source of noise due to AC-Stark shift associated with off-resonant levels Fleischhauer et al. (2000); however, this source of noise is absent for an isolated transition such as the one considered here. is given in terms of the number $`N_\gamma `$ of transmitted photons per unit time, the slope of rotation with respect to $`B`$, and the measurement time $`t`$:

$`(\delta B)^1`$ $`=2{\displaystyle \frac{d\phi }{dB}}\sqrt{N_\gamma t}`$
$`={\displaystyle \frac{g\mu \omega }{\mathrm{}c}}{\displaystyle \frac{d\phi }{d\mathrm{\Omega }}}\sqrt{{\displaystyle \frac{2At}{\pi }}\gamma \kappa (\mathrm{})}`$
$`{\displaystyle \frac{g\mu \omega }{\mathrm{}c}}\sqrt{{\displaystyle \frac{2At}{\pi }}{\displaystyle \frac{\kappa _0\frac{3}{2}\alpha _0\mathrm{}}{\gamma }}}\mathrm{ln}\left({\displaystyle \frac{\kappa _0}{\kappa _03\alpha _0\mathrm{}/2}}\right)`$
$`={\displaystyle \frac{g\mu \omega }{\mathrm{}c}}\sqrt{{\displaystyle \frac{3}{\pi }}{\displaystyle \frac{At\alpha _0\mathrm{}}{\gamma }}}\sqrt{q1}\mathrm{ln}\left({\displaystyle \frac{q}{q1}}\right),`$ (17)
where $`A`$ is the cross-sectional area of the light beam and we have made the change of variables $`\kappa _0=3q\alpha _0\mathrm{}/2`$. The factor

$`\sqrt{q1}\mathrm{ln}\left({\displaystyle \frac{q}{q1}}\right)`$
reaches a maximum of $`0.8`$ at $`q1.2`$ so we see that for media of sufficient thickness, i.e. where $`\alpha _0\mathrm{}3`$, and for the optimum initial saturation parameter $`\kappa _01.8\alpha _0\mathrm{}`$, we have from Eq. (14) that $`\kappa (x)1`$ for all $`x`$. For the optimum light intensity, then,

$`(\delta B)^1`$ $`0.8{\displaystyle \frac{g\mu \omega }{\mathrm{}c}}\sqrt{{\displaystyle \frac{At\alpha _0\mathrm{}}{\gamma }}}`$
$`1.1{\displaystyle \frac{g\mu }{\mathrm{}}}\sqrt{{\displaystyle \frac{A\mathrm{}nt}{\gamma }}}.`$ (18)
This result is consistent with the general observation Budker et al. (2000a) that optimized NMOR provides a method for measuring a spin-system at the standard quantum limit (SQL) given by spin-projection noise 22endnote: 2The use of so-called spin-squeezed quantum states of light (see, e.g., Ref. Ulam-Orgikh and Kitagawa (2001) and references therein) can, in principle, allow one to overcome the SQL. We consider only non-squeezed states of atoms and light here.. The sensitivity is then expected to scale as the square root of the product of the number of available atoms and the spin-relaxation time (see, e.g., Ref. Khriplovich and Lamoreaux (1997), Sec. 3.1.1), which is indeed the result (18).

## IV The Doppler-broadened case

Now we consider the case where the Doppler width $`\mathrm{\Gamma }`$ is much greater than the power-broadened line width. In this case, we need to average over the over the atomic velocity distribution, which is equivalent to averaging over the Doppler-free spectral profiles. On resonance with the Doppler broadened transition, the absorption coefficient is given by

$`\alpha _{DB}`$ $`{\displaystyle \frac{1}{\sqrt{\pi }\mathrm{\Gamma }}}{\displaystyle _{\mathrm{}}^{\mathrm{}}}\alpha (\mathrm{\Delta })𝑑\mathrm{\Delta }`$
$`{\displaystyle \frac{\alpha _0}{\sqrt{\pi }\mathrm{\Gamma }}}{\displaystyle _{\mathrm{}}^{\mathrm{}}}{\displaystyle \frac{d\mathrm{\Delta }}{\left(2\mathrm{\Delta }/\gamma _0\right)^2+2\kappa /3}}`$
$`=\sqrt{{\displaystyle \frac{3\pi }{8}}}{\displaystyle \frac{\gamma _0}{\mathrm{\Gamma }}}{\displaystyle \frac{\alpha _0}{\sqrt{\kappa }}}.`$ (19)
In this case the Doppler-broadened unsaturated absorbtion coefficient is given in terms of the Doppler-free unsaturated absorbtion coefficient by

$`\alpha _0|_{DB}={\displaystyle \frac{\sqrt{\pi }}{2}}{\displaystyle \frac{\gamma _0}{\mathrm{\Gamma }}}\alpha _0.`$ (20)
Comparing Eqs. (19) and (7), we see that we have reproduced a well-known result (see, e.g., Ref. Demtröder (1996), Sec. 7.2.1) that resonant absorption falls as $`1/\kappa `$ for Doppler-free media and as $`1/\sqrt{\kappa }`$ for Doppler-broadened media when $`\kappa 1`$. The change in $`\kappa `$ per unit length is

$`{\displaystyle \frac{d\kappa (x)}{dx}}`$ $`=\alpha _{DB}\kappa (x)`$
$`{\displaystyle \frac{\sqrt{6\pi }}{4}}{\displaystyle \frac{\gamma _0}{\mathrm{\Gamma }}}\alpha _0\sqrt{\kappa (x)},`$ (21)
and solving for $`\kappa `$ as a function of position,

$`\kappa (x)\left(\sqrt{\kappa _0}{\displaystyle \frac{\sqrt{6\pi }}{8}}{\displaystyle \frac{\gamma _0}{\mathrm{\Gamma }}}\alpha _0x\right)^2.`$ (22)
Note that the behavior of the saturation parameter as a function of distance is different in the Doppler-broadened case compared to the Doppler-free case \[cf. Eq. (14)\] where the saturation parameter falls approximately linearly with distance.

Taking the average of small-field rotation over the Doppler distribution gives

$`{\displaystyle \frac{d\phi }{d\mathrm{\Omega }dx}}|_{DB}`$ $`{\displaystyle \frac{1}{\sqrt{\pi }\mathrm{\Gamma }}}{\displaystyle _{\mathrm{}}^{\mathrm{}}}{\displaystyle \frac{d\phi }{d\mathrm{\Omega }dx}}\left(\mathrm{\Delta }\right)𝑑\mathrm{\Delta }`$
$`{\displaystyle \frac{\alpha _0}{\sqrt{\pi }\gamma \mathrm{\Gamma }}}{\displaystyle _{\mathrm{}}^{\mathrm{}}}{\displaystyle \frac{1}{\left(2\mathrm{\Delta }/\gamma _0\right)^2+2\kappa /3}}𝑑\mathrm{\Delta }`$
$`{\displaystyle \frac{\sqrt{6\pi }}{4}}{\displaystyle \frac{\gamma _0}{\mathrm{\Gamma }}}{\displaystyle \frac{\alpha _0}{\gamma \sqrt{\kappa }}}.`$ (23)
We see that the rotation per unit length scales as $`1/\sqrt{\kappa }`$, in contrast to the $`1/\kappa `$ scaling for the Doppler-free case \[similar to the situation with absorption, Eqs. (7,19)\]. This is because the number of atoms producing the effect is, in a sense, not fixed; with increasing light power, a larger fraction of the Doppler distribution is involved.

Substituting (22) into Eq. (23) to find the contribution to the slope as a function of position gives

$`{\displaystyle \frac{d\phi }{d\mathrm{\Omega }dx}}(x)|_{DB}`$ $`{\displaystyle \frac{2}{\gamma }}{\displaystyle \frac{\alpha _0}{\frac{8}{\sqrt{6\pi }}\frac{\mathrm{\Gamma }}{\gamma _0}\sqrt{\kappa _0}\alpha _0x}}.`$ (24)
It is interesting to note that while the light power and the rotation slope per unit length behave differently in the Doppler-broadened case compared to the Doppler-free case, Eq. (24) has the same functional form as for the Doppler-free case \[Eq. (15)\].

Integrating over the length of the medium, we obtain

$`{\displaystyle \frac{d\phi }{d\mathrm{\Omega }}}|_{DB}`$ $`{\displaystyle _0^{\mathrm{}}}{\displaystyle \frac{d\phi }{d\mathrm{\Omega }dx}}(x)|_{DB}dx`$
$`{\displaystyle \frac{2}{\gamma }}{\displaystyle _0^{\mathrm{}}}{\displaystyle \frac{\alpha _0}{\frac{8}{\sqrt{6\pi }}\frac{\mathrm{\Gamma }}{\gamma }\sqrt{\kappa _0}\alpha _0x}}𝑑x`$
$`={\displaystyle \frac{2}{\gamma }}\mathrm{ln}\left({\displaystyle \frac{\sqrt{\kappa _0}}{\sqrt{\kappa _0}\frac{\sqrt{6\pi }}{8}\frac{\gamma _0}{\mathrm{\Gamma }}\alpha _0\mathrm{}}}\right).`$ (25)
The behavior of (25) as a function of $`\alpha _0\mathrm{}`$ is qualitatively similar to that of Eq. (16) shown in Fig. 2. However, the dependence on $`\kappa _0`$ is different.

The magnetometric sensitivity is given by

$`(\delta B)_{DB}^1`$
$`={\displaystyle \frac{g\mu \omega }{\mathrm{}c}}{\displaystyle \frac{d\phi }{d\mathrm{\Omega }}}|_{DB}\sqrt{{\displaystyle \frac{2At}{\pi }}\gamma \kappa (\mathrm{})}`$
$`{\displaystyle \frac{g\mu \omega }{\mathrm{}c}}\sqrt{{\displaystyle \frac{8At}{\pi \gamma }}}\left(\sqrt{\kappa _0}{\displaystyle \frac{\sqrt{6\pi }}{8}}{\displaystyle \frac{\gamma _0}{\mathrm{\Gamma }}}\alpha _0\mathrm{}\right)\mathrm{ln}\left({\displaystyle \frac{\sqrt{\kappa _0}}{\sqrt{\kappa _0}\frac{\sqrt{6\pi }}{8}\frac{\gamma _0}{\mathrm{\Gamma }}\alpha _0\mathrm{}}}\right)`$
$`={\displaystyle \frac{g\mu \omega }{\mathrm{}c}}{\displaystyle \frac{\sqrt{3}}{2}}\sqrt{{\displaystyle \frac{At}{\gamma }}}{\displaystyle \frac{\gamma _0}{\mathrm{\Gamma }}}\alpha _0\mathrm{}\left(p1\right)\mathrm{ln}\left({\displaystyle \frac{p}{p1}}\right),`$ (26)
with the change of variables

$`\kappa _0={\displaystyle \frac{3\pi }{64}}\left(p{\displaystyle \frac{\gamma _0}{\mathrm{\Gamma }}}\alpha _0\mathrm{}\right)^2.`$ (27)
The factor

$`\left(p1\right)\mathrm{ln}\left({\displaystyle \frac{p}{p1}}\right)`$
goes to unity as $`p`$ goes to infinity; it is approximately 0.9 at $`p=5`$. (Further gain from increased power is minimal, and if the power becomes too high, the approximation of the Doppler width being much larger than the power-broadened width breaks down.) Thus, in this case,

$`(\delta B)_{DB}^1`$ $`0.8{\displaystyle \frac{g\mu \omega }{\mathrm{}c}}\sqrt{{\displaystyle \frac{At}{\gamma }}}{\displaystyle \frac{\gamma _0}{\mathrm{\Gamma }}}\alpha _0\mathrm{}`$
$`0.3{\displaystyle \frac{g\mu }{\mathrm{}}}\lambda \mathrm{}n{\displaystyle \frac{\gamma _0}{\mathrm{\Gamma }}}\sqrt{{\displaystyle \frac{At}{\gamma }}},`$ (28)
for sufficiently high $`\kappa _0`$.

This result, where sensitivity increases linearly with optical thickness, holds for the case where the power-broadened width $`\gamma _0\sqrt{\kappa (x)}`$ is smaller than the Doppler width for all $`x`$ within the sample, i.e.

$`1`$ $`\kappa _0\left(\mathrm{\Gamma }/\gamma _0\right)^2\mathrm{and}`$
$`\mathrm{\Gamma }/\gamma _0`$ $`\alpha _0\mathrm{}\left(\mathrm{\Gamma }/\gamma _0\right)^2,`$ (29)
since $`\alpha _0\mathrm{}`$ is related to $`\kappa _0`$ by Eq. (27). As power and optical thickness are increased beyond this range, i.e.

$`\kappa _0`$ $`\left(\mathrm{\Gamma }/\gamma _0\right)^2\mathrm{and}`$
$`\alpha _0\mathrm{}`$ $`\left(\mathrm{\Gamma }/\gamma _0\right)^2,`$ (30)
we obtain the Doppler-free case, where sensitivity increases as the square root of the thickness \[Eq. (18)\].

## V The general case

A numerical result can be obtained for the general case where the restrictions (29,30) on $`\kappa `$ and $`\alpha _0\mathrm{}`$ are removed. In a typical experiment, light power is $`1\mathrm{mW}`$, the laser beam diameter is $`0.1\mathrm{cm}`$, $`\lambda 800\mathrm{nm}`$, and $`\mathrm{\Gamma }/\gamma _060`$. Thus the effective ground state relaxation rate due to the transit of atoms through the laser beam is $`\gamma 2\pi 50\mathrm{kHz}`$ and the initial saturation parameter is $`\kappa _04\times 10^3`$. Here, as in a typical experimental procedure, the optical depth is varied (by changing atomic density) while the laser power is kept constant. Normalized transmission, differential small-field rotation, total small-field rotation, and magnetometric sensitivity are plotted in Fig. 3 as a function of optical depth.

For small optical depth, $`\kappa >\left(\mathrm{\Gamma }/\gamma _0\right)^2`$ and the medium is effectively Doppler-free. Transmission ($`\kappa `$) falls linearly until the transition to the Doppler-broadened case is made (dashed line). Then transmission falls quadratically until the linear regime is reached (dotted line), after which it falls exponentially. Differential small-field rotation ($`d\phi /\left[d\mathrm{\Omega }d\left(\alpha _0\mathrm{}\right)\right]\gamma `$) initially rises, as $`\kappa `$ falls and power broadening is reduced, until non-linear effects begin to turn off. (Linear optical rotation is neglected in this plot.) Since magnetometric sensitivity depends both on total optical rotation and transmission, an intermediate value for the optical depth produces the greatest sensitivity (solid line). Multiplying the normalized inverse sensitivity $`d\phi /d\mathrm{\Omega }\gamma \sqrt{\kappa }`$ by

$`2\sqrt{2\pi }{\displaystyle \frac{g\mu }{\mathrm{}}}{\displaystyle \frac{\sqrt{A}}{\lambda \sqrt{\gamma }}}10^8\left(\mathrm{G}/\sqrt{\mathrm{Hz}}\right)^1`$ (31)
gives the absolute magnitude of sensitivity, $`2\times 10^{10}\mathrm{G}/\sqrt{\mathrm{Hz}}`$ (we assume $`g=1`$). Although this sensitivity is not as high as could be achieved with low atomic density paraffin-coated cells ($`3\times 10^{12}\mathrm{G}/\sqrt{\mathrm{Hz}}`$ Budker et al. (2000a)), it is, nevertheless, sufficiently high to be of interest in practical applications Novikova and Welch (2002). In particular, the power-broadening of the magnetic field dependence of optical rotation at high light power provides an increased dynamic range for magnetometry over the low-power case. There are, however, techniques for shifting the narrow resonance obtained with a paraffin-coated cell to higher magnetic fields, e.g., frequency modulation of the laser light Budker et al. (2002).

## VI Large-field optical rotation

The maximum optical rotation with respect to magnetic field can be determined with a numerical calculation. Maximum rotation is plotted as a function of optical depth in Fig. 4, for the same parameters as used for Fig. 3. The rotation initially rises linearly, the same behavior seen in the experiment Novikova et al. (2001) mentioned in the introduction. For large optical depth, however, rotation in Fig. 4 then begins to rise more quickly, and finally saturates, whereas experimentally a slower increase and then a decrease in rotation is seen. This is evidence for an additional relaxation mechanism not accounted for in the present theory; in Ref. Novikova et al. (2001) it is attributed to the effect of radiation trapping.

## VII Conclusion

In conclusion, we have analyzed magnetometric sensitivity of NMOR measurements optimized with respect to light intensity in the case of negligible Doppler broadening, and in the case of large Doppler broadening. In the former case, we find that the sensitivity improves as the square root of optical density, while in the latter, it improves linearly. In the present discussion, we have neglected the effect of velocity-changing collisions, which makes this analysis not directly applicable to buffer-gas and anti-relaxation-coated cells. However, since there is full mixing between velocity components in these cells, one can expect that the sensitivity should scale as square root of optical density (if this quantity can be varied independently of the ground state relaxation rate).

The authors are grateful to D. F. Kimball, I. Novikova, V. V. Yashchuk, and M. Zolotorev for helpful discussions. This work has been supported by the Office of Naval Research (grant N00014-97-1-0214).


format: text

text:
# The first detection of weak gravitational shear in infrared observations: Abell 1689

## 1 Introduction

A cluster acting as a weak gravitational lens distorts the shapes of background galaxies by virtue of its tidal gravitational field. The seminal paper of Kaiser & Squires (1993) describes how to use this information to obtain a parameter-free reconstruction of a mass distribution. Their technique, and variants thereof, has been applied to optical observations of a number of clusters (e.g. Clowe et al. 2000; Hoekstra et al. 2000).

Prior to this work, there has been no published report on the detection of gravitational shear in infrared observations. For a given integration time, the attainable number density of background galaxies for use in weak lensing studies is usually much lower than for optical observations. In addition, infrared detectors currently lag behind optical detectors in their field-of-view making them even less efficient for wide-field studies. Nevertheless, with developments in detector technology, and new wide-field infrared cameras becoming available, this is an avenue to be explored. Gray et al. (2000) have considered the magnification induced depletion effect in near-infrared CIRSI observations of the cluster Abell 2219, in order to obtain best-fit SIS (singular isothermal sphere) and NFW (Navarro, Frenk & White 1996) models.

Here we report on J, H and Ks infrared observations in which we detect the gravitational shear signature of the cluster Abell 1689 ($`z=0.182`$), manifest in the distorted images of background galaxies. We start by outlining the observations, and the strategy for obtaining a catalogue of background galaxies for use in the weak lensing analysis. Then we describe how this catalogue was used to perform a mass reconstruction of the cluster. The $`M_{\mathrm{ap}}`$ statistic of Schneider (1996) is then calculated, along with a quantitative measure of its significance. Next we use Ks selected cluster members to determine the luminosity of the lens, and obtain its mass-to-light ratio by comparison with the mass derived from fitting an SIS profile to the weak lensing data. We conclude with a brief discussion of the results, and the implications for future work. For quantitative estimates, we assume a cosmology with $`\mathrm{\Omega }=0.3`$, $`\mathrm{\Lambda }=0.7`$ and $`H_0=70\mathrm{km}\mathrm{s}^1\mathrm{Mpc}^1`$, and that the background galaxy population is at $`z=1.0`$. The conversion between angular and linear scales at $`z=0.182`$ is $`1.^{}00.18\mathrm{Mpc}`$.

## 2 Observations and data reduction

Our data were obtained on the nights 2000 April 29th and 30th, using the SofI camera (Moorwood et al. 1998) on the ESO-NTT telescope at Cerro La Silla Observatory, Chile. The instrument was used in “Large Field Imaging” mode, where the field of view is $`4.^{}9`$ on a side, and the pixel scale is $`0.^{\prime \prime }29`$/pixel. Abell 1689 was observed in three filters, J, H and Ks. Individual integrations lasted between 10 and $`30\mathrm{s}`$ and several of these were co-added to form an individual image. Between images, the telescope was moved by up to $`2.^{}0`$, while always keeping the target within the field of view. This offset is somewhat larger than is typical in such observations; however, it was necessary given the high density of galaxies in the centre of A1689. The field of view is effectively increased to $`7.^{}0`$ (1440 pixels). During the observations the seeing was typically $`0.^{\prime \prime }8`$ and conditions were perfectly photometric, except during the early hours of the first night.

The data were reduced in the standard way. From each image, the zero-level offset was removed, the flat field correction was applied, and then an estimate of the sky from other images in the sequence was subtracted. The images were combined so that the $`\frac{S}{N}`$ (signal-to-noise) ratio of the final image was maximised and with linear registration relative to the first frame in each sequence.

The photometric calibration was done by observing several stars from the list of HST NICMOS photometric standards (Persson et al. 1998). Although there was some thin cirrus early in the first night, we were able to observe a sufficiently large number of standards to derive accurate instrumental zero points and atmospheric extinction coefficients in each of the three filters. The total integration times ($`t`$), the final image quality (IQ), the sensitivity limit (1-$`\sigma `$ noise limit in 1 square arcsecond) ($`SL`$) and extinction coefficients (A) of the combined images in the J, H and Ks bands are given in Table 1. The instrumental zero points were very stable over the two nights and the scatter in the fit of the instrumental magnitude versus airmass is less than 0.01 magnitudes in all three filters.

We did not correct for colour terms between the SofI instrumental system and that used for the HST NICMOS standards. Preliminary estimates suggest that these terms are all less than 0.02 magnitude. Nor did we correct the flat fields for illumination effects which cause the photometry to vary differentially across the field by up to 0.03 magnitudes. Thus, the accuracy of the relative and absolute photometry is around 0.03 magnitudes for both.

## 3 Obtaining a catalogue of background objects

The goal described in this section is the acquisition of a catalogue of background object positions and ellipticities, that can be used in the weak lensing analysis.

Our strategy was to start by using SExtractor (Bertin & Arnouts 1996) to detect objects on each of the J, H and Ks images separately. For all filters, an object detection and analysis threshold of at least 3 contiguous pixels above $`1\sigma `$ was required. Since the J image contained most detections, the H and Ks images were transformed to the J reference frame, using IRAF tasks (see Tody 1993 and references therein). SExtractor was then run in “double-image mode”, using the J image for detection (with an object detection and analysis threshold of at least 2 contiguous pixels above $`0.7\sigma `$) and each of the J, H and Ks images for measurements. For each object, magnitudes were measured using both a 15 pixel ($`4.^{\prime \prime }4`$) diameter aperture, and an outer isophote at 70% of the sky noise (J: 24.07, H: 23.11, Ks: 22.55 mag/sq. arcsec.). Aperture magnitudes are used for colour determination, and isophotal magnitudes are used for all other purposes. The LDAC tools111Freeware available from ftp://ftp.strw.leidenuniv.nl/pub/ldac/ were used to obtain a single catalogue of these sources.

Ellipticities for the objects were determined from the second moments of their surface brightness using a modified version of Nick Kaiser’s IMCAT222http://www.ifa.hawaii.edu/$``$kaiser/imcat software. Bright (but unsaturated) stars were selected on the basis of their half-light radii and were used to determine the point spread function (PSF) as a function of position in the images. A bi-cubic polynomial was fit to the stellar ellipticities (for 25 stars) and used to correct the ellipticities of the galaxies using the technique detailed in Kaiser et al. (1995). The smearing of the galaxy ellipticities by the now circular PSF was removed using the basic method outlined in Luppino and Kaiser (1997), but where the pre-seeing shear polarizability tensor was fit as a 5th order polynomial function of $`r_g`$, the Gaussian smoothing radius at which the galaxy achieves maximum significance from the background, and the post-seeing ellipticity. The resulting shear estimate has been shown to be good to a few percent accuracy for well behaved PSFs such as those in these images (Erben et al. 2001).

Stellar objects were filtered out, based on their locus in a magnitude-radius plot. Next, making use of the small dispersion in colour at a particular magnitude for cluster galaxies, likely cluster members were excised to exclude even the faintest cluster members, where the dispersion in the cluster sequence becomes larger.

Further, only galaxies fainter than $`\mathrm{m}_\mathrm{J}=16.6`$ and having a $`\frac{S}{N}>5`$ in J were retained. The $`\frac{S}{N}`$ limit corresponds to m$`{}_{\mathrm{J}}{}^{}=23.2`$, for the smallest objects of similar size to the PSF. Finally, only galaxies redder than the cluster sequence galaxies were kept in the catalogue (they are most probably at higher redshifts than the cluster; see for example the colour evolution plots in Poggianti 1997).

The selection process for the catalogue as outlined above is the most conservative of our trials, and the final number of sources used in the weak lensing analysis is 157, just over 3/arcmin2. In Fig. 1, these objects are indicated (open stars) on a colour-magnitude diagram for the observations. We experimented with various cuts, which induce only a marginal change in the results of the lensing analysis described in the next Section.

## 4 Weak Lensing Analysis

For the weak lensing analysis, we used the mean of the shear estimates determined from the separate J, H and Ks images, weighted by the $`\frac{S}{N}`$ of the object detections. The best-fit SIS model to the radial shear profile (centred on the brightest cluster galaxy) has a velocity dispersion $`\sigma _{1\mathrm{D}}=1030_{80}^{+70}\mathrm{km}\mathrm{s}^1`$, with a significance of 6.1 relative to the null model ($`\sigma _{1\mathrm{D}}=0\mathrm{km}\mathrm{s}^1`$). This value is consistent with the best-fit SIS to WFI optical data ($`\sigma _{1\mathrm{D}}=998_{42}^{+33}\mathrm{km}\mathrm{s}^1`$) presented in King et al. (2001). The best-fit NFW model has $`r_{200}=1.94\mathrm{Mpc}`$ and $`c=5.7`$, the significance and quality of fit being the same as the SIS. We checked that fitting a radial shear profile to the objects bluer than the cluster sequence results in an insignificant weak lensing signal.

### 4.1 Mass reconstruction

We use the non-parametric finite-field mass reconstruction algorithm detailed in Seitz & Schneider (2001). The reconstruction was performed with a smoothing scale of $`0.^{}54`$. The result is shown as solid contours in Fig. 2, normalised (necessary because of the mass-sheet degeneracy) so that the mean surface mass density on a $`0.^{}097`$ wide border is $`\overline{\kappa }=0.1`$. Note that the mass peak is consistent with that of the luminous infrared image.

Since the noise properties for non-parametric reconstructions are not easily interpreted, multiple reconstructions were performed during which the catalogue object positions were kept fixed and ellipticities were assigned at random (with replacement) from the ellipticity ensemble. This enabled us to get a handle on the noise properties, and confirms that the central mass concentration is not a noise artefact.

### 4.2 Aperture mass statistic

Next we consider the $`M_{\mathrm{ap}}`$ statistic, which measures the surface mass density convolved with a compensated Mexican Hat filter function $`𝒰`$. Equivalently, the suitably weighted ellipticities give a practical estimator of $`M_{\mathrm{ap}}`$

$$M_{\mathrm{ap}}(\vartheta )=\frac{\pi \eta ^2}{N}\underset{i}{}ϵ_{\mathrm{t}i}(\vartheta )𝒬(|\theta _i\vartheta |),$$
(1)
where the sum extends over $`N`$ galaxies, located at positions $`\theta _i`$, inside a (filter) radius $`\eta `$ centred on position $`\vartheta `$. The quantity $`ϵ_\mathrm{t}(\vartheta )`$ is the tangential component of the ellipticity with respect to position $`\vartheta `$, and $`𝒬`$ is a filter function related to $`𝒰`$. We use the simplest form of the filter functions from Schneider et al. (1998).

An asset of this statistic is that the noise properties are well understood, and the significance of a $`M_{\mathrm{ap}}`$ detection is readily calculated by comparing its value with that obtained from the same set of galaxies with random position angles. The signal to noise ratio $`\frac{S}{N}`$ for $`M_{\mathrm{ap}}`$$`(\vartheta )`$ is $`\frac{M_{\mathrm{ap}}}{\sqrt{\sigma ^2}}(\vartheta )`$ where $`\sigma ^2=M_{\mathrm{ap}}^2_{\mathrm{randomisations}}(\vartheta )`$. The $`rms`$ dispersion of $`M_{\mathrm{ap}}`$ is calculated in the absence of lensing, since the deviation from this approximation is very small in the weak lensing regime ; it gives an upper limit on $`\sigma `$. The $`\frac{S}{N}`$ of $`M_{\mathrm{ap}}`$ calculated for a filter scale $`\eta =1.^{}95`$ and for 1 000 randomisations, is shown as dashed contours in Fig. 2. $`\frac{S}{N}5`$ for the peak of the detection; its location is consistent with the peak of the non-parametric mass reconstruction and with the brightest cluster galaxy on the infrared image. Also, note that the $`M_{\mathrm{ap}}`$ detection is highly significant over a large area.

The calculation of $`\frac{S}{N}`$ in the previous paragraph relies on the fact that $`M_{\mathrm{ap}}`$($`\vartheta `$) has a Gaussian probability distribution. In addition, $`M_{\mathrm{ap}}`$ was calculated around the peak of the detection for $`4\times 10^7`$ randomisations of the background galaxy orientations, and a higher $`\frac{S}{N}`$ peak was obtained for only 18 of these cases. The probability that such a random distribution of galaxies would give a $`\frac{S}{N}`$ in excess of the observed peak value is therefore $`4.5\times 10^7`$. As another check, the position angles of the background galaxies were rotated through $`45^\mathrm{o}`$, and no highly significant $`M_{\mathrm{ap}}`$ detection was obtained.

## 5 Mass-to-light ratio

Selection of cluster members was performed in the Ks-band data, down to isophotal m$`{}_{\mathrm{Ks}}{}^{}=18.5`$ and using their tight colour-magnitude and colour-colour loci. Further, we restricted selection to objects within 500 pixels ($`2.^{}43`$) of the brightest cluster galaxy. Outside of this radius, the sky noise increases dramatically due to dithering causing a decrease in exposure time. We indicate the objects selected as cluster members on Fig. 1. The apparent magnitudes of these galaxies were converted to absolute magnitudes, and subsequently to luminosities using the evolutionary and k-corrections obtained from the single burst, solar metallicity, Scalo (1986) IMF spectral energy distribution contained in HyperZ (Bolzonella et al. 2000) (amounting to $`0.53\mathrm{mag}`$) and the extinction correction from our data ($`0.09\mathrm{mag}`$).

The differential K-band luminosity function of the cluster members is shown in Fig. 3, with error bars corresponding to the Poisson variance in each bin, along with the best-fit Schechter (1976) luminosity function

$$\mathrm{\Phi }(L)dL=\mathrm{\Phi }^{}\left(\frac{L}{L_{}}\right)^\alpha \mathrm{exp}\left(\frac{L}{L_{}}\right)\frac{dL}{L_{}}.$$
(2)
The nonlinear least-squares Marquardt-Levenberg algorithm was used to obtain the best-fit parameters for the faint slope $`\alpha =1.01\pm 0.14`$ and characteristic absolute magnitude $`\mathrm{M}_{}=24.7\pm 0.40`$ (i.e. $`\mathrm{m}_{}=14.6`$), with a reduced $`\chi ^2`$ of 1.3. Taking into account the difference in adopted cosmology, this value of $`\mathrm{M}_{}`$ is very similar to that obtained by Barger et al. (1996), for a sample of distant clusters.

The integrated luminosity (extrapolated to zero luminosity) for the best-fit Schechter function is 8.0$`\times 10^{12}L_{}`$ and the direct sum of galaxy luminosities (down to $`L_{\mathrm{min}}`$) is 8.8$`\times 10^{12}L_{}`$. These are consistent, given the errors of the fit parameters. We use the direct sum in the mass-to-light ratio estimate, since we reach well down the luminosity function: the fraction of light observed is approximately $`\mathrm{\Gamma }(\alpha +2,\frac{L_{\mathrm{min}}}{L_{}})/\mathrm{\Gamma }(\alpha +2)97\%`$. In order to estimate the error on the luminosity, first consider the error arising from the various cuts in the colour-magnitude and colour-colour planes, given that each object has a photometric uncertainty. We started with the original uncut catalogue, and created 100 000 mock catalogues, assigning the $`i^{\mathrm{th}}`$ object a flux Fi=Fi(mea)\+ $`g\mathrm{F}_{i(\mathrm{err})}`$, where $`g`$ is a uniform gaussian random variable, and Fi(mea), Fi(err) are the measured flux and its error, respectively. The magnitude and colour cuts were then applied to these catalogues, and the dispersion in the recalculated luminosity determined to be only $`0.35\%`$ of the total. Another very small error is the Poissonian error on the total flux, which amounts to $`0.05\%`$. The dominant error is the 2% error on the luminosity because of uncertainty in the Ks zero point. Hence for the cluster luminosity we take $`8.8\pm 0.18\times 10^{12}L_{}`$ (which is a lower limit because $`3\%`$ of the luminosity is expected to be below the Ks magnitude limit, and light outside the isophotal radii of galaxies is not accounted for).

The projected mass of the best-fit SIS to the infrared data within the same $`2.^{}43`$ radius is $`M`$=3.44$`{}_{0.51}{}^{+0.49}\times 10^{14}M_{}`$. This gives an upper-limit K-band mass-to-light ratio of $`39.1\pm 5.7M_{}/L_{}`$. The best-fit NFW model has a mass at this radius which is $``$10% higher.

## 6 Discussion and conclusions

We find that it is possible to fit a radial shear profile (the best-fit SIS model has $`\sigma _{1\mathrm{D}}=1030_{80}^{+70}\mathrm{km}\mathrm{s}^1`$), perform a mass reconstruction, and obtain a $`M_{\mathrm{ap}}`$$`\frac{S}{N}5`$ detection using only the infrared data on a relatively small field. The locations of the peaks of the reconstruction and of the $`M_{\mathrm{ap}}`$ detection are consistent with that of the brightest cluster galaxy in the infrared image.

The K-band luminosity function of the cluster galaxies is well described by a Schechter function, with a faint end slope $`\alpha =1.01\pm 0.14`$ and characteristic absolute magnitude $`\mathrm{M}_{}=24.7\pm 0.40`$. In addition, we find a K-band mass-to-light ratio of $`40M_{}/L_{}`$ in the cosmology adopted, in accord with the study of rich Abell clusters made by Uson & Boughn (1991).

The capability to perform routine wide-field infrared imaging is just around the corner with projects such as the UKIRT WFCAM, the CFHT WIRCAM and VISTA. Weak gravitational lensing observations of distorted galaxies in the infrared will provide another window through which to probe the mass distribution on cluster scales. Although the number density of background sources achievable is roughly an order of magnitude lower at infrared wavelengths, there are a few points worth noting:

* Colours of galaxies in the near-infrared are dominated by light from old evolved stars, so the well defined cluster sequence in the colour-magnitude diagram contains both ellipticals and (unlike the optical cluster sequence) spirals. This makes excising potential cluster members easier using infrared colours, as anything redder than the sequence is likely to be at a higher redshift, and anything bluer is probably at a lower redshift. The infrared selection of background objects can then be used as a consistency check of the optical selection, to assess the contamination from cluster members and foreground objects.
* With optical observations, measurement of the ellipticities of high redshift galaxies can be hampered by the presence of blue knots of star formation.
* The total K-band luminosity of the cluster members can be compared with a weak lensing mass, in order to obtain a mass-to-light ratio. Further, Broadhurst et al. (1992) underlined that the K-band luminosity is a good measure of the underlying stellar mass irrespective of how it assembled. If a high enough resolution mass map became available, we could investigate the correlation between any peaks in the total mass and stellar mass content.

Although infrared observations are not preferable to optical ones if the goal is simply to do a lensing analysis, this work highlights the complementary nature of this passband. Our next step is to combine the infrared and optical information. On a galaxy-by-galaxy basis, the selection of background galaxies can be examined. In addition, we can compare the ellipticities of the lensed galaxies determined from each data set. Finally, we will consider constraints from strong lensing features to obtain the most precise description of the mass distribution of this spectacular lens.

###### Acknowledgements.

This work was supported by the TMR Network “Gravitational Lensing: New Constraints on Cosmology and the Distribution of Dark Matter” of the EC under contract No. ERBFMRX-CT97-0172, and by the DFG and the CNRS. Thanks to Joan-Marc Miralles, Roser Pelló, Tom Broadhurst, Michael Hilker, Martin Altmann and Neil Trentham for useful discussions, and to Leonardo Vanzi for providing us with flat field images from the SofI calibration plan. Thanks also to the referee for very helpful comments on the manuscript.