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

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 the assumption that radiation and matter are decoupled, and use in addition a static model for the chromosphere. Aims. We test the validity of these two assumptions in refining the physics included in existing 1D models, and we propose guides for further improvement. 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 ionization 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 dynamical heating of the accretion flow and the chromosphere, leading to a possible unburrial of the whole column. These two last conclusions are in agreement with the computed monochromatic flux. Both effects (acoustic heating and radiation coupling) have an influence on the amplitude and temporal variations of the net X ray luminosity.


Accretion onto CTTSs
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 † deceased. 1This shock disappears rapidly during the development of the reverse shock, and forms back during the collapse phase of QPOs (one may compare Figure 9 at t = 97 s and t = 140 s for instance).Most of the time, according to simulations, the accretion shock is replaced by a contact discontinuity.The remaining shock, i.e. the reverse shock, is sometimes confusingly called accretion shock in the literature.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.A&A proofs: manuscript no.DE_SA-arXiv

X-ray observations and simulations
The main signature of 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;Stelzer & Schmitt 2004 for TW Hya, Schmitt et al. 2005 for BP Tau, Günther et al. 2006for V4046 Sgr, Argiroffi et al. 2007, 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 (Muzerolle et al. 1998).
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 of the stream 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.Drake et al. (2009) studied thoroughly soft X-ray emission from TW Hydrae and found no periodicity in the range 0.0001-6.811Hz.Günther et al. (2010) completed this study with optical and UV emission, and they come 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 (Orlando et al. 2013), 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 (Sacco et al. 2010;Bonito et al. 2014;Colombo et al. 2016;Costa et al. 2017).The global QPO behavior is therefore suppressed and the measured X-ray emission of the post-shock structure is then systematically low compared to the one derived from simulations (Curran et al. 2011).

Radiation effects
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, and 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 nonphysical 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.In a second stage, we analyze the role of radiation absorption/emission on the column, and in the emerging spectra.In Section 2, we present the radiation hydrodynamics model and the numeric tool we used 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.

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 characterized by a (volumetric mass) density ρ, a velocity u, a thermal pressure p and a volumetric internal energy density e.The plasma evolution is modeled by solving the hydrodynamics equations, written in the conservative form: The gas source terms2 (s e and s m ) include the contributions of thermal conduction (q C , Spitzer & Härm 1953;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 number density of free particles (neutrals, electrons and ions, see Section 2.1.2),and T represents their kinetic temperature 3 .The effect of collisional ionization and 3-body recombination processes (see Section 2.1.2) on the gas energy density is included through their corresponding rates in the thermochemistry term q χ .

Collisional-radiative ionization
The accretion shock forms where the ram pressure is balanced by the local thermal pressure: the stellar chromosphere, that needs to be modeled.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 from the solar case: therefore, we have chosen to use solar parameters in our simulations 4 , and the chemical composition (solar abundances) is then taken from Grevesse & Sauval (1998).
To estimate the total free electron density n e , 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 .For the evaluation of n e , most simulations are performed using a modified Saha equilibrium model (Brown 1973); the last simulation presented in this paper (referred to as the Hybrid setup) uses a time-dependent collisional-radiative ionization model with: collisional ionization rates given by Voronov (1997); radiative recombination rates computed by Verner & Ferland (1996); helium dielectronic recombination rate proposed by Hui & Gnedin (1997); photoionization rates derived from Spitzer (1998) and Yan et al. (1998) cross-sections, and the local radiation energy density.Ions and neutral number volumetric densities n are then computed by a conservative set of equations (see e.g. ( 1)).The electron number volumetric density is then derived from the neutrality conservation: n e = n H ii + n He ii + 2 n He iii .This calculation is performed independently from the opacity computation (see Appendix A), that uses a more refined version of the chemical composition (Grevesse & Sauval 1998).

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/emission coefficients κ ν (also called monochromatic opacity) 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 (Oxenius 1986).
The main issue in performing such calculations is an intricate coupling between the kinetic equilibrium equations (which would be easy to solve had the radiation field is known), and the radiative transfer equation (which would be easy if the atomic level populations, and hence the absorption and emission coefficients, are known).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 was satisfactorily solved in the case of stationary stellar atmospheres (see, e.g.Hubeny & Mihalas 2014), using efficient iterative methods.However, this remains difficult in the case of a non-stationary plasma, where the equations of hydrodynamic 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 flux, energy 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 band 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 take part in the computation of the spectrum emerging from this structure.This is usually done by a post processing of the hydrodynamic results by more detailed spectral synthesis tools, as detailed in Section 2.3.3.

Momenta equations
The radiation field is described here by the momenta equations (see e.g.Mihalas & Mihalas 1984) for the frequency-integrated radiation energy (E r ) and momentum or flux (M r = F r /c 2 ) volumetric densities, written in the comoving frame (Lowrie et al. 2001): 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: 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.
Depending on the expression of the radiation source terms, these equations can continuously model optically thin to thick propagation media (see e.g.Mihalas & Mihalas 1984).

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 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 radiation flux varies strongly from the X rays to the infrared.We 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 the reverse process.In LTE and regimes close to LTE, the radiation energy and momentum source terms are defined by (see e.g.Mihalas & Mihalas 1984): Two radiation-matter coupling factors appear here (in cm 2 g −1 ).The Planck mean opacity κ P is based on the frequency-integrated absorption and emission coefficients κ ν 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 (Mihalas & Mihalas 1984): In these frequency averages, the Planck mean is dominated by the contribution of strong (lines) absorption 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 (Hubeny & Mihalas 2014).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 ).The radiation source/sink terms become then: 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).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.Lequeux 2005, equation 3.66): ρ and κ R values are taken at the photon emission position.The characteristic length L c is here taken as the accretion column mean radius (i.e.1000 km, see Section 3.1).Radiation source terms become then (see Appendix C for further details): the star ( * ) and dagger ( †) denoting respectively the LTE (Eq.( 4)) and the coronal (Eq.( 6)) expressions.

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

SYNSPEC -a spectrum synthesizer
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).

Strategy and common parameters
We 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 K6 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 luminosity gets through the domain from the inner to the outer boundaries.

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 ionization 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 atmosphere is modeled 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 A&A proofs: manuscript no.DE_SA-arXiv   n e (cm −3 ) 10 11 -10 12 10 11 -10 11.5 T max (K) 10 6.5 10 6.5 to the coronal regime.Below a threshold temperature 7 , 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. MP Muscae parameters & "window" vs. chromospheric heating function), the results are in good agreement with each other. 7i.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.

X-ray luminosity
Although this case is relatively academic, it can be interesting to compute the instantaneous X-ray luminosity F Λ and its time average FΛ , to compare them with the values obtained in the two more realistic cases discussed in the next sections: 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, we compare F Λ and FΛ to the incoming mechanical energy flux: where the gravitation energy potential origin is set at the mean accretion shock position (r m 10 3 km), and the velocity is counted positively.
Figure 6 shows the time variation of F Λ .As expected, this quantity increases during the propagation of the reverse shock and decreases during the collapse.The time averaged flux FΛ is equal to 1.5 × 10 9 erg cm −2 s −1 , i.e. 36% of the incoming mechanical energy flux F M .

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 transparent8 Lagrangian interface.The outer zone is described as before, i.e. with modified Saha ionization and optically thin radiative cooling (coronal regime).However, the inner zone is now described by our chromospheric model (see Appendix B).Ionization is still described by the modified Saha equation, but we use the LTE radiation source terms as given in Equation (4).To get a dynamically heated chromosphere, we first compute a hydrostatic equilibrium, with the outer Fig. 6.Time variation of the radiative flux F Λ emerging from the hot slab for the reference (blue), the dynamical chromosphere (green) and the hybrid (red) cases.This quantity is computed assuming an optically thin coronal plasma.To allow comparison, the time is reported in reduced units of t/τ cycle and the radiative flux is normalized 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.
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 lead to small perturbations in the post-shock density (e.g. 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 thick- ness 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 ("new forward") reverse shock which propagates backwards inside the slab.This behavior 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 is therefore overlapped by 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 is slightly reduced (from 400 s for the W-Λ model to 350 s here) when using solar parameters.Since CTTSs' atmospheres have a stronger activity than the Sun's (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.However, 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 their X-ray luminosity and the overall amplitude of the variations (contrast) is reduced compared the the reference case.QPO observations may thus require both higher time resolution and improved sensitivity.
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.

Setup
In this Section, the plasma model includes collisional-radiative ionization (see Section 2.1.2).The radiation-matter coupling is described within the intermediate regime (see Section 2.2.3.3).The goal of this last setup (see Figure 10) is to inspect the net effect of the matter-radiation coupling.We have therefore chosen not to consider any chromospheric activity (i.e. the acoustic heating is kept switched off).

Ionization model
In addition, we test in this setup the effect of the time-dependent ionization through radiative ionization/recombination and collisional ionization with a time-dependent formulation (see Section 2.1.2for more details).The use of collisional and radiative rates with a calculation of the electron density at equilibrium, i.e. the equilibrium value of n e at each time step, brings differences in the transition between the (almost) neutral medium and the fully ionized 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 in Figure 9 are based on this equilibrium calculation of n e .
The main difference brought by a time-dependent calculation of the electron density is a tiny ionization delay behind the reverse shock front, as shown in Figure 11.At the shock front, the kinetic energy is converted into thermal energy, and then a part of this thermal energy is used to ionize the post-shock material with a time scale connected to the ionization rates; the affected gas layer is up to 0.2 km thick, and thus negligible compared to the whole structure (that is at least 4 km thick, see Table 1).This justifies the use of the (time-independent) Saha-Brown model for ionization for the two other cases (W-Λ and Chr-Λ).Günther et al. (2007) and Sacco et al. (2008) obtain the same conclusion from different approaches.

Atmospheric heating and cycle reduction
The first cycle is presented in Figure 9; it shows the time variations over 160 s of the gas temperature and mass density, and of the photon escape (ζ) and absorption (1 − ζ) probabilities (see Section 2.2.3.3) for the same snapshots.Although 1 − ζ shows strong variations, its net value beyond the accretion shock remains small compared to unity, and the post-shock dynamics is roughly the same as in the test case (W-Λ, see Section 3.2).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.2.
The global behavior follows the trends of the two previous models.However, several effects must be highlighted: a global motion of the whole post-shock structure, the reduction of the oscillations period and of the post-shock extension.These effects are discussed below.

QPO reduction
Changing the radiation feedback on the plasma directly affects its cooling efficiency, and thus the cycle duration along with its maximum spatial extension.This simulation shows that a weak coupling between radiation and matter 9 is able to significantly change the dynamics of QPOs: the cycle is reduced to 150 s and at most half the maximum extension of the previous simulations.

Chromospheric beating
Another effect is the heating of the upper chromosphere by the radiating post-shock plasma at the basis of the accretion column.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 12).As a consequence, the whole post-shock structure is pushed upwards from 875 km to 3150 km, thus out of the unperturbed chromosphere (by about 2000 km, see e.g.Vernazza et al. 1973).At the end of the cycle, the chromosphere is not heated any more and the slab buries back into the atmosphere.The expected behavior 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.Such effect, overestimated in 1D, has to be treated in 3D MHD as the topology and the strength of the magnetic field is known to influence the basis of the column.However the burial/unburial of the column into the chromosphere will certainly have observational consequences in terms of X-ray signatures.We note for similar reasons, that the pre-shock is also heated up to ∼ 8500K (at t =93 s).We believe that such preheating of the accretion flow, which was already pointed out in Calvet & Gullbring (1998); Costa et al. (2017), is attributed to the radiating hot slab.
These effects are quantified in Figure 13, which reports the time variations of the chromosphere/slab interface position and velocity, as also hot slab/accretion flow interface and velocity.The temperatures of the heated chromosphere and pre-shock are also reported.

Monochromatic radiative flux
As this simulation is performed using only one group of radiation frequencies, it is interesting to analyze more precisely the details of the previous radiative heating via its feedback on the monochromatic radiation flux.
To this purpose, the hydro structures has been post-processed with the SYNSPEC code (Section 2.3.3).To be consistent, we take the atomic data already provided for the calculation of the average opacities used previously (see Section 2.2.3 and Appendix A).We thus estimate the Eddington flux H λ , i.e. the first moment of the specific intensity .Since the line profile behavior is not investigated here, velocity effects are neglected.
It is important to recall that a quantitative comparison of this synthetic flux with observations, especially in the X rays (see e.g.Güdel et al. 2007;Robrade & Schmitt 2007;Drake et al. 2009) would require NLTE and 3D radiative transfer post-processing.However, 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 (Brickhouse et al. 2010;Ingleby et al. 2013;Calvet & Gullbring 1998).
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 Fig- ures 9 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 flux that emerges from this layer presents three characteristic spectral bands: in the range 1-100 Å (X-rays), the bump is attributed to the hot plasma of the column, with intense lines up to 10 12 erg/cm 2 /s/Å; in the range 100-900 Å (EUV), radiation is efficiently absorbed by the inflow and will lead to its heating as observed in our hydro simulation; in the range 900-10 000 Å (UV+Vis+IR), the second bump is attributed to the stellar chromosphere and photosphere.The strong absorption of the EUV radiation is due to the huge optical depth of the accretion flow 10 .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 flux emerging right after the reverse shock front, at r = 8.3 × 10 3 km.This figure shows that this absorption also af-10 at ρ = 10 −13 g cm −3 and T ∼ 5 × 10 3 -8 × 10 3 K. fects, 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).If these synthetic spectra allow to explain the general trends of the radiation per frequency intervals, as already mentioned, they can't be use for quantitative comparison with observations.Let us for instance consider the wavelength-integrated physical flux in the X-rays, F X in the interval 2-27 Å, observed by Chandra (Brickhouse et al. 2010)), and its time average over the cycle FX : The time variation of F X (see Figure 17, blue curve) reveals four main phases: the building of the post-shock slab (0 − 5 s), the slow radiative cooling (5 − 70 s), the collapse (70 − 90 s) and the chaotic collapse (90 − 160 s).The peak at 90 s is due to the apparition of numerous strong lines between 10 Å and 15 Å.During the chaotic collapse, the accretion shock appears and disappears several times, as well as for the reverse shock, but not in phase: F X presents therefore strong time variations, as shown in Figure 17.
The time averages of F X and F Λ (Eq.( 9)) are expected to remain smaller than the incoming mechanical energy flux F M (Eq.( 10)).This condition is satisfied for FΛ ( 30% of F M ) 11 .However, FX is about 30 times higher than this limit, which is not physical.This discrepancy is imputable to the LTE approximation used in SYNSPEC: the blue curve in Figure 17 is up to 2.5 decades higher than the red curve.
From the mechanical energy flux, and assuming a stellar radius equal to 1 solar radius, we derive from the luminosity measured by Brickhouse et al. (2010) (1.3 × 10 30 erg s −1 ) a filling factor of 2.1%.To recover this luminosity from FX , the filling factor would be even smaller.Using the optically thin approximation ( FΛ ), we calculate a more realistic value of the filling factor (6.9%).

11
FΛ should be slightly overestimated due to the larger wavelength interval used for its estimation.

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 ionization (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 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 post-shock 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 behavior 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 flowis 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 ionization 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.And 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 modeled 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 proved to have a non-negligible impact on the typical characteristics of the accretion dynamics and on its observability.This study could be completed with a 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(Orlando et al. , 2013;;Matsakos et al. 2013Matsakos et al. , 2014)).

Fig. 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 .
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 accretion shock and flow.

LFig. 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 .
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

Fig. 8 .Fig. 9 .
Fig. 8. Snapshots of the density temperature (red) 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 0 and 354 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 (thermal instability, t = 168 s), start of the collapse of the post-shock structure (t = 336 s), impact of the collapsing material on the chromosphere (t = 354 s), launch Photosphere

Fig. 11 .
Fig. 11.Electron density (dark green) and temperature (red) profiles zoomed on the reverse shock front in the early QPO cycle.

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

Fig. 13 .
Fig. 13.Position of the reverse shock (r rev , purple) and the corresponding velocity (u 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. 16 .
Fig. 16.Radiative Eddington flux H λ emerging from the reverse shock front during the QPO cycle of the Hybrid model (see Figure 14).

Fig. 17 .
Fig. 17.Time-variation of the 2-27 Å (LTE) integrated physical flux (blue) and of the optically thin post-shock emission (red) during a QPO cycle for the Hybrid setup.The dashed line represents the incoming kinetic energy flux.

Table 1 .
Brown (1973)ics 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).