Global Alfvénic modes excitation in ohmic tokamak plasmas following magnetic reconnection events

A possible triggering mechanism of Alfvén waves (AWs) in tokamak plasmas, based on localized perturbations induced by magnetic reconnection events, is discussed in the framework of nonlinear viscoresistive 3D magnetohydrodynamics (MHD) modeling. Numerical simulations are performed with the SpeCyl code (Cappello and Biskamp 1996 Nucl. Fusion 36 571) that solves the equations of the viscoresistive MHD model in cylindrical geometry. We investigate a ohmic tokamak configuration where the m = 1, n = 1 internal kink mode (m is the poloidal mode number and n is the toroidal mode number) undergoes a complete reconnection process. An in-depth investigation of the process shows a spatio-temporal correlation between the velocity perturbations associated with the reconnection and the excitation of the shear AW in the core region and the global Alfvén eigenmodes, both with dominant m = 1, n = 0 periodicity. In particular they are observed to emanate from the outflow cones of the reconnection layer associated with the internal kink. The excitation mechanism described in this paper could explain the observations of Alfvénic fluctuations in the absence of energetic ions in several tokamak experiments documented in the literature and could contribute to AWs excitation in general, even in the presence of fast particles. This result shares similarities with analogous study in reversed-field pinch (RFP) configuration (Kryzhanovskyy et al 2022 Nucl. Fusion 62 086019) where AWs were found to be excited by the RFP sawtoothing.


Introduction
The study of the properties of Alfvén waves (AWs) in modern toroidal devices is a crucial contribution to reactor-relevant physics [1,2].AWs are commonly observed in the presence of energetic ions in the MeV energy range (produced by neutral beam injection, ion cyclotron resonance heating, or even fusion-born alpha particles) which can satisfy the conditions of effective resonance and energy exchange with such waves [3].However, modes with Alfvénic frequency scaling have also been detected in ohmic discharges (i.e. in the absence of fast particles) in numerous tokamak experiments, including JET [4], EAST [5], TFTR [6], FTU [7], TEXTOR [8], ASDEX Upgrade [9], TUNAM-3M [10], COMPASS [11], MAST [12], Globus-M2 [13] and SUNIST [14].In most of the above-mentioned experiments (EAST, TFTR, TUMAN-3M, COMPASS, MAST, Globus-M2 and SUNIST), AWs were detected following relatively long-timescale MHD events in the plasma, such as sawteeth, internal reconnection events (IREs) and edge localized modes (ELMs).In some cases (in JET [4], FTU [7,15] and TEXTOR [8]) the AWs excitation was observed in presence of magnetic islands due to selfinduced low frequency MHD modes (for example the m = 2, n = 1 tearing mode, where m and n are the poloidal and the toroidal mode numbers respectively) or applied resonant magnetic perturbations.The dominant toroidal mode number for AWs in ohmic plasmas in JET, TFTR, MAST and Globus-M2 was generally found to be n = 0, which excludes gap modes (e.g.toroidal Alfvén eigenmodes, TAEs) for those discharges.The possible mechanisms of AWs excitation in ohmic plasmas have not been comprehensively investigated yet as currently the research on this subject is mostly centered around the excitation of AWs by fast particles (either energetic ions or runaway electrons).Nevertheless, some possible processes for the excitation of AWs in tokamaks in the ohmic discharges without energetic particle injection have been investigated in [16], where a correlation between high-frequency mode activity and relatively long-timescale MHD events in the plasma, such as IREs or ELMs, was proposed.Furthermore, the possible destabilization of Alfvénic modes in the presence of magnetic islands, even without fast particles, was investigated in [17,18].Finally, we mention that the presence and role of AWs during tokamak disruptions has been studied through nonlinear modeling in ohmic plasmas [19].
In this paper, by means of nonlinear 3D MHD numerical simulations, we describe a mechanism of AWs excitation by magnetic reconnection events, associated with the dynamics of the internal m = 1, n = 1 kink mode.Such a process could be responsible for experimental observations of AWs in several ohmic tokamak discharges, and could provide an additional excitation mechanism in more general cases where fast particles are expected to play a major role.The present paper extends to the tokamak case the results of an analogous study performed on the reversed-field pinch (RFP) configuration [20].In that work AWs were found to be excited following quasiperiodic magnetic reconnection events associated with the RFP sawtoothing activity [21].The present paper investigates AWs following the MHD activity of tokamak sawtooth instability, and extends the previous one by showing in detail the mechanism (intense local flow patterns) that gives rise to Alfvénic excitation.The generality of the mechanism suggests that in principle it could apply to any plasma instability that involves magnetic reconnection process.
The paper is organized as follows.In section 2 the employed MHD model and numerical setup are described.In section 3 numerical results showing AWs excitation by magnetic reconnection in the tokamak are presented.In section 4 we investigate in more detail the AWs excitation mechanism in our simulations and we provide some additional analysis regarding the choice of the dissipative parameters and the spectral resolution.A summary and final remarks are given in section 5.

MHD model and numerical setup
The simulations reported in this paper are performed in cylindrical geometry with the nonlinear visco-resistive 3D MHD code SpeCyl [22], solving the following set of MHD equations in the approximation of negligible pressure and fixed density: Lengths are normalized to the cylinder minor radius a, the density ρ to the on-axis ion mass density ρ 0 , the magnetic field B to the initial on-axis magnetic field B 0 , the velocity v to the on-axis Alfvén velocity v A = B 0 / √ µ 0 ρ 0 and time to the Alfvén time τ A = a/v A .In these units, the resistivity η is the inverse Lundquist number, η = τ A /τ R ≡ S −1 , where τ R ≡ µ 0 a 2 /η, and the kinematic viscosity ν corresponds to the inverse viscous Lundquist number, ν = τ A /τ V ≡ M −1 , with τ R and τ V resistive and viscous time scales.The above equations are solved in periodic cylindrical coordinates (r, θ, z) using a finite difference staggered mesh in the radial coordinate r and adopting a spectral formulation in the periodic coordinates θ and z, with the azimuthal (poloidal-like) angle θ going from 0 to 2π and axial (toroidallike) coordinate z going from 0 to 2π R 0 /a, where R 0 and a are the major and minor radius of the corresponding torus, respectively.In this approximation, the azimuthal and axial mode numbers m and n correspond to the poloidal and toroidal mode numbers in a torus.Time advancement is performed by a predictor-corrector type scheme.To sustain the plasma current a uniform induction electric field E = E 0 ẑ is imposed.The plasma boundary conditions at r = a are chosen to represent an ideal, i.e. perfectly conducting, shell.The SpeCyl code was subject of a nonlinear verification benchmark with the finite volume extended-MHD PIXIE3D code [23], as reported in [24].The main tokamak simulation investigated in this paper starts from an axisymmetric circular static ohmic equilibrium with central safety factor q 0 = 0.8 and bell-shaped density profile, as shown in figure 1.This equilibrium is unstable to the (m, n) = (1, 1) internal kink mode responsible for the typical sawtoothing in tokamak plasmas.Furthermore, it is marginally stable to the 2,1 tearing mode.To properly capture the details of the magnetic reconnection event, 32 harmonics of the 1,1 internal kink mode (modes with same geometric helicity h = n/m = 1) are retained in the computation and are initialized with up-down asymmetric perturbations to allow for a complete reconnection as discussed in [25], where the adequacy of 32 harmonics with h = 1 for the study of magnetic reconnections at the q = 1 rational surface with high Lundquist number (S = 10 8 ) was demonstrated.In addition to the harmonics with h = 1, we add to the spectrum a number of harmonics with other helicities including the 2,1, the 4,2 and several harmonics with n = 0.The complete set of 3D Fourier components present in the simulation is shown in red in figure 2. Both the 1,1 and 2,1 modes are initialized with finite amplitude, which, by the three-wave coupling, activates a fully 3D sawtoothing dynamics.We remark here that in a general toroidal geometry the initialization of the 2,1 mode (corresponding in the present simulation to a magnetic island width of 10% of the minor radius) would not be necessary to activate a full 3D dynamics.Indeed, the toroidal equilibrium deformation of the plasma column (Shafranov shift) naturally entails a 1,0 Fourier modulation interacting with the unstable 1,1 mode, that would lead to a 3D dynamics similar to the one described in this paper.
The viscosity ν is assumed to be constant and uniform, while the resistivity η has a radial profile consistent with the above-mentioned initial 1D equilibrium.The dissipation parameters are S = 10 8 and M = 10 6 .The employed radial mesh has 256 points and the aspect ratio of the periodic cylinder is R 0 /a = 4.The time step is 10 −4 τ A , and fields are saved every 0.1 τ A .
We are going to focus in particular on the AWs with m = 1, n = 0 space periodicity.The corresponding expected frequency-spectrum is shown in the last panel of figure 1.It was analytically evaluated in cylindrical geometry from the linearized ideal MHD model (η = ν = 0), cold plasma approximation (pressure p = 0), and perfectly conducting shell conditions.The shear AW (SAW) continuous frequency, also called Alfvén continuum in non-uniform plasmas, is given by the dispersion relation: where k θ = m/r and k z = n/R 0 .As a typical fusion plasma, the m = 1, n = 0 Alfvén continuum frequency is maximum in the core and has a minimum around the edge of the plasma.The global Alfvén eigenmode (GAE) discrete frequencies can be obtained in the local or Wentzel-Kramers-Brillouin approximation from the following quantization condition due to the finite size of the plasma (minor radius a) (for full derivation see [26]) with j = 1, 2, 3 . . .and where is the frequency separation between the GAEs and the continuum, with the following constraint on ω 2 j values In the considered simulation case, a single GAE solution is found below the SAW minimum.

Excitation of AWs by magnetic reconnection in nonlinear tokamak simulation
The dynamics of the simulation introduced above is shown in figure 3. The magnetic and kinetic energies of the Fourier modes included in the calculation are plotted in panels (a) and (b).The dynamics is characterized by the exponential growth of the h = m/n = 1 helicity modes, with 1,1 mode highlighted in black, and a collapse phase during which a complete magnetic reconnection event takes place and brings the central q above 1 making the mode stable.The excitation of the harmonics with geometric helicity m/n ̸ = 1 is achieved due to the three-wave coupling with the 1,1 mode through the marginally stable 2,1 mode (plotted in green).The internal kink collapse is followed by the flattening of the current density and safety factor profiles in the core region (displayed with dashed curves in figure 1).A signature of reconnection events is a spike in the intensity of the kinetic energy, as part of the magnetic energy contained in the plasma is converted into it, while the rest is transported outside by the Poynting vector and dissipated by viscosity and resistivity.This can be seen in panel (c), where the variation in total magnetic and kinetic energies are plotted.
To investigate whether the reconnection process could excite AWs in tokamak ohmic configuration, in particular the 1,0 mode, we analyze its frequency spectrum by computing the continuous wavelet transform (CWT) at fixed radius r/a = 0.3.In the result, plotted in figure 3(d), we can indeed observe a sharp excitation in the spectrum, which coincides with the sawtooth crash, and distinguish at least two stable frequency signals which persist longer than the sawtooth relaxation event itself, undergoing only a weak damping due to the finite resistivity and viscosity in the SpeCyl code.We can further highlight the close temporal correlation between the AWs excitation and the internal kink dynamics by overplotting the spectrum with the characteristic temporal shape of the central q value (pink curve).To get an estimate of the contribution to the spectrum by various Alfvén modes, in the last panel (e) we plot the amplitude of the CWT analysis of the radial velocity field at the two fixed frequencies, corresponding to the Alfvén modes in panel (d).
To identify the AWs excited during the magnetic reconnection event, the fast Fourier transform of the magnetic and velocity components with n = 0 spatial periodicity was computed during the last 500 τ A of the simulation.This time window was chosen in order to avoid most of the intense low frequency fluctuations generated by the kink/tearing modes.The resulting frequency spectrum, displayed in figure 4, was compared with  the expected analytical one for the 1,0 mode (displayed in the last panel of figure 1), leading to the identification of the two main signals as the SAW and the GAE associated to the m = 1, n = 0 harmonic.The SAW, in addition to the visco-resistive damping, is affected by the phase mixing phenomena, due to which it tends to decay very quickly at the radial positions with stronger frequency gradient.This is why the SAW signal can be clearly seen only for r/a ⩽ 0.4, that is inside the q = 1 surface, due to the flattening of its frequency following the magnetic reconnection event (as exemplified in the last panel of the figure 1).On the other hand, the GAE (a global mode just below the SAW continuum minimum) is only affected by the visco-resistive damping and therefore can be observed for longer time even at r/a > 0.4.The gradually weaker signals with increasing frequency in the core correspond to the SAW frequencies of the 2,0 and 3,0 modes.
In order to understand the origin of AWs excitation, we look into the details of the evolution of the internal kink mode and of the related reconnection process.Figure 5 shows the axial component of the current density (first row) and the poloidal component of the velocity field (second row) at five selected time snapshots covering the growth of the kink mode and the final relaxed state.The contour levels (black lines) of the helical flux function with h = 1 helicity (depicting the magnetic surfaces associated with the internal kink mode, as described in [25]) are also shown.The magnetic island associated with the 1,1 internal kink mode grows on the right-hand side of the plots, while a quite elongated current sheet develops on the left-hand side around the X-point (panel (a)).This creates a topology with two magnetic axes, one corresponding to the unperturbed axis-symmetric one and the other related to the O-point of the 1,1 magnetic island (clearly visible in panel (b)).The strong negative current sheet at the X-point, and the positive current in the island O-point are both intuitively understandable as inductive response of the plasma to this shift of current profile [27].The elongated current sheet breaks up to form a number of secondary islands (panel (b)), called plasmoids [28].Such topological structures have been the subject of several studies related to high current/low q tokamak configurations and they are considered as a possible explanation for fast magnetic reconnection in tokamaks [29][30][31][32][33]. Later the plasmoids coalesce (panel (c)) into a single larger secondary island (dubbed 'monster plasmoid' in [33,34]), according to a typical competition between granulation and coalescence occurring in a current sheet (as shown in [35]).For plasmoids to develop the aspect ratio of the current sheet (i.e. the ratio between its longitudinal and transverse extent) must exceed a certain threshold.Since the aspect ratio of the current sheet increases with Lundquist number S, this translates into a threshold for the Lundquist number itself.For viscoresistive MHD simulations in cylindrical geometry, such a threshold was previously reported to be around S = 10 7 [31][32][33].The employed asymmetric up-down initial perturbation leads to complete reconnection for high-S dynamics (S ⩾ 10 7 ) as shown in [25] and seen earlier in [31].Due to the small asymmetry, the monster plasmoid forms slightly displaced toward the top, and the displacement increases quite rapidly (panel (d)) until the island is absorbed into the main 1,1 island.As a result, the complete magnetic reconnection occurs leading to a final state with nested circular magnetic surfaces and q > 1 everywhere in the plasma (panel (e)).The negative current halo (blue color in panel (e)) forms in order to conserve the total current [36].Such a configuration will diffuse on a global resistive time scale and will eventually become unstable again to the internal kink mode, resulting in a quasi-periodic sawtoothing dynamics.
The poloidal velocity perturbation in the second row of figure 5 can be first observed as a small symmetrical burst (in panel g) at the boundaries of the current sheet.A second stronger perturbation is observed later in time (panel (i)) at the Y-points between the original plasma core and the monster plasmoid, it peaks when the latter is completely absorbed into the main 1,1 island and lingers for some time until it is completely dissipated.These two bursts can also be recognized in figure 3(b) as the two peaks of the E k at t ∼ 17 400 τ A and t ∼ 17 700 τ A .These velocity perturbations, together with the corresponding current density perturbations shown in the first row, are at the origin of AWs excitation in our simulations, as detailed in the following section.

Discussion
In this section, we aim to provide a more comprehensive explanation of the underlying mechanism responsible for the self-consistent excitation of AWs in our simulation.Additionally, we will address the impact of the chosen viscoresistive parameters and spectral resolution on the results we presented so far.

Excitation mechanism
To thoroughly investigate the excitation of AWs by reconnection processes and understand the underlying mechanism, we perform a comparative analysis of the temporal evolution of the poloidal velocity component and the corresponding frequency spectrum, as function of poloidal angle and radius.Our approach involves placing 'virtual probes' at different locations within the plasma to measure local fluctuations over time.
To achieve this, we employ the CWT to analyze the v θ field component at specific mesh points following the q = 1 surface (i.e.all points with a fixed radius r/a = 0.4) and the mesh points between the plasma core and its edge (all points with a fixed poloidal angle θ = 250 • ).These locations are schematically depicted with purple dashed lines in panel (f ) of figure 5, where the solid dot, located approximately at one of the two boundaries of the current sheet, represents a single 'virtual probe' used as a reference.The choice of the fixed radius is intended to encompass the current sheet with the initial symmetrical velocity perturbation burst (panel g in figure 5), while the fixed poloidal angle is selected to include one of these velocity bursts.In this analysis, we specifically track the temporal evolution of the frequency corresponding to the SAW (ωτ A = 0.24 as shown in panel (d) of figure 3).The results of this analysis are presented in figure 6. Panel (a) illustrates the time evolution of the v θ field component as a function of the poloidal angle θ, depicting the two symmetrical burst perturbations, visible between t ∼ (1.65-1.75)× 10 4 τ A (note the nonlinear color scale used in this plot for better visualization of the v θ temporal evolution compared to the linear scale in figure 5).The solid purple line represents a single 'virtual probe' with coordinates (r, θ) = (0.4,110 • ).Panel (b) displays the evolution of the SAW intensity with respect to the poloidal angle θ.It is evident that the SAW is initially excited from two poloidal locations, θ ≃ 110 • and θ ≃ 250 • , during the formation of the current sheet (indicated by the blue arrow, which references the first time snapshot in figure 5).Subsequently, it is excited across the entire poloidal angle in ∼ 400 τ A , as the velocity perturbation covers the entire θ angle (time marked by the red arrow).Comparing this result with panel (a) indicates that the symmetrical velocity perturbation bursts at the boundaries of the current sheet serve as the starting points for the excitation of AWs.
We can further analyze how the AWs are excited in the radial direction from one of their starting points (θ = 250 • ), by examining the time evolution of the SAW frequency as a function of radius in panel (d) of figure 6.The velocity perturbation at θ = 250 • begins exciting the SAW at r/a = 0.4, corresponding to the q = 1 magnetic surface.After ∼ 400 τ A , the SAW is excited everywhere inside the q = 1 surface (r/a ⩽ 0.4), as the original plasma core is absorbed into the main 1,1 island, generating strong velocity perturbations within the entire reconnecting region.The excitation of the SAW beyond the q = 1 surface (r/a > 0.4) follows the outward expansion of the poloidal velocity perturbation over time, as depicted in panel (c), ceasing at r/a ∼ = 0.63.The dashed black curve that marks this expansion is defined by the edge of the flattening region in the q profile, that is where the q profile deviates from the initial one (flattens out) during the reconnection process.This region is also illustrated in the third panel of figure 1, where the dashed curve (q profile after reconnection) deviates from the solid line (initial equilibrium q profile) at r/a ∼ = 0.63.The outward expansion of the perturbation starts at r/a = 0.4 (q = 1 surface) and stops at r/a ∼ = 0.63 after the original core has been completely reconnected (for a more detailed illustration of the flattening of the q profile see supplementary material S2).
In conclusion, the aforementioned analysis reveals that the starting points (r, θ) for the excitation of AWs are located at (0.4, 110 • ) and (0.4, 250 • ).These starting points align with the velocity perturbations occurring at the boundaries of the current sheet, which develops as a consequence of the reconnecting process that characterizes the sawtooth dynamics.The same holds true for the more intense reconnection occurring ∼ 700 τ A later.Therefore, the observed Alfvénic fluctuations in our simulations can be interpreted as secondary modes excited by the magnetic reconnection events and the related velocity perturbations.

Additional analysis
Let us now discuss the role of the the visco-resistive parameters η and ν for the AWs excitation.We have obtained qualitatively similar results to those presented so far in ohmic tokamak simulations with S ≳ 10 7 and M ≳ 10 6 .To illustrate the impact of the visco-resistive parameters on the excitation of AWs during the sawtooth crash, we provide a counter-example in figure 7, where we present a simulation with the same setup as the main one discussed in section 2, but with a lower Lundquist number of S = 10 6 (corresponding to higher resistivity η) and a lower viscous Lundquist number of M = 10 4 (corresponding to higher viscosity ν).The magnetic and kinetic energies of the Fourier components are plotted in panels (a) and (b), respectively.The dynamics observed here are similar to those shown in the corresponding panels of figure 3, but with a significantly shorter time scale due to the higher resistivity, resulting in quicker exponential growth of the 1,1 mode (black curve) and two sawtooth crashes occurring in the same time frame as one in the previous simulation.An important distinction is the evolution of the 2,1 mode, which now remains stable and decays over time, leading to reduced coupling of the 1,1 mode with other harmonics (such as the 1,0 mode) as the simulation progresses.Consequently, the kinetic energy spike in panel (c) is approximately 20 times smaller than in the previous simulation case, and diminishes significantly with subsequent crashes.This much lower kinetic energy spike, coupled with the higher values of η and ν, results in a very weak broadband excitation and negligible Alfvénic activity, even in the core of the plasma (r/a = 0.3), where the SAW is not subject to the phase mixing as previously described.In panel (e), we provide a comparison between this simulation (solid lines) and the main simulation with high S (dashed lines).Both the peak amplitude and dissipation rate are greatly influenced by the choice of visco-resistive parameters.Furthermore, in this simulation, we do not observe the formation of plasmoids during the reconnection processes, in line with previously cited literature.This analysis highlights the necessity of a high Lundquist number (S ⩾ 10 7 ) for a meaningful self-consistent excitation of AWs following reconnection processes.
Lastly, let us discuss the appropriateness of the spectral resolution employed in our main simulation (depicted in red in figure 2).We present an analogous simulation with identical parameters but using a significantly broader range of modes (additional modes displayed in black in the same figure) beyond the h = 1 harmonics.In figure 8, we present the outcome of this simulation which starts from t = 1.5 × 10 4 τ A of the original one.For a straightforward comparison, we employ the same ranges and scales as in the main simulation displayed in figure 3. The temporal evolution of magnetic and kinetic energies exhibits remarkable similarity to the main simulation, with the 2,1 mode marginally stable, providing a coupling between the 1,0 mode and the 1,1 tearing mode.Both the maximum values and their temporal trends demonstrate close qualitative and quantitative agreement with the main simulation.The occurrence of the sawtooth crash remains unchanged, as expected due to the unaltered dissipation parameters and subsequent appearance of plasmoids.The total kinetic energy spike (panel (c)) is only slightly lower than in the previous case, and the excitation of Alfvénic fluctuations (panel (d)) and their dissipation rate (panel (e)) align with the main simulation.The excitation of the SAW is slightly reduced compared to before, while the GAE remains exactly the same (the weak frequency signal near the SAW, preceding the reconnection spike, arises from the stabilization of the additional Fourier modes added during the exponential growth phase).This comparison proves the adequacy of the chosen spectral resolution, as employing a broader resolution does not significantly alter the final results of AWs excitation or the dynamics of the sawtooth crash.

Summary and final remarks
In this paper, we have discussed a number of numerical simulation of the nonlinear resistive 1,1 kink mode using a pressureless viscoresistive MHD model for a cylindrical tokamak geometry.To include multiple helicities into the dynamics a marginally stable 2,1 mode was initialized too.Starting from an axisymmetric ohmic equilibrium with central safety factor below 1, the internal kink mode grows exponentially in time until it saturates and collapses due to magnetic reconnection events as in typical sawtooth dynamics.The latter induces current density and velocity perturbations at the q = 1 surface, which in turn excite Alfvénic fluctuations for S ≳ 10 7 and M ≳ 10 6 that persist beyond the relaxation event itself, being affected only by the relatively weak visco-resistive damping.An in-depth investigation shows a spatio-temporal correlation between the velocity perturbations associated with the reconnection and the excitation of the SAW in the core region and most notably a global mode, namely the GAE, with m = 1, n = 0 periodicity.In particular they are found to emanate from the outflow cones of the reconnection layer associated with the internal kink.The results obtained in this paper are also very robust with respect to the considered coupling approach.In fact, qualitatively analogous results are obtained for the 2,1 island amplitude ranging from 10% (the initialization value here considered) to 1% of the minor radius.Furthermore, qualitatively analogous results are also obtained if the marginally stable 3,2 mode is considered instead of the 2,1 for the three wave coupling, consistently with the mode coupling process discussed in [37,38].The excitation of the 1,0 GAE in our simulations is in line with a previous study carried out in [16] where excitation of the n = 0 GAE by IREs or ELMs was suggested.
The numerical results we have obtained, showing that AWs can be self-consistently excited by magnetic reconnection events in ohmic tokamak configuration characterized by internal kink/tearing dynamics, provides a possible mechanism to explain the Alfvénic activity that has been detected in the absence of superthermal ions during discharges in many conventional and spherical tokamaks, including EAST [5], TFTR [6], TUNAM-3M [10], COMPASS [11], MAST [12], Globus-M2 [13] and SUNIST [14], where oscillations at Alfvénic frequencies were observed following slow MHD events, all characterized by magnetic reconnection events.Furthermore, albeit our proposed mechanism may not be sufficient to explain the details of the observations in all of the referenced papers, it can provide a frame of reference for the interpretation of different regimes where the reconnection processes play an important role, including during the development of large magnetic islands where Alfvénic activity was observed in ohmic discharges in JET [4], FTU [7] and TEXTOR [8].On the other hand, the proposed mechanism is not limited to ohmic discharges and could, in principle, play a relevant role also in more fusion-relevant plasmas with additional heating sources (and possibly in the presence of fast particles) whenever a significant MHD activity is present.We also mention that an analogous mechanism was recently numerically observed in the reversed-field pinch configuration and found to provide a possible interpretation of experimental observations in the RFX-mod device [20].
The simulation results presented here are based on the viscoresistive MHD model in cylindrical geometry, and hence do not include toroidal, finite beta, two-fluids and kinetic effects.The inclusion of those effects will be addressed in the future and will enable us to study the proposed excitation mechanism in a more realistic modelling, which would take into account a wider variety of Alfvénic modes (e.g.toroidal or slow compressional Alfvén eigenmodes) that can be affected by such mechanism and would lead to more precise numerical predictions to be directly compared with experimental results.

Figure 1 .
Figure 1.Radial profiles corresponding to the initial axisymmetric equilibrium (solid lines) and to the configuration after the magnetic reconnection event (dotted lines) of the ohmic tokamak simulation.From left to right: magnetic field components, current density components, safety factor profile, density profile and the resulting analytical solutions of the low-frequency Alfvén waves spectrum with m = 1, n = 0 periodicity.

Figure 2 .
Figure 2. 3D Fourier modes of the tokamak simulations discussed in this article.We show in red the modes of the main simulation discussed in section 3 (indicating also the diagonal modes with helicity h = 1), and in black the wider set of modes considered for the additional analysis in section 4.2.

Figure 3 .
Figure 3. Temporal evolution of the fully 3D ohmic tokamak simulation (S = 10 8 , M = 10 6 ).Panels (a) and (b) show magnetic and kinetic energies of Fourier harmonics around the sawtooth relaxation event.The Fourier harmonics specifically discussed in the paper are shown with colors, while the other harmonics in the spectrum are plotted in gray.(c) Evolution of the total magnetic and kinetic energies, where ∆Em = Em − E f m with E f m = 386 627 the final magnetic energy at the end of the simulation.(d) The corresponding evolution of the frequency spectrum of the normalized velocity field component v 1,0 r , obtained with continuous wavelet transform (CWT), at a fixed radius of r/a = 0.3.The cone of influence (COI) is plotted in dashed black curve and is defined as the region of the wavelet spectrum (underneath the dashed black curve) affected by edge effects which may distort the resulting spectrum.The evolution of the central safety factor is shown in pink.(e) Amplitude of v 1,0 r at fixed frequencies ωτ A corresponding to observed Alfvén modes.

Figure 4 .
Figure 4. Numerical frequency spectra of the ohmic tokamak simulation (S = 10 8 , M = 10 6 ) after the magnetic reconnection event.The frequency spectra of v n=0 and B n=0 normalized fields components are shown as a function of radius at θ = 0 for the time window t = (29 500-30 000)τ A in figure 3(d).In the last panel the frequency spectrum is superimposed with the expected analytical solutions for the (m, n) = (1, 0) SAW and GAE.

Figure 5 .
Figure 5. Evolution of the internal kink mode in the cylindrical tokamak (S = 10 8 , M = 10 6 ) around the sawtooth relaxation event in figure 3. First row: axial component of the plasma current density together with the contour levels of the helical flux function with 1,1 helicity (black lines) are shown for different time snapshots.Plasmoid formation is observed at the tearing unstable current sheet on the q = 1 surface.Second row: poloidal component of the velocity field with the contour levels of the helical flux function.For the corresponding movie of the poloidal velocity component, see supplementary material S1.

Figure 6 .
Figure 6.Temporal evolution of the v θ component and the corresponding wave intensity at fixed frequency ωτ A = 0.24 , at either fixed radius r/a = 0.4 (a)-(b) or fixed poloidal angle θ = 250 • (c)-(d).The solid purple line correspond to the purple dot with coordinates (r, θ) = (0.4,110 • ) in panel (f ) in figure 5, while the arrows mark the snapshot times of the same figure.COI is plotted in black vertical dashed lines.The dashed black curve in panel (c) defines the edge of the flattening region of the q profile.

Figure 7 .
Figure 7. Temporal evolution of the fully 3D ohmic tokamak simulation (S = 10 6 , M = 10 4 ).Panels (a) and (b) show magnetic and kinetic energies of Fourier harmonics.(c) Total magnetic and kinetic energies, where ∆Em = Em − E f m with E f m = 38, 658.(d) Frequency spectrum of the normalized velocity field component v 1,0 r .COI is plotted in dashed black curve.The evolution of the central safety factor is shown in pink.(e) Amplitude of v 1,0 r at fixed frequencies ωτ A corresponding to observed Alfvén modes (solid lines) compared to the same plot (dashed lines) of the main simulation in panel (e) in figure 3.

Figure 8 .
Figure 8. Temporal evolution of the fully 3D ohmic tokamak simulation (S = 10 8 , M = 10 6 ) with higher spectral resolution.Panels (a) and (b) show magnetic and kinetic energies of Fourier harmonics.(c) Total magnetic and kinetic energies, where ∆Em = Em − E f m with E f m = 38 662.(d) Frequency spectrum of the normalized velocity field component v 1,0 r .COI is plotted in dashed black curve.The evolution of the central safety factor is shown in pink.(e) Amplitude of v 1,0 r at fixed frequencies ωτ A corresponding to observed Alfvén modes.