text:
# Stochastic dynamics of the magneto-optical trap

## 1 Introduction

Magneto-optical traps (MOT) produce clouds of cooled atoms at temperature as low as the $`\mu `$K. The achievement of such clouds opened many perspectives, not only in the field of fundamental atomic physics, as e.g. in the domain of the atomic dynamics or the quantum chaosQC , but also leads to several potential applications, in particular the improvement by several orders of magnitude of the atomic clockclocks . The MOT is also the first stage to produce lower temperatures, in particular to obtain Bose-Einstein condensates nobel . Although the use of MOTs is relatively well mastered, some details of the experimental setup remain empirical, because of the existence in the cloud of instabilities that are not well understood. Some studies showed that when the trapping beams are misaligned, the cloud may be spatially altered and become unstable. In particular, ring-shaped clouds and chaos have been observed, and attributed to a vertex force shadow ; mis . However, the main instabilities encountered in the experiments concern well-aligned MOTs. In that case, the cloud, which is usually more or less ellipsoidal, has a complex irregular shape, with an inhomogeneous atomic density (fig 1). This shape may not be stable, and changes as a function of time: a typical example of these instabilities is shown in fig. 1.

There is a deep interest to identify the nature of these instabilities, in the aim to control them, and possibly to take advantage of them. To illustrate these points, let us remember that instabilities may originate in many mechanisms, which can be classified in two families: stochastic or deterministic. In the latter, called also deterministic chaos, the dynamics is described by a set of deterministic equations. When the number of equations is small (low-dimensional deterministic chaos), it has been shown that the dynamics can be controlled to reach various states that are not accessible otherwise, as unstable statescontrolstab  or periodic behaviorscontrolper . It is also well known that the description of such a dynamics gives access to numerous parameters, which are not measurable when the behavior is stationaryparamaccess . The cost for these informations is that there is no other way to reduce the instabilities. On the contrary, if the origin of the instabilities is stochastic, i.e. due to external noise, a reduction of instabilities is obtained by reducing the noise, but the instabilities have no physical meaning, in that sense that they are not intrinsic to the physics of the system.

In nousprl , a first type of instabilities, arising for low beam intensities, has been depicted. It was showed that these so-called “stochastic instabilities” are induced by the absorption of the trapping beams in the cloud. A model taking into account this so-called shadow effect showed that from a dynamical point of view, instabilities arise through a stochastic resonancelike phenomenon, namely the coherent resonance, linked to a Hopf bifurcation in the stationary solutions of the MOT. It is well known that such a bifurcation usually leads to periodic instabilities, and indeed, a recent study evidenced such instabilities, purely deterministic InstDet . However, a correct description of these self-oscillations required to modify the model in nousprl , in particular to take into account the spatial distribution of the cloud.

In the present paper, the results presented in nousprl  are detailed, and extended to a larger range of parameters, where new regimes appear, in particular a behavior without resonance frequency, contrary to the instabilities described in nousprl . Then the modified model introduced in InstDet  is accurately described, and the relative domains of appearance of stochastic and deterministic instabilities are discussed. The difference appears to be the result of the preeminent role of noise in some parameter range, in particular in the vicinity of the bifurcations. The stochastic instabilities predicted by the model are studied, and they are shown to be in agreement with those experimentally observed. In particular, an interpretation of the two types of stochastic instabilities, with and without resonance frequency, is found.

The paper is organized as follows. After this introduction, the section 2 describes the experimental setup. Then, a detailed analysis of the experimental observations of the stochastic instabilities is presented (section 3), showing in particular the existence of two types of instabilities. Section 4 is devoted to the construction of a simple 1D model, already presented in InstDet  in a less detailed way. In section 5, the stationary solutions of the model are discussed. Finally, in section 6, the noise induced dynamics instabilities is studied and compared to the experimental observations.

## 2 Experimental set-up

We work with a Cesium-atom MOT in the usual $`\sigma _+\sigma _{}`$ configuration, with three arms of two counter-propagating beams obtained from the same laser diode. The waist $`w_T`$ of the trap beams may be varied from typically $`3`$ to $`10`$ mm. Two configurations are possible: in the first one, all six beams are independent, by opposition to the second configuration, where counter-propagating beams result from the reflection of the three forward beams. In the last case, the intensity asymmetry resulting from the absorption of the forward beam by the cloud, generates a center-of-mass motion, while in the first case, instabilities are characterized by symmetrical bursts on the cloud shape, much more difficult to measure. However, as the nonlinearities involved in both cases are the same, we expect that the dynamics will be fundamentally of the same nature, and thus we choose the configuration with retro-reflected beams.

A full description of the unstable dynamics of the atomic cloud will be presented in the next section. However, to make easier the understanding of this paragraph, let us depict them briefly. As shown in the introduction (fig. 1), instabilities consist in a deformation of the spatial atomic distribution, leading to fluctuations of the shape of the cloud. Therefore, the relevant dynamical variables allowing us to describe instabilities, could be the shape of the cloud (i.e. for example the local velocities and atomic densities in the cloud). This type of description corresponds to a high dimensional model, associated with partial differential equations. Here, for the sake of simplicity, we choose to limit our description to the center of mass (CM) location $`𝐫`$, and the total number of atoms $`n`$ in the atomic cloud. This allows us to model the system with ordinary differential equations, and reduces the dimension to seven, and even three in a 1D model. As it is shown in the following, the use of this description appears to be sufficient to understand the main mechanisms of the instabilities.

To measure $`𝐫`$ and $`n`$, we used two 4-quadrant photodiodes (4QP) forming an orthogonal dihedral and measuring the fluorescence of the cloud. The differential signal of the 4QP allows us to monitor the motion of the CM. Using only two 4QPs does not allow us to reconstruct the actual 3D motion, but gives access to projections of this motion on two different axes. This prevents the measure from line-of-sight effects due to the optical thickness of the cloud. We checked that whatever the type of dynamical behavior, the motion components $`r`$ recorded by both 4QP have the same properties and are qualitatively identical. In the following, the “CM motion” refers to a component $`r`$ of this motion recorded by one of the 4QPs. The second dynamical variable, namely the number of atoms inside the cloud, is deduced from the total signal received from the 4QPs. In addition to the 4QPs, two video cameras monitor the shape of the cloud. Because of their poor resolution as compared with the 4QPs, the video cameras are not used to record the dynamical variables. They have been essentially used in the first stages of the experiment, to control that there is no discrepancy between the shape dynamics and the CM dynamics.

Instabilities depend on the MOT parameters. Among these parameters, some are easily controllable, and will be referred in the following as control parameters. They are the detuning $`\mathrm{\Delta }_0`$ of the MOT, the magnetic field gradient $`G`$, the MOT beam intensities $`I_1`$ and the repumper laser intensity $`I_{rep}`$. Other parameters cannot be considered as control parameters, because they are not easily controllable or measurable in our experimental setup. Among these parameters, the alignment of the MOT beams appeared to be crucial in the experiment. We limit the present study to the case where beams are aligned. When misaligned, most of the dynamical characteristics of the cloud change qualitatively mis . The MOT beam waists play also a main role on the dynamics, but in practice, they cannot be changed easily independently from the other parameters, in particular $`I_1`$. However, two different values have been used in the experiments, and their effects on the dynamics are well understood. Finally, the vapor pressure in the cell has probably also a large influence on the dynamics. Unfortunately, this parameter is not easily measurable in our experimental setup. Moreover, as it is shown later, it does not appear explicitly in the model, but its impact on the dynamics may be estimated through the equilibrium population of the cloud. The parameter ranges explored in the present experiment are summarized in Tab. 1.

## 3 Experimental results

Instabilities have been described in nousprl  for a given set of parameters. In the following, we extend this description for the whole range of parameters where stochastic instabilities appear. A first fundamental control parameter is the MOT beam intensity, whose value determines the type of observed instabilities. At low MOT beam intensities, typically less than $`10I_S`$ ($`I_S=1.1`$ mW/cm2 is the saturation intensity), the cloud exhibits S instabilities (S stands for Stochastic). For MOT beam intensities larger than $`10I_S`$, C (for Cyclic) instabilities appear InstDet . As the aim of the present paper is to discuss about the stochastic instabilities, we keep a detailed presentation of the deterministic C instabilities for another paper. However, the two types of instabilities cannot be completely separated, and the domains of appearance of both types of instabilities will be discussed.

As discussed above, S instabilities appear at low MOT beam intensities, and are characterized by large fluctuations of the shape of the cloud appearing in a limited range of the parameters. From the experimental point of view, this last point is essential, because it is at the origin of the introduction of the concept of instabilities of the MOT. Indeed, if the noisy dynamics is the same whatever the parameter values, it is clear that the problem becomes as trivial as the reduction of technical noise in an experiment. On the contrary, the motion of the MOT grows suddenly in a narrow range of the parameters, as it is illustrated in fig. 2, where the amplitude $`\mathrm{\Delta }r`$ of the $`r`$ fluctuations is represented as a function of the control parameter $`\mathrm{\Delta }_0`$. Far from the unstable area, e.g. at large $`\left|\mathrm{\Delta }_0\right|`$, the cloud is stable in shape and density. When the control parameter is tuned to the unstable area, instabilities appear progressively, and the amplitude grows until a maximum $`\mathrm{\Delta }r_{\mathrm{max}}`$ at $`\mathrm{\Delta }_{0\mathrm{max}}1.95`$ ($`\mathrm{\Delta }_0`$ is expressed in units of the natural width $`\mathrm{\Gamma }`$ of the atomic transition). There is no abrupt boundary between the stable and unstable areas: the limits given below corresponds to $`\mathrm{\Delta }r=\mathrm{\Delta }r_{\mathrm{max}}/10`$. In this case, for the set of parameters of fig. 2, instabilities appear for $`2.5<\mathrm{\Delta }_0<1.7`$ (for $`\mathrm{\Delta }_0>1.7`$, the population in the cloud vanishes). The unstable range depends on the other parameters, such as the beam intensities, the vapor pressure or the magnetic field gradient. For example, for a less populated cloud, due e.g. to a different vapor pressure in the cell, all other parameters being the same as previously, instabilities will appear at smaller detuning. However, whatever the parameters used in the experiments, in the range given in Tab. 1, we have $`3\mathrm{\Delta }_{0\mathrm{max}}1`$, and instabilities never occur on a range larger than $`1`$.

Fig. 3 shows a typical unstable behavior of $`r`$. It appears as an erratic signal, with a flat spectrum (fig. 4). Low frequency significant components appear typically for values smaller than $`\nu _n2`$ Hz (fig. 4b), and the dynamics is essentially along the first bisector of the three forward beams. The behavior of $`r`$ and $`n`$ are similar, with a cross correlation coefficient larger than 0.8. To determine if these instabilities have a deterministic origin, several tools are offered through the nonlinear dynamical analysis of the time series. As the MOT is dissipative, a deterministic dynamics should have an attractor, which can be reconstructed easily for a low dimensional dynamics. The result, not presented here, appears as a set of randomly distributed points: in particular, it does not present any fine structure. Poincaré section and 1D maps confirm this absence of order in the dynamics. This could be due to a lack of resolution of the measures, but the general shape of the trajectories rather suggests that the behavior is stochastic (or chaotic with a high dimensional dynamics). Because this behavior appears to be a stochastic dynamics with only low frequency components, it will be referred in the following as SL instabilities.

In some situations corresponding to given ranges of parameters nousprl , the $`𝐫`$ dynamics is altered by the appearance of spontaneous large amplitude oscillation-like bursts. As illustrated in fig. 5, the signal inside the bursts is not periodic, although it is clearly dominated by a given frequency. These bursts are relatively scarce, and the global shape of the signal remains that of Fig. 3. The bursts appear in fact as the most spectacular consequence of a deeper change of the dynamics, which is the appearance of a second characteristic time. This appears clearly in the spectrum of $`r`$ (Fig. 6a) as a peak centered at a frequency $`\nu _r`$, which depends on the parameter values, but ranges typically between 10 and 100 Hz. The other characteristics of the dynamics are not modified by the existence of bursts. In particular, the spectrum still exhibits the low frequency component below $`\nu _n`$, and the bursts do not appear on the $`n`$ behavior (Fig. 5b and 6b). The low frequency component dynamics of $`r`$ and the dynamics of $`n`$ remain correlated. This behavior will be referred in the following as SH instabilities.

When the trap beam intensity $`I_1`$ is increased, S instabilities still exist, but they are progressively superseded by C instabilities. These instabilities differ drastically from S onesInstDet : they can be either periodic or erratic, but in both cases, they are cyclic, and the motion amplitude is much larger. However, as for S instabilities, C instabilities exist in a limited range of $`\mathrm{\Delta }_0`$, which value depends on the other parameters. But it is systematically between $`\mathrm{\Delta }_0=3`$ and resonance, and the maximum detuning range is of the order of $`1`$.

As $`I_1`$ is increased, the disappearance of S instabilities occurs progressively, in favor of C instabilities. For intermediate values of $`I_1`$, both types of instabilities exist. Their typical distribution versus $`\mathrm{\Delta }_0`$ is illustrated in Fig. 7: far from resonance, the cloud is stable; as the resonance is approached, S instabilities appear for a detuning $`\mathrm{\Delta }_0=\mathrm{\Delta }_1`$. Then C instabilities appear in $`\mathrm{\Delta }_2>\mathrm{\Delta }_1`$. If the detuning is still increased, C instabilities disappear in $`\mathrm{\Delta }_3`$ at the benefit of a stable behavior. Finally, the cloud vanishes in $`\mathrm{\Delta }_4`$. As $`I_1`$ is increased, the width $`\delta _{12}=\mathrm{\Delta }_2\mathrm{\Delta }_1`$ decreases in favor of the interval $`\delta _{23}=\mathrm{\Delta }_3\mathrm{\Delta }_2`$, while the total unstable interval $`\delta _{13}=\mathrm{\Delta }_3\mathrm{\Delta }_1`$ remains more or less constant. When C instabilities merge for $`I_1=4I_S`$, they appear on a narrow interval $`\delta _{23}0`$. This interval increases rapidly until $`I_1=7.5I_S`$ and $`\delta _{23}=0.8`$. For $`I_1>7.5I_S`$, $`\delta _{23}`$ increases more slowly, to reach the value of $`\delta _{23}=1`$ in $`I_1=20I_S`$, where S instabilities have completely disappeared.

To conclude this section, let us summarize the main properties of the instabilities: they appear to be stochastic at a time scale larger than 0.5 s. They may exhibit a second dynamical component with time scale smaller than $`10^2`$ s and acting only on $`𝐫`$. The low frequency component corresponds to a CM motion along the first bisector of the three forward beams. Finally, they appear only in a limited range of the parameters, suggesting a nonlinear origin of the instabilities. This last point was also suggested by different tests showing that beam phase fluctuations are not at the origin of instabilities nousprl .

## 4 Model

To understand the dynamics observed in the experiments, we build a phenomenological model taking into account the shadow effect. The aim here is not to model as finely as possible the experimental system, but on the contrary to make a model as simple as possible, enlightening the fundamental mechanisms leading to the instabilities.

The model that we used here has already been described in InstDet . It is based on the shadow effect induced by the intensity gradients produced by the absorption of the trapping laser beams in the cloud shadow ; dalibard . The shadow effect generates a force compressing the cloud. In the case where the counter-propagating beams result from the reflection of the three forward beams, this compression is accompanied by a displacement of the CM along the first bisector of the three forward beams: indeed, the backward beams are less intense than the forward ones because of the absorption in the cloud, and so the latter literally push the cloud off its equilibrium location. The main role of the shadow effect in the instabilities is suggested in the experiments (i) by the observed correlation between the $`n`$ and slow CM dynamics and (ii) by the fact that the slow dynamical component of $`r`$ displaces the cloud along the first bisector of the three forward beams. The fundamental role of the shadow effect in coupling the population and the CM position of cold atom clouds was already reported in bistab . In that experiment, bistability in the population dynamics was observed by perturbing the atomic cloud with a highly-focalized laser beam, and shadow effect is shown to be the dominant nonlinearity.

We built a 1D model taking into account the shadow effect in a MOT where counter-propagating beams result from the reflection of the three forward beams. As we are interested here in the collective motion, we use as variables the location $`z`$ of the center of mass of the cloud along the unique axis $`z`$ of the system, and the number of atoms $`n`$ inside the cloud. The origin of $`z`$ coincides with the “trap center”, that is, the zero of the magnetic field. Thus, the motion of $`z`$ may be described by the equation:

$$\frac{d^2z}{dt^2}=\frac{1}{M}F_T$$
(1)
where $`M`$ is the mass of the cloud and $`F_T`$ the global force exerted on the atoms by the two counterpropagating beams. To evaluate $`F_T`$, we assume a multiple scattering regime, i.e. a constant atomic density $`\rho `$ in the cloud. The atoms are distributed between $`z_{\mathrm{min}}`$ and $`z_{\mathrm{max}}`$, and the quantity $`\mathrm{\Delta }z=z_{\mathrm{max}}z_{\mathrm{min}}`$ is the longitudinal size of the cloud. We introduce a constant cross section $`S`$, allowing us to connect $`\mathrm{\Delta }z`$ to $`n`$:

$$n=\rho S\mathrm{\Delta }z$$
(2)
The repulsive force induced by the multiple scattering is an internal force, and thus does not affect directly the center of mass motion. It is also the case for the attractive force induced by the shadow effect, which acts on $`F_T`$ only through the asymmetry between the forward and backward beams. The forward beam is polarized $`\sigma _+`$ and has an input intensity $`I_1`$. After the crossing of the cloud, its output intensity is $`I_2`$. This is also the input intensity of the backward beam, which is polarized $`\sigma _{}`$. After the crossing of the atomic cloud, the output intensity of the backward beam is $`I_3`$. We have obviously $`I_1>I_2>I_3`$. $`F_T`$ is proportional to the number of photons absorbed per unit time, i.e. $`S(I_1I_2)/\mathrm{}\omega _L`$ for the forward beam and $`S(I_2I_3)/\mathrm{}\omega _L`$ for the backward beam ($`\mathrm{}\omega _L`$ is the energy of one photon, $`\omega _L`$ being the laser angular frequency). $`F_T`$ is obtained by multiplying the difference between both quantities by $`\mathrm{}k_L`$

$$F_T=\frac{S}{c}(I_12I_2+I_3)$$
(3)
To evaluate $`I_2`$ and $`I_3`$, we need to solve the propagation equation of light inside the cloud. To simplify the calculations, we choose to model a $`F_g=0F_e=1`$ transition: this is the simplest one leading to a magneto-Doppler effect. This choice does not allow us to reproduce the sub-Doppler effects, neither the photon redistribution between the two waves, but as shown in InstDet , it allows us to explain most of the experimental behaviors.

The $`z`$ axis is chosen as quantification axis, and we introduce the ground state $`|g`$ and the three Zeeman sub-levels of the excited state $`|e_0`$ et $`|e_\pm `$. Because of the kinetic momentum conservation, $`|e_0`$ is not coupled to the light field, and thus, the model is reduced to a three level system, $`\{|g,|e_{},|e_+\}`$. The interaction with light is governed by the Rabi frequency $`\mathrm{\Omega }_\pm `$ and the effective detunings $`\mathrm{\Delta }_\pm `$, taking into account the Doppler et Zeeman variations:

$$\mathrm{\Delta }_\pm =\frac{1}{\mathrm{\Gamma }}(\mathrm{\Delta }_0kv\omega _B^{}z)=\frac{1}{\mathrm{\Gamma }}(\mathrm{\Delta }_0\pm 2\delta )$$
(4)
where $`k`$ is the wave vector of light and $`\omega _B^{}`$ the Zeeman shift, measured in angular frequency by unit of length. The value of the atomic variables is given by the stationary solution of the optical Bloch equations for the density matrix $`\sigma `$:

$$\dot{\sigma }=\frac{1}{i\mathrm{}}[H,\sigma ]+\dot{\sigma }_{relax}$$
(5)
The hamiltonian $`H`$ of coupling with light may be written as a function of $`\mathrm{\Delta }_\pm `$ and $`\mathrm{\Omega }_\pm `$ in the $`\{|g,|e_{},|e_+\}`$ basis:

$$H=\mathrm{}\left(\begin{array}{ccc}0& \mathrm{\Omega }_{}/2& \mathrm{\Omega }_+/2\\ \mathrm{\Omega }_{}^{}/2& \mathrm{\Delta }_{}& 0\\ \mathrm{\Omega }_+^{}/2& 0& \mathrm{\Delta }_+\end{array}\right)$$
(6)
The relaxation rates $`\dot{\sigma }_{relax}`$ are different for the excited states, the optical coherences and the ground state:

$`(\dot{\sigma _{ee}})_{relax}`$ $`=`$ $`\mathrm{\Gamma }\sigma _{ee}`$ (7a)
$`(\dot{\sigma _{eg}})_{relax}`$ $`=`$ $`{\displaystyle \frac{\mathrm{\Gamma }}{2}}\sigma _{eg}`$ (7b)
$`(\dot{\sigma _{gg}})_{relax}`$ $`=`$ $`\mathrm{\Gamma }(\sigma _{++}+\sigma _{})`$ (7c)
The time scale of the internal variables is $`\mathrm{\Gamma }^1`$, and so is much shorter than that of the external variables $`z`$ and $`v`$. Thus the stationary solution of Eq. (5) can be used. The solution for the populations of the excited states $`\mathrm{\Pi }_\pm =\sigma _{\pm \pm }`$ can be obtained analytically. On the other hand, as the redistribution of photons between the beams is forbidden, $`\mathrm{\Pi }_\pm `$ is directly proportional to the absorption rate of the beams, and the intensity gradient can be written:
$$\frac{dI_\pm }{dz}=\mathrm{\Gamma }\mathrm{}\omega _L\rho \mathrm{\Pi }_\pm $$
(8)
By injecting the analytical solution of $`\mathrm{\Pi }_\pm `$ in this equation, we obtain the following equations of propagation:

$`{\displaystyle \frac{dI_+}{dz}}`$ $`=`$ $`\mathrm{\Gamma }\mathrm{}\omega \rho {\displaystyle \frac{\theta \omega _+\omega _{}+C_{}\omega _+}{1A_+\omega _+A_{}\omega _{}+3\theta \omega _+\omega _{}}}`$ (9a)
$`{\displaystyle \frac{dI_{}}{dz}}`$ $`=`$ $`\mathrm{\Gamma }\mathrm{}\omega \rho {\displaystyle \frac{\theta \omega _+\omega _{}+C_+\omega _{}}{1A_+\omega _+A_{}\omega _{}+3\theta \omega _+\omega _{}}}`$ (9b)
where:
$`\omega _\pm `$ $`=`$ $`{\displaystyle \frac{I_\pm }{2\gamma }}`$ (10a)
$`\gamma `$ $`=`$ $`{\displaystyle \frac{1}{4}}\left[\left(\alpha +\beta _++\beta _{}\right)^2+\left({\displaystyle \frac{\beta _+\beta _{}}{\mathrm{\Delta }_0}}\right)^2\right]`$ (10b)
$`\alpha `$ $`=`$ $`{\displaystyle \frac{1+2\mu }{4}}+\mathrm{\Delta }_+\mathrm{\Delta }_{}+{\displaystyle \frac{\delta ^2}{2}}\left(\mu 1\right)`$ (10c)
$`\beta _\pm `$ $`=`$ $`{\displaystyle \frac{\alpha }{2}}\pm \left(\mathrm{\Delta }_0\pm \delta \right)\left({\displaystyle \frac{\delta }{2}}+\mathrm{\Delta }_0\left(\mu _+\mu _{}\right)\right)`$ (10d)
$`\mu _\pm `$ $`=`$ $`{\displaystyle \frac{I_\pm }{4\left(1+\delta ^2\right)}}`$ (10e)
$`\mu `$ $`=`$ $`\mu _++\mu _{}`$ (10f)
$`\theta `$ $`=`$ $`\beta _+\beta _{}+\alpha \left(\beta _+\mu _++\beta _{}\mu _{}\right)`$ (10g)
$`A_\pm `$ $`=`$ $`2\beta _{}+\alpha \left(2\mu _\pm \mu _{}\right)`$ (10h)
$`C_\pm `$ $`=`$ $`\beta _\pm \pm \alpha \left(\mu _+\mu _{}\right)`$ (10i)
These equations appear as the ratio of two polynomials of high order, and an intuitive interpretation of this result is difficult. However, a numerical resolution of these equations leads to a rigorous evaluation of $`F_T`$.
The cloud population dynamics is modeled by a “feed-loss” rate equation feedloss :

$$\frac{dn}{dt}=B\left(n_en\right)$$
(11)
where we have introduce the population relaxation $`B`$ and the atom number in the cloud at equilibrium $`n_e`$, which is linked to the loading rate $`L`$ by $`L=Bn_e`$. $`n_e`$ is assumed to depend on the CM location, to take into account the losses variation when the cloud moves from the trap center. We do not know the exact form of this dependence, because this variation may have several origins. However, we can suppose that the main contribution comes from the transverse distribution of the trap laser beam, which is gaussian. For sake of simplicity, we keep only the first terms from its Taylor’s series, i.e. the quadratic term. In fact, we checked that the exact dependence of $`n_e`$ does not change drastically the results given below. So we write:

$$n_e=n_0\left[1\left(\frac{z}{z_0}\right)^2\right]$$
(12)
where $`n_0`$ is the equilibrium cloud population at the trap center and $`z_0`$ a characteristic length. Atoms in $`\left|z\right|>\left|z_0\right|`$ are considered as lost: this is taken into consideration when $`n`$ is deduced from Eq. 2.

Finally, we introduce the reduced variables $`Z=z/z_0`$, $`V=v/v_r`$ and $`N=n/n_0`$, where $`v_r`$ is the recoil velocity ($`v_r=\mathrm{}k/m`$), and we obtain the following autonomous system of equations:

$`{\displaystyle \frac{dZ}{dt}}`$ $`=`$ $`V{\displaystyle \frac{v_r}{z_0}}`$ (13a)
$`{\displaystyle \frac{dV}{dt}}`$ $`=`$ $`{\displaystyle \frac{1}{Mv_r}}F_T`$ (13b)
$`{\displaystyle \frac{dN}{dt}}`$ $`=`$ $`B\left(1Z^2N\right)`$ (13c)
Most of the theoretical parameters are the exact counterpart of the experimental parameters, as e.g. the magnetic field gradient or the beam intensities. In this case, we used in the model the same values as those of Table 1. It is not the case for all parameters, either because of the simplicity of the model or because they cannot be measured easily in the experiment. For example, $`n_0`$ has not a simple experimental counterpart, but depends on several experimental parameters, as e.g. the repumping laser intensity or the vapor pressure in the cell. For this reason, a large interval of $`n_0`$ values has been used for the simulations (Table 2). The density $`\rho `$ and the cross sectional area $`S`$, which play the same role, have a meaning only in the context of a 1D model, while they correspond to variables in the experiments. They are fixed in the simulations at experimental averaged values, and they have been varied on a wide range to check that their value is not critical. Finally, the parameter $`z_0`$ has not exact experimental counterpart, as it is linked to both the trap beam waist and intensities. Indeed, the relevant size for the atoms is not the beam waist, corresponding to an intensity decreased by a ratio $`e^2`$ as compared to the center of the beam, but rather the location where the local beam intensity decreases under $`I_S`$. This value is much larger than $`w_0`$ for intense beams.

In order to check if the model can be more simplified, we tried to reproduce the unstable dynamics of the cloud with the first terms of the Taylor’s series of Eqs. (9a) et (9b), but we needed to keep several orders and did not obtain simpler equations. Another possible approximation concerns the different terms appearing in $`\mathrm{\Delta }_\pm `$, in Eq. 4. Using the values given in table II, one sees easily that the $`\mathrm{\Delta }_0`$ term and the Zeeman shift are of the same order of magnitude, while the Doppler shift goes to zero at equilibrium, but can be the largest term out of equilibrium. In these conditions, no approximation is possible. Therefore, we use the equations (13) in the calculations.

To perform the comparison between the experiments and the present model, we need to study the behavior of the system when noise is added. Noise can be added in Eqs 13 on any parameter, and as we were not able in our experiments to identify clearly the main source of noise, we tried theoretically several parameters, as the beam intensity $`I_1`$ or the equilibrium population $`n_0`$. In this case, the stochastic model corresponds to Eqs 13 where the parameter $`I_1`$ (resp. $`n_0`$) is replaced by $`\left(1+\zeta \right)I_1`$ (resp. $`\left(1+\zeta \right)n_0`$), where $`\zeta (t)`$ is the noise component. We also simulated several types of noise: gaussian white noise, but also colored noise with different distributions of the frequencies. All configurations give identical results: on the one hand, the choice of the noisy parameter is not critical; on the other hand, the spectrum of noise does not alter the response spectrum, except obviously that it depends on the relative weight of the different frequencies in the noise. Thus, for sake of clarity, all the results reported in the following have been obtained by applying with gaussian white noise on $`I_1`$.

## 5 Stationary solutions

In the previous section, we built a relatively simple model where the cloud is described by a set of three equations that do not depend explicitly on the time. Such a system is known to be able to exhibit a complex dynamics, including chaos, which could explain the experimental instabilities. To know if such a complex dynamics occur in our conditions, the first step is to evaluate the stability of the stationary solutions. The stationary solutions ($`V_s`$, $`Z_s`$, $`N_s`$) of the model with the shadow effect are easily deduced from Eq.13 when the left side is put to zero. Eq. 13a gives immediately $`V_s=0`$. From Eq. 13c, one finds that the stationary solutions $`Z_s`$ of the CM location and $`N_s`$ of the population are linked by the simple expression:

$$N_s=1Z_s^2$$
(14)
Therefore, the discussion is reduced to the analysis of $`Z_s`$. $`Z_s`$ is the solution of $`F_T=0`$ in Eqs 13b. This equation can be resolved numerically: its global shape is illustrated in Fig. 8, where it is plotted as a function of $`\mathrm{\Delta }_0`$ and $`n_0`$.

A first characteristic of the $`Z_s`$ diagram is that $`Z_s`$ goes to $`1`$ (i.e. $`N_s`$ goes to zero) at resonance. This vanishing of the cloud, also observed in the experiments, is a well-known consequence of the inefficiency of the Doppler cooling close to resonance. It is here enhanced by the shadow effect and the displacement of the cloud. The disappearance occurs suddenly for small $`n_0`$, and becomes softer as $`n_0`$ increases. An interesting point is that for small $`n_0`$, the abrupt increase of $`Z_s`$ is linked to a very narrow bistable cycle. As $`n_0`$ is increased, the bistable cycle shifts towards smaller $`Z_s`$ and smaller $`\mathrm{\Delta }_0`$, and the vanishing of $`Z_s`$ becomes progressive.

The main characteristic of the $`Z_s`$ diagram is the presence of several abrupt slope changes in the stationary solutions, leading to a fold in the parameter space. The shape of the fold depends on the parameters, in particular on $`n_0`$. Fig. 9 shows four examples corresponding to situations leading to basically different atomic dynamics. For $`n_0=0.5\times 10^8`$ (fig. 9a), $`Z_s`$ increases smoothly with $`\mathrm{\Delta }_0`$ (i.e. $`N_s`$ decreases slowly). The vanishing of the cloud through a narrow bistable cycle is not visible on the graph, as it occurs closer from resonance. As $`n_0`$ increases, the bistable cycle appears for smaller $`Z_s`$ (and thus larger $`N_s`$), and becomes physically significant. Fig. 9b shows $`Z_s`$ for $`n_0`$ $`=2.5\times 10^8`$ and a bistable cycle for $`0.3\mathrm{\Delta }_00.25`$. If $`n_0`$ is further increased, the bistable cycle disappears, but it remains a fold corresponding to two abrupt slope changes of $`Z_s`$ versus $`\mathrm{\Delta }_0`$ (fig. 9c, $`n_0`$ $`=3.4\times 10^8`$). If $`n_0`$ is still increased, the fold remains, but it becomes smoother (Fig. 9d for $`n_0=`$ $`4\times 10^8`$).

The dynamics of the cloud is determined by the stability of these stationary solutions. In particular, if no stationary solutions are stable, complex dynamics could be obtained. The stability of the above stationary solutions is evaluated through a linear stability analysis, which associates to each stationary solution its three eigenvalues, corresponding to the stability following its three eigendirections in the 3D phase space of our model. The real part of the eigenvalue, corresponding to a damping rate, determines the stability (stable if negative). The imaginary part, when different from zero, is associated to an angular eigenfrequency, also called relaxation frequency, which play a main role in the dynamics. A pleasant – and simple – way to describe the stationary solutions is to use their phase space representation, where each stationary solution corresponds to a fixed point with its properties depending on its eigenvalues. Let us remember that the standard terminology distinguishes the stable node (all eigenvalues real and negative), the stable focus (all real parts negative, two eigenvalues complex conjugate), the saddle node (all eigenvalues real, at least one positive) and the saddle focus (at least one real part positive, two eigenvalues complex conjugate). For sake of simplicity, this terminology will be used in the following.

The set of points of vertical tangency in $`Z_s`$ determines a line which delimits the unstable stationary solution of the medium branch of the bistability cycle. Linear stability analysis shows that the fixed point in that area is a saddle-node with real eigenvalues, two being positive and one negative. But the unstable stationary solutions extend beyond this area. In particular, the stationary solutions can also be unstable on a part of the upper branch of the bistable cycle, and even outside the bistable cycle, when the stationary solution is unique. This is illustrated on fig. 8 where the nature of the stationary solution is indicated by a level of gray and a code. In the SN zone, the fixed point associated with the stationary solution is a saddle node, i.e. the stationary solution is unstable, and its three eigenvalues are real, two being positive. In the SF zone, the fixed point is a saddle focus, i.e. the stationary solution is unstable, and two of its eigenvalues are complex with a positive real part, and the third one is real negative. In the F zone, the fixed point is a stable focus, i.e. the stationary solution is stable, and two of its eigenvalues are complex with a negative real part, the third one is real negative. Finally, in the N zone, the fixed point is a stable node, i.e. the stationary solution is stable and all eigenvalues are real negative. The dashed line indicates the location of the bifurcation, i.e. the transition from stable to unstable stationary solutions. In most cases, it occurs from stable focus to saddle focus, through a super-critical Hopf bifurcation.

As the detuning is varied, four typical situations may occur, already illustrated in Fig. 9. For small $`n_0`$ (Fig. 9a), $`Z_s`$ is always stable, and its dependence versus $`\mathrm{\Delta }_0`$ is almost flat (except very close from resonance): we expect a stationary cloud slightly moving with the detuning. In the bistability area (Fig. 9b), the central branch of the bistability cycle is unstable, as usual in such a situation, but the upper branch is also partly unstable, reducing drastically the parameter range where the system is effectively bistable. However, a narrow bistable zone should be observed. As $`n_0`$ is still increased, only one stationary solution remains (Figs 9c), which is unstable on the fold.

Finally, for large $`n_0`$ (Fig. 9d), the solution is always stable (F zone), except very close from resonance. Thus we expect to observe a stationary cloud moving with $`\mathrm{\Delta }_0`$. However, because of the fold, this moving is non linear, and on the fold, for $`0.45<`$ $`\mathrm{\Delta }_0<0.35`$, $`Z_S`$ is more sensitive to the $`\mathrm{\Delta }_0`$ value, because the slope here is larger. On the other hand, a closer look to the eigenvalues (Fig. 10) shows that at the level of the fold, the real part $`\lambda `$ of the complex eigenvalues $`\lambda \pm i\omega `$ is close to zero, while the eigenfrequency crosses a minimum. This is more evident on Fig. 11, where the same situation is illustrated for a different set of parameters. The vanishing damping rate $`\left|\lambda \right|`$ associated to an increased sensitivity to the parameter values means that a small perturbation will cause a large undamped reaction of the system, with relaxation oscillations at a frequency $`\omega `$ crossing a minimum while $`\mathrm{\Delta }_0`$ is varied on the fold. In the fold inflection point, where the slope is maximum, the eigenfrequency and damping rate are minimum: the effects of perturbations is expected to be maximum at this point. Note that the evolution of $`\lambda `$ and $`\omega `$ around the inflection point is asymmetric: the damping rate and relaxation frequency increase rapidly for detunings smaller than the inflection point, while they remain of the same order of magnitude on the other side. This asymmetry is expected to have consequences on the dynamics.

Deterministic instabilities are expected to occur for parameters where all stationary solutions are unstable, i.e. in the SF monostable zoneInstDet . In all other areas, the stationary solutions are stable, and therefore, deterministic instabilities cannot occur. However, the presence of the fold and the proximity of a Hopf bifurcation generate a type of stochastic behavior similar to instabilitiesnousprl .

Indeed, on the one hand, the proximity of a Hopf bifurcation is known to favor the appearance of coherence resonanceSRlike , as discussed in nousprl . Let us recall that coherence resonance is the counterpart of stochastic resonance in autonomous systemsSRRevue . Stochastic resonance appears in some forced systems: it may be seen as an amplification by noise of the system response to a forcing. In other terms, the signal to noise ratio of the output periodic signal resulting from the modulation exhibits a maximum when the noise amplitude increases. In stable autonomous systems, the behavior is not periodic, but a phenomenon similar to stochastic resonance can lead to an amplification of an internal resonance: it is the internal stochastic resonance, or coherent resonance CohRes ; CohResPrecursor . Coherent resonance is also known to generate noise induced coherent oscillations, which explain some of our experimental observations, as shown in nousprl .

On the other hand, as shown above, the particular configuration of the eigenvalues on the stable fold makes the system very sensitive to small perturbations, i.e. to noise, and thus we expect a large increasing of the response to noise on the fold. The resulting behavior looks like instabilities, whereas it is simply noise amplification. The behavior obtained on the stable fold when noise is taken into account is detailed in the next section.

## 6 The stable fold: stochastic instabilities

In InstDet , it has been shown that the model is able to reproduce the C instabilities, in the monostable SF zone. On the contrary, S instabilities do not appear in this situation. In fact, it has been shown in nousprl  that S instabilities are not instabilities in the usual significance, but the result of the amplification of the system intrinsic noise, and thus it is necessary to add noise in the simulations to observe S instabilities. This dynamics could appear when the stationary solution of the MOT is not stable, but would be difficult to analyze, because in this case the resulting dynamics would be the superimposition of the deterministic instability with the stochastic motion, the latter being masked by the former. For this reason, the present section is devoted to the study of the influence of noise on the stable stationary solutions of the model, in particular in the vicinity of the fold. As the model here is slightly different from that used in nousprl , we could expect different results from those obtained in nousprl . It is shown in the following that the modifications introduced in the model do not alter the previous conclusions, i.e. that S instabilities are produced by noise amplification on the fold.

To evaluate the response of the system to noise in the vicinity of the fold, we have plotted the amplitude of the motion perturbed by noise, versus the detuning, across the fold (Fig. 12). The amplitude of the motion is measured through the standard deviation $`\sigma `$ of $`Z`$. Fig. 12 has been obtained for a gaussian white noise applied on $`I_1`$ with an amplitude of $`10^3`$, in a situation where the stationary solutions are stable on the fold. Noise amplification appears clearly on the fold, and mimics instabilities, as in the experiments (Fig. 2). The maximum of the motion amplitude, obtained for $`\mathrm{\Delta }_0=0.23`$, corresponds to a standard deviation of $`1.1\times 10^3`$, i.e. $`33`$ $`\mu `$m, in good agreement with the experimental value. By comparison, the motion amplitude is $`0.036`$ for $`\mathrm{\Delta }_0=3`$, i.e. the effect of noise is $`30`$ times larger on the fold than in $`\mathrm{\Delta }_0=3`$. Note that, as similar motion amplitudes are obtained in the experiment and in the model for a noise level $`\zeta =10^3`$, this gives us in return an evaluation of the experimental noise level.

Thus, as for the model developed in nousprl , the fold plays the role of a noise amplifier, leading to an unstable cloud over a limited range of the parameters. The maximum of the motion amplitude appears at the fold inflection point. The amplitude decreases progressively and quasi symmetrically on each side of this point, in spite of the asymmetry of the eigenvalues discussed in the previous section. The range where the instabilities can be observed, appears to be of the order of the “width” of the fold, that we could define as the interval between the two abrupt slope changes delimiting the fold (a more precise definition could be given, but is not useful for the following). In the case discussed here, the range in detuning is 0.02, i.e. much smaller than the experimental range, which was typically of the order of 1. However, the width of the slope is very dependent on the MOT parameters: for example, on Fig. 10, where $`I_1`$ and $`n_0`$ are different, the fold is more than twice larger than in the present case. Thus there is no doubt that it is possible to find a parameter set giving a correct fold width. However, because of the extreme simplicity of the model as compared to the experiments, such a search has no meaning: the aim here is just to show that very few ingredients, including the shadow effect and the noise, are able to reproduce the global experimental behavior. We have now demonstrated that the noise response can effectively have the appearance of instabilities for some parameters. We will now examine in detail the resulting dynamics, and compare it to the experimental observations.

Fig. 13 shows the time evolution of $`Z`$ and $`N`$ in $`\mathrm{\Delta }_0=0.23`$, where the amplification is maximum. The signal is of course stochastic, but several differences appear, on the one hand between the cloud response and the applied noise, and on the other hand between the $`Z`$ and $`N`$ evolutions. Indeed, while the applied noise is white, a frequency dominates in the response, in particular for $`Z`$. This is clearly visible on the signal, and of course on the spectrum (Fig. 14), where a peak appears at a frequency $`\omega _S38`$ Hz. For $`N`$, while the frequency $`\omega _S`$ still appears, low frequency fluctuations dominates, as it is also shown by the spectrum (Fig. 14). This behavior with two different characteristic times is very similar to the experimental SH behavior illustrated in the Fig. 6: in both cases, the higher frequency component dominates the motion of the atoms, while the number of atoms in the cloud is mainly driven by the low frequency component.

In the experiments, the origin of the higher frequency was not identified, as it does not correspond to a known experimental characteristic frequency. However, we know now that the stable stationary solution on the fold has complex eigenvalues, and thus is associated with a relaxation frequency. As noise is known to be able to excite such non linear resonance eigenfrequenciescelet , it is interesting to compare in the simulations the frequency appearing in the dynamics with the eigenfrequency of the stationary solution. This is illustrated on Fig. 15 for different values of the detuning. The full and dashed lines concern the stationary solutions: they shows the evolution of $`\omega `$ and $`Z_S`$ versus $`\mathrm{\Delta }_0`$. They are a reproduction of Fig. 11, and are recalled just for comparison. The squares give the resonance frequencies of the noisy dynamics. They are obtained by fitting the calculated spectra of $`Z`$, as illustrated on Fig. 14a, with a lorentzian function. Similar fits on the $`N`$ spectra give the same results. Fig. 15 shows that for detunings larger than the fold inflection point, the noisy resonance frequency and the eigenfrequency correspond exactly. This confirms that the frequency appearing in the dynamics is the relaxation frequency appearing in the eigenvalues associated with the stationary solution. Thus “instabilities” appear in fact as a noisy amplification of the relaxation frequency of the MOT, through a phenomenon already observed in many other systems, as e.g. in lasers celet : the small damping rate allows noise to excite the relaxation frequency, altering the frequency distribution of the system response to noise.

For detunings smaller than the fold inflection point, Fig. 15 exhibits differences between the eigenfrequencies and the resonance frequencies. A close look at the dynamics of the cloud shows that the main point is not this difference, but a dramatic broadening of the noise resonance, leading to an indetermination of the resonant frequency. This is illustrated on Fig. 15, where the bars associated with each point have a height of $`\mathrm{\Delta }\nu `$, the width of the lorentzian peak obtained by fitting the $`Z`$ spectrum. On the inflection point, the resonance is narrow, as it appears on Fig. 14. On the right of the inflection point, the resonance remains narrow, reaching a maximum of $`\mathrm{\Delta }\nu =14`$ Hz in $`\mathrm{\Delta }_0=0.2`$. On the contrary, for detunings smaller than the inflection point, $`\mathrm{\Delta }\nu `$ increases rapidly, reaching already the value $`\mathrm{\Delta }\nu =30`$ Hz in $`\mathrm{\Delta }_0=0.235`$ (the inflection point is in $`\mathrm{\Delta }_0=0.23`$). It is clear on Fig. 15 that the discrepancy between the noise resonance and the relaxation frequencies is linked to the width of the resonance: as the resonance broadens, it also flattens, and the central frequency becomes irrelevant. This is the reason why the resonant frequencies for detunings smaller than $`\mathrm{\Delta }_0=0.26`$ are not reported on the figure.

This asymmetry in the dynamics around the fold inflection point is of course linked to the asymmetry in the eigenvalues associated with the stationary solution. For detunings larger than the inflection point, the damping rate remains small (less than 45 s-1), and thus the noise excitation of the relaxation oscillations remains efficient. On the contrary, for detunings smaller than the inflection point, the damping rate increases rapidly, decreasing the efficiency of the excitation, and leading to a flat resonance. The consequence on the dynamics is illustrated on Fig. 16, where the power spectra of the $`Z`$ dynamics are represented for two different values of the detuning, located symmetrically with respect to the inflection point. In (a), for $`\mathrm{\Delta }_0=0.22`$, i.e. for a larger detuning than the inflection point, the resonance, centered on 50 Hz, remains narrow ($`\mathrm{\Delta }\nu =11`$ Hz), and thus significant. In (b), for $`\mathrm{\Delta }_0=0.24`$, i.e. on the other side of the inflection point, the resonance, centered on 100 Hz, is already 60 Hz wide, and also five times lower (note the different vertical scales in (a) and (b)). Such a resonance is no more significant from an experimental point of view: indeed, the experimental dynamics is analyzed from time series with a necessary limited number of points, leading to spectra with a resolution much smaller than in the simulations. It is clear that a resonance as those observed on the left of the inflection point could not be detected in the experiments, and thus one expects that in this situation, experiments deliver non resolved spectra. It was effectively the case for SL instabilities (Fig. 4), and thus we can conclude that the SL dynamics observed in the experiments is also explained by the present model.

The results obtained in this section can be summarized in the following way: when the parameters of the cloud are such that the unique stationary solution is stable, and therefore that no deterministic instabilities can occur, the action of noise mimics instabilities. This originates in a fold of the stationary solution, which plays the role of a noise amplifier. Moreover, the particular properties of the eigenvalues associated with the stationary solutions on the fold – small damping rate and relaxation frequency – lead to the existence of the stochastic regime: for detuning larger than the fold inflection point, a resonance is excited by the noise, leading to a dynamics with a dominating frequency; for detuning smaller than the inflection point, the damping rate increases, and the resonance vanishes. This allows us to interpret all the regimes observed experimentally.

The particular properties of the eigenvalues discussed above are in fact linked to the proximity of a Hopf bifurcation. Another known consequence of such a situation is to favor the existence of coherence resonance, and such a phenomenon exists effectively in our case nousprl . Now, an interesting question is to determine what are the exact connections between all these phenomena, and in particular if the existence of coherence resonance has a determining role in that of instabilities. The answer is obviously: no! Indeed, the essential ingredient for noise amplification is the fold; the properties of the eigenvalues only determine the time characteristics of the dynamics. In particular, we have checked that with the present model, no coherence resonance occurs when the MOT exhibits SL instabilities.

## 7 Conclusion

It has been shown recently that the cloud of cold atoms obtained from a magneto-optical trap may exhibit a complex dynamics in the regime of high atomic densities. The observed behaviors can be essentially separated in two different types, depending of their nature: stochastic instabilities nousprl  or deterministic instabilitiesInstDet . The aim of this paper was to describe extensively the experimentally observed stochastic dynamics, and to understand its mechanisms, through a model briefly presented in InstDet  and detailed here.

A detailed analysis of the experimental results in the stochastic regime shows a variety of dynamical behaviors, which differ by the frequency components appearing in the dynamics. Indeed, some instabilities exhibit only low frequency components, while in other cases, a second time scale, corresponding to a higher frequency, appears in the motion of the center of mass of the cloud.

The simple stochastic 1D-model that we use here allows us to retrieve and interpret these experimental dynamics. The model shows that the existence of instabilities is linked to folded stationary solutions where noise response is enhanced. Moreover, the proximity of a Hopf bifurcation and the resulting conditions on the stability of the stationary solutions – small damping rate and existence of a relaxation frequency – explains the existence of several types of regimes: indeed, depending on the parameters, noise is sometimes able to excite the relaxation frequency, leading to the appearance of the second time scale in the dynamics. Globally, the agreement between this over-simplified 1D model and the 3D experiments is surprisingly good. Not only it allows us to make a qualitative interpretation of the experimentally observed dynamics, but also gives quantitative results compatible with the experimental measures, in particular concerning the amplitude of the instabilities. However, it is clear that a 3D model is necessary for an accurate quantitative description of all the experimental observations.

The model also emphasized the close relations between the stochastic and deterministic instabilities. Indeed, both types of dynamics appear to be link to the same factors: the effective dynamics depends mainly on the distance between the working point and the bistable cycles (or the Hopf bifurcation). In fact, the dynamics described in the present paper appears as the first sign of the deterministic instabilities described in InstDet , and thus may be considered as a noisy precursor to deterministic instabilitiesparamaccess ; precursor . A study of the phenomenon from this point of view should put in evidence new properties of these instabilities.

## 8 Acknowledgments

The author thanks M. Fauquembergue and A. di Stefano for their participation in the elaboration of the model, D. Wilkowski for its contribution to the first stages of the experiments, and Ph. Verkerk for useful discussions and his comments about this paper. The Laboratoire de Physique des Lasers, Atomes et Molécules is “Unité Mixte de Recherche de l’Université de Lille 1 et du CNRS” (UMR 8523). The Centre d’Etudes et de Recherches Lasers et Applications (CERLA) is supported by the Ministère chargé de la Recherche, the Région Nord-Pas de Calais and the Fonds Européen de Développement Economique des Régions.


format: text

text:
# The Lack of BLR in Low Accretion Rate AGN as Evidence of their Origin in the Accretion Disk

## 1 Introduction

Seyfert 2 galaxies (where only narrow emission lines are visible) are commonly believed to be intrinsically the same as Seyfert 1 galaxies (where both narrow and broad emission lines are visible), the difference being due to orientation. According to the widely accepted Unification Model for AGNs (Antonucci, 1993) type 2 AGNs are seen edge-on, through large columns of circumnuclear obscuring material that prevents the direct view of the nucleus, including the Broad Line Regions (BLRs). This scenario was first proposed by Antonucci & Miller (1985) to explain the presence of polarized broad lines in the archetypical Seyfert 2, NGC 1068, and is now supported other than by spectropolarimetric observations of Hidden Broad Line Regions (HBLRs) in several other sources, also by X–ray observations, which demonstrate that Seyfert 2s have usually absorption columns largely exceeding the Galactic ones.

Despite observations do generally support orientation based unification models for AGNs, exceptions do exist. Only about 50 % of the brightest Seyfert 2s show the presence of HBLRs (Tran, 2001) in their polarized optical spectra, while the remaining half do not. It has now been convincingly shown that the presence or absence of HBLRs in Seyfert 2s depends on the AGN luminosity, with the HBLR sources having on average larger luminosities (Lumsden & Alexander 2001; Gu & Huang 2002; Martocchia & Matt 2002; Tran 2001, 2003). While Lumsden & Alexander (2001) explained this finding still in the framework of an orientation model, Tran (2001) proposed the existence of a population of galactic nuclei whose activity is powered by starburst rather than accretion onto a supermassive black hole and in which therefore the BLRs simply do not exist.

In this paper we present evidence that suggests that the absence or presence of HBLRs is regulated by the ratio between the X-ray luminosity and the Eddington luminosity which, in the accretion–power scenario, is a measure of the rate at which matter accretes onto the central supermassive black hole. Our explanation is based on the model proposed by Nicastro (2000, hereinafter N00), in which the BLRs are formed by accretion disk instabilities occurring in proximity of the critical radius at which the disk changes from gas pressure dominated to radiation pressure dominated. This radius diminishes with decreasing $`\dot{m}`$; for low enough accretion rates (and therefore luminosities), the critical radius becomes smaller than the innermost stable orbit, and BLRs cannot form. Under the Keplerian assumption, the model naturally predicts that AGNs that are accreting close to the lowest possible $`\dot{m}`$ must show the broadest possible emission lines in their optical spectra, either hidden (i.e. in polarimetric light), if the nucleus is obscured, or not, if the nucleus is not obscured. An analogous model has been proposed recently by Laor (2003). In both Nicastro and Laor’s models the existence of BLRs in AGN is related to the breadth of the broad emission lines (BELs), and based on the observed upper limit of $`FWHM25000`$ for the BELs. However, while in Nicastro’s (2000) model this a consequence of the lines been produced in clouds of gas located at a distance from the source depending on the accretion rate in Eddington units (so the accretion rate is the physical driver of the observed correlation), in Laor’s model the driving parameter is the width of the line itself, and no physical origin is proposed (other than the observed correlation between AGN luminosity and line width).

To check this hypothesis and test the above models, we extracted from the Tran (2001) spectropolarimetric sample, the subsample of type 2 AGNs that also have good quality X-ray spectra available (see Martocchia & Matt 2002). Most of the sources in our sample have \[O iii\] luminosities much lower than $`10^{43}`$ erg s-1. At such relatively low luminosities, stellar light from the galaxy may be a strong contaminant in the optical band. X–ray luminosities, instead, are little or not contaminated by stellar components and so are possibly more reliable indicators of the nuclear, accretion-powered activity. For each source of our sample, we then estimate its 2-10 keV intrinsic (i.e. unabsorbed) luminosity, and use it as a reliable measure of the nuclear activity. We then estimate the mass of the black hole, using the relation between the mass and the bulge luminosity (Magorrian et al. 1998; Ferrarese & Merrit 2000). Comparing the black hole mass and the X-ray luminosity, the accretion rate is eventually derived.

The paper is organized as follows: in Sec. 2 we define and discuss the adopted sample and the methods used to derive the various parameters. In Sec. 3 the results of our analysis are presented, and discussed in Sec. 4.

Throughout this paper a value of $`H_0`$=70 km s-1 Mpc-1 is adopted.

## 2 The sample

We extracted our primary sample from the spectropolarimetric survey of Tran (2001; see also Tran 2003) of all, optically classified, Seyfert 2 galaxies in the CfA and 12$`\mu `$ samples. Optical polarized spectra of the sources in this survey are rather homogeneous in signal-to-noise (Tran, 2001), so the lack of polarized Broad Lines in the optical spectra of about half of the sources in the sample should not be an artifact of dramatic differences in data quality. First we searched for the X-ray properties of the sample \[namely (a) the intrinsic 2-10 keV flux, and (b) the equivalent Hydrogen column density $`N_H`$\], and selected only those sources that had been observed at least once with imaging X-ray satellites. This allowed us to minimize confusion problems which, at the typical flux level of our sample, may be relevant (Georgantopoulos & Zezas, 2003). In practice, we used data only from ASCA, BeppoSAX, Chandra and XMM–Newton, either taken from literature, when available and reliable, or (re)analyzed by ourselves. Finally, we applied two further selections, discarding (a) all sources suspected to be Compton–thick, and (b) sources known to be strongly variable, in X-rays, on time scales of years (e.g. NGC 7172, Dadina et al. 2001). For Compton–thick sources no direct measurements of the intrinsic luminosity is available, while highly variable sources cannot be assigned univocally a X-ray luminosity. Additional details on our primary sample selection can be found in Martocchia & Matt (2002). We ended up with a small but reliable sample consisting of 10 HBLR and 6 non-HBLR sources. Previous works (Gu & Huang, 2002, Martocchia & Matt, 2002, Tran, 2003) have shown that all HBLR sources have 2-10 keV luminosities larger than $`3\times 10^{42}`$ erg s-1, while all non-HBLR sources have luminosities smaller than this value. No correlation with the column density of the absorber is, instead, found. The X–ray luminosity is clearly correlated with the \[O iii\] luminosity (even if with a large scatter), confirming the reliability of the estimates.

For the 16 sources in our primary sample, we then tried to estimate the nuclear accretion rate, in units of Eddington. To calculate $`\dot{m}`$, we needed to evaluate the mass of the central black hole, which we derived using the Ferrarese & Merritt (2000) empirical relationship between nuclear mass and bulge luminosity (we used the numerical coefficients they derived from their sample A). The bulge luminosity $`L_{bulge}`$ was derived from the empirical correlation between the galaxy’s morphological type (i.e. T-type) and the bulge to total luminosity ratio given by Simien & de Vaucouleurs (1986) 111 The $`M_{BH}`$-$`\sigma `$ correlation is more tight than the $`M_{BH}`$-$`L_{bulge}`$ correlation, as discussed in detail by Ferrarese & Merritt (2000). However, for only 4 sources (three of them HBLR) in our sample we could find measured nuclear velocity dispersions (only 2 more can be added, both HBLR, considering all sources in the Gu & Huang 2002 compilation). Moreover, sometimes the reported values differ significantly, with dramatic effects on the estimate of the mass, even of an order of magnitude. For these reasons we were forced to use the $`M_{BH}`$-$`L_{bulge}`$ relationship, despite its larger scattering. . The (corrected for extinction) total luminosity and T-type were taken from the RC3 catalog (de Vaucouleurs et al. 1991). The complete information was available only for a subset of our primary sample. This further selection reduced then our final primary sample to a total of 10 sources: 6 HBLRs and 4 non-HBLRs.

For completeness, and given the limited size of our final primary sample, we also considered a secondary sample, derived applying the same criteria and rules used to select our primary sample to the Gu & Huang (2002) compilation. This compilation contains all Seyfert 2s with published (between 1995 and 2002) spectropolarimetric information, and therefore is not as homogeneous as the Tran (2001) sample in the quality of the polarized optical spectra. This search added 3 more HBLRs and 2 more non-HBLRs to our final sample.

For all sources in our final sample we calculated the accretion rate in Eddington units (actually the ratio between the nuclear bolometric luminosity $`L`$ and the Eddington luminosity $`L_E`$) assuming a factor 10 correction between the 2-10 keV and bolometric luminosities (see footnote 5 in §3). Results are summarized in Table 1 (for the sources of our secondary sample luminosities are slightly different from those reported by Gu & Huang, 2002, because of the different choice of H0).

## 3 Results

It has been recently shown (Gu & Huang, 2002, Martocchia & Matt, 2002, Tran, 2003) that the intrinsic AGN luminosity of the sources of the Tran (2001) sample can be used to clearly separate the two classes of HBLR and non-HBLR sources. Sources with 2-10 keV luminosities larger than the threshold value of $`L_X^{thres}3\times 10^{42}`$ erg s-1 do show HBLs in their optical polarized spectra, while sources with $`L_X<L_X^{thres}`$ do not (see Figure 1 of Martocchia & Matt, 2002). Here we provide evidence that suggests that this separation is due to differences in nuclear accretion rate, from HBLR to non-HBLR sources.

Figure 1 shows fractional luminosities in units of Eddington luminosities, versus black hole masses for all the 15 sources of our primary (circles) and secondary (squares) samples. Open symbols in this plot represent HBLR sources, while filled symbols represent non-HBLR sources. We first note that a very broad range of accretion rates is spanned by the sources of our sample (more than three orders of magnitude), which are otherwise powered by central black holes with rather homogeneous masses (only a factor of about 15 across the entire sample). Most importantly, Figure 1 clearly shows that HBLR sources are accreting at much faster rates compared to non-HBLR sources. 222We note that this result is independent on the particular value chosen to convert 2-10 keV into bolometric luminosities: a different value would only shift the threshold accretion rate value that separates HBLR from non-HBLR sources. Of course a potentially more serious problem may arise if the sources in our sample have Spectral Energy Distributions that dramatically differ from each other. However, such a random effect would likely destroy rather than create the correlation that we find. .

The threshold value of $`\dot{m}_{thres}10^3`$ divides up HBLR from non-HBLR sources in the $`M_{BH}`$ vs $`\dot{m}`$ plane (dashed vertical line in Fig. 1). The only exceptions are NGC 3081 and NGC 3281, both sources from our “secondary” sample. NGC 3081 is a HBLR source, with an accretion rate of $`6\times 10^4`$, a factor of $`2`$ below the threshold value, while NGC 3281 is a non-HBLR source accreting at a rate of $`2\times 10^3`$, again a factor of 2, but above the threshold. We looked at these two exceptions in somewhat higher detail. The flux of NGC 3081 reported in Gu & Huang (2002) is taken from Maiolino et al. (1998; BeppoSAX data). We reanalyzed those data and found very similar results. However, the spectral parameters are very poorly constrained. If the power law index is fixed to 2 (instead of the value of 1.7 chosen by Maiolino et al. 1998), an almost equally good fit is found, but the inferred luminosity (and then $`\dot{m}`$) is twice as large 333This large uncertainty is due to the low S/N of the observation coupled with the large value of the $`N_H`$. , moving the source to the right side of our $`\dot{m}`$ vs $`M_{BH}`$ plane, and so supporting our main conclusion.

The X–ray luminosity reported by Gu & Huang (2002) for the second exception in our sample, the non-HBLR source NGC 3281, is based on an ASCA observation (Bassani et al. 1999), and is underestimated by a factor $``$3. A subsequent BeppoSAX observation, in fact, demonstrated that the source is actually moderately Compton–thick (Vignali & Comastri, 2002), which implies a higher intrinsic luminosity. This, in turn, implies a higher accretion rate compared to that reported in Table 1, and then amplify, rather than moderate, the discrepancy found. However, for this source we could not find in the literature details on the spectropolarimetric observation, so it is impossible to judge how significant is the upper limit on the presence of polarized broad lines.

## 4 Discussion

We have presented evidence for a correlation between the presence of HBLs in the polarized optical spectra of nearby Seyfert 2s, and their nuclear accretion rate. Virtually all sources with HBLRs are found above the threshold value of $`\dot{m}=10^3`$, in the $`\dot{m}`$ vs $`M_{BH}`$ plane, and viceversa. Our sample is admittedly small, but both the spectropolarimetric observations (at least for the sources of our primary sample) and the X–ray luminosities are quite reliable. It should be also noted that the many sources of uncertainties in our $`\dot{m}`$ estimates \[namely (a) the bulge-to-total luminosity ratio vs T type, (b) the black hole mass to bulge luminosity correlation, and (c) the bolometric correction\], would likely destroy, certainly not artificially create, the observed correlation.

Our findings fit nicely with predictions of the model of N00, in which the BLRs originate from the accretion disk at the transition radius between the gas pressure and radiation pressure dominated regions. This suggests that the accretion rate (N00), rather than the “line-width” (Laor, 2003) is indeed the parameter physically responsible for the presence or not of BLRs in AGN.

It is remarkable that the threshold value of $`\dot{m}_{thres}10^3`$ that we find is so close to what predicted in the framework of the N00 model (i.e. about 4$`\times 10^3`$ for a black hole mass of 108M). The threshold value $`\dot{m}_{thres}`$, however, is almost independent of the spin of the black hole, and so, unfortunately, this model cannot be used to discriminate between spinning and non-spinning black holes. This is because in the N00 model the critical radius at which the BLR should form decreases with increasing radiative efficiency in the disk and this decrease is almost exactly balanced by the decrease of the radius of the last stable orbit.

Of course, if our explanation is correct, a fraction (set by the relative number of type 1 to type 2 Seyferts in the nearby Universe) of the non-HBLRs should be actually not obscured. In the Gu & Huang (2002) compilation, 4 out of the 23 non-HBLR sources with X-ray information available (i.e. 17 %) show little ($`N_H<10^{21}`$ cm-2) cold X-ray absorption. In the Martocchia & Matt (2002) sample (i.e. our primary sample, before the bulge luminosity estimate selection) 1 out of the 6 non-HBLR sources (again, 17 %) has an estimated column density only a factor of about 3 larger than the Galactic column along that line of sight. In our final sample only one object (NGC 7590, extracted from the Gu & Huang, 2002, compilation) is virtually unabsorbed. Our estimate of the accretion rate for this objects could then in principle be affected by the AGN light in this source, which would artificially increase the total galaxy luminosity (the observed parameter we use), and so the bulge one. However, the 2-10 keV luminosity of this source is very low (see table 1) and so presumably is its B band luminosity. The effect on the estimate of $`\dot{m}`$, in this case, should therefore be small.

Though suggestive of the existence of a population of unobscured (i.e. type 1) non-HBLR sources, the current statistics is too poor to draw any firm conclusion. Moreover current estimates of equivalent H column densities are based on X-ray observations taken with satellites with poor spatial resolutions (mostly ROSAT, ASCA and BeppoSAX). At such low luminosities, a single spatially unresolved luminous X-ray binary, or a population of those, could lead to underestimate the column of matter actually absorbing the nuclear light (Georgantopoulos, & Zezas, 2003). Better (i.e. higher spatial resolution and signal-to-noise) and more numerous X-ray data are required to disentangle this issue. Enlarging the sample of low luminosity supposedly type 2 AGNs with both good quality optical spectropolarimetric and X-ray (i.e. Chandra and/or XMM–Newton) data, would allow us to eventually find a sub-population of non-HBLR sources, with none or little cold X-ray absorption. These would be the genuine type 1 counterparts of the low accretion rate non-HBLR “true” type 2 sources of our current sample. This population should fall on the left side ($`\dot{m}<\dot{m}_{thres}`$) of the $`\dot{m}`$ vs $`M_{BH}`$ plane (see also Laor, 2003).

We note that, recently, few examples of this class of low-accretion rate type 1 AGNs may have actually been discovered in X-rays (i.e. Pappa et al., 2001, Georgantopoulos, & Zezas, 2003, Boller et al., 2002). The ASCA and Chandra spectra of the optically classified Seyfert 2 NGC 4698, shows little or no absorption (Pappa et al., 2001, Georgantopoulos, & Zezas, 2003). The estimated accretion rate for the central AGN in this galaxy is lower than $`\dot{m}\genfrac{}{}{0pt}{}{_<}{^{}}10^4`$ (Georgantopoulos, & Zezas, 2003), consistent with our proposed scenario. Similarly ROSAT and Chandra spectra of 1ES 1927+654 shows an unobscured, highly variable steep power law continuum, typical of Narrow Line Seyfert 1s (Boller et al., 2002). This object has been recently re-classified as a Seyfert 2 based on its optical spectrum (Bauer et al., 2000). 1ES 1927+654 has an optical (B-band) luminosity of $`10^{43}`$ erg s$`1`$, but no mass estimate is available for its central black hole, and so for the accretion rate in Eddington units. However, given the overall X-ray properties of this object (i.e. steep spectrum and large amplitude variability), we speculate that 1ES 1927+654 is actually accreting at very high rates, in terms of Eddington, and belongs then to the opposite extreme of the line-width vs accretion rate correlation in the framework of the N00 model. If this is the case, 1ES 1927+654 would actually be an optically mis-identified Narrow Line Seyfert 1. The two objects described above have no spectropolarimetric observation, so we do not know whether they have a HBLR or not.

Our findings are based on X-ray luminosities and black hole masses derived from the mass-to-bulge luminosity relation. It is interesting to note that Gu & Huang (2002), also find a similar accretion-rate versus HBLR-existence correlation (see their Fig. 6), using \[O iii\] (rather than X-ray) luminosities and black hole masses obtained from the MBH-$`\sigma `$ correlation.

The results presented here support the N00 model, and suggest therefore that the accretion rate is the physical driver of the observed bimodal distribution of HBLR and non-HBLR sources with luminosities. While other explanations (see e.g. Lumsden & Alexander 2001; Tran 2001, 2003, Martocchia & Matt, 2002, Gu & Huang, 2002) cannot at present be ruled out, the results presented here are very encouraging, and worth checking more throughly by enlarging the sample with new spectropolarimetric and/or X–ray observations.

Part of this work has been done during a visit of GM at CfA, whose hospitality he gratefully acknowledges. We thanks Andreas Zezas for enlightening discussions, and an anonymous referee for useful comments that helped improved the paper. We acknowledge financial support from Chandra grant GO2-3122A (FN), MIUR under grant cofin-00-02-36 (AM and GM), and CNES (AM).


text:
# The Relation between Black Hole Mass, Bulge Mass, and Near-Infrared Luminosity

## 1 Introduction

Central massive black holes ($`\mathrm{BH}`$s) are now thought to reside in virtually all galaxies with a hot spheroidal stellar component (hereafter bulge). Such $`\mathrm{BH}`$s seem to be a relic of past quasar activity (e.g., Sołtan 1982; Marconi & Salvati 2002; Yu & Tremaine 2002; Aller & Richstone 2002) and related to host galaxy properties, with the implication that $`\mathrm{BH}`$ and galaxy formation processes are closely linked. Previous work has shown that $`\mathrm{BH}`$ mass $`M_{\mathrm{BH}}`$ is correlated with both blue luminosity $`L_{\mathrm{B},\mathrm{bul}}`$ and bulge mass $`M_{\mathrm{bul}}`$, although with considerable intrinsic scatter ($`rms0.5`$ in log$`M_{\mathrm{BH}}`$; Kormendy & Richstone 1995). However, $`M_{\mathrm{BH}}`$ and the bulge effective stellar velocity dispersion $`\sigma _e`$ correlate more tightly ($`rms0.3`$) than $`M_{\mathrm{BH}}`$-$`L_{\mathrm{B},\mathrm{bul}}`$ (Ferrarese & Merritt, 2000; Gebhardt et al., 2000). The smaller scatter of the $`M_{\mathrm{BH}}`$-$`\sigma _e`$ correlation suggests that the bulge dynamics (or mass), rather than luminosity, is the agent of the correlation. But the smaller spread relative to $`M_{\mathrm{BH}}`$-$`L_{\mathrm{bul}}`$ appears to be an artefact of the manipulations necessary to derive $`L_{\mathrm{bul}}`$. Indeed, recent work has shown that when bulge parameters are measured with more accuracy \[e.g. profile fitting rather than average correction for disk light (Simien & de Vaucouleurs, 1986)\], the resulting scatter is comparable to that of $`M_{\mathrm{BH}}`$-$`\sigma _e`$ (McLure & Dunlop, 2002; Erwin et al., 2003). The correlation between $`M_{\mathrm{BH}}`$ and bulge light concentration also has a comparably low scatter (Graham et al., 2001). Nevertheless, there are strong indications that $`L_{\mathrm{B},\mathrm{bul}}`$ of the brightest elliptical galaxies, for which decomposition issues are unimportant, deviate significantly from the $`M_{\mathrm{BH}}`$-$`L_{\mathrm{bul}}`$ relation (Ferrarese, 2002). Hence, longer wavelengths may also be necessary to better define the intrinsic scatter in $`M_{\mathrm{BH}}`$-$`L_{\mathrm{bul}}`$ compared to that of $`M_{\mathrm{BH}}`$-$`\sigma _e`$.

In this paper, we reexamine the $`M_{\mathrm{BH}}`$-$`L_{\mathrm{bul}}`$ correlation by accurately measuring the bulge luminosity in the near-infrared (NIR) for all galaxies with a well-determined $`M_{\mathrm{BH}}`$. All previous studies have used optical light (B or R) to test the $`M_{\mathrm{BH}}`$-$`L_{\mathrm{bul}}`$ relation, but NIR light provides a clear advantage over the optical: it is a better tracer of stellar mass and less subject to the effects of extinction. If the physical correlation is between the $`\mathrm{BH}`$ mass and bulge mass, the NIR correlations $`M_{\mathrm{BH}}`$-$`L_{\mathrm{bul}}`$ should be tighter than those in the optical, because of the smaller variation of $`M/L`$ ratio $`\mathrm{{\rm Y}}`$ with mass (e.g., Gavazzi 1993). Moreover, we use a two-dimensional (2D) bulge/disk decomposition to determine bulge parameters, an improvement on earlier work which applied 1D fits only. Here we construct the largest possible sample, by considering all galaxies which have been used for the $`M_{\mathrm{BH}}`$-$`\sigma _e`$ and $`M_{\mathrm{BH}}`$-$`L_{\mathrm{B},\mathrm{bul}}`$ correlations. In §2 we present the sample of galaxies with direct dynamical $`\mathrm{BH}`$ mass measurements, and in §3 describe the images and the 2D bulge/disk decomposition applied to them. Finally, in §4 we discuss the results of the analysis.

## 2  The Sample

To date, there are 37 galaxies with direct gas kinematical or stellar dynamical determination of the central $`\mathrm{BH}`$ mass. These galaxies have been compiled and made into a uniform sample (e.g., for distances) by a number of authors (e.g., Merritt & Ferrarese 2002; Tremaine et al. 2002; hereafter MF02 and T02, respectively). We adopt the data from the recent paper by T02 with some modifications and additions. The data in Columns 1-5 and 9 of Table 1 are from the compilation by T02 and the reader can refer to that paper for more details. Differently from T02, when galaxy distances from surface brightness fluctuations (Tonry et al. 2000) are not available, we use recession velocities corrected for Virgocentric infall from the LEDA database (http://leda.univ-lyon1.fr/) with $`H_0=70\mathrm{km}\mathrm{s}^1\mathrm{Mpc}^1`$. In a few cases, we also consider $`\mathrm{BH}`$ mass estimates from different papers than those used by T02; thus, in Col. 6, we indicate the appropriate references. With respect to the 31 galaxies considered by T02 we add: Cygnus A (Tadhunter et al., 2003), M81 (Devereux et al., 2003), M84 (Bower et al., 1998), NGC 4594 (Kormendy, 1988), Centaurus A (Marconi et al., 2001) and NGC 5252 (Capetti et al., 2003).

Following MF02, we divide the galaxies into two groups. In the first group, we place all the galaxies which have a secure $`\mathrm{BH}`$ mass measurement and an accurate determination of the bulge NIR luminosity. We consider ‘secure’ those $`\mathrm{BH}`$ masses for which the black hole sphere of influence, $`R_{\mathrm{BH}}=GM_{\mathrm{BH}}/\sigma _e^2`$ (column 7 of Table 1), has been clearly resolved, i.e., $`N_{res}=2R_{\mathrm{BH}}/R_{res}>1`$, where $`R_{res}`$ is the spatial resolution of the observations. Additional reasons for placing galaxies in Group 2 are given in Table 1.

## 3  Image Analysis

We have constructed a homogeneous set of NIR images of the galaxies presented in Table 1 (except for the Milky Way and M31) by retrieving J, H, and K atlas images from the 2-Micron All Sky Survey (2MASS; http://www.ipac.caltech.edu/2mass). When a single atlas image contained only a portion of the galaxy, we also retrieved adjacent tiles and mosaicked the images after subtracting the sky background and rescaling for the different zero points. The 2MASS images are photometrically calibrated with a typical accuracy of a few percent. More details can be found in Hunt & Marconi (2003; hereafter Paper II).

We performed a 2D bulge/disk decomposition of the images using the program GALFIT (Peng et al., 2002) which is made publicly available by the authors. This code allows the fitting of several components with different functional shapes (e.g., generalized exponential (Sersic) and simple exponential laws); the best fit parameters are determined by minimizing $`\chi ^2`$. More details on GALFIT can be found in Peng et al. (2002). We fit separately the J, H and K images. Each fit was started by fitting a single Sersic component and constant background. When necessary (e.g., for spiral galaxies), an additional component (usually an exponential disk) was added. In many cases these initial fits left large residuals and we thus increased the number of components (see also Peng et al. 2002). The fits are described in detail in Paper II. In Table 1 we present the J, H and K bulge magnitudes, effective bulge radii $`R_e`$ in the J band, and their uncertainties. The J, H and K magnitudes were corrected for Galactic extinction using the data by Schlegel, Finkbeiner, & Davis (1998). We used the J band to determine $`R_e`$ because the images tend to be flatter, and thus the background is better determined.

## 4  Results and Discussion

In Fig. 1 we plot, from left to right, $`M_{\mathrm{BH}}`$ vs $`L_{\mathrm{K},\mathrm{bul}}`$, $`M_{\mathrm{BH}}`$ vs $`M_{\mathrm{bul}}`$, and the residuals of $`M_{\mathrm{BH}}`$-$`\sigma _e`$ vs $`R_e`$ (based on the fit from T02). Only Group 1 galaxies are shown. $`M_{\mathrm{bul}}`$ is the virial bulge mass given by $`kR_e\sigma _e^2/G`$; if bulges behave as isothermal spheres, $`k=8/3`$. However, comparing our virial estimates $`M_{\mathrm{bul}}`$ with those $`M_{\mathrm{dyn}}`$, obtained from dynamical modeling (Magorrian et al., 1998; Gebhardt et al., 2003), shows that $`M_{\mathrm{bul}}`$ and $`M_{\mathrm{dyn}}`$ are well correlated ($`r=0.88`$); setting $`k=3`$ (rather than $`8/3`$) gives an average ratio of unity. Therefore, we have used $`k=3`$ in the above formula. Considering the uncertainties of both mass estimates, the scatter of the ratio $`M_{\mathrm{bul}}`$/$`M_{\mathrm{dyn}}`$ is 0.21 dex. We fit the data with the bisector linear regression from Akritas & Bershady (1996) which allows for uncertainties on both variables and intrinsic dispersion. The FITEXY routine (Press et al., 1992) used by T02 gives consistent results (see Fig. 1). Fit results of $`M_{\mathrm{BH}}`$ vs galaxy properties for Group 1 and the combined samples are summarized in Table 2. The intrinsic dispersion of the residuals ($`rms`$) has been estimated with a maximum likelihood method assuming normally-distributed values. Inspection of Fig. 1 and Table 2 show that $`L_{\mathrm{K},\mathrm{bul}}`$ and $`M_{\mathrm{bul}}`$ correlate well with the $`\mathrm{BH}`$ mass. The correlation between $`M_{\mathrm{BH}}`$ and $`M_{\mathrm{bul}}`$ is equivalent to that between the radius of the $`\mathrm{BH}`$ sphere of influence $`R_{\mathrm{BH}}`$ ($`=GM_{\mathrm{BH}}/\sigma _e^2`$) and $`R_e`$.

### 4.1 Intrinsic Dispersion of the Correlations

To compare the scatter of $`M_{\mathrm{BH}}`$-$`L_{\mathrm{bul}}`$ for different wavebands, we have also analyzed the B-band bulge luminosities for our sample. The upper limit of the intrinsic dispersion of the $`M_{\mathrm{BH}}`$-$`L_{\mathrm{bul}}`$ correlations goes from $`0.5`$ dex in $`\mathrm{log}M_{\mathrm{BH}}`$ when considering all galaxies, to $`0.3`$ when considering only those of Group 1. Hence, for galaxies with reliable $`M_{\mathrm{BH}}`$ and $`L_{\mathrm{bul}}`$, the scatter of $`M_{\mathrm{BH}}`$-$`L_{\mathrm{bul}}`$ correlations is $`0.3`$, independently of the spectral band used (B or JHK), comparable to that of $`M_{\mathrm{BH}}`$-$`\sigma _e`$. This scatter would be smaller if measurement errors are underestimated. McLure & Dunlop (2002) and Erwin et al. (2003) reached a similar conclusion using R-band $`L_{\mathrm{bul}}`$ but on smaller samples. The correlation between R-band bulge light concentration and $`M_{\mathrm{BH}}`$ has a comparable scatter (Graham et al., 2001).

Since $`M_{\mathrm{BH}}`$-$`L_{\mathrm{B},\mathrm{bul}}`$ and $`M_{\mathrm{BH}}`$-$`L_{\mathrm{NIR},\mathrm{bul}}`$ have comparable dispersions, the rough bulge/disk decomposition (§1), the larger reddening and stellar population effects do not apparently compromise the correlation. Most of the galaxies in the sample are early types, and thus may be less sensitive to the above effects. However, the scatter in the $`M_{\mathrm{BH}}`$-$`L_{\mathrm{bul}}`$ correlations does not decrease significantly when considering only elliptical galaxies.

The correlation between $`M_{\mathrm{BH}}`$ and $`M_{\mathrm{bul}}`$ has a slightly lower dispersion (0.25 versus 0.3) than $`M_{\mathrm{BH}}`$-$`L_{\mathrm{bul}}`$. If the scatter of $`M_{\mathrm{bul}}`$-$`M_{\mathrm{dyn}}`$ (0.21 dex) is an indication of the additional uncertainties on our virial estimates, then the intrinsic scatter of $`M_{\mathrm{BH}}`$-$`M_{\mathrm{bul}}`$ drops to $`0.15`$ dex. Judging from the present data where only secure $`M_{\mathrm{BH}}`$ and $`L_{\mathrm{bul}}`$ are included, $`\sigma _e`$, $`M_{\mathrm{bul}}`$, $`L_{\mathrm{B},\mathrm{bul}}`$ and $`L_{\mathrm{NIR},\mathrm{bul}}`$ provide equally good $`M_{\mathrm{BH}}`$ estimates to within a factor of $`2`$.

### 4.2 Correlation slopes

All the slopes are roughly unity, but those of $`M_{\mathrm{BH}}`$-$`L_{\mathrm{bul}}`$ are systematically larger than that of $`M_{\mathrm{BH}}`$-$`M_{\mathrm{bul}}`$. This is expected if $`L_{\mathrm{bul}}`$ correlates with $`M_{\mathrm{BH}}`$ because of its dependence on $`M_{\mathrm{bul}}`$ through the stellar $`M/L`$ ratio. From our $`M_{\mathrm{bul}}`$-$`L_{\mathrm{bul}}`$ relation, we find that $`\mathrm{log}\mathrm{{\rm Y}}_K=\mathrm{\hspace{0.17em}0.18}\mathrm{log}L_{\mathrm{K},\mathrm{bul}}2.1`$; the weak dependence of $`\mathrm{{\rm Y}}`$ on $`L`$ fully accounts for the different slopes of $`M_{\mathrm{BH}}`$-$`L_{\mathrm{K},\mathrm{bul}}`$ and $`M_{\mathrm{BH}}`$-$`M_{\mathrm{bul}}`$, and the same applies to the J and H bands.

All correlations are thus consistent with a direct proportionality between $`M_{\mathrm{BH}}`$ and bulge mass. This contrasts with previous claims of a non-linearity of the $`M_{\mathrm{BH}}`$-$`M_{\mathrm{bul}}`$ relation (Laor, 2001) but is in agreement with McLure & Dunlop (2002). A partial correlation analysis of log$`M_{\mathrm{BH}}`$ (variable $`x_1`$), log$`\sigma _e`$ ($`x_2`$), and log$`R_e`$ ($`x_3`$) shows that $`M_{\mathrm{BH}}`$ is separately significantly correlated both with $`\sigma _e`$ and $`R_e`$. The Pearson partial correlation coefficients, in which the known dependence of $`\sigma _e`$ and $`R_e`$ is eliminated, are $`r_{12}=0.83`$, $`r_{13}=0.65`$, with a significance of $`>99.9\%`$. This is shown graphically in Fig. 1c, where the residuals of the T02 $`M_{\mathrm{BH}}`$-$`\sigma _e`$ correlation are plotted against $`R_e`$; there is a weak, but significant, correlation of these residuals with $`R_e`$. Consequently, when galaxy structural parameters are measured carefully from 2D image analysis, the additional, weaker, dependence of $`M_{\mathrm{BH}}`$ on $`R_e`$ is uncovered. Thus, a combination of both $`\sigma _e`$ and $`R_e`$ is necessary to drive the correlations between $`M_{\mathrm{BH}}`$ and other bulge properties. This fundamental plane of black holes will be further investigated elsewhere.

The average $`\mathrm{log}M_{\mathrm{BH}}/M_{\mathrm{bul}}`$ can be estimated assuming a log-normal distribution with normally distributed observational errors. With maximum likelihood we find $`\mathrm{log}M_{\mathrm{BH}}/M_{\mathrm{bul}}=2.63`$ with an intrinsic dispersion of 0.27 dex (-2.79 and 0.49 dex for all galaxies). Adopting the method of Merritt & Ferrarese (2001), we find $`\mathrm{log}M_{\mathrm{BH}}/M_{\mathrm{bul}}=2.81`$ and $`rms=0.36`$ (-2.86 and 0.44 for all galaxies) consistently with their result of -2.9 and 0.45 dex (see also McLure & Dunlop 2002).

We thank A. Capetti, E. Emsellem, W. Maciejewski, C. Norman, and E. Oliva for useful discussions. This paper was partially supported by ASI (I/R/112/01 and I/R/048/02) and MIUR (Cofin01-02-02). This publication uses the LEDA database, the NASA/IPAC Extragalactic Database operated by JPL, CalTech, under contract with NASA, and data products from 2MASS, a joint project of the University of Massachusetts and the IPAC/CalTech, funded by NASA and NSF.