X-ray emitting MHD accretion shocks in classical T Tauri stars

S. Orlando; G. G. Sacco; C. Argiroffi; F. Reale; G. Peres; A. Maggio

doi:10.1051/0004-6361/200913565

Free Access

Issue		A&A Volume 510, February 2010


Article Number		A71
Number of page(s)		12
Section		Interstellar and circumstellar matter
DOI		https://doi.org/10.1051/0004-6361/200913565
Published online		12 February 2010

A&A 510, A71 (2010)

Case for moderate to high plasma- $\beta$ values

S. Orlando¹ - G. G. Sacco^2,1 - C. Argiroffi^3,1 - F. Reale^3,1 - G. Peres^3,1 - A. Maggio¹

1 - INAF-Osservatorio Astronomico di Palermo, Piazza del Parlamento 1, 90134 Palermo, Italy
2 - Chester F. Carlson Center for Imaging Science, Rochester Inst. of Technology, 54 Lomb Memorial Dr., Rochester, NY, 14623, USA
3 - DSFA-Università degli Studi di Palermo, Piazza del Parlamento 1, 90134 Palermo, Italy

Received 28 October 2009 / Accepted 9 December 2009

Abstract
Context. Plasma accreting onto classical T Tauri stars (CTTS) is believed to impact the stellar surface at free-fall velocities, generating a shock. Current time-dependent models describing accretion shocks in CTTSs are one-dimensional, assuming that the plasma moves and transports energy only along magnetic field lines ( $\beta \ll 1$ ).
Aims. We investigate the stability and dynamics of accretion shocks in CTTSs, considering the case of $\beta \ga 1$ in the post-shock region. In these cases the 1D approximation is not valid and a multi-dimensional MHD approach is necessary.
Methods. We model an accretion stream propagating through the atmosphere of a CTTS and impacting onto its chromosphere by performing 2D axisymmetric MHD simulations. The model takes into account the stellar magnetic field, the gravity, the radiative cooling, and the thermal conduction (including the effects of heat flux saturation).
Results. The dynamics and stability of the accretion shock strongly depend on the plasma $\beta$ . In the case of shocks with $\beta > 10$ , violent outflows of shock-heated material (and possibly MHD waves) are generated at the base of the accretion column and intensely perturb the surrounding stellar atmosphere and the accretion column itself (therefore modifying the dynamics of the shock). In shocks with $\beta \approx 1$ , the post-shock region is efficiently confined by the magnetic field. The shock oscillations induced by cooling instability are strongly influenced by $\beta$ : for $\beta > 10$ , the oscillations may be rapidly dumped by the magnetic field, approaching a quasi-stationary state, or may be chaotic with no obvious periodicity due to perturbation of the stream induced by the post-shock plasma itself; for $\beta \approx 1$ the oscillations are quasi-periodic, although their amplitude is smaller and the frequency higher than those predicted by 1D models.

Key words: accretion, accretion disks - magnetohydrodynamics (MHD) - shock waves - stars: pre-main sequence - X-rays: stars - instabilities

1 Introduction

In the last few years, high resolution (R=600) X-ray observations of several young stars accreting material from their circumstellar disk (TW Hya, BP Tau, V4046 Sgr, MP Mus and RU Lupi) suggested the presence of X-ray emission at a temperature T=2-5 MK from plasma denser than $n_{\rm H} = 10^{11}$ cm^-3 (Kastner et al. 2002; Robrade & Schmitt 2007; Argiroffi et al. 2007; Günther et al. 2006; Argiroffi et al. 2009; Schmitt et al. 2005). The emitting plasma is too dense to be enclosed inside coronal loop structures (Testa et al. 2004). Several authors therefore suggested that this soft X-ray emission component could be produced by the material accreting onto the star surface, flowing along the magnetic field lines of the nearly dipolar stellar magnetosphere, and heated to temperatures of few MK by a shock at the base of the accretion column (Lamzin 1998; Calvet & Gullbring 1998).

The idea that an accreting flow which impacts onto a stellar surface leads to a shocked slab of material heated at millions degrees is not new. For instance, in the context of degenerate stars, several authors have studied the dynamics and energetic of accretion shocks and the effects of radiation on the formation, structure, and stability of the shocks (Imamura 1985; Chanmugam et al. 1985; Chevalier & Imamura 1982; Langer et al. 1981,1982). Hujeirat & Papaloizou (1998) have also investigated the role of the magnetic field in confining the post-shock accreting material and in determining the evolution of the accretion shock. These studies have shown that the role of thermal and dynamical instabilities is critical to our ability to model radiative shock waves.

The first attempt to interpret the evidence of soft X-ray emission from dense plasma in classical T Tauri stars (CTTSs) in terms of accretion shocks has quantitatively demonstrated that some non-coronal features of the X-ray observations of TW Hya (Günther et al. 2007) and MP Mus (Argiroffi et al. 2007) could be explained through a simplified steady-flow shock model. But, steady models are known to be an unsuitable approximation for radiative shocks, since they neglect the important effects of local thermal instabilities as well as global shock oscillations induced by radiative cooling (Sutherland & Dopita 1993; Dopita & Sutherland 1996; Sutherland et al. 2003b,a; Safier 1998; Mignone 2005).

In fact, the first time-dependent 1D models of radiative accretion shocks in CTTSs have shown quasi-periodic oscillations of the shock position induced by cooling (Sacco et al. 2008; Koldoba et al. 2008). In particular, Sacco et al. (2008) have developed a 1D hydrodynamic model including, among other important physical effects, a well-tested detailed description of the stellar chromosphere, and they also investigated the role of the chromosphere in determining the position and the thickness of the shocked region. For hydrodynamical simulations based on the parameters of MP Mus, Sacco et al. (2008) synthesized the high resolution X-ray spectrum, as it would be observed with the Reflection Grating Spectrometers (RGS) on board the XMM-Newton satellite, and found an excellent agreement with the observations.

Up to date time-dependent models of accretion shocks in CTTSs have been one dimensional, assuming that the plasma moves and transports energy only along magnetic field lines. This hypothesis is justified if the plasma $\beta \ll 1$ (where $\beta =\rm gas$ pressure/magnetic pressure) in the shock-heated material. The photospheric magnetic field magnitude of CTTSs is believed to be around 1 kG (e.g. Johns-Krull et al. 1999). Such a powerful stellar field is enough to efficiently confine accretion shocks with the particle density below 10¹³ cm^-3 and the temperature around 5 MK ( $\beta < 0.3$ ). However, recent polarimetric measurements indicate that in some cases the photospheric magnetic field strength could be less than 200 G (Valenti & Johns-Krull 2004), and the plasma $\beta$ in the slab may be around 1 or even larger. In these cases, the magnetic field configuration in the post-shock region may change, influencing the physical structure of the material emitting in X-rays and the stability of the accretion shock.

The low- $\beta$ approximation was challenged by recent findings of Drake et al. (2009). In the case of $\beta \ll 1$ , the radiative shock instability is expected to lead to a detectable periodic modulation of the X-ray emission from the shock-heated plasma, if the density and velocity of the accretion stream do not change over the time interval considered (in agreement with predictions of 1D models). But, the analysis of soft X-ray observations of the CTTS TW Hya (whose emission is believed to arise predominantly from accretion shocks) produced no evident periodic variations (Drake et al. 2009). These authors therefore concluded that 1D models may be too simple to describe the multi-dimensional structure of the shock, and that the magnetic field may play an important role through the generation and damping of MHD waves.

In this paper, we investigate the stability and dynamics of accretion shocks for cases where the low- $\beta$ approximation cannot be applied. We analyze the role of the stellar magnetic field in the dynamics and confinement of the slab of shock-heated material. To this end, we model an accretion stream which propagates through the magnetized atmosphere of a CTTS and impacts onto its chromosphere, using 2D axisymmetric MHD simulations and, therefore, an explicit description of the ambient magnetic field. We investigate cases of $\beta \ga 1$ in the post-shock region for which the deviations from 1D models are expected to be the largest. In Sect. 2 we describe the MHD model and the numerical setup; in Sect. 3 we describe the results; in Sect. 4 we discuss the implications of our results and finally we draw our conclusions in Sect. 5.

2 MHD modeling

The model describes an accretion stream impacting onto the surface of a CTTS. We assume that the accretion occurs along magnetic field lines that link the circumstellar disk to the star and consider only the portion of the stellar atmosphere close to the star.

The fluid is assumed to be fully ionized with a ratio of specific heats $\gamma = 5/3$ . The model takes into account the stellar magnetic field, the gravity, the radiative cooling, and the thermal conduction (including the effects of heat flux saturation). Since the magnetic Reynolds number $\gg$ 1 considering the typical velocity (10⁷ cm s^-1) and length scale (10⁹ cm) of the system, the flow is treated as an ideal MHD plasma. The impact of the accretion stream is modeled by numerically solving the time-dependent MHD equations (written in non-dimensional conservative form):

$\begin{displaymath}\frac{\partial \rho}{\partial t} + \nabla \cdot (\rho \vec{u}) = 0, \end{displaymath}$

(1)

$\begin{displaymath}\frac{\partial \rho \vec{u}}{\partial t} + \nabla \cdot (\rho \vec{u}\vec{u}-\vec{B}\vec{B}) + \nabla P_* = \rho \vec{g}, \end{displaymath}$

(2)

$\displaystyle \frac{\partial \rho E}{\partial t} +\nabla\cdot [\vec{u}(\rho E+P... ...c{u}\cdot \vec{g} -\nabla\cdot \vec{F}_{\rm c} -n_{\rm e} n_{\rm H} \Lambda(T),$

(3)

$\begin{displaymath}\frac{\partial \vec{B}}{\partial t} +\nabla \cdot(\vec{u}\vec{B}-\vec{B}\vec{u}) = 0, \end{displaymath}$

(4)

where

$\begin{displaymath}P_* = P + \frac{B^2}{2}~,~~~~~~~~~~~~~ E = \epsilon +\frac{1}{2} u^2+\frac{1}{2}\frac{B^2}{\rho}, \end{displaymath}$

are the total pressure, and the total gas energy (internal energy, $\epsilon$ , kinetic energy, and magnetic energy) respectively, tis the time, $\rho = \mu {\rm m}_{\rm H} n_{\rm H}$ is the mass density, $\mu = 1.28$ is the mean atomic mass (assuming metal abundances of 0.5 of the solar values; Anders & Grevesse 1989), $m_{\rm H}$ is the mass of the hydrogen atom, $n_{\rm H}$ is the hydrogen number density, $\vec{u}$ is the gas velocity, g is the gravity, T is the temperature, $\vec{B}$ is the magnetic field, $\vec{F}_{\rm c}$ is the conductive flux, and $\Lambda(T)$ represents the optically thin radiative losses per unit emission measure (see Fig. 1) derived with the PINTofALE spectral code (Kashyap & Drake 2000) with the APED V1.3 atomic line database (Smith et al. 2001), assuming the same metal abundances as before (as deduced from X-ray observations of CTTSs; Telleschi et al. 2007). We use the ideal gas law, $P=(\gamma-1) \rho \epsilon$ .

$\begin{figure} \par\includegraphics[width=8.2cm,clip]{13565fg1.ps} \end{figure}$	Figure 1: Radiative losses for an optically thin plasma from the APED V1.3 atomic line database (Smith et al. 2001), assuming the metal abundances of 0.5 of the solar values (Anders & Grevesse 1989).
Open with DEXTER

The thermal conductivity in an organized magnetic field is known to be highly anisotropic, and it can be highly reduced in the direction transverse to the field. The thermal flux is therefore locally split into two components, along and across the magnetic field lines, $\vec{F}_{\rm c} = F_{\parallel}~\vec{i}+F_{\perp}~\vec{j}$ , where (see, for instance, Orlando et al. 2008)

$\begin{displaymath}\begin{array}{l}\displaystyle F_{\parallel} = \left(\frac{1}{... ...rp}}+ \frac{1}{[q_{\rm sat}]_{\perp}}\right)^{-1}, \end{array}\end{displaymath}$

(5)

to allow for a smooth transition between the classical and saturated conduction regime. In Eq. (5), $[q_{\rm spi}]_{\parallel}$ and $[q_{\rm spi}]_{\perp}$ represent the classical conductive flux along and across the magnetic field lines (Spitzer 1962)

$\begin{displaymath}\begin{array}{l}\displaystyle [q_{\rm spi}]_{\parallel} = -\k... ...\frac{n^2_{\rm H}}{T^{1/2}B^2}~ [\nabla T]_{\perp}, \end{array}\end{displaymath}$

(6)

where $[\nabla T]_{\parallel}$ and $[\nabla T]_{\perp}$ are the thermal gradients along and across the magnetic field, and $\kappa_{\parallel}$ and $\kappa_{\perp}$ (in units of erg s^-1 K^-1 cm^-1) are the thermal conduction coefficients along and across the magnetic field lines, respectively. The saturated flux along and across the magnetic field lines, $[q_{\rm sat}]_{\parallel}$ and $[q_{\rm sat}]_{\perp}$ , are (Cowie & McKee 1977)

$\begin{displaymath}\begin{array}{l}\displaystyle [q_{\rm sat}]_{\parallel} = -\m... ...[\nabla T]_{\perp}\right)~ 5\phi \rho c_{\rm s}^3, \end{array}\end{displaymath}$

(7)

where $c_{\rm s}$ is the isothermal sound speed, and $\phi$ is a number of the order of unity; we set $\phi = 1$ according to the values suggested for stellar coronae (Borkowski et al. 1989; Fadeyev et al. 2002; Giuliani 1984, and references therein).

$\begin{figure} \par\includegraphics[width=8.cm,clip]{13565fg2.ps} \end{figure}$	Figure 2: Initial geometry of the system in cylindrical coordinates. The stellar surface lies on the r axis, and the unperturbed stellar magnetic field is uniform and oriented along the z axis (vertical lines). The accretion stream propagates downwards through the stellar corona with a velocity $u_{\rm str}$ . The z-axis is the symmetry axis of the problem.
Open with DEXTER

We solve the MHD equations using cylindrical coordinates in the plane (r,z), assuming axisymmetry and the stellar surface lying on the r axis (see Fig. 2). The initial unperturbed stellar atmosphere is assumed to be magneto-static and to consist of a hot (maximum temperature $\approx$ 10⁶ K) and tenuous ( $n_{\rm H} \approx 2\times 10^8$ cm^-3) corona linked through a steep transition region to an isothermal chromosphere that is uniformly at temperature 10⁴ K and is $8.5\times 10^8$ cm thick (see Fig. 3). The choice of an isothermal chromosphere is adopted for ease of modeling, and different choices of more detailed chromospheric models have been shown not to lead to significantly different results (Sacco et al. 2010, in preparation). The unperturbed stellar magnetic field is uniform, oriented along the z axis and perpendicular to the stellar surface. The gravity is calculated assuming the star mass $M=1.2 ~M_{\odot}$ and the star radius $R=1.3~R_{\odot}$ , which is appropriate for the CTTS MP Mus (see Argiroffi et al. 2007). Different choices of stellar parameters should not lead to different results.

$\begin{figure} \par\includegraphics[width=8.5cm,clip]{13565fg3.ps} \end{figure}$	Figure 3: Initial hydrogen number density (dashed line) and temperature (dotted line) as a function of height above the stellar surface z for the unperturbed stellar atmosphere.
Open with DEXTER

The accretion stream is assumed to be constant and propagates along the z axis. Initially the stream is in pressure equilibrium with the stellar corona and has a circular cross-section with a radius $r_{\rm str} = 5\times 10^{9}$ cm; its radial density distribution is given by

$\begin{displaymath}n_{\rm str}(r) = n_{\rm cor}+\frac{n_{\rm str0}-n_{\rm cor}} {\cosh\left[\sigma\left(r/r_{\rm str}\right)^{\sigma}\right]}, \end{displaymath}$

(8)

where $n_{\rm str0}$ and $n_{\rm cor}$ are the hydrogen number density at the stream center and in the corona, and $\sigma=20$ is a steepness parameter. The above distribution describes a thin shear layer ( $\sim$ $0.1~r_{\rm str}$ ) around the stream that smooths the stream density to the value of the surrounding medium and maintains numerical stability, allowing for a smooth transition of the Alfvén speed at the interface between the stream and the stellar atmosphere. The stream temperature is determined by the pressure balance across the stream boundary.

The velocity of the stream along the z axis also has a radial profile:

$\begin{displaymath}u_{\rm str}(r) = \frac{u_{\rm str0}} {\cosh\left[\sigma\left(r/r_{\rm str}\right)^{\sigma}\right]}, \end{displaymath}$

(9)

where $u_{\rm str0}$ is the free fall velocity at the center of the stream. The smooth shear layer of the velocity avoids the growth of random perturbations with wavelengths on the order of the grid size (e.g. Bodo et al. 1994).

The simulations presented here describe the evolution of the system for a time interval of about 3000 s. We used the accretion parameters (velocity and density) that match the soft X-ray emission of MP Mus (Argiroffi et al. 2007), namely $n_{\rm str0} = 10^{11}$ cm^-3and $u_{\rm str0} = -500$ km s^-1 at a height $z =2.4\times 10^{10}$ cm above the stellar surface. The initial field strengths of $\vert\vec{B}_{\rm0}\vert = 1,~10,~50$ G in the unperturbed stellar atmosphere (see Table 1) maintain $\beta \ga 1$ in the post-shock region, where $\beta=P/(B^2/8\pi)$ is the ratio of thermal to magnetic pressure. For instance, $\vert\vec{B}_{\rm0}\vert = 10$ G corresponds to $\beta$ ranging between $\approx$ 10⁴, at the base of the chromosphere, and $\approx$ 10^-2 up in corona. We also performed a 1D hydrodynamic simulation to provide a baseline for the 2D calculations. This simulation corresponds to the strong magnetic field limit ( $\beta \ll 1$ ), in which the plasma moves and transports energy exclusively along the magnetic field lines (see Sacco et al. 2008).

Table 1: Relevant parameters of the simulations.

The calculations described in this paper were performed using PLUTO (Mignone et al. 2007), a modular, Godunov-type code for astrophysical plasmas. The code provides a multiphysics, multialgorithm modular environment particularly oriented toward the treatment of astrophysical flows in the presence of discontinuities as in the case treated here. The code was designed to make efficient use of massive parallel computers using the message-passing interface (MPI) library for interprocessor communications. The MHD equations are solved using the MHD module available in PLUTO which is based on the Harten-Lax-van Leer Discontinuities (HLLD) approximate Riemann solver (Miyoshi & Kusano 2005). Miyoshi & Kusano (2005) have shown that the HLLD algorithm can exactly solve isolated discontinuities formed in the MHD system; consequently the adopted scheme is particularly appropriate to describe accretion shocks. The evolution of the magnetic field is carried out using the constrained transport method of Balsara & Spicer (1999) that maintains the solenoidal condition at machine accuracy. PLUTO includes optically thin radiative losses in a fractional step formalism (Mignone et al. 2007), which preserves the $2{\rm nd}$ time accuracy, as the advection and source steps are at least of the $2{\rm nd}$ order accurate; the radiative losses $\Lambda$ values are computed at the temperature of interest using a table lookup/interpolation method. The thermal conduction is solved with an explicit scheme that adds the parabolic contributions to the upwind fluxes computed at cell interfaces (Mignone et al. 2007). Such a scheme is subject to a rather restrictive stability condition (i.e. $\Delta t < (\Delta x)^2/(2\eta)$ , where $\eta$ is the maximum diffusion coefficient), as the thermal conduction timescale is generally shorter than the dynamical one (e.g. Orlando et al. 2008; Hujeirat & Camenzind 2000; Orlando et al. 2005; Hujeirat 2005).

The symmetry of the problem allows us to solve the MHD equations in half of the spatial domain with the stream axis coincident with the z axis. The 2D cylindrical (r,z) mesh extends between 0 and $1.2 \times 10^{10}$ cm in the r direction and between 0 and $2.4 \times 10^{10}$ cm in the z direction and consists of a uniform grid with $512\times 1024$ grid points. Additional runs were done with setups identical to those used for runs By-01, By-10, and By-50, but with a higher resolution ( $1024\times 2048$ grid points; see Table 1) to evaluate the effect of the spatial resolution (see Sect. 3.4); the simulations with higher resolution cover a time interval of about 1000 s.

We use fixed boundary conditions at the lower (z = 0) boundary, imposing zero material and heat flux across the boundary. With this condition, matter may progressively accumulate at the base of the chromosphere; we have estimated that this effect may become significant on timescales a factor of 100 longer than than those explored by our simulations. Axisymmetric boundary conditions are imposed at r = 0 (i.e. along the symmetry axis of the problem), a constant inflow at the upper boundary ( $z =2.4\times 10^{10}$ cm) for $r \leq r_{\rm str}$ , and free outflow elsewhere.

3 Results

3.1 One-dimensional reference model

The results of our 1D reference hydrodynamic simulation HD-1D are analogous to those discussed by Sacco et al. (2008): the impact of the accretion stream onto the stellar chromosphere generates a reverse shock, which propagates through the accretion column and produces a hot slab. According to Sacco et al., the expected temperature $T_{\rm\rm slab}$ and the maximum thickness $D_{\rm\rm slab}$ of the slab are in the strong shock limit (Zel'Dovich & Raizer 1967):

$\begin{displaymath}T_{\rm slab} \approx \frac{3}{32}\frac{\mu {\rm m}_{\rm H}}{k_{\rm B}}u_{\rm str0}^2\approx 5~\mbox{MK}, \end{displaymath}$

(10)

$\begin{displaymath}D_{\rm slab} \approx 4.2\times 10^2 \frac{u_{\rm str0}}{n_{\rm str0}} T_{\rm slab}^{3/2}\approx 2.5\times 10^9~\mbox{cm}, \end{displaymath}$

(11)

where $k_{\rm B}$ is the Boltzmann constant.

Figure 4 shows the time-space plot of the temperature evolution for run HD-1D. The base of the hot slab penetrates the chromosphere (the dashed line in Fig. 4 marks the initial position of the transition region between the chromosphere and the corona) down to the position at which the ram pressure, $P_{\rm ram} = \rho_{\rm str0} u_{\rm str0}^2$ , of the post-shock plasma equals the thermal pressure of the chromosphere (Sacco et al. 2008). As is evident from the figure, the shock front is not steady, and the amplitude and period of the pulses rapidly reach stationary values as the initial transient disappears. The shock position oscillates with a period of $\approx$ 600 s due to intense radiative cooling at the base of the slab, which robs the post-shock plasma of pressure support, causing the material above the cooled layer to collapse back before the slab expands again (see Sacco et al. 2008, for a detailed description of the system evolution). The maximum thickness of the slab is $D_{\rm slab}\approx 2.5\times 10^9$ cm, in agreement with the prediction (see Eq. (11)). The post-shock plasma reaches the temperature $T_{\rm slab} \approx 5{-}7$ MK during the expansion phase, and $T_{\rm slab} \approx 2{-}4$ MK during the cooling phase.

As already mentioned, the evolution of the accretion shock described by the 1D reference model is appropriate if the magnetic field lines can be considered to be rigid to any force exerted by the accreting plasma flowing along the lines ( $\beta \ll 1$ ). In case of plasma $\beta$ values around 1 or even higher, the hot slab is expected to be only partially confined by the magnetic field, and flow can occur sideways because of the strong pressure of the post-shock plasma and may ultimately perturb the dynamics of the shock itself. Moreover, 2D or 3D structures leading to fast cooling zones are expected to develop in the post-shock plasma if the latter is characterized by $\beta \ga 1$ . In this case, Sutherland et al. (2003a) showed that the evolution of 1D and 2D radiative shocks may be significantly different.

$\begin{figure} \par\includegraphics[width=8.5cm,clip]{13565fg4.eps} \end{figure}$	Figure 4: Time-space plot of the temperature evolution for run HD-1D. The spatial extent of the shock lies in the vertical direction at any instant in time. The dashed line marks the initial position of the transition region between the chromosphere and the corona.
Open with DEXTER

3.2 Two-dimensional radiative MHD shock model

The cases of $\beta \ga 1$ are described by our 2D model. The main differences of the evolution from the 1D model are expected to occur at the border of the stream where the shock-heated plasma may not be efficiently confined by the magnetic field; so we focus our attention there.

Figure 5 shows the spatial distribution of temperature (on the left) and plasma $\beta$ (on the right) for runs By-01, By-10, and By-50 at time t=530 s (at early stage of evolution). Movies showing the complete evolution of 2D spatial distributions of mass density (on the left) and temperature (on the right) in log scale for runs By-01, By-10, and By-50 are provided as on-line material. The accretion flow follows the magnetic field lines and impacts onto the chromosphere, forming a hot slab at the base of the stream with temperatures $\approx$ 5 MK and $\beta \ga 1$ . In all cases, the slab is rooted in the chromosphere, where the thermal pressure equals the ram pressure, and part ( $\approx$ 1/3) of the shock column is buried under a column of optically thick material and may suffer significant absorption. In runs By-01 and By-10, the dense hot plasma behind the shock front causes a pressure-driven flow parallel to the stellar surface, expelling accreted material sideways (see upper and middle panels in Fig. 5; see also on-line movies for runs By-01 and By-10). This outflow strongly perturbs the shock dynamics and is absent in the 1D shock model. As expected, this feature determines the main differences between these 2D and the 1D simulations.

$\begin{figure} \par\includegraphics[width=16.5cm,clip]{13565fg5.ps} \end{figure}$	Figure 5: Temperature ( left panels) and plasma $\beta$ ( right panels) distributions in the (r,z) plane in log scale in the simulations By-01, By-10, and By-50 at a time t=530 s. The initial position of the transition region between the chromosphere and the corona is at z=0. The magnetic field is initially oriented along the z axis; the white lines mark magnetic field lines.
Open with DEXTER

In run By-01, the magnetic field is too weak to confine the post-shock plasma ( $\beta \approx 10^4$ at the border of the slab), and a conspicuous amount of material continuously escapes from the border of the accretion column at its base, where the flow impacts onto the stellar surface. The maximum escape velocity is comparable to the free-fall velocity $u_{\rm str0}$ , and this outflow acts as an additional cooling mechanism. The resulting outflow advects and stretches the magnetic field lines (see upper panels in Fig. 5), taking the material away from the accretion column and strongly perturbing the stellar atmosphere even at several stream radii. As a result of the outflow, a large component of $\vec{B}$ perpendicular to the stream velocity ( $B_{\rm r} \approx 0.3 \vert\vec{B}_{\rm0}\vert$ ) appears in the post-shock region (see upper panels in Fig. 5). As discussed by Toth & Draine (1993), the presence of even a small magnetic field perpendicular to the flow can stabilize the overstable oscillations. In fact, at variance with the force due to gas pressure, the Lorentz force is not affected by cooling processes (see also Hujeirat & Papaloizou 1998), and this mechanism may contribute to stabilize the shock oscillations (see the on-line movie).

$\begin{figure} \par\includegraphics[width=15cm,clip]{13565fg6.ps} \end{figure}$	Figure 6: Temperature ( left panel) and mass density ( right panel) distributions in the (r,z) plane in log scale in the simulation By-10, at a time t=1800 s. The white lines mark magnetic field lines.
Open with DEXTER

In the intermediate run By-10, the magnetic field is trapped at the head of the escaped material, leading to a continuous increase of the magnetic pressure and field tension there. As a result, the escaped material is kept close to the accretion column by the magnetic field (at variance with the By-01 case) and may eventually perturb the stream itself. In fact, the escaped material accumulates around the accretion column, forming a growing sheath of turbulent material gradually enveloping the stream. The magnetic pressure and field tension increase at the interface between the ejected material and the surrounding medium, pushing on the material and forcing it to plunge into the stream after $\approx$ 1.3 ks. Figure 6 shows the temperature and mass density distributions at a time t=1.8 ks, when the expelled material has already entered into the stream and has deeply perturbed the accretion column (see also the on-line movie to follow the complete evolution). As a result of the stream perturbation, the hot slab may temporarily disappear altogether, as shown in Fig. 6. In this phase, the region of impact of the stream onto the chromosphere is characterized by a rather complex structure with knots and filaments of material and may involve possible mixing of plasma of the surrounding corona with accretion material.

In the By-50 case, the post-shock plasma is confined efficiently by the magnetic field, and no outflow of accreted material forms (see lower panels in Fig. 5). In particular, the plasma $\beta$ is $\approx$ 10^-4 and the magnetic pressure is $P_{\rm B} \approx 100 P_{\rm ram}$ at the stream border, where $P_{\rm ram} = \rho u^2$ is the ram pressure. The 2D shock therefore evolves similarly to the 1D overstable shock simulation, with alternating phases of expansion and collapse of the post-shock region (see the on-line movie). The maximum thickness of the slab is $D_{\rm slab} \approx 1.4\times 10^9$ cm, i.e. less than in the 1D case by a factor of $\approx$ 1.8. This result is analogous to that described by Sutherland et al. (2003a) for the case of 2D hydrodynamic radiative shocks and is explained there as due to the formation of denser knots of more rapidly cooling gas. Note that, at variance with our MHD simulations, the hydrodynamic model of Sutherland et al. (2003a) does not include the thermal conduction. In our simulations, the thermal conduction acts as an additional cooling mechanism of the hot slab, draining energy from the shock-heated plasma to the chromosphere, and partially contrasts the radiative cooling. Depending on the temperature of the shock and on the plasma $\beta$ , the thermal conduction may therefore lead to significant differences between our results and those of Sutherland et al. (2003a).

Figure 7 shows snapshots of the evolution of temperature distribution (in linear scale) in run By-50. The post-shock region gets bigger (t=1.05 ks), reaches the maximum extension (t= 1.15 ks), and collapses (t= 1.25 ks); the cycle then repeats and the slab begins to get larger again (t= 1.35 ks). Although the magnetic field is powerful enough to confine the post-shock plasma, the plasma $\beta$ is slightly larger than 1 inside the slab, and the 1D approximation is not valid there. As in the 2D hydrodynamic simulations (Sutherland et al. 2003a) therefore, complex 2D cooling structures, including knots and filaments of dense material, form there due to the thermal instability of the post-shock plasma (see t=1.15 ks in Fig. 7). As discussed by Sutherland et al. (2003a), these 2D complex structures lead to the zones cooling more efficiently than those in 1D models (by virtue of the increased cold-hot gas boundary) and consequently the amplitude of the oscillations is expected to be reduced. It is worth noting that the similarity of our results with those of the hydrodynamic models is due to the fact that, because $\beta \ga 1$ in the slab in run By-50, 2D cooling structures can form. In the case of $\beta \ll 1$ everywhere in the slab, we expect that the stream can be considered as a bundle of independent fibrils (each of them describable in terms of 1D models), and the 2D MHD simulations would produce the same results as 1D models.

$\begin{figure} \par\includegraphics[width=8.5cm,clip]{13565fg7.ps} \end{figure}$	Figure 7: Temperature distribution in the (r,z) plane in linear scale in the simulation By-50 at the labeled times. The white lines mark magnetic field lines.
Open with DEXTER

We also studied the global time evolution of 2D shocks by deriving time-space plots of the temperature evolution analogous to that derived for the 1D reference simulation (to compare directly our 1D and 2D results). From the 2D spatial distributions of temperature and mass density, we first derive profiles of temperature along the z-axis by averaging the emission-measure-weighted temperature along the r-axis for $r< 4\times 10^9$ cm (i.e. 80% of the stream radius); then we derive the time-space plots of the average temperature evolution from these profiles.

$\begin{figure} \par\includegraphics[width=8cm,clip]{13565fg8.eps} \end{figure}$	Figure 8: As in Fig. 4 for the 2D simulations By-01 ( top), By-10 ( middle) and By-50 ( bottom). The blue curve in the bottom panel marks the shock position derived in the 1D reference model, HD-1D.
Open with DEXTER

Figure 8 shows the results for runs By-01, By-10, and By-50, and can be compared with Fig. 4. In all the cases, the amplitude of the shock oscillations is smaller than the one observed in the 1D reference model. In runs By-01 and By-10, the evolution of the 2D shock is markedly different from that of the 1D shock (compare Fig. 4 with top and middle panels in Fig. 8): in run By-01, the oscillations of the shock are stabilized after 2.5 ks and the solution approaches a quasi-stationary state; in run By-10, the oscillations appear chaotic without an evident periodicity (at least in the time lapse explored here). As discussed before, the stabilization of the shock oscillations in the former case is due to a magnetic field component perpendicular to the flow (Toth & Draine 1993), whereas the chaotic oscillations in the latter case are due to the continuous stream perturbation by the material ejected sideways at the stream base (see Fig. 6). In the run By-50, the shock evolves somewhat similarly to the 1D overstable shock simulation, showing several complete cycles with expansion and collapse of the post-shock region (see lower panel in Fig. 8). However, as already discussed before, some important differences arise: by comparing By-50 with HD-1D, the oscillations are reduced in amplitude by a factor of $\sim$ 1.8 and occur at higher frequency (period $P_{\rm sh} \approx 300$ s).

3.3 Effect of stream parameters

The details of the shock evolution described in this paper depend on the model parameters adopted. In particular, the temperature of the shock-heated plasma is determined by the free-fall velocity with which the plasma impacts onto the star (Eq. (10)); the stand-off height of the hot slab generated by the impact depends on the velocity and density of the stream (see Eq. (11)). Assuming a typical value for the free-fall velocity of 400-500 km s^-1 (leading to temperatures of $\approx$ 3-5 MK, as deduced from observations), the maximum thickness of the slab is determined only by the density of the stream: the heavier the stream, the thinner is the expected slab. Also the sinking of the stream in the chromosphere depends on the stream density and velocity at impact. In fact, we find that the slab penetrates the chromosphere to the position at which the ram pressure $\rho_{\rm str0} u_{\rm str0}^2$ of the post-shock plasma equals the thermal pressure of the chromosphere (see also Sacco et al. 2008). As a result, heavier streams are expected to sink more deeply into the chromosphere.

In the case of $\beta > 1$ , the shock evolution is expected to be influenced also by the stream radius $r_{\rm str}$ . In fact, as shown by our simulations, the complex plasma dynamics close to the border of the stream (e.g. the generation of outflows of accreted plasma; see runs By-01 and By-10) may affect the shock evolution in the inner portion of the slab. The shock evolution is expected to be modified in the whole slab if the oscillation period of the shock, $P_{\rm sh}$ , is larger than the dynamical response time of the post-shock region, which can be approximated as the sound-crossing time of half slab: $\tau_{\rm dyn} \sim r_{\rm str}/c_{\rm s}$ , where $c_{\rm s}$ is the isothermal sound speed. In the cases discussed in this paper, the oscillation period is $P_{\rm sh} \approx 300$ s, and the sound-crossing time of the slab is $\tau_{\rm dyn}\approx 200$ s; as shown by our simulations, the shock evolution is modified in the whole slab. Assuming an accretion stream with twice the radius considered here, we find $\tau_{\rm dyn}\approx 400$ s, and the region at the center of the slab should not be affected by the plasma dynamics at the stream border. In this case, we expect quasi-periodic shock oscillations (analogous to those described by run By-50) at the center of the stream and a strongly perturbed shock at the stream border.

3.4 Effect of spatial resolution

In problems which involve radiative cooling, the spatial resolution and the numerical diffusion play an important role to determine the accuracy with which the dynamics of the system is described. In particular, in the case of radiative shocks, we expect that the details of the plasma radiative cooling depend on the numerical resolution: a higher resolution may lead to different peak density and hence influence the cooling efficiency of the gas, and therefore also the amplitude and frequency of shock oscillations. In the simulations presented here, the thermal conduction partially contrasts the radiative cooling and consequently alleviates the problem of numerical resolution (see, for instance, Orlando et al. 2008).

To check if our adopted resolution is sufficient to capture the basic shock evolution over the time interval considered, we repeated simulations By-01, By-10, and By-50, but with twice the spatial resolution (runs By-01-HR, By-10-HR, and By-50-HR). The efficiency of radiative cooling is expected to be the largest in the model with $\vert\vec{B}_{\rm0}\vert=50$ G, because there is no loss of material through sideway outflows. In fact, quasi-periodic oscillations of the accretion shock occur (see Fig. 8) and make this case adequate for a comparison of different spatial resolutions. Since this case is one of the most demanding for resolution, it can be considered a worst case comparison of convergence.

$\begin{figure} \par\includegraphics[width=8.5cm,clip]{13565fg9.ps} \end{figure}$	Figure 9: Evolution of the averaged shock-front location for runs By-50 (solid blue line) and By-50-HR (dashed red line). The dotted line marks the initial position of the transition region between the chromosphere and the corona.
Open with DEXTER

Figure 9 compares the evolution of the averaged shock-front location for runs By-50 and By-50-HR. The main differences between the simulations appear in the first bounce, which starts earlier and is larger by a factor of $\approx$ 2 in By-50-HR than in By-50. The first expansion of the hot slab occurs after that the stream has penetrated the chromosphere to the position where the ram pressure of the shock-heated material equals the thermal pressure of the chromosphere. The first bounce is a transient, therefore, and in fact its amplitude and width are different from those of subsequent bounces in all the simulations examined. On the other hand, Fig. 9 clearly shows that apart from the first bounce the results of the two simulations agree in general quite well, showing differences of < $10\%$ . We expect therefore that, with an increase of the spatial resolution, the results of our simulations may slightly change quantitatively (e.g. the amplitude of oscillations), but not qualitatively.

4 Discussion

4.1 Distribution of emission measure vs. temperature of the shock-heated plasma

As discussed in the introduction, there is a growing consensus that two distinct plasma components contribute to the X-ray emission of CTTSs: the stellar corona and the accretion shocks. This idea has been challenged recently by Argiroffi et al. (2009), who compared the distributions of emission measure EM(T) of two CTTSs with evidence of X-ray emitting dense plasma (MP Mus and TW Hya) with that of a star (TWA 5) with no evidence of accretion (i.e. only the coronal component is present) and with the EM(T) derived from a 1D hydrodynamic model of accretion shocks (i.e. only the shock-heated plasma component is present; Sacco et al. 2008). They proved that the EM(T) of MP Mus and TW Hya can be naturally interpreted as due to a coronal component, that dominates at temperatures of T> 5 MK, plus a shock-heated plasma component, that dominates at T < 5 MK.

In case of shocks with $\beta > 1$ , our simulations show that the distributions of temperature and density at the base of the accretion column can be rather complex (see, for instance, Fig. 7). It is interesting therefore to investigate how the EM(T) of the shock-heated plasma changes with $\beta$ and whether it is possible to derive a diagnostic of the plasma- $\beta$ in the post-shock region.

In order to derive the EM(T) distribution of the accretion region from the models, we first recover the 3D spatial distributions of density and temperature by rotating the corresponding 2D distributions around the symmetry z axis (r=0). The emission measure in the jth domain cell is calculated as ${\rm em}_{j} = n_{{\rm H}j}^2 V_{j}$ , where $n_{{\rm H}j}^2$ is the particle number density in the cell, and V_j is the cell volume. The EM(T) distribution is then derived by binning the emission measure values into slots of temperature; the range of temperature [ $5 < \log T (\mbox{K}) < 8$ ] is divided into 15 bins, all equal on a logarithmic scale ( $\Delta \log T = 0.2$ ).

Figure 10 shows the EM(T) distributions averaged over 3 ks for runs By-01, By-10, and By-50 together with the average EM(T)derived from our 1D reference model HD-1D (blue dashed lines). The figure shows the EM(T) distributions of plasma with a density of $n_{\rm H} > 10^{11}$ cm^-3 (black) and $10^{10} < n_{\rm H} < 10^{11}$ cm^-3 (red). The 2D MHD models have been normalized to have an X-ray luminosity $L_{\rm X}\approx 4\times 10^{29}$ erg in the band [0.5-8.0] keV, in agreement with the luminosity derived from the low-temperature ( $\log T < 6.7$ ) portion of the EM(T) distribution of MP Mus (Argiroffi et al. 2009). This is obtained assuming that $\approx$ 10 accretion streams similar to that modeled here are present simultaneously. The 1D model has been normalized to match the EM peak in By-50. Inspecting Fig. 10 we note that: i) in all the cases, the EM(T) has a peak at $\sim$ 5 MK and a shape compatible with those observed in MP Mus and TW Hya and attributed to shock-heated material (e.g. Argiroffi et al. 2009); ii) most of the X-ray emission arises from shock-heated plasma with a density of >10¹¹ cm^-3 regardless of the $\beta$ (i.e. the most dense component of the post-shock region dominates); iii) the slope of the ascending branch of the EM(T) distribution is comparable in runs By-50 and HD-1D, and gets steeper for decreasing values of $\beta$ . The time-averaged X-ray luminosity derived from these EM(T) distributions ranges between $5\times 10^{29}$ erg (run By-10) and $8\times 10^{29}$ erg (run By-01), with only a weak dependence on the plasma $\beta$ .

$\begin{figure} \par\includegraphics[width=8.3cm,clip]{13565f10.ps} \end{figure}$	Figure 10: Distributions of emission measure vs. temperature, averaged over 3 ks, for runs By-01, By-10, and By-50. Black (red) lines mark the average EM(T) distributions of plasma with density $n_{\rm H} > 10^{11}$ cm^-3 ( $10^{10} < n_{\rm H} < 10^{11}$ cm^-3); dashed blue lines mark the average EM(T) distribution for the 1D reference model HD-1D.
Open with DEXTER

Our model consequently supports the idea that dense shock-heated plasma may contribute significantly to the low temperature portion of the EM(T) distributions of CTTSs regardless of the $\beta$ value. We also suggest that the shape of the EM(T) could be used as a diagnostic of $\beta$ in the post-shock region, if the coronal contribution to the low temperature tail can be neglected.

4.2 Mass accretion rates

Time-dependent models of radiative accretion shocks provide a convincing theoretical support to the hypothesis that soft X-ray emission from CTTSs arises from shocks due to the impact of the accretion columns onto the stellar surface. But, several points remain still unclear. Among these, the most puzzling is probably the fact that the mass accretion rates derived from X-rays, M $_{\rm X}$ , are consistently lower by one or more orders of magnitude than the corresponding M values derived from UV/optical/NIR observations (e.g. Günther et al. 2007; Argiroffi et al. 2009; Schmitt et al. 2005; Drake 2005, Curran et al. 2010, in preparation). We have here the opportunity to discuss the problem of the accretion rate, in particular focusing on the cases with $\beta \ga 1$ in the post-shock region.

The model parameters adopted in this paper describe a stream with an accretion rate of $\approx$ $9\times 10^{-12}~ M_{\rm\odot}$ yr^-1. According with the discussion in Sect. 4.1, $\approx$ 10 streams are needed to match the soft X-ray luminosity of MP Mus. In this case the accretion rate is $\approx$ $9\times 10^{-11}~ M_{\rm\odot}$ yr^-1, which agrees nicely with that deduced from observations, namely M $_{\rm X} = 5{-}8\times 10^{-11}~ M_{\rm\odot}$ yr^-1(Argiroffi et al. 2007). Alternatively, the observed M $_{\rm X}$ may be reproduced by our model if the accretion stream has a larger cross section (with a radius of $r_{\rm str} \approx 10^{10}$ cm). In this case, we expect some changes to the dynamics of the shock-heated plasma as described in Sect. 3.3.

On the other hand, the mass accretion rate of MP Mus, as deduced from optical observations, is M $_{\rm opt} \approx 3\times 10^{-9}~M_{\rm\odot}$ yr^-1 (Argiroffi et al. 2009), which exceeds by more than one order of magnitude the value obtained from X-rays. Similar discrepancies are found in all CTTSs for which it is possible to derive M $_{\rm X}$ (see also Curran et al. 2010, in preparation). The discrepancy might be reconciled if M $_{\rm X}$ values are underestimated due for instance to absorption from optically thick plasma. But, as we explain below, even assuming that the absorption can account for the observed M discrepancy, the idea that the same streams determine both M $_{\rm opt}$ and M $_{\rm X}$ has to be discarded. Indeed, assuming the accretion parameters adopted here, we derive that $\approx$ 300 streams must be present to match M $_{\rm opt}$ or the stream should have a cross section $\approx$ $2\times 10^{22}$ cm², implying that in both cases $\approx$ 20% of the stellar surface would have to be involved in accretion. None of these hypotheses is realistic, as the surface filling factor of hot spots due to accretion is up to a few percent (e.g. Gullbring 1994; Hartmann & Kenyon 1990; Bouvier & Bertout 1989; Bouvier et al. 1993; Drake et al. 2009; Kenyon et al. 1994 and references therein).

The above inconsistency can be removed if the mass density of accretion streams is higher than assumed here; for instance, for a density of the stream $n_{\rm str0} = 5\times 10^{12}$ cm^-3 (a factor 50 higher than modeled), the mass accretion rate is $\approx$ $5\times10^{-10}~M_{\rm\odot}$ yr^-1, and few streams are needed to match M $_{\rm opt}$ . In this case, however, the model cannot explain the lower values of density derived from X-ray observations ( $n_{\rm str0} \approx 10^{11}$ cm^-3; Argiroffi et al. 2007).

A possible solution to remove the discrepancy in the M values (and in the $n_{\rm str0}$ values) is that few accretion streams characterized by different mass density values are present at the same time: those with a density of $\approx$ 10¹¹ cm^-3 would produce shocks that are visible in the X-ray band, leading to the observed M $_{\rm X}$ ; those with higher densities would produce shocks not visible in the X-ray band and leading to the observed M $_{\rm opt}$ (being the dominant component in the UV/optical/NIR bands). The reason of the different visibility could be due to local absorption of the X-ray emission, as explained below.

Our simulations show that the accretion stream penetrates the chromosphere to the position at which the ram pressure of the post-shock plasma equals the thermal pressure of the chromosphere (see also Sacco et al. 2008). Part of the shock-heated plasma is buried in the chromosphere and is expected to be obscured by significant absorption from optically thick plasma (see also Sacco et al. 2010, in preparation). In the simulations presented here, this portion is 1/3 of the hot slab ( $h_{\rm sink} \approx 5\times 10^8$ cm; see Fig. 8), and most of the post-shock plasma (above the chromosphere) is expected to be visible with minimum absorption. Instead, in denser streams, the shock column is buried more deeply in the chromosphere, due to the larger ram pressure, and its maximum thickness is smaller, according to Eq. (11) (see Sect. 3.3). Assuming the accretion parameters adopted in this paper, but with $n_{\rm str0} = 5\times 10^{12}$ cm^-3, Eq. (11) gives $D_{\rm slab} \approx 5\times 10^7$ cm, that is much smaller than the expected sinking of the stream in the chromosphere ( $h_{\rm sink} \approx 7\times 10^8$ cm). In this case, the post-shock column is buried under a hydrogen column density $N_{\rm H} = n_{\rm H} h_{\rm sink} \approx 2\times 10^{22}$ cm^-2, and the photoelectric absorption of the O VII triplet (i.e. the lines commonly used to trace the accretion in the X-ray band) is given by $\exp[-\sigma_{\rm OVII}N_{\rm H}]\approx 10^{-6}$ . We conclude accordingly that heavy streams may produce X-ray emitting shocks that are still hardly visible in X-rays, because they are buried too deeply in the chromosphere.

4.3 Variability of X-ray emission from shock-heated plasma

Another important point in the study of accretion shocks in CTTSs is the periodic variability of X-ray emission due to the quasi-periodic oscillations of the shock position induced by cooling, predicted by time-dependent 1D models. Still in the only case analyzed up to date, namely TW Hya, no evidence of periodic variations of soft X-ray emission (thought to arise predominantly in an accretion shock) has been found (Drake et al. 2009). This result apparently contradicts the prediction of current 1D models, and Drake et al. suggested that these models might be too simple to explain the 3D shock structure.

On the other hand, quasi-periodic shock oscillations are expected if the accretion stream is homogeneous and constant (no variations of mass density and velocity). But, there is substantial observational evidence that the streams are clumped and inhomogeneous (e.g. Bouvier et al. 2007,2003; Safier 1998; Gullbring et al. 1996). In these conditions, periodic shock oscillations are expected to be hardly observable.

The 2D MHD simulations presented here show that, even assuming constant stream parameters, periodic oscillations are not expected if $\beta \gg 1$ in the post-shock region (runs By-01 and By-10). In these cases the time-space plots of temperature evolution (top and middle panels in Fig. 8) actually predict that the oscillations may be rapidly dumped, approaching a quasi-stationary state with no significant variations of the shock position, or the variability may be chaotic (with no obvious periodicity) due to an intense perturbation of the stream by the accreted material ejected sideways. In none of these cases we therefore expect to observe periodic modulation in the X-ray emission.

At the other extreme, for shocks with $\beta \ll 1$ , we expect that the single accretion stream is structured in several fibrils, each independent on the others due to the strong magnetic field which prevents mass and energy exchange across magnetic field lines. Time-dependent 1D models describe one of these fibrils. As they are independent from each other, the fibrils can be characterized by both a slightly different mass density and velocity, which would result in different instability periods, as also by random phases of the oscillations. Drake et al. (2009), assuming that 1D models describe a single stream, derived the number of streams needed to account for the absence of a periodic variability in TW Hya and concluded that this number contrasts with the presence of conspicuous rotationally modulated surface flux with a small filling factor. Following Drake et al. (2009), but considering the fibrils instead of the streams, we may argue that an accretion stream consisting of 200-300 different fibrils with a radius of $r_{\rm fibr}\approx 10^8$ cm and with different instability periods and random phases would produce a signal pulsed at a level of less than 5% as measured in TW Hya.

The intermediate situation is for shocks with $\beta$ around 1. In this case, our run By-50 shows quasi-periodic oscillations of the shock position (see bottom panel in Fig. 8) and predicts periodic variations of the X-ray emission arising from a single stream. In fact, the magnetic field is strong enough in this case to confine the shock-heated plasma, but it is too weak to consider as valid the 1D approximation inside the slab that cannot be described as a bundle of fibrils. As a result, the shock oscillates coherently in the slab. As noted by Drake et al. (2009), in this case it is not possible to reproduce the absence of the periodic modulation of X-ray emission observed in TW Hya, and we conclude that shocks with $\beta \approx 1$ do not occur in this star.

It is worth to emphasize that further investigation is needed to understand to which degree the intermediate case described by run By-50 is frequent and observable: an intensive simulation campaign is needed to assess the range of $\beta$ values leading to streams with quasi-periodic oscillations; a systematic analysis of the variability of soft X-ray emission should also be performed, considering a complete sample of CTTSs with evidence of X-ray emitting accretion shocks, to assess if and when periodic variability is observed.

4.4 Effects of accretion shocks on the surrounding stellar atmosphere

In the case of shocks with $\beta \gg 1$ , our simulations predict that the stellar atmosphere around the region of impact of the stream can be heavily perturbed by the impact, leading to the generation of MHD waves and plasma motion parallel to the stellar surface: the larger the plasma $\beta$ in the post-shock region, the larger the perturbation of the atmosphere around the shocked slab. The resulting ejected flow also advects the weak magnetic field in such a way that the conditions for ideal MHD may break down, magnetic reconnection may be possible and eventually a release of the stored energy from the magnetic field may occur (which is not described by our model however, as it does not include resistivity effects).

A possible effect of the perturbation of the stellar atmosphere is that shocks with $\beta \gg 1$ may contribute to the stellar outflow. In fact, Cranmer (2008) suggested that the MHD waves and the material ejected from the stream (as in our runs By-01 and By-10) may trigger stellar outflow and proposed a theoretical model of accretion-driven winds in CTTSs (see, also, Cranmer 2009). His model originates from a description of the coronal heating and wind acceleration in the Sun and includes a source of wave energy driven by the impact of accretion streams onto the stellar surface (in addition to the convection-driven MHD turbulence which dominates in the solar case). The author found that this added energy seems to be enough to produce T Tauri-like mass loss rates. It would be interesting to assess how the different plasma- $\beta$ cases discussed in this paper contribute to the added energy.

4.5 Limits of the model

Our simulations were carried out in 2D cylindrical geometry, implying that all quantities are cyclic on the coordinate $\phi$ . This choice is expected to affect some details of the simulations, but not our main conclusions. In particular, adopting a 3D Cartesian geometry, the simulations would provide an additional degree of freedom for hydrodynamic and thermal instabilities, increasing the complexity of cooling zones and thus the cooling efficiency of the plasma in the post-shock region. As a result the amplitude of the shock oscillations might be slightly smaller in 3D simulations and the frequency slightly higher than that observed in 2D simulations.

As discussed in Sect. 3.3, some details of our simulations depend on the choice of the model parameters. For instance, the temperature of shock-heated plasma, the stand-off height of the hot slab and also the sinking of the stream in the chromosphere depend on the stream density and the velocity at impact. The accretion parameters adopted in this paper originate from the values derived by Argiroffi et al. (2007) for MP Mus; the cases presented here are representative of a regime in which the shock-heated plasma has a temperature of $\approx$ 5 MK and part of the slab is above the chromosphere, which renders it observable with minimum absorption. The range of the magnetic field strength considered has been chosen in order to explore shocks with $\beta \ga 1$ . Nevertheless, our results undoubtedly show that the stability and dynamics of accretion shocks strongly depends on $\beta$ , and that a variety of phenomena not described by 1D models arise, including the generation of plasma motion parallel to the stellar surface and MHD waves.

We wish to stress that our model does not include magnetic resistivity effects. The complex shock evolution described by our simulations shows that violent plasma flows may advect the magnetic field in a way that the conditions for ideal MHD may break down and magnetic reconnection may possibly occur (not described by our simulations). In particular, field lines at the stream border are squeezed close together and may reconnect, resulting in a change of the magnetic topology as well as in a release of magnetic energy. Since most of these events are expected to occur at the stream border, magnetic reconnection may play an important role in the dynamics and energetic of the ejected accreted plasma. For instance, the magnetic reconnection may trigger flares at the stream border outside the accretion column. This scenario does not contradict Reale (2009), who showed that it is unlikely that flares can be triggered in an accreting flux tube. This issue deserves further investigation in future studies.

5 Conclusion

We investigated the stability and dynamics of accretion shocks in CTTSs, considering the case of $\beta \ga 1$ in the post-shock region, through numerical MHD simulations. To our knowledge, the simulations presented here represent the first attempt to model 2D accretion shocks that simultaneously include magnetic fields, radiative cooling, and magnetic-field-oriented thermal conduction. Our findings lead to several conclusions:

1.: In all the cases, a hot slab of shock-heated material is generated at the base of the accretion column due to the impact of the stream with the chromosphere. In the case of shocks with $\beta > 10$ , violent outflows of shock-heated material and possibly of MHD waves are generated at the border of the hot slab, and they may perturb the surrounding stellar atmosphere. For shocks with $\beta \approx 1$ , the shock-heated plasma is efficiently confined by the magnetic field, and no outflow forms.
2.: If the magnetic field is too weak to confine the shock-heated plasma, but is powerful enough to keep it close to the stream, the escaped accreted material may strongly perturb the accretion column, modifying the dynamics and stability of the shock itself.
3.: The accretion shocks are susceptible to radiative shock instability. The resulting shock oscillations strongly depend on the plasma $\beta$ : for $\beta > 10$ , the oscillations may be rapidly dumped by the magnetic field, approaching a quasi-stationary state, or may be chaotic with no obvious periodicity due to perturbation of the stream induced by the post-shock plasma itself; for $\beta$ around 1 the oscillations are quasi-periodic with an amplitude smaller and a frequency higher than those predicted by 1D models.

We also suggest that the M discrepancy observed in CTTSs may be solved if few streams with a significantly different mass density are present: heavy streams are not visible in X-rays due to heavy absorption and determine the M value deduced from observations in UV/optical/NIR bands; light streams can be observed in X-rays and determine the M value deduced from X-ray observations. As for the absence of periodic X-ray modulation due to shock oscillations found in TW Hya by Drake et al. (2009), we find that no periodic modulation of X-rays is expected in cases of shocks with either $\beta \ll 1$ or $\beta \gg 1$ , whereas periodic modulation is found for $\beta$ around 1, indeed a rather limited set of possible cases. We interpret the absence of periodic X-ray modulation in TW Hya as evidence of the fact that in this star $\beta$ differs significantly from 1 in the accretion shocks, or that otherwise the streams are inhomogeneous.

Acknowledgements

We thank the referee for constructive and helpful criticism. We are grateful to Andrea Mignone and Titos Matsakos for their support in using the PLUTO code. PLUTO is developed at the Turin Astronomical Observatory in collaboration with the Department of General Physics of the Turin University. The simulations have been executed at CINECA (Bologna, Italy), and at the HPC facility (SCAN) of the Osservatorio Astronomico di Palermo. This work was supported in part by the EU Marie Curie Transfer of Knowledge program PHOENIX under contract No. MTKD-CT-2005-029768 and by Agenzia Spaziale Italiana under contract No. ASI-INAF I/088/06/0.

References

Anders, E., & Grevesse, N. 1989, Geochim. Cosmochim. Acta, 53, 197 [NASA ADS] [CrossRef] [Google Scholar]
Argiroffi, C., Maggio, A., & Peres, G. 2007, A&A, 465, L5 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
Argiroffi, C., Maggio, A., Peres, G., et al. 2009, A&A, 507, 939 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
Balsara, D. S., & Spicer, D. S. 1999, J. Comp. Phys., 149, 270 [NASA ADS] [CrossRef] [MathSciNet] [Google Scholar]
Bodo, G., Massaglia, S., Ferrari, A., & Trussoni, E. 1994, A&A, 283, 655 [NASA ADS] [Google Scholar]
Borkowski, K. J., Shull, J. M., & McKee, C. F. 1989, ApJ, 336, 979 [NASA ADS] [CrossRef] [Google Scholar]
Bouvier, J., & Bertout, C. 1989, A&A, 211, 99 [NASA ADS] [Google Scholar]
Bouvier, J., Cabrit, S., Fernandez, M., Martin, E. L., & Matthews, J. M. 1993, A&A, 272, 176 [NASA ADS] [Google Scholar]
Bouvier, J., Grankin, K. N., Alencar, S. H. P., et al. 2003, A&A, 409, 169 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
Bouvier, J., Alencar, S. H. P., Boutelier, T., et al. 2007, A&A, 463, 1017 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
Calvet, N., & Gullbring, E. 1998, ApJ, 509, 802 [NASA ADS] [CrossRef] [Google Scholar]
Chanmugam, G., Langer, S. H., & Shaviv, G. 1985, ApJ, 299, L87 [NASA ADS] [CrossRef] [Google Scholar]
Chevalier, R. A., & Imamura, J. N. 1982, ApJ, 261, 543 [NASA ADS] [CrossRef] [Google Scholar]
Cowie, L. L., & McKee, C. F. 1977, ApJ, 211, 135 [NASA ADS] [CrossRef] [Google Scholar]
Cranmer, S. R. 2008, ApJ, 689, 316 [NASA ADS] [CrossRef] [Google Scholar]
Cranmer, S. R. 2009, ApJ, 706, 824 [NASA ADS] [CrossRef] [Google Scholar]
Dopita, M. A., & Sutherland, R. S. 1996, ApJS, 102, 161 [NASA ADS] [CrossRef] [Google Scholar]
Drake, J. J. 2005, in 13th Cambridge Workshop on Cool Stars, Stellar Systems and the Sun, ed. F. Favata, G. A. J. Hussain, & B. Battrick, ESA Spec. Pub., 560, 519 [Google Scholar]
Drake, J. J., Ratzlaff, P. W., Laming, J. M., & Raymond, J. 2009, ApJ, 703, 1224 [NASA ADS] [CrossRef] [Google Scholar]
Fadeyev, Y. A., Le Coroller, H., & Gillet, D. 2002, A&A, 392, 735 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
Giuliani, J. L. 1984, ApJ, 277, 605 [NASA ADS] [CrossRef] [Google Scholar]
Gullbring, E. 1994, A&A, 287, 131 [NASA ADS] [Google Scholar]
Gullbring, E., Barwig, H., Chen, P. S., Gahm, G. F., & Bao, M. X. 1996, A&A, 307, 791 [NASA ADS] [Google Scholar]
Günther, H. M., Liefke, C., Schmitt, J. H. M. M., Robrade, J., & Ness, J.-U. 2006, A&A, 459, L29 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
Günther, H. M., Schmitt, J. H. M. M., Robrade, J., & Liefke, C. 2007, A&A, 466, 1111 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
Hartmann, L. W., & Kenyon, S. J. 1990, ApJ, 349, 190 [NASA ADS] [CrossRef] [Google Scholar]
Hujeirat, A. 2005, Comput. Phys. Commun., 168, 1 [NASA ADS] [CrossRef] [Google Scholar]
Hujeirat, A., & Camenzind, M. 2000, A&A, 362, L41 [NASA ADS] [Google Scholar]
Hujeirat, A., & Papaloizou, J. C. B. 1998, A&A, 340, 593 [NASA ADS] [Google Scholar]
Imamura, J. N. 1985, ApJ, 296, 128 [NASA ADS] [CrossRef] [Google Scholar]
Johns-Krull, C. M., Valenti, J. A., Hatzes, A. P., & Kanaan, A. 1999, ApJ, 510, L41 [NASA ADS] [CrossRef] [Google Scholar]
Kashyap, V., & Drake, J. J. 2000, Bull. Astron. Soc. India, 28, 475 [NASA ADS] [EDP Sciences] [Google Scholar]
Kastner, J. H., Huenemoerder, D. P., Schulz, N. S., Canizares, C. R., & Weintraub, D. A. 2002, ApJ, 567, 434 [NASA ADS] [CrossRef] [Google Scholar]
Kenyon, S. J., Hartmann, L., Hewett, R., et al. 1994, AJ, 107, 2153 [NASA ADS] [CrossRef] [Google Scholar]
Koldoba, A. V., Ustyugova, G. V., Romanova, M. M., & Lovelace, R. V. E. 2008, MNRAS, 388, 357 [NASA ADS] [CrossRef] [Google Scholar]
Lamzin, S. A. 1998, Astron. Rep., 42, 322 [NASA ADS] [Google Scholar]
Langer, S. H., Chanmugam, G., & Shaviv, G. 1981, ApJ, 245, L23 [NASA ADS] [CrossRef] [Google Scholar]
Langer, S. H., Chanmugam, C., & Shaviv, G. 1982, ApJ, 258, 289 [NASA ADS] [CrossRef] [Google Scholar]
Mignone, A. 2005, ApJ, 626, 373 [NASA ADS] [CrossRef] [Google Scholar]
Mignone, A., Bodo, G., Massaglia, S., et al. 2007, ApJS, 170, 228 [NASA ADS] [CrossRef] [Google Scholar]
Miyoshi, T., & Kusano, K. 2005, J. Comp. Phys., 208, 315 [NASA ADS] [CrossRef] [Google Scholar]
Orlando, S., Peres, G., Reale, F., et al. 2005, A&A, 444, 505 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
Orlando, S., Bocchino, F., Reale, F., Peres, G., & Pagano, P. 2008, ApJ, 678, 274 [NASA ADS] [CrossRef] [Google Scholar]
Reale, F. 2009, in Protostellar Jets in Context, ed. K. Tsinganos, T. Ray, & M. Stute, 179 [Google Scholar]
Robrade, J., & Schmitt, J. H. M. M. 2007, A&A, 473, 229 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
Sacco, G. G., Argiroffi, C., Orlando, S., et al. 2008, A&A, 491, L17 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
Safier, P. N. 1998, ApJ, 494, 336 [NASA ADS] [CrossRef] [Google Scholar]
Schmitt, J. H. M. M., Robrade, J., Ness, J.-U., Favata, F., & Stelzer, B. 2005, A&A, 432, L35 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
Smith, R. K., Brickhouse, N. S., Liedahl, D. A., & Raymond, J. C. 2001, ApJ, 556, L91 [NASA ADS] [CrossRef] [Google Scholar]
Spitzer, L. 1962, Physics of Fully Ionized Gases (New York: Interscience) [Google Scholar]
Sutherland, R. S., & Dopita, M. A. 1993, ApJS, 88, 253 [NASA ADS] [CrossRef] [Google Scholar]
Sutherland, R. S., Bicknell, G. V., & Dopita, M. A. 2003a, ApJ, 591, 238 [NASA ADS] [CrossRef] [Google Scholar]
Sutherland, R. S., Bisset, D. K., & Bicknell, G. V. 2003b, ApJS, 147, 187 [NASA ADS] [CrossRef] [Google Scholar]
Telleschi, A., Güdel, M., Briggs, K. R., Audard, M., & Scelsi, L. 2007, A&A, 468, 443 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
Testa, P., Drake, J. J., & Peres, G. 2004, ApJ, 617, 508 [NASA ADS] [CrossRef] [Google Scholar]
Toth, G., & Draine, B. T. 1993, ApJ, 413, 176 [NASA ADS] [CrossRef] [Google Scholar]
Valenti, J. A., & Johns-Krull, C. M. 2004, Ap&SS, 292, 619 [NASA ADS] [CrossRef] [Google Scholar]
Zel'Dovich, Y. B., & Raizer, Y. P. 1967, Physics of shock waves and high-temperature hydrodynamic phenomena, ed. W. D. Hayes, & R. F. Probstein (New York: Academic Press) [Google Scholar]

Online Material

13565mv1.gif

13565mv2.gif

13565mv3.gif

Footnotes

... stars: Three movies are only available in electronic form at http://www.aanda.org
... measure: Note that the radiative losses are dominated by emission lines in the temperature regime common to CTTSs and are set equal to zero for T< 10⁴ K.
... chromosphere: Note that the radiative losses are set to zero in the chromosphere to keep it in equilibrium.
... surface: Note that the stream velocity is $u_{\rm str0} = -580$ km s^-1 when the stream impacts onto the chromosphere, due to gravity.
... conditions: Variables are symmetrized across the boundary, and both radial and angular $\phi$ components of vector fields ( $\vec{u}, \vec{B}$ ) change their sign.
... outflow: Set zero gradients across the boundary.
... resolution: Note that in the case of simulations with higher spatial resolution (very CPU time consuming), the time interval covered was $\sim$ 1 ks instead of 3 ks to reduce the computational cost.
... luminosity: The synthetic X-ray spectra in the band [0.5-8.0] keV have been derived with the PINTofALE spectral code (Kashyap & Drake 2000) with the APED V1.3 atomic line database (Smith et al. 2001), assuming the metal abundances adopted in the whole paper, namely 0.5.

All Tables

Table 1: Relevant parameters of the simulations.

All Figures

$\begin{figure} \par\includegraphics[width=8.2cm,clip]{13565fg1.ps} \end{figure}$	Figure 1: Radiative losses for an optically thin plasma from the APED V1.3 atomic line database (Smith et al. 2001), assuming the metal abundances of 0.5 of the solar values (Anders & Grevesse 1989).
Open with DEXTER
In the text

$\begin{figure} \par\includegraphics[width=8.cm,clip]{13565fg2.ps} \end{figure}$	Figure 2: Initial geometry of the system in cylindrical coordinates. The stellar surface lies on the r axis, and the unperturbed stellar magnetic field is uniform and oriented along the z axis (vertical lines). The accretion stream propagates downwards through the stellar corona with a velocity $u_{\rm str}$ . The z-axis is the symmetry axis of the problem.
Open with DEXTER
In the text

$\begin{figure} \par\includegraphics[width=8.5cm,clip]{13565fg3.ps} \end{figure}$	Figure 3: Initial hydrogen number density (dashed line) and temperature (dotted line) as a function of height above the stellar surface z for the unperturbed stellar atmosphere.
Open with DEXTER
In the text

$\begin{figure} \par\includegraphics[width=8.5cm,clip]{13565fg4.eps} \end{figure}$	Figure 4: Time-space plot of the temperature evolution for run HD-1D. The spatial extent of the shock lies in the vertical direction at any instant in time. The dashed line marks the initial position of the transition region between the chromosphere and the corona.
Open with DEXTER
In the text

$\begin{figure} \par\includegraphics[width=16.5cm,clip]{13565fg5.ps} \end{figure}$	Figure 5: Temperature ( left panels) and plasma $\beta$ ( right panels) distributions in the (r,z) plane in log scale in the simulations By-01, By-10, and By-50 at a time t=530 s. The initial position of the transition region between the chromosphere and the corona is at z=0. The magnetic field is initially oriented along the z axis; the white lines mark magnetic field lines.
Open with DEXTER
In the text

$\begin{figure} \par\includegraphics[width=15cm,clip]{13565fg6.ps} \end{figure}$	Figure 6: Temperature ( left panel) and mass density ( right panel) distributions in the (r,z) plane in log scale in the simulation By-10, at a time t=1800 s. The white lines mark magnetic field lines.
Open with DEXTER
In the text

$\begin{figure} \par\includegraphics[width=8.5cm,clip]{13565fg7.ps} \end{figure}$	Figure 7: Temperature distribution in the (r,z) plane in linear scale in the simulation By-50 at the labeled times. The white lines mark magnetic field lines.
Open with DEXTER
In the text

$\begin{figure} \par\includegraphics[width=8cm,clip]{13565fg8.eps} \end{figure}$	Figure 8: As in Fig. 4 for the 2D simulations By-01 ( top), By-10 ( middle) and By-50 ( bottom). The blue curve in the bottom panel marks the shock position derived in the 1D reference model, HD-1D.
Open with DEXTER
In the text

$\begin{figure} \par\includegraphics[width=8.5cm,clip]{13565fg9.ps} \end{figure}$	Figure 9: Evolution of the averaged shock-front location for runs By-50 (solid blue line) and By-50-HR (dashed red line). The dotted line marks the initial position of the transition region between the chromosphere and the corona.
Open with DEXTER
In the text

$\begin{figure} \par\includegraphics[width=8.3cm,clip]{13565f10.ps} \end{figure}$	Figure 10: Distributions of emission measure vs. temperature, averaged over 3 ks, for runs By-01, By-10, and By-50. Black (red) lines mark the average EM(T) distributions of plasma with density $n_{\rm H} > 10^{11}$ cm^-3 ( $10^{10} < n_{\rm H} < 10^{11}$ cm^-3); dashed blue lines mark the average EM(T) distribution for the 1D reference model HD-1D.
Open with DEXTER
In the text

Current usage metrics show cumulative count of Article Views (full-text article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.

Data correspond to usage on the plateform after 2015. The current usage metrics is available 48-96 hours after online publication and is updated daily on week days.

Initial download of the metrics may take a while.