Turbulent field fluctuations in gyrokinetic and fluid plasmas

A key uncertainty in the design and development of magnetic confinement fusion energy reactors is predicting edge plasma turbulence. An essential step in overcoming this uncertainty is the validation in accuracy of reduced turbulent transport models. Drift-reduced Braginskii two-fluid theory is one such set of reduced equations that has for decades simulated boundary plasmas in experiment, but significant questions exist regarding its predictive ability. To this end, using a novel physics-informed deep learning framework, we demonstrate the first ever direct quantitative comparisons of turbulent field fluctuations between electrostatic two-fluid theory and electromagnetic gyrokinetic modelling with good overall agreement found in magnetized helical plasmas at low normalized pressure. This framework is readily adaptable to experimental and astrophysical environments, and presents a new technique for the numerical validation and discovery of reduced global plasma turbulence models.

A key uncertainty in the design and development of magnetic confinement fusion energy reactors is predicting edge plasma turbulence. An essential step in overcoming this uncertainty is the validation in accuracy of reduced turbulent transport models. Drift-reduced Braginskii two-fluid theory is one such set of reduced equations that has for decades simulated boundary plasmas in experiment, but significant questions exist regarding its predictive ability. To this end, using a novel physics-informed deep learning framework, we demonstrate the first ever direct quantitative comparisons of turbulent field fluctuations between electrostatic two-fluid theory and electromagnetic gyrokinetic modelling with good overall agreement found in magnetized helical plasmas at low normalized pressure. This framework presents a new technique for the numerical validation and discovery of reduced global plasma turbulence models.

I. INTRODUCTION
A hallmark of a turbulence theory is the nonlinear connection that it sets between dynamical variables. For ionized gases, perturbations in the plasma are intrinsically linked with perturbations of the electromagnetic field. Accordingly, the electric potential response concomitant with electron pressure fluctuations constitutes one of the key relationships demarcating a plasma turbulence theory and presents a significant consistency test en route to its validation. But to meaningfully evaluate this dynamical interplay across theories or experiment using standard numerical methods is an intricate task. Namely, when comparing global plasma simulations to other numerical models or diagnostics, one must precisely align initial states, boundary conditions, and sources of particles, momentum, and energy (e.g. plasma heating, sheath effects, atomic and molecular processes) 1,2 . These descriptions of physics can be increasingly difficult to reconcile when applying different representations-for example, fluid or gyrokinetic-or when measurements characterizing the plasma are sparse as commonly encountered in thermonuclear fusion devices. Additionally, imperfect conservation properties (both in theory and numerics due to discretization) and the limited accuracy of time-integration schema can introduce errors further misaligning comparative modelling efforts. Such difficulties exist even when examining the same nominal turbulence theory across numerical codes [3][4][5][6] and become magnified for cross-model comparisons 1 . To overcome these classical limitations, we apply a custom physics-informed deep learning framework 7 with a continuous spatiotemporal domain and partial electron pressure observations to demonstrate the first direct comparisons of instantaneous turbulent fields between two distinct global plasma turbulence models: electrostatic drift-reduced Braginskii theory and electromagnetic long-wavelength gyrokinetics. This is significant since validating the accuracy of reduced turbulence theories is amongst the greatest challenges to the advancement of plasma modelling. In the fusion community, both comprehensive full-f global gyrokinetic [8][9][10][11][12][13][14][15] and twofluid [16][17][18][19][20] direct numerical simulations are under active development to improve predictive modelling capabilities and the design of future reactors, but no single code can currently capture all the dynamics present in the edge with sufficient numerical realism. It is therefore of paramount importance to recognize when certain numerical or analytical simplifications make no discernible difference, or when a given theoretical assumption is too risky-to clearly identify the limits of present simulations 1 . Validating reduced plasma turbulence models is accordingly essential. Adaptations of the drift-reduced Braginskii equations have recently investigated several important edge phenomena [21][22][23][24] , but precise quantitative agreement with experiment is generally lacking on a wide scale and their accuracy in fusion reactors is uncertain, especially since kinetic effects may not be negligible as edge plasma dynamics approach the ion Larmor radius and collision frequencies are comparable to turbulent timescales 25,26 . Yet these fluid codes remain in widespread usage for their reduced computational cost, which offers the ability to perform parameter scans with direct numerical simulations that can help uncover underlying physics and optimize reactor engineering designs. Nevertheless, the Vlasov-Maxwell system of equations is the gold standard in plasma theory, but such modelling is prohibitively expensive when wholly simulating the extreme spatiotemporal scales extant in fusion reactors. For the purposes of more tractable yet high fidelity numerical simulations, this system is averaged over the fast gyro-motion of particles to formulate 5D gyrokinetics. Fluid-gyrokinetic comparisons therefore present a way to analyze the robustness of reduced modelling while also illuminating when collisional fluid theories may be insufficient. Quantifying potential shortcomings is thus vital for the testing and development of improved reduced models 27 .
A particularly important quantity for testing predictive modelling is edge electric field fluctuations since turbulent E × B drifts strongly influence edge plasma stability and transport of blob-like structures which can account for over 50% of particle losses during standard operation of tokamaks 26,[28][29][30][31][32] . Resultant turbulent fluxes impacting the walls can cause sputtering, erosion, and impurity injection which adversely affect safe operation and confinement of fusion plasmas 26,33,34 . While explicit comparisons of global nonlinear dynamics across turbulence models were previously impractical using standard numerical simulations, our physics-informed deep learning 35 technique enables direct quantitative cross-examination of plasma turbulence theories by uncovering the turbulent electric field expected in twofluid and gyrokinetic models based upon an electron pressure fluctuation 7 . This validation framework is potentially transferable to magnetized plasma turbulence in contexts beyond fusion energy such as space propulsion systems and astronomical settings as well. Overall, these physics-informed neural networks present a new paradigm in the advancement of reduced plasma turbulence models.
To demonstrate these results, we proceed with a description of gyrokinetic modelling based upon Mandell et al. 8 in Section II, outline a custom physics-informed machine learning architecture from Mathews et al. 7 suited for the analysis of two-fluid theory in Section III, present the first direct comparisons of instantaneous turbulent fields in section IV, and conclude with a summary and future outlook in Section V.

II. GYROKINETIC MODELLING
Following the full treatment outlined in Mandell 36 and Mandell et al. 8 , the long-wavelength gyrokinetic modelling analyzed constitutes nonlinear global electromagnetic turbulence simulations on open field lines using the Gkeyll computational plasma framework. The full-f electromagnetic gyrokinetic equation is solved in the symplectic formulation 37 , which describes the evolution of the gyrocenter distribution function f s (Z,t) = f s (R, v , µ,t) for each species (s = {e, i}), where Z is a phase-space coordinate composed of the guiding center position R = (x, y, z), parallel velocity v , and magnetic moment µ = m s v 2 ⊥ /2B. In terms of the gyrocentre Hamiltonian and the Poisson bracket in gyrocentre coordinates with collisions C[ f s ] and sources S s , the gyrokinetic equation is or, equivalently, where the gyrokinetic Poisson bracket is defined as and the gyrocentre Hamiltonian is The nonlinear phase-space characteristics are given bẏ Here, B * =b·B * is the parallel component of the effective ∇A ×b is the perturbed magnetic field (assuming that the equilibrium magnetic field varies on spatial scales larger than perturbations so that A ∇ ×b can be neglected). Higher-order parallel compressional fluctuations of the magnetic field are neglected so that δ B = δ B ⊥ and the electric field is given by The species charge and mass are q s and m s , respectively. In (6),v has been separated into a term coming from the Hamiltonian,v H = {v , H s }, and another term proportional to the inductive component of the parallel electric field, (q s /m s )∂ A /∂t. In the absence of collisions and sources, (1) can be identified as a Liouville equation demonstrating that the distribution function is conserved along the nonlinear phase space characteristics. Accordingly, the gyrokinetic equation can be recast in the following conservative form, where J = B * is the Jacobian of the gyrocentre coordinates andb·∇ ×b ≈ 0 so that B * ≈ B. The symplectic formulation of electromagnetic gyrokinetics is utilized with the parallel velocity as an independent variable instead of the parallel canonical momentum p in the Hamiltonian formulation 37,38 . We use this notation for convenience, and to explicitly display the time derivative of A , which is characteristic of the symplectic formulation of electromagnetic gyrokinetics. The electrostatic potential, φ , is determined by quasi-neutrality, with the guiding centre charge density (neglecting gyroaveraging in the long-wavelength limit) given by Here, dw = 2π m −1 s dv dµ = m −1 s dv dµ dα is the gyrocentre velocity-space volume element (dv = m −1 s dv dµ dα J ) with the gyroangle α integrated away and the Jacobian factored out (formally, J should also be included in dw). The polarization vector is then where ∇ ⊥ = ∇ −b(b·∇) is the gradient orthogonal to the background magnetic field. A linearized time-independent polarization density, n 0 , is assumed which is consistent with neglecting a second-order E × B energy term in the Hamiltonian. Such an approximation in the scrape-off layer (SOL) is questionable due to the presence of large density fluctuations, although a linearized polarization density is commonly used in full-f gyrokinetic simulations for computational efficiency and reflective of common numerical modelling practices 8,39,40 . Adding the nonlinear polarization density along with the second-order E × B energy term in the Hamiltonian are improvements kept for future work. Consequently, the quasi-neutrality condition can be rewritten as the long-wavelength gyrokinetic Poisson equation, where, even in the long-wavelength limit with no gyroaveraging, the first-order polarization charge density on the left-hand side of (11) incorporates some finite Larmor radius (FLR) or k ⊥ effects in its calculation. It is worth emphasizing that this "long-wavelength" limit is a valid limit of the full-f gyrokinetic derivation since care was taken to include the guidingcenter components of the field perturbations at O(1). Further, although one may think of this as a drift-kinetic limit, the presence of the linearized ion polarization term in the quasineutrality equation distinguishes the long-wavelength gyrokinetic model from versions of drift-kinetics that, for example, include the polarization drift in the equations of motion. The parallel magnetic vector potential, A , is determined by the parallel Ampère equation, Note that we can also take the time derivative of this equation to get a generalized Ohm's law which can be solved directly for ∂ A /∂t, the inductive component of the parallel electric field E 41-43 : Writing the gyrokinetic equation as where ∂ (J f s ) /∂t denotes all terms in the gyrokinetic equation (including sources and collisions) except ∂ A /∂t term, Ohm's law can be rewritten after an integration by parts as To model collisions, we use a conservative Lenard-Bernstein (or Dougherty) operator [44][45][46] , where n s u 2 = dwJ v 2 f s , 3n s v 2 th,s = 2 dwJ µB f s /m s , n s u = dwJ v f s , n s = dwJ f s , and T s = m s v 2 th,s . This collision operator contains the effects of drag and pitch-angle scattering, and it conserves number, momentum, and energy density. Consistent with our present long-wavelength treatment of the gyrokinetic system, FLR effects are ignored. In this work, both like-species and cross-species collisions among electrons and ions are included. The collision frequency ν is kept velocity-independent, i.e. ν = ν(v). Further details about this collision operator, including its conservation properties and discretization, can be found in 46 .
To clarify the approximations undertaken in deriving the gyrokinetic model formulated above and its consequent effects on turbulent fields, the key assumptions are reviewed: The orderings in gyrokinetic theory that effectively reduce the full phase space's dimensionality are ω/Ω s 1 and k ⊥ /k 1. These orderings express that the charged particle gyrofrequency (Ω s = q s B/m s ) in a magnetic field is far greater than the characteristic frequencies of interest (ω) with perpendicular wavenumbers (k ⊥ ) of Fourier modes being far larger than parallel wavenumbers (k ). Such properties are generally expected and observed for drift-wave turbulence in magnetically confined fusion plasmas where ω Ω s 47 . An additional "weak-flow" ordering 48-50 is applied where v E×B /v th,s ≈ k ⊥ ρ s q s φ /T s 1 which allows for large amplitude perturbations. This approximation is also generally valid for electrons and deuterium ions in edge tokamak plasmas 51 as it assumes E × B flows are far smaller than the thermal velocity, v th,s , as observed in experiment 26,52 . By constraining gradients of φ instead of φ itself, this weak-flow ordering simultaneously allows perturbations of order unity (q s φ /T s ∼ 1) at long wavelengths (k ⊥ ρ s 1) and small perturbations (q s φ /T s 1) at short wavelengths (k ⊥ ρ s ∼ 1) along with perturbations at intermediate scales. Alternatively, this approximation can be intuitively viewed as the potential energy variation around a gyro-orbit being small compared to the kinetic energy, Here ρ s is the gyroradius vector which points from the center of the gyro-orbit R to the particle location x. To ensure consistency in the gyrokinetic model at higher order (although the guidingcentre limit is eventually taken in the continuum simulations, i.e. k ⊥ ρ s 1), a strong-gradient ordering is also employed which assumes the background magnetic field varies slowly relative to edge profiles 47 . As noted above, we have taken the long-wavelength limit and neglected variations of φ on the scale of the gyroradius. This yields guiding-center equations of motion which do not contain gyroaveraging operations. While extensions to a more complete gyrokinetic model are in progress, these contemporary modelling limitations are worth noting for the scope of the present results.
In accordance with Mandell et al. 8 , the gyrokinetic tur-bulence is simulated on helical, open field lines as a rough model of the tokamak scrape-off layer at NSTX-like parameters. A non-orthogonal, field-aligned geometry 53 is employed for numerical modelling whereby x is the radial coordinate, z is the coordinate parallel to the field lines, and y is the binormal coordinate which labels field lines at constant x and z. These coordinates map to physical cylindrical coordinates The field-line pitch sin θ = B v /B is taken to be constant, with B v the vertical component of the magnetic field (analogous to the poloidal field in tokamaks), and B the total magnitude of the background magnetic field. The open field lines strike material walls at the top and bottom of the domain consistent with the simple magnetized torus configuration studied experimentally via devices such as the Helimak 54 and TORPEX 55 . This system without magnetic shear contains unfavorable magnetic curvature producing the interchange instability that drives edge turbulence. There is no good curvature region to produce conventional ballooning-mode structure in the current setup. Further, R c = R 0 + a is the radius of curvature at the centre of the simulation domain, with R 0 the major radius and a the minor radius. The curvature operator is where B = B axis (R 0 /R)ê z in the last step and B axis is the magnetic field strength at the magnetic axis. This geometry represents a flux-tube-like domain on the outboard strictly bad curvature side that wraps helically around the torus and terminates on conducting plates at each end in z. Note that although the simulation is on a flux-tube-like domain, it is not performed in the local limit commonly applied in δ f gyrokinetic codes; instead, the simulations are effectively global as they include radial variation of the magnetic field and kinetic profiles. The simulation box is centred at (x, y, z) = (R c , 0, 0) with dimensions L x = 56ρ i0 ≈ 16.6 cm, L y = 100ρ i0 ≈ 29.1 cm, and L z = L p / sin θ = 8.0 m, where ρ i0 = √ m i T i0 /q i B 0 and L p = 2.4 m approximates the vertical height of the SOL. The velocity-space grid has extents −4v th,s0 ≤ v ≤ 4v th,s0 and 0 (32,64,16,10,5) the number of cells in each dimension. At high-β , due to the increased computational cost, the resolution is lowered to (N x , N y , N z , N v , N µ ) = (16,32,14,10,5). For p = 1, one should double these numbers to obtain the equivalent number of grid points for comparison with standard grid-based gyrokinetic codes.
The radial boundary conditions model conducting walls at the radial ends of the domain, given by the Dirichlet boundary condition φ = A = 0. The condition φ = 0 effectively prevents E × B flows into walls, while A = 0 makes it so that (perturbed) field lines never intersect the walls. For the latter, one can think of image currents in the conducting wall that mirror currents in the domain, resulting in exact cancellation of the perpendicular magnetic fluctuations at the wall. Also, the magnetic curvature and ∇B drifts do not have a radial component in this helical magnetic geometry. These boundary conditions on the fields are thus sufficient to guarantee that there is no flux of the distribution function into the radial walls. Conducting-sheath boundary conditions are applied in the z-direction 40,56 to model the Debye sheath (the dynamics of which is beyond the gyrokinetic ordering), which partially reflects one species (typically electrons) and fully absorbs the other species depending on the sign of the sheath potential. This involves solving the gyrokinetic Poisson equation for φ at the z-boundary (i.e. sheath entrance), and using the resulting sheath potential to determine a cutoff velocity below which particles (typically low energy electrons) are reflected by the sheath. Notably, this boundary condition allows current fluctuations in and out of the sheath. This differs from the standard logical sheath boundary condition 57 which imposes zero net current to the sheath by assuming ion and electron currents at the sheath entrance are equal at all times. The fields do not require a boundary condition in the z-direction since only perpendicular derivatives appear in the field equations. The simulations are carried out in a sheathlimited regime but there can be electrical disconnection from the plasma sheath if the Alfvén speed is slow enough. Periodic boundary conditions are used in the y-direction. The simulation parameters roughly approximate an H-mode deuterium plasma in the NSTX SOL: B axis = 0.5 T, R 0 = 0.85 m, a = 0.5 m, T e0 = T i0 = 40 eV. To model particles and heat from the core crossing the separatrix, we apply a non-drifting Maxwellian source of ions and electrons, with source temperature T S = 70 eV for both species and v 2 = v 2 + 2µB/m s . The source density is given by so that x S − 3λ S < x < x S + 3λ S delimits the source region, and we take x S = 1.3 m and λ S = 0.005 m. This localized particle source structure in z-space results in plasma ballooning out largely in the centre of the magnetic field line. The source particle rate S 0 is chosen so that the total (ion plus electron) source power matches the desired power into the simulation domain, P src . Since we only simulate a flux-tube-like fraction of the whole SOL domain, P src is related to the total SOL power, P SOL , by P src = P SOL L y L z /(2πR c L pol ) ≈ 0.115P SOL . Amplitudes are adjusted to approximate P SOL = 5.4 MW crossing into these open flux surfaces at low-β conditions relevant to NSTX 52 . An artificially elevated density case with P SOL = 54.0 MW is also tested to study edge turbulence at high-β . The collision frequency is comparable in magnitude to the inverse autocorrelation time of electron density fluctuations at low-β . For the high-β case, ν is found to be about 10× larger, and the plasma thus sits in a strongly collisional regime. Simulations reach a quasi-steady state with the sources balanced by end losses to the sheath, although there is no neutral recycling included yet which is a focus of ongo-ing work 58 . Unlike in Shi et al. 40 , no numerical heating nor source floors are applied in the algorithm to ensure positivity. In summary, the full-f nonlinear electromagnetic gyrokinetic turbulence simulations of the NSTX plasma boundary region employ the lowest-order, i.e. guiding-center or driftkinetic limit, of the system. Implementing gyroaveraging effects given by the next order terms in advanced geometries is the focus of future work. These present approximations in modern full-f global gyrokinetic simulations should be kept in mind when attributing any similarities or differences to two-fluid theory in the next sections since the gyrokinetic formulation can itself be improved. The turbulence simulations presented are thus a reflection of the current forefront of numerical modelling. A full exposition of the derivation and benchmarking including the energy-conserving discontinuous Galerkin scheme applied in Gkeyll for the discretization of the gyrokinetic system in 5-dimensional phase space along with explicit time-stepping and avoidance of the Ampère cancellation problem can be found in Mandell 36 .

III. MACHINE LEARNING FLUID THEORY
A fundamental goal in computational plasma physics is determining the minimal complexity necessary in a theory to sufficiently represent (anomalous) observations of interest for predictive fusion reactor simulations. Convectiondiffusion equations with effective mean transport coefficients are widely utilized 59,60 but insufficient in capturing edge plasma turbulence where scale separation between the equilibrium and fluctuations is not justified 36,61 . Other reduced turbulence models such as classical magnetohydrodynamics are unable to resolve electron and ion temperature dynamics with differing energy transport mechanisms in the edge [62][63][64] . Following the framework set forth in Mathews et al. 7 , we consider here the widely used first-principles-based twofluid drift-reduced Braginskii equations 25,65 in the electrostatic limit relevant to low-β conditions in the edge of fusion experiments 52 for comparison to gyrokinetic modelling 8 . This is also a full-f 66-68 nonlinear model in the sense that the evolution of the equilibrium and fluctuating components of the solution are merged and relative perturbation amplitudes can be of order unity but evaluated in the fluid limit. This twofluid theory assumes the plasma is magnetized (Ω i ∂ ∂t ), collisional (ν ei ∂ ∂t ), and quasineutral (∇ · j ≈ 0) with the reduced perpendicular fluid velocity given by E × B, diamagnetic, and ion polarization drifts. After neglecting collisional drifts and terms of order m e /m i , one arrives at the following set of equations (in Gaussian units) governing the plasma's density (n ≈ n e ≈ n i ), vorticity (ω), parallel electron velocity (v e ), parallel ion velocity (v i ), electron temperature (T e ), and ion temperature (T i ) 1,65 whereby the field-aligned electric current density is j = en v i − v e , the stress tensor's gyroviscous terms with some FLR effects contain G s = η s 0 2∇ v s + c C (φ ) +C (p s ) /(q s n) , and η s 0 is the viscosity. The convective derivatives are d s f /dt = ∂ t f + (c/B) [φ , f ] + v s ∇ f with [F, G] = b 0 × ∇F · ∇G and b 0 representing a constant unit vector parallel to the unperturbed magnetic field. Consistent with the gyrokinetic modelling, the field's magnitude, B, decreases over the major radius of the torus (B ∝ 1/R), and κ = −R/R. Curvature drifts are encoded via diamagnetism in this fluid representation 69 with the curvature operator, C ( f ) = b 0 × κ · ∇ f , following (17) but with an unperturbed field vector distinguishing it from electromagnetic theory. The constant coefficients κ s and η s correspond to parallel thermal conductivity and electrical resistivity, respectively. The 3-dimensional shearless field-aligned coordinate system over which the drift-reduced Braginskii theory is formulated in the physics-informed machine learning framework exactly matches the gyrokinetic code's geometry. The two-fluid model similarly consists of deuterium ions and electrons with real masses (i.e. m i = 3.34 × 10 −27 kg and m e = 9.11 × 10 −31 kg).
To consistently calculate the electric field response and integrate (20)- (25) in time, classically one would compute the turbulent φ in drift-reduced Braginskii theory by directly invoking quasineutrality and numerically solving the boundary value problem given by the divergence-free condition of the electric current density 70,71 . For the purposes of comparison with gyrokinetic modelling, we apply a novel technique for computing the turbulent electric field consistent with the driftreduced Braginskii equations without explicitly evaluating ∇ · j = 0 nor directly applying the Boussinesq approximation. Namely, this work applies a validated physics-informed deep learning framework 7 to infer the gauge-invariant φ directly from (20) and (24) for direct analysis of electron pressure and electric field fluctuations in nonlinear global electromagnetic gyrokinetic simulations and electrostatic two-fluid theory. To review the pertinent aspects of the machine learning methodology 7 , every dynamical variable in equations (20)-(25) is approximated by its own fully-connected neural network with local spatiotemporal points (x, y,t) from a reduced 2-dimensional domain as the only inputs to the initial layer. The dynamical variable being represented by each network is the sole output. In the middle of the architecture, every network consists of 8 hidden layers with 50 neurons per hidden layer and hyperbolic tangent activation functions all utilizing Xavier initialization 72 . The turbulent electric field consistent with drift-reduced Braginskii theory is then learnt by optimizing the networks against both n e and T e measurements along with the physical theory encoded by (20) and (24). To be precise, as visualized in Figure 3 of Mathews et al. 7 , learning proceeds by training the n e and T e networks against the average L 2 -norm of their respective errors where correspond to the set of observed 2-dimensional data and the variables n * e and T * e symbolize the electron density and temperature predicted by the networks, respectively. Properties of all other dynamical variables in the 6-field turbulence model are unknown. The only existing requirement on the observational data in this deep learning framework is that the measurements' spatial and temporal resolutions should be be finer than the autocorrelation length and time, respectively, of the turbulence structures being analyzed. This scale condition is well-satisfied for the low-β case and marginally-satisfied for the high-β data analyzed. Model loss functions are synchronously learnt by recasting (20) and (24) in the following implicit form and then further normalized into dimensionless form 1,7 . These unitless evolution equations of n e and T e from the two-fluid theory are jointly optimized using where denote the set of collocation points, and f * n e and f * T e are the null partial differential equations prescribed by (28) and (29) in normalized form directly evaluated by the neural networks. The set of collocation points over which the partial differential equations are evaluated correspond to the positions of the observed electron pressure data, i.e.
. This physicsinformed framework is consequently a truly deep approximation of drift-reduced Braginskii theory since every dynamical variable is represented by its own fully-connected deep network which is then further differentiated to build the cumulative computation graph reflecting every single term in the evolution equations. Loss functions are optimized with mini-batch sampling where N 0 = N f = 1000 using L-BFGSa quasi-Newton optimization algorithm 73 -for 20 hours over 32 cores on Intel Haswell-EP processors. We highlight that the only locally observed dynamical quantities in these equations are 2-dimensional views of n e and T e (to emulate gas puff imaging 74 ) without any explicit information about boundary conditions nor initializations nor ion temperature dynamics nor parallel flows. This machine learning framework's collocation grid consists of a continuous spatiotemporal domain without time-stepping nor finite difference schema in contrast with standard numerical codes. All analytic terms encoded in these equations including high-order operators are computed exactly by the neural networks without discretization. In that sense, it is a potentially higher fidelity continuous representation of the continuum equations. While the linearized polarization density-analogous to the Boussinesq approximation-is employed in the gyrokinetic simulations, no such approximations are explicitly applied by our physicsinformed neural networks.
With sparse availability of measurements in fusion experiments, designing diagnostic techniques for validating turbulence theories with limited information is crucial. On this point, we note that our framework can potentially be adapted to experimental measurements of electron pressure [74][75][76][77] . To handle the particular case of 2-dimensional turbulence data, we essentially assume slow variation of dynamics in the zcoordinate and effectively set all parallel derivatives to zero. In computational theory comparisons where no such limitations exist, training on n e and T e in 3-dimensional space over long macroscopic timescales can be easily performed via segmentation of the domain and parallelization, but a limited spatial view away from numerical boundaries with reduced dimensionality is taken to mirror experimental conditions for field-aligned observations 78 and fundamentally test what information is indispensable to compare kinetic and fluid turbulence. To compare the gyrokinetic and two-fluid theories as directly as possible, we analyze the toroidal simulations at z = L z /3 where no applied sources are present. We further remark that, beyond the inclusion of appropriate collisional drifts and sources, this technique is quite robust and generalizable to boundary plasmas with multiple ions and impurities present due to the quasi-neutrality assumptions underlying the two-fluid theory 79 . This turbulence diagnostic analysis technique is hence easily transferable, which permits its systematic application across magnetic confinement fusion experiments where the underlying physical model governing turbulent transport is consistent. The framework sketched can thus be extended to different settings especially in the interdisciplinary study (both numerical and experimental) of magnetized collisional plasmas in space propulsion systems and fusion energy and astrophysical environments if sufficient observations exist. A full description and testing of the overall deep learning technique is provided in Mathews et al. 7 .

IV. QUANTIFYING PLASMA TURBULENCE CONSISTENCY
A defining characteristic of a nonlinear theory is how it mathematically connects dynamical variables. The focus of this work is quantitatively examining the nonlinear relationship extant in plasma turbulence models between electron pressure and electric field fluctuations. As outlined in Section III, using a custom physics-informed deep learning framework whereby the drift-reduced Braginskii equations are embedded in the neural networks via constraints in the form of implicit partial differential equations, we demonstrate for the first time direct comparisons of instantaneous turbulent fields between electrostatic two-fluid theory and electromagnetic long-wavelength gyrokinetic modelling in low-β helical plasmas with results visualized in Figure 1. This multinetwork physics-informed deep learning framework enables direct comparison of drift-reduced Braginskii theory with gyrokinetics and bypasses demanding limitations in standard forward modelling numerical codes which must precisely align initial conditions, physical sources, numerical diffusion, and boundary constraints (e.g. particle and heat fluxes into walls) in both gyrokinetic and fluid representations when classically attempting comparisons of turbulence simulations and statistics. Further, the theoretical and numerical conservation properties of these simulations ordinarily need to be evalu-ated which can be a significant challenge especially when employing disparate numerical methods with differing discretization schema and integration accuracy altogether 1 . We are able to overcome these hurdles by training on partial electron pressure observations simultaneously with the plasma theory sought for comparison. We specifically prove that the turbulent electric field predicted by drift-reduced Braginskii two-fluid theory is largely consistent with long-wavelength gyrokinetic modelling in low-β helical plasmas. This is also evident if analyzing y-averaged radial electric field fluctuations and accounting for the inherent scatter (σ PINN ) from the stochastic optimization employed as displayed in Figure  2. To clarify the origins of σ PINN , every time our physicsinformed deep learning framework is trained from scratch against the electron pressure observations, the learned turbulent electric field will be slightly different due to the random initialization of weights and biases for the networks along with mini-batch sampling during training. To account for this stochasticity, we have trained our framework anew 100 times while collecting the predicted turbulent E r after each optimization. the quantity σ PINN corresponds to the standard deviation from this collection. The two-fluid model's results displayed in Figure 2 is thus based upon computing E r y from 100 independently-trained physics-informed machine learning frameworks. Their turbulent outputs are averaged together to produce the mean, while the scatter associated with the 100 realizations-which approximately follows a normal distribution-comprises the shaded uncertainty interval spanning 2 standard deviations. As visualized, the E r y profiles predicted by the the electromagnetic gyrokinetic simulation and electrostatic drift-reduced Braginskii model are generally in agreement within error bounds at low-β .
In high-β conditions where the particle source is artificially increased by 10×, electromagnetic effects become important and the electrostatic two-fluid theory is markedly errant juxtaposed with electromagnetic gyrokinetic simulations when considering the comparison in Figure 3. As remarked above, multiple realizations are conducted to analyze the sample statistics of the learned turbulent fields consistent with driftreduced Braginskii two-fluid theory based solely upon the intrinsic scatter during training to account for the stochastic nature of the optimization. By collecting 100 independentlytrained realizations, the inherent uncertainty linked to this intrinsic scatter can be evaluated as demonstrated in Figure  4. These discrepancies indicate that the fluctuations in B ⊥ , which are evaluated by solving the parallel Ampère equation (or, equivalently, generalized Ohm's law), cannot be neglected when considering plasma transport across the inhomogeneous background magnetic field as in electrostatic theory. While the fluid approximation is generally expected to be increasingly accurate at high density due to strong coupling between electrons and ions, these results underline the importance of electromagnetic effects even in shear-free high-β plasmas as found in planetary magnetospheres and fusion experiments (e.g. dipole confinement 80 ), and for the first time enables the degree of error between instantaneous fluctuations to be precisely quantified across turbulence models. Uncertainty estimates stemming from the stochastic framework in both regimes are reflected in Figure 5.
One should note that there are novel and different levels of errors to be mindful of in this evaluation. For example, poor convergence arising from nonuniqueness of the turbulent E r found during optimization against the drift-reduced Braginskii equations, or L f ne and L f Te remaining non-zero (and not below machine precision) even after training 7 . These potential errors exist on top of standard approximations in the discontinuous Galerkin numerical scheme representing the underlying gyrokinetic theory such as the implemented Dougherty collision operator. Notwithstanding, when comparing the electrostatic drift-reduced Braginskii theory to electromagnetic long-wavelength gyrokinetic simulations at low-β , the results represent good consistency in the turbulent electric field and all observed discrepancies are mostly within the stochastic optimization's expected underlying scatter. Alternatively, when analyzing high-β conditions, we observe that the electrostatic two-fluid model cannot explain the turbulent electric field in the electromagnetic gyrokinetic simulations.
In particular, ∆E r 15σ PINN in the bottom plot of Figure 5, where ∆E r is the difference in the instantaneous E r predicted by two-fluid theory and gyrokinetics. This signals that the two models' turbulent E r fluctuations are incompatible at high-β conditions while quantifying the separation.
Our multi-network physics-informed deep learning frame- work demonstrates the suitability of electrostatic two-fluid theory as a good approximation of turbulent electric fields in modern gyrokinetic simulations for low-β helical plasmas with sufficient initial and boundary conditions. Conversely, the electrostatic turbulence model is demonstrably insufficient for high-β plasmas. This finding is indicative of the importance of including electromagnetic effects such as magnetic flutter in determining cross-field transport even at β e ∼ 2%. But field line perturbations and the reduced numerical representation of gyrokinetic theory are not the only effects at play causing mismatch at high-β : due to the artificial nature of the strong localized particle source in z-space to produce high-β conditions, parallel dynamics including electron flows along field lines and Ohmic heating effects become increasingly important. These variables should now be observed or learnt first from the 2-dimensional electron pressure measurements to accurately reconstruct the turbulent electric field which signifies a departure from the expected low-β edge of tokamaks. Going forward, consideration of advanced magnetic geometries with squeezing and shearing of fusion plasmas near null-points, which may even couple-decouple upstream turbulence in tokamaks, will be important in the global validation of reduced turbulence models with realistic shaping [81][82][83][84][85] . Also, while there is generally good convergence in the lowβ gyrokinetic simulations at parameters relevant to the edge of NSTX, numerical convergence in the artificially elevated high-β case is currently questionable and a potential source of discrepancy. Running the high-β gyrokinetic simulation with proper collision frequency at an improved resolution of (N x , N y , N z , N v , N µ ) = (48, 96, 18, 10, 5) would cost roughly 4 million CPU-hours to check and we must leave such investigations for future analysis.
As for distinctiveness, we point out that the techniques used for computing the displayed turbulent electric fields in the two cases are markedly different. In particular, the longwavelength gyrokinetic Poisson equation, which is originally derived from the divergence-free condition of the electric current density, is employed in the gyrokinetic simulations. In contrast, simply the electron fluid evolution equations are used to infer the unknown turbulent field fluctuations consistent with drift-reduced Braginskii theory 7 . A principle underlying these models is quasineutrality, but this condition is not sufficient on its own. If one were to apply equilibrium models such as the Boltzmann relation or simple ion pressure balance as expected neoclassically, the turbulent electric field estimates for these nonequilibrium plasmas with nontrivial crossfield transport would be highly inaccurate 7 . Further, no external knowledge of boundary conditions such as sheath effects are explicitly provided to the physics-informed deep learning framework, but this information implicitly propagates from the walls into the observations of electron pressure. This novel approach resultantly permits using limited 2D measurements to compare a global 3D fluid turbulence model directly against 5D gyrokinetics beyond statistical considerations for characterizing non-diffusive intermittent edge transport 61 . All in all, the agreement between drift-reduced Braginskii theory and gyrokinetics supports its usage for predicting turbulent transport in low-β shearless toroidal plasmas, but an analysis of all dynamical variables is required for the full validation of drift-reduced Braginskii theory. Alternatively, when 2dimensional experimental electron pressure measurements are available 74,86 , this technique can be used to infer E r and the resulting structure of turbulent fluxes heading downstream. And, if independently diagnosed measurements of the laboratory plasma's turbulent E r exist, the overall turbulence theory can be expressly tested in experiment.

V. CONCLUSION
To probe the fundamental question of how similar two distinct turbulence models truly are, using a novel technique to analyze electron pressure fluctuations, we have directly demonstrated that there is good agreement in the turbulent electric fields predicted by electrostatic drift-reduced Braginskii theory and electromagnetic long-wavelength gyrokinetic simulations in low-β helical plasmas. As β e is increased to approximately 2%, the 2-dimensional electrostatic nature of the utilized fluid theory becomes insufficient to explain the microinstability-induced particle and heat transport. Overall, by analyzing the interconnection between dynamical variables in these global full-f models, physics-informed deep learning can quantitatively examine this defining nonlinear characteris-FIG. 5. The relative error in the instantaneous radial electric field fluctuations (∆E r /σ PINN ) between electrostatic drift-reduced Braginskii theory and electromagnetic long-wavelength gyrokinetic modelling is displayed in low-β (top) and high-β (bottom) conditions. While all errors in the low-β scenario are generally within approximately 3-4 standard deviations and representative of mostly good quantitative agreement, one must go to over 15σ PINN to fully account for the turbulent fields using an electrostatic two-fluid theory at β e ∼ 2%. The results are based upon collecting 100 independentlytrained physics-informed neural networks to compute the turbulent E r and the intrinsic scatter in these predictions. tic of turbulent physics. In particular, we can now unambiguously discern when agreement exists between multi-field turbulence theories and identify disagreement when incompatibilities exist with just 2-dimensional electron pressure measurements. This machine learning tool can therefore act as a necessary condition to diagnose when reduced turbulence models are unsuitable, or, conversely, iteratively construct and test theories till agreement is found with observations. While we focus on the electric field response to electron pressure in this work, extending the analysis to all dynamical variables (e.g. T i , v e , v i ) for full validation of reduced multifield turbulence models in a variety of regimes (e.g. collisionality, β , closed flux surfaces with sheared magnetic field lines) using electromagnetic fluid theory is the subject of future work. Also, since plasma fluctuations can now be directly compared across models, as gyrokinetic codes begin including fully kinetic neutrals, this optimization technique can help validate reduced source models to accurately account for atomic and molecular interactions with plasma turbulence in the edge of fusion reactors 23,87 since these processes (e.g. ionization, recombination) affect the local electric field. Further progress in the gyrokinetic simulations such as the improved treatment of gyro-averaging effects 37 , collision operators 46 , and advanced geometries 36 will enable better testing and discovery of hybrid reduced theories as well 88 . For example, in diverted reactor configurations, electromagnetic effects become increasingly important for transport near X-points where β p → ∞. A breakdown of Alfvén's theorem in these regions can also arise due to the impact of Coulomb collisions and magnetic shear contributing to an enhanced perpendicular resistivity 89 which presents an important test case of non-ideal effects within reduced turbulence models. While our work supports the usage of electrostatic two-fluid modelling, with adequate initial and boundary conditions, over long-wavelength gyrokinetics for low-β magnetized plasmas without magnetic shear, a comparison of all dynamical variables beyond the turbulent electric field is required for a full validation of the reduced model. Further investigations into reactor conditions may suggest the modularization of individually validated fluid-kinetic turbulence models across different regions in integrated global device modelling efforts 90,91 . This task can now be efficiently tackled through pathways in deep learning as demonstrated by this new class of validation techniques. In addition, precisely understanding the fundamental factors-both physical and numerical-determining the prediction interval, σ PINN , is the subject of ongoing research in analyzing the nature (e.g. uniqueness, smoothness) of chaotic solutions to fluid turbulence equations and the chaotic convergence properties of physics-informed neural networks.