New insight on Young Stellar Objects accretion shocks -- a claim for NLTE opacities
Lionel de Sá, Jean-Pierre Chièze, Chantal Stehlé, Ivan Hubeny, Thierry Lanz, Véronique Cayatte
AAstronomy & Astrophysics manuscript no. DE_SA-arXiv c (cid:13)
ESO 2019August 13, 2019
New insight on Young Stellar Objects accretion shocks – a claim for NLTE opacities –
L. de Sá , , J.-P. Chièze † , , C. Stehlé , I. Hubeny , T. Lanz , and V. Cayatte LERMA, Sorbonne Universités, Observatoire de Paris, PSL Research University, CNRS, F-75252, Paris, Francee-mail: [email protected] CEA / IRFU / SAp, CEA Saclay, Orme des Merisiers, 91191 Gif-sur-Yvette, France Steward Observatory, University of Arizona, 933 North Cherry Avenue, Tucson, AZ 85721, USA Observatoire de la Côte d’Azur (OCA), 06304 Nice, France LUTh, Observatoire de Paris, PSL University, CNRS, Paris-Diderot University, Meudon, FranceStarted September 17, 2015
ABSTRACT
Context.
Accreted material onto CTTSs is expected to form a hot quasi-periodic plasma structure that radiates in X-rays. Simulationsof this phenomenon only partly match with observations. They all rely on a static model for the chromosphere model and on theassumption that radiation and matter are decoupled.
Aims.
We explore the e ff ects on the structure and on the dynamics of the accretion flow of both a shock-heated chromosphere and ofthe coupling between radiation and hydrodynamics. Methods.
We simulate accretion columns falling onto a stellar chromosphere using the 1D ALE code AstroLabE. This code solvesthe hydrodynamics equations along with the two first momenta equations for radiation transfer, with the help of a dedicated opacitytable for the coupling between matter and radiation. We derive the total electron and ions densities from collisional-radiative NLTEionisation equilibrium.
Results.
The chromospheric acoustic heating has an impact on the duration of the cycle and on the structure of the heated slab. Inaddition, the coupling between radiation and hydrodynamics leads to a heating of the accretion flow and the chromosphere, inducinga possible unburial of the whole column. These two last conclusions are in agreement with the computed monochromatic intensity.Both e ff ects (acoustic heating and radiation coupling) have an influence on the amplitude and temporal variations of the net X-rayluminosity, which varies between 30 and 94% of the incoming mechanical energy flux, depending on the model considered. Key words.
Stars: pre-main sequence – Accretion, accretion disk – Methods: numerical – Hydrodynamics – Radiative transfer –Opacity
1. Introduction
Classical T Tauri Stars (CTTSs) are solar-type pre-main se-quence stars surrounded by a thick disk composed of gas anddust (see e.g. Feigelson & Montmerle 1999). Disk material fol-lows a near-Keplerian infall down to the truncation radius , atwhich thermal and magnetic pressures balance. Free-falling ma-terial flows then from the inner disk down to the stellar surfacein magnetically confined accretion columns (Calvet & Gullbring1998). Hot spots observations (Gullbring et al. 2000) suggestfilling factors of up to 1% (Bouvier et al. 1995).Accreted gas is stopped where the flow ram pressure and thethermal pressure of the stellar chromosphere balance: a forwardshock forms and the post-shock material accumulates at the basisof the column. The hot slab of post-shock material is separatedfrom the accretion flow by a reverse shock . A typical simulatedstructure of an accretion shock can be found e.g. in Orlando et al.(2010) and is sketched in Figure 1.One of the most direct probes for the accretion processcomes from the X-rays emitted by the dense ( n e > cm − )and hot ( T e (cid:39) † deceased. The reverse shock is sometimes called accretion shock in the litera-ture. F o r w a r d s h o c k R e v e r s e s h o c k Chromosphere Hotslab Accretion flow
Fig. 1.
Sketch of the basis of an accretion column and its three dis-tinctive zones: the chromosphere ( left, dark grey ), the accretionflow ( right, mid-grey ) and the zone in between ( middle, lightgrey ) hereafter called hot slab or post-shock medium . et al. 2002 and Stelzer & Schmitt 2004 for TW Hya, Schmittet al. 2005 for BP Tau, Günther et al. 2006 for V4046 Sgr,Argiro ffi et al. 2007, 2009 for MP Muscae, Robrade & Schmitt2007 for RU Lup and Huenemoerder et al. 2007 for Hen 3-600).Another signature is the UV-optical veiling, which is attributedto the post shock medium, the heated atmosphere and the Article number, page 1 of 18 a r X i v : . [ a s t r o - ph . S R ] A ug & A proofs: manuscript no. DE_SA-arXiv pre-shock medium (Calvet & Gullbring 1998). In addition,Doppler profiles of several emission lines trace the high velocityin the funneled flow (up to 500 km s − , according to Muzerolleet al. 1998).1D hydrodynamical models (Sacco et al. 2008, 2010) pre-dict Quasi-Periodic Oscillations (QPOs) of the post-shock slabwith periods ranging from 0.01 to 1000 s, depending on the in-flow density, metallicity, velocity and inclination with respectto the stellar surface. For a typical free-fall radial velocity of400 km s − , Sacco et al. (2010) found for instance a period of160 s at 10 cm − . These oscillations are triggered by the cool-ing instability (for further details, see e.g. Chevalier & Imamura1982; Walder & Folini 1996; Mignone 2005).Although plasma characteristics derived from X-ray obser-vations are consistent with the density and the temperature pre-dicted by these numerical studies, there is no obvious observa-tional evidence for such periodicity. Drake et al. (2009) studiedthoroughly soft X-ray emission from TW Hydrae and found noperiodicity in the range 0.0001–6.811 Hz. Günther et al. (2010)completed this study with optical and UV emission, and theycame to the same conclusion in the range 0.02–50 Hz. However,a recent photometric study of TW Hya based on MOST satel-lite observations reports possible oscillations with a period of650–1200 s, which could be assigned to post-shock plasma os-cillations (Siwak et al. 2018).Observations thus raise the question of the existence of anoscillating hot slab in the accretion context. Several numericalstudies explored multi-dimensional magnetic e ff ects, like leaksat the basis of the column (Orlando et al. 2010), the tapering ofthe magnetic field (Orlando et al. 2013), or perturbations in theflow (Matsakos et al. 2013). Although QPOs are still obtainedin these numerical studies, the accretion funnel basis is eitherfragmented in out-of-phase fibrils, or buried under a cooler anddenser gas layer that strongly absorbs X-rays. The observationof global synchronous QPOs becomes therefore very challeng-ing (Curran et al. 2011; Bonito et al. 2014; Colombo et al.2016; Costa et al. 2017). The e ff ect of the slab burial into thechromosphere has also been explored in several 1D simulations(Drake 2005; Sacco et al. 2010). Depending on the depth of theburial, the radiation may only escape the post-shock structurefrom its upper part, leading to a significant reduction of theX-ray luminosity.In these numerical works, the accretion is supposed to takeplace on a quiet medium (an isothermal atmosphere in the bestcases). Moreover, the post-shock medium is assumed to be op-tically thin, and the coupling between radiation and matter isreduced to a gas cooling function (see e.g. Kirienko 1993, re-ported in Figure 3). Although this assumption can be justifiedto model the infalling gas and the post-shock plasma, it is in-consistent with any stellar atmosphere model. The energy bal-ance between radiation and gas in the lower stellar atmosphereis then replaced by a non-physical tuning (heating function, o ff threshold, . . . ). Such an assumption may a ff ect the burial of thepost-shock structure as well as the accretion structure itself.In this work, we focus and refine the physics encompassed inexisting 1D models. We first explore the e ff ect of chromosphericshocks perturbations on the accretion dynamics. We analyse thenhow radiation may a ff ect the chromospheric, post-shock and ac-creted plasmas as well as the QPO duration and the hot slabburial; we also synthesise and discuss the accretion signature inthe emerging radiative spectra. In Section 2, we present the radi-ation hydrodynamics model and the numerical tools we use for the hydrodynamics and the spectra synthesis. We detail in Sec-tion 2.2.3 the two extreme radiative regimes encountered in thiscontext, and a simple model for intermediate radiative regimes.Section 3 is dedicated to accretion simulations and to the cor-responding discussions. The last section (Section 4) presentscaveats and possible improvements to this work.
2. Physical and numerical models
We consider a star of radius R (cid:63) and mass M (cid:63) . The accreted andstellar atmospheric plasmas at position r ( r = (cid:107) r (cid:107) ), hereaftertaken from the stellar surface , are characterised by a (volumetricmass) density ρ , a velocity u , a thermal pressure p and a volumet-ric internal energy density e . The plasma evolution is modelledby solving the hydrodynamics equations, written in the conser-vative form: ∂ t ρ + ∇ · ( ρ u ) = ∂ 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 (cid:63) ρ/ ( R (cid:63) + r ) r / r .The gas source terms ( s e and s m ) include the contributionsof thermal conduction ( q C , Spitzer & Härm 1953; Vidal et al.1995), gravity ( g ( r )), artificial viscosity ( p vis and q vis , von Neu-mann & Richtmyer 1950) and the coupling with radiation ( s M r and s E r , see Section 2.2.3). The closure relation for this systemof equations – the equation of state – is adapted from the idealgas law : p = n tot k T ⇔ e = / p , where n tot stands for the totalvolumetric number density of free particles (neutrals, electronsand ions), and T represents their kinetic temperature . The con-tribution of ionisation / recombination on the gas energy densityis included in the thermochemistry term q χ , and is discussed inthe subsequent section (2.1.2). The forward shock forms where the ram pressure is balanced bythe 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 thedynamics of this chromosphere, and in absence of any reliableinformation, our chromospheric model (see Appendix B) isinspired by the solar case: therefore, we have chosen to use solarparameters in our simulations, and the chemical composition(solar abundances ) 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 iscompleted by a "catch-all" metal "M" .Most simulations are performed using time-independent ion-isation models, for instance the modified Saha equilibrium of 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 (Fitz-patrick 1996). However, this phenomenon is not included in this study. with a number abundance of 0.12%, and a mass (averaged over abun-dances) of 17 u.Article number, page 2 of 18. de Sá et al.: New insight on Young Stellar Objects accretion shocks Brown (1973) (see e.g. Sacco et al. 2008) or a detailed collisionalionisation calculation (e.g. Günther et al. 2007). To estimate thetotal free electron density n e in the two first setups, we use themodified Saha model (for which q χ = Hybrid setup) uses a time-dependent collisional-radiative ioni-sation 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 Spitzer (1998) andYan et al. (1998) cross-sections, and the local radiation en-ergy 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 derivedfrom the neutrality conservation: n e = n H ii + n He ii + n He iii . Fi-nally, the thermochemistry term q χ sums all these contributions,weighted by the corresponding gained / lost energy.These calculations are performed independently from theopacity computation (see Appendix A), that uses a more refinedversion of the chemical composition (Grevesse & Sauval 1998). The coupling between radiation and matter enters at di ff erentscales in astrophysical plasmas. At a microscopic scale, radia-tion a ff ects the thermodynamical state of the matter through itscontribution to the populations of the electronic energy levelsof each plasma ion. The computation of these populations isbased on a large set of kinetic equilibrium equations that takeinto account excitation and de-excitation processes due to col-lisions (interactions with massive particles, mostly electrons) aswell as radiative processes (interactions with photons). This stepallows to derive also the monochromatic absorption and emis-sion coe ffi cients, resp. κ ν (also called monochromatic opacity ,in cm g − ) and η ν (in erg cm − s − ), which in turn are used tocompute the local radiation intensity by solving the equations ofradiative transfer. Two limiting (and simplifying) cases are ex-pected: at large electron densities, one recovers the Local Ther-modynamic Equilibrium (LTE), whereas at low density and foran optically thin medium, the coronal limit is reached (Oxenius1986).The main issue in performing such calculations is an intri-cate coupling between the kinetic equilibrium equations (eas-ily solved given the radiation field), and the radiative transferequation (simple to calculate knowing the atomic level popula-tions, and hence the absorption and emission coe ffi cients). Sincea mean free path of photons is typically much larger than themean free path of massive particles, an explicit treatment of theradiation transport necessarily involves a significant non-localityof the problem. This issue is satisfactorily solved in the caseof stationary stellar atmospheres (see, e.g. Hubeny & Mihalas2014), using e ffi cient iterative methods. However, this remainsdi ffi cult in the case of a non-stationary plasma, where the equa-tions of hydrodynamics need to be coupled, at each time, withthe equations for the radiative transfer.Therefore, the previous kinetic equations have to be solvedsimultaneously with the monochromatic radiative transfer equa- tions. This allows computing the frequency-averaged local radi-ation energy, flux and pressure, and helps including these quan-tities in the hydrodynamics equations (Eq. (1)). In practice, thisexact description would require extensive numerical resources:the di ffi culty is commonly reduced by averaging the radiationquantities by frequency bands. In the multi-groups approxima-tion, the absorption and emission coe ffi cients are averaged overseveral frequency bands using adapted weighting functions: thelarger the number of groups, the better the precision of the com-putation. The simplest and most commonly used approach is the monogroup approximation, which means that the radiation quan-tities are averaged over the whole frequency domain covered.Besides these delicate issues, radiative transfer takes part inthe computation of the spectrum emerging from this structure.This is usually done by the post-processing of the hydrodynamicresults by more detailed spectral synthesis tools, as detailed inSection 2.3.3. The radiation field is described here by the momenta equations (see e.g. Mihalas & Mihalas 1984) for the frequency-integratedradiation energy volumetric density ( E r , in erg cm − ) and mo-mentum ( M r , in erg cm − s) or flux ( F r = c M r , in erg cm − s − ),written in the comoving frame (Lowrie et al. 2001): ∂ t E r + u · ∂ t M r + c ∇ · M r + ( P r : ∇ ) · u + ∇ · ( E r u ) = s E r ∂ t M r + u · ∂ t P r / c + ∇ · P r + ( M r · ∇ ) u + ∇ ( M r · u ) = s M r (2)The (monogroup) radiation quantities are integrated from 1 to10 Å. The M1 closure relation allows then to derive the radia-tion pressure P r from the radiation energy density: P r = D E r . D and χ are respectively the Eddington tensor and factor ( D ≡ χ in1D) and are defined as follows: D = − χ I + χ − i ⊗ i , χ = + f + (cid:112) − f (3)with the reduced flux f = F r / ( c E r ) (and f = (cid:107) f (cid:107) ), the flux di-rection i = f / f = F r / F r and I the second-order identity tensor.As a drawback, the M1 radiation transfer may not properlymodel the radiation field in structures that involve more than onemain radiation source (see e.g. Jiang et al. 2014a,b; S ˛adowskiet al. 2014). Moreover, contrarily to the radiation energy, thecontribution of the radiation flux to the hydrodynamics is notstraightforward to interpret ; both are presented and discussedwith 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 thickpropagation media (see e.g. Mihalas & Mihalas 1984). This work aims at describing in a consistent way the systemcomposed of three zones, which are coupled together throughradiation but in di ff erent thermodynamical states (Figure 1):the dense and optically thick near-LTE chromosphere (Section2.2.3.1) on the one hand, the optically thin coronal hot accretionslab and cold accretion flow (Section 2.2.3.2) on the other hand. In our 1D hydrodynamics simulations, we only consider the compo-nent of vector quantities collinear to the accretion column. For instance, in the case of an isotropic radiation, F r = whereas theradiation energy can be important. Article number, page 3 of 18 & A proofs: manuscript no. DE_SA-arXiv
We also expect, according e.g. to Calvet & Gullbring (1998),that the frequency distribution of the measured radiation variesstrongly from the X-rays to the infrared. We have decided towork step by step, using a model which makes a continuoustransition between the optically thick LTE approximation andthe coronal limit, as described in Section 2.2.3.3.2.2.3.1. Optically thick limitThe deep stellar atmosphere is optically thick and can be consid-ered at LTE, i.e. each microphysics process is counter-balancedby its reverse process. In LTE and regimes close to LTE, themonochromatic absorption and emission coe ffi cients are linkedthrough the Planck distribution function: η ν = κ ν ρ c B ν . The ra-diation energy and momentum source terms are then defined by(see e.g. Mihalas & Mihalas 1984): s ∗ E r = κ P ρ c (cid:16) a R T − E r (cid:17) and s ∗ M r = − κ R ρ c M r (4)where a R is the radiation constant. Two radiation-matter cou-pling factors appear here (in cm g − ). The Planck mean opac-ity κ P is based on the frequency-integrated absorption coe ffi cient κ ν weighted by the Planck distribution function B ν , while theRosseland mean opacity κ R is the harmonic mean of κ ν weightedby the temperature derivative of the Planck function ∂ T B ν , asfollows (Mihalas & Mihalas 1984): κ P = (cid:82) κ ν B ν d ν (cid:82) B ν d ν and κ − = (cid:82) κ − ν ∂ T B ν d ν (cid:82) ∂ T B ν d ν (5)In these frequency averages, the Planck mean is dominated bystrong absorption features (typically lines), whereas the Rosse-land mean is dominated by the regions in the spectrum of low-est monochromatic opacity. As a consequence, at large opticaldepths, κ P correctly describes the energy exchange between par-ticles and photons, while κ R gives the correct total radiative flux(Hubeny & Mihalas 2014). − − − − − ρ (cid:16) g cm − (cid:17) l og T ( K ) log κ P (cid:16) cm g − (cid:17) − − − − − ρ (cid:16) g cm − (cid:17) log κ R (cid:16) cm g − (cid:17) − . − . − . − . − . − . . . . . Fig. 2.
Planck ( κ P , left ) and Rosseland ( κ R , right ) opacities with re-spect to gas density and temperature, in log scale (cf. Appendix A). The black curve represents typical conditions met with chromosphere,accretion shock and flow. Several opacity tables are available for a variety of chemicalcompositions. However, they all fail to cover the full ( ρ, T ) do-main explored in our simulations (see solid black line in Figure2). 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. 2.2.3.2. Optically thin limitDue to its very low density ( ρ (cid:39) − g cm − ), the accretedplasma can be described by the limit regime where the gas den-sity tends towards zero: the coronal regime . The coupling be-tween radiation and matter boils down in this case to an opticallythin radiative cooling function Λ ( T ) (in erg cm s − ). In Eq. (2),the radiation source / sink terms become then: s † E r = n e n H Λ ( T ) and s † M r = (6)The first quantity represents the net radiation power emitted byunit volume in all directions (4 π sr) by a hot optically thin plasma(in erg cm − s − sr − ). The term s † M r is set to zero since there isno coupling between radiation and matter in this regime (see Ap-pendix C for more details). Fig. 3.
Optically thin radiative cooling (in erg cm s − ) for di ff erentmetallicities Z, versus gas temperature (K), adapted from Kirienko(1993). The present work is based on the cooling function providedby Kirienko (1993), reproduced in Figure 3, with Z / Z (cid:12) = ζ = − exp( − τ e )3 τ e , τ e = κ P ρ L c (7) ρ and κ P values are taken at the photon emission position. Thecharacteristic length L c is here taken as the accretion columnmean radius (i.e. 1000 km, see Section 3.1). Radiation sourceterms become then (see Appendix C for further details): 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. Article number, page 4 of 18. de Sá et al.: New insight on Young Stellar Objects accretion shocks
Observations indicate that, in general, the ambient magnetic fieldis of the order of 1 kG (Johns-Krull et al. 1999; Johns-Krull2007). The resulting Larmor radius (1 mm) is very small, i.e. theplasma follows the magnetic field lines. Moreover, the Alfvénvelocity reaches 3% of the speed of light and the magnetic wavesbehave thus like usual light waves. Therefore, focusing on theheart of an accretion column in strong magnetic field case, wecan model the accreted material along one field line, that will beassumed to be radial relative to the stellar center. Since the ac-cretion process is expected to involve strong shocks, we chose anumerical tool able to achieve very high spatial resolution.
The present work is based on numerical studies performed withthe 1D Arbitrary-Lagrangian-Eulerian (ALE) code AstroLabE(see e.g. de Sá et al. 2012; Chièze et al. 2012). It is based onthe Raphson-Newton solver (Numerical Recipes, Section 9) anda 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 chemistryequations (see Sections 2.1 and 2.2), the equations describingthe behaviour of the grid points. The space discretisation canfollow an Eulerian or a Lagrangian description. Moreover, thegrid can freely adapt to hydrodynamics situations (the arbitrary description, Dorfi & Drury 1987): this helps us reach high reso-lution around shocks with fixed cardinality ( δ r / r max (cid:39) − with150–300 grid points).Beside its application to stellar accretion (de Sá 2014), Astro-LabE has been used in several astrophysical situations such asthe interstellar medium (Lesa ff re 2002; Lesa ff re et al. 2004),experimental radiative shocks (Stehlé & Chièze 2002; Bouquetet al. 2004) or type Ia supernovae (Charignon & Chièze 2013)studies. For the computation of the opacities and of the emerging spec-tra, we used the public 1D spectrum synthesis code SYNSPEC(Hubeny & Lanz 2017). It is a multi-purpose code that can ei-ther construct a detailed synthetic spectrum for a given modelatmosphere or disks, or generate LTE opacity tables. In this pa-per, we used SYNSPEC both for generating opacity tables (seeSection 2.2.3.1 and Appendix A), and for the snapshots spectrapresented in Section 3.4.4.The resulting synthetic spectrum reflects the quality of the in-put astrophysical model; using an LTE model results in an LTEspectrum, while using a NLTE model results in a NLTE spec-trum. The snapshots of our hydrodynamic simulations providetemperature and density as a function of position; it is thereforestraightforward to compute LTE spectra for such structures. Itwould be in principle possible to construct approximate NLTEspectra, keeping temperature and density fixed from the hy-drodynamic simulations (the so-called "restricted NLTE prob-lem"). This could be done for instance by the computer programTLUSTY (Hubeny & Lanz 1995, 2017), which would provideNLTE level populations that can be communicated to SYNSPECto produce detailed spectra. However, as previously mentioned,such a study is computationally very demanding and is well be-yond the scope of the present paper. Nevertheless, since NLTE e ff ects may be important, this will be done in a future paper. Itwill allow to inspect the e ff ect of the LTE approximation on ourresults.This synthetic spectrum, computed at di ff erent altitudes ofthe accretion column, will reveal the role played by the di ff erentparts of the spectrum, from X-ray to Visible (1–10 Å). However,it is important to note that, as the accretion column is limited indiameter, some e ff ects, like the absorption by the coldest partsare only pertinent for an observation along or near the directionof the accretion column. A 3D radiative transfer post-processingwould then be more suitable to the geometry of the system (Ibguiet al. 2013).
3. Accretion basis simulations
We have simulated for this study several physical situations inorder to check the net e ff ect on the QPOs of the chromosphericmodel on one side and of the matter-radiation coupling on theother side. We present first the reference case: a gas flow hits afixed, rigid and non-porous interface (W– Λ case, Section 3.2).We check then the e ff ect of a dynamically heated chromosphereon the accretion process (Chr– Λ case, Section 3.3) and wefinally check the e ff ect of the radiation feedback on matter(Hybrid case, Section 3.4). The conditions and main results ofeach simulation are resumed in Table 1.The simulations presented in this paper share few parame-ters: – the computational domain size is r out = km (the outerboundary limit); – the column / fibril radius is set to L c = × − ; – for the gravity magnitude, we use R (cid:63) = R (cid:12) and M (cid:63) = M (cid:12) ; – the accreted gas enters the computational domain through theouter boundary with ρ acc = − g cm − , T acc = and v acc =
400 km s − .The velocity of the accreted gas is derived from the free-fall ve-locity at r = r out above the stellar surface, considering a nullradial velocity at the truncation radius R tr = . (cid:12) (taken herefrom the center of the star).When the M1 radiation transfer is used (either near-LTEtransfer or intermediate regime), one solar surface luminosity( L (cid:12) = . × erg cm − s − ) enters from the inner boundary,and c × E outr / , with E outr beingthe radiation energy density of the last computational cell. Λ ) In the reference case, we simulate the accretion stream using thesame physics and assumptions than in previous models (see e.g.Sacco et al. 2008; Koldoba et al. 2008). The matter-radiation The ratio of the lateral to the longitudinal extension (in terms of typ-ical radiative mean free path) of the column should be ideally large tojustify 1D approximation for the computation of the e ff ect of the radia-tive transfer throughout the system. In the Hybrid case, the temperature of the accretion flow is radiativelyheated by the chromosphere up to 5730 K, before the accretion processstarts. This expression is derived from the flux radiated outwards by an op-tically thin medium containing the radiation energy density E outr .Article number, page 5 of 18 & A proofs: manuscript no. DE_SA-arXiv
Table 1.
Characteristics of 3 simulations used in this work and their main results. "
W(cid:21) Λ " corresponds to our reference case. Name
Atmos. Chromos. Radiation ionisation H max τ cycle Section Fig.heating source terms model ( × km) (s) W– Λ "Window"* – Λ * Modified Saha* 20 400 3.2 5Equilibriumatmosphere L (cid:12) * & LTE (chromos.)& Λ * (acc. flow) Chr– Λ acoustic Modified Saha* 17 350 3.3 8heating Hybrid
Equilibrium L (cid:12) * Intermediate Time-dependent 9 160 3.4 11atmosphere (transition: ζ ) collisional radiative Notes. H max : maximum extension reached by the post-shock medium; τ cycle : cycle duration; "Window" : fixed rigid non-porous transparent inter-face; Λ : optically thin radiative cooling; L (cid:12) : one solar luminosity enters the simulation box from the inner boundary; Modified Saha : Brown(1973). coupling is then described by the coronal radiative cooling (Sec-tion 2.2.3.2) and the plasma ionisation is computed with themodified Saha equation (Section 2.1.2). In order to simplifythe discussion, we focus on the post-shock structure and on theglobal dynamics. The stellar chromosphere is modelled in thesimplest way, hereafter called the "window" model. It consistsin a fixed rigid non-porous transparent interface. The main pa-rameters are resumed in Figure 4. W a ll ( c h r o m o s . ) Accretion flow v acc T acc ρ acc fixedrigidnon-porousModified Saha ionisationRadiation: coronal regime Fig. 4. "W– Λ " simulation setup and boundary conditions. Besides the fact that matter accumulates on the left (inner) rigidboundary interface, the system is found to be perfectly periodic.Figure 5 presents five snapshots of density, temperature and ve-locity profiles during a QPO cycle far from the initial stages. Theaccreted gas falls from right to left. A hot slab of shocked mate-rial builds first ( t = , the fast,quasi-isochoric, cooling of the slab basis causes the collapse ofthe post-shock structure ( t = t = ff erences (Sun vs.MP Muscæ parameters & "window" vs. chromospheric heatingfunction), the results are in good agreement with each other. An X-ray radiative power of 1.3 × erg s − was measured inthe range 2–27 Å by Brickhouse et al. (2010) for TW Hydræ i.e. the temperature at which the thermal instability is triggered( ∼ × K) as expected from the optically thin radiative cooling vari-ations with respect to temperature, see Section 2.2.3.2 and referencestherein for further details.
Table 2.
Comparison between our reference case ("W– Λ ") and resultsobtained by Sacco et al. (2008). Parameters
Sacco et al. "W– Λ " & quantities (2008)Object MP Muscæ Sun
Atmosphere
Heating function "Window"
Radiation
Λ Λ
Ionisation
Modified Saha Modified Saha ρ acc (g / cm ) 10 − − u acc (km / s) 450 400 T acc (K) 10 × τ cycle (s) 400 400 H max (Mm) 18 20 n e (cm − ) 10 –10 –10 T max (K) 10 and an accretion flow velocity estimated at 500 km s − . We com-pute therefore the instantaneous X-ray surface luminosity L Λ (inerg cm − s − ) and its time average ¯ L Λ : L Λ = (cid:90) slab n e n H Λ ( T ) dr & ¯ L Λ = τ cycle (cid:90) τ cycle L Λ dt (9)to compare them with the values obtained with the di ff erent mod-els presented in the subsequent sections and with the observa-tional work of Brickhouse et al.. These quantities are commonlycompared to the incoming kinetic energy flux. However, sincethe flow accelerates in its free-fall from the outer boundary downto the reverse shock, the plasma velocity and density may changebetween the outer boundary of the simulation box and the loca-tion of the reverse shock. To get round this issue, one must con-sider the mechanical energy flux. This flux is calculated at anyposition r by: F M = ρ v ( r ) + (cid:90) rr G M (cid:63) ρ v z dz (10)where the origin of the gravitational energy potential is set at themean forward shock position ( r (cid:39) km). The conservation The "slab" is here defined as the plasma at temperature above 10
K.Article number, page 6 of 18. de Sá et al.: New insight on Young Stellar Objects accretion shocks r ( × km) − − − − l og ρ (cid:16) g c m − (cid:17) r ( × km) 0 10 20 r ( × km) 0 10 20 r ( × km) 0 10 20 r ( × km) cycle3 . . . . l og T ( K ) | | |− − − − l og ρ (cid:16) g c m − (cid:17) t = | | | t = | | | t = | | | t = | | | t = . . . . l og T ( K ) r (km) − . − . − . . v (cid:16) × k m / s (cid:17) r (km) 2 3 4log r (km) 2 3 4log r (km) 2 3 4log r (km) 3 . . . . l og T ( K ) 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 de Sá et al. 2014). From left to right : beginning of a new cycle ( t = t = t = t = t = of the mechanical energy induces that F M does not dependon the position r . The value derived from our simulations is F M = . × erg cm − s − . . . . . . . t /τ cycle − . − . − . − . . . l og L Λ / F M W − Λ Chr − Λ Hybrid
Fig. 6.
Time variation of the surface luminosity L Λ emitted by the hotslab for the reference ( blue ) , the dynamical chromosphere ( green ) andthe hybrid ( red ) cases. These quantities are computed assuming an op-tically thin coronal plasma. To allow comparison, the time is reportedin reduced units of t /τ c y cle and the luminosity is normalised to the in-coming mechanical energy flux F M defined in Eq. (10). The values ofthe cycle duration for each setup is reported in Table 1. Figure 6 shows the time variation of L Λ . As expected, thisquantity increases during the propagation of the reverse shockand decreases during the collapse. The time-averaged luminosity¯ L Λ is equal to 1.5 × erg cm − s − , i.e. 36% of the incomingmechanical energy flux F M . Λ ) In this second setup (see Figure 7), we aim at studying the ef-fect of a dynamically heated chromosphere on the phenomenondescribed in the previous Section. To achieve this, we "divide" the computational domain into two zones separated by a trans-parent Lagrangian interface.The outer zone is described as before, i.e. with modified Sahaionisation and optically thin radiative cooling (coronal regime).However, the inner zone is now described by our chromosphericmodel (see Appendix B). Ionisation is still described by the mod-ified Saha equation, but we use the LTE radiation source termsas given in Eq. (4). To get a dynamically heated chromosphere,we first compute a radiative-hydrostatic equilibrium, with theouter zone inactivated, and with one solar luminosity crossingthe entire domain (no e ff ect on the outer zone). Acoustic energyis then injected in the form of monochromatic sinusoidal motionof the first interface (a "window") with a 60 s period to mimicsolar granulation. Several snapshots of temperature profiles arepresented in Figure B.1. Once the shock-heated chromospherereaches its stationary regime, the accretion process is launched(in the outer zone). P ho t o s ph e r e cE outr / v acc T acc ρ acc L (cid:12) u g = A sin(2 π t /τ )rigidnon-porous < radiation > LTE coronal regimeModified Saha ionisation ( ou t e r z on e ) A cc r e ti on fl o w ( i nn e r z on e ) C h r o m o s ph e r e Transparent Lagrangian interface
Fig. 7. "Chr– Λ " simulation setup and boundary conditions. A = . − and τ =
60 s. Although the column plasma is expected to be at coronal regime, LTEradiation transfer is needed to build the chromosphere layer. It is there-fore essential to allow radiation to escape from the first zone throughthe second (optically thin) one. Article number, page 7 of 18 & A proofs: manuscript no. DE_SA-arXiv
Figure 8 shows seven snapshots of density and temperature pro-files during the first QPO cycle (1–354 s). They are followedin the second line by 5 snapshots of the second QPO cycle(354–415 s). The second cycle di ff ers from the first one only dur-ing the slab building (354–397 s). The sixth snapshot (415 s) isvery close to the snapshot of the first cycle at t =
71 s. The (un-changed) end of the second cycle is then not reported.During the installation phase (1–336 s) of the reverse shock,the post-shock structure follows more or less the same scenariothan for the reference case (W– Λ ). After several periods of theacoustic waves, small di ff erences occur. The transmission ofthese waves / shocks to the accretion column depends on the leapof the acoustic impedance between the upper chromosphere andthe hot slab, which results in reflection / transmission of thesewaves / shocks at this interface. The smallest leap is reached atthe end of the collapse, near 336 s, leading to a transmission in-crease, which however remains still low. Their e ff ect leads tosmall perturbations in the post-shock density (as it can alreadybe seen at 168 s). Table 3.
Position of the old and new reverse shocks between t =
358 sand t =
397 s (see Figure 8).
Time (s)
358 380 386 397 r old (km) r new (km) After this time, the transmitted waves start to feed with mat-ter 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, Figure5). This structure collapses and hits at 354 s the dense chromo-sphere, leading to a secondary reverse shock which propagatesbackwards inside the slab. This behaviour is confirmed by thevelocity variations shown in grey in Figure 8. The two reverseshocks pass then each other: the positions of the new shock (orcontact discontinuity) and the previous ( old ) one are resumed inTable 3. The end of one cycle therefore overlaps the beginningof a new one.
This model implies two main observational consequences.First, compared to the reference case, the QPO cycle periodis modified by the acoustic heating. The question of possibleresonance is pointless regarding multi-mode acoustic heating byout-of-phase waves emitted in di ff erent locations. The period τ cycle is slightly reduced (from 400 s for the W– Λ model to 350 shere, Table 1) when using solar chromospheric parameters.Since CTTSs’ atmospheres have a stronger activity than theSun’s one (that we use for the chromospheric model), the e ff ectis expected to be enhanced in CTTSs.The second e ff ect deals with the X-ray luminosity variation dur-ing a cycle, as reported in green in Figure 6. The growth phase iscomparable with the W– Λ setup, but the acoustic perturbationsfrom the chromosphere induce strong di ff erences in the collapsephase. Moreover, the overlapping of the beginning and end ofthe cycles a ff ect the X-ray luminosity and the overall amplitudeof the variations (contrast) is reduced compared to the referencecase. QPO observations may thus require both higher timeresolution and improved sensitivity. The time averaged surface luminosity (Eq. (9)) is here equal to 4.0 × erg cm − s − , i.e.94% of the mechanical energy flux F M .These results show that, compared to the reference case, thedynamical heating of the chromosphere impacts the duration ofthe QPO period and its observability. Of course, a more realisticdescription of the chromospheric heating would require at leasta 2D MHD picture. For instance, we know that chromosphericperturbations may lead – inside the column – to the developmentof fibrils (see e.g. Matsakos et al. 2013, ChrFlx ff ectedby the chromospheric perturbations. In this Section, the plasma model includes collisional-radiativeionisation (see Section 2.1.2). The radiation-matter coupling isdescribed within the intermediate regime (see Section 2.2.3.3)and the outer radiation flux is set to c × E outr /
4. The goal of thislast setup (see Figure 9) is to inspect the net e ff ect of the matter-radiation coupling. We have therefore chosen not to consider anychromospheric activity. Following the preliminary process of theprevious setup (see Section 3.3.1), the outer zone is first inacti-vated, and the radiative hydrostatic equilibrium is computed inthe inner zone; once the stationary regime is reached, the ac-cretion process is launched. A key advantage of this process isthat nothing is needed to maintain the chromospheric structure,which can therefore freely evolve depending on the physical pro-cesses in play only. We have tested in this setup the e ff ect of the time-dependentionisation through radiative ionisation / recombination andcollisional ionisation with a time-dependent formulation (seeSection 2.1.2 for more details).The main di ff erence brought by a time-dependent calcula-tion of the electron density is a tiny ionisation delay behind thereverse shock front, as shown in Figure 10.At the shock front, the kinetic energy is converted into thermalenergy, and then a part of this thermal energy is used to ionisethe post-shock material with a time scale connected to theionisation rates; the a ff ected gas layer is up to 0.2 km thick,and thus negligible compared to the whole structure (that isat least 10 km thick, see Table 1). This justifies the use of atime-independent model for ionisation in the previous setups(W– Λ and Chr– Λ ). Günther et al. (2007) and Sacco et al. (2008)obtain the same conclusion from di ff erent approaches.However, compared to the Saha-Brown equilibrium calcula-tions, the use of collisional and radiative rates to derive the equi-librium electron density brings di ff erences in the transition be-tween the (almost) neutral medium and the fully ionised plasma.This transition lays between 10 and 10 K. However, suchtemperatures are only reach by the accreted gas during the cool-ing instability. Its overall e ff ect is hence negligible. The resultspresented for the Hybrid case (see Figure 11) are thus based onthis collisional-radiative equilibrium calculation of n e . Article number, page 8 of 18. de Sá et al.: New insight on Young Stellar Objects accretion shocks r (km) − − − − l og ρ (cid:16) g c m − (cid:17) t =
21 s r (km) t =
46 s r (km) t =
71 s r (km) t =
168 s r (km) t =
336 s r (km) t =
351 s r (km) t =
353 s . . . . l og T ( K ) | | |− − − − l og ρ (cid:16) g c m − (cid:17) t =
354 s | | | t =
358 s | | | t =
380 s | | | t =
386 s | | | t =
397 s | | | t =
415 s . . . . l og T ( K ) cycle2 3 4log r (km) − . − . − . . v (cid:16) × k m / s (cid:17) r (km) 2 3 4log r (km) 2 3 4log r (km) 2 3 4log r (km) 2 3 4log r (km) 3 . . . . l og T ( K ) Fig. 8.
Snapshots of the density ( green ), temperature ( red ) and gas velocity ( grey ) profiles of the first QPO cycle with the "Chr– Λ " setup; theaccreted 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 correspondto 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. Atypical sequence is: growth of a hot slab of shocked material ( t =
21 s), quasi-isochoric cooling at the slab basis (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 of a newshock 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 ( t =
386 s) andgrowth of the new slab ( t =
415 s). P ho t o s ph e r e cE outr / v acc T acc ρ acc L (cid:12) fixedrigidnon-porous Collisional-radiative ionisationRadiation: intermediate A cc r e ti on fl o w C h r o m o s ph e r e Transparent Lagrangian interface
Fig. 9. "Hybrid" simulation setup and boundary conditions. − . − . − . . . . r (km)6 . . . . . . l og T ( K ) − . − . − . − . − . . l og ξ Fig. 10.
Electron ionisation rate ( ξ = n e n H + n He + n M , light green ) andtemperature ( red ) profiles zoomed on the reverse shock front in theearly QPO cycle. The first cycle is presented in Figure 11; it shows the time vari-ations over 160 s of the gas temperature and mass density, ofthe photon escape ( ζ ) and absorption (1 − ζ ) probabilities and ofthe radiation energy volumetric density and flux ( E r and F r ) (seeSection 2.2.3.3) for the same snapshots. The next cycles only di ff er from this first one by the position of the interface betweenthe slab and the chromosphere, as discussed in Section 3.4.3.4.The global behaviour follows the trends of the two previousmodels. However, several e ff ects must be highlighted: a heatingof the chromosphere and of the accretion flow, already pointedout by Calvet & Gullbring (1998) and Costa et al. (2017), andthe reduction of the oscillation period and of the post-shock ex-tension. These e ff ects are discussed below.3.4.3.1. QPO cycle reductionAlthough 1 − ζ shows strong variations, its net value beyond theforward shock remains negligible, and the post-shock material isin the coronal limit (as in Section 3.2). The temperature behindthe reverse shock is here equal to 3.1 × K, to be comparedto 4 × K in the reference case. In addition, the compressionis enhanced from 4 (W– Λ case) to 4.4. As a consequence, thecooling is more e ffi cient: the cooling time is reduced from 400 sdown to 220 s, which is compatible with the duration of the cy-cles. This e ff ect is due to the ionisation / recombination energycost ( q χ ), which is included in the gas energy equation for theHybrid case, but not for reference case ( q χ = Λ case, cf.Eq. (1) and Section 2.1.2).3.4.3.2. Radiation energy and fluxThe radiation energy density increases between 9 and 100 s,which corresponds to the growth phase of the hot slab. Thisincrease is however correlated to the upper chromosphere heatedup to 12 000 K (discussed in Section 3.4.3.3 and presented Fig-ure 12). E r remains almost flat in the optically thin post-shockmedium, with a value driven by the heated upper chromosphere.In the accretion flow, during the growth of the hot slab, there is Article number, page 9 of 18 & A proofs: manuscript no. DE_SA-arXiv ||| − − − − − log ρ (cid:16) gcm − (cid:17) t = s ||| t = s ||| t = s ||| t = s ||| t = s ||| t = s ||| t = s ||| t = s . . . . . log T (K) ||| − . − . − . . . v (cid:16) × km / s (cid:17) ||||||||||||||||||||| . . . . . log n e (cm − ) ||| − − − log ζ ||||||||||||||||||||| − − − log (1 − ζ ) l og r ( k m ) . . . . . log E r (cid:16) erg / cm (cid:17) l og r ( k m ) l og r ( k m ) l og r ( k m ) l og r ( k m ) l og r ( k m ) l og r ( k m ) l og r ( k m ) − . . . . . F r (cid:16) × erg / cm / s (cid:17) F i g . . S n a p s ho t s o f t h e m a ss d e n s it y ( green ) , g a s t e m p e r a t u r e ( red ) , e s c ap e ζ ( darkblue ) a nd ab s o r p ti onp r obab iliti e s − ζ ( cyan , s ee S ec ti on2 . . . ) , v e l o c it y ( grey ) , e l ec t r ond e n s it y ( lightgreen ) , r a d i a ti on e n e r gyd e n s it y ( magenta ) a nd fl ux ( orange ) p r o fi l e s o f t h e fi r s t Q P O c y c l e w it h t h e " H yb r i d " s e t up . T h eacc r e t e dg a s f a ll s fr o m t h e r i gh tt o t h e l e f t on a n e qu ili b r i u m a t m o s ph e r e . Article number, page 10 of 18. de Sá et al.: New insight on Young Stellar Objects accretion shocks a tiny decrease due to the absorption by the accreted material upto 0 .
5% at 70 s (see Section 3.4.3.4).The radiative properties of the inner chromosphere is welldescribed by the di ff usive limit: E r (cid:39) a R T and F r (cid:39) cE r /
4. Themost peculiar feature of the radiative flux is its linear growththrough the post-shock slab. Such a pattern is characteristic of avolume emission by an optically thin medium. As discussed be-fore, the radiative energy in the hot slab is somehow imposed bythe heated upper chromosphere: as a consequence, the outgoingradiation flux is cE r / (cid:39) c a R T /
4, with T chr the temperature ofthe upper chromosphere (cf. Figure 12). Then, the radiation fluxpropagates through the accretion flow with negligible changes.Since F r is counted negatively towards the star, the flux emittedby the slab is o ff set by the one produced by the chromosphere.the net radiation flux rises then back in the chromosphere.However, the variations of F r within the slab comes from the in-terweaving of several radiation sources (the chromosphere, theslab itself and the accretion flow): due to the limitations of theM1 radiation transfer (cf. Section 2.2.2), these variations mustbe interpreted with care.3.4.3.3. Chromospheric heating and beating | | | | | | | | | r r e v ( × k m ) − − − − − v r e v ( × k m / s ) t (s)5678910111213 ( × K ) T c h r & T acc Fig. 12.
Position of the reverse shock ( r rev , purple ) and the correspond-ing velocity ( v rev , grey ) together with the accreted plasma temperaturebefore the hot slab ( T acc , dark gold ) and the upper chromosphere tem-perature ( T chr , red ) for the first QPO cycle in the Hybrid case. The upper chromosphere is heated by the radiating post-shock plasma up to 12 000 K (Figure 12). For instance, between9 and 70 s, its temperature varies from 7000 to 10 800 K at800 km, and the pressure increases from 800 to 2600 dyn / cm at this location (Figure 13). As a consequence, the whole post-shock structure is pushed upwards from 875 km to 3150 km, thusout 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 anymore and the slab buries back into the atmosphere. The expectedbehaviour is an oscillation of the slab burial with the same pe-riodicity as QPOs, since it originates from the hot post-shockplasma radiation. At the end of this first cycle, the chromospheredoes not recover its initial thickness: this e ff ect does not af- | | | |− l og p (cid:16) dyn / c m (cid:17) t = . . . . . . . l og T ( K ) r (km) − l og p (cid:16) dyn / c m (cid:17) t =
70 s . . . . . . . l og T ( K ) Fig. 13.
Snapshots of the temperature ( red ) and pressure ( blue ) at 9( top ) and 70 s ( bottom ) for the Hybrid case. fect the post-shock dynamics and cycle characteristics. All thesee ff ects are overestimated in a 1D model. However, this studyshows that the general question of the (un)burial, which is im-portant for X-ray observations, can only be addressed within amodel that takes into account the radiative heating of the chro-mosphere by the hot slab.3.4.3.4. Accretion flow pre-heatingWhile reaching the hydro-radiative steady state of the chromo-sphere, the flow has been homogeneously heated from 3000 Kto 5730 K before the start of the accretion process. During thecycle, the accretion leads to an additive heating of the flow upto ∼ t =
93 s). These e ff ects are quantified in Figure12, which reports the time variations of the position and veloc-ity of the interface between the hot slab and the accretion flow,as well as the temperatures of the heated chromosphere and ofthe pre-shock material. The use of the escape probability formal-ism (see Section 2.2.3.3) induces a dependence of the absorptionby the accretion flow with the section of the column; changingthis section from 1000 km to 10 000 km for instance will varythe parameter ζ from 1 − × − to 1 − × − , increasing theabsorption and thus the radiative heating of the pre-shock flow.Such preheating has already been pointed out by other au-thors (Calvet & Gullbring 1998; Costa et al. 2017). In theseworks, this heating is induced by radiation coming from thehot slab through photo-ionisation. Although radiative cooling ofthe accretion flow may be included in some cases, the radiationtransfer is not taken into account. Depending on the conditions,the pre-shock temperature may reach from 20 000 K (in CG98)up to 10 K (in Co17) close to the reverse shock (up to 10 km).In the latter, this precursor is preceded by a flatter and cooler( ∼ K) zone with an extension of 10 km, thus smaller thanours ( > km).Our simulation shows that part of the heating is a conse-quence of the chromospheric radiation already in play before thestart of the accretion. The analysis of the variation of the radia-tive energy indicates that an additional heating operates duringthe development of the hot slab. However, as we do not includeany dependence with the wavelength, it remains very di ffi cult todiscriminate in details the role played by the radiation emitted by Article number, page 11 of 18 & A proofs: manuscript no. DE_SA-arXiv the hot slab (X-rays) and from the (heated) chromosphere (UV-visible). Complementary information will be given by the syn-thetic spectra computed as a post-process of the hydrodynamicsstructures (Section 3.4.4).3.4.3.5. X-ray luminosityThe X-ray luminosity of the system is computed following themethod described in Section 3.2.3. Its time variation (in unit of τ cycle =
160 s) is reported in Figure 6 for comparison with thetwo previous cases. Compared to the reference case, in additionwith a shortening of the period, this case presents a more pro-nounced radiative collapse (70–90 s), followed by a chaotic col-lapse (90–160 s). The time average of the radiative surface lumi-nosity is here equal to 1.2 × erg cm − s − , which represents30% of the mechanical energy flux (cf. Figure 6). As this simulation is performed using only one group of radia-tion frequencies, it is interesting to analyse more precisely thedetails of the previous radiative heating via its feedback on themonochromatic emergent intensity.To this purpose, the hydrodynamic structures has been post-processed with the SYNSPEC code (Section 2.3.3). For con-sistency purpose, we take the atomic data already used forthe calculation of the average opacities (see Section 2.2.3 andAppendix A). We thus estimate the specific intensity I (cid:113) λ (inerg / cm / s / Å / sr) along the direction of the column. Since the lineprofile behaviour is not investigated here, velocity e ff ects are ne-glected.It is important to recall that a quantitative comparison of thissynthetic spectrum with observations, especially in the X-rays(see e.g. Güdel et al. 2007; Robrade & Schmitt 2007; Drakeet al. 2009) would require NLTE and 3D radiative transfer post-processing. Nonetheless, using 1D radiative transfer and the LTEapproximation is here interesting as it corroborates or not thegeneral accepted trends, e.g. a strong X-ray emission and anexcess of luminosity in the UV-VIS range (Calvet & Gullbring1998; Brickhouse et al. 2010; Ingleby et al. 2013). r ( × km) − − − − l og ρ (cid:16) g c m − (cid:17) t =
70 s 4567 l og T ( K ) Fig. 14.
Snapshot of the density ( green ) and temperature ( red ) profilesat t =
70 s with the "Hybrid" setup, post-processed hereafter.
A typical spectrum emerging from 4.6 × km (locatedwithin the accretion flow) is reported in Figure 15). It is com-puted from a snapshot ( t =
70 s) of the
Hybrid model (see Fig-ures 11 and 14). At this stage, the chromosphere extends up to1.4 × km, the hot plasma from 1.4 × km to 8.3 × kmand the accretion flow from 8.3 × km to 1 × km. The in-tensity that emerges from this layer presents three characteristicspectral bands: λ (Å)4681012 l og I q λ ( e r g / c m / s / Å / s r) Fig. 15.
Specific intensity I (cid:113) λ parallel to the column during the QPOcycle of the Hybrid model (at t =
70 s, see Figure 14). – in the range 1–100 Å (X-rays), the bump is attributedto the hot post-shock plasma, with intense lines up to10 erg / cm / s / Å; – in the range 100–900 Å (EUV), radiation is e ffi ciently ab-sorbed by the inflow; – in the range 900–10 000 Å (UV + Vis + IR), the second bumpis attributed to the heated stellar chromosphere and photo-sphere, i.e. a black body at T (cid:39)
11 300 K (cf. Figure 12). λ (Å)4681012 l og I q λ ( e r g / c m / s / Å / s r) Fig. 16.
Specific intensity I (cid:113) λ emerging from the reverse shock front dur-ing the QPO cycle of the Hybrid model (at t =
70 s, see Figure 14).
The strong absorption of the EUV radiation is due to the hugeoptical depth of the accretion flow . This e ff ect may then be at-tenuated in the case of a bent column or when the observationis performed side-on and not along the column. This absorptione ff ect on the spectrum is illustrated in Figure 16, which presentsthe intensity emerging right after the reverse shock front , at r = . × km. This figure shows that this absorption also af-fects, to a lesser degree, the visible spectrum originating from thechromosphere. This must be considered when interpreting theUV excess (see e.g. Calvet & Gullbring 1998; Hartmann et al.2016; Colombo et al. 2019). Note that a pre-heating of the accre-tion flow is expected as a result of the EUV absorption. A pre-heating is also obtained independently by AstroLabE (Section3.4.3.4); however, a one-to-one correspondence would require amulti-group description of the radiation field in AstroLabE.We compute from I (cid:113) λ (Figure 14) the net X-ray outgoing in-tensity ( I (cid:113) X ) and the corresponding coronal quantity ( I Λ ): I (cid:113) X = (cid:90) I (cid:113) λ d λ and I Λ = L Λ π (11)The time variations of these two quantities, reported in Figure17, present similar characteristics. However, the values derivedby SYNSPEC are higher by about two to three orders of mag-nitude. This discrepancy is either imputable to the LTE approxi-mation or to the assumed 1D plane-parallel geometry. Thus our at ρ = − g cm − and T (cid:39) t (s)67891011 l og I ( e r g / c m / s / s r) Fig. 17.
Time-variation of the 2–27 Å (LTE) integrated X-ray outgoingintensity ( I (cid:113) X , blue ) and of the optically thin post-shock emission ( I Λ red ) during a QPO cycle for the Hybrid setup. synthetic spectra can’t be used for quantitative comparison withobservations.
4. Refining the models
It should be pointed out that this study uses a solar model for thechromosphere with acoustic heating. Compared to the descrip-tion of this heating, a more important improvement would be toconsider a realistic T Tauri chromospheric model, which is todaynot very well known. This may a ff ect the ionisation (and thengas pressure with another chemical abundances) as well as slabcharacteristics (through gravity) and radiation e ff ects (throughopacities and incoming luminosity). Our results are then to beconsidered qualitatively and not quantitatively. We use in this work radiation momenta equations with the M1closure relation. Although this is already a strong improvementcompared to other approaches like the di ff usion model, it couldbe improved by using radiation half-fluxes (i.e. the inward anoutward components of the radiation flux). This should disen-tangle the radiation flux coming from the star and from the post-shock structure.The M1 closure relation allows the radiation field to reach atmost one direction of anisotropy; half-fluxes can extend it to two,i.e. the maximum number of anisotropy directions reachable in1D. Half-fluxes (along with M1) would then be equivalent to themomenta equations with the M2 closure relation (Feugeas 2004),without its prohibitive numerical cost.The M1 model and its limits have been thoroughly studied (seee.g. Levermore 1996; Dubroca & Feugeas 1999; Feugeas 2004).The behaviour of this model along with half-fluxes needs how-ever to be examined.More important is the approximation made with the monogroup approach used in this work. The whole spectrum isthen approximated as a black body providing the adequate opac-ity averages. However, our computed spectra emerging from ac-cretion structures are expected to present three discernible fre-quency groups: – up to the visible domain, the spectrum is dominated by theblack body emerging from the stellar photosphere; – the EUV band is expected to be depleted due to high absorp-tion by the accreted gas; – the X-ray band is thought to be optically thin and to have thehot slab signature on it. Although the multi-group approach is numerically heavier, it willimprove the study of the consequences of the radiation absorp-tion by the surrounding medium. A consequence of the X-rayand EUV absorption by the cold accretion flow is the presenceof a radiative precursor. Such a phenomenon cannot be obtainedthrough a monogroup approach. Moreover, a 3 groups approachwill provide a better description of the feedback of the hot slabon the stellar chromosphere. Two other points may be improved. First, the transition model( ζ ) remains qualitative and may need to be extended to the ion-isation calculation. The work done by Carlsson & Leenaarts(2012) o ff ers paths to reach such consistency and may need tobe investigated further. A better model of both the LTE trans-fer, line cooling and intermediate regimes may demand ded-icated NLTE opacities, namely plasma emissivity (equivalentto n e n H Λ ), radiation energy absorption ( κ P ) and radiation fluxsinking ( κ R ). Moreover, all these quantities, computed with aradiative-collisionnal model, have to be averaged over adequateweighting functions. Due to recent progresses in this topic (Ro-driguez et al. 2018), new results are expected in a near future.Independently, a NLTE description should be used to computethe emerging spectra: this work is already in progress usingTLUSTY code.
5. Conclusion
In this study, we used 1D simulations with detailed physics tocheck the validity of the two following common assumptions inaccretion shocks simulations: the stellar atmosphere can be ei-ther modelled by a hydrostatic or a steady hydrodynamic struc-ture, and the dynamics of accretion shocks is governed by op-tically thin radiation transfer. We checked first that we are ableto recover previous results (Sacco et al. 2008, W– Λ case, Sec-tion 3.2) and tested independently each of these assumptions(Chr– Λ case, Section 3.3, and Hybrid case, Section 3.4). Eachof them proves to have a non-negligible impact on the typicalcharacteristics of the accretion dynamics and on the estimationof its X-ray surface luminosity. This one varies between 1 and4 × erg cm − s − . Taking as a reference the radiative power ofBrickhouse et al. (2010), we derive a section of the accretion spotfrom 3 × to 1 × cm , corresponding to a filling factorof the solar disk between 2 and 8%, i.e. a stream composed of ∼ fibrils of radius L c (cf. Section 3.1) or a column of radius100 × L c , supposing that the global dynamics of the system is notinfluenced by this larger section of the column through radiativee ff ects.In the case of the chromosphere which is heated by acousticperturbations that degenerate into small shock waves, we haveshown that these perturbations do not strongly modify the cycleperiod compared to the reference case. However, the cycle be-comes chaotic due to the generation of secondary shock waves.As a result, the relative duration of the hot phase in the cycleremains longer, and thus the variability in the X-rays is lesspronounced than for the reference case. To be detected, it wouldrequire a better sensitivity of the photometric measurements.In the case of an initially steady atmosphere at radiative equi-librium, the coupling between the radiation and the hydrodynam-ics leads to: Article number, page 13 of 18 & A proofs: manuscript no. DE_SA-arXiv – a radiative feedback (heating) of the atmosphere, which suc-cessively expands and retracts, inducing in particular an un-burial of the column, which is favorable to the lateral escapeof the X-ray emitted from the hot slab; – a chaotic radiative collapse, with an impact on the time vari-ation of the X-ray flux (Figure 6); – a radiative pre-heating of the incoming flow, over the lengthof the simulation box.Moreover the inclusion of ionization in the energy balance leadsto important e ff ects in the post-shock temperature that modifythe cooling e ffi ciency and therefore the cycle duration.In this hybrid case, we computed at LTE the radiative inten-sity 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 ofemission and absorption lines, –3– a EUV component which isvery strong at the position of the reverse shock and disappears atthe outer boundary, due to the importance of the absorption.This study could be completed with a more complete sim-ulation that would include both a dynamically-heated chromo-sphere and the hybrid setup. However, it appears at this stagemore important to take into account a NLTE radiative descrip-tion based on adapted opacities and radiative power losses. An-other necessary improvement will be through a multi-group ra-diation transfer to catch at least the e ff ect of EUV absorptionand X-ray radiative losses on the structure of the column, andto analyse the possibility of a radiative precursor which couldpre-heat the incoming flow. The study is also to be extended tomulti-dimensional simulations in order to check the e ff ects ofboth radiation and magnetic field closer to the real picture (Or-lando et al. 2010, 2013; Matsakos et al. 2013, 2014). Acknowledgements.
We express our gratitude to Jason Ferguson for providingus with the molecular LTE opacity tables used in this work and to Franck De-lahaye for his contribution to the atomic LTE ionisation data used for both theopacity tables and for the radiative transfer post processing. We thank Rafael Ro-driguez for sharing with us helpful preliminary results about NLTE microscopiccollisional-radiative data, Salvatore Orlando for reading the manuscript, ZianeIzri for its contribution at the beginning of this project and Christophe Sauty forhelpful discussions.I.H. thanks the Physics Department of Sorbonne Université for his visiting pro-fessorship.This work was supported by the french ANR StarShock and LabEx Plas@Parprojects (resp. ANR–08–BLAN–0263–07 and ANR–11–IDEX–0004–02), PICS6838, Programme National de Physique Stellaire of CNRS / INSU and Observa-toire de Paris.
References
Alfvén, H. & Lindblad, B. 1947, MNRAS, 107, 211Argiro ffi , C., Maggio, A., & Peres, G. 2007, A&A, 465, L5Argiro ffi , C., Maggio, A., Peres, G., et al. 2009, A&A, 507, 939Auer, L. H. 2003, in Stellar Atmosphere Modeling, Vol. 288 (ASPC), 3–15Ayres, T. R. 1979, ApJ, 228, 509Batalha, C. C. & Basri, G. 1993, ApJ, 412, 363Biermann, L. 1946, Naturwissenschaften, 33, 118Bonito, R., Orlando, S., Argiro ffi , C., et al. 2014, ApJ, 795, L34Bouquet, S., Stehlé, C., Koenig, M., et al. 2004, Phys. Rev. Lett., 92, 225001Bouvier, J., Covino, E., Kovo, O., et al. 1995, A&A, 299, 89Brickhouse, N. S., Cranmer, S. R., Dupree, A. K., Luna, G. J. M., & Wolk, S. J.2010, ApJ, 710, 1835Brown, J. C. 1973, Sol. Phys., 29, 421Calvet, N. 1983, Rev. Mexicana Astron. Astrofis., 7, 169Calvet, N., Basri, G., & Kuhi, L. V. 1984, ApJ, 277, 725Calvet, N. & Gullbring, E. 1998, ApJ, 509, 802Carlsson, M. & Leenaarts, J. 2012, A&A, 539, A39Charignon, C. & Chièze, J.-P. 2013, A&A, 550, A105Chevalier, R. A. & Imamura, J. N. 1982, ApJ, 261, 543Chièze, J.-P., de Sá, L., & Stehlé, C. 2012, EAS Publications Series, 58, 143 Colombo, S., Ibgui, L., Orlando, S., et al. 2019, A&A, submittedColombo, S., Orlando, S., Peres, G., Argiro ffi , C., & Reale, F. 2016, A&A, 594,A93Costa, G., Orlando, S., Peres, G., Argiro ffi , C., & Bonito, R. 2017, A&A, 597,A1Cram, L. E. 1979, ApJ, 234, 949Curran, R. L., Argiro ffi , C., Sacco, G. G., et al. 2011, A&A, 526, A104de Sá, L. 2014, PhD thesis, Université de Paris VIde Sá, L., Chièze, J.-P., Stehlé, C., et al. 2012, in SF2A-2012: Proceedings ofthe Annual meeting of the French Society of Astronomy and Astrophysics,ed. S. Boissier, P. de Laverny, N. Nardetto, R. Samadi, D. Valls-Gabaud, &H. Wozniak, 309–312de Sá, L., Chièze, J.-P., Stehlé, C., et al. 2014, in European Physical Journal Webof Conferences, Vol. 64, 04002Dorfi, E. A. & Drury, L. O. 1987, JCoPh, 69, 175Drake, J. J. 2005, in 13th Cambridge Workshop on Cool Stars, Stellar Systemsand the Sun, 519Drake, J. J., Ratzla ff , P. W., Laming, J. M., & Raymond, J. C. 2009, ApJ, 703,1224Dubroca, B. & Feugeas, J.-L. 1999, CRAS Paris Série 1, 329, 915Dumont, S., Heidmann, N., Kuhi, L. V., & Thomas, R. N. 1973, A&A, 29, 199Feigelson, E. D. & Montmerle, T. 1999, ARA&A, 37, 363Ferguson, J. W., Alexander, D. R., Allard, F., et al. 2005, ApJ, 623, 585Feugeas, J.-L. 2004, LPB, 22, 121Fitzpatrick, E. L. 1996, ApJ, 473, L55Fritsch, F. N. & Butland, J. 1984, SIAM J. Sci. and Stat. Comput., 5, 300Grevesse, N. & Sauval, A. J. 1998, Space Sci. Rev., 85, 161Güdel, M., Skinner, S. L., Mel’nikov, S. Y., et al. 2007, A&A, 468, 529Gullbring, E., Calvet, N., Muzerolle, J., & Hartmann, L. W. 2000, ApJ, 544, 927Günther, H. M., Lewandowska, N., Hundertmark, M. P. G., et al. 2010, A&A,518, A54Günther, H. M., Liefke, C., Schmitt, J. H. M. M., Robrade, J., & Ness, J. U. 2006,A&A, 459, L29Günther, H. M., Schmitt, J. H. M. M., Robrade, J., & Liefke, C. 2007, A&A,466, 1111Hartmann, L. W., Herczeg, G. J., & Calvet, N. 2016, ARA&A, 54, 135Hubeny, I. & Lanz, T. 1995, ApJ, 439, 875Hubeny, I. & Lanz, T. 2017, eprint arXiv:1706.01859,Hubeny, I. & Mihalas, D. 2014, Theory of Stellar Atmospheres An Introductionto Astrophysical Non-equilibrium Quantitative Spectroscopic Analysis, 1stedn. (Princeton University Press)Huenemoerder, D. P., Kastner, J. H., Testa, P., Schulz, N. S., & Weintraub, D. A.2007, ApJ, 671, 592Hui, L. & Gnedin, N. Y. 1997, MNRAS, 292, 27Ibgui, L., Hubeny, I., Lanz, T., & Stehlé, C. 2013, A&A, 549, A126Ingleby, L., Calvet, N., Herczeg, G. J., et al. 2013, ApJ, 767, 112Jess, D. B., Morton, R. J., Verth, G., et al. 2015, Space Sci. Rev., 190, 103Jiang, Y.-F., Stone, J. M., & Davis, S. W. 2014a, ApJ, 796, 106Jiang, Y.-F., Stone, J. M., & Davis, S. W. 2014b, ApJS, 213, 7Johns-Krull, C. M. 2007, ApJ, 664, 975Johns-Krull, C. M., Valenti, J. A., Hatzes, A. P., & Kanaan, A. 1999, ApJ, 510,L41Judge, P. 2006, in Solar MHD Theory and Observations: A High Spatial Reso-lution Perspective, Vol. 354 (Astronomical Society of the Pacific ConferenceSeries), 259Kalkofen, W. 2007, ApJ, 671, 2154Kastner, J. H., Huenemoerder, D. P., Schulz, N. S., Canizares, C. R., & Wein-traub, D. A. 2002, ApJ, 567, 434Kirienko, A. B. 1993, AstL, 19, 11Koldoba, A. V., Ustyugova, G. V., Romanova, M. M., & Lovelace, R. V. E. 2008,MNRAS, 388, 357Lequeux, J. 2005, The Interstellar Medium, Astronomy and Astrophysics Li-brary (Berlin / Heidelberg: Springer-Verlag)Lesa ff re, P. 2002, PhD thesis, Université Paris VIILesa ff re, P., Chièze, J.-P., Cabrit, S., & Pineau des Forêts, G. 2004, A&A, 427,147Levermore, C. D. 1996, JSP, 83, 1021Lowrie, R. B., Mihalas, D., & Morel, J. E. 2001, JQSRT, 69, 291Matsakos, T., Chièze, J.-P., Stehlé, C., et al. 2013, A&A, 557, A69Matsakos, T., Chièze, J.-P., Stehlé, C., et al. 2014, Proceedings of the Interna-tional Astronomical Union, 9, 66Mignone, A. 2005, ApJ, 626, 373Mihalas, D. & Mihalas, B. W. 1984, Foundations of radiation hydrodynamics(New York, Oxford University Press)Muzerolle, J., Calvet, N., & Hartmann, L. W. 1998, ApJ, 492, 743Opacity Project Team. 1995, The Opacity Project, Vol. 1 (Institute of PhysicsPublications, Bristol, UK)Orlando, S., Bonito, R., Argiro ffi , C., et al. 2013, A&A, 559, A127Orlando, S., Sacco, G. G., Argiro ffi , C., et al. 2010, A&A, 510, A71 Article number, page 14 of 18. de Sá et al.: New insight on Young Stellar Objects accretion shocks
Oxenius, J. 1986, Kinetic theory of particles and photons. Theoretical founda-tions of Non-LTE plasma spectroscopy (Springer Series in Electrophysics,Berlin: Springer)Peres, G., Rosner, R., Serio, S., & Vaiana, G. S. 1982, ApJ, 252, 791Press, W. H., Teukolsky, S. A., Vetterling, W. T., & Flannery, B. P. 1994, inFortran Numerical Recipes (Cambridge: University Press)Rammacher, W. & Ulmschneider, P. 1992, A&A, 253, 586Robrade, J. & Schmitt, J. H. M. M. 2007, A&A, 473, 229Rodriguez, R., Espinosa, G., & Miguel Gil, J. 2018, Physical Review E, 98,033213Sacco, G. G., Argiro ffi , C., Orlando, S., et al. 2008, ApJ, 491, L17Sacco, G. G., Orlando, S., Argiro ffi , C., et al. 2010, A&A, 522, A55S˛adowski, A., Narayan, R., McKinney, J. C., & Tchekhovskoy, A. 2014, MN-RAS, 439, 503Schmitt, J. H. M. M., Robrade, J., Ness, J. U., Favata, F., & Stelzer, B. 2005,A&A, 432, L35Schwarzschild, M. 1948, ApJ, 107, 1Siwak, M., Ogloza, W., Mo ff at, A. F. J., et al. 2018, MNRAS, 478, 758Sobotka, M., Heinzel, P., Švanda, M., et al. 2016, ApJ, 826, 49Spitzer, L. 1998, Physical Processes in the Interstellar Medium (Wiley-VCH)Spitzer, L. & Härm, R. 1953, PhRv, 89, 977Stehlé, C. & Chièze, J.-P. 2002, in SF2A-2002: Semaine de l’AstrophysiqueFrancaise, ed. F. Combes & D. Barret, 493Stelzer, B. & Schmitt, J. H. M. M. 2004, A&A, 418, 687Ulmschneider, P., Rammacher, W., Musielak, Z. E., & Kalkofen, W. 2005, ApJ,631, L155van Leer, B. 1973, in Proceedings of the Third International Conference on Nu-merical Methods in Fluid Mechanics, Vol. 1 (Springer, New York), 163–168Vernazza, J. E., Avrett, E. H., & Loeser, R. 1973, ApJ, 184, 605Verner, D. A. & Ferland, G. J. 1996, ApJS, 103, 467Vidal, F., Matte, J. P., Casanova, M., & Larroche, O. 1995, Phys. Plasmas, 2,1412von Neumann, J. & Richtmyer, R. D. 1950, J. Appl. Phys., 21, 232Voronov, G. S. 1997, Atomic Data and Nuclear Data Tables, 65, 1Walder, R. & Folini, D. 1996, A&A, 315, 265Yan, M., Sadeghpour, H. R., & Dalgarno, A. 1998, ApJ, 496, 1044 Article number, page 15 of 18 & A proofs: manuscript no. DE_SA-arXiv
Appendix A: Opacity tables
The specificity of the accretion shocks study led us to work ondedicated opacity tables. We expose in the appendix the reasonsbehind this choice and the creation process. The resulting opac-ity table is accessible upon request.
Appendix A.1: Motivation − − − ρ (cid:16) g cm − (cid:17) l og T ( K ) log κ P − − − ρ (cid:16) g cm − (cid:17) log κ R − . − . − . . . . . . Fig. A.1.
Planck ( left ) and Rosseland ( right ) opacities (in cm g − )with respect to gas density and temperature – any in log scale, as pro-vided by the Opacity Project (Opacity Project Team 1995). The blackcurve is a typical characteristic of an accretion column. Most available opacity tables are defined on a slanted ( ρ, T ) or( n e , T ) domain (see e.g. Figure A.1). However, one peculiarityof accretion shock structures is the presence of a low densityhot post-shock plasma (black curve vertex in Figure A.1) thatexplores a domain uncovered by publicly available tables. Morecomplete tables are thus mandatory for the present study. Appendix A.2: Choice of primary tables l og T ( K ) − − − − − ρ (cid:16) g cm − (cid:17) l og T ( K ) − − − − − ρ (cid:16) g cm − (cid:17) − − − − Fig. A.2.
SYNSPEC ( top ) and Ferguson ( bottom ) Planck ( left ) andRosseland ( right ) opacities (in cm g − ) with respect to gas densityand temperature – any in log scale. The dotted lines show transitiontemperatures chosen for each table (see Annexe A.4.1). To cover the density and temperature range correspondingto our conditions, we implement in the code SYNSPEC (see Section 2.3.3), initially dedicated to stellar atmospheres, mod-ules allowing to generate LTE monochormatic opacities at agiven density and temperature. These monochromatic opacitieswere then averaged with the proper weighting functions togenerate the adequate Rosseland and Planck mean opacitiestables (hereafter called "SYNSPEC tables", see Figure A.2, toppanels). These opacities are consistent with Opacity Project (seee.g. Opacity Project Team 1995) data, that we use as reference,for T between 10 K and 10
K. The advantage of SYNSPECcomes from the high number of atomic species considered, sincea very detailed chemical composition is necessary to model theradiation properties of a plasma at high temperatures.However, below 10
K, the molecular chemistry cannotbe neglected, but is not included in this work on SYNSPEC.We completed thus the SYNSPEC tables with low temperaturemolecular opacities provided by Ferguson et al. (2005) between10 K and 10 K ("Ferguson tables", see Figure A.2, bottom pan-els), that show excellent agreement with Opacity Project at up-per temperatures. To facilitate the merging process, we obtainedfrom the authors tables with compatible density and temper-ature grid (Ferguson, priv. comm.): mesh points from Fergu-son and SYNSPEC tables are identical in the common domain(10 –10 K and 10 − –10 − g cm − ). Appendix A.3: Preliminary study
Appendix A.3.1: Analysis of primary tables
Considering opacity variations as well as temperature and den-sity ranges, we decided to work with the logarithm of all thesequantities. As first derivatives, we use then: ∂ lT l κ = ∂ log κ∂ log T = T κ ∂κ∂ T ∂ l ρ l κ = ∂ log κ∂ log ρ = ρκ ∂κ∂ρ (A.1)where κ stands for κ P or κ R . l og T ( K ) -14 -12 -10log ρ (cid:16) g cm − (cid:17) l og T ( K ) − − − − . . . Fig. A.3. ∂ lT l κ R fromSYNSPEC table ( top )and ∂ l ρ l κ P from Fergusontable ( bottom ) with re-spect to the logarithm ofgas density and tempera-ture. The dotted lines show transition temper-atures chosen for eachtable. Some anomaliesare revealed, especiallyaround 10 K and highdensities: this zone is cut-ted during the mergingprocess.
Preliminary analysis of SYNSPEC and Ferguson tables re-vealed local aberrations, especially looking at the temperature ordensity derivatives (see Figure A.3). We may use the mergingprocess to smooth most aberrations.
Article number, page 16 of 18. de Sá et al.: New insight on Young Stellar Objects accretion shocks
Appendix A.3.2: Physical and numerical constraints
In order to get a satisfying merging, several numerical and phys-ical constraints must be respected: – as far as possible, opacities must be of class C (values andfirst derivatives must be continuous); – the transition region should be as narrow as possible; – the transition region must encompass anomalies encounteredin both primary tables.Such a table is composed of a limited number of discretepoints: the first constraint can be reported to the interpolationmethod as far as opacity values in the transition present smoothvariations.To ensure a smooth transition between the molecular and theatomic (primary) tables, the transition must not take into con-sideration the values within the transition. The transition valuesloose then any physical meaning, and must be as few as possi-ble . Appendix A.4: Merging process
Appendix A.4.1: Method
In this Section, the index "A" refers to values taken at the lowertransition temperature, as the index "B" for the upper ones. Thetransition temperatures chosen to merge SYNSPEC and Fergu-son tables are: – T A = . K and T B = . K for κ P ( ∼ – T A = . K and T B = . K for κ R ( ∼ Appendix A.4.2: Merging along temperature
To satisfy the class C constraint, we combined (see for instanceAuer 2003 and Ibgui et al. 2013 ): – piecewise cubic Hermite polynomials, which ensure continu-ity of values ( κ A , κ B ) and derivatives ( ∂ lT l κ A , ∂ lT l κ B ) at eachtransition limit; – van Leer (1973) slopes to compute ∂ lT l κ A and ∂ lT l κ B , so asto prevent the apparition of spurious extrema in forcing theirlocation 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:log κ ( T i ; ρ j ) = u i (3 − u i ) log κ B + ( u i − u i h ∂ lT l κ B + ( u i − (2 u i +
1) log κ A + ( u i − u i h ∂ lT l κ A (A.2)with h = log T B − log T A and u i = (log T i − log T A ) / h . Thisexpression can be rewritten as a 3 rd degree polynomial in u i . We note that the Ferguson tables showed opacity discontinuities intheir hottest and densest part, as SYNSPEC tables in their coolest anddensest part ( ∗ ) (see Figure A.3). Since the values within the transitionregion are ignored, we use it to artificially remove anomalies: as far aspossible, the transition region must be chosen so that it covers most ofthem. ( ∗ ) Few anomalies remains in regions that are not explored in our simu-lations (see Figure A.4); this problem is postponed for now. Fritsch & Butland (1984) derivatives are used in these papers; theygeneralise van Leer slopes to non-regular grids.
Appendix A.4.3: Density correction
At this stage, we reached class C along temperature, but there isno guarantee of continuity along density. However, in practice,it was C , except for few mesh temperatures T i ∗ .Since the dependency in density is held by the 3 rd degreepolynomial coe ffi cients, we look at the behaviour of each ofthem with respect to density. Every coe ffi cient showed spuriousvariations nowhere but at densities ρ j ∗ . We apply then piecewisecubic Hermite polynomials along with van Leer slopes (densityderivatives) to estimate these coe ffi cients for each ρ j ∗ . These newcoe ffi cients are then used to reestimate opacity values along thetemperature for the ρ j ∗ . Appendix A.4.4: Final tables – interpolation process − − − − . − . − . . . . . . l og T ( K ) ∂ lT l κ P ∂ lT l κ R -14 -12 -10 -8 -6log ρ (cid:16) g cm − (cid:17) l og T ( K ) ∂ l ρ l κ P -14 -12 -10 -8 -6log ρ (cid:16) g cm − (cid:17) ∂ l ρ l κ R Fig. A.4.
Merged table Planck ( left ) and Rosseland ( right ) opacitytemperature ( top ) and density ( bottom ) first derivatives in the ( ρ , T )plane, any in log scale; the grey shape represents the transition region. We checked smoothness of the result by looking at the firstderivatives. Figure A.4 shows no anomaly within the transitiontemperature range [ T A , T B ] (grey shape). The remaining anoma-lies are not reached in our simulations.The interpolation process is copied from the mergingmethod, i.e. piecewise cubic Hermite polynomials along withvan Leer slopes, since it satisfies criteria described in SectionA.3.2. Interpolation is first performed along temperature at the2 × Article number, page 17 of 18 & A proofs: manuscript no. DE_SA-arXiv is arguably due to stronger variations of opacities (especiallyPlanck opacity) with respect to temperature.
Appendix B: Chromospheric model
One of our objectives is to describe the dynamics of the columnand its impact on the chromosphere, as well as the feedback ofthe chromosphere on the column. This requires then to includean adequate description of the physical mechanism leading tothe chromospheric heating. This appendix presents the simplebut 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 accretioncolumn lies in the stellar chromosphere: the dynamics and ob-servability of the column base may then depend on its structureand dynamics. Moreover, the chromosphere may be heated lo-cally by the accretion process. The inner heating mechanism inthe chromosphere is still subject of debates: it is mainly thoughtto originate either from acoustic waves dissipation (Biermann1946; Schwarzschild 1948; or more recently Sobotka et al. 2016)or from MHD waves dissipation (Alfvén & Lindblad 1947; Jesset al. 2015).Most accretion simulations model the stellar atmosphere –when it is modelled – as a hydrostatic plasma layer "tuned up"with ad-hoc sources to recover both temperature and pressureprofiles (see e.g. the heating function empirically introduced byPeres et al. 1982). Although this must work for a static structure,it is delicate to predict the dynamic behaviour of such a struc-ture facing the continuous perturbation from an infalling plasmaflow: such solution is not adapted to studies involving (in a self-consistent way) the dynamics of a perturbed atmosphere, like inthe context of accretion.We do not pretend to develop a "state of the art" model inthis paper: we only aim at using a reasonable model that is bothdynamic and self-consistent with our radiation hydrodynamicsmodel. In our 1D model, we do not consider any magnetic ef-fect but a very e ff ective confinement of the accretion flow alongthe field lines. To allow fast qualitative comparison between ourmodel and theoretical models & observations (see Figure B.1),we only used solar parameters (i.e. abundances, luminosity, massand radius). Appendix B.2: Acoustic waves and shocks
Acoustic waves are generated by photospheric granulation (seee.g. Judge 2006). These waves propagate upwards up to theheight where their velocity overcome the local sound speed, anddegenerate then into shocks. The nature of this mechanism israndom: two di ff erent locations at the stellar surface will becrossed over by acoustic shocks that ought to be out of phaseone with each other.In our simulations, acoustic energy is supplied in the form ofa monochromatic sinusoidal motion of the first Lagrangian inter-face ( T =
60 s and f acc = erg cm − s − , see e.g. Rammacher& Ulmschneider 1992; Ulmschneider et al. 2005; Kalkofen2007). Resulting acoustic waves propagate and degenerate intoshocks. Figure B.1 shows several temperature snapshots of suchsimulation along with the chromospheric model from Vernazzaet al. (1973). Below 300 km, acoustic waves are damped andhardly appear on snapshots. Above 500 km, waves are fully de- generated into shocks: their strength is then governed by the bal-ance between steepening in the pressure gradient and dissipation.Since the corona and the upper chromosphere (above 10 km) arereadily crushed by the accretion flow, the heating of these areasis not considered in our model. r (10 km)246810 T (cid:16) K (cid:17) Fig. B.1.
Successive snapshots of acoustic waves propagation ( thinlines ) and mean chromospheric temperature ( thick red line , Ver-nazza et al. 1973). The simulation setup is described in Section 3.3.1; r = Appendix B.3: From solar to stellar chromosphere
Observations of the solar chromosphere provide time and spaceaverages of thermodynamics quantities ( ρ , T , p , . . . ). Detailedobservation of CTTS chromospheres would demand higherspace and time resolution than the ones permitted by current ob-servational technologies. Most works on this field rely then onscaling laws (see e.g. Ayres 1979; Calvet 1983) or ad hoc fittingsto 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 Hybridmodel
This work encompasses several radiation regimes, from opticallythick LTE radiation transfer (Section 2.2.3.1) to optically thincoronal NLTE regime (Section 2.2.3.2). The momenta equations(Section 2.2) can handle all of them, assuming the proper radia-tion source terms are provided.In the LTE case, both radiation energy and momentumsource terms are well defined (Eq. (4)). In coronal regime, this isnot the case. Gas and radiation are decoupled in such a regime.Radiation only acts then as a gas energy sink: the radiation en-ergy source term (the gas sink) boils down to a cooling function(Kirienko 1993, see e.g.). Computing the radiation flux is irrel-evant in such regime and then no radiation momentum sourceterm is provided. That is why we set s † M r to .In the "Hybrid" setup, we aim at modelling radiative con-ditions that are neither LTE nor coronal regimes but somethingin between. To determine if the situation is closer to one or theother, and how close, we choose to look at the probability for aphoton to escape the accretion column (see Eq. (7)). We use itas a weighting factor to average the source terms, as shown inSection 2.2.3.3.The process is straightforward for the radiation energy sourceterm, but not for the radiation momentum source term since s † M r remains unknown. We assume then that the coronal Rosselandmean opacity may not significantly di ff er from its LTE value.This intuition is reinforced by preliminary calculations concern-ing NLTE radiative collisional opacities (Pérez, priv. com.).er from its LTE value.This intuition is reinforced by preliminary calculations concern-ing NLTE radiative collisional opacities (Pérez, priv. com.).