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 modeling 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 modeling efforts. Such difficulties exist even when examining the same nominal turbulence theory across numerical codes^{3–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 among the greatest challenges to the advancement of plasma modeling. In the fusion community, both comprehensive full-*f* global gyrokinetic^{8–15} and two-fluid^{16–20} direct numerical simulations are under active development to improve predictive modeling 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–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 modeling 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 modeling 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 modeling is edge electric field fluctuations since turbulent $E\xd7B$ 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–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 a direct quantitative cross-examination of plasma turbulence theories by uncovering the turbulent electric field expected in two-fluid 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 modeling based upon Mandell *et al.*^{8} in Sec. II, outline a custom physics-informed machine learning architecture from Mathews *et al.*^{7} suited for the analysis of two-fluid theory in Sec. III, present the first direct comparisons of instantaneous turbulent fields in Sec. IV, and conclude with a summary and future outlook in Sec. V.

## II. GYROKINETIC MODELING

Following the full treatment outlined in Mandell^{36} and Mandell *et al.,*^{8} the long-wavelength gyrokinetic modeling 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 $fs(Z,t)=fs(R,v\u2225,\mu ,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\u2225$, and magnetic moment $\mu =msv\u22a52/2B$. In terms of the gyrocenter Hamiltonian and the Poisson bracket in gyrocenter coordinates with collisions $C[fs]$ and sources *S _{s}*, the gyrokinetic equation is

or, equivalently,

where the gyrokinetic Poisson bracket is defined as

and the gyrocenter Hamiltonian is

The nonlinear phase-space characteristics are given by

Here, $B\u2225*=b\u0302\xb7B*$ is the parallel component of the effective magnetic field $B*=B+(msv\u2225/qs)\u2207\xd7b\u0302+\delta B$, where $B=Bb\u0302$ is the equilibrium magnetic field and $\delta B=\u2207\xd7(A\u2225b\u0302)\u2248\u2207A\u2225\xd7b\u0302$ is the perturbed magnetic field (assuming that the equilibrium magnetic field varies on spatial scales larger than perturbations so that $A\u2225\u2207\xd7b\u0302$ can be neglected). Higher-order parallel compressional fluctuations of the magnetic field are neglected so that $\delta B=\delta B\u22a5$ and the electric field is given by $E=\u2212\u2207\varphi \u2212(\u2202A\u2225/\u2202t)b\u0302$. The species charge and mass are *q _{s}* and

*m*, respectively. In (6), $v\u0307\u2225$ has been separated into a term coming from the Hamiltonian, $v\u0307\u2225H={v\u2225,Hs}$, and another term proportional to the inductive component of the parallel electric field, $(qs/ms)\u2202A\u2225/\u2202t$. 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,

_{s}where $J=B\u2225*$ is the Jacobian of the gyrocenter coordinates and $b\u0302\xb7\u2207\xd7b\u0302\u22480$ so that $B\u2225*\u2248B$. The symplectic formulation of electromagnetic gyrokinetics is utilized with the parallel velocity as an independent variable instead of the parallel canonical momentum $p\u2225$ in the Hamiltonian formulation.^{37,38} We use this notation for convenience, and to explicitly display the time derivative of $A\u2225$, which is a characteristic of the symplectic formulation of electromagnetic gyrokinetics. The electrostatic potential, $\varphi $, is determined by quasi-neutrality,

with the guiding center charge density (neglecting gyroaveraging in the long-wavelength limit) given by

Here, $dw=2\pi \u2009ms\u22121dv\u2225\u2009d\mu =ms\u22121dv\u2225\u2009d\mu \u222bd\alpha $ is the gyrocenter velocity-space volume element $(dv=ms\u22121dv\u2225\u2009d\mu \u2009d\alpha \u2009J)$ with the gyroangle *α* integrated away and the Jacobian factored out (formally, $J$ should also be included in $dw$). The polarization vector is then

where $\u2207\u22a5=\u2207\u2212b\u0302(b\u0302\xb7\u2207)$ 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\xd7B$ 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 modeling practices.^{8,39,40} Adding the nonlinear polarization density along with the second-order $E\xd7B$ 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\u22a5$ 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 guiding-center 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\u2225$, 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 $\u2202A\u2225/\u2202t$, the inductive component of the parallel electric field $E\u2225$,^{41–43}

Writing the gyrokinetic equation as

where $\u2202(Jfs)\u22c6/\u2202t$ denotes all terms in the gyrokinetic equation (including sources and collisions) except $\u2202A\u2225/\u2202t$ 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–46}

where $nsu\u22252=\u222bdwJv\u22252fs,\u20093nsvth,s2=2\u222bdwJ\mu Bfs/ms,\u2009nsu\u2225=\u222bdwJv\u2225fs,\u2009ns=\u222bdwJfs$, and $Ts=msvth,s2$. 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., $\nu \u2260\nu (v)$. Further details about this collision operator, including its conservation properties and discretization, can be found in Ref. 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 $\omega /\Omega s\u226a1$ and $k\u22a5/k\u2225\u226b1$. These orderings express that the charged particle gyrofrequency ($\Omega s=qsB/ms$) in a magnetic field is far greater than the characteristic frequencies of interest (*ω*) with perpendicular wavenumbers ($k\u22a5$) of Fourier modes being far larger than parallel wavenumbers ($k\u2225$). Such properties are generally expected and observed for drift-wave turbulence in magnetically confined fusion plasmas where $\omega \u226a\Omega s$.^{47} An additional “weak-flow” ordering^{48–50} is applied where $vE\xd7B/vth,s\u2248k\u22a5\rho sqs\varphi /Ts\u226a1$, 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\xd7B$ flows are far smaller than the thermal velocity, $vth,s$, as observed in experiment.^{26,52} By constraining gradients of $\varphi $ instead of $\varphi $ itself, this weak-flow ordering simultaneously allows perturbations of order unity ($qs\varphi /Ts\u223c1$) at long wavelengths ($k\u22a5\rho s\u226a1$) and small perturbations ($qs\varphi /Ts\u226a1$) at short wavelengths ($k\u22a5\rho s\u223c1$) 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, $qs\varphi (R+\rho s)\u2212qs\varphi (R)\u2248qs\rho s\xb7\u2207\u22a5\varphi \u226aTs$.^{36} Here $\rho 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 guiding-center limit is eventually taken in the continuum simulations, i.e., $k\u22a5\rho s\u226a1$), 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 $\varphi $ 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 modeling limitations are worth noting for the scope of the present results.

In accordance with Mandell *et al.,*^{8} the gyrokinetic turbulence 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 modeling 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 ($R,\phi ,Z)$ via *R *=* x*, $\phi =(y/\u2009sin\u2009\theta +z\u2009cos\u2009\theta )/Rc,\u2009Z=z\u2009sin\u2009\theta $. The field-line pitch $sin\u2009\theta =Bv/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 an unfavorable magnetic curvature producing the interchange instability that drives edge turbulence. There is no good curvature region to produce a conventional ballooning-mode structure in the current setup. Further, $Rc=R0+a$ is the radius of curvature at the center of the simulation domain, with

*R*

_{0}the major radius and

*a*the minor radius. The curvature operator is

where $B=Baxis(R0/R)e\u0302z$ in the last step and $Baxis$ 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 centered at $(x,y,z)=(Rc,0,0)$ with dimensions $Lx=56\rho i0\u224816.6$ cm, $Ly=100\rho i0\u224829.1$ cm, and $Lz=Lp/sin\u2009\theta =8.0$ m, where $\rho i0=miTi0/qiB0$ and $Lp=2.4$ m approximate the vertical height of the SOL. The velocity-space grid has extents $\u22124vth,s0\u2264v\u2225\u22644vth,s0$ and $0\u2264\mu \u22646Ts0/B0$, where $vth,s0=Ts0/ms$ and $B0=BaxisR0/Rc$. The low-*β* simulations presented here use piecewise-linear (*p *=* *1) basis functions, with $(Nx,Ny,Nz,Nv\u2225,N\mu )=(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 $(Nx,Ny,Nz,Nv\u2225,N\mu )=(16,32,14,10,5)$. For *p *=* *1, one should double these numbers to obtain the equivalent number of grid points for comparison with the 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 $\varphi =A\u2225=0$. The condition $\varphi =0$ effectively prevents $E\xd7B$ flows into walls, while $A\u2225=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 cancelation of the perpendicular magnetic fluctuations at the wall. Also, the magnetic curvature and $\u2207B$ 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 $\varphi $ 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 sheath-limited 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: $Baxis=0.5$ T, $R0=0.85$ m, *a *=* *0.5 m, $Te0=Ti0=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 $v2=v\u22252+2\mu B/ms$. The source density is given by

so that $xS\u22123\lambda S<x<xS+3\lambda S$ delimits the source region, and we take $xS=1.3$ m and $\lambda S=0.005$ m. This localized particle source structure in *z*-space results in plasma ballooning out largely in the center 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, $Psrc$. Since we only simulate a flux-tube-like fraction of the whole SOL domain, $Psrc$ is related to the total SOL power, $PSOL$, by $Psrc=PSOLLyLz/(2\pi RcLpol)\u22480.115PSOL$. Amplitudes are adjusted to approximate $PSOL=5.4$ MW crossing into these open flux surfaces at low-*β* conditions relevant to NSTX.^{52} An artificially elevated density case with $PSOL=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\xd7$ 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 ongoing 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 drift-kinetic 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 Sec. III. Since the gyrokinetic formulation can itself be improved. The turbulence simulations presented are thus a reflection of the current forefront of numerical modeling. 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 five-dimensional phase space along with explicit time-stepping and avoidance of the Ampère cancelation 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. Convection–diffusion 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–64} Following the framework set forth in Mathews *et al.,*^{7} we consider here the widely used first-principles-based two-fluid 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 modeling.^{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 two-fluid theory assumes the plasma is magnetized ($\Omega i\u226b\u2202\u2202t$), collisional ($\nu ei\u226b\u2202\u2202t$), and quasineutral ($\u2207\xb7j\u22480$) with the reduced perpendicular fluid velocity given by $E\xd7B$, diamagnetic, and ion polarization drifts. After neglecting collisional drifts and terms of order $me/mi$, one arrives at the following set of equations (in Gaussian units) governing the plasma's density ($n\u2248ne\u2248ni$), vorticity (*ω*), parallel electron velocity ($v\u2225e$), parallel ion velocity ($v\u2225i$), electron temperature (*T _{e}*), and ion temperature (

*T*),

_{i}^{1,65}

whereby the field-aligned electric current density is $j\u2225=en(v\u2225i\u2212v\u2225e)$, the stress tensor's gyroviscous terms with some FLR effects contain $Gs=\eta 0s{2\u2207\u2225v\u2225s+c[C(\varphi )+C(ps)/(qsn)]}$, and $\eta 0s$ is the viscosity. The convective derivatives are $dsf/dt=\u2202tf+(c/B)[\varphi ,f]+v\u2225s\u2207\u2225f$ with $[F,G]=b0\xd7\u2207F\xb7\u2207G$ and $b0$ representing a constant unit vector parallel to the unperturbed magnetic field. Consistent with the gyrokinetic modeling, the field's magnitude, *B*, decreases over the major radius of the torus ($B\u221d1/R$), and $\kappa =\u2212R\u0302/R$. Curvature drifts are encoded via diamagnetism in this fluid representation^{69} with the curvature operator, $C(f)=b0\xd7\kappa \xb7\u2207f$, following (17) but with an unperturbed field vector distinguishing it from electromagnetic theory. The constant coefficients $\kappa \u2225s$ and $\eta \u2225s$ correspond to parallel thermal conductivity and electrical resistivity, respectively. The three-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., $mi=3.34\xd710\u221227\u2009kg$ and $me=9.11\xd710\u221231\u2009kg$).

To consistently calculate the electric field response and integrate (20)–(25) in time, classically one would compute the turbulent $\varphi $ 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 modeling, we apply a novel technique for computing the turbulent electric field consistent with the drift-reduced Braginskii equations without explicitly evaluating $\u2207\xb7j=0$ nor directly applying the Boussinesq approximation. Namely, this work applies a validated physics-informed deep learning framework^{7} to infer the gauge-invariant $\varphi $ 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 Eqs. (20)–(25) is approximated by its own fully connected neural network with local spatiotemporal points (*x*, *y*, *t*) from a reduced two-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*measurements along with the physical theory encoded by (20) and (24). To be precise, as visualized in Fig. 3 of Mathews

_{e}*et al.,*

^{7}learning proceeds by training the

*n*and

_{e}*T*networks against the average $l2$-norm of their respective errors,

_{e}where ${x0i,y0i,z0i,t0i,ne,0i,Te,0i}i=1N0$ correspond to the set of observed 2-dimensional data and the variables $ne*$ and $Te*$ symbolize the electron density and temperature predicted by the networks, respectively. Properties of all other dynamical variables in the six-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 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*from the two-fluid theory are jointly optimized using

_{e}where ${xfi,yfi,zfi,tfi}i=1Nf$ denote the set of collocation points, and $fne*$ and $fTe*$ 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., ${x0i,y0i,z0i,t0i}i=1N0={xfi,yfi,zfi,tfi}i=1Nf$. This physics-informed 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 $N0=Nf=1000$ using L-BFGS—a quasi-Newton optimization algorithm^{73}—for 20 h 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*(to emulate gas puff imaging

_{e}^{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 physics-informed 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–77} To handle the particular case of two-dimensional turbulence data, we essentially assume a slow variation of dynamics in the *z*-coordinate and effectively set all parallel derivatives to zero. In computational theory comparisons where no such limitations exist, training on *n _{e}* and

*T*in three-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

_{e}^{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=Lz/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 are 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 the electron pressure and electric field fluctuations. As outlined in Sec. 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 modeling in low-*β* helical plasmas with the results visualized in Fig. 1. This multi-network physics-informed deep learning framework enables direct comparison of drift-reduced Braginskii theory with gyrokinetics and bypasses demanding limitations in standard forward modeling 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 evaluated, 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 modeling 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 Fig. 2. To clarify the origins of

*σ*, every time our physics-informed 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

_{PINN}*E*after each optimization. the quantity

_{r}*σ*corresponds to the standard deviation from this collection. The two-fluid model's results displayed in Fig. 2 is thus based upon computing $\u27e8Er\u27e9y$ 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 $\u27e8Er\u27e9y$ profiles predicted by the electromagnetic gyrokinetic simulation and electrostatic drift-reduced Braginskii model are generally in agreement within error bounds at low-

_{PINN}*β*.

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 Fig. 3. As remarked above, multiple realizations are conducted to analyze the sample statistics of the learned turbulent fields consistent with drift-reduced Braginskii two-fluid theory based solely upon the intrinsic scatter during training to account for the stochastic nature of the optimization. By collecting 100 independently trained realizations, the inherent uncertainty linked to this intrinsic scatter can be evaluated as demonstrated in Fig. 4. These discrepancies indicate that the fluctuations in $B\u22a5$, 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 Fig. 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 $lfne$ and $lfTe$ 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, $\Delta Er\u227315\sigma PINN$ in the bottom plot of Fig. 5, where $\Delta Er$ is the difference in the instantaneous

*E*predicted by two-fluid theory and gyrokinetics. This signals that the two models' turbulent

_{r}*E*fluctuations are incompatible at high-

_{r}*β*conditions while quantifying the separation.

Our multi-network physics-informed deep learning framework 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 the electromagnetic effects, such as the magnetic flutter, in determining cross field transport even at $\beta e\u223c2%$, 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 two-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–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 $(Nx,Ny,Nz,Nv\u2225,N\mu )=(48,96,18,10,5)$ would cost roughly 4 × 10^{6} 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 long-wavelength 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 cross field 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 two-dimensional 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. If independently diagnosed measurements of the laboratory plasma's turbulent

*E*exist, the overall turbulence theory can be expressly tested in experiment.

_{r}## 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 two-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 characteristic 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 two-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., $Ti,v\u2225e,v\u2225i$) for full validation of reduced multi-field 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 $\beta p\u2192\u221e$. 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 modeling, 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 modeling 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.

## ACKNOWLEDGMENTS

The authors wish to thank A. Q. Kuang, D. R. Hatch, and G. W. Hammett for insights shared and helpful discussions. Gyrokinetic simulations with the Gkeyll code (https://gkeyll.readthedocs.io/en/latest/) were performed on the Perseus cluster at Princeton University and the Cori cluster at NERSC. All further analysis and codes run were conducted on MIT's Engaging cluster, and we are grateful for assistance with computing resources provided by all teams. The work is supported by the Natural Sciences and Engineering Research Council of Canada (NSERC) by the doctoral postgraduate scholarship (PGS D), Manson Benedict Fellowship, and the U.S. Department of Energy (DOE) Office of Science under the Fusion Energy Sciences program by Contract Nos. DE-SC0014264, DE-SC0014664, and DE-AC02-09CH11466 via the Scientific Discovery through Advanced Computing (SciDAC) program, and No. DE-AR0001263 via the Advanced Research Projects Agency-Energy (ARPA-E) program.

## AUTHOR DECLARATIONS

### Conflict of Interest

The authors have no conflicts to disclose.

## DATA AVAILABILITY

The data that support the findings this study are openly available in Github at https://github.com/AbhilashMathews/PlasmaPINNtheory and https://github.com/ammarhakim/gkyl-paper-inp/tree/master/2021_PoP_GK_Fluid_PINN, Refs. 92 and 93.

## References

*)*(PMLR,