# Results on heavy ion physics at LHCb

Giacomo Graziani  on behalf of the LHCb Collaboration  INFN, Sezione di
Firenze<EMAIL_ADDRESS>
###### Abstract

In the last years, the LHCb experiment established itself as an important
contributor to heavy ion physics by exploiting some of its specific features.
Production of particles, notably heavy flavour states, can be studied in
$p$-$p$, $p$-Pb and Pb-Pb collisions at LHC energies in the forward rapidity
region (pseudorapidity between 2 and 5), providing measurements which are
highly complementary to the other LHC experiments. Moreover, owing to its
forward geometry, the detector is also well suited to study fixed-target
collisions, obtained by impinging the LHC beams on gas targets with different
mass numbers. In this configuration, $p$-A collisions can be studied at the
relatively unexplored scale of $\sqrt{s_{\scriptscriptstyle\rm NN}}$ $\sim$
100 GeV, also providing valuable inputs to cosmic ray physics. An overview of
the measurements obtained so far by the LHCb ion program is presented.

##  1 Introduction

The LHCb experiment [1] has been conceived with the main goal of studying
heavy flavour physics in $p$-$p$ collisions at the LHC, exploiting the
unprecedented yield of $b$-quark pairs, which are mainly produced at small
angles with respect to the direction of the colliding proton beams. The
detector is therefore designed as a forward spectrometer covering the
pseudorapidity region $2<\eta<5$ and providing excellent vertexing, tracking
and particle identification capabilities for the reconstruction of heavy
flavour decays. Another key feature is the online selection system, consisting
of a hardware level with high output bandwidth (up to 1 MHz), followed by a
software level providing high flexibility.

Though heavy ion physics was not among the original motivations for the
experiment, the detector capabilities offer some unique possibilities also for
this field:

  * •

the forward acceptance, with a fully instrumented detector, highly
complementary to the other LHC experiments;

  * •

the excellent reconstruction performance, unrivaled at the LHC, for exclusive
heavy flavour states down to null transverse momentum ($p_{\mathrm{T}}$),
disentangling charmed particles produced promptly from those coming from
$b$-hadron decays;

  * •

the possibility to operate the detector in fixed target mode, by injecting
small amount of noble gas (He, Ne and Ar) in the LHC vacuum [2] and studying
beam-gas collisions, for which the forward geometry of the detector is very
well suited.

On the other hand, the most central Pb-Pb collisions can’t be properly
reconstructed due to the high track density in the forward region. LHCb is
therefore more suited for smaller collision systems like $p$-Pb, but can
contribute also to study peripheral Pb-Pb collisions. It has to be noted that
studies in small collision systems have become increasingly important to the
interpretation of relativistic heavy ion collisions after the discovery,
performed at the LHC, of sizeable collective behaviour even in $p$-$p$
collisions [3].

A comparison among the kinematic reaches of LHCb and the other LHC experiments
in $p$-A collisions, in terms of the Bjorken-$x$ value of the nucleon in the
nuclear target and the squared parton-parton invariant mass $Q^{2}$, is shown
in Figure 1. Two regions are covered in $p$-Pb collisions, depending on the
orientation of the proton and lead beams. The so-called forward configuration,
when the proton beam points toward the detector (i.e., it enters the detector
region from its vertex detector), corresponds to values of $x$ down to
$10^{-5}$, where gluon saturation is expected to occur. In the backward
configuration, when the Pb beam points toward the detector, measurements are
sensitive to the anti-shadowing region up to $x\sim 0.1$.

Figure 1: Kinematic reach corresponding to the acceptance of the four LHC
experiments in $p$-Pb collisions. The kinematic regions accessible in the
fixed-target configuration at LHCb and in $e$-$p$ collisions at HERA are also
shown.

In fixed target collisions, the nucleon-nucleon centre-of-mass (c.m.) energy
reaches 110 GeV for proton beams of 6.5 TeV energy. Measurements at this
energy scale, intermediate between SPS fixed target experiments and beam-beam
collisions at LHC, provides an additional test bed to explore the energy
evolution in the dynamics of nuclear matter. The detector acceptance
corresponds to mid and backward rapidities in the c.m. frame:
$-2.8<y^{*}<0.2$. As also depicted in Figure 1, this gives access to large
values of $x$ in the target nucleon, where nuclear PDFs are modified by the
EMC effect and where the contribution of a possible intrinsic heavy quark
content in the nucleon could be substantiated.

##  2 $p$-Pb collisions

The experiment collected a first dataset of $p$-Pb collisions at
$\sqrt{s_{\scriptscriptstyle\rm NN}}=5.02$ TeV in 2013, corresponding to
integrated luminosities of 1.1 and 0.5 $\mbox{\,nb}^{-1}$ in the forward and
backward configuration, respectively. After this succesfull experience, larger
luminosities were delivered to LHCb during the 2016 $p$-Pb run at
$\sqrt{s_{\scriptscriptstyle\rm NN}}=8.16$ TeV, with 13.6 and 30.8
$\mbox{\,nb}^{-1}$ recorded. The results published so far by LHCb focus on
production of heavy flavour states, whose modification with respect to $p$-$p$
collisions constitutes a major probe for the hot and dense matter, known as
quark-gluon plasma (QGP), produced in Pb-Pb collisions at LHC. In $p$-Pb
collisions QGP is not expected to be formed and the system is considered to be
a reference for the understanding of cold nuclear matter (CNM) effects.

Figure 2: Distributions of (left) reconstructed mass and (middle) pseudo-
proper time $t_{z}\equiv(z_{{{J\mskip-3.0mu/\mskip-2.0mu\psi\mskip
2.0mu}}}-z_{PV})\times(M/p_{z})_{{{J\mskip-3.0mu/\mskip-2.0mu\psi\mskip
2.0mu}}}$ for the ${{J\mskip-3.0mu/\mskip-2.0mu\psi\mskip
2.0mu}}\rightarrow{\mu^{+}\mu^{-}}$ decay candidates from the 8 TeV Pb-$p$
sample [4] in the rapidity bin $-4.0<y^{*}<-3.5$
($z_{{{J\mskip-3.0mu/\mskip-2.0mu\psi\mskip 2.0mu}}}$ and $z_{PV}$ are the
reconstructed longitudinal position of the
${J\mskip-3.0mu/\mskip-2.0mu\psi\mskip 2.0mu}$ decay vertex and the
interaction primary vertex, respectively). The result of a fit to determine
prompt signal, from-$b$ signal and background fractions is overdrawn. In the
right plot. mass distribution for the $\mathchar 28935\relax{(nS)}$ candidates
from the same sample in $-5.0<y^{*}<-2.5$ [5].

Figure 2 illustrates the ability to distinguish promptly produced
${{J\mskip-3.0mu/\mskip-2.0mu\psi\mskip 2.0mu}}$ mesons from those produced in
$b$-hadron decays in reconstructed ${{J\mskip-3.0mu/\mskip-2.0mu\psi\mskip
2.0mu}}\rightarrow{\mu^{+}\mu^{-}}$ decays. The nuclear modification factor
$R_{pPb}\equiv\sigma_{p\text{Pb}}/(A_{\text{Pb}}~{}\sigma_{pp})$, where
$\sigma_{p\text{Pb}}$ and $\sigma_{pp}$ are the production cross-sections in
the two collision systems and $A_{\text{Pb}}=208$ is the lead mass number, can
be measured separately for the two components. The result for the prompt
component obtained from the 8 TeV sample [4] is shown in Figure 3 as a
function of $p_{\mathrm{T}}$ and the c.m. rapidity $y^{*}$. A clear
suppression with respect to $p$-$p$ collisions is observed at forward rapidity
and low $p_{\mathrm{T}}$, compatible with the expected effect from nuclear PDF
modifications (shadowing) as computed in the framework of NRQCD factorisation
using several collinear nuclear PDF sets with the HELAC-Onia package [6, 7].
However, the rapidity dependence can also be well explained by the coherent
energy loss model [8]. The result at forward rapidity is also compatible with
the latest calculations based on the Color Glass Condensate model [9]. More
measurements, notably on Drell-Yan production, are needed to disentangle these
physical effects [10].

Figure 3: Nuclear modification factor $R_{pPb}$ for prompt
${J\mskip-3.0mu/\mskip-2.0mu\psi\mskip 2.0mu}$ production as a function of
(left plot) $y^{*}$ and (right plot) $p_{\mathrm{T}}$ for the forward
configuration [4].

If final-state effects in CNM can’t be excluded from the
${J\mskip-3.0mu/\mskip-2.0mu\psi\mskip 2.0mu}$ result, they are definetely
needed to explain the different modifications of quarkonia states. The first
measurement of $\psi{(2S)}$ production in the 5 TeV sample [11] indicated a
larger suppression with respect to ${J\mskip-3.0mu/\mskip-2.0mu\psi\mskip
2.0mu}$, in agreement with results from the other LHC experiments [12, 13, 14,
15]. A more recent result [5] shows evidence for different suppression among
the three $\mathchar 28935\relax{(nS)}$ states, which are cleanly observed in
the 8 TeV sample (see Figure 2).

Figure 4: Ratio between $\mathchar 28935\relax$ nuclear modification factors
as a function of rapidity for (left plot) $\mathchar
28935\relax{(2S)}$/$\mathchar 28935\relax{(1S)}$ and (right plot) $\mathchar
28935\relax{(3S)}$/$\mathchar 28935\relax{(1S)}$ [5].

The result, shown in Figure 4, is nicely described in the framework of the
“comovers” model [16], where dissociation of the quarkonia states is
attributed to interaction with final-state particles which are close in phase-
space. This relatively large effect needs to be taken into account in the
interpretation of the spectacular suppression of $\Upsilon$ states recently
observed by CMS in Pb-Pb collisions [17].

Production of open charm states, ${D}^{0}$ mesons [18] and ${\mathchar
28931\relax}^{+}_{c}$ baryons [19], has also been measured with large
statistics already in the 5 TeV sample. Data are precise enough to constrain
nuclear PDFs, assuming that initial-state effects dominate the observed
nuclear modification. The ${\mathchar 28931\relax}^{+}_{c}$/${D}^{0}$ ratio is
an important input to the hadronisation phenomenology, since baryon
enhancement in heavy ion collisions, notably at low $p_{\mathrm{T}}$, is
expected from production via coalescence and is affected by the thermal
properties of the nuclear medium. LHCb measures a ratio between the
${\mathchar 28931\relax}^{+}_{c}$ and ${D}^{0}$ prompt production around 0.3,
with no evidence of strong dependence on the rapidity or transverse momentum
(see Figure 5). Results are in substantial agreement with HELAC-onia
computations where the only nuclear effects are due to PDF modifications and
largely cancel in the ratio.

Figure 5: Ratio between ${\mathchar 28931\relax}^{+}_{c}$ and ${D}^{0}$
production in $p$-Pb collisions at 5 TeV, as a function of (left plot) $y^{*}$
and (right plot) $p_{\mathrm{T}}$ for the forward configuration [19]. The
result from ALICE [20] at mid rapidities is also shown for comparison.

Figure 6: In the upper plots, reconstructed mass distribution of candidates
for the (left) ${{{B}^{+}}}\rightarrow{{\kern
1.99997pt\overline{\kern-1.99997ptD}{}}{}^{0}}{{\pi}^{+}}$ and (right)
${{\mathchar 28931\relax}^{0}_{b}}\rightarrow{{\mathchar
28931\relax}^{+}_{c}}{{\pi}^{-}}$ decays in the $p$-Pb forward configuration
data at 8 TeV. In the lower plots, the resulting nuclear modification factor
$R_{pPb}$ for ${{B}^{+}}$ meson production is plotted as a function of (left
plot) $y^{*}$ and (right plot) $p_{\mathrm{T}}$ for the forward configuration
[21]. Results are compared with the measurement for
${J\mskip-3.0mu/\mskip-2.0mu\psi\mskip 2.0mu}$ from $b$-hadron decays [4] and
with HELAC-onia calculations using different nuclear PDF sets.

For the first time, exclusive $b$-hadron decays have been reconstructed in
nuclear collisions [21]. Clean samples of four decay modes
(${{{B}^{+}}}\rightarrow{{\kern
1.99997pt\overline{\kern-1.99997ptD}{}}{}^{0}}{{\pi}^{+}},{{{B}^{+}}}\rightarrow{{J\mskip-3.0mu/\mskip-2.0mu\psi\mskip
2.0mu}}{{K}^{+}},{{B}^{0}}\rightarrow{{D}^{-}}{{\pi}^{+}},{{\mathchar
28931\relax}^{0}_{b}}\rightarrow{{\mathchar 28931\relax}^{+}_{c}}{{\pi}^{-}}$)
have been obtained from the 8 TeV $p$-Pb dataset, down to $p_{\mathrm{T}}$
values below the hadron mass. The ${{B}^{+}}$ modes have been used to measure
the nuclear modification factor as a function of $y^{*}$ and $p_{\mathrm{T}}$.
The results, shown in Figure 6, confirm the suppression pattern observed with
detached ${J\mskip-3.0mu/\mskip-2.0mu\psi\mskip 2.0mu}$, consistently with
nuclear PDF effects. The ratio between prompt production of ${\mathchar
28931\relax}^{0}_{b}$ baryons and ${B}^{0}$ mesons is also measured and found
consistent with the $p$-$p$ measurement.

##  3 Pb-Pb collisions

A first small sample of Pb-Pb collisions was recorded by LHCb during the 2015
run, corresponding to about $10~{}\mu\text{b}^{-1}$. The tracking detector
performance in this challenging environment was found to be satisfactory for
events of centrality above 50%.

Figure 7:  Candidates for ${J\mskip-3.0mu/\mskip-2.0mu\psi\mskip 2.0mu}$
exclusive production in Pb-Pb collisions at 5 TeV [22], reconstructed from
${{J\mskip-3.0mu/\mskip-2.0mu\psi\mskip 2.0mu}}\rightarrow{\mu^{+}\mu^{-}}$
decays. In the upper plot, the reconstructed mass distribution is shown in an
extended mass range where, beside the clean
${J\mskip-3.0mu/\mskip-2.0mu\psi\mskip 2.0mu}$ peak, a small signal for
$\psi{(2S)}$ can also be seen. In the lower plot, the
$\log(p_{\mathrm{T}}^{2})$ distribution in the
${J\mskip-3.0mu/\mskip-2.0mu\psi\mskip 2.0mu}$ mass range ($3096.9\pm
65~{}\mathrm{\,Me\kern-1.00006ptV}$) is fitted with templates for coherent and
incoherent production.

The first preliminary physics result [22] has been obtained from
ultraperipheral collisions (UPC), where hadron photoproduction is enhanced by
the large photon flux from the lead nuclei. The observation of photoproduction
of heavy flavour states, providing a hard scale for perturbative QCD
calculations, is particularly interesting to explore the gluon density down to
the saturation region at $x\sim 10^{-5}$. The exclusive production of
${J\mskip-3.0mu/\mskip-2.0mu\psi\mskip 2.0mu}$ is cleanly observed (see Figure
7). The excellent $p_{\mathrm{T}}$ resolution allows to distinguish coherent
and uncoherent production, whose $p_{\mathrm{T}}$ distributions are found to
be well described by templates obtained with the STARlight generator [23].

The accuracy of the result is limited by the size of the data sample, but this
analysis demonstrates the LHCb potential for physics in Pb-Pb UPC. During the
run performed in november 2018, an integrated luminosity of
$210~{}\mu\text{b}^{-1}$ was collected, providing the possibility for a
precision measurement of exclusive ${J\mskip-3.0mu/\mskip-2.0mu\psi\mskip
2.0mu}$ and $\psi{(2S)}$ production, and possibly to access rarer channels
like $\mathchar 28935\relax$ production and light-by-light scattering.

##  4 Fixed-target collisions

Figure 8: Summary of collected fixed-target physics samples at LHCb. The data
size is given in terms of delivered protons (ions) on target (POT). For a
nominal target pressure of $2\times 10^{-7}$ mbar, $10^{22}$ POT correspond to
an integrated luminosity of about $5$ nb-1 per meter of gas, though the actual
target pressure and the data taking efficiency vary among samples.

Samples of beam-gas collisions with proton and lead beams of different energy,
and with the three possible targets (He, Ne and Ar) have been collected by
LHCb during the LHC Run 2, as summarized in Figure 8. The largest sample,
$p$-Ne collisions at $\sqrt{s_{\scriptscriptstyle\rm
NN}}=69~{}\mathrm{\,Ge\kern-1.00006ptV}$, corresponds to an integrated
luminosity of about 100 $\mbox{\,nb}^{-1}$. A sample of Pb-Ne collisions at
the same energy has been recently collected during the 2018 Pb-Pb run.

The first physics results, demonstrating the potential of this novel program,
have been obtained from some of the first samples collected in 2015 ane 2016.
The first measurement of charm production [24] is based on the $7.6\pm
0.5~{}\mbox{\,nb}^{-1}$ of $p$-He collisions at
$87~{}\mathrm{\,Ge\kern-1.00006ptV}$, and a few $\mbox{\,nb}^{-1}$ of $p$-Ar
data at $110~{}\mathrm{\,Ge\kern-1.00006ptV}$. Order of $10^{3}$ ${D}^{0}$
mesons and $10^{2}$ ${J\mskip-3.0mu/\mskip-2.0mu\psi\mskip 2.0mu}$ mesons are
reconstructed from both samples. These data provide the first determination of
the ${c}{\overline{c}}$ cross-section at this relatively unexplored energy
scale (see Figure 9). The rapidity dependence of the production is found to
agree with predictions based on collinear factorisation not including a
contribution from intrinsic charm, as shown in Figure 9. Therefore, these
measurements, which are sensitive to $x$ values up to $\sim 0.5$, do not
favour the large intrinsic charm contributions predicted by some models in
this kinematic range[25, 26]. The larger available datasets, and possibly the
use of a hydrogen target in the future, will allow for more accurate
constraints.

Figure 9:  In the left plot, the result for the ${c}{\overline{c}}$ production
cross-section per nucleon, obtained from $p$-He collisions at
$69~{}\mathrm{\,Ge\kern-1.00006ptV}$ [24], is compared with other experimental
results and with a NRQCD calculation [27]. The right plot shows the rapidity
dependence of the ${D}^{0}$ production cross-section in the same sample. The
result is compared with HELAC-onia calculations, not accounting for intrinsic
charm, which are rescaled to reproduce the integrated value of the measured
cross-section.

The antiproton production in the $p$-He sample at
$110~{}\mathrm{\,Ge\kern-1.00006ptV}$ [28] has also been measured. This study
is motivated by the recent precision measurements performed in space, notably
by AMS-02 [29], of the antiproton content in cosmic rays, which is sensitive
to possible exotic contributions like dark matter annihilation. For
antiprotons above 10 GeV, the largest uncertainty on the expected flux of
antiprotons from known sources, namely production in collisions between
primary cosmic rays and the interstellar medium, is due to the limited
knowledge of the corresponding production cross-sections. LHCb performed the
first $\overline{p}$ production measurement in $p$-He collisions, which are
responsible for about 40% of the expected cosmic $\overline{p}$ flux, in the
range of $\overline{p}$ momentum between $12$ and
$110~{}\mathrm{\,Ge\kern-1.00006ptV}$. The results, shown in Figure 10, are
significantly more precise than the spread among the predictions from
different phenomenological models, and are contributing to improve the models
for cosmic secondary $\overline{p}$ and the resulting constraints on dark
matter contributions [30, 31].

Figure 10: On the left plot, result for antiproton production in $p$-He
collisions at $110~{}\mathrm{\,Ge\kern-1.00006ptV}$. The differential
production cross-section $d\sigma/dp$ is shown as a function of the momentum
$p$ for different ranges of $p_{\mathrm{T}}$ and is compared with several
generators included in the CRMC package [32]. The normalization is obtained
from the $p$-$e^{-}$elastic scattering events from the same sample,
illustrated in the right plots. The upper plot shows the $p_{\mathrm{T}}$
distribution for single scattered electron candidates with negative and
positive charge, the latter sample being used to subtract the background from
hadronic collisions. After subtraction, the distribution is found to agree
well with the expected one from simulated $p$-$e^{-}$scattering, as shown in
the lower plot.

For both results, the integrated luminosity of the sample is estimated from
the yield of elastically scattered electrons from the target atoms. The yield
is determined by selecting events with a single low-pt track in the detector,
identified as an $e^{\pm}$. The background to this normalization channel is
due to soft diffractive collisions with a single track reconstructed in the
detector and is expected to be symmetric in charge, so that it can be
estimated and subtracted using the $e^{+}$ candidates, as illustrated in
Figure 10. The integrated luminosity is determined using this method with an
accuracy of 6% [28], dominated by the systematic uncertainty on the electron
reconstruction efficiency.

##  5 Conclusions and Prospects

The results obtained so far by the LHCb collaboration from heavy ion
collisions demonstrate the capability of the experiment to provide unique
contributions to this field, exploiting a variety of collision systems and
energy scales, where exclusive final states, notably in the heavy flavour
sector, can be reconstructed with high efficiency and purity.

The rich datasets collected so far provide many additional possibilities that
are being investigated. In $p$-Pb collisions, studies on direct photons, both
inclusive and in gamma-jet events, are underway, with unique sensitivity to
the saturation region [33]. The size of the sample collected in 2018 at 8 TeV
is expected to provide access to new channels, as $\chi_{c}$ and $\eta_{c}$
quarkonia states and Drell-Yan dimuon events. Studies of flow and correlations
with identified particles in the unique forward acceptance region covered by
LHCb are also planned.

Substantial development of the heavy ion program in LHCb is also expected
during the future LHC runs, taking profit of the currently ongoing detector
upgrade [34] which includes a new vertex detector with improved granularity.
The proposed plans for future heavy ion running, discussed in detail in
Reference [35], foresee an increase in integrated luminosity for $p$-Pb and
Pb-Pb collisions by more than an order of magnitude during the LHC Runs 3 and
4 (2021-2029). This will open novel possibilities, as precision studies using
Drell-Yan events and correlations in heavy flavour production. A proposal for
a second detector upgrade for the LHC Run 5 (starting 2031) has been put
forward [36, 37]. The granularity of such detector, conceived to take profit
of the full potential of LHC luminosity in $p$-$p$ collisions, would make it
possible to reconstruct even the most central Pb-Pb collisions.

The fixed-target program will also be considerably developed, thanks to the
installation of a new target device [38]. This will consist of a storage cell
containing the injected gas into a 20 cm long region located just upstream the
LHCb vertex detector. An increase in fixed-target luminosity by up to two
orders of magnitude with the same injected gas flow is anticipated.
Furthermore, it will be possible to inject more gas species, notably hydrogen,
providing $p$-$p$ reference for the other fixed-target samples. The physics
potential of such program, discussed in more detail in Reference [39], covers
a wealth of measurements that won’t be possible at any other facility in the
coming years, including heavy flavour and Drell-Yan production giving access
to nuclear PDFs at large $x$, measurements of quarkonia suppression in
different systems at $\sqrt{s_{\scriptscriptstyle\rm NN}}\sim
100~{}\mathrm{\,Ge\kern-1.00006ptV}$, and studies of particle production of
great interest for cosmic ray physics.

## References

## References

  * [1] Alves Jr A A et al. (LHCb collaboration) 2008 JINST 3 S08005
  * [2] Barschel C 2014 Ph.D. thesis RWTH Aachen U. CERN-THESIS-2013-301
  * [3] Khachatryan V et al. (CMS collaboration) 2010 JHEP 09 091 (Preprint 1009.4122)
  * [4] Aaij R et al. (LHCb collaboration) 2017 Phys. Lett. B774 159 (Preprint 1706.07122)
  * [5] Aaij R et al. (LHCb collaboration) 2018 JHEP 11 194 (Preprint 1810.07655)
  * [6] Shao H S 2016 Comput. Phys. Commun. 198 238–259 (Preprint 1507.03435)
  * [7] Lansberg J P and Shao H S 2017 Eur. Phys. J. C77 1 (Preprint 1610.05382)
  * [8] Arleo F and Peigné S 2013 JHEP 03 122 (Preprint 1212.0434)
  * [9] Ducloué B, Lappi T and Mäntysaari H 2016 Phys. Rev. D94 074031 (Preprint 1605.05680)
  * [10] Arleo F and Peigné S 2017 Phys. Rev. D95 011502 (Preprint 1512.01794)
  * [11] Aaij R et al. (LHCb collaboration) 2016 JHEP 03 133 (Preprint 1601.07878)
  * [12] Abelev B B et al. (ALICE collaboration) 2014 JHEP 12 073 (Preprint 1405.3796)
  * [13] Adam J et al. (ALICE collaboration) 2016 JHEP 06 050 (Preprint 1603.02816)
  * [14] Aaboud M et al. (ATLAS collaboration) 2018 Eur. Phys. J. C78 171 (Preprint 1709.03089)
  * [15] Sirunyan A M et al. (CMS collaboration) 2019 Phys. Lett. B790 509–532 (Preprint 1805.02248)
  * [16] Ferreiro E G and Lansberg J P 2018 JHEP 10 094 [Erratum: JHEP03,063(2019)] (Preprint 1804.04474)
  * [17] Sirunyan A M et al. (CMS collaboration) 2019 Phys. Lett. B790 270–293 (Preprint 1805.09215)
  * [18] Aaij R et al. (LHCb collaboration) 2017 JHEP 10 090 (Preprint 1707.02750)
  * [19] Aaij R et al. (LHCb collaboration) 2019 JHEP 02 102 (Preprint 1809.01404)
  * [20] Acharya S et al. (ALICE collaboration) 2018 JHEP 04 108 (Preprint 1712.09581)
  * [21] Aaij R et al. (LHCb collaboration) 2019 Phys. Rev. D99 052011 (Preprint 1902.05599)
  * [22] LHCb collaboration 2018 LHCb–CONF–2018–003
  * [23] Klein S R, Nystrand J, Seger J, Gorbunov Y and Butterworth J 2017 Comput. Phys. Commun. 212 258–268 (Preprint 1607.03838)
  * [24] Aaij R et al. (LHCb collaboration) 2018 Phys. Rev. Lett. 122 132002 (Preprint 1810.07907)
  * [25] Pumplin J, Lai H L and Tung W K 2007 Phys. Rev. D75 054029 (Preprint hep-ph/0701220)
  * [26] Brodsky S J, Kusina A, Lyonnet F, Schienbein I, Spiesberger H and Vogt R 2015 Adv. High Energy Phys. 2015 231547 (Preprint 1504.06287)
  * [27] Maltoni F et al. 2006 Phys. Lett. B638 202–208 (Preprint hep-ph/0601203)
  * [28] Aaij R et al. (LHCb collaboration) 2018 Phys. Rev. Lett. 121 222001 (Preprint 1808.06127)
  * [29] Kounine A 2012 Int. J. Mod. Phys. E21 1230005
  * [30] Reinert A and Winkler M W 2018 JCAP 1801 055 (Preprint 1712.00002)
  * [31] Korsmeier M, Donato F and Di Mauro M 2018 Phys. Rev. D 97(10) 103019 (Preprint 1802.03030)
  * [32] Pierog T, Baus C and Ulrich R CRMC (Cosmic Ray Monte Carlo package) web.ikp.kit.edu/rulrich/crmc.html
  * [33] Da Silva C (LHCb collaboration) 2018 Contribution to the Quark Matter 2018 conference URL https://indico.cern.ch/event/656452/contributions/2859624
  * [34] LHCb collaboration 2012
  * [35] Citron Z et al. 2018 HL/HE-LHC Workshop Report (Preprint 1812.06772)
  * [36] LHCb collaboration 2017 CERN–LHCC–2017–003
  * [37] LHCb collaboration 2018 CERN–LHCC–2018–027 LHCb–PUB–2018–009 (Preprint 1808.08865)
  * [38] Di Nezza P et al. 2018 Tech. Rep. CERN-PBC-Notes-2018-007
  * [39] Graziani G et al. 2018 Tech. Rep. LHCb-PUB-2018-015


11institutetext:  University of Zagreb, Faculty of Geodesy, Hvar Observatory,
Zagreb, Croatia
11email<EMAIL_ADDRESS>22institutetext: University of Zagreb, Faculty of
transport and traffic sciences, Zagreb, Croatia 33institutetext: Swedish
Institute of Space Physics, Uppsala, Sweden 44institutetext: Christian-
Albrechts University in Kiel, Department of Extraterrestrial Physics, Kiel,
Germany

# A new method of measuring Forbush decreases

M. Dumbović 11 L. Kramarić 1122 I. Benko 1133 B. Heber 44 B. Vršnak 11

###### Abstract

Context. Forbush decreases (FDs) are short-term depressions in the galactic
cosmic ray (GCR) flux and one of the common signatures of coronal mass
ejections (CMEs) in the heliosphere. They often show a two-step profile, the
second one associated with the CME’s magnetic structure. This second step can
be described by the recently developed analytical FD model for flux ropes
(FRs) –ForbMod.

Aims. The aim of this study is to utilise ForbMod to develop a best-fit
procedure to be applied on FR-related FDs as a convenient measurement tool.
Our motivation is to develop a best-fit procedure that can be applied to a
data series from an arbitrary detector. Thus, the basic procedure would
facilitate measurement estimation of the magnitude of the FR-related FD, with
the possibility of being adapted for the energy response of a specific
detector for a more advanced analysis.

Methods.  The non-linear fitting was performed by calculating all possible
ForbMod curves constrained within the FR borders to the designated dataset and
minimising the mean square error (MSE). In order to evaluate the performance
of the ForbMod best-fit procedure, we used synthetic measurements produced by
calculating the theoretical ForbMod curve for a specific example CME and then
applying various effects to the data to mimic the imperfection of the real
measurements. We also tested the ForbMod best-fit function on the real data,
measured by detector F of the SOHO/EPHIN instrument on a sample containing 30
events, all of which have a distinct FD corresponding to the magnetic obstacle
(MO). The extraction of FD profiles (from the onset to the end) was performed
manually by an observer, whereby we applied two different versions of border
selection and assigned a quality index to each event.

Results. We do not find notable differences between events marked by a
different quality index. For events with a selection of two different borders,
we find that the best fit applied on extended interplanetary coronal mass
ejection (ICME) structure borders results in a slightly larger MSE and
differences compared to the traditional method due to a larger scatter of the
data points. We find that the best-fit results can visually be categorised
into six different FD profile types. Although some profiles do not show a
visually pleasing FD, the ForbMod best-fit function still manages to find a
solution with a relatively small MSE.

Conclusions. Overall, we find that the ForbMod best-fit procedure performs
similar to the traditional algorithm-based observational method, but with
slightly smaller values for the FD amplitude, as it’s taking into account the
noise in the data. Furthermore, we find that the best-fit procedure has an
advantage compared to the traditional method as it can estimate the FD
amplitude even when there is a data gap at the onset of the FD.

###### Key Words.:

Sun: coronal mass ejections (CMEs) – solar-terrestrial relations – Sun:
heliosphere – galactic cosmic rays

##  1 Introduction

Forbush decreases (FDs) are short-term depressions in the galactic cosmic ray
(GCR) flux that result from the interaction of GCRs with coronal mass
ejections (CMEs). They are thus one of the common signatures of CMEs in the
heliosphere (see overviews by e.g. Lockwood, 1971; Cane, 2000; Belov, 2009).
Detecting CMEs at various different locations in the heliosphere is crucial
for the analysis of their evolution, as well as the validation of the CME
propagation and evolution models. They are typically observed by in situ
detectors that measure magnetic field and plasma properties (e.g. Zurbuchen &
Richardson, 2006; Kilpua et al., 2017); however, these are not available as
frequently as particle detectors that can detect FDs. Therefore, a better
understanding of FDs can lead to CME analysis where other types of observation
are not possible.

Figure 1: Sketch of the ForbMod model for FR-related FDs. The GCRs diffuse
into the self-similarly expanding FR (r/a=const.), which is initially empty at
its centre (r/a=0), where a(t) is the FR radius and r(t) the radial distance
from the centre of the FR to its borders. As a result, at an arbitrary time
after the FR eruption from the Sun (i.e. at an arbitrary heliocentric
distance), a distribution of GCRs inside the FR is observed, given by the
Bessel function. The relative GCR density (FD[%]) at the borders of the FR is
0%, whereas in the centre of the FR it is given by the amplitude of the
depression ($A_{\mathrm{FD,max}}$). For model details, see the main text.

FDs caused by interplanetary coronal mass ejections (ICMEs) often show a two-
step profile, one associated with the turbulent and compressed shock-sheath
region ahead of the CME magnetic structure and the other with the CME magnetic
structure itself (Cane, 2000). Because of their different physical properties,
the mechanism through which these two regions produce the depression is
different, and thus should be modelled differently. The current paradigm is
that the CME magnetic structure is a flux rope (FR), a structure with field
lines helicoidally winding around a central axis (e.g. Zhang et al., 2021). It
was first proposed by Morrison (1956) that FDs might be caused by clouds of
closed magnetised plasma, whereby they are initially empty of particles and
slowly fill as they propagate through the interplanetary space. This approach
has since been utilised to explain FR-associated FDs, where the particles are
allowed to enter via “cross-field” transport; in other words, perpendicular
diffusion (e.g. Cane et al., 1995; Quenby et al., 2008). In addition to the
particle diffusion, another important contribution comes from the expansion of
the FR. Namely, as they propagate through the IP space, CMEs expand due to the
pressure imbalance (e.g. Klein & Burlaga, 1982; Démoulin & Dasso, 2009),
where, consequently, as the size of the magnetic structure increases, its
magnetic field weakens (Bothmer & Schwenn, 1998; Démoulin et al., 2008; Vršnak
et al., 2019). It was first proposed by Laster et al. (1962) and Singer et al.
(1962) that expansion is needed to explain the observational properties of
FDs. More recently, the diffusion-expansion approach was applied by Munakata
et al. (2006) and Kuwabara et al. (2009) in a numerical model and by
Subramanian et al. (2009) and Arunbabu et al. (2013) in an analytical model,
considering isotropic self-similar expansion. In the recent FD analytical
model ForbMod (Dumbović et al., 2018), different types of expansion are
considered for the expansion-diffusion approach, constrained by the
observational ranges of the expansion factors (Gulisano et al., 2012). ForbMod
was advanced to take into account the energy dependence of the detector when
comparing the model results to measurements (Dumbović et al., 2020), since FDs
are energy-dependent phenomena (i.e. their amplitudes depend on the energy
range and cut-off of the detector which detects them Cane, 2000). This
facilitates direct comparison of modelled and observed FDs, as was shown by,
for example, Freiherr von Forstner et al. (2021).

It should be noted that there are other approaches to FD modelling in the FR,
such as perpendicular transport by guiding centre drifts (e.g. Krittinatham &
Ruffolo, 2009; Tortermpun et al., 2018), full trajectory integration using CME
FR-type models (Petukhova et al., 2019), or CME magnetic field reconstructions
from in situ measurements (Benella et al., 2020), as well as describing FDs
via the change in the single GCR spectrum modulation parameter attributed with
a CME (Guo et al., 2020). However, a recent study by Laitinen & Dalla (2021),
in which they performed full-orbit particle simulations, with the interface
between the external interplanetary magnetic field and the isolated FR field
lines modelled analytically, revealed that the transport around the FR border
is fast compared to diffusive radial propagation within the FR. As a result,
the propagation of GCRs into the FR can be modelled as diffusion into a
cylinder. Moreover, the diffusive approach is further supported by most recent
observations (see e.g. recent studies by Janvier et al., 2021; Davies et al.,
2023).

The aim of this study is to utilise ForbMod to develop a best-fit procedure
that can be applied to FR-related FDs as a convenient measurement tool. Our
motivation is to develop a best-fit procedure that can be applied to a data
series from an arbitrary detector. Thus, the basic procedure would facilitate
measurement estimation of the onset, duration, and magnitude of the FR-related
FD, with the possibility of being adapted for the energy response of a
specific detector for a more advanced analysis. It should be noted that in
many cases the onset of the FR-related FD does not necessarily correspond to
the onset of the FD itself, as the latter might correspond to the preceding
shock-sheath region. In such cases we refer to the onset of the second
decrease as the onset of FD (or FR-related FD). In section 2 we describe the
mathematical formalism of the best-fit procedure as well as two samples used
in the study to test its performance. We first used synthetic measurements to
test the performance of the best-fit procedure under controlled conditions.
Next we used a real sample of events detected by the F-detector of the
Electron Proton Helium INstrument instrument aboard the SOHO spacecraft
(SOHO/EPHIN Müller-Mellin et al., 1995). The results are presented and
discussed in section 3.

##  2 Methods and data

###  2.1 ForbMod best-fit procedure

ForbMod is an analytical diffusion-expansion model that assumes that the CME
magnetic structure (i.e. FR) can locally be represented by a straight
cylinder, which is initially empty at its centre. The FR is assumed to be
propagating at a constant speed, expanding self-similarly, and the particles
can only enter it via perpendicular diffusion. It is assumed that GCRs are
mostly protons. The time-dependent transport equation with these assumptions
simplifies to a time-dependent radial diffusion of protons into an infinite
straight cylinder. The corresponding solution is a combination of the radial-
dependent solution, which is given by the cylindrical Bessel function of the
zero-order, and the time-dependent solution, which is given by the exponential
function (for details on the model see Dumbović et al., 2018). The amplitude
of the FD, defined as normalised to border conditions, depends on the
diffusion rate of particles into the FR. The diffusion rate is governed by the
behaviour of the magnetic field, which is generally decreasing due to
expansion. Thus, the FD amplitude is expected to decrease with diffusion time
(i.e. helisopheric distance). The sketch visualising ForbMod is given in
Figure 1. The diffusion rate of particles also depends on the particle energy,
as higher-energy particles enter the FR more easily. Therefore, it is expected
that the FD amplitude is energy-dependent and that it should appear
differently in detectors with a different energy response (for details see
Dumbović et al., 2018, 2020).

Figure 2: Effects applied to the data to mimic the “imperfection” of the real
measurements: Y-noise (upper panels), X-noise (middle panels), and asymmetry
(lower panels) for three different levels of noise and asymmetry,
respectively. The original theoretical curve is given in black, whereas the
data with the applied effect are coloured blue. For details, see the main
text.

Figure 3: Effects applied to the data to mimic the “imperfection” of the real
measurements: data gaps at certain positions (upper panels), data gaps of a
certain size (middle panels), and different data resolutions (lower panels).
The original theoretical curve is given in black, whereas the data with the
applied effect are coloured blue. For details, see the main text.

Figure 4: In situ measurements of four different events assigned with four
different quality indices, QI4–QI1, from top left to bottom right. Top-to-
bottom panels for each event show: 1) magnetic field strength (black, in nT)
and its fluctuations (grey, in nT); 2) magnetic field GSE components, x, y,
and z, coloured red, blue, and green, respectively (in nT); 3) plasma density
(black, in $\mathrm{cm}^{-3}$), temperature (red, in 1e5 K), and expected
temperature (blue, in 1e5 K); 4) plasma flow speed (black, in km/s) and beta
(grey, non-dimensional); 5) the SOHO/EPHIN F-detector particle count (in %).
The shaded region represents the magnetic obstacle, where the borders were
selected according to the low-beta condition.

Figure 5: Two different versions of border selection: inner structure (yellow
shading) and extended ICME signatures (grey shading). For details, see the
main text.

Figure 6: Relative difference of the best-fit FD amplitude to the FD amplitude
of the theoretical FD vs goodness of fit for different levels of data
imperfection effects: a) Y-noise, b) X-noise, c) asymmetry, d) data gap
position, e) data gap size, and f) resolution (for details, see the main
text).

We developed the ForbMod best-fit procedure by applying a non-linear fitting
to the data. For that purpose we first used synthetic measurements, as is
described in Section 2.2. In general, the non-linear fitting was performed by
calculating all possible ForbMod curves constrained within FR borders to the
designated dataset and minimising the mean square error, $MSE$:

$MSE(y,\hat{y})=\frac{1}{n_{\mathrm{SAMPLES}}}\sum_{i=1}^{n_{\mathrm{SAMPLES}}-1}(y_{i}-\hat{y}_{i})^{2}\,,$
(1)

where $\hat{y}_{i}$ is the predicted value of the $i$-th sample, $y_{i}$ is
the corresponding true value, and there are $n_{\mathrm{SAMPLES}}$ samples.
The ForbMod best-fit procedure actually performs non-linear fitting of the
Bessel function of the zero order (i.e. the radial part of the ForbMod
solution) to the designated time series. The input for the procedure is the
designated dataset with defined FR borders (i.e. the dataset translated into
the FR coordinate system radial direction, with the centre of the FR at
$\hat{r}=0$ and FR borders at $\hat{r}=-1$ and $\hat{r}=1$, respectively;
$\hat{r}=r/a$, where $a$ is the radius of the FR). The output of the procedure
is the best-fitted amplitude of the measured FD, $A_{\mathrm{FR,max}}$, and
the corresponding goodness of fit, that is, $MSE$ (as defined by equation 1).
The FD profiles needed to be translated from a timescale into the spatial
scale of the FR in order for us to apply the ForbMod best-fit. For that
purpose, we assumed that the spacecraft crosses the full diameter of the FR,
which has a circular cross-section. Furthermore, we assumed that the
spacecraft trajectory is oriented vertically with respect to the FR axis. We
note that with this assumption the model might not realistically capture FRs
with orientations deviating significantly from low or high inclinations.
Nevertheless, with this assumption, due to self-similar expansion, the
spacecraft observes deceleration of plasma flow across the FR. Thus, the
diameter can be given by the equation for accelerated linear motion,
$s=0.5at^{2}+v_{0}t$. The borders of the FR were set at [-1,1], per the
definition of the ForbMod best-fit function.

###  2.2 Synthetic Forbush decrease measurements

In order to evaluate the performance of the ForbMod best-fit procedure, we
used synthetic measurements. Synthetic measurements were produced by
calculating the theoretical ForbMod curve for a specific example CME and then
applying various effects to the data to mimic the “imperfection” of the real
measurements.

The arbitrary CME example used to produce synthetic measurements (loosely
based on event analysed by Freiherr von Forstner et al., 2021) had the
following inputs:

  * •

detector response function for perfect detector: $Y(E)=1$

  * •

detector cutoff: $E_{\mathrm{cutoff}}=0.05$GeV

  * •

event month and year: April 2020

  * •

ICME transit time in hours (i.e. diffusion time), $TT=123.38$h

  * •

starting distance from the Sun, $R_{0}=9.54\mathrm{R_{SUN}}$

  * •

final distance from the Sun, $R_{fin}=1.005$au

  * •

initial FR radius, $a_{0}=1.81\mathrm{R_{SUN}}$

  * •

expansion factor for FR radius, $n_{a}=1$

  * •

in situ magnetic field strength, $B=16.2$nT

  * •

expansion factor for mag. field, $n_{B}=2$

The effects applied to the data to mimic the “imperfection” of the real
measurements were:

  * •

Y-noise: \- random noise in the y direction

  * •

X-noise: \- random noise in the x direction

  * •

asymmetry: \- systematic left-shift of the FD minimum in the x direction (to
simulate the asymmetric FD profile)

  * •

data gaps: \- simulating a data gap by removing certain data points around the
selected position within the FD

  * •

low data resolution: \- reducing the number of data points in the FD profile

The Y-noise value was calculated and added to each data point of the ForbMod
solution in the following way:

${y_{i,\mathrm{SYNTH}}=y_{i,\mathrm{ForbMod}}+\mathrm{d}y_{i}\,\,\,;\,\,\,\mathrm{d}y_{i}=A_{\mathrm{norm}}\cdot\delta_{n}\,,}$
(2)

where $y_{i,\mathrm{SYNTH}}$ is the synthetic measurement,
$y_{i,\mathrm{ForbMod}}$ is the ForbMod solution data point, $\delta y_{i}$ is
the noise, and $A_{\mathrm{norm}}$ is the randomly selected value from the
normal distribution defined by the ForbMod solution amplitude as the mean.
$A_{\mathrm{norm}}$ was determined using the Python numpy package function
random.normal. $\delta_{n}$ was an arbitrary percentage applied to the
randomly selected $A_{\mathrm{norm}}$. The parameter $\delta_{n}$ could be
varied in order to obtain different percentages (i.e. different levels) of
noise. The X-noise value was calculated and applied in the same manner as the
Y-noise, only in the x direction. The Y-noise and X-noise of different
percentages are shown in the upper and middle panels of Figure 2,
respectively.

The asymmetry was introduced by incrementally shifting positive $r$ values
corresponding to data points of the ForbMod solution:

$\delta r_{i,\mathrm{pos}}=r_{i}+i\cdot\delta_{a}\,,$ (3)

where $r_{i}$ is the positive $r$ value corresponding to the $i$-th data point
of the ForbMod solution (counted for positive $r$ values only) and
$\delta_{a}$ is a shift percentage between two iterations, which could be
varied in order to analyse different percentages (i.e. levels) of asymmetry.
After applying the asymmetry, the dataset was re-scaled to be defined on a
[-1,1] interval using the sklearn package function minmax_scale. The asymmetry
of different percentages is shown in the lower panels of Figure 2.

The data gaps were introduced by removing a certain number of data points from
a certain location within the FD profile. Thus, we could vary both the
location of the data gap, as well as its size. This is shown in the example in
Figure 3. In the upper panels, a gap 10% in size (i.e. 200 data points out of
2000) is consecutively shifted from the start of the FD to its minimum (five
different positions). In the middle panels of Figure 3 various gap sizes are
shown, where the gap is centred around the middle of the FD main phase (i.e.
the decrease phase). Finally, as the last effect applied to the data to mimic
the “imperfection” of the real measurements, we considered data resolution, as
is shown in the bottom panels of Figure 3.

###  2.3 Event sample and measurements

A sample of events was created by associating SOHO/EPHIN detected FDs with in
situ detected ICMEs and remotely observed CMEs. We used a list of SOHO/EPHIN
detected FDs compiled by (Belov et al., 2021) and a catalogue of near-Earth
ICMEs compiled by (Richardson & Cane, 2010) and regularly updated at
CALTECH111http://www.srl.caltech.edu/ACE/ASC/DATA/level3/icmetable2.htm. Since
we aimed to use the sample for the follow-up analysis, where a 3D
reconstruction using stereoscopic data to obtain CME remote properties would
be performed, we constrained our analysis to the time period when at least one
STEREO spacecraft was operational, where a reliable CME-ICME association could
be made, and where CME signatures in white light coronagraphs were significant
enough to perform a 3D reconstruction. The sample contains 30 events in the
time period 2007–2019. In this study, we only used the part of the sample
containing the interplanetary in situ signatures of ICMEs and the
corresponding FD, and thus we only provided a description of the ICME-FD
sample in detail, whereas the description of the CME-part of the sample will
be given in a future study.

Interplanetary coronal mass ejections are roughly simultaneous with FD events
(Cane et al., 1996; Dumbović et al., 2011), and therefore this association was
trivial and performed based on the timing provided by Belov et al. (2021) or
Richardson & Cane (2010)’s catalogues. To analyse the in situ properties, we
used one-minute plasma and magnetic field data in the geocentric solar
ecliptic (GSE) system provided by the
OMNIWeb222https://omniweb.gsfc.nasa.gov/html/ow_data.html#norm_pla database
(King & Papitashvili, 2005). We identified ICME signatures using criteria
described by (e.g. Zurbuchen & Richardson, 2006; Kilpua et al., 2017),
considering the notification of (Rouillard, 2011), according to which an ICME
consists of the compressed, heated, and turbulent shock-sheath region,
followed by a magnetic cloud or ejecta. The magnetic cloud is characterised by
a low-temperature, low-plasma beta parameter, an enhanced smoothly rotating
magnetic field, and an expanding flow speed profile, whereas the term magnetic
ejecta is used when these criteria are only partially met. We note that both
magnetic cloud and ejecta are nowadays broadly considered as magnetic
obstacles (MO, Jian et al., 2006; Nieves-Chinchilla et al., 2018), and
therefore we use this more general term.

Example events are given in Figure 4. For each event we analysed: the magnetic
field strength, $B$, and its fluctuations, $dBrms$333the root mean square
variation in the vector of the interplanetary magnetic field; the magnetic
field components, $B_{x}$, $B_{y}$, and $B_{z}$; the plasma density,
temperature, and expected temperature444empirical formula based on the speed
of the solar wind calculated according to Lopez (1987) and Richardson & Cane
(1995); the plasma flow speed and plasma beta parameter; and the SOHO/EPHIN
F-detector particle count, $CR\,\,count$555normalised to the pre-decrease
value for each event. We used hourly averaged detector F data of the
SOHO/EPHIN, as they were shown to be suitable for a FD analysis, especially
the MO part (Heber et al., 2015; Dumbović et al., 2020; Belov et al., 2021).
The SOHO/EPHIN data were time-shifted to the Earth’s bow shock using the
OMNIWeb formula for time-shift:

${\Delta t=\frac{X}{V}\cdot\frac{1+\frac{Y\cdot W}{X}}{1-\frac{V_{e}\cdot
W}{V}},}$ (4)

where $X$ and $Y$ are GSE X and Y components of the spacecraft position vector
(we used daily positions of SOHO), $V$ is the observed solar wind speed
(assumed to be radial), $V_{e}$ is the speed of the Earth’s orbital motion (30
km/s), and $W$=tan(0.5atan ($V$/428)) is the geometry factor (for details, see
the OMNIWeb webpage666https://omniweb.gsfc.nasa.gov/html/ow_data.html.

In order to achieve better visibility of smaller FD variations, SEP events
were removed from the data using the $<600$ counts threshold on the
B-detector, which is sensitive to lower energy particles, that is, sharp
increases due to SEPs. We visually inspected each event and discarded events
where either no MO or corresponding FD was observed. In addition, we
classified events according to the quality index (QI), subjectively determined
by the observer, in the range 1–4, with 1 being the worst quality and 4 being
the best quality (see Figure 4). The QI is based on criteria regarding the
clarity of MO signatures and borders, as well as the clarity of FD signatures
and borders. Both are less clear for events designated as having a lower QI.
We note that in QI determination, some more weight was given to the clarity of
an onset corresponding to both MO and FD, than to one corresponding to other
MO/FD properties, due to its importance for the best-fit application. We also
note that although we analysed the clarity of MO borders visually, we did not
actually determine or measure MO borders, but FD borders.

It can be seen in Figure 4 that the MO signatures are weakest for the QI1
event (bottom right). The magnetic field is increased and the $B_{z}$ is
rotating, although the magnetic field is not very smooth. Plasma beta is low,
although the temperature is only partially lower than expected and the
expansing profile of the flow speed is missing. The onset of the FD is clear;
however, the recovery of the FD seems to be missing. This is likely related to
another structure that can be observed trailing the MO. In the QI2 event
(bottom left) the MO signatures are clearer. The magnetic field is again
increased and the $B_{y}$ is rotating smoothly. Plasma beta is low, but the
temperature is again only partially lower than expected. The flow speed seems
to show a globally expanding profile, but with local irregularities. The
corresponding FD has a clear onset; however, the onset is slightly shifted in
time with respect to the MO onset (according to the low-beta condition).
Moreover, the FD shows a double depression without the full recovery to the
pre-decrease level. We note that the substructuring of the FD corresponds to
the substructuring of the MO. The first depression seems to correspond to the
rotation of the Bz component from positive to negative values and a small
density increase, whereas the second depression seems to correspond to the
rotation of the By component from negative to positive values and the low
density region. This substructuring might indicate a complex inner structure
of the MO, or the passage of the spacecraft under some oblique angle. However,
although they both show substructuring, it appears the global structure may
correspond to that of a single MO, as was found for example by Nieves-
Chinchilla et al. (2018). The QI3 event shows clear MO signatures as, in
addition to other properties that the QI2 event shows, the flow speed in the
QI3 event shows a clearly expanding profile. Although a depression can be seen
to start with the shock-sheath, the second depression of the FD has a clear
onset that corresponds to the MO with a clear tendency to return to the pre-
decrease level, although it does not (most likely due to the high-speed stream
that follows). Moreover, the FD seems to have different recovery rate at the
back (which again corresponds to the substructuring of the MO, similar to the
QI2 event, but much less pronounced). Finally, the QI4 event shows clear MO
signatures and a clear, symmetric FD with borders that roughly match the FR
borders. However, we do note that although globally this event shows clear MO
signatures, again some substructuring can be observed. To summarise, in
general the best-quality events have nicely observed FD and MO signatures and
clear borders. The worst-quality events have unclear borders, pronounced
substructures, and/or an unclear recovery. Out of the 30 events in the sample,
6, 12, 8, and 4 have QIs 1, 2, 3, and 4, respectively. The full sample of
events is given in Table 2 in the Appendix A.

###  2.4 Extraction of FD profiles: The border selection problem

In order to apply the best-fit procedure, the FD profiles needed to be
extracted, that is, the borders of the FD needed to be selected. The onset of
FD is usually defined as the point where the GCR count starts to drop, whereas
the end of FD, in a textbook example, would be the point where it returns to
its pre-decrease level (e.g. Dumbović et al., 2011). In textbook examples, the
recovery of two-step FD is prolonged after the passage of the ICME due to the
“shading effect” of the shock-sheath region, whereas for a FR, FD (which can
be self-standing or a second decrease in a two-step FD), the end of FD
corresponds to the end of the FR (Dumbović et al., 2020). However, in many
non-textbook cases, the determination of the onset and end of FD is not a
trivial procedure that can be defined by a simple algorithm applicable evenly
to all events, as can clearly be seen from Figure 4. Therefore, the extraction
of FD profiles (from the onset to the end) was performed manually by an
observer. This may result in a certain subjectivity when determining the
borders, especially the end of FD. Therefore, we applied two different
versions of border selection, as is shown in Figure 5:

  * •

extended ICME borders \- these encompass typical ICME signatures such as, for
example, low temperature or low plasma beta, but with borders set based on the
GCR depression onset or end and/or based on the total magnetic field increase
onset or end. They thus overlap with typical MO signatures only in a subset of
events and they often contain external disturbed parts of the FR.

  * •

inner structure borders \- we find that many events show clearer MO signatures
in a sub-interval of what would normally be selected as the global MO
structure according to, for example, the low-beta condition. Therefore, first
we selected borders to encompass this inner substructure of the CME, which
shows clear MO signatures. However, we did this only for events where
corresponding clear FD substructure borders could also be found.

We find that the inner structure can be identified in 24 events. The inner
structure is usually small compared to the typical ICME duration, which is
roughly 24 hours (Richardson & Cane, 2010). Consequently, these events have
fewer data points to be fitted, but also a more compact profile, without
substructures. We did not further refine the QI based on the different border
selection, but instead treated the border selection and QI as independent
classifications. In some events the inner structure could not be identified,
whereas in some other events the whole duration of the event could be
considered as the inner structure. We found four such events, which we in turn
regard as not having extended ICME borders. The best-fit procedure is applied
to the GCR measurements within the selected borders, and thus slightly
different results of best-fit FD amplitude, FD_bf, can be expected for
different border selections. In addition, for each event we determined the FD
amplitude using a traditional algorithm-based observational method, FD_obs,
whereby FD was measured from the onset to the minimum. Therefore, this
observed FD, FD_obs, may also be slightly different for different border
selections. The borders corresponding to the inner structure and extended ICME
are given in Table 2 in the Appendix A for each event, along with the
corresponding observed and best-fitted FD.

##  3 Results

###  3.1 ForbMod best-fit function testing through synthetic measurements

Figure 7: Relative difference of the best-fit FD amplitude to the FD amplitude
of the theoretical FD vs $MSE$ for a mixture of different levels of data
imperfection effects. In the left plot, no data refinement is included,
whereas in the right plot data with extreme deviations from the identity line
have been removed (for details, see the main text). The identity line
($y=100*x$) is marked as the solid black line in both figures, whereas the
solid red line shows the linear regression line (linear regression equation
and squared Pearson’s coefficient are given in red in both plots).

We applied the ForbMod best-fit procedure to the synthetic data in order to
check how the goodness of fit ($MSE$) behaves for each of the different
effects, but also for different levels of a certain effect. For each best-fit
to the synthetic data we calculated $MSE$ and the relative difference of the
best-fit curve minimum to the minimum of the theoretical curve (i.e. the
relative difference of the best-fit FD amplitude to the FD amplitude of the
theoretical FD, $\mathrm{d}A/A$). We first applied the ForbMod best-fit
procedure to synthetic data with different levels of Y-noise and X-noise. The
results are shown in Figure 6 a and b. We can see that, for the same level of
Y-noise, $\mathrm{d}A/A$ and $MSE$ are correlated, whereas for the same level
of X-noise the $\mathrm{d}A/A$ remains constant. Therefore, by increasing the
X-noise we increase both $MSE$ and $\mathrm{d}A/A$. However, by increasing the
Y-noise we only increase $MSE$, and $\mathrm{d}A/A$ does not change
drastically. This is related to the fact that Y-noise is applied in the same
direction in which the $MSE$ is calculated. Simply put, the random noise up
and down around the curve will contribute to the $MSE$, but the average curve
will stay the same as the original one.

We next applied the ForbMod best-fit procedure to synthetic data with
different levels of X-noise and asymmetry. This is shown in Figure 6c. We can
see that for lower values of asymmetry ($<0.15\%$) asymmetry works in the same
direction as the X-noise: the larger the $MSE$, the larger the
$\mathrm{d}A/A$. However, for large asymmetry there seems to be a saturation
of $\mathrm{d}A/A$, which remains scattered around 60% regardless of the
$MSE$. We note that combining three different data imperfection effects
(Y-noise, X-noise, and asymmetry) very quickly leads to large $\mathrm{d}A/A$,
as the data does not resemble the FD anymore. Therefore, it is already
impossible to apply the best-fit procedure to data with Y-noise, X-noise, and
asymmetry greater than 15%, 15%, and 0.15%, respectively.

We next applied the data gap to the synthetic measurements and repeated the
best-fit procedure. We applied a data gap 10% in size to five different
positions, starting from the onset to the minimum, as is described in Section
2.2. The results are shown in Figure 6d), where we find that for a data gap
position closer to the FD minimum, the increase in $MSE$ results in an
increase in $\mathrm{d}A/A$ greater than for data gaps closer to the onset of
FD. In other words, we might expect to fit worse when our dataset has a gap
close to the FD minimum. We also checked how the data size influences the
best-fit procedure (Figure 6e). We find that for gap sizes of up to roughly
35% the scatter of the $MSE$ for the same $\mathrm{d}A/A$ increases -
depending on where the data gap is and what other noise the data has, the data
gap of the same size can lead to a range of different $\mathrm{d}A/A$, and
this range increases with the gap size. However, for gap sizes $>35\%$ we
observe the opposite trend. This is related to the fact that with an increase
in the data gap we removed data and with it the other data imperfection
effects.

Finally, we applied different combinations of data imperfection effects using
three different resolutions of the original theoretical curve data (2000, 200,
and 20 data points). In Figure 6f we can see that for smaller resolutions we
expect a larger scatter of $MSE$ for the same $\mathrm{d}A/A$.

Figure 8: Synthetic data and the best-fit results for several selected
examples. Different examples correspond to different combinations of data
imperfection effects, where the set of values of the data imperfection effects
is different for different rows and varies from high to low and top to bottom,
respectively: 1) top: Y-noise=30%, X-noise=30%, asym=0.15%, gap position=3,
gap size=30%, data resolution=20pts; 2) middle: Y-noise=20%, X-noise=20%,
asym=0.10%, gap position=2, gap size=20%, data resolution=20pts; 3) bottom:
Y-noise=5%, X-noise=5%, asym=0.05%, gap position=2, gap size=10%, data
resolution=20pts. The right-most plot in each row shows the scatter of
different $\mathrm{d}A/A$ and $MSE$ values obtained through 1000 runs for the
same maximum data imperfection effect levels. Different combinations of data
imperfection effects may lead to different $MSE$ and $\mathrm{d}A/A$
behaviour. The three examples shown left of the corresponding scatter plot are
highlighted in a colour corresponding to the colour of the synthetic data for
the given event (green for the left-most, yellow for the middle, and red for
the right-most examples). Blue data points correspond to the original
theoretical curve, and the best-fit function is given by the black curve.

In Figure 7 we show the relative difference of the best-fit FD amplitude to
the FD amplitude of the theoretical FD versus the $MSE$ for a mixture of
different levels of data imperfection effects. Combinations of different
imperfection effects can quickly lead to very high $MSE$ and $\mathrm{d}A/A$,
and therefore we constrained combinations to the following levels:
Y-noise=15%, X-noise=15%, asym=0.15%, and gap size=30%. However, we allowed
parameters to go beyond these levels when one or more of the data imperfection
parameters was kept at level 0. We can see in Figure 7 that the linear fit is
not far away from the identity line ($y=100*x$); however, the scatter of the
data is quite large.

The linear fit moves even closer to the identity line when we refine the
sample by removing extremely outlying points, that is, those that have
$\mathrm{d}A/A>20\%$ away from the identity line and extremely low or high
corresponding $MSE$, that is, $<0.2$/$>0.3$, respectively. This is shown in
the right plot of Figure 7. Therefore, we next analysed visually the synthetic
data and the best-fit function, as well as the $MSE$ and $\mathrm{d}A/A$ for
cases which are that to the identity line and those far away. This is shown
for nine selected examples in Figure 8. The examples are given in the nine
panels on the left of the Figure, whereas the right-most plots show the
scatter plots of $\mathrm{d}A/A$ vs. $MSE$ for three different sets of
randomly produced synthetic data. The scatter plots differ slightly due to
randomness of producing synthetic data, but in general show similar behaviour.
The three examples shown left of the corresponding scatter plot are
highlighted in a colour corresponding to the colour of the synthetic data for
the given event (green for the left-most, yellow for the middle, and red for
the right-most examples). For synthetic data of the same colour, top-to-bottom
panels show different levels of data imperfection effects, ranging from
highest (upper plots) to lowest (bottom plots).

We can immediately notice that the level of data imperfection affects the
visual representation of the FD. The top plots showing very high levels of
data imperfection effects have a barely identifiable FD, whereas in the bottom
plots the data depict an almost ideal symmetric FD profile. This visual
representation is in principal reflected by the $MSE$ – lower plots have much
smaller $MSE$ compared to the upper plots. The middle plots present
intermediate cases between extreme examples. Looking at the difference between
the theoretical and best-fit curve, $\mathrm{d}A/A$, we can see that it is
negligible for the bottom left-most plot, whereas in the upper left-most plot
it almost reaches 10%. However, for the top-to-bottom right-most examples it
is reversed – $\mathrm{d}A/A$ is small for the top right example, whereas in
the bottom right-most example it almost reaches 10%. This reflects the
previously observed fact that for some data imperfection effects $MSE$ is not
correlated with $\mathrm{d}A/A$. These are the cases that are far away from
the identity line. We note that in the example of low levels of the data
imperfection effects (bottom plots) $MSE$ is low ($<0.03$) in all three cases
and $\mathrm{d}A/A$ is within 10%. Therefore, we do not regard any of the
bottom plots cases as being far away from the identity line.

From this analysis we see that, if we exclude events for which visual
representation is not clear enough to identify FD, we are left with events
close to the identity line. If we further restrict ourselves only to events
with $MSE<0.1$, we are left with events with $\mathrm{d}A/A<$10%. Therefore we
can conclude that $MSE<0.1$ can serve as a goodness-of-fit measure resulting
in $<10\%$ error of the fitted amplitude under the following assumptions: 1)
the real data behave as synthetic data; 2) the data in general depict the
behaviour of the theoretical ForbMod FD (with some levels of data noise); 3)
FD is clearly visible in the data and not masked by noise.

###  3.2 ForbMod best-fit function application to SOHO/EPHIN measurements

Figure 9: Six types of ForbMod fits to EPHIN data based on visual analysis: a)
a good fit (event1, inner structure borders); b) an asymmetric fit (event1,
extended ICME borders); c) a few-data-points fit (event2, inner structure
borders); d) a large-scatter fit (event2, extended ICME borders); e) an under-
recovery fit (event4, extended ICME borders); and f) a substructuring fit
(event17, extended ICME borders). FD_obs was FD amplitude-determined using a
traditional algorithm-based observational method, whereby FD was measured from
the onset to the minimum. FD_bf is the FD amplitude obtained from the ForbMod
best-fit procedure and MSE is the mean standard error.

We next performed the ForbMod best fit on the sample of the events, as is
described in Section 2.3. We first analysed the fits visually and found we
could divide the fits into six types:

  * •

a good fit \- where there was visually a good match between data points and
the best-fit function

  * •

an asymmetric fit \- where the fitted function was shifted along the x axis
compared to the data points

  * •

a few-data-points fit \- where the fit was performed on only a small number of
data points

  * •

a large-scatter fit \- where the fit was performed on data points with a large
scatter

  * •

an under-recovery fit \- where the data points did not return to the pre-
decrease level

  * •

a substructuring fit \- where a single fit was performed while the data points
indicated possible substructuring

Examples of fit types are given in Figures 9a-f. However, we note that some
events may fall into more than one category; for example, an event might show
under-recovery as well as substructuring. As can be seen in Table 2 all but
one good fit events have inner structure borders. They also on average have
the smallest MSE (3.9e-5). A few-data-points fit is another type present only
for inner structure borders events. On the other hand, the large-scatter fit
and the substructuring fit seem to be categories present predominantly for
extended ICME borders and on average have the largest MSEs (1.8e-4 and 1.5e-4,
respectively). We note that for all of the events $MSE<0.01$. Therefore,
although some profiles do not show a visually pleasing FD, the ForbMod best-
fit function still manages to find a solution with a relatively small mean
standard error (MSE).

We next analysed best-fit results for four different QIs, as is described in
Section 2.3 (see Table 1). We do not find notable differences between events
marked by different QIs, in either the difference between FD_bf and FD_obs,
the values of MSE, or the spread of MSE. Surprisingly, we find that MSE on
average takes the smallest values for QI 1, i.e. the least clear events. This
is likely related to the fact that these events also on average have the
smallest FD amplitudes, both FD_bf and FD_obs (1.1% and 1.2%, respectively).

Table 1: Statistics for the best-fit amplitude, FD_bf, the amplitude obtained using a traditional algorithm-based observational method, FD_obs, and the MSE for four different QIs (as is described in Section 2.3). | QI 1 | QI 2 | QI 3 | QI 4
---|---|---|---|---
No. of events | 6 | 12 | 8 | 4
inner structure
No. of events | 4 | 10 | 7 | 3
FD_obs mean | 1,2 | 3,0 | 4,0 | 3,5
FD_bf mean | 1,1 | 2,8 | 3,6 | 3,1
MSE mean | 7,8E-06 | 5,0E-05 | 5,7E-05 | 2,7E-05
MSE stdev | 4,7E-06 | 6,0E-05 | 4,2E-05 | 3,3E-05
extended ICME structure
No. of events | 6 | 10 | 6 | 3
FD_obs mean | 1,8 | 3,6 | 4,0 | 5,3
FD_bf mean | 1,7 | 2,2 | 2,6 | 5,2
MSE mean | 4,8E-05 | 1,1E-04 | 1,3E-04 | 2,0E-04
MSE stdev | 4,7E-05 | 1,1E-04 | 9,8E-05 | 2,7E-04

Figure 10: Linear regression between the best-fitted FD amplitude, FD_bf, and
the FD amplitude measured using the traditional algorithm-based observational
method (from the onset to the minimum), FD_obs, for two different border
selections (inner structure borders and extended ICME borders)

We next compared the calculated best-fit amplitude, FD_bf, with the amplitude
obtained using a traditional algorithm-based observational method, FD_obs,
whereby FD_obs was measured from the onset to the minimum. For the inner
structure borders, we find the average FD_bf of 2.8%, which is somewhat
smaller than the average FD_obs (3.1%). We note that this difference is
expected, because the traditional algorithm-based observational method
measures maximum amplitude, without taking into account possible scatter of
the data. For extended ICME borders we find an even larger difference between
the two, with an average FD_bf of 2.5% and an average FD_obs of 3.5%. This
reflects the fact that the scatter of the data that was fitted is on average
larger for extended ICME borders compared to inner structure borders, as can
be seen from the average values of MSE being 1.1e-4 and 4.2e-5, respectively.

We also performed a linear regression and calculated the Pearsons correlation
coefficient for FD_bf versus FD_obs separately for inner structure borders and
extended ICME borders. The results are presented in Figure 10 and confirm the
findings above. We observe a negative intercept in both linear regressions,
confirming a systematically smaller FD_bf compared to FD_obs. The larger
negative intercept for extended ICME borders indicates a larger difference
between FD_bf and FD_obs. In addition, we can see that the Pearsons
correlation coefficient (i.e. its squared value given by $R^{2}$) is smaller
for extended ICME borders related to the larger scatter. However, we do note
that the correlation is high for both border selection cases ($R>0.8$). Thus,
we can conclude that the ForbMod best-fit procedure performs similarly to the
traditional algorithm-based observational method, but with slightly smaller
values for the FD amplitude, as it’s taking into account the noise in the
data. We note that in two events FD_obs could not be determined because of the
data gap at the onset of the FD (see Table 2). In those cases, FD_bf could
still be determined, which is an obvious advantage of the best-fit procedure.

##  4 Summary and conclusions

We developed a new method to measure the FD amplitude, which is based on the
analytical diffusion-expansion FD model for FRs (ForbMod). The ForbMod best-
fit procedure we developed performs a non-linear fitting of the Bessel
function of the zero order (i.e. the radial part of the ForbMod solution) to
the designated time series. The non-linear fitting was performed by
calculating all possible ForbMod curves constrained within FR borders to the
designated dataset and minimising the MSE. The input for the procedure is the
designated dataset with defined FR borders (i.e. the dataset translated into
the FR coordinate system radial direction, with the centre of the FR at 0 and
FR borders at +/-1, respectively). The output of the procedure is the best-
fitted amplitude of the measured FD and MSE.

In order to evaluate the performance of the ForbMod best-fit procedure, we
used synthetic measurements. Synthetic measurements were produced by
calculating the theoretical ForbMod curve for a specific example CME and then
applying various effects to the data to mimic the imperfection of the real
measurements. Using synthetic measurements we could test how far away the
best-fitted function is from the real one for different cases. We found that
if the FD is clearly visible in the data the reliable ForbMod best-fit
function should result in the MSE being ¡0.1 with an expected ¡10% relative
error in the fitted amplitude, roughly along the identity line.

We next tested the ForbMod best-fit function on the real data, measured by
detector F of the SOHO/EPHIN instrument. We used a list of SOHO/EPHIN detected
FDs from the time period 1995–2015 and a catalogue of near-Earth ICMEs to
compile a list of events suitable for analysis. The sample contains 30 events,
all of which have a distinct FD corresponding to the magnetic obstacle (MO).
We categorised the events according to different QIs, which mark different
levels of clarity of the events. We extracted FD profiles for each event in
the sample from the onset to the end of FD and prepared them for the
application of the ForbMod best-fit function. The extraction of FD profiles
(from the onset to the end) was performed manually by an observer, whereby we
applied two different versions of border selection: 1) the inner substructure
of the CME, which shows clear MO signatures, and 2) the extended ICME
structure, which may also encompass outer, disturbed layers of the ICME.

We do not find notable differences between events marked by different QIs. We
compared the calculated best-fit amplitude with the amplitude obtained using a
traditional algorithm-based observational method (measured from the onset to
the minimum) and found a good match between the two. The ForbMod best-fit
procedure performs similarly to the traditional algorithm-based observational
method, but with slightly smaller values for the FD amplitude, as it takes
into account the noise in the data. For events with two different versions of
border selection we find that the best-fit applied on extended ICME structure
borders results in slightly larger MSE and differences compared to the
traditional method due to the larger scatter of the data points. We find that
the best-fit results can visually be categorised into six different FD profile
types. Although some profiles do not show a visually pleasing FD, the ForbMod
best-fit function still manages to find a solution with a relatively small
MSE.

Finally, we find that the best-fit procedure has an advantage compared to the
traditional method as it can estimate the FD amplitude even when there is a
data gap at the onset of the FD. The method could be further advanced by
taking into account best-fitting for multiple-border selection, and thus
reducing the subjective impact of the observer. With further improvements in
the direction of border selection, the method would not only advance the
estimation of the magnitude of the FD, but also the determination of its
onset, end, and consequently duration. Such an application would provide the
means for a quick, objective, reliable, and robust way of measuring FDs and
may eventually facilitate the automation of FD measurements.

###### Acknowledgements.

This research was performed in the scope of the ESA funded project “Forbush
decrease analysis using model fitting to SOHO/EPHIN data (ForbMod)” (ESA
contract No. 4000135501 lz1lNL/SC/hm). M.D. acknowledges support by the
Croatian Science Foundation under the project IP-2020-02-9893 (ICOHOSS). We
acknowledge the OMNIWeb database as the source of data used, as well as the
SOHO/EPHIN. SOHO/EPHIN is supported by the Ministry of Economics via DLR grant
50OG1702. The OMNI data were obtained from the GSFC/SPDF OMNIWeb interface at
https://omniweb.gsfc.nasa.gov.

## References

  * Arunbabu et al. (2013) Arunbabu, K. P., Antia, H. M., Dugad, S. R., et al. 2013, Astron. Astrophys., 555, A139
  * Belov et al. (2021) Belov, A., Papaioannou, A., Abunina, M., et al. 2021, Astrophys. J., 908, 5
  * Belov (2009) Belov, A. V. 2009, in IAU Symp., Vol. 257, IAU Symp., ed. N. Gopalswamy & D. F. Webb, 439–450
  * Benella et al. (2020) Benella, S., Laurenza, M., Vainio, R., et al. 2020, Astrophys. J., 901, 21
  * Bothmer & Schwenn (1998) Bothmer, V. & Schwenn, R. 1998, Ann. Geophys., 16, 1
  * Cane (2000) Cane, H. V. 2000, Space Sci. Rev., 93, 55
  * Cane et al. (1996) Cane, H. V., Richardson, I. G., & von Rosenvinge, T. T. 1996, J. Geophys. Res., 101, 21561
  * Cane et al. (1995) Cane, H. V., Richardson, I. G., & Wibberenz, G. 1995, Int. Cosmic Ray Conf., 4, 377
  * Davies et al. (2023) Davies, E. E., Winslow, R. M., & Lawrence, D. J. 2023, Astrophys. J., 943, 83
  * Démoulin & Dasso (2009) Démoulin, P. & Dasso, S. 2009, Astron. Astrophys., 498, 551
  * Démoulin et al. (2008) Démoulin, P., Nakwacki, M. S., Dasso, S., & Mandrini, C. H. 2008, Solar Phys., 250, 347
  * Dumbović et al. (2018) Dumbović, M., Heber, B., Vršnak, B., Temmer, M., & Kirin, A. 2018, Astrophys. J., 860, 71
  * Dumbović et al. (2011) Dumbović, M., Vršnak, B., Čalogović, J., & Karlica, M. 2011, Astron. Astrophys., 531, A91+
  * Dumbović et al. (2020) Dumbović, M., Vršnak, B., Guo, J., et al. 2020, Solar Phys., 295, 104
  * Freiherr von Forstner et al. (2021) Freiherr von Forstner, J. L., Dumbović, M., Möstl, C., et al. 2021, Astron. Astrophys., 656, A1
  * Gulisano et al. (2012) Gulisano, A. M., Démoulin, P., Dasso, S., & Rodriguez, L. 2012, Astron. Astrophys., 543, A107
  * Guo et al. (2020) Guo, J., Wimmer-Schweingruber, R. F., Dumbović, M., Heber, B., & Wang, Y. 2020, Earth and Planetary Physics, 4, 62
  * Heber et al. (2015) Heber, B., Wallmann, C., Galsdorf, D., et al. 2015, Central European Astrophysical Bulletin, 39, 75
  * Janvier et al. (2021) Janvier, M., Démoulin, P., Guo, J., et al. 2021, Astrophys. J., 922, 216
  * Jian et al. (2006) Jian, L., Russell, C. T., Luhmann, J. G., & Skoug, R. M. 2006, Solar Phys., 239, 393
  * Kilpua et al. (2017) Kilpua, E., Koskinen, H. E. J., & Pulkkinen, T. I. 2017, Living Reviews in Solar Physics, 14, 5
  * King & Papitashvili (2005) King, J. H. & Papitashvili, N. E. 2005, J. Geophys. Res., 110, A02104
  * Klein & Burlaga (1982) Klein, L. W. & Burlaga, L. F. 1982, J. Geophys. Res., 87, 613
  * Krittinatham & Ruffolo (2009) Krittinatham, W. & Ruffolo, D. 2009, Astrophys. J., 704, 831
  * Kuwabara et al. (2009) Kuwabara, T., Bieber, J. W., Evenson, P., et al. 2009, J. Geophys. Res., 114, A05109
  * Laitinen & Dalla (2021) Laitinen, T. & Dalla, S. 2021, Astrophys. J., 906, 9
  * Laster et al. (1962) Laster, H., Lenchek, A. M., & Singer, S. F. 1962, J. Geophys. Res., 67, 2639
  * Lockwood (1971) Lockwood, J. A. 1971, Space Sci. Rev., 12, 658, 10.1007/BF00173346
  * Lopez (1987) Lopez, R. E. 1987, J. Geophys. Res., 92, 11189
  * Morrison (1956) Morrison, P. 1956, Phys. Rev., 101, 1397
  * Müller-Mellin et al. (1995) Müller-Mellin, R., Kunow, H., Fleißner, V., et al. 1995, Solar Phys., 162, 483
  * Munakata et al. (2006) Munakata, K., Yasue, S., Kato, C., et al. 2006, Advances in Geosciences, Volume 2: Solar Terrestrial (ST), 2, 115
  * Nieves-Chinchilla et al. (2018) Nieves-Chinchilla, T., Vourlidas, A., Raymond, J. C., et al. 2018, Solar Phys., 293, 25
  * Petukhova et al. (2019) Petukhova, A. S., Petukhov, I. S., & Petukhov, S. I. 2019, Astrophys. J., 880, 17
  * Quenby et al. (2008) Quenby, J. J., Mulligan, T., Blake, J. B., Mazur, J. E., & Shaul, D. 2008, J. Geophys. Res., 113, A10102
  * Richardson & Cane (1995) Richardson, I. G. & Cane, H. V. 1995, J. Geophys. Res., 100, 23397
  * Richardson & Cane (2010) Richardson, I. G. & Cane, H. V. 2010, Solar Phys., 264, 189
  * Rouillard (2011) Rouillard, A. P. 2011, J. Atmos. Sol. Terr. Phys., 73, 1201
  * Singer et al. (1962) Singer, S. F., Lenchek, A. M., & Laster, H. 1962, Journal of the Physical Society of Japan Supplement, 17, 583
  * Subramanian et al. (2009) Subramanian, P., Antia, H. M., Dugad, S. R., et al. 2009, Astron. Astrophys., 494, 1107
  * Tortermpun et al. (2018) Tortermpun, U., Ruffolo, D., & Bieber, J. W. 2018, Astrophys. J. Lett., 852, L26
  * Vršnak et al. (2019) Vršnak, B., Amerstorfer, T., Dumbović, M., et al. 2019, Astrophys. J., 877, 77
  * Zhang et al. (2021) Zhang, J., Temmer, M., Gopalswamy, N., et al. 2021, Progress in Earth and Planetary Science, 8, 56
  * Zurbuchen & Richardson (2006) Zurbuchen, T. H. & Richardson, I. G. 2006, Space Sci. Rev., 123, 31

##  Appendix A The sample

Table 2: Full sample of analysed ICMEs and FDs event | ICME start | QI | inner structure | extended ICME
---|---|---|---|---
No. | date |  | DOY start | DOY end | FD_obs | FD_bf | MSE | fit type | DOY start | DOY end | FD_obs | FD_bf | MSE | fit type
1 | 2008-09-17 | 4 | 261,2 | 262,1 | 2,7 | 2,3 | 2,9E-06 | 1 | 261,1 | 262,6 | 2,7 | 2,0 | 2,4E-05 | 2
2 | 2009-09-30 | 1 | 273,3 | 273,6 | 0,8 | 0,3 | 2,5E-06 | 3 | 273,1 | 274,4 | 1,0 | 0,8 | 6,3E-06 | 4
3 | 2009-11-14 | 2 | 318,5 | 318,9 | 0,5 | 0,6 | 1,9E-06 | 3 | 317,9 | 319,2 | 2,6 | 2,5 | 2,1E-05 | 5, 6
4 | 2010-02-07 | 4 | 38,8 | 40,0 | 5,1 | 4,6 | 6,4E-05 | 1 | 38,7 | 40,4 | 4,7 | 5,1 | 5,4E-05 | 5
5 | 2010-04-05 | 4 | N/A | N/A | N/A | N/A | N/A | N/A | 95,4 | 96,8 | 8,5 | 8,4 | 5,1E-04 | 2, 5, 6
6 | 2010-05-28 | 2 | 148,8 | 149,1 | 5,0 | 4,7 | 1,6E-04 | 3 | 148,6 | 149,9 | 7,5 | 5,5 | 3,0E-04 | 5, 6
7 | 2010-06-21 | 2 | 172,3 | 172,6 | 2,2 | 1,9 | 1,1E-05 | 3, 5 | 172,1 | 173,8 | 3,0 | 1,3 | 3,3E-05 | 6
8 | 2010-12-20 | 1 | 354,4 | 355,0 | 0,9 | 0,9 | 5,5E-06 | 3 | 353,7 | 355 | 1,6 | 1,5 | 4,1E-05 | 5, 6
9 | 2011-03-30 | 2 | 88,9 | 90,0 | 5,9 | 5,4 | 1,4E-04 | 5 | 88,9 | 91,7 | 5,7 | 3,1 | 2,8E-04 | 2
10 | 2011-09-17 | 3 | 260,6 | 262,0 | 6,1 | 6,3 | 8,0E-05 | 1 | N/A | N/A | N/A | N/A | N/A | N/A
11 | 2012-10-31 | 1 | N/A | N/A | N/A | N/A | N/A | N/A | 305,9 | 307,3 | N/A | 3,0 | 1,4E-04 | 2
12 | 2012-11-26 | 2 | 331,6 | 332,2 | 1,7 | 1,4 | 3,0E-06 | 1 | 331,5 | 332,5 | 1,6 | 0,8 | 4,1E-05 | 4
13 | 2013-02-16 | 1 | 48,3 | 48,9 | 1,3 | 1,2 | 1,3E-05 | 2 | 47,8 | 49,8 | 2,1 | 1,4 | 1,8E-05 | 4, 6
14 | 2013-04-30 | 1 | N/A | N/A | N/A | N/A | N/A | N/A | 120,5 | 121,3 | 2,0 | 2,0 | 4,3E-05 | 4, 5
15 | 2013-06-06 | 3 | 157,6 | 158,4 | 4,6 | 4,2 | 4,7E-05 | 5 | 157,6 | 159 | 5,3 | 5,0 | 2,7E-04 | 5, 6
16 | 2013-06-27 | 3 | 179,1 | 179,6 | 5,3 | 4,2 | 1,2E-04 | 3 | 179,1 | 180,3 | 5,3 | 2,4 | 1,9E-04 | 6
17 | 2013-07-05 | 1 | 186,7 | 187,3 | 1,8 | 2,0 | 1,0E-05 | 1 | 186,7 | 188,7 | 2,1 | 1,6 | 3,8E-05 | 6
18 | 2013-07-12 | 3 | 194,1 | 195,2 | 4,7 | 4,3 | 9,5E-05 | 2, 5 | 194,1 | 196,0 | 4,7 | 3,1 | 2,0E-04 | 2, 5, 6
19 | 2013-11-11 | 2 | 315,7 | 316,2 | 1,5 | 1,3 | 1,5E-05 | 3 | 315,5 | 316,5 | 3,9 | 1,3 | 1,2E-04 | 4
20 | 2013-12-15 | 2 | 349,7 | 350,2 | 6,0 | 6,0 | 1,0E-04 | 1 | N/A | N/A | N/A | N/A | N/A | N/A
21 | 2014-04-29 | 3 | 119,9 | 120,7 | 2,2 | 2,4 | 1,6E-05 | 1 | 119,8 | 120,9 | 2,6 | 2,1 | 3,7E-05 | 2
22 | 2014-05-30 | 4 | 150,6 | 152,3 | 2,8 | 2,3 | 1,3E-05 | 5 | N/A | N/A | N/A | N/A | N/A | N/A
23 | 2014-08-19 | 3 | 231,7 | 232,2 | 1,8 | 1,5 | 2,6E-05 | 3 | 231,7 | 233,8 | 1,8 | 0,4 | 3,9E-05 | 4
24 | 2016-04-14 | 3 | 105,3 | 106,0 | 3,4 | 2,5 | 1,4E-05 | 5 | N/A | N/A | N/A | N/A | N/A | 6
25 | 2016-07-24 | 2 | N/A | N/A | N/A | N/A | N/A | N/A | 206,7 | 207,7 | N/A | 1,6 | 6,8E-06 | 1
26 | 2016-10-12 | 3 | N/A | N/A | N/A | N/A | N/A | N/A | 287 | 288,7 | 4,3 | 2,7 | 7,2E-05 | 4
27 | 2016-11-09 | 2 | N/A | N/A | N/A | N/A | N/A | N/A | 314,7 | 315,7 | 2,2 | 2,0 | 4,2E-05 | 6
28 | 2017-04-13 | 2 | 104 | 104,6 | 2,3 | 1,7 | 2,1E-05 | 2 | N/A | N/A | N/A | N/A | N/A | N/A
29 | 2018-08-25 | 2 | 237,5 | 238,2 | 2,4 | 2,3 | 3,0E-05 | 6 | 237,5 | 238,4 | 2,4 | 1,6 | 1,0E-04 | 2
30 | 2019-05-13 | 2 | 134,2 | 134,7 | 2,5 | 2,4 | 1,9E-05 | 5 | 134 | 136,0 | 3,7 | 2,0 | 1,5E-04 | 6
The ICME start date corresponds to the first detected ICME signature, which
may be shock, sheath, pile-up, or MO (depending on the event).
The QI marks the clarity of MO and FD signatures and borders and is
subjectively determined by the observer, as is explained in Section 2.3.
Measurements are separated for the inner structure and extended ICME, as is
explained in Section 2.4.
For each, we determined the onset and end of the FD using the day of the year
(DOY), DOY start, and DOY end, respectively;
FD amplitude measured using the traditional method, FD_obs; best-fitted FD
amplitude, FD_bf; MSE;
and fit type as visually determined by the observer (for details see Section
3.2).