Velocity-space distribution function of fast ions in a sawtoothing plasma

This study explores the influence of sawtooth oscillations on the velocity space distribution of fast ions in tokamak plasma discharges. The relevant Fokker–Planck equation for fast ions is solved analytically. Two distinct effects arising from the temperature drop associated with a sawtooth crash and their impact on the distribution function of fast ions are considered. The first effect involves the modulation of the fusion alpha particle source on the timescale of the sawtooth period, linked to the drop in fusion yield resulting from the sawtooth temperature relaxations. The second effect is tied to the increase of the slowing-down time during the sawtooth ramp, causing particles born later in the sawtooth cycle to experience reduced slowing down compared to those born right after the crash, creating an accumulation-like mechanism at higher energies. In regimes where the sawtooth period is shorter than the fast ion slowing-down time, the combined influence of these effects gives rise to fast ion distribution functions that transiently exhibit positive slopes in velocity space.


Introduction
The next generation of magnetic confinement fusion devices will focus its effort towards the production of a burning plasma, in which strong self-heating is provided by the alpha particles generated with an energy of 3.5 Mev by fusion reactions in D-T plasmas.The presence of alpha particles and other energetic ions, produced via auxiliary heating, strongly influences the magnetohydrodynamic (MHD) stability of the plasma discharge.Wave-particle resonant interactions can induce various MHD instabilities, ranging from energetic particle modes (EPMs) to the entire class of Alfvén Eigenmodes [1].Conversely, fast ion populations may also be strongly affected by the plasma MHD activity.
An important example are sawtooth oscillations [2], a well known type of MHD activity leading to periodic relaxations of the pressure, density and q-profile in the core region of the plasma.On one hand, sawtooth oscillations can be stabilized by the presence of fast ions [3], leading to so-called monster sawteeth, characterized by long sawtooth periods and large crashes [4].On the other hand, the loss of core plasma confinement and the redistribution of the pressure towards more peripheral plasma regions affects the fast ions.Therefore the radial redistribution of energetic ions due to sawtooth crashes has been the subject of several studies [5][6][7].It has been shown that fast ions with MeV energies may be very weakly redistributed in space during sawtooth crashes as discussed by Kolesnichenko et al [8][9][10] and recently confirmed by kinetic-MHD simulations by Bierwage et al [11].
Here, we will assume that fusion alpha particles with MeV energies will not undergo any radial redistribution during the sawtooth dynamic.This work focus on the changes of the distribution function in velocity-space due to sawtooth oscillation for the energetic tail of isotropic fast particle distributions.Indeed, an important element in characterizing the waveparticle interaction is the velocity-space gradient of the fast ion distribution function.This becomes even more crucial for modes characterized by toroidal mode number n = 0, like Energetic particle induced geodesic acoustic modes [12] and vertical displacement oscillatory modes destabilized by fast ions [13], for which the sign of the velocity-space gradient may determine whether wave-particle resonance leads to the mode destabilization or damping.
Two effects associated with the core temperature drop caused by sawtooth crashes can modify the alpha population velocity distribution function: the first is the drop of fusion yield after the core temperature crash, effectively modulating the source of fast ions on the sawtooth period, (τ saw ), timescale.Secondly, the slowing-down time of fast ions, (τ sd ), grows during the sawtooth ramp, inducing an accumulation at higher velocities of the particles generated in the later stages of the ramp.Considering these effects in a simplified Fokker-Plank equation, it is shown that, instead of the steady-state slowing down distribution function, characterized by the 1/v 3 behavior, a time dependent distribution function, periodic with period τ saw , can be obtained.Depending on the ratio between the slowing-down time and the sawtooth period, τ s /τ saw , the derivative over velocity of the newly obtained distribution function can change sign.
The paper is organized as follows: in section 2 the plasma pressure evolution is described with a fully analytical model considering triangular sawtooth oscillations and taking into account the Kadomtsev reconnection model [14].The consequences of the periodic oscillations of the plasma pressure on the alpha source and slowing-down time are also discussed in section 2 .In section 3, the standard simplified Fokker-Plank equation considered is reported, as well as the main assumptions required in order to proceed analytically in our analysis.This section is then divided in four subsections: section 3.1 briefly describes the standard slowing down distribution.In sections 3.2 and 3.3 the distribution functions obtained considering the source modulation and the slowingdown time evolution, are discussed independently.Finally, section 3.4 describes the distribution function obtained considering both effects together.Results and concluding remarks follow in section 4. Appendix A shows the analytic procedure followed to describe the plasma pressure evolution during sawtooth oscillations, while in appendix B the full derivation of the complete time-dependent distribution function is reported.

Sawtooth-relaxed core pressure profile
Sawtooth oscillations [2] are relaxation oscillations of the plasma profiles in tokamak plasma discharges in the core region near the magnetic axis where the safety factor, q, Figure 1.Precrash (Blue) and relaxed (Orange) q profiles up to the mixing radius r mix = √ 2rs.The q = 1 surface of the precrash profile is at rs/a = 0.3 (dashed vertical line) and its value on axis is q 0 = 0.9.drops below unity.Two main stages of this quasi-periodic phenomenon can be identified: a slow temperature rise, determined by heating and transport, called the sawtooth ramp, and a much faster phase, called the sawtooth crash, triggered by the onset of a kink instability dominated by toroidal and poloidal mode numbers n = 1 and m = 1.This second phase describes the loss of core confinement on a timescale much faster than the transport one, and the consequent flattening of the profiles up to a radius called mixing radius r mix , which is larger than the q = 1 resonant surface radius, r s .Different theoretical models have been proposed for the determination of the relaxed profiles, i.e. the profiles after the crash [14,15].Here we consider the well known Kadomtsev relaxation model [14], that, under certain conditions, allows for analytic work.This relaxation model prescribes that the relaxed q profile is a monotonically increasing function of radius, with the relaxed on-axis value q 0,rel = 1.Assuming the following analytic profile for the inverse of q before the crash, the relaxed q and pressure profiles up to the mixing radius r mix can be obtained analytically where ∆q = 1 − q 0 and q(r s ) = 1.Following the Kadomtsev relaxation model, as explicitly described in appendix A, the resulting mixing radius is r mix = √ 2r s , and the relaxed inverse q profile reads: Figure 1 shows an example of precrash and relaxed profiles up to the mixing radius as a function of normalized radius r/a, a being the plasma minor radius.
The relaxed pressure can be determined through energy conservation arguments.The pre-crash pressure profile in the plasma core is assumed to be well approximated by a parabolic pressure: with p 0 on-axis value and r p > r mix (notice that, with r p = r mix , the pressure would go to 0 at r = r mix ).For the sake of simplicity, a constant density profile not affected by the sawtooth dynamics is assumed.Due to this profile choice, as detailed in appendix A, the relaxed profile is uniform up to the mixing radius.The sawtooth ramp phase is assumed to be linear with time, corresponding to so-called triangular sawteeth.The temperature evolution within the mixing radius during the first cycle (0 ⩽ t ⩽ τ saw ) then follows, for r < r mix : where As shown in figure 2, due to the choice of q and pre-crash pressure profiles of equations ( 1) and (3), the relaxed temperature profile flattens to a value T 0,rel up to the mixing radius.Just before the next crash, i.e. at t = τ saw , the central temperature peaks at ] .In the region beyond the mixing radius, the profile is assumed to be unchanged with The Kadomtsev relaxation model introduces discontinuities in the q and T profiles at r = r mix after the crash.These discontinuities decay rapidly during the early part of the subsequent sawtooth ramp as a consequence of current and temperature diffusion.
For standard sawtooth oscillations, the crash time τ crash is much smaller than the sawtooth period.In the following the crash is assumed to be instantaneous (τ crash → 0) with respect to the timescales of interest.This allows the temperature evolution during the sawtooth dynamic to be described with: where ⌊t/τ saw ⌋ is the integer part of the quantity t/τ saw .It is useful to define a piece-wise temperature evolution during each sawtoooth cycle: where H(x) is the Heaviside step function; the integer n labels the sawtooth cycles and is associated to the integer part of t/τ saw.An alternative expression for the temperature dynamics during multiple sawtooth cycles is: where n max labels the last sawtooth cycle.Notice that the expressions in equations ( 5) and ( 7) are completely equivalent for positive t, up to t = n max τ saw .
A time-dependent, monochromatic source of alpha particles with birth-velocity v α is assumed: where δ(x) is a Dirac delta function, and the parameter K is related to the total number of alphas born per m 3 per second.For a D − T plasma with n D = n T = n e /2, K can be obtained integrating the source of equation ( 8) over the velocity: The reaction rate σv for D − T reactions averaged over Maxwellian distribution can be expressed as a function of T only [16] for the typical temperature range of tokamak experiments: with σv| DT in m 3 s −1 and T in keV.Taking into account the temperature evolution of equation ( 4), the source of alpha particles will be time dependent, and periodic with period τ saw : where Alpha particles generated by D − T fusion reactions will exchange their energy with the thermal plasma via collisions, leading to their slowing down.For energies larger than a critical value E c ≈ 41T e , collisions with electrons are dominant.
In the regime characterized by E > E c , which is the one relevant fusion born alphas, the slowing-down time required by the alphas to reach the critical velocity via electron collisions is: with E α = 3.5 MeV is birth energy of the alpha particles and the Spitzer time defined as: with c τ = 1.77 • 10 −2 , T in keV, n e in 10 20 m −3 , and τ se in s [16].We neglect the weak temperature dependence of E c in the logarithm of equation (12).Following equation ( 4), also the slowing-down time will depend on time due to the sawtooth dynamics.
Considering a density n e = 1 • 10 20 m −3 , and a value of the relaxed temperature of 8 keV, the on-axis peak temperature is T 0 ≈ 10.7 keV (for r p = 0.6, r s = 0.3).The corresponding slowing-down time will be varying between 0.33 < τ sd [s] < 0.50.

Fokker Plank equation and resulting distribution functions
We restrict the analysis to energetic particles with E > E c .In this regime, the alpha particles collide more importantly with the thermal electrons and the collision operator in the Fokker-Plank equation can be simplified.During the sawtooth ramp, the relevant equation governing the time evolution of the distribution function f α can expressed as [17,18]: where S α is an isotropic source term, and a loss term, absent in references [17,18], is included, with ν α a constant loss frequency.The introduction of a constant loss term is considered for simplicity.While this is a very crude approximation, used for instance in reference [19], it serves as a first step in describing charge exchange losses.The radial derivative in the Fokker-Planck equation has been neglected, following the approach adopted e.g. in [17,18].This choice is justified by two main reasons: firstly, our focus lies on the process of fast ions slowing down during the sawtooth ramp, assuming no radial redistribution of fast ions due to the sawtooth crash.Secondly, in the regime where orbits are thin compared to the q = 1 radius, the slowing down in velocity space is mainly a local phenomenon.However, it is worth noting that for orbit widths that are not negligibly small, the radial derivative becomes relevant, as discussed in, for instance, reference [20].In equation ( 14), we neglect also corrections of order (v c /v) 3 , where v c is the critical velocity.
Defining F α = v 3 f α , equation ( 14) can be solved for F α using standard methods, as detailed in appendix B. In the following sections, the solution of equation ( 14) is discussed under different assumptions for the source term and the slowing-down time.

Standard slowing-down distribution function and the effect of particle losses
Firstly, the solution of the Fokker-Plank equation is recalled for the standard case with constant temperature and constant loss frequency [18,21].Under these assumptions, the slowingdown time is constant, τ sd0 , while the source term is simply The source is considered to be switched on at t = 0.The time-dependent solution for the distribution function under these assumptions is: where l α = ν α τ sd0 /3.Equation ( 15) time dependence comes only from the last Heaviside function.The steady state solution is This solution is illustrated in figure 3.
When the loss term exceeds the value ν α = 3/τ sd0 , i.e. l α > 1, most particles are lost before their complete slow down, and the distribution function develops a positive slope as function of time, ∂f α /∂v > 0, as shown in figure 3(b).

Source temperature modulation at constant slowing-down time
In this section, we consider the solution of equation (14) where, for the sake of simplicity, the slowing-down time is kept constant (to an average value τ sd0 ) but the source term varies in time following the temperature evolution of equation (7).Thanks to the assumption of instantaneous crashes, the solution can be studied during a single cycle, i.e. for a fixed n.The time dependent source term during a cycle can be written as: where S rel = c S n 2 e T 2 0,rel /v 2 α is the source term associated with the relaxed (flat profile) temperature.As detailed in appendix B, the time dependent solution of equation ( 14) for particles born during the nth sawtooth cycle is: Figure 4 shows the solution of equation (19) as a function of time, assuming that the source term is turned on for a single sawtooth cycle.In other words, this figure shows the evolution in velocity space of a group of alpha particles created during a single sawtooth cycle as they slow-down in time.The resulting slope is reduced compared to the standard v −3 behavior, since the source is increasing with time during the sawtooth cycle.
The full solution for the relevant case where the source term is on over several cycles can be obtained by summing all the contributions coming from each n: After one slowing-down time, i.e. the time necessary for the particles born during the first (n = 0) cycle to reach the critical velocity, the solution of equation ( 20) becomes periodic, with period τ saw .In figure 5, a snapshot at t = t * = 5.5τ saw of the solution with τ sd0 /τ saw = 3 is plotted.At this time, all the particles born during the first sawtooth period (n = 0) already reached v c , and only 4 steps are visible associated with Thus, the source modulation has an important effect on the slope of the distribution function in the velocity space when the slowing-down time is longer than the sawtooth period.If the ratio τ sd,0 < τ saw is of order of unity or smaller, this effect becomes negligible.
The situation described in this section can be considered as representative also for the case of modulated NBI sources with a constant background plasma temperature.In this case, the solution of equation ( 19) describes a pitch-angle averaged distribution function for NBI ions injected following a sawtoothlike waveform signal, in contrast with the pulses considered in [22].The effect on the fast ion population of a modulation of the NBI on timescales comparable with the slowing down  time has been successfully proven experimentally on DIII-D in [23].

Time-dependent slowing-down time and constant source
In this sub-section, the effect of a non-constant slowingdown time is analyzed next.In this sub-section, the source is assumed time independent, S α = S 0 − v α ), while the slowing-down time depends on time through temperature: where τ sd,rel is the slowing-down time in the relaxed, postcrash phase, where the on-axis temperature is lowest and the slowing-down time is shortest.Also in this case we can follows fast particles generated during different sawtooth cycles and eventually combine their contributions.However, due to the discontinuities in τ sd associated with the temperature variation at the crashes, particular care is required, as discussed extensively in appendix B. In order to simplify the discussion, losses are neglect in this sub-section, ν α = 0.The distribution function for particles born during the nth sawtooth cycle is: Here, tn (r, t, v) is a rather complicated function proportional to τ saw , defined in equation (B.36) of appendix B. As shown for the source modulation case in figure 4, particles born in each sawtooth cycle are associated with time-dependent minimum and maximum velocities.The same is true for varying τ sd , and these velocities are vn,min and vn,max , defined in equations (B.37) and (B.38).Similarly to the case with constant τ sd , the full solution can be obtained summing the contributions of each cycle: A snapshot of the full solution at time t = t * = 5.5τ saw is plotted in figure 7 for τ sd,rel /τ saw = 3.Also in this situation, the distribution function becomes periodic, with period τ saw , after the particles born during the first (n = 0) cycle have reached the critical velocity.
During each sawtooth ramp, inside the q = 1 radius the temperature linearly increases.Particles born early-on experience stronger slowing down due to electron collisions with the colder plasma, with respect to those born later that are less affected by collisions.An accumulation at high energy of particles born in the later stages of the ramp is thus obtained.Depending on the parameter τ sd,rel /τ saw , this effect can lead to a sign change of the derivative ∂f α (r, t, v)/∂v.However, similarly to the source modulation case, this effect alone is not able to invert the slope of the distribution function for values τ sd,rel /τ saw < 10.
For comparison, the standard slowing-down solution for a constant slowing-down time corresponding to its value in the relaxed state is shown by the dashed curve in figure 7. Note that the two solutions (blue and dashed line) coincide at velocities related to relaxation times, and that the blue curve corresponds to a larger population of fast particles, since with increasing slowing-down time the fast particle take longer to reach the sink at v = v c .
The analogy with NBI sources introduced at the end of section 3.2, can be extended to this section considering the solution of equation ( 22) as a pitch-angle averaged distribution.It is then appropriate to describe the fast ion population injected with a constant source into a sawtoothing plasma.

Time-dependent τ sd and time-dependent source
The two effects discussed in sections 3.2 and 3.3 are considered together in this section, which discusses the relevant case for the alpha particles generated by D-T fusion reactions in a sawtoothing plasma.The effect of the constant loss term is neglected, i.e. we assume ν α = 0.The source of particles follows equation (18), and the slowing-down time varies according to equation (21).As in section 3.2, we first focus on the single nth cycle solution.Following the procedure detailed in appendix B, the resulting distribution function is: The full solution is obtained again as the combination of the solutions associated to different sawtooth cycles, each numbered by ns.Summing over all cycles, we obtain: Figure 8 shows how the combination of the two effects enhances the slope change of the distribution function in velocity space, arising due to the periodic temperature evolution.
A situation with piece-wise positive slope can be achieved with values of τ sd,rel /τ saw ∼ 3, relevant for present day tokamak experiments.In the limit of τ saw → ∞, corresponding to a constant temperature, the standard slowing-down solution of equation ( 15), corresponding to the black dashed curve in figure 8, is recovered.

Discussion and conclusions
In this work, distribution functions describing the velocity space population of fast ions in sawtoothing plasmas have been derived analytically.The case of particles generated by an isotropic source in the pitch angle variable, appropriate for alpha particles generated by fusion reactions, has been considered.We have restricted our derivation to particles with birth velocity much larger than the critical velocity, v α ≫ v c , allowing us to neglect corrections of the order of (v c /v) 3 and to consider the simplified form of the Fokker-Plank equation given by equation (14).Furthermore, we have assumed that the particle energy is sufficiently high for the fast ion radial profile to remain peaked after a sawtooth crash, allowing us to neglect fast ion redistribution in space.The validity of this assumption relies on the conservation of the third adiabatic invariant of fast ion motion during the crash.The condition for the conservation of the third adiabatic invariant is described in terms of the precession frequency of the fast ions compared to the sawtooth crash duration.[8,9].For alpha particles with MeV energies, this condition is satisfied, as discussed in [11].Thus, we have assumed in this article that the radial profile of the fast ions is unaffected during a sawtooth crash.It is important to note that, depending on plasma parameters, the condition for the conservation of the third adiabatic invariant may be satisfied by fast ions with magnetically trapped (i.e.banana) orbits, but not by fast ions with circulating orbits, potentially introducing velocity-space anisotropy not considered in this work.The temperature evolution during the sawtooth dynamics has been described in accordance with the Kadomtsev reconnection model, while sawtooth crashes have been assumed instantaneous, being their characteristic timescale much faster than the ones of interest in our work.
Time-dependent distribution functions can be derived analytically by solving the simplified Fokker-Planck equation.
Following an initial transient, corresponding to a slowingdown time, these distribution functions become periodic on the sawtooth period timescale.
Two main effects determining the slope of the distribution function in velocity have been identified.Firstly, sawtooth oscillations induce a modulation of the alpha particle source on the sawtooth period, τ saw .Secondly, the periodic relaxation of the temperature profile gives rise to a time-dependent slowing-down time.As a consequence, an accumulation-like mechanism of energetic particles born in the later stages of a sawtooth ramp is obtained.
For both effects, the ratio τ sd,rel /τ saw is the main parameter related to the changes of the distribution function in velocity space.Here, τ sd,rel represents the slowing-down time associated with the plasma temperature on axis in the relaxed state.In a sawtooth cycle, the lowest temperature on axis occurs just after the crash, and therefore τ sd,rel reaches its minimum value in the relaxed state.When the sawtooth period is short compared with the slowing-down time, both effects may lead to a distribution function with positive ∂f α (r, t, v)/∂v.
If the two effects are considered separately, as in sections 3.2 and 3.3, distribution functions with regions where ∂f α (r, t, v)/∂v > 0 are only obtained for values of τ sd,rel /τ saw ∼ 10.On the other hand, the combination of source modulation and accumulation-like mechanisms, discussed in section 3.4, results in a distribution function with piecewise regions of positive slope in velocity space for values of τ sd,rel /τ saw ≈ 3.In present-day tokamaks like JET, scenarios where τ sd,rel /τ saw > 1 can arise.Interestingly for our work, this condition has been observed alongside the emergence of fastion driven n = 0 modes, suggesting the existence of a distribution function in velocity space with ∂f α /∂v > 0, [24].
The parameter τ sd,rel /τ saw plays an important role in our analysis.The temperature recovery time after a sawtooth crash is of the same order as the sawtooth period, τ saw , if standard sawteeth with triangular traces are considered.However, in cases involving different types of sawteeth, such as the so-called 'monster sawteeth' [4], characterized by robust fast ion stabilization that prevents crashes, the relevant time scale, to be compared with the fast ion slowing-down time, is the temperature recovery time immediately following a sawtooth crash, and not the duration of the monster sawtooth.Our analysis suggests that fast ion distribution functions with a positive slope may appear transiently, following a monster sawtooth crash, if the temperature recovery time is shorter than the fast ion slowing-down time.
In addition to the two effects associated with temperature oscillations during the sawtooth dynamics, particle losses also impact the velocity space distribution function.In sections 3.1 and 3.2, we discussed how particle losses on a timescale comparable to the slowing-down time can facilitate the sign change of the distribution function slope.
The analytic distribution functions obtained in this work focused on the case of alpha particles generated by fusion reactions.However, the source modulation effect described in section 3.2 can be obtained for energetic ions injected with NBI modulated in time, as discussed theoretically in [22] and experimentally proven in [23].In this context, the distribution function obtained in section 3.2 qualitatively represents a pitch-angle averaged distribution function of NBI fast ions when their source follows a sawtooth-like waveform signal.The second effect, related to the time-dependent slowingdown time, is not limited to the case of alpha particles generated by fusion reaction, but affects fast ions injected with any source in a sawtoothing plasma.
The present work represents a proof-of-principle that positive slopes in the distribution function of fast ions may result because of sawtooth activity.Clearly, sawtooth oscillations can appear in different forms, going from the more standard, triangular sawteeth [2], to monster sawteeth [4], compound sawteeth [25], humpback sawtooth oscillations [26], and many other types.Also, while Kadomtsev model for sawtooth relaxations appears to be a fair approximation for smaller and more resistive tokamaks, large-size tokamaks such as JET have entered regimes where semi-collisional effects do play an important role, giving rise to incomplete magnetic reconnection (see, e.g.[15,27], and other references quoted therein).Each different type of sawtooth phenomenology would require specific considerations with regards to the magnitude of the effect discussed in the present work leading to a positive slope in the fast ion distribution function.Nevertheless, the main mechanism is present in different types of sawtooth activity, in that they all have in common that temperature profiles flatten rapidly during a sawtooth crash and recover slowly after the sawtooth crash.For example, temperature flattening occurs for both complete and incomplete sawtooth reconnection models, as discussed in [15,27].The analysis presented in this article suggests that, for any type of sawtooth activity, if the temperature recovery time after a sawtooth crash is fast as compared with the fast ion slowing-down time, then a positive slope in the fast ion distribution function may be produced, at least transiently.For the more standard, triangular sawteeth, the temperature recovery time is of the order of the sawtooth period, which is the specific case considered in this article.
The regions of positive slope in velocity space described in this work may be the cause of transient destabilization of MHD modes with toroidal mode number n = 0, as observed in recent JET experiments [24,28].Thus, the distribution function obtained analytically here may be used as a model for the analytical assessment of the kinetic-MHD stability of EPMs with toroidal mode number n = 0 in plasmas characterized by fast sawtooth activity.
In order to determine the relaxed pressure profile, we consider that, during the reconnection, the thermal energy of the plasma involved in the process is conserved, i.e. we impose the conservation of the thermal energy integral within the cross-sectional area between reconnecting surfaces.After simple manipulations, we can write and, differentiating with respect to r 2 k : If we assume that the pre-crash pressure profile inside the mixing radius r mix = √ 2r s is parabolic, where p 0 is the onaxis pressure, and r p > r mix , then the relaxed pressure profile is constant: Note that a different assumption for the pre-crash pressure profile in general leads to a non-constant relaxed pressure.The pressure profile is assumed to depend on time during sawtooth ramps according to: and of the parameter β(t) = 1 − p rel /p 1 (t).The time dependent parameter p 1 (t) can be determined describing the temperature evolution during the ramp.Here, we assume that the onaxis temperature increases linearly in time from p rel to p pre (0).In this way one obtains: With β(t) = 1 − p rel /p 1 (t), the pressure profile evolution reads:

Appendix B. Solution of the Fokker Planck equation for fast ions
In this appendix we detail the solution of equation ( 14), under the different assumptions considered in the main body of this article.Defining F α = v 3 f α , one must solve It is possible to restrict the analysis to the solution of the differential equation for particles born during the generic nth sawtooth cycle (with nτ saw < t < (n + 1)τ saw ).Considering a constant slowing-down time and a varying source term, as in section 3.2, we can solve equation (B.2) using the method of characteristics.The characteristics equations are: v = c 1 e −t/τ sd0 (B.3) These are used to solve the differential equation: To solve equation (B.5) we consider: and thus the integral equation: The time-dependent source term can be expressed as with S rel = c S n 2 e T 2 0,rel /v 2 α being the source term associated with the relaxed temperature after a sawtooth crash.Applying the characteristic relation of (B.3) in equation(B.8)one obtains: Integrating along the characteristics and exploiting the properties of the delta function in equation (B.9), now written in terms of c 1 and time as δ(c 1 e −t/τ sd0 − v α ), C reads: where t 0 = τ sd0 log (c 1 /v α ).Finally, considering the second characteristic relation (B.4), the following expression for C(t, r, v) is obtained: Lastly, from expression (B.6) for F α,n , and manipulating the argumet of the Heaviside function, one obtains the relevant distribution function describing particles born during the nth sawtooth cycle: where l α = ν α τ sd0 /3.The same technique can be considered in order to solve equation (B.2) with constant source term and varying slowingdown time, as in section 3.3.However, more care is required due to the discontinuities in τ sd introduced by the sawtooth crashes.For simplicity, we neglect the loss term, i.e. we set ν α = 0. Due to the discontinuities in τ sd , we consider again particles born during the nth sawtooth cycle, but we solve the differential equation (B.2) separately for each successive sawtooth cycle.During the nth sawtooth cycle, the slowing-down time varies according to: We can proceed using the method of characteristics, with characteristic equations: Integrating in time, Using equations (B.16) and (B.17), one obtains the solution valid for the nth sawtooth period: where We can rewrite the Heaviside function in terms of velocities: with: We now look for the solution associated to particles born during the nth cycle during the (n+1)th sawtooth cycle.In this case the source term is zero and we have to solve the following differential equation: It is still possible to proceed considering the characteristic equation, but during this cycle the slowing-down time follows τ sd (r, t) = τ sd,rel h n+1 (r, t)  The same procedure can be followed in order to derive the solution for later sawtooth periods (n + 2,n + 3,…).A recursive expression for the distribution function can be obtained for the evolution of the particles born during the nth cycle for all times: Remembering the definitions of F α,n and h n (t, r), the distribution function reads: × τ sd,rel In the situation with both the source term and slowingdown time depend on time through the plasma temperature, as in section 3.4, the analytic procedure is again the one just described for the case of non-constant slowing-down time of section 3.3.Since the source term and the slowing-down time vary in the same way, following the temperature evolution, the only additional effect of the varying source with respect to the previous case lies in the exponent of the varying temperature function, h n (r, t).The full distribution function can then be expressed as:

Figure 2 .
Figure 2. Temperature profile normalized with the relaxed on-axis value T 0,rel , at different times during the sawtooth period.The q = 1 radius at rs/a = 0.3 and the mixing radius at r mix /a = 0.3 √ 2 are shown with dashed vertical lines.The pressure scale length is rp/a = 0.6, while outside the mixing radius the profile is time-independent.The relaxed temperature profile exhibits a jump at r = r mix that is rapidly smoothed out by temperature diffusion in the early part of the following sawtooth cycle.

Figure 3 .
Figure 3. Steady state distribution functions as function of normalized velocity, for constant τ sd and source term, plotted for v ⩾ vc (the latter velocity indicated by the vertical dotted line).Panel (a): the standard slowing down distribution with no losses (lα = 0).Panel (b): the slowing-down solution for lα = 2.

Figure 4 .
Figure 4. Distribution function for particles born during the nth cycle with time dependent source and constant τ sd , plotted for different times as a function of normalized velocity for r = 0.The critical velocity vc is showed with vertical dotted line.The minimum and maximum velocities reached by the nth cycle particles at different times are highlighted.The loss term is lα = 0 and the ratio between slowing down and sawtooth period is τ sd0 /τsaw = 3.

Figure 5 .
Figure 5. Blue lines: snapshots, at time t * = 5.5τsaw, of the time-dependent distribution functions as functions of normalized velocity, obtained for constant τ sd and time-dependent source term.Black dashed lines: the corresponding, standard slowing-down distributions functions obtained for constant slowing-down time and source time.In (a), no losses are assumed, lα = 0.In (b), lα = 0.5.For all curves, τ sd0 /τsaw = 3.Only values of v ⩾ vc are considered (vc is indicated by the vertical dotted line).

Figure 6 .
Figure 6.Blue lines: distribution function with non-constant source plotted for different times during the sawtooth ramp.The critical velocity vc is showed with vertical dotted line.The loss exponent considered is lα = 0.5 and τ sd0 /τsaw = 3.The times displayed are t * = 5.1τsaw (a), t * = 5.5τsaw(b) and t * = 5.9τsaw (c).The dashed curves correspond to the standard slowing-down solution with losses (f ∝ v −3(1−lα) ).

Figure 7 .
Figure 7. Time dependent distribution function as a function of normalized velocity at time t * = 5.5τsaw, obtained for varying τ sd and constant source term.The distribution function is plotted up to the critical velocity showed with vertical dotted line.The new time dependent distribution function (solid blue) is plotted for ratio τ sd,rel /τsaw = 3, together with the slowing-down obtained for a constant slowing-down time τ sd,rel (dashed black).

Figure 8 .
Figure 8. Snapshot, at time t * = 5.5τsaw, of the time-dependent distribution function as a function of normalized velocity, obtained varying in time both τ sd and the source term.The distribution function is plotted up to the critical velocity showed by the vertical dotted line.The new time dependent distribution function (solid blue) is shown for ratio τ sd0 /τsaw = 3, together with the standard slowing-down distribution obtained for a constant τ sd,rel and constant source S rel (dashed black).The combination of the two effects leads to a distribution function with piecewise ∂f/∂v > 0.

Thermal energy conservation determines the value of p rel , ˆr2 mix 0 dr 2 p
(r, t) = p rel r 2 mix , (A.11)