Development of Bounce-Time-Based Orbit-Following Monte-Carlo Code

Abstract

We developed a bounce-time (BT)-based orbit-following Monte-Carlo code as an extension of the OFMC code in QST. In the BT-based method, we take a bounce time as a time step of the orbit following. The time step is ~ 100 times longer than the gyro period which is a typical time step for the conventional guiding-center (GC) method. In the BT-based method, an accurate and simple estimation of a poloidal projection of the bounce orbit and a staying time are essential. An expression for the orbit gives us an orbit shape by a small calculation with the difference of less than 1% of the minor radius, compared with the GC method with the same fast ion parameters. And an approximate expression for the staying time also gives us the staying time with a good accuracy for our purpose. We can see a good agreement between calculation results for the BT-based method and those for the GC method in an axisymmetric condition. The BT-based method is 70–140 times faster than the GC method, depending on the slowing-down time.

1.  Introduction

Heating, current drive, and torque input distribution by fast ions such as alpha particles or generated by neutral particle beams (NB) are important for an operational scenario design of a fusion reactor and a quantitative interpretation of experimental data. At QST, we have determined the fast ion distribution required for these distributions by tracking the trajectories of many test particles or marker particles using a Monte-Carlo numerical calculation code called OFMC [1]. In the OFMC, the beam source position and direction of the NB are arranged to match the actual NB injector, so that the initial distribution of birth fast ions, including the pitch angle, can simulate the actual situation or experimental condition. In addition, the code can easily handle realistic magnetic perturbation and first-wall shape. Due to this implementation, the code well reproduced the experimental results such as wall heat-load in JT-60U [2, 3], JFT-2M [4], JET and DIII-D [5] tokamaks in various conditions and contributed to ITER design activities [6–9].

In the many standard orbit-following codes, when tracking particle trajectories, a guiding center (GC) trajectory calculation which uses a gyro period as a typical time step unit is used instead of tracking “particle” trajectories in order to shorten calculation time. The guiding center calculation, or a drift-orbit calculation, can be applied to a condition with electric and magnetic field variations where the gyro radius can be considered closed during the gyro cycle. This condition is well met in most magnetic confinement devices, even in a low field-strength tokamak such as MAST [10].

In an operation scenario design, it is necessary to try many NB injection patterns, thus numerical calculations with shorter calculation time or lower calculation load are desired. In addition, a local velocity distribution function in phase space is required when trying to simulate fast ion diagnostics such as neutron distribution measurement, or when discussing velocity space instability such as an ion cyclotron emission (ICE). Determining the local velocity distribution function using the test-particle method requires a large number of test particles, which requires a large computational resource, in order to reduce a Mote-Carlo noise. A computational method with a low computational load is desired from these backgrounds.

We consider the calculation-time reduction using the same idea as in introducing guide-center calculation. A particle in an axis-symmetric tokamak has a closed orbit projection on a poloidal cross-section, characterized by a toroidal canonical angular momentum. The time period of the closed orbit is a bounce time. The axis-symmetry is well realized in most regions in tokamak devices. Thus, the assumption of the axis-symmetry is a reasonable assumption which we are based on. The bounce period is the next shortest period after the gyro period for a particle orbit in the tokamak. Note the bounce time is used for many different periodic motions in plasma physics. In this paper, the bounce time is defined as the period of a periodic orbit in a poloidal cross section, and is expressed as

  

$${\tau }_b={\oint }_{pol}^{}\frac{d{l}_{pol}}{v_{pol}},$$(1)

where $d{l}_{pol}$ is the line element along the trajectory of the guide center projected onto the poloidal cross section, and $v_{pol}$ is the velocity of the guiding center projected on the poloidal cross section. The integration is performed for one round of a closed periodic orbit on a poloidal cross section. This periodic orbit is a banana shape for trapped particles, or a shape which is like a shifted poloidal flux surface in the radial direction for passing particles (Fig. 1). This bounce period is more than 100 times longer than the period of the gyromotion. If we can define this closed orbit, we can utilize this bounce period for a time step of an “orbit-following” calculation. Then, it is expected that we can reduce the calculation time.

Fig. 1.  A poloidal projection of an orbit for a co-passing particle (a), for a counter-passing particle (b), and for a trapped particle (c). The particle species is deuteron. The points obtained by Eq. (10) are shown in red closed circles. The blue curve represents an orbit obtained by the GC calculation. The first “wall” for the calculation in the “Results” section is indicated by an arrow in the panel (b). The R_h, R_l, R_ctr, and R_co are the major radial positions which appear in Fig. 2. Each of them corresponds to the major radial position of an orbit on the midplane.

Here the word “orbit-following” should be understood in the sense of following the collisional motion of test particles that represent entire GC drift orbits in the reduced 3-D phase space of the ($P_{\phi }$, $E$, $\mathrm{\Lambda }$), which comprise the complete set of constants of motion (CoM) in an axisymmetric system, with $P_{\phi }$ being the toroidal canonical angular momentum, $\mathrm{\Lambda }$ the pitch parameter, $E$ the kinetic energy of the particle.

This idea itself is not new and has been used in a framework of a bounce-averaged Fokker-Plank approximation, in which, a Fokker-Plank equation is numerically solved. Our approach differs from the usual Fokker-Planck solvers in two ways. First, instead of evolving the distribution function directly, we use the framework of the OFMC code; that is, the Monte-Carlo method. Second, instead of entirely eliminating the ($R$, $z$)-dependence by bounce-averaging the equations of motion, we retain it in the form of an analytical approximation of the orbit contour and parallel motion. Thus, our approach can be seen as option to analytically short-cut the conventional orbit-following calculation in OFMC. It is thus straightforward to switch between the conventional orbit-following solver and the accelerated bounce-time (BT)-method when necessary, and they may even be run side-by-side for different regions of phase space (for instance, in the case when the analytical approximations are not sufficiently accurate for certain orbits as shown later). In addition, our approach fully accounts for the magnetic drifts known as finite orbit width (FOW) effect, since we utilize the orbit shape projection on the poloidal cross-section, while conventional Fokker-Plank codes ignore the FOW effect or use the model in which an orbit deviation from a poloidal flux surface is assumed to be small. (Note some Fokker-Plank codes with the capability of fully taking into account the FOW have been developed; e.g. ATEP-3D code [11, 12].) However, the orbit deviation is finite. It reaches 1/4 of the minor radius in some cases and the banana width reaches 1/2 of the minor radius. and it can be about 5% of the minor radius for the alpha particles in the ITER for the parameter in Ref. [13]. The orbit deviation or the finite orbit width affects a heating and current drive profile. Our method can handle the orbit width effect.

This paper is organized as follows. We describe the method and its implementation in Sec. 2. The key of the implementation is how we well and simply evaluate an orbit shape. We describe the method and show a comparison between the evaluated orbit and the orbit calculated by the motion equation of the GC. In the connection with a slowdown process and a heating distribution, we also describe the implementation on a collision between a bulk plasma and a fast ion. In Sec. 3, the results are compared with the GC calculation. The summary is given in Sec. 4.

2.  Method and Implementation

2.1  Abstract of fast ion transport and its interaction with bulk plasma

Here, we describe the abstract of the implementation. In an axisymmetric system or an ideal tokamak, a fast ion orbit is characterized by the labels of ($P_{\phi }$, $E$, $\mathrm{\Lambda }$, $sign\left(v_{\parallel }\right)$), where $sign\left(v_{\parallel }\right)$ is the sign of the velocity of the particle against the plasma current. A set of the labels represents a closed GC orbit on a poloidal cross-section. The summation of orbits gives the fast ion distribution in real space. The orbit is not changed unless the set of the labels is changed; that is, the simulation particle does not “move” in the 3-D CoM space. The mechanism to change the labels or to move the particle is a collision. The collision is mainly the interaction with bulk plasma. The collision introduces heating of the bulk plasma. Since we take into account the finite-orbit effect or know the orbit shape for a set of ($P_{\phi }$, $E$, $\mathrm{\Lambda }$, $sign\left(v_{\parallel }\right)$) in the real space as described below, we can evaluate a heating distribution accurately, compared with conventional Fokker-Plank solvers. The collision also changes the quantities of ($P_{\phi }$, $E$, $\mathrm{\Lambda }$, $sign\left(v_{\parallel }\right)$), namely transports the fast ion in a phase space, depending on the plasma parameters at the collision. In this way, how the orbit in real space is well and simply evaluated is the key of this approach. In the following, we show how we evaluate the orbit and show a comparison between the evaluated orbit and the orbit calculated by the motion equation of the GC.

2.2  Description of orbit in the real space

The key is to evaluate the orbit characterized by ($P_{\phi }$, $E$, $\mathrm{\Lambda }$, $sign\left(v_{\parallel }\right)$) as accurately and simply as possible.

The definition of these values are

  

$$P_{\phi }=q{\psi }_p-MR{v}_{\phi },$$(2)

  

$$E=\frac{1}{2}M{v}^2=\frac{1}{2}M(v_{\parallel }^2+v_{\perp }^2),$$(3)

  

$$\mathrm{\Lambda }=\frac{\mu B_0}{E}.$$(4)

Here, $q$ is the electric charge of a fast ion, ${\psi }_p$ is the poloidal flux function, $M$ is the mass of the fast ion, $R$ is the major radial position at the fast ion position, $B_0$ is the magnetic field strength at the magnetic axis, and $v_{\phi }$ is the velocity of the fast ion in the toroidal direction, respectively. The $v_{\parallel }$ and $v_{\perp }$ are the velocity of the fast ion in the direction parallel and perpendicular to the magnetic field, respectively. The $\mu$ is the magnetic moment defined as $\mu \equiv M{v}_{\perp }^2\mathrm{\text{/}}2B$. We use

  

$${\hat{P}}_{\phi }\equiv P_{\phi }\mathrm{\text{/}}q={\psi }_p-\frac{M}{q}R{v}_{\phi },$$(5)

instead of $P_{\phi }$. This form easily relates to the poloidal flux function ${\psi }_p$ and the second term represents the finite orbit effect or the orbit deviation from a flux surface. (This form ${\hat{P}}_{\phi }$ is also known as ${\psi }_p^{*}$.) We assume the toroidal field strength is inversely proportional to $R$ or $B_{\phi }R=B_0{R}_0$ = const., ignoring the diamagnetic effect and we also assume $B_{\phi }\sim B$, where $R_0$ is the major radius value at the magnetic axis, $B$ is the total magnetic field strength. Then, we approximate ${\hat{P}}_{\phi }$ as follows;

  

$$\begin{array}{cc}{\hat{P}}_{\phi }& ={\psi }_p-\frac{M}{q}R\frac{B_{\phi }}{B}{v}_{\parallel }\hfill \\ & \sim {\psi }_p-\frac{M}{q}R{v}_{\parallel }.\hfill \end{array}$$(6)

We can also have

  

$$\begin{array}{cc}E& =\frac{1}{2}M{v}_{\parallel }^2+\mu B\hfill \\ & \sim \frac{1}{2}M{v}_{\parallel }^2+\mu \frac{B_0{R}_0}{R}.\hfill \end{array}$$(7)

Now we define $R_{tip}$ as

  

$$R_{tip}\equiv \mathrm{\Lambda }{R}_0.$$(8)

In the case of a trapped particle, this $R_{tip}$ represents the major radius position of the banana tip because the $v_{\parallel }=0$ at the banana tip of $R=R_{tip}$ in Eq. (7).

By using Eqs. (7) and (8), we can write the parallel velocity at the position of $R$ as

  

$$v_{\parallel }=\mathrm{\text{sign}}\left(v_{\parallel }\right)\sqrt{2E\mathrm{\text{/}}M}\sqrt{1-R_{tip}\mathrm{\text{/}}R}.$$(9)

Then, Eq. (6) leads to

  

$${\psi }_p={\hat{P}}_{\phi }+\mathrm{\text{sign}}\left(v_{\parallel }\right)\frac{\sqrt{2EM}}{q}\sqrt{R(R-R_{tip})}.$$(10)

From this equation, we can have the poloidal projection of an orbit for a given set of ($P_{\phi }$, $E$, $\mathrm{\Lambda }$, $sign\left(v_{\parallel }\right)$). When we specify a radial position $R$, the right-hand-side of the Eq. (10) is determined and we can specify a poloidal flux. Then we can determine two vertical positions using a mapping of ${\psi }_p\leftrightarrow (R,Z)$ since a specific poloidal flux passes two vertical position at a specific radial position inside a separatrix.

In Fig. 1, we compare the poloidal projection of a orbit between the one obtained by the Eq. (10) and the one obtained by the GC calculation for an equilibrium of E039672 in the JT-60U, with $B_0$/ Ip/ ${\psi }_{p0}$ = 1.2T/0.6MA/−0.26, where Ip and ${\psi }_{p0}$ are the plasma current and the value of the poloidal flux function at the magnetic axis, respectively. We define the plasma surface with ${\psi }_p=0$. The direction of the plasma current and the toroidal field are both in the clock-wise direction from the top of the torus. Figure 1(a) is the case of (${\hat{P}}_{\phi }$, $E$, $\mathrm{\Lambda }$, $sign\left(v_{\parallel }\right)$) = (−0.166, 85.9 keV, 0.77, +). This is a co-passing orbit, where co- means the particle is going in the plasma current direction. The points obtained by Eq. (10) are shown in red closed circles. And the blue curve is obtained by the GC calculation. Figure 1(b) is the case of (${\hat{P}}_{\phi }$, $E$, $\mathrm{\Lambda }$, $sign\left(v_{\parallel }\right)$) = (0.035, 83.5 keV, 0.147, −). This is a counter (ctr)-passing orbit. Figure 1(c) is the case of (${\hat{P}}_{\phi }$, $E$, $\mathrm{\Lambda }$) = (−0.222, 85.9 keV, 0.993), which is a trapped orbit. We can see a good agreement. The normalized difference $\mathrm{\Delta }r/a$ to the minor radius is less than 1% to GC. This is good enough for our purpose. Note that the dots in Fig. 1 do not represent the time steps of an “orbit following” in real space. The number of dots just indicates the spatial resolution for the evaluation of distributions shown in Fig. 8 or collision point which affects the heating distributions shown in Fig. 9.

2.3  Collision

The mechanism to change the orbit is a collision in an axisymmetric system. At the collision point, the energy of a fast ion is transferred to an electron or ion.

We use the Trubnikov model to describe the collision [1, 14]. In the model, the effect of a collision, which is a change in velocity, depends on a specified interval time between collisions as follows. In the model, the change of the velocity in the direction to the original movement or the longitudinal direction, $\mathrm{\Delta }{v}_L$, consists of the averaged change, $\mathrm{\Delta }{\bar{v}}_L$, and the statistical deviation from the averaged change, $\sqrt{\mathrm{\Delta }{\overline{v_L}}^2}$. The average change $\mathrm{\Delta }{\bar{v}}_L$ is proportional to the interval time between collisions while the deviation $\sqrt{\mathrm{\Delta }{\overline{v_L}}^2}$ is proportional to the square root of the collision interval time. The velocity change in the transverse direction to the original movement, $\mathrm{\Delta }{v}_T$, consists of two components in the same manner as $\mathrm{\Delta }{v}_L$, however the average component is always zero. Thus, the changes in both longitudinal and transverse directions get larger as the collision interval is longer.

The orbit is characterized by ($P_{\phi }$, $E$, $\mathrm{\Lambda }$, $sign\left(v_{\parallel }\right)$). A collision changes these quantities. In order reasonably to describe a transport of fast ions, it is preferable that the variations ($\mathrm{\Delta }{P}_{\phi }/{\psi }_{p0}$, $\mathrm{\Delta }E$, $\mathrm{\Delta }\mathrm{\Lambda }$) are small at “a collision” since the variations are the step-size of the evolution in the phase space. For the collision interval of ${\tau }_b$, the variations are $\mathrm{\Delta }{P}_{\phi }/{\psi }_{p0}$ ~ 1 × 10⁻⁴, $\mathrm{\Delta }E$ ~ 4 × 10⁻² keV, $\mathrm{\Delta }\mathrm{\Lambda }$ ~ 3 × 10⁻⁴ for the beam energy of $40\times T_{e0}$ or 84 keV, $\mathrm{\Delta }{P}_{\phi }/{\psi }_{p0}$ ~ 7 × 10⁻⁴, $\mathrm{\Delta }E$ ~ 8 × 10⁻² keV, $\mathrm{\Delta }\mathrm{\Lambda }$ ~ 1 × 10⁻³ for the beam energy of $10\times T_{e0}$, and $\mathrm{\Delta }{P}_{\phi }/{\psi }_{p0}$ ~ 6 × 10⁻⁴, $\mathrm{\Delta }E$ ~ 2 × 10⁻² keV, $\mathrm{\Delta }\mathrm{\Lambda }$ ~ 3 × 10⁻³ for the beam energy of $5\times T_{e0}$ for the plasma parameters in the Sec. 3, where $T_{e0}$ is the electron temperature at the magnetic axis. Thus, the variations are reasonably small for the evaluation of the transport, though the resolution gets worse for lower energy. Namely, it is reasonable to choose the bounce period, ${\tau }_b$, as a collision interval. The usage of a long collision-interval of this order was also validated by the GC method. The fast ion distribution for the collision-interval of 5 × 10⁻⁴ of the slowing down time, which is comparable to a typical bounce time, was almost identical to those for the collision-interval of 1 × 10⁻⁶ of the slowing down time in the GC calculations.

The collision position, $R_{col}$, is randomly determined with taking into account the staying time as shown below. The plasma parameters for a collision are given in a coordinate of a poloidal flux. The poloidal flux is determined at $R=R_{col}$ by using Eq. (10). Here, this must be done at each time step, because we follow MC marker particles in CoM space, so they may travel into regions of varying collisionality. In the case of ATEP-3D code [11, 12], which operates on a fixed CoM mesh, the collision operator needs to be computed for each mesh point only once in a stationary system (and may be updated occasionally when the background evolves slowly). Currently, ATEP-3D performs the average procedure along numerically computed GC orbits.

2.4  Staying time and bounce time

When we choose a collision point in a specific orbit, it is reasonable to give a higher possibility to a point where a staying time or traversal time of a particle is longer. The staying time is not constant at each position in an orbit since the velocity projected on a poloidal cross-section depends on its position. It is an important factor to evaluate the staying time as accurately and simply as possible. We use the following approximation for the staying time $\mathrm{\Delta }t$ at the major radius $R$.

  

$$\begin{array}{cc}\mathrm{\Delta }t& =\left|\frac{\mathrm{\Delta }{l}_{pol}}{v_{pol}}\right|=\left|\frac{\mathrm{\Delta }l}{v_{GC}}\right|\sim \left|\frac{\mathrm{\Delta }{l}_{\parallel }}{v_{\parallel }}\right|\hfill \\ & =\left|\frac{1}{v_{\parallel }}\frac{B}{B_R}\mathrm{\Delta }R\right|,\hfill \end{array}$$(11)

where $\mathrm{\Delta }{l}_{pol}$ is a line segment of a GC orbit projected on a poloidal cross-section, $v_{pol}$ is the velocity of the GC motion projected on the poloidal cross-section, $\mathrm{\Delta }l$ is a line segment along the GC orbit, $\mathrm{\Delta }{l}_{\parallel }$ is its parallel component to the magnetic field, $B_R$ is the radial component of the magnetic field strength. When $R$ is given, we determine ${\psi }_p$ from Eq. (10). Then, we can calculate all the components of the right-hand-side; $v_{\parallel }\left(R\right)$ using Eq. (9), $B(R,{\psi }_p)$, and $B_R(R,{\psi }_p)$.

In the code, we use a bounce time, ${\tau }_b$, as an interval time for a collision. The bounce time is evaluated by the following expression;

  

$$\begin{array}{cc}{\tau }_b& ={\oint }_{pol}^{}\frac{d{l}_{pol}}{v_{pol}}\hfill \\ & ={\int }_0^{{\tau }_b}dt\hfill \\ & \sim {\int }_{R_l}^{R_h}{\left|\frac{1}{v_{\parallel }}\frac{B}{B_R}dR\right|}_{\mathrm{\text{upper half}}}\hfill \\ & +{\int }_{R_h}^{R_l}{\left|\frac{1}{v_{\parallel }}\frac{B}{B_R}dR\right|}_{\mathrm{\text{lower half}}},\hfill \end{array}$$(12)

for a passing particle, where $R_h$ and $R_l$ are the smallest and the largest radial positions in a specific orbit, respectively as shown in Fig. 1. The h and l denote the high- and the low- field-sides. The first term in the last line is the integration for the upper-half part of the orbit and the second term is that for the lower-half part of the orbit. For a trapped particle,

  

$$\begin{array}{cc}{\tau }_b& \sim {\int }_{R_{co}}^{R_{tip}}{\left|\frac{1}{v_{\parallel }}\frac{B}{B_R}dR\right|}_{\mathrm{\text{upper half}}}\hfill \\ & +{\int }_{R_{tip}}^{R_{ctr}}{\left|\frac{1}{v_{\parallel }}\frac{B}{B_R}dR\right|}_{\mathrm{\text{upper half}}}\hfill \\ & +{\int }_{R_{ctr}}^{R_{tip}}{\left|\frac{1}{v_{\parallel }}\frac{B}{B_R}dR\right|}_{\mathrm{\text{lower half}}}\hfill \\ & +{\int }_{R_{tip}}^{R_{co}}{\left|\frac{1}{v_{\parallel }}\frac{B}{B_R}dR\right|}_{\mathrm{\text{lower half}}},\hfill \end{array}$$(13)

where $R_{co}$ is the largest radial position in a co-passing movement to the plasma current, and $R_{ctr}$ the largest radial position in a counter(ctr)-passing movement as shown in Fig. 1(c). The first and second terms are the integrations for the upper-half part of the orbit and the third and fourth terms are those for the lower-half part of the orbit.

To calculate the bounce time formulated in the above, we need to know the values of the $R_h$, $R_l$, $R_{co}$, and $R_{ctr}$ in addition to the $R_{tip}$. These four values are obtained by evaluating the intersection of two curves shown in Fig. 2. The blue-solid curve in Fig. 2 is the poloidal flux function, ${\psi }_p$, on the mid-plane. The other red curve is obtained by the Eq. (10), which represents the relation between a poloidal flux and the radial coordinate for an orbit.

Fig. 2.  Plots which show that the positions of $R_h$, $R_l$, $R_{ctr}$, and $R_{co}$. These positions are the intersection of two curves. The red curve is given by Eq. (10). The blue curve represents a poloidal flux on the mid-plane. (a) is for a co-passing particle, (b) for a counter-passing particle, and (c) for a trapped particle (c).

In the panels (a) in Figs. 3–5, we compare the staying time, $\mathrm{\Delta }t$, between the one obtained by Eq. (12) or Eq. (13) and the one obtained by the GC calculation for the same equilibrium as Fig. 1. In addition, the time coordinate along an orbit, which is defined in

Fig. 3.  The comparison of the staying time (a) and time coordinate (b) at a position along an orbit. The horizontal axis is described in the normalized poloidal flux, ${\psi }_{pN}$. The red dots represent the one by Eq. (11) or (14), and the blue dots represent the one obtained by a GC calculation. The comparison is carried out for the upper-half part of an orbit. The parameters of ($P_{\phi }$, $E$, $\mathrm{\Lambda }$, $sign\left(v_{\parallel }\right)$) which characterizes the orbit are those in Fig. 1(a). The orbit type is the co-passing orbit.

Fig. 4.  The comparison of the staying time (a) and time coordinate (b) for the parameters in Fig. 1(b). The orbit type is the ctr-passing orbit.

Fig. 5.  The comparison of the staying time (a) and time coordinate (b) for the parameters in Fig. 1(c). The orbit type is the banana orbit.

  

$$t\left(l\right)={\int }_{}^{}\frac{d{l}_{\parallel }}{v_{\parallel }}={\int }_{}^{}\left|\frac{1}{v_{\parallel }}\frac{B}{B_R}dR\right|,$$(14)

is also compared in the panels (b) in these figures. Figure 3 is the case of the co-passing orbit with ($P_{\phi }$, $E$, $\mathrm{\Lambda }$, $sign\left(v_{\parallel }\right)$) in Fig. 1(a). Figure 4 is the case of the ctr-passing orbit in Fig. 1(b). And Fig. 5 is the case of the trapped orbit in Fig. 1(c). We can see a good agreement in the staying time and time coordinate. This is good enough for our purpose. The coordinate, $l$, is resented by the normalized poloidal flux, ${\psi }_{pN}$, in the figures. (The relation between ${\psi }_{pN}$ and ${\psi }_p$ is ${\psi }_{pN}\equiv 1-{\psi }_p/{\psi }_{p0}$.) For this coordinate, the integration region in Eq. (14) is from $R_l$ to $R$ between $R_h$ and $R_l$ in Fig. 3, is from $R_h$ to $R$ between $R_h$ and $R_l$ in Fig. 4, and is from $R_{co}$ to $R$ between $R_{co}$ and $R_{tip}$ and then between $R_{tip}$ and $R_{ctr}$ for the upper-half part of the orbit. The relation between $R$ and ${\psi }_{pN}$ can be obtained through Eq. (10). For the full bounce motion, it is also evaluated in the lower-half part of the orbit in the same way.

2.5  Loss orbit

A particle cannot travel when the orbit intersects a loss boundary. The loss boundary is a first wall structure in a realistic case. However, a complicated procedure is required to know whether the orbit intersects the realistic first wall. Here, we set a loss boundary at the plasma surface or separatrix. We consider a particle is lost when ${\psi }_p$ is larger than 0 since the separatrix lies at ${\psi }_p=0$ in our calculation. To implement a realistic loss boundary is a future work.

2.6  Evaluation of fast ion density and driven current

Here, we describe how to evaluate profiles related to fast ions. Each orbit in the real space is divided into many segments. The “staying time”, $\mathrm{\Delta }t$, at a segment is multiplied with the weight, $w_i$, of the test particle. In this calculation, the weight is the same for all the particles. Instead, the number of particles are proportional to the value of$P_{in{j}_k}/E_{in{j}_k}$, where $P_{in{j}_k}$is the injected beam power, $E_{in{j}_k}$ the injected energy from k^th injector. The value $w_i\mathrm{\Delta }t$ is summed up at the position of ${\psi }_p$ where the segment lies. In this way, the fast ion density profile is obtained. The fast ion energy density is obtained when the particle energy is multiplied in addition at the segment.

To evaluate the NB current drive (NBCD) profile, the value of ${⟨j_{\parallel }B⟩}_{{\psi }_p}$ is used [15, 16], where ${⟨X⟩}_{{\psi }_p}$ is the flux averaged value of $X$. The value is approximately obtained as follows;

  

$$\begin{array}{c}⟨j_{\parallel }B{⟩}_{{\psi }_p}\hfill \\ =\left(q\mathrm{\text{/}}\mathrm{\Delta }{V}_{{\psi }_p}\right)\sum _{i,k_i}^{}\left(v_{\parallel }B\right)\left(w_i\mathrm{\Delta }t\right){|}_{\mathrm{\text{for a segment with}}{\psi }_p\mathrm{\text{in}}{k}_i^{th}\mathrm{\text{orbit}}}\hfill \\ =\left(q\mathrm{\text{/}}\mathrm{\Delta }{V}_{{\psi }_p}\right)\sum _{i,k_i}^{}{w}_i\frac{B^2}{B_R}\mathrm{\Delta }R{|}_{\mathrm{\text{for a segment with}}{\psi }_p\mathrm{\text{in}}{k}_i^{th}\mathrm{\text{orbit}}},\hfill \end{array}$$(15)

where $i$ represents the index for a test particle, $w_i$ the weight of the test particle, $\mathrm{\Delta }{V}_{{\psi }_p}$ is the volume of a shell specified by a ${\psi }_p$. Namely, the value $q{w}_i\left(B^2/B_R\right)\mathrm{\Delta }R$ in a segment with ${\psi }_p$ in the $k_i^{\mathit{th}}$ bounce orbit is summed over all test particles till the slowdown. The NBCD profile is obtained by multiplying $⟨j_{\parallel }B{⟩}_{{\psi }_p}$ with a geometrical coefficient and a coefficient for an electron shielding [15, 16].

3.  Comparison between BT and GC Method

Here, we show results of the BT method by comparing with those of the GC method in an axisymmetric condition.

The BT method is implemented as a module in the framework of the OFMC code. Thus, we start the calculation with the same birth profile of beam ions.

The plasma parameters are as follows. The equilibrium is shown in Fig. 1. The profiles of electron density, electron temperature, and ion temperature are depicted in Fig. 6. The effective charge number, Z_eff, is uniform and its value is 2.1. The impurity species is carbon alone. The beams are three deuteron beams of ~ 85 keV. The birth profile of the induced ions are shown in Fig. 7. One beam is injected in the co-parallel direction to the plasma current, and one is injected in the counter-parallel direction to the plasma current, and the last beam is injected nearly into the perpendicular direction to the current. The number of test particles is 20,000.

Fig. 6.  Profiles of electron density (a), electron temperature (b), and ion temperature (c).

Fig. 7.  The birth distribution of NB injected ions. It is viewed from the top of the torus. The direction of plasma current and toroidal field is the clock-wise direction.

The results by the BT method are compared with those by the GC method in Fig. 8. Profiles of fast ion energy density, $⟨j_{\parallel }B{⟩}_{{\psi }_p}$, and NBCD are compared in Figs. 8(a), (b), and (c), respectively. In each panel, a black solid curve is obtained by the GC method and three non-solid curves (dotted, dashed, and dash-dot) are obtained by the BT method. The different results for the BT method are obtained with different seeds for random numbers for the Monte-Carlo calculation. Thus, the scatter indicates the Monte-Carlo noise level. (We also observed the same level of the Monte-Carlo noise in the GC method.) We can see a good agreement in the level of Monte-Carlo noise except for a difference of < 20% near magnetic axis in the fast ion energy density. The possible cause of the difference will be discussed below.

Fig. 8.  The comparison of profiles of the fast ion energy density (a), $⟨j_{\parallel }B{⟩}_{{\psi }_p}$ (b), neutral beam current density (c). The solid curve obtained by the GC method and three non-solid curves are obtained by the BT method.

The profiles of electron, bulk deuteron, and impurity (carbon) heating are also compared in Fig. 9. We can also see a good agreement except for the region near magnetic axis.

Fig. 9.  The comparison of profiles of deposit power to bulk deuterons (a), deposit power to bulk electrons (b), deposit power to impurity carbon ions (c). The solid curve obtained by the GC method and three non-solid curves are obtained by the BT method.

The calculation time has gotten shorter by introducing the BT-based method. The BT method is more than 70 times faster than the GC method, or 85 minutes for the GC method and 1.2 minutes for the BT method in our environment. The speed up ratio depends on the slowing-down time. The ratio was ~ 70 times for the slowing-down time of 0.3s and it reached ~ 140 times for the slowing-down time of 1.4s.

Here, we shortly discuss the difference observed in the fast ion energy density in Fig. 8(a). The possible cause could be the small difference of the orbit shape of “shallowly” trapped particles. One of them is shown in Fig. 10. We accept a small difference in the orbit shape which is $\mathrm{\Delta }r/a<0.01$ as described in the above, however its relative difference of ~ 10%, shown by the shaded box in the Fig. 10(c), in an inner leg part of the trapped orbit can affect the population near magnetic axis, where such inner legs of “shallowly” trapped particles can lie. Here, the inner leg indicates the part of a trapped orbit which lies on the side of the magnetic axis or has its major radius between $R_{tip}$ and $R_{ctr}$. This difference is induced by the approximation used in Eq. (10). Because of the approximation of $B_{\phi }\sim B$, the second term on the right-hand-side in Eq. (10) is affected. Then, the red curve in Fig. 2(c) by Eq. (10) tends to be narrower in the vertical direction than the corresponding curve evaluated by the GC method. This results in the situation that the inner leg tends to be longer in our method.

Fig. 10.  A plot for a “shallowly” trapped obit (a) and its expanded plot (b). The staying time is compared in the panel (c). The red makers are calculated by the BT method and the blue markers/line are calculated by the GC method. The shaded region which consists of red markers in the panel (c) is the contribution of a “larger” part of the inner leg in the BT method.

The absence of blue dots in the shaded region means that the orbit by the GC calculation is absence in the region. This leads to the fact that the orbit by the GC calculation stays for a shorter time near the magnetic axis than the orbit by Eq. (10) by the amount of the sum of $\mathrm{\Delta }t$ in the shaded region. Then, the orbit for the BT method stays for a somewhat longer time than that for the GC method near the magnetic axis. In addition, the staying time evaluation in Eq. (11) can be larger since B_R is small in the region. These could lead to the larger population in the region.

4.  Summary

We developed a bounce-time-based orbit-following Monte-Carlo code as an extension of the OFMC code in QST. In the BT-based method, we take a bounce time as a time step of the orbit-following which is the “orbit-following” in the 3D phase space of ($P_{\phi }$, $E$, $\mathrm{\Lambda }$). The time step is ~ 100 times longer than the gyro period which is the characteristic time for the GC method. In the BT-based method, an accurate and simple estimation of a poloidal projection of the bounce orbit and a staying time are essential. The simple formulation Eq. (10) gives us an orbit shape with the difference of less than 1% of the minor radius, compared with the GC method with the same fast ion parameters. And the simple formulation Eq. (11) also gives us a staying time with a good accuracy for our purpose.

We can see a good agreement between calculation results for the BT-based method and those for the GC method in an axisymmetric condition. The BT-based method is 70–140 times faster than the GC method, depending on the slowdown time. It is expected that the BT-based method will contribute to speeding-up the scenario development.

The next step is the development of the hybrid code in which the GC method is used only in given conditions.

The first target for the condition is the “shallowly” trapped orbit. The GC method will be used for “shallowly” trapped orbits. We expect a reduction of differences near magnetic axis.

The second target is the particle which intersects the separatrix. To handle the realistic first wall as the loss boundary, we will switch the orbit-following method to the GC method when the orbit intersects the separatrix of the plasma, making use of our implementation in the OFMC framework.

The third target is to handle a magnetic perturbation. We know a resonant condition is important for transport. The resonant condition is evaluated by using orbit helicity [17]. By switching the orbit-following methods depending on the orbit helicity, we would like to handle a non-axisymmetric magnetic field such as toroidal field ripple or ELM-control perturbation.

 Acknowledgements

This work is partially supported by the Grant-in-Aid for Specially Promoted Research (No. 21H04973).

References

• [1]   K.  Tani et al., J. Phys. Soc. Japan 50, 1726 (1981).
• [2]   K.  Tobita et al., Nucl. Fusion 35, 1585 (1995).
• [3]   K.  Shinohara et al., Nucl. Fusion 47, 997 (2007).
• [4]   K.  Shinohara et al., Nucl. Fusion 43, 586 (2003).
• [5]   G.J.  Kramer et al., Nucl. Fusion 53, 123018 (2013).
• [6]   K.  Shinohara et al., Fusion Eng. Des. 84, 24 (2009).
• [7]   K.  Shinohara et al., Nucl. Fusion 51, 063028 (2011).
• [8]   K.  Tani et al., Nucl. Fusion 52, 013012 (2012).
• [9]   K.  Shinohara et al., Nucl. Fusion 52, 094008 (2012).
• [10]   K.G.  McClements et al., Plasma Phys. Control. Fusion 60, 095005 (2018).
• [11]   Ph.  Lauber et al., Nucl. Fusion 64, 096010 (2024).
• [12]   G.  Meng et al., Nucl. Fusion 64, 096009 (2024).
• [13]   A.  Fasoli et al., Nucl. Fusion 48, S264 (2007).
• [14]   B.A.  Trubnikov, Rev. Plasma Phys. 1, 105 (1965).
• [15]   K.  Tani et al., J. Comp. Phys. 98, 332 (1992).
• [16]   K.  Tani et al., Nucl. Fusion 48, 085001 (2008).
• [17]   K.  Shinohara et al., Nucl. Fusion 60, 096032 (2020). The orbit helicity is referred to as orbit pitch in the paper.