text:
# Shake-up Processes in a Low-Density Two-Dimensional Electron Gas: Spin-Dependent Transitions to Higher Hole Landau Levels

## I  Introduction

Photoluminescence (PL) and photoluminescence excitation (PLE) spectroscopy of a two-dimensional electron gas (2DEG) in a perpendicular magnetic field proved effective in studying few- and many-body correlation effects (see, e.g., Refs. Kukushkin, ; Fink96, ; Yak97, ; Ossau2001, ; trions, ; Yusa2001, ; Schueller2002, ; Broocks2002, ; PRL2002, ; Sanvitto,  and references therein). In the 2DEG of density $`n_e`$ in quantizing magnetic fields with electron filling factor $`\nu _e=2\pi l_B^2n_e<1`$, the lowest energy absorption channel corresponds to creation of a bound $`e`$-$`h`$ pair, a neutral magnetoexciton in zero Landau levels (LL’s), $`\mathrm{𝑝ℎ𝑜𝑡𝑜𝑛}X_{00}`$; $`l_B=(\mathrm{}c/eB)^{1/2}`$ is the magnetic length. Transitions to higher LL’s, $`\mathrm{𝑝ℎ𝑜𝑡𝑜𝑛}X_{N_eN_h}`$, are also allowed when $`N_e=N_h`$. All these transitions are strong and gain strength with increasing magnetic field $`B`$. With increasing electron density $`n_e`$, these transitions become broadened by exciton-electron interactions, and optical processes associated with three-particle bound $`2e`$-$`h`$ states of charged excitons, or trions, $`X^{}`$, start to dominate the absorption Yak97 ; Ossau2001  and PL trions  spectra. Correlation effects in electron LL’s manifest themselves in the rich PL spectrum in strong magnetic fields. Yusa2001 ; Schueller2002 ; Broocks2002 ; PRL2002 ; Sanvitto ; Ashkinadze  The various theoretical approaches were used for studying three-particle bound trion $`X^{}`$ states Palacios96 ; Whit97 ; PRL2000 ; Wojs2000 ; Riva2001 ; PRB2002  and magneto-PL of the 2DEG Cooper97 ; Hawrylak97 ; Rashba2000  in strong magnetic fields. Despite a large amount of experimental and theoretical work, the spectra are not fully understood even in the low-density 2DEG limit. The role of electron-electron interactions, Sanvitto  disorder and localization, PRL2000  and the interplay between them Ashkinadze  are yet to be fully assessed.

Other optical manifestation of many-body effects in the interacting 2DEG are shake-up processes Fink96  in magneto-PL: After the recombination of an $`e`$-$`h`$ pair, one electron is excited to one of the higher LL’s. A closely related phenomenon in absorption (PLE) was also identified Yak97 ; Ossau2001  in low-density 2DEG systems: Here, an incident photon creates an $`e`$-$`h`$ pair and simultaneously excites one electron to higher LL’s. This process was calledYak97  “combined exciton-cyclotron resonance” (ExCR) and may be considered to be a shake-up process (or a radiative Auger process) in magneto-photoabsorption of the 2DEG. The shake-up phenomena in absorption and magneto-PL and the relation between them remain only partially understood. Fink96 ; Yak97 ; Ossau2001 ; PRB2001

At finite filling factors $`\nu _e`$, magneto-optical absorption involves many-particle “$`N`$-electron” initial states and “$`(N+1)`$ electrons + one valence band hole” final states. Generally, one has to deal with a many-body problem that is difficult to solve because the 2DEG in lowest LL’s is a strongly correlated system with no obvious small parameters. In the low-density limit, however, when filling factor $`\nu _e=2\pi l_B^2n_e1`$, distance between electrons is large (compared to $`l_B`$), and absorption can be considered to be an optical resonance $`e_0^{}+\mathrm{𝑝ℎ𝑜𝑡𝑜𝑛}2e`$-$`h`$, which involves an interacting charged three-particle system of two electrons and one hole in the final state. Yak97 ; Ossau2001 ; PRB2001  A similar situation occurs in the vicinity of integer filling factors. Cooper97 ; Hawrylak97 ; Rashba2000 ; PRL2002  In this limit, one needs (i) to study the spectra of charged three-particle Coulomb $`2e`$-$`h`$ problem in a magnetic field, and (ii) to find selection rules and dipole transition matrix elements. The techniques PRL2000 ; PRB2002 ; JETPL99  exploiting an exact symmetry, magnetic translations, magS  were applied in our previous paper PRB2001  to describe the spectra of charged $`e`$-$`h`$ states in higher LL’s and ExCR transitions in the high-field and low-density limit. In particular, a double-peak structure for the transitions to the first electron LL was predicted. PRB2001  The origin of the second peak in ExCR are transitions to the continuum formed by optically inactive excitons $`X_{N_eN_h}`$ ($`N_eN_h`$) which, somewhat unexpectedly, may have large ExCR oscillator strengths. Some experimental indications of the double-peak structure of the ExCR transitions to the first electron LL were reported recently. Ossau2001

In this work, which is an extension of Ref. PRB2001, , we study magneto-optical absorption in the low-density 2DEG in which the hole—in the presence of excess electrons—is excited to higher hole LL’s. In what follows, such shake-up absorption processes will be called hole-ExCR. The theory developed in Sect. II and Sect. III of this paper predicts the shape and fine structure of the high-energy tail of the main magneto-absorption peak of the 2DEG. No detailed experimental studies of this energy region (with sufficiently high spectral resolution) are known to the author.

The results of the present theory have also some relevance for the studies of a two-dimensional hole gas (2DHG). Experimentally, Coulomb correlation effects in magneto-PL spectra of the 2DHG were studied in a number of papers. Butov94 ; Kulik95 ; Ponomarev96 ; Volkov97 ; GlasbergPRB2001  Note that 2DHG magneto-PL spectra are often featureless and characterized by broad lines associated with the Zeeman split levels. This is mainly because of a large difference in electron and hole effective masses. For example, in GaAs quantum wells the electron and hole cyclotron energies are $`\mathrm{}\omega _{\mathrm{ce}}=1.7`$ meV/T and $`\mathrm{}\omega _{\mathrm{ch}}=0.2`$ meV/T, respectively, and differ by nearly an order of magnitude. GlasbergPRB2001  Recently, however, by examining the low energy tail of the magneto-PL spectra of the low-density 2DHG Glasberg et al. GlasbergPRB2001  succeeded in resolving several groups of recombination lines. Each group consists of several equidistant peaks directly related to the valence band LL’s (see also Refs. Butov94, ; Kulik95, ; Ponomarev96, ; Volkov97, ). In particular, in the dilute 2DHG the shake-up recombination lines of the positively charged exciton $`X^+`$, a complex consisting of two holes and one electron ($`2h`$-$`e`$), were observed. Quantitatively, the present theory cannot be directly applied for a description of the 2DHG magneto-optics. However, as is shown in Sec. IV, there are qualitative similarities between the 2DHG and 2DEG magneto-absorption spectra in the strong magnetic field and low-density limit. We also discuss shake-ups in magneto-PL of positive $`X^+`$ and negative $`X^{}`$ trions in Sec. V. A short summary is given in Sec. VI.

## II  Theory

The paper examines the hole-ExCR transitions $`e_0^{}+\mathrm{𝑝ℎ𝑜𝑡𝑜𝑛}2e`$-$`h`$ in the spin-polarized 2DEG in zero electron LL $`N_e=0`$$``$ with final three-particle states $`2e`$-$`h`$ that belong to the first hole LL $`N_h=1`$ with electrons still in zero LL’s (see Figs. 1, 2). The classification of states according to electron and hole LL’s is valid in sufficiently strong magnetic fields $`B`$ such that

$$\mathrm{}\omega _{\mathrm{ce}},\mathrm{}\omega _{\mathrm{ch}}>E_0=\sqrt{\frac{\pi }{2}}\frac{e^2}{ϵl_B},$$
(1)
where $`E_0`$ is the characteristic energy of the Coulomb interactions. Correlated electron-hole states can be obtained in this regime as an expansion in electron and hole LL’s. Before describing this procedure, some general features of the spectra that hold at arbitrary fields are discussed below.

Charged $`2e`$-$`h`$ states form families of degenerate states; each family is labeled by the index $`\nu `$ that plays a role of the principal quantum number and can be discrete (bound $`X^{}`$ states) or continuous (unbound $`X+e^{}`$ states forming a continuum). PRL2000  There is a macroscopic number of degenerate states in each family labeled by the discrete exact oscillator quantum number $`k`$ ($`k=0,1,\mathrm{}`$). This quantum number is associated with magnetic translations magS ; PRL2000 ; PRB2002  and describes the position of the center of the “cyclotron” orbit of the composite charged complex in $`B`$. Each family starts with its Parent State $`\mathrm{\Psi }_{k=0M_zS_eS_h\nu }^{}`$ that has $`k=0`$ and, roughly speaking, rotates in a magnetic field about the origin. All other daughter states, $`\mathrm{\Psi }_{kM_zkS_eS_h\nu }^{}`$ (with $`k=1,2,\mathrm{}`$) can be constructed iteratively PRL2000  out of the Parent State; note the $`M_zk`$ values. They have exactly the same physical properties and, because of this, will mostly be neglected in the discussion that follows. Accordingly, the quantum number $`k=0`$ of the Parent States will also be suppressed in notations. Importantly, the Parent State has the largest possible in the family value of the total angular momentum projection $`M_z`$. It is this value that enters JETPL99 ; PRL2000  into the exact optical selection, see Eqs. (3) and (12) below. In addition to orbital quantum numbers $`k`$ and $`M_z`$, the three-particle states are characterized by the total spin of two electrons, either $`S_e=0`$ (electron spin-singlet states, $`s`$) or $`S_e=1`$ (electron spin-triplet states, $`t`$) and by the spin state of the heavy hole $`S_h=\pm 3/2`$ ($``$, $``$). The heavy holes are described in this work in the effective mass approximation and the states of light holes are neglected.

The exact magneto-optical selection rules, in particular for ExCR, can be derived PRB2001  as follows: Interband transitions with $`e`$-$`h`$ pair creation are described by the luminescence operator $`\widehat{}_{\mathrm{PL}}=p_{\mathrm{cv}}𝑑𝒓\widehat{\mathrm{\Psi }}_e^{}(𝒓)\widehat{\mathrm{\Psi }}_h^{}(𝒓)`$, where $`p_{\mathrm{cv}}`$ is the interband momentum matrix element. When an $`e`$-$`h`$ pair is photocreated in the presence of the low-density 2DEG in $`N_e=0`$ LL, the dipole transition matrix element to the $`\nu `$-th family can be written as

$$D(\nu )=\mathrm{\Psi }_{M_zS_eS_h\nu }^{}|\widehat{}_{\mathrm{PL}}|e_0^{},$$
(2)
where only the final state three-particle correlations are taken into account. The combination of the two exact dipole selection rules, (i) conservation of the oscillator quantum number, $`\mathrm{\Delta }k=0`$ (the centers-of-rotation of charged complexes in the initial and final states are conserved) and (ii) no change in the total angular momentum for the envelope functions, $`\mathrm{\Delta }M_z=0`$, leads to a very simple but powerful result:

$$D(\nu )\delta _{N_e=0,M_z},$$
(3)
where $`M_z`$ is the angular momentum projection of the Parent State in the $`\nu `$-th family. PRL2000 ; PRL2002  The meaning of Eq. (3) is that in the magneto-absorption processes involving electrons from zero LL $`N_e=0`$, the achievable final $`2e`$$`h`$ states must have $`M_z=0`$. If the 2DEG is spin-up polarized ($``$, $`S_{\mathrm{ez}}=1/2`$), the photons of $`\sigma ^{}`$ and $`\sigma ^+`$ circular polarizations produce the electron triplet $`t`$ and singlet $`s`$ final states, respectively (Figs. 1, 2). The spin-dependence follows from the form of the luminescence operator in two circular polarizations,

$`\widehat{}_{\mathrm{PL}}(\sigma ^+)`$ $`=`$ $`p_{\mathrm{cv}}{\displaystyle 𝑑𝒓\widehat{\mathrm{\Psi }}_e^{}(𝒓)\widehat{\mathrm{\Psi }}_h^{}(𝒓)},`$ (4)
$`\widehat{}_{\mathrm{PL}}(\sigma ^{})`$ $`=`$ $`p_{\mathrm{cv}}{\displaystyle 𝑑𝒓\widehat{\mathrm{\Psi }}_e^{}(𝒓)\widehat{\mathrm{\Psi }}_h^{}(𝒓)}.`$ (5)
Interband transitions are described by polarization operators $`\pi _t(\omega )`$ and $`\pi _s(\omega )`$ that involve a single-particle electron propagator in $`N_e=0`$$``$ LL and three-particle propagators, which describe correlated negatively charged states $`\mathrm{\Psi }_{M_zS_eS_h\nu }^{}`$ with $`M_z=0`$ (see Figs. 1c, 2c). The electron-photon vertices, $`D`$, are determined by Eq. (2), and absorption at frequency $`\omega `$ is given by $`2\mathrm{I}\mathrm{m}\pi _{t(s)}(\omega )|D|^2\nu _e`$, where the dependence on electron filling factor $`\nu _e`$ comes from the electron propagator.

Let us now turn to the high-field regime (1) in the strictly-2D geometry. The final three-particle $`2e`$-$`h`$ states in ($`N_eN_h`$) electron and hole LL’s are found using an expansion in the complete orthonormal basis of squeezed oscillator states PRB2002  which simultaneously preserves both axial and magnetic translational symmetries

$$|N_RN_rN_h;kml.$$
(6)
The conserved oscillator quantum number is fixed automatically in basis (6) and equals $`k`$, while the total angular momentum projection must be fixed by imposing the condition $`M_z=N_R+N_rN_hkm+l`$. Here $`N_h`$ is the hole LL number, $`N_r`$ and $`N_R`$ are the LL numbers of the electrons relative and center-of-mass motions, respectively ($`N_r+N_R=N_e`$). The electron permutational symmetry requires that the two-electron relative angular momentum projection $`m_z=N_rm`$ must be even (odd) for electron singlet $`S_e=0`$ (triplet $`S_e=1`$) states; see Ref. PRB2002,  for more details. The quantum numbers $`k`$, $`M_z`$, $`S_e`$, and $`S_h`$ are exact so that the total Hamiltonian is block-diagonal in them. A weak exchange $`e`$-$`h`$ Coulomb interaction, which mixes different electron and hole spin states, is neglected in this work.

Because of the degeneracy in $`k`$, in what follows only Parent States $`\mathrm{\Psi }_{M_zs(t)\nu }^{}`$ are considered. The quantum numbers $`k=0`$ and $`S_h`$ are omitted for brevity and $`s`$ ($`t`$) denotes the singlet $`S_e=0`$ (triplet $`S_e=1`$) electron spin state. Neglecting mixing between LL’s (the high-field limit), the triplet $`2e`$-$`h`$ states in zero electron and first hole ($`N_eN_h`$)=(01) LL’s, $`\mathrm{\Psi }_{M_zt\nu }^{(01)}`$, can be obtained as the following expansion in (6)

$$\mathrm{\Psi }_{M_zt\nu }^{(01)}=\underset{l=0}{\overset{\mathrm{}}{}}\underset{m=0}{\overset{\mathrm{}}{}}\alpha _{M_zt\nu }(2m,l)|\mathrm{0\hspace{0.17em}0\hspace{0.17em}1};\mathrm{0\hspace{0.17em}\hspace{0.17em}2}ml,$$
(7)
where expansion coefficients $`\alpha _{M_zt\nu }(2m,l)\delta _{M_z,l2m1}`$ reflect the block-diagonal structure of the Hamiltonian. These are found by solving the corresponding secular equation. The Coulomb matrix elements of the total Hamiltonian are calculated analytically PRB2002  and the eigenspectra are obtained by numerical diagonalization of finite matrices of order $`5\times 10^2`$. Such finite-size calculations do not break either magnetic translational or axial symmetries and provide very high accuracy for bound $`X^{}`$ states. It is important to stress that such calculations are also capable of reproducing the structure of the three-particle continuum in some detail. PRB2002  The singlet states $`\mathrm{\Psi }_{M_zs\nu }^{(01)}`$ are found using the expansion

$$\mathrm{\Psi }_{M_zs\nu }^{(01)}=\underset{l=0}{\overset{\mathrm{}}{}}\underset{m=0}{\overset{\mathrm{}}{}}\beta _{M_zt\nu }(2m+1,l)|\mathrm{0\hspace{0.17em}0\hspace{0.17em}1};\mathrm{0\hspace{0.17em}\hspace{0.17em}2}m+1l$$
(8)
similar to (7) with expansion coefficients $`\beta _{M_zs\nu }(2m+1,l)\delta _{M_z,l2m2}`$.

The ExCR transitions are absent in the high-field limit (when no LL mixing allowed). PRB2001  The lowest order non-vanishing contribution to the dipole transition matrix element (2) appears when the triplet $`2e`$-$`h`$ states in zero LL’s are admixed to (7):

$$\delta \mathrm{\Psi }_{M_zt\nu }^{(01)}=\underset{l=0}{\overset{\mathrm{}}{}}\underset{m=0}{\overset{\mathrm{}}{}}\gamma _{M_zt\nu }(2m,l)|\mathrm{0\hspace{0.17em}0\hspace{0.17em}0};\mathrm{0\hspace{0.17em}\hspace{0.17em}2}ml.$$
(9)
The expansion coefficients $`\gamma _{M_zt\nu }(2m,l)\delta _{M_z,l2m}`$ are non-zero only because of the Coulomb mixing of different LL’s. Therefore, their characteristic smallness in high fields is $`\gamma _{M_zt\nu }E_0/\mathrm{}\omega _{\mathrm{ch}}B^{1/2}`$. A similar procedure is used to find $`\delta \mathrm{\Psi }_{M_zs\nu }^{(01)}`$ for the singlet states. Technical details of calculations of the dipole transition matrix elements have been described elsewhere. PRB2001

## III  Results and discussion

The triplet $`\mathrm{\Psi }_{M_zt\nu }^{(01)}`$ and singlet $`\mathrm{\Psi }_{M_zs\nu }^{(01)}`$ eigenspectra in the ($`N_eN_h`$)=(01) LL’s are presented in Fig. 3 and Fig. 4, respectively. The eigenstates were calculated in the strictly-2D high-field limit and their Coulomb energies are shown relative to free LL’s, $`\mathrm{}(\omega _{\mathrm{ce}}+3\omega _{\mathrm{ch}})/2`$. The shaded areas of width $`0.574E_0`$ in Figs. 3, 4 correspond to the three-particle continuum formed by the neutral magnetoexciton $`X_{01}`$ ($`e`$ in its zero and the hole $`h`$ in the first LL) with the second electron in a scattering state in zero LL. The lower continuum edge lies at the $`X_{01}`$ ground state energy $`0.574E_0`$. The density of $`X_{01}`$ states exhibits at this energy an inverse square-root van Hove singularity PRB2002  typical for 2D excitations whose dispersion has an extremum at a finite momentum. The presence of the singularity may lead to peculiarities in the $`e^{}`$ scattering on the neutral exciton $`X_{01}`$ and one might expect, in particular, formation of quasibound three-particle states and optical resonances (see below).

Filled dots in Figs. 3, 4 show positions of the excited bound three-particle states, denoted as $`(2e)`$-$`h`$, in which the hole is bound to the two-electron pair. These states originate from the excited states of two electrons that are always bound in 2D because of the confining effect of the magnetic field. FQHE  These truly bound three-particle states form discrete spectra and lie above free LL’s. They appear with relatively large positive values of $`M_z`$, when the hole can be at a sufficiently large distance from the two-electron pair. These states are not relevant for the considered optical transitions and will not be discussed further.

There are also many bound states of charged triplet $`X_{t01}^{}`$ and singlet $`X_{s01}^{}`$ magnetoexcitons in the first hole LL lying below the exciton continuum (open dots in Figs. 3, 4). This is in contrast to the situation in the first electron and zero hole LL’s, ($`N_eN_h`$)=(10), where only one bound state, the strongly bound triplet $`X_{t10}`$, exists in the 2D high-field limit. PRL2000  This difference is explained by the fact that in ($`N_eN_h`$)=(10) LL’s a strongly bound neutral magnetoexciton $`X_{00}`$ is formed. As a result, the three-particle continuum of $`X_{00}+e_1^{}`$ states is rather wide (width $`E_0`$), and essentially one low-lying bound triplet state $`X_{t10}^{}`$ can only be supported outside the continuum.

According to selection rule Eq. (3), the states that can optically be excited from zero electron LL $`N_e=0`$ must have $`M_z=0`$. There is one triplet and one singlet such states in the ($`N_eN_h`$)=(01) LL’s. These ExCR-active states are denoted by the filled triangles in Figs. 3, 4. The binding energies of the ExCR-active singlet $`X_{s01}^{}`$ and triplet $`X_{t01}^{}`$ states are $`0.009E_0`$ and $`0.024E_0`$, respectively. (Note that these are not the ground states and there are other low-lying triplet $`X_{t01}^{}`$ and singlet $`X_{s01}^{}`$ states that have larger binding energies.) Therefore, the ExCR transitions in the 2DEG to the first hole LL can terminate both in the bound states of charged magnetoexcitons, $`e_0^{}+\mathrm{𝑝ℎ𝑜𝑡𝑜𝑛}X_{s(t)01}^{}`$, and, also, in the continuum of scattering $`X_{01}+e_0^{}`$ states (see inset to Fig. 5). This is in contrast to the situation with transitions to the first electron LL, where transitions to the continuum are only allowed. PRB2001

The calculated dipole matrix elements of the hole-ExCR transitions in two circular polarizations $`\sigma ^+`$ and $`\sigma ^{}`$ that terminate, respectively, in the final singlet $`\mathrm{\Psi }_{M_z=0s\nu }^{(01)}`$ and triplet $`\mathrm{\Psi }_{M_z=0t\nu }^{(01)}`$ states are shown in Fig. 5 for a representative case $`\gamma =1`$. The dimensionless parameter $`\gamma =E_0/\mathrm{}\omega _{\mathrm{ch}}B^{1/2}`$ characterizes the magnetic field strength. The $`\sigma ^+`$ spectra at three other magnetic fields are shown in Fig. 6. All energies are given in units of (the magnetic-field dependent value) $`E_0`$ taken at $`\gamma =1`$. The spectra have been convoluted with a Gaussian of the $`0.015E_0`$ width that mimics inhomogeneous broadening. Note that both figures depict only the Coulomb interaction energies. The shown ExCR transitions require the extra photon energy $`\mathrm{}\omega _{\mathrm{ch}}`$ (blue shifted) relative to the fundamental band-gap absorption $`E_{\mathrm{gap}}(B)=E_g(0)+\mathrm{}(\omega _{\mathrm{ce}}+\omega _{\mathrm{ch}})/2`$ (with final states in the lowest LL’s) and they are also red-shifted $`0.5E_0B^{1/2}`$ due to the Coulomb interactions. The ExCR lines are also split in $`\sigma ^+`$ and $`\sigma ^{}`$ polarizations by the Zeeman energies $`\mu _B(g_eS_{\mathrm{ez}}+g_hS_{\mathrm{hz}})BB`$. The latter is schematically shown in the inset to Fig. 5.

It should be stressed that the ExCR transitions only arise because of LL mixing. PRB2001  The dependence of the dipole transition matrix elements on magnetic field is approximately $`|D|^2B^1`$ (Fig. 6), as expected. The intensity of transitions depends also on electron filling factor $`\nu _e`$ (see above). Therefore, the hole-ExCR transitions are suppressed in strong fields as $`\nu _e|D|^2n_el_B^2(E_0/\mathrm{}\omega _{\mathrm{ch}})^2B^2`$.

The ExCR-active trions $`X_{s01}^{}`$ and $`X_{t01}^{}`$ have rather small binding energies. Because of this, transitions to these bound states merge, in the presence of moderate broadening, with the strong transitions that terminate near the lower continuum edge. These are the main peaks in the hole-ExCR spectra (transitions 1 and 2 in Fig. 5). These peaks have intrinsic finite linewidths, in high fields $`0.1E_0`$, and have asymmetric lineshapes with high-energy tails. This is because the three-particle continuum is mainly responsible for these transitions. Generally, the hole-ExCR absorption band becomes wider and weaker with increasing magnetic field (Fig. 6).

The present theory predicts in high fields another strong hole-ExCR feature, the prominent and rather narrow higher-lying peak that can only be observed in the $`\sigma ^+`$ polarization (transition 3 in Fig. 5 and second peak in Fig. 6). This resonance does not show any appreciable change in its width with increasing $`B`$. This peak might be associated with formation of a quasi-bound three-particle resonance $`X_{s01}^{}`$, which lies within the three-particle continuum (see the inset to Fig. 5). The existence of the $`2e`$-$`h`$ resonances is physically plausible in higher LL’s because of the van Hove singularities in the density of states of the neutral 2D mangnetoexcitons (see above). A probability amplitude of finding all three particles in the same region of real space is large for a well-defined resonance. This may give large $`e`$-$`h`$ overlap PRB2001  and lead to unexpectedly large oscillator strengths of ExCR transitions to the continuum originating from optically inactive neutral magnetoexcitons $`X_{N_eN_h}`$ ($`N_eN_h`$).

The inset to Fig. 6 presents the data with no Gaussian convolution and shows some asymmetry in the shape of the optical resonance. This is due to the interaction between the quasi-discrete level and the continuum and is typical for a Fano-resonance. In the considered case of $`2e`$-$`h`$ complexes in the first hole LL, the interaction with the continuum appears to be rather weak and only leads to a slight asymmetry in the resonance shape (inset to Fig. 6). This should be contrasted with the situation in the first electron LL, where the interaction is much stronger, and the quasi-discrete state may show up as an anti-resonance in the intraband absorption spectra (see Fig. 2 in Ref. PRB2002, ). Details of formation of the three-particle resonances and the evolution of such features in the spectra with decreasing magnetic field are not currently understood.

The polarization dependence of the spectra may be explained by different $`e`$-$`e`$ correlations in the final spin-singlet ($`\sigma ^+`$ polarization) and spin-triplet ($`\sigma ^{}`$ polarization) three-particle $`2e`$-$`h`$ states; the former are characterized by much stronger $`e`$-$`e`$ repulsion. In the case considered this results in formation of the higher lying ExCR-active resonance $`X_{s01}^{}`$. Note also redistribution of the oscillator strength to the higher-energy peak 3 at the expense of the main peak 2 that occurs in $`\sigma ^+`$ polarization (see Fig. 5).

Note that one might expect some qualitative or even quantitative similarities between the spectra of the hole-ExCR interband transitions in the low-density 2DEG and the spectra of intraband internal transitions of isolated negative trions $`X^{}`$. This is because these transitions may have similar or exactly the same final three-particle correlated states in higher LL’s. Consider, for example, internal transitions from the ground singlet $`X_s^{}`$ state that has $`M_z=0`$. trions ; PRL2002 ; Whit97  When transitions to the first hole LL are induced by a FIR-photon of $`\sigma ^{}`$ polarization, the final states are $`\mathrm{\Psi }_{M_zs\nu }^{(01)}`$ with $`M_z=1`$. Note that in contrast to interband hole-ExCR, which is suppressed in strong magnetic fields, such internal transitions are hole-CR active, strong and gain strength JETPL99  with increasing $`B`$. This may be useful for experimental studies of the 2DEG using optically-detected cyclotron resonance techniques. PRL2002

## IV  Relation to shake-ups in 2DHG

In this section, magneto-optics of an interacting 2DHG is discussed on a qualitative level in the approximation neglecting the complex valence band structure. Generally, a difference in electron and hole effective masses determines different slopes of electron and hole LL’s and, at finite fields, leads to different interparticle correlations and different magneto-optical spectra in the interacting 2DEG and 2DHG.

In extremely strong magnetic fields $`\mathrm{}\omega _{\mathrm{ce}},\mathrm{}\omega _{\mathrm{ch}}E_0`$, however, LL mixing is suppressed and effective masses no longer determine interparticle correlations. In this limit one can find a very simple relation between the 2DEG and 2DHG magneto-optical spectra.

Consider, for example, a creation of an $`e`$-$`h`$ pair in the presence of low-density spin-polarized 2DHG in zero LL $`N_h=0`$$``$. The exact optical selection rule involving the initial hole state in $`N_h`$-th LL is (Fig. 7)

$$D(\nu )=\mathrm{\Psi }_{M_zS_eS_h\nu }^+|\widehat{}_{\mathrm{PL}}|h_{N_h}^+\delta _{N_h,M_z}.$$
(10)
Therefore, when the initial hole state belongs to $`N_h=0`$ LL, the final positively charged $`2h`$-$`e`$ states $`\mathrm{\Psi }_{M_zS_eS_h\nu }^+`$ must have $`M_z=0`$. In $`\sigma ^{}`$ polarization these are hole triplets ($``$$``$); a short-hand notation $`\mathrm{\Psi }_{M_z=0t\nu }^+`$ will be used in what follows. Notice now that the negatively charged electron triplet state $`\mathrm{\Psi }_{M_z=0t\nu }^{}`$ has—in the high field limit—exactly the same Coulomb interaction energy (see Appendix A). Also, the state $`\mathrm{\Psi }_{M_z=0t\nu }^{}`$ can be reached in transitions from low-density spin-polarized 2DEG in $`N_e=0`$$``$ LL in same $`\sigma ^{}`$ polarization. Therefore, the ExCR spectra (and, more generally, absorption spectra) of the spin-polarized ($``$) 2DHG can be obtained from the spin-polarized ($``$) 2DEG spectra in the same polarization by simply changing to the hole cyclotron energies $`\mathrm{}\omega _{\mathrm{ce}}\mathrm{}\omega _{\mathrm{ch}}`$ and making the appropriate changes in the Zeeman energies.

At finite fields, LL mixing (especially for the hole) breaks this electron-hole duality. In sufficiently strong fields, however, such that $`\mathrm{}\omega _{\mathrm{ch}}E_0`$, one may expect the indicated quantitative similarities between the spectra of the magnetically quantized 2DEG and 2DHG.

## V  Shake-ups in magneto-PL of isolated trions? Prohibited.

As we have seen, the ExCR (shake-up) processes in magneto-photoabsorption are allowed in the low-density 2DEG and 2DHG. Contrary to this, the exact selection rules prohibit shake-up processes in magneto-PL of isolated negative $`X^{}`$ and positive $`X^+`$ charged trions. Indeed, consider, for example, magneto-photoabsorption in a low-density 2DEG with an isolated electron in an arbitrary $`N_e`$-th LL in the initial state. According to the exact optical selection rule

$$D(\nu )=\mathrm{\Psi }_{M_zS_eS_h\nu }^{}|\widehat{}_{\mathrm{PL}}|e_{N_e}^{}\delta _{M_z,N_e},$$
(11)
the final negatively charged three-particle state $`\mathrm{\Psi }_{M_zS_e\nu }^{}`$ must have $`M_z=N_e`$ and may belong to different LL’s; this is generalization PRL2000  of Eq. (3).

Consider now magneto-PL from an isolated charged trion state $`\mathrm{\Psi }_{M_zS_eS_h\nu }^{}`$ that is not necessarily bound ($`X^{}`$), but may also belong to the continuum of scattering states ($`X+e^{}`$). The final state in magneto-PL is now a single electron (Fig. 8), and the exact selection rule

$$D(\nu )^{}=e_{N_e}^{}|\widehat{}_{\mathrm{PL}}^{}|\mathrm{\Psi }_{M_zS_eS_h\nu }^{}\delta _{N_e,M_z}$$
(12)
now states that the electron must reside in a single and given LL with number $`N_e=M_z`$. Therefore, no various LL’s are achievable in magneto-PL of a given trion state $`\mathrm{\Psi }_{M_zS_e\nu }^{}`$. This result is easily extended to positively charged trions $`X^+`$ and holds at arbitrary fields in systems with translational invariance and complex valence band structure (in the axial approximation). PRL2000 ; PRB2002

However, the shake-up processes in magneto-PL attributed to negatively Fink96  $`X^{}`$ and positively GlasbergPRB2001  $`X^+`$ charged trions are commonly observed in experiments in the dilute limit. This may be interpreted as an indication toward the relevance of the various scattering processes such as by disorder PRL2000 ; Ashkinadze  and/or the remaining 2DEG for $`X^{}`$ (or 2DHG for $`X^+`$). It was suggested recently Sanvitto  that the latter process, the scattering on low-density 2DEG, may also give rise to appreciable oscillator strength of the so-called “dark” triplet $`X^{}`$ state, Palacios96 ; PRL2000 ; Wojs2000  which is also commonly observed in experiments. trions ; Yusa2001 ; Schueller2002 ; Broocks2002

## VI  Summary

This paper theoretically considered hole-ExCR in a low-density 2DEG, an optical resonance in the 2DEG photoabsorption in which the hole is excited to higher hole Landau levels. This resonance may be observed in high magnetic fields as a fine structure in the high-energy tail of the main magneto-photoabsorption peak of the 2DEG. It has been shown, in particular, that the high-field hole-ExCR has different absorption patterns in two different circular polarizations $`\sigma ^\pm `$. The duality between correlated states in the 2DEG and 2DHG in the high magnetic field limit was indicated. This may be useful for magneto-optical studies of spin-dependent electron-hole correlations in 2D systems in strong magnetic fields.

###### Acknowledgements.

This work is supported in part by NSF grants DMR-0203560 and DMR-0224225 and by a research grant of CSUB.

*

## Appendix A  Duality between negatively and positively charged states

In the high-field limit (when no LL mixing is allowed) few-particle interacting states can be classified according to electron and hole LL’s ($`N_eN_h`$). In this limit, the Coulomb interaction energies of the positively charged ($`2h`$-$`e`$) hole triplet states $`\mathrm{\Psi }_{M_zt\nu }^{(N_eN_h)+}`$ and of the negatively charged ($`2e`$-$`h`$) electron triplet states $`\mathrm{\Psi }_{M_zt\nu }^{(N_hN_e)}`$ are the same; note the $`M_z`$ values and the order of LL labels. This holds in the effective mass approximation neglecting the complex valence band structure. In this approximation, the orbital (envelope) wavefunctions of the $`2e`$-$`h`$ and $`2h`$-$`e`$ states are related to each other by time reversal (are complex conjugate to each other). As a result, their Coulomb energies coincide. JETPL97  The difference in their total energies comes only from different cyclotron and spin Zeeman energies, which can easily be taken into account. The same, of course, holds for positively $`\mathrm{\Psi }_{M_zs\nu }^{(N_eN_h)+}`$ and negatively charged $`\mathrm{\Psi }_{M_zs\nu }^{(N_hN_e)}`$ singlet states whose interaction energies also coincide.

This duality between $`2h`$-$`e`$ and $`2e`$-$`h`$ states can be extended in high fields to other (charged and neutral) complexes consisting of other number of particles. JETPL97  It is known, for example, that dispersions of neutral magnetoexcitons $`X_{N_eN_h}`$ and $`X_{N_hN_e}`$ coincide in the high-field limit. L&L80  Another closely related example is the duality between the 2D correlated electron states that occur at electron filling factors $`\nu _e`$ and $`1\nu _e`$; this result is known in the theory of the fractional quantum Hall effect FQHE  as the electron-hole symmetry in the lowest LL.


format: text

text:
# STUDY OF THE ANGULAR RESOLUTION OF THE ARGO-YBJ DETECTOR

## 1 Introduction

The ARGO-YBJ detector, currently under construction at the Yangbajing Laboratory (P.R. China, 4300 m a.s.l.), is a full coverage array of dimensions $`74\times 78m^2`$ realized with a single layer of RPCs. The area surrounding this central detector (carpet), up to $``$ 100 $`\times `$ 110 $`m^2`$, is partially ($`50\%`$) instrumented with RPCs. The basic element is the logical pad ($`56\times 62cm^2`$) which defines the time and space granularity of the detector.

In a search for cosmic point-like sources the main problem is the background due to charged cosmic rays, so a good angular resolution (measurement of the arrival direction) is necessary. In order to determine the primary direction, the shower front has to be reconstructed from measurement of the time and position of the fired detectors. This calculation is usually performed for internal events only. In fact, the direction of showers with the core outside the detector is in general badly reconstructed due to the cone-like shape of the shower front and the unknown core position.

In Fig. 2 is shown an exaggerated shower with two scenarios: case A where the shower core strikes at the center of the ARGO-YBJ detector, and case B where the core hits outside the carpet. The primary energy in case B must be larger than that in case A to be able to trigger the array with the same detector response (number of fired pads). In case A, the direction obtained from the plane fitting is close to the true shower direction. On the other hand, in case B, the carpet sees a tilted portion of the same air shower, and in turn the fitted direction acquires a pointing error and deviates from the true shower direction.

In this paper we present different algorithms to reconstruct the shower direction. The calculations have been performed for different primary energies. The effect of a lead converter sheet on the detector is also investigated. The performance of a conical correction to the shower front with a suitable fixed cone slope is discussed.

## 2 Reconstruction procedure of the shower direction

The arrival time of the earliest particles provides the information concerning the time profile of the shower front. From MC simulations (see Fig. 15 of ref.) it appears that the shower front assumes a parabolic shape more pronounced for the electron component than for the photon component. Near the core charged particles form, at a first approximation, a plane disk.

The thickness of an EAS is defined by the distribution of arrival times of the shower particles which are delayed relative to the prompt particles which make up the shower front. Therefore, the shower thickness is determined by the lower-energy particles in a shower. This means that the shower thickness is the smallest near the core and increases with the increasing distance from the shower axis position. Hence, even in the case of vertical incidence the particles arrivals are not simultaneous.

In addition to the curvature, the combined effects of the finite size and time resolution of the detector make the shower apparently non-planar: since the shower plane has a thickness and the arrival time of the first particle is measured by the TDCs, the larger the number of particles, the earlier the measurement (on average).

The usual method for the reconstruction of the shower direction is a $`\chi ^2`$ fit to the recorded arrival times $`t_i`$ by minimization of

$$\chi ^2=\underset{i}{}w(ft_i)^2$$
(1)
where the sum includes all pads with a time signal. Usually the function $`f`$ describes a plane, a cone with a fixed cone slope, or a plane with curvature corrections depending on $`r`$ and $`N`$. The weights $`w`$ are generally chosen to be an empirical function of the number $`N`$ of particles registered in a detector.

An improvement to this scheme can be achieved by excluding from the analysis the time values belonging to the non-gaussian tails of the arrival time distributions with successive $`\chi ^2`$ minimizations for each shower.

In this paper we will investigate the performance of the $`\chi ^2`$ fit reconstructing the primary direction of the showers sampled by the ARGO-YBJ carpet by means of the following iterative procedure:

* Planar fit
In the first step we recover the shower direction cosines $`\{l_p,m_p\}`$ by means of an unweighted planar fit performed minimizing the function

$$\chi ^2(ns^2)=\frac{1}{c^2}\underset{i}{}\{lx_i+my_i+nz_i+c(t_it_0)\}^2$$
(2)
The sum includes all pads with a time signal $`t_i`$, $`c`$ is the velocity of light, $`(x_i,y_i,z_i)`$ are the coordinates of the central position of the $`i`$-th pad. The parameters of the fit are the time offset $`t_0`$ and the direction cosines $`l,m`$. After each minimization the time signals which deviate more than $`K\sigma `$ from the fitting function are excluded from further analysis and the fit is iterated until all times do not verify this condition or the remaining hits number is $`5`$ (in this case the event is discarded). Here $`\sigma `$ is the standard deviation of the time distribution around the fitted plane (i.e., $`\sigma =\sqrt{\frac{\chi ^2}{N3}}`$). In principle, the value of $`K`$ could depend on the features of the reconstructed showers as well as on the experimental condition (pad dimension, shaping of the signals, multihit capability, etc.).
* Conical correction
With this first determination of the arrival direction $`\{l_p,m_p\}`$ we calculate the conical correction

$$\delta t_i=\frac{\alpha }{c}R_i$$
(3)
where

$$R_i=\sqrt{(\mathrm{\Delta }x_i^2+\mathrm{\Delta }y_i^2)(\mathrm{\Delta }x_il_p+\mathrm{\Delta }y_im_p)^2}$$
(4)
and $`\mathrm{\Delta }x_i=x_ix_c`$, $`\mathrm{\Delta }y_i=y_iy_c`$ are the pad distances from the shower core position $`\{x_c,y_c\}`$. Then we correct the experimental time values $`t_i`$

$$t_it_i^{}=t_i\delta t_i(l_p,m_p)$$
(5)
In this approach the slope parameter $`\alpha `$ is not a fit parameter but it is fixed to a value that improves the angular resolution.
* New planar fit
With the time values $`t_i^{}`$ we reconstruct a new shower direction by means of a further planar fit.

In Fig. 2 the reduction of the non-gaussian tailes as a function of the successive plane fit iterations is shown. After few iterations the time residuals distribution assume a symmetrical shape (dotted histograms).

In order to optimize this procedure we have studied the resulting angular resolution as a function of the number of planar fit iterations and of the values for $`K`$. To estimate the performance of the algorithms we have used the $`\psi _{70}`$ parameter, defined as the value in the distribution of the angle between the true and the reconstructed directions within which 71.5 $`\%`$ of the events are contained. The angular resolution $`\sigma _\theta `$ is given by $`\sigma _\theta =1.58\psi _{70}`$.

In Fig. 4 the behaviour of the opening angle $`\psi _{70}`$ as a function of the number of iterations is shown, for different choices of $`K`$. The upper plot refers to 1 TeV $`\gamma `$-induced showers, the lower to 500 GeV $`\gamma `$-induced showers. In these calculations, all showers have been sampled with the core at the center of the carpet and we have selected events with a pad multiplicity $`N_{hit}`$ in the range 20 - 50. The same plots for higher multiplicities are shown is Fig. 4. No converter has been added on the RPCs.

From these figures it results:

* The angular resolution improves substantially with successive iterations. Indeed, the successive minimizations, as expected, reject the farthest temporal hits from the “theoretical plane” (see Fig. 2). This behaviour is qualitatively independent of the primary energy, the pad multiplicity and the values of $`K`$.
* After a few iterations the improvement stops and, for high multiplicities, the angular resolution worsens (see Fig. 4). This feature is due to the large fraction of discarded temporal hits lying far from the shower core. As a consequence, we must fix a suitable maximum number of iterations.

A systematic study performed for different energies and zenithal angles, sampling the shower cores over a wide area, led us to the following best tuning of the planar fit procedure: $`K=1.52`$ with the maximum number of iterations in the range 4 - 6.

A further improvement of the angular resolution can be obtained by fixing the cone slope of the successive conical correction to the value $`\alpha =0.03`$ $`ns/m`$. As an example, in Fig. 6 we show the parameter $`\psi _{70}`$ as a function of pad multiplicity for different values of the cone slope $`\alpha `$. The calculations refer to vertical 1 TeV $`\gamma `$-induced showers with the core at the center of the carpet. As can be seen, the value $`\alpha =0.03`$ $`ns/m`$ improves the angular resolution for all the investigated multiplicities (see Fig. 6).

We also investigated the performance of a conical fit in which the cone slope is a free parameter of the minimization. The angular resolution calculated with fixed $`\alpha `$ gives better results because the minimization tends to adapte the cone slope to the temporal hits, thus introducing a larger error in the direction recontruction for events with the core outside the carpet or in its boundaries.

### 2.1 The effect of a lead converter on the angular resolution

The consequences of placing a thin sheet of converter above the detector are, qualitatively: (1) absorbtion of low energy electrons (and photons) which no longer contribute to the time signal; (2) multiplication process of high-energy electrons and photons which produces a signal enhancement. The lead converter thus reduces the temporal fluctuations: the contributions gained are concentrated near the ideal time profile since the high energy particles travel near the shower front while those lost tend to lag far behind.

In order to study the effect of the converter on the angular resolution we have simulated, via a GEANT3-based code, a 0.5 cm lead sheet located 5 cm above the RPCs. The results discussed in the previous section are qualitatively reproduced also when the converter is added for what concerns the dependence of the angular resolution on the number of iterations of the planar fit. As expected, the improvement of the showers temporal profile implies a fewer number of iterations in the iterative procedure: the best tuning of the planar fit is achieved after only 3 iterations, with K = 1.5 .

In Fig. 8 we show the effect of a 0.5 cm lead sheet on the angular resolution by comparing the $`\psi _{70}`$ parameter as a function of the fired pads with and without the lead. The calculations refer to 1 TeV $`\gamma `$-induced showers with a fixed zenithal angle $`\theta =20^{}`$ and with the core sampled on a 50 $`\times `$ 50 m2 area. In the upper plot we show the values calculated after the best tuned iterative plane fit: 5 iterations without lead and 3 iterations with the converter. In the lower part of the figure $`\psi _{70}`$ is calculated after the conical correction. The improvement is clearly evident.

## 3 Angular resolution of the ARGO-YBJ detector

In oder to study the angular resolution of the ARGO-YBJ detector we have simulated, via the Corsika code, $`\gamma `$-induced showers with a Crab-like spectrum ($`E^{2.5}`$) and proton events with $`E^{2.75}`$, both ranging from 100 GeV to 50 TeV. The $`\gamma `$-rays have been simulated for different zenith angles ($`<40^{}`$), following the daily path of the source in the sky. The detector response has been simulated via a GEANT3-based code. A 0.5 cm lead sheet have been considered 5 cm above the RPCs. The core positions have been randomly sampled in an area, energy-dependent, large up to 800 $`\times `$ 800 m2 and centered on the detector.

In Fig. 8 we compare the dependence of the opening angle $`\psi _{70}`$ on pad multiplicity for the different algorithms discussed above. The $`\gamma `$-induced showers are reconstructed inside a fiducial area similar to the carpet dimensions ($`A_{fid}=80\times 80`$ $`m^2`$). The shower core positions are reconstructed by means of the Maximum Likelihood Method applied on the RPCs (the so-called ’LLF2 method’). The conical correction needs a good shower core position resolution, so for very low multiplicities the planar fit is more performant because of the high contamination of external events. However, when suitable “cuts” are applied and a high fraction of external showers is rejected, the conical correction improves the planar fit.

In Fig. 10 we show the parameter $`\psi _{70}`$ as a function of pad multiplicity for $`\gamma `$ showers reconstructed inside $`A_{fid}`$ with and without the application of the cuts. For comparison, the $`\psi _{70}`$ values for the events with the sampled core truly inside $`A_{fid}`$ are also plotted (circles). As can be seen, for low multiplicities ($`N_{hit}<100`$), the rejection of external events is crucial to improve the angular resolution.

In Fig. 10 a comparison between the angular resolutions obtained after the cuts for gamma and proton showers is shown. Due to the larger lateral spread of particles in proton showers, the angular resolution for protons is worse than that for photons.

## 4 Conclusions

In this paper we described an algorithm optimized to reconstruct the direction of the showers detected by the ARGO-YBJ experiment. A conical correction with a cone slope parameter fixed to the value $`\alpha =0.03`$ $`ns/m`$ is proposed and discussed. The resulting opening angle $`\psi _{70}`$ is better than $`0.3^{}`$ for $`\gamma `$-induced events which fire more than 600 pads on the carpet.

For $`N_{hit}>100`$ (when the dependence on the fiducial area is reduced) the angular resolution of the ARGO-YBJ experiment, for $`\gamma `$-induced events with a Crab-like energy spectra, can be described by the equation $`\psi _{70}=0.046+6.666/\sqrt{N_{hit}}`$.


text:
# Chondrule Formation and Protoplanetary Disk Heating by Current Sheets in Non-Ideal Magnetohydrodynamic Turbulence

## 1 Introduction

Chondrules are millimeter-sized silicate spheres with igneous textures embedded in the most common type of meteorites, comprising 50–80 % of their mass. Their formation is a long-standing mystery in solar system formation. Radioisotope datings place them among the first solids to have condensed out of the slowly cooling solar nebula. As such, they represent the remains of early materials that later formed the terrestrial planets. Their nearly spherical shapes and glassy compositions indicate that they cooled within a few hours after being molten, but not in the few seconds that objects that size would cool in free space. The associated heating event, whatever the origin, must have been transient and localized as evidenced by the inferred rapid heating and cooling rates, and the observed retention of volatiles such as sulfur (which vaporizes by prolonged heating above 650 K). To explain the high abundance of chondrules, it must have been a widespread phenomenon. Hence, encrypted in chondrules are important clues to physical conditions that prevailed in the early stages of solar system formation. Some even argue that taking the form of millimeter-sized molten spheres is a necessary intermediate step for initially micron-sized dust grains to grow to kilometer-sized objects (Hewins, 1997).

A successful theory of chondrule formation should at least meet the following constraints: (1) the extremely short heating timescale of roughly seconds to minutes (Connolly & Love, 1998), (2) the cooling rate of $`10^2`$$`10^3`$ K hr-1 (Radomsky & Hewins, 1990; Lofgren & Lanier, 1990), and (3) the narrow range of sizes (0.1–1$`\mathrm{mm}`$), inferred from laboratory analyses of textures and compositions as well as from attempts to reproduce chondrules. In addition, many chondrules show evidences of multiple heating events, so the heating event must have occurred more than once. Since the initial discovery 200 years ago (Howard, 1802), their formation process has received considerable attention in both meteoritics and solar system astronomy, yet it remains elusive (Wood, 1996). Proposed mechanisms in the past (see the summary by Boss, 1996, for the pros and cons) include nebula shocks, lightning, and an X-wind. None of these mechanisms have really been widely accepted (Hewins, 1996), although recent work on shock heating by Desch & Connolly (2002) has gained wide attention.

Infrared observations and subsequent theoretical modeling of T Tauri stars revealed the existence of circumstellar (protoplanetary) disks (Mendoza, 1968; Shu, Adams, & Lizano, 1987; Adams, Lada, & Shu, 1987; Beckwith et al., 1990), many of which are actively accreting (see, e.g. the review by Calvet, Hartmann, & Strom, 2000). These systems are generally thought to be analogs of the early solar nebula. On the theoretical side, magnetorotational instability (MRI) provided a solution to the long-standing problem of identifying the cause of anomalous viscosity, which allows mass accretion of the disk (Balbus & Hawley, 1991; Hawley & Balbus, 1991, collectively referred to as BH91 hereafter). When applied to partially ionized disks, it remains an efficient mechanism of generating turbulence as long as ions are sufficiently well-coupled to neutrals, i.e., if ion-neutral collision rate far exceeds the orbital rotation rate: $`\nu _{in}\mathrm{\Omega }`$ (Blaes & Balbus, 1994; Mac Low et al., 1995; Hawley & Stone, 1998). In the case of protoplanetary disks, however, a large portion, especially near the disk midplane, is associated with relatively low ionization fraction ($`10^{12}`$), and hence thought to be magnetically dead (Gammie, 1996; Glassgold, Najita, & Igea, 1997).

In partially ionized disks, ambipolar diffusion arises because magnetic forces act on the charged component, i.e. ions, only (Mestel & Spitzer, 1956). As the ions drift through neutrals, the magnetic field lines drag the ions with them, steepening field gradients near magnetic nulls stretched by shear flows (Brandenburg & Zweibel, 1994, hereafter BZ94). Consequently, sheets of high electric current form, as was confirmed by numerical simulations (Mac Low et al., 1995; Suzuki & Sakai, 1996). BZ94 originally proposed this as a mechanism for rapid magnetic reconnection in the interstellar medium. Ambipolar diffusion is important in another context: as a heating source (Scalo, 1977; Goldsmith & Langer, 1978; Zweibel & Josafatsson, 1983; Elmegreen, 1985). Recently, using three-dimensional simulations of turbulent, magnetized molecular clouds, Padoan, Zweibel, & Nordlund (2000) pointed out that ambipolar diffusion heating is a significant source of heating in molecular clouds and that it can be larger than the cosmic ray heating. In this paper, we examine the possibility that similar physics operates in protoplanetary disks, such as in the solar nebula that eventually became our solar system.

The purpose of this paper is twofold. First, we show that current sheets are generated in magnetically active regions within protoplanetary disks. We study the properties of these current sheets and determine when and where the associated heating becomes significant. Second, we apply this formalism to the following two unresolved issues: the heating mechanism that melted chondrules in the early solar nebula and the near-infrared excesses observed in spectral energy distributions (SEDs) of young stellar objects.

We later demonstrate that the heating of dust grains by current sheets is ineffective for a standard model of the early solar nebula in which the dust-to-gas ratio $`\zeta 0.01`$. However, at lower gas densities and much higher dust-to-gas ratios than those usually assumed in the standard model, these current sheets are found to reach high enough temperatures to melt millimeter-sized dust grains. We propose that these unusual conditions occur near the midplane in late stages of solar nebula evolution, when the gas surface density has decreased due to various mechanisms discussed below, and dust grains have gravitationally settled. The surfaces of this dust-enriched midplane layer are then sufficiently dense and ionized to host thin ($`l10^2\mathrm{km}`$), hot ($`T>1600\mathrm{K}`$) current sheets. Circulation of dust to the surface by turbulence within the dust layer is invoked to explain the observed high abundance of chondrules.

The basic features of our current sheet model are set forth in § 2. In § 3.1 we use this model to investigate a scenario in which chondrules are heated inside current sheets, solving numerically the coupled equations of energy transfer for gas and dust particles. The results are presented in § 3.2. We argue, in § 3.3, that in late stages of solar nebula evolution current sheets capable of melting most of the millimeter-sized dust particles form in the dust-rich midplane. In § 3.4 we discuss possible observational signatures of current sheets. Finally, in § 4 we list limitations of our model.

## 2 Current Sheet Model

### 2.1 Basic Equations

The magnetic field evolves according to the induction equation,

$$\frac{𝑩}{t}=\mathbf{}\times (𝒗_e\times 𝑩\eta \mathbf{}\times 𝑩),$$
(1)
where $`𝒗_e`$ denotes the electron velocity and $`\eta `$ is the associated resistivity given by

$$\eta =\frac{c^2}{4\pi \sigma _e}=\frac{234}{\chi _e}T_g^{1/2}\mathrm{cm}^2\mathrm{s}^1$$
(2)
(Krall & Trivelpiece, 1973). Here, $`c`$ is the speed of light, $`\chi _e=n_e/n_n`$ is the electron fraction, and $`T_g`$ is the ambient gas temperature. For the last equality in Eqn. (2), we adopt the value of the electrical conductivity, $`\sigma _e`$, from Draine, Roberge, & Dalgarno (1983) and Blaes & Balbus (1994). Subscripts $`e`$, $`i`$, and $`n`$ denote quantities pertaining to electron, ion, and neutral components, respectively. Other quantities of relevance in this paper are the mean mass per particle of $`m_n=2.33m_p`$, where $`m_p`$ is the proton mass, and the mean ion mass of $`m_i=39m_p`$ at densities where ionized potassium, K+, is generally the dominant ion species (Balbus & Terquem, 2001).

In the ideal MHD limit, the electron velocity is indistinguishable from the center-of-mass velocity of the fluid. This approximation breaks down in two ways as $`\chi _e`$ decreases. First, ions slip through neutrals, and for low enough $`\chi _e`$, when ion pressure and inertia can be neglected, the strong coupling approximation can be used to compute the ambipolar drift velocity by balancing the Lorentz force $`𝑱\times 𝑩/c`$ acting on the ions against the frictional drag $`\rho _i\nu _{in}(𝒗_i𝒗_n)`$ arising from collisions with neutrals. Second, when ions and electrons move along a direction perpendicular to a magnetic field, they experience the Lorentz force acting normal to both directions but they drift away in two opposite directions due to the different signs of their electric charges. This causes the Hall effect. The electron velocity in a weakly ionized medium is then

$$𝒗_e=𝒗+(𝒗_i𝒗)+(𝒗_e𝒗_i)=𝒗+\frac{𝑱\times 𝑩}{c\rho _i\nu _{in}}\frac{𝑱}{n_ee},$$
(3)
where $`\rho _i=m_in_i`$ and $`\nu _{in}=\gamma \rho _n`$ are the ion mass density and ion-neutral collision rate. The value of the drag coefficient, $`\gamma =2.75\times 10^{13}`$ cm3 s-1 g-1 (Balbus & Terquem, 2001). Since $`\chi _e1`$ in our case (see § 2.2), the neutral velocity is the same as the center-of-mass velocity, $`𝒗`$. The current density has the usual form $`𝑱=(c/4\pi )\mathbf{}\times 𝑩`$. We plug Eqn. (3) into Eqn. (1) and consider the case in which the currents and fields are perpendicular ($`𝑱𝑩=0`$), so the induction equation becomes

$$\frac{𝑩}{t}=\mathbf{}\times (𝒗\times 𝑩)\mathbf{}\times \left[\left(\frac{B^2}{4\pi \rho _i\nu _{in}}+\eta \right)\frac{4\pi 𝑱}{c}\right].$$
(4)
The Hall term does not affect the time evolution of $`𝑩`$, since $`\mathbf{}\times (𝑱\times 𝑩)\mathbf{}\times \mathbf{}B^2`$ vanishes. Thus, ambipolar diffusion acts as an additional field-dependent resistivity. We define the ambipolar diffusion coefficient as

$`\lambda _{AD}`$ $``$ $`{\displaystyle \frac{B^2}{4\pi \rho _i\nu _{in}}}`$
$``$ $`1.1\times 10^{17}B_G^2\left({\displaystyle \frac{n_n}{10^{13}\mathrm{cm}^3}}\right)^2\left({\displaystyle \frac{\chi _i}{10^{12}}}\right)^1\mathrm{cm}^2\mathrm{s}^1,`$
where $`B_G`$ means $`B`$ in gauss.

For our model of the protoplanetary disk, we take the minimum-mass solar nebula of Hayashi (1981). It is characterized by the surface mass density $`\mathrm{\Sigma }=\mathrm{\Sigma }_0(R/3\mathrm{AU})^{3/2}`$ with $`\mathrm{\Sigma }_0=3.3\times 10^2`$ g cm-2, where $`R`$ is the radial distance from the central star, and the temperature profile given by $`T_g=1.6\times 10^2(R/3\mathrm{AU})^{1/2}`$. We expect $`\mathrm{\Sigma }_0`$ to decrease in time as the nebular gas clears out in the late stage of protoplanetary disk evolution. Hence we use the form $`\mathrm{\Sigma }_0=3.3\times 10^2f_\mathrm{\Sigma }`$ g cm-2, where $`f_\mathrm{\Sigma }`$ is an arbitrary parameter less than or equal to unity. If we take the gas to be isothermal, the vertical density profile in hydrostatic equilibrium has a Gaussian form (Kenyon & Hartmann, 1987; Chiang & Goldreich, 1997, hereafter CG97)

$$n_g=n_0\mathrm{exp}\left(\frac{z^2}{2h^2}\right),$$
(6)
where $`n_0=\mathrm{\Sigma }/(\sqrt{2\pi }m_nh)=1.4\times 10^{13}f_\mathrm{\Sigma }(R/3\mathrm{AU})^{11/4}`$ cm-3 is the number density of gas in the disk midplane, $`z`$ is the vertical distance measured from the midplane, and $`hc_s/\mathrm{\Omega }`$ is the scale height. Here, as usual, $`c_s=(kT_g/m_n)^{1/2}`$ is the sound speed and $`\mathrm{\Omega }`$ is the Keplerian rotation rate. With these values, we can roughly estimate the importance of ambipolar diffusion by taking the ratio of the two coefficients in Eqn. (4):

$`{\displaystyle \frac{\lambda _{AD}}{\eta }}`$ $`=`$ $`{\displaystyle \frac{B^2}{4\pi \eta \rho _i\nu _{in}}}`$
$``$ $`20f_\mathrm{\Sigma }^2B_G^2\left({\displaystyle \frac{R}{3\mathrm{AU}}}\right)^{23/4}\mathrm{exp}\left({\displaystyle \frac{z^2}{h^2}}\right).`$
We note that (1) with a magnetic field sufficiently amplified (up to $``$1 G), ambipolar diffusion is at least as important as Ohmic dissipation and (2) the ratio has strong spatial dependence within the disk. Remnant magnetism in meteorites suggests that $`B`$ 1 G was typical in the solar nebula. If an efficient dynamo is in action, the ambipolar diffusion term would generally be important in the outer ($`R`$ 1 AU) region of the disk. In the inner region, for a similar value of $`B`$ and $`f_\mathrm{\Sigma }=1`$, it is nonnegligible only a few scale heights above and below the midplane, where gas density is much lower. It generally becomes increasingly important as $`f_\mathrm{\Sigma }`$ decreases.

### 2.2 Magnetorotational Instability in Protoplanetary Disks

Compared to other astrophysical systems, protoplanetary disks are unusually cold and dense and are therefore associated with very low ionization fractions. Ionization fraction is a key parameter in determining whether a disk will be turbulent or not, since the MRI is confined to regions where the ionization fraction is high enough to ensure good coupling between ions and neutrals. This occurs when the ion-neutral collision rate ($`\nu _{in}`$) is significantly greater (by 3–100 at least) than the rotation rate ($`\mathrm{\Omega }`$) (Blaes & Balbus, 1994; Mac Low et al., 1995; Hawley & Stone, 1998). The critical ionization fraction is calculated to be

$$\chi _{crit}=\frac{\mathrm{\Omega }}{m_in_n\gamma }=2\times 10^{12}f_\mathrm{\Sigma }^1\left(\frac{R}{3\mathrm{AU}}\right)^{5/4}\mathrm{exp}\left(\frac{z^2}{2h^2}\right)$$
(8)
(e.g., Balbus & Hawley, 1998). Interestingly, this value is only slightly higher than the values of $`10^{12}`$$`10^{13}`$ estimated in the bulk of protoplanetary disks (Gammie, 1996; Glassgold, Najita, & Igea, 1997; Sano et al., 2000).

Many numerical simulations have demonstrated that MRI-induced velocity shears amplify the magnetic field (see, e.g., Hawley & Balbus, 1991; Brandenburg et al., 1995; Hawley, Gammie, & Balbus, 1995). In our model, perturbations in the magnetic field are also assumed to be generated by the MRI. Considering the low values of $`\chi _i`$ estimated in protoplanetary disks, it is initially unclear whether this assumption is valid or not. Therefore we ask: which region of the solar nebula, if any, is prone to the MRI? (Previous numerical simulations assumed ionization fraction $`\chi _i\chi _{crit}`$.) This question has been investigated by several authors (Umebayashi, 1983; Sano et al., 2000; Fromang, Terquem, & Balbus, 2002). Their results are as follows. The ionization fraction is sufficiently high in the innermost region of the disk up to a few tenths of an AU, where the ambient gas temperature is high enough ($`T_g10^3`$ K) for collisional ionization to be effective (where it dominates over cosmic ray and radioactive element ionizations; Umebayashi, 1983). Also, cosmic rays and stellar X-ray radiation ionize the surface layers (with thickness given by the stopping depth of $`\mathrm{\Sigma }_a10^2`$ g/cm2) at $`R10`$ AU and the entire disk at $`R10`$ AU, where the surface density of the disk falls below $`\mathrm{\Sigma }_a`$. These, then, are the regions where the disk is turbulent, hence effectively transferring angular momentum and also actively accreting. On the other hand, the region near the disk midplane, where terrestrial planets and the asteroid belt now lie, is generally thought to contain no significant heating source and be quiescent, i.e. magnetically dead (radioactive heating is too inefficient to maintain $`\chi _i>\chi _{crit}`$; see Gammie, 1996). These considerations led to the layered accretion model of Gammie (1996; see also Glassgold, Najita, & Igea, 1997).

There are two length scales, other than the scale height $`h`$, of relevance to this discussion. The wavelength of perturbation most unstable to the MRI is only slightly larger than the critical wavelength at which the instability sets in (BH91) and is given by

$$\lambda _{BH}=\frac{2\pi }{\sqrt{3}}\frac{v_A}{\mathrm{\Omega }}0.23f_\mathrm{\Sigma }^{1/2}B_G\left(\frac{R}{3\mathrm{AU}}\right)^{23/8}\mathrm{exp}\left(\frac{z^2}{4h^2}\right)\mathrm{AU}$$
(9)
(BH91). Other processes affect the growth of perturbations: in dense parts of the disk, collisions with neutrals hinder the charged components from drifting with the magnetic field. In this case, Ohmic dissipation acts to stabilize wavelengths shorter than the resistive diffusion scale $`\lambda _\mathrm{\Omega }2\pi \eta /v_A`$, where $`v_A=B/(4\pi \rho _n)^{1/2}`$ is the Alfven speed. In other words, the critical wavelength at which MRI sets in is given by $`\lambda _{crit}=max(\lambda _{BH},\lambda _\mathrm{\Omega })`$ (Blaes & Balbus, 1994; Sano et al., 2000). When $`\lambda _\mathrm{\Omega }>\lambda _{BH}`$, the perturbation growth becomes much slower. An additional constraint enters when $`\lambda _{crit}`$ becomes comparable to the scale height $`h`$; at this point, saturation of the MRI is thought to occur (Balbus & Hawley, 1998).

Sano et al. (2000) investigated where the MRI operates in protoplanetary disks by including a detailed chemical reaction scheme and the effect of the recombination of ions and electrons on grain surfaces. They showed that for the minimum-mass solar nebula ($`f_\mathrm{\Sigma }=1`$) the 2–3 AU region of the disk would be stabilized by Ohmic dissipation. However, as the surface density parameter, $`f_\mathrm{\Sigma }`$, is lowered, the unstable region expands (see their Figure 8). Also, as the disk evolves, dust grains settle to a thin layer in the midplane, and so the grain abundance in the bulk of the nebula gradually decreases. For this reason, Sano et al. (2000) introduced a dust depletion factor, $`f_d`$, which measures the abundance of dust grains relative to that in molecular clouds. Thus a low value of $`f_d`$ may correspond to an old, evolved disk. They demonstrated that as $`f_d`$ decreases for a given surface density, the stable region (the “dead zone”) shrank in size, and for $`f_d10^2`$, the entire disk becomes unstable for $`R>3`$ AU (see their Figure 11). This happens because the recombination rate on grain surfaces decreases. They found a similar trend as grain growth occurs. From these results, we conclude that the region at $`R3`$ AU will become magnetically active at a sufficiently late stage of the evolution. We use and discuss this result further in § 3. Figure 1 compares the three relevant length scales at an early (a) and a late stage (b) of solar nebula evolution. The region of interest where magnetic fields amplify due to the MRI extends from 2 AU to 4 AU in (a), while with a reduced gas density, the active region moves inward to 1–2 AU (b).

### 2.3 Formation of Current Sheets

In this and the following sections, we consider a simple, one-dimensional geometry illustrated in Figure 2. Also shown are the directions of various terms in Eqns. (3) and (4). Consider first a case in which a magnetic null is present. Near the null, a magnetic pressure gradient, $`\mathbf{}(B^2/8\pi )`$, acts on ions and pushes the field lines into the null point in both directions. The ion-neutral drift is then directed toward the null. Along a gradient in the magnetic field, this nonlinear term produces faster field diffusion in regions of stronger field, driving field into regions of weaker field and steepening the gradient, rather than spreading it out as the linear Ohmic term does. Ultimately, the gradient sharpens towards a singularity, marked by a sheet of high electric current. This steepening of magnetic field profile (BZ94) occurs within a timescale of $`\tau _{AD}L_{BH}^2/(2\lambda _{AD})=(\pi ^2/\mathrm{\hspace{0.17em}6})(\nu _{in}/\mathrm{\Omega }^2)`$, calculated to be $`\pi /12`$ of the orbital period when $`\chi _i=\chi _{crit}`$. The factor of 1/2 in $`\tau _{AD}`$ takes into account the fact that the characteristic length decreases as steepening proceeds.

Remarkably, this sharpening has been shown to occur even in the absence of nulls due to a variety of processes including a rotating eddy (BZ94; Zweibel & Brandenburg, 1997), accretion disks sustaining the MRI (Mac Low et al., 1995), large amplitude Alfven waves (Suzuki & Sakai, 1996), and the nonlinear stage of Wardle C-shocks (Mac Low & Smith, 1997). As previously mentioned, we take the example of magnetohydrodynamic turbulence driven by the MRI. In this case, velocity shears produce magnetic nulls and current singularities . Stretched by shear flows, magnetic nulls develop where the magnetic field changes directions. We therefore crudely estimate that the separation between two adjacent magnetic nulls $`L_{BH}=\lambda _{BH}/2`$.

We note that ambipolar diffusion allows the field to slip through the neutrals and tends to weaken the MRI. For strong enough diffusion, i.e. when the ion-neutral collision rate ($`\nu _{in}`$) is much higher than the growth rate of the MRI ($`\mathrm{\Omega }`$), the instability is entirely prevented (Mac Low et al., 1995; Hawley & Stone, 1998). On the other hand, ambipolar diffusion is required to generate current sheets. For our model, we use parameters on the interface between the two regimes. This choice is justified by two-dimensional simulations performed by Mac Low et al. (1995), which showed that for drag coefficients $`\gamma `$ a few times the minimum value that allows the instability to develop, $`\gamma _{min}\mathrm{\Omega }/(m_in_n\chi _{crit})`$, both the instability and ambipolar diffusion sharpening occur. In other words, we require $`1\gamma /\gamma _{min}<100`$. Alternatively, magnetic field amplification could have preceeded the epoch of current sheet formation.

We now make quantitative estimates of the properties of current sheets. BZ94 showed that, in the absence of resistivity, the magnetic field profile relaxes to $`Bx^{1/3}`$, and the current density $`JB/xx^{2/3}`$ becomes singular at the origin. In reality, both the finite resistivity and the ion pressure act to remove this singularity (Brandenburg & Zweibel, 1995). Under the physical conditions of protoplanetary disks, ion pressure can be neglected (Heitsch & Zweibel, 2003). The induction term in Eqn. (4), $`\mathbf{}\times (𝒗\times 𝑩)`$, becomes progressively less important as the field steepens; $`I/OvL/\eta `$ where $`I`$ and $`O`$ denote the magnitudes of the induction and Ohmic dissipation term respectively, and is therefore ignored. The quasi-steady-state magnetic field profile (“Quasi” since the magnetic field profile lasts only as long as the velocity shear is present, for about a dynamical time) is calculated by balancing the Ohmic dissipation term against the ambipolar diffusion term in Eqn. (4):

$$(B_0^2+B_y^2)\frac{d^2B_y}{dx^2}+2B_y\left(\frac{dB_y}{dx}\right)^2=0.$$
(10)
Here, for simplicity, we use the definition

$$B_0(4\pi \eta \rho _i\nu _{in})^{1/2}0.26n_{g13}T_{g3}^{1/4}\mathrm{G},$$
(11)
which measures the relative importance of Ohmic dissipation to ambipolar diffusion. We used $`n_{g13}`$ and $`T_{g3}`$ to denote $`n_g/(10^{13}\mathrm{cm}^3)`$ and $`T_g/(10^3\mathrm{K})`$, respectively. Ohmic dissipation dominates where $`B<B_0`$; ambipolar diffusion wins where $`B>B_0`$ (see Eqn. 2.1). Remarkably, $`B_0`$ is independent of the ionization fraction; so is the magnetic field profile. Solving Eqn. (10), we obtain

$$\frac{1}{3}B_y^3+B_0^2B_y=\mathrm{\Gamma }x,$$
(12)
assuming that the field vanishes at the origin. An integration constant, $`\mathrm{\Gamma }`$, is determined by the boundary condition $`B_y(x=L_{BH})=B_{max}`$, which leads to $`\mathrm{\Gamma }=B_{max}(B_0^2+B_{max}^2/3)/L_{BH}`$. The amplitude of the magnetic field, $`B_{max}`$, is determined from the remnant magnetic field strengths recorded in the matrix of chondrites (0.1–1 G) and in chondrules themselves (1–10 G) (Desch & Connolly, 2002, and references therein). From Eqn. (12), we see that $`B_0^2B_y^2/3`$ near the origin, and the effect of resistivity yields $`B_y\mathrm{\Gamma }x/B_0^2`$. Note the current density singularity is now removed; the maximum current density near the origin is given by $`J(x=0)=c\mathrm{\Gamma }/4\pi B_0^2`$, which is finite. Sufficiently far from the origin where $`B_0^2B_y^2/3`$ is satisfied, $`B_y(3\mathrm{\Gamma }x)^{1/3}`$ holds, and we retain the result in the non-resistive limit of BZ94. There, the current density declines steeply ($`Jx^{2/3}`$). We therefore set $`B_0^2=B_y^2/3`$ and compute the width of the current sheet to be

$`l_{cs}`$ $`=`$ $`{\displaystyle \frac{2\sqrt{3}B_0^3}{\mathrm{\Gamma }}}`$
$``$ $`1.5\times 10^5n_{g13}^{5/2}T_{g3}^{3/4}B_G\left({\displaystyle \frac{B_{max}}{3\mathrm{G}}}\right)^3\left({\displaystyle \frac{R}{3\mathrm{AU}}}\right)^{3/2}\mathrm{km}.`$
This expression is valid for $`B_{max}>\sqrt{3}B_0`$, which is satisfied under the conditions of interest.

### 2.4 Heating Mechanisms

Inside current sheets, two heating mechanisms act. First, electrons and ions are accelerated by the electric field present, and then collide with neutrals and heat up the entire gas, leading to familiar Ohmic dissipation. The corresponding heating rate per unit volume essentially follows from the magnetic field profile:

$$\epsilon _\mathrm{\Omega }=\frac{j^2}{\sigma _e}=\frac{\eta }{4\pi }\left(\frac{dB_y}{dx}\right)^2.$$
(14)
There is additional heating due to friction when ions, accelerated by the Lorentz force, diffuse through neutrals (see, e.g., § 2.1 of this paper; Padoan, Zweibel, & Nordlund, 2000). In this case of ambipolar diffusion heating,

$$\epsilon _{AD}=\rho _i\nu _{in}|𝒗_i𝒗_n|^2=\frac{\eta B_y^2}{\pi B_0^2}\left(\frac{dB_y}{dx}\right)^2$$
(15)
yields the volumetric heating rate. Note that their ratio is given simply by $`\epsilon _{AD}/\epsilon _\mathrm{\Omega }=B_y^2/4B_0^2`$. Within the range of parameters we are interested in, the ambipolar diffusion heating integrated through the current sheet along $`x`$ exceeds (by a factor of $``$10 in our fiducial model) the contribution from Ohmic dissipation.

The resulting heating profiles (log $`\epsilon `$ vs. log $`x`$) as a function of $`B_0`$, defined in Eqn. (11), are displayed in Figure 3. They are naturally symmetric around $`x=0`$. In terms of shape, these double peaks are similar to what has been observed in numerical simulations of Padoan, Zweibel, & Nordlund (2000), although their sheets are produced by supersonic shocks and density compressions. Ambipolar diffusion sharpening becomes more efficient at lower gas densities (see Eqn. 2.3), and the current densities become close to singular. This implies stronger and more concentrated heating far above and below the disk midplane. Alternatively, if the disk becomes thinner with time (Ruden, 1999), the heating rate also increases. For example, soon after thermonuclear reactions began in the Sun, a powerful solar wind is thought to have blown away the gas $``$0.1 to a few million years after the solar nebula formed. Accretion onto the central star or, in local regions, planetary accretion as gap formation occurs would have the same effect. A density drop caused by any of these mechanisms introduces a regime where chondrule formation is favored. As the overall gas density $`n_g`$ drops in time, $`B_0n_g`$, so the peak heating rate rises.

## 3 Applications

### 3.1 Energy Equations and Chondrule Formation

We investigate chondrule heating in current sheets by solving numerically the coupled equations of energy transfer for gas and dust particles, and compare the results with the experimental constraints on chondrule formation.

Dust particles are usually assumed to possess a power-law distribution of sizes with upper and lower size cutoffs (e.g. Chiang et al., 2001). For simplicity, we use just two types of dust particles: small, micron-sized dust grains with radius $`a_d=1\mu \mathrm{m}`$ and bigger, millimeter-sized chondrules with $`a_c=1\mathrm{mm}`$. Also, it is assumed that the same amount of mass is contained in each type. Near the dust melting temperature, small grains radiate in the near-infrared most efficiently. Hence their thermal evolution is of considerable importance to the overall development of current sheets.

To study how these current sheets affect the temperatures of dust particles that emit the bulk of the infrared radiation, we use the one-dimensional slab geometry considered in Hood & Horanyi (1993). Dust and gas travel through a current sheet with relative velocity of $`v_r`$, where $`v_r`$ is a free parameter of the model but assumed to be comparable to the Alfven speed, $`v_A=B/(4\pi \rho _n)^{1/2}`$, and smaller than 17 km s-1, the Keplerian orbital speed at $`R=3`$ AU. For our fiducial model, $`v_r=1`$ km s-1. We discretize the region for computation so that each zone has a thickness of 1 km and we use $`10^3`$ zones. Note that the lateral extent of our computational domain exceeds that of our fiducial current sheet ($`l_{cs}10^2`$ km). The gas and dust temperatures are uniform in each zone but in general differ from each other. The temperatures evolve in time and are assumed to vary only in one spatial direction ($`x`$).

Entering a current sheet, gas molecules gain energy from the two heating mechanisms mentioned in § 2.4 and lose it through collisions with dust grains and chondrules. The collisional heat transfer rate per unit surface area between gas and dust, denoted by $`q_d`$, is defined as $`q_d=(3/4\sqrt{\pi })(2k/m_n)^{3/2}\rho _g(T_gT_d)`$ and $`q_c`$ is defined similarly for chondrules (Hood & Horanyi, 1993). Then the gas temperature evolves in time according to

$$n_gk\frac{dT_g}{dt}=\epsilon \underset{i=d,c}{}4\pi a_i^2n_iq_i+\rho _g\kappa _{IR}\left(J_r4\sigma T_g^4\right),$$
(16)
where $`k`$ is the Boltzmann constant, $`n_i`$ is the number density of either dust grains ($`i=d`$) or chondrules ($`i=c`$), and $`\epsilon =\epsilon _\mathrm{\Omega }+\epsilon _{AD}`$. The opacity of the gas, $`\kappa _{IR}`$, is mainly due to dust associated with it, and is set to 1.0 cm2 g-1, appropriate for infrared wavelengths (Chiang et al., 2001). The intensity of radiation, $`J_r`$ (integrated over all wavelengths) at $`x=x_0`$ from plane-parallel, temperature-stratified atmospheres has the following form (Mihalas, 1978; Desch & Connolly, 2002):

$$J_r(x_0)=2\sigma _{\mathrm{}}^+\mathrm{}(\rho _g\kappa _{IR}T_g^4+\underset{i=d,c}{}n_i\pi a_i^2ϵ_iT_i^4)E_1(\beta |xx_0|)𝑑x$$
(17)
with $`\beta \rho _g\kappa _{IR}T_g^4+_{i=d,c}n_i\pi a_i^2ϵ_i`$. Here, $`ϵ_i`$ is the emissivity of a grain of type $`i`$ averaged over the Planck function at temperature $`T_i(x)`$, which we simply set at 1.0 (for comparison, Desch & Connolly (2002) uses values of 0.8 and 0.4 for 1$`\mathrm{mm}`$\- and 1$`\mu \mathrm{m}`$-sized grains, respectively), and $`E_1`$ is the exponential integral of the first type. Dust grains and chondrules, on the other hand, are heated both by colliding with hotter gas molecules and by absorbing radiation from nearby smaller, micron-sized, grains which reach high temperatures quickly. They lose energy by radiating an approximate blackbody spectrum at $`T_d`$ and $`T_c`$, with emission efficiency of $`ϵ_d`$ and $`ϵ_c`$ for dust and chondrules, respectively. Therefore,

$$m_iC_i\frac{dT_i}{dt}=\pi a_i^2\left[4q_i+ϵ_i\left(J_r4\sigma T_i^4\right)\right].$$
(18)
As before, the index $`i`$ can be either $`d`$ (for dust) or $`c`$ (for chondrule). The last term in Eqn. (18) represents the incident flux on a given particle radiated by dust grains in the other zones minus its own emission.

The experimental constraints on chondrule formation are (1) the grains are completely melted for a few minutes; (2) the cooling timescale of gas and dust is approximately an hour, consistent with the oft-quoted value of $`10^2`$$`10^3`$ K hr-1; and (3) small grains ($``$ 0.1 mm) are absent, possibly evaporated or destroyed completely after passing through a current sheet, leading to the observed narrow size distribution centered at $`a_d`$ 1 mm.

From the properties of current sheets described in § 2.32.4, we infer that the above constraints on chondrule formation will be satisfied only if the following conditions are met:

(1) In order to obtain high enough heating rate, the gas density has to be low ($`n_g10^{12}`$ cm-3), because the heating rate increases steeply as $`n_g`$ decreases.

(2) To have ”localized” heating events, spatially averaged chondrule/dust densities should be high ($`n_d>1`$ cm-3). Otherwise, the mean free path for IR photons will exceed a few hundred kilometers and the heating timescale will be much more than “a few minutes.”

(3) To ensure good thermal coupling between gas and dust, either the gas density has to be high ($`n_g10^{12}`$ cm-3) because $`q_d`$, $`q_cn_g`$ which cannot be satisfied because of the first condition, or dust grains need to be very small ($`a_d1\mu \mathrm{m}`$) so they approach $`T_g`$ quickly.

In order to satisfy conditions (1)–(3) simultaneously, the dust-to-gas mass ratio, $`\zeta `$, should be very large, about 50. Usually $`\zeta 0.01`$ is assumed for the minimum mass solar nebula, but after a few million years the gas density is expected to have decreased due to viscous diffusion (Ruden & Lin, 1986) as well as the various mechanisms mentioned in the last paragraph of § 2.4. Dust settles into a thin, dense layer on a timescale $`\tau _s(a/1\mu \mathrm{m})^1f_\mathrm{\Sigma }`$ Myr, independent of $`R`$ (Ruden, 1999, CG97), comparable to the lifetimes of protoplanetary disks. There will still be some dust grains, particularly submicron-sized ones, that remain well above the midplane mixed with the gas. But at late times the dust density outside the thin, dense dust layer is likely to be quite small, decreasing the dust-to-gas ratio there. Inside the dust layer, however, solids may dominate the local density (Cuzzi, Dobrovolskis, & Champney, 1993). Hence, $`\zeta 50`$ may be reached in the midplane, where most solids reside.

We solve Eqns. (16)–(18) simultaneously for $`T_g(x,t)`$, $`T_d(x,t)`$, and $`T_c(x,t)`$ using a fifth-order Runge-Kutta scheme appropriate for a stiff set of equations with adaptive stepsize control (Press et al., 1992). Both gas and dust particles begin at the same temperature, $`T_{init}=500`$ K. The center of a current sheet is initially placed at $`x=1000`$ km and moves to the left at $`v_r`$. To simplify the problem, we ignore the time evolution of $`\epsilon `$ and $`\chi _i`$ (see § 4). For our fiducial model, $`n_g=10^{12}`$ cm-3 and both $`n_d`$ and $`n_c`$ are set to yield desired values of the dust-to-gas ratio $`\zeta `$.

### 3.2 Results

Figure 4 shows the time evolution of gas, dust, and chondrule temperatures at a fixed position as a current sheet moves through, for our model runs with $`\zeta =`$ (5, 50). We also ran the code for $`\zeta =0.05`$ but do not plot the result because both $`T_d`$ and $`T_c`$ remain practically unchanged at $`T_{init}`$ for over $`10^3`$ s. We find that for $`\zeta =5`$ the gas temperature rises appreciably (above 1000 K) but, due to poor thermal coupling between gas and dust, neither $`T_d`$ nor $`T_c`$ changes significantly, reaching only $`700`$ K after 103 s. In the run with $`\zeta =`$ 50, however, due to the increase in optical depth, dust grains follow the gas temperature closely and the radiation field also steadily increases in strength. The grains melt after about 570 seconds and evaporate after about 630 seconds, at which point we stop the simulation, as it does not take into account the evaporation of particles. We changed the boundary condition between flat at $`T_{init}`$ and a floating boundary condition but found the result remained qualitatively the same. In an improved model that accounts for evaporation and a range of dust grain sizes, we expect the smallest grains to evaporate first, decreasing the optical depth, thus widening the heated region and lowering the net volumetric heating rate. The larger grains will therefore probably not evaporate. This is interesting since it could naturally explain the peaking of the size distribution of chondrules around 1 mm.

The fact that chondrules get melted only for high values of $`\zeta `$ is consistent with the FeO record in chondrules, which seems to indicate oxygen fugacities (relative abundance of oxygen to hydrogen) at formation that are higher than solar values by at least two or three orders of magnitude (Wood, 1967). This is most often interpreted to mean that chondrules were heated in regions with extremely high dust-to-gas ratios. The dust is taken to be oxides of Si, Mg, Ca, Al, etc. in chondritic proportions which bring in a large amount of oxygen. For this reason also, it is advantageous to evaporate the smaller grains first because then oxides will be incorporated into the gas and the partial pressures of condensable elements will increase locally. This helps to make liquids stable at high temperature, and also suppresses the evaporation of silicate liquids which would have caused heavy isotope enrichment of resulting chondrules; such enrichment is not observed (Galy, Young, Ash, & O’Nions, 2000). Ebel & Grossman (2000) gives a full treatment of what dust enrichments are necessary to make FeO-bearing silicates and also to make chondrule liquids stable at low pressure relative to solids, as well as a discussion of dust-enrichment scenarios with high oxygen fugacities. The idea that the midplane is a better place to consider chondrule formation is supported again in this context.

The value of the gas density used in Figure 4 is $`n_g=10^{12}`$ cm-3 corresponding to $`f_\mathrm{\Sigma }=10^1`$. Since we expect the gas density in the disk to decrease to this level after a time of order a few million years, this argues for chondrule formation in a late, evolved disk, and is consistent with the age difference between CAIs (refractory Ca-Al-rich inclusions) and chondrules, estimated from the decay of 26Al (with a half-life of 0.7 Myr) to be 2–3 Myr (Hutcheon et al., 1994; Russell, 1996; Swindle et al., 1996).

### 3.3 Current Sheets in the Disk Midplane

From the results in the previous section, we conclude that in a standard model of the early solar nebula in which dust and gas are well-mixed, the disk will have $`\zeta 0.01`$ everywhere and hence heating of dust grains by current sheets will not be effective. This conclusion no longer holds in the late stage of the disk evolution. As argued by Sano et al. (2000), the bulk of the disk becomes unstable to the MRI, therefore susceptible to current sheet formation, for low enough values of the surface density factor $`f_\mathrm{\Sigma }`$ or of the dust depletion factor $`f_d`$ because the recombination rate on grain surfaces decreases and thus $`\chi _i`$ rises except in a thin layer around the midplane. Current sheets capable of melting mm-sized dust grains form only when $`n_g<10^{12}`$ cm-3 and $`\zeta `$ is quite high, i.e. near the disk midplane. For these reasons, we henceforth consider a late stage of solar nebula evolution in which the gaseous disk has depleted by a factor of 10 ($`f_\mathrm{\Sigma }=10^1`$) and settling of dust grains onto the midplane has proceeded significantly ($`f_d10^2`$).

Let us first estimate the thickness of this dust-rich layer. We assume that $`\zeta `$ has increased from the well-mixed value of $`10^2`$ to $`50`$. Taking into account the reduction in gas density, the vertical scale height of dust grains is expected to decrease by a factor of 500 compared to the gas scale height of 0.1 AU at $`R=2`$ AU. Thus a thickness of $`3\times 10^4`$ km is estimated. This thin dust-rich layer around the midplane will be mostly magnetically dead because of the extremely high dust density there, and the resulting high recombination rate on grain surfaces (Sano et al., 2000). Yet, the thin surfaces of this dust layer will be sufficiently ionized by cosmic rays and stellar X-ray radiation to support the MRI and generate current sheets, considering that there are now few dust grains suspended in the upper atmosphere. Current sheets that form in these thin layers will ionize nearby regions within the dust layer, into optical depth of order unity. Dust grains are heated after passing through these current sheets, sometimes up to their melting temperature ($`T_{melt}1600`$ K) as shown in § 3.1, and then radiate strongly in the near-infrared, subsequently raising the temperature of dust grains in the surrounding environment. Collisions with hot dust grains may then ionize111The dominant ion species in protoplanetary disks are alkali metals (Na+, K+). Grains are dominant charge carriers at $`T_g10^3`$ K, but at higher $`T_g`$ desorption of charged ions from the surface of dust grains becomes effective (Umebayashi, 1983). Thus, our neglect of the effect of charged dust grains is partially justified, at least inside current sheets. the gas in an initially magnetically dead area above the critical value, $`\chi _{crit}`$.

However, the integrated optical depth of the protosolar nebula at $`R=2\mathrm{AU}`$ to the near-infrared radiation in the vertical direction is estimated to be $`\tau =\kappa _{IR}\mathrm{\Sigma }2\times 10^2`$ (D’Alessio, Calvet, & Hartmann, 2001; Chiang et al., 2001). Even if we assume that some of the smallest, submicron-sized grains have been removed from the nebula with the gas, the magnetically active region affected by current sheets would still be only a few percent of the dust-rich layer. To melt the majority of the dust grains that reside in the midplane and form chondrules, the entire layer must be in motion so that all dust particles are exposed in the magnetically active surface layer at some point.

When most of the solids have settled to the midplane, vertical shear between the gas and the differentially rotating, dense particle layer can drive turbulence within a 104–105 km thick boundary layer (Cuzzi, Dobrovolskis, & Champney, 1993; Cuzzi, Dobrovolskis, & Hogan, 1996). Although vertical velocity gradients vanish by symmetry very near the midplane, a turbulence model that self-consistently deals with generation, transport, and damping shows that the turbulence diffuses inward and persists even very near the midplane (Cuzzi, Dobrovolskis, & Champney, 1993, and references therein)222But see Sekiya (1998) and Youdin & Shu (2002) for claims that turbulent mixing becomes ineffective if the dust-to-gas surface mass ratio is significantly enhanced.. We take a typical velocity of $`U10^2`$ cm s-1 \[ see Eqn. (65) of the previous reference \] and the thickness of the dust layer, $`L_d3\times 10^4`$ km, and estimate the eddy turnover time to be $`L_d/U1`$ yr. Keeping in mind that the fraction of current sheets covering the dust surface is quite small ($`l_{cs}/L_{BH}3\times 10^5`$; see the next subsection), we expect an average dust grain to be exposed to current sheet heating about once every $`3\times 10^4`$ yrs. After the brief heating, now spherical, mm-sized chondrules will settle down to the midplane, collecting rims of dust in the process. Such rims are observed (Metzler & Bischoff, 1996, and references therein).

### 3.4 Near-Infrared Excesses

The higher than predicted flux at 2–8 $`\mu \mathrm{m}`$ observed in the spectra of some young protoplanetary disks has been known for a decade (Hartmann, Kenyon, & Calvet, 1993; Natta et al., 1999; Chiang et al., 2001). Using a disk model that included an accounting of grain size distributions and grain compositions, and numerical solution of the equations of radiative and hydrostatic equilibrium under the two-layer approximation of CG97, Chiang et al. (2001) rediscovered this in 4 of their 5 Herbig Ae and T Tauri stars, and found a relative dearth of emission shortward of 10$`\mu \mathrm{m}`$ arising partly from low absorption efficiencies of silicate grains. It is noteworthy that these sources also showed evidence for dust settling; parameters appropriate for gas and dust that are well-mixed in interstellar proportions did not fit the overall level of infrared excess at $`\lambda 100\mu \mathrm{m}`$ for 4 of them. They interpreted this to mean that dust in disk surface layers has settled vertically towards the midplane.

Dullemond, Dominik, & Natta (2001) proposed a disk model in which an inner rim located at the dust evaporation radius forms a vertical wall due to direct heating by the central star and, as a result, puffs up and intersects more starlight. Its reemitted flux produces a bump that peaks at 2–3$`\mu \mathrm{m}`$ in the SED. They note that this modification would be significant only for Herbig Ae stars. This mechanism works only if optical depth due to gas opacity in the inner hole is negligible, but in fact it is quite large, of order 103, when a standard surface mass density profile $`\mathrm{\Sigma }R^{3/2}`$ is assumed, although the surface density of the inner disk is not constrained by the data (Natta et al., 2001).

We suggest an alternative mechanism by noting that current sheets at the dust melting temperature would radiate at just the right wavelengths to fill in the deficient flux at 2–8 $`\mu `$m, perhaps explaining the excess near-IR fluxes. A simple estimate can be made as follows. We take the mean separation between current sheets (produced by the Brandenburg-Zweibel mechanism) to be $`L_{BH}=\lambda _{BH}/2`$, roughly the distance between two adjacent magnetic nulls produced by the MRI. The total emission from dust grains inside current sheets is given by (Chiang et al., 2001)

$`L_\lambda `$ $``$ $`{\displaystyle \frac{l_{cs}}{L_{BH}}}\left[\mathrm{\hspace{0.17em}8}\pi ^2\lambda {\displaystyle _{1\mathrm{AU}}^{2\mathrm{AU}}}2B_\lambda (T_{melt})R𝑑R\right]`$ (19)
$``$ $`10^{32}\left({\displaystyle \frac{\lambda _\mu ^4[\mathrm{exp}(9/\lambda _\mu )1]}{2\times 10^3}}\right)^1\mathrm{erg}\mathrm{s}^1,`$
for $`n_g=10^{12}`$ cm-3 at $`R=3\mathrm{AU}`$. Here, $`\lambda _\mu `$ is the wavelength of radiation in microns. Taking $`l_{cs}10^2\mathrm{km}`$ and $`L_{BH}0.03\mathrm{AU}`$, we obtain $`l_{cs}/L_{BH}3\times 10^5`$ for the fraction covered by current sheets in the region between 1 and 2 AU, the site of their formation in an evolved disk thinned by $`f_\mathrm{\Sigma }=10^1`$ (see Figure 1b). This crude estimate for the excess near-IR emission from current sheets may partly explain the observed excesses in the range of 1032–1033 erg s-1. Uncertainties in modeling of protoplanetary disks at present are admittedly large. Our mechanism involving current sheets could be taken as an example of producing near-IR excesses when magnetohydrodynamical effects are considered.

Ultimately, these ideas can be tested by searching for distinguishable observable signatures using techniques that have recently emerged. Interferometric imaging of LkH$`\alpha `$ in the H and K bands by Tuthill, Monnier, & Danchi (2001) has shown a central gap (or cavity) and a hot inner edge, observed at 3.4 AU from the central star when a distance of 160 pc to LkH$`\alpha `$ is assumed. Earlier theoretical predictions for its radius to fit the near-IR inflections in the SED ($`10R_{}`$; Hillenbrand et al., 1992) are too small by around an order of magnitude. If we can divide the disk into several annular regions and obtain photometric measurements of each annulus, we suggest that there would be a local maximum of near-IR radiation between 1 and 2 AU, for systems like our primordial solar nebula (although for other young stellar objects, the exact position will depend on stellar and disk parameters and should be recalculated). In standard protoplanetary disk models such as CG97, the mentioned region reaches temperatures of only $`100400`$ K and is expected to emit poorly in the near-IR. In the near future, a multi-band, adaptive optics coronagraph under development at the American Museum of Natural History (Oppenheimer, 2000) will be used to directly view protoplanetary disks and will present another possibility to test these theories.

## 4 Caveats

Our conclusions suffer from several caveats. All models of the current sheets that have been done to date have assumed a steady state flow (BZ94; Brandenburg & Zweibel, 1995; Zweibel & Brandenburg, 1997; Heitsch & Zweibel, 2003). However, the heating and ionization that occur in the current sheet are not accounted for self-consistently in these models. Both will increase the ion pressure in the center of the sheet, either saturating the singularity or even blowing the sheet apart. Also, we invoke turbulence induced by the MRI when calculating magnetic field fluctuations but ignore its effect on the evolution of current sheets, e.g. the converging motion of the neutrals in ambipolar diffusion. In our simplified chondrule formation model, we let gas and dust travel through a current sheet at $`v_r1`$ km s-1, but it is not clear how well this approximates the actual condition near a current sheet or if a current sheet would survive in that environment. To study these effects in detail a dynamical model will be necessary. This requires modeling the interplay between at least four physical processes: pressure imbalance, ionization evolution, temperature evolution, and steepening of the magnetic field. These four processes are inherently connected and feed back on one another. That means they must be solved simultaneously in a self-consistent fashion. Radiative transfer and a power-law size distribution of dust grains must be incorporated in a subsequent step to understand chondrule formation in these dynamical models. The large range of timescales and length scales involved presents another challenge for modeling of this complicated set of processes.

We are grateful to Axel Brandenburg, Eugene Chiang, Fabian Heitsch, and Ellen Zweibel for helpful discussions, as well as an anonymous reviewer of an NSF grant proposal for pointing out the importance of high oxygen fugacities derived for chondrules. M. K. R. J. was supported by an AMNH Graduate Student Fellowship. M.-M. M. L. acknowledges support by NSF Career grant AST99-85392, NSF grant AST03-07793, and NASA ATP grant NAG5-10103. D. S. E. acknowledges support from NASA Cosmochemistry grant NAG5-12855.