Drift of ablated material after pellet injection in a tokamak

Abstract

Pellet injection is used for fuelling and controlling discharges in tokamaks, and it is foreseen in ITER. During pellet injection, a movement of the ablated material towards the low-field side (or outward major radius direction) occurs because of the inhomogeneity of the magnetic field. Due to the complexity of the theoretical models, computer codes developed to simulate the cross-field drift are computationally expensive. Here, we present a one-dimensional semi-analytical model for the radial displacement of ablated material after pellet injection, taking into account both the Alfvén and ohmic currents which shortcircuit the charge separation creating the drift. The model is suitable for rapid calculation of the radial drift displacement, and can be useful for e.g. modelling of disruption mitigation via pellet injection.

1. Introduction

Pellet injection is an effective tool for modifying the density profile in fusion devices, and can be used for both fuelling and plasma control (Pégourié 2007). It has also been employed successfully to mitigate transient events in tokamaks, e.g. edge localised modes (Lang et al. 2015) and disruptions (Reux et al. 2021). The use of pellets to control such events is also planned for ITER (Baylor et al. 2009; Hollmann et al. 2015; Lehnen et al. 2018).

In order to assess the performance of pellet injection schemes for future tokamaks, such as ITER, it is important that accurate estimates of the modified density profile created by the pellets are included in the modelling tools used to simulate such events. This can only be achieved through an understanding of the underlying physics of the mass deposition after pellet injection.

When a pellet is injected into a hot, magnetically confined plasma, it travels through the plasma in solid form while the outer layers are continuously ablated by the energy flux from the hot background plasma, resulting in material being deposited along the pellet trajectory. The cloud of ablated material initially has a cold dense structure – a plasmoid – which drifts towards the low-field side of the torus. This is caused by the charge separation that takes place due to electron and ion drifts in the inhomogeneous magnetic field, leading to the build-up of a vertical electric field, and the resulting $\boldsymbol {E}\times \boldsymbol {B}$-drift moves the ablated material across magnetic field lines, where $\boldsymbol {E}$ and $\boldsymbol {B}$ denote the electric and magnetic field, respectively (Parks, Sessions & Baylor 2000; Rozhansky et al. 2004; Pégourié et al. 2006).

The strength of the electric field, and hence the drift velocity, is determined by the mechanisms which can shortcircuit the charge separation inside the plasmoid. The dominant ones are the emission of Alfvén waves from the two ends of the plasmoid (Parks et al. 2000) and the flow of ohmic current parallel to the field lines (Pégourié et al. 2006). Mathematically, the evolution of the pellet cloud is governed by a vorticity equation similar to that used to describe so-called blob transport in the plasma scrape-off layer (Krasheninnikov, D'Ippolito & Myra 2008).

There is a wealth of experimental evidence for radial cross-field drift following pellet injection in current tokamaks and stellarators, e.g. in DIII-D (Baylor et al. 2007), ASDEX Upgrade (Lang et al. 1997; Müller et al. 1999), FTU (Terranova et al. 2007), MAST (Garzotti et al. 2010) and W7-X (Baldzuhn et al. 2019). However, these studies consider small fuelling pellets, and there is much less experimental data on radial drifts from strongly perturbing pellets used for disruption mitigation, although there are recent indications that radial drifts may be important also in such cases at DIII-D and JET (Kong et al. 2022; Lvovskiy et al. 2022). Due to the complexity of the theoretical models, computer codes developed to simulate the cross-field drift are computationally expensive (Strauss & Park 1998, 2000; Aiba et al. 2004; Ishizaki & Nakajima 2011; Samulyak et al. 2021). Therefore, simplified scaling laws, based on current experimental observations, are often used (Baylor et al. 2007; Koechl et al. 2018). Such expressions are of limited use for modelling ITER plasmas, which will have much higher temperatures and magnetic fields. In many cases, e.g. in the currently used disruption mitigation models, the radial drift of the pellet cloud is neglected altogether, for simplicity (Vallhagen et al. 2022). This is particularly problematic in the case of pure hydrogen pellets (Matsuyama 2022), as their clouds can reach significant over-pressure due to negligible radiative energy losses, thus their drifts can be large and therefore affect the pellet penetration and material deposition substantially.

The purpose of this paper is to develop a semi-analytical model for the cross-field drift motion of the ionised plasmoid, taking into account both the Alfvén and ohmic currents. Our aim is to extract the key physical mechanisms described by the codes mentioned above and condense the result into a computationally efficient model. We consider current conservation directly, rather than formulating a vorticity equation for the system, generalising the description of the parallel connection of the ohmic current, and clarifying elements present in the existing literature. Factors such as the assumed shape of the plasmoid and our neglect of its structure along the magnetic field will quantitatively affect the plasmoid dynamics, but will not affect the qualitative nature of the results presented here.

2. Physical model

The motion of the plasmoid arises because of an $\boldsymbol {E}\times \boldsymbol {B}$-drift in the direction of the major radius; the electric field builds up due to the current from the magnetic (curvature + $\boldsymbol {\nabla } B$) drift of the particles, while the time variation of this electric field gives rise to a partially cancelling polarisation drift current. The total radial shift is determined by the drift velocity reached and its duration, which is approximately the time it takes for the cloud to expand one connection length along the field lines ($t\sim {\rm \pi}R_{m} q/c_s$, where $R_{m}$ is the major radius, $c_s$ is the sound speed and $q$ is the safety factor or inverse of the rotational transform of the magnetic field). At this time, magnetic drift currents in the outboard and inboard portions of the cloud cancel out (analogously to a tokamak equilibrium).

In order to mathematically describe the pellet dynamics, we formulate the current-conservation equation for the system, describing the balance between the divergent parts of the currents necessary to maintain quasineutrality. Working within a single-fluid formalism, we introduce the mass density $\rho$, the mass flow velocity $\boldsymbol {v}$, which appears in the total time derivative $d_t=\partial _t+\boldsymbol {v}\boldsymbol {\cdot }\boldsymbol {\nabla }$, the total pressure including the electron and ion pressure components $p=p_e+p_i$, as well as the current density and the magnetic field vectors, $\boldsymbol {j}$ and $\boldsymbol {B}$. In addition, $\boldsymbol {B}=\boldsymbol {b}B$ with the unit vector $\boldsymbol {b}$, the curvature vector of the field lines is $\boldsymbol {\kappa }=\boldsymbol {b}\boldsymbol {\cdot } \boldsymbol {\nabla } \boldsymbol {b}$, and $\mu _0$ denotes the vacuum permeability.

The pellet cloud has higher pressure than the surrounding plasma since it is continuously heated by hot electrons from the latter (Parks & Turnbull 1978). A current perpendicular to the magnetic field lines arises in response to this excess pressure, but we note that the dynamics involved in the drift of the plasmoid is slower than the timescale of compressional Alfvén waves, so that the largest terms in the magnetohydrodynamic (MHD) force balance equation

(2.1)\begin{equation} \rho \frac{{\rm d} {\boldsymbol{v}}}{{\rm d}t} =\boldsymbol{j}\times \boldsymbol{B}-\boldsymbol{\nabla} p, \end{equation}

describe an approximately static force balance between the plasma pressure and the magnetic field. The total current takes the form

(2.2)\begin{equation} \boldsymbol{j} = j_\|\boldsymbol{b} + \frac{\boldsymbol{B}\times \boldsymbol{\nabla} p}{B^2} + \frac{\rho}{B}\boldsymbol{b}\times\frac{{\rm d} {\boldsymbol{v}}}{{\rm d}t}, \end{equation}

and the divergence of the diamagnetic current, the second term on the right, describes the charge accumulation due to the magnetic drifts (driven by field line curvature and field strength inhomogeneity). This is approximately given by

(2.3)\begin{equation} \boldsymbol{\nabla} \boldsymbol{\cdot} \left(\frac{\boldsymbol{B}\times \boldsymbol{\nabla} p}{B^2}\right) = \boldsymbol{\nabla} \boldsymbol{\cdot} \left[ p \boldsymbol{\nabla} \times \left(\frac{\boldsymbol{B}}{B^2}\right)\right] \approx \boldsymbol{\nabla} \boldsymbol{\cdot} \left( 2p\frac{\boldsymbol{b} \times \boldsymbol{\nabla} B}{B^2}\right) \equiv \boldsymbol{\nabla} \boldsymbol{\cdot} \boldsymbol{j}_{\boldsymbol{\nabla} B}. \end{equation}

Using $(\boldsymbol {b} \times \boldsymbol {\nabla } B)/B = \boldsymbol {b} \times \boldsymbol {\kappa }$ we can write the expression for current conservation in the form

(2.4)\begin{equation} 0 = \boldsymbol{\nabla} \boldsymbol{\cdot} \boldsymbol{j} \approx \boldsymbol{\nabla} \boldsymbol{\cdot} \left[j_\|\boldsymbol{b} + \frac{\boldsymbol{b}}{B} \times\left( 2p \boldsymbol{\kappa} +\rho \frac{{\rm d} {\boldsymbol{v}}}{{\rm d}t} \right)\right]. \end{equation}

The time-dependent term in (2.4) is the current due to the polarisation drift

(2.5)\begin{equation} \boldsymbol{j}_{\dot{\boldsymbol{E} }}=\rho \frac{\boldsymbol{b}}{B}\times \frac{{\rm d} \boldsymbol{v}}{{\rm d}t}. \end{equation}

The resistive-MHD Ohm's law $\boldsymbol {E}+(\boldsymbol {v}\times \boldsymbol {B})=\eta \boldsymbol {j}$ implies that, in the limit of modest resistivity $\eta$, the perpendicular mass flow $\boldsymbol {v}_\perp$ is dominated by $\boldsymbol {E}\times \boldsymbol {B}$ motion. For the low-frequency process of interest inside the pellet cloud, the electric field is electrostatic $\boldsymbol {E}=-\boldsymbol {\nabla } \phi$, with the electrostatic potential $\phi$, and thus we write the cross-field velocity as $\boldsymbol {v}_\perp \approx (\boldsymbol {b}\times \boldsymbol {\nabla } \phi )/B$.

The parallel current ($j_\|$) must adjust to make the total current divergence free; that is

(2.6)\begin{equation} \boldsymbol{\nabla}\boldsymbol{\cdot} \boldsymbol{j}=\boldsymbol{\nabla} \boldsymbol{\cdot}\left[j_\| \boldsymbol{b}+\boldsymbol{j}_{\boldsymbol{\nabla} B}+\boldsymbol{j}_{\dot{\boldsymbol{E} }}\right]=0. \end{equation}

Eliminating the contribution which describes the balance of the diamagnetic and parallel currents in the background equilibrium plasma, we are left with the perturbation of the current continuity equation driven by the excess pressure of the plasmoid.

In the very early phase of plasmoid acceleration, $\boldsymbol {j}_{\boldsymbol {\nabla } B}$ is approximately balanced by $\boldsymbol {j}_{\dot {\boldsymbol {E} }}$. At such short times, the length of the pellet cloud is much shorter than the distance around the torus, $t\ll R_m q/c_s$, and the plasmoid is thus poloidally and toroidally localised. If the aspect ratio of the torus is large, the curvature vector of the magnetic field is approximately $\boldsymbol {\kappa }=-\hat {R}/R_m$, where $\hat {R}$ is the unit vector in the direction of increasing major radius. For convenience we introduce the unit vector $\hat {Y}=\boldsymbol {b}\times \hat {R}$, so the direction of $\boldsymbol {j}_{\boldsymbol {\nabla } B}$ is $-\hat {Y}$, which is nearly vertical. As the electric field rises in this early stage, $\boldsymbol {j}_{\dot {\boldsymbol {E} }}$ evolves to point in the $\hat {Y}$ direction everywhere in the cloud. Later the $j_\|$ term starts to dominate over $\boldsymbol {j}_{\dot {\boldsymbol {E} }}$ in balancing $\boldsymbol {j}_{\boldsymbol {\nabla } B}$, setting the quasi-steady speed of the plasmoid.

We may integrate (2.4) over some convenient volume $V$ with boundary $\partial V$, and apply the divergence theorem to obtain

(2.7)\begin{equation} 0 = \int _{\partial V} \left[ \left(\frac{\rho}{B} \frac{{\rm d}\boldsymbol{v}}{{\rm d}t}-\frac{2p}{BR_m}\hat{Y}\right)+ j_{\|}\boldsymbol{b}\right]\boldsymbol{\cdot} \hat{n}\,{\rm d}S, \end{equation}

where $\hat {n}$ is a unit vector pointing outwards from $V$. We align the integration volume $V$ with the cloud by choosing it to be a magnetic flux tube extending along the length of the cloud. Since the magnetic field lines are curved, the end faces of the flux tube, which we denote by $\delta S$, are not quite parallel. We choose the flux tube to have rectangular cross-section with the lower boundary running through the middle of the cloud, separating the upper, blue and lower, red, parts of the cloud shown in figure 1, where the integration volume $V$ is sketched. The length of the cloud along the field line is $L_{\rm cld}$, and the upper boundary of the domain is located just above the cloud.

Figure 1. Schematic views of the ablation cloud and the field lines connecting the various parts of it from different perspectives; the green lines indicate the boundaries of the integration volume $V$: (a) parallel currents and magnetic drift currents indicated in the $y$–$z$ plane, (b) from the side looking in the toroidal direction, (c) from the top and (d) with unwrapped field lines (black dashed), connecting different parts of the cloud after a distance $L$. The cloud expands at the speed of sound $c_s$ in both directions, so that $L_{\mathrm {cld}}=2 c_s t$. We assume the pellet ablation cloud to be symmetric in $z$ (and $y$) with respect to the $y$ ($z$) axis in (a). The pellet is indicated in (b,c) by the black dot, from which the cloud diverges.

For simplicity, we assume that the pellet is injected in the horizontal midplane and therefore (by symmetry) is always located in the middle of the cloud in the direction along the magnetic field and in the vertical direction. The surface normal $\hat {y}$ of the lower surface of $V$ coincides with $\hat {Y}$ in the poloidal plane that contains the pellet, and rotates in the poloidal plane as one follows the field line along the flux tube $V$. The relation between $\hat {y}$, $\hat {Y}$ and $\hat {R}$ is

(2.8)\begin{equation} \hat{y} = \cos{\theta} \,\hat{Y} + \sin{\theta}\, \hat{R}, \end{equation}

where $\theta \approx \varphi /q \approx z/qR_{m}$ is the poloidal angle, $\varphi$ is the toroidal angle and $z$ is the coordinate along the magnetic field lines; we take $z=0$ in the poloidal plane of the pellet. The dimensions of the integration volume in the $\hat {R}$ and $\hat {y}$ directions are $\Delta R$ and $\Delta y$, respectively.

The contribution from the first term, $\boldsymbol {j}_{\dot {\boldsymbol {E} }}$, to (2.7) thus becomes

(2.9)\begin{equation} I_{\dot{\boldsymbol{E} }} = \int _{{-}L_{\mathrm{cld}}/2}^{L_{\mathrm{cld}}/2} \int _0^{\Delta R} \frac{\rho}{B^2} \frac{{\rm d}E_y}{{\rm d}t}\hat{y}\boldsymbol{\cdot} \hat{y}\,{\rm d}R\,{\rm d}z = \frac{\bar{n}\langle m_i \rangle\Delta R}{(1+\langle Z \rangle)B^2} \frac{{\rm d}E_y}{{\rm d}t}, \end{equation}

where we have noted that the field-line-integrated mass density is $\bar {n}\langle m_i\rangle /(1+\langle Z \rangle )$ (neglecting the mass of the electrons), $\bar {n} = \sum _{i} \bar {n}_i + \bar {n}_e = \sum _i \bar {n}_i(1+\langle Z \rangle )$ is the field-line- integrated total density of all species (including electrons) inside the cloud (with $\bar {n}_i$ and $n_e$ denoting the field-line-integrated density of ion species $i$ and electrons, respectively), $\langle m_i \rangle$ is the average ion mass inside the cloud and $\langle Z \rangle$ is the average ion charge inside the cloud.

Considering the second term, $\boldsymbol {j}_{\boldsymbol {\nabla } B}$, we assume that the pressure is constant along the field lines inside the cloud, with equal electron and ion temperatures, denoted by $T$. The contribution from the second term of (2.7) then becomes

(2.10)\begin{align} I_{\boldsymbol{\nabla} B} & = \int _{{-}L_{\mathrm{cld}}/2}^{L_{\mathrm{cld}}/2} \int _0^{\Delta R} -\frac{2(p-p_{\mathrm{bg}})}{BR_{m}} \hat{Y}\boldsymbol{\cdot}\hat{y}\,{\rm d}R\,{\rm d}z = \int _{{-}L_{\mathrm{cld}}/2}^{L_{\mathrm{cld}}/2} -\frac{2(p-p_{\mathrm{bg}})\Delta R}{BR_{m}}\cos{\left(\frac{z}{qR_{m}}\right)}\,{\rm d}z \nonumber\\ & ={-}\frac{4(p-p_{\mathrm{bg}})\Delta Rq}{B}\sin{\left(\frac{L_{\mathrm{cld}}}{2qR_{m}}\right)} ={-}\frac{4(\bar{n}T-L_{\mathrm{cld}}n_{\mathrm{bg}}T_{\mathrm{bg}})\Delta Rq}{BL_{\mathrm{cld}}}\sin{\left(\frac{L_{\mathrm{cld}}}{2qR_{m}}\right)}. \end{align}

The background pressure $p_{\mathrm {bg}}$ enters via the contribution from the upper surface of the integration volume. We see, as noted in the introduction, that the assumptions simplifying the parallel structure of the cloud will quantitatively affect the final results, but accounting for parallel structure will not affect the essential qualitative description of the plasmoid motion.

The key to calculating how the parallel current contributes to the drift motion is to find the relation between the parallel current $j_{\|}$, which flows through the background plasma (beyond the ends of the cloud), and the electric field responsible for $\boldsymbol {E}\times \boldsymbol {B}$ motion, which are related via the electrostatic potential $\phi$ along the plasmoid length. As the pellet flies through the plasma, it undergoes continuous ablation and thus generates a sequence of ablation clouds residing on different field lines. Each of these clouds expands along the magnetic field whilst drifting across it. It is important to note that the cloud drift velocity exceeds the speed of the pellet. We can thus regard the pellet as stationary, which simplifies our discussion.

With these facts in mind, we now study the evolution of the electrostatic potential along each field line. We fix our attention on one particular field line and denote by $\tau$ the time that has elapsed since pellet material first arrived there. This time is in general different from the time $t$ that has passed since this material was originally ablated from the pellet. (Alternatively, in the limit of very high electrical conductivity, it is possible to regard the field lines as ‘frozen into’ the pellet cloud, in which case it is better to consider a field line moving with the pellet cloud. In this case $t=\tau$.)

It is convenient to introduce $L$, the distance along a field line, outside the cloud, which connects the two ends of the cloud; note that $L$ depends on the coordinates identifying a field line and may be different for different field lines in our integration volume $V$. In our large-aspect-ratio approximation, the value of $L = 2 {\rm \pi}R_{m} N$ is equal to the circumference of the torus, $2{\rm \pi} R_{m}$, times the number of turns, $N$, after which the field line connects the two end caps of $V$. This number will in general vary over the cross-section of the flux tube.

The evolution of the electrostatic potential along a field line connecting the oppositely charged parts of the cloud after a length $L$ is illustrated in figure 2. The physical picture of the evolution of this potential is the following: the interface between the end of the plasmoid and the background plasma represents an evolving perturbation, expanding along the field lines at the local sound speed $c_s$ of the pellet material inside the plasmoid. The potential difference between the cloud and the background plasma, along with the plasmoid drift, excite shear Alfvén waves, which are emitted from these interfaces and propagate away from the plasmoid, along field lines through the background plasma, at the local Alfvén speed, $C_A$. For $\tau \ll L/(2 C_A)$, the potential perturbations associated with the Alfvén waves will not have reached each other yet. Thus, the current carried away from the ends of the cloud is determined by the polarisation current resulting from the time-varying potential at the wave fronts, giving rise to the Alfvén current (Scholer 1970). When $\tau =L/(2 C_A)$, the waves emerging from the opposite sides of the cloud meet and interfere with each other. Eventually, a steady state, without propagating waves, is reached when $\tau \gg L/(2 C_A)$. At this stage, the parallel current is instead determined by Ohm's law.

Figure 2. Sketch of the electrostatic potential $\phi (z)$ along a field line connecting the two ends of the cloud, at different values of $y$, characterised by potentials $\phi _A$ and $\phi _B$. We show three representative times: at $\tau _1< L/(2C_A)$ potential perturbations propagating out from the ends of the cloud at the Alfvén speed have not yet met along the field line (solid black line). The perturbations meet at $\tau _2=L/(2C_A)$ (dashed blue). After a long time (compared with Alfvén time scales), $\tau _3\gg L/(2C_A)$, the potential has reached a quasi-steady state where an ohmic current flows between the connected ends of the cloud (dash-dotted green). Note that the cloud length $L_{{\rm cld}}$ is exaggerated in the figure; in reality it is much shorter than the distance along the field line between the connected ends of the cloud.

Thus, the dominant contribution to the $j_{\|}$ current, in the initial phase, is associated with the Alfvén wave propagating from the ends of the drifting cloud (Parks et al. 2000). It is proportional to the electric field inside the cloud, as outlined below, and can be described by the so-called Alfvén conductivity, $\varSigma _A = 1/ R_A = 1/(\mu _0 C_A)$. In the later stages, the ohmic current along the field lines connecting the oppositely charged parts of the cloud (Pégourié et al. 2006) becomes dominant. There is also a contribution to the current caused by the drift resulting from the cloud viscosity, which has been shown by Rozhansky et al. (2004) to be less significant and will be neglected here.

2.1. Parallel current

When calculating the contribution of $j_{\|}$ to the integral (2.7), only the end caps (area $\delta S$) of this flux tube will contribute, as otherwise $\boldsymbol {b}\boldsymbol {\cdot } \hat {n} = 0$. Consider first the contribution from a smaller flux tube, whose end caps have area $\partial s_<$, that only contains field lines for which $C_A t \ll L/2$, that is, for which Alfvén waves propagating from the ends of the cloud have not had time to meet. As the parallel electric field $E_{\|}$ is small in the established hot background plasma outside the cloud, except at the wave front, we can express $E_{\|}$ in Fourier space as

(2.11)\begin{equation} E_{\|} ={-}{\rm i}k_{\|} \phi + {\rm i}\omega A_{\|} \approx 0, \end{equation}

and relate the electrostatic and vector potential via the Alfvén speed

(2.12)\begin{equation} A_{\|} = \phi/(\omega/k_{\|}) = \phi/C_A. \end{equation}

Using Ampère's law we can relate $j_{\|}$ to $A_{\|}$ and thence to $\phi$ as

(2.13)\begin{equation} j_{\|} ={-}\frac{\boldsymbol{\nabla} _\perp^2 A_{\|}}{\mu _0} ={-}\frac{\boldsymbol{\nabla} _\perp^2 \phi}{\mu _0C_A}. \end{equation}

Assuming that the whole cloud moves at the same radial velocity, the electric field $E_y = -\partial \phi /\partial y$ must be constant inside the cloud, i.e.

(2.14)\begin{equation} \boldsymbol{\nabla} _\perp^2 \phi ={-}E_y \left[ \delta(y-\Delta y) - \delta(y+\Delta y) \right], \end{equation}

where $\delta$ denotes the Dirac delta function. If we set $\partial s_<$ to the part of $\partial S$ for which $C_A t < L/2$, the contribution from the Alfvén part of the parallel current becomes

(2.15)\begin{align} I_{\|,A} & = 2\int _{\partial s_<} \frac{E_y\delta(y-\Delta y)}{\mu _0C_A}\,{\rm d}y\,{\rm d}R \nonumber\\ & = 2\int _0^{\Delta R}\int_0^{\Delta y} \varTheta(\partial s_<;y,R)\frac{E_y\delta(y-\Delta y)}{\mu _0C_A}\,{\rm d}y\,{\rm d}R \nonumber\\ & = 2 P_A \Delta R \frac{E_y}{\mu _0C_A} = 2 P_A \Delta R \frac{E_y}{R_A} , \end{align}

where the function $\varTheta (\partial s;y,R)$ is $1$ for the $y$ and $R$ values corresponding to field lines crossing the surface $\partial s$ where $C_A t < L/2$ is satisfied, and zero otherwise; and $P_A$ is the fraction $\partial s_< /\partial S$.

Now consider the field lines crossing the area $\partial s_>$, i.e. the field lines for which $C_A t\gg L/2$. On these field lines, the Alfvén waves emanating from either side of the cloud have already met and decayed, and there is no longer any polarisation current. Only the ohmic current $j_{\|}$ remains and, being divergence free, it must be constant along the field in the large-aspect-ratio limit. This current is related to the parallel electric field by Ohm's law

(2.16)\begin{equation} j_{\|} = \sigma_{\|} E_{\|} ={-}\sigma_{\|}\boldsymbol{\nabla} _{\|}\phi. \end{equation}

As $j_{\|}$ is constant along the field lines, so is $\boldsymbol {\nabla } \phi$, which means that

(2.17)\begin{equation} j_{\|} = \sigma_{\|} E_{\|} ={-}\sigma_{\|} \frac{\phi _B -\phi}{L} ={-}\sigma_{\|} \frac{E_y (y_B-y)}{L}, \end{equation}

where $y$ denotes the vertical coordinate at which the field line emanates from one end of the cloud and $y_B$ that where it hits the other end. The electric field $E_y$ has been assumed to be constant along the field line.

Let us now denote by $\partial s_i$ the subset of $\delta s_>$ containing only field lines connecting to the opposite side of the cloud after a distance $L = 2{\rm \pi} R_{m} i$, i.e. connecting after exactly $i$ toroidal turns. If $i \gg 1$, the connection is essentially random, so that the values of $y$ and $y_B$ are uncorrelated and $\int y_B\,{\rm d}y = 0$. The total ohmic current flowing along field lines in $\delta s_i$ thus becomes

(2.18)\begin{align} I_{\|,\mathrm{ohm}}^{(i)}(\tau\gg L/(2C_A))& ={-}2\int _{\partial s_i} \sigma _{\|} \frac{E_y(y_B - y)}{L}\,{\rm d}y\,{\rm d}R\nonumber\\ & ={-}2P_i \int _0^{\Delta y} \int _0^{\Delta R} \varTheta(\partial s_>;y,R) \sigma _{\|} \frac{E_y(y_B - y)}{2{\rm \pi} R_{m} i}\,{\rm d}y\,{\rm d}R \nonumber\\ & = P_i\sigma _{\|} \frac{E_y\Delta y^2\Delta R}{2{\rm \pi} R_{m} i}, \end{align}

where $P_i = \delta s_i / \partial S$ is the fraction of the cloud connecting to the opposite side after $i$ toroidal turns. This result is similar to the corresponding expression, equation (2), in Commaux et al. (2010), up to an order-unity factor accounting for the finite electron collision time. The total ohmic current is obtained by summing over all values of $i$.

For $\tau \gtrsim L/(2C_A)$, the current will make a transition from 0 to $I_{\|,\mathrm {ohm}}^{(i)}(\tau \gg L/(2C_A))$¹ over a time scale similar to $L/C_A$, so that we can write

(2.19)\begin{equation} I_{\|, \mathrm{ohm}} = \sum _{i=1}^\infty f\left(\frac{\tau}{L/(2C_A)}\right)I_{\|,\mathrm{ohm}}^{(i)}(\tau\gg L/(2C_A)), \end{equation}

where $f(0)=0$ and $f\rightarrow 1$ for large arguments. The detailed form of $f$ is determined by the interaction of the Alfvén waves propagating from opposite sides of the cloud, which is outside the scope of the present work. Here, we instead make the approximation that $f = \theta ({\tau }/({L/(2C_A)})-1)$, where $\theta$ is the Heaviside step function. This is also used in Pégourié et al. (2006). Note that this assumption on $f$ underestimates the time until the onset of the ohmic current, thus overestimating the importance of the ohmic current contribution. With this assumption for $f$, we can write the total ohmic current as

(2.20)\begin{equation} I_{\|, \mathrm{ohm}} = \sum _{i=1}^N I_{\|,\mathrm{ohm}}^{(i)}(\tau> L/(2C_A)), \end{equation}

where $N = \lfloor 2C_A \tau /(2{\rm \pi} R_{m}) \rfloor = \lfloor \tau /t_0 \rfloor$ is the maximum number of toroidal turns the Alfvén wave front has had time to make, with $t_0$ the time for the Alfvén wave to propagate one turn around the torus, accounting for the fact that emission is from both ends of the cloud. The notation $\lfloor x\rfloor$ gives the greatest integer less than or equal to $x$.

2.1.1. Fraction of the cloudlet cross-section connected to the opposite side

We will now calculate the fraction $P_i$ of the cloudlet cross-section that connects to the opposite side during the $i$th turn, assuming an irrational safety factor on the flux surface. We note though that rational flux surfaces have been shown to affect the ablation process – owing to the smaller reservoir of hot electrons that can enter the pellet cloud – and can also produce a drift braking effect, caused by the effective shortcircuiting of the potential variation along field lines (Commaux et al. 2010; Sakamoto et al. 2013). Our model does not account for this, and as such it may provide a conservative upper estimate of the drift distance. As we shall see, most of the contribution to the current comes from terms with $i \gg 1$, i.e. from field lines that encircle the torus many times before connecting the two ends of the cloud. According to Weyl's lemma (Helander 2014), whether a given field line starting from one side of the cloud connects to the other side in a large number of turns is essentially random. We can thus speak of the probability of such a connection, and this probability depends on the fraction of the poloidal cross-section that the cloudlet covers, which is $\Delta y/({\rm \pi} r)$, where $r$ is the characteristic minor radius at the cloudlet position. Therefore, the total connected fraction $P_{\mathrm {con}}^{\mathrm {tot}}$ increases between turn $N$ and $N+1$ in the following way:

(2.21)\begin{equation} P_{\mathrm{con}}^{\mathrm{tot}}(N+1)-P_{\mathrm{con}}^{\mathrm{tot}}(N) = \frac{\Delta y}{{\rm \pi} r}\left(1-P_{\mathrm{con}}^{\mathrm{tot}}(N)\right). \end{equation}

The solution of this difference equation is

(2.22)\begin{equation} P_{\mathrm{con}}^{\mathrm{tot}}(N) = 1 - \left(1-\frac{\Delta y}{{\rm \pi} r}\right)^N, \end{equation}

and we can now express

(2.23)\begin{equation} P_i = P_{\mathrm{con}}^{\mathrm{tot}}(i+1)-P_{\mathrm{con}}^{\mathrm{tot}}(i) = \frac{\Delta y}{{\rm \pi} r}\left(1-\frac{\Delta y}{{\rm \pi} r}\right)^i. \end{equation}

This estimate is consistent with figure 3 in Pégourié et al. (2006). We can also express the fraction $P_A$ (determining the size of the Alfvén current) as $P_A = 1 - P_{\mathrm {con}}^{\mathrm {tot}}$.

Combining (2.18), (2.20) and (2.23), the ohmic current contribution can now be expressed as

(2.24)\begin{equation} I_{\|, \mathrm{ohm}} = \sum _{i=1}^N P_i\sigma _{\|} \frac{E_y\Delta y^2\Delta R}{2{\rm \pi} R_{m} i} = \frac{E_y\Delta R}{R_{\mathrm{eff}}}, \end{equation}

with the inverse effective resistivity $1/R_{\mathrm {eff}}$ given by

(2.25)\begin{equation} \frac{1}{R_{\mathrm{eff}}} = \sum _{i=1}^N P_i\sigma _{\|} \frac{\Delta y^2}{2{\rm \pi} R_{m} i} = \sigma _{\|}\frac{\Delta y^3}{2{\rm \pi} R_{m}{\rm \pi} r }\sum _{i=1}^N \frac{1}{i}\left(1-\frac{\Delta y}{{\rm \pi} r}\right)^i. \end{equation}

For $N\rightarrow \infty$ we may use $\sum _{i=1}^\infty (1-x)^i/i=-\ln {x}$, giving

(2.26)\begin{equation} \frac{1}{R_{\mathrm{eff}}}=\sigma _{\|}\frac{\Delta y^3}{2{\rm \pi}^2 R_{m} r} \ln{\frac{{\rm \pi} r}{\Delta y}}. \end{equation}

Concerning when the $N\rightarrow \infty$ limit is meaningful to take, we must appreciate that, depending on the resistivity of the cloud, the cloud may or may not be frozen into the magnetic field, which determines whether field lines are dragged along with the cloud, or the field lines slip with respect to the cloud.² It is worth re-iterating that the generation of Alfvén waves by the propagating potential perturbation – and so the existence of Alfvén resistivity – does not require the field lines to be frozen in on the drift time scale.

As we will see later, the number of connected turns $N$ becomes large during the drift motion, so that taking $N\rightarrow \infty$ is a valid approximation for the majority of the drift motion, as long as the magnetic field diffusion is slow enough (i.e. the cloud temperature is high enough) that the cloud does not become disconnected from the field lines where the electrostatic potential has been set up.

The picture becomes more complicated if the magnetic field diffusion time scale is fast compared with the drift motion. This is typically the case for low cloud temperatures (e.g. pellets doped with highly radiating impurities), where the conductivity in the cloud is low and the resistive diffusion coefficient is large. In this case, the potential along a given field line will not only be determined by the local cloud properties, but will be affected by all material which has drifted past the field line under consideration. When pellet material first arrives, the Alfvén current dominates. On the other hand, long after the ablation flow started to cross a given field line, the potential along this field line will reach a quasi-stationary profile similar to the case when the field line remains frozen into the cloud for a long time, and thus $N\rightarrow \infty$ also in this case. As the pellet motion is typically slow compared with the other processes of interest, the latter limit should dominate for the majority of the ablated material in most cases even for low cloud temperatures.

The fraction of connected field lines, $P_A$, converges somewhat slower than the effective resistivity $R_{\mathrm {eff}}$. We therefore keep $N$ finite in the expression for $P_A$ for hot clouds. For cold clouds, for the first material drifting past a new part of the background plasma $N$ remains equal to zero, and $P_A=1$. However, as the potential reaches its quasi-stationary value, $P_A \rightarrow 0$ for the whole drift motion (i.e. the parallel current will be dominated by the ohmic component).

2.2. Current balance

We are now finally ready to sum up the various contributions to the current balance and obtain an equation for $E_y$ in terms of the parameters characterising the pellet cloud and the background plasma. From (2.7) we have

(2.27)\begin{align} 0 & = \frac{I_{\boldsymbol{\nabla} B} + I_{\dot{\boldsymbol{E} }} + I_{\|, A} + I_{\|, \mathrm{ohm}}}{\Delta R} \nonumber\\ & ={-}\frac{4(\bar{n}T-L_{\mathrm{cld}}n_{\mathrm{bg}}T_{\mathrm{bg}})q}{BL_{\mathrm{cld}}}\sin{\left(\frac{L_{\mathrm{cld}}}{2qR_{m}}\right)} + \frac{\bar{n}\langle m_i \rangle}{(1+\langle Z \rangle)B^2} \frac{{\rm d}E_y}{{\rm d}t} + 2P_A \frac{E_y}{R_A} + \frac{E_y}{R_{\mathrm{eff}}}. \end{align}

Note that the factor $\sin {({L_{\mathrm {cld}}}/{qR_{m}})}$ will start to oscillate when $t\sim q R_{m}/c_s$, as $L_{\mathrm {cld}} \sim c_s t$, and the amplitude of the term in which this appears in (2.27) decreases as $1/L_{\mathrm {cld}}\propto 1/t$; this oscillation, together with the pressure equilibration (which occurs when $\bar {n}T=L_{\mathrm {cld}}n_{\mathrm {bg}}T_{\mathrm {bg}}$), effectively sets the time scale of the drift duration and eventually leads to a finite displacement for the drift. Also note that $c_{s}t_0/qR_{m}\sim c_{s}/C_A$ is small for typical fusion plasma parameters, meaning that $N$ becomes large during the drift duration, motivating us to take the upper limit of the sum in (2.25) to be infinite, when calculating $R_{\mathrm {eff}}$.

If the plasmoid and background-plasma properties do not depend on $E_y$, (2.27) becomes a linear first-order ordinary differential equation in $E_y$, which can be written in the form

(2.28)\begin{equation} \frac{{\rm d}E_y}{{\rm d}t} + g(t) E_y = f(t), \end{equation}

with

(2.29)\begin{equation} g(t) = \frac{(1+\langle Z \rangle)B^2}{\bar{n}\langle m_i \rangle}\left(2P_A \frac{1}{R_A} + \frac{1}{R_{\mathrm{eff}}}\right), \end{equation}

and

(2.30)\begin{equation} f(t) = \frac{4(1+\langle Z \rangle)B}{\langle m_i \rangle L_{\mathrm{cld}}}\left(T-\frac{L_{\mathrm{cld}}n_{\mathrm{bg}}}{\bar{n}}T_{\mathrm{bg}}\right)q\sin{\left(\frac{L_{\mathrm{cld}}}{2qR_{m}}\right)}. \end{equation}

This equation can be solved by using an integrating factor $e^{G(t)}$, so

(2.31)\begin{equation} E_y = e^{{-}G(t)}\left(E_{y0} + \int _0^t e^{G(t)}f(t)\,{\rm d}t\right), \end{equation}

where $E_{y0}=E_y(t = 0)$ and $G(t) = \int _0^t g(t)\,{\rm d}t$. For a hot cloud, we have

(2.32)\begin{align} G(t) & = \frac{(1+\langle Z \rangle)B^2}{\bar{n}\langle m_i \rangle}\left(2\left[1-\left( 1-\frac{\Delta y }{{\rm \pi} r}\right)^{N+1}\right]\frac{{\rm \pi} r}{\Delta y} \frac{t_0}{R_A} + \frac{t}{R_{\mathrm{eff}}}\right)\nonumber\\ & = 2\left[1-\left( 1-\frac{\Delta y }{{\rm \pi} r}\right)^{N+1}\right]\frac{{\rm \pi} r}{\Delta y} \frac{R_{\mathrm{eff}}}{R_A}\frac{t_0}{t_{\mathrm{acc}}} + \frac{t}{t_{\mathrm{acc}}}, \end{align}

where we have defined

(2.33)\begin{equation} t_{\mathrm{acc}} = \frac{\bar{n}\langle m_i \rangle R_{\mathrm{eff}}}{(1+\langle Z \rangle)B^2}. \end{equation}

This is the characteristic acceleration time scale if the ohmic current dominates over the Alfvén current; if $R_{\mathrm {eff}}/R_A$ is small (corresponding to a hot background plasma), or in the case of a cold cloud long after the ablation flow started to cross the local field line, the expression (2.32) reduces to

(2.34)\begin{equation} G(t) = \frac{t}{t_{\mathrm{acc}}}. \end{equation}

For a cold cloud shortly after the ablation flow started at the local field line, where $P_A = 1$, (2.32) reduces to the same expression but with $R_{\mathrm {eff}}$ replaced with $R_A$ in the expression for $t_{\mathrm {acc}}$.

Finally, as the radially outward drift velocity of the cloudlet is due to the $\boldsymbol {E}\times \boldsymbol {B}$ motion, it can be estimated as $E_y /B$. Time integration leads to an expression for the net radial displacement

(2.35)\begin{equation} \Delta r = \frac{1}{B}\int _{0}^{\infty} E_y \,\mathrm{d}t. \end{equation}

3. Parallel expansion and the final drift displacement

In this section, we complete the description of the pellet cloud by defining the density source resulting from pellet ablation. We then evaluate the drift of the pellet cloud, demonstrating its dependence on pellet composition and background-plasma temperature.

3.1. Model for the line-integrated density and cloud expansion

The line-integrated density can be determined based on an estimate of how many particles the cloud contains when it detaches from the pellet source. The latter can be obtained as the product of the ablation rate and the time during which the pellet source is ablating inside the cloud.

A widely used estimate for the mass ablation rate is given by

(3.1)\begin{equation} \mathcal{G}=\lambda(X)\left(\frac{T_{\mathrm{keV}}}{2}\right)^{5/3}\left(\frac{r_p}{r_{p0}}\right)^{4/3}n_{e20}^{1/3}, \end{equation}

where $\lambda (X)=[27.1+\tan {(1.48X)}]/1000 \,{\rm kg}\,{\rm s}^{-1}$, $X$ is the molecular fraction of deuterium in the pellet, $T_{\mathrm {keV}}$ is the background electron temperature in keV, $r_p$ is the pellet radius, $r_{p0}=2 \,{\rm mm}$ and $n_{e20}$ is the background electron density in units of $10^{20}\,{\rm m}^{-3}$ (Parks 2017). This expression is based on a version of the neutral gas shielding model (Parks & Turnbull 1978) that allows the pellet material to have both hydrogenic and noble gas components.

To determine the average detachment time (during which the pellet source contributes to the cloud), we estimate the initial acceleration $\dot {v}_0=\dot {v}(t=0)=\dot {E_y}(t=0)/B$ by balancing the first two terms in the current balance equation (2.27). The last two terms in (2.27) can be neglected, since in the initial phase, $E_y$ is small. The time derivative of the electric field then becomes

(3.2)\begin{equation} \frac{{\rm d}E_y}{{\rm d}t}=\frac{2B(1+\langle Z \rangle)}{\bar{n}\langle m_i \rangle R_{m}}\left(\bar{n} T-L_{\mathrm{cld}} n_{\mathrm{bg}}T_{\mathrm{bg}}\right), \end{equation}

so that the initial acceleration is

(3.3)\begin{equation} \dot{v}_0=\frac{1}{B}\frac{{\rm d}E_y}{{\rm d}t}=\frac{2(1+\langle Z \rangle)}{\langle m_i \rangle R_{m}}\left(T_0 - \frac{n_{\mathrm{bg}}}{n_0}T_{\mathrm{bg}}\right), \end{equation}

where $n_0=\bar {n}/L_{\mathrm {cld}}$ is the initial cloud density and $T_0$ is the initial temperature.

Initially, the pellet cloud is neutral, and it expands radially, but as soon as the particles are ionised, the expansion will continue along the magnetic field lines. The initial parallel expansion takes place at the speed of sound at a temperature of approximately $T_0$ (which is of the order of a few eV), and starts from a spherical cloud of cross-sectional area ${\rm \pi} \Delta y^2$. We can therefore estimate the density from mass conservation according to

(3.4)\begin{equation} \mathcal{G} = 2n_0\langle m_i \rangle c_s(T_0){\rm \pi} \Delta y^2 \Rightarrow n_0 = \frac{\mathcal{G}}{2\langle m_i \rangle c_s(T_0){\rm \pi} \Delta y^2}. \end{equation}

The average distance the ablated material must drift before it exits the initial expansion tube around the pellet is $\Delta y$. Assuming that the initial motion has a constant acceleration we find

(3.5)\begin{equation} \Delta y=v_0 t_{\mathrm{det}} + \dot{v}_0 t_{\mathrm{det}}^2/2, \end{equation}

and the average detachment time thus becomes

(3.6)\begin{equation} t_{\mathrm{det}} ={-}\frac{v_0}{\dot{v}_0}+\sqrt{\left(\frac{v_0}{\dot{v}_0}\right)^2+\frac{2\Delta y}{\dot{v}_0}}, \end{equation}

where $v_0=v_p$ is the initial cloud velocity relative to the pellet, which is equal in magnitude but opposite in sign to the pellet velocity (assuming the cloud would be frozen in to the field lines where it was ablated in the absence of the $\boldsymbol{E}\times \boldsymbol{B}$ acceleration). During the detachment time $t_{\mathrm {det}}$ the cloud expands to a length

(3.7)\begin{equation} L_c = 2c_s(T_0)t_{\mathrm{det}}, \end{equation}

which serves as an initial condition for the cloud length for the remainder of the drift motion. The line-integrated density is thus

(3.8)\begin{equation} \bar{n} = \frac{\mathcal{G}t_{\mathrm{det}}}{\langle m_i \rangle{\rm \pi} \Delta y^2}. \end{equation}

When the cloud detaches from the pellet, the temperature initially rises quickly due to the heating from hot electrons in the background plasma, but the details depend on the density and composition of the cloud. At low densities, the mean free path of the background-plasma electrons in the cloud is longer than the cloud itself. These electrons thus pass through the cloud and heat it relatively uniformly (Aleynikov et al. 2019; Arnold, Aleynikov & Helander 2021; Runov et al. 2021). Most of the literature, however, considers the opposite limit of a dense cloud, where the stopping power is so great that the hot electrons cannot easily pass through it. We only consider this case but note that it becomes inapplicable at low cloud densities and high background-plasma temperatures. The heating also depends on the pellet composition; if the pellet contains even a small amount of a high-Z radiative component, the radiation from the pellet cloud quickly reaches a balance with the heating from the background plasma, and therefore the temperature rises far more slowly (Matsuyama 2022). For pure hydrogen pellets, on the other hand, the radiation is too weak to have a major impact on the energy balance, and then the cloud temperature will relatively quickly increase to several tens of eV.

The dependence of the cloud temperature on the background-plasma temperature will be rather weak, as a higher background-plasma temperature means both an increased heating and an increased ablation rate, giving more particles to absorb and, in the case of a high-Z-doped pellet, radiate away the energy. The heat flux scales as $q_{\mathrm {bg}}\sim T_{\mathrm {bg}}^{3/2}$ (neglecting any scaling of the cloud cross-sectional area with the temperature), and the ablation rate scales as $\mathcal {G}\sim T_{\mathrm {bg}}^{5/3}$, so that the cloud temperature scales as $T\sim q_{\mathrm {bg}}/\mathcal {G}\sim T_{\mathrm {bg}}^{-1/6}$, i.e. a very weak scaling. Typical values for the cloud temperature, based on the results presented in Matsuyama (2022), are $T=5\,{\rm eV}$ for neon-doped pellets and $T=30\,{\rm eV}$ for pure hydrogen pellets. In the following we will assume the cloud temperature is constant during the drift motion and is independent of the background-plasma temperature. This approximation is, of course, quite crude but not more so than other simplifications we have employed.

Finally, as long as the cloud pressure is much higher than the background-plasma pressure, the cloud will expand by approximately the speed of sound inside the cloud, $c_s\approx \sqrt {(\gamma _e\langle z \rangle + \gamma _i)T/\langle m_i \rangle }$, with $\gamma _e = 1$ and $\gamma _i = 3$, and will slow down when the cloud pressure becomes comparable to the background-plasma pressure. Here, we assume that the expansion speed is equal to $c_s$ as long as the cloud pressure is higher than the background-plasma pressure, and then stops immediately when the cloud pressure becomes equal to the background-plasma pressure, i.e.

(3.9)\begin{equation} L_{\mathrm{cld}}\approx L_c + 2c_s\min(t,t_{\mathrm{pe}}), \end{equation}

where the pressure equilibration time is

(3.10)\begin{equation} t_{\mathrm{pe}} = \frac{T\bar{n}}{2c_s n_{\mathrm{bg}}T_{\mathrm{bg}}}. \end{equation}

With the parallel dynamics model presented here, we have all the details needed to evaluate the electric field inside the cloud, and the drift displacement can be calculated by evaluating the integral in (2.35). Analytical expressions for the drift displacement in various limits are given in Appendix A.

Before evaluating the above expressions for the drift displacement for an ITER-like scenario, as a validation exercise, we evaluate the drift distance for parameters representative of a JET shattered pellet injection (SPI) experiment studied by Kong et al. (2022) using the JOREK code; discharge #96874. Here, the drift displacement was accounted for by imposing a fixed prescribed shift between the ablating pellet shard and the location where the ablated material is eventually deposited. A shift of $\Delta r = 30\,{\rm cm}$ was found to yield the best agreement to the experimental density evolution. We choose this case due to the availability of an estimate of the total drift displacement under experimental conditions as close as possible to the ITER-like scenario studied below in § 3.2.

The injected pellet consisted of $1.6\times 10^{23}$ deuterium atoms and was shattered into $\sim$300 shards, which, using the Parks size distribution model (Parks 2016), gives a characteristic shard radius of $r_p=0.6\,{\rm mm}$. However, as shards of larger volume contribute more to the density build-up – that was matched to the experiment – we consider a representative pellet shard radius of $r_p = 1\,{\rm mm}$ in our estimate. Lacking values for $\Delta y$ and $T$ for the specific case, we set $\Delta y = 1.25$ cm and $T=30$ eV based on simulation results by Matsuyama (2022) of the same ITER-like scenario as the one studied below in § 3.2. Finally, using the representative geometrical and background plasma parameters $v_0 = 300\,{\rm m}\,{\rm s}^{-1}$, $q=1.5$, $B = 3.45\,{\rm T}$, $R_{m} = 3.5\,{\rm m}$, $r = 0.5\,{\rm m}$, $n_{\mathrm {bg}} = 8.5\times 10^{19}\,{\rm m}^{-3}$ and $T_{\mathrm {bg}}=7\,{\rm keV}$ (Kong et al. 2022), we arrive at an estimated drift displacement of $\Delta r = 28\,{\rm cm}$, in good agreement with the value found to match the experimental data. Although this estimate may be altered by a factor $\sim$2 within reasonable ranges of the relevant parameters, this result suggests that the present model is sufficiently accurate for order-of-magnitude estimates and qualitative studies, such as those performed for an ITER-like scenario in the next subsection.

3.2. Calculation of the drift distance in an ITER-like scenario

We now evaluate the above expressions for the drift displacement for parameters of interest in an ITER-like scenario, similar to that studied by Matsuyama (2022). In this scenario, the drifting pellet cloud is ablated from a pellet shard with radius $r_{p} = 2\,{\rm mm}$ located at major radius $R_{m} = 5\,{\rm m}$ and travelling with a speed of $v_0 = 500\,{\rm m}\,{\rm s}^{-1}$ towards the high-field side (i.e. the injection is from the low-field side). We also assume that the cloud is initially stationary in the laboratory frame, so that $E_{y0} = 0$. The background plasma has a free electron density of $n_{\mathrm {bg}} = 10^{20}\,{\rm m}^{-3}$ and the magnetic field strength is $B = 5\,{\rm T}$. Moreover, we set $q = 1$, $\Delta y = 1.25\,{\rm cm}$, (based on simulation results by Matsuyama 2022) and the average charge for the neon is approximately $\langle Z_{\mathrm {Ne}} \rangle \approx 2$ at 5 eV. The background-plasma temperature $T_{\mathrm {bg}}$ and the pellet composition will be varied.

Matsuyama (2022) uses a model similar to that used by Pégourié (2007), adapted to mixed neon–deuterium pellets, including a neutral gas and plasma shielding model for the pellet ablation and a volume-averaged single-cell Lagrangian model for the parallel expansion. However, Matsuyama (2022) only considers the early stages of the drift motion during the first $130\,\mathrm {\mu }{\rm s}$ after the cloud has detached from the pellet, for a single isolated cloud, and does therefore not include the effect of ohmic currents and rotational transform. Thus, the model by Matsuyama (2022) accounts for the same physical mechanisms concerning the drift motion as ours in the case of a cold cloud shortly after the ablation flow has started to cross the local field lines.³ He concluded that the drift displacement is likely to be substantial compared with the plasma minor radius for pure hydrogen pellets, but will be strongly reduced in the presence of even a small amount of neon. Here, we attempt to reproduce this result in the corresponding limit, and then extend it by calculating the drift displacement after a long time, including the effect of ohmic currents.

Figure 3 shows the drift displacement for cold clouds (30 eV for pure hydrogen, 5 eV otherwise). This is calculated by integrating (A2) (leading to (A4) if we integrate up to infinity), as a function of the background plasma temperature and pellet composition, with different integration times and assumptions regarding the ohmic currents. In panel (a) we consider the case when the ablation flow has just started to cross the local field lines, i.e. with the parallel current consisting only of the Alfvén current, and panel (b) shows the results for long after the ablation flow started to cross the local field lines, i.e. with the parallel current being purely ohmic.

Figure 3. Drift displacement as a function of background-plasma temperature and pellet composition for cold clouds (30 eV for pure hydrogen, 5 eV otherwise), with different integration times and assumptions for the parallel current. In panel (a) the parallel current is assumed to be purely Alfvénic (corresponding to when the ablation flow has just started to cross the local field lines), and in panel (b) the parallel current is assumed to be purely ohmic (corresponding to long after the ablation flow started to cross the local field lines). The solid lines correspond to performing the time integral of the drift velocity to $t=\infty$, as in (2.35), the dashed lines are obtained by integrating only to $130\,\mathrm {\mu }{\rm s}$.

The dashed lines in panel (a) are calculated with the assumption that the parallel current is purely Alfvénic, as was assumed by Matsuyama (2022), and the results are similar to those shown in figure 11 in Matsuyama (2022) within an order-unity factor, especially at high background-plasma temperatures. The variation with both the background temperature and pellet composition agrees reasonably well. We see, however, that when we extend the integration time to infinity (solid lines), the drift displacement increases significantly at high background plasma temperatures, so that even clouds with 100 % neon would drift several metres in the absence of ohmic currents, although the drift displacement is not strongly affected for temperatures ${\lesssim }1\,{\rm keV}$. This can be understood by considering that the pressure equilibration time becomes longer at high background-plasma temperatures (see (3.10)), so that the cloud can drift a significant distance after the first $130\,\mathrm {\mu }{\rm s}$. Moreover, in the absence of ohmic currents, the acceleration time scale is typically longer than $130\,\mathrm {\mu }{\rm s}$, so that the cloud continues to gain speed even after this time frame. For low background plasma temperatures, on the other hand, the pressure equilibration time becomes shorter than $130\,\mathrm {\mu }{\rm s}$ so the cloud does not drift significantly after this time.

In panel (b), where the parallel current is purely ohmic, we see that the drift displacement is reduced by approximately one order of magnitude when integrating up to $130\,\mathrm {\mu }{\rm s}$ (compare with panel a), and approximately two orders of magnitude when integrating to infinity. The scaling with the background-plasma temperature is also weaker, as anticipated above, because the resistivity determining the parallel current now scales with the background plasma temperature as $R_{\mathrm {eff}}\sim T_{\mathrm {bg}}^{-3/2}$, which mostly cancels the temperature scaling of the ablation rate $\mathcal {G}\sim T_{\mathrm {bg}}^{5/3}$ (there is some dependence on the background temperature left at lower background temperatures where the ratio of the cloud pressure and the background pressure is lower). Moreover, the effect of increasing the integration time beyond $130\,\mathrm {\mu }{\rm s}$ is now much smaller than in the absence of ohmic currents. This follows as the acceleration time scale $t_{\mathrm {acc}}$ is much shorter, so that the cloud decelerates rather than accelerates after the first $130\,\mathrm {\mu }{\rm s}$.

For neon-doped pellets, the drift displacement now ranges from a few cm up to ${\sim }20\,{\rm cm}$ at the highest relevant temperatures, which is small compared with both the plasma minor radius and the plume of shards in case of a SPI in an ITER-like scenario. The pure deuterium pellets, on the other hand, still have a drift displacement of tens of cm, which is a sizeable fraction of the plasma minor radius and comparable to the radial extent of the shard plume in case of an SPI. This result corroborates the conclusion made by Matsuyama (2022).

We now compare the results for the same plasma scenario as above using the expressions obtained with different limits and model assumptions. As we have seen in § 2.1.1, for hot clouds (e.g. pure deuterium pellets), the $N\rightarrow \infty$ limit of $R_{\mathrm {eff}}$ can be used while we keep $N$ finite in the expression of $P_A$. For cold clouds (e.g. neon-doped pellets), in the long-time limit (as the potential reaches its quasi-stationary value), the Alfvén part of the current can be neglected ($P_A=0$).

In figure 4, the full solution, which contains both the $I_{\|,A}$ and $I_{\|,\mathrm {ohm}}$ contributions obtained by numerically integrating (A1), is shown by a black curve for a pure deuterium pellet (panel a) and a $2\,\%$ neon-doped one (panel b). We also consider the cases representing the long and short-time limits, in terms of the time passed after the ablation flow first started to cross the local field lines. In the short-time limit (green long-dashed curve) $I_{\|,\mathrm {ohm}}$ is neglected, and it is calculated by replacing $R_{\mathrm {eff}}$ by $R_A$ in (A4). The long-time limit (blue dashed curve) physically means that $I_{\|,A}$ is neglected, and it is calculated using (A4). (Note that in the case of a cold cloud with a fast magnetic field diffusion time scale compared with the drift motion, in the long-time limit, the $I_{\|,A}=0$ limit is expected to be accurate, as discussed at the end of § 2.1.1.) In addition, we also show results calculated using the simplified expression (A6) (red dash-dotted), that represents the high-background-temperature asymptotic behaviour of the long-time limit.

Figure 4. Comparison of the drift displacement obtained with different limits and model assumptions, for a pellet consisting of (a) $100\,\%$ deuterium and (b) a mixture with $98\,\%$ neon and $2\,\%$ deuterium. Solid black: $I_{\|,A}+I_{\|,\mathrm {ohm}}$, numerical integration of (A1). Dashed blue: $I_{\|,A}=0$, using (A4). Dash-dotted red: $I_{\|,A}=0$ and taking $T_{\rm bg}\rightarrow \infty$ asymptotic behaviour, using (A6). Long dashed green: $I_{\|,\mathrm {ohm}}=0$, using (A4), but with $R_{\rm eff}$ replaced by $R_{A}$.

We see that for both the pure deuterium and the neon-doped pellet, figure 4(a,b), the long-time limit gives similar drift displacement to the general expression (compare dashed and solid), especially at high background-plasma temperatures. There is a discrepancy of $\lesssim 50\,\%$ at background temperatures of $T_{\rm bg}\sim 100\,{\rm eV}$ where the ohmic conductivity is rather low, but at these temperatures the displacement, and therefore the discrepancy, remains moderate. The overall good agreement reflects that the number of connections $N$ continuously increases with time in a hot cloud, so that the Alfvén conductivity is replaced by ohmic conductivity over a short period of time compared with the total drift time.

In the case of a pure deuterium pellet, figure 4(a), we see that the high-background-temperature asymptotic form of the long-time limit (dash-dotted) approaches the more accurate expression (A4) at $T_{\rm bg}\gtrsim 1\,{\rm keV}$, but the approach is much slower in the doped-pellet case, figure 4(b). This difference is due to the higher cloud temperature for a pure deuterium cloud, leading to a longer pressure equilibration time $t_{\mathrm {pe}}$ while the acceleration time $t_{\mathrm {acc}}$ remains only weakly affected by the background temperature, making the approximation $t_{\mathrm {acc}}/t_{\mathrm {pe}}\approx 0$ accurate at lower temperatures.

Finally, we find that the short-time limit (long-dashed curves in figure 4) typically gives unphysically large drift displacements, unlike the general expression and the long-time limit. Only at $T_{\mathrm {bg}}\lesssim 100\,{\rm eV}$ does the short-time-limit expression become comparable to or smaller than the long-time limit; then the ohmic conductivity of the background plasma becomes so low that the Alfvén conductivity starts to dominate. We note that at sufficiently low values of $T_{\mathrm {bg}}$, the short-time limit result starts to asymptotically approach the general expression (black curve), but that happens at very small, inconsequential, values of the drift displacement $\Delta r$.

4. Discussion and conclusion

We have derived a semi-analytical model for the cross-field drift of an ionised cloud following a pellet injection in a tokamak. The model gives the radial drift velocity in terms of the background plasma and cloud properties, assuming the latter to be constant along the field lines inside the cloud. The main phenomena included in the model are the $\boldsymbol {\nabla } B$ current causing the charge separation inside the cloud and the resulting $\boldsymbol{E}\times \boldsymbol{B}$ drift, the rotational transform, pressure equilibration and the currents limiting the charge separation; the latter including the polarisation current and the currents exiting through the ends of the cloud parallel to the field lines, consisting of an Alfvénic and an ohmic contribution. In particular, we have developed a statistical model for the length of the field lines connecting the two ends of the cloud, and the corresponding effective resistivity for the Ohmic current flowing along those field lines.

We then derive semi-analytical expressions for the final drift displacement, combining our model for the cross-field drift with a simple analytical model for the cloud properties. We evaluate the resulting expressions in an ITER-like scenario similar to those studied by Matsuyama (2022), including a wide range of background-plasma temperatures and different neon–deuterium mixtures for the pellet composition. Our results are in reasonable agreement with those obtained by Matsuyama (2022) in the corresponding limit, integrating only up to $130\,\mathrm {\mu }{\rm s}$ after the cloud is detached from the pellet source and neglecting the ohmic part of the parallel current (corresponding to a cold cloud shortly after the pellet material has started to flow across a given field line). We then investigate the effect of adding the ohmic part of the parallel current and integrating to longer times. Without ohmic currents, the final drift displacement becomes unreasonably long, up to several tens (or even hundreds) of metres, while adding the ohmic current reduces the drift displacement by typically 1–2 orders of magnitude.

Our results suggest that a pure deuterium pellet injection in an ITER-like scenario is likely to be significantly affected by the radial drift displacement, and that a substantial part of the injected material may be expelled from the plasma. On the other hand, a neon-doped pellet injection will likely be significantly less affected by the drift displacement. This result corroborates the conclusion made by Matsuyama (2022).

Note, however, that even a relatively small drift displacement can have a significant effect on the ablation and density profile Vallhagen (2021). The reason is that even a small drift means that the pellet will not feel its own cooling effect on the background plasma, which otherwise provides a self-regulating feedback mechanism which decreases the ablation rate. Even a small drift therefore makes the pellet, or pellet shards, ablate faster, so that they deposit more of their material earlier along their trajectories. This applies especially to injections from the low field side, as in that case the drift will displace the ablated material behind the ablating source. On the other hand, an injection from the high-field side will displace the ablated material in front of the pellet or pellet shard, so that it feels the effect of its own cooling along its trajectory. The importance of this effect also depends on the magnetic field strength, which regulates the transverse dimension of the pellet cloud, and on the injection velocity of the pellet. The effect in question is most important if the field is strong and the pellet velocity is small.

In the case of an SPI in an ITER-like scenario, the plume of shards typically extends over several decimetres. Thus, in the case of a neon-doped pellet, our results indicate that the shards will still feel the cooling of the background plasma from most shards ahead of them, even for an injection from the low-field side. For a pure deuterium SPI, on the other hand, the drift displacement will likely be longer than the extent of the plume of shards, which might increase the ablation significantly, especially for an injection from the low-field side. A quantitative assessment of the effect of the drift displacements calculated by the model presented here would require coupling to a model for the full injection dynamics and response of the background plasma, which is outside the scope of the present work.

The accuracy of the results presented in this paper is also limited by a number of simplifications, primarily in the model for the parallel expansion and cloud properties. In particular, the cloud properties are assumed to be constant along the field lines inside the cloud, and the energy balance and temperature evolution is modelled using only a constant, representative value for the cloud temperature. While the cloud temperature remains rather low and constant for a neon-doped pellet due to the high radiated power, the temperature will vary significantly during the drift motion for a pure deuterium pellet; indeed, the discrepancy compared with the results obtained by Matsuyama (2022) is larger for a pure deuterium pellet. The quantitative accuracy of the present model could therefore be significantly improved by combining the present model for the cross-field drift with a more advanced model for the cloud properties, which is outside the scope of the present work.