text:
# Semi-classical theory of collisional depolarization of spectral lines by atomic hydrogen. I: application to 𝑝 states of neutral atoms.

## 1 Introduction

Observations of linear polarization at the limb of the Sun produce rich structures. This polarization spectrum (the so-called “second solar spectrum”) is an entirely new spectrum to be explored and interpreted (Stenfo & Keller Stenflo (1997)). This polarization is due to scattering of the incident anisotropic solar radiation field. The radiative anisotropy is due to the structure of the semi-infinite or layer atmosphere in the vicinity of the surface, and is revealed by the center-to-limb darkening of the solar brightness.

In order to interpret quantitatively the Stokes parameters of the observed lines, a non-optically-thin and non-LTE model in the presence of a magnetic field has to be used. In fact, collision induced transitions between Zeeman sublevels can play a major role in depolarizing the levels (and thus the lines). Depolarizing collisions of the radiating and absorbing atoms with neutral hydrogen atoms of the medium have to be included in the statistical equilibrium equations. Very few depolarization rates have currently been computed, they have been obtained through sophisticated quantum chemistry methods which are accurate but cumbersome. It would be useful to develop methods capable of giving reasonable results for many levels rapidly. This is the object of the present work.

During the 1990’s Anstee, Barklem and O’Mara (Anstee Anstee2 (1992); Anstee & O’Mara Anstee1 (1991), Anstee3 (1995); Anstee et al. Anstee4 (1997); Barklem Barklem3 (1998); Barklem & O’Mara Barklem1 (1997); Barklem et al. Barklem2 (1998)) developed a new semi-classical theory, in order to obtain spectral line widths due to collisions with atomic hydrogen that were better than those obtained by the Van der Waals approximation. In this theory, which will be referred as ABO in the following, the scattering $`S`$ matrix is obtained by solving semi-classical coupled differential equations. The hydrogen-atom interaction potential is derived from time-independent second-order perturbation theory without exchange, allowing the Lindholm-Foley average over $`m`$ states (Brueckner Brueckner (1971); O’Mara O'Mara (1976)) to be removed. Whereas the Van der Waals potential underestimates the broadening by typically a factor two or more, the ABO method gives rather good agreement ($`<`$ 20 % for the spectral line widths) with the results obtained from quantum chemistry calculations. In addition, these authors used their ABO method by applying it to a selection of strong solar lines. The derived abundances are consistent with meteoritic ones.

In the present paper, we extend the ABO theory to the calculation of collisional depolarization rates which enter the statistical equilibrium equations for the polarized atomic density matrix. This paper is focused on the method and is applied to $`p`$ states of neutral atoms. We will present and compare our results with those recently obtained by quantum chemistry methods. We also compare the results to the depolarization rates calculated with the Van der Waals interaction. A discussion of our results is given, studying the region of the potential curves which play the major role in the variation of the depolarization cross sections. In conclusion, we show that our method, owing to the simplicity of the processes taken into account, is a useful alternative to the accurate but time consuming quantum chemistry method for obtaining most of the depolarization rates to better than 20 % accuracy. It should now be possible to rapidly obtain the large amount of data needed for the interpretation of the second solar spectrum, in particular for the heavy atoms that are inaccessible to present quantum chemistry calculations. An extension to $`d`$ and $`f`$ states is in progress. An extension to ions is in progress and this will be the subject of further papers.

## 2 Formulation of the problem

### 2.1 The statistical equilibrium equations (SEE) for polarization studies

The calculation of the Stokes parameters of a spectral line observed at the solar limb is a non-LTE (Local Thermodynamical Equlibrium) problem. It needs the solution of the coupling between the polarized radiative transfer and the statistical equilibrium equations (SEE) for the populations and coherences of the polarized atomic density matrix. In fact the atomic levels are only “aligned”, which means that the observed lines are globally linearly polarized. This linear polarization is due to anisotropic scattering of the incident solar radiation field which contributes to excitation of the atomic levels, the anisotropy being due to limb darkening, even at zero altitude above the limb.

Sahal-Bréchot (Sahal2 (1977)) developed the SEE for an “aligned” multilevel atom, i.e. for the Zeeman sublevels $`|\alpha JM_J`$. She took into account excitation by the anisotropic solar radiation field, deexcitation by spontaneous and induced emission, and excitation, deexcitation, polarization transfer and depolarization by isotropic collisions of the particles of the medium. She included the effect of a “strong” magnetic field (strong case of the Hanle effect) and thus coherences between Zeeman sublevels could be ignored. Only the populations of the Zeeman sublevels were taken into account.

She wrote these equations both in the dyadic basis $`|\alpha JM_J\alpha ^{}J^{}M_J^{}|`$ built on the standard basis of the atomic states $`|\alpha JM_J`$ and in the irreducible tensorial operator basis $`{}_{}{}^{\alpha J\alpha ^{}J^{}}T_{q}^{k}`$. The expansion of the $`{}_{}{}^{\alpha J\alpha ^{}J^{}}T_{q}^{k}`$ over the dyadic basis is given by Fano & Racah (Fano1 (1959)) (cf. also Messiah Messiah (1961); Fano Fano2 (1963)):
$`{}_{}{}^{\alpha J\alpha ^{}J^{}}T_{q}^{k}=_{M_JM_J^{}}(1)^{J^{}M_J^{}}`$

$`\times |\alpha JM_J\alpha ^{}J^{}M_J^{}|JJ^{}M_JM_J^{}|kq,`$ (1)
where $`JJ^{}M_JM_J^{}|kq`$ is a Clebsch-Gordan coefficient. She demonstrated the interest and power of this second basis, which is especially well adapted to problems with spherical and/or cylindrical symmetry processes. The equations are simplified, and the physics is extracted and highlighted. In particular, the odd $`k`$-terms can be eliminated from the SEE when the incident radiation field is anisotropic but unpolarized. In fact the natural radiation field transports $`k=0`$ (population) and $`k=2`$ (alignment) terms only. Orientation terms ($`k=1`$) vanish in the present astrophysical problem.

Sahal-Bréchot (Sahal2 (1977)) gave the general expressions for the various transition probabilities due to isotropic and anisotropic radiation and for isotropic collisions, i.e. expressions for population, orientation and alignment transfer between the levels, creation and destruction of population, orientation and alignment (depolarization rates in particular). At the same time, Omont (Omont2 (1977)) developed a general theory of relaxation for a two-level atom. In addition he introduced coherences.

In fact, for isotropic processes, the coupling terms in the SEE implying transfer of coherence $`qq^{}`$ are zero, and due to the isotropy of the collisions with hydrogen, the depolarization rates are $`q`$-independent. This is one of the advantages of the irreducible tensorial operator basis. Therefore, as we will consider isotropic collisions only in the following, it will be sufficient to calculate collisional cross sections between the sublevels, as for the no-coherence case (Sahal-Bréchot Sahal2 (1977)), even if coherences play a role in the SEE. We recall that coherences can appear when the magnetic field is weak and not directed along the solar radius (Hanle effect). We also notice that in the present paper we are interested in hydrogen-atom collisions and not in electron or proton-ion collisions as in Sahal-Bréchot (Sahal2 (1977)). After integration of the cross sections over a Maxwell distribution of velocities and multiplying by the local hydrogen density, the depolarization rates obtained (and eventually polarization transfer rates between levels) can be used in the statistical equilibrium equations.

We focus on the collisional depolarization rates in the present paper. We will use both the dyadic and the irreducible tensorial operator bases. The latter one is in fact more convenient.

An important preliminary point to notice is the fact that the radiating atom can have a nuclear spin and thus hyperfine structure. The importance of the effect of the hyperfine structure on the Stokes parameters of the line studied was discussed for the example of the hydrogen coronal line Ly$`\alpha `$ by Bommier & Sahal-Bréchot (Bommier2 (1982)). They showed that the SEE must be solved for the hyperfine levels when the inverse of the lifetime of the upper level is smaller than the hyperfine splitting, i.e. the hyperfine levels are separated.

This arises from the Heisenberg uncertainty principle. On the contrary, when the hyperfine splitting is small compared to the inverse of the lifetime, the hyperfine levels are mixed and degenerate, the hyperfine structure can be ignored and the SEE can be solved for the $`|\alpha J`$ fine structure levels. This was the case for the hydrogen coronal line where the $`F`$ levels are degenerate and the $`J`$ levels are well separated.

In the second solar spectrum, where metallic lines are observed, the inverse of the lifetime (9.7 MHz for the Na D lines) (Gaupp et al Gaupp (1982)) is often small compared to the hyperfine splitting (if it exists) (cf. Table 1). Then the SEE must be solved for the $`|\alpha JF`$ hyperfine levels (Bommier & Molodij Bommier3 (2002); Kerkeni & Bommier Kerkeni3 (2002)) and collisional rates between the hyperfine sublevels must be calculated.

### 2.2 Dyadic basis for the depolarization rates

The depolarization arises from isotropic collisions between the Zeeman sublevels of a level $`|\alpha J`$: the collisions reestablish thermodynamical equlibrium between the sublevels and thus atomic polarization vanishes. We have to calculate (Sahal-Bréchot Sahal2 (1977)) for each sublevel $`|\alpha JM_J`$ (or $`|\alpha JFM_F`$ if there is an hyperfine structure):

$`D(\alpha JM_J,T)={\displaystyle \underset{M_J^{}M_J}{}}\zeta (\alpha JM_J\alpha JM_J^{},T)`$ (2)
where $`\zeta (\alpha JM_J\alpha JM_J^{},T)`$ is the collisional transition rate between the sublevels $`|\alpha JM_J|\alpha JM_J^{}`$. It can be written as a function of the local temperature $`T`$ and the hydrogen perturber local density $`n_H`$:

$`\zeta (\alpha JM_J\alpha JM_J^{},T)=n_H`$ $`{\displaystyle _0^{\mathrm{}}}\sigma (\alpha JM_J\alpha JM_J^{},v)`$ (3)
$`\times vf(v)dv,`$
$`f(v)`$ being the distribution function for the atom-perturber relative velocity $`v`$, and $`\sigma (\alpha JM_J\alpha JM_J^{},v)`$ is the cross section. For the Maxwellian distribution,

$`f(v)=\sqrt{{\displaystyle \frac{2}{\pi }}}({\displaystyle \frac{\mu }{kT}})^{\frac{3}{2}}v^2\text{exp}({\displaystyle \frac{\mu v^2}{2kT}}),`$ (4)
where $`\mu `$ is the reduced mass of the system.

The relaxation term for the density matrix (diagonal term of the system of SEE) due to collisions is given by:

$`({\displaystyle \frac{d^{\alpha J}\rho _{M_JM_J}}{dt}})_{coll}`$ $`=`$ $`D(\alpha JM_J,T)^{\alpha J}\rho _{M_JM_J}`$
$`{\displaystyle \underset{\alpha ^{}J^{}\alpha J}{}}{\displaystyle \underset{M_J^{}}{}}\zeta (\alpha JM_J\alpha J^{}M_J^{},T)`$
$`\times ^{\alpha J,\alpha ^{}J^{}}\rho _{M_JM_J^{}}`$
The terms added to the depolarization term $`D(\alpha JM_J,T)`$ are the so-called quenching terms.

The dyadic basis is not suitable because the depolarization terms depend on the chosen dyadic basis. The $`{}_{}{}^{\alpha J,\alpha ^{}J^{}}T_{q}^{k}`$ basis is well adapted because the depolarization terms do not depend on the basis and have a physical interpretation.

### 2.3 Basis of irreducible tensorial operators

In this basis, the SEE relaxation term in the non-coherence case ($`q=0`$) (Sahal-Bréchot Sahal2 (1977); Omont Omont2 (1977)) becomes

$`({\displaystyle \frac{d^{\alpha J}\rho _0^k}{dt}})_{coll}=D^k(\alpha J,T)^{\alpha J\alpha J}\rho _0^k+\text{quenching term,}`$ (6)
where $`{}_{}{}^{\alpha J}\rho _{0}^{k}`$ is the multipole $`k`$-moment of the density matrix for the level $`|\alpha J`$. For isotropic collisions $`D^k(\alpha J,T)`$ is characterized by $`(2J+1)`$ constants ($`0k2J`$) and $`D^0(\alpha J,T)`$ is called the destruction rate of population, which is zero for the no-quenching approximation. $`D^1(\alpha J,T)`$ is the destruction rate of orientation (circular atomic polarization). $`D^2(\alpha J,T)`$ is the destruction rate of alignment which is of interest in astrophysics (linear polarization), and the depolarization rate $`D^k(\alpha J,T)`$ is given by:

$`D^k(\alpha J,T)=n_H{\displaystyle _0^{\mathrm{}}}\sigma ^k(\alpha J,v)vf(v)𝑑v`$ (7)
Using our present notation, formula (20) in Sahal-Bréchot (Sahal2 (1977)) with $`k=k^{}`$ can be written as:

$`\sigma ^k(\alpha J,v)`$ $`=`$ $`(2k+1){\displaystyle \underset{M_J^{}}{}}{\displaystyle \underset{M_JM_J^{}}{}}(1)^{2J2M_J}`$
$`\times \left(\begin{array}{ccc}J& k& J\\ M_J& 0& M_J\end{array}\right)^2\sigma (\alpha JM_J\alpha JM_J^{},v)`$
$`(2k+1){\displaystyle \underset{M_J}{}}{\displaystyle \underset{M_J^{}M_J}{}}(1)^{2JM_JM_J^{}}`$
$`\times \left(\begin{array}{ccc}J& k& J\\ M_J& 0& M_J\end{array}\right)\left(\begin{array}{ccc}J& k& J\\ M_J^{}& 0& M_J^{}\end{array}\right)`$
$`\times \sigma (\alpha JM_J\alpha JM_J^{},v)`$
The expressions between parentheses denote $`3j`$-coefficients (Messiah Messiah (1961)). After some calculations which are not detailed here, this equation can also be written as:

$`\sigma ^k(\alpha J,v)`$ $`=`$ $`\sigma (\alpha J,v)(2k+1){\displaystyle \underset{M_J,M_J^{}}{}}(1)^{2JM_JM_J^{}}`$
$`\times \left(\begin{array}{ccc}J& k& J\\ M_J& 0& M_J\end{array}\right)\left(\begin{array}{ccc}J& k& J\\ M_J^{}& 0& M_J^{}\end{array}\right)`$
$`\times \sigma (\alpha JM_J\alpha JM_J^{},v).`$
Where $`\sigma (\alpha J,v)`$ is the elastic cross section for the level $`|\alpha J`$ and the relative velocity $`v`$:

$`\sigma (\alpha J,v)={\displaystyle \underset{M_J^{}}{}}\sigma (\alpha JM_J\alpha JM_J^{},v)`$ (22)
This cross section can also be written as an average over the initial states $`M_J`$ as usual:

$`\sigma (\alpha J,v)={\displaystyle \frac{1}{2J+1}}{\displaystyle \underset{M_J,M_J^{}}{}}\sigma (\alpha JM_J\alpha JM_J^{},v)`$ (23)
where the angular average over all the directions of the collisions is not necessary here because the averaging over $`M_J^{}`$ is equivalent to the angular average over the directions of the collision. This property will be used in Sect. 6, first line in equation (6).

Since the depolarization cross section $`\sigma ^k(\alpha J,v)`$ can be written as a linear combination of the $`\sigma (\alpha JM_J\alpha JM_J^{},v)`$, the depolarization rate $`D^k(\alpha J,v)`$ is a linear combination of the collisional rates between the Zeeman sublevels $`\zeta (\alpha JM\alpha JM^{})`$ (Sahal-Bréchot Sahal2 (1977)). In particular we have (Sahal-Bréchot Sahal2 (1977)), for $`J=1`$:

$`D^1(\alpha 1)`$ $`=`$ $`\zeta (101\pm 1)+2\zeta (111\pm 1)`$
$`D^2(\alpha 1)`$ $`=`$ $`3\zeta (101\pm 1),`$
and for $`J=3/2`$ (Gordeyev et al Gordeyev (1969); Roueff & Suzor Roueff2 (1974)):

$`D^1(\alpha {\displaystyle \frac{3}{2}})`$ $`=`$ $`{\displaystyle \frac{1}{5}}\zeta ({\displaystyle \frac{3}{2}}{\displaystyle \frac{1}{2}}{\displaystyle \frac{3}{2}}{\displaystyle \frac{1}{2}})+{\displaystyle \frac{2}{5}}\zeta ({\displaystyle \frac{3}{2}}{\displaystyle \frac{3}{2}}{\displaystyle \frac{3}{2}}{\displaystyle \frac{1}{2}})`$
$`+`$ $`{\displaystyle \frac{8}{5}}\zeta ({\displaystyle \frac{3}{2}}{\displaystyle \frac{3}{2}}{\displaystyle \frac{3}{2}}{\displaystyle \frac{1}{2}})+{\displaystyle \frac{9}{5}}\zeta ({\displaystyle \frac{3}{2}}{\displaystyle \frac{3}{2}}{\displaystyle \frac{3}{2}}{\displaystyle \frac{3}{2}})`$
$`D^2(\alpha {\displaystyle \frac{3}{2}})`$ $`=`$ $`2\zeta ({\displaystyle \frac{3}{2}}{\displaystyle \frac{3}{2}}{\displaystyle \frac{3}{2}}{\displaystyle \frac{1}{2}})+2\zeta ({\displaystyle \frac{3}{2}}{\displaystyle \frac{1}{2}}{\displaystyle \frac{3}{2}}{\displaystyle \frac{1}{2}})`$
$`+`$ $`2\zeta ({\displaystyle \frac{3}{2}}{\displaystyle \frac{3}{2}}{\displaystyle \frac{3}{2}}{\displaystyle \frac{1}{2}})`$
$`D^3(\alpha {\displaystyle \frac{3}{2}})`$ $`=`$ $`{\displaystyle \frac{9}{5}}\zeta ({\displaystyle \frac{3}{2}}{\displaystyle \frac{1}{2}}{\displaystyle \frac{3}{2}}{\displaystyle \frac{1}{2}})+{\displaystyle \frac{8}{5}}\zeta ({\displaystyle \frac{3}{2}}{\displaystyle \frac{3}{2}}{\displaystyle \frac{3}{2}}{\displaystyle \frac{1}{2}})`$
$`+`$ $`{\displaystyle \frac{2}{5}}\zeta ({\displaystyle \frac{3}{2}}{\displaystyle \frac{3}{2}}{\displaystyle \frac{3}{2}}{\displaystyle \frac{1}{2}})+{\displaystyle \frac{1}{5}}\zeta ({\displaystyle \frac{3}{2}}{\displaystyle \frac{3}{2}}{\displaystyle \frac{3}{2}}{\displaystyle \frac{3}{2}})`$
## 3 Calculation of the cross sections

We will use a semi-classical treatment for the collision: we will assume that the hydrogen perturbing atom moves along a straight path, unperturbed by the collision process and characterized by an impact-parameter $`b`$.

In the dyadic basis $`|\alpha JM_J\alpha ^{}J^{}M_J^{}|`$ (or $`|\alpha JFM_F\alpha ^{}J^{}F^{}M_F^{}|`$ if there is an hyperfine structure), the cross section follows from an integration over the impact-parameter $`b`$ and averaging the transition probability $`P(\alpha JM_J\alpha JM_J^{},b,v)`$ over all the collision directions:

$`\sigma (\alpha JM_J\alpha JM_J^{},v)`$ $`=`$ $`{\displaystyle _0^{\mathrm{}}}2\pi bdb({\displaystyle \frac{1}{8\pi ^2}}`$
$`\times `$ $`{\displaystyle }P(\alpha JM_J\alpha JM_J^{},b,v)d\mathrm{\Omega })`$
and then the collisional transition rate is given by equation (3). The transition probability can be written in terms of the $`S`$-matrix elements,

$`P(\alpha JM_J\alpha JM_J^{},b,v)`$ $`=`$
$`|\alpha JM_J|IS(b,v)|\alpha JM_J^{}|^2`$
where $`I`$ is the unit matrix and $`T=IS`$ is the so-called transition matrix.

### 3.1 Calculation of the S-matrix

The Hamiltonian of the system atom-perturber can be written as:

$`H=H_A+H_P+V,`$ (26)
where $`H_A`$ is the atomic Hamiltonian, $`H_P`$ that of the perturber and $`V`$ is the interaction energy. $`H_A`$ can be written as:

$`H_A=H_0+H_{FS}+H_{HFS},`$ (27)
where $`H_0`$ contains the electrostatic part, $`H_{FS}`$ the fine structure and $`H_{HFS}`$ the hyperfine structure if it exists.

It can be shown that the hyperfine splitting $`\mathrm{\Delta }E_{HFS}`$ is always small when compared to $`V`$ and thus $`J`$ and $`I`$ (the nuclear spin) can be decoupled during the collision process. This arises from the Heisenberg uncertainty principle: if $`\tau `$ is a typical collision time, one can ignore the hyperfine structure if

$`\tau \mathrm{\Delta }E_{HFS}1,`$ (28)
where $`\tau {\displaystyle \frac{b}{v}}`$, $`b`$ being a typical impact-parameter (of the order of $`10a_0`$, $`a_0`$ being the Bohr radius) and $`v4.8\times 10^3\text{a.u.}`$ for a temperature of about 5000 K so that $`\tau 2.08\times 10^3\text{a.u.}`$. In the case of the Na atom the nuclear spin is $`I=3/2`$, which yields two hyperfine levels in the ground state $`3s^2S_{1/2}`$, $`F=1`$ and $`F=2`$. For the excited $`3p^2P_{3/2}`$ state there are four hyperfine levels: $`F=0`$, $`F=1`$, $`F=2`$ and $`F=3`$. The energy differences between these levels are given in Table 1. In the case of the ground state $`3s^2S_{1/2}`$, $`\mathrm{\Delta }E_{HFS}2.694\times 10^7\text{a.u.}`$ $`(1772\text{MHz})`$ (Ackermann Ackermann (1966)) and then $`\tau \mathrm{\Delta }E_{HFS}5.6\times 10^41`$. For the excited level the most important hyperfine splitting is, of course, between the levels $`F=0`$ and $`F=3`$ where $`\mathrm{\Delta }E_{HFS}1.653\times 10^8\text{a.u.}`$ $`(108.72\text{MHz})`$ (Gangrsky et al Gangrsky (1998)) and then $`\tau \mathrm{\Delta }E_{HFS}3.4\times 10^51`$.

The hyperfine splitting can always be ignored in the treatment of the collision under these conditions. We have $`F=I+J`$ and $`M_F=M_I+M_J`$, and

$`|\alpha JFM_F={\displaystyle \underset{M_I,M_J}{}}IJM_IM_J|FM_F|JM_J|IM_I,`$ (29)
where $`IJM_IM_I|FM_F`$ is a Clebsch-Gordan coefficient. Then the matrix element

$`\alpha JFM_F|S|\alpha ^{}J^{}F^{}M_F^{}={\displaystyle \underset{M_I,M_J,M_I^{},M_J^{}}{}}IJM_IM_J|FM_F`$
$`\times IJ^{}M_I^{}M_J^{}|F^{}M_F^{}\alpha JM_J|IM_I|S|IM_I^{}|\alpha ^{}J^{}M_J^{}.`$ (30)
If the hyperfine structure can be neglected during the collision, the $`S`$-matrix is diagonal in $`I`$ and its elements do not depend on $`M_I`$ ($`I`$ is conserved). Thus $`M_I^{}=M_I`$ and

$`\alpha JFM_F|S|\alpha ^{}J^{}F^{}M_F^{}={\displaystyle \underset{M_I,M_J,M_J^{}}{}}IJM_IM_J|FM_F`$
$`\times IJ^{}M_IM_J^{}|F^{}M_F^{}\alpha JM_J|S|\alpha ^{}J^{}M_J^{}`$ (31)
This relation holds for each impact-parameter and each relative velocity.

We recall the relation between Clebsch-Gordan coefficients and $`3j`$-coefficients:
$`j_1j_2m_1m_2|jm=\sqrt{2j+1}(1)^{j_1j_2+m}\left(\begin{array}{ccc}j_1& j_2& j\\ m_1& m_2& m\end{array}\right)`$
We use the convention of Messiah (Messiah (1961)) for the phase factors.

### 3.2 Calculation of the $`S`$-matrix elements between fine sublevels

The fine structure splitting can generally be neglected during the collision problem at the temperature of interest in the second solar spectrum ($`5000K`$). The validity of the condition is the same as above, except that the fine structure splitting is larger than the hyperfine splitting. For example in the case of the Na atom, at $`5000K`$, $`\mathrm{\Delta }E_{FS}0.78\times 10^4\text{a.u.}`$ and $`\tau 2.08\times 10^3\text{a.u}`$ and so $`\tau \mathrm{\Delta }E_{FS}0.16<1`$. This has been discussed in earlier papers focussed on collisional line broadening, see for example Roueff (Roueff2 (1974)) and related papers. She showed that the fine structure splitting is of importance at low temperatures (200-500 K) for the sodium $`D`$ lines and negligible at higher temperatures for line broadening studies (500 K and higher). The $`D_1`$ and $`D_2`$ lines have in fact the same collisional width for hydrogen perturbers at solar temperatures.

If the fine structure is neglected, where $`J=S+l`$, and $`M_J=M_S+m_l`$, we have:

$`\alpha lJM_J|S|\alpha ^{}l^{}J^{}M_J^{}={\displaystyle \underset{M_S,m_l,m_l^{}}{}}SlM_Sm_l|JM_J`$ (32)
$`\times Sl^{}M_Sm_l^{}|J^{}M_J^{}\alpha lm_l|S|\alpha ^{}l^{}m_l^{},`$
cf. equation (3.1).

However, we notice that the parallelism between line broadening and depolarization cross sections must be used carefully in the peculiar following case. Depolarization of hyperfine components of spherically symmetric $`J=J^{}=1/2`$ levels (ground state of sodium for instance), needs the matrix elements $`M_J=\pm \frac{1}{2}M_J^{}=\pm \frac{1}{2}`$. If the fine structure is neglected, they are zero and the depolarization and polarization transfer rates of the hyperfine components of the ground states are zero. Kerkeni et al (Kerkeni1 (2000)) calculated the hyperfine collisional depolarization rates for the ground state of sodium with a quantum chemistry method which takes into account the fine structure splitting during the interaction processes. They showed that these rates are of the same order of magnitude as those of upper levels for which the fine structure splitting plays a minor role in the collisional processes. In fact this remark does not concern $`l1`$ states and is very specific. However, we consider that it is important to mention it.

## 4 The model of Anstee, Barklem and O’Mara (ABO)

The present semi-classical model used for the depolarization is the same one as that already developed for line broadening theory applications in the 1990’s by Anstee, Barklem and O’Mara (Anstee Anstee2 (1992); Anstee & O’Mara Anstee1 (1991), Anstee3 (1995); Anstee et al. Anstee4 (1997); Barklem Barklem3 (1998); Barklem & O’Mara Barklem1 (1997), Barklem et al. Barklem2 (1998)).

The theory is described in these papers, and thus only a brief review is presented in the following:

1. The model is of semi-classical nature, the internal states of the two atoms are treated within the framework of quantum mechanics and the collision is modelled with the straight path classical approximation: it is assumed that the interparticle interaction is too weak to deflect the perturbed particle from the straight classical path.

2. The interaction energy between the two atoms is smaller than the energy eigenvalues of the isolated atoms, so that the Rayleigh-Schrödinger perturbation theory can be used for obtaining perturbed eigenvalues. Ionic and exchange interactions are neglected. The hyperfine and fine structure interactions are also neglected.

3. Time-independent perturbation theory is applied to second order. It is assumed that the energy denominator in the second order term can be replaced by a suitable average energy $`E_p`$. This is the Unsöld approximation (Unsöld Unsold1 (1927); Unsöld Unsold2 (1955)). The value $`E_p=`$ -4/9, in atomic units, is adopted. This approximation works when the energy level separations of the perturbed atom are small compared with those of hydrogen. For neutral atoms this condition is generally satisfied. This is particularly valid when considering interactions with hydrogen in the ground state, like in this work, which is well separated from any other levels (Anstee & O’Mara Anstee1 (1991)).

4. Unperturbed atomic eigenfunctions are used. The Coulomb approximation with quantum defect is used.

## 5 Treatment of the collision

Each collision is characterized by an impact-parameter vector $`b`$ and a relative velocity vector $`v`$. Since the hyperfine and fine structure interactions are neglected, equation (26) becomes:

$`H=H_0+H_P+V`$ (33)
where $`V`$ is the interaction energy which is assumed to be totally electrostatic. Only $`V`$ depends on time. In the coordinate system of Fig. 1, $`V`$ is given, in atomic units, by (Anstee & O’Mara Anstee1 (1991)):

$`V={\displaystyle \frac{1}{R}}+{\displaystyle \frac{1}{r_{12}}}{\displaystyle \frac{1}{r_2}}{\displaystyle \frac{1}{p_1}}`$ (34)
and atomic units are used hereafter.

Since the hydrogen perturber remains in its ground state $`1s`$ during the collision, the wave function $`|\psi `$ of the system (atom+perturber), is the product of the wave functions $`|\psi _{(A)}`$ of the perturbed atom and that of hydrogen in its ground state $`|100_{(H)}`$:

$`|\psi =|\psi _{(A)}|100_{(H)}.`$ (35)
The eigenenergy of the hydrogen perturber in its ground state is denoted by $`E_H^0`$. Following Anstee & O’Mara (Anstee1 (1991)), we use second-order time-independent perturbation theory for the eigenenergy of the system.

We also use the unperturbed atomic wave functions defined by:

$`H_0|nlm_l=E_{nlm_l}^0|nlm_l;|\psi _{(A)}=|nlm_l.`$ (36)
They form a complete set of eigenfunctions $`|nlm_l`$ and eigenvalues $`E_{nlm_l}^0`$. The perturbed eigenvalues of the system of Hamiltonian $`H`$ can be written as:

$`E_i=E_i^0+i|V|i+{\displaystyle \underset{ji}{}}{\displaystyle \frac{i|V|jj|V|i}{E_i^0E_j^0}},`$ (37)
the summation being over all product states $`|j=|\psi _{(A)}|\psi _{(H)}`$ except $`|i=|nlm_l_{(A)}|100_{(H)}`$, and where $`E_i^0=E_{nlm_l}^0+E_H^0`$.

Quenching is neglected and thus we consider only the subspace $`nl`$ (2$`l`$+1 states) and we denote the product state of the two separated atom at $`R=\mathrm{}`$ by $`|M_l`$ (Anstee & O’Mara Anstee1 (1991)),

$`|M_l=|i=|nlm_l_{(A)}|100_{(H)}.`$ (38)
We introduce an effective potential $`V_{eff}`$, the diagonal elements of which are defined as

$`M_l|V_{eff}|M_l=E_{M_l}E_{M_l}^0`$ (39)
The evaluation of the second order contribution to the above expression (37) is very difficult because it requires a summation over the complete set of eigenstates of the two isolated atoms. However, the expression can be greatly simplified by using the Unsöld approximation, i.e. by replacing $`E_i^0E_j^0`$ by the constant value $`E_p`$. With this approximation, the expression for $`V_{eff}`$ becomes

$`M_l\left|V_{eff}\right|M_l`$ $`=`$ $`M_l\left|V\right|M_l{\displaystyle \frac{1}{E_p}}(M_l\left|V\right|M_l)^2`$ (40)
$`+{\displaystyle \frac{1}{E_p}}M_l\left|V^2\right|M_l,`$
where $`M_l|V_{eff}|M_l`$ is the so-called Rayleigh-Schrödinger-Unsöld (RSU) potential which is used throughout the ABO model. Then, for computing the RSU potential, explicit expressions for the wave functions are needed.

The perturbed atom is modelled as a single valence electron outside a spherical ionic core. The chosen radial functions are of Coulomb type (Seaton Seaton (1958)). The quantum defect and thus the effective quantum number $`n^{}=[2(E_{\mathrm{}}E_{nl})]^{1/2}`$ is a fundamental ingredient of the calculation of the interaction potential used in the ABO model. We refer to the ABO papers for more details concerning the interaction potential.

The $`S`$-matrix requires the time-dependent Schrödinger equation to be solved, i.e.

$`(H_0+H_P+V_{eff})|\psi (t)=\text{i}{\displaystyle \frac{d|\psi (t)}{dt}}`$ (41)
where the interaction potential $`V`$ is now replaced by $`V_{eff}`$. The $`S`$-scattering matrix is defined by:

$`S=U(+\mathrm{},\mathrm{})`$ (42)
where $`U(+\mathrm{},\mathrm{})`$ is the evolution operator in the interaction representation, where

$`\stackrel{~}{\psi }(t+s)=U^{}(t+s,t)\psi (t)=\text{e}^{\text{i}(H_0+H_P)s}\psi (t),`$ (43)
is the corresponding wave function. The wave function is expressed in terms of the basis formed by the eigenvectors $`|M_l`$ of $`(H_0+H_P)`$ given by (38) so that:

$`|\psi (t)={\displaystyle \underset{M_l}{}}a_{M_l}(t)\text{e}^{\text{i}E_{M_l}^0t}|M_l;`$
$`|\stackrel{~}{\psi }(t)={\displaystyle \underset{M_l}{}}a_{M_l}(t)|M_l`$
and we obtain a system of $`(2l+1)`$ coupled linear differential equations of first order:

$`\text{i}{\displaystyle \frac{a_{M_l}(t)}{t}}={\displaystyle \underset{M_k=l}{\overset{l}{}}}a_{M_k}(t)M_l\left|V_{eff}\right|M_k.`$ (45)
We have to solve this system for the initial conditions:

$`a_{M_l}(\mathrm{})=\delta (M_i,M_l),`$ (46)
where $`\delta `$ is the Kronecker symbol and

$`M_l|S|M_i=M_i|S|M_l=a_{M_l}(+\mathrm{}).`$ (47)
Thus the system has to be solved $`(2l+1)`$ times for obtaining all the elements of the $`S`$-matrix for the $`(nl)`$ subspace.

For the calculations, two frames of reference are used. The first one is the atomic frame in which the perturbed atom is stationary at the origin, the quantization axis is taken as perpendicular to the collision plane $`(`$$`b`$,$`v`$$`)`$ and the second one is the rotating or molecular frame in which the interatomic axis is taken as the quantization axis. These reference frames are shown in Fig. 2. This choice of frames gives the simplest relationship between states quantized relative to the atomic and the rotating frames respectively (Roueff Roueff1 (1974)).

The interaction potential $`V`$ has a cylindrical symmetry about the internuclear axis $`(OZ)`$ and so $`V_{eff}`$ is diagonal in the rotating frame.

A transformation of the angular momentum states from the atomic basis to the rotating frame is achieved by using the rotations matrices (Messiah Messiah (1961)). The Euler angles are given by $`(\beta ,\frac{\pi }{2},\frac{\pi }{2})`$ and equation (45) becomes

$`\text{i}{\displaystyle \frac{a_{M_l}(t)}{t}}={\displaystyle \underset{M_kM_k^{}}{}}𝒟_{M_lM_k^{}}^{(l)}𝒟_{M_lM_k^{}}^{(l)}a_{M_k}(t)M_l\left|V_{eff}\right|M_l`$ (48)
when $`𝒟`$ is the rotation operator.

For a given quantum number $`l`$ we obtain $`(2l+1)`$ coupled differential equations. These coupled equations must be solved for obtaining the $`(2l+1)^2S`$-matrix elements.

The interaction matrix is Hermitian and the $`S`$-matrix is unitary and symmetric. If the quantization axis is perpendicular to the collision plane the selection rules are: $`\mathrm{\Delta }m_l=\pm 2,\pm 4\mathrm{}`$. This rule will disappear when we take an average over all possible orientations of the collision plane, since the collisions are isotropic.

## 6 Depolarization calculations for $`l=1`$

Let us consider collisions between an atom in an $`l=1`$ state and hydrogen in its ground state $`1s`$. The state of the system during the collision becomes explicitly:

$`|\psi (t)`$ $`=`$ $`a_1(t)|11|00\text{e}^{\text{i}E_1^0t}+a_0(t)|10|00\text{e}^{\text{i}E_0^0t}`$ (49)
$`+`$ $`a_1(t)|11|00\text{e}^{\text{i}E_1^0t}`$
and the coupled differential equations become (Anstee & O’Mara Anstee1 (1991)):

$`\text{i}{\displaystyle \frac{a_1(t)}{t}}`$ $`=`$ $`a_1(t)V_++\text{e}^{2\text{i}\beta }a_1(t)V_{}`$
$`\text{i}{\displaystyle \frac{a_0(t)}{t}}`$ $`=`$ $`a_0(t)1|V_{eff}|1`$ (50)
$`\text{i}{\displaystyle \frac{a_1(t)}{t}}`$ $`=`$ $`a_1(t)V_++\text{e}^{2\text{i}\beta }a_1(t)V_{}`$
where:

$`V_\pm ={\displaystyle \frac{1}{2}}(1|V_{eff}|1\pm 0|V_{eff}|0).`$ (51)
The equations are integrated with a Runge-Kutta-Merson routine and then we obtain the transition matrix elements in the $`|\alpha lm_l`$ basis for a given velocity and impact-parameter.

We recall that we need, for depolarization calculations, the $`T`$-matrix elements in the $`|\alpha JM_J`$ basis. Since the spin is neglected, the $`S`$-matrix is diagonal in $`S`$ and its elements do not depend on $`M_S`$. Following equation (32) we obtain:

$`\alpha JM_J|T|\alpha JM_J^{}`$ $`=`$ $`{\displaystyle \underset{m_l,m_l^{},M_S}{}}(1)^{2S2l+M_J+M_J^{}}(2J+1)`$ (54)
$`\times `$ $`\left(\begin{array}{ccc}S& l& J\\ M_S& m_l& M_J\end{array}\right)`$
$`\times `$ $`\left(\begin{array}{ccc}S& l& J\\ M_S& m_l^{}& M_J^{}\end{array}\right)\alpha lm_l|T|\alpha lm_l^{}`$ (57)
The matrix elements $`\alpha JM_J|T|\alpha JM_J^{}`$ are obtained in a basis where the quantization axis is perpendicular to the collision plane. We have now to perform an average of the transitions probabilities over all possible orientations since collisions are isotropic. This angular average is denoted by $`_{av}`$.

Using formula (16) by Sahal-Bréchot (Sahal1 (1973)), it can be shown that, in the irreducible tensorial operator basis, the angular average of the depolarization transition probability is given by:

$`P^k(\alpha J,b,v)_{av}={\displaystyle \frac{1}{2J+1}}{\displaystyle \underset{\mu ,\mu ^{}}{}}|\alpha J\mu |T|\alpha J\mu ^{}|^2`$
$`{\displaystyle \underset{\mu ,\mu ^{},\nu ,\nu ^{}}{}}\alpha J\mu |T|\alpha J\mu ^{}\alpha J\nu |T|\alpha J\nu ^{}^{}`$ (58)
$`\times {\displaystyle \underset{\chi }{}}(1)^{2J+k+\mu \mu ^{}}\left(\begin{array}{ccc}J& J& k\\ \nu ^{}& \mu ^{}& \chi \end{array}\right)\left(\begin{array}{ccc}J& J& k\\ \nu & \mu & \chi \end{array}\right)`$ (63)
Owing to the selection rules on the $`3j`$-coefficients, the summation over $`\chi `$ is reduced to a single term, since $`\chi =(\mu ^{}\nu ^{})=(\mu \nu )`$.

The depolarization cross section may then be computed from:

$`\sigma ^k(\alpha J,v)=2\pi {\displaystyle _0^{\mathrm{}}}P^k(\alpha J,b,v)_{av}b𝑑b`$ (64)
At small impact-parameters, the $`P^k(\alpha J,b,v)_{av}`$ values become highly oscillatory, and it is impractical to integrate numerically. In addition, at these separations the perturbation theory breaks down and the calculated RSU potentials become invalid. Thus, an impact-parameter cutoff is used (Anstee & O’Mara Anstee1 (1991)):

$`\sigma ^k(\alpha J,v)\pi b_0^2+2\pi {\displaystyle _{b_0}^{\mathrm{}}}P^k(\alpha J,b,v)_{av}bdb`$ (65)
Integration over the velocity distribution for a temperature $`T`$ can be performed to obtain the depolarization rate which is given by:

$`D^k(\alpha J,T)`$ $``$ $`n_H{\displaystyle _0^{\mathrm{}}}vf(v)dv(\pi b_0^2`$
$`+2\pi {\displaystyle _{b_0}^{\mathrm{}}}P^k(\alpha J,b,v)_{av}bdb)`$
where $`b_0`$ is the cutoff impact-parameter and here $`b_0=3a_0`$ as in Anstee & O’Mara (Anstee1 (1991)) (see Sect. 8).

## 7 Results

### 7.1 Our results for $`{}_{}{}^{1}P_{1}^{}`$ and $`{}_{}{}^{2}P_{3/2}^{}`$ levels

The perturbed atom is modelled as a single valence electron outside a spherical ionic core. If the energy of the state $`|\alpha l`$ of the valence electron is denoted by $`E_{\alpha l}`$ and the binding energy of the ground state by $`E_{\mathrm{}}`$, the binding energy of the level $`|\alpha l`$ is ($`E_{\mathrm{}}E_{\alpha l}`$) and is related to the effective principal quantum number by $`n^{}=[2(E_{\mathrm{}}E_{nl})]^{1/2}`$.

For example: consider the case of the neutral calcium line $`4227\AA `$ (transition $`4s^24s4p`$). The energy of the $`4p`$ level and the binding energy are $`0.1077\text{a.u.}(23652\text{cm}^1)`$ and $`0.2245\text{a.u.}(49305.72\text{cm}^1)`$ respectively (Wiese et al Wiese (1969)). Hence the effective principal quantum number is $`n^{}2.077`$.

Note that if the parent configuration of the level is not the same as that of the ground state, this should be accounted for by using the appropriate series limit for the excited parent configuration instead of the ground state binding energy. A particular example was explained in detail by Barklem, Anstee & O’Mara (Barklem3p (1998)) (see also Barklem, Piskunov & O’Mara Barklem4 (2000)). As in Anstee & O’Mara (Anstee1 (1991)), we can calculate the cross sections as a function of $`n^{}`$. Table 2 gives $`\sigma ^k(\alpha 1,v)`$ and $`\sigma ^k(\alpha \frac{3}{2},v)`$ as a function of $`n^{}`$ for a typical relative velocity $`v=10\text{km}\text{s}^1`$.

We thus can obtain the depolarization cross section for any level of any atom by interpolation of the values given in Table 2 for the appropriate $`n^{}`$ value for the state of interest. We notice that the dependence of cross sections on velocity is very close to

$`\sigma ^k(\alpha J,v)=\sigma ^k(\alpha J,v_0)({\displaystyle \frac{v}{v_0}})^{\lambda ^k(\alpha J)},`$ (67)
where $`v_0`$ is the velocity at which the cross section is known ($`10\text{km}\text{s}^1`$). Table 3 shows the so-called velocity exponent $`\lambda ^k(\alpha J)`$ (Anstee & O’Mara Anstee1 (1991)). In order to obtain the depolarization cross section for any velocity $`v`$ from Tables 2 and 3 equation (67) can be used.

Figures 3 and 4 display the depolarization cross sections (computed directly) and the corresponding depolarization rates for the CaI $`4p^1P_1(n^{}2.077)`$, MgI $`3p^1P_1(n^{}2.030)`$ and the NaI $`3p^2P_{3/2}(n^{}2.117)`$ states.

To obtain the depolarization rates we must integrate over the Maxwell distribution of velocities.

Using equation (6), it can be shown that for a mean velocity $`\overline{v}=\sqrt{\frac{8kT}{\pi \mu }}`$, the depolarization rates can be expressed as (Anstee & O’Mara Anstee1 (1991))

$`D^k(\alpha J,T)`$ $`=`$ $`({\displaystyle \frac{4}{\pi }})^{(\frac{1}{2}\lambda ^k(\alpha J))}\mathrm{\Gamma }(2{\displaystyle \frac{1}{2}}\lambda ^k(\alpha J))v_0\sigma ^k(\alpha J,v_0)`$ (68)
$`\times ({\displaystyle \frac{\overline{v}}{v_0}})^{1\lambda ^k(\alpha J)}`$
Since $`f(v)`$ depends on the reduced mass $`\mu `$, the depolarization rate is specific to a given atom. In order to obtain depolarization rates for any level of any atom by simple interpolation we applied the approximation: $`\mu =m_H,m_H`$ being the hydrogen mass. We verified that this approximation introduces a less than $`4\%`$ error in the depolarization rates.

Table 4 gives the alignment depolarization rates as a function of $`n^{}`$ for the temperatures of 5000 K and 6000 K.

### 7.2 Comparison with other methods

In Figure 5 we compare our alignment depolarization rates with quantum chemistry depolarization rates (Kerkeni Kerkeni2 (2002)) and the alignment depolarization rates obtained by replacing the RSU potential by the Van der Waals potential, $`V={\displaystyle \frac{C_6}{R^6}}`$, where $`C_6`$ is the Van der Waals constant.

We display only the $`k=2`$ case which is related to the linear depolarization (alignment). Our results show rather good agreement with quantum chemistry calculations. However, the Van der Waals potential underestimates the depolarization cross section. To explain this we now examine the sensitivity of the depolarization cross section to the RSU potential.

## 8 Dependence of depolarization cross sections on potentials

As a first check on the sensitivity of the results to the nature of the potential we introduce a localised perturbation to the interaction potential. We have chosen the case CaI 4$`p`$ for the check. The interaction potential is multiplied by a Gaussian magnification factor of the form (Anstee & O’Mara Anstee1 (1991))

$`G(R)=1+\mathrm{exp}((RR_0)^2).`$ (69)
The $`\sigma `$-symmetry interaction $`(1s,4p\sigma )`$ is more important than the $`\pi `$-symmetry interaction $`(1s,4p\pi )`$. Therefore only the $`(1s,4p\sigma )`$ interaction has been altered. $`R_0`$ is the position of the maximum perturbation.

Fig. 6 shows the effect of the $`R_0`$ variation on the depolarization cross section for three different velocities. We remark that the enhancement effect on the cross section shows a peak in the intermediate range ($`10a_0R_018a_0`$). We conclude that the depolarization cross section depends most strongly on the intermediate region of the potential. It is for this reason that the Van der Waals interaction potential, which is inaccurate in this intermediate region, underestimates the depolarization cross sections.

To study the sensitivity of the cross section to close collisions, a second check is made. We have calculated the alignment cross sections by varying the cutoff $`b_0`$ in equation (65). Fig. 7 shows the effect of the cutoff variation at $`v=10\text{km}\text{s}^1`$ and $`v=30\text{km}\text{s}^1`$. For $`b_0<8a_0`$ the cutoff variation does not have any effect on the cross section calculations.

Therefore it can be concluded that short-range collisions which involve ionic and exchange interactions do not influence the depolarization cross sections and hence the depolarization rates. The principal differences between the RSU potentials and those from quantum chemistry occur at small interatomic separations. Since the important impact-parameters are at medium interatomic distances, this explains why we obtain a rather good agreement between the ABO and quantum chemistry results. However, a more detailed study would be desirable but is outside the scope of the present paper because the exact quantum chemistry potentials are needed for quantitative comparison.

## 9 Conclusion

In conclusion, we have shown that the ABO theory is a powerful tool for quickly calculating a great number of depolarization cross sections for collisions with atomic hydrogen in the physical conditions of the solar atmosphere. The present paper is limited to $`p`$ states for testing the method. Comparison with quantum chemistry results gives us the conviction that it would be useful to extend our calculations to higher $`l`$-values which are also of importance for solar line polarization. Although less accurate than quantum chemistry calculations, an extension to heavy atoms and higher $`l`$-values is easier and faster within the present semi-classical method. We will also extend our semi-classical theory to depolarization cross sections by collisions with singly ionized atoms. This will be the subject of further papers.

###### Acknowledgements.

This paper is dedicated to Jim O’Mara’s memory, who followed with interest the progress of the present work and suddenly passed away on April 27, 2002. The authors thank the referee for valuable suggestions which improved the presentation of this work. We are indebted to B. Kerkeni who has communicated a number of her results before publication.


format: text

text:
# Accretion in Young Stellar/Substellar Objects

## 1 Introduction

The collapse ($`10^5`$ yr; Adams, Lada, & Shu, 1987) of a slowly rotating molecular cloud core leads to formation of a protostar and surrounding circumstellar accretion disk which can last up to $`10^7`$ yr in at least some cases (Muzerolle et al., 2000a; Alencar & Batalha, 2002). Observable signatures of disks and disk accretion include far- and mid-infrared emission from dust at a range of temperatures, near-infrared emission from hot dust and gas in the inner disk, and a variety of emission lines and UV/optical continuum excess due to accretion of material from the inner disk directly onto the star (see, e.g. Mannings, Boss, & Russell 2000 for reviews).

Historically, accretion from the disk onto the star was thought to occur through a hot equatorial “boundary layer” (e.g. Lynden-Bell & Pringle 1974; Bertout, Basri, & Bouvier 1988; Basri & Bertout 1989). Over the past decade, however, a new model of magnetospheric accretion has been widely adopted for classical T Tauri stars (hereafter, CTTSs). This model was originally discussed in the context of neutron stars by Ghosh & Lamb (1979) and first applied to CTTSs by Uchida & Shibata (1985), Bertout et al. (1988), and Konigl (1991). In this scenario, the inner regions of the circumstellar disk are disrupted by the stellar magnetosphere, which channels viscously accreting material out of the disk plane and towards the star along magnetic field lines; an accretion shock is formed as the free-falling material impacts the stellar surface. This model can also generate magnetocentrifugal outflows, which are launched from a region just exterior to the point of disk truncation (Shu et al, 1994), or from near the poles of the star (Hirose et al. 1997). Several of the detailed observational characteristics of CTTSs, such as blue continuum excesses including the Paschen and Balmer continuum slopes (e.g. Calvet & Gullbring, 1998), relatively small near-infrared excesses (e.g. Meyer, Calvet & Hillenbrand, 1997), permitted emission-line profiles (e.g. Muzerolle, Calvet, & Hartmann, 2001a; hereafter, MCH), magnetic field strengths (Johns-Krull et al. 1999; Johns-Krull & Valenti 2000), and slow rotation periods (e.g. Edwards et al. 1993), are successfully explained by such a model.

The presence of inverse P Cygni line profiles in the spectra of many CTTSs motivated Calvet & Hartmann (1992) and then Hartmann, Hewett, & Calvet (1994) to first investigate magnetospheric accretion flows as the source of the line emission. Radiative transfer calculations of magnetospheric accretion predict a characteristic blue asymmetry – with or without a redshifted absorption feature – in both hydrogen (Balmer, Paschen, and Brackett lines) and metallic (notably NaD, OI, and the CaII triplet) line profiles. These calculations successfully reproduce both the line fluxes and the line profiles (Hartmann, Hewett, & Calvet 1994; Muzerolle et al. 1998ab; MCH) observed in most CTTSs in the Taurus molecular cloud.

Results thus far, however, have been limited to T Tauri stars of restricted range in mass and age ($`M0.5M_{}`$; $`\tau <`$1-3 Myr ). Few high-resolution observations even exist of stars with masses $`<0.2M_{}`$. In order to develop a broader understanding of accretion disk activity, in this contribution we extend the analysis of accretion to a mass range more representative of the stellar mass function, focusing on young brown dwarfs and stars in the mass range 0.05-0.2 M (hereafter referred to as very low mass, or “VLM”, objects). Our results for a single object among this sample, V410 Anon 13, appeared in Muzerolle et al. (2000b), the first measurement of accretion in a candidate sub-stellar mass young T Tauri object. The derived $`\dot{M}`$ of $`5\times 10^{12}\mathrm{M}_{}\mathrm{yr}^1`$ remains the lowest measured accretion rate, even when the larger sample presented here is considered. We have also recently extended the analysis in age to study $``$10 Myr old stars via an investigation of disk accretion in the TW Hya association (Muzerolle et al 2000a; Muzerolle et al. 2001b).

We describe our sample in section 2 and our observations, moderate- to high-resolution optical spectra, in section 3. Section 4 presents results from the spectra, including H$`\alpha `$line profiles, measurements of continuum veiling, and model-derived mass accretion rates, as well as radial and rotation velocities for part of the sample. Finally, we discuss the implications of our results for accretion in young VLM objects in section 5.

## 2 Sample

The sample was selected from previously identified young low mass objects in several different star forming regions and young clusters, including Taurus, IC 348, $`\sigma `$ Ori, Upper Sco, and $`\rho `$ Oph (see references in Table 1). Our observations were designed in early 2000 to include every known object with spectral type M5 or later and age $``$10 Myr or younger. Since our observations were obtained, many more such objects have been discovered in star forming regions and young clusters (e.g. Briceño et al. 2002) and thus our sample is not complete, though it is likely to be representative of these very low masses. The final observed sample encompasses 45 young stars and brown dwarfs with spectral types M4-M8.5. At the 0.1-8 Myr ages of these regions, the masses inferred from pre-main sequence evolutionary tracks (e.g. D’Antona & Mazzitelli, 1998; Chabrier et al. 2000) span the substellar boundary, ranging from 0.02-0.14 M.

Table 1 contains basic information on our sample, including spectral types found in the literature from which, in combination with literature photometry, the other parameters (logT, logL, extinction, H-K excess) have been calculated. Masses and ages were determined from the D’Antona & Mazzitelli (1997, including unpublished updates provided in 1998) pre-main sequence evolution tracks, adopting the transformations used in Hillenbrand & Carpenter (2000). Our observations do not include enough dwarf and giant spectral templates to re-examine with any great accuracy the previously published spectral types. We note that one particular object, LkH$`\alpha `$358, has types from Cohen & Kuhi (1979; M5.5) and Kenyon et al. (1998; M5-M6) that are clearly inconsistent with our spectrum, which appears to be more like K8. This star also displays spatially extended H$`\alpha `$ and \[SII\] doublet emission, likely originating in a jet. Our spectral type estimate of another object, Haro 6-28, is also discrepant with the previously reported spectral type of M5 (Cohen & Kuhi 1979); we estimate M2.5. Both of these objects exhibit accretion signatures, including continuum veiling. Figure 1 shows the HR diagram for our sample, and illustrates the estimated mass and age ranges spanned in our study of accretion diagnostics.

## 3 Observations

We obtained spectra at Keck II with the ESI echellette spectrograph (Sheinis et al. 2002) on Jan. 25-28, 2001, and at Keck I with the HIRES echelle spectrograph (Vogt et al. 1994) on Jan. 10-11 2000 and June 30, 2001. Extraction of one dimensional spectra from the two-dimensional images was performed using the MAKEE software written by T. Barlow.

ESI provides wavelength coverage from 3900 - 11000 Å, but with only moderate resolution, $`R8000`$ (37.5 km s-1), while HIRES has a smaller wavelength range, for our setup, 6320 - 8730 Å with gaps between the orders and typical echelle resolution of $`R34,000`$ (8.8 km s-1). All sets of data cover the important H$`\alpha `$emission line, our primary diagnostic of accretion. Despite lower resolution, the ESI data still resolve this line when accretion is present, since in that case the line widths are typically several hundred km s-1; however, chromospheric emission line profiles and photospheric atomic absorption features generally are not resolved. The trade-off for lower resolution is the benefit of increased wavelength coverage into the blue which yields a larger baseline for measuring continuum veiling (although in most of our objects, the S/N is too low to make measurements below $`5000`$ Å). Line profiles are more easily resolved with HIRES, but gaps between the orders prevent simultaneous measurement of all emission lines of interest for accretion studies.

## 4 Results

### 4.1 Radial and Rotational Velocities

We determined $`v\mathrm{sin}i`$ values for the HIRES sample only - the resolution of ESI is too low for any meaningful limits given typical values for low mass young stars (e.g., Hartmann et al. 1986; White & Basri 2003, hereafter WB03). The rotation velocities were determined by cross-correlating photospheric absorption features in small portions of each spectrum (see Tonry & Davis 1979; Hartmann et al. 1986). We selected intervals of $`20`$ Å in several different orders containing relatively strong absorption features with good S/N. Values of $`v\mathrm{sin}i`$ were derived from the widths of the cross-correlation peaks, which were calibrated using template spectra of Gl 447 ($`v\mathrm{sin}i<2`$ km s-1, unresolved at HIRES resolution; Delfosse et al. 1998) convolved with a range of rotation velocities. Results for the different orders were then averaged together, with the standard deviation taken as the final uncertainty. We present these values in Table 2. The rotation velocities range from 8 - 35 km s-1, with most being $`<`$ 20 km s-1. The two fastest rotators, FP Tau and LkCa 1, both at 34 km s-1, show evidence of line doubling in the Ca II triplet emission lines; thus, the absorption line widths (which do not show obvious signs of doubling) may be due to binarity rather than rapid rotation.

We determined heliocentric radial velocities through similar cross-correlation techniques for both the HIRES and ESI samples (see Table 2). For the HIRES spectra, we used the dwarf standards Gl 402 and Gl 447 ($`V_r=2`$ and -31 km s-1, respectively; Delfosse et al. 1998) as templates. For the ESI sample, we used Gl 83.1 and Gl 406 ($`V_r=28.567`$ and 19.482 km s-1, respectively; Nidever et al. 2002). Due to the relatively low resolution of the ESI data (R = 8,000), the cross-correlation power is dominated by broad molecular features which forced us to be selective in the wavelength regime and range employed. As for the $`v\mathrm{sin}i`$ measurements, results for different wavelength ranges were averaged together, with the errors taken from the standard deviation.

### 4.2 Membership Information

Many of the stars we have studied with HIRES and ESI are candidate (as opposed to confirmed) members of the Taurus, IC 348, Sigma Ori, Sco-Cen, and Ophiuchus stellar associations. In these cases, membership is supported by consistency in a color-magnitude diagram with the expectations for a low-mass member at the distance and age of the cluster, combined with a late spectral type derived from low dispersion spectroscopy. Our higher dispersion data contain a variety of spectroscopic indicators that may further inform membership status. At the 1-10 Myr ages of the associations contained in our sample, high lithium abundance and low surface gravity are expected (we list Li equivalent widths in Table 2). Further, we have checked consistency of measured radial velocities with those previously published for known cluster members.

In Taurus, we find that CIDA-13 lacks substantial Li I, although its radial velocity is consistent with the average for Taurus members; it is likely a foreground dwarf star. The Sigma Ori cluster has two candidate members (SOri-45 and SOri-46) with K I lines that are more similar to those of M dwarfs than of low-gravity pre-main sequence stars, as well as radial velocities which are substantially different from the other Sigma Ori stars. In the Upper Sco association, we find three stars (UScoCTIO-85, UScoCTIO-121, UScoCTIO-132) which both lack substantial Li I and have dwarf-like K I absorption.

As suggested in the tablenotes to Table 1, we question the membership of these 6 stars in their respective associations. Many of these stars are among the weakest H$`\alpha `$ emitters (Table 3), and the emission line characteristics in general are more consistent with active main sequence dwarfs rather than young TTSs/brown dwarfs. We do not include these objects in any of the subsequent statistics or analysis.

### 4.3 Emission Lines

Emission line measurements for our entire sample are listed in Table 3, including equivalent widths for a variety of permitted and forbidden lines, as well as the H$`\alpha `$line widths at 10% of peak intensity (a good discriminator of accretion versus chromospheric emission; Muzerolle et al. 2000a, WB03). Measurements for two main sequence late-M dwarfs (Gl 406 and LHS 2351) are also given as a comparison. The emission line spectra of the VLM young objects are very similar to those of higher-mass young stars, including both CTTSs and “weak line” TTSs (WTTSs). Balmer and He I emission are virtually ubiquitous, and are typically stronger than in the main sequence dwarfs. The equivalent widths are larger than seen in higher-mass WTTSs, but this is probably a result of the cooler photospheres, rather than intrinsically stronger emission. Ca II triplet emission is detected in nearly half of the sample, while forbidden emission such as \[OI\] $`\lambda `$6300 is detected in only $``$25% of the sample (though this smaller fraction may be due in part to the drop-off in S/N of many of the spectra around this wavelength).

Most of the emission line profiles are narrow ($`200`$ km s-1) and symmetric, typical of WTTSs and dMe stars. However, 13 objects show broader profiles at H$`\alpha `$(see Fig. 2), including many with blueshifted asymmetries, and several with blueshifted absorption components (e.g., FP Tau, MHO-6, IC348-415). Such features are more reminiscent of CTTSs. The shapes of the upper Balmer lines generally mirror H$`\alpha `$, while somewhat narrower, and blueshifted absorption components decrease in strength as one moves up the series. Broad emission components in the Ca II and He I lines are seen only in the earliest-type object in the sample, LkH$`\alpha `$358. This star is also the only one that exhibits a high-velocity component at \[OI\]. Overall, the most active of the $`0.050.2M_{}`$ VLM young objects have emission characteristics similar to the least active among the well-studied $`0.5M_{}`$ CTTSs.

### 4.4 Accretion Signatures

We interpret our data in the context of accretion disks and accretion-driven outflows associated with these young VLM objects, following our previous work on higher mass TTS. Hot UV/optical continuum veiling of photospheric absorption features, as well as narrow emission components in Ca II and He I, provide clues to the density, temperature structure and filling factor of the accretion shock; these parameters may differ from those characteristic of the well-studied $`0.5M_{}`$ stars, as $`M_{}/R_{}`$ varies. Observations of broad Balmer lines provide detailed information on the nature of accreting magnetospheres. Forbidden emission lines such as \[OI\] and \[SII\], and blueshifted absorption components in permitted lines, are typically associated with accretion disk-driven outflows (e.g. Hartigan, Edwards, & Ghandour 1995; Calvet 1997). We use all of these features to identify objects which are likely undergoing active accretion, but focus primarily on veiling and H$`\alpha `$emission line profiles in our analysis.

#### 4.4.1 Veiling

One familiar manifestation of accretion in T Tauri stars is the veiling of photospheric absorption features by excess continuum emission produced in the accretion shock. Veiling is typically defined as the ratio of the excess flux to the photospheric flux at a given wavelength: $`r_\lambda =F_{\lambda ,ex}/F_{\lambda ,cont}`$. Consistent with the historical suggestions that accretion signatures are weaker in lower mass stars, we find very little veiling in our VLM sample. We have measured or placed upper limits on the veiling for each of the 13 likely accretors (as defined from emission line indicators; see below), using comparison spectra of similar spectral type (e.g., Hartigan et al. 1989). In order to avoid mismatches due to gravity-sensitive features (see WB03 for a more detailed discussion), we opted to use the non-accreting subsample, defined as those which clearly lack infalling gas signatures at H$`\alpha `$, as veiling templates. Since veiling at red optical wavelengths is much less sensitive a diagnostic of accretion than H$`\alpha `$, these templates should not have any intrinsic veiling of their own. We note that WB03 measured an upper limit $`r0.22`$ at 6500 Å for one of our templates, V410 X-ray 3, using a combination of dwarf and giant standard spectra.

Veiling measurements at several wavelength regions are given in Tables 4 and 5. Due to the gaps in the wavelength coverage of the HIRES spectra, we could not measure veiling at exactly the same wavelength regions as ESI, hence we list these separately. The wavelength regions chosen for the ESI sample contain broad molecular absorption features; we focused on these rather than narrow atomic features in order to avoid line blending problems due to the lower resolution of the ESI spectrograph. For the two objects earlier than $``$M3, Haro 6-28 and LkH$`\alpha `$358, the derived veiling is more uncertain since the atomic features dominate. These problems are much less significant for the HIRES sample; however, the grid of templates is less extensive, and slight mismatches due to differing rotation rates result in greater uncertainties. Note especially the results for the V410 Anon 13 and CIDA 14 HIRES spectra - the closest template in spectral type is about 0.5 subclasses earlier, and rotates faster, which could result in slight underestimates of the true veiling (although the results for V410 Anon 13 are consistent with those of the ESI spectrum).

Most of the objects at M4 or later exhibit no significant veiling (i.e., $`2\sigma `$ of the typical uncertainty of $`\pm 0.1`$) at any of the wavelengths. One low-mass object, IC348-415 does have significant veiling with $`r0.5`$ at 6200 Å, though this value is highly uncertain due to low S/N; this object also has the strongest H$`\alpha `$emission of the entire sample, and our best profile model indicates the highest $`\dot{M}`$ of the $`>`$M4 sample. Nevertheless, we consider this detection to be marginal at best. Only in the 2 earliest-type objects in the sample, Haro 6-28 and LkH$`\alpha `$358, do we measure robust levels of veiling, indicating higher levels of accretion than in the lowest-mass sample as a whole. The veiling levels are consistent with the general emission line characteristics of the stars, such as the presence or absence of broad Ca II triplet emission components, which were noted by Muzerolle, Hartmann, & Calvet (1998) to have a threshold behavior whereby the broad component does not appear when the mass accretion rate is $`10^8\mathrm{M}_{}\mathrm{yr}^1`$.

#### 4.4.2 H$`\alpha `$Profiles

Continuum-normalized H$`\alpha `$line profiles are shown in Figure 2. Traditionally, the threshold between CTTSs and WTTSs has been set at a canonical boundary of H$`\alpha `$$`EW=10\mathrm{\AA }`$, which divides accreting from non-accreting TTSs around the peak of the mass function in Taurus, $`0.5M_{}`$. However, this threshold is inappropriate for much lower mass objects, where the decreased optical photospheric continuum results in typically larger equivalent widths despite comparable emission fluxes from chromospheric activity. Recently, WB03 compared H$`\alpha `$10% line widths to equivalent widths for TTSs of a wide range of masses, and calibrated a spectral type-dependent equivalent width boundary between chromospheric and accretion emission at H$`\alpha `$, down to late-M objects. Figure 3 shows the distribution of H$`\alpha `$10% line widths, as well as a comparison with equivalent widths, for our bona-fide young object sample. Here, we also focus on the fact that the line profile shape is a much more robust discriminator of chromospheric vs. accretion emission. The broad ($`>`$200 km s-1), asymmetric profiles produced in infalling accretion flows are easily identified even in the lower-resolution ESI spectra. H$`\alpha `$emission from WTTSs, on the other hand, is similar to that seen in chromospherically active main-sequence dwarfs (e.g., Worden, Schneeberger, & Giampapa 1981), with relatively small equivalent widths and narrow ($`<`$200 km s-1), symmetric line profiles. Our accretion models (see below) indicate that accretion profiles in VLM objects should be detectable down to very low accretion rates ($`10^{12}\mathrm{M}_{}\mathrm{yr}^1`$), much lower than the sensitivity limits of any other optical diagnostics, including continuum veiling111There is a possibility that at least some of the broad H$`\alpha `$profiles could be due to flares. We believe this is unlikely for several reasons. Previously published profiles of M-dwarf flares sometimes show a very weak broad gaussian component ($`FWHM200`$ km s-1), often redshifted, superposed on the (stronger) chromospheric profile, and the equivalent widths are typically $`20\mathrm{\AA }`$ (e.g., Abdul-Aziz et al. 1995; Abranin et al. 1998; Martín 1999). V927 Tau, which shows a similar profile, could be such a case. In general, however, these flare profiles are not consistent with our observations. Also, our repeat observations of V410 Anon 13, separated by 9 months, show a nearly identical H$`\alpha `$profile, indicating a relatively stable source for the emission, thus favoring accretion over flare activity..

We have identified likely accretors on the basis of line profile shape, namely whether the H$`\alpha `$lines are broad ($`200`$ km s-1) and asymmetric. A total of 10 objects meet both of these criteria, and we consider them to be strong candidates for accretion. Another 4 objects (V927 Tau, IC348-205, IC348-382, and UScoCTIO-75) meet the line width criterion, but are relatively symmetric. Yet another 3 objects (MHO-4, MHO-5, and UScoCTIO-128) meet neither of the H$`\alpha `$accretion criteria, but show other potential signposts of accretion, such as forbidden line emission or Ca II triplet emission fluxes stronger than typical dwarf/WTTS levels. Of these 7 objects with ambiguous accretion indicators, we believe at least 3 are compelling cases for the presence of active accretion. IC348-205 and IC348-382 have considerably stronger Balmer and Ca II triplet emission than the VLM WTTSs, and, in both cases, H$`\alpha `$is quite broad, if relatively symmetric, with two of the largest EWs in the entire sample. MHO-5 exhibits relatively narrow and symmetric H$`\alpha `$emission; however, the equivalent width is larger than typical of VLM WTTSs of similar spectral type (e.g., WB03). The profile shape of H$`\alpha `$in this object is probably due to a pole-on geometry for the accretion flow, which produces narrower, more symmetric emission profiles since most of the high-velocity infalling gas is moving perpendicular to the line of sight (§5.2; MCH). MHO-5 also exhibits surprisingly strong \[OI\] emission, which is never seen in higher-mass WTTSs. Briceño et al. 1998 detected similarly strong H$`\alpha `$and \[OI\] emission in MHO-5, so we discount variability due to flaring. For our subsequent analysis, we consider these 3 objects, plus the original 10 H$`\alpha `$-selected strong candidates, as our total “accretor” sample (see Tables 4 and 5). The distribution of H$`\alpha `$10% line widths and equivalent widths for this subsample are delineated in Figure 3. The remaining 22 objects exhibit emission characteristics consistent with chromospheric activity alone – narrow, symmetric H$`\alpha `$, no forbidden emission, and little or no Ca II emission.

Yet another ambiguous case of interest is V927 Tau, whose H$`\alpha `$profile includes a narrow, symmetric component with a small self-absorption reversal at line center, all hallmarks of chromospheric emission. However, there is also a very weak, slightly asymmetric broad component that could be due to accretion. The equivalent width of the broad component is roughly 2 Å, more than an order of magnitude lower than is typical of the accretor sample. These characteristics are reminiscent of flare activity, and we cannot reliably distinguish between that or very low-level accretion without repeat observations, so we do not include this object in the accretor sample.

## 5 Modeling Accretion Flows

### 5.1 Continuum Excess Models

Following the procedures of Calvet & Gullbring (1998, CG98), we calculate accretion shock emission and estimate the amount of veiling expected for a given mass accretion rate, In short, infalling magnetospheric material accretes onto the stellar photosphere through a shock just above the photosphere. The material is assumed to be in plane-parallel columns. In each column, soft X-ray emission from the shock heats both the pre-shock region above the shock and the stellar photosphere below it. The shock flux $`F_{shock}`$ is essentially the sum of the spectra from the heated photosphere and the preshock region. The spectral shape of the shock flux depends on the energy flux carried by each accreting column, $``$, while the level of the emission depends on the surface coverage of accreting columns, $`f4\pi R_{}^2`$ (CG98). The total luminosity of the shock is

$$L_{shock}=f4\pi R_{}^2,$$
(1)
neglecting the intrinsic stellar flux, which is at least an order of magnitude lower than the flux carried by the column in these cool stars. At the same time,

$$L_{shock}=\zeta L_{acc}=\zeta \frac{GM_{}\dot{M}}{R_{}},$$
(2)
where $`\zeta =11/R_{mag}`$, with $`R_{mag}`$ being the radius at which magnetospheric infall begins.

We have calculated models for two generic stars with parameters corresponding to our sample: an M5 star ($`T_{eff}=2900`$ K, $`M_{}=0.12\mathrm{M}_{}`$, $`R_{}=1.2\mathrm{R}_{}`$) and an M6 star ($`T_{eff}=2800`$ K, $`M_{}=0.05\mathrm{M}_{}`$, $`R_{}=0.5\mathrm{R}_{}`$). Models were calculated for log $``$ = 10, 11, and 12, which correspond to the range of energy fluxes that provide better fits to the optical and ultraviolet excess fluxes in higher mass T Tauri stars (CG98; Gullbring et al. 2000). The emission from the shock was calculated for $`\dot{M}=10^{10}\mathrm{M}_{}\mathrm{yr}^1,10^9\mathrm{M}_{}\mathrm{yr}^1`$ and $`10^8\mathrm{M}_{}\mathrm{yr}^1`$, using eqs. (1) and (2) to obtain the filling factors $`f`$ for given values of $`\dot{M}`$ and $``$. We take $`R_{mag}=3`$, in accord with the H$`\alpha `$emission line models (§6).

To calculate the veiling, we compare the flux arising from the accretion shock with the intrinsic photospheric flux. Here, we calculated the photospheric flux using bolometric corrections from Luhman (1999) to convert from model luminosities to absolute J magnitudes, and then used standard colors to obtain corresponding V, Rc, and Ic magnitudes. Colors were taken from Leggett (1992) for types later than M4, and from Kenyon & Hartmann (1995) for earlier types. Fluxes at the effective wavelength of each filter were obtained using the zero points of Bessell & Brett (1988). Fluxes were then scaled to the assumed model radii, and then the veiling was calculated at each filter wavelength by

$$r=\frac{F_{excess}}{F_{photosphere}}=\frac{fF_{shock}}{(1f)F_{phot}}.$$
(3)
In Figure 4, we show the veiling expected for the two model stars for log $`\dot{M}`$ = -10, -9, and -8, and log $``$ = 10, 11, and 12 (note the logarithmic scale). For a given $`\dot{M}`$, the veiling at 5500 Å does not change much with $``$, though the model with lower $``$ has a larger veiling at longer wavelengths. This arises from the fact that at optical wavelengths, the shock flux is dominated by quasi-blackbody emission from the heated photosphere, which has a temperature $`T_{hp}(\mathrm{\hspace{0.17em}3}/4\sigma )^{1/4}`$ (CG98). Thus, models with lower $``$ emit more at redder wavelengths.

It is apparent from Figure 4 that the negligible veiling measurements for most of the accretors (Tables 4 and 5) indicate $`\dot{M}10^{10}\mathrm{M}_{}\mathrm{yr}^1`$. In Figure 5, we compare the detected veiling for 4 accretors with predictions from shock models. Energy fluxes that provide the best fit to the observations, as well as the estimated mass accretion rate are also shown. The filling factors of the shock, given in the caption of Figure 5, are rather small, $`<<1`$%, similar to the case of the more massive CTTSs (CG98). The accretion rates are in excellent agreement with those determined from H$`\alpha `$models, as described in the next section.

### 5.2 H$`\alpha `$Emission Models

Previously, we have developed radiative transfer models of Balmer emission line profiles for $`0.5\mathrm{M}_{}`$ CTTS (MCH). These models treat line emission produced by magnetospheric accretion flows, where gas accreting through a circumstellar disk is channeled onto the star by the stellar magnetic field. The model assumes an ideal dipolar field geometry, and a specified gas temperature distribution and mass accretion rate; the gas density distribution is determined directly from the geometry and the gas velocity (a function of the mass and radius of the star). Further details can be found in MCH. Typical results show a characteristic profile shape with blue asymmetry, due to occultation effects and some absorption of the red wing where the receding part of the flow is seen in projection against the star and accretion shock. Frequently, outright redshifted absorption components are seen in the upper Balmer lines, but rarely in H$`\alpha `$because of thermalization and opacity broadening effects. The qualitative similarity of most of the asymmetric H$`\alpha `$profiles of the sample presented herein to these general model characteristics lends strong support for our application of these models to the lower-mass objects.

Here, we employ the same model, using two values of the stellar mass (0.05 and 0.15 $`\mathrm{M}_{}`$) and radius (0.5 and 1 $`\mathrm{R}_{}`$) that span the range appropriate for our VLM sample (see Table 1). The mass range will change somewhat depending on which evolutionary tracks are used; however, this should not affect our modeling applications since the model gas velocities, and, hence, line widths, are only sensitive to $`\sqrt{M}`$, and the derived masses are unlikely to be off by more than a factor of two. The radius is dominated by uncertainties in the spectral type - $`T_{eff}`$ transformation; a 200 K error in $`T_{eff}`$ translates to $`30\%`$ uncertainty in the radius. Since the model gas density does not depend very strongly on the radius (see MCH, and references therein), our model constraints are insensitive to this level of uncertainty.

We use the same temperature constraints as those found to reproduce the 0.5 $`\mathrm{M}_{}`$ CTTS profiles, which should be a reasonable assumption since the gas heating is likely mechanical and independent of the stellar $`T_{eff}`$. These constraints result in an inverse proportionality between the gas temperature and density. Since the accretion rates of all the objects we are considering here are $`10^9\mathrm{M}_{}\mathrm{yr}^1`$ (as inferred from the lack of continuum veiling), the accretion flow gas temperatures must be $``$ 10,000 K to obtain sufficient H$`\alpha `$emission, thus the gas is almost completely ionized. As an independent check of the adopted temperature constraints, we show several H$`\alpha `$models as a function of both temperature and $`\dot{M}`$ in Figure 6. The models at the higher value of $`\dot{M}`$ ($`10^9\mathrm{M}_{}\mathrm{yr}^1`$) exhibit broader line widths than those at the lower $`\dot{M}`$, as opacity broadening effects kick in due to the higher density; the profile is still broader at a lower temperature of $`T_{max}=8000`$ K (temperatures lower than this are not applicable at this value of $`\dot{M}`$ as the line flux becomes much lower than observed). At the lower value of $`\dot{M}`$ ($`10^{10}\mathrm{M}_{}\mathrm{yr}^1`$), the density is too low to produce opacity broadening in the wings, and the profile shape is nearly identical for $`T_{max}=10,000`$ and 12,000 K. Thus, the model temperature can be constrained enough to obtain more accurate constraints on the mass accretion rate.

As for the flow geometry, we have selected the fiducial magnetosphere inner and outer radius values from MCH (2.2-3 R). A larger magnetosphere would produce higher line fluxes and more strongly centrally-peaked profiles for a given accretion rate, while a smaller magnetosphere would produce redward asymmetries in the line profiles due to stellar occultation effects. The fiducial size is also consistent with inner hole sizes inferred from observed ($`HK`$) and ($`KL`$) excesses (Table 1; Meyer et al. 1997; Liu et al. 2002). Varying the magnetospheric size by roughly 20% may still produce an adequate match to the observed profiles, and still be consistent with infrared excesses. This results in an uncertainty in our $`\dot{M}`$ constraints of about a factor of 2-3 (see MCH for a detailed comparison of magnetospheric size, accretion rate, and line flux).

Finally, we have not included the effects of rotation in the models presented here. MCH showed that for $`V_{rot}20`$ km s-1, consistent with most of our $`v\mathrm{sin}i`$ measurements (§4.4), rotation does not have a significant effect on the line profile. A few of the profiles we model below show evidence for faster rotation, but since we cannot measure $`v\mathrm{sin}i`$, we have elected not to include rotation effects for these. $`\dot{M}`$ constraints are not affected by this omission, since the model gas density is virtually unchanged for the plausible range of rotation rates in the low mass objects.

For each accretor H$`\alpha `$profile, we constructed a model by adopting the above parameter constraints, varying only the mass accretion rate and inclination angle to find the best match to the observed profile shape and flux. A range of accretion rates were explored, $`10^{12}<\dot{M}<10^9\mathrm{M}_{}\mathrm{yr}^1`$, with the lower limit set where accretion emission is produced at an undetectable level, and the upper limit set by the veiling constraints. The inclination angle was varied from pole-on to edge-on orientations of the star/disk system. Our best model matches are shown in Figure 7, with parameters listed in Table 6. We have not modeled the 3 accretors earlier than M5, due to the temperature parameter degeneracy problem mentioned above (we did derive accretion rates from the veiling for two of these, Haro 6-28 and LkH$`\alpha `$358, with $`\dot{M}>10^9\mathrm{M}_{}\mathrm{yr}^1`$ in both cases; see Figure 5). We also do not show a model for V410 Anon 13, the subject of an earlier paper (Muzerolle et al. 2000b). Note that in two objects shown in Figure 7, MHO-6 and IC348-415, the H$`\alpha `$profiles show blueshifted absorption components. Such features are produced in accretion-powered winds, which are not treated in the model calculations. In these cases, we have tried to match only the red and high-velocity blue wings.

The overall agreement between the observed and model profiles is remarkably good, especially considering the idealized model geometry. In reality, the accretion flows are probably not axisymmetric, and in many cases may be broken up into discrete “streams” (e.g., Muzerolle et al. 2000b). The discrepancies near line center for some of the IC 348 objects (higher predicted flux than observed) may be a result of neglecting rotation in the model calculations; we do not have rotation velocities for those objects, which could conceivably be rotating much faster than typical of the Taurus objects in the HIRES subsample. Previous models we have generated for higher-mass stars with rapid rotation generally produce somewhat broader, less-strongly peaked profiles with doubling at line center. In addition, the assumption of a large velocity gradient implicit in our use of the Sobolev approximation in the model calculations for the line source function break down near line center, which mainly probes the outermost region of the magnetospheric flow where velocities are relatively small and densities are correspondingly large. Thus, we may be overestimating the amount of flux that can escape from the outer regions of the magnetosphere. Nevertheless, these models provide quantitative estimates of mass accretion rates (which we show in Tables 4 and 5) that are accurate to within a factor of $``$3-5, given the preceding caveats. Note that these H$`\alpha `$-derived accretion rates are consistent with the values and upper limits derived from the veiling, as shown in the previous section. Although our $`\dot{M}`$ estimates are not as accurate as UV excess measurements, they are likely the best achievable constraints for VLM objects at present, given the generally small accretion rates, which yield only a small amount of excess emission in the UV.

## 6 Discussion

The clear signatures of accretion exhibited in roughly 30% of our late-type sample indirectly reveal the presence of circumstellar accretion disks, and indicate that such disks are relatively common in young VLM objects. Observations of infrared excesses yield further evidence of disk accretion, and several recent studies along these lines have shown that the fraction of VLM young stars and brown dwarfs with disks is indeed – and not surprisingly – similar to that of the higher mass TTSs (e.g., Muench et al. 2001; Liu, Najita & Tokunaga 2002). Examination of Table 1 shows that a fairly large fraction of our sample exhibits near-infrared excess, indicative of circumstellar disks. The level of excess at these wavelengths suggests circumstellar disks with fairly small inner holes, relatively consistent with the inner magnetospheric radius of our H$`\alpha `$models. However, most of the excesses are fairly small ($`\mathrm{\Delta }(HK)0.2`$), and we cannot definitively rule out spots, spectral type mismatches, or cooler, unresolved companions as possible sources of apparent excess.

In any case, the general agreement between the observed and model H$`\alpha `$profiles we present here shows that magnetically mediated accretion from a disk is applicable to VLM young objects. Thus, the disk accretion paradigm, which is a common (if not universal) aspect of early stellar evolution, appears to work over a very large range of masses, from substellar objects with $`M<0.08\mathrm{M}_{}`$ to at least the highest-mass TTS at $`M2\mathrm{M}_{}`$ (and possibly to even higher masses, including the Herbig Ae/Be stars, some of which exhibit magnetospheric infall signatures; Muzerolle et al. in preparation). A key parameter in the physical description of accretion is the mass accretion rate, which guides both the growth of a young star/brown dwarf and the evolution of its disk. Mass accretion rates have been measured extensively at typical TTS masses of $`0.51.5\mathrm{M}_{}`$ (e.g. Valenti et al. 1993; Hartigan, Edwards, & Ghandour 1995; Gullbring et al. 1998; Hartmann et al. 1998; White & Ghez 2001); however, in most cases the mass ranges probed have not been sufficiently broad or complete to effectively study any mass dependences. Rebull et al. (2000; 2002) did estimate $`\dot{M}`$ for $`M0.251\mathrm{M}_{}`$ stars in Orion and NGC 2264 using $`U`$-band excesses, and found a loose correlation with $`\dot{M}M`$. However, their photometrically-selected samples are incomplete at average or lower values of $`\dot{M}`$, and rather large extinction uncertainties result in a more uncertain $`U`$-excess calibration.

More recently, WB03 attempted to measure accretion rates from continuum veiling in a sample of VLM objects in Taurus (including four objects in common with our sample). Most of their sample lacked measurable veiling (as does ours), yielding upper limits on accretion rates of the order $`10^{10}\mathrm{M}_{}\mathrm{yr}^1`$. They did find two objects 222Both of these, GM Tau and CIDA-1, were previously identified from low-resolution spectra as continuum objects, with indeterminate spectral types. with $`M<0.2\mathrm{M}_{}`$, including one probable brown dwarf, that exhibited considerable veiling, yielding accretion rates of approximately $`23\times 10^9\mathrm{M}_{}\mathrm{yr}^1`$. Together with our previously reported accretion rate for VLM object V410 Anonymous 13 (Muzerolle et al. 2000b), and values for higher-mass Taurus TTS, WB03 also find an approximate correlation of $`\dot{M}M`$. However, the number of $`\dot{M}`$ values at the lowest masses that are not upper limits is small.

With H$`\alpha `$modeling, we can probe to accretion rates several orders of magnitude lower than possible from veiling/UV excesses. Figure 8 shows our derived $`\dot{M}`$ values as a function of mass, along with the WB03 results, and values for higher-mass TTS from Gullbring et al. (1998) and White & Ghez (2001). Note that WB03 and White & Ghez (2001) used the Baraffe et al. (1998) pre-main sequence evolution tracks to determine masses; we have scaled those masses, as well as their values of $`\dot{M}`$ (since they used $`\dot{M}=L_{acc}R_{}/GM_{}`$), to the D’Antona & Mazzitelli (1998) evolution tracks in order to be consistent with our results. There is a clear correlation between mass accretion rate and stellar/substellar mass, with an approximate proportionality of $`\dot{M}M^2`$. Such a dependence is steeper than in previous investigations.

There are several caveats that must be emphasized when interpreting the above result. One is that the derived (sub)stellar masses are dependent on which evolutionary tracks are used; for example, the Baraffe et al. (1998) tracks yield systematically larger masses, by a factor of $``$ 1.5-2. This also affects accretion rates derived from veiling or the UV excess, since those methods require the stellar mass and radius in order to derive $`\dot{M}`$ from the accretion luminosity. However, we checked the mass-$`\dot{M}`$ comparison using Baraffe et al. (1998) masses, and find the same proportionality.

Another caveat is that objects accreting at levels low enough that H$`\alpha `$emission from the accretion flow (or UV emission from the accretion shock) is undetectable, would be missed in all surveys conducted to date. Such an effect is probably more relevant to the VLM objects since the minimum $`\dot{M}`$ at $`0.1\mathrm{M}_{}`$ needed to produce observable H$`\alpha `$emission, as determined from our models, is within about a factor of 2 of the lowest measured value in our VLM sample. At $`0.5\mathrm{M}_{}`$, the models indicate a higher minimum $`\dot{M}`$ (since the higher mass of the central object results in higher infall velocities, which in turn results in smaller magnetospheric gas densities for the same accretion rate), but one which is over one order of magnitude lower than the lowest measured value at that mass. Thus, at least at the higher-mass TTSs, the observed lower limit to $`\dot{M}`$ is probably real333This is also somewhat dependent on age, since older objects at the same mass have smaller radii, which results in somewhat higher densities at the same $`\dot{M}`$. Nevertheless, all of the objects plotted in Figure 8 are at roughly the same age, 1-3 Myr. See Muzerolle et al. (2000a) for results at 10 Myr..

A final caveat is that VLM objects accreting at rates higher than those measured thus far could be missed due to either veiling completely filling in photospheric features, or very high extinctions. However, the number of so-called “continuum stars” is fairly small (6 in Taurus, as listed in Kenyon & Hartmann 1995), and not necessarily all are VLM objects. One of three Taurus continuum stars observed by WB03, GN Tau, turned out to be M2.5; the other three in Hartmann & Kenyon (1995) have bolometric luminosities higher than likely for accreting VLM objects. Also, there is no evidence for any correlation between extinction and accretion activity in TTSs (not including embedded “Class I” objects). Furthermore, in the extinction-limited sample of Briceño et al. (2002), no correlation between extinction and (sub)stellar mass up to $`A_V4`$ was seen, thus there is no evidence that young brown dwarfs are selectively more extincted than higher-mass TTSs.

The physical basis for a mass vs. $`\dot{M}`$ relation is not entirely clear. In the standard model of viscous disk accretion (e.g., Lynden-Bell & Pringle 1974), the disk mass accretion rate is proportional to the disk surface density times the viscosity. The currently favored viscosity mechanism, the Balbus-Hawley instability (e.g., Balbus & Hawley 1998), requires sufficient ionization of the disk gas to operate efficiently. A likely source of ionization of T Tauri disks is X-ray radiation produced by magnetic activity (Glassgold, Najita, & Igea 1997; Igea & Glassgold 1999). Interestingly, recent studies have revealed a fairly strong correlation between X-ray luminosity and mass in young stellar objects (Feigelson et al. 2003; Flaccomio et al. 2003; Mokler & Stelzer 2002; Preibisch & Zinnecker 2002), where $`L_X`$ on average decreases by a factor of $`50100`$ from 1 to 0.1 M. It is possible that the lower X-ray flux produced by less massive young stars could result in smaller accretion rates by way of decreased disk ionization fractions, driving the correlation we find between $`\dot{M}`$ and mass. However, it is not clear whether the difference in $`L_X`$ over the mass range represented in Figure 8 can create enough of a change in disk ionization levels to produce a relation as steep as $`\dot{M}M^2`$. On the other hand, the correlation we have determined may be artificially steep, due to the possible selection effects described above. In any case, our results should provide interesting new motivation for continued theoretical work on disk accretion and evolution in young stars and VLM objects.

It should be noted that the mass-$`\dot{M}`$ relation we find is for objects with a relatively narrow range of ages, $`13`$ Myr. Thus, we cannot make any conclusions about whether such a relation holds at younger ages, i.e., at the embedded proto(sub)stellar phase. Most of the mass is probably accreted during the first few $`\times 10^5`$ years, either in the initial cloud collapse, or via FU Orionis- type outbursts, which may be different from typical viscous accretion and hence may be unrelated to the scenario we suggest in the previous paragraph. The accretion rates measured here are indicative of a more quiescent, steady accretion phase, and do not have much impact on the subsequent mass evolution of these objects.

We can comment on the fraction of accretors as a function of (sub)stellar age and mass. We have found 13 low-mass accretors with $`\dot{M}>10^{12}\mathrm{M}_{}\mathrm{yr}^1`$, all of which reside in the young (1-3 Myr) IC 348 or Taurus associations. Conversely, we have found no low mass accretors in the somewhat older (3-8 Myr) Sigma Ori or Sco Cen associations. Consistent with previous results on the evolution of disk accretion in higher mass, $`0.5M_{}`$, stars (e.g. Muzerolle et al 2000a; Muzerolle et al. 2001b), we find here an accretion timescale on the order of several Myr for objects near the substellar boundary. At mass $`0.5M_{}`$ and age $`2`$ Myr, the fraction of stars identified as accretors is about 50% (e.g., Kenyon & Hartmann 1995 for Taurus). Here, although we have found only about 30% of our total $`<0.2M_{}`$ sample to be accreting, $``$50% of the Taurus sample is accreting, albeit at lower rates.

We thank Russel White, Michael Meyer, and Kevin Luhman for helpful discussions and advice, and Suzan Edwards for a critical reading of the manuscript. This work was supported in part by NASA Origins of Solar Systems grant NAG5-9670 to NC and LH. CB received partial support from grant S1-2001001144 of the Fondo Nacional de Ciencia y Tecnología (FONACYT) of Venezuela.