# Impacts of energetic particles from T Tauri flares on inner protoplanetary discs

V. Brunn,<sup>1★</sup> Ch. Rab<sup>1,2,3</sup> A. Marcowith<sup>1</sup> C. Sauty,<sup>1,4</sup> M. Padovani<sup>1,5</sup> and C. Meskini<sup>1</sup>

Accepted 2024 April 17. Received 2024 April 17; in original form 2023 December 28

## ABSTRACT

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 ionization 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 ionization rate expected in a radius range 0.08–0.6 au from the central star. We confirm that energetic particles produced by the flares dominate the ionization of the disc up to column densities of $10^{25}$ cm$^{-2}$. We further study the main consequences of this additional source of ionization on the viscosity, the accretion rate, the volumetric heating rate, and the chemical complexity of inner protoplanetary discs.

**Key words:** acceleration of particles – accretion, accretion discs – magnetic reconnection – stars: flare.

## 1 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 (Venuti et al. 2014; Hartmann, Herczeg & Calvet 2016; Manara et al. 2016). All mechanisms that describe accretion–ejection processes involve the interaction between magnetic fields and ionized disc matter (e.g. Shakura & Sunyaev 1973; Blandford & Payne 1982).

The magnetic field geometry of these stars is much more complex than those of non-accreting cool stars (Donati et al. 2007; Gregory et al. 2008). 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 (Zanni & Ferreira 2013). This process possibly leads to very strong flaring events observed in X-rays (Feigelson & Montmerle 1999; Getman et al. 2008; Getman & Feigelson 2021). 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 (Shibata & Magara 2011; Benz 2017; Oka et al. 2018). T Tauri flares appear to be thousands of times more energetic than solar

flares (Getman & Feigelson 2021). They are expected to produce particle fluxes that are five orders of magnitude higher than those produced by the Sun (Feigelson, Garmire & Pravdo 2002; Rab et al. 2017). The high-energy particle flux produced by T Tauri stars strongly impacts the surrounding environment (Rab et al. 2017; Rodgers-Lee et al. 2017; Offner, Gaches & Holdship 2019; Waterfall et al. 2020; Brunn et al. 2023). In particular, these particles can ionize the gas in the surrounding disc (Brunn et al. 2023). Disc ionization 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 (Cleeves, Adams & Bergin 2013; Offner et al. 2019; Waggoner & Cleeves 2022) and jets (Jacquemin-Ide, Ferreira & Lesur 2019; Sauty 2019; Ray & Ferreira 2021). Magnetic fields can also trigger instabilities in the gas that can drive accretion of disc material on to the star (Hartmann et al. 2016; Lesur 2021a). Besides, ionization drives chemistry in the disc, leading to the formation of molecules and complex organic compounds like amino acids (Semenov, Wiebe & Henning 2004), which are the building blocks of pre-biotic chemistry. Ionization also plays a role in the heating and cooling of discs (Glassgold, Galli & Padovani 2012), a process known to be essential in the launching of jets (Ray & Ferreira 2021). Understanding the mechanisms that drive disc ionization appears then necessary for understanding the formation and evolution of stars and the dynamics of accretion–ejection processes.

In the first article (Brunn et al. 2023, hereafter B23), we studied in stationary framework the ionization rate produced by magnetic

<footnote n="1"><sup>1</sup>Laboratoire Univers et Particules de Montpellier, Université de Montpellier/CNRS, place E. Bataillon, cc072, F-34095 Montpellier, France</footnote>

<footnote n="2"><sup>2</sup>Universitäts-Sternwarte, Fakultät für Physik, Ludwig-Maximilians-Universität München, Scheinerstr 1, D-81679 München, Germany</footnote>

<footnote n="3"><sup>3</sup>Max-Planck-Institut für extraterrestrische Physik, Center for Astrochemical Studies, Giessenbachstrasse 1, D-85748 Garching, Germany</footnote>

<footnote n="4"><sup>4</sup>Laboratoire Univers et Théories, Observatoire de Paris, Université PSL, Université Paris Cité, CNRS, F-92190 Meudon, France</footnote>

<footnote n="5"><sup>5</sup>INAF-Osservatorio Astrofisico di Arcetri, Largo E. Fermi 5, I-50125 Firenze, Italy</footnote>

<footnote>*E-mail: valentin.brunn@umontpellier.fr</footnote>
reconnection events occurring above the inner disc. We call inner disc the region at a radial distance between $R = 0.08$ and $0.6$ au from the central star. We estimated ionization rates produced by protons, electrons, and secondary electrons accelerated from these flares. Our model indicates that in the close environment of the flare, ionization by non-thermal particles is the dominant ionization source in the inner disc of T Tauri stars. It produces ionization 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 ionization rates produced by magnetic reconnection are dominant. The results suggest that energetic particles (EPs) accelerated above the disc during magnetic reconnection events enhances the ionization rates in the inner disc region of T Tauri stars. They are a new source of ionization important to be included in magnetohydrodynamics (MHD) or astrochemical codes.

Studying ionization in protoplanetary discs with a non-stationary approach allows to calculate the mean ionization fraction in the disc. Ionization 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 constraints. Another important aspect is to couple the ionization 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 ionization processes at play in these systems, we shall have a better understanding of the disc chemical evolution leading to the formation of planetesimals and the synthesis 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 ionization 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 Getman & Feigelson (2021) and Getman, Feigelson & Garmire (2021).

The layout of the paper is as follows. In Section 2 we summarize the stationary model as presented in B23, then present the physical model of flares, the chemical model, and the model to derive the viscous $\alpha$ parameter in the disc from the ionization fraction. In Section 3, 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 ionization over disc chemistry, the development of the magnetorotational instability (MRI), and the heat deposition in the disc and at the base of the jet. In Section 4, we propose a general discussion about theoretical perspectives of our work and its limitations before concluding in Section 5.

# 2 METHOD

## 2.1 Summary of the stationary model in B23

Elaborating on the work by Orlando et al. (2011) and Colombo et al. (2019), B23 proposed an uncovered source of ionization in the inner discs of T Tauri stars, the ionization by flares triggered by magnetic reconnection. The idea was also inspired by the results of the simulation of magnetospheric ejection by Zanni & Ferreira (2013) 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 $\lesssim 1$ GeV. Given the geometry of the field lines above the disc, we expect that this geometry favour 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 $\lesssim 1$ GeV, particles mainly lose their energy through ionization losses. In other words, they are very efficient at ionizing the disc. This supports the idea that our model could be a significant source of ionization in the inner disc. B23 estimate the ionization 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 *Chandra* Orion Ultradeep Project (COUP) point source catalogue (Favata et al. 2005; Getman et al. 2008). Using these data, B23 detail how to estimate the flare size and the electron density of the flare based on 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 with 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}$ to $10^{25}$ cm$^{-2}$ in order to treat the propagation of particles in the CSDA regime (Padovani et al. 2018). 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 provide 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 PRODIMO (see Section 2.3.1 for a further description). B23 focus 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 per cent 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 Dalgarno, Yan & Liu (1999) at low energy ($< 1$ keV), which are based on theoretical calculations. At high energy ($> 1$ keV) B23 use
the National Institute of Standards and Technology data base,<fnref n="1" /> see Padovani et al. (2022) 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 ionization and at high energy by radiative losses. The data of the Stopping and Range of Ion in Matter (Ziegler, Ziegler & Biersack 2010) are used to construct the proton loss functions. At high energies, above the threshold $E^\pi = 280$ MeV, we add energy losses due to pion production, as reported by Padovani et al. (2018). At very low energy ($< 1–10$ eV) Coulomb losses are dominant, we use the expression from Schlickeiser (2013).

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 ionization of the medium by primary particles (Ivlev et al. 2015).

Finally, B23 show that, in a stationary configuration, EPs are a powerful source of local ionization, with ionization 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 1 MK flare at 0.1 au propagating along a vertical magnetic field line, B23 find an ionization rate $\zeta = 10^{-9} \text{ s}^{-1}$ at column densities of $10^{25} \text{ cm}^{-2}$, while the ionization rate at this depth due to GCR or X-rays produced in a 1 MK stellar flare is of the order of $10^{-17} \text{ s}^{-1}$, while it is of the order of $10^{-18} \text{ 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 B23 suggest that their assumptions might result in an overestimation of the ionization rate, EPs generated by reconnection events are identified as a potentially dominant source of ionization 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 ionization rates produced are dominant over other ionization 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 \propto E^{-\delta}$ for $\delta < 6$, (ii) flares with temperatures above 1 MK, and (iii) particles propagating along the field line with a ratio of the toroidal component to the poloidal component $b_g = B_\phi/B_{pol} < 1$.

In this article, we have updated the method for calculating the non-thermal particle density based on temperature of B23 (see Appendix A). This revised approach is applied to determine the injected flux of non-thermal particles. We note that the ionization 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 ionization rate in the disc as it determines both the normalization and the injection energy of the non-thermal component of the particle flux. Fig. 1 shows the ionization rates generated by flares with temperatures ranging from $10^5$ to $10^9$ K in the fiducial configuration of B23. This configuration assumes a power law for the energy

<figure><img src="image_1.png" /><figcaption>**Figure 1.** Non-thermal particle ionization rate as a function of column density for a range of flare temperature between $10^5$ and $10^9$ K. The results are obtained assuming the fiducial parameters as described in B23. The flare is assumed to occur at a radius of 0.1 au from the star.</figcaption></figure>

distribution function $F(E) \propto E^{-\delta}$, and $\delta = 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 the simulations of Zanni & Ferreira (2013) and Mattia & Fendt (2020a, b, 2022), 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 Jacquemin-Ide, Lesur & Ferreira (2021) 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 (Ferreira 1997). However, self-similar models including MRI (Jannaud, Zanni & Ferreira 2023) 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.

## 2.2 Non-stationary flare model

As stated above, the model in B23 is stationary, and it estimates the ionization rate below the flare, at its maximum luminosity. Consequently, the resulting ionization rate is not representative of the ionization 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 ionization rate due to flares using the X-ray luminosity as an input. Therefore, it is necessary 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 Getman et al. (2008) 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_{\text{obs}}$ linking the peak X-ray luminosity in the *Chandra* detection energy band to the temperature deduced in this band $T_{\text{obs}}$. Secondly, we deduce the relationship

<footnote n="1">¹https://physics.nist.gov/PhysRefData/Star/Text/ESTAR.html</footnote>
$\log T_f - \log T_{\text{obs}}$, linking the bolometric and X-ray flare temperature. Getman et al. (2008) indicate that because the peak X-ray emission of flares with temperature $T_{\text{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_{\text{obs}} < 200$ MK.

In Fig. 2, we show the scattering plots of *Chandra* data points along with their corresponding linear regression. In Fig. 2(a), we plot the relation between the peak X-ray luminosity and the observed temperature of the flare $L_{X,\text{pk}}(\text{erg s}^{-1}) = 9.63 \times 10^{29} (T_{\text{obs}}/\text{MK})^{0.75}$. In Fig. 2(b), we plot the relation between the observed temperature in the *Chandra* band and the bolometric temperature of the flare $T_f = 3.35 \times T_{\text{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 \times \left( \frac{L_{X,\text{pk}}}{10^{30} \text{ erg s}^{-1}} \right)^{0.73} \text{MK.} \quad (1)$$

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 ($\approx 0.3$ yr).

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 Malzac (1999, chapter 6).

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 (Aschwanden & McTiernan 2010), enables 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* (Getman & Feigelson 2021). In addition, the data deduced from the X-ray flares of these young stellar objects provide the distribution of the flares duration.

### 2.2.1 Waiting times distribution

Observations of waiting times between solar flares reveal power-law distributions. Aschwanden & McTiernan (2010) and Aschwanden, Johnson & Nurhan (2021) propose a type II Pareto distribution based on data collected over 60 yr,

$$n_{\text{wt}}(\tau) = \lambda_0(1 + \lambda_0 \tau)^{-m}, \quad (2)$$

where $\tau$ is the time between two flares, and $1/\lambda_0$ is the breaking point of the power law. We use the value $1/\lambda_0 = 0.80 \pm 0.14$ h, 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,

$$\bar{\tau} = \int_{0}^{\infty} \tau n_{\text{wt}}(\tau) d\tau \approx 5 \text{ ks.} \quad (3)$$

### 2.2.2 Luminosity distribution

The exhaustive study by Getman & Feigelson (2021) 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} \text{ ergs}^{-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_{\text{flare}} = 10^{24} \text{ erg}$) up to superflare energies of solar-type stars ($E_{\text{flare}} = 10^{35} \text{ erg}$) (Notsu et al. 2019; Okamoto et al. 2021). This justifies our extrapolation of the power-law distribution observed by Getman & Feigelson (2021) to luminosity lower than $10^{32} \text{ erg s}^{-1}$.

The normalized X-ray luminosity distribution can be written as

$$n_X(L) = \frac{a-1}{L_c} \left(\frac{L}{L_c}\right)^{-a}, \quad (4)$$

where $a = 2.11$ is the power-law index measured by Getman & Feigelson (2021) and $L_c$ is the low-luminosity cut-off of the distribution. $L_c$ is adjusted in order to reproduce the occurrence rate $N_{\text{MF}}$ of mega flares, flares with peak X-ray luminosity $L_{X,\text{pk}} > L_0 = 10^{32.5}$, observed by Getman & Feigelson (2021), which is $N_{\text{MF}} = 1.7$ mega flares per year per star.

In the previous paragraph we have calculated the mean waiting time between two flares $\bar{\tau} = 5$ ks. We can estimate the total number of flares per year, $N_{\text{tot}} = 6300$ flares yr$^{-1}$ star$^{-1}$.

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

$$N_{\text{tot}} = \frac{N_{\text{MF}}}{P(L_{X,\text{pk}} > L_0)}, \quad (5)$$

where $P(L_{X,\text{pk}} > L_0)$ is the probability of occurrence of a megaflare, i.e.

$$P(L_{X,\text{pk}} > L_0) = \int_{L_0}^{\infty} n_X(L) dL = \left(\frac{L_0}{L_c}\right)^{1-a}, \quad a > 1. \quad (6)$$

From the two above equations, we get

$$L_c = L_0 \left( \frac{N_{\text{MF}}}{N_{\text{tot}}} \right)^{\frac{1}{a-1}}. \quad (7)$$

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

$$L_c \approx 1.9 \times 10^{29} \text{ erg s}^{-1}. \quad (8)$$

### 2.2.3 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 (Benz 2002). It is characterized by an exponential fast growth and gradual decay. There are three parameters in this model. The first parameter is the
<figure><img src="image_2.png" /><figcaption>**Figure 2.** Panel (a) shows the flares peak X-ray luminosity $L_{X,pk}$ dependence as a function of the flare observed temperature $T_{obs}$ in the Chandra band. Panel (b) shows the flares bolometric temperature $T_f$ dependence as a function of the flare observed temperature $T_{obs}$ in the Chandra band. The fits are plotted in solid orange line, while the data from Getman et al. (2008) are plotted as blue dots.</figcaption></figure>

characteristic rise time ($\tau_r$), the second is the characteristic decay time ($\tau_d$), and the last is 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) = \begin{cases} L_{X,pk} e^{\frac{(t-t_0)}{\tau_f}} & \text{if } t < t_0 \\ L_{X,pk} e^{-\frac{(t-t_0)}{\tau_d}} & \text{otherwise.} \end{cases} \qquad (9)$$

The peak luminosity is randomly chosen according to the distribution in the previous section. Fig. 6 in Getman & Feigelson (2021) 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(\tau_r) = \lambda \exp(-\lambda \tau_r), \qquad (10)$$

where $1/\lambda = 26$ ks is the average rise time of the flares from the observed sample by Getman & Feigelson (2021).

Fig. 6 in Getman & Feigelson (2021) shows a correlation between flare growth and decay times (both in ks units):

$$\tau_d = 5.28 \tau_r^{0.60} \qquad (11)$$

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 (Benz 2017).

The time-dependent ionization rate, represented by $\zeta(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 ionization rate matches the phase of increasing soft X-ray luminosity.

$$\zeta(t, L_{X,pk}) = \begin{cases} \zeta_{pk}(L_{X,pk}) e^{\frac{(t-t_0)}{\tau_r}} & \text{if } t < t_0 \\ 0 & \text{otherwise,} \end{cases} \qquad (12)$$

where $\zeta_{pk}$ is the ionization 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 \times 10^{30}$ erg s$^{-1}$, which is close to the typical value of $2 \times 10^{30}$ erg s$^{-1}$ of the X-ray continuum observed on PMS (Flaccomio et al. 2003; Wolk et al. 2005).

### 2.2.4 Flare geometry model

Flares, characterized by sudden energy discharges within a tenuous plasma confined to a 'magnetic bottle', have been extensively explored in various works (Reale 2007; Benz 2017). The energy dissipation in these scenarios predominantly occurs via optically thin radiation and efficient heat conduction directed towards the chromosphere of the star (Reale 2014). 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 (Serio et al. 1991; Reale 2007). McCleary & Wolk (2011) proposed an improvement to the original Serio et al. (1991) relationship, expressed as follows:

$$L_f = \frac{\tau_d T^{1/2}}{\gamma F(\epsilon)}, \qquad (13)$$

where $\gamma \simeq 3.7 \times 10^{-4}$ is generally valid for solar-type stars and the empirical function $F(\epsilon) = \frac{0.63}{\epsilon^{-0.32}} + 1.41 > 1$. The introduction of the function $F(\epsilon)$ includes the effect of a progressive heating in the rise phase of the flare, whereas Serio et al. (1991) assumed an impulsive heating with $F(\epsilon) = 1$. The parameter $\epsilon$ depends on the star. It varies for a given star from one flare to the other. McCleary & Wolk (2011) derived the mean value of this parameter for a sample of flares in Classes I, II, and III young stellar objects. In this study
<figure><img src="image_3.png" /><figcaption>**Figure 3.** Sketch of a flaring flux tube reproduced from Waterfall et al. (2020). Magnetic reconnection takes place at the marked X-point, in the interaction region between the star and disc magnetic field lines. Non-thermal particles move in a circular motion from this point, travelling along the field lines towards the surface of the star emitting X-rays. According to this model particles are entering the disc at $R \approx L_f$.</figcaption></figure>

<figure><img src="image_4.png" /><figcaption>**Figure 4.** Possible geometry of a T Tauri flare reproduced from Hamaguchi et al. (2012). Differential rotation between the star and the disc shears the stellar bipolar magnetic fields, and the magnetic fields twist and continuously reconnect. Matter, accelerated by the magnetic reconnection, collides with the stellar surface, thermalizes, and emits hard X-rays. According to this model particles are entering the disc at $R \approx 2L_f$.</figcaption></figure>

we are interested in Class II young stellar object and for this class $\langle\epsilon\rangle = 1.43$ and $\langle F(\epsilon) \rangle \approx 1.98$, we hence assume $F(\epsilon) = 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 Getman et al. (2021) where the two footpoints of the flare are anchored on to the central star and the flare loop-top is in the disc, see Fig. 3. Waterfall et al. (2020) 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 (Orlando et al. 2011; Zanni & Ferreira 2013; Colombo et al. 2019) and by the observations of Favata et al. (2005) with Chandra and Hamaguchi et al. (2012) 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. 4.

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 the disc. Assuming the first model (SSF), with the two footpoints anchored into the star, the particles are penetrating at $R \approx 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 \approx 2L_f$. The second model injects particles farther in the disc. We show that this has a major impact on the ionization spatial distribution in the inner disc, thus potentially to the extent of the MRI active region.

From equation (1), we can assign to each flare a half loop length given its peak X-ray luminosity $L_{X, pk}$ and its decay time $\tau_d$. These values are picked from the distributions obtained in equations (6) and (11) and allow to derive a radial position $R$.

Getman & Feigelson (2021) also made an estimation of $\chi$, 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 $\chi$ as

$$A_f = \pi \chi^2 L_f^2. \qquad (14)$$

In solar and stellar flares as well as in the PMS flares observed by Getman & Feigelson (2021), the typical value of $\chi$ 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 a uniform distribution of angles.

## 2.3 Chemical model

### 2.3.1 PRODIMO

To evaluate the column density explored by EPs, we use a model generated by the radiation thermochemical code PRODIMO (PROto-planetary DIsc MOdel,<fnref n="2" /> Woitke, Kamp & Thi 2009; Kamp et al. 2010; Thi, Woitke & Kamp 2011; Woitke et al. 2016). This model performs a wavelength-dependent radiative transfer calculation, including X-rays (Rab et al. 2018), 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 (Kamp et al. 2017; Rab et al. 2017). Here, we use a fiducial T Tauri disc model described in Rab et al. (2017). It was also used in Brunn et al. (2023) (see their table 1 for the main physical parameters). Here, we want to model in a more realistic way the thermal and ionization structure in the innermost region of the disc ($R < 1$ au). Thus, we additionally include accretion heating following D’Alessio et al. 1998, see also Oberg et al. 2022 appendix A). We use a total mass accretion rate on to the star of $M_{accr} = 10^{-8} M_{\odot} yr^{-1}$. We added a simple treatment of the collisional ionization to the chemical network following Desch & Turner (2015). For further details, see Appendix B.

### 2.3.2 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. None the less, 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 equation (7) by Fromang, Terquem & Balbus (2002), 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:

$$\frac{dn_e}{dt} = \zeta n_n - \beta n_e^2. \qquad (15)$$

In the no-flare case, assuming steady state, the initial electron density and neutral species density, denoted, respectively, by $n_e, 0$

<footnote n="2">²https://prodimo.iwf.oew.ac.at revision: 66efbd75 2023/06/27</footnote>
and $n_{n,0}$, are computed by the PRODIMO 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}. \quad (16)$$

In equation (15), the parameter $\beta = \beta(T)$ is the recombination rate from equation (9) of Fromang et al. (2002). We note that $\zeta$, the total ionization rate, is the sum of two components, $\zeta = \zeta_f + \zeta_{\text{Prod}}$, with the stationary ionization rate, $\zeta_{\text{Prod}}$, computed by PRODIMO in the scenario without flare, ensuring a stationary $n_{e,0}$. On the other hand, $\zeta_f$, the ionization rate due to flares, is derived from equation (12).

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

## 2.4 Viscosity model

As already stated, in order to develop, MRI requires a sufficiently high degree of ionization for the gas to interact with the magnetic field. The instability can be characterized by the viscous parameter $\alpha$.

Wardle & Salmeron (2012) pointed out that Hall diffusion could, under certain circumstances, counteract the resistive damping effects of MRI. However, numerical simulations indicate a more complex scenario (Lesur, Kunz & Fromang 2014). Given these uncertainties, Thi et al. (2019) chose to exclude the effects of Hall diffusivity in estimating the effective viscosity parameter. To compute the viscous $\alpha$ parameter, we rely on the model proposed by Thi et al. (2019), which introduces an empirical non-ideal MHD MRI-driven equation for $\alpha$ suitable for physico-chemical protoplanetary disc codes such as PRODIMO. The adopted expression for $\alpha$ is given by

$$\alpha = \left(\frac{2}{\beta_{\text{mag}}}\right)^{1/2} \min(1, \Lambda_{\text{Ohm}}) \left[ \left(\frac{50}{\text{Am}^{1.2}}\right)^2 + \left(\frac{8}{\text{Am}^{0.8}} + 1\right)^2 \right]^{-1/2}, \quad (17)$$

if $\sqrt{\beta_{\text{mag}}} \Lambda_{\text{Ohm}} > 1$ and $\alpha \approx 0$ otherwise.

$Am$ represents the frequency at which neutral particles collides with ions normalized to the Keplerian orbital frequency,

$$Am = \frac{v_{\text{in}}}{\Omega_K} = \frac{\beta_{\text{in}} n_{\text{charge}}}{\Omega_K} = \frac{x_e \beta_{\text{in}} n_{\text{tot}}}{\Omega_K}, \quad (18)$$

where $n_{\text{charge}} = x_e n_{\text{tot}}$ is the total number density of charged species and $\beta_{\text{in}} = 2 \times 10^{-9} \text{ cm}^3 \text{s}^{-1}$ is the collisional rate coefficient for ions to distribute their momentum to neutrals. Thus,

$$Am \approx 1 \left(\frac{x_e}{10^{-8}}\right) \left(\frac{n_n}{10^{10} \text{ cm}^{-3}}\right) \left(\frac{R}{1 \text{ au}}\right)^{3/2}. \quad (19)$$

The Ohmic diffusivity is quantified by the dimensionless Elsasser number $\Lambda_{\text{Ohm}}$, defined as the ratio of the Lorentz force to the Coriolis force,

$$\Lambda_{\text{Ohm}} = \frac{B_z^2}{4\pi\rho\eta_O\Omega_K} = \frac{v_A^2}{\eta_O\Omega_K} = \left(\frac{4\pi\sigma_O}{\Omega_K}\right)\left(\frac{v_A}{c}\right)^2,$$

where $B_z$ is the vertical component of the magnetic field, $v_A$ is the Alfvén speed, and $\rho$ is the mass density of the plasma and the Ohmic resistivity is given by $\eta_O = \frac{c^2}{4\pi\sigma_O}$. Thi et al. (2019) derived an approximation for the Ohmic Elsasser number in disc regions relevant to our study, where $\sigma_O \approx \sigma_{e,0}$ and $x_e > 10^{-13}$:

$$\Lambda_{\text{Ohm}} \approx 1 \left(\frac{T}{100 \text{ K}}\right)^{1/2} \left(\frac{10^4}{\beta_{\text{mag}}}\right) \left(\frac{R}{1 \text{ au}}\right)^{3/2} \left(\frac{x_e}{10^{-9}}\right).$$

<figure><img src="image_5.png" /><figcaption>**Figure 5.** Plot of the electron density versus time in a range of column densities between $10^{19}$ and $10^{25} \text{ cm}^{-2}$. The flare temperature is 10 MK flare and the particles are penetrating into the disc at 0.1 au from the central star. The EP ionization rate follows the profile of equation (12) with the rising characteristic time $\tau_r = 26$ ks being the average value observed by Getman & Feigelson (2021).</figcaption></figure>

The magnetic term $\beta_{\text{mag}}(R, z)$ is the ratio of the thermal pressure to the magnetic pressure, given by (Thi et al. 2019)

$$\beta_{\text{mag}}(R, Z) = \beta_{\text{mid}} \frac{P_{\text{th}}(R, Z)}{P_{\text{th}}(R, 0)} = \beta_{\text{mid}} \frac{n(R, Z)T(R, Z)}{n(R, 0)T(R, 0)}, \quad (20)$$

where $\beta_{\text{mid}} = \beta_{\text{mag}}(R, 0)$ is the value of $\beta_{\text{mag}}$ on the mid-plane. $\beta_{\text{mid}}$ is assumed to be independent of the radius so that it can be determined from the thermal structure and chemical abundances of PRODIMO. The parameter $\beta_{\text{mid}}$ is typically in the range $10^4-10^6$. We take $\beta_{\text{mid}} = 10^4$ as reference value.

# 3 RESULTS

## 3.1 Single flare effects

### 3.1.1 Ionization rates

We plot in Fig. 5 the solution of equation (15) for the number density of electrons. The solution is obtained by taking the time-dependent ionization 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}$ to $10^{25} \text{ cm}^{-2}$. It illustrates that for column densities below $10^{24} \text{ 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 $\tau_{\text{rec}} = (\beta(T)n_e)^{-1} \sim 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} \text{ cm}^{-2}$, despite the low temperature ($\approx 500$K), the electron density is very low so the recombination times are longer, $\tau_{\text{rec}} \gtrsim 10^5 \text{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.
<figure><img src="image_6.png" /><figcaption>**Figure 6.** Ionization rate as a function of the column density. The solid blue line shows the effective ionization rate from a 10 MK flare averaged over a Keplerian time. In comparison, the solid orange line shows the X-ray ionization rate of a T Tauri star as computed by PRODiMO.</figcaption></figure>

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

In B23, we argued that the ionization rates obtained in the stationary state were overestimated and in order to have more realistic description of the ionization processes, temporal and spatial averaging should be performed. We propose here a simple estimation of the spatial and temporal average of the ionization 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 $\tau_r = 26$ ks corresponding to the average rising time observed by Getman & Feigelson (2021). We define an effective time-averaged ionization rate $\zeta_{\text{eff}}$ from the average ionization fraction during a period of time P,

$$ \zeta_{P,\text{eff}} = \frac{1}{P} \langle x_e \rangle_P = \frac{1}{P} \sqrt{\frac{\langle \zeta \rangle_P}{\beta n_n}}, \quad (21) $$

where $\beta$ is the recombination rate, $x_e$ is the ionization fraction, and $\langle \cdot \rangle_P$ is the time averaging over a period of time P. The temporal averaging of the ionization rate produced by a single flare over a period P is

$$ \langle \zeta_P \rangle = \frac{1}{P} \int_0^P \zeta(t) dt. \quad (22) $$

Given the temporal profile of equation (12), the average ionization rate produced by a flare occurring at $t = t_0 \ll P$ is

$$ \zeta_P = \zeta_{pk}(L_{x,pk}) \frac{\tau_r}{P}. \quad (23) $$

We compute $\zeta_{P,\text{eff}}$ over a Keplerian period at 0.1 au, so $P \approx 10$ Ms We then define the spatially averaged ionization rate as

$$ \zeta_{\text{eff}} = \frac{r_f^2}{(R + r_f)^2 - R^2} \zeta_{P,\text{eff}}, \quad (24) $$

where R is the radial position of flare and $r_f$ the radius of the flare. $\zeta_{\text{eff}}/\zeta_{P,\text{eff}}$ is the ratio between the area illuminated by particles and the disc ring of width $2r_f$ at radius R. Here, we take $R = 0.1$ au and $r_f/R = 0.1$. This effective ionization rate is illustrative for comparison with other sources of ionization in the inner disc. The effective ionization rate is plotted in Fig. 6 and compared to the X-ray ionization rate computed by PRODIMO. It was anticipated that calculating spatial and temporal averaged ionization rates would result in lower values

than the ionization rate estimated at a flare maximum in B23. This hypothesis is validated here. Considering that observations span over extended periods of time, the averaged ionization rates calculated in this study offer a better basis for comparison with observational data than the ionization rate obtained in B23. Furthermore, these averaged ionization rates also provide more relevant values to be used in numerical models of stationary chemistry in discs such as PRODIMO. As seen in Fig. 6, the effective ionization rate produced by flares at low column densities, $N \lesssim 10^{22} \text{cm}^{-2}$, is similar to the ionization rate produced by X-rays. However, at high column densities, $N \gtrsim 10^{22} \text{cm}^{-2}$, the effective ionization rate produced by flares becomes one to two orders of magnitude higher than that produced by X-rays. The increase in ionization rate by EPs compared to X-rays at high column densities is an expected result as has been demonstrated in Rab et al. (2017), Padovani et al. (2018), and Offner et al. (2019).

### 3.1.2 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_{\text{nt}}$, exerted by a distribution of electrons and protons ($i = e, p$) can be expressed as

$$ P_{\text{nt}} = \frac{1}{3} \int_{E_c}^{\infty} p(E) v(E) F_{i,\text{nt}}(E) dE, \quad (25) $$

where E is the kinetic energy, $E_c$ is the injection energy of the non-thermal particles, and $F_{i,\text{nt}}(E)$ is the energy distribution of these non-thermal particles i,

$$ F_{i,\text{nt}}(E) = N_{i,\text{nt}} \left(\frac{E}{E_c}\right)^{-\delta}. \quad (26) $$

We assume equipartition of energy between non-thermal protons and electrons, i.e $F_{e,\text{nt}} = F_{p,\text{nt}}$.

We consider both species with the same injection energy $E_c$. Thus, $P_{e,\text{nt}} = P_{p,\text{nt}} = P_{\text{nt}}$, with

$$ P_{\text{nt}} \approx \frac{2}{3} \frac{\delta - 1}{\delta - 2} E_c n_{\text{nt}}, \quad (27) $$

where $n_{\text{nt}}$ is the non-thermal particle density (equation A13) and $\delta \neq 2$. We assume that the injection energy is proportional to the thermal energy, $E_c = \theta E_{\text{th}}$, where $E_{\text{th}} = \frac{3}{2} k_B T$; here $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_{\text{nt}} = \frac{3\theta^2}{\delta - 2} \sqrt{\frac{3\theta}{2\pi}} \exp\left(-\frac{3\theta}{2}\right) P_{\text{th}}. \quad (28) $$

The thermal pressure can be computed as a function of the flare temperature from $P_{\text{th}} = n_{\text{th}}(T) k_B T$, where $n_{\text{th}}(T)$ is the plasma thermal density constrained from X-ray observation, equation (A11).

Assuming the fiducial values of B23, $\theta = 3$ and $\delta = 3$,

$$ P_{\text{th}} = 26 \left( \frac{T}{1 \text{MK}} \right)^{0.31} \text{erg cm}^{-3} \quad \text{and} \quad P_{\text{nt}} = 9 \left( \frac{T}{1 \text{MK}} \right)^{0.31} \text{erg cm}^{-3}. \quad (29) $$

The total pressure of non-thermal particles accounting for both electrons and protons is $P_{\text{tot,nt}} \approx 2P_{\text{nt}}$ and the total thermal pressure is $P_{\text{tot,th}} \approx 2P_{\text{th}}$.

In Fig. 7, 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–0 G (Orlando et al. 2011; Colombo
<figure><img src="image_7.png" /><figcaption>**Figure 7.** Different components of the pressure in the disc as a function of the flare temperature in the range 1–100 MK: pressure of thermal particle distribution produced by flares (solid blue line); pressure of the suprathelial component of the particle distribution (solid orange line); magnetic pressure resulting from a 27 G magnetic field, as used in the simulation of the thermal effects of flares by Colombo et al. (2019) (dashed green line); thermal pressure at the disc surface, with $T = 10^3$ K and $n = 10^{10}$ cm$^{-3}$ (dashed red line). The grey shaded area shows the range of magnetic pressure exerted by a magnetic field with intensities between 10 and 50 G.</figcaption></figure>

et al. 2019). 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 (Orlando et al. 2011; Colombo et al. 2019). 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. 7. Colombo et al. (2019) 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 equation (29), 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.

## 3.2 Ionization 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. Waggoner & Cleeves (2022) observed 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 time-scales 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 time-scales for these changes vary depending on the molecule in question and the prevailing local conditions. Ilgner & Nelson (2006) also argued that for flares, the chemistry of ions in the disc should be modelled using time-averaged ionization 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 ionization rate. We use a sample of flares that are randomly drawn using the method described in Section 2.2, over a time period $P$.

First, we need to deal with the temporal averaging of the ionization rate produced by a single flare over a period $P$, as in equation (22).

Given the temporal profile of equation (12), the average ionization rate produced by a flare occurring at $t = t_0$ is

$$ \zeta_P = \zeta_{pk}(L_{X,pk}) \frac{\tau_r}{P} \left(1 - e^{-\frac{P-t_0}{\tau_r}}\right) \quad (30) $$

We now calculate the spatial average of the ionization rate produced by the whole sample of flares. The spatial average is computed in disc rings with a radius $R$ and width $\Delta R$. The area of such a ring is $A_R = \pi(R^2 - (R - \Delta 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 equation (14). Flares produce a time-averaged ionization rate $\zeta_{P,i}(N)$ in a volume of section $A_i$.

The spatial average ionization rate $\zeta_{P,R}$ in a disc ring with a large radius $R$ is given by

$$ \zeta_{P,R}(N) = \frac{1}{A_R} \sum_i A_i \zeta_{P,i}(N). \quad (31) $$

We have calculated the average ionization 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 ionization rate plotted in Fig. 8(a) is obtained by taking $\Delta R \ll R$ and averaging over 1000 yr.

From the numerical results we propose a parametrization of the average ionization rate due to the EPs produced by flares of the form

$$ \zeta(R, N) = \zeta_0 \left(1 + \frac{N}{N_0}\right)^{-a} e^{-b(R-R_0)}. \quad (32) $$

The parameters for the SSF and SDF models are presented in Table 1. This expression provides a fit matching the numerical results with an average relative deviation of less than 10 per cent for column densities ranging from $10^{19}$ to $10^{25}$ cm$^{-2}$, for radii 0.1–0.3 au and 0.1–0.6 au for SSF and SDF, respectively.

## 3.3 Feedback on the disc

### 3.3.1 Effects on the disc chemistry

In this section, we explore how the ionization 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 ionization rate generated by both SSF and SDF represented in Fig. 8. PRODIMO then recalculates the entire disc chemistry, accounting for this new ionization source. Our model computes ionization rate for column densities $N < 10^{25}$ cm$^{-2}$. At higher column density the ionization rate is extrapolated using linear interpolation in log space. The revised disc structures are displayed in Fig. 9 where we specifically examine the ionization fraction. In Fig. 10, 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.
<figure><img src="image_8.png" /><figcaption>**Figure 8.** Time-averaged ionization rate in the range $10^{19}-10^{25}$ cm$^{-2}$ over a period of 1000 yr. Panel (a) shows the results obtained from the SSF model. Panel (b) shows the results from the SDF model.</figcaption></figure>

<table><caption>**Table 1.** Parameters of the ionization rate expression equation (32).</caption><thead><tr><th>Parameters</th><th>ζ<sub>0</sub> [s<sup>-1</sup>]</th><th>N<sub>0</sub> [cm<sup>-2</sup>]</th><th>a</th><th>b [au<sup>-1</sup>]</th><th>R<sub>0</sub> [au]</th></tr></thead><tbody><tr><td>SSF</td><td>1.5 × 10<sup>-4</sup></td><td>1.2 × 10<sup>22</sup></td><td>2.35</td><td>26.5</td><td>0.075</td></tr><tr><td>SDF</td><td>4.6 × 10<sup>-5</sup></td><td>7.9 × 10<sup>21</sup></td><td>2.28</td><td>11.3</td><td>0.15</td></tr></tbody></table>

In addition to the ionization fraction, in Fig. 9, the solid black line shows the delimitation of the MRI-inactive region, where $\alpha = 0$ (i.e. where $\sqrt{\beta_{mag} \Lambda_{Ohm}} < 1$). We assume $\beta_{mid} = 10^4$ (see equation 20). The method to determine the MRI-active region is explained in Section 2.4, and further discussed in Section 3.3.2.

The overall ionization fraction within the disc is determined by the equilibrium between ionization 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. 9 reveal that, in disc structures influenced by flaring activity, the ionization 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 ionization fraction extends farther radially, see Fig. 9. Specifically, the region in which the ionization 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 ionization fraction distribution at 0.3 and 0.6 au in Fig. 9(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. 10 shows the distribution of the most abundant ion species in the inner disc without flares. At the disc surface, the ionizing flux arising from the Far Ultraviolet (FUV) photons is predominant, and the gas layer essentially acts as a photodissociation region (PDR). In this region, atomic ions like $H^+$, $C^+$, $S^+$, and $Mg^+$ are the most abundant cations.

In the UV- and X-ray-shielded inner disc, collisions due to high temperature ionize atoms having low ionization 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 ionization 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 an HCNH$^+$-dominated layer that is absent in the model without flare.

An important aspect for the chemistry in the context of variable flares is the time-scales attached to chemical reactions. In case of the EP-induced ionization the picture is twofold. The ionization due to the particles does not affect only the upper UV-dominated layers, but also deeper regions of the disc (towards the disc mid-plane). These two regions can roughly be separated by the location of the $A_V = 1$ layer (see e.g. Fig. 11). As a consequence the region affected by the particles also spans a huge range in chemical time-scales. In the UV-dominated layers those time-scales a very short (about 1 d), whereas towards the mid-plane they span from $\approx 10^3$ to $\gtrsim 10^6$ yr (see Woitke et al. 2009, figs. 12, 13 and equation 117 for the definition of the chemical time-scale in our model).

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^- \rightarrow HCN + H$ (e.g. Long et al. 2021). This process is the most important in deeper layers (below the $A_V = 1$ limit) of the disc with longer chemical time-scales and is a consequence of the ion-molecular formation pathway as described in Walsh, Nomura & van Dishoeck (2015, see e.g. their fig. 5), which also affects other molecules such as CN. Woitke et al. (2024) pointed out that X-ray ionization 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 ionization rates deeper in the disc and hence those formation pathways become even more important.

In Fig. 13 we show the HCN column density for models with and without flares. This result is given by the PRODIMO steady-state chemistry model assuming averaged particle ionization rates from the flares. As the chemical time-scales 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
<figure><img src="image_9.png" /><figcaption>**Figure 9.** Expected distribution of the ionization fraction as a function of the disc radius, *R*, and the normalized height, *Z/R*, in the case of no flare (panel a), SSF (panel b), and SDF (panel c). The hatched areas show the MRI-inactive regions.</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>**Figure 10.** Expected distribution of the most abundant cation as a function of the disc radius, *R*, and the normalized height, *Z/R*, in the case of no flare (panel a), SSF (panel b), and SDF flare (panel c). The white areas are regions where other cations are dominant like $H_3O^+$, $NH_4^+$, or $SiOH^+$.</figcaption></figure>

<figure><img src="image_11.png" /><figcaption>**Figure 11.** From left to right: disc models with no flares, SSF, and SDF. Top row: total (X-ray + SP + cosmic rays) ionization rate per molecular hydrogen. Bottom panel: gas temperature. The red solid contour shows the visual extinction of unity, and the orange solid contour shows the total vertical column density of $10^{25} \text{ cm}^{-2}$.</figcaption></figure>
<figure><img src="image_12.png" /><figcaption>**Figure 12.** Accretion rate in the inner disc assuming a disc model impacted by the EPs produced by SDF (orange), SSF (blue), and in the absence of flare (green).</figcaption></figure>

<figure><img src="image_13.png" /><figcaption>**Figure 13.** Vertically integrated column densities for HCN in the inner one au, for models without flares, with SDF and SSF.</figcaption></figure>

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. First, 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 pre-biotic molecules. For instance, amino acids are commonly believed to originate from such processes. An aldehyde interacts with HCN and NH₃, resulting either in an aminonitrile, which is then hydrolysed to produce the corresponding amino acid or in a reaction between an aldehyde, HCN, NH₃, and CO₂ that leads to the production of amino acids through hydrolysis (Taillades et al. 1998). Our model, accounting for the ionization 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.

### 3.3.2 Effects on MRI and accretion rates

Fig. 14 displays the distribution of the effective viscosity parameter α for a β<sub>mid</sub> value of 10<sup>4</sup>. This parameter traces the turbulence

efficiency, see Section 2.4. 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 β<sub>mag</sub><sup>1/2</sup> Λ<sub>Ohm</sub> < 1, and the portion of the disc where β<sub>mag</sub> < 1 in the disc atmospheres leading to a drop of α < 0.01. The decrease in α, as one moves outwards the disc, is linked to the reduction in gas pressure and consequently, β<sub>mag</sub>.

In regions with higher ionization rates, the coupling between the gas and the magnetic field is stronger, leading to more efficient MRI-driven turbulence and angular momentum transport, see Section 2.4 and Turner, Sano & Dziourkevitch (2007).

We focus now on the effects of additional ionization due to flares. In the region undergoing an increased ionization rate, at R < R<sub>i</sub>, where R<sub>i</sub> = 0.3 au and R<sub>i</sub> = 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. 14 panels (b) and (c), we see that the size of the MRI-active region increases 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, ionization processes regulate the mass accretion rate on to 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 (Perez-Becker & Chiang 2011),

$$ \dot{M} = 3\pi \nu \Sigma = 2 \times 3\pi \mu m_p N_g \alpha c_s^2 \quad (33) $$

$$ = 2.5 \times 10^{-9} \left( \frac{\alpha}{0.01} \right) \left( \frac{N_g}{10^{24} \text{ cm}^{-2}} \right) \left( \frac{T}{10^3 \text{ K}} \right) M_{\odot} \text{yr}^{-1}, \quad (34) $$

where the effective viscosity is given by ν = αc<sub>s</sub>h and μ is the mean molecular weight. c<sub>s</sub> and h are the sound speed and the height scale of the disc, respectively, defined in Section 2.4. The factor of 2 accounts for accretion along the upper and lower disc surfaces; N<sub>g</sub> is the gas column density above the mid-plane and T is the disc temperature. This formula is an approximation in our case as it assumes very efficient angular momentum transport, namely l<sub>in</sub> ≪ l, where l = R<sup>2</sup>Ω 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 Orlando et al. (2011). In the region undergoing the increased ionization rate, i.e. at R < R<sub>i</sub>, where R<sub>i</sub> = 0.3 au and R<sub>i</sub> = 0.6 au for SSF and SDF, respectively, over 1000 yr, the averaged accretion rate increases by a factor of 7–10 compared to the model without flares, see Fig. 12.

### 3.3.3 Disc heating by EPs

Ionization processes in T Tauri discs contribute to the heating of the disc material, which affects the disc thermal structure (D'Alessio et al. 1998; Ercolano, Clarke & Drake 2009). Various ionization 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 (Igea & Glassgold 1999; Rab et al. 2017). In this section, we compute the volumetric heating rate due to the ionization 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.
<figure><img src="image_14.png" /><figcaption>**Figure 14.** Expected distribution of the effective viscosity parameter (equation 17) as a function of the disc radius, R, and the normalized height, Z/R, in the case of no flare (panel a), SSF (panel b), and SDF (panel c). The white area delimits the MRI-inactive region, inside which $\sqrt{\beta_{mag}}\Delta_{Ohm} < 1$.</figcaption></figure>

The volumetric heating rate due to EPs, $\Gamma$, is computed assuming that the ionization by each atomic hydrogen releases an amount of energy $Q_H = 4.3$ eV and ionization by each molecular hydrogen releases $Q_{H_2} = 18$ eV (Glassgold et al. 2012). $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 (Glassgold et al. 2012). So,

$$ \Gamma(R, N) = (n_H(R, N)Q_H + n_{H_2}(R, N)Q_{H_2})\zeta(R, N), \quad (35) $$

where $\zeta$ is the ionization rate per H nuclei. The volumetric heating rate in the disc as calculated in PRODIMO is shown in Fig. 15.

Previous works have indicated that heating the surface of the disc can lead to an increase in the mass of the outflow (Casse & Ferreira 2000; Bai et al. 2016). Casse & Ferreira (2000) 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 ionization 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 Alexander et al. (2014). 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 preparation. 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 (Wang, Bai & Goodman 2019; Gressel et al. 2020; Rodenkirch et al. 2020; Lesur 2021b), we derive a heat source based on microphysical consideration. We compute $\Gamma_s$, the average volumetric heating rate in the surface layer, at column densities between $10^{-19}$ and $10^{-20}$ cm$^{-2}$. Fig. 16 shows the trends of $\Gamma_s$, which is parametrized by the power law

$$ \Gamma_s \approx a \times 10^{-6} \left( \frac{R}{0.1 \text{ au}} \right)^{-b} \text{erg cm}^{-3} \text{s}^{-1}, \quad (36) $$

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 equation (35). Thermal MHD simulations show that strong heating rates as the one computed here lead to uncollimated outflows, see Meskini et al. (in preparation).

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 Hollenbach et al. (1994) 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 energized. 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 $\approx 1$ au (e.g. Ecolano & Picogna 2022) there is also a non-zero flow inside the launching radius (e.g. Dullemond et al. 2007).

We show here that the ionizing 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 (e.g. Bai et al. 2016).

# 4 DISCUSSION

## 4.1 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 $\Gamma_{acc}$. For a thin disc, this last quantity can be expressed as (Pringle 1981),

$$ \Gamma_{\text{acc}} = \frac{9}{4} \alpha \rho c_s^2, \quad (37) $$

where $\rho$ is the gas mass density and $c_s$ is the sound speed. In this framework, the parameter $\alpha$ links the viscous heating to the gas pressure ($P_g = \rho 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 $\alpha$ for determining thermal structure. Former studies, such as those by D'Alessio et al. (1998, 1999) and Glassgold, Najita & Igea (2004), have investigated the impact of different values of $\alpha$ on the thermal properties of the gas, commonly based on models with $\alpha \approx 0.01$. In the aforementioned models, X-
<figure><img src="image_15.png" /><figcaption>**Figure 15.** Expected distribution of the volumetric heating rate as a function of the disc radius, R, and the normalized height, Z/R, as calculated in PRODIMO in the case of SSF (panel a), and SDF (panel b). The white solid contour lines show the total vertical column density.</figcaption></figure>

<figure><img src="image_16.png" /><figcaption>**Figure 16.** Average volumetric heating rate at the surface layer in the range $10^{19}-10^{20}$ cm$^{-2}$. We consider the SSF model (green) and the SDF model (blue). The solid lines plot the fits for the volumetric heating rate (equation 36).</figcaption></figure>

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. First, heating is released through the ionization of the disc atoms and molecules (see Fig. 15). Secondly, these EPs alter the chemical and ionization structure, thereby increasing the disc viscosity (see Fig. 14). This increased viscosity triggers an increase of the accretion rate, leading to enhanced accretion heating (see equation 37).

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 ionization levels for triggering the MRI.

## 4.2 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 \sim 1 \left( \frac{M_*}{M_{\odot}} \right)^{-1/2} \left( \frac{R}{1 \text{ au}} \right)^{3/2} \text{yr.} \quad (38)$$

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_{\text{diff}} = \frac{R^2}{v}$, is

$$t_{\text{diff}}(R) \sim 10^4 \left( \frac{0.01}{\alpha} \right) \left( \frac{0.1}{h/R} \right)^2 t_K(R). \quad (39)$$

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

Moreover, the thermal time-scale defined from the sound speed $c_s = \alpha^{-1}v_{\text{turb}}$ falls between the dynamical and diffusion times,

$$t_{\text{th}}(R) \sim 10^2 \left( \frac{0.01}{\alpha} \right) t_K(R), \quad (40)$$

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 ionization over periods larger or similar to $t_{\text{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 time-scales. Accretion triggered by the heating pulse of a flare, as simulated by Orlando et al. (2011) and Colombo et al. (2019), may play a role on the accretion variability on time-scale of months, which is the typical time-scale of superflares as observed by Getman & Feigelson (2021).

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 time-scales of months.
## 4.3 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 photodissociation and photoionization by various forms of radiation, i.e. FUV, EUV, and X-rays. We have shown that in this surface layer, EPs produced by flares may be an important additional ionizing and heating source.

Above the disc, magnetized 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 photoheating is less effective, as shown by various studies (e.g. Garcia et al. 2001; Panoglou et al. 2012; Wang et al. 2019). Based on the results presented in this article, which focuses on discs, we emphasize 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.

## 4.4 Dependence of the results to the model parameters

### 4.4.1 Waiting-time distribution

In this section, we explore how the distribution of waiting intervals between two flares affects the ionization rate. The aim of this analysis is to assess the possible variation in ionization 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 Section 2.2.1. 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 due to observations of flares on young stars made by Getman & Feigelson (2021) and Getman et al. (2021), 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 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 ionization 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 ionization rate.

### 4.4.2 EP 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 normalized the number of supra-thermal particles entering the disc. The uncertainty in our approach is that we are not certain whether 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 three orders of magnitude compared to that estimated in B23, the ionization rate from particles produced by flares would be negligible compared to the ionization rate by the star X-ray radiation. In their study, Fleishman et al. (2022) 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 Section 4.5.1, we propose a numerical approach that could allow for a more accurate estimation of the density of supra-thermal particles penetrating the disc.

## 4.5 Model limitations

### 4.5.1 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^{25}$ cm$^{-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 (Padovani, Galli & Glassgold 2009). Diffusive processes need to be considered. For a more detailed discussion, see Padovani et al. (2018).

Secondly, 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 on to the stellar surface, and propagate into the outflows.

Lastly, 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 Marcowith et al. (2020) and Ji et al. (2022) 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, dimin-
ishes 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 ionization rates.

While the literature discusses the possibility whether superflares are sufficiently long to connect the star with the disc (Favata et al. 2005; Hamaguchi et al. 2012) or not (Getman et al. 2008; Getman & Feigelson 2021), we highlight the potential for particles to be produced by magnetic reconnection events resulting from Keplerian shearing. Indeed, Fatuzzo et al. (2023) 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 characterized 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 Fatuzzo et al. (2023) to draw analogies between magnetic reconnection processes in inner circumstellar discs and solar flares, a significant assumption of our study.

To summarize, 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 an MHD code in the test-particle limit or beyond, i.e. including a back reaction term.

### 4.5.2 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 PRODIMO. For a more complete picture, it would be interesting to examine the impact of the particles by varying the structure of PRODIMO 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 ionization rate from particles. The PRODIMO 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 Kamp et al. (2017). Therefore, the rise in the abundance of HCNH⁺ might not be physical but rather a dead-end in the reaction network calculated by PRODIMO. 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 remain 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 EPs from flares are included in the model.

Another limitation concerns the calculation of the viscous parameter $\alpha$ inside the disc. Although $\alpha$ 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 $\alpha$ 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 ionization by flares. From the $\alpha$ calculation, we have also estimated the mass accretion rate, taking ionization 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 Shakura & Sunyaev (1973), applicable far from the disc inner edge, a criterion not necessarily met in our study.

## 5 CONCLUSIONS

The ionization of protoplanetary discs is crucial for the triggering of MHD instabilities, responsible for the accretion of matter on to 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 (B23), we examined the effects of an additional ionization source based on microphysical and observational considerations. A stationary state of the particle emission was assumed leading to an overestimation of the ionization rate. In this paper, we address this issue by incorporating time-dependent components into our model.

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

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

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

*Impact on chemistry:* The ionization fraction increases by at least an order of magnitude up to $N \sim 10^{25} \text{ cm}^{-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 \approx 0.12$. The recombination of this cation leads to an increase of a factor of 2–3 in the column density of hydrogen cyanide, HCN. This molecule is crucial in pre-biotic chemistry to understand the synthesis 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 EPs produced by flares penetrate.
**Impact on the volumetric heating rate:** We determined the volumetric heating rate resulting from the ionization of disc matter by EPs 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 (equation 36).

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 PRODIMO. However, multiple aspects related to synthetic *JWST* spectra production are currently under development within PRODIMO. 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 EPs produced by turbulence-induced reconnection events.

Additionally, until now, our focus has been on the influence of EPs on the ionization rates, which subsequently affect the disc dynamics and chemistry. However, EPs can also collide with stable isotopes, creating radioactive isotopes, such as $p + ^{25}Al \rightarrow ^{26}Al$. The increased abundance of these radioactive species acts as a heat source within the disc. Previous work (e.g Gaches et al. 2020) 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). ChR 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.

## REFERENCES

Alexander R., Pascucci I., Andrews S., Armitage P., Cieza L., 2014, in Beuther H., Klessen R. S., Dullemond C. P., Henning T., eds, Protostars and Planets VI. p. 475

Aschwanden M. J., McTiernan J. M., 2010, ApJ, 717, 683

Aschwanden M. J., Johnson J. R., Nurhan Y. I., 2021, ApJ, 921, 166

Bai X.-N., Ye J., Goodman J., Yuan F., 2016, ApJ, 818, 152

Benz A. O., 2002, Phys. Plasmas, 5, 68

Benz A. O., 2017, Living Rev. Sol. Phys., 14, 2

Blandford R., Payne D., 1982, MNRAS, 199, 883

Brunn V., Marcowith A., Sauty C., Padovani M., Rab C., Meskini C., 2023, MNRAS, 519, 5673

Casse F., Ferreira J., 2000, A&A, 361, 1178

Cleeves L. I., Adams F. C., Bergin E. A., 2013, ApJ, 772, 5

Colombo S., Orlando S., Peres G., Reale F., Argiroffi C., Bonito R., Ibgui L., Stehlé C., 2019, A&A, 624, A50

D'Alessio P., Cantò J., Calvet N., Lizano S., 1998, ApJ, 500, 411

D'Alessio P., Calvet N., Hartmann L., Lizano S., Cantò J., 1999, ApJ, 527, 893

Dalgarno A., Yan M., Liu W., 1999, ApJS, 125, 237

Desch S. J., Turner N. J., 2015, ApJ, 811, 156

Donati J.-F. et al., 2007, MNRAS, 380, 1297

Dullemond C. P., Hollenbach D., Kamp I., D'Alessio P., 2007, in Reipurth B., Jewitt D., Keil K., eds, Protostars and Planets V. University of Arizona Press, Tucson, p. 555

Emslie A. et al., 2012, ApJ, 759, 71

Ercolano B., Picogna G., 2022, Eur. Phys. J. Plus, 137, 1357

Ercolano B., Clarke C. J., Drake J. J., 2009, ApJ, 699, 1639

Fatuzzo M., Adams F. C., Feinstein A. D., Seligman D. Z., 2023, ApJ, 954, 15

Favata F., Flaccomio E., Reale F., Micela G., Sciortino S., Shang H., Stassun K. G., Feigelson E. D., 2005, ApJS, 160, 469

Feigelson E. D., Montmerle T., 1999, ARA&A, 37, 363

Feigelson E. D., Garmire G. P., Pravdo S. H., 2002, ApJ, 572, 335

Ferreira J., 1997, A&A, 319, 340

Flaccomio E., Damiani F., Micela G., Sciortino S., F. R. Harnden J., Murray S. S., Wolk S. J., 2003, ApJ, 582, 398

Fleishman G. D., Nita G. M., Chen B., Yu S., Gary D. E., 2022, New Astron., 606, 674

Fromang S., Terquem C., Balbus S. A., 2002, MNRAS, 329, 18

Gaches B. A. L., Walch S., Offner S. S. R., Münker C., 2020, ApJ, 898, 79

Garcia P. J., Ferreira J., Cabrit S., Binette L., 2001, A&A, 377, 589

Getman K. V., Feigelson E. D., 2021, ApJ, 916, 32

Getman K. V., Feigelson E. D., Broos P. S., Micela G., Garmire G. P., 2008, ApJ, 688, 418

Getman K. V., Feigelson E. D., Garmire G. P., 2021, ApJ, 920, 154

Glassgold A. E., Najita J., Igea J., 2004, ApJ, 615, 972

Glassgold A. E., Galli D., Padovani M., 2012, ApJ, 756, 157

Gregory S., Matt S., Donati J.-F., Jardine M., 2008, MNRAS, 389, 1839

Gressel O., Ramsey J. P., Brinch C., Nelson R. P., Turner N. J., Bruderer S., 2020, ApJ, 896, 126

Hamaguchi K., Grosso N., Kastner J. H., Weintraub D. A., Richmond M., Petre R., Teets W. K., Principe D., 2012, ApJ, 754, 32

Hartmann L., Herczeg G., Calvet N., 2016, ARA&A, 54, 135

Hollenbach D., Johnstone D., Lizano S., Shu F., 1994, ApJ, 428, 654

Igea J., Glassgold A. E., 1999, ApJ, 518, 848

Ilgner M., Nelson R. P., 2006, A&A, 455, 731

Ivlev A., Padovani M., Galli D., Caselli P., 2015, ApJ, 812, 135

Jacquemin-Ide J., Ferreira J., Lesur G., 2019, MNRAS, 490, 3112

Jacquemin-Ide J., Lesur G., Ferreira J., 2021, A&A, 647, A192

Jannaud T., Zanni C., Ferreira J., 2023, A&A, 669, A159

Ji H., Daughton W., Jara-Almonte J., Le A., Stanier A., Yoo J., 2022, Nat. Rev. Phys., 4, 263

Kamp I., Tilling I., Woitke P., Thi W. F., Hogerheijde M., 2010, A&A, 510, A18

Kamp I., Thi W. F., Woitke P., Rab C., Bouma S., Ménard F., 2017, A&A, 607, A41

Lesur G., 2021a, J. Plasma Phys., 87, 205870101

Lesur G. R., 2021b, A&A, 650, A35

Lesur G., Kunz M. W., Fromang S., 2014, A&A, 566, A56

Long F. et al., 2021, A&A, 647, A118

Malzac J., 1999, Theses, Université Paul Sabatier – Toulouse III

Manara C. F. et al., 2016, A&A, 591, L3

Marcowith A., Ferrand G., Grech M., Meliani Z., Plotnikov I., Walder R., 2020, Living Rev. Comput. Astrophys., 6, 1

Matthews S. et al., 2021, Exp. Astron., 54, 1

Mattia G., Fendt C., 2020a, ApJ, 900, 59

Mattia G., Fendt C., 2020b, ApJ, 900, 60

Mattia G., Fendt C., 2022, ApJ, 935, 22

McCleary J. E., Wolk S. J., 2011, AJ, 141, 201

Notsu Y. et al., 2019, ApJ, 876, 58

Oberg N., Kamp I., Cazaux S., Woitke P., Thi W. F., 2022, A&A, 667, A95
Offner S. S. R., Gaches B. A. L., Holdship J. R., 2019, ApJ, 883, 121

Oka M. et al., 2018, Space Sci. Rev., 214, 82

Okamoto S., Notsu Y., Maehara H., Namekata K., Honda S., Ikuta K., Nogami D., Shibata K., 2021, ApJ, 906, 72

Orlando S., Reale F., Peres G., Mignone A., 2011, MNRAS, 415, 3380

Padovani M., Galli D., Glassgold A. E., 2009, A&A, 501, 619

Padovani M., Ivlev A. V., Galli D., Caselli P., 2018, A&A, 614, A111

Padovani M. et al., 2022, A&A, 658, A189

Panoglou D., Cabrit S., Des Forêts G. P., Garcia P., Ferreira J., Casse F., 2012, A&A, 538, A2

Perez-Becker D., Chiang E., 2011, ApJ, 735, 8

Pringle J. E., 1981, ARA&A, 19, 137

Rab C., Güdel M., Padovani M., Kamp I., Thi W. F., Woitke P., Aresu G., 2017, A&A, 603, A96

Rab C., Güdel M., Woitke P., Kamp I., Thi W. F., Min M., Aresu G., Meijerink R., 2018, A&A, 609, A91

Ray T. P., Ferreira J., 2021, New Astron. Rev., 93, 101615

Reale F., 2007, A&A, 471, 271

Reale F., 2014, Living Rev. Sol. Phys., 11, 4

Ripperda B., Porth O., Xia C., Keppens R., 2017, MNRAS, 471, 3465

Rodenkirch P. J., Klahr H., Fendt C., Dullemond C. P., 2020, A&A, 633, A21

Rodgers-Lee D., Taylor A. M., Ray T. P., Downes T. P., 2017, MNRAS, 472, 26

Sauty C., 2019, JET Simulations, Experiments, and Theory. Springer 59–65

Schlickeiser R., 2013, Cosmic ray astrophysics. Springer Science and Business Media 105–143

Semenov D., Wiebe D., Henning T., 2004, A&A, 417, 93

Serio S., Reale F., Jakimiec J., Sylwester B., Sylwester J., 1991, A&A, 241, 197

Shakura N. I., Sunyaev R. A., 1973, A&A, 24, 337

Shibata K., Magara T., 2011, Living Rev. Sol. Phys., 8, 6

Taillades J., Beuzelin I., Garrel L., Tabacik V., Bied C., Commeyras A., 1998, *Orig. Life Evol. Biosph.*, 28, 61

Thi W. F., Woitke P., Kamp I., 2011, MNRAS, 412, 711

Thi W. F., Lesur G., Woitke P., Kamp I., Rab C., Carmona A., 2019, A&A, 632, A44

Turner N. J., Sano T., Dziourkevitch N., 2007, ApJ, 659, 729

Venuti L. et al., 2014, A&A, 570, A82

Waggoner A. R., Cleeves L. I., 2022, ApJ, 928, 46

Walsh C., Nomura H., van Dishoeck E., 2015, A&A, 582, A88

Wang L., Bai X.-N., Goodman J., 2019, ApJ, 874, 90

Wardle M., Salmeron R., 2012, MNRAS, 422, 2737

Waterfall C., Browning P., Fuller G., Gordovskyy M., Orlando S., Reale F., 2020, MNRAS, 496, 2715

Woitke P., Kamp I., Thi W.-F., 2009, A&A, 501, 383

Woitke P. et al., 2016, A&A, 586, A103

Woitke P., Thi W. F., Arabhavi A. M., Kamp I., Kospal A., Abraham P., 2024, A&A, 683, A219

Wolk S. J., F. R. Harnden J., Flaccomio E., Micela G., Favata F., Shang H., Feigelson E. D., 2005, ApJS, 160, 423

Zanni C., Ferreira J., 2013, A&A, 550, A99

Ziegler J. F., Ziegler M. D., Biersack J. P., 2010, Nucl. Instrum. Methods Phys. Res., Sect. B, 268, 1818

## APPENDIX A: NORMALIZING 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) \frac{v(E)}{4\pi}, \qquad (A1)$$

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 normalization $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_{\text{th}}(E) = n_{\text{th}} \frac{2}{\sqrt{\pi}} \sqrt{\frac{E}{kT}} \frac{1}{kT} \exp(-E/kT). \qquad (A2)$$

Here, the peak value is at $E = 3/2kT \equiv E_{\text{th}}$.

As it is usually assumed in reconnection models for solar flares (Ripperda et al. 2017) and supported by observation of solar flares (Emslie et al. 2012; Matthews et al. 2021) 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,\text{nt}} = U_{p,\text{nt}} = \int_{E_c}^{\infty} E F_{p}(E) dE = \int_{E_c}^{\infty} E F_{e}(E) dE. \qquad (A3)$$

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_{\text{nt}}(E). \qquad (A4)$$

As a result of the magnetic reconnection process, we introduce the non-thermal component of the energy distribution function $F_{\text{nt}}$. The non-thermal energy distribution is determined by its normalization $N_{\text{nt}}$. In the following text, 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 $\delta$. 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 = \theta E_{\text{th}}$. Most Tauri flare models exhibit a single power law. We assume an exponential cut-off beyond $E_U$,

$$F_{\text{nt}}(E) = N_{\text{nt}} \left(\frac{E}{E_c}\right)^{-\delta} \exp\left(-\frac{E}{E_U}\right), \qquad (A5)$$

The normalization factor $N_{\text{nt}}$ can be expressed in terms of the non-thermal particle density,

$$n_{\text{nt}} = \int_{E_c}^{\infty} \left(\frac{E}{E_c}\right)^{-\delta} \exp\left(-\frac{E}{E_U}\right) dE \approx N_{\text{nt}} \frac{E_c}{\delta - 1}. \qquad (A6)$$

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

$$n_e(T) = 20.5 \times 10^{10} \left(\frac{T}{1 \text{ MK}}\right)^{-0.69} \text{cm}^{-3} \qquad (A7)$$

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

$$n_{\text{nt}} = \frac{3\theta}{\delta - 1} \sqrt{\frac{3\theta}{2\pi}} \exp\left(-\frac{3}{2}\theta\right) n_{\text{th}}. \qquad (A8)$$

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 equation (A7) to deduce the thermal and non-thermal densities. Given that $n = n_{\text{nt}} + n_{\text{th}}$, we get

$$n = n_{\text{th}} \left( 1 + \frac{3\theta}{\delta - 1} \sqrt{\frac{3\theta}{2\pi}} \exp\left(-\frac{3}{2}\theta\right) \right). \qquad (A9)$$
From equations (A9) and (A7), the thermal and non-thermal particle densities can be expressed as a function of the observed temperature,

$$n_{\text{th}} = 20.5 \times 10^{10} \frac{\left(\frac{T}{1 \text{MK}}\right)^{-0.69}}{1 + \frac{3\theta}{\delta-1} \sqrt{\frac{3\theta}{2\pi}} \exp\left(-\frac{3\theta}{2}\right)} \text{cm}^{-3}. \quad (\text{A10})$$

Assuming the fiducial values $\theta = 3$ and $\delta = 3$,

$$n_{\text{th}} = 19 \times 10^{10} \left(\frac{T}{1 \text{MK}}\right)^{-0.69} \text{cm}^{-3}. \quad (\text{A11})$$

The non-thermal particle density $n_{\text{nt}}$ is then recovered from equations (A10) and (A8),

$$n_{\text{nt}} = 20.5 \times 10^{10} \frac{\left(\frac{T}{1 \text{MK}}\right)^{-0.69}}{1 + \frac{\delta-1}{3} \sqrt{\frac{2\pi}{3\theta}} \exp\left(\frac{3}{2\theta}\right)} \text{cm}^{-3} \quad (\text{A12})$$

Assuming the fiducial values $\theta = 3$ and $\delta = 3$,

$$n_{\text{nt}} = 10^{10} \left(\frac{T}{1 \text{MK}}\right)^{-0.69} \text{cm}^{-3}. \quad (\text{A13})$$

We highlight here the difference between equation (A9) and equation (15) of B23. In that article, equation (15) is non-homogeneous. There is obviously a confusion between the normalization factor of the non-thermal flux and the non-thermal particle density. We corrected the issue in this section. The wrong estimation of the normalization factor $N_{\text{B23}}$ lead to an overestimation of the flux in B23. We compare the normalization $N_{\text{B23}}$ with the corrected normalization factor $N_{\text{nt}}$,

$$\frac{N_{\text{B23}}}{N_{\text{nt}}} = \frac{E_c}{kT(\delta-1)} = \frac{3\theta}{2(\delta-1)}. \quad (\text{A14})$$

Taking the fiducial values of B23, $\delta = 3$ and $\theta = 3$, $N_{\text{B23}} \approx 2N_{\text{nt}}$, the fluxes and ionization rates are reduced by a factor of 2 compared to B23.

## APPENDIX B: COLLISIONAL IONIZATION

In the chemical network used in Rab et al. (2017), collisional ionization, which can dominate the ionization fraction in the innermost high density and high temperature ($T > 1000$ K) (e.g. Fromang et al. 2002; Desch & Turner 2015) 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 + \text{Na} \rightarrow H_2 + \text{Na}^+ + \text{e}^- \quad (\text{B1})$$

Usually, gas-phase potassium K is used because of its low ionization potential (see Desch & Turner 2015). However, K is not included in

<figure><img src="image_17.png" /><figcaption>**Figure B17.** Comparison of the ionization fraction for different models and (semi-)analytical expressions in the disc mid-plane. The coloured solid lines show the different thermo-chemical models. In blue, the model from Rab et al. (2017), in green the same model but including accretion heating and in green including accretion heating (AH) and our approximation for collisional ionization (CI). The red dashed line shows the mid-plane gas temperature considering accretion heating. The dotted and dash-dotted grey lines show the ionization fraction using equation (13) of Fromang et al. (2002) for the no-metals and metal-dominated cases. The solid black line shows the results of solving the Saha equation. The vertical dotted line indicates $T = 1000$ K.</figcaption></figure>

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 ionization potential of K to calculate the rate coefficient using equation (6) of Desch & Turner (2015).

To verify our approximation, we compare our models to various analytical expressions in Fig. B17. By simply including accretion heating, which increases the temperature in the mid-plane, the ionization fraction increases by two orders of magnitude at $R \approx 0.1$ au. The inclusion of collisional ionization increases by two orders of magnitude and expands the high ionization region radially. This last model is in good agreement with the ionization fraction solely considering the Saha equation (using Na as the main element with the ionization potential of K). Close to the inner radius, our approximation underpredicts the ionization 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 ionization, the disc regions with $T = 1000$ K now have a high enough ionization rate to become MRI-active, which was not the case in the model of Rab et al. (2017) (see Section 3.3.2).

This paper has been typeset from a TeX/Latex file prepared by the author.