Synthesizing impurity clustering in the edge plasma of tokamaks using neural networks

This work investigates the behavior of impurities in edge plasma of tokamaks using high-resolution numerical simulations based on Hasegawa--Wakatani equations. Specifically, it focuses on the behavior of inertial particles, which has not been extensively studied in the field of plasma physics. Our simulations utilize one-way coupling of a large number of inertial point particles, which model plasma impurities. We observe that with Stokes number ($St$) which characterizes the inertia of particles being much less than one, such light impurities closely track the fluid flow without pronounced clustering. For intermediate $St$ values, distinct clustering appears, with larger Stokes values, {\it i.e.} heavy impurities even generating more substantial clusters. When $St$ is significantly large, very heavy impurities tend to detach from the flow and maintain their trajectory, resulting in fewer observable clusters and corresponding to random motion. A core component of this work involves machine learning techniques. Applying three different neural networks - Autoencoder, U-Net, and Generative Adversarial Network (GAN) - to synthesize preferential concentration fields of impurities, we use vorticity as input and predict impurity number density fields. GAN outperforms the two others by aligning closely with direct numerical simulation data in terms of probability density functions of the particle distribution and energy spectra. This machine learning technique holds the potential to reduce computational costs by eliminating the need to track millions of particles modeling impurities in simulations.


I. INTRODUCTION
The pursuit of a sustainable and abundant energy source has sparked greater attention towards nuclear fusion, which is the process responsible for powering celestial bodies such as stars.Nuclear fusion is distinct from nuclear fission as it involves the fusion of light atomic nuclei to create heavier ones, resulting in the release of a significant amount of energy.Fusion is highly regarded as the ultimate goal in energy production due to its capacity to generate significant power while causing minimal harm to the environment.
Fusion reactors confine and heat a plasma, consisting of ions and electrons, to enable fusion reactions at high temperatures.However, the task of confining plasma poses a considerable challenge due to the presence of multiple instabilities that can result in energy loss.Impurity accumulation in the plasma core can result in heat loss through radiation, leading to a decrease in confinement quality.Therefore, it is crucial to conduct a thorough investigation of impurity in fusion plasma.
The dynamics of magnetically confined plasma flow are influenced by drift-wave turbulence and zonal flows.The numerical high-resolution simulations for modeling the plasma flow carried out within this work are based on the Hasegawa-Wakatani model, which governs cross-field transport by electrostatic drift waves in magnetically confined plasmas.A modified version of the model 1 is likewise investigated to take into account zonal flows.The electric field perpendicular to magnetic field lines is particularly significant because it strongly drives cross-field fluxes, impacting edge pressure profiles and stability. 2,35][6] The resulting fluxes on walls can cause sputtering, erosion, and impurity injection, further degrading confinement and safe operation. 5,7,80][11] This approach involves the analysis of trajectories of numerous tracers.Numerical simulations involve solving equations to determine the trajectories of test particles within a specific velocity field, such as the E × B velocity field.The Lagrangian approach highlights the significant impact of coherent structures on transport.][17] We increase the complexity of these models and add the effect of inertia, which has not been done so far, to the best of our knowledge.Therewith we can study the behavior of heavy and light atoms and see the impact on the impurity distribution.A crucial aspect of understanding inertial impurities is their self-organization, such as clustering and void formation, which can be quantified mathematically by deviation from Poissonian statistics. 18Their dynamics can then be analyzed using tesselation based techniques. 19here is a growing interest in studying spatial patterns in turbulence exploring machine learning. 20,21Our study uses machine learning to estimate mesoscale inertial particle clustering in turbulence, relying on flow field data rather than individual particle information.We apply machine learning tools that have been successfully implemented in the field of hydrodynamics, see Ref. 22 and 23.Three neural networks -Autoencoder, U-Net, and Generative Adversarial Network (GAN) -are trained to predict impurity particle densities from vorticity field snapshots.The results are visually and statistically compared with actual DNS particle data.This approach reduces the cost for predicting particle density distributions.
The paper is structured as follows: Sec.II outlines the theoretical basis used for simulation, introducing the Hasegawa-Wakatani edge plasma model and the motion of inertial impurity particles.Sec.III details the numerical simulations and results.In Sec.IV, we explore machine learning models such as Autoencoders, U-net, and GAN to synthesize impurity concentration fields.Finally, Sec.V summarizes our findings and suggests potential directions for future investigations.

A. Hasegawa-Wakatani model
The numerical simulations are based on Hasegawa-Wakatani (HW) equations for plasma edge turbulence driven by drift-wave instability. 24Our focus is on the twodimensional, slab geometry of the HW model, see e.g.Ref. 10. Fig. 1 depicts a representation of the flow configuration.The magnetic field lines are assumed to be constant, straight and perpendicular to the slab.Ions are considered cold, and the effects of the temperature gradient are ignored.
FIG. 1. Illustration of the two-dimensional slab geometry in the tokamak edge used in the Hasegawa-Wakatani system.The radial direction is represented by x, while y is the poloidal direction.There is an imposed mean plasma density gradient ∇n 0 in the radial direction.The 2D flow is computed within a domain whose size is equivalent to 64 times the Larmor radius (ρ s ).This computational domain is marked with a green square.On the right, a vorticity field for the classical Hasegawa-Wakatani (cHW) is shown.The figure is adapted from Ref. 11. .
The HW model consists of two partial differential equations that describe the time evolution of the plasma electrostatic potential φ and fluctuating plasma density n: All quantities are dimensionless and normalized as described in Refs.11 and 25.The constant D and ν are crossfield diffusion coefficient of plasma density fluctuations and kinematic viscosity, respectively.The constant κ defined as κ ≡ −∂ x ln n 0 , is a measure of the density gradient.The Poisson bracket is defined as: In this system, the electrostatic potential φ acts as a stream-function for the E × B flow velocity, denoted by u = ∇ ⊥ φ , where ∇ ⊥ = (−∂ y , ∂ x ).Thus, we have u x = −∂ φ /∂ y and u y = ∂ φ /∂ x.The vorticity is defined as ω = ∇ 2 φ .The turnover time of turbulent eddies, τ η , is defined as 1/ √ 2Z ms , with Z ms denoting the mean-square vorticity.The adiabaticity parameter c is: where η is the electron resistivity.Here k ∥ is the effective parallel wavenumber.The parameter c controls the phase difference between the electrostatic potential and the plasma density fluctuations.The case c ≫ 1 is known as the adiabatic limit.In this limit, the Hasegawa-Wakatani model reduces to the Hasegawa-Mima equation, and the electrons in the plasma follow a Boltzmann distribution.In that case, there is no phase shift between φ and n.Conversely, when c ≪ 1 (the hydrodynamic limit), the system reduces to a form that is analogous to a two-dimensional Navier-Stokes equation.In this regime, the density fluctuations are passively advected by the E × B flow. 10 The most physically interesting situation in fusion plasma is when c ∼ 1, known as the quasi-adiabatic regime.In this regime, there is a phase shift between φ and n.This phase shift permits the presence of particle transport.
In the context of the classical Hasegawa-Wakatani model (cHW), as detailed above, zonal flows are not observed.To specifically obtain zonal flows, a modified Hasegawa-Wakatani model (mHW) can be considered.This modification was introduced in Ref. 1 where the coupling term c(φ − n) is set to zero for modes k y = 0.

B. Impurity particles model
Previous studies [15][16][17] have assumed that impurity particles perfectly follow the fluid flow without deviating from the flow path.However, this assumption may not always be valid, particularly when the particles have significant mass, resulting in observable inertial effects.

Stokes number
To account for the inertial effects, the concept of the Stokes number (St) is introduced. 26The Stokes number is a dimensionless parameter used to quantify the inertia of particles in a fluid flow.The Stokes number is defined as: where τ p is the impurity particle's relaxation time.This is a measure of how quickly an impurity particle can respond to changes in the fluid flow.It is determined by factors like the particle's size, mass, as well as the viscosity of the fluid.The turbulent eddy turnover time is τ η .This is a measure of how quickly the fluid flow is changing.A high Stokes number suggests that the particle's motion is dominated by its inertia, meaning it will continue moving in its current direction even when the fluid flow changes.Note that the Stokes number is a simplification, and real-world particle-fluid interactions may involve other complex phenomena, like particle-particle interactions, two-way coupling, finite particle-size effects, etc.Nevertheless, the Stokes number provides a valuable prediction of particle behavior in a turbulent flow.

Motion of inertial impurity particles
We assume the behavior of impurities does not influence the plasma's overall dynamics.This study is focused on tracking individual particles using Lagrangian mechanics.The equations governing the motion of these particles are: The particles are denoted by j ( j = 1, . . ., N, where N is the total number of impurity particles in the system).Here, v imp, j is the velocity vector of the j-th impurity particle, and u imp, j represents the fluid velocity vector at the location of the particle, x imp, j .The force resisting the motion of these particles is assumed to be in the form of Stokes drag, which is directly proportional to the relative velocity between the particle and the fluid.8][29] In the initial phase of our study, we excluded electromagnetic forces to focus on establishing a fundamental understanding without introducing additional complexity.The effects of electromagnetic forces on the particles will be integrated incrementally in future studies.

III. NUMERICAL SIMULATIONS A. Flow configurations and physical parameters
The simulation was conducted within a domain with periodic boundary condition that spans an area of 64 2 and it was discretized using a resolution of 1024 2 grid points.The time step was 5 × 10 −4 , while values for diffusivity of plasma density D and kinematic viscosity ν were both set to 5 × 10 −3 .The value for κ was set to 1. Simulations for the flow begin with Gaussian random initial conditions and run until achieving a saturated, fully developed turbulent flow (see Ref. 10).After this phase,10 6 uniformly distributed inertial impurity particles are injected.

B. Numerical schemes
The systems of HW equations, specifically Eq. (1) and Eq. ( 2) are solved using a Fourier pseudospectral method.It is particularly well-suited for problems with periodic boundary conditions and provides excellent accuracy for smooth solutions. 30The equations governing the motion of the impurity particles, namely Eq. ( 5) and Eq. ( 6), are solved using a second-order Runge-Kutta (RK2) scheme.This method provides a balance between accuracy and computational efficiency.Simulating a large ensemble of particles presents a considerable computational challenge.Our solution leverages High-Performance Computing (HPC) to overcome this complexity, employing parallel computing with Message Passing Interface (MPI).For further details on the numerical method and its implementation, we refer to Ref. 10.

C. Results
Fig. 2 shows the vorticity field (without impurity particles) in a fully developed turbulence regime of classical Hasegawa-Wakatani model (cHW) for different c values and modified Hasegawa-Wakatani model (mHW) for c = 2. Larger vortices were exhibited in hydrodynamic cases (cHW, c = 0.01), while smaller ones were observed in quasi-adiabatic (cHW, c = 0.7) and adiabatic cases (cHW, c = 2).For the modified Hasegawa-Wakatani model, with c = 2, zonal flows are presented.To quantify the distribution of kinetic energy of the flow across different scales of the turbulent flow, we compute the Fourier energy spectra: 31 where u(k) = u(x) exp(−i2πk • x) dx is the 2D Fourier transform of the flow velocity, i = √ −1 and k = (k x , k y ) the wave vector.The summation in Eq. ( 7) is carried out over concentric shells in wave number space.In Fig. 3 the Fourier energy spectra E(k) are shown.We observe a peak at smaller wavenumbers for the cases c = 0.01 (cHW and mHW) and c = 2 (mHW), indicating the presence of larger structures.Furthermore, the spectra display a power-law scaling, with an exponent nearing -3.5 for the case c = 0.7.
In this paper, we focus on the quasi-adiabatic regime (c = 0.7), which is the relevant case for edge plasma of tokamaks.For readers interested in zonal flow, relevant discussions can be found in Appendix A. Fig. 4 shows the visualization of vorticity fields and 10 4 superimposed impurity particles for various Stokes numbers in quasi-adiabatic regime.From Fig. 4 we can observe clustering of impurity particles for the case of St = 0.25 and St = 1, while for St = 0 and St = 50, the impurity density field homogeneously fills the entire physical domain and does not show significant spatial correlation.When the particles are very light, the response time of the particles is much smaller than the turnover time of the turbulent eddies, i.e., τ p ≪ τ η (St ∼ 0).Their response time is so quick that they can adjust their velocity to match any changes in the flow immediately.6][17] When τ p ≫ τ η (St ≫ 1), the inertia of the impurity particles is so high that they are almost unaffected by the eddies.Particles are moving ballistically and are randomly distributed, resulting in less clustering and more dispersion.When the response time of the particles is comparable to the eddy turnover time, i.e., τ p ∼ τ η (St ∼ 1), particles are subject to centrifugal effect of turbulent eddies.Eddies can accelerate the particles to a degree where they are thrown out of the eddy due to their inertia.This ejection and the subsequent movement towards areas of lower vorticity result in particle clustering.The clustering of impurity particles results in inhomogeneous poloidal distributions (y direction in the domain) of the impurity density, which can modify radial collisional transport. 32 analyze and statistically interpret the distribution of impurity particles, we increase the number of particles from 10 4 to 10 6 .Displaying the position of each particle in the flow is impractical due to this large quantity.Instead, we visualize the particle density.We achieve this by dividing the entire domain into 512 × 512 boxes and defining the impurity density n imp as the number of particles in each box, i.e. we compute a histogram.Fig. 5 shows the impurity density field for the case of c = 0.7 considering different Stokes numbers.When St = 0, the impurity density is distributed uniformly throughout the flow.As St increases, variations in impurity density become apparent, with some areas exhibiting high density and others low.This uneven distribution indicates the formation of voids and clustering of particles.When St = 50, the impurity density returns to a more uniform distribution.The impurity density fields exhibit patterns similar to those seen in the earlier visualization involving 10 4 impurity particles.
To analyze the impurity density distribution at different scales, we compute the Fourier energy spectra similar to Eq. ( 7) replacing the flow velocity u by the impurity density distribution n imp .The spectral characteristics of impurity density field are shown in Fig. 6 in terms of the Fourier energy spectra.The variance of the impurity density is equivalent to the cumulative sum of the corresponding spectral energy.This is visually represented as the area under the energy spectra curve.As observed, the variance initially increases with the Stokes number up to St = 1, indicating stronger fluctua-

IV. NEURAL NETWORKS FOR SYNTHESIZING PREFERENTIAL CONCENTRATION OF PARTICLES
Simulating flow without impurity particles to obtain the vorticity field is not computationally expensive, but including and tracking 10 6 impurity particles increases costs.Since a specific statistically stationary distribution of impurity density is linked to a specific vorticity field, we propose creating a neural network that can directly predict the distribution of impurity density based on the vorticity field.This network uses as input the vorticity field (easily obtained from DNS) and produces as output the corresponding impurity density field for a specific Stokes number.With this neural network, we can run DNS to get the vorticity field, use it as input for the network, and efficiently predict the impurity density distribution for a specific Stokes number.This approach eliminates the need for costly simulations involving millions of particles.
We will evaluate and compare three different neural network architectures that have previously been utilized in particle-laden hydrodynamic turbulence studies: 22,23 the Autoencoder, U-Net, and Generative Adversarial Network (GAN).To determine the most effective neural network architecture for accurately predicting the impurity density fields, a thorough analysis of the statistical attributes of the impurity density fields generated by these neural networks is carried out.In the following, we will consider the cHW model for the quasi-adiabatic case (c = 0.7).
A. Different machine learning models: Autoencoder, U-Net, and GAN The models used in this study are based on those utilized in Ref. 22 and 23.Fig. 7 displays the architectures tailored from the standard structure, which we have implemented using TensorFlow.
• Autoencoder: Autoencoder 33 is the basic model in our study, comprising of an encoder that compresses the data, and a decoder that reconstructs the output from this compressed version.
• U-Net: U-Net 34 is essentially an Autoencoder but with added skip connections.These connections facilitate non-sequential connections between layers which helps in the preservation of information at different scales throughout the network.
• GAN: GAN (Generative Adversarial Network) consists of two neural networks, 35 a Generator and a Discriminator, trained simultaneously through adversarial processes.The Generator attempts to produce synthetic data, while the Discriminator tries to distinguish between real and synthetic data.We use U-Net as the generator. 36More details can be found in Ref. 22 and  23.In this study, we use the data of one long simulation in the statistically stationary regime.This simulation has a temporal span equivalent to 220 eddy turnover times.Within this simulation, we have captured 400 snapshots of the vorticity field.For each snapshot of the vorticity field, there are 5 corresponding impurity density fields (St = 0.01, 0.05, 0.25, 1, 5).The neural network is then trained to learn the relation between the vorticity field and the corresponding impurity density field.70% of the dataset is for training the model and 30% for evaluating its performance.Vorticity data are normalized between -1 and 1. Impurity density data are normalized between 0 and 1 to enable the use of binary cross-entropy.
The Adam optimizer 37 was used to train the neural networks.Training on a single sample from the training set is termed a 'step'.The Autoencoder and U-net models were trained for 28,000 training steps.The GAN model required more steps for convergence and was trained for 56,000 steps The Python codes were run on an NVIDIA graphics card Tesla V100 32GB.

B. Visual Results
Fig. 8 shows the vorticity field and corresponding impurity density field (St = 1) obtained from DNS for the case of c = 0.7 (cHW).Fig. 9 illustrates impurity density field (St = 1) predicted from Autoencoder, Unet and GAN for c = 0.7 (cHW) using vorticity field as input.When comparing DNS impurity density fields with those predicted by the Autoencoder and U-Net, we find that both DNS and predicted fields exhibit similar structures, but the clustering in the predicted fields is slightly fuzzy.The GAN model successfully predicts void areas and filament-like structures to some degree of success, but it sometimes struggles with predicting the details of regions with impurities.

C. Statistical analysis
To better understand the performance of our models, we analyze the results using the probability density functions (PDFs) and energy spectra of the impurity number density.Fig. 10 shows the PDFs of DNS data and predicted data using three different architectures for St = 1.From Fig. 10 we observe that the GAN model excels in predicting the density distribution, as its predictions closely align with DNS.In contrast, U-Net and Autoencoder exhibit shortcomings, particularly in predicting higher densities.The deviations in the PDFs of U-Net and the Autoencoder suggest that these models may struggle to accurately replicate the intricate characteristics present in high-density regions.
Fig. 11 shows the energy spectra of DNS data and predicted data using three different architectures for St = 1.In Fig. 11 we observe that the GAN's energy spectra are nearly identical to that of the DNS data.This demonstrates the remarkable accuracy of the GAN model in predicting the energy spectra of the impurity density.On the contrary, both U-Net and the Autoencoder exhibit a significant drop at high wave numbers, which correspond to fine structures in the data.This indicates FIG. 7. The neural network architectures are illustrated as: Autoencoder in purple; U-Net in purple and red; GAN in purple, red, and green.that U-Net and the Autoencoder are unable to accurately capture these fine-scale structures and may struggle to represent the intricate details present in the impurity density.
Overall, the analysis of the PDFs and energy spectra reaffirms that GAN outperforms the other two for St = 1.Thus, we utilize GAN in the following to investigate more Stokes numbers values.From the PDF of impurity density as shown in Fig. 12, we can see that for St = 0.25, GAN predicts lower impurity density when compared with DNS.For small St (St = 0.01, 0.05), GAN predicts more particles than DNS.When dealing with a large Stokes number (St = 5), the PDF is similar to the DNS results.In the case of St = 1, the GAN predictions align perfectly with DNS.For energy spectra, as shown in Fig. 13, it is noticeable that the GAN predictions align very well with the DNS.In the GAN predicted energy spectra, notable spikes were identified at frequencies corresponding to the edges and corners of the domain.This could be because the initial encoding of domain periodicity was inadequate.To correct this, adjustments in the code are needed for true periodic convolution.Unfortunately, as of the time of this writing, TensorFlow does not yet provide native support for periodic padding.
In our study, we train the AI model specifically for c = 0.7, which is a typical value for edge plasma.This model predicts impurity density fields for St = 0.01, 0.05, 0.25, 1 and 5 using the vorticity field at c = 0.7 as input.The model's validity is specific to c = 0.7 and the aforementioned Stokes numbers.
Its extrapolability is shown and discussed in Appendix B. Now we know that GAN is the optimal model.Next, we can conduct simulations to generate more vorticity fields without particles, a process which is computationally efficient.We then input these vorticity fields into GAN to predict the impurity density fields.Subsequently, we can perform detailed analyses of transport characteristics and clustering, using methods like multiresolution analysis and Voronoi tessellation.There are several possible ways to enhance the model.One method is to adjust architectural features like the number of layers, filter sizes, and channel counts.Another effective approach is to integrate physical constraints into the model, ensuring that predictions align with established physical laws, such as mass or energy conservation, for more accurate and plausible results.

V. CONCLUSIONS
This study focuses on inertial impurity within fusion plasmas, which is different from previous literature where the focus has typically been on 'passive tracers' following the flow.We used high-resolution numerical simulations based on Hasegawa-Wakatani equations with impurity particles having different inertia.We studied the quasi-adiabatic regime (c = 0.7), typical for Tokamak edge plasma.Our findings show that the inertia of particles, quantified by the Stokes

Vorticity
Impurity density obtained from DNS 1.0 0.5 0.0 0.5 1.0 1.5 2.0 2.5 3.0 0.00 0.04 0.08 0.12 0.16 0.20 As simulating 10 6 or more impurity particles is a computationally intensive task, we utilized machine learning techniques to create a surrogate model.Three neural networks -Autoencoder, U-Net, and GAN -were developed to synthesize impurity densities, using vorticity fields as input.All models produced visually comparable results to simulations, with varying degrees of success.The GAN closely matches DNS values in terms of probability density functions and energy spectra of density.This confirms GAN's effectiveness in modeling both the general distribution and finer structures.Based on the results, future work will involve: Improving GAN by adjusting the architectures or involving physical constraints directly into the model; running simulations for more vorticity fields, then using the GAN to predict impurity fields for analyzing clustering and transport properties; Including electromagnetic interactions for more accurate impurity particle models.

DATA AVAILABILITY STATEMENT
The data that support the findings of this study are available from the corresponding author upon reasonable request.In the modified Hasegawa-Wakatani model at c = 2 with 10 4 impurity particles as shown in Fig. 14, the zonal flow is very strong and the vorticity values are weaker compared to cHW.The particles appear to be distributed evenly across the fluid, and in proximity to the zonal flow, it seems that clustering is not evident To gain deeper insights, we examined the impurity density fields obtained from 10 6 impurity particles.As depicted in Fig. 15, there is some preferential concentration of impurities in the shear layer region for St = 0.25, 1 and 5.The impurity particles are ejected by the zonal flow, and for St=1 the effect seems to be the strongest.For St = 0 and 50 the impurity particles are randomly distributed.

FIG. 4 .
FIG. 4. Visualization of vorticity fields in fully developed turbulence regime and superimposed impurity particles (10 4 impurity particles) for various Stokes numbers in the case of c = 0.7 (cHW).

FIG. 10 .FIG. 11 .
FIG.10.PDFs of impurity density of DNS data and predicted data using three different architectures for St = 1.

FIG. 15 .
FIG. 14. Visualization of vorticity fields in fully developed turbulence regime and 10 4 superimposed impurity particles for various Stokes numbers in the modified Hasegawa-Wakatani model (c = 2).