Experimental evaluation of avalanche type of electron heat transport in magnetic confinement plasmas

Turbulent transport is undoubtedly important in magnetic confinement plasmas. It has been reported that a lot of transport dynamics are not satisfying the local diffusive models. Here, we report the quantitative measurement of electron heat flux associated with ballistic propagating long-range transport events, which is considered to be a component of avalanches. In addition, we show the first observations of the substantial impact of avalanche-driven transport on profile resilience (or profile stiffness) observed in JT-60U. We found that, in the channel of the electron heat flux, the ratio between the increment of the avalanche-driven component to that of the total plasma heating becomes dominant (∼80%) in the case of the high-heating limit. This suggests a possible role for avalanche-driven transport to induce profile resilience, which has been evidenced by flux-driven simulations.


Introduction
Turbulent transport is ubiquitous in nature and often a main contributor of cross-field transport in magnetic confinement plasmas. Since plasma particles are subjected to cyclotron motion, the characteristic scale length of turbulence is thought to scale with the cyclotron radius, which is much smaller than the system size. Therefore, the turbulence-driven fluxes a Current address: Institute of Advanced Energy, Kyoto University, Uji 611-0011, Japan * Author to whom any correspondence should be addressed.
Original content from this work may be used under the terms of the Creative Commons Attribution 4.0 licence. Any further distribution of this work must maintain attribution to the author(s) and the title of the work, journal citation and DOI. have been conventionally estimated by applying the local mixing length theory in order to yield the diffusive transport relations, where the coefficients in flux-gradient relations are expressed by the local parameters [1]. However, experiments have shown that the transport phenomena do not always follow the local diffusive nature (i.e. so-called non-local transport problems) [2]. For example, the following observations have been reported: the spontaneous rise of core temperature by edge cooling [3][4][5], the rapid and ballistic response of transport fluxes [6,7] and hysteresis in the flux-gradient relations, which is not explained by the local parameters [8,9].
Stimulated by the non-local transport phenomena, numerous possible explanations have been proposed and discussions have taken place in the magnetic fusion community [10][11][12]. One of the main concepts to understand non-locality is known as an avalanche type of transport, which has been proposed since 1995 [13][14][15] and is still developing [16][17][18]. Avalanches are intermittent transport events that propagate sequentially to neighbors through local critical excitations, which can induce long-range transport [19]. A simple example of an avalanche is provided by the sandpile model, which is a typical concept of the self-organized criticality (SOC) [20]. Avalanches with intermittent and mesoscopic nature have been observed in numerous numerical simulations [21][22][23][24][25] and in some experiments [26][27][28]. In particular, from the fluxdriven simulations, it has been reported that avalanches play a dominant role in forming global plasma profiles, e.g. resilient pressure profiles [21,22] and E × B staircases [25]. Here, we note that [21,22] suggest that the resilient pressure profile or stiffness profile, is maintained through avalanches exhibiting non-diffusive features, and not by the diffusive transport conventionally explained as the critical gradient transport model [29]. In this way, avalanches are expected to induce substantial impact on the plasma profiles. However, the amount of avalanche-driven transport has never been measured definitively in the experiment. This is because avalanches can only be measured transiently with various magnitudes, with the result that the discrimination of these transport fluxes from the background components has not been easy. Quantitative measurements of avalanche-driven transport are required to validate the formation of global plasma profiles, which has been discussed in the flux-driven simulations.
In this article, we have evaluated an element of an avalanche type of electron heat flux quantitatively and show the impact of this transport on the resilient profile formations in the tokamak device, JT-60U. We observed that the characteristic bursts of density fluctuation intensity measured in the reflectometer are synchronized with electron temperature perturbations. By referencing those bursts, an electron temperature perturbation is quantitatively extracted by conditional averaging. The extracted electron temperature perturbation traverses about half (or more) of the plasma radius at a velocity of a few tenths of diamagnetic drift velocity in the manner of joint reflection symmetry [13]. Since the time scale of this perturbation is in the range of 1/ f power-law scaling, which is one of the features of the SOC, we consider the perturbation to be one of the components of electron heat avalanches. The amount of electron heat flux driven by such a long-range avalanche is substantial, and the ratio between the increment of avalanche heat flux and that of mean heat flux shows a sudden increase, when the mean heating power becomes high. In the high-heating limit, the avalanche-driven fluxes reach 80% of the incremented heating power, which can sustain the plasma profiles resilient to the fixed gradients. This discovery shows the essential role of avalanches in global profile formations in tokamak plasmas.

Experimental setup
JT-60U is a large tokamak device for magnetic confinement experiments [30]. The deuterium plasma is confined to a magnetic field of single-null divertor configuration with a toroidal magnetic field B ∼ 3.72 T, a plasma major radius R ∼ 3.3 m and a plasma minor radius a ∼ 0.83 m. In this experiment, plasmas in reversed-magnetic shear (RS) configurations are sustained under neutral beam (NB) heating, and the power of NBs is changed by discharges. The discharge waveforms are shown in figures 1(a)-(e). At the beginning of plasma current ramp-up, tangential NB injections (4 MW) start at t = 3.45 s with no-momentum input (so-called balance injection), as shown in figures 1(a) and (b). The RS configurations are achieved by starting the plasma heating at the current ramp-up phase to reduce the current diffusion to the core region. At t = 4.55 s, heating power scans are provided by perpendicular NBs, and total heating powers reach 8, 10 and 11 MW with respect to discharges of E048162, E048163 and E048166. Here, the control parameters, except the power of perpendicular NBs, are same in every discharge. As NB power increases, the diamagnetic stored energy (figure 1(c)) is saturated, while D αline intensity from the divertor (figure 1(d)) is increased. This suggests increases in energy loss with respect to NB power. In addition, since the D α -line intensity is relatively steady, there are no H-mode transitions that can strongly affect the core plasma transport [31]. Because the electron density is low (n e < 1 m −3 ), the threshold power of the H-mode requires higher power injection [30]. Here, the diamagnetic stored energy and the D α -line intensity are increased and decreased at t ∼ 6 s, respectively. This relates to the transient increases of temperature when the q min is close to the rational surface  [32,33].In these discharges, the q min across 4 around t ∼ 6 s. Thus, we focused on the analysis before this transition phase.
Plasma profiles at t = 5.2 s are shown in figures 2(a)-(c). Here, ρ indicates the normalized minor radius defined by the plasma volume average. As shown in figure 2(a), a safety factor (q) measured by motional Stark effect (MSE) spectroscopy shows RS configurations. Electron density (n e ) and electron temperature (T e ) profiles are measured by the Thomson scattering (TS) system, as can be seen in figures 2(b) and (c). Regardless of the increases in NB power, the plasma profiles are quite similar in every discharge. The dotted line in figure 2(c) indicates the T e profiles averaged in every discharge at t = 4.5 s (4 MW heating), just before the perpendicular NB injection. Since T e increases from 4to 8 MW, the profile resilience (or profile stiffness) [34][35][36] clearly appears at 10 and 11 MW discharges.
In this experiment, a conventional O-mode reflectometer with proving frequency fixed at 34 GHz is provided for analysis [37]. The corresponding cut-off electron density at 1.4 × 10 19 m −3 is shown as a dashed line in figure 2(b). The probing wave is launched to the plasma almost perpendicular to the cut-off layer However, there is a finite incident angle of ∼3 degrees, which is estimated by the ray-tracing code. In addition, since the beam radius of the probing waves ∼15 cm is large, our reflectometer system is dominated by the multiple scattering [38]. As discussed in [38,39], the reflectometry is predominantly sensitive to the long-wavelength fluctuations (k ⊥ ρ i < 1), which is larger than the typical turbulence scale. Amplitude and phase signals from the reflectometer are obtained through a quadrature detector as in-phase and quadrature-phase (I/Q) signals. The I/Q signals are measured at a sampling rate of 1 MHz.
The electron cyclotron emission (ECE) diagnostic is equipped to obtain a temporal evolution of local electron temperature. The detection of ECE is provided by a heterodyne radiometer system, which provides good S/N ratio [40]. The sampling rate of the ECE diagnostic is 50 kHz and the bandwidth of the intermediate frequency (IF) is 500 MHz. The measurement locations of ECE are limited at ρ ∼ 0-0.75. The ECE data are calibrated by the TS data.

Extraction of electron temperature perturbations by referencing the reflectometer signal
The characteristic nature of density fluctuations is observed by the O-mode reflectometer. Since the electron density profiles of each discharge are similar, the cut-off locations of each discharge are ρ ∼ 0.4 at t ∼ 4.8 s, as shown in figure 3(a). Short-time-Fourier transform with a time window of 500 μs is performed on the I/Q signals, as shown in figures 3(b)-(d). At around t ∼ 4.8 s, a spikey increase in auto-power spectrum density (APSD) in negative frequency ranges (the circle in the figure) appeared especially in the 10 and 11 MW discharges. We call this spontaneous enhancement of intensities of fluctuation in wide frequency ranges BF. The BFs can be defined as An important feature of BFs is that they are synchronized to the electron temperature perturbations. Figure 4(a) shows an envelope of I/Q signals band-pass filtered in −250 k < f < −20 kHz (black) and integrated APSD for −250 k < f < −20 kHz (red). Note that both the signals represent the BFs. The ECE diagnostic provides normalized electron temperature fluctuations,T e /T e = (T e − T e )/T e , where T e indicates the time-averaged electron temperature including the trend less than 50 Hz components. Here,T e /T e is arranged asT e /T e + ρ, where ρ shows the normalized minor radius of each measurement location. As can be seen in figure 4(b), temperature fluctuations at around 0.4 < ρ < 0.7 are somewhat synchronized with the BFs observed at the cut-off locations of ρ ∼ 0.4. These dynamics indicate global features. However, the magnetic fluctuations measured by the saddle loop coil are not correlated to the BFs, as can be seen in figure 4(c). Thus, the T e fluctuations are not due to the magnetohydrodynamic activities.
To investigate the characteristics of T e fluctuations, Fourier analysis is performed on the ECE signals. The APSD of electron temperature fluctuations at ρ ∼ 0.67 is shown in figure 5. Here, the APSD of ECE shows two characteristics, i.e. one is a low-frequency region ( f 2 kHz) that scales with 1/ f , and the other is a high-frequency region ( f 2 kHz) that is at constant power with some spikes. This is because the high-frequency region suffers from thermal noise and electronic noise. The electronic noise is estimated by calculating the APSD before the plasma breakdown (t = 1.0-1.7 s), which is shown as a gray line in figure 6. Similar spikes seen in the spectra of the ECE signal with and without plasma indicate the inherent electronic noise. Since the thermal noise can be reduced by correlating independent ECE signals [41], we also calculate the cross-power spectrum density (CPSD). The orange line in figure 6 shows the CPSD ofT e between ρ ∼ 0.67 and Indeed, the constant power of CPSD for f > 2 kHz is reduced by an order of magnitude compared to that of the APSD. Thus, the constant power (or white noise) seen in f 2 kHz is due to the thermal noise. Since these noises are not affecting the low-frequency region, the 1/ f power-law scaling in the range of 0.05 k f 2 kHz seen in the APSD and CPSD is reasonable. The power-law dependence ofT e suggests the existence of avalanches that appear in the SOC system [20]. The 1/ f power-law dependence of theT e spectrum indicates the self-similarity of theT e perturbations, i.e. the perturbations occur similarly in every scale. This is a similar characteristic to that previously observed in electron heat avalanches in tokamak plasmas [26,27]. As can be seen from figures 4(a) and (b), the time scale ofT e /T e synchronizing to the BFs is ∼10 ms, which is in the range of the 1/ f power law.
Since the BFs synchronize with T e perturbations, the BFs can be used as a reference signal for quantifying the spatiotemporal evolution of electron temperature fluctuations. The conditional averaging is performed on the ECE signals by referencing the BF emergence (see appendix A.1). The analysis is performed at t = 5.0-5.7 s, where the plasma is relatively steady. Here, the number of signals sampled from t = 5.0-5.7 s is adjusted to ∼70. This indicates that conditional averaging is performed for theT e with a characteristic time scale of ∼10 ms. The result for the 11 MW discharge is shown in figure 6, which indicates the conditionally averaged temporal evolution of the envelope of I/Q signals at −250 k < f < −20 kHz andT e /T e . Here, τ = 0 represents the time at the BF emergence, determined by the threshold of its time derivative. In figures 6(b)-(i), the green line indicates the low-pass 0.5 kHz components of the conditionally averagedT e /T e . When BF occurs (τ > 0), theT e /T e decreases (voids) at ρ 0.49, while it increases (bumps) at ρ 0.54. Here, the pulse width of the voids and bumps are almost constant. The delay times of the voids and bumps with respect to the signal at ρ = 0.49 and ρ = 0.54 are shown by blue markers in figures 6(b)-(i). These delay times are estimated by cross-correlation analysis. The nearly constant delay times with respect to the distance indicates the ballistic radial propagation of the voids and bumps to the inward and outward direction, respectively. The bi-directional propagation of voids and bumps expected in avalanches is known as joint reflection symmetry [13]. From the linear fittings, the propagation velocity of the voids and bumps is estimated to be −166 ± 4 m s −1 (radially inward) and 108 ± 7 m s −1 (radially outward), which are a few tenths of the diamagnetic drift velocity. The estimated propagation speed is much faster than the typical diffusive time speed of χ e /a ∼ 1 m s −1 , where χ e ∼ 1 m 2 s −1 and a ∼ 0.83 m are the diffusive coefficient of electron heat estimated by power balance analysis and plasma minor radius, respectively. The radial propagation of avalanches close to the diamagnetic drift velocity is anticipated from the numerical simulations [23]. In addition, the amplitude of the voids and bumps seems constant, which should be decreased in the framework of local diffusive transport. This can be clearly seen from the Lissajous diagrams against the BFs andT e /T e shown in figure 7. The rotating direction of Lissajous for clockwise and counterclockwise indicate the voids and bumps, respectively. The amplitude of the voids and bumps is shown as a width of the y-axis in the diagrams. For voids (ρ 0.49), the amplitude is not decreased but even increased with propagation. Similarly, the amplitude of the bumps (ρ 0.54) is increased or even constant with propagation. The tilting of Lissajous also indicates the propagation of the voids and bumps. The finite tilting of the ellipse seen at ρ = 0.49 (negative correlation direction) and ρ = 0.54 (positive correlation direction) is gradually degraded towards the inside and outside radial direction. This indicates the inward/outward propagation of voids/bumps.
In addition to the above observations, we have seen the 1/ f power law at 0.05 k f 2 kHz in the ECE spectrum. Here, the extracted T e perturbation has a specific time scale of ∼10 ms. Thus, we have only extracted a component of T e perturbation (relatively large scale) that belongs to the 1/ f power law. From this point of view, the extracted T e perturbation is one of the elements of avalanches that is composed of various scales with self-similar nature [19]. Next, we discuss the heating power dependence of electron temperature perturbations obtained by conditional averaging.  Note that in the 11 MW case, the boundary of the voids and bumps is almost coincident to the peak of R/L T e , which may suggest that the electron heat avalanches relax the local peak of electron temperature gradients. In addition, the radial scale of T e perturbations elongates beyond the peak of R/L T e , where there is the free energy source of turbulence. Such a long-radial transport is not explained by the local diffusion model.
The increase in electron heat avalanches with respect to the NB-heating power can also be seen in the power spectra of the ECE. Figure 9(a) shows the APSD ofT e /T e at ρ ∼ 0.67. As can be seen from figure 9(a), the power at the f ∼ 0.1 kHz component, which is coincident to the time scale of conditionally averaged T e perturbations, is increased with NB heating. In addition, the power-law scaling at 0.05 k f 2 kHz seems to be changed by NB heating. The exponent of power spectrum (P( f )) is defined as, P( f ) ∝ f α , where α can be obtained by a linear fitting of a logarithmic plot. Figure 9(b) shows the radial profiles of α with regard to the NB heating. Note that the power spectrum analysis for ρ > 0.4 is limited by the noise, whereas the noise is suppressed by the conditional averaging because the S/N ratio increases with 1/N, where N is the ensemble number. As can be seen in figure 9(b), the powerlaw exponent is close to −1 with heating power increases, which suggests that large-sized avalanches are activated. This is consistent with the result obtained byconditional averaging.

Characteristics of BFs
So far, we have discussed the characteristics of T e perturbations and found that they are large-scale components of electron heat avalanches. Now, we focus on the characteristics of BFs, which are used as a reference signal of conditional averaging to extract the T e perturbations.
First, we discuss the power spectrum of the BFs. The APSD of I/Q signals in the negative frequency region and that of BF signals defined by the envelope (i.e. envelope of I/Q signals ranging from −250 k < f < −20 kHz) are shown in On the other hand, the spectra of envelopes at 0.05 k < f < 2 kHz showe power-law scaling (∝ f −0.6 ) at 8 MW, while it disappears and the power is increased at 10 MW and further increased at 11 MW with a quasi-coherent peak. Note that the quasi-coherent peak at f ∼ 0.1 kHz is not observed in the I/Q signals. Since phase components of I/Q signals are steady, the BF is the result of envelope modulations of carrier waves at −250 k < f < −20 kHz. The appearance of quasi-periodic BFs could be a feature of the SOC system, including the strong diffusion process [42,43], which will be discussed again in section 5.
The characteristic variation of BFs against the NB heating is also observed in the radial correlation. In this experiment, in addition to the 34 GHz reflectometer, the O-mode reflectometer with probing frequency of 36 GHz is simultaneously launched to the plasma [37]. The two-point crosscorrelation analysis is performed on the fluctuations constituting BFs. Figures 11(a) and (b) show the auto-correlation function (ACF) and cross-correlation function (CCF) of bandpass filtered (−250 k < f < −20 kHz) I/Q signals between the 34 and 36 GHz reflectometer. As the NB power increases, the absolute value of the CCF becomes larger, and those at 10 and 11 MW are above the noise level. The time duration of the finite CCF is comparable to that of the ACF. Here, the radial distance of the cut-off layers of two probing waves are ∼4 cm, which is longer than the Airy width ∼1.2 cm. Thus, the radial correlation length of fluctuations constituting BFs could be increased with NB power, whereas the amplitude saturates at 10 MW, which is shown in figure 10(a).
Since the BF is constituted by the negative frequency components (−250 k < f < −20 kHz), we consider the effect of Doppler shifts. Figure 12(a) shows temporal evolution of conditionally averaged APSD of the I/Q signals (−250 k < f < −20 kHz and 20 k < f < 250 kHz components) and Doppler shift ( f D ) that is calculated by a weighted average of the spectrum. When the BF occurs, f D increases from −20 to −70 kHz. The spectra before (τ = −1 ms) and after (τ = 1 ms) the BF are shown in figure 12(b). Here, the reflectometer is sensitive to a lot of scattering and reflecting components (m = −2, −1, 0, 1, 2, . . . ) because the probing beam radius is large. Thus, basically, it is difficult to determine the sensitivity of the reflectometer used in this experiment. The observation of negative Doppler shift even before the burst suggests that the finite tilting angle of incident beams, which is observed in the ray-tracing code, causes the difference in sensitivity of m = −1 and 1 scattered waves.
Finally, the possible mechanism of Doppler shift is mentioned. The Bragg backscattering condition gives the Doppler shift as 2π f D = v · k = v r k r + v θ k θ + v ϕ k ϕ , where r, θ and ϕ represent the radial, poloidal and toroidal direction, respectively [44]. Here, we neglect the v ϕ k ϕ term because k ϕ k r , k θ can be expected. Since avalanches can produce finite radial velocity, here we consider the poloidal and radial components, respectively. With regard to the poloidal components, scattering waves of poloidal wavenumber are estimated to be k θ = 2k i sin α ∼ 0.74 cm −1 , where wavenumber of incident beam k i ≈ 7.1 cm −1 and tilting angle α ≈ 3 degree are  [38,39] and observations of poloidal velocity with an order of v θ ∼ v E×B ∼ 1 km s −1 [45]. On the other hand, the radial wavenumber of scattering waves can be obtained as k r = 2k i cos α ∼ 14.3 cm −1 . The radial velocity is estimated to be v r = 2πΔ f D /k r ∼ −220 m s −1 , which is coincident to the measurement that voids propagating with −166 m s −1 (inward direction) at the cutoff locations. However, k r ∼ 14.3 cm −1 is above the maximum value of the radial wavenumber evaluated from Airy width (∼1.2 cm) as ∼5.2 cm −1 . Thus, based on simple estimates, the radial Doppler shift is unexpected. During the BFs, the possibility of Doppler shift in the poloidal direction can be the result of the changes in the fluctuations of wavenumber and/or poloidal velocity. The relationship between avalanches, turbulence and momentum transport is an interesting subject to discuss with regard to intrinsic rotation [46], and detailed experiments are required for future study.

Avalanche-driven electron heat flux and its contribution to profile resilience
In this section, we have evaluated the electron heat flux of an element of avalanches, which is extracted by conditional averaging. Since conditional averaging reduces the noise level of ECE signals, the electron heat flux (q e ) can be evaluated directly from the energy conservation equation. Details  of the analysis are provided in appendices A.2 and A.3. Figures 13(a)-(c) show the time evolution of the conditionally averaged BF and electron density normalized heat flux (q e /n e ) with respect to NB power. The changes in BF and q e /n e becomes strongly synchronized when the heating power increases. In particular, a bursty increase in q e /n e is clearly observed in the 11 MW case. Note that the decreases in q e /n e at τ < 0 indicate the increase in the T e gradient at ρ ∼ 0.5 before the burst, which can be seen in figure 6. Considering how the number of extracted avalanche components affect the total transport, the time-averaged quantities of the electron heat flux are considered. Figure 14(a) indicates the radial profiles of mean electron heat flux (q e,mean ) and time-averaged effective electron heat flux driven by the bursty components (q eff e,burst ), which indicates avalanche-driven fluxes occupied in the q e,mean (detailed definitions are provided in appendix A.2). The ratio of q eff e,burst /q e,mean is shown in figure 14(b). Here, the error bars are estimated by considering error propagation of several uncertainties (time derivative of electron temperature, electron density and electron heating source profiles). It is found that the extracted avalanchedriven transport amounts at most to ∼10% of the total electron heat transport in the case of 11 MW heating. Since the avalanche-driven transport increases abruptly when the NB power reaches 11 MW, the increments of the q eff e,burst are compared with the increments of the heating power. Figure 14(c) shows the ratio of the incremented quantities (from 8/10-10/11 MW) of effective to mean electron heat flux (Δq eff e,burst /Δq e,mean ). The incremented avalanche-driven transport is almost equivalent (∼80%) to the increased NB power from 10-11 MW. The ratio of incremented quantities is radially averaged at ρ = 0.4-0.6, as shown in figure 15. The occupancy of the extracted avalanche components to the incremented heating power is increased with increases in mean heat flux. This shows that the profile resilience at the strong heating limit (11 MW) can be sustained by the avalanching transport.

Discussion
In this study, we have extracted the electron heat avalanches that are correlated with the BFs, which is necessary to discuss here. In the model of avalanching transport, turbulence also spreads to other radial regions [19,47,48], which can enhance the local fluctuation amplitude at the linearly stable zone [49,50]. In addition, the correlation length of fluctuations can be modified by the enhancement of the turbulence spreading [19], which is consistent with the tendency observed from figure 11. However, we cannot conclude that the BFs are the result of turbulence spreading because the reflectometer system used in this experiment is basically sensitive to the long-wavelength fluctuations, which are larger than the typical scale of turbulence.
The changes in the power spectra of BFs from powerlaw scaling at 8 MW discharge to a quasi-periodic feature at 11 MW discharge should also be mentioned. If the BFs are composed of density avalanches described in the classical SOC systems, they are not expected to appear regularly. The situation could be related to the SOC system including the diffusion, where the 1/ f type of avalanches can be replaced by the quasi-periodic type with enhanced diffusion [42,43]. Indeed, in some flux-driven simulations, quasi-periodic avalanches occur with increased heat conductivity [51], and the quasiregular bursts are observed co-existing with the 1/ f type of avalanches [52]. However, once again, we cannot conclude that the BFs are density avalanches because the measurement points are limited, and thus their radial profiles and propagations cannot be discussed. Another possibility is that the BFs can be considered to be long-range fluctuations nonlinearly coupling to the turbulence [53]. Detailed measurements of the BFs are necessary for future works.
It is necessary to mention again that the evaluated electron heat flux is only a component of electron heat avalanches. As shown in the power spectrum of ECE, avalanches contain 1/ f scaling, indicating scale invariant nature. In our study, conditionally averaged avalanches with a characteristic time scale of ∼10 ms are provided for estimating the radial heat flux, which reaches ∼10% of total transport (without this process, we could not distinguish avalanches from noise). Since the 1/ f region in the power spectrum is wide, it can be said that the total amount of avalanche inducing transport is more than this estimation, which is predicted in the global gyrokinetic simulations [54]. After all, the avalanche-driven transport is not negligible, substantial and important, especially in hightemperature core plasmas. The evaluation of total avalanchedriven fluxes is necessary for future study for comparison to the local diffusive transport components.

Summary
We have shown experimental observations of the actual role of avalanche-type electron heat transport in the resilient profile observed in JT-60U. Since avalanches are synchronized to the burst-like increases of density fluctuations, we can conduct conditional averaging to reduce the noise level of electron temperature perturbations and successfully extract the specific time scale of electron temperature perturbations quantitatively. The conditionally averaged T e perturbations are (i) extracted from the fluctuations including the 1/ f -type of power-law scaling, (ii) characterized by bi-directional propagation of voids and bumps satisfying the joint reflection symmetry, (iii) propagating with constant amplitude, pulse width and velocity with a few tenths of diamagnetic drift velocity and (iv) arriving at the long-radial distance that is out of the peak of the T e gradient. Based on these observations, the conditionally averaged T e perturbations are considered to be one of the components of avalanches. We found that these types of avalanches increase with increased NB power. The electron heat flux driven by the element of avalanches is first quantitatively evaluated, and we found that the amount of this avalanche-driven transport reached 10% of the total electron heat transport at the high-heating limit. This value is almost equivalent to the incremented NB-heating power from 10-11 MW. It has been argued from flux-driven simulations [21,22] that the avalanches could cause the profile resilience observed in high-temperature core plasmas. The results shed light on the transport process of magnetic confinement plasmas, which has been discussed with unresolved enigmas for a long time. Due to the conditional averaging, the noise level of ECE signals is dramatically reduced. However, the residual noise could still make the heat flux analysis difficult. We performed the low-pass filtering with f < 0.5 kHz on the conditionally averaged ECE signals (green line in figures 7(b)-(i)). The root mean squared errors (RMSEs) are defined from the residuals between the conditionally averaged original signals and the low-pass filtered signals. Therefore, analysis using the conditionally-averaged ECE signals is performed through the low-pass filtered signals with uncertainties defined by the RMSE.

A.2. Estimation of avalanche-driven electron heat flux
The electron heat flux is calculated at t = 5.0-5.7 s, where the plasma is relatively steady. The spatio-temporal evolution of electron heat flux (q e (r, t)) can be obtained from the energy conservation equation [55] as q e (r, t) = 1 S(r) r 0 P r − 3 2 ∂(n e (r ,t)T e (r ,t)) ∂t dV, where P(r) is the electron heating source profiles and V and S are the plasma volume and plasma surface area, respectively. The electron heating source includes NB heating, ohmic current heating, energy transfer between electron-ion collisional equipartition and radiation loss. The NB-heating deposition profiles are calculated by the orbit-following Monte Carlo code (OFMC) [56], and the other heating sources and losses are calculated through the transport analysis code TOPICS [57], every 20 ms at t = 5.0-5.7 s. The time-averaged P(r) and their standard errors are used for the calculations. The standard deviation of the calculated P(r) at t = 5.0-5.7 s is ∼24%, which is comparable to the error bar of P(r) from one time slice. The ensemble of N = 35 is used to calculate the standard error of P(r) to be 24/ √ N ∼ 4%. Here, we have assumed the steady state of heating sources and losses, which is discussed in appendix A.3. The time derivative terms of electron temperature are estimated from the low-pass filtered conditionally averaged data, and their uncertainties are evaluated from the RMSE introduced in the conditional-averaging technique. Here, the level of electron density fluctuations, which is evaluated using the conditionally averaged interferometer data, is much smaller than that of the electron temperature. Thus, electron density is assumed to be constant in time, n e (r, t) = n e (r), which is estimated by the Gaussian process regression of TS data [58]. The error bars of q e (r, t) are calculated by considering the error propagation of uncertainties in P(r), n e (r) and ∂T e (r,t) ∂t . Time-averaged quantities of electron heat flux are also introduced. The long-time-averaged mean electron heat flux (q e,mean ) is described as q e,mean (r) = lim where T = T bg+burst + T bg , and T bg+burst and T bg are the duration of time windows for bursty q e components and background q e components, respectively. The q e,bg+burst and q e,bg are time-averaged q e in the respective time windows. A schematic view of the analysis is shown in figure A1. By introducing the net bursty components of electron heat flux, i.e. q e,burst = q e,bg+burst − q e,bg , the q e,mean is simply described as, q e,mean (r) = T bg+burst T q e,burst (r) + q e,bg (r) ≡ q eff e,burst (r) + q e,bg (r).
Here, q eff e,burst (r) = T bg+burst T q e,burst (r) indicates the bursty components of effective electron heat flux, which occupies q e,mean . The q eff e,burst is simply obtained as q eff e,burst = q e,mean − q e,bg . In the analysis, conditionally averaged electron temperature perturbation provides q e (τ ). The time windows of bursty q e (τ ) components are set to 0 < τ < 3 ms, which is equivalent to the full width at half maximum of BF signals. The time windows of background q e (τ ) components are set to −4 < τ < 0 ms and 3 < τ < 7 ms, and thus, T bg+burst = 3 ms and T bg = 8 ms. The error bars of time-averaged quantities are estimated from the error of q e (τ ).

A.3. Remarks on heating source and losses
The assumption of neglecting the time variation of heating sources and losses in the heat flux analysis is validated. Here, the actual values of the electron heating sources and losses are as follows (11 MW discharge, ρ = 0.50): NB heating is 48 kW m −2 , ohmic current heating is 6 kW m −2 , equipartition of collisional energy transfer between ion and electron is 4 kW m −2 and radiation loss is −6 kW m −2 . Thus, the ratio of each heating source and losst to the total power (52 kW m −2 ) is 92%, 11%, 8% and 11% for NB heating, ohmic current heating, equipartition of ion-electron collisional energy transfer and radiation loss, respectively. These values are larger or comparable to the avalanche-driven electron heat flux ∼10%.
The time scale of the calculations of the heating sources and losses through numerical codes are restricted by the measurement of TS, which measures T e and n e every 20 ms. Thus, the heating sources and losses are also evaluated with 20 ms time steps, which is larger than the time scale of the avalanches ∼10 ms. However, due to the following reasons, it is reasonable to assume that these heating sources and losses are constant in time during the calculation of the avalanche-driven electron heat flux.
Considering heating sources, the time scale of NB heating is ∼200 ms, which is the relaxation time estimated from the OFMC calculation. The time scale of ohmic current heating is dependent on the current diffusion time, which is ∼1 s. The equipartition of ion−electron collisional energy transfer time is ∼300 ms. Thus, the time scale of NB heating, ohmic current heating and the equipartition term are much larger than that of the bursts.
The radiation loss from the plasma core (ρ = 0.50) is dominated by the cyclotron radiation and Bremsstrahlung. The power of the cyclotron radiation (P cy ) and the Bremsstrahlung (P Brems ) are described as P cy ∝ B 2 n e T e and P Brems ∝ n e √ T e Z 2 n z , where Z and n z indicate the charge and density of impurities, respectively. During the bursts, the changes in n e fluctuation can be neglected because it is much smaller than theT e /T e ∼ 0.01. Therefore, the P cy can be changed during the burst on the order of 0.01, which is equivalent to ∼−0.06 kW m −2 (0.11%). This is much smaller than the avalanche-driven electron heat flux (10%). On the other hand, the variation in P Brems is evaluated from the charge exchange spectroscopy (CXRS) measurements. CXRS measures the carbon impurity (C VI, 529.05 nm) and the Bremsstrahlung for 2.5 ms time steps. To improve the time resolution of the measurements, we used the conditional reconstruction technique [59]. By taking the relative time difference of the measurements with respect to the BFs, the effective time resolution can be improved by using many repeated BFs in the steady-state phase. As a result, we found thatP Brems /P Brems T e /T e . Thus, changes in the P Brems are less than the order of 0.01, which amounts to ∼−0.06 kW m −2 (0.11%). Thus, the variation of radiation loss of Bremsstrahlung is also negligible compared to the avalanche-driven transport.