A&A 491, L17-L20 (2008)
DOI: 10.1051/0004-6361:200810753

LETTER TO THE EDITOR

X-ray emission from dense plasma in classical T Tauri stars: hydrodynamic modeling of the accretion shock

G. G. Sacco1,2 - C. Argiroffi3,2 - S. Orlando2,1 - A. Maggio2 - G. Peres3,2,1 - F. Reale3,2,1


1 - Consorzio COMETA, via S. Sofia 64, 95123 Catania, Italy
2 - INAF - Osservatorio Astronomico di Palermo, Piazza del Parlamento 1, 90134 Palermo, Italy
3 - DSFA - Università degli Studi di Palermo, Piazza del Parlamento 1, 90134 Palermo, Italy

Received 5 August 2008 / Accepted 1 October 2008

Abstract
Context. High spectral resolution X-ray observations of classical T Tauri stars (CTTSs) demonstrate the presence of plasma at temperature $T\sim 2{-}3\times 10^6$ K and density $n_{\rm e}\sim 10^{11}{-}10^{13}$ cm-3, which are unobserved in non-accreting stars. Stationary models suggest that this emission is due to shock-heated accreting material, but do not allow us to analyze the stability of the material and its position in the stellar atmosphere.
Aims. We investigate the dynamics and stability of shock-heated accreting material in classical T Tauri stars and the role of the stellar chromosphere in determining the position and thickness of the shocked region.
Methods. We perform one-dimensional hydrodynamic simulations of the impact of an accretion flow on the chromosphere of a CTTS, including the effects of gravity, radiative losses from optically thin plasma, thermal conduction and a well tested detailed model of the stellar chromosphere. We present the results of a simulation based on the parameters of the CTTS MP Mus.
Results. We find that the accretion shock generates an hot slab of material above the chromosphere with a maximum thickness of  $1.8 \times 10^9$ cm, density $n_{\rm e}\sim 10^{11}{-}10^{12}$ cm-3, temperature $T\sim 3\times 10^6$ K, and uniform pressure equal to the ram pressure of the accretion flow ($\sim$450 dyn cm-2). The base of the shocked region penetrates the chromosphere and remains at a position at which the ram pressure is equal to the thermal pressure. The system evolves with quasi-periodic instabilities of the material in the slab leading to cyclic disappearance and re-formation of the slab. For an accretion rate of $\sim$ $10^{-10}~M_{\odot}$ yr-1, the shocked region emits a time-averaged X-ray luminosity of $L_{\rm X}\approx 7\times 10^{29}$ erg s-1, which is comparable with the X-ray luminosity observed in CTTSs of identical mass. Furthermore, the X-ray spectrum synthesized from the simulation reproduces in detail all the main features of the O VIII and O VII lines of the star MP Mus.

Key words: X-rays: stars - stars: formation - accretion, accretion disks - hydrodynamics - shock waves - methods: numerical

1 Introduction

Pre-main sequence stars are strong X-ray emitters, with X-ray luminosities of up to 103 times the solar value ( $L_{\rm X}\sim 10^{27}$ erg s-1). In a similar way to main-sequence stars, the X-ray emission is probably generated by low density ( $n_{\rm e} \sim10^{10}$ cm-3) plasma enclosed in coronal loop structures and heated to temperatures of $T\sim10^6{-}10^7$ K (Feigelson & Montmerle 1999). During the last five years, high resolution ($R\sim 600$) X-ray observations of some classical T Tauri stars (CTTSs) (TW Hya, BP Tau, V4046 Sgr, MP Mus and RU Lupi) have indicated the presence of X-ray plasma at $T\sim 2{-}3\times 10^6$ K and denser than 1011 cm-3 (Günther et al. 2006; Kastner et al. 2002; Robrade & Schmitt 2007; Schmitt et al. 2005; Argiroffi et al. 2007), which suggests an origin that differs from being coronal.

Calvet & Gullbring (1998) and Lamzin (1998) proposed that X-ray emission from CTTSs could also be produced by the accreting material. In fact the material heated by the accretion shock at the base of the accretion column could reach a temperature of $T \sim 10^6$ K. By assuming that a stationary, strong accretion shock is present and by solving the conservation laws for the shock-heated region, Günther et al. (2007) demonstrated quantitatively that the accretion shock model describes some non coronal features of the X-ray observations of TW Hya. However, these models were based on several approximations and did not include a detailed model of the stellar atmosphere, which determines the shock position and could influence the profile of pressure and density and, therefore, the thickness of the shocked region. They assumed stationary conditions, but different studies have proven the existence of thermal instabilities in an accreting flow impacting the stellar surface (Langer et al. 1981; Chevalier & Imamura 1982; Koldoba et al. 2008). Considering that the physical structure of the shock-heated region and its temporal evolution determines the fraction of the accreted energy re-emitted in the X-ray band and its spectral behavior, further investigations of the issues discussed above are required to understand the precise nature of the X-ray emission from shock-heated plasma, to check its observability and, then, to derive the physical properties of the accretion process from high resolution X-ray observations.

In this paper, we address these issues with the aid of a time-dependent hydrodynamic numerical model describing the impact of an accretion stream on the chromosphere of a CTTS. As a first application, we have compared the X-ray spectra of the CTTS MP Mus (Argiroffi et al. 2007) with those synthesized from the results of a simulation tuned on this star.

2 The model

We assume that accretion occurs along a magnetic flux tube linking the circumstellar disk to the star. We consider a flux tube that is analogous to closed coronal loops observed on the Sun; therefore, as done for standard coronal loop models (e.g. Peres et al. 1982; Betta et al. 1997), we assume that the plasma moves and transports energy exclusively along the magnetic field lines. The argument that this hypothesis is appropriate is supported by the low value of the plasma parameter $\beta = P/(B^2/8\pi) \sim 10^{-2}$, where $P\sim 4\times 10^2$ dyn/cm2 is the expected post-shock zone pressure and $B \sim 10^3$ G is the typical magnetic field at the surface of a CTTS (Johns-Krull 2007). We note that both the plasma density and velocity are expected to vary across the section of an accretion stream (e.g. Romanova et al. 2003). However, in the case of $\beta \ll 1$, the stream can be considered to be a bundle of ``elementary'' streams, each characterized by different values of density and velocity. Our 1-D model describes one of these elementary streams. We limit our analysis to the impact of the accretion stream on the chromosphere and consider the portion of the flux tube close to the star. We assume that the magnetic field lines are perpendicular to the stellar surface and a plane-parallel geometry.

The impact of the accretion stream on the chromosphere is modeled by numerically solving the time-dependent fluid equations of mass, momentum, and energy conservation, taking into account the gravity stratification, the thermal conduction (including the effects of heat flux saturation) and the radiative losses from an optically thin plasma:

 \begin{displaymath}\frac{\partial \rho}{\partial t} + \frac{\partial \rho v}{\partial s} = 0, \end{displaymath} (1)
 \begin{displaymath}\frac{\partial \rho v}{\partial t} +\frac{\partial (P+\rho v^2)}{\partial s} = \rho g, \end{displaymath} (2)
 \begin{displaymath}\frac{\partial \rho E}{\partial t} +\frac{\partial (\rho E+P)... ...frac{\partial q}{\partial s} - n_{\rm e} n_{\rm H} \Lambda(T), \end{displaymath} (3)
\begin{displaymath}\epsilon =\frac{P}{\rho (\gamma-1)},\hspace{0.8cm} P=(1+\beta(\rho,T))\frac{\rho k_{\rm b} T}{m_{\rm A}} \end{displaymath}

where t is the time; s is the coordinate along the magnetic field lines; $\rho$ is the mass density; v is the plasma velocity; P is the thermal pressure; g(s) is the gravity; T is the plasma temperature; $E = \epsilon + v^2/2$ is the total gas energy per unit mass; $\epsilon$ is the internal energy per unit mass; q is the heat flux in the formulation of Spitzer (1962) corrected for the effect of flux saturation (Cowie & McKee 1977); $n_{\rm e}$ and $n_{\rm H}$ are the electron number density and hydrogen number density, respectively; $\Lambda(T)$is the radiative losses per unit emission measure; $\beta (\rho,T)$ is the fractional ionization ( $n_{\rm e}/n_{\rm H}$) derived from the modified Saha equation (see Brown 1973, for the solar chromosphere condition); $\gamma= 5/3$ is the ratio between the specific heats, $m_{\rm A} = 2.14 \times 10^{-24}$ g is the average atomic mass (assuming heavy elements abundances of 0.5 the solar values; Anders & Grevesse 1989), and $E_{\rm H}$ is a parametrized chromospheric heating function ( $E_{\rm H} = 0$ for $T > 8\times10^3$ K) defined, as in Peres et al. (1982), to maintain the unperturbed chromosphere in stable equilibrium. The radiative losses per unit emission measure are derived by the PINTofALE spectral code (Kashyap & Drake 2000) with the APED V1.3 atomic line database (Smith et al. 2001), and by assuming heavy elements abundances of 0.5 the solar values (Anders & Grevesse 1989), as obtained from X-ray observations of CTTSs by Telleschi et al. (2007). The equations are solved numerically using the FLASH code (Fryxell et al. 2000), an adaptive mesh refinement multiphysics code for astrophysical plasmas. The code has been extended with additional computational modules to simulate the radiative losses, thermal conduction, and the evolution of fractional ionization ( $n_{\rm e}/n_{\rm H}$).

The computational domain extends over a distance $D = 4.34 \times 10^9$ cm above the stellar surface. We allow for 5 levels of refinement in the adaptive mesh algorithm of FLASH ( PARAMESH; MacNeice et al. 2000), with resolution increasing twice at each refinement level: $\Delta s_{\rm max}=2.1 \times 10^6$ cm at the coarsest resolution, and $\Delta s_{\rm min}=1.3 \times 10^5$ cm at the finest level, which corresponds to a uniform mesh of $\sim$30 000 grid points. We analyzed the effect of spatial resolution on our results by considering two additional simulations that adopt an identical setup to the one discussed here, but with either 4 or 6 levels of refinement. We found that the adopted resolution is the most appropriate compromise between accuracy and computational cost and that the system evolution is described in detail accurately.

The simulation presented here covers a time interval of about 2000 s. We used the accretion parameters (velocity and density) derived by Argiroffi et al. (2007) to reproduce the soft X-ray emission of MP Mus. We calculated the gravity considering the star mass $M=1.2~M_{\odot}$and the star radius $R=1.3~R_{\odot}$ used by Argiroffi et al. (2007).

The external part of the initial configuration, extending from $2.75 \times 10^8$ to $4 \times 10^9$ cm, consists of an accretion stream constant in density ( $n_{\rm e}=10^{11}$ cm-3), temperature (T= 103 K) and velocity (v=450 km s-1). The inner part of the initial configuration consists of a static chromosphere. We reproduce the pressure gradient of a young stellar chromosphere by considering the temperature profile prescribed by the solar chromosphere models of Vernazza et al. (1973), scaled to reproduce a pressure of $\sim$ $7\times 10^4$ dyn cm-2at the base of the chromosphere. As boundary conditions, we consider fixed values both at the top ( $n_{\rm acc}=10^{11}$ cm-3, $T_{\rm acc}= 10^3$ K, $v_{\rm acc}=450$ km s-1) and at the base ( $n_{\rm chr} =10^{17}$ cm-3, $T_{\rm chr} =4.44 \times 10^3$ K, $v_{\rm chr} =0.0$ km s-1) of the computational domain. In principle these boundary conditions produce an accumulation of matter at the base of the chromosphere. We have estimated, however, that this effect is significant only for timescales longer by a factor 200 than those explored by our simulations. In addition, we confirmed that for s<108 cm the chromosphere remains virtually unperturbed (with variations in mass density below 1%) during the timescale considered. We note that we neglect the heating of the chromosphere (in particular, at the lower boundary) due to the X-ray emission originating from hot plasma. In the case of MP Mus (effective temperature $T_{\rm eff} \sim 5000$ K), this approximation is justified by the low ratio of the energy flux originating in the accretion flow to that of the photospheric emission, $(\rho_{\rm acc}v_{\rm acc}^3)/(4\sigma T_{\rm eff}^4) \approx 0.1$, where $\sigma$ is the Stefan-Boltzmann constant. The stability of the chromosphere was tested by dedicated simulations longer than 100 ks, some of which considering also strong transient heating.

   
3 Results


  \begin{figure} \par\includegraphics[width=8.5cm,clip]{0753fig1.ps} \end{figure} Figure 1: Evolution of plasma temperature  A), density  B), pressure  C), and velocity  D) distributions along the flux tube from the chromosphere to the unperturbed accretion stream, sampled every 60 s from 0 to 300 s ( left panels) and every 30 s from 330 to 420 s ( right panels). The figure shows the inner portion of the spatial domain, including the chromosphere and the hot slab.
Open with DEXTER

The impact of the accretion stream on the stellar chromosphere generates a transmitted (into the chromosphere) and a reverse (into the accretion column) shocks. The latter propagates through the accretion column producing an hot slab. As pointed out by Argiroffi et al. (2007), in the strong shock limit (Zel'Dovich & Raizer 1967), the expected temperature of the slab is $\sim$ $(3/16)(\mu m_{\rm H}/k_{\rm b})v_{\rm acc}^2\approx 3\times 10^6$ K and, since it is subject to the radiative cooling, the expected maximum thickness is $\sim$ $v_{\rm ps}\tau_{\rm rad}$, where $v_{\rm ps} = v_{\rm acc}/4 \approx 100$ km s-1is the post-shock plasma velocity in the slab, and $\tau_{\rm rad}$ is the cooling time for the shocked gas defined as

\begin{displaymath}\tau_{\rm rad}=\frac{1}{\gamma-1}\frac{P}{n_{\rm e}n_{\rm H}\Lambda(T)}\sim 6.7\times10^{3}\frac{T^{3/2}}{n_{\rm e}}, \end{displaymath} (4)

where, for temperatures characteristic of the slab ( $T \approx 10^5{-}10^7$ K), we have approximated the cooling function as $\Lambda(T) \approx 6.0\times 10^{-20}~T^{-1/2}$ erg s-1 cm3. For typical values of Tand $n_{\rm e}$ in the hot slab ( $T\approx 3\times 10^6$ K and $n_{\rm e} = 4~n_{\rm acc} \sim 4\times 10^{11}$ cm-3), $\tau_{\rm rad}\approx 90$ s and the thickness of the slab is expected to be $D_{\rm slab} \approx 10^9$ cm. We note that deviations from equilibrium of ionization may be present beyond the shock surface due to the rapid growth in temperature (e.g. Lamzin 1998). However, these deviations are expected within a distance $d_{\rm NEI} \la v_{\rm ps}\tau_{\rm eq}\sim 10^7$ cm from the shock, where $\tau_{\rm eq}$, the ionization equilibrium timescale, has been estimated for the characteristic density and temperature of the slab. Since $d_{\rm NEI} \ll D_{\rm slab}$, non-equilibrium ionization effects are irrelevant for the evolution of the hot slab.

The evolution of temperature, density, pressure, and velocity during the first part of our simulation is shown in Fig. 1. During the first 300 s of evolution, the reverse shock gradually heats up the accretion column to a temperature of $T\sim 3\times 10^6$ K and the inflow accumulates material creating a hot slab of density $n_{\rm e} = 3{-}7\times$1011 cm-3 and constant pressure equal to the accretion flow ram pressure $P_{\rm ram}= \rho_{\rm acc} v_{\rm acc}^2 = 450$ dyn cm-2 (see left panels of Fig. 1). The base of the shocked region initially penetrates the chromosphere and comes at rest at the point at which the ram pressure equals the thermal pressure. The ram pressure, therefore, determines the position of the slab inside the chromosphere, while the chromosphere below the shock region remains virtually unperturbed. Due to accumulation of material at the base of the shocked column, the radiative losses gradually increase there. At the end of this phase, the extension of the hot slab is $\approx$ $1.8 \times 10^9$ cm.

In the subsequent phase of evolution, the strong radiative cooling of the high density plasma that has accumulated at the base of the hot slab triggers thermal instabilities there (Field 1965; see right panels of Fig. 1): the plasma temperature and pressure both decrease by more than two orders of magnitude in few seconds, leading to the collapse of the upper layers of the accretion column. As a consequence, the residual hot slab falls down and the reverse shock moves downward to the chromosphere. The additional compression due to the collapse leads to the further increase of plasma density and radiative cooling at the base of the hot slab which gradually cools down within $\sim$100 s.

After the accretion column is completely cooled, a new hot slab is generated (see the dashed-dotted lines in the right panels of Fig. 1) and the system starts a quasi-periodic evolution with alternate phases of heating and cooling lasting $\sim$400 s each, linked to the up and down displacement of the reverse shock. Our simulation follows the evolution of the system for five periods and finds no significant differences between the results for each period. We therefore presume that our results can be extrapolated to describe events on longer time intervals. We note that the fluctuation period depends on the parameters of the inflow and in part those of the chromosphere. The assumption of steady flow is appropriate for the limited time explored. The whole issue should be revisited by simulations covering timescales comparable to those of typical variations due to stellar rotation, magnetic field evolution, etc.

We note that Langer et al. (1981) discovered periodic variations in the accretion shock position and, therefore, in the thickness of the hot slab in the context of compact objects. As in our case, the cycle consisted of time intervals in which the shock moved upward and the gas accumulated behind and phases during which the shock moved downward as the hot gas cooled radiatively. In the context of CTTSs, Koldoba et al. (2008) found oscillations in the accretion shock position with periods of $\sim$ 0.02-0.2 s, due to the same mechanism described here. Short periods are primarily due to the values of accretion flow parameters (velocity and density) and heavy elements abundances adopted by these authors. In addition, our model includes the effects of thermal conduction (important in the energy budget) and the stellar chromosphere (which determines the shock position).

We checked the observability of the X-ray emission produced by the shock-heated material by synthesizing the spectrum in the energy range [0.5-8.0] KeV from the simulation results, and by calculating the time-averaged luminosity over the time interval covered by the simulation[*]. Since our model is one-dimensional, we have to assume a value of flux tube cross section. Considering the velocity of the accretion stream v=450 km s-1, its density $n_{\rm e}=10^{11}$ cm-3, and an accretion rate $\sim$ $10^{-10}~M_{\odot}$ yr-1, the accretion stream cross section is $\sim$ $6.5\times 10^{20}$cm2. Using this cross section, we derive an X-ray luminosity varying from $\sim$ $8.4\times 10^{21}$ to  $\sim1.6 \times 10^{30}$ erg s-1 with a time-averaged value of $L_{\rm X} = 7.2 \times 10^{29}$ erg s-1. Since this value is comparable with the typical overall luminosity observed in young T Tauri stars of the same mass (Preibisch et al. 2005), our model has demonstrated that the shock-heated material can contribute to the X-ray emission of the CTTSs.

In Fig. 2, we compare the X-ray spectrum synthesized from the simulation (averaged over the time interval covered by the simulation) with the spectrum of the star MP Mus observed with the Reflection Grating Spectrometers (RGS) on board of the XMM-Newton satellite (Argiroffi et al. 2007). The synthesis of the X-ray spectrum accounts for the instrumental response of XMM-RGS and the interstellar absorption ( $N_{\rm H} =5\times 10^{20}$ cm-2) derived from the observations, and assumes that the emitting source is at the same distance of MP Mus (86 pc). We also assume that the time-averaged X-ray luminosity produced by the shock-heated slab in the 18-23 Å wavelength range is equal to the observed MP Mus luminosity in the same range. With the above assumptions, it turns out that the accretion rate is $\sim$ $7.7\times 10^{-11}~M_{\odot}$ yr-1 and the accretion stream cross section is $\sim$ $5\times 10^{20}$ cm2 (i.e. a surface filling factor of about 0.5%). These values are in agreement with those ($\sim$ $5\times 10^{-11}~M_{\odot}$ yr-1 and $\sim$ $3\times 10^{20}$ cm2) derived by Argiroffi et al. (2007) from the analysis of the observations.


  \begin{figure} \par\includegraphics[width=8.5cm,clip]{0753fig2.ps}\end{figure} Figure 2: Observed X-ray spectrum of the star MP Mus (gray line) with the synthetic spectrum derived from the simulation (black line).
Open with DEXTER

Figure 2 shows O VIII Ly$\alpha$ and the O VII triplet. The ratio of the O VIII Ly$\alpha$ to O VII resonance line is a tracer of plasma temperature and the ratio of the O VII forbidden to intercombination line is a tracer of plasma density. The good agreement between the observed and synthetic line profiles demonstrates that our model is capable of explaining the origin of the entire dense X-ray plasma observed in CTTSs. Furthermore, our hypothesis on heavy elements abundances is supported by the agreement with the ratio of lines to continuum intensity.

4 Conclusions

We have modeled the impact of an accretion flow on the chromosphere of a young T Tauri star. Our main results are:

A large set of simulations, exploring the domain of the physical parameters of the system, will be studied in a future paper.

Acknowledgements

We thank Jeremy Drake for useful discussions. This work was supported in part by the Italian Ministry of University and Research (MIUR) and by Istituto Nazionale di Astrofisica (INAF). The software used in this work was in part developed by the DOE-supported ASC/Alliance Center for Astrophysical Thermonuclear Flashes at the University of Chicago. This work makes use of results produced by the PI2S2 Project managed by the Consorzio COMETA, a project co-funded by the Italian Ministry of University and Research (MIUR) within the Piano Operativo Nazionale ``Ricerca Scientifica, Sviluppo Tecnologico, Alta Formazione'' (PON 2000-2006). More information is available at http://www.pi2s2.it and http://www.consorzio-cometa.it.

References

 
Copyright ESO 2008