3D ion gyro-orbit heat load predictions for NSTX-U

High power tokamaks operate with divertor heat loads capable of destroying the plasma facing components (PFCs). High fidelity heat load predictions are necessary to ascertain the PFC state for design and during operation. Typical heat flux calculations are 2D, time invariant, and assume that power flows directly along the magnetic field lines (the optical approximation). These assumptions neglect the complex 3D geometries employed to protect the PFCs, the time varying nature of the plasma and PFC thermal state, and the helical trajectories of ions with finite Larmor radii (the gyro-orbit approximation). An integrated software framework, the heat flux engineering analysis toolkit (HEAT), was developed to generate time varying optical heat loads applied to real engineering computer aided design (CAD) (Looby et al 2022 Fusion Sci. Technol. 78 10–27). Recently, an ion-gyro orbit module has been added to HEAT. This module calculates the helical trajectories of ions as they gyrate about the magnetic field lines using kinetic theory macro-particles to accelerate the calculation. First, the new gyro-orbit module will be presented. Next, a comparison to existing research is performed. Finally, an analysis of the gyro-orbit heat loads for NSTX-U is presented for diverted discharges using the engineering CAD models utilized for PFC fabrication. Including these gyro-orbit effects can enhance the PFC performance by ‘smearing’ out the magnetic shadows associated with the castellated fish-scaled geometry. Simultaneously, the helical trajectories can degrade performance when they load narrow regions on edges and corners with high heat fluxes. Analysis of the trade-offs between these competing effects is included, and regions for further investigation are identified.


Introduction and background
The heat fluxes in tokamaks are some of the most extreme in any engineering domain. For example, the perpendicular heat flux incident upon the nose cone of a missile is approximately * Author to whom any correspondence should be addressed.
Original content from this work may be used under the terms of the Creative Commons Attribution 4.0 licence. Any further distribution of this work must maintain attribution to the author(s) and the title of the work, journal citation and DOI.
1 MW m −2 [2], whereas the perpendicular heat flux incident upon a modern tokamak divertor can exceed 300 MW m −2 [3]. Designing plasma facing components (PFCs) that can withstand the harsh environment in the tokamak divertor is a significant challenge, and ultimately the PFC design will constrain tokamak operations. High precision three-dimensional (3D) heat flux predictions are necessary to determine the sensitivity of the tokamak operational boundary to 3D divertor topology. More specifically, the engineering computer aided design (CAD) models used to fabricate PFCs should be used in the heat loading (heat flux, temperature, stress, and other collective effects) analysis and subsequent operational scenario design.
Recently, an open source Python suite, the heat flux engineering analysis toolkit (HEAT), has been developed to predict heat loads using 3D engineering CAD geometries and axisymmetric MHD equilibria [1]. Up until now, HEAT has operated under the assumption that power flows directly along magnetic field lines, which is called the 'optical approximation'. This assumption neglects the helical trajectories of charged particles as they gyrate about the magnetic field lines at finite Larmor radii. This manuscript presents results from a new HEAT module that calculates the heat loads incident upon 3D engineering CAD that arise from these helical trajectories.
The manuscript is organized as follows. Section 1 lays out background information, such as the relevant physics, mathematical formulations, and an introduction to the NSTX-U PFCs. Section 2 describes the HEAT numerical algorithms that simulate charged particle trajectories and the accompanying heat loads on the PFCs. Section 3 provides some gyro-orbit power exhaust predictions for NSTX-U. Section 4 provides a synopsis and describes future work.

The optical approximation
The simplest 3D heat flux calculations assume that power flows directly along the magnetic field lines, which is called the 'optical approximation'. The cyclotron motion of charged particles as they gyrate about the magnetic field lines is ignored. This can be stated mathematically by constraining the velocity to be parallel to the magnetic field, where v opt is the particle velocity vector under the optical approximation, v is the component of velocity parallel to the magnetic field, v is the scalar particle speed, and B is the magnetic field vector. The divertor incident heat flux profile can be calculated using the equation provided in [1] for the optical approximation, where • q target is the surface heat flux at the target location; • q 0 is a scaling coefficient for the parallel heat flux; • ψ is the poloidal flux at the target location; •q (ψ) is the normalized user defined heat flux profile evaluated at ψ; • B target is the magnetic field at the target location; • B origin is the magnetic field at the outboard midplane (OMP); •b is the normalized magnetic field vector at the target; •n is the target surface normal vector; • ¬ is the logical negation (bitwise inversion) operator; • Shadow is a binary variable that is 1 if face is magnetically shadowed.
The optical approximation is useful for fast heat load predictions, but it neglects the true cyclotron motion of particles as they gyrate about the magnetic field lines. Including the cyclotron motion in the heat flux calculation is referred to as the 'gyro-orbit approximation', and it will generate heat flux profiles that are different from the optical approximation. This difference is especially pronounced where the helical trajectories of the particles resonate with the PFC geometry, and can enhance or degrade the PFC performance depending upon the plasma operational scenario.

The gyro-orbit approximation
For identical non-relativistic energies, the Larmor radius of an electron is approximately √ 1836 times smaller than an ion due to the mass difference. For numerical heat load calculations on 3D meshes with surface mesh element lengths ranging from hundreds of microns to several millimeters, it is reasonable to assume that electrons obey the optical approximation because the Larmor radii are small enough that the optical approximation and the gyro-orbit approximation will generate the same heat flux profiles on the mesh. The ions, however, will have Larmor radii that can range from sub-millimeter to several millimeters for relevant divertor energies (i.e. 10 eV-1 keV). The helical trajectories associated with these gyrating ions can cause them to strike the PFC mesh in locations that are magnetically shadowed under the optical approximation. It is therefore desireable to calculate the helical trajectories of ions when calculating the 3D divertor power exhaust profile, as the optical approximation may predict an incorrect heat flux footprint.
The Lorentz force governs a particle's trajectory under electromagnetic forces, where F is the electromagnetic force, q is the particle charge, E is the electric field, v is the particle's velocity vector, and B is the magnetic field. For the results presented herein, the HEAT gyro-orbit module calculates the gyro-orbit trajectories neglecting electric fields ( E = 0) in equation (3). The validity of neglecting electric fields has been explored in research for ITER [4][5][6][7]. Generally, the validity of this assumption is highly objective dependent. If the goal is to have predictions that capture the physics entirely, then ignoring electric fields is undesireable. If the goal is to obtain a fast prediction with a tolerable amount of error then ignoring the electric fields is allowable, as it eliminates the need to run self-consistent PIC simulations and still captures the characteristic features of particle cyclotron motion. For many simulation objectives, such as PFC design scoping studies or inter-shot physics model validation, speed is a priority and the consequences of ignoring electric fields are minor. Future work will aim to incorporate electric fields into the HEAT gyro-orbit calculations. The field line tracing engine that is used in HEAT is called the manifold and footprint tracer (MAFOT) [8]. MAFOT already has the capability to incorporate electric fields for 2D and 3D plasmas in the plasma and sheath, which should simplify including these fields into HEAT analysis. That being said, all HEAT results included herein ignore electric fields in equation (3).
The Lorentz force equation can be solved and integrated to yield the equations of motion for a particle in a magnetic field with no electric field, assuming the magnetic field is aligned to the z-axis.ẍ where |q|B/m can be identified as the cyclotron frequency (ω c ),v i represents a velocity component in the i direction, and the top sign is for ions. These equations of motion can be integrated to yield the particle trajectories, where t represents time, (x 0 , y 0 , z 0 ) represents the guiding center location at t = 0, α represents the gyro orbit phase angle, the ∓ accounts for diamagnetism (− for ions), and v ⊥ can be defined by the following energy sharing relationships, where v represents the particle instantaneous speed, v ⊥ represents the component of velocity perpendicular to the magnetic field, v represents the component of the velocity parallel to the magnetic field, and β is called the velocity phase angle. Section 2.2 describes how these equations are used to map gyro-orbit power onto a 3D CAD mesh in HEAT.

The NSTX-U divertor PFCs
The NSTX-U divertor PFCs have been optimized to handle large heat and particle fluxes, specifically by employing a castellated design [9,10]. The castellations are cubes of graphite approximately 28 mm in width, and each castellation is separated by a 0.5 mm gap designed to arrest thermal diffusion across the PFC bulk. An additional benefit to the gaps is that each castellation can be utilized as a calorimeter to measure the shot integrated deposited energy and the heat flux decay width [11,12]. The height of the castellations varies in the toroidal direction due to a bevel, resulting in magnetic shadows that protect the castellation leading edges from near normal angle of incidence heat fluxes [1]. The toroidal bevel (also called a 'fish-scale') is approximately 1 • , corresponding to a bevel height of 0.44 mm. The heat flux profile on this complex castellated geometry is 3D, and as such it requires a tool like HEAT to generate accurate 3D heat flux predictions. Figure 1 provides a CAD rendering of the NSTX-U divertor, as well as enlarged views of one of the PFC tiles and an individual castellation. The numerous rectangular boxes that can be observed on the PFCs in the CAD figure are each individual castellations. The divertor regions are labeled in the figure: center stack angled surface, inboard divertor vertical, inboard divertor horizontal (IBDH), and outboard divertor. For the results presented in this manuscript, analysis will be performed for individual castellations and for the IBDH tile. Future work will aim to analyze other regions of the NSTX-U divertor with respect to ion gyro-orbit heat loads.
The aforementioned castellated design of the high heat flux PFCs in NSTX-U is inherently comprised of many sharp (i.e. 90 • ) corners and edges. As will be discussed in section 3, these edges and corners can act as gyro-orbit heat load 'amplifiers' because the gyrating particles can dip into the inter-castellation gap and load the corners from multiple sides. The merit of including so many sharp features in the design will be investigated in this manuscript. The geometry presented herein is similar to ITER, which uses tungsten monoblocks that also have fish-scales (bevels) and sharp edges [13]. Interestingly, ITER edge localized mode (ELM) scenarios described in [4] exhibit gyro-radii of similar size to steady state NSTX-U scenarios when accounting for the differences in fuel mass and machine magnetic field, which implies that NSTX-U operation may be able to provide insight about gyro-orbit heat loading for ITER full field ELMs.

HEAT gyro-orbit algorithms
This section outlines the HEAT algorithms that calculate gyroorbit trajectories and calculate the power incident upon the PFC surface meshes. In HEAT the user supplies a gyro-source mesh in the CAD (usually an arbitrarily placed plane hovering just over the divertor), and macro-particles are traced from the mesh onto the PFCs of interest. A macro-particle represents a solution to equations (7)- (11), and can be described as a helical gyro-orbit trajectory that a single particle would follow. The magnetic field line that intersects the center of each triangular mesh element on the gyro-source plane serves as the guiding center for a family of macro-particles that the user will specify. The gyro-source plane must be oriented such that the magnetic field lines intersect it (b ·n = 0), and it must be wide enough to provide total coverage of the PFCs of interest (in flux coordinates) with a margin of several gyro-radii on the edges.
Once the gyro-source has been constructed, the optical heat flux is calculated on its surface, per equation (2), neglecting any magnetic shadows cast upon it from other objects in the CAD (but including self-shadows for non-plane objects). Because the gyro-source is discretized into mesh elements, the heat flux is assumed to be constant across each mesh element and its value is calculated at the mesh center. Each mesh element therefore can be thought of as a flux tube. The power deposited onto a mesh element, j, is calculated by multiplying the heat flux at the mesh center by the mesh element area. The power, P j , deposited onto each mesh element on the gyrosource will be divided and redistributed to the divertor mesh elements via the helical trajectories described in section 1.2 and the distribution functions described in the next section. The power carried by each macro-particle, is therefore a fractional component of the optical power carried at its guiding center.
One reason why we redistribute power (W) instead of energy flux (W m −2 ) is because it provides a self-consistent way of normalizing the heat flux between disparate meshes. The mesh elements on the gyro-source are not aligned to the same flux surfaces as the mesh elements of the divertor PFCs. The resolution and mesh tolerances of each mesh is independently controlled by the user, and is therefore dynamic and adjustable. This implies that there is not necessarily a 1:1 mapping of gyro-source mesh elements to divertor mesh elements; multiple gyro-source mesh elements can strike a single divertor mesh element.
When we calculate the optical heat flux profile on a mesh object in HEAT, we are not mapping heat flux between meshes, but rather mapping each mesh element to a continuous function of poloidal flux. By calculating the poloidal flux coordinate, ψ N , at each divertor mesh element center, finding the corresponding heat flux value at the OMP, q optical (ψ N ), and then using the principle of conservation of flux in a flux tube to scale the heat flux at the divertor mesh element, we avoid the problem of calculating power transfer between meshes. But this is not the case when we are calculating the gyro-orbit heat flux and mapping heat fluxes between meshes, so we must normalize the fluxes by area to account for potential geometric overlap.
Furthermore, if we are to redistribute the optical heat flux that arrives at the gyro-source plane to the divertor PFCs via helical trajectories, we need to not only scale the parallel heat flux to account for the flux expansion of the flux tubes (as is done by the ratio of magnetic field magnitudes in equation (2)), but also scale the flux to account for the 'helically wetted area' that arises from each individual helical trajectory. For the optical approximation, we scale the parallel heat flux using the principle of conservation of magnetic flux in a flux tube. We cannot use this principle for the gyro-orbit approximation as the helical trajectories do not perfectly follow the flux surfaces. Without scaling the heat flux appropriately between the gyrosource plane and the target meshes, we run the risk of violating the conservation of energy.
At each step in the calculation, confirming the conservation of energy represents the most important criteria for a valid HEAT result. The power that is deposited across an infinitely wide gyro-source plane's surface should be equal to the power that crosses into the separatrix, P SOL . Then, the power that falls upon the divertor should again equal P SOL . Power balance is the primary diagnostic that informs the HEAT user about the state of the simulation. It is for these reasons that we redistribute power in the gyro-orbit calculations.

Sampling algorithms
If the macro-particle's guiding center trajectory is known, and time is the independent simulation variable, equations (7)-(11) contain three unknown variables: (α, β, v) which, when chosen, constrain the family of potential trajectory solutions to a single trajectory. HEAT follows a kinetic theory approach, where these variables are sampled according to distribution functions that represent the macroscopic physics. The user specifies the number of discrete samples that are to be taken from each of these distributions, where NgP is the number of gyro phase angle (α) samples, NvP is the number of velocity phase angle (β) samples, NvS is the number of speed (v) samples. These numbers of samples are written in a shorthand notation, (NgP, NvP, NvS). A simulation that includes (NgP, NvP, NvS) samples corresponds to N macro-particles, where N = NgP · NvP · NvS. Each distribution function, (α, β, v), is discussed below.
Far upstream from the PFCs, the distribution of particle speeds is assumed to be Maxwellian. HEAT uses the Maxwellian speed distribution function, where m is the mass of the ion, k is the Boltzmann constant, T is the temperature, and v is the particle speed. This function is derived by integrating the velocity vector Maxwellian distribution over solid angle, based on the assumption that the parallel and perpendicular components of velocity are both independently Maxwellian. This represents a simplification of the true situation, however, because as ions approach the target the parallel component of velocity will become non-Maxwellian, as described in [14][15][16]. For the initial development of the HEAT gyro-orbit algorithms it is assumed that the velocity and speed distributions are purely Maxwellian, but this could be modified in future work by providing independent distribution functions for the perpendicular and parallel components of velocity. In the code this would amount to changing a few lines of Python. We understand that there are many ways to sample velocity, energy, and flux, and we do not claim our method to be the only one. We expect that many distribution functions and sampling techniques will be explored in future work. Currently in HEAT, the Maxwellian speed distribution is sampled such that each speed has approximately equal probability. This is achieved by splitting the cumulative distribution function (CDF) of equation (13) into NvS bins, and making each bin equally probable. Mathematically, the CDF bin boundaries can be described as the CDF values, The bin boundaries are thus spaced at regular intervals throughout the CDF, which has a total range of (0, 1). The velocity slices (HEAT velocity sample locations) are the velocities that correspond to the center point of these CDF bins. The parallel heat flux to the divertor target can be expressed in spherical coordinates using the aforementioned variables (α, β, v), as an adapted form of equation (2.30) from Stangeby [14], This equation is used in HEAT to weight the fraction of energy carried by each macro-particle, as each of the (α, β, v) variables can be separated and sampled independently. The mathematical definition of the energy fraction for each speed (also called a 'velocity slice') is calculated by converting the speed integral in equation (15) to energy and then converting to a bin fraction, where vSFrac is the velocity slice fraction for a sample from the Maxwellian distribution function, and E a and E b are the energy bin boundaries for a specific vSlice. To obtain E a and E b we convert the velocity bin boundaries to kinetic energy using E = 1/2mv 2 . Note that the normalization factors have been omitted from equation (16) because they cancel out. The velocity phase angle, β, defines the relationship between the parallel and perpendicular velocity components, v and v ⊥ . For a given total velocity magnitude, v, there are infinite combinations of v and v ⊥ . In order to constrain v and v ⊥ , a unique combination of these components that satisfies equation (10) must be chosen. This constraint represents the energy sharing between the parallel and perpendicular velocity components. Neglecting particles that flow backwards up the magnetic field line (v > 0) and using the fact that the plasma is diamagnetic (v ⊥ > 0), leaves only the first quadrant of velocity space. Therefore the velocity phase angle ranges from (0, π/2), as observed in equation (15). The PDF of the velocity phase angle can be extracted from equation (15), The sample locations of β are taken to be the center of each bin, where the bin boundaries are defined uniformly in the CDF of equation (17) CDF bins, vP = x : The fraction of energy carried by a macro-particle with velocity phase angle, β, is (19) where β a represents the lower bin boundary, β b represents the upper bin boundary, and we have once again omitted the normalization constants.
The gyro phase angle, α, represents the orientation of the perpendicular velocity vector in the plane orthogonal to the magnetic field line. The gyro phase angle ranges from (0, 2π), and samples are taken such that they are uniformly distributed between these bounds. The fraction of energy carried by each α sample is taken to be

Gyro-orbit power
HEAT samples the (α, β, v) distribution functions to define macro-particles, then follows the trajectories of these macro-particles looking for intersections with mesh elements. Each macro-particle (originating at an arbitrary source mesh defined by the user) is traced and the mesh element index of the intersection location is recorded via a matrix called intersectRecord. After all macro-particles have been traced, the gyro-source plane optical power calculated at the macro-particle's guiding center is reassigned to the corresponding intersection mesh element with an appropriate weighting factor to reflect the numerous macro-particles.
The total optical power deposited across the gyro-source plane can be described by where j is used to iterate over the mesh elements in the gyrosource plane, A represents the gyro-source plane area, and A j represents the area of mesh element j. Once the fractions for gyro phase angle (gPFrac), velocity phase angle (vPFrac), and velocity slice (vSFrac), have been determined, the total power carried by each macro-particle can be calculated by multiplying the power at its guiding center, P j , by the fractions. The power transported from the gyro-source to the target at location y is the superposition of all the macro-particle trajectories that strike y where x is the location of a macro-particle's guiding center on the gyro-source, σ is a function that returns the strike location of the helical trajectory ( x , α, β, v) on the target, the Dirac function is used to only include trajectories that strike the target location y , and h ion is the fraction of power carried by ions as opposed to electrons. This can be discretized to form the formula that is used in HEAT to calculate power redeposition. Explicitly, the ion gyro orbit power mapped from the gyro-source mesh elements, j, to divertor mesh element, i, is described by the following equation, where • i is the mesh element gyro power is being reassigned to • j is the mesh element pseudo-optical power is being reassigned from via macro-particle klm • k corresponds to the macro-particle gyro orbit phase angle, α • l corresponds to the macro-particle velocity phase angle, β • m corresponds to the macro-particle velocity slice (speed, v) from the Maxwellian • P j,optical is the power that would have been assigned to element j using the optical approximation • ionFrac is the fraction of power carried by ions (as opposed to electrons) • gPFrac j,k is the fraction of power prescribed to the k gyro phase angle for distribution j. Equation (20). • vPFrac j,l is the fraction of power prescribed to the l velocity phase angle for distribution j. Equation (19).
• vSFrac j,m is the fraction of power prescribed to the m velocity slice for distribution j. Equation (16). • intersectRecord i, jklm is a 5D binary matrix that is 1 if macro-particle jklm strikes mesh element i, and corresponds to the delta function in equation (22) The portion of equation (23) inside the summation represents the power carried by macro-particle (k, l, m). The summation serves to calculate the power superposition of all of the macro-particles launched from all the source mesh elements.

Gyro-orbit results for NSTX-U
This section provides the results for HEAT gyro-orbit simulations. A simple block geometry is used to compare the HEAT output to existing gyro-orbit research. Convergence studies are performed on the block geometry to ensure that increasing the number of macroparticles converges to a single solution. Then, the engineering CAD models for NSTX-U are used in a series of analyses that generate heat fluxes and temperatures arising from the gyro-orbit approximation.

NSTX-U plasma scenario
Recent analysis has been performed to predict the divertor plasma temperature for the upcoming campaigns on NSTX-U. Predictions indicate that ion temperature will be approximately equal to the electron temperature, and range from 10 eV to 30 eV [17]. The NSTX-U plasma MHD equilibrium used in this analysis was taken from EFIT reconstructions of the 2016 experimental campaign on NSTX-U [18,19]. The onaxis magnetic field is aligned to the negative toroidal direction (counter clockwise as viewed from above the tokamak), and the plasma current flows in the positive toroidal direction. The simulations presented here are for CAD objects that are located directly underneath the strike point, where the heat flux magnitude is largest. When calculated from the equilibrium, the poloidal field magnitude is approximately 0.09 T, and the toroidal field magnitude is approximately 1.1 T, at the strike point on the lower outer divertor.
The heat flux profile for all the simulations contained in this manuscript follows an Eich profile [20] where the heat flux width, λ q , is calculated to be 3.79 mm from the plasma equilibrium assuming 7 MW of P SOL . The Gaussian spreading term, S, for the Eich profile is also calculated from the equilibrium to be 2.58 mm using figure 6(a) from Makowski [21]. More information about this scenario can be found in section 3.4 and figure 10 from [1].
Using the aforementioned magnetic field strengths and deuterium mass, one can perform a napkin calculation to help comprehend the typical sizes of gyro-radii for the NSTX-U case presented in this section. A mean ion gyro-radius at the lower divertor outer strike point can be approximated as 400 microns for a 10 eV (average temperature) plasma and 1.3 mm for a 100 eV plasma, if the energy is shared equally between perpendicular and parallel components of velocity. This is a poor representation of the true situation, however, as there will be distributions of particles with varying gyroradii and helical slope. There will be a family of gyro-radii for a given mean temperature, corresponding to the family of solutions that satisfy equations (7)- (12). For each mean temperature, a user defined number of macroparticles will be sampled from (1) a Maxwellian distribution function, (2) a velocity phase distribution function, and (3) a gyro phase distribution function, as described in section 2.
Ions from the high energy tail of the Maxwellian speed distribution function can have gyro-radii of several millimeters. These large gyro-radii enable particles to access regions that would be shadowed to an ion carrying the mean energy of the distribution. While from a statistical standpoint these high energy ions have lower probability, they can load shadowed regions and therefore are included in the calculation of the gyro-orbit heat load. Additionally, for a fixed total energy (Maxwellian sample), ions whose perpendicular component of velocity carries the majority of the energy will also have larger gyro-radii. For example, a 10 eV NSTX-U deuterium ion's gyro-radius can vary from approximately 10 microns to nearly 600 microns, depending upon how the total energy is divided between the perpendicular and parallel components of velocity (equations (10)-(12)). Sampling from the distributions as described in section 2 results in a family of macro-particles and a corresponding family of gyro-radii that are all traced along each flux tube during a HEAT simulation.

Comparison to existing research
This section benchmarks the HEAT output against research performed in support of ITER for gyro-orbit heat fluxes. Gunn has published several examples of using the ion gyro-orbit approximation to predict heat loads for the ITER tungsten monoblocks [4,5]. The Gunn calculations ignore electric field effects, and as a means of validating this assumption a series of self consistent particle in cell (PIC) analyses that included electric fields were performed by Komm [6]. Experimental validation of the gyro-orbit approximation was performed by Dejarnac [22]. Using the gyro-orbit approximation instead of the optical approximation can alter the surface heat flux profile, predominantly for larger gyro-radii that are of similar length scales to the geometric features of the PFCs. For gyro-radii that are much smaller than the geometric feature sizes, or for simulation mesh representations that are much larger than the gyro-radius, the difference between the optical and gyro-orbit approximations will be negligible. For cases when the gyroradii are of similar length scales to the geometric features, the optical approximation power can be redistributed around the PFC edges by the helical trajectories of the ions, but the integrated power remains approximately the same. The peak heat flux can be greatly increased in the gyro-orbit approximation, when power that is transported around the PFC corner loads a narrow region.
To validate the HEAT algorithms against the poloidal profiles from Gunn [5] and the toroidal profiles from Komm [6], a test scenario was developed in HEAT. To make the problem geometrically simple and computationally fast, a simple block geometry was used. The simple block was created in a CAD model and its surface was placed flush with the top of an NSTX-U castellation in the lower divertor. Figure 2 provides a visualization of this geometry relative to the NSTX-U IBDH tile, as well as close up views of the entire test scenario and the PFC of interest. The PFC of interest is the 3D CAD object where the heat fluxes will be calculated. The PFC of interest extends 5 mm in the poloidal direction and 10 mm in the toroidal direction. These dimensions were chosen because they are sufficiently large to quantify the gyro orbit power redistribution, but small enough to run quickly on a desktop workstation. In the toroidal direction, the PFC of interest protrudes from the upstream PFC by 0.3 mm, which matches the step height in Komm [6]. The adjacent poloidal PFCs are spaced 0.5 mm apart, yet their surfaces are perfectly aligned with the PFC of interest as in Gunn's analyses [4]. Because this test case matches the geometry toroidally in Komm's PIC analysis and poloidally in Gunn's ion-orbit analysis, a comparison to both can be performed from a single set of HEAT results. The average divertor plasma temperature for the HEAT runs was set to 10 eV, which produces slightly smaller average gyro-radii than Gunn (10 eV) and Komm (approx. 5 eV) due to the magnetic field difference, although the HEAT algorithms sample uniformly in velocity phase space, which will result in a variety of gyro-radii. A second HEAT run was then performed with the average divertor plasma temperature set to 100 eV.
Because HEAT predicts the heat loads over the entire 3D PFC surface, the 3D heat fluxes must be collapsed to a 1D profile in order to compare with the results shown in figure 19 from [4] and figures 7-9 from [6]. To this end, an algorithm for generating heat flux profiles along a 1D 'chord' of the PFCs has been developed in HEAT. First, a plane is constructed that will define the chord. For these comparisons the planes are taken to be the centerlines of the PFCs in the longitudinal and transverse directions, as shown in figure 2. Next, all mesh points that are further than a user specified threshold (here +/−0.25 mm) from the plane are omitted. The remaining points are then projected orthogonally down to the plane. Figure 2 shows the points that are included in the heat flux profile for the toroidal (green) and poloidal (magenta) directions. The 1D profile (chord) is constructed by defining the PFC corner to be S pol,tor = 0, similarly to figure 13 in [4], and then plotting the heat fluxes as one walks the PFC surface away from the corner. Due to the non-uniform mesh, there is considerable noise in this method, but it is generally sufficient to compare against the aforementioned results for ITER. Figure 3 shows the 1D heat flux profiles as a function of the poloidal and toroidal chords. The poloidal plot shows an optical leading edge (fuscia dotted trace) with heat fluxes at grazing angles of incidence for S pol < 0 and with a uniform heat flux across the PFC surface. Two ion gyro-orbit profiles are shown: one for a plasma with 10 eV mean ion temperature (blue solid trace), and another for a plasma with 100 eV mean ion temperature (green dotted trace). Both of the ion gyro-orbit profiles peak on both corners of the PFC, where the helical trajectories of ions enables them to strike the corners. On the right corner, ions that descend into the toroidally running gap between adjacent PFCs are scraped off by the optically shadowed edge, resulting in the right peak. It appears that the peak magnitude is bigger for the 100 eV trace, but this could be a mesh artifact as it only represents a single mesh element. On the left corner, ions whose guiding center trajectories descend to just above the PFC surface wrap around the corner and strike the leading edge of the PFC, resulting in the left peak. In both cases, the interaction between the inter-PFC gap width and the gyro radius size is what defines the magnitude of the power that strikes the PFC edges, and the width of the edge loading.
These results match the Gunn results shown in figure 29 from [4], although the plasma scenarios are different (ITER vs NSTX-U). A difference between the Gunn analysis and the HEAT analysis is that HEAT samples from distribution functions to generate a family of macro-particles each with different v ⊥ and v , as described in section 2. This results in a variety of combinations of helical radius and pitch at a given mean plasma temperature, T i . This is in contrast to the Gunn results, where the authors constrain the value of v ⊥ /v to be a single value.
The results from the top plot in figure 3 can also be compared to work performed in support of ITER on COM-PASS [22], where both optical and gyro-orbit heat fluxes on a toroidally-running gap (TG) were predicted with simulation and experimentally observed. In [22] the PFC of interest was a center stack limited PFC with toroidal bevel, and the plasmalimiter contact point was aligned to be approximately on the TG. This configuration constrains the power flow to only be in the positive toroidal direction, yet also allows for changes to the plasma helicity which are observed as changes in heat flux deposition on both sides of the gap.
While figure 1 from [22] does well to capture the heat flux physics for a TG placed on a center stack limiter at the plasma contact point, it must be taken as a local representation of the interaction between magnetic field and surface geometry. The magnetic field vectors plotted in the upper right of each plot in figure 1 from [22] are local to a single spatial location. As one progresses along the TG the values of the magnetic field components will change, albeit very slightly for the case presented therein. For situations in which the TG on a PFC is not on the center stack limiter, the change in magnetic field components can be substantial.
Specifically, for TGs on PFCs in the divertor (all the NSTX-U simulations presented herein), the situation is more complex. Straight TGs running through PFCs in the divertor floor will have inner and outer surfaces whose normals can align with a radial axis, and therefore can have a sign change in (b ·n ) as a function of toroidal angle. This sign change is due to the interaction of a flat surface with the curved magnetic field structure. For these TGs, the side that gets loaded by the parallel heat flux is a function of toroidal angle, can change along a TG, and is periodic with a toroidal mode number defined by the PFC geometry.
Returning to the analysis of the top plot of figure 3 in this manuscript, both the optical and gyro-orbit heat fluxes land on the left corner, similarly to figures 1(a) and (d) in [22]. An optical heat flux magnitude of approximately 5 MW m −2 can be observed on the upstream (left) side of the PFC. If one were to move this plot to a different toroidal location on the PFC, the magnitude and spatial extent would change, as discussed in the preceding paragraphs. The downstream (right) side of the PFC is magnetically shadowed, and there is no optical heat flux. The gyro-orbit heat flux in the top plot in figure 3 lands primarily on the left corner of the PFC, as this surface is aligned such that it intercepts the left-handed helicity of the gyro-orbit helical traces. However, one can also observe a gyro-orbit peak on the right corner in the top plot in figure 3 for both temperatures. The origin of this peak is under investigation, but it could be a result of a family of macro-particles traveling in the TG adjacent the right corner have gyro-radii smaller than the gap width, and can load the optically shadowed face despite it not being the favored side for left handed helicity.
The toroidal plot (bottom) in figure 3 can be directly compared to figures 7-9 from [6], which describes the physics of plasma deposition around a leading edge. In the top plot of figure 3, the PFC of interest is elevated 0.3 mm above the toroidally upstream PFC, resulting in a leading edge that is exposed at near normal angle of incidence to the magnetic field lines. This results in an extreme optical heat flux (fuscia dotted) on the leading edge (where S tor < 0) that exceeds 100 MW m −2 . On the top surface the optical heat flux is uniform. On the exposed leading edge the gyro-orbit heat fluxes are much smaller than the optical heat flux due to the helical pitch of the ion trajectories. The helical motion enables ions with certain gyro phase angles to gyrate up onto the top surface of the PFC, even when their guiding center strikes the leading edge. As S tor increases along the PFC surface, the gyro-orbit heat load decays as particles that have 'gyrated up' are scraped off. The noise in the gyro-orbit traces in figure 3 are a result of the nonuniform 3D CAD mesh being collapsed down to a 1D chord, as well as the discrete sampling from physically continuous gyro-orbit algorithm distribution functions.
Increasing the average plasma temperature increases the ion gyro-radius, and this can also be observed in figure 3. The ions in the 100 eV plasma have a gyro-radius and helical pitch that is an average of √ 10 times the size of the ion radii/pitch in the 10 eV plasma, and this results in a higher heat flux magnitude on the PFC top surface. This elevated flux magnitude is because the 100 eV particles have longer helical pitches, so they can traverse almost the entire toroidal length of the test scenario PFC in a single gyro-period, enabling particles whose guiding centers strike the leading edge to load the entire top surface of the PFC. In contrast, the 10 eV particles with smaller helical pitches are all scraped off within approximately 3 mm of the toroidal leading edge, and the flux normalizes to the optical approximation flux magnitude thereafter. It should be re-emphasized that the 10 eV and 100 eV gyro-orbit values represent the mean plasma temperature, and that HEAT traces many macro-particles that are sampled from the distributions described in section 2. So the fluxes observed in figure 3 represent the contributions of macro-particles with varying helical radii and pitch.
Another difference between the ions in the 100 eV plasma and the 10 eV plasma can be observed on the corners of figure 3. For the poloidal case, the heat loading in the vertical gap (S pol < 0) is across a wider area for the 100 eV plasma when compared to the 10 eV plasma, as a result of the larger gyro-radii enabling the particles to access regions deeper in the gap. Non-zero heat fluxes as far as 1 mm down into the gap (S pol = −1 mm) can be observed for the 100 eV plasma on the left corner of the poloidal plot, which corresponds to the mean gyro-radius of 1.3 mm. For the 10 eV plasma, the mean gyro-radius is approximately 400 um, and in the figure the heat flux trends to 0 at approximately 0.5 mm (S pol = −0.5 mm) below the left corner. For the toroidal case, both the 100 eV and 10 eV plasmas load the 0.3 mm exposed leading edge, but the redistribution of power onto the top surface extends the entire length of the PFC for the 100 eV plasma, as opposed to approximately 3 mm for the 10 eV plasma. This is because the larger helical pitch of the 100 eV plasma ions enables them to deposit energy further downstream than the 10 eV plasma ions. Future work will aim to investigate how plasma temperature changes the PFC heat flux profile, as this example is solely a demonstration of the capabilities on a simple CAD geometry for comparison against existing research, and further study is warranted.
A comparison between the bottom plot in figures 3 and 7 from [6] shows both qualitative and quantitative agreement. The traces in our figure 3 consist of only the ion contribution to the respective heat flux profiles. In [6], there is no ion optical trace, but there is a total optical trace (q e + q i ), as well as an electron optical trace (q e ). Subtraction of the electron optical trace from the total optical trace in figure 7 from [6] would result in a trace with a peak at roughly 15 MW m −2 . This peak is approximately two times the ion gyro-orbit trace peak, q i (PIC), which is similar to the 10 eV plasma in our figure 3. It should be noted that the ion gyro-orbit trace, q i (PIC), includes electric fields whereas our code does not. If a trace of the ballistic ion gyro-orbit profile were included in figure 8 [6], that would provide a more direct comparison.
One way to validate the gyro-orbit output in HEAT is to use energy balance in a manner identical to the optical approximation power balance calculation described in [1], where the power that lands upon each CAD mesh element is summed to calculate the spatially integrated power across the entire PFC surface. The spatially integrated power will only equal the theoretical power if the gyro-orbit trajectories and associated shadows have been correctly calculated. For a magnetic field that does not vary significantly across a PFC, the ion optical and ion gyro-orbit integrated power can be identical but the heat flux profiles can be very different. The caption to figure 3 provides the optical and gyro orbit integrated heat flux values along the 1D chords. For the poloidal (top) plot, the integral of the 10 eV gyro-orbit trace is within 3% of the integral of the optical trace. For the 100 eV plasma the integral of the poloidal trace is approximately 30% higher than the optical approximation, because the larger helical pitches distribute leading edge power along the entire toroidal length of the PFC, as was already discussed. For the toroidal (bottom) plot, the integral of the gyro-orbit traces are both within 4% of the integral of the optical trace, implying that power has simply been redistributed around the PFC surface in different ways for each trace.
While collapsing the 3D HEAT output to a 1D plot is useful for validating the HEAT output against existing research, it does not capitalize on the inherent 3D time varying capabilities of HEAT analysis. HEAT is capable of predicting the time varying heat fluxes on fully featured CAD geometry, as it uses time varying magnetic equilibria and engineering CAD as its inputs. Figure 4 provides a visualization of the 10 eV plasma gyro-orbit heat flux applied to the test scenario PFC of interest The geometry and magnetic field configuration determine the difference between the 3D optical heat flux profile and the 3D gyro-orbit heat flux profile. Geometries that employ many edges, toroidal shaping, and small features, tend to have a stronger difference between the optical and gyro-orbit profiles than smooth axisymmetric surfaces. Future work will investigate if the ion gyro-orbit heat flux profile can be predicted by convoluting a 3D kernel with the optical approximation heat flux profile, rather than tracing macro-particles, but for the work presented herein all results were generated by tracing macro-particles.

Convergence studies
To validate the computational integrity of the HEAT gyro-orbit model, various convergence studies were performed. Each of the studies performed demonstrates convergence, and taken together provide confidence in the HEAT output. That being said, convergence does not necessitate physical reality, and future work will aim to experimentally validate the HEAT ion gyro-orbit module on various machines.
HEAT uses macro-particles to represent the behavior of a chunk of particles from the velocity and spatial domains without having to simulate every single particle. The number of macro-particle simulations is a user input in HEAT, and the first convergence study was performed by incrementally increasing the number of macro-particles while observing the change to the output heat flux profile.
A visual example of the change in heat flux profile as a function of the number of macro-particles is provided in figure 5 for the test scenario CAD model described in section 3.2. Each image in the top figure has a set of numbers in parentheses (NgP, NvP, NvS), corresponding to the number of samples from each dimension of the macro-particle space (gyro phase angle, velocity phase angle, Maxwellian speed distribution). Increasing the number of macro-particles from 1 to 64 (in the figure (1, 1, 1) to (4,4,4)) greatly changes the heat flux footprint. Increasing the number of macro-particles beyond (3,3,3) does little to qualitatively change the heat flux profile on the PFC surface, but does generate a smoother profile.
It should be noted that these results are for a 10 eV deuterium plasma with 0.63 T on axis field (B t0 ), and a maximum mesh element length of 0.2 mm. The mesh size must be considered when analyzing the convergence of the gyro orbit heat flux profiles, as it can greatly change the convergence rate. For example, if the combination of plasma temperature and magnetic field strength yield a gyro radius of 1 mm, then a 20 mm PFC mesh will not be able to resolve the differences between optical and gyro orbit profiles. The mesh size must be chosen so that for the smallest gyro radius in the kinetic domain, aliasing will not occur. From a computational speed vantage point, tracing the trajectories of as few macro-particles as possible is desireable. A useful analysis is to determine how many macroparticles are necessary for convergence of the gyro-orbit predictions with respect to a 'ground truth'. The lower plot in figure 5 displays a histogram for the (NgP, NvP, NvS) HEAT predictions associated with the top figure heat flux profiles. The (5,5,5) run was taken to be the ground truth. Each datapoint, i, corresponds to the difference between the heat flux assigned to mesh element i in the (N, N, N) run compared to the (5, 5, 5) ground truth run. As expected, the results converge as the number of macro-particles is increased. More specifically, increasing the number of macro-particles reduces the data spread (interquartile range) and skew. Other studies were performed that increased the number of samples up to 12 (i.e. 12gP12vP12vS) and the convergence increases for these  figure 9. The combination of the helical trajectory and the magnetic field line curvature create a situation where power can load the 'helically unfavorable' side of the toroidally running gap. Including the full 3D trajectory is necessary to determine the complex interaction between the magnetic field and the 3D geometry. Figure 7 provides heat flux profiles on the castellation of interest. cases as well, although they are not shown here to preserve clarity of the figures. An interesting observation is that when a single macro-particle is used to represent the entire gyro-orbit domain (α, β, v), skewed peaks in the histogram appear, which correspond to resonances between that particular helical pitch and the CAD geometry.

Single castellation analysis
This section provides gyro-orbit approximation predictions for a steady state lower single null diverted discharge for the NSTX-U castellated PFCs, and compares these predictions to the optical approximation. The results provided here are for the lower divertor, specifically the IBDH tile, which employs a castellated graphite design as described in section 1.3 and illustrated in figure 1. Each individual castellation is protected from near normal angle of incidence heat loads on the leading edges by a toroidal fish-scale 0.44 mm in height. There is no poloidal fish-scale, as the poloidally-facing edges do not see angles of incidence as high as the toroidally-facing edges. A complete analysis of the effectiveness of poloidal shaping can be found in [5].
An example visualization of a single gyro-orbit trajectory as it descends into the toroidally-running gap can be observed in figure 6. The figure depicts the CAD geometry (grey), a magnetic field line trace (yellow), and the helical trajectory of a deuterium ion (green). The poloidal coordinate S pol is shown on the castellation of interest, and corresponds to the negative radial direction. The combination of magnetic field line curvature, straight PFC geometry, and the helical trajectory, enables the particle to strike the helically "unfavorable" side of the toroidally running gap. This complex interaction between the geometry and the particle trajectory require a fully 3D model, such as HEAT, to quantify the power deposition. A HEAT gyro-orbit simulation will consist of many such trajectories, as the gyro-orbit distribution functions are sampled per the discussion in the preceding sections.
For a 10 eV (average ion temperature) steady state deuterium plasma, a side by side comparison between the optical and gyro orbit approximation is provided in figure 7. The castellation of interest described in figure 1 and annotated in figure 6 is the CAD that heat fluxes are overlaid onto in figure 7. Four frames are provided in the figure, where the left column is the optical approximation and the right column is the gyro orbit approximation. The top row of the figure provides a view of the castellation from a downstream vantage point, while the bottom row provides a view of the castellation from the upstream vantage point. Each castellation surface is labeled with two coordinates, S tor which extends in the toroidal direction, and S pol which extends in the poloidal direction. The S variables begin at S = 0 on the top surface corner at each respective location. The colorbar has been capped at 15 MW m −2 .
The left column in figure 7 demonstrates some of the characteristic effects of the optical approximation. The top surface consists of a magnetic shadow on the toroidally upstream leading edge, caused by the toroidal fishscale on the upstream (not-pictured) castellation. The heat loaded region on the top of the castellation has a relatively uniform heat flux. The magnetic shadow has a sharp edge, which is characteristic of the optical approximation, and is caused by the relatively 'straight' magnetic field lines when compared to the small 0.5 mm inter-castellation gaps. This shadow reduces the heat loaded area by approximately 15%, which increases the peak heat flux on the uniform region of the PFC. Additionally, an optical hot spot occurs on a poloidally-facing leading edge. This hot spot occurs where the rectangular PFC surface and the curved magnetic field line interact to generate a faceted heat flux on the downstream PFC corner. More specifically, because the optical power does not have a gyrating component, it can descend into narrow gaps between castellations without getting scraped off. The curvature of the magnetic field line in these narrow gaps results in a changingb ·n (from equation (2)) along the flat PFC side surface (i.e. faceting). The changingb ·n causes the peak of the heat flux facet to occur at the PFC corner.
The gyro-orbit approximation heat flux predictions (right column) in figure 7 differ from the optical approximation predictions in the shadows and on leading edges. The optical shadow described in the preceding paragraph is almost completely smeared out by the gyro-orbit heat fluxes. This is because the helical trajectories of ions enable them to dip into the magnetic shadow and strike the optically shadowed surface. A family of gyro-radii and gyro-periods sampled from the distribution functions will have helical slopes that enable them to access these shadowed regions. The heat fluxes in this shadow are generally less than 5 MW m −2 , but they do increase the plasma wetted area of the castellation top surface, thereby reducing the peak heat flux on the uniform region of the PFC.
On the other hand, the corners in the gyro-orbit predictions are hot. The bottom right frame of six illustrates that the heat flux on the castellation poloidally-facing corner is nearly 25 MW m −2 , which is the highest heat load on the entire PFC surface. This hot corner is caused by ions whose helical trajectories enable them to wrap around the castellation corner and strike the vertical edge of the castellation, as well as the horizontal top surface, creating an amplified heat load on the corner.
The optical hot spot that occurs as a result of the faceting effect is nearly absent in the gyro-orbit prediction. Optical power that strictly follows the magnetic field line is able to descend into the gap between castellations and strike the flat castellation edge. Gyro-orbit power, however, is scraped off on the top corners of castellations before descending down into the gaps. A certain subset of macroparticles, with specific helix pitch and gyro-radii, will be able to descend into the gap. In this manner, the gap widths can be interpreted as a bandpass particle filter, but only a small fraction of the gyroorbit distribution of particles will be in the passband, resulting in the disappearance of the optical hot spot. This phenomenon could also be used as a diagnostic, in a manner similar to a Katsumata probe [23].
While viewing the heat flux over an entire castellation is useful for gaining intuition into the gyro-orbit heat fluxes, analyzing a 1D chord in both the toroidal and poloidal directions, as described in section 3.2, is useful for quantitatively providing insight into the heat loads. Analysis in 1D also enables the NSTX-U castellated graphite design to be compared to similar work performed for ITER by Gunn [4] and Komm [6]. The 1D profiles used for the castellation analysis are depicted in figure 7 as S tor and S pol , corresponding to the toroidal and poloidal directions, respectively. A range of mesh elements on either side of the S tor and S pol coordinate lines in figure 7 are 'collapsed' into a single 1D coordinate. The heat flux profiles can then be plotted in 1D, as depicted in figures 8 and 9. Figure 8 illustrates the heat fluxes along the 1D coordinate, S tor , depicted in figure 7. S tor = 0 corresponds to the top surface corner and the castellation top surface is called out with an arrow in purple. Optically shadowed regions are greyed out. At S tor = 0, the optical heat flux (fuscia dotted trace) is zero, due to the magnetic shadow cast by an upstream castellation fish-scale. At approximately S tor = 4 mm, the magnetic shadow ends and the optical heat flux maintains approximately 12 MW m −2 across the remainder of the castellation top surface. The gyro-orbit heat flux (blue trace), however, can access regions that are optically shadowed. For this example the gyro-orbit heat fluxes begin on the S tor < 0 surface, which is designed to be protected by a toroidal fish-scale. In the center of the castellation, the heat fluxes are similar between the optical and gyro-orbit approximations. On the downstream edge of the castellation, at S tor ≈ 28 mm, there is a peak in the gyro-orbit profile. This peak appears because the magnetic field at this location has a negative z component, which enables the helical trajectories of some helices to intercept the downstream optically shadowed face. Figure 9 illustrates the heat fluxes along the 1D coordinate, S pol , depicted in figure 7. S pol = 0 corresponds to the top surface corner and the castellation top surface is called out with an arrow in purple. Optically shadowed regions are greyed out. The heat fluxes in the center of the castellation are similar between the optical and gyro-orbit approximations for the poloidal chord. The decay of the Eich heat flux profile [20] can be observed on the top surface, where the heat flux magnitude changes as a function of poloidal coordinate for both profiles. The primary difference between the two profiles is on the corners, where S pol = 0 and S pol = 28 mm. The gyro-orbit heat flux on the corners is far higher than the optical approximation, because the helical trajectories of ions enables them to wrap around the castellation corners. This spike in heat flux is nearly two times the optical approximation for the poloidally upstream corner, and peaks out at nearly 25 MW m −2 . This high heat flux on the castellation poloidally-facing corners is a characteristic feature of the gyro-orbit heat flux model. Figure 6 illustrates the particle loading on the downstream corner.
In order to quantify how the gyro-orbit heat fluxes will impact the PFCs with respect to their engineering limits, a thermal analysis has been performed using the heat fluxes depicted in figure 7. Figure 10 provides a visualization of the 3D temperature profile for the optical and gyro orbit approximations. The NSTX-U PFCs are temperature limited rather than stress limited, so temperature is the quantity of interest for this discussion. Figure 10 provides six snapshots of the castellation temperature. Each column of snapshots corresponds to a unique point in time, and three time points have been chosen to illustrate the temporal evolution of the castellation temperature. The top row in each column uses the optical approximation heat flux as a boundary condition for the finite volume mesh, and the bottom row in each column uses the gyro-orbit approximation heat flux.  The discharge simulated in figure 10 lasts for 5 s. At the onset, the PFC is assumed to be in thermal equilibrium at 300 K. Once the discharge commences, heat begins to load the top surface and castellation edges as illustrated in figure 7. At 30 ms, which corresponds to the left column in figure 10, the optical approximation temperature begins to peak on the corner that is associated with the faceted optical hot spot. Simultaneously, the gyro orbit temperature peaks on the narrow poloidally-facing leading edge that is exposed to the helical ion gyro orbits. In the center column of figure 10, at 3200 ms, the castellation surfaces begin to show temperatures in excess of 1873 K, which is the thermal limit for the NSTX-U divertor. A large continuous region of the optical castellation surface is above this temperature at 3200 ms, whereas the gyro castellation has two distinct regions of high temperature (on the edge and in the center). The optical hot spot that arises from the faceting effect in the optical approximation drives a faster temperature increase than the narrow helical leading edge observed in the optical approximation, despite the helical leading edge having a higher peak heat flux. At the final timestep in figure 10 (right column), both the heat flux models have temperature profiles that have exceeded 2500 K. The optical hot spot in the optical approximation reaches 2720 K, and the gyro-orbit helical leading edge reaches 2581 K. Both of these simulations exceeded the engineering limit of NSTX-U graphite, which is 1873 K.
In figure 11, the peak temperatures over the entire castellation for the optical (green dotted trace) and gyro-orbit (blue solid trace) approximations are plotted. These traces represent the maximum temperature in any volume mesh cell over the entire castellation at each timestep. The engineering temperature limit has also been overlaid onto the plot at 1873 K. The simulated discharge consists of a 5 s pulse followed by a 500 ms cooldown period. At approximately 2.21 s, the optical approximation temperature exceeds the engineering limit. At approximately 2.44 s, the gyro orbit trace also exceeds the engineering limit. At 5.01 s, both traces reach their maximum temperature for the entire discharge. The optical heat flux generates a temperature of 2720 K, and the gyro-orbit heat flux generates a temperature of 2581 K. Both traces are within 5% of each other at the peak.
It is somewhat surprising that the gyro-orbit heat flux does not generate a higher peak temperature than the optical heat flux, given that the gyro-orbit heat flux is nearly 25 MW m −2 on the poloidally-facing leading edge. The optical heat flux generates an optical hot spot on the vertical faceted corner, but it has a far lower magnitude of 5 MW m −2 . The poloidallyfacing leading edge that is loaded by the gyro-orbit heat flux is extremely narrow, extending approximately 0.2 mm on either side of the corner. While the peak of this heat flux is 25 MW m −2 , its spatial extent is small, resulting in a small spatially integrated deposited energy. The optical heat flux extends over 3 mm down the side of the castellation, so although the optical hot spot heat flux is smaller in magnitude, the spatially integrated deposited energy is higher. The approximate area of the hot 25 MW m −2 leading edge in the gyro-orbit profile is ≈5.0 mm 2 , while the area of the optical hot spot is ≈33.1 mm 2 . The bulk material can diffuse the energy from the narrow gyro-orbit leading edge more easily than the energy from the wider optical hot spot. Additionally, the magnetic shadow cast onto the toroidally-facing leading edge by the upstream castellation fish-scale reduces the heat flux wetted area by approximately 15% for the optical approximation. For the gyro-orbit approximation, however, this shadow gets filled in by ions with helical trajectories that can load the surface inside the magnetic shadow, thereby increasing the plasma wetted area and decreasing the peak heat flux on the uniform region on the top surface. The combination of these effects results in a lower peak temperature for the ion gyro-orbit approximation.
Simulating the gyro-orbit heat fluxes takes far longer than simulating the optical heat fluxes. For situations when speed is the primary motive, such as physics operation or when generating a large number of HEAT results for a database, using the gyro-orbit approximation may not be feasible. In these cases, the optical approximation calculation can be completed quickly (usually within a few minutes for medium resolutions), and can be utilized as a conservative estimate of the heat flux. Because the peak temperature using the optical approximation will generally be higher than the gyro-orbit approximation,  Temperatures correspond to figure 10. For both heat flux models, the peak temperature is reached at the end of the discharge when t = 5 s. The peak temperature is 2720 K for the optical approximation and 2581 K for the gyro-orbit approximation. Optical trace crosses the temperature limit (1873 K) at 2.21 s and gyro-orbit trace crosses the limit at 2.44 s. the engineering limits will be reached earlier when using the optical heat flux. That being said, there are certain situations when the gyro-orbit calculation is imperative. These include simulations that have geometric features that are resonant with gyro-orbits, or when the gyro orbit trajectories load components that are not designed to tolerate heat fluxes. It is up to the HEAT user to intelligently determine when the gyro-orbit calculation can be omitted, and when it is mandatory.

Full tile analysis
The previous section outlined some of the gyro-orbit effects that arise for a single castellation. While this is useful to ascertain the temperatures and fluxes over a single castellation, it neglects the macroscopic effects of gyro-orbit heat fluxes on an entire PFC. It also neglects the full power of HEAT, which can analyze heat loading on entire 3D CAD geometries. This section is designed to accompany the previous section, and describe the macroscopic effects of the ion gyro-orbit heat fluxes on entire PFCs. The NSTX-U IBDH tile has been chosen for this example, as it will see some of the highest heat fluxes of the entire machine. Section 1.3 provides an introduction to the NSTX-U PFCs and figure 1 provides a visualization of the IBDH tile and where it resides in the NSTX-U divertor. Figure 12 provides a visual comparison between the optical approximation and two gyro-orbit predictions for the entire IBDH tile in the lower NSTX-U divertor. HEAT requires the user to input the average ion temperature (eV) at the PFCs (see equation (13)), and uses this temperature to generate a Maxwellian speed distribution as described in section 2. Changing the input temperature will change the gyro-radii of particles, and increasing the temperature increases the average distribution of gyro-radii. This dependence of gyro-radii on temperature is important when calculating the gyro-orbit heat flux profiles, as the gyro-radius of an ion determines the amount it can 'dip' into an optical magnetic shadow. Heat fluxes for the optical approximation (left) and gyro-orbit approximation (center and right) for the entire NSTX-U IBDH tile in the lower divertor. With the optical approximation, shadows have sharp edges and cover a large fraction of the PFC. One large inter-PFC shadow can be observed on the left edge of the tile, and small inter-castellation shadows can be observed on the left edge of each castellation. Turning on the gyro-orbit approximation with average ion temperature set to 10 eV fills in the inter-castellation shadows, but has a minimal effect on the inter-PFC shadows. Increasing the average ion temperature to 100 eV increases the gyro-radii of ions, which results in total smearing of the optical approximation shadows. HEAT can perform gyro-orbit analyses on entire PFC geometries directly from engineering CAD. Figure 12 shows the heat fluxes corresponding to the optical approximation and two plasma ion input temperatures.
The optical approximation (left frame) in figure 12 has a heat flux footprint that is characterized by alternating areas of heat flux and shadows, which correspond to the geometric features of the IBDH tile. The narrow bands of shadows that appear at the left side of every column of castellations is due to the inter-castellation fish-scale. The inter-castellation fish-scale is small, so these shadows are narrow. The wide shadow at the left side of the IBDH tile is a result of the inter-PFC fish-scale, which is designed to protect the leading edges on adjacent IBDH tiles. As discussed in the preceding sections, the helical trajectories of ions when the gyro-orbits are accounted for enables them to access optically shadowed regions. This can be observed in the center and right frames in figure 12, which correspond to average plasma temperatures of 10 eV and 100 eV, respectively.
For the 10 eV case (center frame), the stark shadow lines of the optical approximation begin to be smeared out by ions with finite Larmor radii. The inter-castellation fish-scale shadows are approximately halfway filled in, while the inter-PFC fish-scale shadows are approximately one quarter filled. At 100 eV average ion temperature (right frame), the intercastellation fish-scale shadows have completely disappeared, and the majority of the inter-PFC fish-scale shadow has disappeared. The distribution of ion gyro-radii at 100 eV smear out almost all of the optical shadows on the PFC, which increases the PFC plasma wetted area. This increase in plasma wetted area serves to enhance the PFC performance, as the same amount of spatially integrated power is distributed over a larger surface area. Figure 12 does not focus on the small gyroorbit effects on the edges (as was performed in the previous sections), but rather on the macroscopic PFC performance. For the macroscopic PFC analysis, the helical trajectories of the ions serves to spread the heat flux across a surface area that is inaccessible via the optical trajectory.

Conclusion and future work
The helical trajectories of ions as they gyrate about the magnetic field lines at finite Larmor radii can result in 3D power exhaust profiles that load geometrical features non-uniformly. Helical trajectories that resonate with specific geometric features, such as edges and corners, can create localized hot spots. Simultaneously, the helical trajectories can reduce the average heat flux on a surface by spreading the power into regions that are magnetically shadowed. These competing effects require advanced simulation frameworks that are capable of analyzing the interaction between particle gyro-orbit trajectories and real engineering CAD geometry. HEAT provides such a framework, and has been utilized to generate power exhaust predictions under the optical and gyro-orbit approximations for NSTX-U.
The HEAT algorithms and mathematical formulations were described, and then several example HEAT outputs were presented. A comparison to existing research on ion gyro-orbit heat fluxes was performed. Convergence studies illustrate that the HEAT gyro-orbit algorithms are valid and theoretically sound, though an experimental validation is still necessary. The IBDH tile in the NSTX-U divertor was analyzed under both the optical approximation and the ion gyro-orbit approximation, and the results were compared. There are generally two effects that are characteristic of the ion gyro-orbit heat fluxes for the NSTX-U castellated PFCs: hot corners and shadow filling.
The corners of the NSTX-U divertor PFCs will be loaded by more intense heat fluxes when the gyro-orbits are included. For NSTX-U, these corners are loaded with heat fluxes of ≈25 MW m −2 , but only on a region of ≈200 μm in width, resulting in a relatively low spatially integrated deposited energy. These edges could potentially represent a problematic situation if the bulk PFC material cannot readily diffuse away the deposited energy, or if the width of these heatloaded regions expands. Ironically, the optical approximation hot spot that appears on the NSTX-U downstream corner as a result of faceting drives the PFC temperature faster than the gyro-orbit leading edge, despite its much lower peak heat flux.
Additionally, some regions of the PFC that are magnetically shadowed will see gyro-orbit heat fluxes, because the helical trajectories of the ions enables particles to 'dip' into the magnetic shadows and strike PFC surfaces. Except for the aforementioned cases when the gyrating particle loads a corner or edge, this effect enhances PFC performance. The performance enhancement is a result of the increased plasma wetted area, which decreases the peak heat flux across the PFC surface. Both the hot corners and the shadow filling must be considered when designing PFCs for high power tokamaks.
Future work will aim to extend this analysis to other regions of the NSTX-U divertor, as well as to limiters. Some work has already been performed on center stack limited discharges, but the extremely shallow (i.e. <1 • ) angles of incidence near the separatrix contact point require a significant number of macro-particles to cover an entire PFC tile. It should also be noted that for limited discharges, there will likely be heat fluxes on the limiter PFCs due to particles confined inside the separatrix that have large (i.e. up to 10 mm) gyro-radii. A true accounting of the gyro-orbit heat fluxes during limited operation would require including a model for these confined particles. Finally, near the separatrix contact point in a limited discharge the electric field likely becomes non-negligible due to the shallow field line angles. It is for these reasons that the simulations to date for limited discharges are considered only partially valid, and are not presented in this manuscript.
To accelerate the calculation and make higher fidelity simulations possible, the gyro-orbit heat flux profile could be generated by other mathematical tricks. One such method is by applying a 3D convolution kernel to the optical heat flux profile, which 'smears' the optical profile and amplifies it on the edges. Another method is to map the PFCs directly into flux space and assign gyro-orbit power using flux-coordinate versions of the distribution functions. A third idea is to perform a frequency domain analysis on the PFCs to find the overlapping resonances between the ion helical trajectories, the magnetic field line trajectories, and the geometry. Any of these methods would enable the gyro-orbit power to be mapped directly to the CAD mesh, eliminating the computationally expensive helical trajectory tracing and associated mesh-triangle intersection calculations, which consume the bulk of the simulation time.