{'fraction_non_alphanumeric': 0.034392181605140154, 'fraction_numerical': 0.013983492782921082, 'mean_word_length': 4.820143884892087, 'pattern_counts': {'":': 0, '<': 9, '<?xml version=': 0, '>': 8, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Global gyrokinetic simulations of ASDEX Upgrade up to the transport time-scale with GENE-Tango

An accurate description of turbulence up to the transport time scale is essential for predicting core plasma profiles and enabling reliable calculations for designing advanced scenarios and future devices. Here, we exploit the gap separation between turbulence and transport time scales and couple the global gyrokinetic code GENE to the transport-solver Tango, including kinetic electrons, collisions, realistic geometries, toroidal rotation and electromagnetic effects for the first time. This approach overcomes gyrokinetic codes' limitations and enables high-fidelity profile calculations in experimentally relevant plasma conditions, significantly reducing the computational cost. We present numerical results of GENE-Tango for two ASDEX Upgrade discharges, one of which exhibits a pronounced peaking of the ion temperature profile not reproduced by TGLF-ASTRA. We show that GENE-Tango can correctly capture the ion temperature peaking observed in the experiment. By retaining different physical effects in the GENE simulations, e.g., collisions, toroidal rotation and electromagnetic effects, we demonstrate that the ion temperature profile's peaking is due to electromagnetic effects of submarginal MHD instability. Based on these results, the expected GENE-Tango speedup for the ITER standard scenario is larger than two orders of magnitude compared to a single gyrokinetic simulation up to the transport time scale, possibly making first-principles ITER simulations feasible on current computing resources.


I. INTRODUCTION
Reliable calculations of plasma profiles are essential to design scenarios with improved confinement and enable accurate performance predictions in present magnetic confinement devices and future devices. This is a challenging task due to the large variety of coexisting effects acting at different scales in time and space, affecting the evolution of plasma profiles. While particle and energy sources (either externally injected or arising from fusion reactions) act on macroscopic scales, microturbulence -driven by plasma instabilities destabilized by gradients in the pressure profile -leads to radial transport of particle and energy on the microscopic scales. A primary goal in fusion research is to develop tools able to describe these mutually interacting effects selfconsistently. The main obstacle in this regard is represented by the large timescale separation between microscopic and macroscopic physics in the core of magnetic confinement devices [1]. Assuming gyro-Bohm scaling of turbulent transport [2], this timescale gap increases as (1/ρ * ) 2 , with ρ * = ρ s /a being the ratio between the ion sound Larmor radius and the minor radius of the device. Therefore, for a device such as ITER (with 1/ρ * ∼ 1000) [3], the microscopic timescales are expected to be six orders of magnitude smaller than the macroscopic timescales at least (while the electron turbulence time scale is expected to be roughly 1µs, the energy confinement time a few seconds) [4]. This timescale gap has been a major limiting factor for high-fidelity gyrokinetic codes [5,6] to simulate the micro-turbulence dynamics up to the transport timescale. This is due to the prohibitive computational resources often required by these codes to simulate turbulence, even only at the microscopic timescales. Typically, ion or electron scale gyrokinetic simulations require at least ∼ 10 4 − 10 6 CPUh at a single radial position for covering only fractions of milliseconds [7]. This cost increases substantially (several order of magnitudes) for multi-scale simulations (see e.g., Ref. [8]).
In the past few decades, a major focus has been to develop accurate reduced turbulence models able to grasp the main signatures of plasma micro-instabilities and turbulence at a considerably reduced computational cost. This effort led to the emergence of a hierarchy of different codes mostly based on quasi-linear theory. Among the most well-known and widely used reduced turbulence codes, we note TGLF (trapped-gyro-Landau-fluid) [9,10] and QuaLiKiz [11][12][13][14]. While TGLF is a gyrofluid electromagnetic code retaining shaped tokamak geometry effects, QuaLiKiz is a reduced gyrokinetic code developed in the electrostatic limit and for shifted circular plasmas.
The great advantage of these models is the tremendous speedup (a few seconds on a single CPU) compared to gyrokinetic codes that enabled dynamic profile calculations. This is routinely done nowadays with integrated modelling tools, such as TGLF-ASTRA [15,16], TGLF-T-GYRO [17], QuaLiKiz-JINTRAC [18], by coupling these reduced turbulence models with transport codes. While the turbulence codes compute the transport coefficients for given plasma profiles, the transport codes evolve the profiles according to sources and turbulent fluxes. These tools are shown to accurately reproduce the experimental measurements in different devices and on a large variety of scenarios, contributing to performance predictions and scenario development.
However, the reduced turbulence models fail to fully describe plasma micro-instabilities and turbulence in strongly nonlinear regimes. In particular, in highly electromagnetic regimes and in cases with substantial external plasma heating, a growing number of studies show that the reduced turbulence codes under-predict turbulence suppression [7,[19][20][21]. This leads to a systematic mismatch of the on-axis temperatures in these regimes, thus affecting plasma performance predictions. This limitation is not expected to affect gyrokinetic codes that are shown to capture correctly nonlinear electromagnetic effects on the microscopic time scales [22][23][24][25][26]. These findings motivate (i) the further development of the reduced turbulence models and (ii) the improvement of the highfidelity models aimed at self-consistent profile evolution.
In recent years, an increasing effort has been spent on extending gyrokinetic codes and allowing simulations from the turbulent to the transport timescales. One of the most promising approaches developed so far is the so-called multiple-timescale approach. It employs the basic scheme used with reduced turbulence models and transport codes and consists of simulating turbulence and transport phenomena only on their natural timescales [27][28][29]. This is achieved e.g., by coupling a gyrokinetic turbulence code to a transport code. More precisely, the turbulence code evaluates the turbulence levels for a given pressure profile over several microscopic or fastscale characteristic times, while the transport code evaluates the new plasma profiles consistent with the given turbulence levels and the experimental sources. These new profiles are transferred back to the turbulence code, and the process is repeated until the volume-averaged energy and particle injected by the external sources is equal to the turbulent fluxes (steady-state solution). This code coupling leads to a reduction in computational expense of several orders of magnitude, making first-principles simulations over the confinement time -even for ITERlike devices -feasible on presently available computing resources.
The applicability of such a multiple-timescale approach has been already demonstrated when using gyrokinetics directly for the surface-averaged turbulent fluxes in local (flux-tube) frameworks with GENE [30], GS2 [31] and Trinity [32] and in simplified global setups (circular plasmas and adiabatic electrons) [33,34]. In particular, this global coupling has been achieved by coupling the radially global version of the gyrokinetic code GENE [35] with the transport-solver Tango. In this paper we further extend this new integrated modelling tool (GENE-Tango) by including kinetic electrons, collisions, realistic plasma geometries, electromagnetic effects and toroidal rotation (not yet evolved with Tango). A radially global approach is most likely to be required to investigate global transport features such as internal transport barrier formation [36][37][38], turbulence spreading [39,40], transport avalanches [41][42][43], and supra-thermal particle modes [44]. In particular, GENE-Tango is here applied to study two ASDEX Upgrade discharges. We first benchmark GENE-Tango against previously published results obtained with the radially local model GENE-Trinity [32] on the ASDEX Upgrade shot #13151 at t = 1.35s. Afterwards, we investigate the shot #31555 at t = 1.45s, which exhibits a pronounced peaking of the ion temperature profile that is not captured by TGLF-ASTRA. We show that the electromagnetic GENE-Tango simulations can correctly reproduce the experimental findings. To assess the physical effects responsible for the ion temperature peaking, we perform GENE-Tango simulations retaining different physical effects. This paper is organized as follows. Section II briefly introduces the underlying basic equations solved in the gyrokinetic code GENE and transport-solver Tango. The numerical scheme employed to couple these codes (socalled LoDestro method [45,46]) is described in Section III. While the numerical setup and grid resolutions are discussed in Section IV, the first results with kinetic electrons and collisions with GENE-Tango are presented in Section V together with a comparison with the profiles computed with ASTRA and GENE-Trinity (local approach). In Section VI, we apply GENE-Tango to study an ASDEX Upgrade H-mode plasma which exhibits a pronounced on-axis peaking not captured by TGLF-ASTRA, showing excellent agreement between the GENE-Tango profiles and the experimental measurements. In Section VII, we perform a detailed study of the impact of different physical effects on the evolution of the plasma profiles for this ASDEX Upgrade plasma, focusing in identifying the mechanism responsible for the ion temperature peaking. In particular, the role of collisions (Sections VII A and VII C), electromagnetic effects (Section VII D) and toroidal plasma rotation (Sections VII B and VII D) -induced by external neutral-beam-injection (NBI) -is analyzed. Section VIII contains a comparison of the overall speed-up obtained by running GENE-Tango compared to the macroscopic timescales. Finally, conclusions are drawn in Section IX.

II. THEORETICAL FRAMEWORK
The high-fidelity model used to study turbulence in magnetic confinement devices for many decades is gyrokinetics [5,6]. It can be derived rigorously starting from the coupled system of Fokker-Planck and Maxwell's equations under the so-called gyrokinetic ordering [47]. This ordering imposes several constraints -motivated by experimental observations -on the amplitudes and spacetime scales of the macro-and micro-physics. In particular, it assumes that the fluctuations are (i) highly anisotropic with parallel correlation length largely exceeding the perpendicular one; (ii) of small amplitude with respect to the background quantities (e.g., magnetic-surface-averaged profiles); and (iii) of low frequency compared to the gyro-frequency with background quantity evolution being significantly slower than the fluctuations. With these assumptions, the particle gyromotion reduces in the underlying equations by a gyroring description. It is therefore particularly convenient to solve the gyrokinetic equations in the so-called guidingcenter coordinate system. This leads to a reduction of the dimensions of the distribution function from six to five and the elimination of irrelevant space-time scales with a significant saving in computational time without neglecting any meaningful physical effects. Many gyrokinetic codes, including the gyrokinetic code GENE, exploit the first ordering assumption -particularly valid in the core of magnetic confinement devices -and further reduce the equation complexity by employing the so-called δf ordering, namely by splitting the distribution function into a background and a small perturbation Under these assumptions, the system evolves on the microscopic scale as follows Here, Γ s = ∇δf 1,s − q s ∇φ 1 /m s v (∂f 0 /∂v ); and C [δf 1,s ] represents the two-particle Coulomb operator, v c the magnetic drift velocity (containing the curvature, the grad-B and the E × B drifts), B 0 the toroidal magnetic field, b 0 the unit vector along the magnetic field, c the speed of light, µ = m s v 2 ⊥ /(2B 0 ) the magnetic moment, q s the charge, T 0,s the temperature and m s the mass of the species s. The overbar on the electrostatic (φ 1 ) and electromagnetic (A 1, ) fields denotes gyroaveraged quantities. The contribution of the parallel magnetic-field fluctuations is, for the moment, neglected in the global version of the code GENE. We refer the reader to Refs. [35] for further details. As mentioned above, this equation is coupled to the Maxwell's equations, which in the particle coordinate system read Here, δf pc 1,s represents the perturbed part of the distribution function f in the particle coordinate system. The solution of Eqs. (1), (2), (3) requires the additional field equation to evaluate the time derivative of A 1, This is solved by replacing the time derivative of the perturbed distribution function with the right-hand-side of Eq. (1) written in particle coordinates. The transformation between the gyrocenter (removed gyro-angle dependence from guiding center coordinates) and particle variable on the perturbed quantities is performed via the pullback operator. For a detailed description and derivation of the gyrokinetic equations we refer the reader to Ref. [5]. This set of equations is particularly suited to study micro-turbulence dynamics at the microscopic timescales. By solving Eqs. (1), (2), (3), (4) we evolve the perturbed distribution function δf 1 on the microturbulence time scale and compute particle and energy losses due to plasma micro-instabilities. These are determined by energy and particles fluxes that are respectively defined for the generic species s and Here, m s represents the mass of the species s, while v E×B the E × B velocity and · S the surface average. The background quantities are considered constant on the micro-turbulence time scale since their evolution occurs on time scales that are orders of magnitude slower than the micro-turbulence scales. However, to aim for a self-consistent description of experimental discharges, it is essential to simulate not only the microscopic scales (a few milliseconds in ITER) but the macroscopic scales as well (a few seconds in ITER). This is particularly important since particle and energy losses will inevitably affect plasma profiles which, in turn, will modify the turbulence drive. In this regard, single gyrokinetic simulations covering these long time scales (at least a confinement time) are particularly rare in literature due to the prohibitive computational resources required. Most of the existing work has been carried out on reduced setups (adiabatic electrons, collisionless and/or electrostatic plasma, simplified geometry, no impurities nor supra-thermal particles) [48][49][50][51]. In recent years, an increasing effort has been spent on extending gyrokinetic codes improving code performances, e.g., via GPU porting [52]; and developing new algorithms exploiting the time separation between micro-and macro-physics, e.g., time-telescoping methods [32,33,53]. One of the most promising approaches developed so far allowing gyrokinetic simulations up to the transport time-scale, consists of coupling the gyrokinetic code to a 1D transport code (multipletimescale approach). In this manuscript, this is done by coupling the global version of the gyrokinetic code GENE to the transport code Tango [33,34]. More precisely, GENE evaluates the turbulence levels for a given pressure profile over several microscopic time steps, while Tango evaluates the new plasma profiles consistent with the given turbulence levels and the experimental sources. For kinetic-electron-setups, the transport code solves the 1D equations for density and pressure (temperature) for each plasma species 3 2 with V ′ = dV /dx differential volume of the flux surfaces, Q s and Γ s , respectively, heat and particle fluxes defined in Eqs. (5), (6), S n and S i,e external particle and heat sources and ν e,i collisional energy exchange frequency [54] defined as Here, Λ e,i = 10, the density is expressed in m −3 , the temperature in eV and the mass in units of proton mass. The collisional energy exchange ν e,i couples Eq. (8) for the ion and electron pressure. Once the transport code evolves the plasma density and temperatures at the next transport time step, these new profiles are transferred back to the turbulence code, which computes the new turbulence levels consistent with the new profiles. This iterative process is repeated until the volume-averaged energy and particle injected by the external sources is equal to the turbulent fluxes. The multiple time-step approach leads to a significant reduction in computational expense. Interestingly, this approach becomes more and more efficient (with more favorable computational saving) as the time scale separation between turbulence and transport time scales grows.
For ITER-like devices, it is expected to speed-up convergence to the steady-state solution by several orders of magnitude, making first-principles simulations over the confinement time feasible on presently available computing resources.
It is worth mentioning that momentum transport and variations to the plasma geometry are not considered in the current version of Tango and neglected for the analyses described within this manuscript. These capabilities will be included in the near future. Moreover, for simplicity, the heat and particle fluxes in Eqs. (7), (8) consist only of the turbulent contribution, namely Eqs. (5), (6).

III. SOLUTION OF THE TRANSPORT EQUATIONS
In this Section, we describe the numerical scheme used to couple the global gyrokinetic code GENE and the transport-solver Tango. This coupling is based on the co-called LoDestro method [45,46]. It was first introduced in Ref. [45] to evolve the plasma density by coupling the two-dimensional Hasegawa-Wakatani equations to Eq. (8). Afterwards, the Hasegawa-Wakatani model was replaced by gyrokinetic simulations in reduced setups involving adiabatic electrons, circular plasmas and simplified shapes for the external sources [33,34].
The LoDestro method prescribes a procedure to solve the system of equations 7, 8 within an implicit timestep advance. In particular, Eqs. (7), (8) are discretized in time using backward differences and solved -for each macroscopic timestep -iteratively. The set of Eqs. (7), (8) is hence written Here, the subscript m denotes the macroscopic time index and the subscript l the iteration index within a timestep ∆t. The species subscript has been removed for simplicity. Moreover, the sum over the index u indicates a sum over all the plasma species s. It is worth mentioning that in this work, although Eqs. (10), (11) are fully time-dependent, we solve them only for its steady-state solution taking a single large timestep. As mentioned above, the ion and electron pressure equations are coupled via the collisional energy exchange. Therefore, the ion and electron pressures are updated together implicitly to achieve stable iteration. The coupled set of Eqs. (10), (11) are solved in the LoDestro method by setting ad-hoc rules for the turbulent fluxes that depend nonlinearly on the plasma pressure and its logarithmic gradient. In particular, particle and heat fluxes are decomposed in a convective and diffusive contribution, which are assumed to fulfill a local dependence with the pressure and its gradient, namely Γ m,l = −D n m,l−1 ∂ x n m,l + c n m,l−1 n m,l .
Here, we have introduced, respectively, the diffusive and convective transport coefficients D m and c m . They are computed at the iteration l − 1 and employed to decompose the turbulent fluxes at the iteration l.
The splitting of the turbulent fluxes in its convective and diffusive contributions is performed in Tango with a user-specified parameter called θ m . We have and Given the mostly diffusive nature of turbulent transport in the core of magnetic confinement devices, we have fixed this parameter to θ m = 0.95 when D m > 0 (c m > 0) and to θ m = −1 for D m < 0 (c m < 0). We note that the convective and diffusive coefficients D and c should not be interpreted as physical diffusive and pinch terms. At convergence, Eqs. (12) -(15) provide the turbulent fluxes. With this splitting of the turbulent fluxes, the discretized transport equations 10, 11 can be written for the generic species s as (16) n m,l − n m−1 ∆t Therefore, the plasma density and pressure at the iteration l depend only on known quantities, e.g., transport coefficients computed at the iteration l − 1, external sources and plasma pressure and density (and their gradients) evaluated at the macroscopic timestep m − 1. For each timestep m, the iterations proceed until the turbulent fluxes -computed with the high-fidelity gyrokinetic code GENE -match the volume-average of the external particle and heat sources. As already discussed in Ref. [46], the high sensitivity of turbulent transport on changes in the plasma pressure gradients might lead to stability issues within the iteration loop. This is avoided in Tango by applying relaxation on the fluxes (and hence transport coefficients) and plasma pressure. In particular, we replace the turbulent fluxes with the "relaxed" fluxes, namely and This relation procedure acts as an average, and the effective turbulent fluxes seen by the transport code Tango at the iteration l − 1 correspond to the fluxes computed by the turbulence code GENE at the iteration l−2 (with the profiles p l−2 ) and a minor correction that is determined by the fluxes at the present iteration (fixed by the amplitude of α). Similarly, the pressure profiles passed by the transport code Tango to the turbulence code GENE at the iteration l − 1, are evaluated as the average of the pressure profiles at the iteration l − 2 and the one at the present iteration. The magnitude of the relaxation coefficients depends on the sensitivity of the turbulent transport on the gradients of the plasma pressure and needs to be properly tuned case by case. The optimal values found for the present studies range from α p = [0.05 − 0.2] and α q = [0.1 − 0.4].

IV. NUMERICAL SETUP
In the following, we summarize the main parameters employed throughout this work. The GENE simulations are performed with deuterium ions and kinetic electrons with realistic ion-to-electron mass ratio. Collisions are retained (except for the analyses of Section VII C) and modelled with a linearized Landau-Boltzmann collision operator with energy and momentum conserving terms [55]. To account for the effect of impurities -neglected as an active species in this study -we include an effective charge Z eff . Neoclassical effects are typically small compared to turbulence and are neglected in the turbulent calculations [56,57]. As mentioned previously, the current version of GENE-Tango does not allow a self-consistent evolution of the magnetic geometry as the plasma pressure profiles evolve within the iterations. Therefore, the magnetic equilibrium is fixed to the one reconstructed via CLISTE for each ASDEX Upgrade discharge at the time-slice of interest.
Different grid sizes and resolutions have been employed in Sections V and VI. While a radial domain of ρ tor = [0.1 − 0.8] is used for the analyses of the ASDEX Upgrade discharge #13151, this is reduced to ρ tor = [0.05 − 0.7] for the discharge #31555. Here, ρ tor is the radial coordinate based on the toroidal flux Φ, ρ tor = Φ/Φ LCF S . The grid resolution employed in the radial (x), bi-normal (y) and field-aligned (z) directions is respectively (nx0 × nky0 × nz0) = (256 × 48 × 32) for the ASDEX Upgrade discharge #13151 and (nx0 × nky0 × nz0) = (225 × 36 × 32) for the discharge #31555. The discretized toroidal mode number is given by n = n 0,min · j with j being integer-valued in the range j = [0, 1, 2, ..., nky0 − 1] and n 0,min = 2. The velocity space grids of species (s) are set from −3.5 to 3.5 for v /v th,s and from 0 to 12 for the magnetic moment µB 0 /T s . Here, v th,s = (2T s /m s ) 1/2 represents the thermal velocity and B 0 the on-axis magnetic field. We employed 32 points along the parallel velocity and 24 along the magnetic moment.
The numerical GENE simulations have been performed running GENE in gradient-driven mode. Therefore, Krook-type particle and heat operators have been applied to keep the plasma profiles (on average) fixed to the ones provided by Tango at each iteration. However, as GENE-Tango approaches the steady-state solution, the amplitude of the Krook heat (γ k ) and particle (γ p ) coefficients can be significantly reduced. In particular, we used, respectively, γ k = 0.01c s /a and γ p = 0.01c s /a. Here, c s = (T e /m i ) 1/2 represents the sound speed, with T e the electron temperature at the reference radial position and m i the bulk ion mass in proton units.
Numerical fourth-order hyperdiffusion is used to damp fluctuations at large toroidal mode numbers due to electron temperature gradient (ETG) modes to reduce the otherwise prohibitive computational cost of these simulations [58]. Moreover, we used buffer regions covering 10% of the GENE simulation radial domain to damp fluctuations to zero near the domain boundaries and enforce consistency with the Dirichlet boundary conditions. In these areas we applied a Krook operator γ b with an amplitude γ b = 1.0c s /a.
In Tango, different approaches are used to treat the turbulent fluxes and sources in these regions to account for the unphysical damping on plasma fluctuations. In the inner buffer region, the physical sources are set to zero (see e.g., Fig. 1, 5). This is typically a minor correction to the volume integral of the injected sources due to the reduced plasma volume close to the magnetic axis. On the other hand, at the outer buffer, we introduce extrapolation regions in Tango (which might differ for temperatures and density). Specified radial domains are selected in areas unaffected by the Krook operator. Here, Tango performs an interpolation of the GENE turbulent fluxes. The turbulent fluxes outside this region (towards the right boundary) are replaced by a linear extrapolation until the end of the GENE-Tango grid. The specific choice of the extrapolation regions has been adjusted accordingly during the GENE-Tango iterations.
Moreover, Tango applies Neumann boundary conditions for the inner boundary and Dirichlet for the outer one.

V. GENE-TANGO BENCHMARK
We begin by applying GENE-Tango to study the profile evolution due to external auxiliary power and microturbulence for the ASDEX Upgrade discharge 13151 at t = 1.35s [59]. It is an H-mode plasma heated with 5MW of neutral-beam-injection (NBI). Although a dated experiment, this discharge was selected since it was previously studied with flux-tube simulations with GENE-Trinity [32], thus providing an interesting test case to compare different (radially local and global) transport and turbulent models. For the comparison, the same heating, particle source and geometry used in GENE-Trinity are employed in this section. The radial profile of the power absorbed by the thermal ions and electrons and particle refueling are extracted from the ITER database [60] and shown in Fig. 1. In this regard, the ohmic heating, the total radiated power and the particle source profile from SOL neutrals were not accessible and hence neglected in the simulations. We note that the amplitude of the external sources has been set to zero in the radial domain covered by the inner buffer regions in GENE to avoid nonphysical steepening of the pressure profile in this area. In particular, due to the effect of the Krook operators in GENE the turbulent fluxes will be artificially damped in this region, forcing Tango to increase the pressure gradients. The magnetic geometry is reconstructed by TRACER and read into GENE via numerical field-line tracing provided by the TRACER-EFIT interface [61]. Impurities have been neglected for this case and the effective charge is set to Z eff = 1. The numerical setup and physical parameters are the same as the ones summarized in Section II. Before starting the coupled GENE-Tango iterations, a GENE standalone simulation is performed to evaluate the turbulence levels corresponding to the initial ion and electron temperatures and density. These profiles are fixed to the ASTRA ones at t = 1.35s (see Fig. 3). This choice sets a constraint over the temperatures and density at the right boundary ρ tor = 0.8. The resulting time-averaged (over the steady-state nonlinear phase) heat and particle fluxes are shown in Fig. 2.
With the initial profiles computed by ASTRA, GENE predicts large turbulent fluxes that are not consistent with the volume integral of the injected heating and particle sources. Excellent agreement between turbulence fluxes and volume integral of injected sources is achieved only after running GENE-Tango for 80 iterations, corresponding to an overall GENE run-time of 41ms. This is illustrated in Fig. 2 where the GENE turbulent heat and particle fluxes -averaged over the last five GENE-Tango iterations -are shown. To speed up the convergence, the relaxation coefficient α p and the GENE run-time per iteration have been initially fixed to α p = 0.1 and t = 250c s /a and then progressively reduced to α p = 0.05 and t = 150c s /a as the GENE-Tango iterations were approaching the final solution. This reduction of α p is required to make the GENE-Tango iterations more stable, as the system is brought closer to marginal stability. The relaxation coefficient acting on the turbulent fluxes has been kept constant to the value α q = 0.4. Looking at steady-state profiles of Fig. 3, we observe excellent agreement between GENE-Tango and the local GENE-Trinity (flux-tube) model for the ion and electron temperature, with a relative error among the profiles of ∼ 5%. The largest deviations are observed on the electron temperature at ρ tor ≈ 0.2. These differences in the profiles might be caused by finite-size effects, with 130 [62,63], which get more pronounced as ρ tor is reduced. These essentially global effects typically lead to an over-prediction of the turbulence levels in the flux-tube GENE simulations when compared to the radially global ones [62,63], thus causing a more pronounced flattening of the bulk profiles to match the fixed injected power. For a medium-size device as ASDEX Upgrade the flux-tube calculation can lead to an over-prediction of the turbulent fluxes up to a factor of two compared to the global code [64]. The plasma density also shows also qualitatively good agreement between GENE-Tango and GENE-Trinity, although larger differences (∼ 7%) are observed compared to the temperature profiles.
Although the GENE-Tango profiles are largely consistent with the GENE-Trinity ones, more pronounced deviations are found with respect to the ion temperature and density computed by ASTRA. In Ref. [32] the lack of flow shear in the GENE-Trinity simulations and the oversimplified magnetic equilibrium used in ASTRA (possibly MHD unstable) were proposed as possible causes of discrepancy. On the other hand, the three transport codes' electron temperature profiles are in close agreement, apart from a deviation by Trinity at small ρ tor . For completeness we add the comparison of the logarithmic temperature and density gradients in Fig. 4.
As noted above, these results reveal that the plasma profiles for this ASDEX Upgrade discharge are particularly stiff, characteristic of turbulence near marginal stability. Relatively small differences in the ion temperature and density logarithmic gradients (Fig. 4) are found to impact the turbulent fluxes considerably, as shown in Fig. 2.

VI. ION TEMPERATURE PEAKING CAPTURED WITH GENE-TANGO AT ASDEX UPGRADE
We apply GENE-Tango to analyze a more recent discharge at ASDEX Upgrade, namely the shot #31555 at t = 1.45s (previously studied in Ref. [21] with TGLF-ASTRA). This experiment exhibits a pronounced peaking of the ion temperature profile not reproduced by the reduced turbulence model TGLF-ASTRA (see, e.g., Fig. 8). This is a known limitation of reduced turbulence models, typically under-estimating the stabilizing effect of electromagnetic effects and supra-thermal particles on ion-scale turbulence, leading to rather flat profiles not consistent with the experimental measurements [7,[19][20][21].
The ASDEX Upgrade shot #31555 at t = 1.45s is an H-mode plasma with 5MW of NBI heating with plasma current I p = 0.6MA, magnetic field B t = 2.8T and relatively low plasma fueling rate Γ n = 0.7 × 10 22 e/s. The power absorbed by the thermal species and particle fueling computed by ASTRA are illustrated in Fig. 5. While the electron heating source is defined as the sum of the ohmic, NBI heating, exchange electrons-ions energy term and the radiated power; the ion heating source as the sum of NBI heating and exchange electrons-ions energy term. The particle source retains the contribution coming from both the NBI source and the SOL neutrals. We note that the amplitude of the external sources has been set to zero in the radial domain covered by the inner buffer regions in GENE to avoid nonphysical steepening of the pressure profile in this area. In particular, due to the effect of the Krook operators in GENE the turbulent fluxes will be artificially damped in this region, forcing Tango to increase the pressure gradients.
On this specific discharge, we apply GENE-Tango simultaneously retaining electromagnetic effects, collisions and toroidal external rotation in the GENE simulations. The numerical setup and physical parameters are the same as the ones summarized in Section II.
As mentioned previously, the current version of GENE-Tango does not evolve the toroidal angular momentum. Therefore, we employed a fixed v tor profile throughout the GENE-Tango iterations. The profile of the toroidal rotation is taken from the measurements of the Charge Exchange diagnostic [65], fitted with a third order polynomial and it is shown in Fig. 6. Positive values of v tor denote a counter-clockwise rotation.
In the global version of GENE, a phase factor is applied on the perturbed part of the distribution function of each species at every time-step to mimic the effect of finite toroidal rotation, i.e., is defined as follows with C y = x 0 /q(x 0 ), x 0 = 0.375 center of the radial box, and a minor radius of the device. We note that (global) GENE does not include -at the moment -the parallel flow shear term. Therefore, no model for the parallel velocity gradient (PVG) instability is retained in our simulations. This might possibly lead to an over-estimation of the impact of a finite toroidal rotation on the turbulent fluxes [66]. Before starting the Tango iterations, we performed a GENE stand-alone global simulation with the ion and electron temperature and density profiles computed by TGLF-ASTRA (see Fig. 8). The TGLF simulations have been performed using the SAT1geo saturation rule [10,67]. The time-averaged fluxes are shown in Fig. 7. Although the GENE standalone results are close to the volume integral of the injected sources for ρ tor ≈ [0.3 − 0.5], an excellent agreement is obtained after running GENE-Tango for 29 iterations, each covering t = 150c s /a. The overall GENE run time is 11.7ms. The Tango relaxation coefficients acting on the plasma profiles and turbulent fluxes are fixed, respectively, to the α p = 0.15 and α q = 0.4.
The steady-state plasma profiles -averaged over the last five GENE-Tango iterations -are compared in Fig. 8 with those obtained with TGLF-ASTRA and the experimental measurements. While the electron temperature was reconstructed by combining two Thomson scattering systems (respectively one in the core while the other in the plasma edge) [68] and the electron cyclotron emission (ECE) system [69], the plasma (electron) density with the lithium beam emission spectroscopy diagnostic [70,71] and the deuterium-cyanide-nitrogen laser interferometer [72]. The charge exchange recombination spectroscopy systems was used for the ion temperature. To keep the separatrix electron temperature at ∼ 100eV (typical value at ASDEX Upgrade) [73], the electron temperature and density were shifted together. The dotted continuous red line in Fig. 12 represents a fit over the experimental data. The fit was performed using a third  order polynomial for the temperatures, and the two-line fit [74] for the density. Further details can be found in Ref. [21]. To better quantify the differences among the plasma profiles and experiment, we add in Fig. 9 the logarithmic temperature and density profiles. A first striking observation from Fig. 8 is the significant peaking of the ion temperature profiles -with respect to the TGLF-ASTRA profiles that remains rather flat -at ρ tor = 0.3. This is seen more clearly when looking at the logarithmic ion temperature gradient, that exhibits a localized increase for ρ tor = [0.2 − 0.3] which is consistent with the experimental measurements. Nevertheless, the steady-state ion temperature profile computed with GENE-Tango is still under-predicted for ρ tor < 0.2. This is likely to be caused by the absence of supra-thermal particles in the GENE simulations. Energetic particles are known to stabilize -when properly optimized -ITG turbulence strongly, thus possibly leading to a further peaking of the ion temperature on-axis [23,25,26]. The radial region ρ tor = [0 − 0.3] is the fraction of the plasma with the largest supra-thermal particle concentration in this ASDEX Upgrade discharge, due to the on-axis NBI heating scheme used in the experiments. The impact of supra-thermal particles on these results will be analyzes in the near future. This required some non trivial changes in the interface between GENE and Tango and will be performed in the near future. Moreover, we observe from Fig. 8 an excellent agreement on the electron temperature between the different models and the experimental measurements. The only minor deviations are observed close to the GENE left buffer region, possibly related to the procedure used in Tango to account for the GENE inner buffer (namely by setting the amplitude of the sources to zero in the inner buffer region).
Despite the excellent agreement obtained on the ion and electron temperature profiles for ρ tor = [0.2−0.7], we notice that GENE-Tango systematically under-predicts the plasma density for this ASDEX Upgrade discharge with respect to the experimental measurements. This might be either related to missing physical effects in the GENE-Tango calculations or to an under-prediction of the plasma density at the GENE-Tango right boundary condition due to the large experimental uncertainties.
A. Improved agreement on the plasma density with boundary density value fixed by IDA In the previous section, we demonstrate that GENE-Tango can reproduce quantitatively the experimental temperature profiles of the ASDEX Upgrade discharge #31555 at t = 1.45s, showing a pronounced peaking of the ion temperature not captured by TGLF-ASTRA. However, the density profile computed by GENE-Tango did not recover the experimental measurements, suggesting that the value of the plasma density at the GENE-Tango outer boundary was not consistent with the GENE turbulent levels. The choice of the plasma density (and temperatures) at the GENE-Tango outer boundary strongly influences the whole density profile at the inner core region. By looking at Fig. 8 we ob-serve large uncertainties on the experimental measurements of the plasma density especially at ρ tor = 0.7, thus making an accurate evaluation of the plasma density at the boundary particularly challenging. To assess whether an improved agreement with the experimental density profile could be achieved with a different choice of the boundary condition, we performed GENE-Tango simulations with a different value of the plasma density at ρ tor = 0.7. The ion and electron temperature values at the GENE-Tango outer boundary are unchanged. The value of the density at ρ tor = 0.7 has been fixed to that of the IDA [75] density profile (red line in Fig. 10). It differs with respect to those employed by TGLF-ASTRA by ≈ 8% (compatible with the experimental error bars), which was obtained with the integrated model based on engineering parameters (IMEP) introduced in Ref. [76].
To reduce the computational cost of the simulation, we initialized the GENE-Tango iterations by running a GENE stand-alone simulation with the steady-state temperature profiles obtained with the TGLF-ASTRA By looking at Fig. 10 we observe not only an excellent agreement on the ion and electron temperature profiles, but also an improved agreement for the plasma density, which well sits inside the experimental measurements in ρ tor = [0.3 − 0.7]. However, we still notice an underprediction of the density profile in ρ tor = [0 − 0.2]. As mentioned above, this result might be due to the absence of supra-thermal particles in the GENE simulations, which are known to strongly suppress heat and particle transport (see e.g., [78]), thus leading to increased core gradients. The impact of fast particles on the evolution of the plasma profiles computed by GENE-Tango will be addressed in the near future.
It is worth mentioning here that the GENE-Tango sim-ulations presented in this and the previous Section lasted, respectively, 24h and 72h on 16 IBM POWER9 AC922 nodes, each with 4 Nvidia Volta V100 GPUs on Mar-coni100. Although the computational cost of the GENE-Tango simulations is greater than reduced models like TGLF-ASTRA (a few hours on a single CPU [76]), it still allows the routine use of GENE-Tango for discharge analyses. A detailed discussion on the GENE-Tango speedup and extrapolations to ITER can be found in Section VIII.

VII. IMPACT OF DIFFERENT PHYSICAL EFFECTS ON PLASMA PRESSURE PROFILES
In the previous Section, an excellent agreement between the GENE-Tango profiles and the experimental measurements was found even in regimes where TGLF-ASTRA under-predicts the peaking of the on-axis ion temperature profile. Given the moderate computa-tional resources required to run GENE-Tango up to the transport time-scale, in this Section we perform several GENE-Tango simulations that retain different physical effects in the turbulence calculations. These analyses are essential to identify the mechanism responsible for the increased on-axis ion temperature and, possibly, provide insights into improving reduced turbulence models.
Therefore, the impact of collisions, toroidal plasma rotation and electromagnetic effects on the evolution of the plasma profiles is studied with GENE-Tango simulations of the ASDEX Upgrade discharge #31555 at t = 1.45s. To compare more easily the GENE-Tango profiles with TGLF-ASTRA, we employ the same boundary values on temperatures and density used in TGLF-ASTRA. However, we note -as shown in Section VI A -that GENE-Tango could achieve an improved agreement with the experiment when fixing the plasma density at ρ tor = 0.7 to the value obtained with IDA. The extrapolation regions used in Tango

A. Electrostatic GENE-Tango simulation with collisions and without external E × B shear
We begin these analyses in the electrostatic limit (β e (x 0 ) = 8πn(x 0 )T e (x 0 )/B 2 0 = 1e − 4 with x 0 = 0.375) with collisions modeled with a linearized Landau-Boltzmann operator. The numerical setup and physical parameters are the same as the ones summarized in Section II. We initialize the iterations with Tango by running a GENE standalone simulation using the TGLF-ASTRA profiles. The results are illustrated in Fig. 11. Large turbulence fluxes are observed in all the different channels, which are not compatible with the volume average of the injected sources, thus suggesting that the TGLF-ASTRA profiles cannot be sustained in the GENE-Tango electrostatic simulations with collisions. Starting from these reference results, an excellent agreement between the GENE fluxes and the experimental power balance is achieved after 20 GENE-Tango iterations, corresponding to an overall GENE run-time of 8.1ms. This is shown in Fig. 11 where the heat and particle fluxes obtained over the last 5 iterations are compared with the volume integral of the injected sources. For these analyses, the Tango relaxation coefficients acting on the plasma profiles and turbulent fluxes are kept constant, respectively, to the values α p = 0.2 and α q = 0.4. The GENE run time (for each iteration) was set to 150c s /a. The temperatures and density computed by Tango over the last 5 GENE-Tango iterations are compared in Fig. 12 with the TGLF-ASTRA profiles and the experimental measurements.
These results show a good agreement between the elec-tron temperature profiles. However, GENE-Tango predicts a more pronounced flattening of the ion temperature for ρ tor < 0.4 compared to TGLF-ASTRA, thus making T i further deviating with respect to the experimental measurements. While the electron temperatures follow the experimental measurements closely, we notice that both the electrostatic GENE-Tango ion temperature and the TGLF-ASTRA one do not exhibit the on-axis peaking observed in the experiment. Moreover, Fig. 12 shows that the density profile computed with GENE-Tango deviates respect to the TGLF-ASTRA profile and the experimental measurements. In particular, GENE-Tango predicts a rather flat profile not in agreement with the experiment and TGLF-ASTRA. This result is consistent with Fig. 8 and provide further evidence indicating that the value of the plasma density on the right boundary of the GENE-Tango simulation (fixed to the TGLF-ASTRA one) is not fully consistent with the transport levels computed by GENE. The impact of the plasma density at the GENE-Tango outer boundary was discussed in detail in Section VI A.
To better quantify variations among the different profiles, we show in Fig. 13 the logarithmic temperature and density gradients. The logarithmic gradient of the experimental data has been computed on the fitted line. These analyses confirm the good agreement between GENE-Tango and TGLF-ASTRA on the electron temperature. More pronounced variations are observed for ω Ti and ω n with reduced gradients predicted by Tango in most of the radial domain considered.

B. Electrostatic GENE-Tango simulation with collisions and external E × B shear
The under-prediction of the on-axis temperature with TGLF-ASTRA is often attributed to electromagnetic and supra-thermal particle effects on plasma turbulence that are not fully captured yet by reduced turbulence models, such as TGLF [7,[19][20][21]. However, these codes have been shown recently to qualitatively recover the experimental peaking on T i by artificially increasing the amplitude of the toroidal rotation (called v tor in the remaining of this Section) above the experimental uncertainties [20]. Nevertheless, high-fidelity gyrokinetic codes exhibit only a weak dependence in the plasma core on the turbulent fluxes with the toroidal rotation (for experimentally relevant amplitudes of the flow shear), despite the strong effect observed on reduced turbulence codes [23,79,80]. The role of the toroidal rotation on the evolution of the plasma profiles computed with GENE-Tango will be addressed in detail in the present Section and in Section VII D.
Toroidal rotation in magnetic confinement devices leads to a radial electric field and hence to a finite E × B plasma flow, typically reducing the turbulence correlation lengths and outward particle and energy fluxes [81,82]. An efficient way of inducing these large-scale plasma flows in experiments is via neutral beam injection. Given the non-negligible NBI heating power used in the ASDEX Upgrade discharge #31555 at t = 1.45s, we investigate here the role of the E × B shear on the T i peaking observed in the experiment with gyrokinetic simulations up to the transport time scale with GENE-Tango. The profile of the toroidal rotation (not evolved in Tango) used in this Section is the same as the one shown in Fig. 6.
To begin the GENE-Tango iterations, we run at first a GENE standalone simulation to convergence using the TGLF-ASTRA temperature and density profiles. The time-averaged turbulent fluxes obtained with GENE are illustrated in Fig. 14 (red line) and compared with the initial GENE simulation without external toroidal rotation (cyan line) (Fig. 11). Fig. 14 shows (with respect to the electrostatic GENE-Tango results without external toroidal rotation) that the inclusion of the external toroidal rotation leads to a turbulence stabilization in all the different channels. While only a mild reduction in heat and particle fluxes is observed for ρ tor < 0.3 (roughly ∼ 15%), we observe a progressive enhancement in the outer regions (roughly 42% at ρ tor = 0.6). This turbulence stabilization brings the initial GENE turbulent fluxes closer to the volume integral of the injected sources (black dots of Fig. 14). However, the dependence of the ion heat flux with v tor is not consistent with the experimental measurements of T i . In particular, the on-axis temperature profile peaking observed in the experiment would require a localized stabilization in the deep core regions. This is in contrast with the effect of v tor on the GENE fluxes observed in Fig. 14 (localized more in the plasma edge).
Starting from the initial fluxes of Fig. 14, we perform 15 GENE-Tango iterations before reaching the steadystate solution where the GENE fluxes match the power balance in all turbulent channels. For each iteration, GENE has been run up to t = 150c s /a, thus leading to an overall GENE run-time of 6.1ms. The Tango relaxation coefficients acting on the plasma profiles and turbulent fluxes are kept constant, respectively, to the value α p = 0.2 and α q = 0.4. The converged temperature and density profiles computed by Tango -averaged over the last five GENE-Tango iterations -are illustrated in Fig. 15 and compared with the TGLF-ASTRA ones, the electrostatic GENE-Tango profiles without external toroidal rotation and the experimental measurements. We also add the logarithmic temperature and density profiles in Fig. 16 to better quantify the differences among the models and the experiment. Fig. 15 shows an excellent agreement between the GENE-Tango and TGLF-ASTRA profiles with only minor variations close to the GENE left buffer region. As discussed already in Section VI, the flattening of the GENE-Tango profiles at the left boundary might be related to the procedure used in Tango to account for the GENE inner buffer (namely by setting the amplitude of the sources to zero in the inner buffer region).
The most considerable effect of the toroidal rotation is observed on the plasma density. While it was underpredicted in the GENE-Tango simulations of Fig. 12, it is fully recovered when the toroidal plasma rotation is retained in the GENE simulations. In particular, we notice a significant peaking of the density profile primarily located in the outer regions, thus improving the agreement of the GENE-Tango density profile with the experimental measurements and TGLF-ASTRA. These findings are in agreement with the flux-tube results of Ref. [83] suggesting that the E × B shear might enhance inward particle pinch.
Consistently with the results of Fig. 14, we observe that the effect of toroidal rotation on the temperature profiles is localized in the outer core regions, i.e., ρ tor = [0.4−0.6]. This is illustrated clearly in Fig. 16 showing the comparison between the logarithmic temperature gradients obtained with and without toroidal rotation.
However, the ion temperature profile still does not exhibit the pronounced on-axis peaking characteristic of the ASDEX Upgrade discharge #31555 at t = 1.45s. In agreement with previous flux-tube results, these findings suggest that the toroidal rotation can hardly explain the T i peaking in gyrokinetic codes in the plasma core [23,79].
As it will be shown and discussed in Section VII D the increase of the ion temperature on-axis experienced in experiments with large external heating is most likely due to electromagnetic and fast particle effects [23,24,26], which are known to be strongly underestimated by reduced turbulence models. This interpretation is consistent with the results of Fig. 15, showing that electrostatic GENE-Tango simulations retaining toroidal rotation are sufficient to reproduce the temperature and density profiles computed by TGLF-ASTRA, suggesting that the model used in TGLF is essentially electrostatic.

C. Electrostatic GENE-Tango simulation without collisions and without external E × B shear
The effect of collisions on the GENE electrostatic fluxes and plasma profiles is investigated within this Section. In particular, we perform electrostatic GENE-Tango simulations neglecting collisions in the turbulence calculations. To initialize the GENE-Tango iterations, we carry out a reference standalone GENE simulation with the profiles computed by TGLF-ASTRA. These initial turbulent fluxes are illustrated in Fig. 17 by the continuous red line. When comparing these reference results with those obtained in the electrostatic GENE-Tango results without toroidal rotation (Fig. 11), we note an overall increase of the ion heat flux in the absence of collisions by roughly 40% at ρ tor = 0.5. On the other hand, the electron heat flux is reduced by 25% at ρ tor = 0.5. The particle flux undergoes the most violent changes, exhibiting large inward fluxes in the whole radial region of interest. This result is well consistent with the theoretical model derived in Ref. [84]. In particular, collisions FIG. 14. Time-averaged radial profile of the a) ion, b) electron heat fluxes in MW and c) particle flux in 1/s corresponding to the stand-alone GENE simulation with the initial TGLF-ASTRA profiles, respectively, with (red) and without (cyan) toroidal rotation and the last 5 GENE-Tango iterations (blue). The shaded blue area represents the fluctuations of the turbulent fluxes over the last five iterations. The gray areas denote the buffer regions and the black circles the volume integral of the injected particle and heat sources. alter the pure convective term [85,86] in the quasi-linear decomposition of the particle flux by adding a positivedefined contribution in ITG turbulence [87,88]. This term typically balances -together with the diagonal diffusion (proportional to the logarithmic density gradient) -the inward thermo-diffusion term in ITG regimes, leading to a net outward flux (as shown in Fig. 11). However, the outward collisional contribution vanishes in collisionless plasmas, resulting in an inward particle flux if the diffusion term is not large enough to balance this negative contribution.
Although the initial profiles are far from matching the volume integral of the injected particle and energy sources, an excellent agreement is obtained in Fig. 17 after running GENE-Tango for 30 iterations, each covering 150c s /a. This corresponds to an overall GENE run-time of 12.7ms. The Tango relaxation coefficients acting on the plasma profiles and turbulent fluxes are kept constant, respectively, to the values α p = 0.1 and α q = 0.4.
The converged temperatures and density computed by Tango over the last 5 GENE-Tango iterations are compared in Fig. 18 with the TGLF-ASTRA profiles, the experimental measurements, and the GENE-Tango profiles obtained retaining the effect of collisions. Consistently with the reference results (Fig. 17) -showing an increase in the ion heat flux in the absence of collisions -, we note that the ion temperature is further flattened respect to profile including collisions to match the power balance. Moreover, we observe a significant peaking in the plasma density. This is a consequence of the large inward particle flux contribution obtained without collisions. In particular, a positive particle flux -fixed by the external sources -can be matched only by enhancing the positive-defined diagonal diffusion term, hence by increasing the logarithmic density gradient. This dependence of the plasma density with collisionality is consistent with the theoretical results of Ref. [87,89] and the experimental evidence observed in different magnetic confinement devices [89][90][91][92][93][94].
Another observation from Fig. 18, is the flattening of the electron temperature profile respect to the results of Fig. 12. This is due to the enhancement of trappedelectron-mode (TEM) turbulence with the increased density gradient. For completeness, we show in Fig. 19 a comparison of the logarithmic temperature and density gradient for the different cases studied within this section. Given the central role played by collisions in setting the plasma density profile, more accurate operators (such as the Sugama model [95] and the exact Fokker-Planck operator [96]) will be employed in the near future.

D. Electromagnetic GENE-Tango simulation with collisions and without external E × B shear
In the previous Sections, we have investigated the effect of toroidal rotation and collisions on the electrostatic turbulent fluxes and bulk profiles. The inclusion of external E × B shear in the GENE simulations leads to an excellent agreement between the temperature and density profiles computed by GENE-Tango and TGLF-ASTRA. Collisions are found to be essential to model the particle flux correctly and hence to set the plasma density. A minor impact has been identified, on the other hand, on the bulk species temperatures. Interestingly, both the TGLF-ASTRA and GENE-Tango (electrostatic) profiles with/without toroidal rotation do not capture the onaxis peaking of the ion temperature observed in the experiment. The computed profiles start to deviate from the measurements for ρ tor < 0.4. This is the plasma region where electromagnetic effects are expected to become more relevant. In such electromagnetic regimes, an increasing number of studies are showing that presently available reduced turbulence models fail to fully capture electromagnetic effects on plasma micro-turbulence, often leading in gyrokinetic codes to a significant turbulence stabilization. In this Section, we assess the role of finite beta on the bulk profiles by performing GENE-Tango electromagnetic simulations with collisions. We begin by running a GENE standalone simulation to convergence, starting from the reference TGLF-ASTRA profiles. The timeaveraged turbulent fluxes are shown in Fig. 20 by the continuous red line. Compared to the reference GENE electrostatic simulation (cyan line), we note that the inclusion of electromagnetic effects strongly affects the turbulent fluxes in all different channels. To begin with, we observe a significant increase of heat and particle fluxes at ρ tor = 0.1, related to the destabilization of an electromagnetic mode with the TGLF-ASTRA profiles (not shown here). Moreover, while no quantitative differences are found on the ion heat flux for ρ tor > 0.4, the electron heat and particle fluxes increase by roughly ∼ 40% at ρ tor = 0.5. This is consistent with previous flux-tube findings showing that the inclusion of electromagnetic effects in gyrokinetic simulations without supra-thermal particles might destabilize the electron heat flux [97]. However, and more importantly, we observe a significant turbulence suppression in the radial region ρ tor = [0.15 − 0.35], which coincides with the radial locations where the ion temperature profile peaks in the experiments but not in the TGLF-ASTRA and electrostatic GENE-Tango profiles.
To further analyze the impact of the toroidal plasma rotation on the GENE-Tango results (which might affect the turbulent fluxes differently in electromagnetic regimes), the turbulent fluxes obtained in the electromagnetic GENE simulation retaining v tor (see Fig. 7) are added to Fig. 20 (green line). While minor differences are observed on the heat and particle fluxes in the deep core regions (i.e., ρ tor < 0.35), a significant turbulence suppression is found in all channels for ρ tor > 0.4 for a finite toroidal rotation. The electron heat flux decreases with respect to the case with toroidal rotation at ρ tor ≈ 0.65. These results are consistent with the ones shown in Fig. 14, showing an increased effect of a finite E × B in the plasma edge.
Starting from these initial profiles, we performed 28 GENE-Tango iterations, each covering a time domain of t = 150c s /a to reach steady-state. The overall GENE run time measures 11.7ms. For these analyses the relaxation coefficients used in Tango for the plasma profiles and fluxes are respectively α p = 0.15 and α q = 0.4. The beta profile is computed self-consistently from the Tango profiles. The results are illustrated in Fig. 20, showing an excellent agreement between the GENE fluxes and the volume integral of the injected sources. The steady-state temperatures and density profiles computed by Tango over the last five iterations are depicted in Fig. 21 and compared with the electrostatic ones (Fig. 12), the ones with electromagnetic and toroidal rotation (Fig. 8), the TGLF-ASTRA profiles and the experimental measurements. To better visualize the differences among these profiles, we show in Fig. 22 the logarithmic gradients as well.
We start by comparing the electrostatic (Fig. 12) with the electromagnetic GENE-Tango profiles. Interestingly, we observe a significant peaking of the ion temperature -with respect to the electrostatic GENE-Tango profile -starting from the radial location ρ tor = 0.3. Its onaxis value increases from T i = 1.6keV to T i = 2.3keV when electromagnetic effects are retained in the GENE simulations. Due to this pronounced peaking the electromagnetic GENE-Tango temperature profile shows an excellent agreement with the experimental measurements for ρ tor > 0.2. This can be observed also by looking at the logarithmic ion temperature gradient, which exhibits the characteristic increase at ρ tor = [0.2 − 0.3] of the experimental measurements. These findings provide numerical evidence showing that the peaking of the ion temperature observed in several discharges and still not fully captured by reduced models is likely to be caused by electromagnetic effects on plasma turbulence. As discussed in Section VI, the under-prediction of T i for ρ tor < 0.2 in GENE-Tango could be caused by the absence of suprathermal particles in the GENE simulations.
The electron temperature profile undergoes a mild flattening for ρ tor > 0.4 and a minor peaking for ρ tor = 0.3 (more evident in Fig. 22). This is consistent with the reference results of Fig. 20. We note the excellent agreement of the electron temperature compared with the experimental measurements.
The density profile, on the other hand, further deviates from the experimental profile for the GENE-Tango case. In particular, its on-axis value reduces in the electromagnetic simulations. This result is in agreement with the theoretical predictions of Ref. [98], showing that electromagnetic fluctuations largely enhance the outward flux contribution of passing electrons in ITG dominated turbulence regimes, thus leading to an increase in the particle flux. To match the power balance on the particle flux -fixed only by the NBI scheme and particle refueling coming from the SOL neutrals -, Tango further flattens the density profile. The impact of finite beta on the density profile is particularly visible in Fig. 22, where the logarithmic density gradient is depicted.
When we compare the electromagnetic GENE-Tango profiles in Fig. 21 obtained retaining and neglecting toroidal rotation, we observe no qualitative differences on the ion and electron temperatures. Both these profiles are found to quantitatively reproduce the experimental measurements up to ρ tor ≈ 0.2. To measure more clearly the (minor) impact of the toroidal plasma rotation of the temperature profiles we refer the reader to Fig. 22.
We notice -in agreement with the electrostatic findings of Section VII C -that the inclusion of finite toroidal rotation leads to a mild peaking in the outer core regions (ρ tor = [0.3−0.6]). Similarly as observed in Section VII C in the electrostatic limit, the largest effect of a finite E×B is observed on the plasma density. As shown in Fig. 21 and Fig. 22, it undergoes a large peaking throughout the radial domain ρ tor = [0.25 − 0.7].
It is worth mentioning that the electromagnetic turbulence stabilization observed within this Section cannot be related to linearly unstable MHD modes. On the contrary, the steady-state profiles keep these modes sub-marginal. This is shown in Fig. 23, where the linear growth rates and frequencies of radially global linear simulations are shown for the final GENE-Tango profiles (see Fig. 21) for different toroidal mode-numbers and values of beta computed at the center of the GENE radial box. The value of the nominal plasma beta is marked in Fig. 21) by the horizontal black line.

VIII. GENE-TANGO SPEEDUPS AT ASDEX UPGRADE AND EXTRAPOLATIONS TO ITER
Particular attention has been paid throughout this work to optimize the overall GENE run time during the GENE-Tango simulations for the ASDEX Upgrade discharge #31555 at t = 1.45s. We found that the GENE-Tango convergence can be improved by applying different relaxation coefficients on the plasma profiles (α p ) and turbulent fluxes (α q ). More precisely, we used α q = 0.4 and α = [0.05−0.2]. A larger value of α q allows Tango to react more promptly to changes on the turbulent fluxes, thus adjusting the pressure profiles accordingly. On the other hand, the relaxation coefficient acting on the profiles had to be reduced to achieve faster and more stable iterations. This is due to the large sensitivity of the turbulent fluxes observed on even minor changes of the plasma pressure. The convergence speed depends strongly on the Tango extrapolation regions. Their location needs to be carefully adjusted to avoid nonphysical changes of the plasma profiles at the GENE buffer regions, which will inevitably affect the plasma pressure in the core. With these settings, the GENE run time for the turbulent calculation for the ASDEX Upgrade discharge #31555 was t = 150c s /a.
The overall run time of the GENE-Tango simulations with these optimal settings provides an estimate of the overall speedup achieved by properly exploiting the time-FIG. 20. Time-averaged radial profile of the a) ion, b) electron heat fluxes in MW and c) particle flux in 1/s corresponding to the stand-alone GENE simulation with the initial TGLF-ASTRA profiles, respectively, with electromagnetic effects and no toroidal rotation (red), in the electrostatic limit and no toroidal rotation (cyan) and with electromagnetic effects and toroidal rotation (green) and the last 5 GENE-Tango iterations (blue). The shaded blue area represents the fluctuations of the turbulent fluxes over the last five iterations. The gray areas denote the buffer regions and the black circles the volume integral of the injected particle and heat sources.  scale separation between micro-and macroscopic physics in the core of magnetic confinement devices. This can be done by comparing the time simulated by GENE with the energy confinement time for this ASDEX Upgrade discharge. In particular, a single GENE simulation letting the plasma profiles evolve due to the effect of physical sources (flux-driven) is expected to converge on timescales comparable to few confinement times. The energy confinement time is defined as the ratio between the plasma stored energy and the overall integral of the injected sources. This measures τ E ≈ 50ms for the ASDEX Upgrade discharge #31555 at t = 1.45s. To compare more easily the speedup achieved for each of the different cases studied in this paper, every GENE-Tango simulation was started from the reference TGLF-ASTRA profiles. We obtained the following speedups with GENE-Tango on 16 IBM POWER9 AC922 nodes each with 4 Nvidia Volta V100 GPUs on Marconi100.
(i) Section VI: electromagnetic with collisions and external E × B shear -the GENE-Tango run time was 11.7ms. Therefore, we achieved a speedup of τ E /τ GENE ≈ 4.3. The simulation lasted approximately 72h.
(ii) Section VII A: electrostatic with collisions -the GENE-Tango run time was 8.1ms. Therefore, we achieved a speedup of τ E /τ GENE ≈ 6.2. The simulation lasted approximately 24h.
(iii) Section VII B: electrostatic with collisions and external E × B shear -the GENE-Tango run time was 6.1ms. Therefore, we achieved a speedup of τ E /τ GENE ≈ 8.2. The simulation lasted approximately 24h.
(iv) Section VII C: electrostatic without collisions -the GENE-Tango run time was 12.7ms. Therefore, we achieved a speedup of τ E /τ GENE ≈ 4. The simulation lasted approximately 24h. for DIII-D, from Ref. [101] for JET and from Ref. [102,103] for ITER (H98(y, 2) ITER scaling). The ASDEX Upgrade GENE-Tango run-time has been set to 12ms, which is the one obtained in the GENE-Tango simulation retaining the full physics of Section VI.
(v) Section VII D: electromagnetic with collisions -the GENE-Tango run time was 11.7ms. Therefore, we achieved a speedup of τ E /τ GENE ≈ 4.3. The simulation lasted approximately 72h.
These results reveal a significant speedup achieved when exploiting the multiple time-step approach of coupling the global version of the code GENE to the transport-solver Tango. At ASDEX Upgrade, we reached the steady-state solution roughly a factor of 4-8 faster than the expected run time of a single (flux-driven) GENE simulation.
These findings are particularly encouraging when performing extrapolations to larger devices, such as JET or ITER, since the time-step separation between turbulence and transport timescales is expected to scale as (1/ρ * ) 2 and thus the GENE-Tango speedup as (1/ρ * ) 3 considering a fixed radial resolution relative to the gyroradius. In particular, assuming GENE-Tango would converge to the steady-state solution (turbulent fluxes match the volume integral of the injected sources) with comparable speed of the ASDEX Upgrade discharge #31555, we can compute a rough estimate for the overall GENE run time for different devices with GENE-Tango and compare this value to the confinement time. This is shown in Fig. 24 by considering to reach convergence with 40 iterations each corresponding to a GENE run time of 150c s /a. Fig. 24 reveals a considerable increase in the GENE-Tango speedup (when compared to the confinement time, hence to a single flux-driven gyrokinetic simulation) as the size of the device increases. While this is negligible (only a factor of ≈ 1.5) for TCV, it becomes significant for JET with an expected speedup of a factor of ≈ 15 − 18. For the ITER standard scenario, this GENE-Tango convergence criterion leads to a speedup of roughly a factor of 200 compared to a single GENE (flux-driven) simulation for the same case and numerical setup. This is likely to be an underestimation and larger speedups are likely to be achieved. As discussed above, a time-dependent gyrokinetic simulations might require 2-3 confinement times to reach the steady-state solution, thus increasing the GENE-Tango speedup of ITER up to a factor of 600. This would make high-fidelity gyrokinetic simulations up to the transport time-scale possible with the current computing resources also for ITER-like devices for the first time.

IX. CONCLUSIONS
Global nonlinear simulations of different ASDEX Upgrade discharges have been performed with GENE-Tango up to the transport time-scale. We started our analyses by benchmarking GENE-Tango against previously published results obtained by coupling the flux-tube version of the GENE code to the transport-solver Trinity. The results refer to the ASDEX Upgrade shot #13151 at t = 1.35s. We found an excellent agreement between the ion and electron temperature profiles computed by GENE-Trinity and GENE-Tango with minor deviations in the plasma density. The analysis on the logarithmic gradients reveal that GENE-Trinity systematically under-predicts the temperature gradients with respect to GENE-Tango. This is likely to be caused by the local flux-tube approximation, which is known to over-predicts turbulent fluxes for medium-size devices such as ASDEX Upgrade [64], thus leading to reduced logarithmic gradients. The local and global turbulent fluxes are expected to converge as 1/ρ * → ∞ [62,63].
Afterwards, we analyzed the ASDEX Upgrade shot #31555 at t = 1.45s with GENE-Tango. This plasma discharge shows a large peaking of the ion temperature profile, which is not captured by TGLF-ASTRA simulations. This is a known limitation for presently available reduced turbulence codes, which have been shown to strongly under-predict electromagnetic and supra-thermal particle effects on turbulence. These effects are typically found to reduce turbulent transport in gyrokinetic simulations consistently with the experimental signatures observed in different tokamak devices. We demonstrated that GENE-Tango can indeed recover the experimental measurements, correctly capturing the ion temperature peaking up to ρ tor ≈ 0.2. The correct description of the deep core regions (ρ tor < 0.2) is likely to require the ex-tension of GENE-Tango including supra-thermal particle physics. To assess the impact of β-stabilization, MHD modes, collisions and toroidal rotation on the evolution of the plasma profiles and T i peaking in the experiment, we performed GENE-Tango simulations retaining different physical effects in the GENE runs.
The main findings can be summarized as follows: (i) In the electrostatic limit, GENE-Tango well reproduces the electron temperature profile, but shows a rather flat ion temperature and plasma density compared to ones computed with TGLF-ASTRA.
(ii) An excellent agreement between the electrostatic GENE-Tango and TGLF-ASTRA profiles is obtained only when toroidal rotation is retained in the GENE turbulent calculations. The inclusion of a finite E × B shear is found to stabilize turbulent transport in all different channels. However, the effect is mostly localized in the outer core regions ρ tor > 0.4. Moreover, the largest effect of the toroidal rotation is observed on the plasma density, which undergoes a significant peaking for ρ tor > 0.4. While these findings suggest that the TGLF model is essentially electrostatic, they hint that the significant on-axis ion temperature peaking observed in the experiment is not related to effects of the E × B shear, as also previously suggested in Ref. [20]. In particular, both the GENE-Tango (with toroidal rotation) and the TGLF-ASTRA profiles do not exhibit this characteristic increase in the ion logarithmic temperature gradient.
(iii) Collisions play a major role in the evolution of the plasma density profile, partially affecting the temperatures (mainly electron). In particular, we found that collisions contribute to the particle flux by adding a positive contribution (namely a radially outward flux) that, together with the (positivedefined) diagonal diffusion, competes with the inward thermo-diffusion term. When collisions are neglected in the GENE-Tango simulations, the diagonal diffusion can no longer balance the thermodiffusion contribution for the nominal profiles, resulting in negative particle flux. As a result, GENE-Tango strongly increases the amplitude of the logarithmic density gradient to enhance the diagonal diffusion until matching the volume integral of the injected particle flux (coming from both NBI and neutrals coming from the SOL). This process leads to a strong density peaking and a corresponding mild flattening of the electron temperature (due to the enhancement of TEM turbulence). These findings are well consistent with previous theoretical and experimental results, showing a favorable dependence of the density peaking with reduced collisionality [89][90][91][92][93][94].
(iv) The peaking of the ion temperature profile observed in the experiment is recovered in the electromag-netic GENE-Tango simulations with collisions. In particular, we observe in the initial GENE standalone simulation a significant turbulence suppression (respect to the electrostatic case) in all channels in the radial regions where the ion logarithmic temperature gradient increases in the experiment. This leads to a pronounced peaking in the ion profile computed by the GENE-Tango electromagnetic simulation. However, the on-axis temperature of the GENE-Tango ion temperature profile is still under-predicted with respect to the experimental measurements. This is likely to be caused by the absence of supra-thermal particles in the GENE simulations, which are well-known to stabilize ITG turbulence strongly when properly optimized. Fast ions will be included in the GENE-Tango workflow in the near future. Moreover, while the electron temperature is only mildly affected by electromagnetic fluctuations, we observe a flattening of the density profile. This is due to the enhancement of the outward flux of passing electrons in ITG turbulence as discussed in Ref. [98]. We note that the steady-state plasma profiles do not excite any electromagnetic modes, which are submarginal, otherwise leading to a large increase in turbulent trans-port and possibly to a flattening of plasma profiles.