# Numerical study of the Transverse Diffusion coefficient for a one component
model of plasma

Lorenzo Valvo  Department of Mathematics, Università degli Studi di Milano,
Via Saldini 50, 20133 Milano, Italy. Email<EMAIL_ADDRESS>Andrea Carati
Department of Mathematics, Università degli Studi di Milano, Via Saldini 50,
20133 Milano, Italy. E-mail<EMAIL_ADDRESS>
###### Abstract

In this paper we discuss the results of some Molecular Dynamics simulations of
a magnetized One Component Plasma, targeted to estimate the diffusion
coefficient $D_{\perp}$ in the plane orthogonal to the magnetic field lines.
We find that there exists a threshold with respect to the magnetic field
strength $|\vec{B}|$: for weak magnetic field the diffusion coefficients
scales as $1/|\vec{B}|^{2}$, while a slower decay appears at high field
strength. The relation of this transition with the different mixing properties
of the microscopic dynamics is investigated by looking at the behavior of the
velocity auto correlation.

The diffusion process is well understood for stochastic motions (see Ref.
[3]), that are supposed to mimic the behaviour of a chaotic dynamical system.
Many questions are instead left open in the study of the diffusive properties
of a system which is in a partially ordered state (see for example Ref. [27]).
A central issue, as regards magnetized plasma confinement, is the diffusion of
charged particles in the direction perpendicular to the magnetic field lines.
A widely accepted law, predicting that the transversal diffusion $D_{\perp}$
coefficient is proportional to the inverse of the square of the magnetic field
strength $|\vec{B}|$, was proposed more than 50 years ago( see Ref. [21]).
Being based on kinetic theory, this law is expected to hold whenever the
microscopic dynamics is chaotic. However, as the magnetic field $|\vec{B}|$ is
increased, a partially ordered state seems to set in (see Refs. [7, 9]), at
least for a pure electron plasma. Our purpose was to investigate the
consequences (if any) of this transition on the diffusion process. So we have
performed Molecular Dynamics simulations of a magnetized one component plasma,
that is a set of mutually interacting electrons subject to a constant external
magnetic field. We estimate the diffusion coefficient $D_{\perp}$ in the plane
orthogonal to the field, for different values of the magnetic field strength
$|\vec{B}|$. We find that the kinetic law holds for low $|\vec{B}|$, when the
microscopic dynamics is chaotic. But as the magnetic field grows the diffusion
coefficient seems to saturate to a plateau, while the microscopic state turns
to a partially ordered one.

##  1  Introduction

In the years Sixties there was a great exchange among research groups of
plasma physics and of dynamical systems, as both were interested in the study
of the $1\tfrac{1}{2}$ Hamiltonian system that can represent the magnetic
field lines; see Ref. [12] for a broad historical overview. In more recent
years, another point of connection have become the study of diffusion process
in the phase space, see Ref. [3]. It has been shown that if the dynamics is
not fully chaotic, then the process of diffusion in the phase space can be
“not normal”, i.e., the mean square displacement doesn’t necessarily grow
linearly in time: the process is called super diffusive if the growth is
faster than linear, or sub diffusive otherwise; for a review see [27]. In both
cases there is no widespread agreement of the correct definition of the
diffusion coefficient to be adopted.

Figure 1:  Diffusion coefficient $D_{\perp}$ perpendicular to the magnetic
field versus $\beta$, computed by MD simulations. Circles are the numerical
results, while the solid line is the plot of the function
$D_{\perp}=\frac{\gamma A}{\gamma^{2}+\beta^{2}}$, with $\gamma=0.168$ and
$A=9.9$. Such values are obtained from the auto correlation of the electron
velocities, as explained in the text. The broken line is the plot of the
function $D_{\perp}\simeq\beta^{-0.7}$, corresponding to the law found in
paper [13].

Now, as it will be shown below, the equations of motion for a single electron
in a one component plasma subject to an external constant magnetic field
$\vec{B}$ can be recasted in a dimensionless form, in which the only parameter
appearing is the ratio $\beta=\sqrt{4\pi}\omega_{c}/\omega_{p}$ among the
cyclotron frequency $\omega_{c}=eB/mc$ and the plasma frequency
$\omega_{p}=\sqrt{4\pi ne^{2}/m}$ (in the c.g.s. system); here $n$ is the
particle density and $e$ the electronic charge. The parameter $\beta$ measures
the relative strength the magnetic Lorentz force acting on a single electron,
with respect to the electrostatic force due to all other electrons. In the
limit of $\beta\to+\infty$ the equations of motions decouple and the system
reduces (formally) to a set of independent electrons in a constant magnetic
field, i.e., to an integrable system.

So, it can be conjectured that for high magnetic field strength, the dynamics
will not be fully chaotic, (see papers [7] and [9]), and that the diffusion
process in phase space may be “anomalous”. Actually, it is impossible for us
to study numerically such process, and we limit ourselves to study the
diffusion of the electrons in the physical space. To this end we study the
diffusion coefficient, defined as usual (see for example Ref. [1]) by

$D_{\perp}=\lim_{t\rightarrow\infty}\frac{\langle\left|\Delta\vec{x}_{\bot}\right|^{2}\rangle}{4t}$
(1)

being $\left|\Delta\vec{x}_{\bot}\right|^{2}$ the mean displacement, in the
plane perpendicular to the magnetic field, of the electrons from their initial
positions, while $\langle\cdot\rangle$ should be the phase average. Instead in
this paper, following a common attitude, averages will be always taken as time
averages along the orbits. We have not investigated the relations between the
two averages.

As regards physical applications, small values of $D_{\perp}$ are important
for the purpose of plasma confinement in fusion devices. This is another
reason to investigate in which regime the diffusion coefficient is small.

Another dimensionless parameter which characterize the state of a plasma is
the coupling parameter $\Gamma=e^{2}/(ak_{B}T)$, where $a$ is the inter
particle spacing (related to the particle density $n$ by the relation
$a=n^{1/3}$), $T$ the plasma temperature and $k_{B}$ as usual the Boltzmann
constant. So defined, $\Gamma$ is the ratio among the mean coulomb energy of a
couple of nearest particles, and the mean kinetic energy. The weak coupling
regime is then defined by $\Gamma<1$ and the strong coupling regime by
$\Gamma\geq 1$. Up to now, because of the reasons explained in the following,
Molecular Dynamics (MD) simulations have been performed mainly in the strongly
coupled case, while the weakly coupled regime have been addressed mostly by
kinetic theory.

In literature it is possible to find different estimates for the diffusion
coefficient. The oldest one (see Ref. [21]) predicting the scaling law
$D_{\perp}\propto\beta^{-2}$, is obtained in the frame of the kinetic theory,
in the weak coupling regime. However, other different theories have been
proposed in the years, each giving a different law for the dependence of the
diffusion coefficient $D_{\perp}$ on the parameter $\beta$. A few of them are
summarized in Table 1 on page 135003-2 of Ref. [22]; another one is
percolation theory (see Refs. [17, 18, 13]), which predicts a scaling like
$\beta^{-0.7}$. This law is in the closest agreement with our results. It was
brought to our attention by an anonymous referee that we warmly thank.

But none of the proposed theories is based on the loss of chaoticity in the
Newtonian microscopic dynamics. Also the numerical works found in literature
fail to address this point. In fact, up to now, the behaviour of the diffusion
coefficient have been investigated by MD mostly in the unmagnetized case, see
for instance Refs. [10, 24, 4]. The magnetized case was studied in Ref. [22],
but only in the strong coupling regime: at the smallest value $\Gamma=1$ a
transition at $\beta\approx 1$ in the behavior of $D_{\perp}$ was observed,
but the origin of such a transition was not discussed. A similar transition
was observed also in two more recent works [2] and [25], where the diffusion
coefficient was studied numerically for $\Gamma$ down to $0.1$. However, those
works were based on a so-called Particle-Particle Particle-Mesh ($P^{3}M$)
code, which is a sort of hybrid between a kinetic and a MD code. We think that
such method is not suited to study the chaoticity of the microscopic dynamics.
More on the connections between our results and the cited paper will be said
in the conclusions.

So in this paper we perform a fully MD simulation for different values of
$\beta$ for a system of $N=4096$ electrons, for a fixed value of $\Gamma=0.1$,
which is the smallest value we were able to manage. The aim is to verify for
what value of $\beta$ the transition in the behavior of $D_{\perp}$ occurs,
and to observe the chaoticity of the dynamics, by computing the time auto
correlation of the transverse particle velocity.

In Section 2 we describe the model and the numeric algorithm used, in Section
3 the numerical results are reported while in Section 4 the conclusions
follow.

##  2 The model and the numerical scheme

The system we are considering is called in the literature a one component
plasma and it consists of a number $N$ of electrons in a cubic box of side $L$
with periodic boundary conditions, the electrons being subject to mutual
Coulomb interactions, and to an external constant magnetic field
$\vec{B}=B\vec{e}_{z}$ ($\vec{e}_{z}$ is the unit vector directed along the
$z$ axis). The density is then defined by $n=N/L^{3}$. This is considered a
model of a plasma as the positive ions are supposed to constitute a uniform
neutralizing background.

If $t$ denotes time and $\vec{x}_{i}$ the position of the $i$-th electron
(with $i=1,\dots,N$), we define dimensionless variables

$\vec{y}_{i}=a^{-1}\vec{x}_{i},\quad\tau=\omega_{c}t\ ,$ (2)

by rescaling distances with the inter particle spacing $a=n^{-1/3}$ and time
with the cyclotron frequency $\omega_{c}$. Using such variables, the equations
of motion for the $i$-th electron read

$\frac{d^{2}{\vec{y}}_{i}}{d\tau^{2}}=\vec{e}_{z}\times\frac{{d\vec{y}}_{i}}{d\tau}+\frac{1}{\beta^{2}}\vec{E}(\vec{y}_{i})$
(3)

where $\vec{E}(\vec{y}_{i})$ is the electric field acting on $i$-th electron
due to all other charges. The electric field of a periodic system of charges
can be computed via the Ewald formula (see for example Ref. [14]),

$\begin{split}&\vec{E}(\vec{y}_{i})=\\\
&\sum_{\vec{l}}\sum_{j=1}^{N}\frac{\vec{y}_{ij\vec{l}}}{|\vec{y}_{ij\vec{l}}|^{3}}\left[\operatorname{erfc}({\alpha}|\vec{y}_{ij\vec{l}}|)+\frac{2\alpha|\vec{y}_{ij\vec{l}}|}{\sqrt{\pi}}\exp(-{\alpha}^{2}{|\vec{y}_{ij\vec{l}}|}^{2})\right]\\\
&+\frac{4\pi}{N}\sum_{\vec{q}\neq
0}\sum_{j=1}^{N}\frac{\vec{q}}{q^{2}}e^{-q^{2}\\!/4{\alpha}^{2}}\sin(\vec{q}\cdot\vec{y}_{ji}),\qquad\alpha=\sqrt{\pi}N^{-1/6}\\\
\end{split}$ (4)

where $\vec{l}$ is a triplet of integers denoting the position of an image
cell, while $\vec{q}$ is a vector in the reciprocal lattice, i.e. is defined
by $\vec{q}=2\pi\vec{n}/L$ with $\vec{n}$ an integer vector, and finally we
have defined $\vec{y}_{ij\vec{l}}=\vec{y}_{i}-\vec{y}_{j}+\vec{l}\sqrt[3]{N}$.
The two series are truncated as to assures a relative error below $10^{-7}$,
which is smaller than the relative error of the energy conservation in a
single numerical step.

Of the two dimensionless parameters of the problem, only $\beta$ appears into
the equations of motion. The parameter $\Gamma$ enters through the choice of
the initial data: in fact, while the positions are extracted from a uniform
distribution, the velocities are taken from a Maxwell distribution, and the
temperature is uniquely determined by $\Gamma$. With this choice, at the
beginning of each simulation, the system is out of equilibrium: so there is a
drift of the kinetic energy, and the system reaches a different, random,
temperature. In order to fix the temperature to the desired value, we operate
in this way: after extracting the initial values, we let the system evolve
until equilibrium is reached, i.e. until the kinetic energy appears to
stabilize. We then generate new velocities again with a Maxwell distribution
at temperature $T$, and repeat the process until the kinetic energy appears to
be constant and close to the chosen value.

The equations (3) were integrated using a symplectic splitting algorithm (see
for example Ref. [16]). The total Hamiltonian

$H=\frac{1}{2}\sum_{k}\Big{(}\vec{p}_{k}+\frac{1}{2}\vec{e}_{z}\times\vec{y}_{k}\Big{)}^{2}+\frac{1}{\beta^{2}}V(\vec{y}_{1},\ldots,\vec{y}_{N})\
,$ (5)

where $V(\vec{y}_{1},\ldots,\vec{y}_{N})$ is the Coulomb potential of the
electrons computed according the Ewald prescription, was split as
$H=H_{1}+H_{2}$ where

$H_{1}\stackrel{{\scriptstyle\mathrm{def}}}{{=}}\frac{1}{2}\sum_{k}\Big{(}\vec{p}_{k}+\frac{1}{2}\vec{e}_{z}\times\vec{y}_{k}\Big{)}^{2}\
,\quad
H_{2}\stackrel{{\scriptstyle\mathrm{def}}}{{=}}V(\vec{y}_{1},\ldots,\vec{y}_{N})\
.$ (6)

Now, denoting by $\Phi^{t}_{1}$ the flow determined in the phase space by the
Hamiltonian $H_{1}$ and by $\Phi^{t}_{2}$ the one due to the Hamiltonian
$H_{2}$, the integration algorithm is obtained simply by composition
$\Phi\stackrel{{\scriptstyle\mathrm{def}}}{{=}}\Phi^{\delta
t/2}_{1}\circ\Phi^{\delta t}_{2}\circ\Phi^{\delta t/2}_{1}$, where $\delta t$
is the time step. Such time step was chosen so that the energy conservation
was better then a part over $10^{4}$ in all the simulations. In particular we
choose $\delta t=2\pi C(\beta)\,10^{-4}$, where the factor $C(\beta)$ is taken
equal to $\beta$ for $\beta<1$, and costant equal to 1 for larger $\beta$. In
figura 2 we report the relative error of the energy conservation as a function
of the number of integration step, for a tipical run. The number of step used
was $4\,10^{6}$ for $\beta\leq 2.5$, four times this number for $\beta$
larger.

Figure 2:  Relative error of the energy conservation as a function of the
number of integration steps. Data refers to a case with $\beta=1$ and $N=4096$
particles.

In MD simulations the choice of the number $N$ of particles is always an
issue, all the more in the context of a weakly coupled Plasma. For a very
rough estimation, we made the following considerations: the fundamental cell
of the simulation should have a side larger than the Debye length
$\lambda_{D}$, which in our units reads $\lambda_{D}=\sqrt{1/\Gamma}$. The
first of (2) implies that $L=N^{1/3}$, so that the requirement $L>\lambda_{D}$
in our units becomes $N>\Gamma^{-3/2}$. As the Coulomb force is long range,
the computational cost scales as a power of $N$. With a clever use of the
Ewald summation formula, see Ref. [23], the computational cost scales
(asymptotically) as111It is possible to conceive algorithms with an even
slower asymptotic growth, but for the order of magnitude of $N$ that we are
considering, the Ewald summation formula works best. $N^{3/2}$. So we cannot
afford to choose a value much bigger than $N=4096$. In any case, for
$\Gamma=0.1$ our constraint reads $N>32$ and so with our choice of $N=4096$ we
have $L\simeq 5\lambda_{D}$.

Another interesting length scale which appears in the problem is the Larmor
radius $r_{l}$, i.e., the gyration radius of the electrons determined by the
magnetic Lorentz force due to the external magnetic field $B$. A simple
computation shows that one has $r_{l}/L=\sqrt{2/\Gamma}/\beta N^{1/3}$. For
the smallest value of $\beta$ used in our computations, i.e. $\beta=0.25$, the
Larmor radius is slightly larger than the side of the simulation cell, because
one gets $r_{l}\simeq 1.1L$. Instead $r_{l}$ turns out to be well below the
Debye length, for the largest value of $\beta=10$.

##  3 Results

We recall that the transverse diffusion coefficient is defined by (1), where

$\left|\Delta\vec{x}_{\bot}\right|^{2}\stackrel{{\scriptstyle\mathrm{def}}}{{=}}\frac{\sum_{k}\Big{(}|x_{k}(t)-x_{k}(t_{0})|^{2}+|y_{k}(t)-y_{k}(t_{0})|^{2}\Big{)}}{N}$
(7)

is the mean particles displacement in the plane orthogonal to the magnetic
field $\vec{B}$ (we recall that the latter was taken directed as the $z$–axis)
and the brackets denote the time average along the orbit. To estimate this
quantity we proceed as follows.

Let $\delta t$ be the integration step and $M_{tot}$ the total number of
integration steps. Let also $j$ be an integer smaller than a fixed fraction
$M^{\prime}$ of $M_{tot}$ (we take one sixteenth), then the time averages of
$\langle\left|\Delta\vec{x}_{\bot}\right|^{2}\rangle$ at time $t_{j}=j\delta
t$ were computed as

$\langle\left|\Delta\vec{x}_{\bot}\right|^{2}\rangle(t_{j})=\frac{1}{M}\sum_{h=1}^{M}\left(\frac{1}{N}\sum_{k=1}^{N}\Big{(}|x_{k}(t_{h+j})-x_{k}(t_{h})|^{2}+|y_{k}(t_{h+j})-y_{k}(t_{h})|^{2}\Big{)}\right)\
,$ (8)

where $M=M_{tot}-M^{\prime}$. After plotting
$\langle\left|\Delta\vec{x}_{\bot}\right|^{2}\rangle$ as a function of time,
we fit the tail with a straight line whose angular coefficient (divided by
four) is an estimate of the diffusion coefficient. The values computed in such
way can be found in Table 1 (second entry).

Table 1:  Values of the diffusion coefficient $D_{\perp}$, toghether with the
values of the parameter $A$, $\omega$ and $\gamma$ appearing in the fit of
velocity auto correlation.

$\beta$ |  $D_{\perp}$$~{}^{\mathrm{a}}$ |  $D_{\perp}$$~{}^{\mathrm{b}}$ | $A$ | $\omega$ | $\gamma$
---|---|---|---|---|---
0.25 | 18.1 | 17.9 | 9.90 | 0.253 | 0.169
0.50 | 5.61 | 5.91 | 9.85 | 0.506 | 0.169
0.75 | 2.68 | 2.83 | 9.80 | 0.754 | 0.167
1.00 | 1.56 | 1.63 | 10.1 | 1.01 | 0.164
1.25 | 1.03 | 1.07 | 9.95 | 1.258 | 0.166
2.00 | 0.428 | 0.434 | 9.90 | 2.022 | 0.165
2.50 | 0.327 | 0.342 | 9.75 | 2.518 | 0.171
4.00 | 0.212 | 0.200 | 10.0 | 4.046 | 0.168
5.00 | 0.188 | 0.184 | 10.35 | 5.035 | 0.164
10.0 |  0.133 |  0.146 | 10.30 |  10.06 |  0.157

$~{}^{\mathrm{a}}$ computed from the linear fit of the tail of function (8).

$~{}^{\mathrm{b}}$ computed from the velocity correlation by formula (9).

The whole set of our plots (in log–log scale) can be found in figure 3; as
usual, they display a diffusive (linear) behaviour only after a certain time
(the so called “ballistic jump”). So we also restricted the sets of points to
the latter time window to exclude small times.

Figure 3:  The average square displacement as a function of time. It is
interesting to observe the initial non diffusive behaviour, which becomes more
evident as $\beta$ is increased.

The results of our computations are summarized in figure 1. There, in
logarithmic scale, the numerically computed values of the coefficient
$D_{\perp}$ are reported (full circles) as a function of $\beta$. All the
simulations were performed at the same value of $\Gamma=0.1$. The solid line
is the plot of the function $D_{\perp}=\frac{\gamma A}{\gamma^{2}+\beta^{2}}$,
with $\gamma=0.168$ and $A=9.9$. It can be checked that this function agrees
very well with the numerical results for $\beta<2$. We remark that, for small
values of $\gamma$, this function decrees as $\beta^{-2}$, i.e., for
$\beta<2$, the computed values of $D_{\perp}$ agrees with the prediction of
the kinetic theory. But by further increasing the magnetic field above a value
about $2$, the plot shows a “knee”: the diffusion coefficient appears to obey
a different law. These results are in agreement with figure 5(a) of Ref. [25].
Notice that for such value of $\beta$, the Larmor radius $r_{l}$ is slightly
smaller than Debye length.

Figure 4:  Velocity auto correlation
$\langle\vec{v}_{\perp}(t)\cdot\vec{v}_{\perp}(0)\rangle$, as a function of
time, for different values $\beta=0.5$, $\beta=1.25$, $\beta=2.5$ and
$\beta=5$, of the magnetic field (circles), together with a best fit by a
damped oscillation $A\cos\omega te^{\gamma t}$ (solid line). The values of the
parameter $A$, $\gamma$ and $\omega$ are reported in each panel. The fit is
very good either below and above $\beta=2$. Notice that $A$, $\gamma$ are
essentially independent from $\beta$, while $\omega$ agrees well with the
cyclotron frequency $\omega_{c}$.

The law $D_{\perp}=\frac{\gamma A}{\gamma^{2}+\beta^{2}}$, was not obtained by
interpolation, but by the use of the following argument. Let us introduce the
velocity auto correlation
$\langle\vec{v}_{\perp}(t)\cdot\vec{v}_{\perp}(0)\rangle$ where
$\vec{v}_{\perp}$ in the component of the mean particles velocity transverse
to the magnetic field, and the brackets denote the time averages along the
orbit. Then the diffusion coefficient can be expressed (see the handbook [26])
as

$D_{\perp}=\frac{1}{2}\int_{0}^{+\infty}\langle\vec{v}_{\perp}(t)\cdot\vec{v}_{\perp}(0)\rangle\,\mathrm{d}t$
(9)

whenever the velocity auto correlation decays at $t\to+\infty$ fast enough.
Let us introduce a function $f(t)$ to describe this decay by setting

$\langle\vec{v}_{\perp}(t)\cdot\vec{v}_{\perp}(0)\rangle\simeq\langle\vec{v}_{\perp}(0)\cdot\vec{v}_{\perp}(0)\rangle\cos(\omega_{c}t)f(t)\
.$ (10)

Figure 5:  Autocorrelation of the magnetization (along the magnetic field
direction) as a function of time, for two different values $\beta=0.5$
(square) and $\beta=5$ (circles) of the magnetic field. Notice that, for a
field below the threshold the auto correlation vanishes very fast, while,
above the threshold, autocorrelation vanishes eventually on a totally
different time scale. Figure 6:  Spectrum of the velocity as a function of
$\omega$ in semi logarithmic scale, for two different value of $\beta=0.5$
(circle) and $\beta=5$ (square). The continuous lines are the Fourier
transform of the damped oscillation $2Ae^{-\gamma t}\cos\omega t$. For
$\beta=1.0$ the slope of the spectrum vanishes at $\omega=0$; for $\beta=5$
the slope remains different from zero. This latter case hints to a decay of
the velocity auto correlation as $1/t^{2}$ for $t\to+\infty$ (see text).

Recalling the electronic dynamics, we expect that, due that the gyration along
the magnetic field lines alone,
$\langle\vec{v}_{\perp}(t)\cdot\vec{v}_{\perp}(0)\rangle$ would oscillate with
the cyclotron frequency $\omega_{c}$. However, the electric interaction
between the electrons instead determines a loss of coherence of the electronic
motion, and thus the decaying to zero of the auto correlation. A common choice
is to consider an exponential decay, and thus $f(t)=2Ae^{-\gamma t}$, where
$\gamma$ is a parameter representing the inverse of the decorrelation time.
Then formula (9) gives exactly $D_{\perp}=\frac{\gamma
A}{\gamma^{2}+\beta^{2}}$.

In figure 4 the velocity auto correlation is plotted as a function of time,
for different values of $\beta$, either above and below the threshold
$\beta=2$. As one can check the law (10) is nicely agreed. From the values
reported in table 1, one finds that the parameters $A$ and $\gamma$ are quite
independent from $\beta$, while $\omega$ turns out to be very close to the
cyclotron frequency $\omega_{c}$ (equal to $\beta$ in our units) as expected.
So, taking average values $A=9.9$ and $\gamma=0.168$, and $\omega=\omega_{c}$
the values of the diffusion coefficient $D_{\perp}$ obtained by numerical
computations can be recovered from formula (9).

Now, a peculiar fact happens. In fact, while formula (10) appears to be in
very good agreement with the velocity auto correlation plots for all the
values of $\beta$ considered, the formula $D_{\perp}=\frac{\gamma
A}{\gamma^{2}+\beta^{2}}$ is valid only below a threshold $\beta_{cr}\simeq
2$. This notwithstanding, if we compute the time integral numerically and
inject it into formula (9), the resulting values of $D_{\perp}$ agree very
well with those found using linear regression, as one can check from the Table
1 comparing the values reported in the second and third column. Notice that,
when computing the integral appearing formula (9), one has to truncate the
integral to an upper limit $t_{max}$ chosen carefully, i.e., not to large,
otherwise the numerical errors introduced in computing the auto correlation
$\langle\vec{v}_{\perp}(t)\cdot\vec{v}_{\perp}(0)\rangle$ spoil the result.

A possible explanation of this transition when passing from a weakly
magnetized to a strongly magnetized regime can be given in the spirit of the
paper [9], were it was surmised that at low $\beta$ the plasma is in a fully
chaotic regime, but as $\beta$ is increased a transition to a partially
ordered regime occurs.

It was shown in that paper how, in such a partially ordered regime, a
perturbation theory could apply by adapting the Hamiltonian perturbation
theory developed for system of interest in statistical mechanisc (see for
example Refs. [5, 8, 15, 11]) to the case of plasma. The idea is that, in the
thermodynamical limit, one cannot controll the adiabatic invariants along each
individual trajectory (in the phase space), but it is instead possible to
controll their auto correlations with respect to an invariant measure, showing
that they vanish exponentially slow in the perturbative parameter. Notice
that, in virtue of the linear response theory, such auto correlations
correspond often to important physical observables.

In the above mentioned paper [9], it was considered the case of the component
$M$ of the magnetization along the magnetic field $\vec{B}$, defined (as
usual) as

$M\stackrel{{\scriptstyle\mathrm{def}}}{{=}}\frac{e}{2mc}\frac{1}{N}\sum_{k}(\vec{v}_{k}\times\vec{x}_{k})_{z}\
.$ (11)

Notice that, the auto correlation
$\mathcal{C}_{M}(t)\stackrel{{\scriptstyle\mathrm{def}}}{{=}}<M(t)M(0)>$ is an
important quantity because, according the linear response theory, its Fourier
transform $\hat{\mathcal{C}}_{M}(\omega)$ gives the magnetic susceptibility
$\chi(\omega)$ at the frequency $\omega$ (see for example Refs. [20, 19] or
the Appendix B of Ref. [6]).

In our case, the behavior the auto correlation is different below and above
the threshold. This is shown in Figure 5: for low magnetic field the auto
correlation relaxes to zero, while for high magnetic field it keeps
oscillating and eventually vanishes on a totally different time scale. So the
magnetizazion could be considered an adiabatic invariant of the system, thus
implying that the dynamics remains partially correlated for long times.

A similar mechanism could be at work also in the case of the velocity auto
correlation, even if from our plots this is not as evident as in the case of
the magnetization. In fact, a clue can be obtained by looking at the spectrum
of the velocity auto correlation, i.e., at their cosine Fourier transform, as
shown in figure 6 (in semi logarithmic scale). Notice that $D_{\perp}$ is
simply half the value of the spectrum at $\omega=0$.

For $\beta=0.5$ the slope of the spectrum seems to vanish at $\omega=0$. As
the spectrum is an even function of the frequency this is coherent with the
behavior of a smooth function. Instead, for $\beta=5$, the spectrum slope
seems to remain (negative and) different from zero at $\omega=0$. Now, it can
be very easily shown that, by denoting with $\hat{f}(\omega)$ the Fourier
cosine transform of the function $f(t)$ then, for $t\to+\infty$, one has
$f(t)=\frac{2}{\pi}\frac{1}{t^{2}}\left(-\hat{f}^{\prime}(0)+o(1)\right)$, if
all the derivatives of $\hat{f}(\omega)$ up to the second one are
integrable.222 This can be inferred by the inverse transform formula
$f(t)=\frac{2}{\pi}\int_{0}^{+\infty}\hat{f}(\omega)\cos\omega
t\,\mathrm{d}\omega$ by integrating two times by part and using the
Riemann–Lebesgue lemma on the remainder. So it seems that in our data on the
velocity auto correlation, there is a small component which decays very slowly
(i.e., as $t^{-2}$) to zero. Nevertheless such a component gives a very big
contribution to the diffusion coefficient (more then the double of the one due
to the exponentially decreasing part $2Ae^{-\gamma t}\cos\omega t$).

In any case, at the moment, it is not clear, for what reason, in a less
chaotic regime, the diffusion coefficient apparently decreases, as a function
of $\beta$, at a slower rate with respect to the fully chaotic regime.

##  4 Final considerations

In this work we performed MD simulations of a magnetized electron gas, also
called a One Component Plasma. We have shown that such a system shows a
transition between two different regimes as the value of the parameter $\beta$
is raised above a critical threshold of about $2$.

The transition occurs both on a macroscopic level, with a change in the
diffusive behavior, and on a microscopic level as well: when the parameter is
raised, the system passes from a chaotic state to a partially ordered one.
This is the main finding that we have pointed out, and is a rather new
phenomenon, only remotely addressed in the literature, up to now. As a matter
of fact a similar result, for what concerns the behavior of the diffusion
coefficient, was published quite recently in Ref. [25]. But the authors tried
to explain this phenomenon in the framework of the kinetic theory, looking at
the behavior of the particles during the “collisions”. We refrain from this
approach, and we try to discuss it according to the ergodic theory of
dynamical system, using its standard tools. In particular, our aim is to
understand if the dynamics is truly chaotic or not, and, in this latter case,
what are the consequence for what concerns the macroscopic quantity
characterizing the plasma.

As regards the direct consequences of our results on plasma physics, a strong
objection may be raised on the dependence, we have not explored, of the
results on the number of particles $N$. In particular, in Ref. [25] it was
claimed that a few hundred particles are sufficient in the low $\beta$ regime,
but for high $\beta$ a huge number of particles (above $10^{6}$) is needed. In
particular, they show that the value of the diffusion coefficient slowly
decreases as $N$ is increased. But it seems unlikely that those values would
ever agree with the kinetic law $D_{\perp}\propto\beta^{-2}$, although a Bohm
relation $D_{\perp}\propto\beta^{-1}$ may finally show up in the high $\beta$
regime, as in Ref. [22]. Also, such a high value of $N$ is going towards the
real number of particles, so at this point even the use of periodic boundary
conditions become questionable. Finally, maybe the problem is that the
diffusive regime is not normal, so that the diffusion coefficient is not well
defined. This problem is discussed at length in Ref. [3].

In any case, we think that there are much more fundamental questions to
address about the portability of our results to real plasmas: above all, the
absence of positive ions in our models. We hope to be able to address such a
big issue in the (near) future.

## Data Availability Statement

The data that support the findings of this study are available from the
corresponding author upon reasonable request

## References

  * [1] M.P. Allen and D.J. Tildeslay.  Computer simulation of liquids.  Clarendon Press, Oxford, 1991.
  * [2] Scott D. Baalrud and Jerome Daligault.  Transport Regimes Spanning Magnetization-Coupling Phase Space.  Physical Review E, 96:043202, 2017.
  * [3] R Balescu.  Aspects of Anomalous Transport in Plasmas, volume 18 of Series in Plasma Physics.  Taylor & Francis, 2005.
  * [4] M. E. Caplan and I. F. Freeman.  Precise Diffusion Coefficients for White Dwarf Astrophysics.  Monthly Notices of the Royal Astronomical Society, 505(1):45–49, 2021.
  * [5] A. Carati.  An averaging theorem for hamiltonian dynamical systems in the thermodynamic limit.  Journal of Statistical Physical, 128, 2007.
  * [6] A. Carati, F. Benfenati, and L. Galgani.  Relaxation properties in classical diamagnetism.  Chaos, 21, 2011.
  * [7] A. Carati, F. Benfenati, A. Maiocchi, M. Zuin, and L. Galgani.  Chaoticity threshold in magnetized plasmas: numerical results in the weak coupling regime.  Chaos, 24:013118, March 2014.
  * [8] A. Carati and A. Maiocchi.  Exponentially long stability times for a nonlinear lattice in the thermodynamical limit.  Communication in Mathematical Physical, 314, 2012.
  * [9] A. Carati, M. Zuin, A. Maiocchi, M. Marino, E. Martines, and L. Galgani.  Transition from order to chaos, and density limit, in magnetized plasmas.  Chaos, 22(3):033124, 2012.
  * [10] Jérôme Daligault.  Diffusion in Ionic Mixtures across Coupling Regimes.  Physical Review Letters, 108:225004, 2012.
  * [11] W. De Roeck and F. Huveneers.  Asymptotic localization of energy in nondisordered oscillator chains.  Pure and Applied Mathematics, 68, 2015.
  * [12] D. F. Escande.  Contributions of plasma physics to chaos and nonlinear dynamics.  Plasma Physics and Controlled Fusion, 58(11):113001, 2016.
  * [13] G. G Zimbardo, P. Veltri, and P. Pommois.  Anomalous, quasilinear, and percolative regimes for magnetic-field-line transport in axially symmetric turbulence.  Physical Review E, 61:1940, 2000.
  * [14] P. Gibbon and G. Sutmann.  Long-Range Interactions in Many-Particle Simulation.  In Quantum Simulations of Complex Many-Body Systems: From Theory to Algorithms, volume 10 of NIC Series, pages 467–506. John von Neumann Institute for Computing, Jülich,, 2002.
  * [15] A. Giorgilli, S. Paleari, and T. Penati.  An extensive adiabatic invariant for the klein-gordon model in the thermodynamic limit.  Annales Henri Poincaré, 16, 2015.
  * [16] E. Haire, C. Lubich, and G. Wanner.  Geometric Numerical Integration.  Springer, Berlin, second edition, 2009.
  * [17] M. B. Isichenko.  Effective plasma heat conductivity in ’braided’ magnetic field-ii. percolation limit.  Plasma Physics and Controlled Fusion, 33:809, 1991.
  * [18] M. B. Isichenko.  Percolation, statistical topography, and transport in random media.  Reviews of Modern Physics, 64:961, 1992.
  * [19] Yu.L. Klimontovich.  Statistical Physics.  Taylor & Francis, New York, 1986.
  * [20] R. Kubo.  Statistical-mechanical theory of irreversible processes. i. general theory and simple applications to magnetic and conduction problems.  Journal of the Physical Society of Japan, 12, 1957.
  * [21] C. L. Longmire and M. N. Rosenbluth.  Diffusion of Charged Particles across a Magnetic Field.  Physical Review, 103(3):507–510, 1956.
  * [22] T. Ott and M. Bonitz.  Diffusion in a Strongly Coupled Magnetized Plasma.  Physical Review Letters, 107(13):135003, 2011.
  * [23] J.W. Perram, H.G. Petersen, and S.W.D. Leeuw.  An algorithm for the simulation of condensed matter which grows as the $n^{3/2}$ power of the number of particles.  Molecular Physics, 65:875–893, 1988.
  * [24] D. Saumon, C. E. Starrett, and J. Daligault.  Diffusion coefficients in white dwarfs.  2014\.  arXiv: 1410.7645.
  * [25] Keith R. Vidal and Scott D. Baalrud.  Extended space and time correlations in strongly magnetized plasmas.  Physics of Plasmas, 28:042103, 2021.
  * [26] Gregory H. Wannier.  Statistical Physics.  Dover Publications, New York, revised edition, 2010.
  * [27] G. M. Zaslavsky.  Chaos, fractional kinetics, and anomalous transport.  Physics Reports, 371(6):461–580, 2002.


# MH370 Burst Frequency Offset Analysis and Implications on Descent Rate at
End-of-Flight

Ian D. Holland
Defence Science and Technology Group, Edinburgh, Australia Copyright © 2018
Commonwealth of Australia.

##  I Introduction

On 7 March 2014 at 16:41:43Z111All times given in this article are in UTC or
“Zulu” time, denoted by “Z”., Malaysian Airlines flight MH370 departed Kuala
Lumpur (KL) International Airport bound for Beijing. Less than an hour later,
following the last recorded radio transmission from MH370 at 17:19:30Z, the
plane’s secondary radar transponder went offline. As evidenced by Malaysian
military radar, the plane (registration number 9M-MRO) then veered off course
unexpectedly, back-tracked across the Malaysian Peninsula, and was then
tracked heading northwest from Penang through the Malacca Straits. After
disappearing from radar at 18:22:12Z, it re-established a satellite
communications (SATCOM) link with the Inmarsat satellite I-3F1 at 18:25:27Z.
By analyzing a series of automated messages exchanged via that satellite
between the plane and an Inmarsat ground station in Perth, Australia, it was
determined that the plane continued to fly for six hours, before finally
ceasing message exchange with the ground station at 00:19:37Z on 8 March 2014.
This article discusses specifically the analysis of burst frequency offset
(BFO) metadata from the SATCOM messages. Importantly, it is shown that the
BFOs corresponding to the last two SATCOM messages from the plane at 00:19:29Z
and 00:19:37Z suggest that flight MH370 was rapidly descending and
accelerating downwards when message exchange with the ground station ceased.

An initial analysis by Inmarsat of the SATCOM metadata for MH370 in the last
six hours of flight suggested that MH370 had flown into the Southern Indian
Ocean before SATCOM was ultimately lost (see [1] for further details). As
summarized in [2] an intensive aerial and surface search was undertaken in the
Southern Indian Ocean by an international search team during March and April
2014, with no MH370 related debris found. On 28 April 2014, the aerial search
concluded and the search transitioned to an underwater phase [2]. The
Australian Transport Safety Bureau (ATSB) took responsibility for the
definition of the underwater search zone. It convened an international flight
path prediction working group bringing together experts in satellite
communications, statistical data processing and aviation, in order to estimate
the most likely final location of flight MH370. The group consisted of
representatives from the Australian Defence Science and Technology (DST) Group
and the other organizations listed in the Acknowledgment section of this
article.

New methods of analyzing the Inmarsat data were developed by the group,
resulting in the release of reports concerning the likely final location of
flight MH370 from the ATSB in August 2014 [3], October 2014 [4], and December
2015 [5]. Inmarsat also published an article regarding their contribution to
the flight path reconstruction effort [1]. The DST Group contribution that
assisted in the definition of an extended priority search area in December
2015 [5] has been detailed in [2]. This demonstrated how Bayesian analysis was
used to identify a high probability region of where the plane was believed to
be at the time of last SATCOM transmission (00:19:37Z 8 March 2014). The DST
Group Bayesian method used a prior probability distribution defined by the
Malaysian military radar, a likelihood function describing the relationship
between SATCOM measurements and the aircraft position and velocity during the
flight, and a model of the aircraft dynamics.

It should be noted that this article does not cover the Bayesian method used
in defining the underwater search area. For details on that method, the reader
is referred to [2]. Instead, this article focuses on:

  1. 1.

A brief review of the statistical analysis of BFOs for several previous
flights of 9M-MRO.

  2. 2.

Examination of the effects of the plane’s vertical velocity on the BFO.

  3. 3.

An analysis of the effect of the track angle of the plane on the BFO towards
the end-of-flight.

  4. 4.

The establishment of a BFO trend throughout the last 6 hours of flight so as
to determine an expected BFO at end-of-flight.

  5. 5.

An analysis of the behavior of the frequency oscillator in the plane’s
satellite data unit (SDU) after power outage events such as those believed to
have occurred twice during flight MH370, and how this affects the BFO.

  6. 6.

An analysis that shows that the final two BFO’s are consistent with MH370
being in a rapid descent and accelerating downwards.

The remainder of this article is structured as follows. Section II presents a
brief timeline of key events during the MH370 accident flight. In Section III,
a review of the SATCOM model is provided, along with a brief review of the BFO
statistics. This serves as a summary of work previously presented in [2].
Section IV describes the effects of aircraft position and velocity on the BFO.
The settling behavior of the SDU’s frequency oscillator for 9M-MRO after
power-up is then detailed in Section V. This is important because the SDU is
believed to have undergone a power outage between 00:11Z and 00:19Z on 8 Mar
2014, immediately preceding the last two SATCOM transmissions from MH370. The
effect of oscillator warm-up on the BFO after power-up can be used in bounding
the descent rates for MH370 at 00:19:29Z and 00:19:37Z. In Section VI an
analysis of the descent rate of MH370 based on the last two BFO values is
presented. This analysis derives lower and upper bounds on the descent rate at
00:19:29Z and 00:19:37Z. Conclusions are presented in Section VII.

##  II Timeline of Events During the Accident Flight

A comprehensive description of the events occurring during the MH370 accident
flight is provided in [6]. A summary of salient events to this article is
provided in Table I.

TABLE I: Timeline of Events During MH370 Date and Timestamp(s) | Description | Comments
---|---|---
7 Mar. 2014 16:42Z | MH370 departs Kuala Lumpa Airport | Normal take-off
Last ACARS transmission | 7 Mar. 2014 17:07Z | Everything normal, on path to Beijing
Loss of secondary radar | 7 Mar. 2014 17:21:13Z | Onboard transponder goes inactive; last civilian radar
7 Mar. 2014 18:22:12Z | Last Malaysian military radar contact | MH370 shown tracking NW through Malacca Straits
7 Mar. 2014 18:25Z-18:28Z | Satellite Data Unit (SDU) log-on sequence | If flying level, MH370 still tracking NW
7 Mar. 2014 18:39-18:41Z | Unanswered ground-to-air phone call | If flying level, MH370 tracking southwards
7 Mar. 2014 19:41Z to 8 Mar. 2014 00:11Z | SATCOM metadata approx. hourly from MH370 |  [1, 2] suggest flight into Southern Indian Ocean
8 Mar. 2014 00:19:29Z-00:19:37Z | Partial SDU log-on sequence | Implies rapid descent

##  III Review of SATCOM Model

The accident aircraft was fitted with a SATCOM terminal that used the Inmarsat
Classic Aero system [1], which uses geosynchronous satellites to relay
messages between aircraft and ground stations. During the flight, messages
were passed between the aircraft and a ground receiving station located in
Perth, Australia, via the Inmarsat-3F1 satellite. Figure 1 illustrates the
SATCOM system in use during the flight. The aircraft is referred to as the
Aircraft Earth Station (AES) and the ground receiving unit is referred to as
the Ground Earth Station (GES). Inmarsat-3F1 is a satellite in geosynchronous
orbit at 64.5∘ East longitude and it was used exclusively for the duration of
the flight.

Figure 1: System model of the satellite communication system [2].

An AES is equipped with an SDU comprising a satellite modem with auxiliary
hardware and software. Transmission of data over the satellite is via bursts
which are scheduled to arrive at the GES at a specified time and a given
frequency. As explained in [2, 1], communications from multiple users are
coordinated by the allocation of different time and frequency slots to each
user. This is done without knowledge of individual AES locations or precise
knowledge of the satellite location. Therefore, messages from a given AES
might not arrive at the GES at exactly the expected time, and generally would
arrive slightly later. The difference between the expected time of arrival
(based on a nominal assumed position for the satellite and the AES) and the
actual time of arrival is referred to as the Burst Timing Offset (BTO). The
BTO is a measure of how far the aircraft is from the sub-satellite
position222The sub-satellite position is the point on the earth directly below
the satellite..

The relative velocity between the satellite and the AES, as well as between
the satellite and the GES, leads to a Doppler frequency offset on the signals
received at the GES. Coupled with small frequency offsets inherent in the
reference frequency oscillators in the AES, satellite and GES, this results in
a net difference between the expected and actual frequency of the signal
presented to the modem in the GES for a given user. Frequency compensations
applied onboard the aircraft (aircraft induced Doppler pre-compensation) and
at the ground station (Enhanced Automatic Frequency Correction, which utilizes
the reference signal transmitted from a reference station in Burum,
Netherlands), [2, 1] serve to reduce the possible difference between the
expected and actual frequency of the messages received from the aircraft. The
residual difference between the expected frequency of each communications
burst and the actual received frequency is referred to as the BFO.

###  III-A Review of BFO Statistics

Based on 20 previous flights of 9M-MRO in the week leading up to the accident
flight (see [2] for further details), a histogram was produced for the
difference between the predicted BFO (based on known details of the plane and
the satellite’s position and velocity) and the measured BFO (based on Inmarsat
ground station logs). This difference (i.e. predicted minus measured) is
referred to as the BFO error. The histogram of the BFO error is shown in Fig.
2, along with a Gaussian distribution fit line. It can be seen that the
distribution is somewhat Gaussian. The standard deviation of the BFO error was
found in [2] to be 4.3 Hz. Whilst it is reasonable to apply bounds on the
possible BFO error based on $\pm 3$ standard deviations as was done for the
approximate analysis described in [7], for the purpose of the descent analysis
presented later in Sec. VI, it is assumed the BFO error is strictly bounded on
the larger interval $\left[-28,+18\right]$ Hz, which corresponds to the bounds
of all 2501 observed valid333One outlier was removed as explained in Sec. V.
in-flight BFO error values available from the preceding 20 flights of 9M-MRO.

Figure 2: Histogram of BFO errors for 20 flights of 9M-MRO prior to MH370
(reproduced from Fig. 5.5 of [2]). The mean and standard deviation for the
distribution are 0.18 Hz and 4.3 Hz, respectively [2].

##  IV Effects of Aircraft Position and Velocity on the Burst Frequency Offset

###  IV-A Mathematical Description of the BFO

In [2], the BFO is defined mathematically at time step $k$ as the sum of a
noiseless component $h_{k}^{\mathsf{BFO}}$ and a scalar $w_{k}^{\mathsf{BFO}}$
that represents the BFO noise. The noiseless component of the BFO is defined
in [2] as the sum:

$\displaystyle h_{k}^{\mathsf{BFO}}\left(\mathbf{x}_{k},\mathbf{s}_{k}\right)$
$\displaystyle=$ $\displaystyle\Delta
F_{k}^{\mathsf{up}}\left(\mathbf{x}_{k},\mathbf{s}_{k}\right)+\Delta
F_{k}^{\mathsf{down}}\left(\mathbf{s}_{k}\right)$ (1) $\displaystyle+\delta
f_{k}^{\mathsf{comp}}\left(\mathbf{x}_{k}\right)+\delta
f_{k}^{\mathsf{sat}}\left(\mathbf{s}_{k}\right)$ $\displaystyle+\delta
f_{k}^{\mathsf{AFC}}\left(\mathbf{s}_{k}\right)+\delta
f_{k}^{\mathsf{bias}}\left(\mathbf{x}_{k},\mathbf{s}_{k}\right)\mbox{,}$

where

  * •

$\mathbf{x}_{k}$ denotes the state vector of the aircraft;

  * •

$\mathbf{s}_{k}$ denotes the state vector of the satellite;

  * •

$\Delta F_{k}^{\mathsf{up}}\left(\mathbf{x}_{k},\mathbf{s}_{k}\right)$ is the
uplink (aircraft to satellite) Doppler shift;

  * •

$\Delta F_{k}^{\mathsf{down}}\left(\mathbf{s}_{k}\right)$ is the downlink
(satellite to ground station) Doppler shift;

  * •

$\delta f_{k}^{\mathsf{comp}}\left(\mathbf{x}_{k}\right)$ is the frequency
compensation applied by the aircraft;

  * •

$\delta f_{k}^{\mathsf{sat}}\left(\mathbf{s}_{k}\right)$ is the variation in
satellite translation frequency;

  * •

$\delta f_{k}^{\mathsf{AFC}}\left(\mathbf{s}_{k}\right)$ is the frequency
compensation applied by the ground station receive chain;

  * •

$\delta f_{k}^{\mathsf{bias}}\left(\mathbf{x}_{k},\mathbf{s}_{k}\right)$ is a
slowly varying bias due to errors in the aircraft and satellite oscillators
and processing in the SDU.

By treating the bias $\delta
f_{k}^{\mathsf{bias}}\left(\mathbf{x}_{k},\mathbf{s}_{k}\right)$ as a constant
determined at the source tarmac for any particular flight, as was done in [2]
for MH370, any small time-varying component of the bias during a particular
flight can be considered as part of the BFO noise (indeed this was done when
compiling the results used to obtain the BFO error histogram shown in Fig. 2).
Details regarding the terms $\delta
f_{k}^{\mathsf{sat}}\left(\mathbf{s}_{k}\right)$ and $\delta
f_{k}^{\mathsf{AFC}}\left(\mathbf{s}_{k}\right)$ are provided in [1].
Tabulated values of the sum of these two terms were provided by Inmarsat to
the MH370 Flight Path Reconstruction Group to use in estimating the likely
trajectory flown. These two terms depend on the satellite state
$\mathbf{s}_{k}$ only, and not on the aircraft state $\mathbf{x}_{k}$.
Moreover, the downlink Doppler $\Delta
F_{k}^{\mathsf{down}}\left(\mathbf{s}_{k}\right)$ does not depend on the
location or velocity of the aircraft, and can be calculated given the known
frequency of the downlink and the known satellite state at any given time.

Equation (1) can be then be simplified to:

$\displaystyle h_{k}^{\mathsf{BFO}}\left(\mathbf{x}_{k},\mathbf{s}_{k}\right)$
$\displaystyle=$ $\displaystyle\Delta
F_{k}^{\mathsf{up}}\left(\mathbf{x}_{k},\mathbf{s}_{k}\right)+\delta
f_{k}^{\mathsf{comp}}\left(\mathbf{x}_{k}\right)+\delta
f_{k}^{\mathsf{det}}\left(\mathbf{s}_{k}\right)\mbox{,}$

where $\delta f_{k}^{\mathsf{det}}\left(\mathbf{s}_{k}\right)$ is effectively
a known deterministic value for any time step $k$. The other terms in
(LABEL:eq:BFO_simple) couple the aircraft state $\mathbf{x}_{k}$ by way of the
aircraft position and velocity to the BFO as per the following equations
adapted from [2]444Note that the sign convention used in [2] is opposite to
that used in (3):

$\displaystyle\Delta
F_{k}^{\mathsf{up}}\left(\mathbf{x}_{k},\mathbf{s}_{k}\right)$
$\displaystyle=$
$\displaystyle\frac{F^{\mathsf{up}}}{c}\frac{\left(\mathbf{v}_{\mathsf{s}}-\mathbf{v}_{\mathsf{x}}\right)^{\mathsf{T}}\left(\mathbf{p}_{\mathsf{x}}-\mathbf{p}_{\mathsf{s}}\right)}{\left|\mathbf{p}_{\mathsf{x}}-\mathbf{p}_{\mathsf{s}}\right|}\mbox{,}$
(3) $\displaystyle\delta f_{k}^{\mathsf{comp}}\left(\mathbf{x}_{k}\right)$
$\displaystyle=$
$\displaystyle\frac{F^{\mathsf{up}}}{c}\frac{\left(\hat{\mathbf{v}}_{\mathsf{x}}\right)^{\mathsf{T}}\left(\hat{\mathbf{p}}_{\mathsf{x}}-\hat{\mathbf{p}}_{\mathsf{s}}\right)}{\left|\hat{\mathbf{p}}_{\mathsf{x}}-\hat{\mathbf{p}}_{\mathsf{s}}\right|}\mbox{,}$
(4)

where the dependence on the time step $k$ on the right hand side of the
equations has been removed for simplicity of notation. In equations (3) and
(4), $\left|\cdot\right|$ is the three dimensional Cartesian distance, and:

  * •

$F^{\mathsf{up}}$ is the uplink carrier frequency;

  * •

$c$ is the speed of light;

  * •

$\mathbf{v}_{\mathsf{s}}$ is the velocity vector of the satellite;

  * •

$\mathbf{v}_{\mathsf{x}}$ is the velocity vector of the plane;

  * •

$\mathbf{p}_{\mathsf{s}}$ is the position vector of the satellite;

  * •

$\mathbf{p}_{\mathsf{x}}$ is the position vector of the plane;

  * •

$\hat{\mathbf{v}}_{\mathsf{x}}$ is the SDU’s estimate of the plane’s velocity
vector, which is obtained using the plane’s track angle and ground speed,
whilst assuming the vertical speed is zero;

  * •

$\hat{\mathbf{p}}_{\mathsf{s}}$ is the SDU’s estimate of the position vector
of the satellite, which assumes the satellite is at its nominal orbital slot
of 0 degrees North and 64.5 degrees East;

  * •

$\hat{\mathbf{p}}_{\mathsf{x}}$ is the SDU’s estimate of the position vector
of the plane, which is obtained using the plane’s latitude and longitude,
whilst assuming the plane is at sea level.

###  IV-B Effect of uncompensated vertical velocity

The vertical speed of the plane is not used in the SDU Doppler compensation.
As such, there is a direct contribution of Doppler due to the proportion of
the vertical velocity vector projected onto the radial direction from the
aircraft to the satellite. It is straightforward to understand that if the
plane was directly below the satellite, the vertical velocity vector would be
fully towards or away from the satellite if the plane was climbing or
descending. The direct contribution of the Doppler to the BFO in that case
would be governed by the following standard Doppler equation.

$\displaystyle\Delta
F_{k}^{\mathsf{up}}\left(\mathbf{x}_{k},\mathbf{s}_{k}\right)$
$\displaystyle=$ $\displaystyle\frac{v_{z}\cdot F^{\mathsf{up}}}{c}\mbox{,}$
(5)

where $F^{\mathsf{up}}$ and $c$ are as previously defined, and $v_{z}$ is the
vertical speed of the plane. Substituting an uplink frequency of 1646.6525 MHz
(the uplink frequency stated in [1]) and a vertical velocity of 100 feet per
minute (fpm), equivalent to 0.508 meters per second, equation (5) implies that
the predicted BFO would increase by 2.8 Hz per 100 fpm of climb rate or
decrease by 2.8 Hz per 100 fpm of descent rate if the plane were directly
below the satellite. This is the maximum possible contribution of the plane’s
climb or descent rate to the BFO. In the more general case, equation (5) is
moderated by the sine of the elevation angle $\theta$ from the aircraft to the
satellite. This is expressed as

$\displaystyle\Delta
F_{k}^{\mathsf{up}}\left(\mathbf{x}_{k},\mathbf{s}_{k}\right)$
$\displaystyle=$ $\displaystyle\frac{v_{z}\cdot
F^{\mathsf{up}}\sin(\theta)}{c}\mbox{.}$ (6)

As such, at 00:19Z (the time at which the plane crosses the $7^{th}$ arc
provided by Fig. 3)555The arcs referred to in this paper and other MH370
literature are segments of the rings in Fig. 3., where the elevation to the
satellite is $38.8$ degrees, the contribution to the BFO of climb or descent
rate is reduced to approximately $+1.7$ Hz or $-1.7$ Hz per 100 fpm
respectively.

###  IV-C Effect of track angle on BFO towards end-of-flight

Another factor that needs to be considered in interpreting the BFOs at 00:19Z
is the track angle of the aircraft. In [2], Fig. 5.6, a set of curves were
shown illustrating the relationship of the BFO error and the aircraft track
angle at 18:39Z. That was done to illustrate that under the assumption of
level flight during the unanswered telephone call period of 18:39Z-18:41Z,
MH370 would have been tracking in a southerly direction. The same model can be
used to examine the relationship between the BFO error and the track angle at
00:11Z, when MH370 crossed the 6th arc, just 8 minutes prior to the last
messages were received from MH370.666Note that the 6th arc crossing is
considered here to isolate effects of ground track angle variation as opposed
to vertical velocity, since the measured BFO at the 6th arc is still
consistent with level flight.

Figure 3: BTO rings during the MH370 flight [3]. Note that the latitude of the
arrows shown for each ring is arbitrary.

To generate the curve of BFO error vs. track angle, it is assumed that the
aircraft crossed this arc at 38.67S, 85.11E. A previous sensitivity analysis
has revealed the curve is relatively insensitive to the actual crossing point.
Also, calculations revealed that (as shown in Fig. 5.6 of [2] at 18:40Z) the
track angle of the aircraft influences the BFO to a greater degree for a
faster assumed ground speed. To assess the maximum possible effect of track
angle on the BFO at 00:19Z, a ground speed of 500 kts is assumed, and the BFO
error vs. track angle is then as shown in Fig. 4. Another curve assuming 450
kts ground speed is shown for comparative purposes. It is noted that the peak-
to-peak variation of the BFO difference is similar for both considered ground-
speeds.

Figure 4: BFO errors as a function of track angle at 00:11Z.

###  IV-D BFO Trend During the MH370 Flight

The measured BFOs from 19:41Z to 00:11Z are shown in Fig. 5 together with a
line-of-best-fit. This line is extended forward to the time of the 00:19Z log-
on resulting in an expected BFO of roughly 254 Hz777Note the BFOs observed at
00:19Z were much lower than the expected value.. With reference to the track
angle curves presented in Fig. 4, the 00:11Z BFO error value for the southern-
most track is roughly 6 Hz, meaning the measured BFO of 252 Hz was 6 Hz lower
than the lowest value it could have been (assuming roughly level flight at
approximately 450 kts). If we assume the same error value at 00:19Z, and use
the extrapolated BFO, assuming a south track with level flight and similar
ground speed, the expected BFO value is 260 Hz.

Figure 5: Measured BFOs during the MH370 flight with trend-line extrapolated
back to 18:40Z and forward to 00:19Z.

The linear trend-line in Fig. 5 is also extended back by 1 hour to the time of
the first unanswered satellite telephone call to MH370. The mean of the
measured BFOs during that call attempt is also shown. Since the mean of the
measured BFOs from the 18:39-18:41Z call attempt are in broad agreement with
the linear trend observed in the BFOs from 19:41Z to 00:11Z (for which the
BTOs themselves were consistent with straight and level flight [2]), this
supports the finding in [2] that there were most-likely no major turns after
the unanswered call attempt (see ([2], Fig. 10.5)).

##  V Effect of SDU Startup on the BFO

The BFOs logged for messages from the plane received at the Perth GES at
00:19:29Z and 00:19:37Z 8 March 2014 were much lower than the expected value
of 260 Hz. In the absence of any other factors affecting the BFOs at these
times, this suggests a large uncompensated negative vertical velocity
component in the BFOs. In determining how to interpret this data, it is
important therefore to carefully consider any known other factors that could
effect these BFOs. One potential factor is an SDU oscillator startup
transient. This section analyses the range of expected effect on the BFO
caused by a period of power outage followed by a restart of the SDU, as was a
possibility between 00:11Z and 00:19Z on 8 March 2014 (see [3] for further
details).

###  V-A Detailed Analysis of SDU Startup Effects on the BFO

At 18:25:27Z 7 March 2014, the Inmarsat GES in Perth received a SATCOM log-on
request from 9M-MRO. A series of messages were exchanged in the following few
minutes as part of a standard log-on sequence. The BFOs over those minutes
displayed somewhat unusual behavior in that (barring the first BFO) the BFO
appeared to be exhibiting a transient settling behavior. This is shown in Fig.
6. It was noted in [1] that the spike in BFO observed at this time888When
viewed on a time-scale of hours as shown in Fig. 9 of [1], this behavior looks
like a spike in BFO between 18:25Z and 18:28Z. was not fully understood and
whilst originally attributed to a possible aircraft maneuver it could also be
related to actual frequency changes occurring during the logon sequence. A
subsequent study [8] by the SATCOM sub-group of the MH370 Flight Path
Reconstruction Group revealed this was most likely due to power-on frequency
drift and subsequent stabilization of the oven-controlled crystal oscillator
(OCXO) in the SDU. In [8], a number of different SDUs were tested to
investigate the effects of power outages on the BFO during SDU power on,
triggering a SATCOM logon such as that which occurred for MH370 at 18:25:27Z.
It was found that whilst the frequency settling behavior was different for
each individual SDU, any given SDU displayed repeatable behavior for a fixed
outage duration, and similar settling characteristics for outages of different
duration.

Figure 6: Measured BFOs following the 18:25:27Z log-on message from MH370.

Based on sequences of Inmarsat BFO logs for 9M-MRO concerning four log-ons
after periods of SATCOM outage in the week leading up to the accident flight,
it was found (in [8]) that the SDU in 9M-MRO resulted in a BFO that was
initially too high at the time of log-on, followed by a simple decay over the
next few minutes to a steady-state value. Regarding the log-on event at
18:25:27Z for MH370, it was noted that there was a non-zero bit error rate
(BER) associated with the log-on request at that time. The associated received
signal level and carrier-to-noise density ratio ($C/N_{0}$) were also
unusually low. As such, the first BFO was deemed untrustworthy.999An
alternative explanation for this outlier BFO could be the unknown conditions
of the aircraft between 17:22Z and 18:25Z, for instance the temperature of the
SDU might have been much lower than the expected ambient temperature. Note
that this is not the only time such a behavior was observed for the SDU in
9M-MRO. In the analysis of BFOs presented in [2], there was an outlier in the
BFO (with an error of -170 Hz) found to have occurred for a SATCOM message
associated with 9M-MRO in one of the 20 previous flights of 9M-MRO. This was
discarded as an outlier and not presented in [2]. It too had a non-zero BER
and a low $C/N_{0}$ (37.6 dBHz compared to surrounding values between 41.5 and
42 dBHz). With the removal of the untrustworthy BFO logged for the first
message in the sequence shown in Fig. 6, the 18:25Z log-on from MH370 was also
determined to follow the simple decay trend observed in other instances.

It was noted by DST Group that if the simple BFO decay trend established after
SDU power-up [8] also occurred during the MH370 SATCOM log-on event beginning
at 00:19:29Z, then it would allow bounds to be established on the BFO that
would have been expected for MH370 if the plane was flying level at that time.
This in turn would allow the determination of bounds on the possible descent
rates of MH370 during that final log-on event. In order to build confidence
that the simple BFO decay trend would hold for the 00:19Z log-on event, DST
Group reviewed additional Inmarsat logs for 9M-MRO, corresponding to the
period 22 February to 28 February 2014. Two additional cases were identified
in which 9M-MRO logged back onto an Inmarsat satellite after a sustained
period of SATCOM outage. The sequence of BFOs observed in these two cases,
along with the five already considered by the SATCOM working group in [8], are
shown in Fig. 7, and details about each of the log-ons are given in Table II.

TABLE II: Details of Log-on Sequences Used for Analysis Identifier | Date and Timestamp of Log-on | Duration of Preceding Power Outage | Any Other Comments
---|---|---|---
Log-on 1 | 23rd Feb. 23:57Z | Between 381 and 442 minutes | After scheduled A1 maintenance check, some non-zero BERs
Log-on 2 | 26th Feb. 14:11Z | Between 295 and 354 minutes |
Log-on 3 | 5th Mar. 03:06Z | Between 35 and 95 minutes |
Log-on 4 | 6th Mar. 13:29Z | Between 43 and 103 minutes |
Log-on 5 | 6th Mar. 15:02Z | Between 35 and 92 minutes |
Log-on 6 | 7th Mar. 12:50Z | Between 228 and 288 minutes | Some non-zero BERs
Log-on 7 | 7th Mar. 18:25Z | Between 20 and 78 minutes | First point untrustworthy

The periods of SATCOM outages followed by log-on events were identified from
Inmarsat-provided ground-station logs by identifying sequences of three or
more unsuccessful log-on interrogations to 9M-MRO (suggesting the SDU was
likely powered off) followed by a log-on to the satellite system initiated
from 9M-MRO at some later time. The exact duration of the power-off period for
the SDU was unable to be determined from the logs, however the timestamps
associated with the unsuccessful log-on interrogations and subsequent 9M-MRO
initiated log-ons were used to determine bounds on the outage time as shown in
the $3^{rd}$ column of Table II.101010Note that the bounds shown are derived
under the assumption that the SATCOM outage was exclusively due to a power
outage. It is also possible that the first part of some outages was due to a
blocked line of sight to the Inmarsat satellites. Additional comments about
each log-on sequence are given in the final column. Whilst some of the BERs
for log-ons 1 and 6 were non-zero, the BFOs did not appear to be outliers, and
were therefore considered valid. The uncertainty window in the length of the
outage durations shown is roughly 60 minutes; this relates to timers in the
Inmarsat satellite ground stations that result in a log-on interrogation being
sent to the plane after roughly 60 minutes of SATCOM inactivity. Indeed, it is
this same timing responsible for the roughly 60 minute separation between the
four BTO arcs from 19:41Z to 22:41Z (e.g. [3, 1]). It is also of note that
based on other evidence related to MH370 (e.g. secondary radar transponder
loss of signal at 17:22Z) that SDU power was probably first lost at 17:22Z, in
which case the SDU power outage duration for log-on 7 would be approximately
63 minutes.

Figure 7: Measured BFOs for 7 log-ons of 9M-MRO.

For log-ons 3 and 5, there were some points available about half an hour after
the log-on in each case (noting the plane was on the tarmac for the whole
time) that demonstrate the BFO appears to have fully settled to a steady-state
value. Log-ons 1 to 5 only lasted for less than 1 minute after the log-on
request message. The results for log-ons 6 and 7, which had log-on sequences
lasting a few minutes suggest a settling behavior of approximately 3 minutes
duration. For all observed log-on sequences, there seems to be a simple decay
of BFO over a few minutes after SDU startup and log-on request (see also Fig.
8). It is therefore possible a similar decay was occurring during the SDU log-
on sequence beginning at 00:19:29Z.

In order to establish what the likely range of steady-state equivalent BFOs
would have been at 00:19Z, it is useful to re-plot the curves from Fig. 7
translated vertically such that the log-on acknowledgment has 0 Hz offset.
This is done in Fig. 8, which has also been zoomed in the time-axis for
clarity.

For log-ons 1 to 6, the aircraft was on the ground at an airport. It has been
established that the most common scenario for an on-ground-at-airport log-on
after sustained power outage of the SDU is for the AES to use a different
self-Doppler compensation method than was in use for log-on 7 and the log-on
at 00:19Z. This is due to a lack of navigational input data from the
aircraft’s inertial navigation unit in the few minutes after power-up, and is
referred to as the closed-loop Doppler compensation approach.111111The use of
closed-loop Doppler compensation in this circumstance was established in
consultation with members of the MH370 Flight Path Reconstruction Group,
including representatives from Thales and Inmarsat.

The standard AES Doppler compensation mode when navigational data is available
to the SDU is the mode used for log-on 7 and the log-on at 00:19Z. This
approach was briefly described in Sec. IV-A. It used the aircraft’s own
estimates of its position and horizontal velocity to predict and correct the
AES contributed Doppler shift. This is referred to as open-loop Doppler
compensation.

For log-ons 1 to 6, we work on the premise that closed-loop Doppler
compensation was used.121212Previous analysis had not considered the
difference in Doppler compensation mode, yet still resulted in similar descent
rates to those presented in this paper. In that approach, the frequency of
broadcast P-channel transmissions received by 9M-MRO were used to determine
the AES Doppler contribution. In this approach, all of the difference between
the received and expected P-channel frequency is attributed by the SDU to be
caused by AES motion induced Doppler. As such, if the OCXO frequency is higher
than its nominal settled value, due to incomplete temperature stabilization,
the received frequency will appear lower than expected by the corresponding
amount. When applying the AES Doppler compensation in this case, the transmit
frequency is increased by the difference between the stable and pre-stabilized
oscillator frequency (scaled to account for the difference between P-channel
frequency and the R-channel transmit frequency). Moreover, the transmit
frequency itself (prior to AES Doppler compensation) will already be too high
by the same amount due to the incomplete temperature stabilization. The net
combined effect to the BFO decay observed due the OCXO settling will therefore
be doubled compared to the case where open-loop AES Doppler pre-compensation
is used. Hence, the BFO decay results for log-ons 1 to 6 shown in Fig. 8 have
been halved for the analysis when compared to those inferred from the raw
measured BFOs presented in Fig. 7.

Figure 8: Decay in measured BFOs for 7 log-ons of 9M-MRO. The decay rates for
log-ons 1 to 6 have been adjusted to account for the different AES Doppler
compensation mode used during these log-ons.

Recall that the two messages at 00:19:29Z and 00:19:37Z are a log-on and log-
on acknowledgment respectively. From the data used to plot Fig. 8, it can be
determined that the log-on acknowledgment BFOs are in the range of [0,6] Hz
lower than the log-on BFOs. Also the maximum difference between the log-on BFO
and the settled BFO value is in the range [17,136] Hz. The only log-on that
doesn’t appear to be approaching a settled value in this range is log-on 1,
for which it is anticipated that if more data points were available several
minutes later, as was the case for log-ons 3, 5, 6 and 7, the settling
behavior would be similar to that for log-ons 6 and 7, bearing in mind the
outage duration for log-on 1 is probably at least 381 minutes, and the BFO
already dropped substantially in the first 30 seconds after log-on.

###  V-B Summary of SDU Startup Effects on the BFO

It is likely that the SDU in 9M-MRO lost power sometime after 00:11Z, and then
regained power (presumably due to Auxiliary Power Unit (APU) startup) prior to
00:19Z, leading to the SDU making the log-on attempt at 00:19:29Z. Previous
such events for 9M-MRO have shown that this results in an initially too-high
BFO for the log-on message, followed by a simple decay characteristic to reach
a steady-state BFO after several minutes. In order to establish the relevance
of the BFOs at 00:19:29Z and 00:19:37Z to the motion (including descent rate)
of MH370 at those times, this section has established bounds for the likely
range of steady-state equivalent BFOs at those times. This was done by
analyzing BFOs from seven previous SDU start-up events for a period of several
minutes. The specific results obtained suggest that:

  1. 1.

The recorded BFO for the 00:19:29Z log-on was between 17 and 136 Hz higher
than it would have been if the OCXO in the SDU was in a steady state at that
time;

  2. 2.

The recorded BFO for the 00:19:37Z log-on acknowledge message was between 17
and 130 Hz higher than it would have been if the OCXO in the SDU was in a
steady state at that time.

##  VI Bounding the Descent Rates of MH370

The results from Sections IV and V can be combined to provide bounds on the
descent rate of MH370 implied by the BFOs from the two last SATCOM messages
for the flight. In this section, it is shown how this is done for two
different possibilities that could explain the attempted SATCOM log-on from
9M-MRO at 00:19Z. In [3] the most likely cause of this log-on was stated to be
a power interruption resulting from insufficient fuel and subsequent engine
flameout. It could also have been due to a temporary software failure, a loss
of systems providing critical input to the SDU, or a loss of the SATCOM link
due to aircraft attitude being such that the line-of-sight to the satellite is
blocked. If it was indeed a power interruption to the SDU caused by loss of
fuel and subsequent reboot using the APU, the SDU would be without power for
about one minute. In this case, the results of Sec. V-B need to be considered
when interpreting the last two BFOs.131313Given the relatively short duration
of the power-loss under Hypothesis 1, it is possible that the extent of BFO
decay due to OCXO warm-up would be less than shown for log-ons 1 to 7. This is
essentially covered within Hypothesis 2, which considers a momentary SATCOM
outage. If on the other hand, the power loss was momentary (resulting in a
reset of the SDU) or if the temporary SATCOM outage leading to the log-on
request was due to one of the other listed reasons, there would be no “warm-up
drift” to consider, so the results of Sec. V would not need to be applied.
Both cases are considered separately in the following two sub-sections, and
then combined overall bounds are presented.

###  VI-A Hypothesis 1: SATCOM outage due to insufficient fuel

In the event that the SDU log-on at 00:19:29Z was due to engine flame-out,
followed by a restart of the SDU using power from the APU, the SDU outage
preceding the log-on would have lasted about one minute. This would result in
some cooling of the OCXO in the SDU. The likely effect of the SDU startup and
consequent ‘warm-up drift’ of the OCXO is summarized in Sec. V-B. Table III
presents the recorded BFOs and bounds on the adjusted BFOs to remove these
effects of warm-up drift. In the last column of the table, the bounds are
extended taking into account the BFO noise bounds of $\left[-28,+18\right]$ Hz
established in Sec. III-A.

TABLE III: The Last Two BFOs for MH370 Under Hypothesis 1 Timestamp | Recorded BFO (Hz) | BFO Range if Start-Up Drift Removed (Hz) | Extended Range Considering BFO Noise (Hz)
---|---|---|---
00:19:29Z | 182 | $\left[46,165\right]$ | $\left[28,193\right]$
00:19:37Z | -2 | $\left[-132,-19\right]$ | $\left[-150,9\right]$

As established in Sec. IV-B, the recorded BFO would be roughly 1.7 Hz lower
for every 100 fpm of descent rate. As such, depending on whether the plane was
tracking South or North (minimum or maximum expected BFOs, respectively, see
Sec. IV-C), bounds on the descent rate of MH370 at the times of transmission
corresponding to the last 2 BFOs can be determined by subtracting the values
given in the rightmost column of Table III from the expected BFO values for
level flight tracking South or North, dividing the result by 1.7 Hz and
multiplying by 100 fpm. The expected BFO for a south track is approximately
260 Hz, whilst for a north track it is close to 280 Hz. Using these numbers,
it is straightforward to obtain the bounds shown in Table IV. Note that the
bounds have been rounded to the nearest 100 fpm. Looking at all values in the
table, it can be concluded that irrespective of ground track angle and for
assumed ground speeds less than approximately 500 kts, under Hypothesis 1,
MH370 would have been descending at between 3,900 and 14,800 fpm at 00:19:29Z,
and just 8 seconds later at between 14,800 and 25,300 fpm. These descent rates
are consistent with simulations of an uncontrolled phugoid descent reported in
[7].

TABLE IV: MH370 Descent Rates at 00:19Z 8 March 2014 Under Hypothesis 1 Timestamp | Min. Desc. Rate, South Track | Min. Desc. Rate, North Track | Max. Desc. Rate, South Track | Max. Desc. Rate, North Track
---|---|---|---|---
00:19:29Z | 3,900 fpm | 5,100 fpm | 13,600 fpm | 14,800 fpm
00:19:37Z | 14,800 fpm | 15,900 fpm | 24,100 fpm | 25,300 fpm

###  VI-B Hypothesis 2: SATCOM outage due to some other reason

It is still of interest to determine bounds on the descent rates under the
alternate hypothesis that something else not associated with a power outage
led to the SDU-initiated log-on event at 00:19Z. In this case, the warm-up
drift would not apply. Therefore, the recorded BFOs can be treated normally
(though still subject to BFO noise). In this case, the BFO bounds are as set
out in Table V. Using these bounds, the lower and upper descent rates can be
calculated again, resulting in the bounds presented in Table VI. It can be
concluded from the table that the recorded BFOs indicate that at 00:19:29Z the
plane was descending at between 2,900 fpm and 6,800 fpm. At the time of the
last SATCOM message, it can be concluded under Hypothesis 2 that the descent
rate would have been between 13,800 fpm and 17,600 fpm.

TABLE V: The Last Two BFOs for MH370 Under Hypothesis 2 Timestamp | Recorded BFO (Hz) | BFO Bounds With Noise (Hz)
---|---|---
00:19:29Z | 182 | $\left[164,210\right]$
00:19:37Z | -2 | $\left[-20,26\right]$
TABLE VI: MH370 Descent Rates at 00:19Z 8 March 2014 Under Hypothesis 2 Timestamp | Min. Desc. Rate, South Track | Min. Desc. Rate, North Track | Max. Desc. Rate, South Track | Max. Desc. Rate, North Track
---|---|---|---|---
00:19:29Z | 2,900 fpm | 4,100 fpm | 5,600 fpm | 6,800 fpm
00:19:37Z | 13,800 fpm | 14,900 fpm | 16,500 fpm | 17,600 fpm

###  VI-C Summary of Descent Rate Bounds

Combining the descent rate bounds for Hypotheses 1 and 2, outer bounds can be
established for the descent rate of MH370 at 00:19:29Z and 00:19:37Z.
Specifically, regardless of ground track angle and for ground speeds less than
approximately 500 kts, accounting for possible BFO noise, and regardless of
whether the SATCOM outage between 00:11Z and 00:19Z was due to an SDU power
outage or another reason, the outer bounds on the possible descent rates at
the times of the last two SATCOM messages from MH370 are given in Table VII.

###  VI-D Estimated Downwards Acceleration

When interpreting the bounds presented in Table VII, the different conditions
under which each bound was derived need to be considered. For instance the
lowest rate at 00:19:29Z was derived assuming no period of SDU outage and a
southwards track, whereas the highest descent rate for 00:19:37Z was derived
assuming a one minute outage of the SDU with the maximum BFO decay observed in
the previous 7 SDU outage events for 9M-MRO. Hence, it would not be reasonable
to say that the plane could have been descending at 2,900 fpm at 00:19:29Z and
25,300 fpm 8 seconds later, which would imply a downwards acceleration on the
order of 2,800 fpm per second. Whilst it is not possible to determine a
precise acceleration value, a rough approximation of the average descent rate
over those 8 seconds could for instance be taken using the mid-points of the
bounds at each time, which would result in an average downwards acceleration
of 10,700 fpm in 8 seconds or around 1,300 fpm per second, which equates to
around 6.7 ms-2 or $0.68g$, where $g$ denotes Earth’s gravitational constant,
which is approximately 9.8 ms-2. It is straightforward to see that other
reasonable methods of estimating the downward acceleration (such as taking the
difference between the two minimum descent rates or the two maximum descent
rates at 00:19:29Z and 00:19:37Z) would yield similar acceleration values.

TABLE VII: Range of Possible MH370 Descent Rates at 00:19Z 8 March 2014 Timestamp | Min. Desc. Rate | Max. Desc. Rate
---|---|---
00:19:29Z | 2,900 fpm | 14,800 fpm
00:19:37Z | 13,800 fpm | 25,300 fpm

##  VII Conclusion

This article has discussed in detail all known factors that could have
influenced the BFOs associated with the last two SATCOM messages received from
MH370 at 00:19Z. Lower and upper bounds on its descent rate at this time were
then derived.141414The conclusions regarding descent rate have been made
taking into account BFO noise bounds, all reasonably feasible ground tracks
and speeds, and possible OCXO warm-up drift. It has been implicitly assumed
that there were no otherwise unknown factors that could have affected the last
two BFOs. The downwards acceleration over the 8 second interval between these
two messages was found to be approximately $0.68g$. The derived bounds and
approximate downwards acceleration rate are consistent with simulations of an
uncontrolled descent near the $7^{th}$ arc as reported in [7]. This suggests
that 9M-MRO should lie relatively close to the $7^{th}$ BTO arc.

At the time of submission of this paper, MH370 had not been located. The
search was in indefinite suspension pending ‘credible new evidence pointing to
a specific location for the plane’. Further efforts in data analysis and in
drift modeling of located MH370 debris have been conducted and reported (e.g.
[9, 10, 7]). Note that the descent analysis from this paper was used in [7]
when considering the width of any potential future search area. Some recent
news reports have suggested a US underwater exploration firm named Ocean
Infinity may soon be launching a renewed search effort for MH370. The reader
should refer to the Australian Transport Safety Bureau website [11] and the
general media for updates.

## Acknowledgment

The author would like to thank his Defence Science and Technology Group
colleagues Dr. Neil Gordon, Dr. Samuel Davey, Dr. Jason Williams, and Dr. Mark
Rutten who worked alongside him as part of the MH370 Flight Path
Reconstruction Group, as well as Mr. Balachander Ramamurthy for his support to
the satellite data analysis. Special thanks goes to Dr. Gerald Bolding, both
for his support to the satellite data analysis and his extensive efforts in
reviewing and improving this article. The author also thanks the many other
members of the MH370 Flight Path Reconstruction Group, which comprised
representatives from Inmarsat, Thales, Boeing, US National Transportation
Safety Board (NTSB), UK Air Accidents Investigation Branch (AAIB), Honeywell,
SED, Square Peg Communications Inc., and Panasonic. The group has been
expertly led by the Australian Transport Safety Bureau (ATSB).

## References

  * [1] C. Ashton, A. Shuster Bruce, G. Colledge, and M. Dickinson, “The search for MH370,” _The Journal of Navigation_ , vol. 68, no. 1, pp. 1–22, Oct. 2014\. [Online]. Available: http://dx.doi.org/10.1017/S037346331400068X
  * [2] S. Davey, N. Gordon, I. Holland, M. Rutten, and J. Williams, _Bayesian Methods in the Search for MH370_. SpringerOpen, Aug. 2016.
  * [3] ATSB, “MH370 Definition of Underwater Search Areas,” Australian Transport Safety Bureau, ATSB Transport Safety Report, External Aviation Investigation AE-2014-054, Aug. 2014.
  * [4] ——, “MH370 flight path analysis update,” Australian Transport Safety Bureau, ATSB Transport Safety Report, External Aviation Investigation AE-2014-054, Oct. 2014.
  * [5] ——, “MH370 search area definition update,” Australian Transport Safety Bureau, ATSB Transport Safety Report, External Aviation Investigation AE-2014-054, Dec. 2015.
  * [6] Safety Investigation Team for MH370, “Factual Information : Safety Investigation for MH370,” Ministry of Transport, Malaysia, Tech. Rep. MH370/01/15, 2015. [Online]. Available: http://mh370.mot.gov.my/download/FactualInformation.pdf
  * [7] ATSB, “MH370 - Search and debris examination update,” Australian Transport Safety Bureau, ATSB Transport Safety Report, External Aviation Investigation AE-2014-054, Nov. 2016.
  * [8] “Internal study regarding SATCOM ground-station logs,” MH370 Flight Path Reconstruction Group - SATCOM Subgroup.
  * [9] D. Griffin, P. Oke, and E. Jones, “The search for MH370 and ocean surface drift,” CSIRO, Tech. Rep. EP167888, Dec. 2016.
  * [10] ——, “The search for MH370 and ocean surface drift - Part II,” CSIRO, Tech. Rep. EP172633, Apr. 2017.
  * [11] “ATSB website.” [Online]. Available: http://www.atsb.gov.au

|  Ian D. Holland received the Bachelor of Electronic and Communication
Engineering in 2000 and a Ph.D. in wireless telecommunications in 2005, both
from Curtin University of Technology, Western Australia. Since then he has
held positions in the University of Western Australia, the Institute for
Telecommunications Research at the University of South Australia, EMS Satcom
Pacific and Lockheed Martin Australia. Since January 2011, Ian has been
working as a Research Scientist in Protected Satellite Communications at the
Defence Science and Technology Group, Australia.
---|---