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. (Reference Commaux, Pégourié, Baylor, Köchl, Parks, Jernigan, Géraud and Nehme2010), 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))$Footnote ¹ 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. (Reference Pégourié, Waller, Nehme, Garzotti and Géraud2006). 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. Reference Commaux, Pégourié, Baylor, Köchl, Parks, Jernigan, Géraud and Nehme2010; Sakamoto et al. Reference Sakamoto, Pégourié, Clairet, Géraud, Gil, Hacquin and Köchl2013). 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 Reference Helander2014), 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. (Reference Pégourié, Waller, Nehme, Garzotti and Géraud2006). 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.Footnote ² 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.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 Reference Parks2017). This expression is based on a version of the neutral gas shielding model (Parks & Turnbull Reference Parks and Turnbull1978) 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. Reference Aleynikov, Breizman, Helander and Turkin2019; Arnold, Aleynikov & Helander Reference Arnold, Aleynikov and Helander2021; Runov et al. Reference Runov, Aleynikov, Arnold, Breizman and Helander2021). 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 Reference Matsuyama2022). 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 (Reference Matsuyama2022), 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. (Reference Kong, Nardon, Hoelzl, Bonfiglio, Hu, Sheikh, Boboc, Carvalho, Hender and Jachmich2022) 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 Reference Parks2016), 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 (Reference Matsuyama2022) 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. Reference Kong, Nardon, Hoelzl, Bonfiglio, Hu, Sheikh, Boboc, Carvalho, Hender and Jachmich2022), 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 (Reference Matsuyama2022). 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 Reference Matsuyama2022) 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 (Reference Matsuyama2022) uses a model similar to that used by Pégourié (Reference 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 (Reference Matsuyama2022) 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 (Reference Matsuyama2022) 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.Footnote ³ 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 (Reference Matsuyama2022), and the results are similar to those shown in figure 11 in Matsuyama (Reference Matsuyama2022) 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 (Reference Matsuyama2022).

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$.