Analysis of Coronal Mass Ejection Flux Rope Signatures Using 3DCORE and Approximate Bayesian Computation

Coronal mass ejections (CMEs) are the most violent and energetic events that occur within our solar system and have a significant impact on the interplanetary magnetic field and planetary magnetospheres (Schwenn et al. 2005; Chen 2011; Webb & Howard 2012). An enormous amount of magnetized plasma is ejected into the interplanetary medium and propagates as an extremely large, and continuously expanding, structure (e.g., Burlaga et al. 1981; Farrugia et al. 1993; Gopalswamy et al. 2000) that can reach the outer planets of our solar system (e.g., Witasse et al. 2017). CMEs also carry along a strong internal magnetic field, believed to be in the form of a magnetic flux rope (MFR, e.g., Marubashi 1986; Burlaga 1988; Lepping et al. 1990), that can manifest itself as a magnetic cloud when measured in situ (Burlaga et al. 1981; Klein & Burlaga 1982; Bothmer & Schwenn 1998) by a spacecraft. This strong magnetic field, given a certain configuration of the CME, can also induce what is known as a geomagnetic storm (e.g., Farrugia et al. 1993; Gonzalez et al. 1994). Geomagnetic storms are associated with a variety of phenomena such as aurorae, geomagnetically induced currents (e.g., Pirjola 1983; Boteler et al. 1998; Bolduc 2002), and disturbances within the ionosphere (e.g., Proelss 1980), which adversely affect high-frequency ground-to-ground radio or spacecraft communication. CMEs, alongside other solar events, can furthermore pose a significant radiation hazard for human space travel (e.g., Zeitlin et al. 2013).

As such, the study of CMEs has been of high interest to the solar physics and space weather community ever since their discovery during the Skylab era (Tousey 1973) and connection with terrestrial phenomena (Gosling et al. 1991). Nonetheless, there remain a number of unresolved issues regarding their generation, interplanetary evolution, and internal structure. Some of these issues are exacerbated due to the fact that most CMEs are observed only by individual satellites, and their global structure therefore remains largely hidden. Solar imagers and coronagraphs are capable of showing the formation and eruption phases of CMEs, but these observations are susceptible to projection effects. They are also not fully representative of the global structure within the interplanetary medium, because CMEs can undergo drastic change due to rotation or deflection (e.g., Vourlidas et al. 2011; Kay et al. 2015; Möstl et al. 2015), or due to interaction with the solar wind (e.g., Riley & Crooker 2004; Manchester et al. 2017; Luhmann et al. 2020).

Further uncertainty exists on the structure of the internal magnetic field, which determines the geoeffectivity. The magnetic field structures have been approximately described using several different flux rope models, such as those based on the analytical cylindrical Lunquist solution (Lundquist 1950; Burlaga 1988; Zhang & Burlaga 1988; Lepping et al. 1990; Owens 2006; Owens et al. 2006; Kay et al. 2017) or the analytical cylindrical Gold–Hoyle solutions (e.g., Gold & Hoyle 1960; Farrugia et al. 1999). The most recent studies have also introduced models for slightly altered geometries, such as torii (e.g., Hidalgo & Nieves-Chinchilla 2012; Vandas & Romashets 2017) or elliptical cylinders (e.g., Hidalgo et al. 2002; Nieves-Chinchilla et al. 2018), in order to better accommodate for distortions further away from the idealized cylindrical flux rope picture. Other models that do not use flux ropes include those that are based on spheromaks (Farrugia et al. 1995; Vandas et al. 1997). Additionally, nonspecified field structures were calculated with the Grad–Shafranov equation (Hu & Sonnerup 2002; Möstl et al. 2009), and they result in a field closely resembling a Gold–Hoyle flux rope (Hu et al. 2015). In many cases multiple different models are able to reproduce the same observed measurements with similar accuracy and are therefore indistinguishable in the absence of independent auxiliary measurements.

A better understanding of the global structure of CMEs and their embedded magnetic field can be gained by observing multipoint events, i.e., CMEs that were observed by multiple spacecraft at distinct positions within the heliosphere. These events are inherently rare due to the low number of operating spacecraft, and only a few dozen have been cataloged so far (e.g., Burlaga et al. 1981; Leitner et al. 2007; Good et al. 2019; Vršnak et al. 2019; Salman et al. 2020). With the advent of next-generation spacecraft dedicated to solar physics, such as the Parker Solar Probe (PSP; Fox et al. 2016) or Solar Orbiter (Müller et al. 2013), along with other missions, such as BepiColombo, that have onboard magnetometers; longer-operating spacecraft, such as STEREO-Ahead (Kaiser et al. 2008); as well as possible interplanetary CubeSats, there will be up to half a dozen spacecraft within 1 au sampling the heliosphere in the upcoming years. Combined with the rise of solar cycle 25, there will likely be a considerable number of additional multipoint CME observations available.

The subsequent study of such multipoint events will require global CME models that can simulate measurements at several distinct positions simultaneously. The general approach would require expensive MHD simulations of the inner heliosphere for the solar wind environment and the evolving CME. Examples of such simulation codes include Enlil (Odstrcil et al. 2004), EUHFORIA (Pomoell & Poedts 2018; Verbeke et al. 2019), and MAS (Török et al. 2018). Unfortunately, the use of these models is limited due to their computational complexity and inherent issues in the generation of the boundary conditions on and near the Sun. On the other hand, the most simple analytical flux rope models that are used for describing the magnetic field locally (Lepping et al. 1990; Leitner et al. 2007) make use of rigid geometries that are incapable of accounting for any of the interplanetary evolution for CME structure that is expected to occur. Due to these obstacles, there has been a recent focus on using semi-empirical models as an alternative (Isavnin 2016; Möstl et al. 2018; Rouillard et al. 2020). Semi-empirical models have the major advantage that they are computationally inexpensive and conceptually easy to implement. It is furthermore possible, with varying degrees of complexity, to include simplified interactions with the solar background wind. These models should be seen as an attempt to bridge the gap between the "expensive" MHD simulations and the geometrically simple but physically accurate analytical models.

In this paper we will showcase a significant update to the 3D coronal rope ejection (3DCORE) model (Möstl et al. 2018), which is such a semi-empirical model. 3DCORE is a forward simulation model that describes a CME as a propagating and self-similarly expanding torus-like structure, influenced by a simple drag model that stays attached to the Sun. The cross section of this torus is extended to also incorporate elliptical shapes (i.e., different aspect ratios). This torus-like structure contains an embedded Gold–Hoyle-like field with a time-invariant twist number.

Our assumption of self-similar expansion for a torus-like structure is not able to describe deformations of the flux rope structure that can occur due to drastic longitudinal solar wind speed gradients. Nonetheless, it can be expected that the approximation is generally adequate for modeling measurements on a small scale. The model can therefore also be used to infer the distance scales at which the approximation breaks down. This can be of use for future studies because it gives an estimate for the required spatial resolution for resolving deformations.

We improve the 3DCORE model implementation so that it allows us to generate extremely large mega-ensembles on the order of 10⁶ runs per second. We use this newly gained efficiency to deploy an approximate Bayesian computation sequential Monte Carlo (ABC-SMC) algorithm for fitting the in situ magnetic flux rope measurements. This ABC-SMC algorithm is a Bayesian inference algorithm that is capable of not only generating global best-fit solutions, but also estimating constraints on the model parameters. The algorithm can be additionally fine-tuned by incorporating priors from either auxiliary measurements or physical considerations. The specific implementation of the fitting algorithm is easily extendable to multipoint constellations, which will facilitate the study of multipoint events.

Because of the novelty of this algorithm in the context of studying CMEs, we will briefly show a numerical test, using synthetically generated measurements, to show that our method delivers the correct results in the best-case scenario when the ground truth is known. Additionally, we apply our model and algorithm to a flux rope observed in situ by PSP on 2018 November 12 at 0.25 au, shortly after its first close flyby of the Sun, to show that our method is also applicable to real-world scenarios.

This study is structured as follows. In Section 2 we introduce our improved 3DCORE model and provide an overview of the basic assumptions that we use. Section 3 describes the motivation for the usage of the ABC-SMC algorithm, gives a very brief introduction to the topic, and goes into detail on our specific implementation. Section 3.4 shows the results of our ABC-SMC algorithm with respect to synthetic measurements. In Section 4 we apply the 3DCORE model and the ABC-SMC algorithm to the flux rope observed by PSP and interpret the results. We discuss the applications of our novel approach in Section 5.

In this section we give an in-depth description of the 3DCORE model and the accompanying data-generating process. We repeat most of the basic concepts of the model, which can also be found in the original paper by Möstl et al. (2018), for the purpose of clarity. The model serves to empirically describe and reproduce the general properties of CME flux rope measurements as we believe them to be. We attach ourselves to the concept of a bent flux tube that is connected on both sides to the Sun, describing the flux rope on a global scale. The self-consistency of this global picture can be tested by applying it to multipoint events, i.e., flux ropes that were measured by multiple spacecraft at significantly different positions. Such a study is well beyond the scope of the current paper, but the methods that we introduce can very easily be extended to multipoint analysis, and we therefore lay the necessary groundwork for the future.

The solar wind environment heavily influences the evolution of any CME during its propagation within the interplanetary medium. As the solar wind is far from uniform in density or speed, specifically in the angular component, it is expected that the global shape of the CME is continuously deformed. Because of kinematic interactions, the flux rope is also expected to be flattened in the direction of propagation (Riley & Crooker 2004). For our purposes we limit ourselves to describing the solar wind background using a single global speed value and associated global density coefficient. We assume that the geometry of our coronal mass ejection expands self-similarly during its slowed-down or accelerated propagation without any deformation effects. The kinematical flattening effect is approximated using a constant elliptical cross section.

The magnetic field is inserted into our flux rope shape using the same procedure as for local analytical magnetic field models. The specific implementation allows us to insert any arbitrary analytical magnetic field model into our chosen shape with only some minor modifications. This allows us greater freedom in choosing the magnetic field models for analysis, and future work will include extensive comparisons between different models.

In summary, we split the 3DCORE model into three different components. These are namely the shape model, which describes the global geometry; the propagation model, which describes the interaction with the background solar wind and the self-similar expansion; and the magnetic field model, which inserts an analytical magnetic field into our chosen shape.

2.1. Shape Model

The global shape of our flux rope is described using a custom curvilinear coordinate system, further denoted as Q, that is defined using the parameterization shown in Equation (1).

$$\begin{eqnarray}&&x^{\prime} =-\left[{\rho }_{0}+r{\rho }_{1}{\left(\sin \displaystyle \frac{\psi }{2}\right)}^{2}\sin \phi \right]\cos \psi +{\rho }_{0}\end{eqnarray} \tag{ 1a }$$

$$\begin{eqnarray}&&y^{\prime} =\left[{\rho }_{0}+r{\rho }_{1}{\left(\sin \displaystyle \frac{\psi }{2}\right)}^{2}\sin \phi \right]\sin \psi \end{eqnarray} \tag{ 1b }$$

$$\begin{eqnarray}&&z^{\prime} =\delta r{\rho }_{1}{\left(\sin \displaystyle \frac{\psi }{2}\right)}^{2}\cos \phi .\end{eqnarray} \tag{ 1c }$$

The mappings $f^{\prime} :Q\mapsto {{\mathbb{R}}}^{3}$ and $g^{\prime} :{{\mathbb{R}}}^{3}\mapsto Q$ define the forward and backward coordinate transformations that are required to transform Q-coordinates into Cartesian coordinates and vice versa. The coordinates in Q are denoted as $r\in [0,\infty ),\psi \in [0,2\pi ),\phi \in [0,2\pi )$. The inverse mapping g' is shown in Equation (2)

$$\begin{eqnarray}&&\tan \psi =-\displaystyle \frac{y^{\prime} }{{\rho }_{0}-x^{\prime} },\end{eqnarray} \tag{ 2a }$$

$$\begin{eqnarray}&&\tan \phi =\displaystyle \frac{z^{\prime} }{\delta \left(\sqrt{{\left({\rho }_{0}-x^{\prime} \right)}^{2}+y{{\prime} }^{2}}-{\rho }_{0}\right)},\end{eqnarray} \tag{ 2b }$$

$$\begin{eqnarray}&&r=\displaystyle \frac{\sqrt{{\left({\rho }_{0}-x^{\prime} \right)}^{2}+y{{\prime} }^{2}}-{\rho }_{0}}{{\left(\sin \tfrac{\psi }{2}\right)}^{2}{\rho }_{0}\cos \phi }.\end{eqnarray} \tag{ 2c }$$

This curvilinear coordinate system is based on the toroidal coordinate system. The parameters ${\rho }_{0}$ and ${\rho }_{1}$ define the major and minor radius of the base torus. The δ parameter controls the ellipticity of the cross section. A value of δ > 1 corresponds to a "flattened" cross section, as one would expect from flux ropes, due to kinematic effects (some other papers use an inverse definition). The flux rope structure itself is defined by the implicit volume r ≤ 1. This leads to a torus-like shape, with an elliptical cross section of variable size that is permanently attached to the point of origin, which acts as the Sun on the thinnest side. The resulting geometry is illustrated in Figure 1.

Zoom In Zoom Out Reset image size

The parameterization is chosen so that the frontal part of the flux rope geometry, which indicates the direction of propagation, always points toward the positive x'-axis. To allow for general propagation in any direction, we introduce a rotation R that is defined by three further parameters. These are the latitude L_1; the longitude L_2; and the inclination, tilt or orientation O. They represent the direction of the propagation of the flux rope structure and its respective orientation. By our definition, an inclination value of 0° (and 180°) lies within the XY plane of the rotated Cartesian coordinate system. The final mappings that we use in order to transform from general Cartesian coordinates, which are ideally coordinates in an inertial reference frame, are defined by $f=R\,\circ \,f^{\prime} ,g=g^{\prime} \,\circ \,{R}^{-1}$. Unless otherwise stated we will be using the Heliocentric Inertial (HCI) coordinate system as default.

The fact that g' can still be given in an analytical form for this particular shape is of great advantage. It allows extremely simple collision detection and evaluation of the magnetic field when given a magnetic field function $B(r,\psi ,\phi )$ in Q-coordinates. This approach leads to a significant computational speed-up and is the most important improvement when compared to the implementation used in the original 3DCORE paper (Möstl et al. 2018), because it opens up the possibility of using more computationally expensive methods for analysis. In the general case, when using an arbitrary parameterization for the shape, this is no longer necessarily true, and numerical approximation schemes must be used.

2.2. Propagation Model

Flux rope propagation is implemented by adding time dependence to the shape parameters ${\rho }_{0}$ and ${\rho }_{1}$. This time dependence takes the following form:

$$\begin{eqnarray}&&{\rho }_{1}(t)=\displaystyle \frac{{D}_{1\mathrm{au}}{R}_{\mathrm{apex}}^{1.14}(t)}{2},\end{eqnarray} \tag{ 3 }$$

$$\begin{eqnarray}&&{\rho }_{0}(t)=\displaystyle \frac{{R}_{\mathrm{apex}}(t)-{\rho }_{1}(t)}{2},\end{eqnarray} \tag{ 4 }$$

where R_apex(t) is the distance of the apex point A to the Sun and ${D}_{1\mathrm{au}}$ is the diameter of the flux rope cross section at the widest point and at 1 au. For the propagation and interaction with the solar background wind, we consider only the kinematics of the apex point. The rest of the flux rope structure expands with the apex point in order to preserve the condition of self-similarity over time. The distance of our apex point is described using a drag-based model, based on Vršnak et al. (2013), with the following analytical form:

$$\begin{eqnarray}\begin{array}{rcl}{R}_{\mathrm{apex}}(t) & = & \pm \displaystyle \frac{1}{\gamma }\mathrm{ln}\left[(1\pm \gamma ({V}_{0}-{V}_{w})t\right]\\ & & +\,{V}_{w}t+{R}_{0}\end{array}\end{eqnarray} \tag{ 5 }$$

where V₀ is the initial CME velocity, V_w is the background solar wind speed, R₀ is the initial apex distance from the Sun, and γ is the solar background wind drag-coefficient. This propagation model is unchanged from Möstl et al. (2018). An important aspect of this particular drag-based model is the fact that there is an analytical expression for R_apex(t). We can therefore directly evaluate ${\rho }_{0}(t)$ and ${\rho }_{1}(t)$ at any point in time without simulating the propagation from the start. The computational time needed for a single 3DCORE simulation therefore scales linearly with respect to the number of chosen time points at which the magnetic field is to be simulated.

2.3. Magnetic Field Model and Embedding

Extensive work has been performed finding analytical magnetic field solutions for various flux rope geometries. The most common cases are fully analytical or approximate solutions for cylindrical (Lundquist 1950; Gold & Hoyle 1960) or toroidal geometries (Hidalgo & Nieves-Chinchilla 2012; Vandas & Romashets 2017). Additional approximate solutions exist for elliptical-cylindrical geometries (Hidalgo et al. 2002; Nieves-Chinchilla et al. 2018). Because no solutions exist for our model shape, we are forced to adapt an existing solution to our needs. In this particular case, we choose to embed the toroidal Gold–Hoyle-like solution described in Vandas & Romashets (2017). In normal toroidal coordinates this magnetic field takes the following form:

$$\begin{eqnarray}&&{B}_{r}=0,\end{eqnarray} \tag{ 6a }$$

$$\begin{eqnarray}&&{B}_{\psi }=\displaystyle \frac{{B}_{0}}{1+{b}^{2}{r}^{2}},\end{eqnarray} \tag{ 6b }$$

$$\begin{eqnarray}&&{B}_{\phi }=\displaystyle \frac{{B}_{0}{br}}{\left(1+{b}^{2}{r}^{2}\right)\left(1+r\tfrac{{\rho }_{1}}{{\rho }_{0}}\cos \phi \right)}\end{eqnarray} \tag{ 6c }$$

where b is the twist parameter. We can relate the twist parameter with the total number of twists in the structure using:

$$\begin{eqnarray}&&b=\tau \displaystyle \frac{{\rho }_{1}}{2\pi {\rho }_{0}}\end{eqnarray} \tag{ 7 }$$

where τ is now the total number of twists along the entire torus. In the case of our slightly different shape, we correct the twist factor using:

$$\begin{eqnarray}&&b=\tau \displaystyle \frac{{\rho }_{1}}{2\pi {\rho }_{0}}E(\delta )\sin {\left(\displaystyle \frac{\psi }{2}\right)}^{2}\end{eqnarray} \tag{ 8 }$$

where $E(\delta )$ is the circumference of an ellipse with an aspect ratio of δ and a minor axis length of one. Because this circumference cannot be calculated analytically, we use a numerical approximation.

Figure 2 shows traced magnetic field lines using our magnetic field model for both a normal torus and our torus-like shape. For each shape we trace two magnetic field lines positioned at coordinates r = 1 and r = 0.5 across the entire structure. This can be used as a visual verification for our magnetic field model as we used integer twist values for τ. Figure 3 shows two in situ magnetic profiles generated by simulated Earth-directed CMEs. In both profiles we see the characteristic rotation of the magnetic field in the B_z component, and evolution of the total magnetic field strength, with a peak that is slightly off-center because of expansion.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

2.4. Noise

In this section we introduce and discuss the numerical approach that we use to fit the in situ magnetic field measurements using the 3DCORE model that we described in the previous section. In particular we state the motivation and then introduce the ABC-SMC algorithm that we use for our subsequent analysis. In Section 3.4 we test the algorithm on a synthetic data set, generated by the 3DCORE itself, in order to verify that our implementation can correctly extract a known ground truth. This represents an analysis in the best-case scenario, in which the model is fully capable of reproducing and describing the measurements. We could furthermore test how our results vary in the synthetic case for different levels of artificial noise and which model parameters can be inferred more accurately.

The most common approach to fitting MFRs with an analytical magnetic field model, cylindrical or toroidal, is to minimize a customly defined error metric using an iterative gradient-descent-based minimization algorithm. A very popular error metric is the rms error, which is based on either the three magnetic field components only, or the components and the total magnetic field strength. This type of approach can be seen in the studies by Lepping et al. (1990), Vandas & Romashets (2017), and Nieves-Chinchilla et al. (2019), or papers that focus on comparing the various different methods (Riley et al. 2004). These algorithms generally guarantee only convergence toward a local minimum, and, while this approach has been shown to work very well for local magnetic field models, they tend to fail for global flux rope models (e.g., Isavnin 2016) because of the increased geometrical complexity. In general, Monte Carlo–based methods are required that are capable of searching the full parameter space and finding global minima. The primary drawback is a significant increase in the computational cost of the fitting algorithm itself.

Flux rope fitting methods that rely on the minimization of an error metric also derive only a single parameter estimate. The estimation of an error on the model parameters themselves requires multiple independent measurements of the same event (i.e., multipoint events with many satellites). Because the majority of events are observed by only a single spacecraft, and in some rare cases by up to two or three, it is hard to perform any statistical analysis for these types of measurements. This error estimate is important in order to indicate the level of confidence in the derived results. Because of the complex 3D structure of the flux ropes that are measured at only one single point, i.e., strongly projected, there may be significant ambiguities, and there is a strong possibility that a large range of different model parameters can reproduce essentially the same result. This ambiguity can be further amplified when attempting to analyze very noisy or strongly distorted flux ropes.

These issues serve as a motivation to explore a different class of inference algorithms with the intent of at least partially alleviating the aforementioned issues. In this paper we will showcase a particular implementation on a Monte Carlo–based Bayesian inference algorithm that is capable of searching the full valid parameter space (i.e., finding global minima) and generating error estimates on the resulting model parameters in the form of probability distributions.

3.1. Bayesian Inference

For our purposes we reformulate the fitting problem using Bayes' Theorem:

$$\begin{eqnarray}&&p(\theta | y)=\displaystyle \frac{{ \mathcal L }(y| \theta )p(\theta )}{p(y)}\end{eqnarray} \tag{ 9 }$$

where θ are the model parameters and y is the data set that we use for the fitting procedure. The goal is to compute the posterior $p(\theta | y)$, which is a multidimensional probability distribution over the parameter space and gives the conditional probability of θ when given the data set y. It is important to keep in mind that the probability distributions encode our belief in which value a parameter takes and that the true underlying value is still fixed.

To compute the posterior, we first need to compute the likelihood ${ \mathcal L }(y| \theta )$ and define a prior p(θ). The p(y) term can be safely ignored for all application because it simply reduces to a normalization factor due to $\int p(\theta | y)=1$.

The prior p(θ) encompasses all information that we know about our model parameters before running any analysis. This can include constraints due to physical considerations or statistics from past events (generated from extensive CME catalogs). In the case where we do not want to include any additional information, we are required to use noninformative priors, which can be constructed in various different ways. Throughout the remaining work we will always use noninformative priors with uniform distributions over a certain range of interest. As we will briefly comment on later, this may introduce significant bias and error.

The second, and more important, component is the expression for the likelihood ${ \mathcal L }(y| \theta )$. The likelihood gives the probability of generating y given model parameters θ. This highlights the need for a stochastic forward simulation model, because ${ \mathcal L }$ would simply reduce to a δ function in the case of a deterministic model. Because we only inserted randomness into our model using random Gaussian noise (see Section 2.4), we are able to write down an analytical expression for ${ \mathcal L }$. For a k-dimensional magnetic field measurement y, a simulation output $x=M(\theta )$, and a noise level of σ, the likelihood takes the following form:

$$\begin{eqnarray}&&{ \mathcal L }(y| \theta )={ \mathcal N }(y-x,\sigma {{\boldsymbol{I}}}_{k})\end{eqnarray} \tag{ 10 }$$

where ${ \mathcal N }$ is the k-dimensional multivariate normal distribution. In the more general case, and for more complicated models, it is impossible or highly impractical to find an expression for the likelihood.

The posterior is then finally computed, or more correctly approximated, by sampling from the product ${ \mathcal L }(y| \theta )p(\theta )$. Various different sampling algorithms have been developed to achieve this, with the most popular class of algorithms using Markov Chain Monte Carlo (MCMC; e.g., Hastings 1970). While an algorithm like MCMC should in theory be applicable in our case, because we are given an analytical expression for the likelihood, we found it to not work well in practice. This is most likely due to our model being capable of generating nonresults, i.e., CMEs that completely miss the observer in space or time. These simulation runs do not allow the evaluation of a likelihood and therefore define invalid regions in the parameter space. These invalid regions require special handling and significantly slowed down our MCMC implementations to the point that they were not able to generate satisfactory results.

An alternative to MCMC are sequential Monte Carlo (SMC) sampling algorithms (del Moral et al. 2006), which are also sometimes referred to as particle filters. Their primary advantage is that these algorithms are generally easily parallelizable, which is not the case for MCMC, leading to a significant speed-up for the sampling process. The SMC sampler can also better handle the invalid regions of the parameter space without some of the major pitfalls that may hamper the convergence of an MCMC, albeit still at a significant computational cost. Despite these advantages, we were still not able to generate satisfactory results using an SMC sampler because, among other things, it converged too slowly. While it should be expected that a more sophisticated SMC algorithm can overcome these issues, the method would still be limited by the requirement for an analytical expression for the likelihood, which may restrict the applicability for future models.

For this reason we decided to make use of a more general class of algorithms known as approximate Bayesian computation (ABC), sometimes also referred to as likelihood-free algorithms. In ABC the likelihood is replaced by a summary statistic. They therefore do not require an expression of the likelihood, which is their primary advantage, and as such greatly expand the realm of models to which they can be applied. Their only requirement is a more or less accurate and, more important, fast-forward simulation that can reproduce the observed measurements.

3.2. Approximate Bayesian Computation

For a more extensive introduction to ABC we refer the interested reader to Sisson et al. (2019). The principal idea behind any ABC algorithm is to bypass the likelihood ${ \mathcal L }$ using a distance metric $\rho (x,y)$, which can be defined in various ways, that measures the difference between data y and simulation outputs $x=M(\theta )$. The posterior is approximated using the ensemble:

$$\begin{eqnarray}&&p(\theta | y)\approx \mathop{\mathrm{lim}}\limits_{\epsilon \to 0}\{{\theta }_{i}| \rho ({x}_{i}=M({\theta }_{i}),y)\lt \epsilon \}\end{eqnarray} \tag{ 11 }$$

where is a threshold value and $\{{x}_{i}\}=\{M({\theta }_{i})| {\theta }_{i}\sim p(\theta )\}$ (Tavaré et al. 1997; Beaumont et al. 2002). The set of ${\theta }_{i}$ used in the ensemble in Equation (11) is initially drawn from the prior p(θ), and only the parameter candidates ${\theta }_{i}$ that satisfy $\rho ({x}_{i}=M({\theta }_{i}),y)\lt \epsilon$ are used to approximate the posterior $p(\theta | y)$. This is the simplest form of an ABC algorithm and is called the rejection algorithm. While it illustrates the core principles of any ABC algorithm, it is of practically no use except for the most simplistic cases or toy models.

In practice the rejection algorithm suffers from multiple problems. For smaller values, or as tends toward zero, the number of rejected parameter candidates ${\theta }_{i}$ becomes extremely large. This generally leads to the existence of a threshold value ${\epsilon }_{\min }$, below which no accepted candidates can be reliably found. Another source of error is the distance metric $\rho (x,y)$. Using extra-large or multidimensional data sets or model outputs for the ρ statistic will additionally increase the rejection rate of the ABC algorithm due to the curse of dimensionality. This problem is commonly circumvented by using a summary statistic S so that $\rho =\rho (S(x),S(y))$, where S reduces the dimensionality and complexity of the data set and model outputs. While the usage of this summary statistic is highly beneficial for computational efficiency, it can be another additional source of bias error (e.g., Prangle 2019). Lastly, sampling from the entire prior p(θ) is highly inefficient because large regions of the prior space can be identified as being of very little interest when using wide priors.

While some of the inherent deficiencies of the basic ABC algorithm are hard to tackle, one can significantly improve the sampling process by which candidates ${\theta }_{i}$ are drawn from p(θ). Various more sophisticated ABC algorithms have been proposed, such as those based on MCMC (Marjoram et al. 2003) and particle filters (SMC, Sisson et al. 2007; Beaumont et al. 2009). For our analysis we opted to use an ABC-SMC algorithm, for which we will explain the implementation in detail in the next section.

3.3. ABC-SMC Implementation

The first choice that we make for our implementation is the definition of the distance metric ρ and the summary static S. For our study we opted to simply use the rms error between the simulation result x and reference data y at K time points t_i:

$$\begin{eqnarray}&&{\rho }_{\mathrm{rms}}(x,y)=\sqrt{\displaystyle \frac{1}{K}\displaystyle \sum _{i=1}^{K}{\left({x}_{i}-{y}_{i}\right)}^{2}}\end{eqnarray} \tag{ 12 }$$

where x_i and y_i are magnetic field vectors. The summary static S reduces the full time-series y to only a handful of measurements at significant time intervals $\{{y}_{0},\ldots ,{y}_{K-1}\}$. We found K values of around a dozen to be sufficient for our analysis, and this generally corresponds to time intervals of one to a few hours if the time points are uniformly spaced out across the entire flux rope. Because of the setup of our model, as was detailed in Section 2.2, we can directly generate a simulation result at any time point. Using smaller K values is therefore significantly faster than using larger ones with a linear relationship for the computational complexity.

This statistic has a significant issue when two signatures with a time-shift are compared. To remedy this issue, we further introduce two control points ${t}_{S},{t}_{E}$, with t_S being set just before the CME arrives, and t_E being set after it has passed. These control points are generally inserted a few hours before and after the start and end of the flux rope. Any sample is rejected if the observer is within the synthetic flux rope at times t_S or t_E. Additionally, we require that any sample generates a valid magnetic field measurement at the K reference time points. This allows only flux rope measurements with a small time-shift and similar duration to be accepted. Assuming that the model returns a null vector when the observer is outside of the flux rope, the final summary statistic can be described as follows:

$$\begin{eqnarray}\rho (x,y)=\left\{\begin{array}{ll}\infty , & {x}_{S}\ne 0\\ & {x}_{E}\ne 0\\ & \exists i:{x}_{i}=0\\ {\rho }_{\mathrm{rms}}(x,y), & \mathrm{otherwise}\end{array}\right.\end{eqnarray} \tag{ 13 }$$

which ensures that any simulations with signatures that vary greatly in duration compared with the reference measurement are rejected. The time discretization used in the summary statistics is most likely a significant source of error when performing inference, and some of the likely effects will be discussed later in the examples.

The basic idea of the ABC-SMC algorithm is to iteratively approximate the posterior using intermediate distributions with larger threshold values . Furthermore each intermediate approximation of the posterior is generated by sampling parameter candidates from the previous approximation instead of the full initial prior. In the first iteration the ABC-SMC algorithm matches the simple rejection algorithm. We sample a set of candidate parameters $\{{\theta }_{i}^{0}\}\sim p(\theta )$ and generate the first intermediary distribution using ${P}^{0}=\{{\theta }_{i}^{0}| \rho (M({\theta }_{i}^{0}),y)\lt {\epsilon }_{0}\}$. The sampling is repeated until the set P⁰ reaches a predetermined size N that depends on the required resolution. Each "particle" within P⁰ is then assigned the same weight ${\omega }_{i}^{0}=\tfrac{1}{N}$. Because of the specific construction of our distance metric ρ, the exact value that we chose for ${\epsilon }_{0}$ is not as important, because most candidates will be rejected due to a significant time-shift when comparing the flux rope signatures. For our implementation we found ${\epsilon }_{0}=\rho (0,y)$ to be an initially high enough starting value.

In each successive iteration we generate a new set of candidate parameters by sampling from the previous intermediary instead of the prior, i.e., $\{{\theta }_{i}^{j}\}\sim {P}^{j-1}$. Because the intermediary P^j−1 is given in only approximate form by $\{{\theta }_{i}^{j-1}\}$, specific methods must be used to correctly draw new candidate parameters. This can be done by randomly picking (accounting for the weights ${\omega }_{i}^{j-1}$) a particle ${\overline{\theta }}_{i}^{j}$ from the intermediary P^j−1 and perturbing the selected particle via a perturbation kernel. Different variants of this method are explained in detail and compared in Filippi et al. (2013). The most common approach is to perturb ${\overline{\theta }}_{i}^{j}$ by a vector drawn from a multivariate normal distribution, which can be constructed by computing the covariance matrix of the overall distribution P^j−1. Larger perturbation kernels are capable of probing larger areas of the parameter space and lessen the probability that the iterative algorithm finds itself stuck in a local minimum. On the other hand, smaller kernels will lead to a higher probability of drawn candidates being accepted, which leads to faster convergence (but not necessarily toward the correct result).

For our case we opted to use a transition kernel based on M-nearest neighbors (Filippi et al. 2013), in which the covariance matrix is computed from only half of all particles within P^j−1 that are closest to ${\overline{\theta }}_{i}^{j}$ (M = N/2). This choice delivered significantly better results compared with the standard method of using the full covariance matrix due to the particular shape of the intermediary distributions. Because of our model, there are various degeneracies that are expected to occur at large values that make the entire sampling process extremely inefficient when using a large transition kernel. The only downside with our choice is that special care needs to be taken that the algorithm converges correctly and does not get stuck in a local minimum due to insufficient exploration of the parameter space.

After generating a new intermediary distribution P^j, we assign weights to each particle according to Beaumont et al. (2009):

$$\begin{eqnarray}&&{\omega }_{i}^{j+1}=\displaystyle \frac{p({\theta }_{i}^{j+1})}{{\displaystyle \sum }_{i^{\prime} =0}^{N}{\omega }_{i^{\prime} }^{j}K({\theta }_{i}^{j}| {\theta }_{i^{\prime} }^{j-1})}\end{eqnarray} \tag{ 14 }$$

where $K(\cdot | \cdot )$ describes the transition probability given by the perturbation kernel.

The final component required for the algorithm is the determination of the threshold values for . For this purpose we use an adaptive scheme based on quantiles (e.g., Lenormand et al. 2013). The ${\epsilon }_{j}$ value is computed as the αth percentile from the set $\{\rho (M({\theta }_{i}^{j}),y)| {\theta }_{i}^{j}\in {P}^{j}\}$. We generally found values of $\alpha \in [0.2,0.5]$ to work well.

Any of the generated intermediary distributions can serve as an approximation for the final posterior. As such we furthermore need a criterion when to abort the iterative algorithm. In practice, as the threshold value continues to decline, the algorithm becomes successively slower as most parameter candidates (and simulation results) are rejected. In our implementation we stop the algorithm when the number of iterations exceeds a predetermined amount, or the accept-reject ratio (the number of drawn samples that are accepted compared with those that are rejected) falls below a predetermined value. Both of these two abort conditions can be modified as required. For any reasonably complex situation, ABC algorithms are not expected to fully converge or converge very slowly. As such the main criteria for choosing the stopping point are computational and time constraints.

3.4. ABC-SMC Synthetic Example

We provide a detailed illustration of the type of results our ABC-SMC algorithm can deliver by applying the algorithm on a synthetic measurement. For this purpose we generate a magnetic field measurement using the 3DCORE model itself, with which we will then perform the analysis. This type of test is important in order to verify that our algorithm delivers correct results under ideal conditions.

Figure 4 shows the synthetic flux rope that we used for this test and the resulting fit from our ABC-SMC algorithm. The smooth solid line shows the underlying signal from the 3DCORE model. The dotted markers show the actual noisy measurement used for the analysis. The resulting fit is the dashed line with a shaded area representing the 1σ and 2σ spread of the ensemble solution.

Zoom In Zoom Out Reset image size

The synthetic measurement used was generated by the 3DCORE model using the model parameters shown in Table 1. The single parameter that is not listed in the table is the initialization time t₀, which is set to 2018 January 1 00:00. To simplify the analysis for the algorithm we furthermore fixed multiple parameters in the inference. We set the parameters ${t}_{0},{R}_{0},{V}_{0},{V}_{\mathrm{sw}},{b}_{s}$ and σ to their true underlying values. With the exception of the σ parameter, these parameters can be estimated using coronagraph or heliospheric imagers for any real event.

Table 1.  Fiducial Initial Parameters Used in Our ABC-SMC Mock Example

L1	L2	O	${D}_{1\mathrm{au}}$	δ	R0
20°	17°	315°	$0.165\,\mathrm{au}$	2.68	$20\,{R}_{\odot }$
V0	τ	${B}_{1\mathrm{au}}$	Vsw	Γ	σ
675 km s−1	3.15	17.2 nT	368 km s−1	0.65	1 nT

Note. Further parameters, not shown in the table, are the initialization distance at 20 solar radii and the initial launch time, which is set to 2018 January 1 00:00 UTC.

Download table as:  ASCIITypeset image

The fitting is performed on a 16 hr interval, with the individual fitting time points and the corresponding noisy measurements shown as dots in Figure 4, using a total of 16 fitting points. The t_S and t_E markers are set to exactly 30 minutes before/after the flux rope is measured. The initial threshold value is calculated as ${\epsilon }_{0}=12.3\,$ nT, which is reduced to the final value of ${\epsilon }_{16}=1.7\,$ nT after 17 iterations. The α hyper parameter was set to 0.25 and the number of particles N per iteration was set to 4096.

Figure 5 shows a so-called scatter-plot matrix of the remaining free parameters. From this result we can see that our algorithm is able to accurately infer the correct range where the ground truth is located. The parameters that control the orientation ${L}_{1},{L}_{2}$ and O are inferred with an extremely unrealistic accuracy of only a few degrees. The two parameters with the lowest accuracy are δ and τ. As can be seen in the δ − τ 2D scatter plot, these two parameters are not independent, i.e., there exists a degeneracy. This is most likely an artifact of our magnetic field model because of the twist being modified by the factor $E(\delta )$. This degeneracy will severely limit the accuracy with which we can infer the twist parameter in flux ropes. Given a fixed δ parameter, possibly determined from auxiliary measurements, the accuracy on τ could be considerably improved.

Zoom In Zoom Out Reset image size

As the ABC-SMC algorithm probes the entire parameter space of our model, the resulting uni-model probability distributions for the parameters show that are no multiple local minima, at least within the achieved threshold of the summary statistic, for this particular synthetic event. While there is no absolute guarantee that this is actually the case, we were able to gather evidence for the synthetic case by running the ABC-SMC algorithm multiple times with different random seeds (that control how the random samples are drawn).

In this section we apply the 3DCORE model and the associated ABC-SMC algorithm to an observed flux rope measurement. For this purpose we select an event that was captured by PSP on 2018 November 11–12 shortly after its first flyby at the Sun at approximately 0.25 au. This event represents the in situ observation of a CME magnetic flux rope at the smallest heliocentric distance in space history. A recent detailed study of this event was also presented in Nieves-Chinchilla et al. (2020).

Figure 6 shows the spacecraft and planetary constellation on 2018 November 12 00:00 UT. PSP was positioned almost exactly at the backside of the Sun as seen from Earth, at 1786 heliospheric longitude, while STEREO-A was close to quadrature, at −1028 longitude from Earth, which is a favorable position for imaging CMEs that are either Earth directed or on the backside as seen from Earth.

Zoom In Zoom Out Reset image size

Figure 7 shows an image from STEREO-A's COR2 coronagraph at November 12 02:09 UT and the corresponding Jplot. The coronagraph shows a small CME structure propagating away from the Sun within the ecliptic at around 10 R_⊙ (left side on the image). Using the Jplot we estimated the CME velocity to be approximately 280 ± 20 km s⁻¹ at 15 R_⊙.

Zoom In Zoom Out Reset image size

Figure 8 presents the in situ measurements for the magnetic field, the proton bulk velocity, and the resulting β_p coefficient for the event around 24 hr after the CME was observed by COR2 on STEREO-A. The magnetic field was measured by the FIELDS instrument (Bale et al. 2016), and the proton bulk speed, density, and temperate by SWEAP (Kasper et al. 2016). For our analysis we used the publicly available L2 data set for the magnetic field and the L3 moments for the proton measurements. In the case of the magnetic field, we further applied smoothing in the form of a Gaussian kernel of width σ = 900 s (solid line).

Zoom In Zoom Out Reset image size

The measurements show clear characteristics of a transient flux rope event, in the form of a significantly enhanced magnetic field and a very small β_p coefficient. The magnetic cloud can be identified by inspecting the β_p coefficient, which describes the ratio in between thermal and magnetic pressure. We use the extremely small β_p value, lower than 0.25 for most of the period in between 23:52 UT and 6:13 UT, as the indicator for the magnetic cloud interval. The coronal mass ejection itself is slow, with an estimated speed of $425\,\mathrm{km}\,{{\rm{s}}}^{-1}$ at the front edge. During its propagation over PSP, the speed decreases to an average of $387\,\mathrm{km}\,{{\rm{s}}}^{-1}$, only slightly faster as the solar background wind speed of approximately $350\,\mathrm{km}\,{{\rm{s}}}^{-1}$ in the hours preceding the event.

As is described in detail in Nieves-Chinchilla et al. (2020), there are indications that the local measurements show two interacting structures. This may explain the strong distortion that is present in the observed flux rope, which undergoes little change from 23:52 until around 04:00, at which point the magnetic field changes very rapidly. Because of the limited scope of our study and the applicability of our model, we will assume that there is only one single structure.

We apply the ABC-SMC algorithm to the PSP magnetic field measurements using the time range 2018 November 12T01:00–2018 November 12T06:00 with 11 equidistant fitting time points at half-hour intervals. We furthermore set the t_S and t_E markers to 2018 November 12T00:00 and 2018 November 12T07:30 respectively, allowing for flux rope solutions with up to 7.5 hr in duration. The CME initialization is set at 2018 November 11T06:00 at a distance 15 R_⊙ from the Sun with a fixed initial velocity of 280 km s⁻¹. These are the only initial conditions that we use for our algorithm, and all other parameters, including the noise level, are described by flat priors within sensible ranges.

The hyperparameters for the algorithm were chosen similarly to the synthetic example in Section 3.4. The particle count was doubled to N = 8192, because a higher number of free parameters requires a higher resolution. The adaptive threshold value α was set to an aggressive value of 0.25. The initial threshold value is computed to be around ${\epsilon }_{0}=80\,$ nT.

Figure 9 shows the resulting ensemble fit of the magnetic field measurements. We find that we are able to generally reproduce the measurements very well within the first half of the flux rope until 04:00 UTC. Beyond that time the fit starts to diverge, especially the B_N component, which only reaches a maximum of 50 nT instead of the measured 80 nT. Furthermore, all members of the ensemble are significantly longer in duration and only end at 07:30, which was the defined t_E marker. These issues with fitting the latter part of the flux rope are clear indicators of the distortion that is present in the measurement.

Zoom In Zoom Out Reset image size

In Figure 10 we show the inferred 1D and 2D posterior distributions in the form of a scatter-plot matrix. The parameters that define the propagation direction and orientation of the CME are given in HCI coordinates. For most parameters the confidence intervals are similarly low as for the synthetic example, but as the flux rope fit itself still contains a significant error in addition to the distortion, these results should be treated with some suspicion. The confidence intervals on the direction and orientation parameters are lower than 10 deg. The results for the latitude L = 8 ± 5 deg and the orientation parameter O = 198 ± 4 deg are particularly interesting, as they can be roughly verified using the coronagraph images from STEREO-A from Figure 7. By comparing our results with the coronagraph images, we find that our inferred latitude is within the same region, as the structure in the coronagraph is at most propagating at an angle of 10 deg to the ecliptic. The inferred orientation parameter indicates that the structure should be tilted at an angle of around 15–20 deg to the ecliptic, which is harder to verify. This type of comparison could be made clearer with results from the graduated cylindrical shell model, if images of the CME existed at higher altitudes. The ${D}_{1\mathrm{au}}=0.26\pm 0.01\,\mathrm{au}$ estimate is most definitely too large by at least 10%–15% due to the 1.5 hr time difference between the actual flux rope event and the fits. The ${B}_{1\mathrm{au}}=11.5\pm 0.8\,$ nT parameter indicates that the CME was not particularly magnetically strong under the assumption that the scaling relations up to 1 au hold as they are defined in the model.

Zoom In Zoom Out Reset image size

Some of these results are put into a better perspective in Figure 11, where we show the structure of the fit with respect to the planetary and satellite positions. The parameters for this single representative 3DCORE run are the medians of the ensemble solution generated by the ABC-SMC algorithm (except for the δ parameter, which was set to 2 for better visualization).

Zoom In Zoom Out Reset image size

An interesting result in this particular analysis are the values for the cross-section aspect ratio δ and the magnetic field twist τ. When defining the prior $p(\delta )$, we limited the maximum value to δ = 8.5. As we can see from Figure 10, the results prefer an extreme aspect ratio. As was already observed in more detail in Section 3.4, there is an inverse relationship between the aspect ratio δ and twist τ. A large δ value is generally combined with a low magnetic field twist. The total number of twists over our torus-like shape is estimated to τ = 1.73 ± 0.5. At 1 au this value corresponds to around 0.6 twists/au, which is rather low. Both the aspect ratio and twist parameters are most likely the least accurate results in the analysis of this event. There is likely to be a strong bias due to the uniform priors that we have used. Because our model is semi-empirical and no physical simulations are used to generate the results, there is no inherent weighting for any of the parameters. This effect may be the most visible for the ${D}_{1\mathrm{au}}$ and the δ parameter (and therefore indirectly also the twist). Large ${D}_{1\mathrm{au}}$ and δ parameters correspond to large CME structures that therefore also have a higher chance of hitting any observer. Therefore when analyzing any event, there will always be a strong bias toward larger structures, because they have a large probability of having generated the measurement when assuming that small structures have the same occurrence rate as larger ones. This is most certainly not true, as one would expect larger CMEs to be rarer than their smaller counterparts. This shows that in the cases where the measurements are not fully conclusive, the priors can have a strong effect on the end result.

The second oddity is the result for solar wind speeds V_SW that are unusually high, with values over 470 km s⁻¹. This value is larger than any speed that is measured locally by PSP at any time just before, during, or just after the CME of interest. This is a strong indicator that our chosen initial CME speed of 280 km s⁻¹ is strongly underestimated. A simple explanation that can lead to such an underestimation is the viewing angle of STEREO-A. The longitudinal position of STEREO-A was approximately 230° on November 12 (HCI); as such the viewing plane in which we measured the initial velocity was at 140°. Assuming that our inferred value of ${L}_{1}=168^\circ$ holds true nonetheless, we would have underestimated the initial speed by 12% (almost 40 km s⁻¹). Proper handling of this issue would require a coupling of the L₁ and V₀ parameter in our fitting algorithm or in the model.

In this paper we have presented an improved version of the 3DCORE model and implemented an ABC-SMC fitting algorithm that we combined with our model to fit magnetic field flux rope measurements. The primary improvement with respect to the original 3DCORE model (Möstl et al. 2018) was a significant increase in the number of simulation runs that can be performed allowing the usage of a more computationally expensive fitting algorithm. This was achieved due to the definition of the geometrical model shape using a separate coordinate system, which allows for efficient collision detection and evaluation of the internal analytical magnetic field. On the implementation side we made heavy use of parallelization and created a data generation pipeline that operates on multiple simulations simultaneously. In general our 3DCORE implementation is fastest when using thousands of simulations at once.

The computational efficiency is important for two reasons. First, it could in the future open up the possibility of real-time analysis, which ideally requires simulations that run on the minute scale. Second, it will allow us to further generalize or make the model more complex without losing its usability. This will be of particular interest in future studies incorporating a more sophisticated shape or propagation model that includes a solar wind model.

We tested our ABC-SMC algorithm on a synthetic data set and found that under ideal conditions, it is capable of extracting the known ground truth. This assumes that we know the exact boundaries of the flux rope signature in time and have a correct statistical model of the noise or fluctuations that occur in the measurements. This algorithm also allows us to estimate the constraints on the model parameters, which means we can determine the accuracy of our results. In the case of the synthetic test results this worked very well, although we still found that there were significant errors on the cross-section aspect ratio δ and magnetic field twist τ parameters. We also found that there is a relationship between the pairs of each (δ, τ) solution, where a large aspect ratio corresponds to a small twist value and vice versa. Whether this is just an artifact of our magnetic field model assumptions or if it still persists for more physically accurate models will have to be discussed in the future. If this degeneracy persists, it shows that estimating any of these two parameters independently will generally fail and auxiliary measurements will be required.

While no big issues appeared in our synthetic test case, the picture changes drastically once we apply our fitting technique to real-world data, as in the PSP event that we presented. The model is no longer fully capable of describing the measurements, particularly the distortion of the flux rope toward the end. The assumption of Gaussian random noise is also no longer adequate. An issue that is not apparent from the results is that the algorithm is very sensitive to how we set up the fitting points, especially the t_S and t_E markers. The algorithm is more efficient, in the sense of the accepted-rejection ratio, for larger intervals defined by t_S and t_E, because fewer simulations are sorted out due to generating signatures that are too long in duration. In general, setting the t_S and t_E markers to the exact positions of the start and end of the flux rope is impractical as the algorithm will be too slow. A larger window will lead to an overestimation of the ${D}_{1\mathrm{au}}$ parameter and may also have a non-negligible effect on others.

In all our analyses we also always used uniform priors out of convenience. In the synthetic case this is no issue because the data are fully self-explanatory and can be modeled completely by our simulation. In the more complicated case, when less information can be extracted from measurements, the priors start having a significant effect on the end result. For these cases, which will most likely include most observed flux ropes, one has to construct more useful priors. In particular there should be a lower weighting for larger CME structures, which are inherently rarer than smaller ones. Such priors could be constructed from CME catalogs that contain well-behaved flux rope signatures (e.g., from HELCATS⁵ ).

The novel component of our work is the first application of Bayesian analysis with respect to CME flux ropes. For this purpose we implemented an ABC-SMC algorithm with which we can approximate the posterior distributions of all model parameters when fitting the model to observations. This approach, if done correctly, can have multiple advantages over the simpler fitting methods for flux ropes that are primarily used today. First and foremost the usage of this class of algorithms allows us to estimate the intrinsic errors on the inferred model parameters. This gives us information on how well we are able to determine certain parameters or if they are determinable within a reasonable accuracy at all, as they may be obscured completely due to projection effects. This can furthermore allow us to compare different models and various techniques in terms of accuracy and not only in the final result, which can be ambiguous as the ground truth is unknown.

Furthermore, our fitting algorithm is easily extendable, without major changes, to the case of multipoint event analysis. In this context it is possible to investigate if and how multiple measurements agree or disagree with one another. It is additionally possible to measure, in absolute terms, the information gain that is achieved when continuously adding independent measurements to the analysis. One could, for example, perform a cost-benefit analysis on how many in situ spacecraft would be required to sufficiently determine the properties of any CME flux rope.

Lastly, the construction of priors allows us to easily incorporate additional information into our analysis. This extra data can come from other measurements, physical considerations, or previous analysis results from our algorithm. It also allows us to continuously update any result whenever new measurements or model improvements become available.

In summary, we have developed an alternative path for fitting CME in situ flux rope measurements and presented the first results. In the near future, we aim to further develop our technique by improving the underlying 3DCORE model, constructing a better noise model and defining better initial priors. The next logical step for our technique is to test its self-consistency when analyzing multipoint events. This will hopefully allow us to gain new insights into both our model and technique and the overall structure of CMEs.

C.M., A.J.W., R.L.B., M.A.R., and T.A. thank the Austrian Science Fund (FWF): P31521-N27, P31659-N27, P31265-N27, J4160-N27. We acknowledge the NASA Parker Solar Probe Mission and FIELDS team, led by S.D. Bale, for use of data. We acknowledge the NASA Parker Solar Probe Mission and SWEAP team, led by Justin Kasper, for use of data. The open source code for the 3DCORE model is available at https://github.com/ajefweiss/py3dcore, and the heliosat package to download and process spacecraft data is accessible at https://github.com/ajefweiss/heliosat. The Python scripts that were used to generate the figures and perform analysis in this paper are partially available as Jupyter notebooks at https://github.com/helioforecast/Papers/tree/master/Weiss2020_3DCORE.