# An Analytic Solution to the Kozai-Lidov Evolution Equations

Stephen H. Lubow1
1Space Telescope Science Institute, 3700 San Martin Drive, Baltimore, MD
21218, USA
Email<EMAIL_ADDRESS>
(Accepted XXX. Received YYY; in original form ZZZ)

###### Abstract

A test particle in a noncoplanar orbit about a member of a binary system can
undergo Kozai-Lidov oscillations in which tilt and eccentricity are exchanged.
An initially circular highly inclined particle orbit can reach high
eccentricity. We consider the nonlinear secular evolution equations previously
obtained in the quadrupole approximation. For the important case that the
initial eccentricity of the particle orbit is zero, we derive an analytic
solution for the particle orbital elements as a function of time that is exact
within the quadrupole approximation. The solution involves only simple
trigonometric and hyperbolic functions. It simplifies in the case that the
initial particle orbit is close to being perpendicular to the binary orbital
plane. The solution also provides an accurate description of particle orbits
with nonzero but sufficiently small initial eccentricity. It is accurate over
a range of initial eccentricity that broadens at higher initial inclinations.
In the case of an initial inclination of $\pi/3$, an error of 1% at maximum
eccentricity occurs for initial eccentricities of about 0.1.

###### keywords:

celestial mechanics – planetary systems – methods: analytic – methods: N-body
simulations – binaries: general

††pubyear: 2021††pagerange: An Analytic Solution to the Kozai-Lidov Evolution
Equations–A

##  1 Introduction

The evolution of a test (massless) particle orbit around a member of a binary
system has been the subject of many studies. Particle orbits that are
initially sufficiently inclined with respect to the orbit of the binary can
undergo oscillations in which inclination and eccentricity are exchanged
(Kozai 1962; Lidov 1962; see review by Shevchenko 2017). An initially highly
inclined circular test particle orbit around a member of a binary can
periodically gain high eccentricity when the inclination decreases. The theory
of such orbits has been extended in many ways to include effects such as
higher order (octupole) corrections to the gravitational effects of the binary
(Naoz et al., 2013; Naoz, 2016), tidal friction within stars (Kiseleva et al.,
1998; Fabrycky & Tremaine, 2007), general relativity (Eggleton & Kiseleva-
Eggleton, 2001; Fabrycky & Tremaine, 2007), and particle mass (Hamers, 2021).

The dynamical effects of such Kozai-Lidov (KL) oscillations have found
application to many areas of astronomy, including asteroids (Kozai, 1962),
artificial satellites (Lidov, 1962; Tremaine & Yavetz, 2014), the
eccentricities of planets (Takeda & Rasio, 2005; Saleh & Rasio, 2009; Dawson &
Chiang, 2014), the production of close binaries and close orbiting (hot
Jupiter) planets (Mazeh & Shaham, 1979; Wu & Murray, 2003; Anderson et al.,
2016; Petrovich & Tremaine, 2016), merging black holes (Blaes et al., 2002;
Miller & Hamilton, 2002; Safarzadeh et al., 2020), blue straggler stars
(Perets & Fabrycky, 2009; Antonini et al., 2016; Fragione & Antonini, 2019),
and the evolution of inclined protostellar disks in binaries (Martin et al.,
2014; Zanazzi & Lai, 2017; Lubow & Ogilvie, 2017).

The analytical studies of KL oscillations are typically carried out by
developing secular orbital evolution equations that average over the binary
and particle orbital periods. The perturbing companion is taken to be far from
the test particle so that its effects can be treated in the quadrupole or
octupole approximation. By making this approximation, the orbit averaging for
secular evolution can be accomplished analytically. The secular equations can
be analyzed by analytic methods to provide insight into the nature of these
oscillations. The analytic studies have emphasized the evolution properties
starting with general initial conditions. Kinoshita & Nakai (2007) provides
analytic solutions in terms of elliptic integrals and an infinite Fourier
series expansion for the longitude of the ascending node. Fabrycky & Tremaine
(2007) determine the accessible states of systems by the taking advantage of
the two constants of motion.

In this paper, we concentrate on the important case of an initially circular
inclined particle orbit within the quadrupole approximation using the secular
equations given in Kiseleva et al. (1998). Section 2 reviews these equations.
Section 3 derives the analytic solution that is expressed in terms of the
simple trigonometric and hyperbolic functions. Section 4 describes the
behavior of the solution at high initial inclination. Section 5 discusses how
the analytic solution extends to low initial inclination. Section 6 discusses
a comparison of the analytic solution with a numerical solution. Section 7
describes some phase portraits. Section 8 discusses the errors made in using
the analytic solution in cases where the initial particle orbit is
noncircular. Section 9 contains the summary.

##  2 Orbital Evolution Equations

Consider a test (massless) particle that orbits a member of a binary with mass
$M_{1}$ and is perturbed by an exterior companion of mass $M_{2}$. The
particle orbit is inclined relative to the binary orbital plane. We adopt the
quadrupole approximation for the secular evolution of the particle orbital
elements. The secular evolution equations for particle eccentricity $e$,
inclination relative to the binary orbital plane $i$, argument of periapsis
$\omega$, and longitude of the ascending node $\Omega$, following equations 2
-5 in Kiseleva et al. (1998), are given by

$\displaystyle\sqrt{1-e^{2}}\frac{de}{dt^{\prime}}$ $\displaystyle=$
$\displaystyle 5e(1-e^{2})\sin^{2}{i}\sin{\omega}\cos{\omega},$ (1)
$\displaystyle\sqrt{1-e^{2}}\frac{di}{dt^{\prime}}$ $\displaystyle=$
$\displaystyle-5e^{2}\sin{i}\cos{i}\sin{\omega}\cos{\omega},$ (2)
$\displaystyle\sqrt{1-e^{2}}\frac{d\omega}{dt^{\prime}}$ $\displaystyle=$
$\displaystyle(5\sin^{2}{\omega}-1)(e^{2}-\sin^{2}{i})$ (3)
$\displaystyle+1-e^{2}+\cos^{2}{i},$
$\displaystyle\sqrt{1-e^{2}}\frac{d\Omega}{dt^{\prime}}$ $\displaystyle=$
$\displaystyle-[1+e^{2}(5\sin^{2}{\omega}-1)]\cos{i},$ (4)

where $t^{\prime}$ is a dimensionless time with $t^{\prime}=t/\tau$ for time
$t$, timescale

$\tau=\frac{2P_{\rm b}^{2}(M_{1}+M_{2})}{3\pi PM_{2}}(1-e_{\rm b}^{2})^{3/2},$
(5)

with binary orbital period $P_{\rm b}$, binary eccentricity $e_{\rm b}$, and
particle orbital period $P$.

In addition, from equations 6 and 7 in Kiseleva et al. (1998), there are two
constants of motion that are

$\displaystyle(1-e^{2})\cos^{2}{i}$ $\displaystyle=$ $\displaystyle{\rm
constant},$ (6)
$\displaystyle\left[5e^{2}\sin^{2}{\omega}+2(1-e^{2})\right]\sin^{2}{i}$
$\displaystyle=$ $\displaystyle{\rm constant}.$ (7)

We take time $t^{\prime}=0$ to be when the eccentricity is maximum. Under the
quadrupole approximation applied in this paper, the oscillation period is
formally infinite with a weak logarithmic singularity in the limit of small
initial eccentricity $e_{\rm init}$. This singularity is removed in octupole
order, as discussed in Kiseleva et al. (1998). The "initial" time
$t^{\prime}_{\rm init}$ is then formally at $t^{\prime}=-\infty$. The particle
orbit is taken to be initially circular, $e_{\rm init}=0$ in the sense of a
limit at $t^{\prime}\rightarrow-\infty$. We are interested in a solution for
$e(t^{\prime})$ for which $e(t^{\prime})>0$ at any finite time. The
eccentricity is then small but nonzero at any early finite time. There is also
a solution for which $e_{\rm init}=0$ at finite initial time. This point is
discussed further in Section 5.

For this case of an initially circular particle orbit these constants reduce
to

$\displaystyle(1-e^{2})\cos^{2}{i}$ $\displaystyle=$
$\displaystyle\cos^{2}{i_{\rm init}},$ (8)
$\displaystyle\left[5e^{2}\sin^{2}{\omega}+2(1-e^{2})\right]\sin^{2}{i}$
$\displaystyle=$ $\displaystyle 2\sin^{2}{i_{\rm init}}.$ (9)

##  3 General Solution for Initially Circular Orbit

Applying Equations (8) and (9) in Equation (1), we obtain a solution for
$e(t^{\prime})$ subject to the initial condition described in the previous
section that can then be used to determine the other orbital elements as a
function of dimensionless time $t^{\prime}$. The result is that

$\displaystyle e\left(t^{\prime}\right)$ $\displaystyle=$
$\displaystyle\frac{q}{\sqrt{6}}\operatorname{sech}{\left(q\,t^{\prime}\right)},$
(10) $\displaystyle\cos^{2}({i(t^{\prime})})$ $\displaystyle=$
$\displaystyle\frac{6\cos^{2}{(i_{\rm
init})}}{6-q^{2}\operatorname{sech}^{2}{(q\,t^{\prime})}},$ (11)
$\displaystyle\sin^{2}({\omega(t^{\prime})})$ $\displaystyle=$
$\displaystyle\frac{4[3(1+\cosh{(2q\,t^{\prime})})-q^{2}]}{12-7q^{2}+3(4+q^{2})\cosh{(2q\,t^{\prime})}},$
(12) $\displaystyle\Omega{(t^{\prime})}$ $\displaystyle=$
$\displaystyle\Omega_{0}-\cos{(i_{\rm
init})}t^{\prime}-\arctan{\left(\frac{q\tanh{(q\,t^{\prime})}}{2\cos{(i_{\rm
init})}}\right)},$ (13)

where $\Omega_{0}$ is the longitude of the ascending node at time
$t^{\prime}=0$ and

$q=\sqrt{1-5\cos{(2i_{\rm init})}}.\\\ $ (14)

Equation (13) for $\Omega$ contains a term that corresponds to a uniform
precession that is present even at low inclinations (as is discussed in
Section 5), as well as a nonuniform precession term that is due to the KL
effect. As discussed in Appendix A, this nonuniform precession term is closely
related to the Gudermannian function and the Mercator projection (see
https://en.wikipedia.org/wiki/Gudermannian_function).

The KL oscillations take place for real values of parameter $q$ defined in
Equation (14). Oscillations occur when $0<q\leq\sqrt{6}$. This constraint is
satisfied over a range of initial tilt angles for which $\cos(2i_{\rm
init})<1/5$ that corresponds to $141^{\circ}\la i_{\rm init}\la 39^{\circ}$
(Kozai, 1962; Lidov, 1962). Equation (11) can also be written as

$\cos^{2}({i(t^{\prime})})=\frac{3(6-q^{2})}{5(6-q^{2}\operatorname{sech}^{2}{(q\,t^{\prime})})},$
(15)

which shows that $\cos^{2}({i(t^{\prime})})\leq 3/5$ for this range of $q$
values. Equation (15) implies that the inclination reaches a minimum value at
$t^{\prime}=0$ that is equal to the Kozai-Lidov critical angle,
$\cos^{2}({i(0}))=3/5$ or equivalently $\cos({2\,i(0}))=1/5$, as is well
known.

For any orbit in the KL regime, notice that by Equation (1), the eccentricity
$e$ is maximum in time ($de/dt^{\prime}=0$) for

$\omega(0)=\pm\frac{\pi}{2}.$ (16)

This value of $\omega(0)$ agrees with the solution given by Equation (12)
which implies $\sin^{2}({\omega{(0)})}=1$. Which of the two values of
$\omega(0)$ is realized depends on initial conditions. (With $\omega=0$, or
$\pi$, derivative $de/dt=0$ can also occur in Equation (1), but corresponds to
an eccentricity minimum because $d^{2}e/dt^{\prime 2}>0$. That follows from
applying Equation (3) to the time derivative of Equation (1).) Notice that
Equation (12) for $\omega$ does not contain an arbitrary constant of
integration, as occurs in Equation (13). The phasing of $\omega(t)$ is
determined at the time of eccentricity maximum ($t^{\prime}=0$), as given by
Equation (16).

From the time derivative of Equation (1) and Equation (16), it follows that
$d\omega/dt(0)>0$. We apply this condition in Equation (12) along with
Equation (16) to determine $\omega(t^{\prime})$. In addition, we determine
$i(t^{\prime})$ from Equation (11) to obtain

$\displaystyle i(t^{\prime})$ $\displaystyle=$
$\displaystyle\arccos{\left(\frac{\sqrt{6}\cos{(i_{\rm
init})}}{\sqrt{6-q^{2}\operatorname{sech}^{2}{(q\,t^{\prime})}}}\right)},$
(17) $\displaystyle\omega(t^{\prime})$ $\displaystyle=$ $\displaystyle\pm\pi
H(\pm t^{\prime})-{\text{sgn}}(t^{\prime})\times$ (18)
$\displaystyle\arcsin{\left(\sqrt{\frac{4[3(1+\cosh{(2q\,t^{\prime})})-q^{2}]}{12-7q^{2}+3(4+q^{2})\cosh{(2q\,t^{\prime})}}}\right)},$

where $H$ is the Heaviside step function. These functions and their
derivatives vary smoothly in time.

At early times ($t^{\prime}\ll 0$) we have that

$\displaystyle e(t^{\prime})$ $\displaystyle=$
$\displaystyle\frac{\sqrt{6}\,q}{3}\exp{(qt^{\prime})},$ (19) $\displaystyle
i(t^{\prime})$ $\displaystyle=$ $\displaystyle i_{\rm init},$ (20)
$\displaystyle\omega(t^{\prime})$ $\displaystyle=$ $\displaystyle-(1\mp
1)\frac{\pi}{2}+$ (21)
$\displaystyle\arcsin{\left(\sqrt{\frac{4}{4+q^{2}}}\right)},$
$\displaystyle\Omega{(t^{\prime})}-\Omega_{0}+\cos{(i_{\rm init})}t^{\prime}$
$\displaystyle=$ $\displaystyle\arctan{\left(\frac{q}{2\cos{(i_{\rm
init})}}\right)}.$ (22)

The KL effect can behave as an instability in eccentricity. A sufficiently
inclined orbit with small initial eccentricity undergoes exponential growth in
time at nearly constant inclination. Equation (19) shows that the eccentricity
grows exponentially at early times with a growth rate $d\ln{e}/dt=q/\tau$ that
agrees with equation 43 of Tremaine & Yavetz (2014) and equation 9 of Lubow &
Ogilvie (2017) 111There is a typo in equation 9 of Lubow & Ogilvie (2017), the
factor $\sqrt{5-3\cos^{2}{i}}$ should be $\sqrt{3-5\cos^{2}{i}}$ . (see also
Tremaine et al. (2009)) The KL effect can also be understood as involving a
resonance, as will be discussed in Section 7.

##  4 Solution at High Initial Inclination

We consider cases where the initial particle orbit is both circular and highly
inclined. We define angle

$\theta=\pi/2-i_{\rm init}$ (23)

and assume that $0<\theta\ll 1$. We consider multiple timescales defined as

$t_{j}=t^{\prime}/\theta^{j}$ (24)

for integer $j$. Equations (10), (17), (18), and (13) respectively
approximately reduce to

$\displaystyle e$ $\displaystyle=$
$\displaystyle\left(1-\frac{5}{6}\theta^{2}\right)\operatorname{sech}{(\sqrt{6}t_{0})},$
(25) $\displaystyle i$ $\displaystyle=$
$\displaystyle\arccos{\left(\sqrt{\frac{1}{\frac{5}{3}+\frac{25}{9}\theta^{2}+6t_{1}^{2}}+\theta^{2}}\right)},$
(26) $\displaystyle\omega$ $\displaystyle=$ $\displaystyle\pm\pi H(\pm
t_{1})-{\text{sgn}}(t_{1})\arcsin{\left(\sqrt{\frac{5+18t_{1}^{2}}{5+45t_{1}^{2}}}\right)},$
(27) $\displaystyle\Omega$ $\displaystyle=$
$\displaystyle\Omega_{0}-t_{-1}-\arctan{\left(3t_{1}\right)}.$ (28)

We have retained each contributing term to at least lowest order in $\theta$
in order to obtain a reasonably accurate solution at all times. We see that
orbital elements $i$, $\omega$, and $\Omega$ undergo large changes on a
relatively short timescale $\theta\,\tau\ll\tau$ from the time of the
eccentricity maximum. The eccentricity varies only on the longer timescale
$\tau$.

Figure 1 plots the orbital elements with $i_{\rm init}=0.45\pi$ and
$\Omega_{0}=0$ using the full analytic expressions given in Section 3 as black
solid lines and the approximate expressions given above by red dotted lines.
The approximate expressions are quite accurate with expansion parameter
$\theta=0.05\pi\simeq 0.16$ in this case. A shown in the figure, $i,\omega,$
and $\Omega$ undergo large changes within a small timescale, while $e$ does
not vary as much on that timescale.

Figure 1: Plot of orbital elements $e$, $i$, $\omega$, and $\Omega$ versus
dimensionless time based on the full analytic solution in black solid lines
given by Equations (10), (17), (18), and (13), respectively, and the
approximate expressions in red dotted lines given by Equations (25) - (28).
The initial inclination is $i_{\rm init}=0.45\pi$ and nodal phase constant
$\Omega_{0}=0$.

##  5 Solution at Low Initial Inclination

For any inclination, there is another solution to Equations (1) - (4) that
begins with $e_{\rm init}=0$, where $t^{\prime}_{\rm init}$ can be a finite or
infinite initial time. For this solution, we have

$\displaystyle e(t^{\prime})$ $\displaystyle=$ $\displaystyle 0,$ (29)
$\displaystyle i(t^{\prime})$ $\displaystyle=$ $\displaystyle i_{\rm init},$
(30) $\displaystyle\Omega{(t^{\prime})}$ $\displaystyle=$
$\displaystyle\Omega_{0}-\cos{(i_{\rm init})}t^{\prime}.$ (31)

Consideration of $\omega(t^{\prime})$ is omitted for this solution, since the
orbit remains circular. For $i_{\rm init}$ within the KL angle range, $q>0$
and the eccentricity given in Section 3 grows exponentially in time for
arbitrarily small but nonzero initial values at finite time (see Equation
(19)). Therefore, this constant zero eccentricity solution is unstable for
$q>0$ and the solution in Section 3 then applies.

Both solutions coincide at $q=0$. That is, the solution in Section 3 and the
above solution are the same for $q=0$. Outside the KL angle range (including
the case of low initial inclination), Equations (29) - (31) provide the
solution. Notice that the solution in Section 3 is valid for all inclination
angles if we take $q$ to be the real part of $q$. In that way, $q=0$ outside
the KL angle range and Equations (10), (11), and (13) reduce to Equations (29)
- (31).

##  6 Numerical Verfification

The analytic solutions for the orbital elements $e$, $i$, $\omega$, and
$\Omega$ given respectively by Equations (10), (17), (18), and (13) were
compared to numerical solutions. The numerical calculations were carried out
using Equations (1) - (4) together with initial conditions provided by the
analytic solutions at time $t^{\prime}=-5$ and $i_{\rm init}=\pi/3$. The
equations were integrated from $t^{\prime}=-5$ to $t^{\prime}=5$. Both
branches of $\omega$ given in Equation (18) were tested. The calculations were
carried out using NDSolve in Mathematica. Over this time interval, the
eccentricity $e$ ranges from about $1.3\times 10^{-4}$ at the endpoints to
about 0.76 at the midpoint. The analytic and numerical results agreed to about
$1\times 10^{-6}$ for all orbital elements throughout this time interval. The
errors are likely due to the limits of the precision in the numerical
integration.

##  7 Phase Portraits

, ,

Figure 2:  Phase portraits of $e\cos{(\omega)}$ versus $e\sin{(\omega)}$ for
three values of the initial orbit inclination $i_{\rm init}$. The heavy black
line is a separatrix and corresponds to the solutions found in Section 3.
Orbits inside (outside) the separatrix all have $\omega_{\rm init}=\pm\pi/2$
(0). The different orbits have different values of $e_{\rm init}$. The green X
marks the maximum eccentricity on the separatrix.

Figure 2 plots some phase portraits of $e\cos{(\omega)}$ versus
$e\sin{(\omega)}$ that result from solving Equations (6) and (7). There is a
separate phase portrait for each of three values of $i_{\rm init}$. For any
point in a plot, the value of the eccentricity $e$ is its distance from the
origin and the argument of periapsis $\omega$ is its polar angle from the
horizontal. The heavy black line in each plot passes through the origin and
therefore corresponds to an orbit that has zero eccentricity at some time. The
time dependence of this orbit is described by the analytic solution given in
Section 3. The heavy line is also a separatrix in each plot.

We consider the initial values with subscript init for each orbit to occur
where the eccentricity is minimum (closest point to the origin along the
plotted orbit). The plotted orbits inside the separatrix are chosen to have
$\omega_{\rm init}=\pm\pi/2$. They are librating orbits, orbits that undergo a
limited range of $\omega$ less that $2\pi$. The orbits outside the separatrix
have $\omega_{\rm init}=0.$ They are circulating orbits, orbits that undergo a
full range of $\omega$ equal to $2\pi$. The different plotted orbits have
different values of $e_{\rm init}$. The plotted orbits inside and outside the
separatrix have nonzero eccentricities at all times.

The existence of the librating orbits reflects the resonant nature of the KL
oscillations (e.g., Malhotra, 2012; Shevchenko, 2017). The libration is
associated with an apsidal-nodal resonance in which the average value of
$d\omega/dt=d\varpi/dt-d\Omega/dt$ is zero over the libration period, where
$\varpi$ is the longitude of the periapsis.

The maximum eccentricity on the separatrix occurs at points marked by a green
X. The orbits inside the separatrix, which all have $\omega_{\rm
init}=\pm\pi/2$, converge to this point with an ever decreasing range of
$\omega$ values.

##  8 Convergence to the Analytic Solution

We examine whether the analytic solutions are stably reached. That is, whether
numerical solutions with different initial conditions would approach the
analytic solutions in Section 3 over time. By doing so, we obtain an estimate
of the accuracy of using the analytic solution given by Section 3 in cases
where the initial eccentricity $e_{\rm init}$ is nonzero.

There are different ways to define the convergence to the analytic solution.
Testing for convergence involves a set of an initial parameters that are
varied and a measure of error at a later time. We consider initial values of
$e_{\rm init}>0$, $i_{\rm init}$, and $\omega_{\rm init}$.

We examine the numerical solutions of Equations (1) - (4) and values of
$e_{\rm init}$ and $\omega_{\rm init}$ that differ from the analytic solution
of Section 3 for a given value of $i_{\rm init}$. Of interest, is whether
those solutions approach the analytic solution over time. In carrying out the
numerical integrations for this test, we do not know in advance the time of
maximum eccentricity. Instead, we take the starting time of the numerical
integrations as $t^{\prime}=0$ and determine the time of peak eccentricity
$t^{\prime}_{\rm peak}$. The comparison with the analytic solutions can then
be made by a simple time shift to redefine time $t^{\prime}_{\rm peak}$ of
numerical solutions to time zero, $t^{\prime}=0$. The analytic solution has an
infinite period, while the numerical solutions have a finite period. We
compare results over the first full period of the numerical solutions centered
at the time when the eccentricity is maximum.

Figure 3:  The plot shows the effect of changing the initial condition on the
argument of periapsis $\omega$ from the value given by the analytic solution
of $\omega_{\rm init}\simeq 0.82$ to $\omega_{\rm init}=0$ (upper plot) and
$\omega_{\rm init}=\pi/2$ (lower plot). The initial inclination is $\pi/3$ in
all cases. The black line shows the analytic solution given by Equation (18).
Nonblack colored lines are the numerical solutions for the argument of
periapsis $\omega$ over the first oscillation in integrating Equations (1) -
(4) as a function of time relative to the time of the eccentricity maximum.
The initial values for eccentricity $e_{\rm init}$ are 0.001 (red), 0.01
(orange), and 0,1 (blue).

###  8.1 Convergence of $\omega(t^{\prime})$

We investigate the convergence of $\omega(t^{\prime})$. We consider cases with
starting values $i_{\rm init}=\pi/3,\omega_{\rm init}=0,\pi/2$, and
$\Omega_{0}=0$. Three values of initial eccentricity are applied that are
$e_{\rm init}=0.001,0.01,$ and $0.1$. Orbits with $\omega_{\rm init}=0$ are
outside the separatrix as shown in Figure 2, while orbits with $\omega_{\rm
init}=\pi/2$ are inside the separatrix. The analytic solution has the same
initial inclination, but with $e_{\rm init}=0$ and $\omega_{\rm
init}=\arcsin{(2\sqrt{2/15})}\simeq 0.82$ at $t^{\prime}=-\infty$ as given by
Equation (21). The results plotted in Figure 3 show that the numerical results
for $\omega(t^{\prime})$ are close to the analytic results midway through the
oscillation, as expected because $\omega(0)=\pm\pi/2$ in all cases as seen in
Equation (16).

$\omega(t^{\prime})$ reaches a value close to the analytic solution and breaks
away at later times due to periodicity. The time interval over which the
numerical results agree well with the analytic results increases with
decreasing $e_{\rm init}$. Therefore, the numerical results for
$\omega(t^{\prime})$ converge towards the analytic result as $e_{\rm init}$
approaches zero, for both values of $\omega_{\rm init}$. However, the infinite
time available as $e_{\rm init}$ approaches zero is an artifact of the
quadrupole approximation. In reality, the convergence would be limited in the
case that the higher order gravitational moments are taken into account.

###  8.2 Convergence of $e(0)$

We apply a measure of error that is the difference between the maximum
eccentricity of an orbit with some set of initial conditions from the maximum
eccentricity of the analytic solution of Section 3 that has the same value of
$i_{\rm init}$. In particular, we want to know how rapidly this eccentricity
difference approaches zero as $e_{\rm init}$ approaches zero. Recall that we
define $t^{\prime}=0$ to be the time at which the maximum eccentricity occurs.
We are then comparing solutions at the same time $t^{\prime}=0$. We know from
Equation (16) that this occurs at $\omega=\pm\pi/2$, as seen in Figure 2.

We denote the analytic solution for maximum eccentricity given by the analytic
solution as

$e_{\rm amax}(i_{\rm init})=\frac{\sqrt{1-5\cos{(2i_{\rm init})}}}{\sqrt{6}}$
(32)

that follows from Equation (10). The maximum eccentricity for an orbit with
$e_{\rm init}$ nonzero is denoted by $e_{\rm nmax}(e_{\rm init},i_{\rm
init},\omega_{\rm init})$. We determine the error

$\Delta e(e_{\rm init},i_{\rm init},\omega_{\rm init})=|e_{\rm nmax}(e_{\rm
init},i_{\rm init},\omega_{\rm init})-e_{\rm amax}(i_{\rm init})|,$ (33)

Error $\Delta e$ can be understood in terms of Figure 2. This function
measures the distance between the green X and the points on the other orbits
at $\omega=\pm\pi/2$. We are interested in how much the orbits in Figure 2
converge from "initial" values near the origin to the values near the green X.
As is evident from these plots, there is significant convergence for orbits
both inside and outside the separatrix. As is also seen in Figure 2, the
convergence is stronger with increasing $i_{\rm init}$.

Figure 4: Plot of function $f$ that is defined by Equation (35).

We determine $\Delta e(e_{\rm init},i_{\rm init},\omega_{\rm init})$
analytically using a series expansion in $e_{\rm init}$. For the orbits
outside the separatrix with $\omega_{\rm init}=0$ as shown in Figure 2, we
have to lowest order that

$\Delta e(e_{\rm init},i_{\rm init},\omega_{\rm init}=0)=f(i_{\rm
init})\,e_{\rm init}^{2}$ (34)

where

$f(i_{\rm init})=\frac{25\sin^{2}{(2i_{\rm init}})}{2\sqrt{6}(1-5\cos{(2i_{\rm
init})})^{3/2}}.$ (35)

The errors decrease rapidly with decreasing initial eccentricity as $e_{\rm
init}^{2}$. Function $f(i_{\rm init})$ is plotted in Figure 4. For $i_{\rm
init}\simeq\pi/2$, we have that

$f(i_{\rm init})=\frac{25}{18}\left(\frac{\pi}{2}-i_{\rm init}\right)^{2},$
(36)

which shows the error drops rapidly as the initial inclination approaches
$\pi/2$. At the KL critical angle, $i_{\rm crit}=0.5\arccos{(0.2)}\simeq
0.217\pi\simeq 39.2^{\circ}$, the denominator vanishes in Equation (35). Just
above this critical angle, we have that

$f(i_{\rm init})=\frac{1}{6^{1/4}4(i_{\rm init}-i_{\rm crit})^{3/2}},$ (37)

which shows that the error grows as $i_{\rm init}$ approaches $i_{\rm crit}$,
as seen in Figure 4. In this regime, $\Delta e\ll 1$ for $e_{\rm
init}\ll(i_{\rm init}-i_{\rm crit})^{3/4}$.

Orbits inside the separatrix with $\omega_{\rm init}=\pi/2$ all have the same
maximum eccentricity, independent of initial eccentricity and inclination and
therefore $\Delta e=0$, as seen in Figure 2. Orbits that begin inside the
separatrix in Figure 2 and have $\omega_{\rm init}\neq\pi/2$ cross that
separatrix. Therefore, the heavy black line is not actually a separatrix for
such orbits. For orbits that begin inside the separatrix in Figure 2 and have
$\omega_{\rm init}\simeq\pi/2$, we obtain in a low order series approximation
that

$\Delta e(e_{\rm init},i_{\rm init},\omega_{\rm init})=f(i_{\rm
init})\,\left(\frac{\pi}{2}-\omega_{\rm init}\right)^{2}\,e_{\rm init}^{2}$
(38)

with $f$ defined by Equation (35).

Figure 5: Plot of fractional error function $\Delta e/e_{\rm amax}$ defined by
Equations (32) and (33) as a function of $e_{\rm init}$ with $i_{\rm
init}=\pi/3$. The solid lines are determined numerically for the five
indicated values of $\omega_{\rm init}$. The red, orange, and green lines are
for orbits that begin outside the separatrix in the middle panel of Figure 2.
The purple and blue lines are for orbits that begin inside that separatrix.
The red dashed line is based on the analytic approximation given by Equation
(34). The purple and blue dashed lines are based on the analytic approximation
given by Equation (38).

Figure 5 plots the fractional error at maximum eccentricity $\Delta e(e_{\rm
init},i_{\rm init}=\pi/3,\omega_{\rm init})/e_{\rm amax}$ given by Equations
(32) and (33) as a function of $e_{\rm init}$ for five values of $\omega_{\rm
init}$. The solid lines plot the numerically determined $\Delta e/e_{\rm
amax}$ values. The dashed lines are based on the analytic approximations
provided by Equations (34) and (38) for orbits that begin outside and inside
the separatrix in the middle panel of Figure 2, respectively. The analytic
approximations agree well with the numerical values of $\Delta e/e_{\rm amax}$
for small $e_{\rm init}$, as expected. The $\Delta e/e_{\rm amax}$ values for
orbits that begin inside the separatrix (purple and blue) are smaller than
those outside (red, orange, and green).

In Figure 5 the values of $\Delta e/e_{\rm amax}$ at fixed $e_{\rm init}$
decrease monotonically with increasing $\omega_{\rm init}$ for
$0\leq\omega_{\rm init}<\pi/2$. Also, $\Delta e/e_{\rm amax}$ is symmetric
about $\omega_{\rm init}=0$ and about $\omega_{\rm init}=\pi/2$. That is,

$\displaystyle\Delta e(e_{\rm init},i_{\rm init},\ \omega_{\rm init})$
$\displaystyle=$ $\displaystyle\Delta e(e_{\rm init},i_{\rm init},-\omega_{\rm
init}),$ (39) $\displaystyle\Delta e(e_{\rm init},i_{\rm init},\ \omega_{\rm
init})$ $\displaystyle=$ $\displaystyle\Delta e(e_{\rm init},i_{\rm
init},\pi-\omega_{\rm init}).$ (40)

Consequently, the results in Figure 5 cover values of $\omega_{\rm init}$ in
all four quadrants, $0\leq\omega_{\rm Init}<2\pi$.

Figure 6: Plot of $e_{\rm init}/\sqrt{\Delta e_{\rm init}/e_{\rm amax}}$ as a
function of $i_{\rm init}$ that is defined by Equation (41).

As seen in Figure 5, the errors in using the analytic solution with nonzero
$e_{\rm init}$ are largest for the orbits with $\omega_{\rm init}=0$.
Combining Equations (32), (34), and (35), we then estimate the $e_{\rm init}$
value that is small enough to reach a given level of fractional error $\Delta
e_{\rm init}/e_{\rm amax}$ as

$e_{\rm init}=\frac{\sqrt{2}(1-5\cos{(2i_{\rm init}))}}{5\sin{(2i_{\rm
init})}}\sqrt{\frac{\Delta e}{e_{\rm amax}}}.$ (41)

Figure 6 plots the ratio $e_{\rm init}/\sqrt{\Delta e_{\rm init}/e_{\rm
amax}}$ as a function of $i_{\rm init}$. To achieve a small fractional error
$\Delta e_{\rm init}/e_{\rm amax}=0.01$ with an initial tilt of $i_{\rm
init}=\pi/4$ requires $e_{\rm init}\la 0.03$, for $i_{\rm init}=\pi/3$
requires $e_{\rm init}\la 0.1$, and for $i_{\rm init}=4\pi/9=80^{\circ}$ only
requires $e_{\rm init}\la 0.5$.

##  9 Summary

This paper considers the case of an initially circular orbit of a test
particle around a member of a binary system. The orbit is significantly
inclined with respect to the binary orbital plane. Such a particle undergoes
Kozai-Lidov oscillations that have been the subject of many previous studies.
The companion is assumed to lie on an orbit that is far outside the orbit of
the test particle. Under this (quadrupole) approximation, Equations (10),
(17), (18), and (13) provide an exact analytic solution to the nonlinear
secular evolution of the particle’s orbital elements $e$, $i$, $\omega$, and
$\Omega$ given by Equations (1) - (4). The solution is expressed in terms of
simple trigonometric and hyperbolic functions of time. In the case that the
particle orbit inclination angle deviates slightly from being perpendicular to
the binary orbital plane by amount $0<\theta\ll 1$, the inclination, argument
of periapsis, and longitude of the ascending node undergo large changes on a
timescale that is shorter than the eccentricity evolution timescale by about a
factor of $\theta$ (Section 4).

The analytic solution extends to all inclinations by taking the real part of
$q$ in Equation (14) (Section 5). As discussed in Section 6, the analytic
solution agrees well with a numerically determined solution. For given initial
inclination value $i_{\rm init}$, the analytic solution determines
$\omega_{\rm init}$ given by Equation (21). Numerical solutions that begin
with other values of $\omega_{\rm init}$ approach the analytic solution for
$\omega(t^{\prime})$ with small values of initial eccentricity $e_{\rm init}$
(see Figure 3).

In Section 8.2 we determine the error in using the maximum eccentricity
provided by the analytic solution for cases with nonzero initial eccentricity
as function of $e_{\rm init},i_{\rm init}$, and $\omega_{\rm init}$. The
errors are quadratic in $e_{\rm init}$ for small $e_{\rm init}$ (see Figure
5). The errors increase as $i_{\rm init}$ approaches the KL critical angle and
drop rapidly as $i_{\rm init}$ approaches $\pi/2$ (see Figure 4). The initial
eccentricity required to reach a given fractional error in maximum
eccentricity depends on the initial inclination (see Figure 6). In the case of
an initial inclination of $\pi/3$, an error of 1% at maximum eccentricity
occurs for initial eccentricities of about 0.1. The analytic solution in
Section 3 provides good accuracy for a range of initial conditions that
broadens at higher initial inclinations.

## Acknowledgements

I thank the referee for suggestions that led to improvements in the analysis
of convergence to the analytic solution. I acknowledge support through NASA
XRP grant 80NSSC19K0443. I thank Rebecca Martin and Dan Romik for useful
discussions.

## Data availability

The data underlying this article will be shared on reasonable request to the
author.

## References

  * Anderson et al. (2016) Anderson K. R., Storch N. I., Lai D., 2016, MNRAS, 456, 3671
  * Antonini et al. (2016) Antonini F., Chatterjee S., Rodriguez C. L., Morscher M., Pattabiraman B., Kalogera V., Rasio F. A., 2016, ApJ, 816, 65
  * Blaes et al. (2002) Blaes O., Lee M. H., Socrates A., 2002, ApJ, 578, 775
  * Dawson & Chiang (2014) Dawson R. I., Chiang E., 2014, Science, 346, 212
  * Eggleton & Kiseleva-Eggleton (2001) Eggleton P. P., Kiseleva-Eggleton L., 2001, ApJ, 562, 1012
  * Fabrycky & Tremaine (2007) Fabrycky D., Tremaine S., 2007, ApJ, 669, 1298
  * Fragione & Antonini (2019) Fragione G., Antonini F., 2019, MNRAS, 488, 728
  * Hamers (2021) Hamers A. S., 2021, MNRAS, 500, 3481
  * Kinoshita & Nakai (2007) Kinoshita H., Nakai H., 2007, Celestial Mechanics and Dynamical Astronomy, 98, 67
  * Kiseleva et al. (1998) Kiseleva L. G., Eggleton P. P., Mikkola S., 1998, MNRAS, 300, 292
  * Kozai (1962) Kozai Y., 1962, AJ, 67, 591
  * Lidov (1962) Lidov M. L., 1962, Planet. Space Sci., 9, 719
  * Lubow & Ogilvie (2017) Lubow S. H., Ogilvie G. I., 2017, MNRAS, 469, 4292
  * Malhotra (2012) Malhotra R., 2012, Encyclopedia of Life Support Systems by UNESCO, 6, 55
  * Martin et al. (2014) Martin R. G., Nixon C., Lubow S. H., Armitage P. J., Price D. J., Doğan S., King A., 2014, ApJL, 792, L33
  * Mazeh & Shaham (1979) Mazeh T., Shaham J., 1979, A&A, 77, 145
  * Miller & Hamilton (2002) Miller M. C., Hamilton D. P., 2002, ApJ, 576, 894
  * Naoz (2016) Naoz S., 2016, ARA&A, 54, 441
  * Naoz et al. (2013) Naoz S., Farr W. M., Lithwick Y., Rasio F. A., Teyssandier J., 2013, MNRAS, 431, 2155
  * Perets & Fabrycky (2009) Perets H. B., Fabrycky D. C., 2009, ApJ, 697, 1048
  * Petrovich & Tremaine (2016) Petrovich C., Tremaine S., 2016, ApJ, 829, 132
  * Safarzadeh et al. (2020) Safarzadeh M., Hamers A. S., Loeb A., Berger E., 2020, ApJ, 888, L3
  * Saleh & Rasio (2009) Saleh L. A., Rasio F. A., 2009, ApJ, 694, 1566
  * Shevchenko (2017) Shevchenko I. I., 2017, The Lidov-Kozai Effect - Applications in Exoplanet Research and Dynamical Astronomy. Astrophysics and Space Science Library Vol. 441, Springer, doi:10.1007/978-3-319-43522-0
  * Takeda & Rasio (2005) Takeda G., Rasio F. A., 2005, ApJ, 627, 1001
  * Tremaine & Yavetz (2014) Tremaine S., Yavetz T. D., 2014, American Journal of Physics, 82, 769
  * Tremaine et al. (2009) Tremaine S., Touma J., Namouni F., 2009, AJ, 137, 3706
  * Wu & Murray (2003) Wu Y., Murray N., 2003, ApJ, 589, 605
  * Zanazzi & Lai (2017) Zanazzi J. J., Lai D., 2017, MNRAS, 467, 1957

##  Appendix A Nonuniform Nodal Precession: The Gudermannian Function and the
Mercator Projection

Figure 7: Plot of the Gudermannian function $\operatorname{gd}(x)$.

We discuss here the nonuniform nodal precession term that is associated with
Kozai-Lidov oscillations and appears in Equation (13)

$\Omega_{\rm{KL}}(t^{\prime})=-\arctan{\left(\frac{q\tanh{(q\,t^{\prime})}}{2\cos{(i_{\rm
init})}}\right)}.$ (42)

The Gudermannian function is defined as

$\operatorname{gd}(x)=\int_{0}^{x}\operatorname{sech}(u)\,du$ (43)

and can be expressed as

$\operatorname{gd}(x)=2\arctan{\left(\tanh{\left(\frac{x}{2}\right)}\right)}.$
(44)

The function is plotted in Figure 7. In general, Equation (42) cannot be
expressed exactly in terms of Equation (44). However, it can be related in an
approximate way by requiring the function derivatives match at $t^{\prime}=0$
and that the function values match at large $|t^{\prime}|$. The result is that

$\Omega_{\rm{KL}}(t^{\prime})\simeq
c_{1}\operatorname{gd}(c_{2}\,t^{\prime}),$ (45)

where

$\displaystyle c_{1}$ $\displaystyle=$
$\displaystyle-\frac{2}{\pi}\arctan{\left(\frac{q}{2\cos{(i_{\rm
init})}}\right)},$ (46) $\displaystyle c_{2}$ $\displaystyle=$
$\displaystyle\frac{\pi q^{2}}{4\cos{(i_{\rm
init})}\arctan{\left(\frac{q}{2\cos{(i_{\rm init})}}\right)}}.$ (47)

The approximation is fairly accurate for typical parameters in the Kozai-Lidov
regime. In particular, Equation (45) is an exact equality for $q=2\cos(i_{\rm
init})$ in Equation (42) that occurs for $i_{\rm
init}=\arccos{(\sqrt{3/7})}\simeq 49.1^{\circ}$.

The Mercator projection is frequently used for constructing maps of the
surface of the Earth. The projection transforms points on the surface with
longitude and latitude $(\lambda,\phi)$ to an $(x,y)$ Cartesian system (the
map). It is based on a (conformal) mapping that preserves angles between
lines. It does not preserve area, resulting in the familiar problem that the
high latitude country Greenland appears to be comparable in size to the entire
continent of Africa. It can be shown that the inverse Mercator projection of
latitude that maps $y$ to $\phi$ is a Gudermannian function of $y$, see
https://en.wikipedia.org/wiki/Mercator_projection. Therefore, the
$\Omega_{\rm{KL}}$ is related to $t^{\prime}$ by an appropriately scaled
Mercator projection with $\Omega_{\rm{KL}}$ playing the role of $\phi$ and
$t^{\prime}$ playing the role of $y$. The stretched $y$-extent of Greenland by
the Mercator projection (the large value of $dy/d\phi$ at high lattitude
$\phi$) is analogous to the stretched time $t^{\prime}$ for larger values of
$\Omega_{\rm{KL}}$ (the large value of $dt^{\prime}/d\Omega_{\rm{KL}}$ at
larger $\Omega_{\rm{KL}}$). This effect is seen in Figure 7 as the flatness of
the Gudermannian function (the inverse Mercator projection) at large $|x|$.


# Filtering Module on Satellite Tracking

1st Wafi, Moh Kamalul  Department of Engineering Physics
Institut Teknologi Sepuluh Nopember
Surabaya, Indonesia
<EMAIL_ADDRESS>
###### Abstract

The scope of satellite has increasingly attained as one of the most
challenging topics due to the attraction of elaborating the outer space. The
satellite, as a means of collecting data and communicating, needs a proper
calculation so as to maintain the movement and its appearance. The concept of
the proposed research lies in the mathematical model along with certain
noises. The mathematical model is started by initial two variable states,
constituting a radius and an angle, with no process noise on it. These two
states then are formulated with certain assumption of noises in terms of the
range and the scaled angle deviations from them in turn. Keep in mind that
those two noises are mutually independent and their covariance are considered.
the model is defined as Algebraic Riccati Equation (ARE) along with Kalman
filter algorithm, from the estimation, the steady-state estimator, the
computational of gain matrix to the stability of the predictor. The findings
show that, as for the two pairs of states, the performance of the estimation
can follow the state with just slight fluctuations in the first a fifth of a
thousand iterations. With respect to the Mean Square Error (MSE), both noises
are around 0.2 for the four states.

###### Index Terms:

Kalman Filter, Satellite Tracing Tracking

##  I Introduction

Tracking an object is becoming more challenging and it has been studying to
get the precise position while tracking it. The object refers to the satellite
and it has increasingly attained as one of the most challenging topics due to
the attraction of elaborating the outer space. The proposed concept of doing
it is to use Kalman filter as conducted by [1] and [2] which presented the
reduced order of the Riccati differential equation, the tractable of the
object in terms of mathematical model, and the the ease in the real
implementation in turn. This also stimulates to upgrade the classic Kalman
filter in order to obtain another best estimation method as done by [3],
comprising the upgrade of gradient decent in terms of error covariance. The
classic Kalman filter [4] is emerged so as to compare the method in [7] in
terms of the mean square estimation error (MSEE) along with its average over
certain iterations.

The basic concept of the classic Kalman is set from [4] while the elaborating
of the basic is well-presented in [6] saying the various possibility of any
engineering and non-engineering being perturbed by any noises including
tracking a moving object. The initial foundation has been inspired by [9], [4]
along with [5] for the concept of filtering.

The classic is matched with the Micro-Kalman Filter developed by [7] which
then was upgraded in the following paper presented in [8]. This algorithm is
the beginning of the possibility of Distributed Kalman Filter (DKF). This
paper proposes the algorithm compared to the classic in terms of the mean
square error and overall performance of the states.

##  II Problem Formulation

This paper presents a mathematical modelling on satellite tracking as shown in
Fig. 1 with the following qualitative objectives:

  1. 1.

The problem is initiated from showing the right state-space in terms of
nominally circular solution by considering small deviations of the state
variables.

  2. 2.

Those small deviations are measured in terms of observations on earth so that
the two possible are suggested, the noisy measurements of both the range
deviations from $R$ and the scaled angle deviations from $\omega t$.

  3. 3.

The states estimation are derived either with or without initial condition
with respect to both noises along with the mean square error (MSE)

##  III Mathematical Model

From the designed objectives, the basis construction of this paper is begun
from the state-space itself with two additional deviations proposed. The
following is to design the suitable Kalman filter with certain parameters
related to the this paper.

###  III-A State-Space Representation

A unit human-made celestial body rotates in the subspace under the influence
of an inverse square law force field. The motion of the mass is governed by a
couple of second order differential equation for the ideal rotational move of
radius $r(t)=R$ at time $t$

$\ddot{r}(t)=r(t)\dot{\theta}^{2}-\frac{G}{r^{2}(t)}$ (1)

$G_{0}$ is a physical constant enabling circular orbits being defined as
$R^{3}\omega^{2}$ and that is also influenced by a certain angle
$\theta(t)=\omega t$, such that:

$\ddot{\theta}(t)=-2\dot{\theta}\frac{\dot{r}(t)}{r(t)}$ (2)

To initiate the linearised equations about the nominally circular solution,
the deviation $r$ from the ideal move has to be small and those are expressed
as $\textbf{x}=\left[x_{1},x_{2},x_{3},x_{4}\right]^{T}$. The state variables,
due to two pairs of $r$ and $\theta$, are defined as $x_{1},\cdots,x_{4}$,
with $x_{1}=r-R$, $x_{2}=\dot{r}$, $x_{3}=R\left(\theta-\omega t\right)$, and
$x_{4}=R(\dot{\theta}-\omega)$ in turn. Keep in mind that the process noise is
ignored from the set linearised equations.

###  III-B The Construction of Noises

It is important to set two sorts of disturbances in measuring those two
parameters $r$ and $\theta$ at sample times $t=kh$ with $k=1,2,\cdots,N$. The
mean of such noise in terms of the range deviations from $R$ (type 1) is
designed as zero whilst its covariance is known with certain parameter of
$\varphi$, such that:

$y_{1}\left(t\right)=x_{1}(t)+v_{1}(t)\rightarrow
v_{1}\sim\mathcal{N}(0,\varphi)$ (3)

Another noisy measurement is the scaled angle deviations from $\omega t$ (type
2) along with certain covariance of $\psi$:

$y_{3}\left(t\right)=x_{3}(t)+v_{3}(t)\rightarrow
v_{3}\sim\mathcal{N}(0,\psi)$ (4)

$v_{1}$ is a Gaussian noise with mean $\hat{v}_{1}(0)=0$ and covariance
$\mathcal{Q}_{v_{1}}(0)=\varphi$ whereas $v_{3}$ is equal to
$\hat{v}_{3}(0)=0$ and $\mathcal{Q}_{v_{3}}=\psi$. Those two types affect the
matrix of $H$ so that the the research is based on the two models leading to
the different simulations.

###  III-C Kalman Estimation

The problems of estimating variables that are not directly available and
without making any assumption on the stationary of the stochastic processes
are addressed. Recalling the states of Kalman filter by taking into account
$n$-vector process $x_{k}$ and $m$-vector process of observation $y_{k}$
generated by the following model, such that:

$\displaystyle\hat{x}(kh)$
$\displaystyle=Fx\left[\left(k-1\right)h\right]+u^{s}+q\left[\left(k-1\right)h\right]$
(5) $\displaystyle\hat{y}_{1,3}$
$\displaystyle=H\hat{x}(kh)+u^{0}+v_{1,3}\left[\left(k-1\right)h\right]$ (6)

$q\sim\mathcal{W}\mathcal{G}\mathcal{N}(0,\delta)$ is a White Gaussian Noise
with mean $\hat{q}_{1}(0)=0$ and covariance $\mathcal{Q}_{q}(0)=\delta$. The
criteria is the same as that of for $y_{1,3}(kh)=Hx(kh)+u^{0}+v_{1,3}$,
constituting
$v_{1,3}\sim\mathcal{W}\mathcal{G}\mathcal{N}(0,\varphi\>\textrm{or}\>\psi)$
with equal $\hat{v}_{1,3}$ and $\mathcal{Q}_{v_{1,3}}$ which is non-singular.
$q(\bullet)$ and $v_{1,3}(\bullet)$ are mutually independent with $x_{0}$
whereas $F,H,c^{T},V_{1}$, and $V_{2}$ are known. As an initial state,
$x_{0}$, with Gaussian random, its mean and covariance are respectively
written as:

$E\left[x_{0}\right]=\hat{x}_{0|0}\quad\textrm{and}\quad\mathcal{Q}_{x}(0)=P_{0|0}$
(7)

A prevalent of the steady-state filter aimed to either time-invariant or time-
varying systems along with non-stationary noise covariance is time-varying
Kalman filter. Bear in mind that the estimate of
$\textbf{x}\left[(k-1)h|(k-1)h\right]$ is described as
$\hat{\textbf{x}}\left[kh|(k-1)h\right]$ obtained from the preceding
measurements with $y\left[(k-1)h\right]$ whereas the updated estimate
according to the last measurement $y(kh)$ is explained as
$\hat{\textbf{x}}\left[(kh|kh)\right]$. Note that the up-to-date estimate is
$\hat{\textbf{x}}\left[(kh|kh)\right]$ and the iteration then predicts the
following time-updated at $(k+1)h$. The measurement update hence employs this
estimation according to the one-step-ahead predictor $y\left[(k+1)h\right]$
which then is used in the error term, the difference between the measured and
estimated, defined by
$y\left[(k+1)h\right]-c\hat{\textbf{x}}\left[(k+1)h|kh)\right]$. The
innovation of gain $K$ is set to be minimum with respect to the covariance of
the estimation error, such that:

$\xi(kh)=\textbf{E}\left[v_{1}(kh)v_{1}(kh)^{T}\right]$ (8)

The expected value of the second type noise is also written as
$\zeta(kh)=\textbf{E}\left[v_{3}(kh)v_{3}(kh)^{T}\right]$.

###  III-D Steady-State Gain Matrices

The equations to each time-step upgrade $K\left[(k+1)h\right]$ and
$P\left[(k+1)h\right]$ entail the inversion of a matrix that, assuming a
$n$-by-$n$ matrix, requires $O(n^{3})$ operation (in general) and hence the
complexity enumeration of each iteration escalates with the cubic power of the
dimension of the state vector. This would be practically very appealing to be
able to replace the time-varying matrices $K\left[(k+1)h\right]$ and
$P\left[(k+1)h\right]$ with constant $K(kh)$ and $P(kh)$ calculated off-line
previously. It would obviously provide a rise to a sub-optimal predictor but
would allow to handle in practice high-dimensional problems, so that:

$\hat{\textbf{x}}\left[(k+1)h|kh\right]=F\textbf{x}\left[kh|kh\right]+Ke(kh)$
(9)

with $e=y-\hat{y}$. The steady-state gain matrix is simply computed as common
Kalman $K=FPH^{T}\left(HPH^{T}+\psi\right)^{-1}$. Moreover, the matrix
$P=P(kh)=P\left[\left(k+1\right)\right]$ is a solution of the Algebraic
Riccati Equation (ARE):

$P=F\left[P-PH^{T}\left(HPH^{T}+\psi\right)^{-1}HP\right]F^{T}+\varphi$ (10)

In case of multiple solutions of the ARE, it is necessary to choose the
positive semi-definite one. Recalling the equations of the following steady-
state Kalman predictor:

$\hat{y}\left[(k+1)h|kh\right]=H\hat{\textbf{x}}\left[(k+1)h|kh\right]$ (11)

If $K$ stabilizes $F-KH$, the solution of the ARE is stabilizing. However,
this does not say anything about the existence of a positive (definite or
semi-definite) solution of ARE. Then, it is worth asking under what conditions
the recursive Riccati equation converges, that is, the ARE has (at least) one
positive semi-definite solution and by recalling that $P$ is the state
prediction error covariance matrix.

###  III-E Micro-Kalman Filter $\mu$KF

The algorithm is based on research conducted by [7] and [8] being called as
information form which was then to examine the possibility of constructing the
Distributed Kalman Filter (DKF). The state of the plant is the same as the
classical, such that:

$\displaystyle x_{k+1}$ $\displaystyle=A_{k}x_{k}+B_{k}\vartheta_{k}$ (12)
$\displaystyle z_{k}$ $\displaystyle=H_{k}x_{k}+\chi_{k}$ (13)

the schemes of $\vartheta_{k}$ and $\chi_{k}$ are set as
$\mathcal{W}\mathcal{G}\mathcal{N}$ with additional properties as in the
following:

$\displaystyle\textbf{E}(\vartheta_{k}\vartheta_{l}^{T})$
$\displaystyle=Q_{k}\bar{\delta}_{kl}$
$\displaystyle\textbf{E}(\chi_{k}\chi_{l}^{T})$
$\displaystyle=R_{k}\bar{\delta}_{kl}$ (14)

initial error covariance $P_{k}$ and an inverse-covariance matrix $\hat{S}$
are combined to obtain matrix $M$, such that:

$\displaystyle S_{k}=H^{T}_{k}R^{-1}_{k}H_{k}$ (15) $\displaystyle
M^{-1}=(P^{-1}+S)$ (16)

whereas matrix of $K$ is influenced by the $M$-matrix as follows:

$K_{k}=M_{k}H^{T}_{k}R^{-1}_{k}$ (17)

The equation of state estimation is then altered from the common
$\hat{x}_{k}=\hat{x}_{k|k-1}+K_{k}(y_{k}-H\hat{x}_{k|k-1})$ to:

$\displaystyle\hat{x}m$
$\displaystyle=\hat{x}_{k|k-1}+M_{k}(H^{T}_{k}R^{-1}_{k}y_{k}-H^{T}_{k}R^{-1}_{k}H\hat{x}_{k|k-1})$
(18) $\displaystyle\hat{x}_{k|k-1}$ $\displaystyle=\bar{x}_{0}$ (19)

from the equation, it can be also written as $\hat{S}$ and a new measurement
$z_{k}=H^{T}_{k}R^{-1}_{k}y_{k}$

$\hat{x}_{m}=\bar{x}+M_{k}(z_{k}-\hat{S}\bar{x})$ (20)

compared to the classic Kalman filter, the prediction of $\hat{x}$ and $P$ are
becoming the update of this $\mu$KF such that:

$\displaystyle P^{+}_{k}$ $\displaystyle=AM_{k}A^{T}+B\delta B^{T}$ (21)
$\displaystyle\bar{x}^{+}$ $\displaystyle=A\hat{x}_{m}$ (22)

The algorithm is another, not the sequence of the classical, method and is
compared to classic (centralized) Kalman filter in terms of estimation error.

Algorithm 1 Filtering Module

Assumption:

$x_{0}\sim\mathcal{N}(\acute{x}_{0},\tau)$;
$q\sim\mathcal{N}(\bar{q},\delta)$;$v_{1}\sim\mathcal{N}(\bar{v}_{1},\varphi)$;
$v_{3}\sim\mathcal{N}(\bar{v}_{3},\psi)$

Initialization:

$x_{0}=\bar{x}_{0}$ and $P_{0}=\tau$

for $j=1$ to $\phi$ do

for $k=1$ to $N$ do

1\. Centralized Kalman Filter $\displaystyle\hat{x}_{k|k-1}$
$\displaystyle=Fx_{k-1}+u^{s}$ $\displaystyle\hat{y}_{k|k-1}$
$\displaystyle=H\hat{x}_{k|k-1}+u^{0}$ $\displaystyle P_{k|k-1}$
$\displaystyle=FP_{k-1|k-1}F^{T}+\delta$ $\displaystyle T$
$\displaystyle=HP_{k|k-1}H^{T}+\textrm{noise}\rightarrow\varphi\>\textrm{or}\>\psi$
$\displaystyle P_{k}$
$\displaystyle=P_{k|k-1}-P_{k|k-1}H^{T}(T)^{-1}HP_{k|k-1}$ $\displaystyle
K_{k}$ $\displaystyle=P_{k|k-1}H^{T}(T)^{-1}$ $\displaystyle\hat{x}_{k}$
$\displaystyle=\hat{x}_{k|k-1}+K_{k}(y_{k}-\hat{y}_{k|k-1})$ $\displaystyle
M_{k}$ $\displaystyle=(P^{-1}+\hat{S}_{k})^{-1}$ $\displaystyle\hat{x}_{m}$
$\displaystyle=\bar{x}_{0}+M_{k}(\hat{z}_{k}-\hat{S}_{k}\bar{x}_{0})$
$\displaystyle P_{k}^{+}$ $\displaystyle=FM_{k}F^{T}+B\delta B^{T}$
$\displaystyle\bar{x}^{+}$ $\displaystyle=F\hat{x}_{m}$
$\displaystyle\beta_{k}$ $\displaystyle=|x_{k}-\hat{x}_{k}|$
$\displaystyle\gamma_{k}$ $\displaystyle=|x_{k}-\bar{x}^{+}_{k}|$

end for

Mean Square Estimation Error (MSEE): $\displaystyle\kappa_{j}$
$\displaystyle=\frac{1}{N}\sum_{1}^{N}\beta^{2}$ $\displaystyle\Gamma_{j}$
$\displaystyle=\frac{1}{N}\sum_{1}^{N}\gamma^{2}$

end for

Averaged MSEE: $\displaystyle\Xi_{\kappa}$
$\displaystyle=\frac{1}{\phi}\sum_{1}^{\phi}\kappa$
$\displaystyle\Xi_{\Gamma}$
$\displaystyle=\frac{1}{\phi}\sum_{1}^{\phi}\Gamma$

TABLE I: The Mean Square Estimation Error (MSEE) $\Gamma_{k}$ with noise type 1 $v_{1}$ State $\bar{x}^{+}$ | Iteration $\phi_{i}$
---|---
1 | 2 | 3 | 4 | 5 | 6 | 7 | 8 | 9 | 10 | Averaged
$x_{1}$ | 0.0012 | 0.0015 | 0.0010 | 0.0043 | 0.0019 | 0.0015 | 0.0014 | 0.0011 | 0.0016 | 0.0012 | 0.0017
$x_{2}$ | 0.0023 | 0.0023 | 0.0038 | 0.0043 | 0.0076 | 0.0048 | 0.0052 | 0.0041 | 0.0070 | 0.0061 | 0.0048
$x_{3}$ | 0.0120 | 0.0032 | 0.0040 | 0.1271 | 0.0101 | 0.0059 | 0.0088 | 0.0039 | 0.0082 | 0.0081 | 0.0191
$x_{4}$ | 0.0022 | 0.0034 | 0.0019 | 0.0075 | 0.0037 | 0.0022 | 0.0042 | 0.0030 | 0.0042 | 0.0037 | 0.0036
TABLE II: The Mean Square Estimation Error (MSEE) $\Gamma_{k}$ with noise type 2 $v_{3}$ State $\bar{x}^{+}$ | Iteration $\phi_{i}$
---|---
1 | 2 | 3 | 4 | 5 | 6 | 7 | 8 | 9 | 10 | Averaged
$x_{1}$ | 0.0207 | 0.0365 | 0.0095 | 0.0202 | 0.0116 | 0.0150 | 0.0085 | 0.0163 | 0.0191 | 0.0117 | 0.0169
$x_{2}$ | 0.0232 | 0.0395 | 0.0065 | 0.0205 | 0.0114 | 0.0157 | 0.0074 | 0.0159 | 0.0208 | 0.0112 | 0.0172
$x_{3}$ | 0.0155 | 0.0180 | 0.0129 | 0.0073 | 0.0090 | 0.0080 | 0.0081 | 0.0100 | 0.0114 | 0.0084 | 0.0109
$x_{4}$ | 0.0514 | 0.0779 | 0.0240 | 0.0389 | 0.0251 | 0.0346 | 0.0183 | 0.0353 | 0.0419 | 0.0240 | 0.0371

##  IV Simulation Results

The initialization of the simulation is presented with the following
parameters. Supposing $R=1,\,\omega=1,\,G_{0}=1$ for $k=1,\cdots,N$ with
$N=1000$ in both types. Moreover, assume that the probability density of the
state at the initial time $t=0$ is described in the following equations, such
that:

$x_{0}=[0.1,0,0,0]\sim\mathcal{N}\left(\acute{x},\tau\right)$ (23)

where mean is $0$ and its covariance is $0.1\times I_{4\times 4}$. Recalling
the matrix of the states $\left[x_{1},x_{2},x_{3},x_{4}\right]^{T}$, the one-
step-ahead from those four states are defined as follows:

$\displaystyle\dot{x}_{1}$ $\displaystyle=\dot{r}\rightarrow x_{2}$
$\displaystyle\dot{x}_{2}$ $\displaystyle=\ddot{r}\rightarrow
3\omega^{2}x_{1}+2\omega x_{4}$ $\displaystyle\dot{x}_{3}$
$\displaystyle=R(\dot{\theta}-\omega)\rightarrow x_{4}$
$\displaystyle\dot{x}_{4}$ $\displaystyle=R\ddot{\theta}\rightarrow-2\omega
x_{2}$

where it can be concluded that the matrix of A would be:

$\begin{bmatrix}\dot{x}_{1}(t)\\\ \dot{x}_{2}(t)\\\ \dot{x}_{3}(t)\\\
\dot{x}_{4}(t)\\\ \end{bmatrix}=\underbrace{\begin{bmatrix}0&1&0&0\\\
3\omega^{2}&0&0&2\omega\\\ 0&0&0&1\\\ 0&-2\omega&0&0\\\
\end{bmatrix}}_{\textbf{A}}\underbrace{\begin{bmatrix}x_{1}(t)\\\ x_{2}(t)\\\
x_{3}(t)\\\ x_{4}(t)\\\ \end{bmatrix}}_{\textbf{x}}$ (24)

from those two types of noise, it is showed the matrices of $y_{1}$ and
$y_{3}$ are $\left[1,0,0,0\right]$ and $\left[0,0,1,0\right]$ respectively
along with divergent noises of $\varphi=0.1$ and a hefty $\psi=0.5$. It then
determines the matrix $F$ and the row vector $c^{T}$ in the equations relating
$x(kh)$ and $y(kh)$ to $x\left[\left(k-1\right)h\right]$ and the appropriate
noise variable. The value of $F$ with $\Delta t=h=0.01$ for either types is
written as follows:

$F=\begin{bmatrix}1.0001&0.0100&0&0.0001\\\ 0.0300&1.0000&0&0.0200\\\
-0.0000&-0.0001&1&0.0100\\\ -0.0003&-0.0200&0&0.9998\\\ \end{bmatrix}$ (25)

while the value of $c^{T}_{1}$ and $c^{T}_{2}$ is then described the same as
that of $H$.

Figure 1(a) highlights the true and the estimated states along with the
measurement of the $\mu$KF for those four states being perturbed by type 1
noisy measurement and their MSEE(s) are illustrated in Table I. The Figure
1(e) presents the error among the states in Figure 1. The comparison on Table
III $(\textbf{type 1})$ says that the almost there is no gap between the
$\mu$KF and the classic Kalman filter which highlights the different just a
slight $10^{-4}\times[1,1,2,0]^{T}$. Whilst, in terms of the bigger noisy
measurement type 2 with a hefty $0.5$, the patterns of result are the same as
those of the preceding graphs with respect to MSEE and AMSEE which are no much
different compared to the CKF shown in Table II and Table III respectively.
Table III $(\textbf{type 2})$ portrays the $10^{-4}\times[4,3,6,10]^{T}$
compared to the classic. The error in Figure (red) 2(e) are the conclusion of
the figure from 2(a) to 2(d) which have the bigger span due to the bigger
error set.

TABLE III: AMSEE comparison of $\Xi_{\kappa}$ and $\Xi_{\Gamma}$ with either noises State | type 1 | type 2
---|---|---
$\Xi_{\kappa}$ | $\Xi_{\Gamma}$ | $\Xi_{\kappa}$ | $\Xi_{\Gamma}$
$x_{1}$ | 0.0017 | 0.0016 | 0.0169 | 0.0165
$x_{2}$ | 0.0048 | 0.0047 | 0.0172 | 0.0169
$x_{3}$ | 0.0191 | 0.0193 | 0.0109 | 0.0103
$x_{4}$ | 0.0036 | 0.0036 | 0.0371 | 0.0361

##  V Conclusion

The initial mathematical model of the classic and the micro has been proposed
along with the construction of two noisy measurements. From the results, it
can be concluded that $\mu$KF Kalman filtering is being able to track the
satellite with mean square estimation error (MSEE) being shown in Table I and
Table II perturbed by either noises. This algorithm results in almost exactly
the same as classic Kalman filter as depicted in Table III.

## Acknowledgment

Thanks to Professor Richard B. Vinter from the Imperial College London who has
taught me in the lecture leading to finishing this paper and to LPDP
(Indonesia Endowment Fund for Education) Scholarship from Indonesia.

## References

  * [1] F. A. Faruqi and R. C. Davis, ”Kalman Filter Design for Target Tracking”, IEEE Transactions on Aerospace and Electronics Systems, Vol. AES-16, No. 4 July 1980
  * [2] B. Esmat, ”An Adaptive Kalman Filter for Tracking Maneuvering Targets”, Journal of Guidance, Control, and Dynamics, September, Vol. 6, No. 5: pp. 414-416
  * [3] W. Yuh-Shyang, S. You and N. Matni, ”Localized Distributed Kalman Filters for Large-Scale Systems”, ELSEVIER - IFAC Papers online, 50-1 (2017) 10742 - 10747
  * [4] R. E. Kalman, ”A New Approach to Linear Filtering And Prediction Problems”, ASME Journal of Basic Engineering, 1960.
  * [5] B. D. O. Anderson and J. B. Moore, Optimal Filtering. New York: Dover, 2005.
  * [6] E. Mansour, ”Book Reviews” on Optimal Filtering, IEEE Transactions on Systems, Man, And Cybernetics, Vol. SMC-12, No. 2, March/April 1982.
  * [7] R. Olfati-Saber, ”Distributed Kalman Filter with Embedded Consensus Kalman Filters”, Proceedings of the 44th IEEE Conference on Decision and Control, and the European Control Conference, Seville, Spain, Dec 12 - 15, 2005.
  * [8] R. Olfati-Saber, ”Distributed Kalman Filtering for Sensor Networks”, Proceedings of the 46th IEEE Conference on Decision and Control, New Orleans, LA, USA, Dec. 12 - 14, 2007.
  * [9] M B. Rhudy, R. A. Salguero and K. Holappa, ”A Kalman Filtering Tutorial for Undergraduate Students”, International Journal of Computer Science and Engineering Survey (IJCSES), Vol.8, No.1, February 2017.

(a)

(b)

(c)

(d)

(e)

Figure 1: These figures present the results of the filtering according to four
different states of position (a), velocity (b), angle (c), and angular
velocity (d) whereas the last figure (e) portrays the error between the true
and the estimated. The charts are influenced by noise type 1 whilst (a) is
given additional information of measurement

(a)

(b)

(c)

(d)

(e)

Figure 2: The following figures present the results of the filtering according
to four different states of position (a), velocity (b), angle (c), and angular
velocity (d) whereas the last figure (e) portrays the error between the true
and the estimated. The charts are influenced by noise type 2 whilst (c) is
given additional information of measurement


# A new technique to provide realistic input to CME forecasting models

Nat Gopalswamy1 Sachiko Akiyama2 Seiji Yashiro2 Hong Xie2 1 Solar Physics
Laboratory, NASA Goddard Space Flight Center, Greenbelt, MD 20771, USA
email<EMAIL_ADDRESS>
2Dept. of Physics, The Catholic University of America, Washington DC 20064,
USA

(2017)

###### Abstract

We report on a technique to construct a flux rope (FR) from eruption data at
the Sun. The technique involves line-of-sight magnetic fields, post-eruption
arcades in the corona, and white-light coronal mass ejections (CMEs) so that
the FR geometric and magnetic properties can be fully defined in addition to
the kinematic properties. We refer to this FR as FRED (Flux Rope from Eruption
Data). We illustrate the FRED construction using the 2012 July 12 eruption and
compare the coronal and interplanetary properties of the FR. The results
indicate that the FRED input should help make realistic predictions of the
components of the FR magnetic field in the heliosphere.

###### keywords:

coronal mass ejection, flux rope, reconnection.

††volume: 335††journal: Space Weather of the Heliosphere: Processes and
Forecasts††editors: C. Foullon & O. Malandraki, eds.

##  1 Introduction

Coronal mass ejections (CMEs) are the most important players in space weather
because they cause the severest of geomagnetic storms and accelerate energetic
particles to GeV energies. The intensity of geomagnetic storms as measured by
the $Dst$ index primarily depends on the magnitude of the south-pointing out-
of-the-ecliptic ($B_{z}$) component of the interplanetary magnetic field and
the speed ($V$) with which the structure impacts Earth’s magnetosphere: $Dst$
= $-$0.01$V$$\mathopen{|}B_{z}\mathclose{|}$$-$32 nT ([WuLepping02, Wu &
Lepping 2002]; [Gopalswamy_etal08, Gopalswamy et al. 2008]). One of the vexing
problems in space weather has been the prediction of southward (negative)
$B_{z}$ of the interplanetary field that reconnects with Earth’s magnetic
field to produce geomagnetic storms. While there has been reasonable progress
in the prediction of the arrival time of CMEs, $B_{z}$ prediction has been
very limited due to the lack of realistic input to global MHD models that
track CMEs into the heliosphere and provide asymptotic values of the CME
parameters including the magnetic field. Since most CMEs arriving at Earth
have a flux rope (FR) structure, it is inevitable that the global MHD models
should use a FR input at the near-Sun boundary. Current models use a pressure
pulse (which has no magnetic content) or ad hoc magnetic structures such as a
spheromak ([OdstrcilPizzo99, Odstrcil & Pizzo 1999]; [Jin_etal17, Jin et al.
2017]; [Shen_etal14, Shen et al. 2014]; [Shiota_etal14, Shiota et al. 2014]).
We construct a “flux rope from eruption data” (FRED) by combining two key
results: (i) the reconnected (RC) flux during an eruption approximately equals
the poloidal flux of the ejected flux rope (Longcope et al. 2007; Qiu et al.
2007; Hu et al. 2014; Gopalswamy et al. 2017a), and (ii) white-light or EUV
coronal mass ejections (CMEs) can be fit to a FR to get its geometrical
properties (see e.g., Temmer et al. 2011). The RC flux is computed from the
area under post-eruption arcades (PEAs, McAllister and Martin 2000) and the
underlying unsigned photospheric magnetic field strength. The poloidal flux of
the FR is known from the RC flux; assuming that the FR is force free
(Lundquist 1951) we can get the axial and azimuthal field components and the
toroidal flux of the flux rope. Thus we have a fully-defined FR, complete with
geometric and magnetic parameters that can be used as input to global MHD
models. The axial and poloidal fields are essentially responsible for
geoeffectiveness in high and low-inclination magnetic clouds at Earth,
respectively and hence the FR input should lead to definite prediction
schemes. The FRED technique is complementary to another one that uses the
source-region relative magnetic helicity along with FR geometrical properties
to estimate the FR magnetic properties (Patsourakos et al. 2016).

Figure 1: (a) PEA area (polygon) marked on an SDO/AIA 193 Å image. (b) SDO/HMI
line-of-sight magnetogram with PEA and the filament magnetic field direction
(arrow) overlaid.

##  2 The Coronal Flux Rope

We illustrate the construction of FRED using the 2012 July 12 eruption that
resulted in a large SEP event and a major geomagnetic storm (Dst = -139 nT).
The event has been studied by many authors (see e.g., [Gopalswamy_etal13,
Gopalswamy et al. 2013]; [Gopalswamy_etal14, Gopalswamy et al. 2014];
[HessZhang14, Hess & Zhang 2014]; [Möstl_etal14, Möstl et al. 2014];
[Hu_etal16, Hu et al. 2016]). Most of these papers were concerned with the
kinematics and Sun-to-Earth propagation of the CME that ended up as a shock-
driving magnetic cloud at Earth. Many of the authors also fit a flux rope
using the graduated cylindrical shell (GCS) model ([Thernisien11, Thernisien
2011]). Here we use the GCS model to get the geometric parameters of the FR.
Geometrical Properties. Fitting the CME observed by the three views provided
by the SOHO and STEREO missions to the GCS model, we get the coordinates of
the flux rope as S12W06, which is slightly different from the flare location
in AR 11520 (S15W01, [Gopalswamy_etal14, Gopalswamy et al. 2014]). The ratio
of the FR radius ($R_{0}$) to the leading edge distance from the Sun center
($R_{tip}$) is 0.26. At $R_{tip}$ = 10 Rs, the $R_{0}$ = 2.6 Rs in the
coronagraph field of view. The face-on and edge-on angular widths of the FR
are 94∘ and 40∘, respectively. The tilt angle of the FR axis at its apex with
respect to the horizontal is about 53∘, indicating a northwest-southeast
orientation of the flux rope axis, consistent with the neutral line as
inferred from filament location.
Magnetic Properties. We combine these geometrical information with the RC flux
to get the magnetic properties of the FR. It was recently shown that the total
RC flux during an eruption can be obtained from a snapshot of the post
eruption arcade (PEA) and the underlying photospheric magnetic field strength
([Gopalswamy_etal17a, Gopalswamy et al. 2017a]). Fig. 1 shows a snapshot of
the PEA associated with the 2012 July 12 eruption as observed by the
Atmospheric Imaging Assembly (AIA, Lemen et al. 2012) on board the Solar
Dynamics Observatory (SDO, Pesnell et al. 2012) at 193 Å. Also shown is the
corresponding line-of-sight photospheric magnetic field from SDO’s
Helioseismic and Magnetic Imager (HMI, Scherrer et al. 2012). The arrow in
Fig. 1(b) points to the direction of the axial magnetic field of the erupted
flux rope based on filament connectivity. The axis is tilted by about 46
degrees to the horizontal. The area of the polygon in Fig. 1(a) is
7.2$\times$1019 cm2, and an average magnetic field strength of $B$ of 392 G.
The RC flux $\phi{{}_{r}}$ is thus 1.42$\times$1022 Mx, which is half the
unsigned flux through the PEA area.

Since $\phi{{}_{r}}$ roughly equals the poloidal flux ($\phi{{}_{p}}$) of the
erupted FR, we can get the axial field ($B_{0}$) of the coronal FR from the
relation, $\phi{{}_{p}}$ = $(L/x_{01})B_{0}R_{0}$, where $L$ is the FR length
and $x_{01}$ (=2.4048) is the first zero of the Bessel function $J_{0}$.
Taking $L$ = 2 $R_{tip}$ = 20 Rs, we get $B_{0}$ = $\phi{{}_{r}}x_{01}/LR_{0}$
= 0.13 G, which is near the higher end of the axial field of coronal FRs
reported in Gopalswamy et al. (2017b). The FR axial field strength is larger
than the typical ambient field strength. Gopalswamy and Yashiro (2011)
obtained the strength of the ambient magnetic field as $B_{amb}$ =
0.329$R$-1.23, which gives $B_{amb}$ = 0.019 G at a heliocentric distance $R$
of 10 Rs. Clearly, the FR axial field strength is about 7 times greater than
the ambient field, consistent with the fact that CME FRs are low-beta plasmas.
The toroidal flux of the FR is $\phi{{}_{t}}$ = $\phi{{}_{p}}$$(2\pi
R_{0}/L)J_{1}(x_{01})$, where $J_{1}$ is the first order Bessel function. With
the above numbers from the GCS fit, we get $\phi{{}_{t}}$ = 0.42$\phi{{}_{p}}$
or 5.96$\times$1021 Mx. The poloidal ($B_{p}$) and toroidal ($B_{t}$) field
strengths at any distance $r$ from the FR axis are given by $B_{p}$ =
$HB_{0}J_{1}(\alpha r)$ with $H$ = $\pm 1$ (helicity sign) and $B_{t}$ =
$B_{0}J_{0}(\alpha r)$. Here, $\alpha=x_{01}/R_{0}$ is the force-free
parameter. In the present case, $H$ = $+$1 because the axial field points
mostly to the south and the azimuthal field goes from the positive to the
negative side (see Fig. 1b). The direction is consistent with a southern
hemispheric eruption having right-handed helicity sign (e.g.,
[BothmerSchwenn98, Bothmer & Schwenn 1998]). Finally, the relative helicity
per unit length ($H_{r}/L$) is also well-defined for a Lundquist FR: $H_{r}/L$
= 0.7$B_{0}$2$R_{0}$3 (e.g., Dasso et al. 2003). For the 2012 July 12 coronal
FR, $H_{r}$ = 9.98 $\times$1043 Mx2, if we take the FR length to be
2$R_{tip}$. Thus, the flux rope is fully defined both geometrically and
magnetically. The three-dimensional CME speed was measured using STEREO-B as
$V$ = 1548 km s-1 ([Gopalswamy_etal13, Gopalswamy et al. 2013]). According to
the empirical relation $V$ = 298$\phi{{}_{r}}$0.75 ([Gopalswamy_etal17b,
Gopalswamy et al. 2017b]), this speed indicates a $\phi{{}_{r}}$ of
9.0$\times$1021 Mx, which is only about 36% lower than the observed value
(1.42$\times$1022 Mx).

Figure 2: Scatter plot between observed $B_{tot}$ and $\phi{{}_{r}}$ for a set
of events from solar cycle 23 reported in Gopalswamy et al. (2017b). Two
outliers are excluded from the correlation. Blue and red symbols correspond to
magnetic clouds (MC) and non-cloud ejecta (EJ), respectively.

##  3 Comparison with 1-AU Observations

A recent study showed that the axial field of the coronal FR derived from
$\phi_{r}$ correlates with the axial field of the 1-AU flux rope (fitted) and
also with the strength of the observed total magnetic field
($B_{tot}$)([Gopalswamy_etal17b, Gopalswamy et al. 2017b]). In that study the
white-light CME data were fit with the elliptical flux rope model of Krall and
St. Cyr (2006). Figure 2 shows a scatter plot between $\phi{{}_{r}}$ and
$B_{tot}$ for a set of interplanetary CME (ICME) events from cycle 23. From
the regression line $B_{tot}$ = 21.9$\log$$\phi{{}_{r}}$$-$453.2 nT, we get
$B_{tot}$ =31.9 nT with $\phi{{}_{r}}$ = 1.42$\times$1022 Mx for the 2012 July
12 eruption. The inferred $B_{tot}$ is very close to the observed value at 1
AU: $\sim$30 nT (see e.g., Hess and Zhang 2014). If the FR expands self-
similarly, the 1-AU FR radius is expected to be 0.26 AU (since
$R_{0}$/$R_{tip}$ = 0.26 with $R_{tip}$ = 1 AU). Therefore, one can get the
1-AU axial field strength from the coronal field strength ($B_{0}$ = 0.13 G)
as: $B_{01AU}$ = $B_{0}$(10/214)2 = 28.3 nT, again very similar to the
observed value. The assumption of self-similar expansion seems to be valid
because the FR radius is consistent with that obtained from Grad-Shafranov
reconstruction using in-situ data (Hu et al. 2016). For the same reasons,
$H_{r}$ in the 1 AU FR is similar to that in the coronal FR. The 1-AU data
also shows that the FR Y-component (azimuthal) of the magnetic field rotates
from west to east, while the Z-component (axial) points southward throughout
the cloud duration (Hess & Zhang 2014; Hu et al. 2016), consistent with the
orientation of the coronal FR. In this event, the $B_{z}$ is simply the axial
field component.

##  4 Discussion and Summary

We have shown that the FRED technique is a viable starting point to infer the
expected asymptotic magnetic structure in the heliosphere. It combines the
line-of-sight photospheric magnetograms and PEA observations (X-ray, EUV,
H-alpha, or microwave) with CME observations in the corona to obtain complete
properties of coronal FRs. We have shown these properties at a distance of 10
Rs for the purpose of illustration, but FRED properties can be obtained at any
initial height where the CME is observed. For example, the inner boundary of
many MHD simulation models is located at $R_{tip}$ $\sim$ 21 Rs, where the
axial field $B_{0}$ = $\phi{{}_{r}}x_{01}/LR_{0}$ will be weaker than the
10-Rs case because both $L$ and $R_{0}$ will be larger. Expressing $L$ and
$R_{0}$ in terms of $R_{tip}$ and noting that $R_{0}$/$R_{tip}$ is invariant
for each FR from the GCS model we get $B_{0}$ =
$\phi{{}_{r}}x_{01}/[pq(R_{tip}$)2], where $p$ and $q$ are constants: $L$ =
$pR_{tip}$ (assumed FR length) and $R_{0}$ = $qR_{tip}$ (from the GCS model).
This expression shows that the axial field falls off as the square of the
heliocentric distance. Note that we have taken $p$ = 2 following Nindos et al.
(2003). The length can be higher if the flux rope assumes the shape of a half
torus, so $L$ = $\pi$$R_{tip}$. Démoulin et al. (2016) recently reported a
statistical value of $p$ = 2.6$\pm 0.3$, which is in between the two values
noted above. The invariant toroidal flux can also be written in a general form
as $\phi{{}_{t}}$ = $\phi{{}_{p}}$$(2\pi q/p)J_{1}(x_{01})$. Global MHD models
that use a realistic input such as FRED should be able to provide a realistic
forecast of what to expect in the heliosphere. It must be noted that
additional effects such as CME deflection and rotation can modify the FR,
which can be accounted for using semi-analytic models (see e.g., Kay et al.
2015). The results presented in this paper support the idea that CME FRs are
formed during the eruption. The poloidal flux of any pre-existing FR is
expected to be a small fraction of that added during the eruption.

Acknowledgement. Work supported by NASA Heliophysics GI and LWS programs.

## References

  * [Bothmer & Schwenn(1998)] Bothmer, V., & Schwenn, R. 1998, Annales Geophysicae, 16, 1
  * [Dasso et al.(2003)] Dasso, S., Mandrini, C. H., DéMoulin, P., & Farrugia, C. J. 2003, JGR, 108, 1362
  * [Démoulin et al.(2016)] Démoulin, P., Janvier, M., & Dasso, S. 2016, SolPhys, 291, 531
  * [Gopalswamy et al.(2008)] Gopalswamy, N., Akiyama, S., Yashiro, S., Michalek, G., & Lepping, R. 2008, JASTP, 70, 245
  * [Gopalswamy & Yashiro(2011)] Gopalswamy, N., & Yashiro, S. 2011, ApJ, 736, L17
  * [Gopalswamy et al.(2013)] Gopalswamy, N., Mäkelä, P., Xie, H., & Yashiro, S. 2013, Space Weather, 11, 661
  * [Gopalswamy et al.(2014)] Gopalswamy, N., Xie, H., Akiyama, S., Mäkelä, P. A., & Yashiro, S. 2014, Earth, Planets, and Space, 66, 104
  * [Gopalswamy et al.(2017a)] Gopalswamy, N., Yashiro, S., Akiyama, S., & Xie, H. 2017a, SolPhys, 292, 65
  * [Gopalswamy et al.(2017b)] Gopalswamy, N., Akiyama, S., Yashiro, S., & Xie, H. 2017b, JASTP, arXiv:1705.08912
  * [Hess & Zhang(2014)] Hess, P., & Zhang, J. 2014, ApJ, 792, 49
  * [Hu et al.(2014)] Hu, Q., Qiu, J., Dasgupta, B., Khare, A., & Webb, G. M. 2014, ApJ, 793, 53
  * [Hu et al.(2016)] Hu, H., Liu, Y. D., Wang, R., Möstl, C., & Yang, Z. 2016, ApJ, 829, 97
  * [Jin et al.(2017)] Jin, M., Manchester, W. B., van der Holst, B., et al. 2017, ApJ, 834, 173
  * [Kay et al.(2015)] Kay, C., Opher, M., & Evans, R. M. 2015, ApJ, 805, 168
  * [Krall & St. Cyr(2006)] Krall, J., & St. Cyr, O. C. 2006, ApJ, 652, 1740
  * [Lemen et al.(2012)] Lemen, J. R., Title, A. M., Akin, D. J., et al. 2012, SolPhys, 275, 17
  * [Longcope et al.(2007)] Longcope, D., Beveridge, C., Qiu, J., et al. 2007, SolPhys, 244, 45
  * [Lundquist (1951)] Lundquist, S. 1951, Phys. Rev., 83, 307
  * [McAllister & Martin (2000)] McAllister, H., Martin, S. F. 2000, AdSpR, 26, 469
  * [Möstl et al.(2014)] Möstl, C., Amla, K., Hall, J. R., et al. 2014, ApJ, 787, 119
  * [Nindos et al.(2003)] Nindos, A., Zhang, J., & Zhang, H. 2003, ApJ, 594, 1033
  * [Odstrčil & Pizzo(1999)] Odstrčil, D., & Pizzo, V. J. 1999, JGR, 104, 483
  * [Patsourakos et al.(2016)] Patsourakos, S., Georgoulis, M. K., Vourlidas, A., et al. 2016, ApJ, 817, 14
  * [Pesnell et al.(2012)] Pesnell, W. D., Thompson, B. J., & Chamberlin, P. C. 2012, SolPhys, 275, 3
  * [Qiu et al.(2007)] Qiu, J., Hu, Q., Howard, T. A., & Yurchyshyn, V. B. 2007,ApJ, 659, 758
  * [Scherrer et al.(2012)] Scherrer, P. H., Schou, J., Bush, R. I., et al. 2012, SolPhys, 275, 207
  * [Shen et al.(2014)] Shen, F., Shen, C., Zhang, J., et al. 2014, JGR, 119, 7128
  * [Shiota et al.(2014)] Shiota, D., Kataoka, R., Miyoshi, Y., et al. 2014, Space Weather, 12, 187
  * [Temmer et al.(2011)] Temmer, M., Veronig, A. M., Gopalswamy, N., & Yashiro, S. 2011, SolPhys, 273, 421
  * [Thernisien(2011)] Thernisien, A. 2011, ApJS, 194, 33
  * [Wu & Lepping(2002)] Wu, C.-C., & Lepping, R. P. 2002, JGR, 107, 1346