Impacts of Energetic Particles from T Tauri Flares on Inner Protoplanetary Discs
Valentin Brunn Christian Rab Alexandre Marcowith Christophe Sauty Marco Padovani Chadi Meskini


firstpage–lastpage
From weak to strong-coupling superconductivity tuned by substrate in TiN films
    Zhi-Rong Lin
    Accepted XXX, Received YYY; in original form ZZZ
==============================================================================





T Tauri stars are known to be magnetically active stars subject to strong flares observed in X-rays. These flares are likely due to intense magnetic reconnection events during which a part of the stored magnetic energy is converted into kinetic energy of supra-thermal particles. Since T Tauri stars are surrounded by an accretion disc, these particles may influence the disc dynamics and chemistry. This work continues on a previous stationary model, which showed that energetic particles accelerated during flares can produce a strong ionisation rate at high column densities in the inner accretion disc. The present model includes non-stationary sequences of flaring events sampled by a Chandra X-ray survey of nearby young stellar objects. We calculate the averaged ionisation rate expected in a radius range from 0.08 to 0.6 au from the central star. We confirm that energetic particles produced by the flares dominate the ionisation of the disc up to column densities of 10^25 cm^-2. We further study the main consequences of this additional source of ionisation on the viscosity, the accretion rate, the volumetric heating rate and the chemical complexity of inner protoplanetary discs.





acceleration of particles – magnetic reconnection  – stars: flare – accretion, accretion discs 








§ INTRODUCTION


T Tauri stars are young stars, typically of a few million years old, still contracting and accreting material from their surrounding circumstellar discs. The way T Tauri stars accrete the surrounding infalling matter is still an open subject in astrophysics <cit.>. All mechanisms that describe accretion-ejection processes involve the interaction between magnetic fields and ionised disc matter (e.g., ). 

The magnetic field geometry of these stars is much more complex than those of non-accreting cool stars <cit.>. Some field lines might become twisted due to the differential rotation between the disc and the star, causing them to bulge outwards and possibly even eject matter <cit.>. This process possibly leads to very strong flaring events observed in X-rays <cit.>. It is thought to be caused by magnetic reconnection events. Magnetic reconnection transforms a fraction of the system magnetic energy into particle kinetic energy, producing bursts of high-energy particles and radiation as observed in our Sun <cit.>. T Tauri flares appear to be thousands of times more energetic than solar flares <cit.>. They are expected to produce particle fluxes that are five orders of magnitude higher than those produced by the Sun <cit.>. The high-energy particle flux produced by T Tauri stars strongly impacts the surrounding environment <cit.>. In particular, these particles can ionise the gas in the surrounding disc <cit.>.
Disc ionisation is a critical process in T Tauri discs, as it couples the disc material to magnetic fields, which affects the dynamics and chemistry of the whole disc <cit.> and jets <cit.>. Magnetic fields can also trigger instabilities in the gas that can drive accretion of disc material onto the star <cit.>. Besides, ionisation drives chemistry in the disc, leading to the formation of molecules and complex organic compounds like amino acids <cit.>, which are the building blocks of prebiotic chemistry. Ionisation also plays a role in the heating and cooling of discs <cit.>, a process known to be essential in the launching of jets <cit.>. Understanding the mechanisms that drive disc ionisation appears then necessary for understanding the formation and evolution of stars and the dynamics of accretion-ejection processes.

In a first article (, hereafter B23), we studied in stationary framework, the ionisation rate produced by magnetic reconnection events occurring above the inner disc. We call inner disc the region at a radial distance between R=0.08 au and 0.6 au from the central star. We estimated ionisation rates produced by protons, electrons, and secondary electrons accelerated from these flares. Our model indicates that in the close environment of the flare, ionisation by non-thermal particles is the dominant ionisation source in the inner disc of T Tauri stars. It produces ionisation rates several orders of magnitude higher than previously proposed sources such as X-rays, galactic Cosmic Rays (GCRs), and radionuclides. B23 also proposed a comparative analysis of parameters to determine the conditions under which ionisation rates produced by magnetic reconnection are dominant. The results suggest that energetic particles (EPs) accelerated above the disc during magnetic reconnection events enhances the ionisation rates in the inner disc region of T Tauri stars. They are a new source of ionisation important to be included in magneto-hydrodynamics (MHD) or astrochemical codes.


Studying ionisation in protoplanetary discs with a non-stationary approach allows to calculate the mean ionisation fraction in the disc. Ionisation is a source of heat and a precise computation is crucial for understanding how jets are launched from the disc. Such results will be of interest for MHD simulations studying jet launching to constrain the heating rate based on physical constrains. Another important aspect is to couple the ionisation models with non-ideal MHD simulations of protoplanetary discs. This will provide insights into how the magnetic fields interact with the material in the disc and how it may shape the system and affect its evolution. Another goal of this work is to improve the astrochemical codes of protoplanetary discs. By developing more accurate models of the ionisation processes at play in these systems, we shall have a better understanding of the disc chemical evolution leading to the formation planetesimals and the synthesisation of the building blocks of live. 

In this paper we study in detail the non-thermal particle injection in the star-disc environment. It is the first of a series in which we shall study ionisation of T Tauri discs due to flares in a non-stationary framework. In this paper, we focus on the effects of the flares occurring close to the star, in the inner disc. For this study, we benefit from the X-ray observations of Chandra presented by <cit.> and <cit.>.

The layout of the paper is as follows. In section <ref> we summarise the stationary model as presented in B23, then present the physical model of flares, the chemical model and the model to derive the viscous α parameter in the disc from the ionisation fraction. In section <ref>, we present our results for the case of one flare, then we extend them to the time-dependent solutions from our multiple flare model. We consider the impact of non-thermal particle ionisation over disc chemistry, the development of the magneto-rotational instability and the heat deposition in the disc and at the base of the jet. In section <ref>, we propose a general discussion about theoretical perspectives of our work and its limitations before concluding in Sec. <ref>.



§ METHOD


 §.§ Summary of the stationary model in B23


Elaborating on the work by <cit.> and <cit.>, B23 propose an uncovered source of ionisation in the inner discs of T Tauri stars, the ionisation by flares triggered by magnetic reconnection. The idea was also inspired by the results of the simulation of magnetospheric ejection by <cit.> where the typical magnetic field topology of magnetic reconnection can clearly be identified in the interaction zone of the magnetic field of the star and the magnetic field of the disc. On the other hand, it has been routinely observed in the Sun and in laboratory experiments that such magnetic reconnection events are very efficient at accelerating particles to energies ≲ 1 GeV. Given the geometry of the field lines above the disc, we expect that this geometry favor the propagation of EPs from the acceleration region down to the disc. 

Once these particles penetrate into the disc, they suffer from energy losses. At energies ≲ 1 GeV, particles mainly lose their energy through ionisation losses. In other words, they are very efficient at ionising the disc. This supports the idea that our model could be a significant source of ionisation in the inner disc. B23 estimate the ionisation rate produced by EPs accelerated by magnetic reconnection events in the region of star-disc magnetic field interaction and proposes a parametric study. The methodology of the B23 model is as follows. 

First, B23 estimate the flux of supra-thermal particles produced by a flare. This estimation is based on X-ray observations from the COUP (Chandra Orion Ultradeep Project) point source catalogue <cit.>. Using these data, B23 detail how to estimate the flare size and the electron density of the flare based on a the observed flare temperature. B23 then develop a model of particle acceleration which, allows to estimate the non-thermal particle flux from the temperature, thermal electron density and size associated to the flare. The acceleration model also relies on the results of numerical simulations of magnetic reconnection to estimate the power-law index of the flux of non-thermal particles. Thus, B23 propose an estimation of the non-thermal particle flux produced by flares of T Tauri star based on micro-physical considerations. This flux takes as parameter the flare temperature, the index of the EPs distribution and the maximum energy of the EPs distribution. 

Subsequently, B23 describe their model for particle transport. This transport model allows for the estimation of the propagated particle flux within the disc based on the initial particle flux produced by magnetic reconnection, as a function of column density. Particle transport is described within the Continuous Slowing Down Approximation (CSDA), which is applicable under the assumption that the pitch angle of particles remains constant and that the energy losses due to collision with the medium are small compared to the kinetic energy of the particle. B23 computed the propagated flux for column densities ranging from 10^19 cm^-2 to 10^25 cm^-2 in order to treat the propagation of particles in the CSDA regime <cit.>. Within this approximation, given that the gyroradius of particles with energy E<1 GeV is small compared to the characteristic size of the disc, the trajectories of the particles are determined by the large-scale magnetic field lines. B23 provides a detailed description of the energy loss processes using the formalism of the energy loss function. The loss function is defined as the energy loss per unit column density and depends both on the incoming particle species (protons or electrons) and the composition of the target medium. 

The composition of the target medium, i.e. the abundance of various chemical species within the disc through which the particles travel, is determined using the radiation thermo-chemical code (see Sect. <ref> for a further description). B23 focuses exclusively on the losses associated with the dominant species: free electrons, H, H_2 and He. We have estimated that the energetic losses due to more complex species accounts for less than 5% of the total energy losses. Consequently, we neglect the contribution of metals and more complex molecules to the energetic losses in order to simplify the transport model. 

In B23, the electron loss function is deduced from the data provided by <cit.> at low energy (<1keV), which are based on theoretical calculations. At high energy (>1keV) B23 use the National Institute of Standards and Technology database[<https://physics.nist.gov/PhysRefData/Star/Text/ESTAR.html>], see <cit.> for an updated description of the electron energy loss function. At low energy, energy losses in molecular gases are controlled by rotational and vibrational excitations, at intermediate energy by ionisation and at high energy by radiative losses. The data of the Stopping and Range of Ion in Matter <cit.> are used to construct the proton loss functions. At high energies, above the threshold E^π = 280  MeV, we add
energy losses due to pion production, as reported by <cit.>. At very low energy (<1-10  eV) Coulomb losses are dominant, we use the expression from <cit.>. 

Whether the cross sections and loss functions are determined experimentally or theoretically, the uncertainties in their determination by the authors cited in B23 are negligible compared to the uncertainties linked to the parameters of the acceleration model. Additionally, B23 include the flux of secondary electrons produced issued from the ionisation of the medium by primary particles <cit.>. 

Finally, B23 show that, in a stationary configuration, EPs are a powerful source of local ionisation, with ionisation rates that exceed X-ray, stellar EPs, and radioactivity contributions in the inner disc by several orders of magnitude. For the reference case of a 1MK flare at 0.1 au propagating along a vertical magnetic field line, B23 find an ionisation rate ζ=10^-9   s^-1 at a column densities of 10^25  cm^-2, while the ionisation rate at this depth due to GCR or X-rays produced in a 1  MK stellar flare is of the order of 10^-17   s^-1, while it is of the order of 10^-18   s^-1 due to radionuclides. B23 then proceed to a parametric study where they consider various spectral indices for the EPs flux, various flare temperatures, different poloidal magnetic configurations as well as toroidal components along which the EPs can propagate, and different distance from the star where EPs are entering into the disc.

Although in conclusion B23 indicate that their assumptions may lead to overestimate the ionisation rate even if this process is shown to be a potentially dominant ionisation processes in the inner disc of T Tauri stars. As there are several parameters in the model that are difficult to constrain either experimentally or observationally, B23 conducted a comparative analysis of these parameters. The aim of this analysis is to define a range of flare parameters so that the ionisation rates produced are dominant over other ionisation sources. B23 anticipate that this will be the case for: (i) a reconnection process that accelerates particles following an injection flux with a power law  in kinetic energy j ∝ E^-δ for δ < 6, (ii) flares with temperatures above 1 MK, (iii) particles propagating along the field line with a ratio of the toroidal component to the poloidal component b_ g=B_ϕ/B_ pol < 1.

In the present article, we have updated the method for calculating the non-thermal particle density based on temperature of B23 (see App. <ref>). This revised approach is applied to determine the injected flux of non-thermal particles. We note that the ionisation rate is affected only by a factor of 2, which is relatively modest regarding the lack of constraints on the injected flux. This adjustment has insignificant implications for the chemistry involved.

According to B23, the temperature strongly affects the ionisation rate in the disc as it determines both the normalisation and the injection energy of the non-thermal component of the particle flux. Figure <ref> shows the ionisation rates generated by flares with temperatures ranging from 10^5 K to 10^9 K in the fiducial configuration of B23. This configuration assumes a power-law for the energy distribution function F(E) ∝ E^-δ, and δ=3 for particles produced by a magnetic reconnection event. It also assumes that magnetic field lines are perpendicular to the disc plane and that the EPs follow helical trajectories around the vertical lines with a homogeneous distribution of pitch angles. 

The assumption that the magnetic field is perpendicular to the disc plane stems from the results from simulations of <cit.> where the authors take into account turbulence and dynamo. They show that, especially when turbulence is saturated, the magnetic field is mostly vertical. Simulations in <cit.> of turbulent flows show similar results with almost vertical magnetic field lines. These results are comparable to what is obtained in self-similar accretion-ejection models <cit.>. However, self-similar models including magnetorotational instability <cit.> may show more complicated bent magnetic structures such as those we explored in B23. This however does not affect strongly the results. Of course, as the disc is axisymmetric, there should be an extra toroidal magnetic field that is acting as a guide field.
.






 §.§ Non-stationary flare model

As stated above, the model in B23 is stationary, the model estimates the ionisation rate below the flare, at its maximum luminosity. Consequently, the resulting ionisation rate is not representative of the ionisation of the whole disc on an extended period of time. This is why a time-dependent model should be considered. 

In this section, we propose a stochastic model to simulate the X-ray emission of a T Tauri star. Our aim is to calculate the resulting ionisation rate due to flares using the X-ray luminosity as an input. Therefore, it is needed to establish a relationship between the peak X-ray luminosity (L_X, pk)  and the corresponding flare temperature T_f.

To proceed, we rely on observations from Chandra. The work of <cit.> includes a table that displays the X-ray luminosity and the temperature of 216 flares in the COUP sample. Linear regressions are applied to the data. First, we deduced from this the relationship log L_X, pk - log T_ obs linking the peak X-ray luminosity in the Chandra detection energy band to the temperature deduced in this band T_ obs. Second, we deduce the relationship log T_f - log T_ obs, linking the bolometric and X-ray flare temperature. <cit.> indicate that because the peak X-ray emission of flares with temperature T_ obs>200 MK is outside the Chandra sensitivity band, the measurements deduced for these flares are highly imprecise. So, to perform the regression, we only consider a restricted sample of 179 flares with T_ obs<200 MK. 


In Fig. <ref>, we show the scattering plots of Chandra data points along with their corresponding linear regression. In Fig. <ref>, we plot the relation between the peak X-ray luminosity and the observed temperature of the flare L_X, pk( erg  s^-1)=9.63×10^29(T_ obs/ MK)^0.75. In Fig. <ref>, we plot the relation between the observed temperature in the Chandra band and the bolometric temperature of the flare T_f=3.35× T_ obs^0.55 (both temperatures are in MK units). Combining these two equations, the relation between the observed X-ray luminosity and the temperature of the flare can be derived as,

    T_f= 3.45 ×(L_X, pk/10^30  erg  s^-1)^0.73   MK.



We present now the stochastic model used for probing the parameter space relevant for T Tauri flares. Our main objective is to simulate the X-ray emission from young stars during a Keplerian period P_K. For solar mass stars, the Keplerian period at 0.1 au, P_K, is approximately 10 Ms (≈  0.3 years).

Our methodology is based on a random sampling of flare distributions, based on three parameters, the waiting time between two events, the X-ray luminosity and the flare duration.
To randomly sample the distributions, we use an inverse function technique, as detailed in <cit.>.

A comprehensive observational analysis spanning six decades of solar activity, recorded by instruments such as ISEE-3/ICE, HXRBS/SMM, WATCH/GRANAT, BATSE/CGRO, and RHESSI <cit.>, enable the identification of the distribution of the waiting times between solar flare events. We assume that we can extrapolate this solar flare temporal distribution to T Tauri star flares. This empirical distribution helps estimating the flare frequency within our designated period of study, while also giving a time of occurrence to each flare event. By relying on solar observations, we can estimate the frequency distribution of flares with luminosity that are lower than those observable on T Tauri stars.

Subsequently, to ascribe a luminosity value to every individual flare, we refer to the luminosity distribution of X-ray flares from young stars, as observed by Chandra <cit.>. In addition, the data deduced from the X-ray flares of these young stellar objects provide the distribution of the flares duration.



  §.§.§ Waiting Times Distribution


Observations of waiting times between solar flares reveal power-law distributions. <cit.> and <cit.> propose a type II Pareto distribution based on data collected over 60 years,

    n_ wt(τ)=λ_0(1+λ_0 τ)^-m   ,


where τ is the time between two flares, 1/λ_0 is the breaking point of the power-law.

We use the value 1/λ_0=0.80 ± 0.14 hour, as suggested by the authors. m is the power-law index of the distribution. Solar observations show that m ranges from 2.1 to 2.4. The average time between two events decreases as the power-law index increases. Since T Tauri stars are expected to be more magnetically active than the Sun, we take the upper bound of the solar observed waiting time power-law index, m=2.4. This gives shorter waiting times. 

From this distribution, we compute the mean waiting time between two flaring events,

    τ̅=∫_0^∞τ n_ wt(τ)  d τ≈ 5   ks.




  §.§.§ Luminosity Distribution


The exhaustive study by <cit.> focused on the energetics of X-ray flares from Pre-Main Sequence (PMS) stars. Their observations revealed that the X-ray peak luminosity distribution of high luminosity flares (L>10^32  erg s^-1) follows a power-law with an index of 2.11. For luminosities below this threshold, the trend of the power-law cannot be deduced because the set of data is incomplete. However, this index of 2.11 corresponds to the value of the power-law index observed over a wide range of luminosities, from optical to extreme-UV (EUV) and X-ray domains, in solar and stellar flares. The range of validity of this power-law starts from solar nanoflare energies (E_ flare=10^24 erg) up to superflare energies of solar-type stars (E_ flare=10^35 erg) <cit.>. This justifies our extrapolation of the power-law distribution observed by <cit.> to luminosity lower than 10^32  erg s^-1.

The normalised X-ray luminosity distribution can be written as,

    n_X (L) = a-1/L_c(L/L_c)^-a,

where a=2.11 is the power-law index measured by <cit.> and L_c is the low luminosity cut-off of the distribution. L_c is adjusted in order to reproduce the occurrence rate N_ MF of mega flares, flares with peak X-ray luminosity L_X, pk > L_0=10^32.5, observed by <cit.>, which is N_ MF= 1.7 mega-flares per year per star.

In the previous paragraph we have calculated the mean waiting time between two flares τ̅=5   ks. We can estimate the total number of flares per year, N_ tot = 6300  flares/year/star.

To find the low luminosity cut-off we start from the equation,

    N_ tot= N_ MF/P(L_X, pk>L_0),

where P(L_X, pk>L_0) is the probability of occurrence of a megaflare, i.e.,

    P(L_X, pk>L_0)=∫_L_0^∞ n_X(L) dL= (L_0/L_c)^1-a, a > 1
        .

From the two above equations, we get,

    L_c= L_0 (N_ MF/N_ tot)^1/a-1.


Taking N_ MF=1.7, N_ tot=6300, a=2.11, L_0=10^32.5   erg  s^-1, the low-luminosity cut-off is,

    L_c ≃ 1.9× 10^29   erg  s^-1 .




  §.§.§ Flares luminosity profile


We consider a simple time evolution model of flare soft X-ray luminosity L_X(t) as typically observed in the Sun <cit.>. It is characterised by an exponential fast growth and gradual decay. There are three parameters in this model. The first parameter is the characteristic rise time (τ_r), the second,the characteristic decay time (τ_d) and, the last, the peak X-ray luminosity of the flare (L_X, pk). The temporal evolution of a flare with luminosity peaking at t=t_0 is given by,

    L_X(t) = {[  L_X, pk e^(t-t_0)/τ_r                  t<t_0; L_X, pk e^-(t-t_0)/τ_d                        ].


The peak luminosity is randomly chosen according to the distribution in the previous section. 
Figure 6 in <cit.> shows no convincing correlation between the luminosity of the flares and their decay time. Thus, we derive a growth time following a Poisson distribution of the form,

    n(τ_r)=λexp(-λτ_r),

where 1/λ=26 ks is the average rise time of the flares from the observed sample by <cit.>. 

Figure 6 in <cit.> shows a correlation between flare growth and decay times (both in ks units),

    τ_d=5.28 τ_r^0.60



From the solar flare waiting time distributions, the peak X-ray luminosity distribution and the rise and decay time distribution of the luminosity, we can model the total X-ray luminosity emitted by a PMS. 

In solar flares, the fast rise followed by a gradual decay of the luminosity describes the soft X-ray component emitted by the thermal plasma. Conversely, the non-thermal particles, accelerated during the flare, produce hard X-rays via Bremsstrahlung and such hard X-rays are only observed during the rising phase of the soft X-ray emission <cit.>. 

The time-dependent ionisation rate, represented by ζ(t), reflects the temporal profile of these non-thermal particles. Besides, the temporal profile of non-thermal particles follows the profile of hard X-rays. Consequently, we assume that the profile of the ionisation rate matches the phase of increasing soft X-ray luminosity.


    ζ(t,L_X, pk) = {[ ζ_ pk(L_X, pk) e^(t-t_0)/τ_r                        t<t_0;                            0                              ].

where ζ_ pk is the ionisation rate at the peak of the luminosity that is calculated using the method described in B23, and t_0 is the time at the luminosity peak. 

The waiting time, the rising time and the decaying time of the flares are of the same order of magnitude. Thus, a luminosity continuum appears. Our model produces a mean luminosity over a year of 4.0× 10^30 erg s^-1, which is close to the typical value of 2× 10^30 erg s^-1 of the X-ray continuum observed on PMS <cit.>.



  §.§.§ Flare geometry Model

Flares, characterised by sudden energy discharges within a tenuous plasma confined to a "magnetic bottle",
have been extensively explored in various works <cit.>. The energy dissipation in these scenarios predominantly occurs via optically thin radiation and efficient heat conduction directed towards the chromosphere of the star <cit.>. Magnetic confinement plays a crucial role in determining the typical flare light curve. It yields a fast rise in soft X-ray emissions due to the rapid plasma heating, followed by a gradual, near-exponential decay resulting from combined effects of thermal conduction and radiation.

This temporal evolution follows a very general trend that the majority of observed flares satisfy. It has turned out to be a key to deduce the half-length (L_f) of the magnetic structure that governs the eruption <cit.>. <cit.> proposed an improvement to the original <cit.> relationship, expressed as follows,


    L_f=τ_d T^1/2/γ F(ϵ) ,


where γ≃ 3.7 10^-4 is generally valid for solar-type stars and the empirical function F(ϵ) = 0.63/ϵ-0.32 +1.41 > 1. The introduction of the function F(ϵ) includes the effect of a progressive heating in the rise phase of the flare whereas <cit.> assumed an impulsive heating with F(ϵ)=1. The parameter ϵ depends on the star. It varies for a given star from one flare to the other. <cit.> derived the mean value of this parameter for a sample of flares in Class I, II and III young stellar objects. In this study we are interested in Class II YSO and for this class ⟨ϵ⟩  = 1.43 and ⟨ F(ϵ) ⟩≃ 1.98, we hence assume F(ϵ) = 2 hereafter.

The half-length of the flare loop is a proxy for the position R where the flux tube is in contact with the disc, thus where the EPs are released in the disc.

We will discuss two models for the estimation of the position R. The first model that we call the Star-Star flare (hereafter SSF) model, is supported by X-rays observations of <cit.> where the two footpoints of the flare are anchored onto the central star and the flare loop-top is in the disc, see Fig <ref>. <cit.> estimated the radio emission of a flare in this configuration.


The second model, that we call Star-Disc flare (hereafter SDF) model, is supported by simulations <cit.> and by the observations of <cit.> with Chandra and <cit.> with XMM-Newton and Suzaku observatories. This model assumes that one flare loop is anchored to the star while the other is anchored to the disc, and the loop-top, where reconnection occurs, is above the disc, in the magneto-ejection region, see Fig <ref>. 

The main impact that each of these models has on our results is the position at which the EPs produced in the reconnection region enters in the disc. Assuming the first model (SSF), with the two foot points anchored into the star, the particles are penetrating at R≈ L_f. On the other hand, assuming the second model (SDF) with the flare loop linking the stellar surface with the disc, particles are penetrating at R≈ 2 L_f. The second model injects particles farther in the disc. We show that this has a major impact on the ionisation spatial distribution in the inner disc, thus potentially to the extend of the magnetorotational instability (MRI) active region.


From Eq. (<ref>), we can assign to each flare a half loop length given its peak X-ray luminosity L_X, pk and its decay time τ_d. These values are picked from the distributions obtained in Eq.(<ref>) and Eq. (<ref>) and allow to derive a radial position R.

<cit.> also made an estimation of χ, the ratio between the flare loop half-length and its cross-sectional radius r. So the area of the flare cross section A_f is expressed in terms of χ as, 

    A_f= πχ^2 L_f^2.

In solar and stellar flares as well as in the PMS flares observed by <cit.>, the typical value of χ is 0.1. 

Each flare is given a radial location in the disc based on its luminosity, while the azimuthal location is randomly picked in an uniform distribution of angles.




 §.§ Chemical Model



  §.§.§  
 
To evaluate the column density explored by EPs, we use a model generated by the radiation thermochemical code ProDiMo <cit.>. This model performs a wavelength-dependent radiative transfer calculation, including X-rays  <cit.>, which determines the gas and dust temperatures and the local radiation field. Chemical abundances are used to balance the heating and cooling processes, and the chemical network contains 235 different chemical species and 3143 chemical reactions <cit.>.
Here, we use a fiducial T Tauri disc model described in <cit.>. It was also used in <cit.> (see their  for the main physical parameters). Here, we want to model in a more realistic way, the thermal and ionisation structure in the innermost region of the disc (R< 1 au). Thus, we additionally include accretion heating following , see also  Appendix A). We use a total mass accretion rate onto the star of Ṁ_accr=10^-8 M_⊙ yr^-1. We added a simple treatment of the collisional ionisation to the chemical network following <cit.>. For further details, see Appendix <ref>.




  §.§.§ Electron density equilibrium

In the framework of EP propagation within the protoplanetary disc, we adopt a model based on the "test-particle" approximation. We assume that these particles exert no influence on the disc chemical structure. Nonetheless, our objective is to quantify their effect on the chemistry of the disc. To this end, we determine the electron density in regions of the disc influenced by flares.

From the differential system Eq. 7 by <cit.>, and assuming global neutrality, where the electron density equals the cation density (n_e=n_m^+) while setting aside the impact of metallicity, we can simplify the system down to a single ordinary differential equation,

    d n_e/d t =ζ n_n - β n_e^2.


In the no-flare case, assuming steady state, the initial electron density and neutral species density, noted respectively n_e,0 and n_n,0, are computed by the code. Upholding particle number conservation, the evolution of the neutral species density is related to the electron density and the initial conditions as,

    n_n(t)=n_e,0-n_e(t)+n_n,0.

In Eq.(<ref>), the parameter β=β(T) is the recombination rate from Eq. 9 of <cit.>. We note ζ, the total ionisation rate, being the sum of two components, ζ=ζ_f+ζ_ Prod, the stationary ionisation rate, ζ_ Prod, computed by in the scenario without flare, ensuring a stationary n_e,0. On the other hand, ζ_f, the ionisation rate due to flares, is derived from Eq. (<ref>). 

A global study of the EP feedback on the disc chemistry via the chemical model of will be studied in Sect. <ref>. Nevertheless, it is interesting to compute analytically the electron density via the above-outlined rudimentary chemical model. Notably, temporal dynamics remain elusive within the framework. This simple model provides insights into the temporal trajectory of electron density modulations induced by flares, an aspect presented in Sect. <ref>.



 §.§ Viscosity model

As already stated, in order to develop, MRI requires a sufficiently high degree of ionisation for the gas to interact with the magnetic field. The instability can be characterised by the viscous parameter α. 

<cit.> pointed out that Hall diffusion could, under certain circumstances, counteract the resistive damping effects of MRI. However, numerical simulations indicate a more complex scenario <cit.>. Given these uncertainties, <cit.> chose to exclude the effects of Hall diffusivity in estimating the effective viscosity parameter. To compute the viscous α parameter, we rely on the model proposed by <cit.>, which introduces an empirical non-ideal MHD MRI-driven equation for α suitable for physico-chemical protoplanetary disc codes such as . The adopted expression for α is given by 

    α= (2/β_ mag)^1/2 min(1,Λ_Ohm) [(50/Am^1.2)^2+(8/Am^0.8+1)^2]^-1/2
        ,


if √(β_mag)Λ_Ohm > 1 and α≈ 0 otherwise. 

Am represents the frequency at which neutral particles collides with ions normalised to the keplerian orbital frequency, 


    Am ≡ν_ in/Ω_K = β_ in n_ charge/Ω_K= x_e β_ in n_ tot/Ω_K ,


where n_ charge = x_En_ tot is the total number density of charged species and β_ in = 2 × 10^-9  cm^3 s^-1 is the collisional rate coefficient for ions to distribute their momentum to neutrals. Thus,

    Am≈ 1(x_e/10^-8)(n_n/10^10  cm^-3) (R/1 au)^3/2 .



The Ohmic diffusivity is quantified by the dimensionless Elsasser number Λ_Ohm, defined as the ratio of the Lorentz force to the Coriolis force,

    Λ_Ohm≡B_z^2/4πρη_O Ω_K≡v_A^2/η_O Ω_K≡( 4πσ_O/Ω_K) ( v_A/c)^2,

where B_z is the vertical component of the magnetic field, v_A is the Alfvèn speed and ρ is the mass density of the plasma and the Ohmic resistivity is given by η_O = c^2/4πσ_O. <cit.> derived an approximation for the Ohmic Elsasser number in disc regions relevant to our study where, σ_O ≈σ_e,O and x_e > 10^-13:

    Λ_Ohm≈ 1 ( T/100  K)^1/2( 10^4/β_mag) ( R/1  au )^3/2( x_e/10^-9),

The magnetic term β_mag(R, z) is the ratio of the thermal pressure to the magnetic pressure, given by <cit.>,

    β_mag(R,Z) = β_midP_th(R,Z)/P_th(R, 0) = β_midn(R,Z)T(R,Z)/n(R,0) T(R,0),

where β_mid= β_mag(R,0) is the value of β_mag on the midplane. β_mid is assumed to be independent of the radius so that it can be determined from the thermal structure and chemical abundances of . 
The parameter β_mid is typically in the range 10^4-10^6. We take β_mid=10^4 as reference value.





§ RESULTS




 §.§ Single flare effects



  §.§.§ Ionisation rates 




We plot in Fig. <ref> the solution of Eq. (<ref>) for the number density of electrons. The solution is obtained by taking the time-dependent ionisation rate produced by the particles propagating along vertical field lines emitted by a single 10 MK flare at R=0.1 au. The electron density is computed for column densities ranging from 10^19  cm^-2 to 10^25  cm^-2. It illustrates that for column densities below 10^24  cm^-2, the system rapidly returns to a steady state. In this range of column densities, due to the high electron density, the recombination time τ_ rec=(β(T)n_e)^-1∼ 1 s. Thus, the increase of the electron density only occurs during the time during which particles are injected. On the other hand, at N>10^24 cm^-2, despite the low temperature (≈ 500 K), the electron density is very low so the recombination times are longer, τ_ rec≳ 10^5 s, thus electrons remain free much longer. As a result, the effects of the flare persist for an extended period at high column density. Given that the equilibrium state discussed in this section is rather basic and solely accounts for the recombination of electrons with H^+, it is necessary to refine these results. This can be done especially by taking into account various other recombination processes, including those occurring on grains or metals.

In the next section, we study the effect of multiple flares using the full chemical model of accounting for all recombination processes.

In B23, we argued that the ionisation rates obtained in the stationary state were overestimated and in order to have more realistic description of the ionisation processes, temporal and spatial averaging should be performed. We propose here a simple estimation of the spatial and temporal average of the ionisation rate due to a single flare and postpone a more precise treatment to next section. We assume the same flare as above, 10 MK at 0.1 au, with a rising time τ_r= 26 ks corresponding to the average rising time observed by <cit.>. We define an effective time-averaged ionisation rate ζ_ eff from the average ionisation fraction during a period of time P, 

    ζ_P, eff=1/P⟨ x_e ⟩ _ P = 1/P√(%s/%s)⟨ζ⟩ _ Pβ n_n  ,

where β is the recombination rate, x_e, the ionisation fraction and ⟨·⟩_P is the time averaging over a period of time P. The temporal averaging of the ionisation rate produced by a single flare over a period P is,

    ⟨ζ_P ⟩ =1/P∫_0^Pζ(t) dt.


Given the temporal profile of Eq. (<ref>), the average ionisation rate produced by a flare occurring at t=t_0≪ P is, 

    ζ_P=ζ_ pk(L_X, pk) τ_r/P.


We compute ζ_P, eff over a Keplerian period at 0.1 au, so P≈ 10 Ms. We then define the spatially averaged ionisation rate as,

    ζ_ eff= r_f^2/(R+r_f)^2-R^2ζ_P, eff,

where R is the radial position of flare and r_f the radius of the flare. ζ_ eff/ζ_P, eff is the ratio between the area illuminated by particles and the disc ring of width 2 r_f at radius R. Here we take R=0.1 au and r_f/R=0.1. This effective ionisation rate is illustrative for comparison with other sources of ionisation in the inner disc. The effective ionisation rate is plotted in Fig. <ref> and compared to the X-ray ionisation rate computed by .
It was anticipated that calculating spatial and temporal averaged ionisation rates would result in lower values than the ionisation rate estimated at a flare maximum in B23. This hypothesis is validated here. Considering that observations span over extended periods of time, the averaged ionisation rates calculated in this study offer a better basis for comparison with observational data than the ionisation rate obtained in B23. Furthermore, these averaged ionisation rates also provide more relevant values to be used in numerical models of stationary chemistry in discs such as .
As seen in Fig. <ref>, the effective ionisation rate produced by flares at low column densities, N≲ 10^22  cm^-2, is similar to the ionisation rate produced by X-rays. However, at high column densities, N≳ 10^22 cm^-2, the effective ionisation rate produced by flares becomes one to two orders of magnitude higher than that produced by X-rays. The increase in ionisation rate by EPs compared to X-rays at high column densities is an expected result as has been demonstrated in <cit.>.






  §.§.§ Non-thermal pressure



In addition to altering the disc chemical state, non-thermal particles also influence its dynamics. To gauge this impact, we compare the pressure exerted by the non-thermal particle distribution to the magnetic and gas pressures. The non-thermal pressure, P_i,nt, exerted by a distribution of of electrons and protons (i=e,p), can be expressed as,

    P_i,nt = 1/3∫_E_c^∞ p(E) v(E) F_i,nt(E)   dE,

where E is the kinetic energy, E_c is the injection energy of the non-thermal particles, and F_i,nt(E) is the energy distribution of these non-thermal particles i,


    F_i,nt(E)= N_i,nt(E/E_c)^-δ.


 
We assume equipartition of energy between non-thermal protons and electrons, i.e F_e,nt = F_p,nt.


We consider both species with the same injection energy E_c. Thus, P_e,nt=P_p,nt=P_nt, with 

    P_nt≃2/3δ-1/δ-2 E_c n_ nt,

where n_ nt is the non-thermal particle density (Eq. <ref>) and δ 2. We assume that the injection energy is proportional to the thermal energy, E_c=θ E_ th, where E_ th=3/2 k_B T, k_B is the Boltzmann constant and T the flare temperature. Then, we express the pressure of the non-thermal distribution in terms of the pressure of the thermal distribution, 

    P_nt=3θ^2/δ-2√(3θ/2π)exp(-3θ/2) P_th.

The thermal pressure can be computed as function of the flare temperature from P_ th= n_ th(T) k_B T, where n_ th(T) is the plasma thermal density constrained from X-ray observation, Eq. (<ref>).  

Assuming the fiducial values of B23, θ=3 and δ=3,

    P_th= 26 (T/1  MK)^0.31  erg   cm^-3  and    P_nt= 9 (T/1  MK)^0.31  erg   cm^-3.


The total pressure of non-thermal particles accounting for both electrons and protons is P_tot,nt≈ 2 P_nt and the total thermal pressure is P_tot,th≈ 2 P_th.

In Fig. <ref>, we plot the total thermal and non-thermal particle pressure alongside the gas and the magnetic pressures. Above the inner disc surface, where a flare is occurring, the magnetic field intensity should be in the range  10-50 G <cit.>. The graph indicates that the pressure of thermal and supra-thermal particles generated by flares is typically comparable to the magnetic pressure at the inner disc surface. Meanwhile, the particle pressure greatly exceeds the disc gas pressure. Therefore, these particles are expected to influence the magnetic and gas structures of the discs. Earlier simulations examining the effect of pressure from flares focused on the pressure contribution from the thermal particles <cit.>. They assumed a magnetic intensity of 27 G in the inner disc, resulting in the magnetic pressure represented by the green dashed line in Fig. <ref>. <cit.> demonstrated that the pressure from the thermal component could greatly alter the thermal structure of the disc. It can even initiate overpressure waves that cross the disc, forming accretion funnels on the opposite side of the disc. 

From Eq. (<ref>), we have shown that the non-thermal component accounts for over a third of the thermal pressure. Thus, it represents a non-negligible portion of the particle pressure exerted by particles from a flare. It appears essential for future flare simulations in T Tauri stars to consider the effect of the supra-thermal component.



 §.§ Ionisation by multiple flares

A single flare appears to cause a temporary change in the disc chemistry, but the system comes back to its original chemical state over time. 

<cit.> observe that the collective effects the X-rays produced by flares can push some chemical species into a "new pseudo-steady state" if their formation or destruction timescales are longer than the intervals between flares. For instance, in their study, if the creation of a particular species is boosted by a byproduct of a flare and its decay time exceeds the time between strong flares, the system will not come back to the original state before the next flare occurs. The specific timescales for these changes vary depending on the molecule in question and the prevailing local conditions. <cit.> also argued that for flares, the chemistry of ions in the disc should be modelled using time-averaged ionisation rates. 

So, in order to have a more complete understanding of the effects of EPs produced by multiple flares, we compute the spatial and temporal averaged ionisation rate. We use a sample of flares that are randomly drawn using the method described in Sect. <ref>, over a time period P. 

First, we need to deal with the temporal averaging of the ionisation rate produced by a single flare over a period P, as in Eq. (<ref>).

Given the temporal profile of Eq. (<ref>), the average ionisation rate produced by a flare occurring at t=t_0 is, 

    ζ_P=ζ_ pk(L_X, pk) τ_r/P(1-e^-P-t_0/τ_r)


We now calculate the spatial average of the ionisation rate produced by the whole sample of flares. The spatial average is computed in disc rings with a radius R and width Δ R. The area of such a ring is A_R= π(R^2-(R-Δ R)^2). 

The flux tube of each flare i from the sample intersects the disc ring with an area A_i, where A_i is smaller than A_f. The flare cross section with the whole disc is given by Eq. (<ref>). 
Flares produce a time averaged ionisation rate ζ_P,i(N) in a volume of section A_i. 

The spatial average ionisation rate ζ_P,R in a disc ring with a large radius R is given by,

    ζ_P,R(N)=1/A_R∑_i A_i ζ_P,i(N).


We have calculated the average ionisation rate by varying the duration P. On shorter time scales, the likelihood of encountering a powerful flare capable of emitting particles at large distances from the central star is relatively low. Consequently, the primary variation introduced by changing the period is the radial extent of the region impacted by flares. 
The spatial average ionisation rate plotted in Fig. <ref> is obtained by taking Δ R ≪ R and averaging over 1000 years.



From the numerical results we propose a parametrisation of the average ionisation rate due to the EPs produced by flares of the form,

    ζ(R,N)= ζ_0 (1+N/N_0)^-a e^-b (R-R_0).

The parameters for the SSF and SDF models are presented in Tab. <ref>. This expression provides a fit matching the numerical results with an average relative deviation of less than 10% for column densities ranging from 10^19 to 10^25 cm^-2, for radii between 0.1-0.3 au and 0.1-0.6 au for SSF and for SDF, respectively.




 §.§ Feedback on the disc




  §.§.§ Effects on the disc chemistry

In this section, we explore how the ionisation rate, resulting from flaring activities, influences the chemical composition of the disc. This impact is quantified using ProDiMo, which takes as inputs the spatial distribution of the ionisation rate generated by both SSF and SDF represented in Fig. <ref>. ProDiMo then recalculates the entire disc chemistry, accounting for this new ionisation source. Our model computes ionisation rate for column densities N< 10^25 cm^-2. At higher column density the ionisation rate is extrapolated using linear interpolation in log space. The revised disc structures are displayed in Fig. <ref> where we specifically examine the ionisation fraction. In Fig. <ref>, we examine the distribution of the most abundant cation in the inner disc. On these figures, we plot the effects of three flare models on disc structures, with the horizontal axis indicating the disc radius R in astronomical units (au) and the vertical axis representing Z/R. We show the reference disc structure without flare contributions and those incorporating effects from SSF and SDF. 

In addition to the ionisation fraction, in Fig. <ref>, the solid black line shows the delimitation of the MRI inactive region, where α=0 (i.e. where √(β_ mag)Λ_ Ohm<1). We assume β_ mid=10^4 (see Eq. (<ref>)). The method to determine the MRI active region is explained in Sec. <ref>, and further discussed in Sec. <ref>.


The overall ionisation fraction within the disc is determined by the equilibrium between ionisation sources, namely UV, X-rays, cosmic rays, and EPs from flares, and recombination as well as charge exchange reactions. The principal charge carriers in this system are free electrons for negative charges and gas-phase cations for positive charges.

Panels (b) and (c) in Fig. <ref> reveal that, in disc structures influenced by flaring activity, the ionisation fraction is higher by more than one order of magnitude compared to structures that do not account for flares (panel (a)) at heights where the Z/R ratio is between 0.05 and 0.2.

Due to the longer loops of SDF compared to SSF, the impact of the former on the ionisation fraction extends farther radially, see  Fig. <ref>. Specifically, the region in which the ionisation fraction increases by more than one order of magnitude compared to the structure without flares, extends up to ∼ 0.6 au in structures that account for SDF. Conversely, it is limited to ∼ 0.3 au in structures accounting for SSF. The sharp profile observed in the ionisation fraction distribution at 0.3 au and 0.6 au in Fig. <ref> (b) and (c), respectively, can be attributed to the maximum distances from the star that particles produced by the flares can reach. 



Panel (a) in Fig. <ref> shows the distribution of the most abundant ion species in the inner disc without flares. At the disc surface, the ionising flux arising from the Far Ultraviolet (FUV) photons is predominant, and the gas layer essentially acts as a photo-dissociation region (PDR). In this region, atomic ions like H^+, C^+, S^+, and Mg^+ are the most abundant cations.

While in the UV and X-ray shielded inner disc, collisions due to high temperature ionise atoms having low ionisation potential like the alkali, producing a region where Na^+ is the dominant cation.

However, injection of non-thermal particles have additional effects on the chemistry. In the region undergoing the increased ionisation rate, at R < R_i, where R_i=0.3 au  and R_i=0.6 au for SSF and SDF, respectively, the chemical structure is similarly modified for both flare models. The disc surface layer is no more a standard PDR, dominated by C^+ and  S^+, but the most abundant ion is H^+. However, a surprising feature appears. There is a HCNH^+ dominated layer, that is absent in the model without flare. 




An important aspect for the chemistry in context of variable flares are the timescales attached to chemical reactions. In case of the EP-induced ionisation the picture is twofold. The ionisation due to the particles does not affect only the upper UV-dominated layers, but also deeper regions of the disc (towards the disc midplane). These two regions can roughly be separated by the location of the A_V=1 layer (see e.g. Fig. <ref>). As a consequence the region affected by the particles also spans a huge range in chemical timescales. In the UV-dominated layers those timescales a very short (about 1 day), whereas towards the midplane they span from ≈ 10^3 to ≳ 10^6 years <cit.>.

In our model the HCNH^+- dominated layer crossed the A_V=1 layer and hence an interpretation of the chemistry in this region is complex. One of the consequences of the increase in the abundance of HCNH^+ is an enhancement of HCN abundances through HCNH^+ + e ⟶ HCN + H <cit.>. This process is most important in deeper layers (below the A_V=1 limit) of the disc with longer chemical timescales and is a consequence of the ion-molecular formation pathway as described in <cit.> which also effects other molecules such as CN. <cit.> pointed out that X-ray ionisation plays an important role for the formation of HCN but also C_2H_2 in UV shielded regions. However, compared to X-rays, an irradiation by EPs produce higher ionisation rates deeper in the disc and hence those formation pathways become even more important.

In Fig. <ref> we show the HCN column density for models with and without flares. This result is given by the steady-state chemistry model assuming averaged particle ionisation rates from the flares. As the chemical timescales are long in the deeper layers of the disc (which dominate the column density), we expect that this result does not depend strongly on the flare frequency, but more importantly on the period of time in which the disc suffers from enhanced particle irradiation due to flares. However, a more detailed analysis of the complex chemistry and the time-dependent behaviour of the disc is required to make more firm conclusions and will be treated in a forthcoming work.

 

Nevertheless, the enhancement of HCN is a promising outcome of our study and interesting in two ways. Firstly, the increase in HCN abundance due to flare suggests that this species may become a candidate tracer for comparing our model with observational data. Secondly, on the theoretical side, further hydrolysis of HCN may produce amino-acids. Indeed, several chemical mechanisms have been suggested for the abiotic formation of prebiotic molecules. For instance, amino acids are commonly believed to originate from such processes. An aldehyde interacts with HCN and NH_3, resulting either in an aminonitrile, which is then hydrolysed to produce the corresponding amino acid or in a reaction between an aldehyde, HCN, NH_3 and CO_2 that leads to the production of amino acids through hydrolysis <cit.>. Our model, accounting for the ionisation by flares, produces a strong increase in the HCN abundance in region at distance R= 0.1 - 1 au from the central star. It may be of interest for the synthesis of the "building blocks of life" on primitive Earth-like systems.



  §.§.§ Effects on MRI and Accretion Rates

Figure <ref> displays the distribution of the effective viscosity parameter α for a β_mid value of 10^4. This parameter traces the turbulence efficiency, see Sec. <ref>. Three zones can be distinguished, the MRI-active region with high values of α, the Ohmic diffusion-limited "dead zone" delimited by the total damping criterion β_mag^1/2Λ_Ohm < 1, and the portion of the disc where β_mag < 1 in the disc atmospheres leading to a drop of α<0.01. The decrease in α, as one moves outward the disc, is linked to the reduction in gas pressure and consequently, β_mag. 


In regions with higher ionisation rates, the coupling between the gas and the magnetic field is stronger, leading to more efficient MRI-driven turbulence and angular momentum transport, see Sec. <ref> and <cit.>. 

We focus now on the effects of additional ionisation due to flares. In the region undergoing an increased ionisation rate, at R < R_i, where R_i=0.3 au and R_i=0.6 au for SSF and SDF, respectively, the values of the viscosity parameters increase from 0 to a finite value ≳ 0.01, activating MRI. On Fig. <ref> panels (b) and (c), we see that the size of the MRI active region increase from Z/R≈0.1 to Z/R≈ 0.06. In this interval, the MRI active region size increases by more than one order of magnitude in terms of column density. 

By controlling the disc viscosity, ionisation processes regulate the mass accretion rate onto the central star. If the mass accretion is driven by turbulence originating from either hydrodynamical or magnetohydrodynamical processes, assuming efficient angular momentum transport via the MRI, the disc mass accretion rate can be approximated in terms of an α-disc model <cit.>,


    Ṁ   =  3πνΣ = 2× 3πμ m_p N_g α c_s^2 
        = 2.5 × 10^-9(α/0.01) (N_g/10^24  cm^-2) (T/10^3  K)    M_⊙/yr,


where the effective viscosity is given by ν = α c_s h and μ is the mean molecular weight. c_s and h are the sound speed and the height scale of the disc, respectively, defined in Sec. <ref>. The factor of 2 accounts for accretion along the upper and lower disc surfaces, N_g is the gas column density above the midplane and T is the disc temperature. 
This formula is an approximation in our case as it assumes very efficient angular momentum transport, namely l_ in≪ l where l=R^2Ω is the specific angular momentum. This is correct far from the inner disc radius, but it may not necessarily be the case here. Nevertheless, the formula provides an order of magnitude estimation, which requires to be further tested using numerical simulations, like in <cit.>. In the region undergoing the increased ionisation rate, i.e. at R < R_i, where R_i=0.3 au  and R_i=0.6 au for SSF and SDF, respectively, over 1000 yr, the averaged accretion rate increases by a factor of 7 to 10 compared to the model without flares, see Fig. <ref>. 





  §.§.§ Disc heating by EPs 


Ionisation processes in T Tauri discs contribute to the heating of the disc material, which affects the disc thermal structure <cit.>. Various ionisation mechanisms, including cosmic rays, X-rays, and UV radiation, deposit energy into the disc gas, which in turn raises the gas temperature and modifies the thermal equilibrium <cit.>.
In this section, we compute the volumetric heating rate due to the ionisation of EPs produced by flares. 

When they are separated from their nucleus, electrons deposit heat in the surrounding medium. 
The amount of energy deposited per unit time and unit volume is called the volumetric heating rate. The volumetric heating rate due to EPs, Γ is computed assuming that the ionisation by each atomic hydrogen releases an amount of energy Q_ H=4.3 eV and ionisation by each molecular hydrogen releases Q_ H_2=18 eV <cit.>. Q_ H and Q_ H_2 are the sum of the heating produced by elastic collisions, rotational excitation, excitation of H_2 vibrational levels, dissociation of H_2, and chemical heating produced by the reactions of the H^+ and H_2^+ ions with other chemical species in the disc. All these heating processes are detailed in <cit.>. So, 

    Γ(R,N)= (n_ H (R,N) Q_ H+ n_ H_2 (R,N) Q_ H_2) ζ(R,N),

where ζ is the ionisation rate per H nuclei.
The volumetric heating rate in the disc as calculated in is shown in Fig.<ref>.




Previous works have indicated that heating the surface of the disc can lead to an increase in the mass of the outflow <cit.>. <cit.> showed that introducing a heating source comparable to emissions from an embedded stellar object could offset magnetic compression, thereby allowing for lower mass loads. Previous works simulating the launching of jets incorporate heating rates at the disc surface. The primary sources of disc ionisation previously considered were high-energy stellar X-ray and UV radiation, which also act as the predominant heating agents for both, the disc surface and its entire wind zone. In the absence of magnetic fields, this level of heating is powerful enough to initiate the photoevaporation of protoplanetary discs. This process is widely considered as the leading mechanism for disc mass loss and dissipation, as highlighted in the review by <cit.>.

The intensity of this volumetric heating rate influences in turn the jet and the wind characteristics, such as the degree of collimation, velocities, mass loss rates, and temperatures, see Meskini et al. in prep.
In these studies however, the microphysical origin of the volumetric heating rate is not discussed. In the context of the recent advances in the development of thermal-MHD outflow models <cit.>, we derive a heat source based on microphysical consideration. 
We compute Γ_S, the average volumetric heating rate in the surface layer, at column densities between 10^-19 and 10^-20 cm^-2. Figure <ref> shows the trends of Γ_S, which is parameterised by the power-law,




    Γ_S≈ a × 10^-6(R/0.1   au)^-b  erg   cm^-3 s^-1,

where a= 3 and b=6.3 for the SSF model, while a=8 and b=5.2 for the SDF model.


The volumetric heating rate computed here can be considered as an upper limit for the volumetric heating rate in outflows. In outflows, the density is lower than the density in the disc. This would lead to lower heating rates in the outflows, see Eq. (<ref>). Thermal MHD simulations show that strong heating rates as the one computed here lead to uncollimated outflows, see Meskini et al. (in prep.).

While MHD winds are generally considered to drive stellar accretion through mass loss, thermal winds have been suggested as a mechanism for dispersing protoplanetary discs. Specifically, these thermal winds, also known as photoevaporative winds, are driven purely by thermal forces and do not contribute to the angular momentum transfer. They were initially proposed by <cit.> to explain the relatively short lifetimes of discs, of the order of a few million years. In these models, the gas at the surface is heated by photons and forms a wind that is thermally energised. The rotating gas has enough energy to escape from the stellar gravitational potential. Although the launching radius of a thermal wind for a T Tauri star is around ≈ 1 au <cit.> there is also a non-zero flow inside the launching radius <cit.>
 
We show here that the ionising particles produced by flares are also a source of heating that should be taken into account in the thermal equilibrium of the disc for the launching of magneto-thermal winds <cit.>.




§ DISCUSSION




 §.§ Viscous heating

The thermal equilibrium in discs should not only integrate chemical and radiative heating and cooling processes but also account for viscous accretion heating Γ_ acc. For a thin disc, this last quantity can be expressed as <cit.>,

    Γ_ acc= 9/4αρ c_s^2,

where ρ is the gas mass density and c_s, the sound speed. In this framework, the parameter α links the viscous heating to the gas pressure (P_g=ρ c_s^2). This equation is fundamentally phenomenological and comes with substantial uncertainties. We use it here just to build an illustrative idea of the physical processes at play. 

Disc models that include such a heating source are particularly dependent on the value of α for determining thermal structure. Former studies, such as those by <cit.>, have investigated the impact of different values of α on the thermal properties of the gas, commonly based on models with α≈ 0.01. In the aforementioned models, X-ray heating dominates in the disc corona, while accretion heating is more important in the inner disc regions.

In our study, we propose that EPs from flares could serve as a significant source of heating, both at the surface and within the inner disc, primarily emanating from two mechanisms. Firstly, heating is released through the ionisation of the disc atoms and molecules (see Fig. <ref>). Secondly, these EPs alter the chemical and ionisation structure, thereby increasing the disc viscosity (see Fig. <ref>). This increased viscosity triggers an increase of the accretion rate, leading to enhanced accretion heating (see Eq. (<ref>)).

These two additional sources of heating could, in turn, alter the disc thermal structure. We will quantify these changes caused by the additional heating sources in an upcoming paper. Meanwhile, we speculate that the temperature increase due to heating by EPs produced by flares may lead to an expansion of the region in which the disc temperature is high enough to maintain sufficient ionisation levels for triggering the MRI.



 §.§ Variability of accretion-ejection processes


Accretion variability depends on the physical conditions within the disc and the processes that generate unstable flows. Additionally, the amplitude of this variability depends on how efficiently the involved processes can extract energy and angular momentum from the disc. 

For events occurring at a disc radius R and influenced by gravitational forces, the local Keplerian time is given by,

    t_K∼ 1 (M_*/M⊙)^-1/2(R/1   au )^3/2 yr.

It ranges from just a few days near the inner edge of the disc, to centuries for material in the outer disc. However, if the diffusion of the disc material is the determinant process as in viscous α-disc models for mass transport, then the diffusion time, t_diff =R^2/ν, is,

    t_diff (R) ∼ 10^4 (0.01/α) (0.1/h/R)^2 t_K(R).

It sets the time scale. Here, h represents the disc scale height, where typically α < 0.01 and h/R < 0.1. Consequently, even in the innermost parts of the disc, viscous diffusion is expected to require a minimum timescale of several years.

Moreover, the thermal timescale defined from the sound speed c_s = α^-1 v_ turb falls between the dynamical and diffusion times, 

    t_th(R)∼ 10^2 (0.01/α) t_K(R)

and corresponds to months near the disc inner edge. 

Considering these two characteristic time scales, we infer that the increase in the accretion rate due to flares can be triggered by a corresponding average increase in ionisation over periods larger or similar to t_diff. Consequently, viscous accretion likely does not play a role in accretion variability over periods shorter than a hundred years in the inner disc. 

On the other hand, flare activity can impact the thermal structure on shorter timescales. Accretion triggered by the heating pulse of a flare, as simulated by <cit.>, may play a role on the accretion variability on timescale of months, which is the typical timescale of superflares as observed by <cit.>.

In addition, as the modification of the thermal structure also affects the launching of winds and jets, the variability in flare activity might play a role in the variability of outflows at timescales of months.



 §.§ Effects of EP in jets and winds launching

In the current understanding, the chemical state and the
thermal energy equilibrium in wind regions is determined by the local flux of stellar high energy UV and X-rays photons. Determining the gas temperature is crucial for accurately calculating pressure gradients, wind flow characteristics, and interpreting line emissions. The gas temperature is influenced by the interplay between various heating and cooling processes. As for the heating, factors like collisions with dust grains, and line emissions from ions, atoms, and molecules come into play. These molecular abundances are also dependent on underlying chemical reactions.

At the surface of the disc, where the wind originates, multiple factors contribute to the gas heating. This includes grain photoelectric heating by FUV photons. It also includes the conversion of binding energy to heating during the formation of molecular hydrogen. Finally, it includes the excess electron energy resulting from photo-dissociation and photo-ionisation by various forms of radiation, i.e. FUV, extreme Ultraviolet and X-rays. We have shown that in this surface layer, EPs produced by flares may be an important additional ionising and heating source.

Above the disc, magnetised jets and winds may also experience heating through non-ideal MHD processes like Alfvén wave damping heating, Ohmic heating and ambipolar diffusion heating, especially at altitudes where photo-heating is less effective, as shown by various studies (e.g., ). Based on the results presented in this article, which focuses on discs, we emphasise the importance of studying the effects of particles generated by flares also on winds and jets. Such a model will be presented in a forthcoming publication.






 §.§ Dependence of the results to the model parameters



  §.§.§ Waiting-Time distribution
 

In this section, we explore how the distribution of waiting intervals between two flares affects the ionisation rate. The aim of this analysis is to assess the possible variation in ionisation rates considering the uncertainty associated with the determination of the parameters of this distribution.
To randomly generate the luminosity of the flares, we have developed this distribution using three parameters: the total annual frequency of flares (N_ tot), the annual frequency of super flares (N_ MF), and the power-law index (a), see Sec. <ref>.
The estimation of N_ tot and a is based on the extrapolation of solar properties, which is a strong assumption and introduces significant uncertainty. Conversely, the estimation of N_ MF is supported by robust statistics thanks to observations of flares on young stars made by <cit.>, leading to short range of acceptable values.

We note that our model is primarily dependent on the frequency of super flares N_ MF. Because according to our flare geometry model, the loops of these high-luminosity flares are the only long enough to send particles into the disc. However, given the short range of acceptable values for this quantity, its variation has minimal impact on the results. On the other hand, although the uncertainty on the parameters N_ tot and a is larger, their limited influence on the ionisation rate means this uncertainty does not significantly affect the results. Varying N_ tot and a with N_ MF constant is equivalent to varying the flare frequency while keeping the energy deposited by the super flares constant.

In conclusion, varying the three parameters of the waiting time distribution within their uncertainty range does not significantly affect the average ionisation rate.



  §.§.§ Energetic particle density entering the disc

The most critical parameter in our model is the density of EPs penetrating the disc. This quantity is notably challenging to be constrained by observations. We have outlined a method in B23 based on the assumption that energetic electrons emit the observed X-ray radiation via bremsstrahlung. From these observations, we have estimated the density of the thermal plasma emitting X-rays. Assuming that this density is that of the reconnection region, we normalised the number of supra-thermal particles entering the disc. The uncertainty in our approach is that we are not certain that the density derived from bremsstrahlung radiation corresponds to that of the reconnection zone. It could actually be the density of the stellar or the disc chromosphere heated by the particles. If this is the case, the density that we assumed would be overestimated. We find that if the density of EPs entering the disc is reduced by a factor greater than about 3 orders of magnitude compared to that estimated in B23, the ionisation rate from particles produced by flares would be negligible compared to the ionisation rate by the star X-ray radiation. In their study, <cit.> have shown the challenges and requisite resolution for determining the density of supra-thermal electrons in solar flares through radio emission observations. This highlights the complexities involved and suggests that the observational determination of supra-thermal particle density in T Tauri flares is not currently feasible. In Sec. <ref>, we propose a numerical approach that could allow for a more accurate estimation of the density of supra-thermal particles penetrating the disc.



 §.§ Model limitations




  §.§.§ Particle model
 
First, there are aspects of the particle propagation model that we have not addressed, especially when it comes to propagation at column densities greater than 10^25cm^-2. Beyond this column density, the dominant losses for protons are through pion production. During pion production, protons lose a significant portion of their energy and change their pitch angle. Due to these factors, particle transport can no longer be approached using the CSDA framework <cit.>. Diffusive processes need to be considered. For a more detailed discussion, see <cit.>.

Second, in our particle model, we must establish what is the fraction of particles accelerated in the reconnection zone of the flare, that penetrate into the disc. This aspect relies on the flare geometry and the magnetic field structure within the star-disc interaction zone. For now, we will restrict ourselves to a parametric study. We scale down the incoming particle flux to estimate the minimum fraction of particles that must enter the disc to significantly impact its chemistry and dynamics. For a quantitative assessment, a test particle distribution would have to be integrated into a reconnection zone of an MHD simulation of the inner disc. This should allow to compute the fraction of particles that penetrate the disc, fall onto the stellar surface, and propagate into the outflows.

Third, we are also aware that the particle acceleration model from B23 is somewhat simple. However, the models and simulations of particle acceleration by magnetic reconnection are an active research field, see <cit.> for reviews. The magnetic structure of the star-disc interaction region does not only influence the location and proportion of particles entering the disc but also, as detailed in B23, the particle acceleration process itself. As shown in B23, the presence of a guide field, the magnetic field component that is not involved in the reconnection process, diminishes the particle acceleration efficiency. We have not accounted for a guide field here, placing us in an optimal configuration for particle acceleration. In a model with a more realistic magnetic structure, featuring a toroidal component that could act as a guide field, the acceleration process would be less effective, thus resulting in lower ionisation rates. 

While the literature discusses the possibility that superflares are sufficiently long to connect the star with the disc <cit.> or not <cit.>, we highlight the potential for particles to be produced by magnetic reconnection events resulting from Keplerian shearing. Indeed, <cit.> conducted simulations to model the occurrence of magnetic reconnection events within circumstellar discs, induced by Keplerian shear. The simulations show the generation of turbulence within the magnetic truncation region, leading to flaring phenomena. The resultant energy and luminosity distributions are characterised by power-law behaviours. This is what we postulate in our work based on observations of soft X-rays. This is now corroborated through numerical support. These distributions exhibit a universal character, consistent under various conditions such as the nature of the reconnection mechanisms and the boundary conditions implemented. Such universality leads <cit.> to draw analogies between magnetic reconnection processes in inner circumstellar discs and solar flares, a significant assumption of our study.

To summarise, the particle model used here provides an upper limit to the impact of particles on discs. A model that describes large-scale particle acceleration by magnetic reconnection in the specific magnetic environment of T Tauri stars is currently not available. Nevertheless, such a model would be crucial for a complete understanding of the effects of particles produced by magnetic reconnection on discs. For a more realistic picture, a first step could be to add particles in a MHD code in the test-particle limit or beyond, i.e. including a back reaction term.



  §.§.§ Disc model


A crucial point for interpreting our results is that all our conclusions about the impact of particles are specific to the chemical structure computed by . 
For a more complete picture, it would be interesting to examine the impact of the particles by varying the structure of discs. For instance, we anticipate that the particle impact would depend on the radial and vertical structure of the disc mass, and this dependence will be evaluated in a forthcoming paper.  

In addition, we would like to caution about the increase in the density of HCNH^+ and HCN due to the additional ionisation rate from particles. The code, although it considers hundreds of chemical species and thousands of reactions, does not take into account molecules with more than four atoms, except for methanol. For details, we refer to ). Therefore, the rise in the abundance of HCNH^+ might not be physical but rather a dead-end in the reaction network calculated by . The recombination of HCNH^+ would then lead to an increase in HCN, which might be overestimated. To assess the increase of the abundance of HCN, we would need to conduct a more detailed analysis of the chemical reactions leading to this increase. Subsequently, we could assess the influence of particles from flares in a chemical model that considers more complex molecules with additional atoms. Nevertheless, we can highlight that our conclusions on the increasing complexity of molecules remains valid. Even if the rise in HCN abundance might be due to a chemical reaction dead-end in the network, the chemistry of more complex species will surely be boosted anyway if energetic particles from flares are included in the model.

Another limitation concerns the calculation of the viscous parameter α inside the disc. Although α was initially treated as a parameter to study accretion processes in discs, this parameter can be physically deduced from the stress tensor of the Navier-Stokes' equations. The physics underlying this parameter is thus inherently of MHD origin. In this paper, we computed the value of this parameter based on the properties of a static disc. Hence, the derived α value is inevitably imperfect. The intent of this computation is more a proof of concept, illustrating that the size of a disc dead zone is reduced when considering ionisation by flares. From the α calculation, we have also estimated the mass accretion rate, taking ionisation due to flares into account. This estimation too, should be viewed as a proof of concepts, given that it relies on the thin disc model of <cit.>, applicable far from the disc inner edge, a criterion not necessarily met in our study.



§ CONCLUSIONS


The ionisation of protoplanetary discs is crucial for the triggering of MHD instabilities, responsible for the accretion of matter onto the central star and is crucial in the disc chemistry, key to understand the formation of planets and the synthesis of the building blocks of life. 

In the first paper of this series <cit.>, we examined the effects of an additional ionisation source based on microphysical and observational considerations. A stationary state of the particle emission was assumed leading to an overestimation of the ionisation rate. In this paper, we address this issue by incorporating time-dependent components into our model.

We started by examining the average ionisation rate generated by a flare to compare it to the ionisation rate from the stellar X-rays. We conclude that both sources produce similar average ionisation rates at low column densities (N < 10^22cm^-2). However, particles generated by the flares produce an ionisation rate between 1 and 3 orders of magnitude higher than the ionisation rate from X-rays at high column densities (N > 10^22cm^-2). 

We used a Monte Carlo analysis to determine the spatial distribution of a temporally averaged ionisation rate. To facilitate the reuse of our numerical results, we proposed a parametric expression for the spatial distribution of the ionisation rate produced by flares (Eq. <ref>). 

Based on the thermo-chemical disc structures computed with , accounting for energetic particles produced by flares, we studied the impact of this additional ionisation source on the chemistry, viscosity, accretion rate, and volumetric heating rate of protoplanetary discs.

Impact on Chemistry: The ionisation fraction increases by at least an order of magnitude up to N∼ 10^25cm^-2 in disc region where particles penetrate. The distribution of cations is also significantly altered, especially the density of H^+ increases in the photodissociation region. Additionally, taking into account particles from flares leads to the formation of a layer where HCNH^+ is the most abundant cation at Z/R≈0.12. The recombination of this cation leads to an increase of a factor of 2 to 3 in the column density of hydrogen cyanide, HCN. This molecule is crucial in prebiotic chemistry to understand the synthetisation of amino-acids, the building blocks of life.

Impact on Viscosity and Accretion: 
The extent of the MRI-stable region or "dead zone", is reduced. We estimated that the mass accretion rate increases by almost one order of magnitude in disc regions where the energetic particles produced by flares penetrate.

Impact on the Volumetric Heating Rate: 
We determined the volumetric heating rate resulting from the ionisation of disc matter by energetic particles produced by flares. We showed that this volumetric heating rate is strong and would lead to collimated outflows. To facilitate the use of our numerical results, we proposed a parametric expression for the volumetric heating rate at the disc surface (Eq. <ref>). 


Given that particles from flares affect both the chemical and the thermal structure of the disc in regions that JWST/MIRI is capable of tracing, it is reasonable to expect that these flares would impact JWST spectra. We plan to refine our flare model using the synthetic spectra production feature of . However, multiple aspects related to synthetic JWST spectra production are currently under development within . So, we postpone this study to a forthcoming publication. 


We have shown that EPs produced by flares have strong impact on the inner disc of T Tauri discs. In a forthcoming work, we will present a particle acceleration model due to magnetic reconnection that could occur further out, at the surface of the discs. This work will extend the study presented in this paper to study the impact of energetic particles produced by turbulence-induced reconnection events. 

Additionally, until now, our focus has been on the influence of energetic particles on the ionisation rates, which subsequently affect the disc dynamics and chemistry. However, energetic particles can also collide with stable isotopes, creating radioactive isotopes, such as p + ^25Al→ ^26Al. The increased abundance of these radioactive species acts as a heat source within the disc. Previous work (e.g ) has shown that planet formation is significantly influenced by this heating source. Although the formation of planets was not the research topic of this paper, our particle acceleration model could be applied to this broader subject.




§ ACKNOWLEDGEMENTS

We acknowledge financial support from "Programme National de Physique Stellaire" (PNPS) and from "Programme National des Hautes Energies" (PNHE) of CNRS/INSU, France. This work has been conducted under the INTERCOS project of CNRS/IN2P3. CS thanks LUPM for hosting him  (délégations CNRS). Ch. Rab is grateful for support from the Max Planck Society and acknowledges funding by the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) - 325594231. The authors thank A. Vidotto and P. Caselli, Z. Meliani, S. Masson, T.P. Downes, J. Ferreira, G. Lesur, J. Morin, V. Guillet, D. Rodgers-Lee, for fruitful discussions. 





§ DATA AVAILABILITY

The data underlying this article will be shared on reasonable request to the corresponding author.







mnras
 



















§ NORMALISING THE PARTICLE INJECTION FLUX - REVISION OF B23 


We present here a revision of the calculation of B23 of the non-thermal particle flux produced by the flare. 

The energy distribution per unit volume, denoted as F(E), or the flux, denoted as j(E), is expressed in terms of the kinetic energy, E, for both electrons and protons, where 

    j(E)=F(E) v(E)/4π,

and v is the particle speed. The flux and distribution consist of a thermal and a non-thermal component, denoted by the subscript "th" and "nt", respectively. It is expressed in units of particles per unit energy, time, area, and solid angle.

The thermal component is defined by its temperature T and its normalisation n_ th, where n_ th is the thermal particle density in cm^-3. For energetic protons and electrons we assume equal temperatures. This gives us the following equation for their energy distribution,


    F_th(E) =  n_ th2/√(π)√(E/kT)1/k Texp(-E/kT)  .


Here, the peak value is at E = 3/2 k T ≡ E_ th.

As it is usually assumed in reconnection models for solar flares  <cit.> and supported by observation of solar flares <cit.> energy equipartition between non-thermal protons and electrons below 1 GeV is a good proxy of non-thermal content energetics. The energy density of non-thermal particles can be represented as,


    U_e,nt= U_p,nt=∫_E_ c^∞ E F_p(E)  E=∫_E_ c^∞ E F_e(E)  E
     ,


If we assume that electrons and protons have the same injection energy E_ c, they also share the same non-thermal energy distribution, 

    F_e(E)=F_p(E)=F_ nt(E).
 

As a result of the magnetic reconnection process, we introduce the non-thermal component of the energy distribution function F_ nt. The non-thermal energy distribution is determined by its normalisation N_ nt. In the following, we show a way to derive this quantity.

We define a maximum energy for the non-thermal distribution, E_ U > E_ c. The distribution is a power-law between these two energies. The index of the power law is called δ. This model is consistent with the non-thermal electron distribution observed in solar flares. 

The injection energy, E_ c, is proportional to the thermal energy, E_c=θ E_ th. Most T Tauri flare models exhibit a single power law. We assume an exponential cut-off beyond E_ U,


    F_ nt(E)=N_ nt(E/E_ c)^-δexp(-E/E_ U),


The normalisation factor N_ nt can be expressed in terms of the non-thermal particle density,

    n_ nt=∫_E_ c^∞(E/Ec)^-δexp(-E/E_ U) dE ≃ N_ ntE_c/δ-1.


The total electron density as function of the flare temperature T derived from X-ray observation (B23), 

    n_e(T)=20.5 × 10^10(T/1   MK)^-0.69  cm^-3


is the sum of the thermal and non-thermal components, n= n_ th+ n_ nt. By setting the non-thermal and the thermal energy distribution to be equal at the injection energy (E=E_c=θ E_ th), we can derive n_ nt,


    n_ nt=3θ/δ-1√(3θ/2 π)exp(-3/2θ) n_ th.

The non-thermal particle density depends on the local density of the thermal plasma. We use the total plasma density n from X-ray observations given by Eq. (<ref>) to deduce the thermal and non-thermal densities.
Given that n=n_ nt+n_ th, we get,


    n=n_ th(1+3θ/δ-1√(3θ/2 π)exp(-3/2θ)) .


From Eqs. (<ref>) and (<ref>), the thermal and non-thermal particle densities can be expressed as function of the observed temperature,

    n_ th= 20.5 × 10^10(T/1   MK)^-0.69/1+3θ/δ-1√(3θ/2 π)exp(-3θ/2)   cm^-3,

Assuming the fiducial values θ=3 and δ=3,

    n_ th=19 × 10^10(T/1   MK)^-0.69  cm^-3.


The non-thermal particle density n_ nt is then recovered from Eqs. (<ref>) and (<ref>),

    n_ nt=20.5 × 10^10(T/1   MK)^-0.69/1+δ-1/3√(2π/3 θ)exp(3/2θ)   cm^-3

Assuming the fiducial values θ=3 and δ=3,

    n_ nt= 10^10(T/1   MK)^-0.69  cm^-3.

  
We highlight here the difference between Eq. (<ref>) and Eq. 15 of B23. In that article, Eq. (15) is non-homogeneous. There is obviously a confusion between the normalisation factor of the non-thermal flux and the non-thermal particle density. We corrected the issue in this section. 
The wrong estimation of the normalisation factor N_ B23 lead to an overestimation of the flux in B23. We compare the normalisation N_ B23 with the corrected normalisation factor N_ nt,

    N_ B23/N_ nt= E_c/k T (δ-1)= 3θ/2(δ-1).


Taking the fiducial values of B23, δ=3 and θ=3, N_ B23≈ 2 N_ nt, the fluxes and ionisation rates are reduced by a factor of two compared to B23.





§ COLLISIONAL IONISATION 

In the chemical network used in <cit.>, collisional ionisation, which can dominate the ionisation fraction in the innermost high density and high temperature (T > 1000 K) <cit.> region of the disc, is not included. To model this process, which is usually described by the Saha equation, we add the following reaction to our network. 


    H_2 + Na→H_2 + Na^+ + e^-.

Usually, gas phase potassium K is used because of its low ionisation potential <cit.>. However, K is not included in our chemical network, and including it would require adding a couple of other chemical reactions (e.g. recombination). We therefore use the element Na as it is already in our network, but assume the ionisation potential of K to calculate the rate coefficient using Eq. 6 of <cit.>. 

To verify our approximation, we compare our models to various analytical expressions in Fig. <ref>. By simply including accretion heating, which increases the temperature in the midplane, the ionisation fraction increases by two orders of magnitude at R≈ 0.1 au. The inclusion of collisional ionisation increase by two orders of magnitude and expands the high ionisation region radially. This last model is in good agreement with the ionisation fraction solely considering the Saha equation (using Na as the main element with the ionisation potential of K). Close to the inner radius, our approximation under-predicts the ionisation rate compared to the Saha equation. This might be because our chemical network also considers other elements and molecules, which are neglected in the Saha equation. Nevertheless, with the addition of the accretion heating and collisional ionisation, the disc regions with T=1000 K now have a high enough ionisation rate to become MRI active, which was not the case in the model of <cit.> (see Sec. <ref>). 







