#  Edge Detection Quantumized: A Novel Quantum Algorithm for Image Processing

Syed Emad Uddin Shubha  Department of Electrical and Computer Engineering,
North South University, Bashundhara, Dhaka  Mir Muzahedul Islam  Department of
Electrical and Computer Engineering, North South University, Bashundhara,
Dhaka  Tanvir Ahahmed Sadi  Department of Electrical and Computer Engineering,
North South University, Bashundhara, Dhaka  Md. Hasibul Hasan Miraz
Department of Electrical and Computer Engineering, North South University,
Bashundhara, Dhaka  M.R.C. Mahdy  Department of Electrical and Computer
Engineering, North South University, Bashundhara, Dhaka

###### Abstract

Quantum image processing is a research field that explores the use of quantum
computing and algorithms for image processing tasks such as image encoding and
edge detection. Although classical edge detection algorithms perform
reasonably well and are quite efficient, they become outright slower when it
comes to large datasets with high-resolution images. Quantum computing
promises to deliver a significant performance boost and breakthroughs in
various sectors. Quantum Hadamard Edge Detection (QHED) algorithm, for
example, works at constant time complexity, and thus detects edges much faster
than any classical algorithm. However, the original QHED algorithm is designed
for Quantum Probability Image Encoding (QPIE) and mainly works for binary
images. This paper presents a novel protocol by combining the Flexible
Representation of Quantum Images (FRQI) encoding and a modified QHED
algorithm. An improved edge outline method has been proposed in this work
resulting in a better object outline output and more accurate edge detection
than the traditional QHED algorithm.

_Keywords_ Quantum Image Representation $\cdot$ Quantum Algorithm $\cdot$
Image Processing, $\cdot$ Edge Detection $\cdot$ Hadamard $\cdot$ Amplitude
Encoding $\cdot$ Angle Encoding $\cdot$ Unitary $\cdot$ Partial Measurement.

##  1 Introduction

††footnotetext: * Corresponding author. E-mail address:
<EMAIL_ADDRESS>(M.R.C. Mahdy).

In modern technology, image processing is an indispensable tool, employing
sophisticated algorithms to extract valuable insights from visual data. Its
applications are vast and diverse, spanning fields such as image segmentation,
enhancement, computer vision, recognition, and medical imaging. With the
emergence of quantum image processing algorithms, significant strides have
been made in critical operations like edge detection and beyond [1, 2].

Quantum computers have the prospect of revolutionizing many fields by solving
problems that are currently computationally expensive. The development of
quantum computers is a rapidly advancing field of research and potential
applications of quantum computers are currently being investigated. Quantum
Image Processing shows a promising aptness by utilizing quantum superposition
and entanglement as a powerful resource [3, 4, 5]. Converting image data into
a quantum state, we can perform a wide range of linear operations, [6, 7, 8],
thus transforming and processing images for various applications. There are
many Quantum Image Representation techniques, such as Flexible Representation
of Quantum Images (FRQI), Quantum Probability Image Encoding (QPIE), and Novel
Enhanced Quantum Representation (NEQR), offering efficient pathways for
translating pixel information into quantum data. Among these two methods, the
FRQI method uses angle encoding, which requires fewer qubits [9, 10]. The QPIE
method requires one qubit less than FRQI. The NEQR method, on the other hand,
requires a much higher number of qubits. [11, 12].

Edge detection, a cornerstone of image analysis and feature extraction,
assumes paramount importance across a multitude of applications, ranging from
fingerprint recognition to medical imaging. While classical edge detection
algorithms, such as, Sobel, Laplacian, Prewitt, and Canny, have proven
effective in many scenarios, they encounter significant scalability challenges
when dealing with larger images, owing to the increased computational demands
[13, 14, 15, 16].

The emergence of quantum algorithms for edge detection presents a promising
avenue for overcoming these challenges. Leveraging the inherent parallelism of
quantum computation, these algorithms enable the simultaneous processing of
multiple pixels, thereby reducing the overall time complexity. Among these
pioneering algorithms, Quantum Hadamard Edge Detection (QHED) stands out for
its constant time complexity and memory efficiency [17, 18]. However, despite
its inherent advantages, QHED’s effectiveness diminishes when confronted with
real-world images, exhibiting inaccuracies and noise due to its reliance on
binary image bases.

In response to these challenges, this paper proposes a series of innovative
modifications to enhance QHED’s accuracy and efficiency. At first, we used the
angle encoding (FRQI) method to convert a classical image into a quantum
state, followed by a partial measurement. After that, we reset the existing
qubit and use it as a resource to implement traditional QHED with an
additional qubit method. Then a classical post processing has been applied to
generate an edge-detected output image. We have also analyzed the procedure
from an information-theoretic perspective. We have designed the protocol to
encompass the development of a versatile quantum circuit adaptable to diverse
image inputs, refinement of object outline methodologies to improve precision
in representation, and significant enhancements in edge detection for complex
images through a combination of classical and quantum modifications applied to
existing QHED algorithms. Through these endeavors, we seek to push the
boundaries of quantum image processing, paving the way for more effective and
robust edge detection methodologies in both theory and practice. Fig. 1 shows
our proposed protocol and Fig. 2 gives a comprehensive overview of our
improved protocol over traditional QHED.

Figure 1: System Diagram of the Proposed Model Figure 2: Comparison between
Traditional QHED and Proposed Model.

##  2 Preliminary Concepts

###  2.1 Quantum Image Processing

An image is given by a matrix $F=(F_{ij})_{M\times N}$. If the image is a
grey-scale image, then $0\leq M_{kl}\leq 1$ represents the pixel value of
$j^{th}$ row and $k^{th}$ column. If the image is a color image, then we have
$M_{lk}=[R_{lk}G_{lk}B_{lk}]$, representing the RGB values of a pixel,
typically these values are 8 bits each, thus having a value between $0$ to
$255$. Some images have an additional component, named channel. Classically
image processing can be given by the Map $G=PFQ$, where $G=(G_{ij})_{M\times
N}$.
In quantum image processing, We first need to encode the image into a quantum
state called Quantum Image Representation. The state is represented by a
normalized column vector $|f\rangle$ that encodes color and position
information. There are several ways to encode an image, for example, in Novel
Enhanced Quantum Representation (NEQR), a $2^{n}\times 2^{n}$ grey scale image
is represented as [19, 20]

$|I\rangle=\frac{1}{2^{n}}\sum_{Y=0}^{2^{n}-1}\sum_{X=0}^{2^{n}-1}|{f(Y,X)}\rangle|{YX}\rangle$
(1)

Here, $\Omega_{YX}=\otimes_{i=0}^{i=q-1}|C_{YX}^{i}\rangle$ encodes pixel
information, where $C_{YX}^{i}\in{0,1}$ and $q=8$ for grey scale image.
Now quantum image processing can be described by linear quantum operations,
given by:

$\epsilon(|f\rangle)=\sum_{i}E_{i}|f\rangle\langle f|E_{i}^{\dagger}$ (2)

Here $\sum_{k}E_{k}^{\dagger}E_{k}=\mathds{1}$. We call $\\{E_{j}\\}$ Kraus
Operators, and the map $\epsilon$ represents a completely positive trace-
preserving map, and $\epsilon(|f\rangle)$ is called a density operator.

###  2.2 Flexible Representation of Quantum Images

The Flexible Representation of Quantum Images encodes the pixel data in
angles, which requires only one qubit to store pixel information, rest of the
qubits stores position information. Given
$\\{\theta_{0},\theta_{1},...,\theta_{4^{n}-1}\\}\quad(\theta_{i}\in[0,\pi/2])$,
the FRQI state is given by,

$|I(\theta)\rangle=\frac{1}{2^{n}}\sum_{i=0}^{4^{n}-1}(cos\theta_{i}|{0}\rangle+sin\theta_{i}|{1}\rangle)\otimes|{i}\rangle$
(3)

We can transform $|0\rangle^{\otimes 2n+1}$ into FRQI state,
$|I(\theta)\rangle$ by applying unitary operator $\mathcal{R}\mathcal{H}$,
with $\mathcal{H}=I\otimes H^{\otimes 2n}$ and
$\mathcal{R}=\prod_{i=0}^{4^{n}-1}R_{i}$. where $R_{i}$ is a unitary matrix
named controlled rotation, given by

$R_{i}=I\otimes\sum_{j=0,j\neq
i}^{4^{n}-1}|{j}\rangle\langle{j}|+R_{y}(2\theta_{i})\otimes|{i}\rangle\langle{i}|=I^{\otimes
2n+1}+(R_{y}(2\theta_{i})-I)\otimes|{i}\rangle\langle{i}\mid$ (4)

We can implement the controlled rotation $R_{i}$’s using $C^{2n}R_{y}$
(controlled $R_{y}$) and $X$ gate. Here,

$C^{2n}R_{y}(2\theta_{i})=I^{\otimes
2n+1}+(R_{y}(2\theta_{i})-I)\otimes|{4^{n}-1}\rangle\langle{4^{n}-1}\mid$ (5)

We can convert $C^{2n}R_{y}$ into $R_{i}$ by applying NOT gates on the left
and right sides. For example, in case of $2\times 2$ grey image, we see,
$R_{3}=C^{2}R_{y}$, and

$\begin{split}R_{2}&=(I\otimes I\otimes X)C^{2}R_{y}(I\otimes I\otimes X)\\\
R_{1}&=(I\otimes X\otimes I)C^{2}R_{y}(I\otimes X\otimes I)\\\
R_{0}&=(I\otimes X\otimes X)C^{2}R_{y}(I\otimes X\otimes X)\\\ \end{split}$
(6)

Therefore controlled rotations $R_{i}$ can be implemented by
$C^{2n}R_{y}(2\theta_{i})$ and NOT operations.

Basically, applying $I\otimes P_{i}$ in the both side of $C^{2n}R_{y}$ to make
the conversion:

$P_{i}|{4^{n}-1}\rangle=|{i}\rangle$ (7)

Here $P_{i}$ is the tensor product of some $I$ and $X$, which is a hermitian.
Fig. 1 refers to the circuit of the FRQI method for an arbitrary $2\times 2$
color image.
Say, one colored pixel has the RGB values $(R_{i},G_{i},B_{i})$, then we
encode the information as follows:

$\theta_{i}=cos^{-1}(\frac{R_{i}}{256}+\frac{G_{i}}{256^{2}}+\frac{B_{i}}{256^{3}})$
(8)

The advantage is, given $\theta_{i}$, it is easier to get correspondent
$R_{i},G_{i},B_{i}$ by expressing $256^{3}cos(\theta_{i})$ as $256$ base
number. Fig. 3 refers to a circuit implementation of a $2\times 2$ image.

Figure 3: FRQI Encoding circuit for a $2\times 2$ sample color image.

###  2.3 Quantum Hadamard Edge Detection

QHED edge detection algorithm is based on QPIE image representation. QPIE will
convert digital image information into a quantum state. To begin with, the
images are represented as normalized probability amplitude corresponding to a
specific quantum state. Given an image $I_{xy}$, it can be expressed as

$|img\rangle=\sum_{x,y}\frac{I_{yx}}{\sqrt{\sum
I_{yx}^{2}}}|xy\rangle=\begin{pmatrix}c_{0}&c_{1}&\ldots&c_{N-1}\end{pmatrix}^{\mathrm{T}}$
(9)

To build the QHED circuit an additional qubit is required. This will help to
calculate all the pixel differences at once. At first, the image information
is initialized to the N number of data qubits using the QPIE technique. Then,
the Hadamard gate is applied to the auxiliary qubit. Finally, we have the
following state that encodes the image:

$\left|img^{\prime}\right\rangle=\frac{1}{\sqrt{2}}\begin{pmatrix}c_{0}&c_{0}&c_{1}&c_{1}&\ldots&c_{N-1}&c_{N-1}\end{pmatrix}^{\mathrm{T}}$
(10)

After that, we performed an amplitude permutation unitary operation
$D_{2^{n+1}}$ to shift down the element of the state vector, by one row. This
operation will simultaneously facilitate the gradient calculation for both
even and odd pixel pairs. The expression of $D_{2^{n+1}}$ and its operation on
the image is given in Eqn. (11) and (12) respectively.

$D_{2^{n+1}}=\left[\begin{array}[]{ccccc}0&1&0&\cdots&0\\\ 0&0&1&\cdots&0\\\
\vdots&\vdots&\vdots&\ddots&\vdots\\\ 0&0&0&0&1\\\
1&0&0&0&0\end{array}\right]$ (11)
$D_{2^{n+1}}\left|img^{\prime}\right\rangle=\frac{1}{\sqrt{2}}\begin{pmatrix}c_{0}&c_{1}&c_{1}&\ldots&c_{N-1}&c_{N-1}&c_{0}\end{pmatrix}^{\mathrm{T}}$
(12)

After the above operation, an H-gate is applied to the auxiliary qubit. The
resultant state would be,

$\left|img^{\prime\prime}\right\rangle=\frac{1}{2}\begin{pmatrix}c_{0}+c_{1}&c_{0}-c_{1}&c_{1}+c_{2}&c_{1}-c_{2}\ldots
c_{N-1}+c_{0}&c_{N-1}-c_{0}\end{pmatrix}^{\mathrm{T}}$ (13)

The quantum circuit of the QHED algorithm for a $4\times 4$ binary image is
shown in Fig. 4.

Figure 4: Quantum Circuit of QHED Algorithm.

Edge detection is the process of detecting the edge of an object from a
picture based on dissimilarity of color intensity. The discontinuation can be
measured from the difference between two neighbor pixels. If we measured the
auxiliary qubit being in state |1⟩ then the necessary information for edge
detection can be retrieved. The whole procedure carries out a horizontal scan
to detect edges in the horizontal direction. For vertical scan, the image is
transposed and the above operations are repeated. Finally, the horizontal and
vertical scan outputs are superimposed into one another using classical post-
processing for the final edge detected output.

##  3 Proposed Edge Detection Algorithm

In this paper, we have modified the edge outline of an object by analyzing the
relative edge detection order and imposing a hyperparameter to filter out the
noise. We have also modified the encoding part using FRQI and partial
measurement. In this section, let’s assume we have an image with size
$2^{n}times2^{n}$.

###  3.1 Proposed Encoding Method

Suppose we have two grey-scale images, Image$-1=\\{I_{i}\\}$ and
Image$-2=\\{I^{\prime}_{i}\\}$, such that $I^{\prime}_{k}=\sqrt{1-I_{k}^{2}}$
($\forall{k}$).
Since $0\leq I_{k}\leq 1$, this map is a bijection with smooth behavior, that
preserves the edge information, as shown in Fig. 5:

Figure 5: (a) image $-1$, (b) image $-2$. We note that the map doesn’t change
the edges.

Now, say, $I_{i}=cos(\theta_{i})$. Therefore the QPIE representations of
Images 1 and 2 are given by,

$|img_{1}\rangle=\sum_{i=0}^{4^{n}-1}\frac{I_{i}}{\sqrt{\sum_{j=0}^{4^{n}-1}I_{j}^{2}}}|i\rangle=\sum_{i=0}^{4^{n}-1}\frac{cos(\theta_{i})}{\sqrt{\sum_{j}{cos^{2}(\theta_{j}})}}|i\rangle$
(14)

$|\text{img}_{2}\rangle=\sum_{i=0}^{4^{n}-1}\frac{I_{i}^{\prime}}{\sqrt{\sum_{j=0}^{4^{n}-1}I_{j}^{\prime
2}}}|i\rangle=\sum_{i=0}^{4^{n}-1}\frac{\sin(\theta_{i})}{\sqrt{\sum_{j}\sin^{2}(\theta_{j})}}|i\rangle$
(15)

Now in our proposed protocol, we first encode our image using FRQI, given by
Eq. (3). Say, the pixel information is stored at ancilla qubit, which will now
be measured in the Z-basis, i.e., the Krauss Operators are given by,
$\\{|0\rangle\langle 0|\otimes\mathds{1}^{\otimes 2n},|1\rangle\langle
1|\otimes\mathds{1}^{\otimes 2n}\\}$. The collapsed state will be given by
either Eq. (14) or Eq. (15), depending on the measurement outcome. Since both
collapsed state resembles the same edge, we will now apply the QHED algorithm
with the additional qubit. The modified algorithm is shown in the diagram
below:

Figure 6: System Diagram of Modified Algorithm.

###  3.2 Edge Detection Improvement

The ancilla qubit will be measured at state $|0\rangle$ or $|1\rangle$, we
will put an $X$ gate if the first case happens. After that a Hadamard gate
will be applied to the ancilla qubit, that is how the state of the ancilla
qubit will now be in $|-\rangle$ to apply conventional QHED. QHED process
compares two neighbor pixels, and based on their dissimilarity, it plots the
edge. The edge detection operation occurs based on the superimposed scan
results of left-to-right and top-to-bottom scans. One of the issues with QHED
edge detection is that, it does not consider object boundaries. If there is a
difference in color intensity between the current and the next pixel, QHED
portrays an edge in the current pixel. QHED defines the edges inside an object
boundary and does not wrap around the object. As a consequence, the resultant
object outline is ambiguous and does not represent the actual object and shape
correctly. We considered the object’s exact position to address this
complication and analyzed the pixel gradient result. The edge outline around
the object borders are then portrayed according to these findings.

Comparison of edge detection between traditional approach (middle column) and
modified approach (far right column) for input image (far left column). It is
evident the conventional approach does not capture the object outline
correctly, and the actual object retrieval might not be easy. In comparison,
our proposed model outlines the object accurately. The main QHED algorithm
cannot detect all the edges properly for complex images and produces a lot of
noise, while our modified approach can reduce noise using a dynamic threshold.
In Fig. 7,8 & 9, we have demonstrated how our edge detection algorithm
provides better visualization.

Figure 7: Scanning result difference between traditional QHED and the proposed
method. Figure 8: Comparison of edge detection between traditional approach
(middle column) and modified approach (far right column) for input image (far
left column). It is evident the traditional approach does not capture the
object outline correctly and the actual object retrieval might not be easy.
Whereas our proposed model outlines the object accurately. Figure 9:
Comparison of edge detection output. The middle column represents the
conventional QHED implementation result, and the improved implementation of
QHED is shown in the far-right column.

Threshold has been used to remove noise in grey-scale images and improve the
edge outline. For any $2^{n}\times 2^{n}$ image, if the represented state is
given by $|\psi\rangle=\sum_{i}c_{i}|i\rangle$, we define the threshold as:

$thr=\frac{|max_{i}\\{c_{i}-c_{i+1}\\}|}{2n}$ (16)

In the post-processing, we will generate a binary image, if the difference is
greater than the threshold of the main picture, we will call it an edge. If
the sign of difference is mismatched with the first edge, the edge will be
shifted one pizel either in the vertical or horizontal direction based on
scanning order. Lastly, noises will be removed if the pixel difference is less
than $(-thr)$. Thus the post-processing gives us an improved edge outline
without noise.

###  3.3 Quantum Thermodynamics Perspective

When we partially measure the state given by Eq. (3), we have either collapsed
the state in Eq.(14) or Eq. (15). We note that, even if we were to use Eq. (8)
to preserve the color information, the measurement will make the state into
the representation of a grey-scale image. We can easily note that by observing
the fact that, for $i=1,2,..,n$, the set $\\{cos(\theta_{i})\\}$ is
independent but $\\{\frac{cos(\theta_{i})}{\sum_{i}{cos^{2}(\theta_{i}}}\\}$
is not. Hence, it can be said that grey-scale representation loses some
information, we can produce the state by partially measuring the color-encoded
state, i.e., measuring the pixel qubit. The Edge detection overall means
losing some information but keeping the edge outline invariant. This also will
be examined thoroughly in the next version of this work.

##  4 Conclusion

In this paper, we provided a generalized implementation scheme for NEQR and
FRQI algorithms. We also have addressed the shortcomings of the QHED algorithm
and proposed a modified version of QHED. Traditional QHED does not work for
all input images, and it struggles with complex images. Our modified algorithm
allowed the encoding and edge-detecting circuits to work with various types of
complex images. We implemented a dynamic circuit for the input image, fixed
the edge outline, and, most importantly, adjusted the threshold value to
remove unnecessary noise. Overall, our project has significantly improved the
QHED algorithm, making it a more efficient and accurate method for edge
detection in quantum image processing. The changes we have made will allow the
QHED algorithm to be applied to a broader range of images and provide more
accurate results. We have also used the FRQI method, along with a partial
measurement, to replace the QPIE method. This inclusion makes our model a more
novel and physics-oriented approach.

## Acknowledgment

The authors would like to thank Research Assistants Deponker Sarkar Depto and
Md Shawmoon from the NSU Optics Lab for their constructive suggestions and
technical help, respectively. M.R.C. Mahdy acknowledges the support of NSU
CTRGC grant 2023-24 (approved by the members of BOT) and NSU internal grant.

## References

  * [1] Junfeng Jing, Shenjuan Liu, Gang Wang, Weichuan Zhang, and Changming Sun.  Recent advances on image edge detection: A comprehensive review.  Neurocomputing, 503:259–271, 2022.
  * [2] Zhaobin Wang, Minzhe Xu, and Yaonan Zhang.  Review of quantum image processing.  Archives of Computational Methods in Engineering, 29(2):737–761, 2022.
  * [3] Quantum image processing on real superconducting and trapped-ion based quantum computers.  Tm-technisches Messen, 2023.
  * [4] Yongquan Cai, Xiaowei Lu, and Nan Jiang.  A survey on quantum image processing.  Chinese Journal of Electronics, 27(4):718–727, 2018.
  * [5] Andrei N Soklakov and Rüdiger Schack.  Efficient state preparation for a register of quantum bits.  Physical review A, 73(1):012307, 2006.
  * [6] Novel Design of Quantum Circuits for Representation of Grayscale Images.  2023\.
  * [7] An improved order-encoded quantum image representation model and its application.  International Journal of Quantum Information, 2023.
  * [8] Yue Ruan, Hanwu Chen, Jianing Tan, and Xi Li.  Quantum computation for large-scale image classification.  Quantum Information Processing, 15:4049–4069, 2016.
  * [9] Quang Phuc Le, Dong Fangyang, Arai Yoshinori, and Hirota Kaoru.  Flexible representation of quantum images and its computational complexity analysis.  In Proceedings of the 25th Fuzzy Systems Symposium of the Japanese Society for Fuzzy Theory and Intelligent Information, pages 185–185. Japanese Society for Fuzzy Theory and Intelligent Information, 2009.
  * [10] Phuc Q Le, Fangyan Dong, and Kaoru Hirota.  A flexible representation of quantum images for polynomial preparation, image compression, and processing operations.  Quantum Information Processing, 10:63–84, 2011.
  * [11] Quantum image processing algorithm using line detection mask based on neqr.  Entropy, 2023.
  * [12] An evolutionary quantum scrambling scheme for medical image with neqr representation.  International Journal of Computing and Digital Systems, 2023.
  * [13] Quantum-inspired edge detection algorithms implementation using new dynamic visual data representation and short-length convolution computation.  arXiv.org, 2022.
  * [14] Rajib Chetia, S. M. B. Boruah, and P. P. Sahu.  Quantum image edge detection using improved sobel mask based on neqr.  Quantum Information Processing, 20:1–25, 2021.
  * [15] Pengao Xu, Zhenxing He, Tianhui Qiu, and Hongyang Ma.  Quantum image processing algorithm using edge extraction based on kirsch operator.  Optics express, 28 9:12508–12517, 2020.
  * [16] Ri gui Zhou, Han Yu, Yu Cheng, and Feng-Xin Li.  Quantum image edge extraction based on improved Prewitt operator.  Quantum Information Processing, 18, 2019.
  * [17] Analisis deteksi tepi citra dengan quantum hadamard edge detection (qhed).  Techno.COM Jurnal, 2022.
  * [18] Xi-Wei Yao, Hengyan Wang, Zeyang Liao, Ming-Cheng Chen, Jian Pan, Jun Li, Kechao Zhang, Xingcheng Lin, Zhehui Wang, Zhihuang Luo, et al.  Quantum image processing and its application to edge detection: theory and experiment.  Physical Review X, 7(3):031041, 2017.
  * [19] Yi Zhang, Kai Lu, Yinghui Gao, and Mo Wang.  Neqr: a novel enhanced quantum representation of digital images.  Quantum information processing, 12:2833–2860, 2013.
  * [20] Realization of quantum mean filters with different sized on neqr quantum images by qft based operations.  International Journal of Information Security Science, 2023.


* [15] L. Coutin and Z. Qian.  Stochastic analysis, rough path analysis and fractional brownian motions.  Probability theory and related fields, 122(1):108–140, 2002.
  * [16] T. Fawcett.  Problems in stochastic analysis: connections between rough paths and non-commutative harmonic analysis.  PhD thesis, University of Oxford, 2003.
  * [17] E. R. Ferrucci.  Rough path perspectives on the Itô-Stratonovich dilemma.  PhD thesis, Imperial College London, 2022.  https://spiral.imperial.ac.uk/handle/10044/1/96036.
  * [18] P. Friz and S. Riedel.  Convergence rates for the full Gaussian rough paths.  Annales de l’Institut Henri Poincaré, Probabilités et Statistiques, 50(1):154 – 194, 2014.
  * [19] P. Friz and N. Victoir.  Differential equations driven by Gaussian signals.  Annales de l’Institut Henri Poincaré, Probabilités et Statistiques, 46(2):369 – 413, 2010.
  * [20] P. K. Friz and A. Shekhar.  General rough integration, Lévy rough paths and a Lévy–Kintchine-type formula.  The Annals of Probability, 45(4):2707 – 2765, 2017.
  * [21] P. K. Friz and N. B. Victoir.  Multidimensional stochastic processes as rough paths, volume 120 of Cambridge Studies in Advanced Mathematics.  Cambridge University Press, Cambridge, 2010.  Theory and applications.
  * [22] B. Hambly and T. Lyons.  Uniqueness for the signature of a path of bounded variation and the reduced path group.  Annals of mathematics, 171(1):109–167, 2010.
  * [23] S. T. Huang and S. Cambanis.  Stochastic and multiple wiener integrals for gaussian processes.  The Annals of probability, 6(4):585–614, 1978.
  * [24] M. Jolis.  On the Wiener integral with respect to the fractional Brownian motion on an interval.  J. Math. Anal. Appl., 330(2):1115–1127, 2007.
  * [25] P. E. Kloeden and E. Platen.  Numerical solution of stochastic differential equations, volume 23 of Applications of Mathematics (New York).  Springer-Verlag, Berlin, 1992.
  * [26] P. Lévy.  Random Functions: General Theory with Special Reference to Laplacian Random Functions.  University of California publications in statistics. University of California Press, 1953.
  * [27] S. Li and H. Ni.  Expected signature of stopped brownian motion on $d$-dimensional ${C}^{2,\alpha}$-domains has finite radius of convergence everywhere: $2\leq d\leq 8$.  Journal of Functional Analysis, 282(12):109447, 2022.
  * [28] T. Lyons, M. Caruana, and T. Lévy.  Differential equations driven by rough paths ecole d’eté de probabilités de saint-flour xxxiv-2004, 2007.
  * [29] T. Lyons and H. Ni.  Expected signature of Brownian motion up to the first exit time from a bounded domain.  The Annals of Probability, 43(5):2729 – 2762, 2015.
  * [30] T. Lyons and Z. Qian.  System control and rough paths.  Oxford mathematical monographs. Clarendon Press, Oxford, 2002.
  * [31] T. Lyons and N. Victoir.  Cubature on Wiener space.  volume 460, pages 169–198. 2004.  Stochastic analysis with applications to mathematical finance.
  * [32] B. B. Mandelbrot and J. W. V. Ness.  Fractional brownian motions, fractional noises and applications.  SIAM review, 10(4):422–437, 1968.
  * [33] D. Marinucci and P. M. Robinson.  Alternative forms of fractional brownian motion.  Journal of Statistical Planning and Inference, 80(1):111 – 122, 1999.
  * [34] I. Nourdin and G. Peccati.  Normal approximations with Malliavin calculus: from Stein’s method to universality.  Cambridge tracts in mathematics; 192. Cambridge University Press, Cambridge [England]; New York, 2012.
  * [35] D. Nualart.  The Malliavin Calculus and Related Topics.  Probability and Its Applications. 2nd ed. edition, 2006.
  * [36] D. Nualart and S. Tindel.  A construction of the rough path above fractional Brownian motion using Volterra’s representation.  The Annals of Probability, 39(3):1061 – 1096, 2011.
  * [37] R. Passeggeri.  On the signature and cubature of the fractional brownian motion for h¿12.  Stochastic Processes and their Applications, 130(3):1226–1257, 2020\.


# The current state and future directions of modeling thermosphere density
enhancements during extreme magnetic storms

Denny M. Oliveira 1,2∗, Eftyhia Zesta 2, Piyush M. Mehta 3, Richard J. Licata
3, Marcin D. Pilinski 4, W. Kent Tobiska 5, and Hisashi Hayakawa 6,7,8,9

###### Abstract

## 1

Satellites, crewed spacecraft and stations in low-Earth orbit (LEO) are very
sensitive to atmospheric drag. A satellite’s lifetime and orbital tracking
become increasingly inaccurate or uncertain during magnetic storms. Given the
planned increase of government and private satellite presence in LEO, the need
for accurate density predictions for collision avoidance and lifetime
optimization, particularly during extreme events, has become an urgent matter
and requires comprehensive international collaboration. Additionally, long-
term solar activity models and historical data suggest that solar activity
will significantly increase in the following years and decades. In this
article, we briefly summarize the main achievements in the research of
thermosphere response to extreme magnetic storms occurring particularly after
the launching of many satellites with state-of-the-art accelerometers from
which high-accuracy density can be determined. We find that the performance of
an empirical model with data assimilation is higher than its performance
without data assimilation during all extreme storm phases. We discuss how
forecasting models can be improved by looking into two directions: first, to
the past, by adapting historical extreme storm datasets for density
predictions, and second, to the future, by facilitating the assimilation of
large-scale thermosphere data sets that will be collected in future events.
Therefore, this topic is relevant to the scientific community, government
agencies that operate satellites, and the private sector with assets operating
in LEO.

##  2 Introduction

During magnetic storms, large amounts of magnetospheric energy enters the
ionosphere-thermosphere system at high latitudes through field-aligned
currents (Prölss, 2011; Emmert, 2015). Since ionospheric currents are
intensified, the interaction of the moving plasma with the local neutral gas
is enhanced through the collision between ions and neutral molecules, thereby
further heating the neutral atmosphere. The atmosphere then expands upwards
and satellites flying in higher regions experience increased atmospheric drag
forces, which in turn intensify orbital drag effects. These effects ultimately
decrease the satellite’s life time and introduce significant errors in orbital
tracking which increase as the storm becomes more intense (Doornbos and
Klinkrad, 2006; Zesta and Huang, 2016; Nwankwo et al., 2021).

The third satellite in the Sputnik series provided arguably the first
observations of storm-induced orbital drag effects. Jacchia (1959) observed a
strong decay of Sputnik 1958$\delta$1 due to increased atmospheric drag
forces. These pioneer observations led to the creation of many thermospheric
empirical models throughout the decades, including the Jacchia (1970) model
and subsequent series; the Mass Spectrometer Incoherent Scatter model series
developed by Hedin (1987) and later improved by the Naval Research Laboratory
to become the Mass Spectrometer Incoherent Scatter Extended model (Picone et
al., 2002), and the Drag Temperature Model developed by Bruinsma (2015). In
this work, we will discuss results provided by the High Accuracy Satellite
Drag Model (HASDM; Storz et al., 2005), and the improved version of the
Jacchia model series described by Bowman et al. (2008), henceforth JB2008.

The understanding of the thermosphere response to magnetic storms and the
accurate capability of predicting subsequent satellite orbital drag effects is
of paramount interest of, e.g., the U.S. Federal Government and the private
sector (National Science and Technology Council, 2015b, a; Cakaj, 2021). The
correct orbital tracking of low-Earth orbit (LEO) satellites particularly in a
time window of 72 hours into the storm is of great interest for taking actions
such as pre-determined maneuvers as a means to avoid a satellite’s collision
with debris in space or even with other satellites (Pardini and Anselmo, 2009;
Wang, 2010; Lewis, 2019). Such actions may be characterized as an important
tool for preventing the occurrence of the Kessler Syndrome. First introduced
by Kessler and Cour-Palais (1978), the Kessler Syndrome suggests that by the
dawn of the 22nd century LEO regions will pose high risks for satellite
traffic due to the high probability of satellite collisions with significantly
increased space debris levels. As we will discuss later, the need for accurate
predictions of satellite orbital track during extreme magnetic storms is
twofold: (i), the number of satellites in LEO has been and will be
significantly increased; and (ii), prediction models and historical data
suggest solar activity will increase in the next years and decades.

In this paper, we will focus on four major LEO spacecraft missions that
carried/carry high-precision accelerometers that can be used to obtain high-
quality density data. These satellites are named CHAllenge Mini-satellite
Payload (CHAMP; Reigber et al., 2002), Gravity Recovery And Climate Experiment
(GRACE; Tapley et al., 2004), Gravity field and steady-state Ocean Circulation
Explorer (GOCE; Drinkwater et al., 2003), and Swarm (Siemes et al., 2016). We
will then pay particular attention to model performance during the only seven
extreme magnetic storms observed by CHAMP and GRACE.

##  3 Data and models

###  3.1 Data

CHAMP was launched in July 2000 at the initial altitude 456 km completing a
full orbit around Earth in 90 minutes with orbit inclination 87.25∘. CHAMP
completed a longitudinal cycle in $\sim$130 days. The acceleration measurement
precision by CHAMP was 3.0$\times 10^{-8}$ m/s2 with 0.1 Hz cadence (Bruinsma
et al., 2004). CHAMP re-entered in September 2010.

The GRACE mission was composed by two twin satellites, named GRACE-A and -B,
that were launched in March 2002. GRACE-B followed GRACE-A within a controlled
distance of $\sim$220 km. Therefore, densities derived by both spacecraft were
generally very similar. We then use data from the first spacecraft, hereafter
termed GRACE data. The initial altitudes of GRACE were around 500 km. The
orbital period of GRACE was 95 minutes with orbit declination 89.5∘. The
longitudinal coverage by GRACE was usually completed within $\sim$160 days.
The GRACE acceleration precision and cadence are 10 times higher in comparison
to CHAMP (Flury et al., 2008). GRACE re-entered in March 2018.

GOCE was launched in March 2009 at the initial altitude 270 km and orbit
inclination 96.5∘. The GOCE accelerometer precision was as high as 1.0$\times
10^{-11}$ m/s2, with measurements being performed within the bandwith
0.005-0.1 Hz (Bruinsma et al., 2014). GOCE re-entered in October 2013.

The Swarm mission is composed of three identical satellites that were launched
after November 2013 at orbits within 480-538 km altitude with orbital
inclination $\sim$88∘. The Swarm acceleration precision is 1.0$\times
10^{-11}$ m/s2 within a bandwith whose upper limit is 0.1 Hz (Siemes et al.,
2016). As a simple choice, we use only data from Swarm A, hereafter Swarm.
Swarm is the only mission analyzed in this study that is still in operation.

###  3.2 Models

The JB2008 empirical model computes thermospheric mass density using as inputs
solar indices and the 3-hour time resolution ap index accounting for low
geomagnetic activity. When geomagnetic activity intensifies (Dst $<$ –75 nT),
the model switches over to Dst as the proxy for geomagnetic activity (Bowman
et al., 2008). The solar indices used by JB2008 map the energy input into the
thermosphere provided by different solar irradiance sources, which in turn
affect density variability in the most important thermospheric layers (Tobiska
et al., 2008).

HASDM dynamically calibrates densities by assimilating the observed drag
effects on a carefully selected set of low-perigee calibration satellites and
it solves for the global thermospheric neutral density with 3-hour cadence by
“correcting” the JB2008 density prediction. Typically, 75 calibration
satellites are used to fit and correct the JB2008 density and the greater the
number of calibration satellites, the higher the accuracy. Because HASDM is
able to mitigate uncertainty in atmospheric effects, by resolving the
differences between calibration satellites and the core JB2008 model at each
epoch (time step), real-time nowcasting specification of the LEO drag
environment has improved significantly (e.g., Storz et al., 2005; Bowman et
al., 2008; Licata et al., 2021b; Calabia and Jin, 2021).

Forecasting is achieved with the JB2008 model, driven by forecasted solar and
geomagnetic drivers, with uncertainties growing significantly a few hours post
epoch. Marcos et al. (2010) demonstrated with detailed comparisons and
validation that JB2008 is the most accurate empirical model. In JB2008, the
1-sigma uncertainties at 400 km for a given epoch dropped from 15% to 8%. The
use of JB2008 in HASDM, and the assimilation of calibration satellite data
into the initial JB2008 density solutions, further dropped the comparable
HASDM density uncertainties to less than 5% for current epoch, lower altitude,
and quiet atmosphere.

Although HASDM is not available for direct use by the general scientific
community, Space Environment Technologies (SET) has made global density
outputs available to the public (Tobiska et al., 2021). The SET HASDM density
data base, available from 2000 to 2019, covers altitudes in the range 175-825
km, with resolution of 15∘ (longitude), 10∘ (latitude), and 25 km (altitude)
(Tobiska et al., 2021; Licata et al., 2021a).

##  4 Satellite coverage in LEO during magnetic storms, including extreme
events

Figure 1a summarizes the altitudes plotted as a function of time of the four
LEO satellites since the launching of CHAMP in May 2001 to December 2021.
Although there has been significant coverage from $\sim$250 to $\sim$540 km
altitude, this coverage is relatively sparse. As noted by Bruinsma et al.
(2021), the most well-covered storm was a mild event that occurred in April
2010, with CHAMP and GOCE around 300 km altitude, and GRACE around 480 km
altitude, all with very limited local solar time coverage.

Figure 1: Top panel: CHAMP, GRACE, GOCE, and Swarm altitudes. Bottom row, left
panel: superposed epoch analysis of satellite observation times plotted as a
function of local times and altitudes. The bin sizes are 5 km and 0.5 hour.
Right panel: the same as plotted in the left panel, but for the 7 extreme
magnetic storms (minimum SYM-H $\leq$ –250 nT) observed by CHAMP and GRACE.

Figure 1b shows the observation time distributions as a function of altitude
and local times for over 320 magnetic storms. The altitude data were collected
within the 72-hour interval after the storm main phase onset as described by
Oliveira et al. (2017). There is significant coverage above 400 km due to
CHAMP, GRACE, and Swarm. Additionally, there is high coverage by GOCE between
$\sim$180 km and $\sim$300 km concentrated around the dawn and dusk sectors.
The coverage between $\sim$300 km and $350$ km is very low and comes from the
late days of CHAMP’s observations. Many advancements in our understanding of
the thermosphere response to magnetic storms in the past two decades come from
these observations, including, e.g., thermosphere heating in polar cusp
regions (Lühr et al., 2004), traveling atmospheric disturbance (TAD)
propagation from high- to low-latitude regions (Bruinsma and Forbes, 2007;
Fujiwara and Miyoshi, 2006), orbital drag effects induced by magnetic storms
(Krauss et al., 2015; Oliveira and Zesta, 2019), high-latitude density
enhancements due to magnetospheric compressions by solar wind pressure pulses
(Connor et al., 2016; Shi et al., 2017), and thermosphere global time response
to magnetic storms (Sutton et al., 2005, 2009; Oliveira et al., 2017).

Altitude coverages of all extreme events (minimum SYM-H $\leq$ –250 nT) ever
observed by LEO spacecraft with high-accuracy accelerometers are shown in
Figure 1c. The figure shows observation time plotted as a function of local
time and altitude for the extreme storms observed by CHAMP (7 events) and
GRACE (6 events). In comparison with all events, the LEO observational
coverage during extreme magnetic storms is remarkably sparser and briefer,
particularly due to the rarity of extreme events and the low number of
concurring LEO missions operating during these events.

##  5 Comparing model performance during the extreme events observed by CHAMP
and GRACE

The left column of Figure 2 shows thermosphere density response to the extreme
storms whose coverages are shown in Figure 1c and dates shown in Table 1 of
Zesta and Oliveira (2019). Density data are superposed in epoch time $\times$
MLAT bins with size 15 minutes $\times$ 3∘. The zero epoch time is taken as
the time of each respective storm main phase onset. The panels show
$\log[\rho/\rho_{0}]$, with $\rho$ being the storm-time observed or modeled
density and $\rho_{0}$ being the background density estimated by JB2008 if
there was no storm.

As described by Oliveira et al. (2017), $\rho_{0}$ is obtained by the
following steps: (i) JB2008 is used to compute the neutral density $\rho_{1}$
during periods of low geomagnetic activity ($|$Dst$|$ $<$ 30 nT) to exclude
effects of high magnetospheric compressions and intense magnetic storms; (ii)
a polynomial expansion (Arlinghaus, 1994) with degree 15 is fitted to the
ratio $\rho/\rho_{1}$ to obtain a calibration function $f(t)$ (this degree
order provided the best polynomial fitting for a few control events, including
moderate and intense storms); and (iii) with density computed by JB2008 with
Dst = 0, i.e., $\rho(Dst=0)$, and f(t) interpolated for the whole dataset, the
background density $\rho_{0}$ is given by:

$\rho_{0}=\rho(Dst=0)\times f(t)$ (1)

Results shown in Figure 2a were published by Zesta and Oliveira (2019), but
here we removed effects caused by double storms, i.e., events in periods when
another CME impacted Earth when the magnetosphere was highly driven. The data
clearly show (i) high-latitude density enhancements due to the impacts of
shocks in the driver leading edges; (ii) TAD propagation effects from high to
low latitudes after storm main phase onset within 1-2 orbits ($\sim$2 hours,
(iii) thermosphere cooling effects at all latitudes $\sim$22 hours during
storm time presumably due to Nitric Oxide (NO) cooling effects (Mlynczak et
al., 2003; Knipp et al., 2017; Zesta and Oliveira, 2019) and (iv) thermosphere
overcooling during storm recovery meaning that density levels are now lower
than density levels during pre-storm periods (Lei et al., 2012; Zhang et al.,
2019). See Roman numerals in panel a for the corresponding effect described
above on density data.

Figure 2: Left column: superposed epoch analysis results of thermospheric
neutral mass density response using the methodology introduced by Oliveira et
al. (2017) with observations (a), JB2008 estimates (b), and HASDM estimates
(c). Right column: density averaged within –50∘ $\leq$ MLAT $\leq$ 50∘ in a
15-minute cadence (d); density estimate errors with respect to observations
computed as $\eta_{i}$ = ($\rho_{i}$ \- $\rho_{obs}$)/$\rho_{obs}$, with $i$ =
(JB2008, HASDM) (e); storm-time drag effects computed for GRACE’s orbit during
the 15 May 2005 extreme storm using the framework provided by Oliveira and
Zesta (2019) (f).

The comparison between JB2008 estimates (b) with observations (a) is quite
remarkable. JB2008 does not capture effects caused by pre-storm shock
compressions, nor does it capture TAD effects at storm main phase onset (white
arrows in panel a). Additionally, JB2008 predicts thermosphere heating at low
and high latitudes at the same time, which coincides with the time TADs take
to reach equatorial latitudes as shown by observations. Although JB2008
predicts a major thermosphere cooling at $\sim$22 hours, the model predicts a
secondary cooling occurring $\sim$8 hours later in comparison to observations,
and the model does not predict thermosphere overcooling effects during the
storm recovery phase. Quite the contrary, densities are slightly higher during
recovery in comparison to pre-storm density levels.

HASDM results (c) are considerably improved in comparison to JB2008 results.
Although HASDM does not clearly show high-latitude thermosphere heating due to
pre-storm shock compressions, the model captures some patterns of TAD
propagation, the extreme latitudinal cooling at t $\sim$ 22 hours, and the
thermosphere overcooling during storm recovery. Panel d shows mid- and low-
latitude densities averaged within – 50∘ $\leq$ MLAT $\leq$ 50∘ every 15
minutes throughout the storms. Our results show that (i) JB2008 and HASDM
densities agree moderately well with observations during pre-storm periods;
(ii) JB2008 underestimates and HASDM overestimates densities during storm main
phase, and (iii) HASDM estimates agree remarkably well with observations
during storm recovery, reproducing cooling effects quite well, but JB2008
density levels always surpass observations during the recovery phase. Licata
et al. (2021b) noted similar trends when comparing CHAMP and GRACE
observations with HASDM and JB2008 results for the Halloween storms. Panel e
shows the relative errors of the model estimates with respect to observations.
The shaded grey area corresponds to the $\pm$5% confidence interval suggested
by the U.S. Space Force for optimizing orbit predictions (Lewis, 2019).
Results show that HASDM errors stay within the confidence interval during all
storm periods, but the model introduces high errors during times of pre-storm
shock compressions. On the other hand, the standalone JB2008 performance is
considerably worse during all storm periods. The model underestimates density
during the main phase, but it considerably overestimates density during the
recovery phase, reaching error levels as high as 35% during late recovery. As
discussed in Oliveira and Zesta (2019) and Zesta and Oliveira (2019), JB2008
should be modified to accommodate effects induced by NO cooling for the
improvement of density predictions particularly during extreme magnetic
storms.

Finally, Figure 2f shows the comparison of observed and modeled drag effects
on GRACE during the 15 May 2005 extreme magnetic storm. The drag effects are
computed according to the framework developed by Oliveira and Zesta (2019).
The nearly straight grey line corresponds to the drag effect caused by the
background density $\rho_{0}$ if there was no storm. Observations (blue curve)
show that GRACE decayed by 55 m 72 hours after storm onset. This is 11 times
more severe in comparison to the orbital decay caused by the background
density. As a result of the errors introduced by JB2008 and HASDM during
different phases of the storm cycle (e), drag effects are underestimated by
the former and overestimated by the latter.

##  6 Discussion and conclusion

The Sputnik satellites opened the doors for human exploration in space
(Launius, 2004). Jacchia (1959) used ephemeris data collected by Sputnik
during an extreme magnetic storm in 1958 to link orbital drag effects with
geomagnetic activity for the first time. Many empirical models were created
and satellite missions were designed to study the dynamic thermosphere
response to magnetic storms throughout the decades. However, orbital drag
effects induced in LEO by extreme storms are yet relatively less understood
partly due to the rarity of extreme storms during the space age. The
International Geophysical Year (IGY) is marked by the adoption of the Dst
index that became a standard trademark in assessing the intensity of magnetic
storms. Interestingly, the largest yearly sunspot number ever recorded
occurred in 1957 (e.g., Clette and Lefèvre, 2016), only a year before Sputnik
launched and observed intense orbital drag effects (Jacchia, 1959). Although
human exploration in space had just begun, the overall solar activity
throughout the solar cycles (SCs) has been steadily decreasing since then.
Since the IGY, relatively few extreme events (Dst $\leq$ – 250 nT) have
occurred. Meng et al. (2019) reported on the occurrence of $\sim 40$ extreme
events, from which only 7 were observed by spacecraft with high-precision
accelerometers for density derivations. The last extreme magnetic storm took
place two SCs ago (SC23) on 15 May 2005.

Given the very low number of extreme magnetic storms amongst hundreds of
events on record, a way to investigate more extreme events is to look in the
past. Recently, many efforts have been undertaken to investigate magnetic
activity represented by Dst-like indices reconstructed from archival materials
generally drafted in the first half of the 20th century. Some of these events
occurred in October/November 1903 (Hayakawa et al., 2020b), September 1909
(Love et al., 2019b; Hayakawa et al., 2019a), May 1921 (Love et al., 2019a),
and March 1946 (Hayakawa et al., 2020a). These works showed that minimum Dst-
like values of these storms went below – 500 nT, which characterizes the event
as a superstorm. Each reference provides the respective Dst-like index data
for further investigations. For example, Oliveira et al. (2020) used the Dst-
like data for the first 3 events mentioned above along with real Dst data for
the March 1989 storm (Allen et al., 1989; Boteler, 2019) to show with JB2008
that, by comparison, long-lasting and less intense superstorms can induce more
severe drag effects than short-lasting and more intense superstorms. Further
investigations are needed here, as we have yet the unexplored events occurring
in September 1859 (Tsurutani et al., 2003; Hayakawa et al., 2019b), February
1872 (Silverman, 2008; Hayakawa et al., 2018), and March 1946 (Hayakawa et
al., 2020a).

Now looking to the future, a possibility to improve our understanding of
extreme orbital drag is to foster international collaboration for the
development of a central database of thermospheric density data. There has
been increased interest of the private sector in the elaboration and
construction of a very large constellation of satellites in LEO. For example,
SpaceX launched in 2018 the first satellite prototypes for the Starlink
project (Space Exploration Holdings, 2016). Starlink’s primary goal is to
create a megaconstellation of over 12,000 LEO satellites bellow 600 km
altitude with latitudinal distribution of 0.005-0.01 spacecraft per square
meter for the formation of a worldwide internet network (McDowell, 2020).
Another private company, OneWeb, intends so launch by the end of 2022 a
648-satellite constellation to provide worldwide internet service as well
(Barnett, 2016). Using an evolutionary model with parameters provided by the
respective operations application (Barnett, 2016; Space Exploration Holdings,
2016), Le May et al. (2018) showed that, within a 5-year operation time, the
probability of occurring a catastrophic collision involving an OneWeb
spacecraft is $\sim$5%, whereas the same for SpaceX is near 50%. However, the
authors did not consider any effects introduced by solar activity, and, if
they did, these figures would have certainly been higher. Therefore, the need
of accurate models for orbital track prediction during magnetic storms,
particularly during extreme events, is an important tool for preventing the
occurrence of the Kessler Syndrome in space.

A future large-scale thermospheric density data base can be used in studies
involving Machine Learning (ML) applications. ML studies have become very
popular in the field of Earth and Space Sciences in the past decade (e.g.,
Keesee et al., 2020; Smith et al., 2020; Bortnik and Camporeale, 2021; Haines
et al., 2021). For example, historical data sets (geomagnetic and solar
indices, sunspot numbers) prior to 2000 can be used for training a model to
predict storm drivers and the subsequent global thermospheric density and
orbital drag of a LEO satellite in a given location (Licata et al., 2020,
2021c). This can be accomplished, e.g., by the use of ML linear regression for
“training” the model to “learn” how to predict these effects (Rong and Bao-
wen, 2018). Such technique is named by Bortnik and Camporeale (2021) as Time
Series and Spatiotemporal Prediction, and is recognized by the authors as an
important step for ML applications in Earth and Space Sciences. As shown in
Figure 2, if HASDM provided remarkable results with the use of only 75
calibrated satellites, such results can be further improved by the
assimilation of large-scale data provided by a few hundreds satellites in LEO,
particularly for extreme events.

Finally, recent studies suggest that solar activity of SC25 will be
approximately the same as the solar activity of SC24 (e.g., Javaraiah, 2017).
However, McIntos et al. (2020) predict that SC25 will not only be stronger
than SC24, but SC25 will rival the magnitude of the strongest SCs on record.
Additionally, Javaraiah (2017) suggests that the transition between SC25 and
SC26 will coincide with the end of the current Gleissberg cycle. The
Gleisseberg cycle is characterized by the periodic occurrences of solar maxima
every 77-88 years or so (Gleissberg, 1967; Feynman and Ruzmaikin, 2014).
Therefore, all these predictions indicate that in the next years and decades
the Sun will shift from its relatively quiet conditions to a much more active
behavior. As a result, very accurate and precise thermospheric density models
will play crucial roles in keeping and guaranteeing the safety of the ever
increasing number of satellites in LEO by improved orbital tracking
particularly during extreme magnetic storms.

## Conflict of Interest Statement

The authors declare that the research was conducted in the absence of any
commercial or financial relationships that could be construed as a potential
conflict of interest.

## Author Contributions

DMO is the primary author of the article and performed the superposed epoch
analyses and drag computations. EZ contributed with the overall discussions
and helped the primary author set up the overall structure of the manuscript.
PMM and RLL performed the computations of the JB2008 and SET HASDM densities
for the extreme events. WKT and MDP provided support for the SET HASDM density
computations. HH provided importance guidance with respect to historical data
set discussions and interpretation of past and future solar activity
discussions. All co-authors provided editorial comments thus contributing to
the article and approved the submitted version.

## Funding

This work was funded by the NASA Space Weather Science Applications Operations
2 Research program, under grant number 20-SWO2R20-2-0014.

## Data Availability Statement

LEO satellite (ephemeris and density) data for CHAMP, GRACE, GOCE, and Swarm
were obtained from the repository provided by the Delft University of
Technology (http://thermosphere.tudelft.nl). Although not explicitly shown,
geomagnetic index data represented by SYM-H (1-minute resolution), were
downloaded from the World Data Center for Geomagnetism, Kyoto et al. (2015).
The JB2008 model is publicly available at the Space Environment Technologies
(SET) website (https://spacewx.com/jb2008/). Dst and solar index data are also
provided by SET with the source codes. The JB2008 and SET HASDM outputs can be
obtained from the repository located at
https://zenodo.org/record/4602380#.YRwSXC2cYkh.

## Acknowledgments

The authors of this manuscript warmly thank Delft University of Technology for
providing and making public the LEO satellite data used in this study. HH
thanks World Data Center for the production, preservation and dissemination of
the international sunspot number (WDC SILSO) and Data Analysis Center for
Geomagnetism and Space Magnetism (WDC for Geomagnetism at Kyoto) for providing
the international sunspot number and the Dst index.

## References

  * Allen et al. (1989) Allen, J., Sauer, H., Frank, L., and Reiff, P. (1989).  Effects of the March 1989 solar activity.  _Eos Transactions AGU_ 70, 1479–1488.  10.1029/89EO00409
  * Arlinghaus (1994) Arlinghaus, S. L. (1994).  _Practical Handbook of Curve Fitting_ (Boca Raton, FL: CRC Press)
  * Barnett (2016) Barnett, R. J. (2016).  _OneWeb Non-geostationary Satellite System: Technical Information to Supplement Schedule S_.  Tech. rep., Attachment to FCC Application SAT-LOI-20160428-00041, Washington, D.C.
  * Bortnik and Camporeale (2021) Bortnik, J. and Camporeale, E. (2021).  Ten Ways to Apply Machine Learning in Earth and Space Sciences.  _Eos Transactions AGU_ 202\.  10.1029/2021EO160257
  * Boteler (2019) Boteler, D. H. (2019).  A Twenty‐First Century View of the March 1989 Magnetic Storm.  _Space Weather_ 17, 1427–1441.  10.1029/2019SW002278
  * Bowman et al. (2008) Bowman, B. R., Tobiska, W. K., Marcos, F. A., Huang, C. Y., Lin, C. S., and Burke, W. J. (2008).  A new empirical thermospheric density model JB2008 using new solar and geomagnetic indices.  In _AIAA/AAS Astrodynamics Specialist Conference, AIAA 2008–6438_ (Honolulu, HI: American Institute of Aeronautics and Astronautics (AIAA) and American Astronautical Society (AAS)), 1–19
  * Bruinsma et al. (2021) Bruinsma, S., Fedrizzi, M., Yue, J., and Siemes, C. (2021).  Charting Satellite Courses in a Crowded Thermosphere.  _Eos Transactions AGU_ 102\.  10.1029/2021EO153475
  * Bruinsma (2015) Bruinsma, S. L. (2015).  The DTM-2013 thermosphere model.  _Journal of Space Weather and Space Climate_ 5\.  10.1051/swsc/2015001
  * Bruinsma et al. (2014) Bruinsma, S. L., Doornbos, E., and Bowman, B. R. (2014).  Validation of GOCE densities and evaluation of thermosphere models.  _Advances in Space Research_ 54, 576–585.  10.1016/j.asr.2014.04.008
  * Bruinsma and Forbes (2007) Bruinsma, S. L. and Forbes, J. M. (2007).  Global observations of traveling atmospheric disturbances (TADs) in the thermosphere.  _Geophysical Research Letters_ 34\.  10.1029/2007GL030243
  * Bruinsma et al. (2004) Bruinsma, S. L., Tamagnan, D., and Biancale, R. (2004).  Atmospheric densities derived from CHAMP/STAR accelerometer observations.  _Planetary and Space Science_ 62, 297–312.  10.1016/j.pss.2003.11.004
  * Cakaj (2021) Cakaj, S. (2021).  The Parameters Comparison of the “Starlink” LEO Satellites Constellation for Different Orbital Shells.  _Frontiers in Astronomy and Space Science_ 2, 7.  10.3389/frcmn.2021.643095
  * Calabia and Jin (2021) Calabia, A. and Jin, S. (2021).  Thermospheric Mass Density Disturbances Due to Magnetospheric Forcing From 2014–2020 CASSIOPE Precise Orbits.  _Journal of Geophysical Research: Space Physics_ 126, e2021JA029540.  10.1029/2021JA029540
  * Clette and Lefèvre (2016) Clette, F. and Lefèvre, L. (2016).  The New Sunspot Number: Assembling All Corrections.  _Solar Physics_ 291, 2629–2651.  10.1007/s11207-016-1014-y
  * Connor et al. (2016) Connor, H. K., Zesta, E., Fedrizzi, M., Shi, Y., Raeder, J., Codrescu, M. V., et al. (2016).  Modeling the ionosphere-thermosphere response to a geomagnetic storm using physics-based magnetospheric energy input: OpenGGCM-CTIM results.  _Journal of Space Weather and Space Climate_ 6, 1–15.  10.1051/swsc/2016019
  * Doornbos and Klinkrad (2006) Doornbos, E. and Klinkrad, H. (2006).  Modelling of space weather effects on satellite drag.  _Advances in Space Research_ 37, 1229–1239.  10.1016/j.asr.2005.04.097
  * Drinkwater et al. (2003) Drinkwater, M. R., Floberghagen, R., Haagmans, R., Muzi, D., and Popescu, A. (2003).  GOCE: ESA’s First Earth Explorer Core Mission.  In _Earth Gravity Field from Space – From Sensors to Earth Sciences_ , eds. G. Beutler, M. R. Drinkwater, R. Rummel, and R. Von Steiger (Dordrecht, The Netherlands: Springer Netherlands), vol. 17, 419–432.  10.1007/978-94-017-1333-7_36
  * Emmert (2015) Emmert, J. T. (2015).  Thermospheric mass density: A review.  _Advances in Space Research_ 56, 773–824.  10.1016/j.asr.2015.05.038
  * Feynman and Ruzmaikin (2014) Feynman, J. and Ruzmaikin, A. (2014).  The Centennial Gleissberg Cycle and its association with extended minima.  _Journal of Geophysical Research: Space Physics_ 119, 6027–6041.  10.1002/2013JA019478
  * Flury et al. (2008) Flury, J., Bettadpur, S., and Tapley, B. D. (2008).  Precise accelerometry onboard the GRACE gravity field satellite mission.  _Advances in Space Research_ 42, 1414–1423.  10.1016/j.asr.2008.05.004
  * Fujiwara and Miyoshi (2006) Fujiwara, H. and Miyoshi, Y. (2006).  Characteristics of the large-scale traveling atmospheric disturbances during geomagnetically quiet and disturbed periods simulated by a whole atmosphere general circulation model.  _Geophysical Research Letters_ 33\.  10.1029/2006GL027103
  * Gleissberg (1967) Gleissberg, W. (1967).  Secularly smoothed data on the minima and maxima of sunspot frequency.  _Solar Physics_ 2, 231–233.  10.1007/BF00155925
  * Haines et al. (2021) Haines, C., Owens, M. J., Barnard, L., Lockwood, M., Ruffenach, A., Boykin, K., et al. (2021).  Forecasting Occurrence and Intensity of Geomagnetic Activity With Pattern-Matching Approaches.  _Space Weather_ 19, e2020SW002624.  10.1029/2020SW002624
  * Hayakawa et al. (2019a) Hayakawa, H., Ebihara, Y., Cliver, E. W., Hattori, K., Toriumi, S., Love, J. J., et al. (2019a).  The extreme space weather event in September 1909.  _Monthly Notices of the Royal Astronomical Society_ 484, 4083–4099.  10.1093/mnras/sty3196
  * Hayakawa et al. (2020a) Hayakawa, H., Ebihara, Y., Pevtsov, A., Bhaskar, A., Karachik, N., and Oliveira, D. M. (2020a).  Intensity and Time Series of Extreme Solar-Terrestrial Storm in March 1946.  _Monthly Notices of the Royal Astronomical Society_ 497, 5507–5517.  10.1093/mnras/staa1508
  * Hayakawa et al. (2018) Hayakawa, H., Ebihara, Y., Willis, D. M., Hattori, K., Giunta, A. S., Wild, M. N., et al. (2018).  The Great Space Weather Event during 1872 February Recorded in East Asia.  _The Astrophysical Journal_ 862, 1–10.  10.3847/1538-4357/aaca40
  * Hayakawa et al. (2019b) Hayakawa, H., Ebihara, Y., Willis, D. M., Toriumi, S., Iju, T., Hattori, K., et al. (2019b).  Temporal and Spatial Evolutions of a Large Sunspot Group and Great Auroral Storms around the Carrington Event in 1859.  _Space Weather_ 17, 1553–1569.  10.1029/2019SW002269
  * Hayakawa et al. (2020b) Hayakawa, H., Ribeiro, P., Vaquero, J. M., Gallego, M. C., Knipp, D. J., Mekhaldi, F., et al. (2020b).  The Extreme Space Weather Event in 1903 October/November: An Outburst from the Quiet Sun.  _The Astrophysical Journal Letters_ 897, L10.  10.3847/2041-8213/ab6a18
  * Hedin (1987) Hedin, A. E. (1987).  MSIS-86 Thermospheric Model.  _Journal of Geophysical Research_ 92, 4649–4662.  10.1029/JA092iA05p04649
  * Jacchia (1959) Jacchia, L. G. (1959).  Corpuscular radiation and the acceleration of artificial satellites.  _Nature_ 183, 1662–1663.  10.1038/1831662a0
  * Jacchia (1970) Jacchia, L. G. (1970).  New static models of the thermosphere and exosphere with empirical temperature profiles.  In _Spec. Rep. 313_ (Cambridge, Massachusetts: Smithson, Astrophys. Obs.), 1–82
  * Javaraiah (2017) Javaraiah, J. (2017).  Will Solar Cycles 25 and 26 Be Weaker than Cycle 24?  _Solar Physics_ 292\.  10.1007/s11207-017-1197-x
  * Keesee et al. (2020) Keesee, A. M., Pinto, V., Coughlan, M., Lennox, C., Mahmud, M. S., and Connor, H. K. (2020).  Comparison of Deep Learning Techniques to Model Connections Between Solar Wind and Ground Magnetic Perturbations.  _Frontiers in Astronomy and Space Science_ 7, 72.  10.3389/fspas.2020.550874
  * Kessler and Cour-Palais (1978) Kessler, D. J. and Cour-Palais, B. G. (1978).  Collision frequency of artificial satellites: The creation of a debris belt.  _Journal of Geophysical Research_ 83, 2637–2646.  10.1029/JA083iA06p02637
  * Knipp et al. (2017) Knipp, D. J., Pette, D. V., Kilcommons, L. M., Isaacs, T. L., Cruz, A. A., Mlynczak, M. G., et al. (2017).  Thermospheric nitric oxide response to shock-led storms.  _Space Weather_ 15, 325–342.  10.1002/2016SW001567
  * Krauss et al. (2015) Krauss, S., Temmer, M., Veronig, A., Baur, O., and Lammer, H. (2015).  Thermosphere and geomagnetic response to interplanetary coronal mass ejections observed by ACE and GRACE: Statistical results.  _Journal of Geophysical Research: Space Physics_ 120, 8848–8860.  10.1002/2015JA021702
  * Launius (2004) Launius, R. D. (2004).  _Frontiers of Space Exploration_ (Westport, Connecticut: Greenwood Press)
  * Le May et al. (2018) Le May, S., Gehly, S., Carter, B. A., and Flegel, S. (2018).  Space debris collision probability analysis for proposed global broadband constellations.  _Acta Astronautica_ 151, 445–455.  10.1016/j.actaastro.2018.06.036
  * Lei et al. (2012) Lei, J., Burns, A. G., Thayer, J. P., Wang, W., Mlynczak, M. G., Hunt, L. A., et al. (2012).  Overcooling in the upper thermosphere during the recovery phase of the 2003 October storms.  _Journal of Geophysical Research_ 117\.  10.1029/2011JA016994
  * Lewis (2019) Lewis, C. D. (2019).  The DARPA Space Environment Exploitation (SEE) Program.  In _Friday AGU Chapman Conference Abstracts_ (Presented at the 2019 AGU Chapman Conference on Scientific Challenges Pertaining to Space Weather Forecasting Including Extremes, Pasadena, CA, 11-15 February), 1
  * Licata et al. (2021a) [Dataset] Licata, R. J., Mehta, P. M., Kent Tobiska, W., Bowman, B. R., and Pilinski, M. D. (2021a).  CHAMP and GRACE-A Density Estimates with Associated HASDM and JB2008 Predictions.  10.5281/zenodo.4602380
  * Licata et al. (2021b) Licata, R. J., Mehta, P. M., Tobiska, W. K., Bowman, B. R., and Pilinski, M. D. (2021b).  Qualitative and Quantitative Assessment of the SET HASDM Database.  _Space Weather_ 19, e2021SW002798.  10.1029/2021SW002798
  * Licata et al. (2021c) Licata, R. J., Mehta, P. M., Weimer, D. R., and Kent Tobiska, W. (2021c).  Improved Neutral Density Predictions through Machine Learning Enabled Exospheric Temperature Model.  _Earth and Space Science Open Archive_ 10.1002/essoar.10507687.1
  * Licata et al. (2020) Licata, R. J., Tobiska, W. K., and Mehta, P. M. (2020).  Benchmarking Forecasting Models for Space Weather Drivers.  _Space Weather_ 18, e2020SW002496.  10.1029/2020SW002496
  * Love et al. (2019a) Love, J. J., Hayakawa, H., and Cliver, E. W. (2019a).  Intensity and impact of the New York Railroad superstorm of May 1921\.  _Space Weather_ 17, 1281–1292.  10.1029/2019SW002250
  * Love et al. (2019b) Love, J. J., Hayakawa, H., and Cliver, E. W. (2019b).  On the Intensity of the Magnetic Superstorm of September 1909.  _Space Weather_ 17, 37–45.  10.1029/2018SW002079
  * Lühr et al. (2004) Lühr, H., Tother, M., Köhler, W., Ritter, P., and Grunwaldt, L. (2004).  Thermospheric up-welling in the cusp region: Evidence from CHAMP observations.  _Geophysical Research Letters_ 31, 1–4.  10.1029/2003GL019314
  * Marcos et al. (2010) Marcos, F. A., Delay, S. H., and Sutton, E. K. (2010).  Toward next level satellite drag modeling.  In _AIAA/AAS Astrodynamics Specialist Conference and Exhibit_ (Toronto, Ontario Canada), 1–20.  10.2514/6.2010-7840
  * McDowell (2020) McDowell, J. C. (2020).  The Low Earth Orbit Satellite Population and Impacts of the SpaceX Starlink Constellation.  _The Astrophysical Journal Letters_ 892, 1–10.  10.3847/2041-8213/ab8016
  * McIntos et al. (2020) McIntos, S. W., Chapman, S., Leamon, R. J., Egeland, R., and Watkins, N. W. (2020).  Overlapping Magnetic Activity Cycles and the Sunspot Number: Forecasting Sunspot Cycle 25 Amplitude.  _Solar Physics_ 295\.  10.1007/s11207-020-01723-y
  * Meng et al. (2019) Meng, X., Tsurutani, B. T., and Mannucci, A. J. (2019).  The Solar and Interplanetary Causes of Superstorms (Minimum Dst $\leq$ –250 nT) During the Space Age.  _Journal of Geophysical Research: Space Physics_ 124, 3926–3948.  10.1029/2018JA026425
  * Mlynczak et al. (2003) Mlynczak, M. G., Martin-Torres, F. J., Russell, J., Beaumont, K., Jacobson, S., Kozyra, J., et al. (2003).  The natural thermostat of nitric oxide emission at 5.3 $\mu$m in the thermosphere observed during the solar storms of April 2002.  _Geophysical Research Letters_ 30\.  10.1029/2003GL017693
  * National Science and Technology Council (2015a) National Science and Technology Council (2015a).  _National Space Weather Action Plan_.  Tech. rep., Executive Office of the President of the United States, Washington, D.C.
  * National Science and Technology Council (2015b) National Science and Technology Council (2015b).  _National Space Weather Strategy_.  Tech. rep., Executive Office of the President of the United States, Washington, D.C.
  * Nwankwo et al. (2021) Nwankwo, V. U. J., Denig, W., Chakrabarti, S. K., Ajakaiye, M. P., Fatokun, J., Akanni, A. W., et al. (2021).  Atmospheric drag effects on modelled low Earth orbit (LEO) satellites during the July 2000 Bastille Day event in contrast to an interval of geomagnetically quiet.  _Annales Geophysicae_ 39, 397–412.  10.5194/angeo-39-397-2021
  * Oliveira and Zesta (2019) Oliveira, D. M. and Zesta, E. (2019).  Satellite Orbital Drag During Magnetic Storms.  _Space Weather_ 17, 1510–1533.  10.1029/2019SW002287
  * Oliveira et al. (2020) Oliveira, D. M., Zesta, E., Hayakawa, H., and Bhaskar, A. (2020).  Estimating satellite orbital drag during historical magnetic superstorms.  _Space Weather_ 18, e2020SW002472.  10.1029/2020SW002472
  * Oliveira et al. (2017) Oliveira, D. M., Zesta, E., Schuck, P. W., and Sutton, E. K. (2017).  Thermosphere global time response to geomagnetic storms caused by coronal mass ejections.  _Journal of Geophysical Research: Space Physics_ 122, 10,762–10,782.  10.1002/2017JA024006
  * Pardini and Anselmo (2009) Pardini, C. and Anselmo, L. (2009).  Assessment of the consequences of the Fengyun-1C breakup in low Earth orbit.  _Advances in Space Research_ 44, 545–557.  10.1016/j.asr.2009.04.014
  * Picone et al. (2002) Picone, J. M., Hedin, A. E., Drob, D. P., and Aikin, A. C. (2002).  NRLMSISE-00 empirical model of the atmosphere: Statistical comparisons and scientific issues.  _Journal of Geophysical Research_ 107, SIA 15–1–SIA 15–16.  10.1029/2002JA009430
  * Prölss (2011) Prölss, G. (2011).  Density perturbations in the upper atmosphere caused by the dissipation of solar wind energy.  _Surveys in Geophysics_ 32, 101–195.  10.1007/s10712-010-9104-0
  * Reigber et al. (2002) Reigber, C., Lühr, H., and Schwintzer, P. (2002).  CHAMP mission status.  _Advances in Space Research_ 30, 129–134.  10.1016/S0273-1177(02)00276-4
  * Rong and Bao-wen (2018) Rong, S. and Bao-wen, Z. (2018).  The research of regression model in machine learning field.  _MATEC Web Conf._ 176, 01033.  10.1051/matecconf/201817601033
  * Shi et al. (2017) Shi, Y., Zesta, E., Connor, H. K., Su, Y.-J., Sutton, E. K., Huang, C. Y., et al. (2017).  High-latitude thermosphere neutral density response to solar wind dynamic pressure enhancement.  _Journal of Geophysical Research: Space Physics_ 122, 11,559–11,578.  10.1002/2017JA023889
  * Siemes et al. (2016) Siemes, C., de Teixeira da Encarnação, J., Doornbos, E., van den IJssel, J., Kraus, J., Peres̆tý, R., et al. (2016).  Swarm accelerometer data processing from raw accelerations to thermospheric neutral densities.  _Earth, Planets and Space_ 68, 1–16.  10.1186/s40623-016-0474-5
  * Silverman (2008) Silverman, S. M. (2008).  Low-latitude auroras: The great aurora of 4 February 1872.  _Journal of Atmospheric and Solar-Terrestrial Physics_ 70, 1301–1308.  10.1016/j.jastp.2008.03.012
  * Smith et al. (2020) Smith, A. W., Rae, J., Forsyth, C., Oliveira, D. M., Freeman, P. M., and Jackson, D. (2020).  Probabilistic Forecasts of Storm Sudden Commencements from Interplanetary Shocks using Machine Learning.  _Space Weather_ 18, e2020SW002603.  10.1029/2020SW002603
  * Space Exploration Holdings (2016) Space Exploration Holdings (2016).  _SpaceX Non-geostationary Satellite System_.  Tech. rep., Attachment to FCC Application SAT-LOA-20161115-00118, Washington, D.C.
  * Storz et al. (2005) Storz, M. F., Bowman, B. R., Branson, J. I., Casalic, S. J., and Tobiska, W. K. (2005).  High accuracy satellite drag model (HASDM).  _Advances in Space Research_ 36, 2497–2505.  10.1016/j.asr.2004.02.020
  * Sutton et al. (2009) Sutton, E. K., Forbes, J. M., and Knipp, D. J. (2009).  Rapid response of the thermosphere to variations in Joule heating.  _Journal of Geophysical Research_ 114\.  10.1029/2008JA013667
  * Sutton et al. (2005) Sutton, E. K., Forbes, J. M., and Nerem, R. S. (2005).  Global thermospheric neutral density and wind response to the severe 2003 geomagnetic storms from CHAMP accelerometer data.  _Journal of Geophysical Research_ 110, 1–10.  10.1029/2004JA010985
  * Tapley et al. (2004) Tapley, B. D., Bettadpur, S., Watkins, M., and Reigber, C. (2004).  The gravity recovery and climate experiment: Mission overview and early results.  _Geophysical Research Letters_ 31\.  10.1029/2004GL019920
  * Tobiska et al. (2008) Tobiska, W. K., Bouwer, D., and R.Bowman, B. (2008).  The development of new solar indices for use in thermospheric density modeling.  _Journal of Atmospheric and Solar-Terrestrial Physics_ 70, 803–819.  10.1016/j.jastp.2007.11.001
  * Tobiska et al. (2021) Tobiska, W. K., Bowman, B. R., Bouwer, D., Cruz, A., Wahl, K., Pilinski, M., et al. (2021).  The SET HASDM density database.  _Space Weather_ 19, e2020SW002682.  10.1029/2020SW002682
  * Tsurutani et al. (2003) Tsurutani, B. T., Gonzalez, W. D., Lakhina, G. S., and Alex, S. (2003).  The extreme magnetic storm of 1–-2 September 1859.  _Journal of Geophysical Research_ 108\.  10.1029/2002JA009504
  * Wang (2010) Wang, T. (2010).  Analysis of Debris from the Collision of the Cosmos 2251 and the Iridium 33 Satellites.  _Science $\&$ Global Security_ 18, 87–118.  10.1080/08929882.2010.493078
  * World Data Center for Geomagnetism, Kyoto et al. (2015) [Dataset] World Data Center for Geomagnetism, Kyoto, Nose, M., Iyemori, T., Sugiura, M., and Kamei, T. (2015).  Geomagnetic AE index.  10.17593/15031-54800
  * Zesta and Huang (2016) Zesta, E. and Huang, C. Y. (2016).  Satellite orbital drag.  In _Space Weather Fundamentals_ , ed. G. V. Khazanov (Boca Raton, FL: CRC Press), chap. 18. 329–351
  * Zesta and Oliveira (2019) Zesta, E. and Oliveira, D. M. (2019).  Thermospheric heating and cooling times during geomagnetic storms, including extreme events.  _Geophysical Research Letters_ 46, 12,739–12,746.  10.1029/2019GL085120
  * Zhang et al. (2019) Zhang, Y., Paxton, L. J., Lu, G., and Yee, S. (2019).  Impact of nitric oxide, solar EUV and particle precipitation on thermospheric density decrease.  _Journal of Atmospheric and Solar-Terrestrial Physics_ 182, 147–154.  10.1016/j.jastp.2018.11.016