\section{Introduction}

Red giants are cool  stars with an extended convective envelope, which can, as in main sequence solar-like stars,  stochastically excite pressure modes
of oscillation.  Although  stochastic  oscillations were already  detected in a few red giants from ground and space observations \citep[e.g.][]{Frandsen02, Joris06, Barban07} it has been only thanks to the photometric  space mission COROT \citep{Baglin02} that an unambiguous detection of radial and non-radial modes in a large number of red-giant stars has been achieved \citep{Joris09, Hekker09, Carrier10}. That confirmation has opened the way to the seismic study of the structure and evolution of these objects that play a fundamental role in fields such as stellar age determination and chemical evolution of galaxies.
The application of \cite{KB95} theoretical scaling laws, which relate basic seismic observables (the large frequency separation -- $\Delta\nu$--  and the frequency at maximum power --$\nu_{\rm max}$) to the stellar global parameters,  allowed \cite{Mosser10}  to estimate the masses of CoRoT and KEPLER red giants (adopting a value for the effective temperature). Moreover, the combination of these scalings  with the predictions of  population synthesis models let \cite{MiglioPop09} and \cite{MiglioPop209} to  characterize the population of  CoRoT and Kepler targets.

These results, showing that a vast amount of information can be extracted from quite   easy-to-access
seismic observables have deeply changed the perception of the predictive capabilities of asteroseismology and have strengthened its interaction with the other fields of astrophysics. However, much more information is contained in the oscillation spectra of these large number of red giants. In this Letter, we analyze  the properties of red giant adiabatic oscillation spectra and relate them with their evolutionary state.

\begin{figure*}[ht!]
\centering
\resizebox{\hsize}{!}{\includegraphics{OC5_2_fig2}\includegraphics{OC5_2_fig3}}
\resizebox{\hsize}{!}{\includegraphics{OC5_2_fig4}\includegraphics{OC5_2_fig5}}
\vspace*{-1.75cm}
\caption{Upper panels: propagation diagrams for two models at almost the same luminosity,  1.5~\msun\  in the RGB (left) and 2.5\msun\  in central He-burning phase (right). Horizontal dotted lines limit the solar-like frequency domain for each model. Solid line is \BV\ frequency, and short- and long-dashed lines correspond to the Lamb frequency for $\ell=1$ and 2 respectively.  Lower panels: corresponding plots of  inertia as a function of frequency for $\ell$=0 (black circles), 1 (grey triangles) and 2 (dark-grey squares) modes in the solar-like oscillation domain (dashed vertical lines).}
\label{label_prop}
\end{figure*}

\section{Stellar Models}
\label{sec_models}
Stellar models were computed with the code ATON3.1 \citep[][and references therein]{Ventura08} for masses between 0.7 and 5.0 \msun\, and  different chemical compositions:  He mass fraction, $Y$=0.25 and 0.278,  and metal mass fraction, $Z$=0.006, 0.01, 0.015, 0.02 and 0.03. The energy transport in the convective regions was modeled with the classic mixing length treatment with  $\alpha_{\rm MLT}=1.6$. For a given chemical composition ($Z=0.02$, $Y=0.278$), models with  $\alpha_{\rm MLT}=1.9$ and  FST treatment of convection \citep{Canutoetal96} were also computed. The evolution of these models was followed from the pre-main sequence until the exhaustion of He in the core for models more massive than 2.3\msun, and until the helium flash for the less massive ones. The core He-burning phase for low-mass ($0.7-2.3$\msun) stars (Red clump stars) has also been followed starting from  zero age horizontal branch  models. Microscopic diffusion was not included but its effects on red giant models \citep{michaud10} are not relevant for the present study.

The evolution rate during  ascending RGB, descending  RGB and core He-burning phases is very different and strongly  depends on stellar mass. For  low-mass stars the  time spent ascending the RGB  may be comparable with that of  core He-burning of more massive stars. As a consequence, observing stars in both evolutionary stages would be equally likely. Concerning the internal structure of these  models, it is worth mentioning that  for a low-mass model (1.5~\msun, for instance) the density contrast ($\rho_c/\langle\rho\rangle$, central to mean density ratio) changes from $10^6$  at the bottom of its RGB to  $3\times 10^9$  at $\log L/L_\odot\sim 2$. Models in the core He-burning phase, on the other hand,  have a much lower value of  $\rho_c/\langle\rho\rangle$ of the order of $2\times 10^7$; moreover,  due to  the high dependence on temperature of the 3$\alpha$ nuclear reactions, they develop  a small convective core. At a given  luminosity, e.g. that of the Red Clump,  the value of $\rho_c/\langle\rho\rangle$ for a 1.5\msun\ RGB model is more than 10 times that of  a He-burning one. So different structures should imply significant effects on the oscillation properties.

\section {Adiabatic oscillation properties}

Adiabatic oscillation frequencies were computed with an Eulerian version of the code LOSC \citep{Scuflaire08} for models from the bottom of the RGB until a maximum luminosity  ($\log L/L_\odot ~ 2.2-3.2$, depending on mass) and as well during the phase  of core  He-burning. In this paper we deal with adiabatic computations and do not consider the problem of excitation and damping of solar-like oscillations in red giants \cite[][]{Dziembowski01, HoudekGough02, Dupret09}. We use the scaling laws \citep[][]{Brownetal91, KB95} to derive  the frequency domain in which solar like oscillations are expected, and the value of the mode inertia as an estimate of the expected mode amplitude \citep[][]{JCD04}.  We search oscillation modes with angular degree $\ell=0$, 1, 2 and 3\footnote{$\ell=3$ only for $Z=0.02$, $Y=0.278$, $\alpha_{\rm MLT}=1.9$}  in the domain of frequencies defined by  an interval around $\nu_{\rm max}$ \citep[Eq. (10) of][]{KB95}. The width of the solar-like frequency domain is taken to be 20\% larger than the  difference between the acoustic cutoff frequency in the stellar atmosphere  and $\nu_{\rm max}$.

\begin{figure}[ht!]
\centering
\resizebox{\hsize}{!}{\includegraphics{Sig1_ls0}}
\resizebox{\hsize}{!}{\includegraphics{Sig1_rho}}
\caption{Mean squared deviation with respect to the main value of $\Delta\nu_{\ell=1}$, for RGB and core He-B models between 0.7\msun\ (darkest symbols) and 5.0\msun\ (lightest symbols), as a function of mean large separation for radial modes (a), and as a function of the density contrast (b). Vertical dotted and dashed lines correspond to 1\msun\ models at $\log L/L_{\odot} \sim 1.68$ in the RGB and core He-B phases respectively.}
\label{label_sigma}
\end{figure}

\begin{figure}[ht!]
\centering
\resizebox{\hsize}{!}{\includegraphics{mosser_theo}}
\caption{ Theoretical oscillation pattern expected for $\ell=0$ (light grey), 1 (black) and 2 (dark-grey) modes from the red giant population in the CoRoT-exofield in the galactic center direction. In ordinate the large frequency separation for radial modes, and in abscissa the frequency of ($n,\ell$) modes are normalized by the large separation.}
\label{label_mosser}
\end{figure}

\begin{figure*}[ht!]
\centering
\resizebox{\hsize}{!}{\includegraphics{d02_LS}\includegraphics{d02N_LS}\includegraphics{d02_d01}}
\resizebox{\hsize}{!}{\includegraphics{d03_LS}\includegraphics{d02N_d03N}\includegraphics{d02_d03}}
\caption{Small frequency separations for models with masses between 0.7 and 5~\msun\ (each grey level corresponds to a different mass, from the darkest for 0.7\msun, to the lightest for 5\msun), $Z$=0.006, 0.015, 0.02 and 0.03, $Y$=0.025 and 0.278, $\alpha_{\rm MLT}=1.6$ and 1.9. Dotted vertical line indicates the minimum $\langle\Delta\nu\rangle$ from the first 34-d of KEPLER data, and the dot-dashed lines correspond to the linear fits for these data from  \cite{Beddingetal10}. a)  $\langle\delta\nu_{02}\rangle$ as a function of $\langle\Delta\nu_0\rangle$, for  low-mass models (0.7-2.3,\msun) from the bottom of the RGB up to  $\log L/L_\odot ~ 2.5-3$, and 2.5-5\msun\ models in the He-burning phase ($Y_c$ between 0.9 and 0.1).  b) Same as a) but for the $\langle\delta\nu_{02}\rangle/\langle\Delta\nu_0\rangle$. c) $\langle\delta\nu_{01}\rangle/\langle\Delta\nu_0\rangle$ as a function of  $\langle\delta\nu_{02}\rangle/\langle\Delta\nu_0\rangle$ for the same models as in a) and b). Note the high concentration of RGB models with small and negative value of $\langle\delta\nu_{01}\rangle/\langle\Delta\nu_0\rangle$. d) $\langle\delta\nu_{03}\rangle/\langle\Delta\nu_0\rangle$ {\it vs.} $\langle\Delta\nu_0\rangle$  for models with $Y=0.278$, $Z=0.02$, and $\alpha_{\rm MLT}=1.9$. e) $\langle\delta\nu_{03}\rangle/\langle\Delta\nu_0\rangle$ {\it vs.} $\langle\delta\nu_{02}\rangle/\langle\Delta\nu_0\rangle$ and models used in d). f) $\langle\delta\nu_{03}\rangle$ as a function of $\langle\delta\nu_{02}\rangle$ and models in d).}
\label{label_d02}
\end{figure*}

The  properties of oscillation modes depend on the behaviour of the \BV\ ($N$) and Lamb ($S_{\ell}$) frequencies. In red-giant models, $N$ reaches huge values in the central regions and therefore the frequency of gravity modes (g-modes)  and their number by frequency interval ($n_g$) increase with respect to main sequence models. On the other hand, the low mean density makes  the frequency of pressure modes (p-modes) to decrease. All that leads to an oscillation spectrum for red-giants where, in addition to   radial modes, one finds a large number of non-radial modes with mixed g-p properties. The dominant character of these non-radial modes depends on the separation between gravity and acoustic cavities, and may be estimated from the value of the normalized  mode inertia ($E$) \citep[see e.g.][ and references therein]{JCD04}. Therefore, some non-radial modes may be well trapped in the acoustic cavity and behave as p-modes presenting a  mode inertia close to that of radial modes, while modes with strong mixed g-p character have larger $E$ value.  Hereafter, we will use the term p-modes in quotation marks to refer to mixed modes with a dominant p-character.

 In Fig.~\ref{label_prop} we present, in top panels, the $\ell=1$, 2  propagation diagrams for a RGB 1.5~\msun\ model (left) and for a core He-burning (He-B) model of 2.5\msun\ (right).  In the bottom panels we plot the variation of the mode inertia with frequency for radial and non-radial  modes ($\ell=1,\,2$). As mentioned above, the RGB model is ten times more centrally condensed than the He-B one. The huge difference in  density between the central region and the convective envelope entails a high potential barrier between the acoustic and the gravity cavities reducing the interaction between p and g modes. As a consequence,  we find  for RGB models that $\ell=1$ modes  with $E_{\ell=1}\approx E_{\ell=0}$  are quite regularly spaced in frequency. For He-B ones,  the coupling between these cavities is more important and $\ell=1$ modes are mixed modes with $E_{\ell=1}> E_{\ell=0}$. Nevertheless, $E_{\ell=1}$ presents  still a minimum value for modes between two consecutive radial ones showing a  somewhat regular pattern. Even if the $E$ value is larger than that corresponding to radial modes  we can still consider those modes, based on the value of $E$, as observable ``p-modes''. For $\ell=2$ modes, the coupling  between the g- and p-cavities is smaller  than for $\ell=1$ and hence the trapping more efficient.  Therefore, independently of the central  condensation of the model,  a spectrum of regularly spaced $\ell=2$ ``p-modes''  with $E_{\ell=2}\approx E_{\ell=0}$ is expected. Finally, note that  the turning points for acoustic modes ($tp_\ell$ defined as the point where $\nu_{\rm max}=S_{\ell}$) are inside the convective envelope for the RGB model and in the radiative region for the He-B one.

In the asymptotic theory for p-modes \citep{vandakurov67, tassoul80, gough86} the frequencies of two modes of same degree and consecutive order are separated by a constant value $\langle\Delta \nu\rangle$  which is  approximately  independent of $\ell$ for low-degree modes, and related to the mean density of the star. Of course, the asymptotic theory is no longer valid for mixed modes and in regions with rapidly varying  physical quantities, nevertheless  the modes partially or well trapped in  the  acoustic cavity  with a dominant p-character (``p-modes'') show such regular pattern. We select then,  as ``p-modes'' of degree $\ell$, the modes with  the  minimum inertia  between  two consecutive radial modes and use them to compute the large and small frequency separations and to analyze their dependence on the stellar parameters and evolutionary state.

\subsection{Large Separation: $\Delta\nu$}

The mean value of the large frequency separation ($\langle\Delta\nu_{\ell}\rangle=\langle\nu_{n,\ell}-\nu_{n-1,\ell}\rangle$ averaged over the radial order $n$) decreases as the star ascends the RGB with a denser and denser core and a more and more diffuse envelope, but it is not possible  to distinguish, on the basis of the $\langle\Delta\nu\rangle$  value alone, among different evolutionary states, ascending RGB, descending RGB, or core He-burning.   Nevertheless, an indirect information about the evolutionary state is provided not by the average value $\langle\Delta\nu\rangle$, but by the deviation of $\Delta\nu$ as a function of frequency with respect to its mean value ($\sigma(\Delta\nu_{\ell})$). Radial and  $\ell=2$ ``p-modes'', as  mentioned above, show a very regular pattern, and  the mean quadratic deviation of $\Delta\nu(\nu)$ with respect to its mean value over the solar-like frequency domain ($\sigma(\Delta\nu_{\ell})$) is always smaller than 5\% for all the evolutionary states and masses considered.  On the contrary, $\sigma(\Delta\nu_1)$ strongly depends on the evolutionary state and  while its value remains small for more concentrated models, it may get values as large as 50\% for core He-burning ones.

Fig.~\ref{label_sigma} shows  how the scatter of $\ell=1$ modes depends on $ \Delta\nu$ (that is on $\langle \rho \rangle$) and on the density contrast. The scatter of $\ell=1$ frequencies decreases as the luminosity increases from the bottom of the RGB (low $\rho_c/\langle\rho\rangle$, or highest $\langle\Delta\nu\rangle$).  As the star goes up the RGB,  the $\ell=1$ modes are better trapped in the acoustic cavity and the spectra of dipole modes are more regular. These results are consistent with the  observational results obtained by \cite{Beddingetal10} for the first 34 days of KEPLER observations:  radial and quadruple modes for a sample of 50 low luminosity red giants ( $< 30$\lsun\ from scaling laws) show a low scatter, while dipole modes present a  significantly larger one.
 The large dispersion shown in Fig.~\ref{label_sigma}  around $\langle \Delta\nu_0\rangle \sim 4\,\mu$Hz and $\rho_c/\langle\rho\rangle \sim 2 \times 10^7 $ corresponds to models  burning He in the core. Following the core expansion, the external convective zone recedes, increasing the coupling  between gravity and acoustic cavities (see. Fig.~\ref{label_prop}) and the mixed character of oscillation modes.  The spectrum of ``trapped'' $\ell=1$ modes is then much less regular. It is worth pointing out  that, for a given mass (1~\msun, for instance) at Red-Clump luminosity, while RGB and core He-burning models have the same value of $\langle \Delta\nu_0\rangle$,  $\rho_c/\langle\rho\rangle$ may differ by two orders of magnitude.

The distributions of $ \Delta\nu$ and $\nu_{\rm max}$ for  the CoRoT-exofield red giant sample show  a single dominant peak located at  $\sim 4\mu$Hz and $\sim 30\mu$Hz respectively \citep[][]{Hekker09} which was interpreted by \cite{MiglioPop09} as being consistent with a population of red giants dominated by Red-Clump stars.  According to population synthesis simulations done with  TRILEGAL \citep[][]{girardietal05}  for  the CoRoT field in the galactic center direction, 70\% of the stars at the Red-clump luminosity are low-mass stars in the core He-burning phase, and therefore, with oscillation spectra significantly different from those of the other 30\% of stars that, at the same luminosity, are in the RGB phase.
In Fig.~\ref{label_mosser} we plot the ``ensemble oscillation pattern''  that we expect for the CoRoT red giant population in the galactic center direction. This diagram was obtained by taking the theoretical adiabatic spectra corresponding to the stellar models whose stellar parameters are closest to those resulting from population synthesis computations. Note the large scatter of $\ell=1$ modes in the region around $\langle \Delta\nu \rangle \sim 4\,\mu$Hz, and the regular pattern of $\ell=0$ and 2 modes.  Both characteristics are in good agreement with the observational results found by \cite{Mosser10b}.

\begin{figure}[ht!]
\centering
\resizebox{\hsize}{!}{\includegraphics{d01_rac}}
\resizebox{\hsize}{!}{\includegraphics{d01Z_rac}}
\caption{a) $\langle\delta\nu_{01}\rangle/\langle\Delta\nu_0\rangle$ as a function of the distance in  acoustic radius between the bottom of the convective envelope $\tau_{\rm BCZ}$  and the $\ell=1$ turning point $\tau_{\rm tp1}$  for models with i $Z=0.02$, $Y=0.278$, and $\alpha_{\rm MLT}=1.9$ and masses between 0.7 and 5~\msun.  b) zoom for low mass models (0.7--2.3\msun) in the RGB. Same grey level criteria as in Fig.4.}
\label{label_d01L}
\end{figure}

\subsection{Small separations: $\delta\nu_{02}$, $\delta\nu_{03}$ and $\delta\nu_{01}$}

According to the asymptotic theory, the  mean small frequency separation ($\langle\delta\nu_{02}\rangle=\langle\nu_{n0}-\nu_{n-1\,2}\rangle$) is  related to the behaviour of the sound speed ($c$) mostly in the central regions, and hence to the stellar evolutionary state. The representation of $\langle\delta\nu_{02}\rangle$ vs. $\langle\Delta\nu\rangle$ is in fact considered as a seismic diagnostic diagram allowing to derive stellar mass and age for main-sequence solar-like pulsators \citep{jcd88}. The corresponding seismic diagram for red-giant models is drawn in Fig.~\ref{label_d02}.a  which shows a linear  dependence of $\langle\delta\nu_{02}\rangle$ on $\langle\Delta\nu_0\rangle$, with a slope that increases as the mass decreases. In Fig.~\ref{label_d02}.b we plot the normalized quantity  $\langle\delta\nu_{02}\rangle/\langle\Delta\nu_0\rangle$ that in main-sequence stars is known to depend mostly on  central physical conditions \citep{RoxburghVorontsov03}.  It is worth noticing that we plotted  $\langle\delta\nu_{02}\rangle$ for models with different chemical compositions and convection treatment, nevertheless, a predominant  dependence on mass and radius appears. For a given mass  $\langle\delta\nu_{02}\rangle/\langle\Delta\nu_0\rangle$ increases with density contrast i.e. as the star expands with smaller values of $\langle\nu_0\rangle$, decreases as the mass increases, and  does not change significantly during the core He-burning phase.

In Fig.~\ref{label_d02}.a and b  a vertical dotted line indicates the lower limit of $\langle\Delta\nu\rangle$ measured from the first 34d of KEPLER mission,  and  the dot-dashed line corresponds to the fit  $\langle\delta\nu_{02}\rangle=0.122\langle\Delta\nu_0\rangle$ proposed by \cite{Beddingetal10} for those observations.  The comparison between our figures and  Fig.4 in \cite{Beddingetal10} indicates that the predictions from  theoretical models of low-mass stars (1-1.5\msun) in the low luminosity part of
the ascending RGB are consistent with  observational data, as the scaling based on $\nu_{\rm max}$ and $\Delta\nu$ also suggests. In that paper, the authors also report  the detection of $\ell=3$ modes, and because of  the scatter of $\ell=1$ mode  frequencies, they suggest to use a new small separation $\delta\nu_{03}(n)=0.5\,(\nu_{0\,n-1}-2\,\nu_{3\,n-2}+\nu_{0\,n})$  instead of the classic $\delta\nu_{13}$.  Note that the observational data provide a ratio $\delta\nu_{03}/\delta\nu_{02}\,\sim\,2.2$ \citep{Beddingetal10} instead of 2 predicted by the asymptotic theory. In the bottom panels of Fig.~\ref{label_d02} we plot the theoretical results for $\delta\nu_{03}$ and its relation with $\Delta\nu_0$ and $\delta\nu_{02}$.  \cite{Beddingetal10}  results are represented by a dash-dotted line while the dotted line corresponds to our fit $\delta\nu_{03}\sim 2.42\, \delta\nu_{02}$. The dependence of $\delta\nu_{03}$ on stellar mass and radius is similar to that of  $\delta\nu_{02}$, therefore no additional information should be expected from    $\delta\nu_{03}$.

The small frequency separation $\delta\nu_{01}(n)=0.5\,(\nu_{0n}-2\,\nu_{1,n}+\nu_{0\,n+1})$ in main-sequence stars is also known  to be  sensitive to the center physical conditions. The asymptotic theory predicts a $\langle\delta\nu_{01}\rangle=1/3\langle\delta\nu_{02}\rangle$ relationship. As it is evident in  Fig.~\ref{label_d02}.c  the ``p-mode'' spectrum for red-giant models does not follow those predictions. In particular, a large number of models have negative or very small values of $\langle\delta\nu_{01}\rangle$  independently of $\langle\delta\nu_{02}\rangle$.  Similar values of $\langle\delta\nu_{01}\rangle$ have also been found in the KEPLER data \citep{Beddingetal10} and in the oscillation spectrum of the CoRoT red-giant HR~7349 \citep{Carrier10}.  We note that the largest concentration of negative/small $\langle\delta\nu_{01}\rangle$ values   correspond to models ascending or descending  the RGB.  $\langle\delta\nu_{01}\rangle/\langle\Delta\nu\rangle$ vs. luminosity confirms this result.

Searching for a common characteristic in the structure of these models, we find that while ascending and descending the RGB the turning points of $\ell=1$ modes are well inside the convective envelope.
 The steady He-burning models have a shallower convective envelope and the turning points of $\ell=1$ modes are inside the radiative region. Fig.~\ref{label_d01L}.a shows the variation of  $\langle\delta\nu_{01}\rangle/\langle\Delta\nu\rangle$  with the distance (in  acoustic radius $\tau(r')=\int_0^{r'} dr/c$) between the bottom of the convective zone (BCZ) and the turning point for a $\ell=1$ mode with frequency close to $\nu_{\rm max}$ ($tp_1$). In Fig.~\ref{label_d01L}.b we highlight the behaviour of $\langle\delta\nu_{01}\rangle$ for low-mass RGB models.  The scatter of $\langle\delta\nu_{01}\rangle/\langle\Delta\nu\rangle$ rapidly decreases as $\tau_{tpl_1}-\tau_{\rm BCZ}$ increases (deep convective envelope) and  $\langle\delta\nu_{01}\rangle/\langle\Delta\nu\rangle$ takes negative values for models in which $tp_1$ is well inside the convective envelope.

\section{Concluding remarks}

We presented the properties of the  theoretical spectrum of solar-like oscillations during the  RGB and core He-burning phases of red giant evolution, and analyzed the behaviour of large and small frequency separations  computed from modes well trapped in the acoustic cavity of these stars. The main results of this global overview are the following:

\begin{itemize}
\item Independently of the evolutionary state, $\ell=2$ and $\ell=3$ modes trapped in the acoustic cavity have an inertia of the same order as that of the corresponding radial mode and behave as ``p-modes'' with frequencies regularly spaced by $\langle\Delta\nu\rangle$. As a consequence, the scatter of $\ell=2$ modes in the folded \'echelle diagrams is rather small.

\item The trapping of $\ell=1$ modes in the acoustic cavity depends on the evolutionary state. While a regular pattern of dipole modes is expected in more centrally  condensed models ascending the RGB , the scatter significantly increases for core He-burning ones. Therefore the regularity  of the $\ell=1$ spectrum  could be used to discriminate between different evolutionary phases.

\item $\langle\delta\nu_{02}\rangle$  and $\langle\delta\nu_{03}\rangle$  depend almost linearly on the large separation, hence on the mean density of the model, with a slope that slightly depends on the stellar mass.

\item $\langle\delta\nu_{01}\rangle$ seems to  reflect the distance between the $\ell=1$ turning point and the bottom of convective envelope.
It takes negative (or small) values if $tp_1$ is well inside the convective envelope, as it occurs in models ascending or descending the RGB.

\item  The properties of our adiabatic spectra are in good agreement with those found in KEPLER data for low luminosity RGB stars  \citep{Beddingetal10};  and also with those for CoRoT-exofield red giants. In this study we predict a large scatter of dipole modes for ``Red-Clump" stars  \citep[$\Delta\nu\sim 4\mu$Hz,][]{Hekker09, Mosser10, MiglioPop09}, as well as a regular pattern of $\ell=0$ and 2 modes for all the red giants. Both features correspond quite well with the properties of the spectra found by \cite{Mosser10b} in the  CoRoT-exofield red-giants (see Fig.~\ref{label_mosser} and their Fig.~3).

\end {itemize}

\acknowledgements
J.M. and A.M. acknowledge financial support from the Prodex-ESA Contract Prodex 8 COROT (C90199) and FNRS.  The authors thank M.A. Dupret for useful suggestions concerning oscillation frequency computations.


\section{Introduction}
\label{sec:intro}

Observational and theoretical studies of galaxy rotation curves
\citep{rubin70}, velocity dispersions of elliptical galaxies \citep{faber76},
baryon fraction in clusters \citep{clowe06,mass07}, gravitational lensing,
structure formation, CMB power spectrum etc. indicate that a significant
fraction of the mass of the Universe has no electromagnetic interaction
\citep{bertone05}. At galactic scale, the main evidence for the existence of
such gravitating matter with very high mass to light ratio comes from
observations of galaxy rotation curves. The gravitational potential derived
from the observed rotation curve can not be explained only by the visible mass
with any reasonable mass to light ratio
\citep{rubin70,rubin80,sofue96,sofue01,spano08}. Though there are modified
theories of the gravitation and other suggestions
\citep{milgrom83,beken84,sand86,fahr90,sanders97,brown06} to explain this
anomaly, the dark matter concept is widely accepted to be a simpler
explanation of these observations. Though observations of galaxy rotation
curves established the existence of the dark matter at galactic scale, there
is no general agreement on the nature and various properties (like mass
distribution) of this major constituent of the Universe \citep{nfw96,blok05}.
Potentially, galaxy rotation curves can also shed light on some aspects of the
dark matter properties and thus help us in deriving a better understanding of
the nature of the dark matter. One such key aspect is the density profile or
the mass distribution of the dark matter at galactic scale.

On the theoretical front, there are various mass distribution models of the
dark matter halo. The isothermal profile, the Navarro-Frenk-White (NFW)
profile and some variant of these two profiles are used extensively in both
theoretical and observational studies. An almost flat observed rotation curve
outside the core of a galaxy led to the dark matter halo model with $1/r^2$
isothermal density profile \citep{bege91,burk95}. The density at the centre
is not finite in this model. However, it is possible to impose appropriate
boundary conditions to derive a non-singular solution with density
$\rho(r)_{NIS} = \rho_0r_c^2/(r_c^2+r^2)$, where $\rho_0$ is the central
density and $r_c$ is the ``core radius''. At large radius, this profile will
be similar to $1/r^2$ singular isothermal profile. A cut-off radius
$r_{max}$ is also required to be imposed to keep the total mass of the halo
finite. However, the non-singular isothermal profile was not very successful
in explaining the observed rotation curves. A significant progress in the
field was the introduction of an alternative density profile, the NFW
density profile, $\rho(r)_{NFW} = \rho_0r_c^3/[r(r_c+r)^2]$, where $\rho_c$
and $r_c$ are the characteristic density and characteristic radius
respectively \citep{nfw96,nfw97,jing00}.

The NFW profile has been quite successful in explaining the observed rotation
curves for many galaxies. Over the time, a class of variant of the isothermal
and the NFW profiles are also introduced to fine tune the agreement of theory
and observations \citep{burk95,zhao96,fuku01}. All these density profiles are
of more or less similar merit, as far as explaining the rotation curve is
concerned. However, detailed analysis shows that the NFW profile has a number
of problems. For example, observations of galaxies and clusters suggest a
central flat density core \citep{blok07,kdn08}, whereas, in NFW profile, dark
matter density has a central cusp with a logarithmic slope $\approx -1$.
Numerical simulations show that the presence of baryon can change the mass
distribution significantly \citep[][and references therein]{ez04,sm10}, and
for isothermal cusp, minor merger and dynamical friction may lead to a
shallower central density slope \citep{rd08}. The NFW profile also seems not
to fit the observed rotation curves of the dark matter dominated low surface
brightness galaxies and low mass dwarf galaxies \citep{blok02,blok05}. A more
serious issue with the NFW profile is that the profile is derived by fitting
analytical function to the density distribution derived from numerical
simulation of structure formation with dark matter. This, in a sense, lacks a
proper physical understanding. Though some physical insight of the NFW profile
have recently come from further high resolution dark matter simulations
\citep[e.g.,][]{taylor01}, the issue is far from being settled. Moreover, even
if these cosmological simulations are very successful on large scale, at
galactic scale there are unsolved issues, like the ``angular momentum
problem'' \citep{sl99,sl01,burk04} or the ``missing satellite problem''
\citep{klypin99,moore99}, which are yet to be addressed.

Here, we have developed a self-consistent model of the dark matter halo
substructure distribution at galactic scale to explain observed NFW-like
rotation curves. Support for the existence of dark matter substructures has
mainly come from numerical simulations \citep{gc08,mdk08,vs08,etws09,lad09}.
But, there are strong observational hints like flux anomalies and time delays
in gravitational lensing \citep{chen09,km09,vk09,xu09}, or enhanced gamma rays
and leptonic cosmic rays \citep{ep09,pinz09}, indicating the presence of
substructures. The present model is based on the assumption of a scale free
nature of the dark matter clustering that leads to a statistically
self-similar distribution of the halo substructures at galactic scale. It is
shown that a simple fractal model of the dark matter halo substructure
predicts an NFW-like rotation curve. Such a model also predicts a scale free
power spectrum of the rotation velocity fluctuations. The model is described
in Section~\ref{sec:model}, and the results are presented in
Section~\ref{sec:result}. Possible limitations of our analysis are discussed
in Section~\ref{sec:disc}. Finally, we summarize and present our conclusions
in Section~\ref{sec:concl}.

\section{The Model}
\label{sec:model}

\begin{figure*}
\begin{center}
\includegraphics[width=6.25cm]{fig0.eps}\includegraphics[trim = 25mm 25mm 25mm 25mm, clip, width=6.25cm]{fig1.eps}\includegraphics[trim = 25mm 25mm 25mm 25mm, clip, width=6.25cm]{fig2.eps}
\caption{\label{fig:dmfrac} {\it Left:} Schematic representation of the halo substructure with three subhalos around any big halo. Note that both the cut-off radius ($r_{max}$) and the halo to subhalo distance ($D$) are scaled down by a factor $f_r$ for each substructure level. {\it Middle and Right} Halo substructure with seven subhalos around any big halo and with $f_r=0.33$. The middle and right panels show the structure with two and four substructure levels respectively.}
\end{center}
\end{figure*}

\subsection{Assumptions, Parameters and Constraints}
\label{sec:param}

Assuming that the dark matter clustering has a scale free nature (i.e., there
exist halo substructures of a wide range of mass), the density profile can be
described as a combination of a smooth radial profile $\overline{\rho(r/r_c)}$
and a stochastic part $\delta\rho$. Here $r_c$ is a characteristic ``core''
radius and $\overline{\rho}$ contains information of the density variation at
scales larger than or comparable to $r_c$. For the purpose of this work, we
have used a simpler model of this density distribution which, however, retains
all relevant key features. In this simplified model, we have assumed that each
dark matter clump has a number of smaller clumps around it. Each of these
smaller clumps are in turn just a scaled down version with even smaller clumps
around them. As a result, the whole structure has an approximate spherical
symmetry and a statistical self-similarity. It is assumed here that all these
clumps have non-singular isothermal density profile where the central density
($\rho_0$), the core radius ($r_c$), the cut-off radius ($r_{max}$) and the
halo to subhalo distance ($D$) are scaled down accordingly. However, this
specific density profile is not a crucial assumption in our model and the
individual substructures may have any non-singular density profile
$\rho(r/r_c)$, where $r_c$ is some characteristic radius. Essentially, this
is a fractal structure with three parameters: (i) $n$, the number of small
clumps around any clump, (ii) $f_r$, spatial scaling factor for core radius,
cut-off radius and distance and (iii) $f_\rho$, central density scaling
factor between any clump and its next smallest clumps. A fractal is a
fragmented and irregular geometrical shape with exact or stochastic
self-similar structures at all scales \citep{mand83}. Independent of the
value of $f_r$ and $n$, the (Hausdorff) fractal dimension of such a
structure is $3$ for $nf_r^3<1$. However, since at each iteration, the linear
size of clumps scales by a factor $f_r$ and the mass scales by a factor $n$,
the local mass dimension for any substructure level is
$D_m = -log(n)/log(f_r)$ over a certain range of scales. A mass dimension of
$D_m$ for a medium implies that the mass enclosed in a sphere of radius $r$ in
such a medium will be $M(r) = kr^{D_m}$. So, for the $N^{\rm th}$ substructure
level, $M_N(r) = k_Nr^{D_m}$ over a range of scales depending on $N$, $n$,
$f_r$, $D$ and $r_{max}$. Note that this range is different for different
substructure levels. Hence, for the complete structure, the total mass $M(r)$,
which is the sum of $M_N(r)$ of all the substructure levels, will not have a
simple power law radial dependence. But, the dark matter halo mass function
will still be a power law, $N(m) \propto m^{-\alpha}$, where the power law
index $\alpha = -log(n)/log(f_r^3 f_{\rho})$ is a more physically motivated
parameter of this model and can be constrained from theoretical and numerical
analysis of dark matter structure formation \citep{gao04,zemp09}. The density
distribution can be written as
\begin{equation}
\rho(r) = \rho_{\rm bg} + \displaystyle\sum_{i=0}^\infty\sum_{j=1}^{n^i}
\rho_s(\rho_{0,i}, r_{c,i}, r_{max,i}, \vec{r}_{i,j})
\label{eq:density}
\end{equation}
where $\rho_{\rm bg}$ is background density and $\rho_s(\rho_{0,i}, r_{c,i},
r_{max,i}, \vec{r}_{i,j})$ is density profile of individual substructure
centred at $\vec{r}_{i,j}$ with central density ($\rho_{0,i}$), core radius
($r_{c,i}$) and cut-off radius ($r_{max,i}$). Considering the self-similarity
of this model,
\begin{eqnarray}
&&\rho_{0,i} = f_\rho \times \rho_{0,i-1} \nonumber\\
&&r_{c,i} = f_r \times r_{c,i-1} \nonumber\\
&&r_{max,i} = f_r \times r_{max,i} \nonumber\\
&&\vec{r}_{i,j} = \vec{r}_{i-1,k} + (f_r \times D_{i-1} \times \vec{d})
\end{eqnarray}
where $k=n^{i-1}$, $\vec{d}$ is a unit length vector with random orientation
and the initial set of parameters $\rho_{0,0}, r_{c,0}, r_{max,0} = \rho_{0},
r_{c}, r_{max}$ is for the largest subhalo centered at the origin. In
principle, $\rho_{\rm bg}$ may be a smooth function of $r$. But, since we are
assuming it to be a small background density threshold, its effect on the final
rotation curve is not very significant. So, for simplicity, we have assumed
$\rho_{\rm bg}$ to be constant over the radius of our interest. This structure
is shown schematically (without any randomness for the shake of clarity) in
the left panel of Figure~\ref{fig:dmfrac}. After introducing randomness in
angular position of the subhalos, one realization of such a structure with
$n=7$ and $f_r=0.33$ is shown with two and four substructure levels in the
middle and right panel respectively. This model can be considered as a
simplified representation of the scale free, clumpy density structure of dark
matter above a small threshold density at the galactic scale.

We note that the parameters for this model are constrained to a good extent by
various physical considerations. Assuming that the structure is extended down
to the infinitely small scale, to avoid a divergence of the total mass, the
quantity $n f_r^3 f_{\rho}$ must be less than unity. Similarly, for a halo
with cut-off radius $r_{max}$, halo to subhalo distance $D$ must be greater
than or equal to $(1-f_r)r_{max}/(1-2f_r)$ to avoid any overlap. Also note
that the density of any subhalo at the cutoff radius $\rho(r_{max})$ scales as
$f_{\rho}$. So, we adopt $f_{\rho} = 1.0$ for further analysis to ensure a
constant density at the cutoff radius for all the substructures. In this
simplified model, we have assumed that all the parameters like $D$, $\rho_0$,
$r_{max}$ etc.~are identical for all the subhalos for a particular
substructure level and $n$, $f_r$ and $f_\rho$ remain constant for all the
levels. In a realistic situation, however, all these parameters may have some
random variation. As a result of these fluctuations, the halo mass function,
the radial mass distribution, and in turn, the rotation curves, are expected
to be somewhat smoother than that are predicted from this analysis.

\subsection{Tidal stability}
\label{sec:tidal}

A much stronger constraint on the parameters comes from the consideration of
the stability of this structures preventing tidal disruption by invoking
self-gravity. Considering a rigid object of mass $m$ and radius $r$ at a
distance $d$ from a bigger object of mass $M$ and radius $R$, from the
standard Roche limit consideration, $D$ should be $\sim r(2M/m)^{1/3}$ to
avoid tidal disruption of the smaller body. In the case of any halo and its
immediate subhalo, $m/M = f_r^3$ and $r/R = f_r$ imply $d \approx 1.25 R$. For
even smaller subhalo structures with $r/R = f_r^k$, the mass is scaled
accordingly $m/M = f_r^{3k}$ and $d \approx 1.25 R$ assures stability. Hence,
a distance
\begin{equation}
D \geq 1.25x\frac{1-f_r}{1-2f_r}r_{max}
\label{eq:roche}
\end{equation}
will make the whole structure stable. Here $x$ is a fudge factor and we use
the value of $x=1.1$ to accommodate non-rigid density clumps. For a given
value of $D$, the number of substructure $n$ is also constrained to be
\begin{equation}
n \leq \frac{4\pi D^2}{\pi (f_rD/1-f_r)^2} = \frac{4(1-f_r)^2}{f_r^2}
\label{eq:nsub}
\end{equation}
to avoid any possible overlap of subhalo structures.

\subsection{Virial stability}
\label{sec:virial}

A detailed virial stability analysis requires numerical simulation of the
dynamics of such a density distribution to get the time-averaged dynamical
quantities. But, a simple ensemble average virial scaling analysis may be used
to constrain the central density $\rho_0$ for a set of model parameter. As the
whole structure is assumed to have an approximate spherical symmetry, average
kinetic and potential energies, $\langle T \rangle$ and $\langle V \rangle$,
for thin spherical shells of radius $r$ and thickness $\delta r$ will be
\begin{eqnarray}
&& \langle T \rangle = \frac{1}{2}m\sigma_{\rm DM}^2 = \frac{1}{2} 4\pi r^2\delta r \rho(r) \sigma_{\rm DM}^2 \nonumber\\
&& \langle V \rangle = -\frac{GM(r)m}{r} = -v_c^2(r) 4\pi r^2\delta r \rho(r)
\end{eqnarray}
where $\sigma_{\rm DM}$ is dark matter velocity dispersion, $M(r)$ is total
mass within radius $r$
and $v_c(r) = (GM/r)^{1/2}$ is the scale dependent
virial velocity equivalent of the ``circular velocity'' for rotating disk.
Since the rotation curve has a roughly constant value $v_0$ at large radius,
the ratio $2\langle T \rangle/|\langle V \rangle|$ will tend to the equilibrium
value of $1$ at large radius if $\sigma_{\rm DM} \approx v_0$. Using the
minimum value for $D$ from equation~(\ref{eq:roche}), the maximum extent of
the structure $R_{max}$ will be $D/(1-f_r)$ and the average density will be
\begin{equation}
\langle\rho\rangle = \frac{3(1-2f_r)^3(s-tan^{-1}s)}{(1.25xs)^3(1-nf_r^3)}\rho_c
\label{eq:vir}
\end{equation}
where $s=r_{max}/r_c$. Note, however, that this is a fractal structure with
significant porosity. So, the average density of any individual clump is
higher than $\langle\rho\rangle$ by a factor of
$(1-nf_r^3)(R_{max}/r_{max})^3$. Now, for global stability of the whole
structure, $R_{max}$ should be less than or equal to the radius within which
the virial equilibrium is maintained. Average density within this virial
radius $r_{vir}$ or $r_{200}$ should be $\approx 200$ times more than the
critical density $\rho_{crit} = 3H^2/8\pi G$. For a choice of model parameters
$n$, $f_r$, $s$ and $x$, this will constrain the lower limit of the central
density $\rho_0$ so that $\langle\rho\rangle \gtrsim 200 \rho_{crit}$. For
individual substructures, both mass and volume scale as $f_r^3$, keeping the
average density constant. This implies that stability for one substructure
level ensures stability for all other levels.

\section{Results}
\label{sec:result}

\subsection{Rotation curve}
\label{sec:results1}

Since the halo density distribution is significantly clumpy, the velocity
field for such a system is also expected to have fluctuations at all scales.
However, due to approximate spherical symmetry of the clump distribution, the
average rotation velocity over a spherical shell at radius $r$ will still be
$\langle v_c(r) \rangle \approx (GM_r/r)^{1/2}$, where $M_r$ is the total mass
within this radius. As pointed out in subsection \ref{sec:virial} using the
virial stability argument, the virial velocity or, equivalently, the
``rotation'' velocity is expected to be approximately same as the local
velocity dispersion. This derived rotation curve for the fractal model is found
to be NFW-like at large radial distance. At small radius, by construction, the
rotation curve is obviously exactly same as that of a non-singular isothermal
halo. This is shown in Figure~\ref{fig:rotcur}. Radial distance is scaled by
$r_c$ and rotation velocity is scaled by $v_0$, which is the rotation velocity
at the maximum radius plotted. The black curve is the predicted rotation curve
for the fractal substructure model with $n=35$, $f_r=0.25$, $f_{\rho} = 1.0$
and $r_{max} = 6.0 r_c$ and the red line is the best fit NFW profile to that.
The background density threshold $\rho_{bg}$ is assumed to be negligible in
this model. Rotation curves for non-singular isothermal sphere with and
without a cutoff (green and blue curves respectively) as well as for NFW halo
(magenta curve) are also shown in Figure~\ref{fig:rotcur} for a comparison.
The velocity scaling ensures that the total mass encompassed by the maximum
radius is same for all models.

\begin{figure}
\begin{center}
\includegraphics[width=9cm]{fig3.eps}
\caption{\label{fig:rotcur} Predicted rotation curve for the fractal substructure model (black) and the best fit NFW profile to that (red). This is for $n=35$, $f_r=0.25$, $f_{\rho} = 1.0$ and $r_{max} = 6.0 r_c$. Rotation curves for non-singular isothermal sphere with and without a cutoff (green and blue curves respectively) as well as for NFW halo (magenta curve) are also shown for a comparison. Radial distance and rotation velocity are scaled by $r_c$ and $v_0$ (rotation velocity at the furthest radial distance) respectively.}
\end{center}
\end{figure}

\begin{figure}
\begin{center}
\includegraphics[width=9cm]{fig4.eps}
\caption{\label{fig:rbgcur} Predicted rotation curve for the fractal substructure model for different background density threshold $\rho_{bg}$. All other parameters are same as in Figure~\ref{fig:rotcur}. Different curves are for $\rho_{bg}/\rho(r_{max}) = 0.00$, $0.03$, $0.10$ and $0.30$ (black, red, green and blue curve respectively).}
\end{center}
\end{figure}

The effect of the background density threshold $\rho_{bg}$ is shown in
Figure~\ref{fig:rbgcur}. Here, the rotation curve is evaluated for different
$\rho_{bg}$ keeping all other parameters same as in Figure~\ref{fig:rotcur}.
For scaling, we have used the density $\rho(r_{max})$ at the cutoff radius
$r_{max} = 6.0 r_c$. Different curves in Figure~\ref{fig:rbgcur} are for
$\rho_{bg}/\rho(r_{max}) = 0.00$, $0.03$, $0.10$ and $0.30$ (black, red, green
and blue curve respectively). The general NFW-like nature of the rotation
curve and the radial fluctuations remain unchanged. But depending on the value
of $\rho_{bg}$, rotation curve at large radius may be rising, flat or
declining. Note that the derived rotation curves shown in
Figure~\ref{fig:rbgcur} are with the simple model of a constant $\rho_{bg}$.
In reality, $\rho_{bg}$ is expected to be decreasing with increasing $r$,
giving rise to a rotation curve somewhat intermediate between the extremes
shown in Figure~\ref{fig:rbgcur}.

\subsection{Velocity fluctuations}
\label{sec:results2}

The predicted rotation curve due to the clustering of dark matter subhalos is
very similar to the observed rotation curve and the empirical NFW rotation
curve. However, unlike the NFW model with a smooth radial density distribution,
the present model predicts a significant fluctuation of rotation velocity in
both angular and radial directions. Since the underlying density field, which
gives rise to this velocity fluctuations, is scale free, the velocity
fluctuation power spectrum is expected to be a power law. Though the fractal
model has many free parameters, the only relevant parameters for the index of
this power law are $n$ and $f_r$ while the rest of them will just introduce
different multiplicative scaling. This prediction can be easily verified from
high spatial and spectral resolution observation of neutral hydrogen of normal
galaxies. Note that part of this fluctuations will cancel out for the
spherically averaged rotation curve and hence it is important to use the full
velocity field to search for such scale free fluctuations. It is also
important to note that fluctuations of the velocity field of the hydrogen gas
will have contributions from the local density perturbations of the disk. But
the scale dependence of these perturbations may makes it possible to
disentangle the scale free perturbations due to the halo substructures. There
is some indication of such a power law scaling of the rotation velocity
fluctuation power spectra from direct observation and analysis of the
H~{\sc i} 21~cm velocity field of some nearby galaxy \citep{mdr10}. We leave a
more detailed treatment of this aspect to a future work.

\section{Discussion}
\label{sec:disc}

Throughout this analysis we have assumed that the clump distribution in this
model has an approximate spherical symmetry. This assumption is most likely
to be untrue in reality. Due to small count, bigger substructures are more
likely to cause departure from spherical symmetry. In this situation, rotation
velocity field will also have strong azimuthal asymmetry. Observationally,
about half of the spiral galaxies show some degree of kinematical lopsidedness
\citep{rs94,hay98,sw99,sofue01,jc09} which may originate from tidal distortion
partly due to deviation from spherical symmetry of the underlying
gravitational potential of the lopsided dark matter halo \citep{sc09,saha09}.
Due to their large number, smaller substructures are likely to have relatively
more symmetrical distribution. So, the velocity fluctuations spectrum is
expected to show a power law scaling at large spatial frequency (large $k$,
i.e., small physical scale), and a departure from power law at small $k$.

We have also assumed in this analysis that the baryons will not significantly
affect the galactic dynamics. However, some recent numerical studies
\citep[e.g.][]{wd08,rd08,rd09} have shown that the presence of baryons have
effects like flattening the central cusp, reducing the halo triaxiality, and
introducing bias in clustering with environment density. Effectively, baryon
dissipation may destroy the similarity between different scales of clustering.
However, \citet{rd10} have found that phenomenon like ``efficient feedback
from stellar evolution and the central supermassive black holes'' may
counterbalance this effect to some extent. On the other hand, \citet{kne10}
have found no significant effect of baryonic physics on properties like shape
and radial alignment of substructures. Numerical simulations also suggest
that, even in the presence of baryons, the subhalo mass function remains a
power law though the power law index changes from $-0.99$ to $-1.13$
\citep{rd10}. It is, hence, important to keep in mind that, in a realistic
situation, baryon dissipation may alter the velocity fluctuations causing a
significant departure from the scale free velocity fluctuations.

We note that the key result of this analysis, that is an NFW-like rotation
curve for the fractal model, is not crucially dependent on the exact density
profile of individual halos. A variety of density profile without any central
singularity (variant of non-singular isothermal sphere) will leads to a similar
NFW-like rotation curve. This strongly suggests that the clustering properties
of the dark matter particles dominantly govern the radial density profile of
the halo. Finally, these assembly of substructures will give rise to a flat
NFW-like rotation curve for normal galaxies. But, for more dark matter
dominated low mass dwarf galaxies and low surface brightness galaxies,
depending on the exact form of the non-singular density profile, the dominant
contribution of the central big halo may make the rotation curve for such
galaxies intrinsically different which is consistent with observational
results \citep{blok02,blok05}.

\section{Conclusions}
\label{sec:concl}

\begin{enumerate}
\item The dark matter halo substructures at galactic scale is found to
significantly affect the rotation curve.
\item A self-consistent model of statistically self-similar hierarchical dark
matter substructures predicts an NFW-like rotation curve, though each clump
has a non-singular isothermal density profile. This NFW-like rotation curve
emerges out of the fractal geometry and is independent of specific density
profile of individual clumps.
\item The model predicts a scale free power spectrum of the rotation velocity
fluctuations which can be observationally verified.
\item The model also provides some plausible explanation of the observed
intrinsic difference between the dark matter halo density profile of normal
galaxies and dark matter dominated low mass dwarf galaxies and low surface
brightness galaxies.
\end{enumerate}

\acknowledgments

We are grateful to Rajaram Nityananda for many useful comments on an earlier
version of this paper. We thank Susmita Chakravorty, Aritra Basu, Sanjay
Bhatnagar, Abhirup Datta, Prasun Dutta, Sanhita Joshi and Chandreyee Sengupta
for helpful discussions. We are also grateful to the anonymous referee for
useful comments and for prompting us into substantially improving this paper.
NR is a Jansky Fellow of the National Radio Astronomy Observatory (NRAO).
This research was supported by the NRAO. The NRAO is a facility of the
National Science Foundation (NSF) operated under cooperative agreement by
Associated Universities, Inc. (AUI).

{\it Disclaimer}: This is an author-created, un-copyedited version of an
article accepted for publication in the Astrophysical Journal. IOP Publishing
Ltd is not responsible for any errors or omissions in this version of the
manuscript or any version derived from it. The definitive publisher
authenticated version will be available online at {\tt http://iopscience.iop.org/}.


\section{INTRODUCTION}
\label{sec:intro}
The VLT Survey Telescope (VST) project is a cooperation between
the European Southern Observatory
(ESO) and the Italian National
Institute of Astrophysics (INAF) for the
design and realization of a wide field Alt-Az
telescope of 2.61-m aperture, specialized
for high quality astronomical imaging,
to be installed and operated at the ESO
Paranal Observatory in Chile close to the
VLT units.

In this paper we briefly describe the optical\cite{2005SPIE.5962..608M} and mechanical\cite{2003SPIE.4837..379M} design of the VST telescope in order to have all the elements to finally describe the alignment procedure
we planned to follow for the mechanical and optical components of the telescope.

\section{VST Design}
The VST is a 2.61 m Alt-az f/5.5 modified Ritchey-Chretien telescope with a corrected Cassegrain Field of View (FoV) diameter of 1.47degrees,
coupled with the OmegaCAM\cite{2002Msngr.110...15K} instrument, a 16k$\times$16k CCD mosaic
camera of 15$\mu$m pixel size and 0.21arcsec/pixel scale.
 In order to cover the wide field of view with the required high image quality, the telescope is not a pure Ritchey-Chretien, but the design includes two different camera correctors respectively with and without dispersing
elements to compensate atmospheric dispersion Also the OmegaCAM dewar window is a low power lens
itself.
One corrector camera is made of two lenses (L1 and L2 hereafter) and operates from U to I bands. The
other one is made with an atmospheric dispersion corrector (ADC) and one different lens (L3) to observe from B to I bands. The two correctors can be exchanged by mean of a large translator stage switching the two systems. The telescope is provided with active optics both for secondary and primary mirrors in order to compensate for telescope
structure deformations.

The telescope opto-mechanical design includes an instrument interface flange, which divides the
telescope in two blocks along the optical axis. The first block includes mirrors and primary mirror cell with correctors and a pick-up
mirror for the active optics (Shack-Hartmann) sensor and the second block includes filters, dewar window lens and the OmegaCAM instrument camera.
The secondary mirror is mounted on a dedicated high resolution positioning system. This system is
a double stage hexapod device which guarantees micro metric and sub-micro metric
displacements and high resolution tilts of the secondary mirror (M2) with respect to the primary (M1).

The final goal of the alignment procedure is to place the M1, M2 and the two correctors aligned along the optical axis and at the proper distances within the range of adjustment offered by the active optics system.

\section{Telescope Alignment}
\label{sec:align}
The alignment procedure was inspired by the one used to align the mirrors of the Large Binocular Telescope (LBT) and the LBT prime focus Camera (LBC)\cite{2004SPIE.5492..513D,2008SPIE.7014E.159G}.
We will use an ad-hoc alignment tool, built by Tomelleri S.r.l., to align the optical components of the VST telescope. This is a pure optical procedure, which does not trust in the precision of any mechanical device, such as in other possible alignment procedure followed for the Telescopio Nazionale Galileo (TNG) and the Very Large Telescope (VLT).
We designed an ad-hoc device (an ``alignment tool'') for this task which will be used for the M1 and M2 optical alignment.
\subsection{Mechanical Alignment}
A pre-requisite for the optical alignment is the mechanical alignment of the altitude and azimuthal axes with respect to the OmegaCAM instrument bearing that we will use as reference for the optical procedure.
In order to have a good starting position that allows to have most of the dynamic range available on the M2 hexapod we plan to shim the spider of M2 to center it with respect to the mechanical axis of the telescopes. The spider is mounted using flanges mounted at $\approx$45degrees (Figure~\ref{fig:fig1}), considering that this shimming will introduce extra Z axis shift of the same order of the decentering introduced.
Similar shimming will be applied to the primary mirror cell, but only for tilt, being installed for the centering of the M1 mirror three translation pads mounted on screws.

   \begin{figure}
   \begin{center}
   \begin{tabular}{c}
   \includegraphics[width=8cm]{Fig1.eps}
   \end{tabular}
   \end{center}
   \caption[fig1]
   { \label{fig:fig1}
The figure shows the M2 cell and the spider.}
   \end{figure}

The reference for the alignment of the M1 cell and secondary mirror spider will be a sighting telescope mounted below the azimuth axis on the basement of the VST dome (the interface is already installed).
In this way we can also define the encoder Zenith value for the elevation movement of the telescope, which places on the same plane the two axes defining the Azimuth and the OmegaCAM rotation bearing. Such as angle will be taken into account in the pointing model of the telescope.
First we will adjust the M2 mounting (using the target mounted on the M2 cell) shimming eventually the spider, then we place the autocollimator target on the instrument bearing and we will adjust the cell tilt mechanically. Such as an operation will be repeated at the very end of the alignment procedure.
The alignment of the bearing with respect to the azimuth axis will prevent field rotation issues because if bearing axis is different with respect to the telescope pointing movements  it will have a different centre of rotation. Please notice that a tilt of the bearing with respect to Azimuth generates a misalignment for the field correction derotation, while a decenter will generated only a different speed for derotation visible in long exposures images.
\subsubsection{Operations}
The sighting telescope will be mounted on the fixed platform below the VST telescope (see Figure~\ref{fig:fig2}) and the target mounted on the azimuth flange. The sighting telescope will be aligned with respect to the azimuth axis by rotating it.
A second target mounted on the OmegaCAM bearing will be adjusted in centering by rotating the bearing itself with respect to the axis defined the bearing itself.
The telescope altitude axis will be moved in order to reach the closest position between the axis of the OmegaCAM bearing and the axis of the telescope azimuth. We will take note of the encoder value of the altitude motion relative to this position because this will define the Zenith position.
A third target will be mounted on the nominal center of the M2 hexapod, with the hexapod at the mid of the available range of motion, in order to measure the decenter of the M2 itself.
This configuration will allow to compute the tilt of the OmegaCAM bearing with respect to the Azimuth motion and to design the shims needed to compensate it. The shim will be inserted below the four flexion bars. This operation implies the demounting and remounting of the cell (without the OmegaCAM dummy mounted). The manufacturing of the shims will be done in the mechanical workshop of the Paranal Observatory.
The measurement relative to M2 target will provide the decentering of the M2 spider, that will be corrected by shimming the spider and/or the hexapod mounting.
After the remounting of the M1 cell and M2 hexapod and spider, the two targets will be realigned, the elevation encoder value measured again, and the decentering and tilt of the cell and the decentering of the M2 measured verified.
The M2 reference center position will be the one taking into account the flexure expected at a mid elevation altitude in order to save range to compensate the decentering of the M2 due to off the Zenith telescope elevations.

   \begin{figure}
   \begin{center}
   \begin{tabular}{c}
   \includegraphics[width=8cm]{Fig2.eps}
   \end{tabular}
   \end{center}
   \caption[fig2]
   { \label{fig:fig2}
The figure shows the sighting telescope ``STaz'' mounted below the VST azimuth axis with a special flange.}
   \end{figure}

After this phase we will check the ``alignment tool'' mounting on board of the telescope: we will test the procedure to mount this tool and the one to dismount it. We strongly recommend to verify this procedure after the mechanical alignment being the last opportunity to have the M1 dummy in the cell. We will check also the wood protection of the M1 mirror to be mounted below it. This device will be attached the alignment tool to the instrument bearing interface.
\subsection{Alignment Procedure}
As before, the reference for the alignment will be the 2$^{\rm nd}$ bearing rotation axis (the instrument one), which precision will set the limit of the final alignment quality.
The ad-hoc mechanical device is mounted on the bearing axis using a flange (see Figure~\ref{fig:fig3}) and, a mechanical vertical extension (see Figure~\ref{fig:fig6}) supports an horizontal rail (see Figure~\ref{fig:fig4}) to the area immediately above the primary mirror. This rail on a side (the short one) supports the laser and on the other, long about as the radius of the M1, supports a 2~inches folding mirror and three 50/50 beam splitters. On this rail the laser is able to fire (remotely commanded) a beam to the secondary or to the primary mirror using these folding mirrors (beam-splitters): the selection of the beam up- or down- warding can be accomplished by a moving cardboard (actually steel made) with selection holes.
Such as device will be mounted on a flange centered with respect to the optical axis, passing through the ADC corrector barrel. Of course the ADC correctors lens and prisms have to be removed from this barrel and they will be mounted only at the end of this procedure.
The lower baffling above the M1 has to be removed. This is doable without the removal of the cover (petals) of the mirror or other elements in the Serrurier truss area. However, in order to place the alignment device, the petals covering the M1 will be removed.
The alignment operation will start with the two lenses of the corrector mounted in the respective barrel, while the ADC prisms and the lens are not mounted in the ADC barrel.
During the alignment the stage mounting the corrector barrels is in the so called ``ADC'' position, in this way there is the physical room to have the ``alignment device", ``T-Tower" hereafter, going through the hole of the M1 mirror. The corrector lenses L1 and L2 can be placed in position without dismounting nothing of the telescope (for example the bearing). The placing in position of the two prisms and of L3 lens is not possible without dismounting the corresponding barrel: therefore the so called ``ADC'' stage will be dismounted from the cell and put back in position using mechanical pins as precision references.
We provide the T-Tower with an insertable detector (``CCD'' in the following) mounted on the focal plane area of the telescope. Such as device will be unavoidably subjected to flexures due to gravitational forces. However, as the telescope is pointing to the local Zenith and the only movement is around the gravity axis (rotation of the bearing) such as flexures will be invariant for the whole alignment operations.

The alignment device is mounted directly to the bearing with a ``Y" flange (see Figure~\ref{fig:fig3}) and with a vertical truss (see Figure~\ref{fig:fig6}) passing through the ADC barrel in order to match the bearing axis. The horizontal arm of the T-Tower will stay $\approx$700mm far from the M1 dish.

   \begin{figure}
   \begin{center}
   \begin{tabular}{c}
   \includegraphics[width=8cm]{Fig3.eps}
   \end{tabular}
   \end{center}
   \caption[fig3]
   { \label{fig:fig3}
The figure shows the flange designed for the T-Tower. It is mounted on the instrument bearing. The hole let the light passing through, either for the adjustable laser beam or for the CCD looking to the light fired by the laser mounted on the horizontal flange above the M1.}
   \end{figure}

   \begin{figure}
   \begin{center}
   \begin{tabular}{c}
   \includegraphics[width=12cm]{Fig4.eps}
   \end{tabular}
   \end{center}
   \caption[fig4]
   { \label{fig:fig4}
   The horizontal arm of the T-Tower, a card-board mechanism will switch between the configurations with beams going to M1 or going to M2. This arm will be mounted on a truss made in way to let the beam coming from the M2.}
   \end{figure}

On top of this, about 700mm far from the M1, an horizontal flange hosts a laser: this fires a beam to the secondary or to the primary mirror. The selection of the up- or down- warding direction is done by using beam splitters and a small folding mirror. The vertical mechanical extension is a truss divided in two parts in order to let beams coming from M2 to reach the insertable CCD at the level of the focal plane.
In this way it is possible to realize both chief and marginal rays. By rotating the bearing with the laser turned-on is possible to illuminate different portions (rings) of the mirrors (M1 and M2) and adjust iteratively centering and tilt of both M1 and M2. The selection of the chief vs. marginal rays can be done through a cardboard shutter possibly remotely controlled.
   \begin{figure}
   \begin{center}
   \begin{tabular}{c}
   \includegraphics[width=8cm]{Fig5.eps}
   \end{tabular}
   \end{center}
   \caption[fig5]
   { \label{fig:fig5}
   The picture shows lower baffle that must be removed  to place in position the T-Tower. It can be removed passing through the Serrurier truss.}
   \end{figure}

   \begin{figure}
   \begin{center}
   \begin{tabular}{c}
   \includegraphics[width=8cm]{Fig6.eps}
   \end{tabular}
   \end{center}
   \caption[fig6]
   { \label{fig:fig6}
   Sketch of the components needed for the ``LBT" like optical alignment. The T-tower vertical support will pass through the empty barrel of the ADC corrector. The two ADC prisms and L3 lens will be inserted in a later stage, after the alignment of M2, M1 and Corrector (L1+L2).}
   \end{figure}

   \begin{figure}
   \begin{center}
   \begin{tabular}{c}
   \includegraphics[width=8cm]{Fig7.eps}
   \end{tabular}
   \end{center}
   \caption[fig7]
   { \label{fig:fig7}
   The figure shows a sketch of the alignment of the M2 mirror.}
   \end{figure}

\subsection{M2 and M1 alignments}
The procedure will start with the M2 alignment, only after this phase M1 will be placed in the correct position. The alignment scheme of M1 and M2 is the following:
\begin{enumerate}
\item the CCD mounted at the bottom of the T-Tower is illuminated with any kind of source fixed with respect to the fixed part of the telescope. By rotation of the device the illumination pattern does apparently rotate onto the CCD and one can individuate the position on the detector that is intercepted by the bearing axis. In the following such a position is called ``center of the CCD".
We will rotate the CCD only using the bearing and illuminating it with a fixed led light: in this way we see a pattern which rotates and we can identify the rotation pivot pixel.
\item 	Using the configuration in Figure~\ref{fig:fig7}, the ``chief" ray is materialized by the laser and illuminates the secondary mirror. It is adjusted in order to match approximately the CCD center. Now the bearing is rotated and the secondary mirror is adjusted with only tilt in order not to move during the rotation;

\item	The ``marginal" ray is materialized and the bearing is rotated and the secondary is adjusted in order to have no apparent movement of the reflected spot with only rotation around the center of curvature of the osculating sphere on the M2;
\item	Steps 2 and 3 are iterated till to reach the no movements of the spot condition on the CCD during the bearing rotation;
\item	NOW the M2 is aligned with the bearing.

   \begin{figure}
   \begin{center}
   \begin{tabular}{c}
   \includegraphics[width=8cm]{Fig8.eps}
   \end{tabular}
   \end{center}
   \caption[fig8]
   { \label{fig:fig8}
   The figure shows a sketch of the alignment of the M1 mirror: the laser beam is reflected by the M1 and M2 before reaching the CCD mounted on the Y shaped flange at the level of the expected focal plane.}
   \end{figure}

   \begin{figure}
   \begin{center}
   \begin{tabular}{c}
   \includegraphics[width=8cm]{Fig9.eps}
   \end{tabular}
   \end{center}
   \caption[fig9]
   { \label{fig:fig9}
   The figure shows a CAD view of the alignment tool system T-Tower. Super-imposed to the image are the optical rays as they actually will fly from the laser reflected towards the M1, back reflected, passing through the holes made on the T-arm and directed to the M2, back reflected again and imaged on the CCD mounted on the bottom of the truss.}
   \end{figure}

\item	The laser beam is now fired toward the primary mirror, reflected by the secondary and finally reaches the CCD mounted on the T-Tower Y-flange. Procedures similar to the steps 2 and 3 are now iterated using the primary mirror adjustments. Using the chief ray we set tilt, and using the marginal one we set decentering. At the end of this level the primary mirror is also co-aligned;
\item After the alignment of M1 we plan to look to a bright star on the sky in order to refine focus and mirror position looking to the pupil image. In this case we will remove the T-tower truss, but not the CCD that will be used both as reference and as detector;
\item Still looking to a star it is possible to make out of focus images and correct coma through M2 rotation around its center of curvature thanks to a trials and errors procedure.
\end{enumerate}
\subsection{Corrector(s) alignment}
The insertion of the corrector lenses, using the ADC stage mechanism, is possible while the laser is in position on the alignment tool. The insertion of the L3 lens and of the two prisms inside the ADC barrel will be done in only on a later stage. We have provision for pins of the corrector mounting mechanics with respect to the M1 cell in order to remove it, insert L3 and prisms in the barrel, and place the mechanics in to its old position. The pinning was already done in Italy before the shipment of the M1 cell and correctors systems.
In principle, since the optical tolerance of the correctors are quite large (of the order of 3-5mm in centering, because of the low optical power of the lenses) there is a not negligible probability that mechanical precision is enough to have the corrector elements already aligned, therefore we planned to point on sky, and to simply place the corrector lenses barrel on axis and to have a first look at the sky in the center of the FoV: if the dimension of the imaged stars on the focal plane will be in the specification the corrector barrel can be considered aligned, otherwise we will follow the following steps.
An adjustable laser is installed on the bearing, using the ``Y" base of the T-Tower, and adjusted iteratively with the help of a close by and a far away CCDs in order to materialize the bearing (and now optical) axis.
The two CCDs for the alignment of the laser will be mounted on ad-hoc mechanical mountings:
\begin{itemize}
\item	 in front M2, mounted on a ad-hoc flange (see Figure~\ref{fig:fig10});
\item on the M1 cover;
\end{itemize}
   \begin{figure}
   \begin{center}
   \begin{tabular}{c}
   \includegraphics[width=8cm]{Fig10.eps}
   \end{tabular}
   \end{center}
   \caption[fig10]
   { \label{fig:fig10}
   The figure shows a CAD view of the interface to be mounted in front of M2 to place the CCD used for the definition of the optical axis.}
   \end{figure}
The barrel with the corrector lenses is place in position using the ADC translation stage and the laser spot refraction and the back reflection are used to adjust the centering and tilting of the lenses. We start with the lower one (L2):
\begin{enumerate}
\item	Tilting of L1 and L2 will be done using ad-hoc 0.1mm shims on the barrel;
\item	Centering is done using the translation correctors stage and moving the whole plate, on which the stage is mounted, on the orthogonal direction (there is a provision for this movement);
\item	In the case of L1, while is not in the L2 case, light back reflected from the first of the two corrector lenses, once this is aligned, is useful, through interference to align the second one.
\item	Please notice that the back reflected beam is observable only through a holed screen and projected there, while the refracted beam can be measured through one of the CCDs used for the laser alignment.
\item	We compute the size of the shims to be inserted on L1 and L2;
\item The corrector barrel will be dismounted and the shims will be inserted on the side of L1 and L2;
\item	The corrector is mounted back;
\item   Adjustable Laser aligned as described at the step 5;
\item	L1 and L2 moved in position using the corrector stage;
\item	alignment will be verified.
\end{enumerate}
Now the two corrector lenses are aligned and the procedure for L3 and prisms alignment starts.
\begin{enumerate}
\item	The ADC -Corrector system will be dismounted the two prisms and the L3 lens will be place inside the ADC barrel and the system mounted back;
\item	The laser will be re aligned, but this time passing through the corrector lenses (not the ADC), when the laser is aligned the beam in principle is not deviated by the L1 and L2, therefore we can use the 2CCDs procedure described before and we will measure, as we planned for the L1 and L2, the position of the reflected beam and we will compute the shim size to be inserted;
\item	The ADC -Corrector system will be dismounted and the shims will be inserted on the side of the L3 and of the two prisms;
\item	The ADC-Corrector is mounted back;
\item	Laser aligned as described at the step 2;
\item	ADC barrel placed in position;
\item	Alignment will be verified for the ADC barrel lenses;
\item	OPTIONAL: Alignment will be verified for the Corrector barrel lenses again.
\end{enumerate}
After this operations the ADC and the corrector are aligned to the telescope optical axis.
\section{Conclusions}
In this paper we described the procedure we plan to follow at the end of this year (according to the last schedule, 2010) for the alignment of the optical and mechanical components of the VST. We designed several opto-mechanical tools to optimally accomplish this goal. The tools are ready and wait to be tested on the telescope. In these last lines we recall that the procedure described above is purely optical and relies only in the precision of the optical axis definition: in our case the instrument bearing of the telescope. However wobble and run-out of the bearing had been already measured with the bearing mount on the M1 cell before the shipment to Paranal Observatory and are fulfilling the (optical) specification we need.

\acknowledgments
Since this paper offers an overview of the final step of the building of the VST telescope with the final positioning of the optical elements, the actual work was made by the INAF-Osservatorio di Capodimonte team, among them we remember M. Capaccioli, P. Schipani, S. d'Orsi, M. Brescia, L. Marty and D. Fierro.
We would like to say thanks to the ESO J. Spyromilio and R. Tamai for useful and meaningful discussions.


\section{}

The study of low temperature molecular collisions and interactions is a rapidly expanding area of research at the interface of chemistry and physics \cite{Carr:2009oz,Krems:2008pccp}.  Experimental and theoretical techniques have been developed for studying a wide variety of phenomena, ranging from quantum threshold scattering \cite {Gilijamse_OH_Xe,Sawyer:2008qq}, inelastic atom - molecule collisions \cite{campbell2009mechanism, Zahzam:2006la,Staanum:2006it}, external field control of dipolar interactions in cold and ultracold molecules \cite{Sawyer:2010mb, Ni:2010xy}, to chemistry at cold \cite{PhysRevLett.100.043203, Cresu:06ag,PhysRevLett.105.033001} and ultracold temperatures \cite{Krems:2008pccp,Ospelkaus:2010sci}. The immense diversity of molecular structure and interactions is the cornerstone of many applications of cold molecules in quantum information science \cite{DeMille:2002xr}, condensed-matter physics \cite{barnett:190401}, precision measurement \cite{Hudson:02ed,Vutha:10jp}, cold controlled chemistry \cite{Softley:09mp}, and astrophysics \cite{Akyilmaz_NO_depletion}.  Achievement of these applications relies on improved methods for cooling and further understanding of low temperature molecular collisions.

Recently, several experimental techniques were developed to study molecular collisions at cold and ultracold temperatures, beyond the reach of traditional molecular beam experiments \cite{Scoles}.  Ultracold ground state molecules with temperatures below 1 $\mu$K can be created via coherent state transfer from magneto-asssociated ultracold atomic gases, producing KRb \cite{Ni:2008eq} and Cs$_2$ \cite{Danzl:2010cs}.  Alternatively, direct cooling techniques such as buffer-gas cooling \cite{Hummon:2008rm} and Stark deceleration \cite{Gilijamse_OH_Xe} provide more general methods for cooling a larger class of molecules, enabling the production of cold molecules starting from room-temperature sources as an input, such as supersonic beams \cite{Gilijamse_OH_Xe,Krems:2008pccp,Carr:2009oz}.  Since the maximum molecular trap densities for these direct cooling techniques are, so far, $10^8-10^9$~cm$^{-3}$,  the collision experiments with these molecular samples typically use a dense rare gas collision partner \cite{campbell2009mechanism, Gilijamse_OH_Xe} or perform collisions at energies $> 10$~cm$^{-1}$ \cite{Gilijamse_OH_Xe,Sawyer:2008qq,Cresu:06ag}.  Further cooling and compression of molecular samples produced via direct cooling techniques is key to study cold molecule-molecule interactions and their applications. One possible approach is to sympathetically cool a molecule via collisions with an atom \cite{Carr:2009oz,Lara_Rb_OH,zuchowski:022701, Tacconi_NH_Rb,Wallis:2009ul,Barker:09nj}. Thus, studies of cold atom-molecule collisions not only uncover the underlying collision physics, but also may lead to important new methods for the production of a variety of ultracold molecules.

A number of theoretical studies have focused on the optimal selection of atomic collision partners for sympathetic cooling of molecules \cite{Lara_Rb_OH,zuchowski:022701, Tacconi_NH_Rb,Wallis:2009ul,Barker:09nj}. Almost all of these studies focused on the alkali-metal atoms, which can be easily prepared at milli-Kelvin temperatures via laser cooling.
However, accurate quantum scattering calculations have shown that collision-induced inelastic relaxation of molecules such as OH \cite{Lara_Rb_OH}, ND$_3$ \cite{zuchowski:022701} and NH \cite{Tacconi_NH_Rb} with the alkali-metal atoms occurs rapidly, causing molecular trap loss and severely limiting the efficiency of sympathetic cooling. A recent theoretical study suggested that sympathetic cooling of NH molecules by collisions with laser-cooled alkaline earth atoms such as Mg might be possible at low magnetic field strengths \cite{Wallis:2009ul}. A related study indicated the possibility of sympathetic cooling of large molecules by collisions with rare-gas atoms in an optical dipole trap \cite{Barker:09nj}. We have recently suggested that cold atomic nitrogen (N) can be used for sympathetic cooling of open-shell molecules, and demonstrated co-trapping of N atoms with NH molecules \cite{Hummon:2008rm}, but the presence of $^3$He buffer gas at 0.5 K precluded measurement of N+NH collisions in that experiment.

In this Letter, we report the observation of cold (570 mK) collisions between magnetically trapped, ground state open-shell atoms and ground-state polar molecules, N and NH, both species spin polarized. We measure a small N+NH trap loss rate coefficient of $k^{(\mathrm{N+NH})}_\mathrm{loss} = 8(4) \times 10^{-13}$~cm$^{3}$s$^{-1}$. To interpret experimental observations, we carry out accurate {\it ab initio} and quantum scattering calculations of N+NH collisions in a magnetic field. Our theoretical results agree with the observations and predict high elastic-to-inelastic ratios ($\gamma$) over a wide range of temperatures (from 1K down to about 1 mK), suggesting the feasibility of sympathetic cooling of NH molecules by collisions with N atoms in a magnetic trap. Recent unpublished theory work by another group has made similar theoretical claims \cite{Zuchowski:2010bh}.

The apparatus used to co-trap atomic nitrogen and NH is similar to that described in our previous work \cite{Hummon:2008rm, tsikata:10nh}.  A pair of super-conducting solenoids produce a 4~T deep spherical quadrupole magnetic trap.  In the bore of the solenoids resides a cryogenic buffer gas cell held at 500~mK.
N and NH are produced in a
molecular beam using a DC glow discharge and enter the trapping region through a 1~cm diameter aperture in the buffer gas cell.
N and NH
thermalize
with $^3$He buffer gas
to 500~mK  and fall into the trap.
The initial buffer gas density at time of loading of $10^{15}$~cm$^{-3}$
is produced using a cryogenic reservoir with fast actuating valve \cite{tsikata:10nh} to inject buffer gas through a 3.8~cm diameter aperture in the cell.  After trap loading, the valve is closed, and the buffer gas exits primarily back out through the larger aperture, yielding background helium densities of $10^{12}$~cm$^{-3}$, corresponding to NH trap lifetimes of several seconds.

Detection of trapped N is performed using two-photon absorption laser induced fluorescence (TALIF)  from the ground $(2p^3)^4$S$_{3/2}$ state to the excited $(3p)^4$S$_{3/2}$ state at 96750 cm$^{-1}$.
NH is detected via laser induced fluorescence, as described in \cite{Hummon:2008rm}.
We initially load NH with densities on the order of $10^8$~cm\ensuremath{^{-3}}. Accurate knowledge of the absolute NH density is unnecessary as the NH density is much smaller than the N density in all our measurements.
To determine the N density, the TALIF signal is calibrated using N+N collisional loss measurements \cite{Tscherbul:2010nn}.
We load $5\times10^{11}$ N atoms into the trap with peak densities of $10^{12}$~cm$^{-3}$.

\begin{figure}
\epsfig{figure=Fig1_Hummon_PRL,width = 3 in}
\caption[NH decay. ]{NH trap decay taken at a trap depth of 3.9~T and cell temperature 570~mK.  The solid line is a fit to the form $n_\mathrm{NH} = a\exp[-\Gamma_\mathrm{NH}t]$.}
\label{fig:nh_decay}
\end{figure}

Figure \ref{fig:nh_decay} shows a typical NH trapping decay. At time $t=0$~s, N atoms and NH molecules are co-loaded into the magnetic trap.  For the first 2~s the NH trap loss is rapid due to the collisions with background helium gas during the pump out of buffer gas from the cell.  By $t=5$~s the NH trap loss reaches a steady rate, with typical $1/e$ lifetimes of about 3~s. Since the excitation laser for the atomic N TALIF detection causes additional loss of the trapped NH,  we measure the trapped atomic N density at  $t=15$~s.

\begin{figure}
\epsfig{figure=Fig2_Hummon_PRL,width = 3in}
\caption[NH loss vs N density]{NH loss rate vs cotrapped N density.  The data are fitted to Eq.~(1), yielding a value of  $k^{(\mathrm{N+NH})}_\mathrm{loss} = 8(4) \times 10^{-13}$~cm$^{3}$s$^{-1}$.}
\label{fig:NHvsN}
\end{figure}

To measure N+NH collisions, we observe the NH trap loss over a range of co-trapped N densities.  The NH trap loss is fit between $t=5$ and 10~s to an exponential decay to extract the total NH trap loss rate, $\Gamma_\mathrm{NH}$.  The cotrapped N density is varied by changing the ratio of molecular beam process gases N$_2$ and H$_2$ between (90\%, 10\%) to (3\%, 97\%).  Figure \ref{fig:NHvsN} shows the total NH loss rate, $\Gamma_\mathrm{NH}$, versus cotrapped nitrogen density.  The solid line in Fig. \ref{fig:NHvsN} is a fit to the equation
\begin{align}
\Gamma_\mathrm{NH} = \frac{k^{(\mathrm{N+NH})}_\mathrm{loss}}{14} n_\mathrm{N} + \Gamma_\mathrm{He}
\end{align}
where $n_\mathrm{N}$ is the peak nitrogen density, $\Gamma_\mathrm{He}$ is the NH loss rate attributable to collisions with background helium gas, and $ k^{(\mathrm{N+NH})}_\mathrm{loss}$ is the N+NH loss rate coefficient.  The factor of 14 arises from averaging the N and NH densities over the volume of the magnetic trap.  We find from this fit that $k^{(\mathrm{N+NH})}_\mathrm{loss} = 8(4) \times 10^{-13}$~cm$^{3}$s$^{-1}$. The uncertainty in $k^{(\mathrm{N+NH})}_\mathrm{loss}$ is dominated by the uncertainty in our nitrogen density calibrations \cite{Tscherbul:2010nn}.  The N+NH loss rate coefficient has contributions from both elastic (evaporative) and inelastic N+NH collisions \cite{Tscherbul:2010nn}, but their individual contributions cannot be determined by a single measurement.  Therefore, our measurement of $ k^{(\mathrm{N+NH})}_\mathrm{loss}$ is an upper limit on the inelastic N+NH rate coefficient and could be lower by up to a factor of 50\%.

To interpret the experimental observations and explore the possibility of sympathetic cooling of NH molecules by collisions with co-trapped N atoms, we performed rigorous quantum scattering calculations of inelastic relaxation in N+NH collisions in an external magnetic field. The Hamiltonian of the atom-molecule collision complex may be written ($\hbar=1$)
\begin{multline}\label{H}
\hat{H} = -\frac{1}{2\mu R}\frac{\partial^2}{\partial R^2}R + \sum_{S=1/2}^{5/2}V_S(R,r,\theta) |SM_S\rangle \langle SM_S| \\ + \frac{\hat{\ell}^2}{2\mu R^2} - \sqrt{\frac{24\pi}{5}} \frac{\alpha^2}{R^3}  \sum_{q}Y^\star_{2q}(\hat{R})[\hat{S}_\text{NH}\otimes\hat{S}_\text{N}]^{(2)}_q + \hat{H}_\text{NH} + \hat{H}_\text{N},
\end{multline}
where $\mu$ is the reduced mass of the complex, $R=|\bm{R}|$ is the N+NH separation vector, $r$ is the internuclear distance in NH, $\theta$ is the angle between the vectors $\bm{R}$ and $\bm{r}$, $V_S(R,r,\theta)$ is the interaction potential of the N+NH collision complex, $\hat{S}$ is the total spin of the complex, and $M_S$ is its projection on the magnetic field axis. The last two terms in Eq. (\ref{H}) describe non-interacting collision partners in the presence of an external magnetic field of strength $B$ \cite{Krems:2004jcp}. The term proportional to $R^{-3}$ represents the magnetic dipole interaction \cite{Tscherbul:2010nn, Krems:2004jcp}.

\begin{figure}[t!]
\epsfig{figure=Fig3_Hummon_PRL,width = 3 in}
\caption[theory]{Contour plot of the {\it ab initio} interaction PES for the fully spin-polarized $S = 5/2$ state of N+NH calculated in this work. The NH bond distance is fixed at its equilibrium value $r=1.958$ $a_0$. Energies are in units of cm$^{-1}$.}
\label{fig:PES}
\end{figure}

The interaction of NH($^3\Sigma$) molecules with N($^4S_{3/2}$) atoms gives rise to three adiabatic potential energy surfaces (PESs) with $S=1/2$, 3/2, and 5/2. Since in our experiments both atoms and molecules are confined in a permanent magnetic trap, the incident collision channel is the maximally spin-stretched Zeeman state $|M_{S_\text{NH}}=1\rangle\otimes|M_{S_\text{N}}=3/2\rangle$ with $S=5/2$. We neglect the weak couplings between the electronic states of different $S$ arising from the fine-structure terms in $\hat{H}_\text{NH}$   and off-diagonal matrix elements of the magnetic dipole interaction  \cite{Tscherbul:2010nn, Krems:2004jcp}.

To evaluate the PES for the $S=5/2$ electronic state of N-NH, we use the state-of-the-art partially spin restricted coupled cluster method with single, double, and non-iterative triple excitations [RCCSD(T)]  \cite{MOLPRO} using quadruple-zeta basis set (aug-cc-pvqz) of Dunning {\it et al.} \cite{Dunning89, Dunning92,Peterson94} augmented with $3s3p2d2f1g$ bond functions placed in the middle of the intermolecular distance. A contour plot of the calculated PES in shown in Fig. 3. The PES has a global minimum
87.83 cm$^{-1}$ deep in the linear N--HN configuration ($R=7.02a_0,\, \theta = 0$) and a secondary minimum of 77.52 cm$^{-1}$ in the HN--N configuration ($R=6.61a_0,\, \theta=180^\text{o}$) separated by a barrier 39.6 cm$^{-1}$ high located at $\theta=92^\circ$. Compared to the PES for the He-NH interaction \cite{Cybulski2005jcp}, our calculated N+NH PES is both deeper and more anisotropic.

To solve the scattering problem, we expand the wave function of the collision complex in a direct-product basis set $|NM_{N}\rangle |S_\text{NH}M_{S_\text{NH}}\rangle |S_\text{N}M_{S_\text{N}}\rangle |\ell m_\ell\rangle$, where $N$ is the rotational angular momentum of NH, and $M_N$, $M_{S_\text{NH}}$, $M_{S_\text{N}}$, and $m_\ell$ are the projections of $\hat{N}$, $\hat{S}_\text{NH}$, $\hat{S}_\text{N}$, and $\hat{\ell}$ on the magnetic field axis. This expansion results in a system of close-coupled (CC) equations parametrized by the matrix elements of the Hamiltonian (\ref{H}), which can be readily evaluated analytically given the direct-product structure of the basis \cite{Krems:2004jcp}. We solve the CC equations numerically for each value of the total angular momentum projection $M=M_N+M_{S_\text{NH}} + M_{S_\text{N}} + m_\ell$ for collision energies between $10^{-4}$ and 1 cm$^{-1}$ and  extract the $S$-matrix elements and scattering cross sections. Large basis sets with $N=0-5$ and $\ell=0-8$ are used to ensure that the results are converged to $<5\%$, resulting in 2906 CC equations for $M=0$.

\begin{figure}[t!]
\epsfig{figure=Fig4_Hummon_PRL,width = 2.65 in}
\caption[theory]{(a) Cross sections for elastic scattering (circles) and inelastic relaxation (triangles) in NH + N collisions calculated as functions of collision energy at a magnetic field of 0.1 T. Also shown is the cross section calculated with the magnetic dipole interaction omitted (dashed line). (b) Thermally averaged ratios of the rate constants for elastic scattering and inelastic relaxation as functions of temperature for $B=0.1$ T (triangles) and $B=1$ T (diamonds).}
\label{fig:theory}
\end{figure}

Figure 4(a) shows the cross sections for elastic energy transfer and inelastic relaxation in N+NH collisions. The cross sections increase with decreasing collision energy before reaching a maximum at $E_C\sim 2$ mK, which we identify as a shape resonance supported by the centrifugal barrier with $\ell=1$ in the incident collision channel. At $E_C< 0.1$ mK the inelastic cross sections assume the characteristic $1/E_C$ dependence on collision energy, and the elastic cross sections become constant, according to the Wigner threshold law for $s$-wave scattering. By thermally averaging the cross sections in Fig. 4(a), we obtain elastic and inelastic rate constants at the experimental temperature of 570 mK of $k_\mathrm{el} = 2.2\times10^{-10}$~cm$^3$s$^{-1}$, $k_\mathrm{in} = 4.1\times10^{-13}$~cm$^3$s$^{-1}$.  The fraction of elastic N-NH collisions that lead to trap loss is on the order of $3\times 10^{-3}$, leading to a total trap loss rate coefficient $k^{(\mathrm{N+NH})}_\mathrm{loss-theory} = 11\times10^{-13}$~cm$^3$s$^{-1}$, in good agreement with the experimental value of  $k^{(\mathrm{N+NH})}_\mathrm{loss-experiment} = (8\pm4)\times10^{-13}$~cm$^3$s$^{-1}$.
The ratio $\gamma = k_\mathrm{el}/k_\mathrm{in}$ is shown in Fig. 4(b) as a function of temperature. The ratio remains large ($\gamma>100$) over the temperature range $\sim$10 mK - 1 K, which indicates that NH molecules can be sympathetically cooled by elastic collisions with spin-polarized N atoms down to the milli-Kelvin regime. As shown in Fig. 4(b), an increasing magnetic field suppresses inelastic relaxation, so that $\gamma$ remains high ($\sim$50) even at $T= 1$ mK. This observation indicates that applying a strong uniform magnetic field of order $1$ T may be used to stabilize spin-polarized atom-molecule mixtures against collisional losses, thereby enhancing the efficiency of sympathetic cooling.

Inelastic collisions of $^3\Sigma$ molecules with open-shell atoms like N can occur due to (i) indirect couplings induced by the anisotropy of the atom-molecule interaction potential and the spin-spin interaction between the $N = 0$ and $N = 2$ rotational states \cite{Krems:2004jcp, campbell2009mechanism}, and (ii) direct couplings between atomic and molecular Zeeman levels induced by the long-range magnetic dipole interaction \ref{H}). The magnetic dipole interaction couples the Zeeman levels directly \cite{Tscherbul:2010nn}, and can thus be more efficient in inducing inelastic relaxation than the indirect couplings. Figure~4(a) shows that omitting the magnetic dipole interaction from scattering calculations reduces the inelastic cross sections by a factor of $\sim$10, confirming that Zeeman transitions in N+NH collisions are indeed driven by the magnetic dipole interaction.  In our experiments, both N and NH collision partners are fully spin-polarized, so the chemical reaction N + NH $\to$ N$_2$ + H  is spin-forbidden and can only proceed via non-adiabatic transitions between different electronic states of the N-NH complex in the entrance reaction channel mediated by the fine-structure and magnetic dipole couplings (see Eq. 2). Our observed value for $k^\text{(N+NH)}_\text{loss-experiment}$ is much smaller than the calculated reaction rate for spin-unpolarized reactants ($k_\text{reaction} \sim 3\times 10^{-11}$ cm$^3$/s at $T=1$ K) \cite{Frankcombe_N_NH}. Since our measurements are consistent with theoretical predictions which do not account for the reaction channel, we conclude that chemical exchange processes in spin-polarized N+NH mixtures occur at a slow rate, and do not contribute to the observed trap loss dynamics. This important finding shows that inelastic relaxation in N+NH collisions occurs via the same mechanism as dipolar relaxation in spin-polarized atomic gases \cite{Tscherbul:2010nn}.

In conclusion, we have measured a small N+NH trap loss rate coefficient of $k^{(\mathrm{N+NH})}_\mathrm{loss} = 8(4) \times 10^{-13}$~cm$^{3}$s$^{-1}$ at a temperature of $~\sim 600$ mK. To interpret experimental observations, we have carried out accurate {\it ab initio} quantum scattering calculations of Zeeman relaxation in N+NH collisions in a magnetic field and find theory and experiment to agree. Our calculations show that the ratio of N+NH elastic to inelastic collisions remains large ($>100$) over the temperature range $\sim$10 mK - 1 K, which indicates that it may be possible to sympathetically cool NH molecules down to the milli-Kelvin regime via elastic collisions with spin-polarized N atoms. It remains to be seen whether this conclusion holds for other paramagnetic molecules of interest such as the highly polar CaH or SrF in their electronic ground states of $^2\Sigma$ symmetry. If it does, it may be possible to create large samples of these molecules via collisional cooling with N atoms in a magnetic trap.

This work was supported by the Department of Energy,  under Grant No. DE-FG02-02ER15316 and by the Air Force Office of Scientific Research, under Grant No. FA9550-07-1-0492.