Dynamics of ultrarelativistic charged particles with strong radiation reaction.
II. Entry into Aristotelian equilibrium

We now summarize Gruzinov’s original arguments for Aristotelian motion [4, 15] using the notation of our paper. Gruzinov considered the case where a particle moves primarily along a PND, i.e.,

	$\displaystyle\gamma$	$\displaystyle\gg 1,$		(15)
	$\displaystyle\sqrt{v_{n}^{2}+v_{k}^{2}}$	$\displaystyle\ll 1.$		(16)

If we approximate the particle motion as a circle of radius equal to the radius of curvature $R=1/\kappa$ of the PND, the power radiated is $(2/3)q^{2}R^{2}\gamma^{4}$. Gruzinov assumed that this “curvature radiation” power is balanced by the Lorentz force power $|q|E_{0}$,

$\displaystyle\frac{2}{3}q^{2}R^{2}\gamma^{4}=|q|E_{0}.$

(17)

Solving for $\gamma$ gives the Gruzinov formula for the Lorentz factor,

$\displaystyle\gamma_{g}$

$\displaystyle=\left(\frac{3E_{0}R^{2}}{2|q|}\right)^{1/4}.$

(18)

Gruzinov obtained an additional condition for this equilibrium based on the idea that it can only occur when the particle has enough space to be accelerated to this terminal Lorentz factor before the PND curves significantly. The curvature radius $R$ sets the scale over which the uniform field approximation breaks down, so the particle must be able to gain $m\gamma_{g}$ energy in a region much smaller than $R$. The typical energy gain over a region of size $D$ is $qE_{0}D$, so we obtain the condition

$\displaystyle m\gamma_{g}\ll qE_{0}R,$

(19)

which may equivalently be written

$\displaystyle\gamma_{g}\ll\frac{E_{0}R}{\mathcal{E}\mathcal{R}}$

(20)

or

$\displaystyle\mathcal{E}^{3/2}\mathcal{R}\ll E_{0}^{3/2}R.$

(21)

In paper I [13] we sought to better understand the Aristotelian equilibrium by studying the fundamental LL equation of charged particle dynamics. We again assumed motion along a PND [Eqs. (15) and (16)]. However, the energy balance condition (17) is inappropriate in this context since (1) it requires the particle motion to be treated as circular, an uncontrolled approximation whose compatibility with the LL equation is not obvious; and (2) the local LL dynamics contains more information than just energy conservation. Instead, to express the idea of energy balance, we assumed that the local change in energy is small compared to the typical value set by the Lorentz force,

$\displaystyle m\left|\frac{d\gamma}{dt}\right|$

$\displaystyle\ll|q|E_{0}.$

(22)

We also assumed that the timescale $T$ for changes in the field is long compared to the lengthscale $L$ for spatial changes in the field,

$\displaystyle T\gg L.$

(23)

We found that Gruzinov’s equilibrium emerged only after a final additional assumption,

$\displaystyle|\iota|$

$\displaystyle\ll\frac{|q|}{m\gamma}\textrm{Max}\{E_{0},|B_{0}|\},$

(24)

where $\iota$ is the torsion of the PND. These three conditions (22), (23) and (24) replace the assumption (17) of global power balance in approximate circular motion. Eq. (22) guarantees approximate power balance locally at the level of the LL equation (as opposed to globally, at the level of total radiated energy), while Eqs. (23) and (22) reflect the approximate circular motion with radius $R$.

Under these five assumptions (15), (16), (22), (23), (24), we found that the field determines the velocities pointwise as

	$\displaystyle\gamma$	$\displaystyle=\left(\frac{9}{4}\frac{R^{2}}{\mathcal{R}^{2}}\frac{E_{0}}{\mathcal{E}}\right)^{1/4}=\gamma_{g}$		(25)
	$\displaystyle v_{n}$	$\displaystyle=-\frac{1+\delta}{\gamma}\sqrt{\frac{E_{0}}{\mathcal{E}}}$		(26)
	$\displaystyle v_{k}$	$\displaystyle=\frac{\delta}{\gamma}\frac{B_{0}}{E_{0}}\sqrt{\frac{E_{0}}{\mathcal{E}}},$		(27)

where we used the definitions (13) and also introduced

$\displaystyle\delta=\frac{E_{0}\mathcal{E}}{E_{0}^{2}+B_{0}^{2}}.$

(28)

We thus reproduced Gruzinov’s Lorentz factor (18) and derived new expressions (26) and (27) for the drift velocities.

The five assumptions (15), (16), and (22)–(24) can be recast as conditions purely on the field configuration by using the results (25)–(27). Presenting the results as five conditions $C_{i}\ll 1$, the conditions are

	$\displaystyle C_{1}$	$\displaystyle=\frac{\mathcal{R}}{R}\sqrt{\frac{\mathcal{E}}{E_{0}}}\ll 1$		(29)
	$\displaystyle C_{2}$	$\displaystyle=\delta\frac{\mathcal{R}}{R}\sqrt{\frac{\mathcal{E}}{E_{0}}}\ll 1$		(30)
	$\displaystyle C_{3}$	$\displaystyle=|\iota|\sqrt{R\mathcal{R}}\left(\frac{E_{0}}{\mathcal{E}}\right)^{1/4}\frac{\mathcal{E}}{\textrm{Max}\{E_{0},|B_{0}|\}}\ll 1$		(31)
	$\displaystyle C_{4}$	$\displaystyle=\eta\sqrt{\frac{\mathcal{R}}{R}}\left(\frac{\mathcal{E}}{E_{0}}\right)^{3/4}\ll 1,$		(32)
	$\displaystyle C_{5}$	$\displaystyle=\frac{L}{T}\ll 1$		(33)

with

$\displaystyle\eta=\left|\vec{v}\cdot\vec{\nabla}R+\frac{R}{2E_{0}}\vec{v}\cdot\vec{\nabla}E_{0}\right|.$

(34)

In this last expression, it is understood that Eqs. (25)–(27) [as well as (7)] are to be used to express $\vec{v}$ in terms of the local field configuration. In obtaining (30) we have used the fact that $\sqrt{v_{n}^{2}+v_{k}^{2}}\approx\sqrt{\delta}/\gamma$ under the assumptions (14).

The names $C_{i}$ are derived from paper I. In this paper, $C_{1}$–$C_{5}$ correspond (respectively) to Eqs. (15), (16), (24), (22), and (23). In particular, $C_{1}\ll 1$ and $C_{2}\ll 1$ express the approximate PND motion (15) and (16), while $C_{4}\ll 1$ expresses the approximate local equilibrium (22).

Our expression for $C_{4}$ differs from that presented in paper I, where we assumed that $\eta\sim R/L$. Generically we will have $\eta\sim 1$. In this paper we also include cases where $\eta$ is very small, since they arise quite naturally in our numerical studies, and help to illustrate the differences between our approach and that originally given by Gruzinov (Sec. II.1).

Notice that $C_{1}$ and $C_{2}$ are related to $C_{4}$ by the equations

	$\displaystyle\eta C_{1}^{3/2}$	$\displaystyle=\frac{\mathcal{R}}{R}C_{4}\ll C_{4}$		(35)
	$\displaystyle\eta C_{2}$	$\displaystyle\leq\sqrt{\frac{\mathcal{R}}{R}}C_{4}\ll C_{4}.$		(36)

In the second line we used the bound $\delta\leq\mathcal{E}/E_{0}$, which follows from the definition (28) of $\delta$. In the generic case that $\eta\sim 1$, we see that both $C_{1}$ and $C_{2}$ are small compared with $C_{4}$. This means that $C_{4}\ll 1$ in fact implies that $C_{1}\ll 1$ and $C_{2}\ll 1$, and generically one can ignore the $C_{1}$ and $C_{2}$ conditions. While $\eta$ can be a small number in special cases (such as the circular fields studied in Sec. VI.1 below), it would have to be less than $\sim\sqrt{\mathcal{R}/R}$ for the $C_{1}$ and $C_{2}$ conditions to become relevant. For macroscopic fields this factor is very small; for example, $\sqrt{\mathcal{R}/R}\approx 5\times 10^{-8}$ if $R$ is one meter. A field configuration with $\eta$ this small would be extraordinarily fine-tuned. We therefore conclude that for the macroscopic fields of relevance in astrophysics, we may always ignore the $C_{1}$ and $C_{2}$ conditions, since they are implied by the $C_{4}$ condition.

Let us therefore focus on the $C_{4}$ condition (32). This condition can be rewritten as

$\displaystyle\eta^{2}\mathcal{R}\mathcal{E}^{3/2}\ll RE_{0}^{3/2},$

(37)

showing that it is in fact equivalent to Gruzinov’s condition (21) in the generic case $\eta\sim 1$. This agreement is interesting since Gruzinov’s condition (21) arose from reasoning about when a particle can enter equilibrium, whereas our condition (37) arose from the assumption (22) that equilibrium had been achieved. We will see numerically in Sec. VI that Eq. (37) is necessary and sufficient for equilibrium, provided the other conditions are satisfied.

Finally we discuss $C_{3}$, which satisfies

$\displaystyle C_{3}$

$\displaystyle=\frac{|\iota|}{\kappa}\eta^{-1}C_{4}\frac{E_{0}}{\textrm{Max}\{E_{0},|B_{0}|\}}.$

(38)

The last factor $E_{0}/\textrm{Max}\{E_{0},|B_{0}|\}$ is always $\lesssim 1$ and will be small for pulsars, where $E_{0}\ll|B_{0}|$. Again assuming the generic case $\eta\sim 1$, Eq. (38) shows that large values of the torsion-to-curvature ratio $|\iota|/\kappa$ are required for the condition $C_{3}\ll 1$ to be violated in a regime where the condition $C_{4}\ll 1$ is satisfied. For fields with $|\iota|/\kappa\lesssim 1$, the $C_{3}$ condition can be ignored, since it is implied by the $C_{4}$ condition. We discuss cases with large $|\iota|/\kappa$ in Sec. IV below.

The conditions $C_{i}\ll 1$ arose in paper I as a minimal set of assumptions under which the LL equation reduced to the simple results (25)–(27). Although we do not attempt any rigorous mathematical proof, it seems clear from the steps of the derivation there that these assumptions were all necessary for the simple results to emerge. We therefore regard the conditions (29)–(33) as necessary conditions for the equilibrium described by (25)–(27). One of the main goals of this paper is to argue that they are also sufficient for equilibrium to occur.

At least for some short time, the motion of a particle can be determined in the approximation that the field is uniform. In this approximation, we can always work in a frame where $\vec{E}$ and $\vec{B}$ are parallel. The full analytic solution to the LL equation was found in this case by [16]. We denote the initial velocity by $(v_{1},v_{2},v_{3})$ with initial Lorentz factor $\gamma_{0}$, such that the initial four-velocity is given by

$\displaystyle u^{\alpha}(0)=\gamma_{0}(1,v_{1},v_{2},v_{3}).$

(57)

Taking the field to be in the $z$ direction, the analytic solution to the LL equation is

	$\displaystyle\gamma(\tau)$	$\displaystyle=\frac{1}{2}\frac{(1+v_{3})+(1-v_{3})e^{-2\tau/\tau_{E}}}{\sqrt{1-v_{3}^{2}-(v_{1}^{2}+v_{2}^{2})e^{-2\tau/(\delta\tau_{E})}}}e^{\tau/\tau_{E}}$		(58)
	$\displaystyle v_{x}(\tau)$	$\displaystyle=A(\tau)e^{-\tau/\tau_{\perp}}\sin(\tau/\tau_{B}+\phi)$		(59)
	$\displaystyle v_{y}(\tau)$	$\displaystyle=-\text{sign}(qB_{0})A(\tau)e^{-\tau/\tau_{\perp}}\cos(\tau/\tau_{B}+\phi)$		(60)
	$\displaystyle v_{z}(\tau)$	$\displaystyle=\text{sign}(q)\frac{(1+v_{3})-(1-v_{3})e^{-2\tau/\tau_{E}}}{(1+v_{3})+(1-v_{3})e^{-2\tau/\tau_{E}}}.$		(61)

We are using Cartesian coordinates with $v_{x}(\tau=0)=v_{1}$, $v_{y}(\tau=0)=v_{2}$, and $v_{z}(\tau=0)=v_{3}$, and $\phi$ is the initial direction of motion in the $xy$ plane ($\tan\phi=v_{2}/v_{1}$). Here $A(\tau)$ is defined by

$\displaystyle A(\tau)$

$\displaystyle=\frac{2\sqrt{v_{1}^{2}+v_{2}^{2}}}{(1+v_{3})+(1-v_{3})e^{-2\tau/\tau_{E}}},$

(62)

and we introduced three time scales

	$\displaystyle\tau_{B}$	$\displaystyle=\frac{m}{|qB_{0}|},$		(63)
	$\displaystyle\tau_{E}$	$\displaystyle=\frac{m}{|q|E_{0}}$		(64)
	$\displaystyle\tau_{\perp}$	$\displaystyle=\frac{m}{|q|E_{0}}\frac{\delta}{\delta+1}.$		(65)

The definition of $\delta$ was given above in Eq. (28).

We see that the particle asymptotically approaches light speed in the $z$ direction, i.e., it eventually moves along the PND. We interpret this to mean that radiation damps only the perpendicular momentum. The particle executes a damped circular motion in the $xy$ plane (around the field direction), in accordance with the low-$E_{0}$ intuition of synchrotron motion and associated radiation damping. The period of oscillation is just the classical synchrotron period $\tau_{B}$ (expressed here in proper time), while the decay time $\tau_{\perp}$ generalizes the classical synchrotron damping result. In particular, we have

$\displaystyle\tau_{\perp}\approx\begin{cases}\tau_{B}\frac{\mathcal{E}}{|B_{0}|}&|B_{0}|\gg E_{0}\\ \tau_{E}&E_{0}\gg|B_{0}|\end{cases},$

(66)

taking into account $E_{0}\ll\mathcal{E}$. The former case agrees in order of magnitude with Eq. (A2) of Ref. [17].

To understand the dynamics it is convenient to expand the Lorentz factor at early and late times,

$\displaystyle\gamma(\tau)\approx\begin{cases}\gamma_{0}\left(1-\tau/\tau_{\rm drop}-\frac{v_{3}}{\gamma_{0}^{2}}\tau/\tau_{E}\right),&\tau\to 0\\ \frac{1}{2}\sqrt{\frac{1+v^{3}}{1-v_{3}^{2}}}e^{\tau/\tau_{E}},&\tau\to\infty.\end{cases}$

(67)

where we introduce yet another timescale

$\displaystyle\tau_{\rm drop}\equiv\frac{\delta}{\gamma_{0}^{2}(v_{1}^{2}+v_{2}^{2})}\tau_{E}.$

(68)

For sufficiently large initial velocity perpendicular to the field, $\gamma_{0}(v_{1}^{2}+v_{2}^{2})\gg\delta$, we have $\tau_{\rm drop}\ll\tau_{\rm E}$. In this case there is a sudden drop in Lorentz factor on timescale $\tau_{\rm drop}$ as the particle loses its perpendicular momentum, followed by a subsequent rise on timescale $\tau_{E}$ as the particle gains parallel momentum. The minimal Lorentz factor is determined by the details of (58); in the special case $\delta\gg 1$, the minimum occurs at approximately $\sqrt{\delta}$.

In a realistic field configuration, the exponential rise in Lorentz factor is ultimately cut off by the effects of non-uniform fields. If the conditions are right for Aristotelian equilibrium, there will be a transition around the equilibrium Lorentz factor. In order to understand this transition, we now study the LL equation near the Aristotelian equilibrium solution.

To study the dynamics near equilibrium, we adopt the same assumptions as of paper I, assuming time derivatives are perturbatively small (instead of dropping them entirely). Following the steps of Sec. IVB of that reference, except retaining time derivatives,³³3We also drop all terms involving $\iota$, a step that was justified in a later section of paper I. When time derivatives are dropped, Eqs. (69), (70), and (71) reduce, respectively, to Eqs.  (43), (54) and (55) of paper I with $\iota=0$. we arrive at three coupled equations for $\gamma$, $v_{n}$ and $v_{k}$

	$\displaystyle\frac{d\gamma}{d\tau}$	$\displaystyle=\frac{|q|\gamma E_{0}}{m}-\frac{|q|\gamma^{3}}{m\mathcal{E}}(E_{0}^{2}+B_{0}^{2})(v_{n}^{2}+v_{k}^{2})$		(69)
	$\displaystyle\frac{dv_{n}}{d\tau}$	$\displaystyle=\frac{|q|B_{0}}{m}v_{k}-\frac{|q|E_{0}}{m}\frac{1+\delta}{\delta}v_{n}-\gamma\kappa$		(70)
	$\displaystyle\frac{dv_{k}}{d\tau}$	$\displaystyle=-\frac{|q|B_{0}}{m}v_{n}-\frac{|q|E_{0}}{m}\frac{1+\delta}{\delta}v_{k}.$		(71)

We now express these quantities as small perturbations of their equilibrium values,

	$\displaystyle\gamma$	$\displaystyle=\Gamma+\bar{\gamma}$		(72)
	$\displaystyle v_{n}$	$\displaystyle=V_{n}+\bar{v}_{n}$		(73)
	$\displaystyle v_{k}$	$\displaystyle=V_{k}+\bar{v}_{k},$		(74)

where $\Gamma,V_{n},V_{k}$ are taken to be Eqs. (25)–(27), respectively. Regarding $E_{0},B_{0},\kappa$ as constants and linearizing in the barred quantities, we find

	$\displaystyle\tau_{E}\kappa\frac{d\bar{\gamma}}{d\tau}$	$\displaystyle=-2\kappa\bar{\gamma}+2\bar{v}_{n}/\tau_{\perp}-\epsilon 2\bar{v}_{k}/\tau_{B}$		(75)
	$\displaystyle\frac{d\bar{v}_{n}}{d\tau}$	$\displaystyle=-\bar{v}_{n}/\tau_{\perp}+\epsilon\bar{v}_{k}/\tau_{B}-\bar{\gamma}\kappa$		(76)
	$\displaystyle\frac{d\bar{v}_{k}}{d\tau}$	$\displaystyle=-\epsilon\bar{v}_{n}/\tau_{B}-\bar{v}_{k}/\tau_{\perp},$		(77)

where the timescales $\tau_{i}$ were introduced in Eqs. (63)–(65) and we have also written $\epsilon=\textrm{sign}(B_{0})$. Defining

	$\displaystyle T$	$\displaystyle=\frac{\tau}{\tau_{E}}$		(78)
	$\displaystyle b$	$\displaystyle=\epsilon\frac{\tau_{E}}{\tau_{B}}=\frac{B_{0}}{E_{0}}$		(79)
	$\displaystyle c$	$\displaystyle=\frac{\tau_{E}}{\tau_{\perp}}=\frac{1+\delta}{\delta}>1,$		(80)

this set of equations may also be written

$\displaystyle\frac{d\bar{X}}{dT}=M\bar{X},$

(81)

with

$\displaystyle\bar{X}=\begin{pmatrix}\bar{\gamma}\kappa\tau_{E}\\ \bar{v}_{n}\\ \bar{v}_{k}\end{pmatrix},\quad M=\begin{pmatrix}-2&2c&-2b\\ -1&-c&b\\ 0&-b&-c\\ \end{pmatrix}.$

(82)

Making the ansatz

$\displaystyle\bar{X}=\bar{X}_{0}e^{\lambda T},$

(83)

presents the eigenvalue problem $M\bar{X}_{0}=\lambda\bar{X}_{0}$ (here $\bar{X}_{0}$ is a constant independent of $T$). If $\lambda^{(i)}$ are the eigenvalues and $\bar{X}^{(i)}_{0}$ the associated eigenvectors (for $i=1,2,3$), then the general solution is

$\displaystyle\bar{X}(T)=\sum_{i}C_{i}\textrm{Re}[\bar{X}_{0}^{(i)}e^{\lambda^{(i)}T}],$

(84)

for three real constants $C_{i}$. The eigenvalues and eigenvectors depend only on $E_{0}$ and $B_{0}$ and can be found numerically for any given values of $E_{0}$ and $B_{0}$.

We now show that the equilibrium is linearly stable, i.e., $\textrm{Re}[\lambda^{(i)}]<0$. The eigenvalues $\lambda^{(i)}$ are the roots of the characteristic polynomial (defined with a minus sign for convenience),

		$\displaystyle f(\lambda)=-\textrm{det}(M-\lambda I)$
		$\displaystyle=\lambda^{3}+2(c+1)\lambda^{2}+(b^{2}+6c+c^{2})\lambda+4(b^{2}+c^{2}).$		(85)

Notice that all the coefficients are real and positive. Thus $f$ is strictly positive for $\lambda\geq 0$, implying that any real root is strictly negative. We have therefore proven stability in the case of real roots only. For complex roots, note that such roots must appear in a complex conjugate pair, since we consider a real polynomial. We may therefore write

$\displaystyle f(\lambda)=(\lambda-\lambda_{r})(\lambda-\alpha-i\beta)(\lambda-\alpha+i\beta),$

(86)

with $\lambda_{r},\alpha,\beta$ all real. Comparing the coefficient of $\lambda^{2}$ with the polynomial (V.2), we find that

$\displaystyle 2\alpha=-\lambda_{r}-2(c+1).$

(87)

However, we find that $f$ is negative at $\lambda=-2(c+1)$:

$\displaystyle f(-2(c+1))=-2\left[(c+3)(c+2)c+b^{2}(c-1)\right]<0,$

(88)

noting from Eq. (80) that $c>1$. Since $f$ is positive for $\lambda\geq 0$, there must be a real root between $\lambda=-2(c+1)$ and $\lambda=0$. In the case (86) we consider, $\lambda_{r}$ is the single real root, so we have

$\displaystyle-2(c+1)<\lambda_{r}<0.$

(89)

It then follows from (87) that $\alpha$ is strictly negative, completing the proof that $\textrm{Re}[\lambda^{(i)}]<0$.

Perturbations away from equilibrium are thus exponentially damped. The qualitative behavior of approach to equilibrium is determined by the mode (or pair of complex conjugate modes) with smallest $|\textrm{Re}[\lambda^{(i)}]|$, which decays the slowest. If this mode has a non-zero imaginary part (i.e., it is part of a complex-conjugate pair), then oscillations will accompany the decay, with frequency equal to $|\textrm{Im}[\lambda^{(i)}]|$. Note also that the properties of the decay are insensitive to the sign of $b$, since the characteristic polynomial (V.2) depends only on $b^{2}$.

It is interesting to ask for what parameter ranges such oscillations will be important. The discriminant of the cubic (V.2) is equal to

	$\displaystyle\Delta$	$\displaystyle=-\frac{4}{27}[(3b^{2}-c^{2}+10c-4)^{3}+$		(90)
		$\displaystyle(c^{3}-15c^{2}+30c+9b^{2}c-45b^{2}-8)^{2}]$		(91)

When $\Delta<0$ there are complex-conjugate roots; otherwise all roots are real. Both signs occur over the parameter ranges $b\in\mathbb{R},c>1$, so both behaviors are possible. If $\Delta\geq 0$, there will be no oscillations. If $\Delta<0$ there will be oscillations, and these will be the dominant behavior if $\alpha<\lambda_{r}$ in the notation of (86).

One tractable case is when $\delta\gg 1$ so that $c\approx 1$. It is easy to see that the discriminant is negative for $c=1$, so that there will be complex-conjugate roots. One can then check that these oscillatory modes decay more slowly than the pure exponential mode provided $|b|\geq 2/\sqrt{3}$, and that they decay with a similar order of magnitude even for $|b|\leq 2/\sqrt{3}$. The oscillations will thus always be important in the case $\delta\gg 1$.

Some further details are helpful for interpreting this case. As a function of $b$, the dominant decay time (real part of eigenvalue with real part closest to zero) ranges between $-1$ at $b=0$ and $-4$ at $|b|\to\infty$, and the frequency of oscillations (imaginary part of complex eigenvalue) ranges between (approximately) $1.32$ at $b=0$ and infinity at large $|b|$, approaching linearly like $|b|$ as $b\to\pm\infty$. Back in the physical time coordinate $\tau=\tau_{E}T$, we see that decay will typically take several to tens of $\tau_{E}$, and the oscillation angular frequency will be at least $1.32\tau_{E}^{-1}$, and will be approximately $|b|\tau_{E}^{-1}=\tau_{B}^{-1}$ at large $|b|$. It is notable that these oscillations exist even in the pure-E case ($\delta\gg 1,b=0$) and are thus physically distinct from synchrotron motion in that limit.

We now show two numerical examples illustrating the features explored in the previous subsections. We consider the case of parallel circular electric and magnetic fields, i.e., Eqs. (52) and (53) with $h=0$. Since the electric and magnetic fields are parallel in the lab frame, no boost is required to relate the parallel-frame uniform field solution (Sec. V.1) to the lab-frame near-equilibrium solution (Sec. V.2). We always place the particle in a region of the circular field configuration where $C_{4}\ll 1$ and hence equilibrium is expected to occur.

First we consider the case where the particle begins with a large momentum perpendicular to the PND. Based on the analysis of the previous subsections, we expect the particle to sharply lose its perpendicular momentum on a timescale $\tau_{\rm drop}$ (68), then gain parallel momentum exponentially on a timescale $\tau_{E}$ (64) until it nears the equilibrium Lorentz factor, and finally approach equilibrium exponentially, with timescale and possible oscillations determined from $\tau_{E},\tau_{B},\tau_{\perp}$ via the eigenvalue problem discussed in Sec. V.2. These expectations are confirmed by our numerical experiments; an example is shown in Fig. 4.

We next consider the case where the particle begins with a large momentum $\gamma\gg\gamma_{g}$ parallel to the PND. Here the uniform field approximation is not useful since the particle quickly reaches the region where the field line bends away from its motion. However, at this stage we can regard the particle as traveling through a new uniform field configuration with some perpendicular momentum, which will be removed by radiation reaction according to the intuition discussed in the first new paragraph below Eq. (65). In this way the particle will lose energy until it either enters equilibrium directly (“from above”) or finds itself in a uniform field configuration that will accelerate it back up to near-equilibrium conditions, after which it enters equilibrium “from below”. An example somewhat intermediate between these two cases is shown in Fig. 5.

Eq. (34) defined a quantity $\eta$ that appears in the conditions for equilibrium. For generic fields, $\eta\approx 1$ and factors of $\eta$ can be dropped. To illustrate a situation where $\eta$ cannot be dropped, consider a purely azimuthal field configuration with constant parallel electric and magnetic fields. The field lines and PNDs are just circles of radius $\rho=\sqrt{x^{2}+y^{2}}$. Since $\vec{\nabla}\rho=\vec{n}$, where $\vec{n}$ is the Frenet-Serret normal direction to the curve, $\eta$ can be calculated as

$\displaystyle\eta=|\vec{v}\cdot\nabla\rho|=|v_{n}|.$

(92)

This quantity is small in equilibrium since $|v_{n}|\ll 1$ by the assumption of motion along a PND. The small size of $\eta$ arises because the field strength and radius of curvature are both constant along the PND direction $\vec{\ell}$, a finely tuned arrangement. From Eqs. (26) and (25), we have

$\displaystyle|v_{n}|^{2}=\frac{2}{3}\frac{\mathcal{R}}{R}\sqrt{\frac{E_{0}}{\mathcal{E}}}(1+\delta)^{2}.$

(93)

The condition (37) then becomes

$\displaystyle R^{2}E_{0}\gg(1+\delta)^{2}\mathcal{R}^{2}\mathcal{E},$

(94)

where we drop a factor of $2/3$.

Eq. (94) is a necessary condition for equilibrium to occur. To check whether it is also sufficient, we simulated a large number of trajectories with different field configuration parameters and particle initial conditions. Specifically, we fixed $\tilde{B}_{0}=0.1$ and chose the particle to begin on the $x$ axis, choosing the other parameters randomly from uniform distributions in the ranges $\log_{10}\tilde{E}_{0}\in(-4,0)$ for the log of the field strength, $\gamma\in(1,\gamma_{g})$ for the initial Lorentz factor, $\theta\in(0,\pi)$ and $\phi\in(0,2\pi)$ for the initial direction of motion (where $\theta,\phi$ are spherical coordinates on the space of velocities), and $\log_{10}x\in(-3,0)$ for the log of the initial position. The initial radius of curvature $\tilde{R}$ is just $x$, so we sample an initial range $\log_{10}\tilde{R}\in(-3,0)$.

For each run of the code we evolve for a maximum time of $\tau_{\rm max}=(6\log\gamma_{g})\tau_{E}$, which is sufficient for the particle to enter equilibrium if it is destined to do so.⁴⁴4According to Eq. (67), the time for a particle to be accelerated to $\gamma_{g}$ is of order $\tau_{E}\log\gamma_{g}$. For the circular field $E_{0}$ and hence $\tau_{E}$ is constant; however, for other field configurations we update the maximum allowed run time after each time step to use the local value of $E_{0}$. We also examined a selection of the trajectories that did not enter equilibrium in this time and evolved them for longer, finding that indeed they never enter equilibrium. At any given time in the trajectory, the particle is considered to have entered equilibrium if $\langle|(\gamma-\gamma_{g})|/\gamma_{g}\rangle<3\%$, where the average $\langle...\rangle$ is calculated over the last $20\%$ of the computed trajectory, using the local value of $\gamma_{g}$. Once a particle enters equilibrium, we record the local electric field $E_{0}$ and curvature radius $R$ and terminate the run. (If the particle never enters equilibrium, we make no corresponding record.) Fig. 6 shows the values of $E$ and $R$ at which equilibrium occurred in our random sample, showing clearly that the condition (94) indeed controls whether equilibrium is obtained. We find that the left-hand-side must be approximately 15 times larger than the right-hand side (red line in the figure) for equilibrium to occur.

When $E_{0}\gg B_{0}$, we have $\delta\approx\mathcal{E}/E_{0}\gg 1$, and this condition reduces to the original Gruzinov condition (21). (This agreement is accidental, due to the specific form of $\eta$ for this field configuration.) By contrast, when $E_{0}\lesssim B_{0}$, the condition (94) remains distinct from the Gruzinov condition (21). This behavior is seen clearly in Fig. 6. The Gruzinov condition is necessary, but not sufficient, in this special case. The correct condition for entry into equilibrium is the modified condition (37).

The circular field is a special case that cleanly illustrates the role of $\eta$ in the condition $C_{4}\ll 1$ (32) for entry into equilibrium. In order to test the condition more generally, we consider two more kinds of field configuration. The “ellipse” field configuration simply changes the shape of the field lines to ellipses instead of circles, while keeping everything else the same. The “separate center” field configuration consists of circular electric and magnetic field lines with equal field strengths $E_{0}=B_{0}$ but with the center of the circles shifted by a distance $d$. We randomly varied the field strength and center-distance while also randomly choosing initial conditions as before.

The results of these surveys are shown along with the circular case in Fig. 7, demonstrating clearly that the condition $C_{4}\ll 1$ is necessary and sufficient for equilibrium to occur. The precise value of $C_{4}$ required depends on the field configuration, ranging from $\sim 0.1$ to $\sim 0.01$ for the cases considered here.