L De Sá

J.-P Chièze

C Stehlé

I Hubeny

T Lanz

V Cayatte

F Delahaye

New insight on Young Stellar Objects accretion shocks -a claim for NLTE opacities

Keywords:

Context. Accreted material onto CTTSs is expected to form a hot quasi-periodic plasma structure that radiates in X-rays. Simulations of this phenomenon only partly match with observations. They all rely on a static model for the chromosphere model and on the assumption that radiation and matter are decoupled. Aims. We explore the effects on the structure and on the dynamics of the accretion flow of both a shock-heated chromosphere and of the coupling between radiation and hydrodynamics. Methods. We simulate accretion columns falling onto a stellar chromosphere using the 1D ALE code AstroLabE. This code solves the hydrodynamics equations along with the two first momenta equations for radiation transfer, with the help of a dedicated opacity table for the coupling between matter and radiation. We derive the total electron and ions densities from collisional-radiative NLTE ionisation equilibrium. Results. The chromospheric acoustic heating has an impact on the duration of the cycle and on the structure of the heated slab. In addition, the coupling between radiation and hydrodynamics leads to a heating of the accretion flow and the chromosphere, inducing a possible unburial of the whole column. These two last conclusions are in agreement with the computed monochromatic intensity. Both effects (acoustic heating and radiation coupling) have an influence on the amplitude and temporal variations of the net X-ray luminosity.

Introduction

Classical T Tauri Stars (CTTSs) are solar-type pre-main sequence stars surrounded by a thick disk composed of gas and dust (see e.g. Feigelson & Montmerle 1999). Disk material follows a near-Keplerian infall down to the truncation radius, at which thermal and magnetic pressures balance. Free-falling material flows then from the inner disk down to the stellar surface in magnetically confined accretion columns (Calvet & Gullbring 1998). Hot spots observations (Gullbring et al. 2000) suggest filling factors of up to 1% (Bouvier et al. 1995).

Accreted gas is stopped where the flow ram pressure and the thermal pressure of the stellar chromosphere balance: an accretion shock forms 1 and the post-shock material accumulates at the basis of the column. The hot slab of post-shock material is separated from the accretion flow by a reverse shock. A typical simulated structure of an accretion shock can be found e.g. in Orlando et al. (2010) and is sketched in Figure 1. † deceased. 1 The reverse shock is sometimes confusingly called accretion shock in the literature. One of the most direct probes for the accretion process comes from the X-rays emitted by the dense (n e > 10 11 cm -3 ) and hot (T e ∼ 2-5 MK) post-shock plasma (see e.g. Kastner et al. 2002 andStelzer &[START_REF] Stelzer | [END_REF] for TW Hya, Schmitt et al. 2005 for BP Tau, Günther et al. 2006 A&A proofs: manuscript no. DE_SA-arXiv Argiroffi et al. 2007Argiroffi et al. , 2009 for MP Muscae, Robrade & Schmitt 2007 for RU Lup and Huenemoerder et al. 2007 for Hen 3-600). Another signature is the UV-optical veiling, which is attributed to the post shock medium, the heated atmosphere and the pre-shock medium (Calvet & Gullbring 1998). In addition, Doppler profiles of several emission lines trace the high velocity in the funneled flow (up to 500 km s -1 , according to [START_REF] Muzerolle | [END_REF].

1D hydrodynamical models (Sacco et al. 2008(Sacco et al. , 2010) ) predict Quasi-Periodic Oscillations (QPOs) of the post-shock slab with periods ranging from 0.01 to 1000 s, depending on the inflow density, metallicity, velocity and inclination with respect to the stellar surface. For a typical free-fall radial velocity of 400 km s -1 , Sacco et al. (2010) found for instance a period of 160 s at 10 11 cm -3 . These oscillations are triggered by the cooling instability (for further details, see e.g. Chevalier & Imamura 1982;Walder & Folini 1996;Mignone 2005).

Although plasma characteristics derived from X-ray observations are consistent with the density and the temperature predicted by these numerical studies, there is no obvious observational evidence for such periodicity. [START_REF] Drake | [END_REF] studied thoroughly soft X-ray emission from TW Hydrae and found no periodicity in the range 0.0001-6.811 Hz. Günther et al. (2010) completed this study with optical and UV emission, and they came to the same conclusion in the range 0.02-50 Hz. However, a recent photometric study of TW Hya based on MOST satellite observations reports possible oscillations with a period of 650-1200 s, which could be assigned to post-shock plasma oscillations (Siwak et al. 2018).

Observations thus raise the question of the existence of an oscillating hot slab in the accretion context. Several numerical studies explored multi-dimensional magnetic effects, like leaks at the basis of the column (Orlando et al. 2010), the tapering of the magnetic field [START_REF] Orlando | [END_REF], or perturbations in the flow (Matsakos et al. 2013). Although QPOs are still obtained in these numerical studies, the accretion funnel basis is either fragmented in out-of-phase fibrils, or buried under a cooler and denser gas layer that strongly absorbs X-rays. The observation of global synchronous QPOs becomes therefore very challenging (Curran et al. 2011;Bonito et al. 2014;Colombo et al. 2016;Costa et al. 2017). The effect of the slab burial into the chromosphere has also been explored in several 1D simulations [START_REF] Drake | 13th Cambridge Workshop on Cool Stars, Stellar Systems and the Sun[END_REF]Sacco et al. 2010). Depending on the depth of the burial, the radiation may only escape the post-shock structure from its upper part, leading to a significant reduction of the X-ray luminosity.

In these numerical works, the accretion is supposed to take place on a quiet medium (an isothermal atmosphere in the best cases). Moreover, the post-shock medium is assumed to be optically thin, and the coupling between radiation and matter is reduced to a gas cooling function (see e.g. Kirienko 1993, reported in Figure 3). Although this assumption can be justified to model the infalling gas and the post-shock plasma, it is inconsistent with any stellar atmosphere model. The energy balance between radiation and gas in the lower stellar atmosphere is then replaced by a non-physical tuning (heating function, off threshold, . . . ). Such an assumption may affect the burial of the post-shock structure as well as the accretion structure itself.

In this work, we focus and refine the physics encompassed in existing 1D models. We first explore the effect of chromospheric shocks perturbations on the accretion dynamics. We analyse then how radiation may affect the chromospheric, post-shock and ac-creted plasmas as well as the QPO duration and the hot slab burial; we also synthesise and discuss the accretion signature in the emerging radiative spectra. In Section 2, we present the radiation hydrodynamics model and the numerical tools we use for the hydrodynamics and the spectra synthesis. We detail in Section 2.2.3 the two extreme radiative regimes encountered in this context, and a simple model for intermediate radiative regimes. Section 3 is dedicated to accretion simulations and to the corresponding discussions. The last section (Section 4) presents caveats and possible improvements to this work.

Physical and numerical models

Hydrodynamics model

Hydrodynamics equations

We consider a star of radius R and mass M . The accreted and stellar atmospheric plasmas at position r (r = r ), hereafter taken from the stellar surface, are characterised by a (volumetric mass) density ρ, a velocity u, a thermal pressure p and a volumetric internal energy density e. The plasma evolution is modelled by solving the hydrodynamics equations, written in the conservative form:

           ∂ t ρ +∇• (ρ u) = 0 ∂ t (ρ u) +∇(ρ u ⊗ u) = s m = -∇ (p + p vis ) + g(r) -s M r ∂ t e +∇• (e u) = s e = -∇ • (p u) + q vis -∇ • q C -s E r -q χ (1) with g(r) = -GM ρ/ (R + r) 2 r/r.
The gas source terms2 (s e and s m ) include the contributions of thermal conduction (q C , [START_REF] Spitzer | [END_REF]Vidal et al. 1995), gravity (g(r)), artificial viscosity (p vis and q vis , von Neumann & Richtmyer 1950) and the coupling with radiation (s M r and s E r , see Section 2.2.3). The closure relation for this system of equations -the equation of state -is adapted from the ideal gas law: p = n tot k T ⇔ e = 3/2 p, where n tot stands for the total volumetric number density of free particles (neutrals, electrons and ions), and T represents their kinetic temperature 3 . The contribution of ionisation/recombination on the gas energy density is included in the thermochemistry term q χ , and is discussed in the subsequent section (2.1.2).

Collisional-radiative ionisation

The accretion shock forms where the ram pressure is balanced by the local thermal pressure, i.e. within the stellar chromosphere, that needs then to be modelled. In contrary to the solar case, there is a very limited information about T Tauri chromospheres. Thus, as our goal is to propose a qualitative description of the dynamics of this chromosphere, and in absence of any reliable information, our chromospheric model (see Appendix B) is inspired by the solar case: therefore, we have chosen to use solar parameters in our simulations, and the chemical composition (solar abundances4 ) is then taken from Grevesse & Sauval (1998). In the hydrodynamics, we only consider hydrogen (H i, H ii) and helium (He i, He ii, He iii); the chemical composition is completed by a "catch-all" metal "M"5 . Most simulations are performed using time-independent ionisation models, for instance the modified Saha equilibrium of Brown (1973) (see e.g. Sacco et al. 2008) or a detailed collisional ionisation calculation (e.g. Günther et al. 2007). To estimate the total free electron density n e in the two first setups, we use the modified Saha model (for which q χ = 0). The last simulation presented in this paper (referred to as the Hybrid setup) uses a time-dependent collisional-radiative ionisation model with:

collisional ionisation rates given by Voronov (1997); radiative recombination rates computed by Verner & Ferland (1996); helium dielectronic recombination rate proposed by Hui & Gnedin (1997); photo-ionisation rates (P) derived from [START_REF] Spitzer | Physical Processes in the Interstellar Medium[END_REF] and Yan et al. (1998) cross-sections, and the local radiation energy density. The time dependent ion and neutral volumetric number densities n are then computed by a conservative set of equations (see e.g. ( 1)). The electron volumetric number density is then derived from the neutrality conservation: n e = n H ii + n He ii + 2 n He iii . Finally, the thermochemistry term q χ sums all these contributions, weighted by the corresponding gained/lost energy.

These calculations are performed independently from the opacity computation (see Appendix A), that uses a more refined version of the chemical composition (Grevesse & Sauval 1998).

Radiation model

Radiation and hydrodynamics

The coupling between radiation and matter enters at different scales in astrophysical plasmas. At a microscopic scale, radiation affects the thermodynamical state of the matter through its contribution to the populations of the electronic energy levels of each plasma ion. The computation of these populations is based on a large set of kinetic equilibrium equations that take into account excitation and de-excitation processes due to collisions (interactions with massive particles, mostly electrons) as well as radiative processes (interactions with photons). This step allows to derive also the monochromatic absorption and emission coefficients, resp. κ ν (also called monochromatic opacity, in cm 2 g -1 ) and η ν (in erg cm -3 s -1 ), which in turn are used to compute the local radiation intensity by solving the equations of radiative transfer. Two limiting (and simplifying) cases are expected: at large electron densities, one recovers the Local Thermodynamic Equilibrium (LTE), whereas at low density and for an optically thin medium, the coronal limit is reached [START_REF] Oxenius | Kinetic theory of particles and photons[END_REF].

The main issue in performing such calculations is an intricate coupling between the kinetic equilibrium equations (easily solved given the radiation field), and the radiative transfer equation (simple to calculate knowing the atomic level populations, and hence the absorption and emission coefficients). Since a mean free path of photons is typically much larger than the mean free path of massive particles, an explicit treatment of the radiation transport necessarily involves a significant non-locality of the problem. This issue is satisfactorily solved in the case of stationary stellar atmospheres (see, e.g. [START_REF] Hubeny | Theory of Stellar Atmospheres An Introduction to Astrophysical Non-equilibrium Quantitative Spectroscopic Analysis[END_REF], using efficient iterative methods. However, this remains difficult in the case of a non-stationary plasma, where the equations of hydrodynamics need to be coupled, at each time, with the equations for the radiative transfer.

Therefore, the previous kinetic equations have to be solved simultaneously with the monochromatic radiative transfer equations. This allows computing the frequency-averaged local radiation energy, flux and pressure, and helps including these quantities in the hydrodynamics equations (Eq. ( 1)). In practice, this exact description would require extensive numerical resources: the difficulty is commonly reduced by averaging the radiation quantities by frequency bands. In the multi-groups approximation, the absorption and emission coefficients are averaged over several frequency bands using adapted weighting functions: the larger the number of groups, the better the precision of the computation. The simplest and most commonly used approach is the monogroup approximation, which means that the radiation quantities are averaged over the whole frequency domain covered.

Besides these delicate issues, radiative transfer takes part in the computation of the spectrum emerging from this structure. This is usually done by the post-processing of the hydrodynamic results by more detailed spectral synthesis tools, as detailed in Section 2.3.3.

Moment equations

The radiation field is described here by the momenta equations (see e.g. [START_REF] Mihalas | Foundations of radiation hydrodynamics[END_REF] for the frequency-integrated radiation energy volumetric density (E r , in erg cm -3 ) and momentum (M r , in erg cm -4 s) or flux6 (F r = c 2 M r , in erg cm -2 s -1 ), written in the comoving frame (Lowrie et al. 2001)

:        ∂ t E r + u•∂ t M r + c 2 ∇• M r + (P r : ∇)• u +∇•(E r u) = s E r ∂ t M r + u•∂ t P r /c 2 + ∇•P r + (M r • ∇) u +∇ (M r •u) = s M r (2)
The (monogroup) radiation quantities are integrated from 1 to 10 4 Å. The M1 closure relation allows then to derive the radiation pressure P r from the radiation energy density: P r = D E r . D and χ are respectively the Eddington tensor and factor (D ≡ χ in 1D) and are defined as follows:

D = 1 -χ 2 I 2 + 3χ -1 2 i ⊗ i , χ = 3 + 4 f 2 5 + 2 4 -3 f 2 (3)
with the reduced flux f = F r /(c E r ) (and f = f ), the flux direction i = f / f = F r /F r and I 2 the second-order identity tensor. Contrarily to the radiation energy, the contribution of the radiation flux to the hydrodynamics is not straightforward to interpret7 . Both are presented and discussed with our last setup (Section 3.4.3.2).

Depending on the expression of the radiation source terms, these equations can continuously model optically thin to thick propagation media (see e.g. [START_REF] Mihalas | Foundations of radiation hydrodynamics[END_REF].

Radiation source terms -opacities & line cooling

This work aims at describing in a consistent way the system composed of three zones, which are coupled together through A&A proofs: manuscript no. DE_SA-arXiv radiation but in different thermodynamical states (Figure 1): the dense and optically thick near-LTE chromosphere (Section 2.2.3.1) on the one hand, the optically thin coronal hot accretion slab and cold accretion flow (Section 2.2.3.2) on the other hand. We also expect, according e.g. to Calvet & Gullbring (1998), that the frequency distribution of the measured radiation varies strongly from the X-rays to the infrared. We have decided to work step by step, using a model which makes a continuous transition between the optically thick LTE approximation and the coronal limit, as described in Section 2.2.3.3.

Optically thick limit

The deep stellar atmosphere is optically thick and can be considered at LTE, i.e. each microphysics process is counter-balanced by its reverse process. In LTE and regimes close to LTE, the monochromatic absorption and emission coefficients are linked through the Planck distribution function: η ν = κ ν ρ c B ν . The radiation energy and momentum source terms are then defined by (see e.g. [START_REF] Mihalas | Foundations of radiation hydrodynamics[END_REF]:

s * E r = κ P ρ c a R T 4 -E r and s * M r = -κ R ρ c M r (4)
where a R is the radiation constant. Two radiation-matter coupling factors appear here (in cm 2 g -1 ). The Planck mean opacity κ P is based on the frequency-integrated absorption coefficient κ ν weighted by the Planck distribution function B ν , while the Rosseland mean opacity κ R is the harmonic mean of κ ν weighted by the temperature derivative of the Planck function ∂ T B ν , as follows [START_REF] Mihalas | Foundations of radiation hydrodynamics[END_REF]:

κ P = κ ν B ν dν B ν dν and κ -1 R = κ -1 ν ∂ T B ν dν ∂ T B ν dν (5)
In these frequency averages, the Planck mean is dominated by strong absorption features (typically lines), whereas the Rosseland mean is dominated by the regions in the spectrum of lowest monochromatic opacity. As a consequence, at large optical depths, κ P correctly describes the energy exchange between particles and photons, while κ R gives the correct total radiative flux [START_REF] Hubeny | Theory of Stellar Atmospheres An Introduction to Astrophysical Non-equilibrium Quantitative Spectroscopic Analysis[END_REF]).

-14 -12 -10 -8 -6 log 10 ρ g cm -3 Several opacity tables are available for a variety of chemical compositions. However, they all fail to cover the full (ρ, T ) domain explored in our simulations (see solid black line in Figure 2). We constructed then with the SYNSPEC code (Section 2.3.3) our own LTE opacity table (see Appendix A for further details), presented in Figure 2. These opacities include atomic (high T ) and molecular (low T ) contributions.

Optically thin limit

Due to its very low density (ρ ∼ 10 -13 g cm -3 ), the accreted plasma can be described by the limit regime where the gas density tends towards zero: the coronal regime. The coupling between radiation and matter boils down in this case to an optically thin radiative cooling function Λ(T ) (in erg cm 3 s -1 ). In Eq. ( 2), the radiation source/sink terms become then:

s † E r = n e n H Λ(T ) and s † M r = 0 (6)
The first quantity represents the net radiation power emitted by unit volume in all directions (4π sr) by a hot optically thin plasma (in erg cm -3 s -1 sr -1 ). The term s † M r is set to zero since there is no coupling between radiation and matter in this regime (see Appendix C for more details).

Fig. 3. Optically thin radiative cooling (in erg cm 3 s -1 ) for different metallicities Z, versus gas temperature (K), adapted from Kirienko (1993).

The present work is based on the cooling function provided by Kirienko (1993), reproduced in Figure 3, with Z/Z = 1 (see Appendix B.1 for the explanation).

Intermediate regimes

The previous source terms describe two well-defined plasma situations. On the one hand, the basis of the stellar chromosphere is optically thick and can be described by the previous LTE radiation source terms. On the other hand, the low density and hot slab is mostly optically thin and can be described in the coronal regime. It is physically expected and numerically compulsory to perform a smooth and continuous transition to encompass intermediate regimes. This could be done using adequate opacities and emissivities, as for instance obtained in a collisional-radiative model, unfortunately not available yet for the whole range of physical conditions of the present study. Thus we have preferred to follow the transition between LTE and coronal regimes with the probability for a photon (emitted from the column center) to escape sideways (see e.g. [START_REF] Lequeux | The Interstellar Medium[END_REF], equation 3.66):

ζ = 1 -exp(-3τ e ) 3τ e , τ e = κ R ρ L c (7)
ρ and κ R values are taken at the photon emission position. The characteristic length L c is here taken as the accretion column

s E r = (1 -ζ) s * E r + ζs † E r and s M r = s * M r (8)
the star ( * ) and dagger ( †) denoting respectively the LTE (Eq. ( 4)) and the coronal (Eq. ( 6)) expressions.

Numerical tools

One-dimensional approach

Observations indicate that, in general, the ambient magnetic field is of the order of 1 kG (Johns-Krull et al. 1999;Johns-Krull 2007). The resulting Larmor radius (1 mm) is very small, i.e. the plasma follows the magnetic field lines. Moreover, the Alfvén velocity reaches 3% of the speed of light and the magnetic waves behave thus like usual light waves. Therefore, focusing on the heart of an accretion column in strong magnetic field case, we can model the accreted material along one field line, that will be assumed to be radial relative to the stellar center. Since the accretion process is expected to involve strong shocks, we chose a numerical tool able to achieve very high spatial resolution.

AstroLabE -an ALE code

The present work is based on numerical studies performed with the 1D Arbitrary-Lagrangian-Eulerian (ALE) code AstroLabE (see e.g. [START_REF] De Sá | Proceedings of the Annual meeting of the French Society of Astronomy and Astrophysics[END_REF]Chièze et al. 2012). It is based on the Raphson-Newton solver (Numerical Recipes, Section 9) and a fully implicit scheme (the CFL condition can then be ignored) to compute primary variables at each time step. This code solves, along with the adequate physics and chemistry equations (see Sections 2.1 and 2.2), the equations describing the behaviour of the grid points. The space discretisation can follow an Eulerian or a Lagrangian description. Moreover, the grid can freely adapt to hydrodynamics situations (the arbitrary description, Dorfi & Drury 1987): this helps us reach high resolution around shocks with fixed cardinality (δr/r max ∼ 10 -7 with 150-300 grid points).

Beside its application to stellar accretion (de Sá 2014), Astro-LabE has been used in several astrophysical situations such as the interstellar medium (Lesaffre 2002;Lesaffre et al. 2004), experimental radiative shocks [START_REF] Stehlé | Semaine de l'Astrophysique Francaise[END_REF]Bouquet et al. 2004) or type Ia supernovae (Charignon & Chièze 2013) studies.

SYNSPEC -a spectrum synthesiser

For the computation of the opacities and of the emerging spectra, we used the public 1D spectrum synthesis code SYNSPEC (Hubeny & Lanz 2017). It is a multi-purpose code that can either construct a detailed synthetic spectrum for a given model atmosphere or disks, or generate LTE opacity tables. In this paper, we used SYNSPEC both for generating opacity tables (see Section 2.2.3.1 and Appendix A), and for the snapshots spectra presented in Section 3.4.4. The resulting synthetic spectrum reflects the quality of the input astrophysical model; using an LTE model results in an LTE spectrum, while using a NLTE model results in a NLTE spectrum. The snapshots of our hydrodynamic simulations provide temperature and density as a function of position; it is therefore straightforward to compute LTE spectra for such structures. It would be in principle possible to construct approximate NLTE spectra, keeping temperature and density fixed from the hydrodynamic simulations (the so-called "restricted NLTE problem"). This could be done for instance by the computer program TLUSTY (Hubeny & Lanz 1995, 2017), which would provide NLTE level populations that can be communicated to SYNSPEC to produce detailed spectra. However, as previously mentioned, such a study is computationally very demanding and is well beyond the scope of the present paper. Nevertheless, since NLTE effects may be important, this will be done in a future paper. It will allow to inspect the effect of the LTE approximation on our results.

This synthetic spectrum, computed at different altitudes of the accretion column, will reveal the role played by the different parts of the spectrum, from X-ray to Visible (1-10 4 Å). However, it is important to note that, as the accretion column is limited in diameter, some effects, like the absorption by the coldest parts are only pertinent for an observation along or near the direction of the accretion column. A 3D radiative transfer post-processing would then be more suitable to the geometry of the system (Ibgui et al. 2013).

Accretion basis simulations

Strategy and common parameters

We have simulated for this study several physical situations in order to check the net effect on the QPOs of the chromospheric model on one side and of the matter-radiation coupling on the other side. We present first the reference case: a gas flow hits a fixed, rigid and non-porous interface (W-Λ case, Section 3.2). We check then the effect of a dynamically heated chromosphere on the accretion process (Chr-Λ case, Section 3.3) and we finally check the effect of the radiation feedback on matter (Hybrid case, Section 3.4). The conditions and main results of each simulation are resumed in Table 1.

The simulations presented in this paper share few parameters:

the computational domain size is r out = 10 5 km (the outer boundary limit); the column radius is set to L c = 1000 km; for the gravity magnitude, we use R = R and M = M ; the accreted gas enters the computational domain through the outer boundary with ρ acc = 10 -13 g cm -3 , T acc = 3000 K8 and v acc = 400 km s -1 . The velocity of the accreted gas is derived from the free-fall velocity at r = r out above the stellar surface, considering a null radial velocity at the truncation radius R tr = 2.2 R (taken here from the center of the star).

When the M1 radiation transfer is used (either near-LTE transfer or intermediate regime), one solar surface luminosity (L = 6.3 × 10 10 erg cm -2 s -1 ) enters from the inner boundary, and c × E out r /4 leaves from the outer boundary9 , with E out r being the radiation energy density of the last computational cell.

A&A proofs: manuscript no. DE_SA-arXiv 3.2. Reference case (W-Λ)

Setup

In the reference case, we simulate the accretion stream using the same physics and assumptions than in previous models (see e.g. Sacco et al. 2008;Koldoba et al. 2008). The matter-radiation coupling is then described by the coronal radiative cooling (Section 2.2.3.2) and the plasma ionisation is computed with the modified Saha equation (Section 2.1.2). In order to simplify the discussion, we focus on the post-shock structure and on the global dynamics. The stellar chromosphere is modelled in the simplest way, hereafter called the "window" model. It consists in a fixed rigid non-porous transparent interface. The main parameters are resumed in Figure 4.

QPO cycle

Besides the fact that matter accumulates on the left (inner) rigid boundary interface, the system is found to be perfectly periodic.

Figure 5 presents five snapshots of density, temperature and velocity profiles during a QPO cycle far from the initial stages. The accreted gas falls from right to left. A hot slab of shocked material builds first (t = 2750 and 2884 s) and cools down according to the coronal regime. Below a threshold temperature10 , the fast, quasi-isochoric, cooling of the slab basis causes the collapse of the post-shock structure (t = 2994 and 3110 s). Just after the full collapse of the slab, since the accretion process is still working, a new slab forms and grows (t = 3156 s). This simulation is to be compared to the ones performed by Sacco et al. (2008); Table 2 resumes the main parameters and results for fast comparison. Despite few key differences (Sun vs. n e (cm -3 ) 10 11 -10 12 10 11 -10 11.5

T max (K) 10 6.5 10 6.5

MP Muscae parameters & "window" vs. chromospheric heating function), the results are in good agreement with each other.

X-ray luminosity

An X-ray radiative power of 1.3 × 10 30 erg s -1 was measured in the range 2-27 Å by Brickhouse et al. (2010) for TW Hydrae and an accretion flow velocity estimated at 500 km s -1 . We compute therefore the instantaneous X-ray surface luminosity L Λ (in erg cm -2 s -1 ) and its time average LΛ , to compare them with the values obtained with the different models presented in the subsequent sections and with the observational work of Brickhouse et al.:

L Λ = slab n e n H Λ(T ) dr & LΛ = 1 τ cycle τ cycle 0 L Λ dt (9)
These quantities are commonly compared to the incoming kinetic energy flux. However, since the flow accelerates in its free-fall from the outer boundary down to the reverse shock, the plasma velocity and density may change between the outer boundary of the simulation box and the location of the reverse shock. To get round this issue, one must consider the mechani- cal energy flux. This flux is calculated at any position r by:

F M = 1 2 ρ v(r) 3 + r r 0 G M ρ v z 2 dz (10)
where the origin of the gravitational energy potential is set at the mean accretion shock position (r 0 10 3 km). The conservation of the mechanical energy induces that F M does not depend on the position r. The value derived from our simulations is

F M = 4.2 × 10 9 erg cm -2 s -1 . 0.0 0.2 0.4 0.6 0.8 1.0 t/τ cycle -2.0 -1.5 -1.0 -0.5 0.0 0.5 log 10 L Λ /F M W -Λ Chr -Λ Hybrid
Fig. 6. Time variation of the surface luminosity L Λ emerging from the hot slab for the reference (blue), the dynamical chromosphere (green) and the hybrid (red) cases. These quantities are computed assuming an optically thin coronal plasma. To allow comparison, the time is reported in reduced units of t/τ cycle and the luminosity is normalised to the incoming mechanical energy flux F M defined in Eq. ( 10). The values of the cycle duration for each setup is reported in Table 1.

Figure 6 shows the time variation of L Λ . As expected, this quantity increases during the propagation of the reverse shock and decreases during the collapse. The time-averaged luminosity LΛ is equal to 1.5 × 10 9 erg cm -2 s -1 , i.e. 36% of the incoming mechanical energy flux F M .

Effect of a dynamical chromosphere (Chr-Λ)

Setup

In this second setup (see Figure 7), we aim at studying the effect of a dynamically heated chromosphere on the phenomenon described in the previous Section. To achieve this, we "divide" the computational domain into two zones separated by a transparent11 Lagrangian interface. The outer zone is described as before, i.e. with modified Saha ionisation and optically thin radiative cooling (coronal regime). However, the inner zone is now described by our chromospheric model (see Appendix B). Ionisation is still described by the modified Saha equation, but we use the LTE radiation source terms as given in Eq. (4). To get a dynamically heated chromosphere, we first compute a hydrostatic equilibrium, with the outer zone inactivated, and with one solar luminosity crossing the entire domain (no effect on the outer zone). Acoustic energy is then injected in the form of monochromatic sinusoidal motion of the first interface (a "window") with a 60 s period to mimic solar granulation. Several snapshots of temperature profiles are presented in Figure B.1. Once the shock-heated chromosphere reaches its stationary regime, the accretion process is launched (in the outer zone). very close to the snapshot of the first cycle at t = 71 s. The (unchanged) end of the second cycle is then not reported.

Acoustic perturbations

During the installation phase (1-336 s) of the reverse shock, the post-shock structure follows more or less the same scenario than for the reference case (W-Λ). After several periods of the acoustic waves, small differences occur. The transmission of these waves/shocks to the accretion column depends on the leap of the acoustic impedance between the upper chromosphere and the hot slab, which results in reflection/transmission of these waves/shocks at this interface. The smallest leap is reached at the end of the collapse, near 336 s, leading to a transmission increase, which however remains still low. Their effect leads to small perturbations in the post-shock density (as it can already be seen at 168 s).

Table 3. Position of the old and new reverse shocks between t = 358 s and t = 397 s (see Figure 8).

Time (s)

358 380 386 397 r old (km) 10 3.40 10 3.25 10 3.10 10 2.95 r new (km) 10 3.15 10 3.70 10 3.80 10 3.90

After this time, the transmitted waves start to feed with matter the hot collapsing layer behind the reverse shock. The thickness of this layer increases, as can be shown in Figure 8 at 351 s, compared for instance with our reference case (3110 s, Figure 5). This structure collapses and hits at 354 s the dense chromosphere, leading to a secondary reverse shock which propagates backwards inside the slab. This behaviour is confirmed by the velocity variations shown in grey in Figure 8. The two reverse shock pass then each other: the positions of the new shock (or contact discontinuity) and the previous (old) one are resumed in Table 3. The end of one cycle therefore overlaps the beginning of a new one.

Observational consequences

This model implies two main observational consequences. First, compared to the reference case, the QPO cycle period is modified by the acoustic heating. The question of possible resonance is pointless regarding multi-mode acoustic heating by out-of-phase waves emitted in different locations. The period τ cycle is slightly reduced (from 400 s for the W-Λ model to 350 s here, Table 1) when using solar chromospheric parameters. Since CTTSs' atmospheres have a stronger activity than the Sun's one (that we use for the chromospheric model), the effect is expected to be enhanced in CTTSs. The second effect deals with the X-ray luminosity variation during a cycle, as reported in green in Figure 6. The growth phase is comparable with the W-Λ setup, but the acoustic perturbations from the chromosphere induce strong differences in the collapse phase. Moreover, the overlapping of the beginning and end of the cycles affect the X-ray luminosity and the overall amplitude of the variations (contrast) is reduced compared to the reference case. QPO observations may thus require both higher time resolution and improved sensitivity. The time averaged surface luminosity (Eq. ( 9)) is here equal to 4.0 × 10 9 erg cm -2 s -1 , i.e. 94% of the mechanical energy flux F M .

These results show that, compared to the reference case, the dynamical heating of the chromosphere impacts the duration of the QPO period and its observability. Of course, a more realistic description of the chromospheric heating would require at least a 2D MHD picture. For instance, we know that chromospheric perturbations may lead -inside the column -to the development of fibrils (see e.g. Matsakos et al. 2013, ChrFlx# models), which is one of the scenarii explaining the absence of observation of QPO. In the acoustic description of the chromospheric heating, these fibrils, evolving out of phase, will also be strongly affected by the chromospheric perturbations.

Radiation effect on accretion (Hybrid)

Setup

In this Section, the plasma model includes collisional-radiative ionisation (see Section 2.1.2). The radiation-matter coupling is described within the intermediate regime (see Section 2.2.3.3) and the outer radiation flux is set to c × E out r /4. The goal of this last setup (see Figure 9) is to inspect the net effect of the matterradiation coupling. We have therefore chosen not to consider any chromospheric activity (i.e. the acoustic heating is kept switched off).

Ionisation model

We have tested in this setup the effect of the time-dependent ionisation through radiative ionisation/recombination and collisional ionisation with a time-dependent formulation (see Section 2.1.2 for more details).

The main difference brought by a time-dependent calculation of the electron density is a tiny ionisation delay behind the reverse shock front, as shown in Figure 10. At the shock front, the kinetic energy is converted into thermal energy, and then a part of this thermal energy is used to ionise the post-shock material with a time scale connected to the ionisation rates; the affected gas layer is up to 0.2 km thick, and thus negligible compared to the whole structure (that is at least 10 4 km thick, see Table 1). This justifies the use of a time-independent model for ionisation in the previous setups (W-Λ and Chr-Λ).

However, compared to the Saha-Brown equilibrium calculations, the use of collisional and radiative rates to derive the equilibrium electron density brings differences in the transition between the (almost) neutral medium and the fully ionised plasma. This transition lays between 10 3.6 and 10 4.2 K. However, such temperatures are only reach by the accreted gas during the cooling instability. Its overall effect is hence negligible. The results presented for the Hybrid case (see Figure 11) are thus based on this collisional-radiative equilibrium calculation of n e .

Radiation and ionisation feedback

The first cycle is presented in Figure 11; it shows the time variations over 160 s of the gas temperature and mass density, of the photon escape (ζ) and absorption (1ζ) probabilities and of the radiation energy volumetric density and flux (E r and F r ) (see Section 2.2.3.3) for the same snapshots. The next cycles only differ from this first one by the position of the interface between the slab and the chromosphere, as discussed in Section 3.4.3.4.

The global behaviour follows the trends of the two previous models. However, several effects must be highlighted: a heating of the chromosphere and of the accretion flow, already pointed out by Calvet & Gullbring (1998) and Costa et al. (2017), and the reduction of the oscillation period and of the post-shock extension. These effects are discussed below.

QPO cycle reduction

Although 1ζ shows strong variations, its net value beyond the accretion shock remains negligible, and the post-shock material is in the coronal limit (as in Section 3.2). The temperature behind the reverse shock is here equal to 3.1 × 10 6 K, to be compared to 4 × 10 6 K in the reference case. In addition, the compression is enhanced from 4 (W-Λ case) to 4.4. As a consequence, the cooling is more efficient: the cooling time is reduced from 400 s down to 220 s, which is compatible with the duration of the cycles. This effect is due to the ionisation/recombination energy cost (q χ ), which is included in the gas energy equation for the Hybrid case, but not for reference case (q χ = 0 in W-Λ case, cf. Eq. (1) and Section 2.1.2).

Radiation energy and flux

The radiation energy density increases between 9 and 100 s, which corresponds to the growth phase of the hot slab. This increase is however correlated to the upper chromosphere heated up to 12 000 K (discussed in Section 3.4.3.3 and presented Figure 12). E r remains almost flat in the optically thin post-shock medium, with a value driven by the heated upper chromosphere. In the accretion flow, during the growth of the hot slab, there is The most peculiar feature of the radiative flux is its linear growth through the post-shock slab. Such a pattern is characteristic of a volume emission by an optically thin medium. As discussed before, the radiative energy in the hot slab is somehow imposed by the heated upper chromosphere: as a consequence, the outgoing radiation flux is cE r /4 c a R T 4 chr /4, with T chr the temperature of the upper chromosphere (cf. Figure 12). Then, the radiation flux propagates through the accretion flow with negligible changes. Since F r is counted negatively towards the star, the flux emitted by the slab is offset by the one produced by the chromosphere. the net radiation flux rises then back in the chromosphere. 0 The upper chromosphere is heated by the radiating postshock plasma up to 12 000 K (Figure 12). For instance, between 9 and 70 s, its temperature varies from 7000 to 10 800 K at 800 km, and the pressure increases from 800 to 2600 dyn/cm 2 at this location (Figure 13). As a consequence, the whole postshock structure is pushed upwards from 875 km to 3150 km, thus out of the unperturbed chromosphere (by about 2000 km, see e.g. [START_REF] Vernazza | [END_REF].

At the end of the cycle, the chromosphere is not heated any more and the slab buries back into the atmosphere. The expected behaviour is an oscillation of the slab burial with the same periodicity as QPOs, since it originates from the hot post-shock plasma radiation. At the end of this first cycle, the chromosphere does not recover its initial thickness: this effect does not affect the post-shock dynamics and cycle characteristics. All these effects are overestimated in a 1D model. However, this study shows that the general question of the (un)burial, which is important for X-ray observations, can only be addressed within a model that takes into account the radiative heating of the chromosphere by the hot slab. While reaching the hydro-radiative steady state of the chromosphere, the flow has been homogeneously heated from 3000 K to 5730 K before the start of the accretion process. During the cycle, the accretion leads to an additive heating of the flow up to ∼ 8500 K (at t = 93 s). These effects are quantified in Figure 12, which reports the time variations of the position and velocity of the interface between the hot slab and the accretion flow, as well as the temperatures of the heated chromosphere and of the pre-shock material. Such preheating has already been pointed out by other authors (Calvet & Gullbring 1998;Costa et al. 2017). In these works, this heating is induced by radiation coming from the hot slab through photo-ionisation. Although radiative cooling of the accretion flow may be included in some cases, the radiation transfer is not taken into account. Depending on the conditions, the pre-shock temperature may reach from 20 000 K (in CG98) up to 10 5 K (in Co17) close to the reverse shock (up to 10 4 km). In the latter, this precursor is preceded by a flatter and cooler (∼ 10 4 K) zone with an extension of 10 5 km, thus smaller than ours (> 10 5 km).

Our simulation shows that part of the heating is a consequence of the chromospheric radiation already in play before the start of the accretion. The analysis of the variation of the radiative energy indicates that an additional heating operates during the development of the hot slab. However, as we do not include any dependence with the wavelength, it remains very difficult to discriminate in details the role played by the radiation emitted by the hot slab (X-rays) and from the (heated) chromosphere (UVvisible). Complementary information will be given by the synthetic spectra computed as a post-process of the hydrodynamics structures (Section 3.4.4).

X-ray luminosity

The X-ray luminosity of the system is computed following the method described in Section 3.2.3. Its time variation (in unit of τ cycle = 160 s) is reported in Figure 6 for comparison with the two previous cases. Compared to the reference case, in addition with a shortening of the period, this case presents a more pro-A&A proofs: manuscript no. DE_SA-arXiv nounced radiative collapse (70-90 s), followed by a chaotic collapse (90-160 s). The time average of the radiative surface luminosity is here equal to 1.2 × 10 9 erg cm -2 s -1 , which represents 30% of the mechanical energy flux (cf. Figure 6).

SYNSPEC monochromatic emergent intensity

As this simulation is performed using only one group of radiation frequencies, it is interesting to analyse more precisely the details of the previous radiative heating via its feedback on the monochromatic emergent intensity.

To this purpose, the hydrodynamic structures has been postprocessed with the SYNSPEC code (Section 2.3.3). For consistency purpose, we take the atomic data already used for the calculation of the average opacities (see Section 2.2.3 and Appendix A). We thus estimate the specific intensity I λ (in erg/cm 2 /s/Å/sr) along the direction of the column. Since the line profile behaviour is not investigated here, velocity effects are neglected.

It is important to recall that a quantitative comparison of this synthetic spectrum with observations, especially in the X-rays (see e.g. Güdel et al. 2007;Robrade & Schmitt 2007;[START_REF] Drake | [END_REF]) would require NLTE and 3D radiative transfer postprocessing. Nonetheless, using 1D radiative transfer and the LTE approximation is here interesting as it corroborates or not the general accepted trends, e.g. a strong X-ray emission and an excess of luminosity in the UV-VIS range (Calvet & Gullbring 1998;Brickhouse et al. 2010;Ingleby et al. 2013). A typical spectrum emerging from 4.6 × 10 4 km (located within the accretion flow) is reported in Figure 15). It is computed from a snapshot (t = 70 s) of the Hybrid model (see Figures 11 and 14). At this stage, the chromosphere extends up to 1.4 × 10 3 km, the hot plasma from 1.4 × 10 3 km to 8.3 × 10 3 km and the accretion flow from 8.3 × 10 3 km to 1 × 10 5 km. The in-tensity that emerges from this layer presents three characteristic spectral bands:

in the range 1-100 Å (X-rays), the bump is attributed to the hot post-shock plasma, with intense lines up to 1012 erg/cm 2 /s/Å; in the range 100-900 Å (EUV), radiation is efficiently absorbed by the inflow; in the range 900-10 000 Å (UV+Vis+IR), the second bump is attributed to the heated stellar chromosphere and photosphere, i.e. a black body at T ∼ 11 300 K (cf. Figure 12). The strong absorption of the EUV radiation is due to the huge optical depth of the accretion flow 12 . This effect may then be attenuated in the case of a bent column or when the observation is performed side-on and not along the column. This absorption effect on the spectrum is illustrated in Figure 16, which presents the intensity emerging right after the reverse shock front, at r = 8.3 × 10 3 km. This figure shows that this absorption also affects, to a lesser degree, the visible spectrum originating from the chromosphere. This must be considered when interpreting the UV excess (see e.g. Calvet & Gullbring 1998;Hartmann et al. 2016;Colombo et al. 2019). Note that a pre-heating of the accretion flow is expected as a result of the EUV absorption. A preheating is also obtained independently by AstroLabE (Section 3.4.3.4); however, a one-to-one correspondence would require a multi-group description of the radiation field in AstroLabE. We compute from I λ (Figure 14) the net X-ray outgoing intensity (I X ) and the corresponding coronal quantity (I Λ ):

I X = 27 2 I λ dλ and I Λ = L Λ 4π (11)
The time variations of these two quantities, reported in Figure 17, present similar characteristics. However, the values derived by SYNSPEC are higher by about two to three orders of magnitude. This discrepancy is either imputable to the LTE approximation or to the assumed 1D plane-parallel geometry assumed in SYNSPEC. Thus our synthetic spectra can't be used for quantitative comparison with observations.

Refining the models

A more realistic chromosphere

It should be pointed out that this study uses a solar model for the chromosphere with acoustic heating. Compared to the description of this heating, a more important improvement would be to consider a realistic T Tauri chromospheric model, which is today not very well known. This may affect the ionisation (and then gas pressure with another chemical abundances) as well as slab characteristics (through gravity) and radiation effects (through opacities and incoming luminosity). Our results are then to be considered qualitatively and not quantitatively.

Improvements of the radiation model

We use in this work radiation momenta equations with the M1 closure relation. Although this is already a strong improvement compared to other approaches like the diffusion model, it could be improved by using radiation half-fluxes (i.e. the inward an outward components of the radiation flux). This should disentangle the radiation flux coming from the star and from the postshock structure. The M1 closure relation allows the radiation field to reach at most one direction of anisotropy; half-fluxes can extend it to two, i.e. the maximum number of anisotropy directions reachable in 1D. Half-fluxes (along with M1) would then be equivalent to the momenta equations with the M2 closure relation (Feugeas 2004), without its prohibitive numerical cost. The M1 model and its limits have been thoroughly studied (see e.g. Levermore 1996;Dubroca Feugeas 1999;Feugeas 2004).

The behaviour of this model along with half-fluxes needs however to be examined.

More important is the approximation made with the monogroup approach used in this work. The whole spectrum is then approximated as a black body providing the adequate opacity averages. However, our computed spectra emerging from accretion structures are expected to present three discernible frequency groups:

up to the visible domain, the spectrum is dominated by the black body emerging from the stellar photosphere; the EUV band is expected to be depleted due to high absorption by the accreted gas; the X-ray band is thought to be optically thin and to have the accretion shock signature on it.

Although the multi-group approach is numerically heavier, it will improve the study of the consequences of the radiation absorption by the surrounding medium. A consequence of the X-ray and EUV absorption by the cold accretion flow is the presence of a radiative precursor. Such a phenomenon cannot be obtained through a monogroup approach. Moreover, a 3 groups approach will provide a better description of the feedback of the hot slab on the stellar chromosphere.

NLTE effects in radiation hydrodynamics and in synthetic spectra

Two other points may be improved. First, the transition model (ζ) remains qualitative and may need to be extended to the ionisation calculation. The work done by Carlsson & Leenaarts (2012) offers paths to reach such consistency and may need to be investigated further. A better model of both the LTE transfer, line cooling and intermediate regimes may demand dedicated NLTE opacities, namely plasma emissivity (equivalent to n e n H Λ), radiation energy absorption (κ P ) and radiation flux sinking (κ R ). Moreover, all these quantities, computed with a radiative-collisionnal model, have to be averaged over adequate weighting functions. Due to recent progresses in this topic (Rodriguez et al. 2018), new results are expected in a near future. Independently, a NLTE description should be used to compute the emerging spectra: this work is already in progress using TLUSTY code.

Conclusion

In this study, we used 1D simulations with detailed physics to check the validity of the two following common assumptions in accretion shock simulations: the stellar atmosphere can be either modelled by a hydrostatic or a steady hydrodynamic structure, and the dynamics of accretion shocks is governed by optically thin radiation transfer. We checked first that we are able to recover previous results (Sacco et al. 2008, W-Λ case, Section 3.2) and tested independently each of these assumptions (Chr-Λ case, Section 3.3, and Hybrid case, Section 3.4). Each of them proves to have a non-negligible impact on the typical characteristics of the accretion dynamics and on the estimation of its X-ray surface luminosity. This one varies between 1 and 4 × 10 9 erg cm -2 s -1 . Taking as a reference the radiative power of Brickhouse et al. (2010), we derive a section of the accretion spot from 3 × 10 20 to 1 × 10 21 cm 2 , corresponding to a filling factor of the solar disk between 2 and 8%.

In the case of the chromosphere which is heated by acoustic perturbations that degenerate into small shock waves, we have shown that these perturbations do not strongly modify the cycle period compared to the reference case. However, the cycle becomes chaotic due to the generation of secondary shock waves. As a result, the relative duration of the hot phase in the cycle remains longer, and thus the variability in the X-rays is less pronounced than for the reference case. To be detected, it would require a better sensitivity of the photometric measurements.

In the case of an initially steady atmosphere at radiative equilibrium, the coupling between the radiation and the hydrodynamics leads to:

a radiative feedback (heating) of the atmosphere, which successively expands and retracts, inducing in particular an unburial of the column, which is favorable to the lateral escape of the X-ray emitted from the hot slab; a chaotic radiative collapse, with an impact on the time variation of the X-ray flux (Figure 6); a radiative pre-heating of the incoming flow, over the length of the simulation box. Moreover the inclusion of ionization in the energy balance leads to important effects in the post shock temperature that modify the cooling efficiency and therefore the cycle duration.

In this hybrid case, we computed at LTE the radiative intensity emerging from the location of the reverse shock (Figure 15) as also from the outer boundary (Figure 16). The flux is char-acterized by -1-a huge number of atomic lines in the X rays, -2-a near blackbody profile in the Visible, with the presence of emission and absorption lines, -3-a EUV component which is very strong at the position of the reverse shock and disappears at the outer boundary, due to the importance of the absorption.

This study could be completed with a more complete simulation that would include both a dynamically-heated chromosphere and the hybrid setup. However, it appears at this stage more important to take into account a NLTE radiative description based on adapted opacities and radiative power losses. Another necessary improvement will be through a multi-group radiation transfer to catch at least the effect of EUV absorption and X-ray radiative losses on the structure of the column, and to analyse the possibility of a radiative precursor which could pre-heat the incoming flow. The study is also to be extended to multi-dimensional simulations in order to check the effects of both radiation and magnetic field closer to the real picture (Orlando et al. 2010[START_REF] Orlando | [END_REF]Matsakos et al. 2013Matsakos et al. , 2014)). Preliminary analysis of SYNSPEC and Ferguson tables revealed local aberrations, especially looking at the temperature or density derivatives (see Figure A.3). We may use the merging process to smooth most aberrations.

Fig. 1 .

 1 Fig.1. Sketch of the basis of an accretion column and its three distinctive zones: the chromosphere (left, dark grey), the accretion flow (right, mid-grey) and the zone in between (middle, light grey) hereafter called hot slab or post-shock medium.

Fig. 2 .

 2 Fig. 2. Planck (κ P , left) and Rosseland (κ R , right) opacities with respect to gas density and temperature, in log scale (cf. Appendix A). The black curve represents typical conditions met with chromosphere, accretion shock and flow.

Fig. 4 .

 4 Fig. 4. "W-Λ" simulation setup and boundary conditions.

Fig. 5 .

 5 Fig.5. Lin-log (top) and log-log (bottom) snapshots of the density (green), temperature (red) and velocity (grey) profiles of a QPO cycle (square brackets) with the "W-Λ" setup; the accreted gas falls from the right to the left (adapted from deSá et al. 2014). From left to right: beginning of a new cycle (t = 2750 s), growth of a hot slab of shocked material (t = 2884 s), quasi-isochoric cooling at the slab basis (thermal instability, t = 2994 s), collapse of the post-shock structure (falling back of the reverse shock, t = 3110 s) and end of the collapse (t = 3156 s).

Figure 8 Fig. 7 .

 87 Figure8shows seven snapshots of density and temperature profiles during the first QPO cycle (1-354 s). They are followed in the second line by 5 snapshots of the second QPO cycle (354-415 s). The second cycle differs from the first one only during the slab building (354-397 s). The sixth snapshot (415 s) is

LFig. 8 .Fig. 9 .

 89 Fig.8. Snapshots of the density (green), temperature (red) and gas velocity (grey) profiles of the first QPO cycle with the "Chr-Λ" setup; the accreted gas falls from the right to the left. The first line (between 1 and 353 s) corresponds to the first cycle. The second and third lines correspond to the beginning of the second cycle. Snapshots at t = 71 s and 415 s are very close: from this time, the cycle behaves like the previous one. A typical sequence is: growth of a hot slab of shocked material (t = 21 s), quasi-isochoric cooling at the slab basis (thermal instability, t = 168 start of the collapse of the post-shock structure (t = 336 s), impact of the collapsing material on the chromosphere (t = 354 s), launch of a new shock before the end of the collapse (t = 358 s), passing of the two shocks (t = 380 s), end of the collapse of the "old" structure = 386 s) and growth of the new slab (t = 415 s).

Fig. 10 .

 10 Fig. 10. Electron ionisation rate (ξ =

F

  Fig. 11. Snapshots of the mass density (green), gas temperature (red), escape ζ (dark blue) and absorption probabilities 1 -ζ (cyan, see Section 2.2.3.3), velocity (grey), electron density (light green), radiation energy density (magenta) and flux (orange) profiles of the first QPO cycle with the "Hybrid" setup. The accreted gas falls from the right to the left on an equilibrium atmosphere.

Fig. 12 .

 12 Fig. 12. Position of the reverse shock (r rev , purple) and the corresponding velocity (v rev , grey) together with the accreted plasma temperature before the hot slab (T acc , dark gold) and the upper chromosphere temperature (T chr , red) for the first QPO cycle in the Hybrid case.

Fig. 13 .

 13 Fig. 13. Snapshots of the temperature (red) and pressure (blue) at 9 (top) and 70 s (bottom) for the Hybrid case.

Fig. 14 .Fig. 15 .

 1415 Fig.14. Snapshot of the density (green) and temperature (red) profiles at t = 70 s with the "Hybrid" setup, post-processed hereafter.

Fig. 16 .

 16 Fig.16. Specific intensity I λ emerging from the reverse shock front during the QPO cycle of the Hybrid model (at t = 70 s, see Figure14).

Fig. 17 .

 17 Fig. 17. Time-variation of the 2-27 Å (LTE) integrated X-ray outgoing intensity (I X , blue) and of the optically thin post-shock emission (I Λ red) during a QPO cycle for the Hybrid setup.

  Fig. A.3. ∂ lT lκ R from SYNSPEC table (top) and ∂ lρ lκ P from Ferguson table (bottom) with respect to the logarithm of gas density and temperature. The dotted show transition temperatures chosen for each table. Some anomalies are revealed, especially around 10 3.7 K and high densities: this zone is cutted during the merging process.

  for V4046 Sgr,

	arXiv:1904.09156v3 [astro-ph.SR] 28 Jun 2019
	Article number, page 1 of 18

Table 1 .

 1 Characteristics of 3 simulations used in this work and their main results. "W Λ" corresponds to our reference case.Notes. H max : maximum extension reached by the post shock medium; τ cycle : cycle duration; "Window": fixed rigid non-porous transparent interface; Λ: optically thin radiative cooling; L : one solar luminosity enters the simulation box from the inner boundary; Modified Saha:Brown (1973).

	Name	Atmos.	Chromos. heating	Radiation source terms	ionisation model	H max (×10 3 km)	τ cycle Section Fig. (s)
	W-Λ "Window"*	-	Λ*	Modified Saha*	20	400	3.2	5
	Chr-Λ	Equilibrium atmosphere	L * & acoustic heating	LTE (chromos.) & Λ* (acc. flow)	Modified Saha*	17	350	3.3	8
	Hybrid	Equilibrium atmosphere	L *	Intermediate (transition: ζ)	Time-dependent collisional radiative	9	160	3.4	11

Table 2 .

 2 Comparison between our reference case ("W-Λ") and results obtained bySacco et al. (2008).

	Parameters & quantities	Sacco et al. (2008)	"W-Λ"
	Object	MP Muscae	Sun
	Atmosphere Heating function	"Window"
	Radiation	Λ	Λ
	Ionisation	Modified Saha	Modified Saha
	ρ acc (g/cm 3 )	10 -13	10 -13
	u acc (km/s)	450	400
	T acc (K) τ cycle (s)	10 3 400	3 × 10 3 400
	H max (Mm)	18	20

  3. ∂ lT lκ R from SYNSPEC table (top) and ∂ lρ lκ P from Ferguson table (bottom) with respect to the logarithm of gas density and temperature. The dotted show transition temperatures chosen for each table. Some anomalies are revealed, especially around 10 3.7 K and high densities: this zone is cutted during the merging process.

A sink is algebraically identified as negative source term.

All particles are assumed here to have the same kinetic temperature, i.e. T neutrals = T ions = T electrons = T .

Accreted material is expected to be depleted in heavy elements(Fitzpatrick 1996). However, this phenomenon is not included in this study. Article number, page 2 of 18 L. de Sá et al.: New insight on Young Stellar Objects accretion shocks

with a number abundance of 0.12%, and a mass (averaged over abundances) of 17 u.

In our 1D hydrodynamics simulations, we only consider the component of vector quantities collinear to the accretion column.

 7 For instance, in the case of an isotropic radiation, F r = 0 whereas the radiation energy can be important.

In the Hybrid case, the temperature of the accretion flow is radiatively heated by the chromosphere up to 5730 K, before the accretion process starts.

This expression is derived from the flux radiated outwards by an optically thin medium containing the radiation energy density E out r .

i.e. the temperature at which the thermal instability is triggered (∼ 8 × 10 5 K) as expected from the optically thin radiative cooling variations with respect to temperature, see Section 2.2.3.2 and references therein for further details.

Although the column plasma is expected to be at coronal regime, LTE radiation transfer is needed to build the chromosphere layer. It is therefore essential to allow radiation to escape from the first zone through the second (optically thin) one.

| | | | | | | | | | | | | | | | | | |

at ρ = 10 -13 g cm -3 and T ∼ 5000-8000 K. Article number, page 12 of 18 L. de Sá et al.: New insight on Young Stellar Objects accretion shocks

We note that the Ferguson tables showed opacity discontinuities in their hottest and densest part, as SYNSPEC tables in their coolest and densest part, as SYNSPEC tables in their coolest and densest part ( * ) (seeFigure A.3). Since the values within the transition region are ignored, this is used to artificially remove anomalies: as far as possible, the transition region must be chosen so that it covers most of them. ( * ) Few anomalies remains in regions that are not explored in our simulations (seeFigure A.4); this problem is postponed for

now.14 Fritsch & Butland (1984) derivatives are used in these papers; they generalise van Leer slopes to non-regular grids.

Article number, page 16 of 18 L. de Sá et al.: New insight on Young Stellar Objects accretion shocks

Acknowledgements. We express our gratitude to Jason Ferguson for providing us with the molecular LTE opacity tables used in this work. We thank Rafael Rodriguez for sharing with us preliminary results about NLTE microscopic collisional-radiative data, Salvatore Orlando for reading the manuscript, Ziane Izri for its contribution at the beginning of this project and Christophe Sauty for helpful discussions. I.H. thanks the Physics Department of Sorbonne Université for his visiting professorship. This work was supported by the french ANR StarShock and LabEx Plas@Par projects (resp. ANR-08-BLAN-0263-07 and ANR-11-IDEX-0004-02), PICS 6838, Programme National de Physique Stellaire of CNRS/INSU and Observatoire de Paris.

is arguably due to stronger variations of opacities (especially Planck opacity) with respect to temperature.

The specificity of the accretion shock study led us to work on dedicated opacity tables. We expose in the appendix the reasons behind this choice and the creation process. The resulting opacity table is accessible upon request. To cover the density and temperature range corresponding to our conditions, we implement in the code SYNSPEC (see Section 2.3.3), initially dedicated to stellar atmospheres, modules allowing to generate LTE monochormatic opacities at a given density and temperature. These monochromatic opacities were then averaged with the proper weighting functions to generate the adequate Rosseland and Planck mean opacities tables (hereafter called "SYNSPEC tables", see Figure A.2, top panels). These opacities are consistent with Opacity Project (see e.g. Opacity Project Team 1995) data, that we use as reference, for T between 10 3.5 K and 10 7.5 K. The advantage of SYNSPEC comes from the high number of atomic species considered, since a very detailed chemical composition is necessary to model the radiation properties of a plasma at high temperatures.

However, below 10 3.5 K, the molecular chemistry cannot be neglected, but is not included in this work on SYNSPEC. We completed thus the SYNSPEC tables with low temperature molecular opacities provided by Ferguson et al. (2005) between 10 3 K and 10 4 K ("Ferguson tables", see Figure A.2, bottom panels), that show excellent agreement with Opacity Project at upper temperatures. To facilitate the merging process, we obtained from the authors tables with compatible density and temperature grid (Ferguson, priv. comm.): mesh points from Ferguson and SYNSPEC tables are identical in the common domain (10 3.5 -10 4 K and 10 -14 -10 -6 g cm -3 ).

Appendix A.3: Preliminary study

Appendix A.3.1: Analysis of primary tables Considering opacity variations as well as temperature and density ranges, we decided to work with the logarithm of all these quantities. As first derivatives, we use then:

Appendix A.3.2: Physical and numerical constraints

In order to get a satisfying merging, several numerical and physical constraints must be respected:

as far as possible, opacities must be of class C 1 (values and first derivatives must be continuous); the transition region should be as narrow as possible; the transition region must encompass anomalies encountered in both primary tables. Such a table is composed of a limited number of discrete points: the first constraint can be reported to the interpolation method as far as opacity values in the transition present smooth variations.

To ensure a smooth transition between the molecular and the atomic (primary) tables, the transition must not take into consideration the values within the transition. The transition values loose then any physical meaning, and must be as few as possible 13 .

Appendix A.4.1: Method

In this Section, the index "A" refers to values taken at the lower transition temperature, as the index "B" for the upper ones. The transition temperatures chosen to merge SYNSPEC and Ferguson tables are:

-T A = 10 3.71 K and T B = 10 3.80 K for κ P (∼1200 K wide); -T A = 10 3.65 K and T B = 10 3.86 K for κ R (∼2800 K wide).

The problem is decoupled in temperature and in density. First, we consider the merging at each mesh density as an isolated problem, and apply a correction -if needed -to improve smoothness along the density. piecewise cubic Hermite polynomials, which ensure continuity of values (κ A , κ B ) and derivatives (∂ lT lκ A , ∂ lT lκ B ) at each transition limit; van Leer (1973) slopes to compute ∂ lT lκ A and ∂ lT lκ B , so as to prevent the apparition of spurious extrema in forcing their location to the estimated closest mesh point.

For each grid density ρ j , opacity at temperature T i ∈ [T A , T B ] is estimated using the formula:

with h = log 10 T B -log 10 T A and u i = (log 10 T i -log 10 T A )/h. This expression can be rewritten as a 3 rd degree polynomial in u i .

Appendix A.4.3: Density correction At this stage, we reached class C 1 along temperature, but there is no guarantee of continuity along density. However, in practice, it was C 1 , except for few mesh temperatures T i * .

Since the dependency in density is held by the 3 rd degree polynomial coefficients, we look at the behaviour of each of them with respect to density. Every coefficient showed spurious variations nowhere but at densities ρ j * . We apply then piecewise cubic Hermite polynomials along with van Leer slopes (density derivatives) to estimate these coefficients for each ρ j * . These new coefficients are then used to reestimate opacity values along the temperature for the ρ j * .

Appendix A.4.4: Final tables -interpolation process We checked smoothness of the result by looking at the first derivatives. The interpolation process is copied from the merging method, i.e. piecewise cubic Hermite polynomials along with van Leer slopes, since it satisfies criteria described in Section A.3.2. Interpolation is first performed along temperature at the 2 × 2 grid densities framing the requested density, so as to calculate van Leer slopes at the requested temperature and interpolate along density. Interpolating along temperature and then density showed to be slightly more accurate than interpolation along density first. This

Appendix B: Chromospheric model

One of our objectives is to describe the dynamics of the column and its impact on the chromosphere, as well as the feedback of the chromosphere on the column. This requires then to include an adequate description of the physical mechanism leading to the chromospheric heating. This appendix presents the simple but self-consistent model of a chromosphere used in this work.

Appendix B.1: Motivations and limits

The study of the solar chromosphere is a tough problem in itself. Its modelling is of interest for us since the base of the accretion column lies in the stellar chromosphere: the dynamics and observability of the column base may then depend on its structure and dynamics. Moreover, the chromosphere may be heated locally by the accretion process. The inner heating mechanism in the chromosphere is still subject of debates: it is mainly thought to originate either from acoustic waves dissipation (Biermann 1946;Schwarzschild 1948; or more recently Sobotka et al. 2016) or from MHD waves dissipation (Alfvén & Lindblad 1947;Jess et al. 2015).

Most accretion simulations model the stellar atmospherewhen it is modelled -as a hydrostatic plasma layer "tuned up" with ad-hoc sources to recover both temperature and pressure profiles (see e.g. the heating function empirically introduced by Peres et al. 1982). Although this must work for a static structure, it is delicate to predict the dynamic behaviour of such a structure facing the continuous perturbation from an infalling plasma flow: such solution is not adapted to studies involving (in a selfconsistent way) the dynamics of a perturbed atmosphere, like in the context of accretion.

We do not pretend to develop a "state of the art" model in this paper: we only aim at using a reasonable model that is both dynamic and self-consistent with our radiation hydrodynamics model. In our 1D model, we do not consider any magnetic effect but a very effective confinement of the accretion flow along the field lines. To allow fast qualitative comparison between our model and theoretical models & observations (see Figure B.1), we only used solar parameters (i.e. abundances, luminosity, mass and radius).

Appendix B.2: Acoustic waves and shocks

Acoustic waves are generated by photospheric granulation (see e.g. [START_REF] Judge | Solar MHD Theory and Observations: A High Spatial Resolution Perspective[END_REF]). These waves propagate upwards up to the height where their velocity overcome the local sound speed, and degenerate then into shocks. The nature of this mechanism is random: two different locations at the stellar surface will be crossed over by acoustic shocks that ought to be out of phase one with each other.

In our simulations, acoustic energy is supplied in the form of a monochromatic sinusoidal motion of the first Lagrangian interface (T = 60 s and f acc = 10 8 erg cm -2 s -1 , see e.g. [START_REF] Rammacher | [END_REF]Ulmschneider et al. 2005;Kalkofen 2007). Resulting acoustic waves propagate and degenerate into shocks. Figure B.1 shows several temperature snapshots of such simulation along with the chromospheric model from [START_REF] Vernazza | [END_REF]. Below 300 km, acoustic waves are damped and hardly appear on snapshots. Above 500 km, waves are fully degenerated into shocks: their strength is then governed by the balance between steepening in the pressure gradient and dissipation. Since the corona and the upper chromosphere (above 10 3 km) are readily crushed by the accretion flow, the heating of these areas is not considered in our model. Appendix B.3: From solar to stellar chromosphere

Observations of the solar chromosphere provide time and space averages of thermodynamics quantities (ρ, T , p, . . .). Detailed observation of CTTS chromospheres would demand higher space and time resolution than the ones permitted by current observational technologies. Most works on this field rely then on scaling laws (see e.g. [START_REF] Ayres | [END_REF]Calvet 1983) or ad hoc fittings to recover specific observational features (see e.g. Dumont et al. 1973;Cram 1979;Calvet et al. 1984;Batalha & Basri 1993).

Appendix C: Radiation source terms in the Hybrid model This work encompasses several radiation regimes, from optically thick LTE radiation transfer (Section 2.2.3.1) to optically thin coronal NLTE regime (Section 2.2.3.2). The momenta equations (Section 2.2) can handle all of them, assuming the proper radiation source terms are provided.

In the LTE case, both radiation energy and momentum source terms are well defined (Eq. ( 4)). In coronal regime, this is not the case. Gas and radiation are decoupled in such a regime. Radiation only acts then as a gas energy sink: the radiation energy source term (the gas sink) boils down to a cooling function (Kirienko 1993, see e.g.).

Computing the radiation flux is irrelevant in such regime and then no radiation momentum source term is provided. That is why we set s † M r to 0.

In the "Hybrid" setup, we aim at modelling radiative conditions that are neither LTE nor coronal regimes but something in between. To determine if the situation is closer to one or the other, and how close, we choose to look at the probability for a photon to escape the accretion column (see Eq. ( 7)). We use it as a weighting factor to average the source terms, as shown in Section 2.2.3.3. The process is straightforward for the radiation energy source term, but not for the radiation momentum source term since s † M r remains unknown. We assume then that the coronal Rosseland mean opacity may not significantly differ from its LTE value. This intuition is reinforced by preliminary calculations concerning NLTE radiative collisional opacities (Pérez, priv. com.).