# Eighth–Order Foldy–Wouthuysen Transformation

Ulrich D. Jentschura  Department of Physics and LAMOR, Missouri University of
Science and Technology, Rolla, Missouri, USA

###### Abstract

The calculation of higher-order binding corrections to bound systems is a
fundamental problem of theoretical physics. For any nonrelativistic expansion,
one needs the Foldy–Wouthuysen Transformation which disentangles the particle
and the antiparticle degrees of freedom. This transformation is carried out
here to eighth order in the momenta, or, to eighth order in the momentum
operators, which is equivalent to the eighth order of the fine-structure
constant. Matrix elements of the eighth-order terms are evaluated for
$F_{5/2}$ and $F_{7/2}$ states in hydrogenlike ions and compared with the
Dirac–Coulomb energy levels.

###### Contents

  1. I Overview
  2. II Direct Calculation
    1. II.1 Unitary Transformation
    2. II.2 General Particle Hamiltonian
  3. III Applications
    1. III.1 Coulomb Field Coupling
    2. III.2 Application: $\maybebm{F_{5/2}}$ States
    3. III.3 Application: $\maybebm{F_{7/2}}$ States
  4. IV Conclusions
  5. A Comparison with the Literature

##  I Overview

The Foldy–Wouthuysen transformation FoWu1950  is one of the most essential
ingredients of the bound-state formalism JeAd2022book . Specifically, the
Foldy–Wouthuysen transformation enables one to disentangle the particle and
antiparticle degrees of freedom, and write separate particle and antiparticle
Hamiltonians for spin-$1/2$ particles coupled to electromagnetic and
gravitational fields JeNo2013pra . the case of a free Dirac particle, the
Foldy–Wouthuysen transformation can be carried out all orders in the momentum
operators, and the transformed Hamiltonian takes the simple form [see Eq.
(11.17) of Ref. JeAd2022book ]

$H_{\mathrm{FW}}=\left(\begin{array}[]{cc}\sqrt{\vec{p}^{\,2}+m^{2}}\cdot\mathbbm{1}_{2\times
2}&0\\\ 0&-\sqrt{\vec{p}^{\,2}+m^{2}}\cdot\mathbbm{1}_{2\times 2}\\\
\end{array}\right)\,,$ (1)

where $\vec{p}$ is the momentum operator. In other cases, where the electron
is subjected to nontrivial couplings to electromagnetic or gravitational
fields JeNo2013pra , it is possible to carry out the Foldy–Wouthuysen
transformation only up to a finite order in a chosen expansion parameter,
which can, in many cases, be chosen as a typical momentum scale of the
physical problem. Here, we are concerned with an electron coupled to a general
electromagnetic field, described by a vector potential $\vec{A}$ and a scalar
potential $\Phi$.

For bound systems, including bound Coulomb systems and electrons bound in a
Penning trap BrGa1986 ; WiMoJe2022 , one can identify the typical momentum
scale as $\alpha\,m\,c$, where $m$ is the electron mass, $c$ is the speed of
light, and $\alpha$ is the fine-structure constant or a generalization
thereof. Henceforth, we shall use natural units with $\hbar=\epsilon_{0}=c=1$,
and $e^{2}=4\pi\alpha$. Specifically, for electrons bound in Penning traps,
one can define a generalized cyclotron fine-structure constant according to
Eq. (35) of Ref. WiMoJe2022 . We define the kinetic momentum

$\vec{\pi}=\vec{p}-e\,\vec{A}\,,$ (2)

where $e$ is the electron charge. For the electric field $\vec{E}$ and the
magnetic field $\vec{B}$, as well as its time derivative
$e\,\partial_{t}\vec{E}$, we assume the following scaling [see Eq. (47) of
Ref. WiMoJe2022 ],

$\vec{\pi}\sim\alpha\,,\qquad e\vec{A}\sim\alpha\,,\qquad
e\,\vec{B}\sim\alpha^{2}\,,\\\ V\equiv e\,\Phi\,,\qquad
e\,\vec{E}=-\vec{\nabla}V\sim\alpha^{3}\,,\qquad
e\,\partial_{t}\vec{E}\sim\alpha^{5}\,.$ (3)

For the fourth-order Foldy–Wouthuysen transformation, a particularly
instructive derivation is presented in Ref. BjDr1964 . The sixth-order
Foldy–Wouthuysen transformation has been considered extensively in the
literature [see Eqs. (36)–(38) of Ref. ZaPa2010 , Ref. HiLePaSo2013 . Eqs.
(15) and (20) of Ref. PaYePa2016 , Eq. (7) of Ref. HaZhKoKa2020 , as well as
Ref. ZhMeQi2023 .

Here, it is our goal to extend the formalism to the eighth order in the fine-
structure constant. We venture to obtain the general Hamiltonian for
electromagnetic coupling in Sec. II, and apply the obtained results to subsets
of bound states in hydrogenlike systems in Sec. III, which leads to a
verification of the results against the analytically known Dirac–Coulomb
energy. Extensive use is made of computer algebra Wo1999 . Conclusions are
reserved in Sec. IV.

##  II Direct Calculation

###  II.1 Unitary Transformation

We start from the well-known Dirac Hamiltonian $H_{\mathrm{D}}$ coupled to a
general electromagnetic field,

$\displaystyle\mathcal{H}_{\rm D}=$
$\displaystyle\;\vec{\alpha}\cdot\vec{\pi}+\beta\,m$ $\displaystyle=$
$\displaystyle\;\left(\begin{array}[]{cc}(m+e\,\Phi)\mathbbm{1}_{2\times
2}&\;\;\;\vec{\sigma}\cdot\vec{\pi}\\\
\vec{\sigma}\cdot\vec{\pi}&\;\;\;(-m+e\,\Phi)\,\mathbbm{1}_{2\times
2}\end{array}\right)\,,$ (6) $\displaystyle\vec{\alpha}=$
$\displaystyle\;\left(\begin{array}[]{cc}0&\vec{\sigma}\\\
\vec{\sigma}&0\end{array}\right)\,,\qquad\beta=\left(\begin{array}[]{cc}\mathbbm{1}_{2\times
2}&0_{2\times 2}\\\ 0_{2\times 2}&-\mathbbm{1}_{2\times
2}\end{array}\right)\,.$ (11)

The Dirac Hamiltonian couples upper and lower components of the Dirac
bispinor. The aim of the (unitary) Foldy–Wouthuysen transformation is to
separate the upper (particle) from the lower (antiparticle) degrees of
freedom, through an iterative procedure. One writes the Hamiltonian as a sum
of even terms ($\mathcal{E}$) and odd terms ($\mathcal{O}$) in bispinor space,

$\displaystyle\mathcal{H}_{\rm D}=$
$\displaystyle\;\mathcal{E}+\mathcal{O}\,,$ (12a) $\displaystyle\mathcal{E}=$
$\displaystyle\;\tfrac{1}{2}\,\left(H+\beta\,H\,\beta\right)\,,\qquad\mathcal{O}=\tfrac{1}{2}\,\left(H-\beta\,H\,\beta\right)\,,$
(12b) $\displaystyle\mathcal{E}=$ $\displaystyle\;\beta m+e\Phi$
$\displaystyle=$
$\displaystyle\;\left(\begin{array}[]{cc}(m+e\,\Phi)\mathbbm{1}_{2\times
2}&\;\;\;0_{2\times 2}\\\ 0_{2\times
2}&\;\;\;(-m+e\,\Phi)\,\mathbbm{1}_{2\times 2}\end{array}\right)\,,$ (12e)
$\displaystyle\mathcal{O}=$
$\displaystyle\;\vec{\alpha}\cdot\vec{\pi}=\left(\begin{array}[]{cc}0&\;\;\;\vec{\sigma}\cdot\vec{\pi}\\\
\vec{\sigma}\cdot\vec{\pi}&0\end{array}\right)\,.$ (12h)

The aim is to eliminate the odd terms $\mathcal{O}$ through unitary
transformations. These transformations, necessarily, in order to preserve the
physical interpretation of the operators, need to conserve parity [see the
remarks following Eq. (7.33) of Ref. ObSiTe2017  and the comprehensive
discussion in Ref. JeNo2014jpa ]. One defines the Hermitian operator $S$ and
the unitary operator $U=\exp(\mathrm{i}S)$ as follows,

$S=-\mathrm{i}\,\beta\,\frac{\mathcal{O}}{2m}\,,\qquad
S=S^{\mathrm{+}}\,,\qquad U=\mathrm{e}^{\mathrm{i}S}\,.$ (13)

The odd operator $\mathcal{O}$ defined in Eq. (12h) is proportional to the
kinetic momentum $\vec{\pi}$, and hence, the Foldy–Wouthuysen eliminates the
odd terms order by order in an expansion in the momenta. The transformed
Hamiltonian is written in terms of nested commutators,

$\displaystyle\mathcal{H}_{\rm FW}=$
$\displaystyle\;\exp(\mathrm{i}S)\,\left(H-\mathrm{i}\,\partial_{t}\right)\exp(-\mathrm{i}S)$
$\displaystyle=$
$\displaystyle\;H+\left[\mathrm{i}\,S,H-\mathrm{i}\,\partial_{t}\right]+\frac{1}{2!}\,\left[\mathrm{i}S,\left[\mathrm{i}\,S,H-\mathrm{i}\,\partial_{t}\right]\right]$
$\displaystyle\;+\frac{1}{3!}\,\left[\mathrm{i}S,\left[\mathrm{i}S,\left[\mathrm{i}\,S,H-\mathrm{i}\,\partial_{t}\right]\right]\right]$
$\displaystyle\;+\frac{1}{4!}\,\left[\mathrm{i}S,\left[\mathrm{i}S,\left[\mathrm{i}S,\left[\mathrm{i}\,S,H-\mathrm{i}\,\partial_{t}\right]\right]\right]\right]+\ldots\,.$
(14)

where one notes the identity
$[\mathrm{i}S,H-\mathrm{i}\partial_{t}]=\mathrm{i}[S,H]-\partial_{t}S$.

Though the intensive use of computer algebra generalized to the symbolic
commutation relations Wo1999  of kinetic momentum operators with the four-
vector and scalar potentials, it is possible to carry out the transformation
through eighth order in the fine-structure constant, under the proviso of the
scaling implied by Eq. (3). The result of the iterative eighth-order
transformation JeNo2014jpa , can be written as follows,

$\mathcal{H}_{\rm
FW}=\mathcal{H}^{[0]}+\mathcal{H}^{[2]}+\mathcal{H}^{[4]}+\mathcal{H}^{[6]}+\mathcal{H}^{[8]}\,.$
(15)

The superscript denotes the power of the coupling parameter at which the term
becomes relevant. The coupling parameter is usually denoted as $\alpha$. From
the zeroth to the third order in $\alpha$, the terms read as follows,

$\displaystyle\mathcal{H}^{[0]}=$ $\displaystyle\;\beta
m\,,\qquad\vec{\Sigma}=\left(\begin{array}[]{cc}\vec{\sigma}&0_{2\times 2}\\\
0_{2\times 2}&\vec{\sigma}\end{array}\right)\,,$ (18)
$\displaystyle\mathcal{H}^{[2]}=$
$\displaystyle\;\frac{\beta}{2m}(\vec{\Sigma}\cdot\vec{\pi})^{2}+V\,.$ (19)

The $\alpha^{4}$ terms can be expressed very compactly, and are found in
agreement with Refs. BjDr1964 ; Pa2005 ,

$\displaystyle\mathcal{H}^{[4]}=$
$\displaystyle\;-\beta\,\frac{1}{8m^{3}}(\vec{\Sigma}\cdot\vec{\pi})^{4}-\frac{\mathrm{i}\,e}{8m^{2}}[\,\vec{\Sigma}\cdot\vec{\pi},\vec{\Sigma}\cdot\vec{E}\,]\,.$
(20)

For the $\alpha^{6}$ terms, we indicate three alternative representations,

$\displaystyle\mathcal{H}^{[6]}=$
$\displaystyle\;\frac{\beta\,(\vec{\Sigma}\cdot\vec{\pi})^{6}}{16m^{5}}-\frac{5\mathrm{i}e}{128m^{4}}[\,\vec{\Sigma}\cdot\vec{\pi},[\,\vec{\Sigma}\cdot\vec{\pi},[\,\vec{\Sigma}\cdot\vec{\pi},\vec{\Sigma}\cdot\vec{E}\,]\,]\,]$
$\displaystyle\;+\frac{\mathrm{i}\,e}{8m^{4}}\;\\{\,(\vec{\Sigma}\cdot\vec{\pi})^{2},[\,\vec{\Sigma}\cdot\vec{\pi},\vec{\Sigma}\cdot\vec{E}\,]\\}+\beta\,\frac{e^{2}\vec{E}^{\,2}}{8m^{3}}$
$\displaystyle=$
$\displaystyle\;\frac{\beta\,(\vec{\Sigma}\cdot\vec{\pi})^{6}}{16m^{5}}+\frac{5\mathrm{i}e}{128m^{4}}\left[(\vec{\Sigma}\cdot\vec{\pi})^{2},\left\\{\vec{\Sigma}\cdot\vec{E},\vec{\Sigma}\cdot\vec{\pi}\right\\}\right]$
$\displaystyle\;+\frac{3\mathrm{i}\,e}{64m^{4}}\;\\{\,(\vec{\Sigma}\cdot\vec{\pi})^{2},[\,\vec{\Sigma}\cdot\vec{\pi},\vec{\Sigma}\cdot\vec{E}\,]\\}+\beta\,\frac{e^{2}\vec{E}^{\,2}}{8m^{3}}\,.$
$\displaystyle=$
$\displaystyle\;\frac{\beta\,(\vec{\Sigma}\cdot\vec{\pi})^{6}}{16m^{5}}+\frac{3\mathrm{i}e}{32m^{4}}\left[\left(\vec{\Sigma}\cdot\vec{\pi}\right)^{3},\vec{\Sigma}\cdot\vec{E}\right]\,.$
$\displaystyle\;-\frac{\mathrm{i}\,e}{128m^{4}}\;[\,(\vec{\Sigma}\cdot\vec{\pi})^{2},\\{\,\vec{\Sigma}\cdot\vec{\pi},\vec{\Sigma}\cdot\vec{E}\,\\}\,]+\beta\,\frac{e^{2}\vec{E}^{\,2}}{8m^{3}}\,.$
(21)

Here, $\\{A,B\\}=A\,B+B\,A$ denotes the anticommutator. The last form of the
sixth-order terms in Eq. (II.1) is in agreement with the sixth-order terms
from Eq. (8) Ref. ZhMeQi2023 . The sixth-order terms are also compatible with
Eqs. (36)–(38) of Ref. ZaPa2010 , with Eq. (7) of Ref. HaZhKoKa2020 , and with
the approach from Ref. HiLePaSo2013 . Alternatively, the result in Eq. (II.2)
can be obtained by applying the unitary transformation outlined in Eq. (19) of
Ref. PaYePa2016  to the Hamiltonian given in Eq. (15) of Ref. PaYePa2016 ,
which is tantamount to the Hamiltonian obtained by adding the terms given in
Eqs. (15) and (20) of Ref. PaYePa2016 .

The eighth-order term are naturally written as a sum of a kinetic term
$\mathcal{K}$, a term $\mathcal{D}$ involving temporal derivatives of the
electric field, terms quadratic in the electric field, denoted by
$\mathcal{Q}$, and linear terms in the electric field, which we denote by
$\mathcal{L}$. The result is as follows,

$\displaystyle\mathcal{H}^{[8]}=$
$\displaystyle\;\mathcal{K}+\mathcal{D}+\mathcal{Q}+\mathcal{L}\,,\qquad\mathcal{K}=-\beta\,\frac{5}{128m^{7}}(\vec{\Sigma}\cdot\vec{\pi})^{8}\,,\qquad\mathcal{D}=-\frac{\mathrm{i}\,e^{2}}{32m^{4}}\,[\,\vec{\Sigma}\cdot\vec{E},\vec{\Sigma}\cdot\partial_{t}\vec{E}\,]+\frac{e}{48m^{5}}\,\beta\,\\{\,(\vec{\Sigma}\cdot\vec{\pi})^{3},\vec{\Sigma}\cdot\partial_{t}\vec{E}\,\\}\,,$
$\displaystyle\mathcal{Q}=$
$\displaystyle\;\frac{7\,\beta\,e^{2}}{192m^{5}}\,[\,\vec{\Sigma}\cdot\vec{\pi},\vec{\Sigma}\cdot\vec{E}\,]\;[\,\vec{\Sigma}\cdot\vec{\pi},\vec{\Sigma}\cdot\vec{E}\,]-\frac{3\,\beta\,e^{2}}{64m^{5}}\,\\{\,\vec{\Sigma}\cdot\vec{\pi},\vec{\Sigma}\cdot\vec{E}\,\\}\;\\{\,\vec{\Sigma}\cdot\vec{\pi},\vec{\Sigma}\cdot\vec{E}\,\\}-\frac{\beta
e^{2}}{24m^{5}}\,[\,\vec{\Sigma}\cdot\vec{\pi},[\,\vec{\Sigma}\cdot\vec{\pi},(\vec{\Sigma}\cdot\vec{E})^{2}\,]\,]\,,$
$\displaystyle\mathcal{L}=$
$\displaystyle\;-\frac{5\,\mathrm{i}\,e}{1024m^{6}}[\,\vec{\Sigma}\cdot\vec{\pi},[\,\vec{\Sigma}\cdot\vec{\pi},[\,\vec{\Sigma}\cdot\vec{\pi},[\,\vec{\Sigma}\cdot\vec{\pi},[\,\vec{\Sigma}\cdot\vec{\pi},\vec{\Sigma}\cdot\vec{E}\,]\,]\,]\,]\,]-\frac{\mathrm{i}\,e}{32m^{6}}\\{\,\vec{\Sigma}\cdot\vec{\pi},\\{\,\vec{\Sigma}\cdot\vec{\pi},[\,\vec{\Sigma}\cdot\vec{\pi},[\,\vec{\Sigma}\cdot\vec{\pi},[\,\vec{\Sigma}\cdot\vec{\pi},\vec{\Sigma}\cdot\vec{E}\,]\,]\,]\\}\\}$
$\displaystyle\;-\frac{\mathrm{i}\,e}{48m^{6}}\,\\{\,\vec{\Sigma}\cdot\vec{\pi},\\{\,\vec{\Sigma}\cdot\vec{\pi},\\{\,\vec{\Sigma}\cdot\vec{\pi},\\{\,\vec{\Sigma}\cdot\vec{\pi},[\,\vec{\Sigma}\cdot\vec{\pi},\vec{\Sigma}\cdot\vec{E}\,]\,\\}\,\\}\,\\}\,\\}\,.$
(22)

For alternative representations, one notes the identities

$\displaystyle\mathcal{Q}=$ $\displaystyle\;-\frac{\beta
e^{2}}{96m^{5}}\,\;\\{\,\vec{\Sigma}\cdot\vec{E},4\,(\vec{\Sigma}\cdot\vec{\pi})^{2}\,\vec{\Sigma}\cdot\vec{E}+4\,\vec{\Sigma}\cdot\vec{E}(\vec{\Sigma}\cdot\vec{\pi})^{2}+\vec{\Sigma}\cdot\vec{\pi}\,\vec{\Sigma}\cdot\vec{E}\,\vec{\Sigma}\cdot\vec{\pi}\,\\}\,,$
(23a) $\displaystyle\mathcal{L}=$
$\displaystyle\;\frac{65\mathrm{i}e}{3072m^{6}}[(\vec{\Sigma}\cdot\vec{\pi})^{4},\\{\,\vec{\Sigma}\cdot\vec{\pi},\,\vec{\Sigma}\cdot\vec{E}\,\\}]-\frac{77\mathrm{i}e}{1536m^{6}}[(\vec{\Sigma}\cdot\vec{\pi})^{5},\vec{\Sigma}\cdot\vec{E}]$
$\displaystyle\;-\frac{43\,\mathrm{i}\,e}{1536m^{6}}\,[\,(\vec{\Sigma}\cdot\vec{\pi})^{3},\,(\vec{\Sigma}\cdot\vec{\pi})^{2}\,\vec{\Sigma}\cdot\vec{E}+\vec{\Sigma}\cdot\vec{E}\,(\vec{\Sigma}\cdot\vec{\pi})^{2}+\vec{\Sigma}\cdot\vec{\pi}\,\vec{\Sigma}\cdot\vec{E}\,\vec{\Sigma}\cdot\vec{\pi}\,]\,.$
(23b)

The particle-antiparticle symmetry implies that the terms are invariant when
the following transformations are simultaneously applied: (i) multiplication
by an overall factor $-1$, (ii) replacement $\beta\to-\beta$, (iii)
replacement $\vec{\pi}\to-\vec{\pi}$, and $\partial_{t}\to-\partial_{t}$ (iv)
replacement $\vec{\Sigma}\to-\vec{\Sigma}$, and (v) replacements $e\to-e$,
$V\to-V$ and $\vec{E}\to-\vec{E}$.

###  II.2 General Particle Hamiltonian

The upper left $2\times 2$ submatrix of $\mathcal{H}_{\rm FW}$ constitutes the
particle Hamiltonian, while the lower left $2\times 2$ submatrix of
$\mathcal{H}_{\rm FW}$ constitutes the antiparticle Hamiltonian. Here, we
concentrate on the particle Hamiltonian. Formally, the particle Hamiltonian
can be found from the results given in Eqs. (19)—(II.1) under the replacements
$\vec{\Sigma}\to\vec{\sigma}$, and $\beta\to\mathbbm{1}_{2\times 2}$, The
general Foldy–Wouthuysen transformed particle Hamiltonian $H$ under the
presence of the external electric and magnetic fields is obtained as

$H_{\rm FW}=\beta m+H^{[2]}+H^{[4]}+H^{[6]}+H^{[8]}\,.$ (24)

One finds

$\displaystyle H^{[2]}=$
$\displaystyle\;\frac{1}{2m}(\vec{\sigma}\cdot\vec{\pi})^{2}+V\,,$ (25)
$\displaystyle H^{[4]}=$
$\displaystyle\;-\,\frac{1}{8m^{3}}(\vec{\sigma}\cdot\vec{\pi})^{4}-\frac{\mathrm{i}\,e}{8m^{2}}[\,\vec{\sigma}\cdot\vec{\pi},\vec{\sigma}\cdot\vec{E}\,]\,,$
(26) $\displaystyle H^{[6]}=$
$\displaystyle\;\frac{(\vec{\sigma}\cdot\vec{\pi})^{6}}{16m^{5}}+\frac{5\mathrm{i}e}{128m^{4}}\left[\left(\vec{\sigma}\cdot\vec{\pi}\right)^{2},\left\\{\vec{\sigma}\cdot\vec{E},\vec{\sigma}\cdot\vec{\pi}\right\\}\right]$
$\displaystyle\;+\frac{3\mathrm{i}\,e}{64m^{4}}\;\\{\,(\vec{\sigma}\cdot\vec{\pi})^{2},[\,\vec{\sigma}\cdot\vec{\pi},\vec{\sigma}\cdot\vec{E}\,]\\}+\,\frac{e^{2}\vec{E}^{\,2}}{8m^{3}}\,.$
(27)

The eighth-order term comprise the kinetic term $K$, a term $D$ involving
temporal derivatives of the electric field, terms quadratic in the electric
field, denoted by $Q$, and linear terms in the electric field, which we denote
by $L$,

$\displaystyle H^{[8]}=$ $\displaystyle\;K+D+Q+L\,,\qquad
K=-\,\frac{5}{128m^{7}}(\vec{\sigma}\cdot\vec{\pi})^{8}\,,\qquad
D=-\frac{\mathrm{i}\,e^{2}}{32m^{4}}\,[\,\vec{\sigma}\cdot\vec{E},\vec{\sigma}\cdot\partial_{t}\vec{E}\,]+\frac{e}{48m^{5}}\,\\{(\vec{\sigma}\cdot\vec{\pi})^{3},\vec{\sigma}\cdot\partial_{t}\vec{E}\\}\,,$
$\displaystyle Q=$
$\displaystyle\;\frac{7\,e^{2}}{192m^{5}}\,\;[\,\vec{\sigma}\cdot\vec{\pi},\vec{\sigma}\cdot\vec{E}\,]\;[\,\vec{\sigma}\cdot\vec{\pi},\vec{\sigma}\cdot\vec{E}\,]-\frac{3\,e^{2}}{64m^{5}}\,\;\\{\,\vec{\sigma}\cdot\vec{\pi},\vec{\sigma}\cdot\vec{E}\,\\}\;\\{\,\vec{\sigma}\cdot\vec{\pi},\vec{\sigma}\cdot\vec{E}\,\\}-\frac{e^{2}}{24m^{5}}\,\;[\,\vec{\sigma}\cdot\vec{\pi},[\,\vec{\sigma}\cdot\vec{\pi},(\vec{\sigma}\cdot\vec{E})^{2}\,]\,]\,,$
$\displaystyle L=$
$\displaystyle\;-\frac{5\mathrm{i}e}{1024m^{6}}[\vec{\sigma}\cdot\vec{\pi},[\,\vec{\sigma}\cdot\vec{\pi},[\,\vec{\sigma}\cdot\vec{\pi},[\,\vec{\sigma}\cdot\vec{\pi},[\,\vec{\sigma}\cdot\vec{\pi},\sigma\cdot\vec{E}\,]\,]\,]\,]]-\frac{\mathrm{i}\,e}{32m^{6}}\,\\{\,\vec{\sigma}\cdot\vec{\pi},\\{\,\vec{\sigma}\cdot\vec{\pi},[\,\vec{\sigma}\cdot\vec{\pi},[\,\vec{\sigma}\cdot\vec{\pi},[\,\vec{\sigma}\cdot\vec{\pi},\sigma\cdot\vec{E}\,]\,]\,]\,\\}\,\\}$
$\displaystyle\;-\frac{\mathrm{i}\,e}{48m^{6}}\,\\{\,\vec{\sigma}\cdot\vec{\pi},\\{\,\vec{\sigma}\cdot\vec{\pi},\\{\,\vec{\sigma}\cdot\vec{\pi},\\{\,\vec{\sigma}\cdot\vec{\pi},[\,\vec{\sigma}\cdot\vec{\pi},\sigma\cdot\vec{E}\,]\,\\}\,\\}\,\\}\,\\}\,.$
(28)

For alternative representations, one notes identities analogous to Eqs. (23a)
and (23),

$\displaystyle Q=$
$\displaystyle\;-\frac{e^{2}}{96m^{5}}\,\;\\{\,\vec{\sigma}\cdot\vec{E},4\,(\vec{\sigma}\cdot\vec{\pi})^{2}\,\vec{\sigma}\cdot\vec{E}+4\,\vec{\sigma}\cdot\vec{E}(\vec{\sigma}\cdot\vec{\pi})^{2}+\vec{\sigma}\cdot\vec{\pi}\,\vec{\sigma}\cdot\vec{E}\,\vec{\sigma}\cdot\vec{\pi}\,\\}\,,$
(29a) $\displaystyle L=$
$\displaystyle\;\frac{65\mathrm{i}e}{3072m^{6}}[(\vec{\sigma}\cdot\vec{\pi})^{4},\\{\,\vec{\sigma}\cdot\vec{\pi},\,\vec{\sigma}\cdot\vec{E}\,\\}]-\frac{77\mathrm{i}e}{1536m^{6}}[(\vec{\sigma}\cdot\vec{\pi})^{5},\vec{\sigma}\cdot\vec{E}]$
$\displaystyle\;-\frac{43\,\mathrm{i}\,e}{1536m^{6}}\,[\,(\vec{\sigma}\cdot\vec{\pi})^{3},\,(\vec{\sigma}\cdot\vec{\pi})^{2}\,\vec{\sigma}\cdot\vec{E}+\vec{\sigma}\cdot\vec{E}\,(\vec{\sigma}\cdot\vec{\pi})^{2}+\vec{\sigma}\cdot\vec{\pi}\,\vec{\sigma}\cdot\vec{E}\,\vec{\sigma}\cdot\vec{\pi}\,]\,.$
(29b)

##  III Applications

###  III.1 Coulomb Field Coupling

One of the important applications of the Hamiltonian (24) concerns Coulombic
bound states, which are relevant to one-electron ions in the central field of
a nucleus of charge number $Z$. In this case, one has the relations

$\displaystyle e\,A^{0}=$
$\displaystyle\;V=-\frac{Z\alpha}{r}\,,\qquad\vec{A}=\vec{0}\,,\qquad\vec{B}=\vec{0}\,,$
(30) $\displaystyle\qquad e\,\vec{E}=$
$\displaystyle\;-e\,\vec{\nabla}A^{0}=-\vec{\nabla}V\,,\qquad\vec{\pi}=\vec{p}\,.$
(31)

The scalar potential and electric field are time-independent in this case. One
ends up with the following leading term, where the subscript $C$ indicates the
relevance for the Coulomb field,

$H^{[2]}_{C}=\frac{\vec{p}^{\,2}}{2m}+V=\frac{\vec{p}^{\,2}}{2m}-\frac{Z\alpha}{r}\,.$
(32)

This is the Schrödinger–Coulomb Hamiltonian in the nonrecoil approximation
(Chap. 4 of Ref. JeAd2022book ). For the evaluation of the fourth-order
corrections, one needs the identities

$\displaystyle\mathrm{i}e\,[\vec{\sigma}\cdot\vec{p},\vec{\sigma}\cdot\vec{E}]=$
$\displaystyle\;-\vec{\nabla}^{2}V-2\,\vec{\sigma}\cdot(\vec{\nabla}V\times\vec{p})\,,$
(33a) $\displaystyle
e\\{\vec{\sigma}\cdot\vec{p},\vec{\sigma}\cdot\vec{E}\\}=$
$\displaystyle\;-\\{\vec{\sigma}\cdot\vec{p},\vec{\sigma}\cdot\vec{\nabla}V\\}=-\mathrm{i}[\vec{p}^{\,2},V]\,.$
(33b)

The well-known leading relativistic corrections in the Foldy–Wouthuysen
Hamiltonian in this case reads as follows,

$\displaystyle H^{[4]}_{C}=$
$\displaystyle\;-\frac{\vec{p}^{\,4}}{8m^{3}}+\frac{1}{8m^{2}}\vec{\nabla}^{2}V+\frac{1}{4m^{2}}\vec{\sigma}\cdot(\vec{\nabla}V\times\vec{p})$
$\displaystyle=$
$\displaystyle\;-\frac{\vec{p}^{\,4}}{8m^{3}}+\frac{\pi\,(Z\alpha)}{2m^{2}}\;\delta^{(3)}(\vec{r})+\frac{Z\alpha}{4m^{2}\,r^{3}}\;\vec{\sigma}\cdot\vec{L}\,,$
(34)

where $\vec{L}$ is the orbital angular momentum operator. The sixth-order
corrections attain the following form,

$\displaystyle H^{[6]}_{C}=$
$\displaystyle\;\frac{\vec{p}^{\,6}}{16m^{5}}-\frac{3\\{\,\vec{p}^{\,2},\vec{\nabla}^{2}V\,\\}}{64m^{4}}\;-\frac{3\\{\,\vec{p}^{\,2},\vec{\sigma}\cdot(\vec{\nabla}V\times\vec{p})\\}}{32m^{4}}\;$
$\displaystyle\;+\frac{5}{128m^{4}}\,\left[\vec{p}^{\,2},[\vec{p}^{\,2},V]\right]+\,\frac{(\vec{\nabla}V)^{\,2}}{8m^{3}}\,.$
(35)

With $\partial_{t}\vec{E}=\vec{0}$, one has for the eighth-order corrections,

$\displaystyle H^{[8]}_{C}=$ $\displaystyle\;K_{C}+D_{C}+Q_{C}+L_{C}\,,\qquad
K_{C}=-\frac{5}{128m^{7}}\vec{p}^{\,8}\,,\qquad D_{C}=0\,,$ (36a)
$\displaystyle Q_{C}=$
$\displaystyle\;\frac{7[\vec{\nabla}^{2}V+2\,\vec{\sigma}\cdot(\vec{\nabla}V\times\vec{p})]\,[\vec{\nabla}^{2}V+2\,\vec{\sigma}\cdot(\vec{\nabla}V\times\vec{p})]}{192m^{5}}+\frac{3\,\;[\vec{p}^{\,2},V]\,[\vec{p}^{\,2},V]}{64m^{5}}-\frac{[\,\vec{\sigma}\cdot\vec{p},[\,\vec{\sigma}\cdot\vec{p},(\vec{\nabla}V)^{2}\,]\,]}{24m^{5}}\,,$
(36b) $\displaystyle L_{C}=$
$\displaystyle\;\frac{5[\vec{\sigma}\cdot\vec{p},[\,\vec{\sigma}\cdot\vec{p},[\,\vec{\sigma}\cdot\vec{p},[\,\vec{\sigma}\cdot\vec{p},\vec{\nabla}^{2}V+2\,\vec{\sigma}\cdot(\vec{\nabla}V\times\vec{p})]\,]\,]\,]}{1024m^{6}}+\frac{\\{\,\vec{\sigma}\cdot\vec{p},\\{\,\vec{\sigma}\cdot\vec{p},[\,\vec{\sigma}\cdot\vec{p},[\,\vec{\sigma}\cdot\vec{p},\vec{\nabla}^{2}V+2\,\vec{\sigma}\cdot(\vec{\nabla}V\times\vec{p})\,]\,]\,\\}\,\\}}{32m^{6}}\,$
$\displaystyle\;+\frac{1}{48m^{6}}\,\\{\vec{\sigma}\cdot\vec{p},\\{\vec{\sigma}\cdot\vec{p},\\{\vec{\sigma}\cdot\vec{p},\\{\vec{\sigma}\cdot\vec{p},\vec{\nabla}^{2}V+2\,\vec{\sigma}\cdot(\vec{\nabla}V\times\vec{p})\\}\\}\\}\\}\,.$
(36c)

###  III.2 Application: $\maybebm{F_{5/2}}$ States

We shall compare the eighth-order corrections to the bound-state energy
obtained from Eq. (36a) to the bound-state energies of the Dirac–Coulomb
problem. It is well known that the Dirac–Coulomb problem can be solved exactly
(Chap. 8 of Ref. JeAd2022book ), with the result

$\displaystyle E_{\mathrm{D}}=$
$\displaystyle\;m\,\left(1+\frac{(Z\alpha)^{2}}{(n_{r}+\gamma)^{2}}\right)^{-\tfrac{1}{2}}\,,$
$\displaystyle n_{r}=$
$\displaystyle\;n-j-\tfrac{1}{2}\,,\qquad\gamma=\sqrt{(j+1/2)^{2}+(Z\alpha)^{2}}\,.$
(37)

For $nF_{5/2}$ states, the presence of large spin-orbit coupling implies the
emergence of nontrivial corrections to the energy from the corresponding
higher-order terms in Eq. (36a). One expands as follows,

$E_{\mathrm{D}}(nF_{5/2})=m-\frac{(Z\alpha)^{2}m}{2n^{2}}+(Z\alpha)^{4}m\left(\frac{3}{8n^{4}}-\frac{1}{6n^{3}}\right)\\\
-(Z\alpha)^{6}m\left(\frac{1}{216n^{3}}+\frac{1}{24n^{4}}-\frac{1}{4n^{5}}+\frac{5}{16n^{6}}\right)\\\
+(Z\alpha)^{8}m\Bigl{(}-\frac{1}{3888n^{3}}-\frac{1}{432n^{4}}-\frac{1}{432n^{5}}\\\
+\frac{5}{48n^{6}}-\frac{5}{16n^{7}}+\frac{35}{128n^{8}}\Bigr{)}+\mathcal{O}(Z\alpha)^{10}\,.$
(38)

Within the Foldy–Wouthuysen method, the eighth-order terms combine several
effects, namely, the combined effect of third-order perturbative effect
$E^{[8]}=E^{[8]}_{c}$ generated by the fourth-order Hamiltonian $H^{[4]}$, the
mixed fourth- and sixth-order $E^{[8]}_{m}$, and diagonal element of eighth
order $E^{[8]}_{d}$,

$E^{[8]}=E^{[8]}_{c}+E^{[8]}_{m}+E^{[8]}_{d}\,.$ (39)

The terms will be further examined in the following. Let $G^{\prime}$ be the
reduced Green function,

$G^{\prime}=\left(\frac{1}{E_{S}-H_{S}}\right)^{\prime}\,,$ (40)

where $E_{S}$ is the Schrödinger–Coulomb energy
$E_{S}=-\frac{(Z\alpha)^{2}m}{2n^{2}}$ and $H_{S}=H^{[2]}_{C}$ is the
Schrödinger–Coulomb Hamiltonian. Then,

$\displaystyle E^{[8]}_{c}=$ $\displaystyle\;\langle
H^{[4]}\,G^{\prime}\,(H^{[4]}-\langle
H^{[4]}\rangle)\,G^{\prime}\,H^{[4]}\rangle\,,$ (41a) $\displaystyle
E^{[8]}_{m}=$ $\displaystyle\;2\langle H^{[4]}\,G^{\prime}\,H^{[6]}\rangle\,,$
(41b) $\displaystyle E^{[8]}_{d}=$ $\displaystyle\;\langle
H^{[8]}\rangle=\langle K_{C}\rangle+\langle Q_{C}\rangle+\langle
L_{C}\rangle\,.$ (41c)

After lengthy algebra, we obtain the results,

$\displaystyle\frac{\langle K_{C}\rangle_{nF_{5/2}}}{(Z\alpha)^{8}m}=$
$\displaystyle\;-\frac{2}{693n^{3}}+\frac{65}{1386n^{5}}-\frac{1549}{5544n^{7}}+\frac{35}{128n^{8}}\,,$
(42a) $\displaystyle\frac{\langle Q_{C}\rangle_{nF_{5/2}}}{(Z\alpha)^{8}m}=$
$\displaystyle\;-\frac{4}{31185n^{3}}+\frac{131}{99792n^{5}}-\frac{23}{13860n^{7}}\,,$
(42b) $\displaystyle\frac{\langle L_{C}\rangle_{nF_{5/2}}}{(Z\alpha)^{8}m}=$
$\displaystyle\;-\frac{281}{249480n^{3}}+\frac{5975}{399168n^{5}}-\frac{599}{13860n^{7}}\,,$
(42c) $\displaystyle\frac{E^{[8]}_{d}(nF_{5/2})}{(Z\alpha)^{8}m}=$
$\displaystyle\;-\frac{1033}{249480n^{3}}+\frac{25219}{399168n^{5}}-\frac{8989}{27720n^{7}}$
$\displaystyle\;+\frac{35}{128n^{8}}\,,$ (42d)

For absolutely clarity, we should emphasize the all matrix elements are
calculated with nonrelativistic Schrödinger–Pauli two-component reference-
state wave functions (Chap. 6 of Ref. JeAd2022book ). Summing up the results
for the combined third-order perturbation theory term, and the mixed term, one
obtains

$\displaystyle\frac{E^{[8]}_{c}(nF_{5/2})}{(Z\alpha)^{8}m}=$
$\displaystyle\;-\frac{169721}{23950080n^{3}}-\frac{19}{1728n^{4}}+\frac{4560727}{41912640n^{5}}$
$\displaystyle\;+\frac{955}{3024n^{6}}-\frac{1052}{693n^{7}}+\frac{21}{16n^{8}}\,,$
(43) $\displaystyle\frac{E^{[8]}_{m}(nF_{5/2})}{(Z\alpha)^{8}m}=$
$\displaystyle\;\frac{262729}{23950080n^{3}}+\frac{5}{576n^{4}}-\frac{3652871}{20956320n^{5}}$
$\displaystyle\;-\frac{40}{189n^{6}}+\frac{28271}{18480n^{7}}-\frac{21}{16n^{8}}\,,$
(44) $\displaystyle\frac{E^{[8]}(nF_{5/2})}{(Z\alpha)^{8}m}=$
$\displaystyle\;-\frac{1}{3888n^{3}}-\frac{1}{432n^{4}}-\frac{1}{432n^{5}}+\frac{5}{48n^{6}}$
$\displaystyle\;-\frac{5}{16n^{7}}+\frac{35}{128n^{8}}\,.$ (45)

The latter terms confirm Eq. (38).

###  III.3 Application: $\maybebm{F_{7/2}}$ States

The calculation proceeds in full analogy with $nF_{5/2}$ states. The
Dirac–Coulomb energy finds the expansion,

$E_{\mathrm{D}}(nF_{7/2})=m-\frac{(Z\alpha)^{2}m}{2n^{2}}+(Z\alpha)^{4}m\left(\frac{3}{8n^{4}}-\frac{1}{8n^{3}}\right)\\\
-(Z\alpha)^{6}m\left(-\frac{1}{512n^{3}}-\frac{3}{128n^{4}}+\frac{3}{16n^{5}}-\frac{5}{16n^{6}}\right)\\\
+(Z\alpha)^{8}m\Bigl{(}-\frac{1}{16384n^{3}}-\frac{3}{4096n^{4}}-\frac{1}{1024n^{5}}\\\
+\frac{15}{256n^{6}}-\frac{15}{64n^{7}}+\frac{35}{128n^{8}}+\mathcal{O}(Z\alpha)^{10}\,.$
(46)

After lengthy algebra, we obtain the results,

$\displaystyle\frac{\langle K_{C}\rangle_{nF_{7/2}}}{(Z\alpha)^{8}m}=$
$\displaystyle\;-\frac{2}{693n^{3}}+\frac{65}{1386n^{5}}-\frac{1549}{5544n^{7}}+\frac{35}{128n^{8}}\,,$
(47a) $\displaystyle\frac{\langle Q_{C}\rangle_{nF_{7/2}}}{(Z\alpha)^{8}m}=$
$\displaystyle\;-\frac{5}{18144n^{3}}+\frac{247}{90720n^{5}}-\frac{1}{315n^{7}}\,,$
(47b) $\displaystyle\frac{\langle L_{C}\rangle_{nF_{7/2}}}{(Z\alpha)^{8}m}=$
$\displaystyle\;\frac{13261}{7983360n^{3}}-\frac{156853}{7983360n^{5}}+\frac{9589}{221760n^{7}}\,,$
(47c) $\displaystyle\frac{E^{[8]}_{d}(nF_{7/2})}{(Z\alpha)^{8}m}=$
$\displaystyle\;-\frac{121}{80640n^{3}}+\frac{2417}{80640n^{5}}-\frac{965}{4032n^{7}}$
$\displaystyle\;+\frac{35}{128n^{8}}\,.$ (47d)

The combined third-order terms $E^{[8]}_{c}(nF_{7/2})$ and mixed terms
$E^{[8]}_{m}(nF_{7/2})$ find the representations,

$\displaystyle\;\frac{E^{[8]}_{c}(nF_{7/2})}{(Z\alpha)^{8}m}=-\frac{548047}{232243200n^{3}}-\frac{253}{61440n^{4}}$
(48)
$\displaystyle\;\qquad+\frac{2636603}{50803200n^{5}}+\frac{1427}{8064n^{6}}-\frac{22847}{20160n^{7}}+\frac{21}{16n^{8}}\,,$
$\displaystyle\;\frac{E^{[8]}_{m}(nF_{7/2})}{(Z\alpha)^{8}m}=\frac{55147}{14515200n^{3}}+\frac{13}{3840n^{4}}-\frac{8417851}{101606400n^{5}}$
$\displaystyle\;\qquad-\frac{1909}{16128n^{6}}+\frac{7649}{6720n^{7}}-\frac{21}{16n^{8}}\,.$
(49)

The sum

$\frac{E^{[8]}(nF_{7/2})}{(Z\alpha)^{8}m}=\frac{E^{[8]}_{c}(nF_{7/2})}{(Z\alpha)^{8}m}+\frac{E^{[8]}_{m}(nF_{7/2})}{(Z\alpha)^{8}m}\\\
+\frac{E^{[8]}_{d}(nF_{7/2})}{(Z\alpha)^{8}m}=-\frac{1}{16384n^{3}}-\frac{3}{4096n^{4}}-\frac{1}{1024n^{5}}\\\
+\frac{15}{256n^{6}}-\frac{15}{64n^{7}}+\frac{35}{128n^{8}}\,,$ (50)

reproduces the terms of order $(Z\alpha)^{8}$ from the Dirac–Coulomb bound-
state energy (46).

##  IV Conclusions

In this article, we have extended the treatment of the Foldy–Wouthuysen
transformation to eighth order, based on the scaling of the operators outlined
in Eq. (3). The results are obtained by a straightforward application of the
elimination of odd operators by repeated unitary transformations of the form
outlined in Eq. (13). The sixth-order terms [Eqs. (II.1) and (II.2)] are
obtained in full agreement with the literature (Refs. ZaPa2010 ; HiLePaSo2013
; PaYePa2016 ; HaZhKoKa2020 ; ZhMeQi2023 ). For the eighth-order terms, we
give results in Eqs. (II.1) and (II.2). We have applied our general results to
the relativistic bound Coulomb problem in Sec. III. An application to
$nF_{5/2}$ and $nF_{7/2}$ states, which present large spin-orbit couplings,
confirms, analytically, that the Dirac–Coulomb bound-state energy can be
obtained, within the Foldy–Wouthuysen formalism, as a sum of combined third-
order perturbative effects generated by the leading relativistic corrections
[Eq. (41a)], mixed fourth-order and sixth-order Hamiltonian terms [Eq. (41b)],
and diagonal elements of the eighth-order Hamiltonian [Eq. (41c)]. The latter
terms are obtained as diagonal elements of our eighth-order terms for the
Coulomb field, evaluated on Schrödinger–Pauli wave functions (Chap. 6 of Ref.
JeAd2022book ). As outlined in Ref. HiLePaSo2013 , the results are important
in a wider context, in view of the fact that the Foldy–Wouthuysen Hamiltonian
determines (part of) the matching coefficients in the Hamiltonian of
Nonrelativistic Quantum Electrodynamics (NRQED).

## Acknowledgements

The author acknowledges helpful conversations with Gregory S. Adkins. This
work has been supported by the National Science Foundation (Grant
PHY–2110294).

##  Appendix A Comparison with the Literature

We compare our results with those of Ref. ZhMeQi2023 . The result for the
eighth-order terms given in Eq. (8) of Ref. ZhMeQi2023  can be written as the
sum
$\mathcal{K}+\mathcal{D}^{\prime}+\mathcal{Q}^{\prime}+\mathcal{L}^{\prime}$,
where

$\displaystyle\mathcal{D}^{\prime}=$
$\displaystyle\;-\frac{\mathrm{i}\,e^{2}}{32m^{4}}\,[\,\vec{\sigma}\cdot\vec{E},\vec{\sigma}\cdot\partial_{t}\vec{E}\,]\,,$
(51a) $\displaystyle\mathcal{Q}^{\prime}=$
$\displaystyle\;-\frac{e^{2}}{32m^{5}}\,\;\\{\,\vec{\sigma}\cdot\vec{E},2\,(\vec{\sigma}\cdot\vec{\pi})^{2}\,\vec{\sigma}\cdot\vec{E}+2\,\vec{\sigma}\cdot\vec{E}(\vec{\sigma}\cdot\vec{\pi})^{2}$
$\displaystyle\;+\vec{\sigma}\cdot\vec{\pi}\,\vec{\sigma}\cdot\vec{E}\,\vec{\sigma}\cdot\vec{\pi}\,\\}\,,$
(51b) $\displaystyle\mathcal{L}^{\prime}=$
$\displaystyle\;\frac{11\mathrm{i}e}{1024m^{6}}[(\vec{\sigma}\cdot\vec{\pi})^{4},\\{\,\vec{\sigma}\cdot\vec{\pi},\,\vec{\sigma}\cdot\vec{E}\,\\}]$
(51c)
$\displaystyle\;-\frac{31\mathrm{i}e}{512m^{6}}[(\vec{\sigma}\cdot\vec{\pi})^{5},\vec{\sigma}\cdot\vec{E}]-\frac{9\,\mathrm{i}\,e}{512m^{6}}\,[\,(\vec{\sigma}\cdot\vec{\pi})^{3},\,$
$\displaystyle\;(\vec{\sigma}\cdot\vec{\pi})^{2}\,\vec{\sigma}\cdot\vec{E}+\vec{\sigma}\cdot\vec{E}\,(\vec{\sigma}\cdot\vec{\pi})^{2}+\vec{\sigma}\cdot\vec{\pi}\,\vec{\sigma}\cdot\vec{E}\,\vec{\sigma}\cdot\vec{\pi}\,]\,.$

Here, $\mathcal{D}^{\prime}$ is the term from Eq. (8) in Ref. ZhMeQi2023
which contains temporal derivatives, $\mathcal{Q}^{\prime}$ is the term from
Eq. (8) in Ref. ZhMeQi2023  which is quadratic in the electric fields, and
$\mathcal{L}^{\prime}$ is the term from Eq. (8) in Ref. ZhMeQi2023  which is
linear in the electric fields, The kinetic term $\mathcal{K}$ of eighth order
agrees with the corresponding term from Eq. (II.1) here.

Specifically, the term $\mathcal{D}^{\prime}$ lacks the terms proportional to
the anticommutator
$\\{(\vec{\sigma}\cdot\vec{\pi})^{3},\vec{\sigma}\cdot\partial_{t}\vec{E}\\}$
in comparison to our result for $\mathcal{D}$ given in Eq. (II.1). The results
for $\mathcal{Q}^{\prime}$ and $\mathcal{L}^{\prime}$ differ from those given
in Eqs. (23a) and (23) in the prefactors of the individual terms.

One can easily specialize the 4 operators given in Eqs. (51) and (51c) to the
case of a Coulomb field, on the basis of Eq. (30). This leads to the operators
$Q^{\prime}_{C}$ and $L^{\prime}_{C}$. The diagonal matrix elements, for
$nF_{5/2}$ evaluate as follows,

$\displaystyle\frac{\langle
Q^{\prime}_{C}\rangle_{nF_{5/2}}}{(Z\alpha)^{8}m}=$
$\displaystyle\;-\frac{37}{748440n^{3}}+\frac{779}{1496880n^{5}}-\frac{29}{41580n^{7}}$
(52a) $\displaystyle\frac{\langle
L^{\prime}_{C}\rangle_{nF_{5/2}}}{(Z\alpha)^{8}m}=$
$\displaystyle\;-\frac{23}{20790n^{3}}+\frac{1399}{95040n^{5}}-\frac{593}{13860n^{7}}\,.$
(52b) For $nF_{7/2}$ states, we have the results, $\displaystyle\frac{\langle
Q^{\prime}_{C}\rangle_{nF_{7/2}}}{(Z\alpha)^{8}m}=$
$\displaystyle\;-\frac{131}{997920n^{3}}+\frac{1301}{997920n^{5}}-\frac{16}{10395n^{7}}\,,$
(52c) $\displaystyle\frac{\langle
L^{\prime}_{C}\rangle_{nF_{7/2}}}{(Z\alpha)^{8}m}=$
$\displaystyle\;\frac{947}{532224n^{3}}-\frac{7921}{380160n^{5}}+\frac{1101}{24640n^{7}}\,.$
(52d)

These results differ, both individually as well as in the sum, from those
given in Eqs. (42) and (47), for the diagonal matrix elements of the operators
$Q_{C}$ and $L_{C}$, for both fine-structure components ($F_{5/2}$ and
$F_{7/2}$) investigated here.

## References

  * (1) L. L. Foldy and S. A. Wouthuysen, On the Dirac Theory of Spin $1/2$ Particles and Its Non-Relativistic Limit, Phys. Rev. 78, 29–36 (1950).
  * (2) U. D. Jentschura and G. S. Adkins, Quantum Electrodynamics: Atoms, Lasers and Gravity (World Scientific, Singapore, 2022).
  * (3) U. D. Jentschura and J. H. Noble, Nonrelativistic limit of the Dirac–Schwarzschild Hamiltonian: Gravitational Zitterbewegung and gravitational spin-orbit coupling, Phys. Rev. A 88, 022121 (2013).
  * (4) L. S. Brown and G. Gabrielse, Geonium theory: Physics of a single electron or ion in a Penning trap, Rev. Mod. Phys. 58, 233–311 (1986).
  * (5) A. Wienczek, C. Moore, and U. D. Jentschura, Foldy–Wouthuysen transformation in strong magnetic fields and relativistic corrections for quantum cyclotron energy levels, Phys. Rev. A 106, 012816 (2022).
  * (6) J. D. Bjorken and S. D. Drell, Relativistic Quantum Mechanics (McGraw-Hill, New York, 1964).
  * (7) J. Zatorski and K. Pachucki, Electrodynamics of finite-size particles with arbitrary spin, Phys. Rev. A 82, 052520 (2010).
  * (8) R. J. Hill, G. Lee, G. Paz, and M. P. Solon, NRQED Lagrangian at order $1/M^{4}$, Phys. Rev. D 87, 053017 (2013).
  * (9) V. Patkos, V. A. Yerokhin, and K. Pachucki, Higher-order recoil corrections for triplet states of the helium atom, Phys. Rev. A 94, 052508 (2016).
  * (10) M. Haidar, Z.-X. Zhong, V. I. Korobov, and J.-P. Karr, Nonrelativistic QED approach to the fine- and hyperfine-structure corrections of order $m\alpha^{6}$ and $m\alpha^{6}(m/M)$: Application to the hydrogen atom, Phys. Rev. A 101, 022501 (2020).
  * (11) W. Zhou, X. Mei, and H. Qiao, The $m\alpha^{8}$-order Foldy–Wouthuysen Hamiltonian and relativistic corrections to Coulomb systems, J. Phys. B 56, 045001 (2023).
  * (12) S. Wolfram, The Mathematica Book, 4 ed. (Cambridge University Press, Cambridge, UK, 1999).
  * (13) Y. N. Obukhov, A. J. Silenko, and O. V. Teryaev, General treatment of quantum and classical spinning particles in external fields, Phys. Rev. D 96, 105005 (2017).
  * (14) U. D. Jentschura and J. H. Noble, Foldy–Wouthuysen transformation, scalar potentials and gravity, J. Phys. A 47, 045402 (2014).
  * (15) K. Pachucki, Higher-order effective Hamiltonian for light atomic systems, Phys. Rev. A 71, 012503 (2005).


# On encounter rates in star clusters

Krisada Rawiraswattana  Division of Physical Science, Faculty of Science,
Prince of Songkla University, 15 Kanjanavanich Road,
Hatyai, Songkhla 90110, Thailand  Simon P. Goodwin Department of Physics and
Astronomy,
University of Sheffield, Sheffield S3 7RH, UK

###### Abstract

Close encounters between stars in star forming regions are important as they
can perturb or destroy protoplanetary discs, young planetary systems, and
stellar multiple systems. We simulate simple, viralised, equal-mass $N$-body
star clusters and find that both the rate and total number of encounters
between stars varies by factors of several in statistically identical clusters
due to the stochastic/chaotic details of orbits and stellar dynamics.
Encounters tend to rapidly ‘saturate’ in the core of a cluster, with stars
there each having many encounters, while more distant stars have none.
However, we find that the fraction of stars that have had at least one
encounter within a particular distance grows in the same way (scaling with
crossing time and half-mass radius) in all clusters, and we present a new
(empirical) way of estimating the fraction of stars that have had at least one
encounter at a particular distance.

methods: numerical — stars: kinematics and dynamics — open clusters and
associations: general

##  1 Introduction

Young stars are commonly found with circumstellar discs (e.g. Hillenbrand et
al., 1998; Lada et al., 2000; Haisch et al., 2000, 2001), and these discs are
thought to be where planet formation occurs. Since most stars are formed in
relatively dense environments (e.g. Lada & Lada, 2003), it is possible for the
discs, and the on going planet formation process within, to be affected by
close encounters between stars.

Simulations have shown that the effect of tidal perturbation from a stellar
fly-by can range from slightly changing the density distribution in the disc
to truncating or even destroying it (e.g. Clarke & Pringle, 1993; Cuello et
al., 2022), depending on how close the encounter is. This dynamical
truncation, as well as photoevaporation (e.g. Concha-Ramírez et al., 2022),
and face-on accretion (Wijnen et al., 2017), can significantly affect the
population of young stars with discs in the early stages of star formation.
Perturbations can also trigger disc instabilities (e.g. Thies et al., 2005,
2010) and may determine the population of planets forming in the disc (Ndugu
et al., 2022). Another interesting effect of encounters on the disc is the
misalignment between the rotational planes of the disc and the host star due
to a non-coplanar encounter (e.g. Heller, 1993; Larwood, 1997). Encounters may
alter already formed planetary systems, changing orbits (e.g. Breslau &
Pfalzner, 2019), or disrupting them (e.g. Parker & Quanz, 2012). And
encounters can similarly alter or destroy multiple stellar systems (e.g.
Goodwin, 2010; Reipurth et al., 2014).

Young stars with masses $\lesssim 1$ M☉ typically have discs with radii of a
few hundreds of au (e.g. Andrews & Williams, 2007). For the discs of those
stars to be significantly perturbed in encounters, the periastron distance
between the encountering stars needs to be less than $\sim 1000$ au.
Therefore, to understand how important encounters are in affecting
discs/planet formation, a key question is how many young stars have encounters
within $1000$ au?

There are two approaches one might take to finding the rates and numbers of
encounters in some star cluster of interest. The first is to perform $N$-body
simulations of a variety of similar systems which is time consuming and
computationally expensive (e.g. Parker & Quanz, 2012; Craig & Krumholz, 2013),
the second would be to have some (ideally analytic) estimate to quickly get at
least a ‘feel’ for the expected values.

In this paper, we examine the encounter rate in a number of $N$-body
simulations of bound star clusters. We show that encounter rates can vary by
up to an order of magnitude between statistically identical clusters, but the
fraction of stars that have had an encounter remains statistically the same.
We present an empirical way of estimating the fraction of stars in a cluster
that have had at least one encounter within a particular distance.

##  2 The encounter rate

The number of encounters per unit time ($\varepsilon$) for a star seems like
it should depend on several factors: the encounter distance of interest, some
average number density of stars, and the typical velocity of stars. The
velocity of stars will affect the encounter rate by both changing how often
stars encounter other stars, and also changing how effective gravitational
focusing is.

###  2.1 The standard method

The most common method of calculating encounter times is based on the
fundamental assumptions that a star is travelling through an effectively
infinite, uniform density medium at a constant speed (see e.g. the derivation
in Binney & Tremaine, 2008, note that these assumptions are perfectly adequate
if one is interested in e.g. the Galactic disc).

The encounter rate for any individual star is typically given by

$\varepsilon=4\sqrt{\pi}n\sigma\left(r_{\rm e}^{2}+\frac{Gm}{\sigma^{2}}r_{\rm
e}\right),$ (1)

where $n$ is the number density of the stars, $\sigma$ is the velocity
dispersion, $r_{\rm e}$ is the closest distance during the encounter, $G$ is
the gravitational constant, and $m$ is the typical mass of the stars (Binney &
Tremaine, 2008). The second term in the brackets is associated with the
gravitational focusing effect which deflects the trajectories and decreases
the distance of closest approach for slow encounters or encounters between
more massive stars.

For an ensemble of $N$ stars (e.g. a cluster), it seems reasonable to assume
that the total encounter rate (total number of encounters per unit time,
$\mathcal{E}$) scales with the total number of stars. Since each encounter
involves two stars, the encounter rate should scale with $N/2$, i.e.
$\mathcal{E}\simeq N\varepsilon/2$.

In convenient units where $r_{\rm e}$ is in au, $\sigma$ in km s-1, $n$ in
pc-3, and $m$ in M☉, the encounter rate in a cluster of $N$ stars is then

$\mathcal{E}=8.5\times 10^{-11}Nn\sigma\left(r_{\rm
e}^{2}+886\frac{m}{\sigma^{2}}r_{\rm e}\right)\ \textrm{Myr}^{-1}.$ (2)

Therefore, it would seem that to calculate the rate of encounters,
$\mathcal{E}$, at a particular distance of interest, $r_{\rm e}$, in an
ensemble of $N$ stars, the correct values of (a) number density, $n$, and (b)
velocity dispersion, $\sigma$, are required. If there is a distribution of
stellar masses, an appropriate value for $m$ must be taken.

While various assumptions that go into this simple calculation are clearly
wrong for star clusters (e.g. moving through an effectively infinite uniform
density medium at a constant speed) one might think that some simple variation
on this approach might work (e.g. taking some appropriate average speed and
density). However, we will show that this approach in star clusters gives an
often wrong, and an always misleading, ‘answer’.

###  2.2 What do we want to know?

It is important to clarify what we want to know about encounters in a cluster.
In most cases, what we would like is an estimate of what fraction of stars
have had a close encounter as this tells us the relative levels of
disc/planetary system/multiple system perturbation/destruction. It is
important to remember that this is not what an estimate of an encounter rate
gives without a further assumption of how the encounters are distributed
between stars.

As an example let us take a cluster that we shall examine in detail later: an
$N=300$, $M=300$ M☉, equal-mass (so $m=1$ M☉) virialised Plummer sphere
cluster with half-mass radius $r_{\rm h}=0.5$ pc. If we want to know the
number of stars that have had an encounter within, e.g. $r_{\rm e}=1000$ au,
we can calculate that $\mathcal{E}\sim 25$ Myr-1 by taking the values for the
velocity dispersion and half-mass density of this cluster and putting them
into equation (2).

If we assume this encounter rate estimate is correct, to calculate how many
stars have had an encounter after some time we need to make a further
assumption that encounters are random so that after $2$ Myr there will be 50
stars that have had an encounter, and after $10$ Myr $250$ stars (ie. $>80$
per cent) will have had an encounter. (One could be somewhat more
sophisticated and estimate as the encounter fraction starts to approach unity
how many stars have had zero, one, two etc. encounters.)

In assuming encounters are random, this calculation ignores that encounters
are much more likely to occur in the core, and that after a few crossing times
some stars in the core are likely to have had multiple encounters, while those
in the halo may have had none. Indeed, when stated like this, this approach
does seem extremely naive and it would be surprising if it gave the correct
answer.

##  3 Simulations

We investigate encounter rates in clusters by performing $N$-body simulations
in which we can record individual encounters, the stars involved, and their
distances of closest approach. The simulations we report here are of the
simplest bound systems: virialised Plummer spheres of equal-mass stars.

We simulate $N=300$ and $N=600$ virialised Plummer spheres (Plummer, 1911)
initialised by the method described in Aarseth et al. (1974) with initial
half-mass radii of $0.5$, $0.75$ and $1$ pc. Simulations are run only with
equal-mass stars so that we can ignore any complicating effects of mass
spectra.

The number of stars ($N$), the stellar mass ($m$), the initial half-mass
radius ($r_{\rm h}$), and the label of simulated clusters are shown in Table
1. Clusters with $r_{\rm h}=0.5$ pc are truncated at $3$ pc, those with
$r_{\rm h}=0.75$ pc are truncated at 5 pc, and $r_{\rm h}=1$ pc clusters are
truncated at 7 pc so that they have approximately the same relative sizes.

Table 1: The properties of cluster ensembles. Each ensemble contains 10 clusters with different random number seeds (labelled A to J) with the same number of stars ($N$), stellar masses ($m$), half-mass radii ($r_{\rm h}$), and crossing time ($t_{\rm cr}$). The ID of a simulation contains information on the initial conditions: N3 or N6 are $N=300$ and $N=600$ respectively, SM stands for single-mass, and after R is the initial half mass radius of the cluster. Cluster IDs | $N$ |  $m$ (M☉) |  $r_{\rm h}$ (pc) |  $t_{\rm cr}$ (Myr)
---|---|---|---|---
N3SMR050-A to J | $300$ | $1$ | $0.5$ | 1.0
N3SMR075-A to J | $0.75$ | 1.9
N3SMR100-A to J | $1$ | 3.1
N6SMR050-A to J | $600$ | $1$ | $0.5$ | 0.74

All simulations are run for $10$ Myr using our own $N$-body code. The code
uses a forth-order Hermite scheme (Makino & Aarseth, 1992) as the integrator.
We keep the energy error well below $10^{-4}$ by employing an adaptive
timestep, i.e. using equation (7) from Makino & Aarseth (1992) with parameter
$\eta=4\times 10^{-4}$ for $N=300$ clusters and $\eta=1\times 10^{-4}$ for
$N=600$ clusters. We also use block timestepping for $N=600$ runs to speed up
the calculations.

The separation between any pair of stars is monitored at every timestep. Once
two stars are closer to each other than $1000$ au, whether they are bound or
unbound, they are considered as having close encounter. During this period,
the closest separation is recorded and taken as the encounter distance once
the stars move away from each other beyond $1000$ au. In binaries multiple
‘encounters’ will occur, if the separation stays below $1000$ au this is only
counted as one ‘encounter’. This is to prevent hard binaries inflating the
close encounter rate, however as we discuss below hard binaries form extremely
rarely in our simulations.

Simulations are run with a gravitational softening length of $0.01$ au to
avoid collisions or computationally expensive very close encounters. This is
only of importance for the details of extremely close encounters at $\ll 10$
au which is much closer than the vast majority of encounters and at a distance
that would completely disrupt discs or planetary systems.

###  3.1 Number density and velocity dispersion

It would seem reasonable to assume that encounter rates and fractions should
depend in some way on both number density and velocity dispersion. Below we go
into some detail on how the number density and velocity dispersion of a
cluster might be quantified. This is important to show that the encounter
rates we measure in simulations disagree with simple calculations because the
assumptions that underlie them are wrong for clusters, rather than that we are
using the ‘wrong’ values for number density or velocity dispersion, or not
accounting for how they change with time. A reader happy to take our word for
this can skip the details below.

####  3.1.1 Number density

The derivation of equation (2) assumes that stars are uniformly distributed
and so $n$ is constant in time and space. This is a reasonable assumption for
e.g. encounters in the Galactic disc, but not for encounters in a cluster (or
any region where the number density varies on short length-scales).

There are a number of ways one could quantify some average number density, and
they can result in very different values for the estimated encounter rate.

The first average number density we use is simply calculated from the half-
mass radius of the cluster ($r_{\rm h}$):

$n_{\rm h}=\frac{3N}{4\pi{}r_{\rm h}^{3}}.$ (3)

The second average number density is the mean number density defined by

$n_{\rm
m}=\frac{\int_{0}^{\infty}nf\mathrm{d}{r}}{\int_{0}^{\infty}f\mathrm{d}{r}}=\int_{0}^{\infty}nf\mathrm{d}{r},$
(4)

where $n$ and $f$ are the number density and the probability density function
of the distance of stars from the centre of mass of the cluster ($r$). In
practice, we can approximate the integral by

$n_{\rm m}\simeq\Delta{}r\sum_{i=1}^{N_{\rm bin}}n_{i}f_{i},$ (5)

where $\Delta{}r$ and $N_{\rm bin}$ are the size and the number of bins of
stellar distances from the centre of mass of the cluster.

Theoretically, the probability density function is defined as
$f=\mathrm{d}{P}/\mathrm{d}{r}$, where $\mathrm{d}{P}$ is the probability of
finding stars at distance between $r$ and $r+\mathrm{d}{r}$ from the centre of
mass of the cluster. But practically, the value of $f_{i}$ in equation (5) may
be obtained from

$f_{i}=\left(\frac{\Delta{}P}{\Delta{}r}\right)_{i}=\frac{N_{i}}{N\Delta{}r},$
(6)

where $N_{i}$ is the number of stars in the $i^{\rm th}$-bin and $N$ is the
number of stars in the cluster. The number density $n_{i}$ in equation (5) is
related to $f_{i}$ via

$n_{i}=\left(\frac{\Delta{}N}{\Delta{}V}\right)_{i}=\frac{N}{4\pi{}r_{i}^{2}}\left(\frac{\Delta{}P}{\Delta{}r}\right)_{i}=\frac{Nf_{i}}{4\pi{}r_{i}^{2}},$
(7)

where $\Delta{}V$ is a spherical volume element containing $\Delta{}N$ stars.
Substituting equations (6) and (7) in (5) gives

$n_{\rm m}\simeq\frac{1}{4\pi{}N\Delta{}r}\sum_{i=1}^{N_{\rm
bin}}\frac{N_{i}^{2}}{r_{i}^{2}}.$ (8)

The half-mass number density is often used as it is simple to calculate, the
more complex mean number density has the advantage of including information on
the full density distribution of the cluster. We also note here that the half-
mass radius is often used as a characteristic radius as it remains roughly
constant over the long-term evolution of a cluster (Aarseth et al., 1974),
however the half-mass radius does fluctuate, especially at early times and so
even the half-mass density changes (sometimes by factors of several).

####  3.1.2 Velocity dispersion

The velocity dispersion ($\sigma$) in equation (1) and (2) comes from the
assumption that the velocity distribution of the stars in the cluster is
Maxwellian. From the Maxwell-Boltzmann distribution, the velocity dispersion
is related to the mode of the distribution ($v_{\rm m}$) by $\sigma=v_{\rm
m}/\sqrt{2}$. The mode can simply be determined by constructing the velocity
histogram and then fitting it with a polynomial regression to find the
position of the peak.

It should be noted that this way of finding the velocity dispersion is only
possible when all (3D) velocities are well known. In any observation the value
of $\sigma$ is either ‘guessed’ by assuming virial equilibrium, or observed in
either 1D (radial velocities) or 2D (proper motions) with usually quite
significant errors and biases (e.g. binary inflation, see Cottaar et al.
(2012)).

####  3.1.3 Time-averaging

There are two ways to calculate the number density and velocity dispersion to
use in equation (2). One is to take the values of $n$ or $\sigma$ calculated
instantaneously at the end of the simulation, the other is to take a time
average.

In a simulation it is possible to calculate full 3D time-averaged values for
various quantities. However, to estimate encounter rates in an observed region
or for a single snapshot only the current values for any quantity can be
calculated, and even they might be uncertain/guesstimated (e.g. no velocity
data is available and only 2D positions).

For later reference, Table 2 shows the initial, time-averaged and final (i.e.
those at 10 Myr) values of $\sigma$, $n_{\rm h}$, and $n_{\rm m}$ for all
simulations.

##  4 Results

In our simulations we follow all encounters at distances of $<1000$ au. We
record when the encounter occurred, which two stars were involved, and the
distance of closest approach. This allows us to find the encounter rate within
a particular distance (i.e. what equation (2) attempts to estimate), and the
number of stars that have had such an encounter – something equation (2)
cannot tell us without further assumptions, but is often what we wish to know.

###  4.1 Comparing an $N=300$ and an $N=600$ cluster

We start by comparing encounter rates at various distances in two clusters
with $N=300$ and $N=600$ equal-mass stars. Both are initially virialised
Plummer Spheres with $r_{\rm h}=0.5$ pc, with stars each of mass $1$ M☉.

Figure 1 shows the cumulative number of encounters over $10$ Myr in the
$N=300$ (top panel) and $N=600$ (bottom panel) clusters. In each panel the
lines from bottom-to-top are the cumulative numbers of encounters at $r_{\rm
e}<50$ (orange), $100$ (green), $500$ (magenta), and $1000$ (blue) au
respectively. At the top left of each sub-figure are three numbers for each of
the encounter distances: the first is the actual encounter rate (Myr-1) as
found in the simulation, the next two are time-averaged estimates that we will
return to later, but note for now that all three numbers are often quite
different. The two simulations shown are N3SMR050-A (top) and N6SMR050-A
(bottom).

Figure 1:  Cumulative numbers of encounters in runs N3SMR050-A (top panel, a)
and N6SMR050-A (bottom panel, b). The four curves in each panels are for
encounter distances $r_{\rm e}<50$ (orange), $100$ (green), $500$ (magenta),
and $1000$ (blue) au, from bottom to top. Numbers in the square brackets in
the top left are the encounter rates $[\mathcal{E}_{\rm sim}/\mathcal{E}_{\rm
est}(n_{\rm h})/\mathcal{E}_{\rm est}(n_{\rm m})]$ (see text).

Figure 1 shows a number of features one might expect.

  1. 1.

The total number of encounters grows roughly linearly with time (in these two
clusters at least).

  2. 2.

The number of encounters at different distances scales very roughly with
$r_{\rm e}^{2}$ (e.g. for $N=300$ after $10$ Myr there have been $317$
encounters at $r_{\rm e}<500$ au, and $27$ at $<50$ au).

  3. 3.

Increasing both $N$ (and therefore also $n$) by a factor of two results in
about $4$ times more encounters (e.g. $317$ when $N=300$, and $1165$ when
$N=600$ at $r_{\rm e}<500$ au).

The second and third numbers in the top left are the estimates of encounter
rate as calculated from equation (2) using the time-averaged values of the
half-mass number density ($n_{\rm h}$) and the mean number density ($n_{\rm
m}$) respectively.

In both cases using $n_{\rm h}$ under-estimates the number of encounters by a
factor of $\sim 2$. Using $n_{\rm m}$ seems better, often giving a reasonable
estimate (but sometimes being off by a factor of $\sim 2$). So at a first
glance at just these two simulations one might consider that using $n_{\rm m}$
often provides a reasonable estimate of the encounter rate in a cluster.

However, we have only compared two simulations which just happened to be those
labelled A in our ensembles. As we show below, when we look at the whole
ensemble the picture becomes much more complicated, and this emphasises the
importance of looking at ensembles of simulations when dealing with $N$-body
systems.

Table 2:  The first three triple columns are the initial (Ini.), time averaged (Avg.) and final (End) values of the velocity dispersion ($\sigma$), the half-mass number density ($n_{\rm h}$) and the mean number density ($n_{\rm m}$) of clusters N3SMR050-A..J (initial conditions in Table 1). In the last triple column are the analytic estimates of the encounter rates at $r_{\rm e}<1000$ au, using the average half-mass number density ($\mathcal{E}_{\rm est}(n_{\rm h})$), and the average mean number density ($\mathcal{E}_{\rm est}(n_{\rm m})$), compared with the actual values of the encounter rate measured in the simulations ($\mathcal{E}_{\rm sim}$). Note that the encounter rate can also be calculated from the initial or final values of $\sigma$, $n_{\rm h}$ and $n_{\rm m}$.  Cluster |  $\sigma$ (km s-1) |  $n_{\rm h}$ (pc-3) |  $n_{\rm m}$ (pc-3) |  $\mathcal{E}$ at $r_{\rm e}<1$ kau (Myr-1)
---|---|---|---|---
Ini. | Avg. | End | Ini. | Avg. | End | Ini. | Avg. | End | $\mathcal{E}_{\rm est}(n_{\rm h})$ | $\mathcal{E}_{\rm est}(n_{\rm m})$ | $\mathcal{E}_{\rm sim}$
N3SMR050-A | $0.57$ | $0.50\pm 0.04$ | $0.55$ | $541$ | $662\pm 144$ | $786$ | $1925$ | $1320\pm 475$ | $1650$ | $38\pm 8$ | $76\pm 28$ | $77$
N3SMR050-B | $0.50$ | $0.51\pm 0.05$ | $0.43$ | $543$ | $634\pm 112$ | $693$ | $1735$ | $1018\pm 482$ | $735$ | $36\pm 7$ | $58\pm 28$ | $58$
N3SMR050-C | $0.55$ | $0.49\pm 0.04$ | $0.43$ | $557$ | $382\pm 66$ | $247$ | $1535$ | $950\pm 780$ | $2402$ | $22\pm 4$ | $56\pm 46$ | $131$
N3SMR050-D | $0.57$ | $0.54\pm 0.04$ | $0.57$ | $571$ | $507\pm 92$ | $474$ | $1877$ | $1396\pm 948$ | $1313$ | $28\pm 5$ | $78\pm 53$ | $143$
N3SMR050-E | $0.57$ | $0.53\pm 0.06$ | $0.49$ | $577$ | $647\pm 120$ | $475$ | $964$ | $1232\pm 413$ | $1813$ | $36\pm 7$ | $69\pm 24$ | $60$
N3SMR050-F | $0.53$ | $0.53\pm 0.05$ | $0.42$ | $571$ | $682\pm 102$ | $532$ | $766$ | $669\pm 174$ | $484$ | $38\pm 6$ | $38\pm 10$ | $30$
N3SMR050-G | $0.49$ | $0.50\pm 0.05$ | $0.50$ | $577$ | $669\pm 130$ | $480$ | $1906$ | $1291\pm 537$ | $1158$ | $39\pm 8$ | $75\pm 31$ | $90$
N3SMR050-H | $0.58$ | $0.51\pm 0.04$ | $0.48$ | $555$ | $582\pm 125$ | $551$ | $1435$ | $1710\pm 723$ | $1100$ | $33\pm 7$ | $98\pm 42$ | $112$
N3SMR050-I | $0.57$ | $0.50\pm 0.04$ | $0.43$ | $549$ | $437\pm 110$ | $228$ | $1029$ | $956\pm 497$ | $531$ | $25\pm 6$ | $55\pm 29$ | $208$
N3SMR050-J | $0.50$ | $0.50\pm 0.05$ | $0.46$ | $577$ | $441\pm 103$ | $375$ | $2359$ | $1583\pm 917$ | $978$ | $26\pm 6$ | $92\pm 53$ | $198$
N3SMR050-A..J | $0.54\pm 0.03$ | $0.51\pm 0.05$ | $0.48\pm 0.05$ | $562\pm 14$ | $564\pm 110$ | $484\pm 175$ | $1553\pm 506$ | $1213\pm 595$ | $1216\pm 603$ | $32\pm 7$ | $69\pm 34$ | $111\pm 60$

###  4.2 An ensemble of statistically identical clusters

We now consider all ten clusters in our $N=300$ equal-mass stars, and a half-
mass radius of $r_{\rm h}=0.5$ pc ensemble. The only difference between these
clusters is the random number seed used to generate the initial positions and
velocities. Therefore one would expect that the encounter rates in each would
be similar - and ideally be close to an analytic estimate.

Figure 2: Top panel: the final encounter rates ($\mathcal{E}_{\rm sim}$)
against the encounter distance ($r_{\rm e}$) from each cluster in the $N=300$,
$r_{\rm h}=0.5$ pc ensemble. Each cluster has a different colour as shown in
the top left. Middle and bottom panels: the ratio of the analytic estimate to
the actual encounter rate using the half-mass density (middle panel), and mean
density (bottom panel).

The top panel of Fig. 2 shows the final encounter rates ($\mathcal{E}_{\rm
sim}$) for each of our identical clusters measured after $10$ Myr from $r_{\rm
e}=0$ to $1000$ au.111The exact values of encounter rates at separations of
less than a few au may be affected by our softening, but a close encounter did
happen, we just might not be able to trust the distance of closest approach
too precisely. The simulation shown in the top panel of Fig. 1 is A which is
the black line towards the middle-bottom of all the lines.

The most important thing to note about the top panel of Fig. 2 is the total
encounter rates after 10 Myr varies very significantly between clusters with a
difference of almost an order of magnitude. There appears to be no ‘typical’
clusters and some outliers – just a seemingly random spread in encounter rates
between clusters that are initially statistically identical.

Interestingly, all of the curves in the top panel of Fig. 2 follow a
distribution that goes roughly as $r_{\rm e}^{2}$ suggesting that the
distribution of encounter distances is what would be expected for unbound
encounters. However, the distribution can slightly deviate from $r_{\rm
e}^{2}$, often due to three-body encounters between a single star and a binary
(which has formed during the simulation). These encounters can cause an
increase in the encounter rate at $r_{\rm e}\sim 1000$ au, as can be most
obviously seen in cluster I (dark red line at the top of the figure).

The other panels of Fig. 2 show the difference between the analytically
estimated encounter rates and the actual encounter rates for each cluster. In
the middle panel the half-mass density ($n_{\rm h}$) is used for the estimate,
and in the bottom panel it is the mean number density ($n_{\rm m}$). A good
match to the analytic estimate is the black dashed line at a ratio of unity.

Using the half-mass density never gives a good estimate, and can be wrong by
an order of magnitude. Using the mean density is slightly better - three or
four clusters stay reasonably close to unity, but most clusters are always
wrong by factors of several.

One obvious explanation would be that different clusters have changed
significantly (some expanding and some contracting?).

In Table 2 we give the initial, time averaged and final values of $\sigma$,
$n_{\rm h}$, and $n_{\rm m}$ for all ten clusters, as well as the final
cumulative encounter rate. The averages and variance of each quantity over all
the clusters are given in the bottom line.

The velocity dispersions ($\sigma$) are very similar between all clusters, but
measures of density can change significantly over time with quite different
time averaged and final values, and between different measures (half mass or
mean).

However, these variations do not seem to correlate with the vastly different
encounter rates. The last three columns show the encounter rates estimated
with $n_{\rm h}$ and $n_{\rm m}$ and then the actual encounter rates from the
simulations. Only once (cluster F) does $n_{\rm h}$ get close to predicting
the actual encounter rate. Estimates using $n_{\rm m}$ are reasonable for 5/10
of the clusters, but very wrong for 5/10. We think this was seen by Craig &
Krumholz (2013) who note that their simulations had far more encounters than
one might expect, but substructure complicated their analysis.

What is particularly interesting is that there is no systematic change in
encounters rates with any measure of density or how they evolve. Cluster F has
the lowest final mean density (484 pc-1) and the lowest encounter rate (30
Myr-1), but cluster I has the second lowest final mean density (531 pc-1), but
the highest encounter rate (208 Myr-1). Clusters B and E have almost the same
encounter rate, but final mean densities that are different by a factor of
over two (735 and 1813 pc-1).

One might think that maybe there was some extreme deviation in density at some
point in time that the time averaged densities do not properly include,
however this is not the case. In Fig. 3 we show the numbers of encounters with
$r_{\rm e}<1000$ au (top panel), half-mass density (middle panel), and mean
density (lower panel) for clusters D and I (cf. Fig. 1).

Cluster D (black line) shows a relatively linear increase in encounter numbers
to end with $\sim 1400$ encounters within 10 Myr. Cluster I (red line) is
similar to cluster D until a period between 6–7 Myr when the encounter rate
increases significantly.

There is no reason to think that the increased encounter rate in cluster I is
due to density variations however. The middle panel shows that both cluster’s
half-mass densities are very similar, and both fairly constant and declining
slightly. The bottom panel shows more variation in the mean densities with
short-lived fluctuations of factors of a few, but both clusters show this
behaviour, and, if anything, cluster D has higher densities. There are
fluctuations in the mean density of cluster I around when the encounter rate
increases significantly, but there are others when it does not.

Examination of the data shows that the large numbers of encounters in cluster
I at 6–7 Myr is due to a few pairs of stars having multiple self-encounters in
weakly bound pairs (cf. Moeckel & Clarke, 2011).

In all the ensembles of statistically identical clusters we have run we find
no systematic relationship between any measure of cluster density and
encounter rates (apart from occasionally in just a few of the clusters, but
these could be chance given that many fluctuations do not correlate).

####  4.2.1 Binaries

All our simulations start with no binaries. An interesting question is how
many binaries can form, and how they might alter the evolution.

Soft binaries are extremely easy to form (see Moeckel & Clarke, 2011), and can
inflate the encounter rate. Any wide binary with periastron below 1000au and
apastron above 1000au will be included as multiple encounters - however such
binaries are extremely soft and short lived in our simulations (this was seen
in cluster I).

Hard binaries, however, are much more difficult to form. The soft binaries we
find are very weakly bound and can appear due to fluctuations in the global
potential. However, to form a hard, long-lived, binary system requires a
three-body encounter as the third body is needed to carry-away the excess
energy (see Goodman & Hut, 1993). A back-of-the-envelope calculation suggests
hard binary formation should be rare in our clusters, and an examination of
the simulations finds only a few hard binaries have managed to form (one every
few simulations, and never more than one in a simulation).

Figure 3: The evolution of encounter rates and density for cluster N3SMR050-D
(black lines) and cluster N3SMR050-I (red lines). Top panel (a): the number of
encounters with time for $r_{\rm e}<1000$ au. Middle panel (b): the half-mass
density $n_{\rm h}$. Bottom panel (c): the mean density $n_{\rm m}$.

###  4.3 The number of stars having had an encounter

We clearly see that the total number of encounters between stars at any
particular distance can be different by maybe an order of magnitude in
initially statistically identical clusters, in a way that cannot be explained
by density fluctuations.

The encounter rates we have shown in Fig. 2 and Table 2 are the number of
times two stars come closer together than a particular distance. However, this
measure does not include information on if a particular star, or a particular
pair of stars, have had multiple encounters.

Star clusters have a density distribution with a high density core and
increasingly lower density as one moves outwards, and a Plummer profile is a
reasonable approximation to young, bound clusters.

Within this density distribution stars can have a variety of orbits (which can
change after encounters). Some stars will spend a significant amount of their
time in the high density core, some will spend most of their time in the low
density halo, and various combinations in between (orbits can be radial or
circular etc. and can change over time).

This means that each individual star will have a unique encounter history.
Those that spend a lot of time in the core may have multiple encounters, while
those in the halo may have none. In addition (as seen above), some stars may
get into loosely bound multiples (cf. Moeckel & Clarke, 2011) and potentially
have numerous self-encounters which can inflate the encounter rate
significantly (see above).

Table 3:  For each cluster in the $N=300$, $r_{\rm h}=0.5$ pc ensemble (column 1) with $N$ stars (column 2), we show the half-mass radius $r_{\rm h}$ (column 3), crossing time $t_{\rm cr}$ (column 4), and the total number $N_{\rm s}$ (column 5) and fraction $f_{\rm s}$ (column 6) of stars that have had an encounter within $1000$ au in $10$ Myr. Half mass radii and crossing times are calculated explicitly for each cluster. The final row shows the means and variances of each quantity over the ensemble.  Cluster | $N$ | $r_{\rm h}$ | $t_{\rm cr}$ | $N_{\rm s}$ | $f_{\rm s}$
---|---|---|---|---|---
N3SMR050-A | $300$ | $0.510$ | $1.06$ | $162$ | $0.54$
N3SMR050-B | $300$ | $0.509$ | $1.06$ | $163$ | $0.54$
N3SMR050-C | $300$ | $0.505$ | $1.05$ | $143$ | $0.48$
N3SMR050-D | $300$ | $0.500$ | $1.03$ | $155$ | $0.52$
N3SMR050-E | $300$ | $0.499$ | $1.03$ | $160$ | $0.53$
N3SMR050-F | $300$ | $0.501$ | $1.03$ | $164$ | $0.55$
N3SMR050-G | $300$ | $0.499$ | $1.03$ | $170$ | $0.57$
N3SMR050-H | $300$ | $0.505$ | $1.05$ | $153$ | $0.51$
N3SMR050-I | $300$ | $0.507$ | $1.05$ | $144$ | $0.48$
N3SMR050-J | $300$ | $0.499$ | $1.03$ | $157$ | $0.52$
N3SMR050-A..J | $300$ | $0.503\pm 0.004$ | $1.04\pm 0.01$ | $157\pm 9$ | $0.52\pm 0.03$

Despite the large variation in encounter rates, When we measure the encounter
fraction in each of our ten clusters we find that this value is statistically
the same. In Table 3 we show the number ($N_{\rm s}$) and fraction ($f_{\rm
s}$) of stars in each of the ten statistically identical clusters from Fig. 2
and Table 2 that have had an encounter within 1000 au after 10 Myr. This
number is between 143 and 170 ($157\pm 9$) – statistically consistent with
being the same number, and a little over half the stars in the cluster at
$f_{\rm s}=0.58\pm 0.03$.

This tells us that in all clusters in this ensemble the same fraction of stars
are having very different numbers of encounters.

We can also examine other ensembles of statistically identical clusters and we
find that the encounter fraction in different clusters in the same ensemble is
statistically the same.

The mean and variances of encounter fractions for each ensemble are given in
Table 4, but to summarise for $r_{\rm e}<1000$ au: for $N=300$ clusters with
$r_{\rm h}=0.75$pc, the encounter fraction is $0.38\pm 0.02$; for $N=300$
clusters with $r_{\rm h}=1$pc, the encounter fraction is $0.28\pm 0.03$; and
for $N=600$ clusters with $r_{\rm h}=0.5$pc, the encounter fraction is
$0.60\pm 0.02$.

We do not present the data here in detail, but the same is true for different
encounter distances: ie. the encounter fraction is lower when the distance is
e.g. 500au, but the encounter fraction is statistically the same within each
ensemble. (It is difficult to say anything about extremely close encounters as
we are into small-$N$ statistics.)

That the encounter fraction is constant is extremely interesting, as the most
useful measure of encounters is often how many stars have had at least one
encounter closer than a particular distance over a particular timescale.

####  4.3.1 Encounter fractions in different clusters

As we saw above, within an ensemble of statistically identical clusters the
fraction of stars that have at least one encounter within 1000 au within 10
Myr is statistically the same, but it is different between different
ensembles.

Figure 4:  Encounter fractions for each ensemble against the absolute time $t$ (top panel), crossing time $t_{\rm cr}$ (middle panel), and the encounter crossing times $t_{\rm cr}^{*}$ (bottom panel, see text). The top panel contains the legend for the IDs of each ensemble. The line is the mean value, and the shaded region shows the variance for $N=600$ clusters with $r_{\rm h}=0.5$ pc (blue), $N=300$ clusters with $r_{\rm h}=0.5$ pc (red), $N=300$ clusters with $r_{\rm h}=0.75$ pc (purple), and $N=300$ clusters with $r_{\rm h}=1$ pc (green).  Table 4:  For each ensemble (column 1) with $N$ stars (column 2), we show the mean and variance of the half-mass radii $r_{\rm h}$ (column 3), crossing times $t_{\rm cr}$ (column 4), and the total number $N_{\rm s}$ (column 5) and fraction $f_{\rm s}$ (column 5) of stars that have had an encounter within $1000$ au in $10$ Myr. Half-mass radii and crossing times are calculated explicitly for each cluster.  Cluster | $N$ | $r_{\rm h}$ | $t_{\rm cr}$ | $N_{\rm s}$ | $f_{\rm s}$
---|---|---|---|---|---
N3SMR050-A..J | $300$ | $0.503\pm 0.004$ | $1.04\pm 0.01$ | $157\pm 9$ | $0.52\pm 0.03$
N3SMR075-A..J | $300$ | $0.751\pm 0.002$ | $1.90\pm 0.01$ | $113\pm 5$ | $0.38\pm 0.02$
N3SMR100-A..J | $300$ | $1.012\pm 0.021$ | $2.98\pm 0.09$ | $84\pm 8$ | $0.28\pm 0.03$
N6SMR050-A..J | $600$ | $0.503\pm 0.004$ | $0.74\pm 0.01$ | $361\pm 13$ | $0.60\pm 0.02$

The top panel of Fig. 4 shows how the encounter fraction, $f_{\rm s}$
increases with (absolute) time for all of our ensembles. The blue line and
shaded region which shows the variance at the top are for the $N=600$ clusters
with $r_{\rm h}=0.5$ pc. The red line and shaded region are for the $N=300$
clusters with $r_{\rm h}=0.5$ pc. The purple line and shaded region are for
the $N=300$ clusters with $r_{\rm h}=0.75$ pc. And at the bottom the green
line and shaded region are for the $N=300$ clusters with $r_{\rm h}=1$ pc.

As can be seen, in each case the encounter fraction within ensembles evolves
in the same general way - rising rapidly and then flattening - but different
ensembles seem to evolve at different rates.

That encounters occur at different rates in these different ensembles should
not be a surprise as each of the clusters have a different internal dynamical
timescale set by their crossing time. In the middle panel of Fig. 4 we show
the encounter fractions by crossing time, rather than by physical time, and
the differences between the different ensembles becomes less pronounced. That
the (green) $N=300$ clusters with $r_{\rm h}=1$ pc have had the fewest
encounters is clearly to a large extent because these clusters are dynamically
much younger.

However, it is clearly not just dynamical age that is important as the lines
are still somewhat different. The reason for this is that the cluster size
plays a role. In all of these simulations we are counting encounters within
$1000$ au which is a more significant fraction of the distances between stars
in an $r_{\rm h}=0.5$ pc cluster than in an $r_{\rm h}=1$ pc cluster. So, we
would expect the encounter timescale to also be sensitive to the relative
impact parameter $(r_{\rm h}/r_{\rm e})^{2}$.

In the bottom panel of Fig. 4, we plot encounter fraction against an
‘encounter crossing time’, $t_{\rm cr}^{*}$, defined as

$t_{\rm cr}^{*}=t_{\rm cr}\left(\frac{r_{\rm h}/{\rm pc}}{r_{\rm e}/1000\,{\rm
au}}\right)^{2}.$ (9)

Now in the bottom panel we appear to have found a timescale on which all
clusters show extremely similar behaviour. For encounters within $1000$ au
there is a sharp rise in $f_{\rm s}$ in the first $5$ $t_{\rm cr}^{*}$ to a
point where roughly a third of all stars have had an encounter. It then takes
another $\sim 50$ $t_{\rm cr}^{*}$ for the next third of stars to have an
encounter.

To test the scaling with $(r_{\rm h}/r_{\rm e})^{2}$, in Fig. 5 we compare
encounters within $500$ au in an $r_{\rm h}=0.5$ pc cluster (red) with
encounters within $1000$ au in an $r_{\rm h}=1$ pc cluster (green). Here
$(r_{\rm h}/r_{\rm e})^{2}$ is the same in both clusters (half the encounter
distance, but half the half-mass radius), therefore we would expect the growth
of the different encounter fractions with just the crossing time to be the
same, which they are as is clear in the figure.

Note that there are various subtleties at play such as different velocity
dispersions causing the effect of gravitational focusing to be different, and
we are dealing with small-$N$ stochastic systems. But overall, the overlap
between evolution between different clusters within ensembles and very
different ensembles is impressive when scaled by crossing time and relative
encounter cross section.

Figure 5:  The encounter fractions against crossing time with $r_{\rm e}<1000$
au for $N=300$, $r_{\rm h}=1$ pc (green shaded region), and with $r_{\rm
e}<500$ au for $N=300$, $r_{\rm h}=0.5$ pc (red shaded region).

####  4.3.2 An empirical relationship

The bottom panel of Fig. 4 provides a way of getting a rough estimate of the
encounter fraction of stars, $f_{\rm s}$, at a particular encounter distance,
$r_{\rm e}$, after some time, $t$, in any cluster if one knows the crossing
time, $t_{\rm cr}$ and half-mass radius, $r_{\rm h}$. Then from equation (9)
one can calculate $t_{\rm cr}^{*}$, and so $t/t_{\rm cr}^{*}$.

It is worth pointing-out that the curve in the bottom panel looks like it
should have a fairly simple function that would provide a fit. However, we
have struggled to find a simple ($2$ or $3$ parameter) function that fits (the
problem is that the initial rise is much steeper than e.g. an exponential will
fit). Therefore we would suggest simply reading-off the value of $f_{\rm s}$
from the bottom panel of Fig. 4 for whatever value of $t/t_{\rm cr}^{*}$.

We stress that this is a rough estimate, however, as rough-and-ready as this
may be, it will still almost certainly provide a much better feeling for how
many stars have had an encounter than any attempt to use equation (2), and
then to extrapolate to an encounter fraction.

####  4.3.3 An example

We can take a roughly Orion Nebula-like cluster with $M=1000$ M☉, $N=2500$, a
half-mass radius $r_{\rm h}=0.7$ pc, and age $3$ Myr and attempt to estimate
what fraction of stars have had an encounter at $<1000$ au. For such a
cluster, $n=1800$ pc-3, and assuming virial equilibrium $\sigma=2.6$ km s-1,
and so $t_{\rm cr}=0.27$ Myr.

From equation (9) if $r_{\rm e}=1000$ au, then $t_{\rm cr}^{*}=0.5t_{\rm
cr}=0.14$ Myr. Therefore this cluster has an ‘encounter age’ of $t/t_{\rm
cr}^{*}\sim 20$. From the bottom panel of Fig. 4 this would suggest around
$40$ per cent of stars will have had an encounter within $1000$ au.

If we use equation (2), we find $\mathcal{E}\sim 1000$ Myr-1. For an age of
$3$ Myr this suggests $3000$ encounters among the $N=2500$ stars. An extremely
naive extrapolation might suggest that therefore all stars had an encounter
within $1000$ au. One can be a little more sophisticated and assume that
encounters are random which finds that $\sim 90$ per cent of stars will have
had an encounter (as each one of the $3000$ encounters involves two stars,
even if they are random most stars will have been involved in one). Even if
the value of $3000$ encounters in $3$ Myr happened, by luck, to be right, the
extension of this encounter number to the number of stars involved in
encounters is certainly not random.

##  5 Conclusions

We have performed $N$-body simulations of small star clusters to investigate
stellar encounters with separations $r_{\rm e}<1000$ au. This is the regime in
which discs, planetary systems, and multiple stellar systems can be
significantly perturbed or destroyed.

We find that the encounter rates vary by up to an order of magnitude between
statistically identical clusters. However, we find that the fraction of stars
that have had an encounter is statistically the same within statistically
identical clusters.

The fraction of stars that have had an encounter increases rapidly at early
dynamical times before flattening significantly once stars in orbits
particularly susceptible to encounters have had at least one encounter. This
depends on both the dynamical timescale of the cluster ($t_{\rm cr}$), and the
relative impact parameter $(r_{\rm h}/r_{\rm e})^{2}$.

We find a consistent, and reasonably tight, relationship between the fraction
of stars that have had an encounter and a modified crossing time $t_{\rm
cr}^{*}\propto t_{\rm cr}(r_{\rm h}/r_{\rm e})^{2}$.

The relationship we have found has a seemingly solid physical basis, but no
detailed theoretical underpinning (we are working on this). However, it
provides a simple way of extracting an estimate of the encounter fraction for
a particular cluster of a particular age from a figure. While this is
empirical, it almost certainly provides a much better estimate of the true
encounter fraction than any attempt to apply standard theory.

SG was partly funded by STFC consolidated grant ST/V000853/1.

## References

  * Aarseth et al. (1974) Aarseth, S. J., Henon, M., & Wielen, R. 1974, A&A, 37, 183
  * Andrews & Williams (2007) Andrews, S. M., & Williams, J. P. 2007, ApJ, 659, 705, doi: 10.1086/511741
  * Binney & Tremaine (2008) Binney, J., & Tremaine, S. 2008, Galactic Dynamics: Second Edition (Princeton University Press)
  * Breslau & Pfalzner (2019) Breslau, A., & Pfalzner, S. 2019, A&A, 621, A101, doi: 10.1051/0004-6361/201833729
  * Clarke & Pringle (1993) Clarke, C. J., & Pringle, J. E. 1993, MNRAS, 261, 190, doi: 10.1093/mnras/261.1.190
  * Concha-Ramírez et al. (2022) Concha-Ramírez, F., Wilhelm, M. J. C., & Zwart, S. P. 2022, MNRAS, doi: 10.1093/mnras/stac1733
  * Cottaar et al. (2012) Cottaar, M., Meyer, M. R., & Parker, R. J. 2012, A&A, 547, A35, doi: 10.1051/0004-6361/201219673
  * Craig & Krumholz (2013) Craig, J., & Krumholz, M. R. 2013, ApJ, 769, 150, doi: 10.1088/0004-637X/769/2/150
  * Cuello et al. (2022) Cuello, N., Ménard, F., & Price, D. J. 2022, arXiv e-prints, arXiv:2207.09752.  https://arxiv.org/abs/2207.09752
  * Goodman & Hut (1993) Goodman, J., & Hut, P. 1993, ApJ, 403, 271, doi: 10.1086/172200
  * Goodwin (2010) Goodwin, S. P. 2010, Philosophical Transactions of the Royal Society of London Series A, 368, 851, doi: 10.1098/rsta.2009.0254
  * Haisch et al. (2000) Haisch, Jr., K. E., Lada, E. A., & Lada, C. J. 2000, AJ, 120, 1396, doi: 10.1086/301521
  * Haisch et al. (2001) —. 2001, ApJ, 553, L153, doi: 10.1086/320685
  * Heller (1993) Heller, C. H. 1993, ApJ, 408, 337, doi: 10.1086/172591
  * Hillenbrand et al. (1998) Hillenbrand, L. A., Strom, S. E., Calvet, N., et al. 1998, AJ, 116, 1816, doi: 10.1086/300536
  * Lada & Lada (2003) Lada, C. J., & Lada, E. A. 2003, ARA&A, 41, 57, doi: 10.1146/annurev.astro.41.011802.094844
  * Lada et al. (2000) Lada, C. J., Muench, A. A., Haisch, Jr., K. E., et al. 2000, AJ, 120, 3162, doi: 10.1086/316848
  * Larwood (1997) Larwood, J. D. 1997, MNRAS, 290, 490, doi: 10.1093/mnras/290.3.490
  * Makino & Aarseth (1992) Makino, J., & Aarseth, S. J. 1992, PASJ, 44, 141
  * Moeckel & Clarke (2011) Moeckel, N., & Clarke, C. J. 2011, MNRAS, 415, 1179, doi: 10.1111/j.1365-2966.2011.18731.x
  * Ndugu et al. (2022) Ndugu, N., Abedigamba, O. P., & Andama, G. 2022, MNRAS, 512, 861, doi: 10.1093/mnras/stac569
  * Parker & Quanz (2012) Parker, R. J., & Quanz, S. P. 2012, MNRAS, 419, 2448, doi: 10.1111/j.1365-2966.2011.19911.x
  * Plummer (1911) Plummer, H. C. 1911, MNRAS, 71, 460, doi: 10.1093/mnras/71.5.460
  * Reipurth et al. (2014) Reipurth, B., Clarke, C. J., Boss, A. P., et al. 2014, in Protostars and Planets VI, ed. H. Beuther, R. S. Klessen, C. P. Dullemond, & T. Henning, 267, doi: 10.2458/azu_uapress_9780816531240-ch012
  * Thies et al. (2010) Thies, I., Kroupa, P., Goodwin, S. P., Stamatellos, D., & Whitworth, A. P. 2010, ApJ, 717, 577, doi: 10.1088/0004-637X/717/1/577
  * Thies et al. (2005) Thies, I., Kroupa, P., & Theis, C. 2005, MNRAS, 364, 961, doi: 10.1111/j.1365-2966.2005.09644.x
  * Wijnen et al. (2017) Wijnen, T. P. G., Pols, O. R., Pelupessy, F. I., & Portegies Zwart, S. 2017, A&A, 604, A91, doi: 10.1051/0004-6361/201731072