© The Authors 2026

1 Introduction

Young, low-mass pre-main-sequence stars in the T Tauri phase, typically a few million years old, continue to contract and to accrete material from their circumstellar discs. How they remove the excess angular momentum to enable accretion remains an open problem (Venuti et al. 2014; Hartmann et al. 2016; Manara et al. 2016). Most viable scenarios invoke magnetic fields coupled to partly ionised gas, through either turbulent torques driven by magneto-hydrodynamic (MHD) instabilities (Shakura & Sunyaev 1973), or laminar torques associated with magnetised winds and jets (Blandford & Payne 1982; Jacquemin-Ide et al. 2019). In this work, we focus on the consequences of MHD turbulence, where present in the disc. We restrict our attention to protoplanetary discs where the magneto-rotational instability (MRI) is likely the primary turbulence driver (Balbus & Hawley 1991; Balbus 2003). Crucially, in any accretion scenario where magnetic stresses regulate angular momentum transport, magnetic energy must be dissipated. In MRI-active layers, current sheets form so that dissipation through magnetic reconnection is a natural by-product of turbulent, magnetised discs (Pucci et al. 2024).

Reconnection can occur in several young stellar object (YSO) environments: in stellar coronae during flares (Benz 2017), in the star–disc interaction region near the truncation radius (Ferreira et al. 2000; Fatuzzo et al. 2023), and, central to this paper, within the disc’s MRI-active layers, where field-line turbulence continually regenerates current sheets (de Gouveia Dal Pino et al. 2010). Besides, sustaining MRI requires a minimum ionisation level. In YSO environments several sources of ionisation have already been investigated: UV and X-ray radiation either coming directly from the star or from the background radiation field (Glassgold et al. 1997), Galactic cosmic rays (GCRs), (Cleeves et al. 2013a; Padovani et al. 2018), radionuclide ionisation (Cleeves et al. 2013b), stellar energetic particles (SEP; Rab et al. 2017; Rodgers-Lee et al. 2017) and in-situ accelerated particles, in shocks (Padovani et al. 2016; Offner et al. 2019) or by magnetic reconnection in flares (Brunn et al. 2023, 2024). In the latter work, reconnection-accelerated EPs were shown to be an exceptionally efficient ionisation source, exceeding all others in their injection region (inner disc, R ≲ 1 au), which motivated the present work.

These sources are necessary to maintain sufficient disc ionisation to trigger MRI. Here we consider an additional, internal channel: non-thermal or energetic particles (EPs hereafter) produced in situ by turbulence-induced reconnection and their impact on the disc’s ionisation balance.

As MHD turbulence and magnetic reconnection develop in accretion discs, a part of the stored magnetic energy is possibly released in the form of non-thermal particles. Magnetic reconnection sites are known to accelerate non-thermal particles (see reviews in different astrophysical or space plasmas contexts by Marcowith et al. 2020; Oka et al. 2023; Guo et al. 2024). Particles can be accelerated in a variety of processes: Fermi-I like acceleration associated with the electric field induced in gas inflows (de Gouveia dal Pino & Lazarian 2005; Drury 2012), or due to curved magnetic fields. Acceleration can be due to the betatron effect in increasing magnetic zones under the conservation of the adiabatic particle moment, in magnetic contracting islands (Oka et al. 2010) or by direct electric field acceleration (Oka et al. 2023). Whatever the exact acceleration process, magnetic reconnection can be seen as a mechanism able to inject non-thermal particles into large scale turbulent flows as magnetic reconnection and turbulence are intimately inter-related.

In this work, we focus on the effect of non-thermal particles locally accelerated by the turbulence-induced magnetic reconnection. Rather than commit to one microphysical route, we adopt an empirical, Fermi-like parametrisation of the non-thermal tail and quantify the resulting ionisation and heating, treating the acceleration physics as an effective source term tied to the local turbulent–reconnection environment. A work in preparation, paper II, will address the consequences of this new ionisation source for the disc chemistry and chemical species abundances.

The paper is organised as follows. Sect. 2 outlines the disc thermo-chemical set-up and links it to turbulence and reconnection scalings. Sect. 3 introduces our EP injection and acceleration prescription as well as the transport set-up. Sect. 4 presents ionisation (and associated heating) maps and identifies where EPs dominate over X-rays and GCRs. Sect. 5 examines implications and limitations. Sect. 6 summarises the main findings and outlines future directions.

2 Physical Framework: from thermochemistry to turbulence-driven reconnection

In this section we introduce the physical framework. We first outline the disc thermo-chemical background, then summarise observational and theoretical constraints on turbulence, and finally present the reconnection model that sets the parameters for EP acceleration. The left panel of Fig. 1 sketches our scenario: a protoplanetary disc with a weakly ionised, MRI-stable mid-plane overlain by a magnetically coupled atmosphere where sub-Alfvénic MRI turbulence bends field lines and builds current sheets. Magnetic reconnection in these atmospheric layers accelerates supra-thermal particles that propagate quasi-ballistically along near-vertical field lines¹ towards the mid-plane. This section presents the quantitative arguments and assumptions that justify this scenario and defines the quantities used in the following model.

2.1 Protoplanetary disc model (PRODIMO)

PRODIMO² is a 2D radiation thermo-chemical code for protoplanetary discs (Woitke et al. 2009; Kamp et al. 2010; Thi et al. 2011; Woitke et al. 2016). It self-consistently couples disc structure, wavelength-dependent radiative transfer, gas and dust chemistry, and thermal balance to predict physical conditions and synthetic spectra over a broad wavelength range. The disc is treated on an axisymmetric grid (R, Z) for R ≈ 0.07–600 AU and Z/R = 0–0.5 with a resolution set to 100 × 100. Frequency-dependent dust continuum radiative transfer yields the local radiation field and dust temperature, including stellar and interstellar UV and stellar X-rays, enabling photochemistry and heating and cooling in the atmosphere. The chemical network spans 235 species and 3143 reactions (Kamp et al. 2017; Rab et al. 2017), with ice freeze-out and grain desorption (Thi et al. 2011). Heating processes include photoelectric emission, X-ray/cosmic-ray ionisation, viscous and collisional heating; cooling is treated via non-local thermal equilibrium line emission (e.g. CO, Fe II, O I), gas–dust coupling, and molecular lines. Iterative coupling ensures a consistent vertical structure (Arabhavi et al. 2022).

We adopted the fiducial T Tauri disc of Rab et al. (2017), with M_⋆ = 0.7M_⊙, T_⋆ = 4000K, L_⋆ = 1.0L_⊙, and an accretion rate of $\dot{M}$ = 10⁻⁸ M_⊙yr⁻¹, as is detailed in their Table 4, for our thermo-chemical background. This set-up has been widely used (Woitke et al. 2016; Kamp et al. 2017; Woitke et al. 2019). From the model, we extracted at each grid cell (R, Z): the total hydrogen nuclei density, n_H,tot = n_H + 2n_H2; the gas temperature, T; the electron fraction, x_e; the local composition (n_H, n_H2, n_He) used to characterise the turbulence and to build mixture-averaged energy-loss functions for electrons and protons; and the vertical column densities, $N(R,Z)\equiv {\int }_{0.5R}^Z{n}_{\mathrm{H},\mathrm{t}\mathrm{o}\mathrm{t}}(R,z)dz$, required for the transport model. These quantities set the injection EP energetics (Sect. 3), determine energy losses during propagation (Sect. 4), and yield the reference X-ray and CR ionisation rates used for comparison.

The top right panel of Fig. 1 displays the disc structure used in this work. As per PRODIMO, it shows a colour map of the total hydrogen number density, n_H,tot. Dashed black contours trace the vertical column density, N, at 10¹⁵, 10¹⁹, 10²¹, and 10²⁴ cm⁻². The dashed red curve marks the H⁺/H transition (ionised–atomic boundary), and the dashed orange curve marks the H/H₂ transition (atomic-molecular boundary). In the following, we adopt the following notation to designate the different regions of the disc: (i) atmosphere: 10¹⁵ ≤ N < 10¹⁹ cm⁻², highly ionised. (ii) Surface: 10¹⁹ ≲ N ≲ 10²¹ cm⁻², predominantly atomic gas. (iii) Intermediate layer: 10²¹ ≲ N ≲ 10²⁴ cm⁻², dominated by molecular gas. (iv) Deep layer: N ≳ 10²⁴ cm⁻², characterised by strong visual extinction (A_V ≳ 10).

2.2 Turbulence in protoplanetary discs

In this section, we summarise the key observational constraints on the amplitude and vertical stratification of turbulence in protoplanetary discs, and outline the parametrisation we adopt to embed MRI-like turbulence into our static thermo-chemical disc model.

Fig. 1 Left: schematic radial cut from R ~ 0.1 to ~50 AU presenting our scenario qualitatively. The disc is vertically stratified into a neutral mid-plane (blue; MRI-stable), where EPs arriving from above travel quasi-ballistically (CSDA-like) along near-vertical paths, and an ionised atmosphere (yellow; low-βp MRI-active), where sub-Alfvénic turbulence (MA < 1) triggers turbulent magnetic reconnection that injects EPs, which then stream downwards (black arrows). In the bottom left part, a patchy variant of the ionised atmospheric layer is sketched, with a lower covering area and possible lateral growth (black arrows). In the ionised turbulent layer, EP transport could be treated as stochastic. The inner star–disc interaction zone treated in Brunn et al. (2023, 2024) is indicated at R ≲ 0.3. Right: radial cut from R ~ 0.07 to ~600 AU displaying the disc chemical and turbulent structure used in this work based on the PRODIMO model. The top panel shows the total hydrogen number density, nH,tot (colour), with vertical column-density contours N = 1015, 1019, 1021, 1024 cm−2 (dashed black). The dashed red line marks the ionised–atomic transition (H+/H) defining the border between the region where H+ dominates and the region where H dominates. Analogously, the dashed orange line delimits the atomic–molecular transition (H/H2). The bottom panel displays the Alfvénic Mach number MA ≡ vturb/VA (colour). The red curve traces MA = 1. The sub-Alfvénic region (MA < 1) identifies where fast, turbulence-enabled reconnection, and thus EP acceleration, is expected. The solid red lines delimit the region where EP are efficiently accelerated, i.e. to energy more than 10 MeV. The white region is the so-called ‘dead zone’, where MRI is inefficient (see criterion Eq. (A.1)), there is no turbulence, MA is set to 0.

2.2.1 Observational constraints on turbulence in discs

Quantifying turbulence in discs relies primarily on (i) spectrally resolved molecular line observations that separate thermal from non-thermal broadening, (ii) the use of multiple tracers and isotopologues that probe different heights, and (iii) indirect constraints from dust settling. High spectral–spatial resolution ALMA data of CO and its isotopologues generally indicate weak turbulence in many systems, particularly in the outer disc (R ≳ 30 AU) and near the mid-plane, with typical non-thermal velocities bounded by v_turb ≲ 0.1c_s at radii of a few tens of astronomical units (e.g. Flaherty et al. 2017; Rosotti 2023). Independent inferences from the thickness of dust layers likewise point to low turbulent stirring. For instance, in HL Tau, dust vertical stratification implies v_turb ~10⁻²c_s at R ~ 100 AU (Pinte et al. 2016). However, a contrasting picture emerges in some sources and at higher altitudes above the mid-plane. Broader CO lines in DM Tau suggest v_turb ~0.2–0.3c_s at R ~ 20–50 AU (Flaherty et al. 2020). Moreover, chemically selective tracers that originate in the upper layers (e.g. CN, C₂H) reveal non-negligible turbulence at Z/R ≳ 0.2–0.3 and large radii in IM Lup (Paneque-Carreño et al. 2024). Recent ALMA CO analyses also in IM Lup, report on v_turb ~0.18–0.30c_s at R ~ 60 AU (Flaherty et al. 2024). Taken together, these measurements suggest a vertical gradient in the turbulent amplitude: turbulence is generally weakest near the mid-plane and stronger in the disc atmosphere, consistent with expectations for more ionised, MRI-active surface layers. This work focuses on a radial domain in the inner disc (R ≃ 0.3 to 30 AU) that is not constrained by ALMA. The turbulence is unknown. However, it is expected that within this region MRI is active (Lesur et al. 2023). The inner boundary reflects the zone of direct stellar magnetospheric influence due to stellar flaring loops, typically extending to R ~ 0.3–0.5 AU (Brunn et al. 2024).

2.2.2 Turbulence modelling in a static disc model

To evaluate the energy distribution of non-thermal particles, we first summarise the key properties of MHD turbulence in the region of interest, focusing on MRI-driven turbulence. MRI requires sufficient ionisation for strong gas–field coupling. The geometric set-up adopted here is sketched on the left part of Fig. 1.

Predicting the amplitude of turbulent velocity fluctuations remains challenging. The classical Shakura & Sunyaev (1973) framework parametrises angular-momentum transport via an effective viscosity, $${\nu }_{\mathrm{e}\mathrm{f}\mathrm{f}}\equiv L{v}_{\text{turb }}={\alpha }_{\mathrm{e}\mathrm{f}\mathrm{f}}{c}_{s}H,$$(1)

where L is a characteristic turbulent length scale, v_turb the turbulent speed,c_s the sound speed, and H the disc scale height. α_eff < 1 encodes transport efficiency.

The velocity amplitude depends both on α_eff and on how this efficiency partitions between spatial and velocity scales, as was discussed in Hughes et al. (2011). If eddies act on scales of ~H, one expects v_turb ~ α_effc_s; if α_eff distributes similarly across length and velocity, then v_turb $\sim \sqrt{{\alpha }_{\text{eff }}}{c}_s$. Empirically, observations and simulations rather support $$v_{\text{turb }}\approx \sqrt{{\alpha }_{\text{eff }}}{c}_s\Rightarrow L\approx \sqrt{{\alpha }_{\text{eff }}}H,$$(2)

which we adopt as the representative eddy size (Cuzzi et al. 2001; Hughes et al. 2011).

We evaluated α_eff using the non-ideal, MRI-based prescription of Thi et al. (2019), described in Appendix A. This yields v_turb ≈ 0.1–0.4c_s in the upper layers and v_turb ≲ 0.01c_s near the mid-plane, consistent with observational constraints, and provides the values of the L parameter encoding the reconnecting current sheet length and the turbulent speed used below in Section 2.3.2 to connect turbulence to reconnection geometry and, ultimately, to the non-thermal particle distribution.

2.3 Magnetic reconnection and turbulence in accretion flows

In this section, we show that the ionised, magnetised disc surface layers make MRI turbulence an efficient driver of current-sheet formation and magnetic reconnection, and we define the reconnection parameters that govern particle acceleration.

2.3.1 Magnetohydrodynamic turbulence and magnetic reconnection

Turbulence by itself does not guarantee magnetic reconnection. However, the atmospheric layers of protoplanetary discs combine physical conditions that make reconnection a natural outcome of the dynamics. Indeed, in this region, irradiation by stellar far UV and X-rays maintains high ionisation fractions (≳0.1, see the H⁺/H transition Fig. 1) ensuring tight magnetic coupling (weak Ohmic and ambipolar slippage). Where coupling is good, the MRI efficiently sustains turbulence and Maxwell stresses. The canonical expectation is that MRI is most efficient for high β_p (>100); MRI in low-β_p require additional conditions to grow. The parameter β_p = P_gas/P_mag is the ratio of the gas pressure to the magnetic pressure. In particular, Kim & Ostriker (2000) showed that reducing the toroidal component of the magnetic field re-opens the MRI window even at low β_p. Recent simulations of disc surface layers are consistent with this picture. Under conditions of strong ionisation and suitable field geometry, MRI-driven turbulence naturally generates magnetic reversals, shears and then thin current sheets, the precursors of reconnection (e.g. Jacquemin-Ide et al. 2021; Rosenberg & Ebrahimi 2021), making this process frequent and energetically relevant in disc atmospheres.

There are only few works linking turbulence and magnetic reconnection in accretion discs. Kadowaki et al. (2018) explore the connection between turbulence and fast magnetic reconnection in accreting systems. The authors conduct 3D ideal MHD simulations in the shearing box approximation. Although ideal, the simulations allow magnetic reconnection thanks to numerical resistivity. Their set-up starts from an initial azimuthal magnetic field and an exponential profile of the gas density over a scale H. The azimuthal magnetic field component dominates the initial radial and vertical components. The initial plasma beta parameter, β_p, is set to 1,10 and 100 in a series of three different runs. The simulations follow the evolution of the magnetic components under the effect of the MRI but also of the Parker-Rayleigh-Taylor instability (PRTI, Parker 1966).

The PRTI together with MRI, sustains an α − Ω dynamo that reaches a quasi-steady state in which the characteristic reconnection speed is V_r ≈ (0.13 ± 0.09)V_A in the disc and (0.17 ± 0.10)V_A in the corona, and the typical current-sheet width is Δ_r ≈ 0.08H (disc) and 0.10H (corona), where V_A is the local Alfvén speed. We computed the scale height using the Keplerian thin disc relation H =c_s/Ω_K³.

We report these numbers to show that a simulation, which self-consistently evolves MRI turbulence and reconnection, within the turbulence regime relevant to our work (v_turb ≈ 0.1–0.4c_s), naturally results in fast magnetic reconnection. Moreover, by applying the static viscous framework adopted here (see Appendix A, Thi et al. 2019) to our PRODIMO thermochemical background, we estimated the local Alfvén speed, $V_A=c_s\sqrt{\frac{2}{{\beta }_p}}$, and the Alfvénic Mach number, $M_A\equiv \frac{v_{\text{turb}}}{V_A}$, wherec_s is the local sound speed and v_turb the turbulent velocity. From this, we derived the reconnection speed, V_r, from Lazarian & Vishniac (1999) scaling, a model particularly suitable in the context of turbulent protoplanetary disc regions. In sub- to trans-Alfvénic turbulence $(M_A\lesssim 1),V_r\sim V_A{M}_A^2$, we predict V_r ~ 0.16V_A at the disc surface at 1 AU from the star, in very good agreement with the Kadowaki et al. (2018) results discussed above. The bottom right panel of Fig. 1 summarises the turbulent state we derived here and used to set the reconnection parameters on the PRODIMO grid. It shows a colour map of M_A, with the dashed red lines marking the M_A = 1 contour. The low-β_p (≲1) sub-Alfvénic domain (M_A < 1), bounded below this line, is the region of interest for EP acceleration, as it is the regime in which fast, turbulence-enabled magnetic reconnection is expected. The white region is the so called ‘dead zone’, where MRI is inefficient, there is no turbulence, and M_A ≡ 0.

The values of {V_r, Δ_r, H, V_A, M_A} computed from PRODIMO thus provide a coherent, simulation-supported model from which we estimate the reconnection–driven EP acceleration parameters, as is detailed in Sect. 3.2.

2.3.2 Multi-scale magnetic reconnection modelling and particle acceleration

At scales larger than kinetic ones, magnetic reconnection is allowed through the breaking of magnetic frozen-in conditions by collisions, even when the plasma is fully ionised and neutrals are not present (see e.g. Sweet 1958; Furth et al. 1963). While at kinetic scales magnetic reconnection is known to be sustained by electron inertia or the electron pressure tensor and to be efficient (i.e. fast) in converting magnetic energy into particle acceleration, fluid models have struggled to provide an explanation for fast energy release. Different models have been proposed in recent years (see e.g. the discussion in Lazarian et al. 2020), of which the tearing instability on thin current sheets is the most promising (see e.g. Loureiro et al. 2007; Bhattacharjee et al. 2009; Pucci & Velli 2014). The survival of laminar current sheets, necessary for the onset of the tearing instability, remains controversial in the context of turbulent protoplanetary discs, and outside the scope of this work. Still, reconnection is observed in 2D and also 3D turbulent simulations even with sufficiently high magnetic Reynolds number (see Dong et al. 2022 and references therein). In the latter work, the formation of reconnection-produced magnetic flux ropes in 3D is shown to exhibit a complex morphology, differing significantly from the simpler island-like morphology in 2D. The authors note the similar size of the flux ropes between the 2D and 3D configurations, though the size of plasmoids is distributed over a wide variety of scales.

The length, L_r, of the largest current sheet forming in a turbulent region is assumed to be approximately equal to the local turbulent injection scale, i.e. L ≈ L_r; the largest current sheets would be destroyed by the turbulence. The thickness, Δ_r, is set by the scale over which a magnetic field line deviates from its original direction (see Lazarian & Vishniac 1999) or, adopting a laminar tearing model, the intrinsic scale at which reconnection becomes fast (Bhattacharjee et al. 2009; Loureiro et al. 2007; Pucci & Velli 2014). Fig. 2 schematically illustrates our working picture: sub-Alfvénic MRI turbulence (M_A < 1) at the injection scale, L (left), generates current sheets where turbulent magnetic reconnection operates (right). The reconnecting current sheet has a typical length, L_r, and thickness, Δ_r, with L_r = L.

In a generic magnetic reconnection geometry, EP acceleration is characterised by three parameters – the acceleration rate, α_acc, the escape time, τ_esc, and the injection time, τ_inj – which depend on the reconnection model and on the characteristic size and lifetime of the reconnecting current sheet (Lazarian & Vishniac 1999; Kowal et al. 2012; Lazarian et al. 2020; Xu & Lazarian 2023), $${\tau }_{\mathrm{i}\mathrm{n}\mathrm{j}}=\frac{L_r}{V_r},{\tau }_{\mathrm{e}\mathrm{s}\mathrm{c}}=\frac{L_r}{V_A}.$$(3)

For the purpose of this paper, we consider the developments of the turbulent-reconnection framework of Zhang et al. (2023) to be appropriate for the estimation of the reconnection parameters. Thus, the turbulent parameters {L, M_A} set the reconnection speed, $V_{\text{rec }}\sim V_A{M}_A^2$ (Lazarian & Vishniac 1999), and, as is shown below, the parameters of the EP acceleration model. The particles are considered to be accelerated multiple times within the turbulent layer, transitioning between different current sheets, and undergo Fermi acceleration through island compression.

The same assumption holds when assuming that the plasmoid instability or any fast tearing instability kicks in. In particular we consider here the ‘ideal’ tearing instability (see Pucci & Velli 2014). Then the acceleration rate is ${\alpha }_{\text{acc}}\sim t_A^{-1}$, where the proportionality constant depends on the topology of the reconnecting magnetic field (Pucci et al. 2018) and other non-ideal effects (see e.g. Pucci et al. 2017; Shi et al. 2019). In this work, we adopt the results of Xu & Lazarian (2023), giving $${\alpha }_{\mathrm{a}\mathrm{c}\mathrm{c}}=\frac{2V_r}{3{\mathrm{\Delta }}_r},{\mathrm{\Delta }}_r=L_r{M}_A^2,V_r=V_A{M}_A^2,$$(4)

and together with Eq. (3), $${\alpha }_{\mathrm{a}\mathrm{c}\mathrm{c}}=\frac{2}{3}{t}_A^{-1},{\tau }_{\mathrm{i}\mathrm{n}\mathrm{j}}=M_A^{-2}{t}_A,{\tau }_{\mathrm{e}\mathrm{s}\mathrm{c}}=t_A,$$(5)

where t_A ≡ L_r/V_A = Δ_r/V_r.

In conclusion, comparing collisional fluid modelling to Zhang et al. (2023) brings correction factors to the energy effectively converted within the turbulent layer in the disc’s atmosphere. The development of detailed particle acceleration models based on the interplay of reconnection and turbulence remains largely phenomenological, and various types of dimensional analysis may be adopted. Not knowing the details of the microscopic acceleration process at this stage, we introduce the efficiency parameter, κ, and defer a full treatment to future work.

Fig. 2 Sketch of MRI-driven turbulence and embedded turbulent magnetic reconnection. Left: eddies at the injection scale L in a sub-Alfvénic regime characterised by the Alfvénic Mach number MA. Right: reconnection layer produced by the turbulence, with length Lr = L and thickness Δr. The turbulence parameters {L, MA} set the reconnection, geometry parameters {Lr, Δr}, speed $V_{\text{rec }}\sim V_A{M}_A^2$, and EP acceleration parameters used in this work.

3 Methodology – energetic particle acceleration by turbulence-induced magnetic reconnection

Low-energy (≲1 GeV) particles can significantly enhance the ionisation rate in the disc, influencing both its chemistry and the coupling between the gas and magnetic fields (Rab et al. 2017; Padovani et al. 2018; Brunn et al. 2023, 2024). In order to study ionisation by non-thermal particles in protoplanetary discs, we start by building a model to determine the energy distribution of the supra-thermal particles produced during magnetic reconnection events.

3.1 Energy distribution from acceleration by a Fermi-like process

Magnetic reconnection can efficiently accelerate a fraction of the thermal particle population to supra-thermal energies through Fermi-like acceleration mechanism. Assuming such an acceleration process is active in the region of interest, the energy gain by a particle is proportional to the particle kinetic energy, $\dot{E}={\alpha }_{\text{acc}}E$, with α_acc being the acceleration rate. The equation for the energy distribution is (Guo et al. 2014, 2015; Oka et al. 2023) $$\frac{\mathrm{\partial }{F}_{\mathrm{n}\mathrm{t}}}{\mathrm{\partial }t}+\frac{\mathrm{\partial }}{\mathrm{\partial }ϵ}\left(\dot{E}{F}_{\mathrm{n}\mathrm{t}}\right)=\frac{F_{\mathrm{i}\mathrm{n}\mathrm{j}}}{{\tau }_{\mathrm{i}\mathrm{n}\mathrm{j}}}-\frac{F_{\mathrm{n}\mathrm{t}}}{{\tau }_{\mathrm{e}\mathrm{s}\mathrm{c}}},$$(6)

where F_nt is the particle energy distribution (erg⁻¹ cm⁻³), $\dot{E}$ is the particle energy gain rate (erg s⁻¹), τ_inj is the timescale for injection of the initial particle distribution, F_inj, into the accelerating region, and τ_esc is the escape timescale from the acceleration region.

Assuming that the initial upstream distribution is a Maxwellian distribution⁴, $$F_{\mathrm{i}\mathrm{n}\mathrm{j}}(ϵ)=\frac{2n_{\mathrm{i}\mathrm{n}\mathrm{j}}}{\sqrt{\pi }}\sqrt{ϵ}\exp (-ϵ),$$(7)

where F_inj is the upstream particle energy distribution in units of thermal energy, ϵ is the energy in units of thermal energy, ϵ = E/E_th, with E_th = 3/2k_BT, where k_B is the Boltzmann constant and T the temperature of the medium, and n_inj is the density of particles entering in the acceleration region.

The solution of Eq. (6) is a non-thermal distribution that can be written in units of thermal energy as $$F_{\mathrm{n}\mathrm{t}}(ϵ,t)=\frac{2n_{\mathrm{i}\mathrm{n}\mathrm{j}}}{\sqrt{\pi }{\alpha }_{\mathrm{a}\mathrm{c}\mathrm{c}}{\tau }_{\mathrm{i}\mathrm{n}\mathrm{j}}}{ϵ}^{-(1+\beta )}[{\mathrm{\Gamma }}_{3/2+\beta }\left(ϵe^{-{\alpha }_{\mathrm{a}\mathrm{c}\mathrm{c}}t}\right)-{\mathrm{\Gamma }}_{3/2+\beta }(ϵ)],$$(8)

where $\beta =\frac{1}{{\alpha }_{\text{acc}}{\tau }_{\text{esc}}}$ and Γ(3/2 + β, ϵ) is the upper incomplete Gamma function of order 3/2 + β.

In the next sections, we evaluate the acceleration parameters assuming a specific turbulent magnetic reconnection model and see how to normalise it based on the power released by the turbulence in the disc.

3.2 Energy distribution parameters in turbulent magnetic reconnection

The shape of the non-thermal energy distribution is set by the three parameters, α_acc, τ_esc, and τ_inj, as defined in Eq. (5), which depend on the characteristic size and lifetime of the reconnecting current sheet.

As reconnection proceeds, particles are continuously injected and accelerated within current sheets. The combination of sustained injection and Fermi-like energisation naturally builds a supra-thermal tail. We take the effective acceleration duration to be the stability time of a reconnecting current sheet, τ_CS, which in a turbulent medium we approximate by the eddy turnover time at the injection scale, ${\tau }_{\mathrm{C}\mathrm{S}}=\frac{L_r}{v_{\text{turb}}}$. In what follows, we evaluate the non-thermal distribution at t = τ_CS. This ‘end-of-burst’ snapshot captures the particle population just before the sheet is disrupted, provides a representative reference for intermittent reconnection in turbulence, and simplifies the analysis for this first, stationary treatment.

After a time, τ_CS, a power-law component emerges if the factor Γ_3/2+β (ϵe^{−α_accτ_CS}) − Γ_3/2+β(ϵ) in Eq. (8) is (approximately) energy-independent. This holds across the supra-thermal range provided ϵe^{−α_accτ_CS} ≲ 1. For ϵe^{−α_accτ_CS} ≳ 1, the spectrum bends away from a pure power law, defining a finite-time cut-off, $${ϵ}_{\mathrm{m}\mathrm{a}\mathrm{x},\mathrm{C}\mathrm{S}}=\exp \left({\alpha }_{\mathrm{a}\mathrm{c}\mathrm{c}}{\tau }_{\mathrm{C}\mathrm{S}}\right)=\exp \left(\frac{2}{3M_A}\right),$$(9)

which directly reflects Fermi-like energisation, with $\dot{ϵ}={\alpha }_{\text{acc}}ϵ$ acting for a duration, τ_CS.

A second limit comes from confinement (Hillas criterion): the Larmor radius must fit within the accelerator, giving $${ϵ}_{max,\mathrm{H}}=\frac{eB{\mathrm{\Delta }}_r}{E_{\mathrm{t}\mathrm{h}}},$$(10)

with B the magnetic field strength. The effective maximum energy is therefore $$E_{max}=min\{{ϵ}_{max,\mathrm{C}\mathrm{S}},{ϵ}_{max,\mathrm{H}}\}{E}_{\mathrm{t}\mathrm{h}}.$$(11)

At the low-energy end, we define the injection threshold, ϵ_inj, by continuity between thermal and non-thermal components, $$F_{\mathrm{i}\mathrm{n}\mathrm{j}}({ϵ}_{\mathrm{i}\mathrm{n}\mathrm{j}})=F_{\mathrm{n}\mathrm{t}}({ϵ}_{\mathrm{i}\mathrm{n}\mathrm{j}}),$$(12)

which we solve numerically at time t = τ_CS. We then refer to the non-thermal population as EPs with normalised energy ϵ ≥ ϵ_inj.

The slope of the supra-thermal tail follows directly from the acceleration and escape parameters, $$p=1+\frac{1}{{\alpha }_{\mathrm{a}\mathrm{c}\mathrm{c}}{\tau }_{\mathrm{e}\mathrm{s}\mathrm{c}}}=2.5.$$(13)

In this work, we take α_acc, τ_esc, and τ_inj to be the same for electrons and protons, so the derived p = 2.5 applies to both species. A species-dependent treatment is deferred to future work.

With the spectral shape fixed by α_acc, τ_inj, and τ_esc (thus setting p, ϵ_inj, and ϵ_max), the remaining task is to determine the normalisation, i.e. the number density of non-thermal particles, from the local energy budget. In the next subsection, we link the particle power to a fraction of the turbulent accretion dissipation to obtain the non-thermal EP density.

3.3 Normalisation of the energy distribution from the available energy of the turbulent accretion

We normalised the non-thermal distribution by assuming that the power sustaining it, ${\dot{U}}_{\mathrm{n}\mathrm{t}}$, is a fixed fraction, κ, of the local viscous dissipation per unit volume, D/H, $${\dot{U}}_{\mathrm{n}\mathrm{t}}\equiv {\int }_{E_{\mathrm{i}\mathrm{n}\mathrm{j}}}^{E_{max}}F(E)\dot{E}\mathrm{d}E=\kappa \frac{D}{H},$$(14)

where D is the viscous energy flux (erg s⁻¹ cm⁻²), H the disc scale height, and κ the (dimensionless) fraction channelled into non-thermal particle acceleration.

As we assumed a Fermi-like acceleration process, $\dot{E}={\alpha }_{\text{acc}}E$, with α_acc independent of energy, and thus $${\int }_{E_{\mathrm{i}\mathrm{n}\mathrm{j}}}^{E_{max}}F(E)E\mathrm{d}E=\kappa \frac{D}{{\alpha }_{\mathrm{a}\mathrm{c}\mathrm{c}}H},$$(15)

so the left-hand side is the non-thermal energy density, U_nt, by definition.

In the Keplerian, thin disc approximation, the viscous dissipation is (Pringle 1981; Martin et al. 2019) $$D(R)=\frac{1}{2}{\nu }_{\mathrm{e}\mathrm{f}\mathrm{f}}\mathrm{\Sigma }{\left(R\frac{\mathrm{d}\mathrm{\Omega }}{\mathrm{d}R}\right)}^2=\frac{9}{8}{\nu }_{\mathrm{e}\mathrm{f}\mathrm{f}}\mathrm{\Sigma }{\mathrm{\Omega }}^2,$$(16)

where ν_eff is defined in Eq. (1), Σ = Hρ is the accreting column density, ρ is the mass density, Ω is the angular velocity, and R is the distance to the star. Hence, $$\frac{D}{H}=\frac{9}{8}\rho c_s^2\frac{L_r{v}_{\text{turb }}}{H^2}=\frac{5}{4}{U}_{\text{th }}\frac{L_r{v}_{\text{turb }}}{H^2},$$(17)

where we usedc_s = HΩ and $U_{\text{th }}=\frac{3}{2}{P}_{\text{th}}$ with $\frac{5}{3}{P}_{\text{th}}=\rho c_s^2$, as the disc is composed of mono-atomic gas in the acceleration region. Substituting into Eq. (15) gives $$\frac{U_{\mathrm{n}\mathrm{t}}}{U_{\mathrm{t}\mathrm{h}}}=\kappa \frac{5}{4{\alpha }_{\mathrm{a}\mathrm{c}\mathrm{c}}}\frac{L_r{v}_{\mathrm{t}\mathrm{u}\mathrm{r}\mathrm{b}}}{H^2}.$$(18)

This condition fixes the normalisation. Hence the injected number density, n_inj, of the distribution of Eq. (8) is $$n_{\mathrm{i}\mathrm{n}\mathrm{j}}\gamma =\kappa \frac{5\sqrt{\pi }}{8}\frac{L_r{\tau }_{\mathrm{i}\mathrm{n}\mathrm{j}}{v}_{\mathrm{t}\mathrm{u}\mathrm{r}\mathrm{b}}}{H^2}{n}_{\mathrm{t}\mathrm{h}},$$(19)

where $\gamma ={\int }_{{ϵ}_{\text{inj }}}^{{ϵ}_{\text{max }}}{ϵ}^{-\beta }[{\mathrm{\Gamma }}_{3/2+\beta }\left(ϵ/{ϵ}_{\text{max }}\right)-{\mathrm{\Gamma }}_{3/2+\beta }(ϵ)]dϵ$ and n_th = U_th/E_th. Using τ_injv_turb = L_r/M_A (Eq. (5)) yields $$\frac{n_{\mathrm{i}\mathrm{n}\mathrm{j}}}{n_{\mathrm{t}\mathrm{h}}}=\kappa \frac{5\sqrt{\pi }}{8}{\gamma }^{-1}{\left(\frac{L_r}{H}\right)}^2{M}_A^{-1},$$(20)

where L_r/H follows Eq. (2), β = 1/(α_accτ_esc) = 1.5 from Eq. (5), M_A is computed locally (Appendix A), ϵ_inj is obtained by numerically solving Eq. (12), γ is likewise solved numerically, and κ remains the free energetics parameter.

Using the injected non-thermal density from Eq. (20), we computed the local non-thermal energy density, U_nt(R, Z). Because the thermochemical background, and thus the acceleration parameters, vary with the position, all quantities were evaluated cell by cell on the PRODIMO (R,Z) grid. Fig. 3 maps the resulting U_nt/U_th and marks, in red, the regions where the local maximum particle energy exceeds 10 MeV (efficient acceleration⁵), while the black contour traces N = 10¹⁹ cm⁻² (the disc surface). We set the non-thermal power to κ = 0.1, by analogy with the canonical ~10–30% efficiency with which supernova-remnant shocks are thought to channel kinetic energy into GCR (Ptuskin et al. 2010).

On Fig. 3, the sharp drop from U_nt > 0 to U_nt = 0 on the disc side tracks the M_A = 1 contour. In our prescription, the maximum non-thermal energy ϵ_max tends to unity for M_A ≳ 1 (Eq. (9)), the supra-thermal tail merges with the thermal distribution. We therefore set U_nt = 0 for M_A ≳ 1, noting that this reflects the disappearance of a distinct non-thermal particle population rather than a suppression of turbulent reconnection itself. Moreover, because the turbulent reconnection model is formulated for strongly magnetised, low-β_p plasmas (Lazarian & Vishniac 1999), and simulations indicate that as β_p increases, the non-thermal particle acceleration efficiency declines (Li et al. 2017). We restricted EP acceleration to low-β_p regions and set U_nt = 0 for β > 10. We note that efficient non-thermal acceleration may occur in hot or relativistic high-β plasmas, yet such regimes are not relevant for protoplanetary discs, and therefore lie outside the scope of this work.

On the atmosphere side, at very low M_A, U_nt also vanishes because the reconnection layer thickness scales as ${\mathrm{\Delta }}_r\sim M_A^2$ (Eq. (4)), becoming extremely thin. The maximum particle energy is then limited by the Hillas criterion. In practice, we set U_nt = 0 when E_max ≲ 10E_th.

Notably, the map shows that our turbulence-driven reconnection model can accelerate particles to mega-electronvolt energies in extended disc region, out to a few tens of astronomical units along the ionised atmospheric layers, i.e. within the intermediate window where M_A < 1 yet not so small that Δ allow E_max ≫ E_th.

Throughout the acceleration region, U_nt is orders of magnitude smaller than the thermal energy density, U_th, and, since β_p ≲ 1, also far below the magnetic energy density. Consequently, EPs are not expected to exert a dynamical back-reaction over the plasma. This validates our treatment in the next section where EP transport is computed in the test-particle limit as a first approximation. In the next section (i.e. Sect. 4.1), the non-thermal energy distribution is computed from this normalisation, together with the acceleration parameters at each PRODIMO grid cell in order to model its propagation through the disc.

Fig. 3 Nonthermal–to–thermal energy density ratio, Unt/Uth, for κ = 0.1. The red contour encloses regions where the local maximum particle energy exceeds 10 MeV (efficient acceleration sites capable of supplying ionising EPs). The black contour marks a vertical column density, N = 1019 cm−2 (approximate disc surface).

Fig. 4 Particle energy distributions at R = 1 AU for κ = 0.1. The dashed curves show the thermal Maxwellian based on PRODIMO density and temperature for Z/R = 0.25–0.35. The solid curves show the non-thermal distribution produced by magnetic-reconnection acceleration.

4 Results: Ionisation from reconnection-accelerated particles

To quantify the ionisation rate in the disc due to EPs generated by turbulence induced magnetic reconnection, we first computed their initial energy distribution, propagation and energy attenuation through the disc material.

4.1 Non-thermal particle energy distributions: Injection and propagation

For each PRODIMO −R,Z grid cell, the local initial (‘injection’) non-thermal energy distribution $F_{\mathrm{n}\mathrm{t}}^{R,Z}$(E), was computed with Eq. (8) at time t = τ_CS, using the normalisation of Eq. (20), together with the acceleration parameters α_acc, τ_esc, and τ_inj. Fig. 4 shows representative energy distributions at R = 1 AU within the height range where U_nt > 0: the thermal Maxwellian (dashed) based on PRODIMO density and temperature, and the associated non-thermal tail (solid) produced by our model of magnetic reconnection in turbulence. We then propagated this local energy distribution through the disc material with the continuous–slowing–down approximation (CSDA) to obtain, for any traversed column density, N, the propagated distributions, $F_{\mathrm{n}\mathrm{t}}^{R,Z}$(E, N), for both electrons and protons, as in Brunn et al. (2023). We assume isotropic injection and approximate trajectories as vertical, which is adequate to first order in the subsonic–to–transonic turbulent regime considered here. A more complete treatment that explicitly includes turbulence, accounting for spatial diffusion along turbulent field lines and for stochastic re-acceleration in MRI turbulence (e.g. Sun & Bai 2021), is deferred to a forthcoming paper. Methodological details on loss functions, cross sections, CSDA validity, and the treatment of primaries and secondaries are provided by Padovani et al. (2009, 2018) and Brunn et al. (2023).

4.2 Ionisation rate

Energetic electrons and protons, produced in the turbulence-induced magnetic reconnection, propagate into the protoplanetary disc and interact with the ambient gas, leading to the ionisation of hydrogen and molecular hydrogen. To compute the ionisation rate, we consider both primary ionisation reactions and secondary processes⁶ triggered by the ejected energetic electrons. The ionisation rate, ζ_R,Z(N), at a vertical column density, N, produced by the EP distribution, $F_{\mathrm{n}\mathrm{t}}^{R,Z}$(E), originating from a pixel at position R, Z, was computed as in Brunn et al. (2023).

The total EP-induced ionisation rate at radius R, ζ_R(N), was computed by vertically averaging the contributions from all reconnection-active pixels (cells with non-zero non-thermal energy density; the coloured pixels in Fig. 3). Each acceleration site at height Z_i within the radial slice contributes by ζ_R,Zi(N) to the total propagated ionisation profile. The total ionisation is $${\zeta }_{\mathrm{R}}(N)=\frac{1}{\mathrm{\Delta }Z}\sum _i{\zeta }_{R,Z_i}(N)\mathrm{\Delta }{Z}_i,$$(21)

where ΔZ_i is the thickness of pixel i and ΔZ = ∑_i ΔZ_i is the total vertical extent of the active layer at that radius (i.e the length of the coloured slice at radius R of Fig. 3).

We adopt a thickness-weighted vertical average because, at a fixed radius, multiple acceleration sites at different heights contribute simultaneously, and not equally, to the ionisation at depth N. Weighting each pixel’s contribution, ζ_R,Zi(N), by its vertical extent, ΔZ_i, and normalising by the total active thickness, ΔZ, approximates the continuous integral, ΔZ⁻¹ ∫ ζ_R,Z(N)dZ. This ensures that layers occupying a larger fraction of the reconnecting zone contribute proportionally more, eliminates biases due to grid resolution, and provides a surrogate for the time and ensemble average of inherently patchy, intermittent reconnection. We ignored EP ionisation from the opposite disc side because the propagated flux is suppressed before reaching the mid-plane at column density N ≳ 10²⁵ cm⁻².

Fig. 5 illustrates the result at R = 1 AU, for κ = 0.1. Coloured curves (shown with the ‘viridis’ colour scale) are the individual ζ_R,Zi(N) for EPs injected at different heights, the thick orange curve is the vertically averaged total from Eq. (21), and the black curve shows, is for comparison, the stellar X-ray ionisation of a typical T Tauri object with an X-ray luminosity of L_X = 10³⁰ erg s⁻¹. While X-rays dominate at low column density (10¹⁶–10²³ cm⁻²), EPs substantially enhance and exceed the ionisation at higher columns (above a few 10²³ cm⁻²). At column density N ≳ 10²⁵ cm⁻² (not shown in Fig. 5), the EP ionisation rate drops, and the GCRs component overtakes close to the disc mid-plain. The column density at which EPs, X-rays or GCR dominate ionisation depends on the distance to the star; we discuss this aspect in the next section.

4.3 Spatial ionisation rate distribution

Using Eq. (21), we computed the total EP–induced ionisation rate at every PRODIMO grid cell (R, Z), yielding the map in Fig. 6. The figure shows ζ_EP on a logarithmic colour scale from 10⁻²⁰ s⁻¹ (dark purple) to 10⁻⁸ s⁻¹ (bright yellow). Dashed black curves are isocontours.

Ionisation is strongest in the inner disc and upper layers, then weakens with increasing radius and approaching the mid-plane, reflecting both attenuation with column density and reduced acceleration efficiency farther from the star. The solid red contour marks where particles are accelerated above 10 MeV, i.e. the portion of the atmosphere that efficiently injects EPs. The white band at the top indicates very low column density (<10¹⁵ cm⁻²) above the disc atmosphere. Outflow propagation in that region is beyond our present scope. The white area beyond ~30 AU corresponds to radii where acceleration is too weak to produce >10 MeV particles, so EPs are stopped at low column density and the ionisation is negligible there.

Importantly, while ζ_EP peaks at the surface, EPs also penetrate to large column densities, sustaining significant ionisation (N ≳ 10²³ cm⁻²) where stellar X-rays are heavily attenuated. This deep reach positions EPs as the dominant ionisation source in otherwise shielded disc layers.

Fig. 5 Ionisation rates as function of the disc vertical column density at R = 1 AU for κ = 0.1. The lines in the viridis colour scale are the contribution to the ionisation rate coming from EPs accelerated in region at different Z/R. The thick orange line is the total ionisation rate corresponding to the weighted sum (Eq. (21)) of each local contribution. The black line is the ionisation rate from the stellar X-rays.

Fig. 6 EP ionisation rate map for κ = 0.1. The solid red contour outlines the region where EPs are accelerated. Dashed lines indicate isocontours of the EP ionisation rate. The upper white region corresponds to column densities below N = 1015 cm−2, indicating regions above the disc atmosphere.

4.3.1 EP-dominated ionisation regions

This section identifies the conditions and locations where EPs dominate over X-rays and GCRs ionisation, mapping EP-dominated domains for κ = 0.1 and 1. In line with Rab et al. (2017), we produce a dominance map, but concentrate on regions closer to the star than in their study.

Fig. 7 shows the fraction of the total ionisation rate produced by EPs, ζ_EP/ζ_tot, for our reference case κ = 0.1, with ζ_tot = ζ_EP + ζ_XR + ζ_GCR. The colour scale is logarithmic. Bright yellow marks regions where EPs contribute most. Solid black contours enclose zones where EPs dominate, i.e. ζ_EP > ζ_XR + ζ_GCR. The dashed black contour marks where EPs supply a non-negligible share of the ionisation rate, at least 5% of ζ_tot. The red contour outlines regions dominated by GCRs. Fig. 7 highlights the efficiency of EP ionisation even when only 10% of the local viscous dissipation is channelled into non-thermal particles. Even with this modest energy budget, EPs form a clear dominant band in the inner disc, between 0.1–2 AU in the intermediate to deep layer at Z/R ~ 0.1. Beyond the strictly dominant zones, large areas exceed the 5% threshold at distances between 0.1–20 AU from the star, in the disc atmosphere and in the intermediate to deep layer between Z/R ~ 0.8–0.15. These zones of EP dominance are largely determined by the EP energy and propagation depth. In particular, GCRs due to their higher energies, propagate at higher column density and dominate the ionisation close the disc mid-plane. While X-rays dominate in the surface to intermediate region, above the EP-dominated zone. The EPs are thus a significant contributor to the ionisation rate over extended regions, relevant for disc chemistry and dynamics.

Fig. 8 shows the spatial distribution of the EP contribution to the total ionisation rate under the extreme assumption that all local viscous dissipation is channelled to non-thermal particles, i.e. κ = 1. The ionisation map shows that EPs dominate in the disc atmosphere from ~0.2–40 AU, and also in intermediate to deep layers from ~0.08–20 AU. This figure highlights that EPs accelerated by turbulence-induced magnetic reconnection can dominate the ionisation in several regions of the disc, and if not dominant, contribute significantly (more than 5%) to the total ionisation rate in most of the disc at radii from 0.1 AU to tens of astronomical units.

The κ = 1 assumption is unrealistic, but we use it to illustrate the upper-limit spatial pattern that our model can produce. We note that a similar ionisation state could also occur in systems with higher accretion rates than adopted here. The value $\dot{M}$ = 10⁻⁸ M_⊙ yr⁻¹ is a conservative choice. A larger $\dot{M}$ increases the viscous dissipation power, and therefore the energy available for EP acceleration. In this sense, the κ = 1 map can be viewed either as an upper bound or as a proxy for cases with κ < 1 but higher $\dot{M}$, which can produce EP-dominated regions of comparable extent and depth. The above results demonstrate that EPs, previously unexplored, may be a major ionisation source in thermochemical models of protoplanetary discs and could be relevant in planetary astrophysics and prebiotic chemistry.

Fig. 7 Ionisation rate due to EPs (ζEP) compared to the total ionisation rate ζtot = ζXR + ζGCR + ζEP, for the reference case κ = 0.1, where ζXR and ζGCR are the ionisation rates due to X-rays and GCRs, respectfully. The solid black contours mark regions where ionisation from EPs dominates over all other ionisation sources, while the dashed black line indicates where the EP ionisation is 5% of the total ionisation rate. The red contour outlines regions dominated by GCRs. The region in between is dominated by X-ray ionisation.

Fig. 8 Same as Fig. 7, but for κ = 1.

4.3.2 Energetic thresholds for the ionisation impact of reconnection-accelerated particles

To test the robustness of our conclusions, we varied the fraction of viscous dissipation power channelled into EP acceleration and examined the resulting ionisation rates. By reducing this fraction κ, we quantified how sensitive EP-driven ionisation is to the available energy budget and identified when EPs remain a significant source:

• When κ < 0.004, i.e., less than 0.4% of the local volumetric viscous dissipation is channelled into non-thermal energy, EPs do not dominate the ionisation in any region of the disc. Under these conditions, X-rays and GCRs remain the prevailing sources of ionisation throughout the disc;

• For κ < 2.5 × 10⁻⁴ (0.025% efficiency), EPs contribute to less than 5% of the total ionisation rate (ζ_EP < 0.05 ζ_tot) across the entire disc. In this case, their influence on the ionisation structure, and consequently on the chemistry, magnetic coupling, and observable signatures, is likely negligible.

These results show that ionisation by EPs from turbulent magnetic reconnection depends sensitively on κ. EP-dominated regions appear for efficiencies ≳0.4%, and even a very small conversion of 0.025% yields local contributions at the level of a few percent. Moreover, although subdominant, such localised enhancements are chemically and dynamically relevant. They can modify molecular abundances, shift the ionisation balance that controls magnetic coupling, and imprint distinct line-emission features. Crucially, this behaviour differs from X-rays and GCRs, whose uncertainties (even by factors of a few) tend to affect the entire disc more uniformly. The locality of EP-driven ionisation therefore offers an observational handle. Standard GCR and X-ray rates may reproduce molecular-ion emission across much of a disc, yet still struggle to model the inner disc ionisation. Indeed, forward modelling of DM Tau (Long et al. 2024) finds that models with X-rays and GCR ionisation matching the observations in the outer disc still under-produce inner-disc HCO⁺, H¹³CO⁺, and N₂H⁺ compared to observations below 50 AU, indicating an additional, more localised ionisation source in the inner regions. Accordingly, EPs accelerated by magnetic reconnection, even under conservative assumptions, should be included in thermochemical and dynamical disc models. Moreover, targeted inner-disc observations of molecular ions (e.g. HCO⁺, H¹³CO⁺, N₂H⁺) can help to discriminate EPs from globally acting X-rays and GCRs via their confined spatial signatures. We also note that this study neglects (i) re-acceleration by MRI-driven turbulence (Sun & Bai 2021) and (ii) other EP acceleration channels (Pinte et al. 2016), which could further strengthen or reshape these localised effects.

Fig. 9 EP volumetric heating rate map for κ = 1. The solid red contour encloses the region where EPs are accelerated. Dashed curves are isocontours of the EP volumetric heating rate. The upper white band marks very low columns (N < 1015 cm−2), i.e. above the disc atmosphere.

4.4 Ionisation volumetric heating rate

Ionisation not only sets chemistry but is also expected to heat the gas. When EPs ionise H or H₂, energy is deposited locally through elastic and inelastic collisions and subsequent chemistry. Following the standard prescription of Glassgold et al. (2012), we took per–ionisation heat deposition Q_H = 4.3 eV and Q_H2 = 18 eV, and computed the EP volumetric heating rate as $${\mathrm{\Gamma }}_{\mathrm{E}\mathrm{P}}(R,Z)=(n_{\mathrm{H}}{Q}_{\mathrm{H}}+2n_{{\mathrm{H}}_2}{Q}_{{\mathrm{H}}_2}){\zeta }_{\mathrm{E}\mathrm{P}},$$(22)

where n_H, n_H2 are the local atomic, molecular hydrogen number density, respectively and ζ_EP is the EP ionisation rate per H nucleus, computed in Eq. (21).

This microphysically-based estimation of the heating term can be directly incorporated into thermochemical models of protoplanetary discs and outflows. Indeed, previous studies have shown that heating at the disc surface significantly affects the launching conditions of photo-evaporative and MHD-driven winds (Casse & Ferreira 2000; Bai et al. 2016; Meskini et al. 2024). While these studies typically include stellar X-ray and UV heating as primary sources, we now argue that they should consider EPs.

Fig. 9 shows the spatial distribution of the EP volumetric heating rate, Γ_EP (eV cm⁻³ s⁻¹), assuming the extreme κ = 1 case. Heating is the strongest in the inner disc (within a few astronomical units) and at intermediate heights (Z/R ~ 0.15–0.3), where EP flux and gas density combine to maximise energy deposition, and it gradually declines toward larger radii and deeper layers. The solid red contour marks the zone where E_max > 10 MeV (efficient EP acceleration), dashed lines are isocontours, and the white region at the top indicates a very low column density (N < 10¹⁵ cm⁻²), above the disc atmosphere. These results indicate that EPs accelerated in the atmosphere channel accretion and magnetic energy into heating the gas deeper in the disc. We compared the volumetric heating from EPs, Γ_EP, to the total heating budget, Γ_tot, computed by PRODIMO (photoelectric, ${\mathrm{H}}_2^{\ast }$ de-excitation, H₂ formation and dissociation, gas–dust accommodation, X-ray Coulomb, viscous, chemical, etc.). Fig. 10 maps the ratio Γ_EP/Γ_tot throughout the disc. Even in the upper-limit case κ = 1, EP heating is negligible near the mid-plane and generally small elsewhere (mostly 10⁻³–10⁻²), becoming significant only in a narrow region at the disc surface in the inner AU and reaches close to unity values in the disc atmosphere beyond ~10 AU. For more realistic efficiencies, κ ~ 0.1, the EP contribution stays below ~1% across the disc, except for a small region at the surface around R ≳ 10 AU where it reaches the few percent level. These fractions are well within typical uncertainties of the total heating budget. Translated into temperature, the impact is expected to be minor compared with other heating sources. Consequently, EP heating appears secondary for thermochemical modelling, and its omission is unlikely to affect global temperature structures. By contrast, the EP ionisation term is the robust, model-relevant effect. It persists at large columns where X-rays are attenuated and can effectively modify chemistry. We therefore emphasise EP ionisation as the key addition, while treating EP heating as optional for extreme or surface-layer wind launching studies.

Fig. 10 Relative EP heating map for κ = 1. Ratio ΓEP/Γtot as a function of radius R (AU) and height Z/R; colours show the fraction on a logarithmic scale from 10−3 to 1. The red contour encloses the zone where EPs are efficiently accelerated (Emax > 10 MeV). The solid black contour marks where EP heating contributes at least to 5% of the total heating (ΓEP/Γtot = 0.05), which is considered to be a significant contribution. The dashed lines delimit the disc surface between column densities of 1019–1021 cm−2.

5 Discussion

In this section we summarise our findings and place them in context. We outline dynamical implications for accretion and wind launching, before closing with the main limitations of our framework and priorities for future work.

5.1 Turbulence-induced ionisation in discs: Interpretation, and thresholds

As we have shown, EPs accelerated in situ by turbulence-induced magnetic reconnection in disc atmospheric layers can substantially reshape the ionisation structure of protoplanetary discs. When the non-thermal power is a fraction κ ≳ 10⁻², EP ionisation overtakes stellar X-rays and GCRs across extended regions, inner radii and surface to intermediate layers, while even κ ~ 10⁻³ yields contributions of a few percent that are still chemically relevant. Because EPs penetrate to large column densities, they maintain significant ionisation where X-rays are attenuated, providing an additional heating channel at depth.

This behaviour follows naturally from expected atmospheric conditions: good magnetic coupling in far UV/X-ray–ionised layers, sustained MRI turbulence with v_turb ~ 0.1–0.3c_s and sub-H coherence, and fast turbulent reconnection with $V_r\sim V_A{M}_A^2$. In such environments, the turbulent cascade continually forms thin current sheets that release magnetic energy and inject supra-thermal particles. Modelling the acceleration with a Fermi-like prescription yields a steady non-thermal tail with p ≃ 2.5, bounded by cut-offs set by geometry and timescales, and normalised by the local energy budget via κ. This adds a disc-internal heating and ionisation channel that complements UV, X-rays, and GCRs, and remains effective at depths where external ionisers fade, with implications for both chemistry and non-ideal MHD.

Practically, for the T Tauri system considered here ($\dot{M}$ = 10⁻⁸ M_⊙yr⁻¹ and L_X = 10³⁰ erg s⁻¹) three regimes emerge:

• κ ≳ 10⁻²: The EPs dominated ionisation in the inner radii in the intermediate/deep layer and in the atmosphere;

• 10⁻³ ≲ κ ≲ 10⁻²: The EPs are subdominant yet chemically significant (at the level of a few percent) in extended regions and can still matter for specific chemical tracers and non-ideal MHD terms;

• κ ≲ 10⁻³: The EP effects on discs can be considered insignificant.

Here κ encapsulates the partition of released magnetic energy into heat versus particles, and electron–proton sharing, and is likely height- and radius-dependent. Here we set κ ≈ 0.1, by analogy with the canonical efficiency with which supernova-remnant shocks are thought to channel kinetic energy into GCR, in a future work, a disc-specific, first-principles estimate of κ will be developed.

The maps are most sensitive to three ingredients: (i) the available non-thermal power, κ; (ii) the ability to reach megaelectronvolt energies (controlled by α_acc, τ_esc, τ_inj, and geometry L_r/H, M_A); and (iii) local gas thermodynamics through the viscous parameter, α_eff, fixing the turbulent speed.

5.2 Implications for chemistry and observables

We expect reconnection-accelerated EPs to measurably alter the disc chemistry because their additional ionisation directly modifies the ion–molecule network and the electron fraction in the layers that they penetrate. As for X-rays and GCRs, EP ionisation of H₂ initiates the ${\mathrm{H}}_2^{+}\to {\mathrm{H}}_3^{+}$ pathways, enhancing the abundance of classical ionisation tracers (e.g. HCO⁺, N₂H⁺) and shifting their emitting layers. Cascades of secondary electrons further generate a CR-like internal UV field (Prasad–Tarafdar effect Prasad & Tarafdar 1983; Padovani et al. 2024), boosting photodissociation, ionisation and grain charging at high column densities. Because the ionisation rate drives the abundance of nitriles and hydrocarbons (e.g. HCN, C₂H₂), introducing ζ_EP(R, Z) where X-rays are attenuated will raise n_e, push ion–molecule chemistry to larger column density, and alter abundances and lines of these tracers, yielding observable signatures. To quantify these effects, the second paper in this series of articles, will post-process our disc backgrounds with ProdiMo, adding to the standard UV/X-ray/CR sources, the EP terms ζ_EP(R, Z) and Γ_EP(R, Z). We will explore a grid in non-thermal power fraction κ, spectral slope p, and acceleration geometry; solve to thermochemical equilibrium with and without ζ_EP; and map changes in gas temperature, electron fraction, and the layers where EPs overtake X-rays/GCRs. We will quantify EP-driven changes in classic ionisation tracers (HCO⁺, N₂H⁺) and, importantly, in mid-IR prebiotic and hydrocarbon species (e.g. HCN, HNC, C₂H₂) that originate in the warm surface layers where EPs effect is strong, exactly the regions probed by JWST/MIRI. We will generate synthetic MIRI diagnostics (e.g. band fluxes, ratios, emitting heights, line-to-continuum ratio) to assess how EPs imprint on those spectra and how to observe them. The outcome will be abundance and temperature maps and practical discriminants that isolate EP-driven ionisation from X-ray/CR scenarios and delineate the (k, p) domains where EP effects remain observable.

5.3 Dynamical implications for accretion and winds

Reconnection-accelerated EPs ionise and raise the electron fraction, n_e, in layers where stellar X-rays are already attenuated. In the simple gas-phase recombination balance, $n_e,\sim \sqrt{\zeta /\left({\beta }_{\text{rec}}{n}_{\mathrm{H}}\right)}$, where β_rec is the reconnection rate (Fromang et al. 2002), so boosting the ionisation rate by a factor, f, increases n_e by~ f^1/2. The EP contributions can therefore significantly increase n_e, with two immediate consequences.

First, on accretion and MRI coupling activation, larger n_e reduces Ohmic resistivity η_O ∝ 1/n_e and, by increasing the ion density, reduces ambipolar diffusivity, η_A ∝ 1/ρ_i. The corresponding Elsasser numbers, ${\mathrm{\Lambda }}_{\mathrm{O}}=v_A^2/\left({\eta }_{\mathrm{O}}\mathrm{\Omega }\right)$ and $Am=v_A^2/({\eta }_{\mathrm{A}}\mathrm{\Omega })$, rise. Layers that were marginal (almost active, Λ_O ~ 1, Am ~ 1) can cross MRI thresholds, shrinking dead zones and pushing the magnetically active region to larger column densities, which will increase the local accretion rate. This creates a bootstrap process: MRI produces turbulent reconnection that accelerates EP, which enhances ionisation triggering further MRI activation. In this picture, EPs produced in an MRI-active region may irradiate adjacent, marginally coupled layers, increasing the EP acceleration region size. A saturation arises from the β–dependence of EP acceleration. Turbulent-reconnection acceleration is most efficient in low-β regions; as the MRI-active layer expands and encroaches on higher-β gas, the reconnection-driven EP production drops. Once the front of the MRI-active zone reaches these high-β layers, it can no longer sustain a supra-thermal EP population able to ionise, so the outward growth of the acceleration region naturally stalls, quenching the bootstrap loop. A related pathway may operate in regions where the field is still sufficiently frozen-in to wind up. Reversals between toroidal field loops can generate thin current sheets even when MRI modes are suppressed. These sheets can undergo bursts of turbulent reconnection, producing rapid EP acceleration and temporary spikes in ionisation, potentially activating MRI and EP acceleration locally and intermittently.

Second, on wind launching and magneto-thermal coupling, magneto-thermal winds require gas that is well coupled to a reasonably coherent large-scale field at the launch base. Disc layers where winds are anchored are expected to be turbulent and only partially ionised, which lowers conductivities and suppresses efficient wind launching solutions. By adding ζ_EP(R, Z) where X-rays fade, EPs can thicken the well-coupled layer and, through their volumetric heating, Γ_EP, aid mass loading. This does not guarantee large-scale field coherence, but it may broaden the parameter space in which a steady magneto-thermal wind can be maintained, shifting the launch region and altering mass-loading and torque.

In a future work, we will post-process this additional EP source on a ProDiMo disc structure to clarify the importance and spatial extent of this bootstrap process and whether this mechanism can patch or reduce the dead zone. However, only dedicated non-ideal MHD simulations, including turbulent reconnection and EP transport self-consistently, can properly assess the impact of this mechanism. Including ζ_EP(R, Z) and Γ_EP(R, Z) in conductivity and energy modules will modify vertical profiles of η_O, η_A (and Hall where relevant), the MRI-active layer thickness, and the wind-launch region. Because EPs penetrate to a larger column density than X-rays, their impact is concentrated exactly where simulations currently struggle to maintain adequate coupling between the disc material and magnetic field. While EPs will not by themselves ‘solve’ the wind-launching problem, they provide an additional, physically motivated ionisation and heating channel that can tip marginal layers into the coupled regime, with measurable consequences for both accretion and outflows.

5.4 Limitations

Our framework entails several simplifying assumptions. First, on the acceleration model. We adopted a Fermi-like gain rate, $\dot{E}={\alpha }_{\text{acc}}E$, and connected α_acc, τ_esc, and τ_inj to a turbulent reconnection geometry. While these estimates are consistent with reconnection-in-turbulence theory and numerical trends, several efficiencies remain poorly constrained and depend on disc height and radius: (i) the partition of released magnetic energy into particles versus heat; (ii) the electron–proton sharing, for which future work will include injection and acceleration efficiencies, η_e, η_p; (iii) the fraction of turbulent power that feeds reconnection, set by the statistics and filling factor of current sheets; and (iv) the role of guide-field geometry and plasma, β_p, in regulating α_acc, τ_esc and pitch-angle scattering. To make these dependencies explicit and to avoid over-constraining uncertain microphysics, we normalised the non-thermal component by the local energy budget through κ. The phenomenological parameter κ thus absorbs uncertainties in the reconnection, particle and heat partition, and species sharing, allowing us to assess how the ionisation maps respond to realistic variations in the acceleration efficiency, while remaining grounded in the turbulent–reconnection geometry. A physically grounded estimation of κ will be provided in a future publication.

Second, on particle transport, we modelled EP propagation using the CSDA, neglecting pitch–angle scattering. CSDA is reliable up to column densities of order ~10²⁵ cm⁻²; beyond this depth (i.e. for protons above the pion–production threshold) a transport treatment including spatial diffusion and discrete losses is required (Padovani et al. 2018). In this set-up we also approximate the magnetic field as locally vertical, omitting its toroidal and poloidal components. A magnetic field configuration, combining toroidal shear and poloidal lines, would induce helical guiding-centre motion, thereby modifying path lengths and redistributing EP flux laterally. Magnetic field wandering would also lengthen particle path–lengths and generally reduce penetration depths (tending to shift EP–dominated regions upwards in Z/R), and may inject particles at higher radii. In addition, field–line channelling could focus fluxes in preferred directions; capturing these effects would demand coupling transport to self-consistent MHD fields. Furthermore, in turbulent flows, non-thermal particles can further gain energy by stochastic acceleration in randomly moving MRI-generated magnetic perturbations (Sun & Bai 2021). In this situation, non-thermal particles can also get (re)accelerated by Fermi second order process and shear acceleration. The particle maximum momentum, p_max, scales as the ratio of the local sound speed to light speed,c_s/c. Because of the small ratioc_s/c, it is clear that the acceleration effect is rather modest. Quantifying such re-acceleration is deferred to future work.

Third, on time dependence and intermittency, we computed the non-thermal distribution at the end of a current-sheet lifetime and then averaged over acceleration sites on a stationary disc background. The accretion rate used is likewise steady and conservative ($\dot{M}$ = 10⁻⁸ M_⊙ yr⁻¹), i.e. it does not include episodic enhancements. In reality, accretion, reconnection, and turbulence are intermittent, and young discs show non-stationary accretion with frequent bursts. During such bursts, a larger fraction of the (temporarily increased) accretion power can be channelled into EPs, raising the ionisation and heating rates (ζ_EP, Γ_EP), pushing EP influence deeper in the disc and driving time-variable chemical responses. Incorporating burst statistics and fully time-dependent transport and thermochemistry is a natural next step, beyond the steady-state framework adopted here.

As a perspective three next steps are appearing. (1) Thermo-chemical feedback and observables: post-process with PRODIMO including ζ_EP(R, Z) and Γ_EP(R, Z) to predict changes in abundances, temperatures, and synthetic observational ALMA/JWST diagnostics, and to identify discriminants between EP-, X-ray-, and CR-dominated regimes. (2) Calibration of efficiencies: constrain α_acc, τ_esc, τ_inj and particle and heat partition from stratified-box simulations that measure reconnection statistics, and confront these with observational diagnostics. (3) Transport: couple EP propagation to 3D MHD fields to include field-line diffusion, and stochastic re-acceleration, and extend beyond the CSDA limit.

6 Conclusion

We investigated the ionisation produced by EPs accelerated in situ by turbulence-induced magnetic reconnection in the atmospheres of protoplanetary discs. Building on a Fermi-like prescription for the energy gain and linking the acceleration, escape and injection timescales to a turbulent-reconnection geometry (Xu & Lazarian 2023), we derived a non-thermal particle distribution with a steady-state spectral index p = 2.5, an energy range bounded by a turbulence– and geometry–set cut-off, and a normalisation constrained by the fraction of accretion energy channelled to EPs through the parameter κ. Using a PRODIMO disc background, we computed EP propagation and the resulting ionisation and heating rate spatial distribution.

Our main conclusions are the following:

• A new leading ionisation source. For κ ≳ 10⁻², EPs dominate the ionisation budget over stellar X-rays and GCRs in the inner disc (sub-astronomical units to a few astronomical units) at normalised height Z/R ≈ 0.1 and in the disc atmosphere out to radii of a few tens of astronomical units (see Figs. 8, 7). Even for κ ~ 10⁻³, EPs contribute at the level of a few percent locally – enough to affect chemistry and magnetic coupling. Below 10⁻³, they are minor;

• Deep penetration. Compared to X-rays, EPs maintain stronger ionisation at large column densities, modifying the ionisation structure in layers that are otherwise shielded from radiation;

• Heating. EP-driven heating is typically modest compared to the total heating budget. For κ ≲ 0.1, Γ_EP stays below ~1% almost everywhere, reaching the percent level only in a narrow atmospheric band at R ≳ 10 AU. Ionisation not heating – is the primary impact of EPs in our framework;

• Observability and locality as a test. Because EPs are generated locally, they imprint spatially confined ionisation patterns distinct from globally distributed GCRs or stellar X-rays. This locality will provide testable signatures in emitting position of ion tracers (HCO⁺, N₂H⁺) and mid-IR species (HCN, C₂H₂) in the same layers that JWST (especially MIRI) probes. Comparing these diagnostics against our maps could offer a way to test our model.

Consequently, EPs from turbulent reconnection should be incorporated into thermochemical and dynamical disc models and for the interpretation of ionisation tracer observations.

In summary, turbulence-induced magnetic reconnection in disc atmospheres provides a physically motivated, in situ source of EPs that can reshape the ionisation structure of protoplanetary discs. Even when energetically subdominant, these particles penetrate and affect layers where other ionisers fade. When reconnection is efficient, they become a leading ionisation channel. Accounting for this component is therefore necessary for more realistic models of disc chemistry, dynamics, and their observational signatures. In turn, including this EP component in disc models will matter in the broader astrophysical context. In astrobiology, ionisation helps building (and destructing) complex organic molecules, sets where prebiotic species can form, and affects the interpretation of possible biosignatures in young planetary systems. More broadly, the same turbulence–reconnection framework can be adapted to accretion onto compact objects, where it may likewise contribute to the accretion–ejection dynamics.

Acknowledgements

V.B. and M.P. acknowledge the INAF grant 2023 MERCATOR (“MultiwavelEngth signatuRes of Cosmic rAys in sTar-fOrming Regions”); M.P. acknowledges the INAF grant 2024 ENERGIA (“ExploriNg low-Energy cosmic Rays throuGh theoretical InvestigAtions at INAF”). CHR acknowledge the support of the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) Research Unit “Transition discs” – 325594231. CHR is grateful for support from the Max Planck Society. C.S. aknowledge support for this work by the “Action Thématiques” “Physique Stellaire” (ATPS), “Phénomènes Extrêmes et Multi-messagers” (AT-PEM) and “Physique Chimie du Milieu Interstellaire (AT-PCMI) of CNRS/INSU PN Astro with INP and IN2P3, co-funded by CEA and CNES”. A.M. thanks the support of IN2P3 through the INTERCOS master project. The authors thank D. Galli, G. Lesur and D. Rodgers-Lee for fruitful discussions.

References

Appendix A Turbulent viscosity model

As already stated, in order to develop, MRI requires a sufficiently high degree of ionisation for the gas to interact with the magnetic field. The instability is characterised by the viscous parameter α_eff. To compute the viscous α_eff parameter, we rely on the model proposed by Thi et al. (2019), which introduces an empirical non-ideal MHD MRI-driven equation for α_eff suitable for physico-chemical protoplanetary disc codes such as Prodimo. The adopted expression for α_eff is given by $${\alpha }_{\mathrm{e}\mathrm{f}\mathrm{f}}={\left(\frac{2}{{\beta }_{\mathrm{m}\mathrm{a}\mathrm{g}}}\right)}^{1/2}min(1,{\mathrm{\Lambda }}_{\mathrm{O}\mathrm{h}\mathrm{m}}){[{\left(\frac{50}{{\mathrm{A}\mathrm{m}}^{1.2}}\right)}^2+{(\frac{8}{{\mathrm{A}\mathrm{m}}^{0.8}}+1)}^2]}^{-1/2},$$(A.1)

if $\sqrt{{\beta }_p}{\mathrm{\Lambda }}_{\text{Ohm}}>1$ and α_eff ≈ 0 otherwise.

Am represents the frequency at which neutral particles collides with ions normalised to the Keplerian orbital frequency, $$Am\equiv \frac{{\nu }_{\mathrm{i}\mathrm{n}}}{{\mathrm{\Omega }}_K}=\frac{{\beta }_{\mathrm{i}\mathrm{n}}{n}_{\mathrm{c}\mathrm{h}\mathrm{a}\mathrm{r}\mathrm{g}\mathrm{e}}}{{\mathrm{\Omega }}_K}=\frac{x_e{\beta }_{\mathrm{i}\mathrm{n}}{n}_{\mathrm{t}\mathrm{o}\mathrm{t}}}{{\mathrm{\Omega }}_K},$$(A.2)

where n_charge = x_En_tot is the total number density of charged species and β_in = 2 × 10⁻⁹ cm³ s⁻¹ is the collisional rate coefficient for ions to distribute their momentum to neutrals. Thus, $$Am\approx {10}^4\left(\frac{x_e}{{10}^{-4}}\right)\left(\frac{n_n}{{10}^{10}{\mathrm{c}\mathrm{m}}^{-3}}\right){\left(\frac{R}{1\mathrm{A}\mathrm{U}}\right)}^{3/2}.$$(A.3)

The Ohmic diffusivity is quantified by the dimensionless Elsasser number Λ_Ohm, defined as the ratio of the Lorentz force to the Coriolis force, $${\mathrm{\Lambda }}_{\mathrm{O}\mathrm{h}\mathrm{m}}\equiv \frac{B_z^2}{4\pi \rho {\eta }_O{\mathrm{\Omega }}_K}\equiv \frac{v_A^2}{{\eta }_O{\mathrm{\Omega }}_K}\equiv \left(\frac{4\pi {\sigma }_O}{{\mathrm{\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 ρ 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, σ_O ≈ σ_e,O and x_e > 10⁻¹³: $${\mathrm{\Lambda }}_{\mathrm{O}\mathrm{h}\mathrm{m}}\approx 1{\left(\frac{T}{100\mathrm{K}}\right)}^{1/2}\left(\frac{{10}^4}{{\beta }_p}\right){\left(\frac{R}{1\mathrm{A}\mathrm{U}}\right)}^{3/2}\left(\frac{x_e}{{10}^{-9}}\right),$$

The magnetic term β_p(R, Z) is the ratio of the thermal pressure to the magnetic pressure, given by (Thi et al. 2019), $${\beta }_p(R,Z)={\beta }_{\mathrm{m}\mathrm{i}\mathrm{d}}\frac{P_{\mathrm{t}\mathrm{h}}(R,Z)}{P_{\mathrm{t}\mathrm{h}}(R,0)}={\beta }_{\mathrm{m}\mathrm{i}\mathrm{d}}\frac{n(R,Z)T(R,Z)}{n(R,0)T(R,0)},$$(A.4)

where β_mid = β_p(R, 0) is the value of β_p on the mid-plane. β_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 β_mid is typically in the range 10⁴–10⁶. We take β_mid = 10⁴ as reference value.

Using the β_p and α_eff fields derived above, we map the turbulence/reconnection parameters across the entire PRODIMO grid. In particular, the local Alfvén speed is $$V_A=c_s\sqrt{\frac{2}{{\beta }_p}},$$(A.5)

and the Alfvénic Mach number is $$M_A\equiv \frac{v_{\text{turb}}}{V_A}$$(A.6)

wherec_s is the local sound speed and v_turb the turbulent velocity computed from Eq. 2.

All Figures

Fig. 1 Left: schematic radial cut from R ~ 0.1 to ~50 AU presenting our scenario qualitatively. The disc is vertically stratified into a neutral mid-plane (blue; MRI-stable), where EPs arriving from above travel quasi-ballistically (CSDA-like) along near-vertical paths, and an ionised atmosphere (yellow; low-βp MRI-active), where sub-Alfvénic turbulence (MA < 1) triggers turbulent magnetic reconnection that injects EPs, which then stream downwards (black arrows). In the bottom left part, a patchy variant of the ionised atmospheric layer is sketched, with a lower covering area and possible lateral growth (black arrows). In the ionised turbulent layer, EP transport could be treated as stochastic. The inner star–disc interaction zone treated in Brunn et al. (2023, 2024) is indicated at R ≲ 0.3. Right: radial cut from R ~ 0.07 to ~600 AU displaying the disc chemical and turbulent structure used in this work based on the PRODIMO model. The top panel shows the total hydrogen number density, nH,tot (colour), with vertical column-density contours N = 1015, 1019, 1021, 1024 cm−2 (dashed black). The dashed red line marks the ionised–atomic transition (H+/H) defining the border between the region where H+ dominates and the region where H dominates. Analogously, the dashed orange line delimits the atomic–molecular transition (H/H2). The bottom panel displays the Alfvénic Mach number MA ≡ vturb/VA (colour). The red curve traces MA = 1. The sub-Alfvénic region (MA < 1) identifies where fast, turbulence-enabled reconnection, and thus EP acceleration, is expected. The solid red lines delimit the region where EP are efficiently accelerated, i.e. to energy more than 10 MeV. The white region is the so-called ‘dead zone’, where MRI is inefficient (see criterion Eq. (A.1)), there is no turbulence, MA is set to 0.

In the text

Fig. 2 Sketch of MRI-driven turbulence and embedded turbulent magnetic reconnection. Left: eddies at the injection scale L in a sub-Alfvénic regime characterised by the Alfvénic Mach number MA. Right: reconnection layer produced by the turbulence, with length Lr = L and thickness Δr. The turbulence parameters {L, MA} set the reconnection, geometry parameters {Lr, Δr}, speed $V_{\text{rec }}\sim V_A{M}_A^2$, and EP acceleration parameters used in this work.

In the text

	Fig. 3 Nonthermal–to–thermal energy density ratio, Unt/Uth, for κ = 0.1. The red contour encloses regions where the local maximum particle energy exceeds 10 MeV (efficient acceleration sites capable of supplying ionising EPs). The black contour marks a vertical column density, N = 1019 cm−2 (approximate disc surface).
In the text

	Fig. 4 Particle energy distributions at R = 1 AU for κ = 0.1. The dashed curves show the thermal Maxwellian based on PRODIMO density and temperature for Z/R = 0.25–0.35. The solid curves show the non-thermal distribution produced by magnetic-reconnection acceleration.
In the text

Fig. 5 Ionisation rates as function of the disc vertical column density at R = 1 AU for κ = 0.1. The lines in the viridis colour scale are the contribution to the ionisation rate coming from EPs accelerated in region at different Z/R. The thick orange line is the total ionisation rate corresponding to the weighted sum (Eq. (21)) of each local contribution. The black line is the ionisation rate from the stellar X-rays.

In the text

	Fig. 6 EP ionisation rate map for κ = 0.1. The solid red contour outlines the region where EPs are accelerated. Dashed lines indicate isocontours of the EP ionisation rate. The upper white region corresponds to column densities below N = 1015 cm−2, indicating regions above the disc atmosphere.
In the text

Fig. 7 Ionisation rate due to EPs (ζEP) compared to the total ionisation rate ζtot = ζXR + ζGCR + ζEP, for the reference case κ = 0.1, where ζXR and ζGCR are the ionisation rates due to X-rays and GCRs, respectfully. The solid black contours mark regions where ionisation from EPs dominates over all other ionisation sources, while the dashed black line indicates where the EP ionisation is 5% of the total ionisation rate. The red contour outlines regions dominated by GCRs. The region in between is dominated by X-ray ionisation.

In the text

	Fig. 8 Same as Fig. 7, but for κ = 1.
In the text

	Fig. 9 EP volumetric heating rate map for κ = 1. The solid red contour encloses the region where EPs are accelerated. Dashed curves are isocontours of the EP volumetric heating rate. The upper white band marks very low columns (N < 1015 cm−2), i.e. above the disc atmosphere.
In the text

Fig. 10 Relative EP heating map for κ = 1. Ratio ΓEP/Γtot as a function of radius R (AU) and height Z/R; colours show the fraction on a logarithmic scale from 10−3 to 1. The red contour encloses the zone where EPs are efficiently accelerated (Emax > 10 MeV). The solid black contour marks where EP heating contributes at least to 5% of the total heating (ΓEP/Γtot = 0.05), which is considered to be a significant contribution. The dashed lines delimit the disc surface between column densities of 1019–1021 cm−2.

In the text