© A. J. Weiss et al. 2021

1. Introduction

Coronal mass ejections (CMEs) are a highly energetic process in which a large amount of magnetized plasma is ejected from the Sun. Interplanetary CMEs (ICMEs), which can be observed and identified using heliospheric imagers, propagate far into the heliosphere and arrive at the inner planets within a few days (Möstl et al. 2009b; Rouillard 2011). On Earth, they are responsible for the most extreme geomagnetic storms that can occur and are therefore the strongest possible driver for the Earth’s space weather system (Gosling et al. 1991; Farrugia et al. 1993; Gonzalez et al. 1999; Schwenn 2006). As such, they represent one of the most extreme events that must be taken into account when relying on modern applications or equipment that are affected by the space weather environment (e.g., Fuller-Rowell et al. 1994; Bolduc 2002; Pulkkinen et al. 2017).

While the geomagnetic effectiveness of an ICME arriving at Earth relies on multiple factors, the two most important indicators are the speed and the north-south magnetic field component B_z. Given a strong southward orientation of the magnetic field, an ICME can facilitate the injection of energy into the magnetosphere, leading to the aforementioned geomagnetic events (Burton et al. 1975; Gonzalez et al. 1994). The in situ magnetic field signatures of ICMEs often exhibit similar structure with a frontal shock and sheath region that is followed by a highly structured rotating magnetic field (Burlaga et al. 1981; Klein & Burlaga 1982; Kilpua et al. 2017). This rotating magnetic field, which is commonly referred to as the magnetic obstacle, can potentially manifest itself as a strongly south orientated magnetic field over a duration of many hours. As such, the understanding of the structure of these magnetic obstacles is key to linking ICMEs with the resulting geomagnetic effects (Bothmer & Schwenn 1998).

The magnetic obstacles are often described as magnetic flux ropes (MFR) structures (Marubashi 1986; Burlaga 1988; Lepping et al. 1990). While alternative pictures exist, based on spheromak models (Vandas & Romashets 2017; Singh et al. 2020), the most common depiction uses some form of a bent flux tube that remains connected to the solar surface while propagating away from the Sun (Dasso et al. 2005). A variety of different successful magnetic flux rope models have been developed for the purpose of describing the in situ signatures of ICMEs, which generally assume a rigid cylindrical or toroidal geometry. They are furthermore divided into two general classes, those based on uniform-twist fields (Gold & Hoyle 1960; Hu et al. 2015; Vandas & Romashets 2017) or constant α fields (Lundquist 1950; Lepping et al. 1990; Hidalgo et al. 2002; Nieves-Chinchilla et al. 2018).

Along with the wide variety of descriptions for the magnetic field, the global shape of ICMEs also remains largely hidden. While the initial structure can be partially observed due to Thomson scattering close to the Sun (Cremades & Bothmer 2004; Thernisien et al. 2006), the in situ magnetic field measurements represent a strong projection of the internal magnetic field that is intertwined with the global or local geometry. This can make it exceptionally difficult to draw conclusions on the shape using just single in situ magnetic field measurements. It is broadly understood that there are various processes that influence the shape and evolution of ICMEs during their propagation through the heliosphere (Manchester et al. 2017). This includes rotations or deflections close to the Sun (Lugaz et al. 2012; Kay et al. 2015) and interactions with the ambient solar wind (Riley & Crooker 2004; Liu et al. 2006; Démoulin & Dasso 2009). More complex interactions are also possible, such as CME-CME collisions (Möstl et al. 2012; Lugaz et al. 2017).

More recently, a focus has been placed on the study of multi-point events, namely, ICMEs that are observed in situ by multiple spacecraft (Möstl et al. 2009a; Kilpua et al. 2009; Good et al. 2019; Davies et al. 2020; Salman et al. 2020) in an attempt to gain a better global understanding of ICME structure and their evolution. Unfortunately, such events are inherently rare due to the limited number of spacecraft with the onboard magnetometers in different solar orbits within similar time frames. The successful modeling of these multi-point events invariably requires global ICME models that include magnetic fields and ICME evolution, which goes beyond the capability of analytical models with rigid cylindrical or toroidal shapes. While it may be tempting to use large-scale magnetohydrodynamic (MHD) models (e.g., Manchester et al. 2008; Török et al. 2018; Verbeke et al. 2019), they are too computationally expensive and unwieldy for analyzing the magnetic field measurements of CME events. For this purpose, multiple semi-empirical models have been developed (Isavnin 2016; Möstl et al. 2018; Kay & Gopalswamy 2018; Rouillard et al. 2020; Weiss et al. 2021), which are aimed at mimicking the general properties that ICMEs are thought to have. These models serve as a bridge between the more complex MHD models and the fully rigid analytical models and can include simple interaction with the solar wind, some form of deformations, and an internal magnetic field while retaining some of the computational simplicity of the simpler analytical models.

In our previous paper (Weiss et al. 2021), henceforth referred to as W21, we showcased the successful application of the semi-empirical 3D coronal rope ejection (3DCORE) model to a single point event detected by Parker Solar Probe. The 3DCORE model is a semi-empirical forward simulation model that describes an ICME as a self-similarly expanding torus-like structure that remains attached to the Sun. It includes ICME evolution factors such as flux rope expansion and deceleration due to slower solar wind which are described using empirical relations. The internal magnetic field structure is based on an analytical toroidal uniform-twist solution described in Vandas & Romashets (2017).

In this paper, we attempt the novel approach of simultaneously fitting multiple magnetic field measurements from different spacecraft at different locations. So far, the majority of multi-point studies have only compared different single point fits from different spacecraft from the same underlying event (Möstl et al. 2012). The results from our analysis are aimed at showing whether it is even possible to describe multiple observations using the global assumptions within our model and, if not, to provide hints regarding its deficiencies. The comparison of the single-point fits and multi-point fits of the same event also allows for a sort of cross-validation that can further the general understanding of the underlying event and the applied model.

Section 2 introduces the two events used for testing our approach. The first event, the April 19 ICME, represents three very nearby, high-quality measurements for which we would expect our approach to have the highest chances of success. The second event, the May 28 ICME, presents a bigger challenge and was primarily chosen to test the limits of our methods. In Sect. 3, we give a brief recap of the 3DCORE model from W21 and also describe the two changes that were made to update the model. Section 4 describes the adapted ABC-SMC fitting algorithm that was used for fitting the multi-point data. Section 5 shows the reconstructed magnetic field measurements for both events and the inferred model parameters from multiple spacecraft combinations. Finally, in Sect. 6, we discuss the results and the conclusions we can draw.

2. Data and events

To demonstrate the applicability of our 3DCORE model and the associated Monte-Carlo based fitting method for modeling the magnetic field of multi-point ICMEs, we choose two recent events that were both observed in situ by the Solar Orbiter (SolO), Wind, and BepiColombo (Bepi) spacecraft.

The first and primary event that we make use of is the 2020 April 19 ICME (henceforth the April event). This event was initially detected in situ by Solar Orbiter at a heliocentric distance of 0.80 AU on 2020 April 19 05:06 UTC, at about 4 degrees longitude (HEEQ) east of Earth. The ICME shock subsequently arrived a day later at Wind and BepiColombo on April 20 at 01:34 UTC and 03:09 UTC, respectively. For an extensive description of this event, we refer to Davies et al. (2021) as here we only focus on the details that are of interest to our analysis.

The second event that we chose is the 2020 May 28 ICME (henceforth the May event). Similarly to the April event, this ICME was again first detected by Solar Orbiter at a significantly smaller heliocentric distance of 0.55 AU on May 28 15:00 UTC. The ICME was again also observed by Wind and BepiColombo at a very large longitudinal separation on May 29 at 15:15 UTC and May 30 at 01:15 UTC respectively. Due to the lack of clear indicators of a shock front, these arrival times represent rough estimates.

Figure 1 shows the in situ measurements from the SolO MAG magnetometer (Horbury et al. 2020), the Wind MFI magnetometer (Lepping et al. 1995), the Bepi MPO-MAG magnetometer (Heyner et al. 2021), and the Wind SWE plasma instrument (Ogilvie et al. 1995) for both events. The magnetic field measurements are shown in their respective local radial-tangential-normal (RTN) coordinate systems. For the plasma measurements from the Wind spacecraft, we only show the proton density, N_p, and the associated computed plasma beta coefficient, β_p. The blue shaded area represents our estimated intervals for the unperturbed magnetic flux ropes (UMFR, see Davies et al. 2021 for the April event), which are the intervals we use for our fitting analysis. The UMFR intervals were estimated using both the magnetic field measurements and the plasma measurements and matched with each other using features in the magnetic field data.

Fig. 1. In situ magnetic field measurements in local radial tangential normal (RTN) coordinates from Solar Orbiter (a/e), BepiColombo (b/f), and Wind (c/g) for the 2020 April 20 (left) and 2020 May 28 (right) event. The two bottom panels (d/h) show the in situ proton density number Np and the associated plasma beta coefficient, βp. The blue shaded area represents our estimated intervals for the unperturbed magnetic flux ropes (UMFR) which we subsequently use in our fitting approach.

Figure 2 shows the different spacecraft positions within the heliosphere on 2020 April 20 00:00 UTC and 2020 May 29 00:00 UTC including other satellites, such as Parker Solar Probe (PSP) and the inner planets. As can be seen from these two figures, the STEREO-A spacecraft, which carries heliospheric imagers (HI, Eyles et al. 2009), was ideally positioned at both times to image and track Earth-directed CMEs. The captured images for both events are briefly discussed in Sect. 2. Unfortunately, PSP was located on the far side of the Sun and therefore incapable of recording any of these two events. Table 1 summarizes the positions of all three spacecraft that observed the events in situ on the same dates in heliocentric Earth equatorial (HEEQ) coordinates.

Fig. 2. Spacecraft positions for Solar Orbiter, BepiColombo, Parker Solar Probe, and STEREO-A on 2020 April 20 00:00 UTC (left) and on 2020 May 29 00:00 UTC (right) in HEEQ coordinates. Furthermore shown are the planets Mercury (brown), Venus (orange), and Earth (green). Colored trails show orbital movements of the spacecraft within a time range of two months around the specifically chosen position. Two open-ended cones show the field of view (FOV) of the HI-1 (20°) and HI-2 (70°) imaging instruments aboard STEREO-A.

The magnetic obstacle of the April event appears as an almost pristine flux rope in all three in situ measurements with a nicely rotating magnetic field. The normal component, B_N, is bipolar, while the transversal component, B_T, unipolar and the radial component, B_R, are negative near zero. According to the standard classification schemes this corresponds to a low inclination left-handed south-east-north flux rope (Bothmer & Schwenn 1998; Davies et al. 2021). Due to the textbook- level structure of this flux rope it represents a unique opportunity to test the validity of our model with multiple nearby in situ measurements. The May event shows significantly more complex structure and was primarily chosen to identify the limits of our modeling approach at extremely large longitudinal separations. While we can still clearly identify a right-handed south-west-north (SWN) flux rope in the SolO MAG measurements, which represents a low inclination flux rope with the opposite axial direction and handedness compared to the April event, the magnetic field is strongly perturbed at the center for Bepi and Wind. Nonetheless, taking only the endpoints of the Bepi and Wind measurement into account the magnetic field still exhibits a magnetic field rotation that is comparable to Solar Orbiter. The differences can be explained due to large spatial distance in between the spacecraft and poses a significant challenge for our modeling efforts. We partially alleviate these issues by applying a strong smoothing onto the magnetic field data from the May event to suppress the smaller scale perturbations and only using the less perturbed endpoints for our analysis.

Figure 3 shows two white light images taken by the HI-1 instrument on 2020 April 16 07:29 UTC (left) and 2020 May 26 15:29 UTC (right). The HI-1 instrument was able to nicely follow the evolution of the April CME over a time range of almost 24 h showing a faint leading edge, an elliptical cavity and a strong density pile up at the trailing edge. The CME cavity, which represents the magnetic obstacle of interest, is seen edge-on and appears to be of elliptical shape with a time-invariant aspect ratio of around 2 (see Davies et al. 2021). These images for the April event appear to be consistent with the approximation of an elliptical cross-section that is one of the central assumptions built into our model. The white-light images for the May event, on the other hand, show a smaller and fainter structure propagating above the ecliptic. The structure also darkens considerably faster making it impossible to observe the evolution of the shape for a longer duration than a mere few hours.

Fig. 3. White-light images, with the stars subtracted using a median filter, from the STEREO-A HI-1 instrument showing the April event on the 16 of April 07:29 UTC (left) and the May event on the 26 of May 15:29 UTC (right). The image of the April event shows a slightly flattened Earth-directed CME with a faint leading edge, an elliptical cavity, and a strong density pile up at the back. The image of the May event shows a structure at a higher latitude that is hard to identify in a single image (marked with the yellow arrow) and does not contain an obvious cavity.

Both the April¹ and May² event are part of the HELCATS catalog. Despite their only being observed from a single vantage point, we can still extract rough estimates for the propagation speed and direction by tracking the CME evolution. Under the assumption of self-similar expansion, the given values for the April event are an average CME speed of 339 ± 11 km s⁻¹ and a propagation direction of −2deg longitude and −6deg latitude in HEEQ coordinates. For the May event, we obtain an average speed of 497 ± 64 km s⁻¹ and a propagation direction of 36deg longitude and 11deg latitude (HEEQ).

3. Model

In this paper, we use an updated version of the 3DCORE model to analyze our two selected events. We give a brief general overview here and then go on focus on the changes to the implementation that we have made in order to facilitate the analysis of multi-point events. The two primary changes that are worth mentioning are the “twist reparametrzitation”, where we substituted the magnetic field twist parameter by using a different related parameter that has beneficial numerical properties and the introduction of a new noise model. For more, we refer to W21 for a detailed description of the original implementation.

3.1. 3DCORE overview

The 3DCORE model is an empirically based forward simulation model which can generate synthetic in situ magnetic field measurements anywhere within the heliosphere. It assumes a self-similarly expanding torus-like MFR structure that propagates along a fixed direction radially away from the Sun. The end points of the torus-like structure stay attached to the Sun. The expansion of the flux rope structure is mimicked using a simple scaling relation. The interaction with the background solar wind is described by using a drag-based model (Vršnak et al. 2013), where the frontal part of the structure is accelerated or decelerated, depending on the relative speed to the ambient solar wind speed. The solar wind speed is assumed to be constant everywhere. The “pancaking” effect, described by Riley & Crooker (2004), is approximated by stretching the cross-section of the torus into an elliptical shape. An example of the 3DCORE model is shown in Figure 4.

Fig. 4. Example of a 3DCORE simulation showing the torus-like structure at near 1 AU with a highly elliptical cross-section and two magnetic field lines at differing radial distances (twist number τ = 5).

The internal magnetic field of the MFR is simulated by inserting an analytical magnetic field model into the structure. We use a model based on the solutions described in (Vandas & Romashets 2017), which we have slightly modified in order to accommodate the elliptical cross-section. We furthermore make the assumption that the total number of twists along the torus stays constant over time, which leads to a decreasing twist per unit length as the structure expands. The decrease of the magnetic field strength as the structure propagates forward and expands is again described by a simple scaling relation (Leitner et al. 2007). While the implemented magnetic field is based on analytically derived equations, it is important to note that it is not fully physical and violates various conservation laws. We assume that the errors that are introduced by this approximation are less severe than those arising from the basic assumptions on the geometry.

One of the key points of our model is that the simulation is extremely fast as it is only based on empirical relations and does not solve any expensive and complicated MHD equations. This efficiency makes it significantly easier to explore the full parameter space of the model and allows for the usage of Monte Carlo-based methods, which require on the order of thousands or millions of simulation runs.

In total the model contains 12 controllable parameters: the initialization time t₀, propagation direction and orientation, lat,  lon,  inc, cross-section width, d_1 au, cross-section aspect ratio δ, initialization distance, r₀, initial speed, v₀, magnetic field strength b_1 au at 1AU, twist number τ, global solar wind speed, v_sw, and drag parameter Γ. Both the diameter of the cross-section and magnetic field strength are modified by scaling relations with d(t) = d_1 aur(t)^1.14 and b(t) = b_1 aur(t)^−1.68, where r(t) is the distance of the CME front from the Sun.

3.2. Twist reparameterization

One of the key findings of W21 was that we were not able to accurately infer the magnetic field twist τ and the aspect ratio of the torus cross-section δ using single-point observations. We found a strong τ − δ degeneracy, meaning that these two parameters are strongly correlated. The dependency relation roughly follows τ ∝ δ/E(δ), where E(δ) is the function for the circumference of an ellipse with a given aspect ratio.

This issue persists as we are not directly capable of measuring the twist of the magnetic field and, instead, we can only measure the local direction of the magnetic field vectors. Depending on the geometry or size – in this case the height of the CME structure given by the aspect-ratio – the total number of twists along a fixed length can change drastically. On top of the issue of not being able to accurately infer these model parameters simultaneously we found this degeneracy to be a numerical issue when using our Monte Carlo-based fitting algorithm, as the degeneracy is strongly non-linear and is inefficiently sampled using multivariate Gaussian transition kernels. The shape of this degeneracy, as seen in a two dimensional histogram, closely resembles a “banana shape” and presents a classical problem for the class of algorithms that we use (Haario et al. 1999; Filippi et al. 2013).

We therefore attempt to at least remove the non-linearity of this degeneracy by defining a new twist parameter T_τ = τE(δ). The introduction of this new parameter significantly reduces the numerical instability of our approach and increases the overall speed of the fitting algorithm as the sampling efficiency is increased. The original τ values remain recoverable in a post-processing step.

3.3. Noise model

In order to estimate the uncertainties within our model with our fitting algorithm, we need to be capable of simulating the measured noise. In the context of our model, we gather all unwanted small scale components of the magnetic field under the general term of “noise”. This includes the expected instrument noise but also other seemingly random phenomena such as MHD turbulence and waves or fluctuations due to undetermined physical processes. In the case of instrument noise, we assume that the data sets are thoroughly cleaned and that the most significant spacecraft influences are correctly removed.

In W21, we made the assumption of a simple additive Gaussian noise model that was described by a single parameter σ, which determined the standard deviation of the additive noise component. As was briefly mentioned in W21, this is a very basic approach that is easy to implement but it fails to fully capture the small scale behaviour of the magnetic field that is seen in the observations. We therefore go on to introduce a more sophisticated additive noise model based on the power spectrum distribution (PSD) P(k). Our implementation for generating the noise is inspired by the analogous process that is used to generate primordial perturbations for dark matter N-body simulations using the matter power spectrum (see Hahn & Abel 2011).

For our purpose, we are interested in using the PSD to reproduce, within our model, the type of fluctuations that are seen in the actual observations. As a first step we require an appropriate PSD from which we generate the noise. The most adequate choice is the PSD from the underlying event which is being analysed as the distributions can vary significantly in between different events. In our case, we ignore cross-correlations between the three different components of the magnetic field, which significantly simplifies the process as we only need to generate three independent noise samples for each component. Given a random noise field μ(t)∼𝒩(0, 1) and its corresponding Fourier transform $\stackrel{\sim }{\mu }(k)=\mathcal{F}(\mu (t)),$ we can generate a random time series:

$$\begin{array}{c}\hfill \delta (t)={\mathcal{F}}^{-1}\left(\sqrt{P(k)}\stackrel{\sim }{\mu }(k)\right),\end{array}$$(1)

so that $\stackrel{\sim }{\delta }(t)\approx P(k)$. This process allows us to generate random and unique noise fields which all have the same underlying PSD statistic. In the final step, we add this random field to our model output to combine the large scale magnetic field rotation with the random small scale fluctuations.

We compute the PSD for each magnetic field component independently using Welch’s algorithm (Welch 1967) from the time interval of interest. In order to remove the large-scale fluctuations that partially describe the magnetic field rotation described by our flux rope model, we pre-process the measurements using a linear de-trending process on an hourly scale.

As the PSD varies only slightly for each component, we generate a single general distribution from the average over all three components. In the next step, we use the computed PSD to generate a random sequence with the same underlying PSD that acts as our noise. As white noise (Gaussian noise) has a flat power spectrum, we can generate the desired sequence by first generating a random Gaussian sequence δ ∼ 𝒩(0, 1) and multiplying the newly generated sequence in Fourier space with P(k) before transforming the sequence back into real space. This process is performed independently for each component of the magnetic field on top of which the noisy δ sequence is then added.

The biggest advantage of this newer noise model is that it eliminates the necessity for a noise level parameter as the severity of noise is directly estimated from the measurements. In theory, this approach also allows us to fit profiles with an extremely high density of fitting points since the point-to-point fluctuations are more accurately modeled (decreasing with increasing proximity). In practice, this is not done as the inaccuracies are dominated by the assumption of the geometry and magnetic field model.

While this noise model is expected to be a large improvement over the model employed in W21 it is still important to note one of its deficiencies. The fluctuations in the magnetic field strength, which represent compressive fluctuations, are generally weaker than those of its components that describe incompressible fluctuations (e.g., Good et al. 2020). This is a strong hint that the underlying fluctuations would be better described using rotations, rather than the additions that we use in our noise model. However, since the summary statistic we use for our fitting algorithm only makes use of the magnetic field components and not the total field, we expect the impact to be negligible.

4. Multi-point fitting algorithm

In this section, we focus on the details of the fitting algorithm that we use to analyze the magnetic field measurements with our 3DCORE model. The underlying algorithm, an approximate Bayesian Computation sequential Monte Carlo (ABC-SMC) fitting algorithm, remains largely unchanged from W21 and requires only minimal adaptions for the multi-point scenario.

The general idea behind the ABC-SMC fitting algorithm is to generate large ensemble simulations from an initial parameter space and classify the goodness of fit for each simulation run according to an error metric, called the summary statistic, and reject all simulations that fall above a certain error threshold. An important step is also the rejection of all simulations that do not generate valid results at all due to the underlying CME geometry not hitting any of the observers or producing magnetic field signatures that are significantly time-shifted due to the simulated CME arriving too late or early.

By repeating this procedure and using the accepted simulations as a representation of the initial parameter space for the next iteration, we can sequentially reduce the error threshold starting from a larger value and arrive at a good estimate of the underlying posterior after a certain number of iterations.

Due to the stochastic nature of the simulations, as we include the noise model introduced in Sect. 3.3, and the correct implementation of an ABC algorithm, we are able to interpret the generate ensembles as probability distributions (in the Bayesian view). This naturally gives us not only best-fit model parameters but also non-Gaussian error estimates and parameter correlations, which can further deepen the understanding of the results and the underlying model.

We discuss the changes that we have made to the summary statistic and the pre-filtering process for determining whether a simulation run can be considered as valid or invalid. We also take a brief look at the inherent biases that may appear in our analysis.

4.1. Summary statistic

For each single observation, we use a normalized root mean square error statistic between model output, x, and reference observation, y, which is defined as followed:

$$\begin{array}{c}\hfill \rho (x,y)=\sqrt{\frac{\sum _{i=1}^K{(x(t_i)-y(t_i))}^2}{\sum _{i=1}^{K}y{(t_i)}^2},}\end{array}$$(2)

where t_i is one of the K time points at which the magnetic field is to be compared.

At this point, it may appear viable to create a combined statistic that condenses multiple statistics into a single number, such as the mean RMSE over multiple observations. In practice, we have found that this approach runs into the issue in which the algorithm can make a hidden trade-off in between improving the fit for one observation at the expense of the others. For this reason, we adapted our method so that it simultaneously optimizes multiple error metrics. While this does not necessarily guarantee optimal convergence it does not allow for a scenario where one observation is optimized at the expense of another. Furthermore, the end result clearly shows if a certain observation cannot be properly fitted.

As a consequence, a simulation output for N observations x = (x¹, ⋯, x^N) is accepted if and only if ρ(x^k, y^k) < ϵ^k for all k ∈ [1, N]. The iterative threshold values e^k are computed in the same way as in W21 with $e_0^k=\rho (0,y^k)$ and $e_{j+1}^k$ is set as a quantile fraction of $e_j^k$.

4.2. Adaptive time markers

An important aspect of the fitting algorithm is to determine whether a simulation run can be considered valid in the first place. In a large ensemble simulation, a majority of runs are not expected to generate any valid results due to a lack of one or many of the defined observers or due to the producing of signatures that are strongly time-shifted when compared with the reference observation that is to be reproduced. In W21, we used two time markers, t_s and t_e, to determine valid simulation runs.

The condition that we initially set was that any simulation that we run is to be considered valid if and only if it produces a valid magnetic field signature at any of the time points, t_i, that are part of the summary statistic and produces an invalid measurement at t_s and t_e. This requires the observer to be within the CME geometry, as assumed by the 3DCORE model, at any time point t_i and to be outside of it at t_s/e. By setting t_s just before t₁ and t_e just after t_K, we are able to closely control the in situ duration of any single valid run.

The problem that arises, specifically in the multi-point scenario, is that it is initially very hard to find simulations that produce signatures of just the right length. As additional constraints are injected by adding observations, the efficiency of the fitting algorithm can suffer substantially. For this reason, we implemented an adaptive time marker scheme, in which the t_s and t_e markers are set at a large temporal distance to the fitting time points, t_i. This distance is iteratively reduced, as is the error threshold, until it reaches a final value after a certain number of iterations. This helps the sampling process considerably in the initial stages when the explored parameter space is still large. From basic tests we have concluded that this initial approximation does not have any significant negative effects on the final results.

4.3. Bias

An inherent property of our 3DCORE model is that it is an empirical model and does not contain any physical calculations. As we show, this can lead to some peculiar results, which require some additional treatment. While our fitting algorithm and error metric guarantees that we only select model runs and the associated model parameters that match the observations well, it can only make use of “valid” runs. Simulation runs that do not generate any results due to the 3DCORE shape neglecting to hit any of the predefined observers or runs that generate magnetic fields that are too strongly time-shifted are automatically rejected prior to the computation of any error metrics.

While the size of the shapes is limited due to the use of bounds for the model parameters, this nonetheless leads to an extremely strong bias towards the largest allowed shapes as they are not affected by this selection bias. The two parameters in our model that control the size of our structure are the d_1 au (torus width) and δ (torus height) parameters. The width parameter d_1 au is constrained by the duration of the in situ measurement and can therefore be more or less estimated by our algorithm but the height is completely unknown due to all observers being largely located within the ecliptic plane. Due to the δ − τ degeneracy that we mentioned previously, we know that we cannot constrain the δ parameter using a single in situ measurement. This generally leads to a linear bias of p(δ)∝δ acting as a structure that has a larger aspect ratio δ and having a larger chance of hitting the observer.

Under certain circumstances, this bias does not occur, such as a CME that is directly propagated towards the observer, which makes this issue hard to handle. In general, studies must be done on a case by case basis in looking at the model simulations that are accepted and rejected due to either a hit or miss of the observers. The conclusion is that our analysis would normally generate extremely large estimates for the δ parameter, while the correct interpretation is that nothing about this parameter is actually known. The primary hope is that this issue would be resolved with the bias eventually disappearing when adding multiple observations.

5. Results

In this section, we present the results of our study, where we fitted the April and May event using our updated 3DCORE model and ABC-SMC fitting algorithm. A constellation of the three different satellites, combined with their in situ measurements, not only allows us to create a single large analysis encompassing all three satellites but also three dual-point analyses in which a combination of two of the three available in situ measurements are used. Counting the single-point results then leads to a total of seven different possible data combinations that can be used. This allows us to study in detail how the results using our method differ when using either the full data set or smaller subsets.

Figure 5 shows the resulting triple-point fit of the 2020 April 19 CME using all three in situ measurements and a dual-point fit of the 2020 May 28 CME for SolO and Wind in local RTN coordinates. While we were capable of producing a triple-point fit for the April event we did not achieve satisfactory results for the May event. The dashed lines show one single representative fit from the generated ensemble. The shaded area around this single fit shows the 1-σ spread of all ensemble members. The vertical lines show the selected fitting time points t_i and the area in between the final t_s and t_e markers is shaded in light blue. For the adaptive time marker scheme, we used single-hour offsets for the single point results, four-hour offsets for dual-point results, and eight hours for the triple point result. In all cases, the offset was reduced by one hour per iteration of the sequential ABC-SMC algorithm.

Fig. 5. Reconstructed synthetic magnetic field measurements (dashed) with 1σ confidence interval (shaded area) using the triple point constellation for the April event and a dual-point constellation for the May event. The reconstructed magnetic field is overlayed on top of the real in situ measurements for Solar Orbiter (a/d), BepiColombo (b/e), and Wind (c/f). The blue shaded area represents the defined UMFR interval and the dotted vertical lines show the used fitting points ti. All magnetic fields are shown in local RTN coordinates.

The model initialization time t₀ is set to a specific date for both events, the initialization distance r₀ is fixed to 35 solar radii, and the initial speed v₀ and background solar wind speed v_sw are bounded within a smaller range. All other parameters are largely unbounded within a sensible range. For the April event, we set t₀ = 2020-04-15T23:00 UTC, v₀ ∈ [350, 750], and v_sw ∈ [275, 375]. For the May event we use t₀= 2020-05-26T14:00 UTC, v₀ ∈ [300, 650], and v_sw ∈ [250, 450]. The initialization times, distances, and CME initial speeds are roughly estimated from the remote observations from STEREO-A. Due to projection effects from the images. we can only estimate the lower bound for the speed. The solar wind speed is taken from in situ data prior to the CME arrivals.

For the April event, we mainly used equidistant fitting points with a one-hour separation. While this was also possible for the SolO measurement for the May event, this was no longer the case at Bepi or Wind due to the strong perturbations at the center. As can be seen in Fig. 5, we used an extremely low number of points for these two measurements that are located towards the endpoints of the flux rope signature.

The results for the April event demonstrate that we were generally capable of reproducing all three measurements simultaneously with our model. The largest discrepancy can be found within the SolO measurement where we were not able to reproduce the sloped radial component, B_R, and instead generate an almost constant field as would be expected from a rigid flux rope structure. We also do not find any issues with fitting the varying arrival times, which is a strong indication that the general assumptions on the propagation are more or less valid within the vicinity of the three spacecraft. The small jumps that can be seen in the representative ensemble members are artifacts that occur due to our magnetic field model when combined with high aspect-ratios.

The results for the May event, despite only being a dual-point fit, are significantly less accurate, as seen just from the spread. Since we only took the SolO and Wind measurements into account here, we can see that the arrival time at Bepi is incorrect, with most of the ensemble runs arriving later than expected. For the magnetic field, we are more or less able to reproduce the positive constant B_T component although its strength is underestimated in both cases, the bipolar B_N component is overestimated, and the radial B_R is far off from the actual measurement at SolO and more or less accurate for Wind.

Figure 6 shows the 3DCORE shapes for representative ensemble members, using the same model parameters as the representative fits in Fig. 5, from a side-top and a frontal viewpoint. The color codes for the plotted inner planets and spacecraft are the same as in Fig. 2. In both cases, we show three different fits, for the April event, these are the individual SolO and Wind fits and the combined triple-point fit. In the case of the May event, we show the individual SolO and Wind fit including the combined SolO + Wind fit.

Fig. 6. Representative 3DCORE geometries for multiple different in situ combinations, including the inner planets and satellites using the same color codes as in Fig. 2. The initial parameters for the representative ensemble members are the same as for the representative in situ signatures in Fig. 5. The two top panels show a view from the top while the bottom panels show a view from the front.

Table 2 shows the inferred model parameters for the April and May event using all seven different spacecraft combinations. The final row in each case shows the final achieved error threshold, e^k, for each observation. From these tabular results, we can see that the general solutions for the geometry of the April event do not change much when using different observations. A larger difference can be found in the parameters that describe the CME kinematics with the single-point observations showing a large range of possible values for v₀, v_sw, and Γ. While the Γ parameter seems to be largely unconstrained in the multi-point results, the variations on the initial CME velocity and solar wind background speed are reduced significantly, demonstrating that these multi-point observations combined with our rigid CME geometry set stronger constraints for the kinematics.

The ensemble solutions for the May event show larger variations, specifically concerning the orientation and propagation direction of the flux rope. Specifically there is a larger difference in the longitudinal propagation parameter that spans a range of over 35 deg. Due to these larger variations, it is harder to find general solutions that can reproduce the observed measurements for two or three spacecraft, as can be seen in the large error values ϵ_final for the dual-point fits. These issues clearly show that our assumed geometry and model has significant issues attempting to describe the May event.

Figure 7 shows a scatter plot for the twist and aspect-ratio solutions from the overall ensembles of the April event. In this plot, we see that the δ − τ degeneracy persist despite the high bias towards higher δ values. As previously explained in Sect. 4.3, this can be explained by way of a purely geometrical effect. The solid lines show the τ ∝ δ/E(δ) trend for each single case. Unfortunately this means that one of the key issues in W21 remains unresolved. We are not able to solve the δ − τ degeneracy using dual or triple point observations. In taking a closer look at the spacecraft positions in the April event in Fig. 6, we are able to see that all three spacecraft are at a similar distance below the inferred flux rope axis. This further reduces the amount of information that can be extracted on the aspect ratio within our magnetic field model.

Fig. 7. Scatter plot showing the τ-δ degeneracy for the April event with the SolO, Bepi, and Solo+Bepi+Wind ensembles. The solid lines show the τ ∝ δ/E(δ) trend.

Lastly we can also compare some of the results to the remote observations shown in Sect. 2. We find that in general the inferred directions match the given values from the HELCATS catalog very well. For the April event we find a HEEQ longitude and latitude of 15 to 28deg and −3 to −10deg where the estimated value is −2/−6 from the images. For the May event we infer a HEEQ longitude and latitude of −23 to 20deg and 1 to 18deg, where the estimated value is +36/+11.

6. Conclusion

As a summary, we have shown in this paper that it is possible to reconstruct multiple in situ flux rope measurements simultaneously using an empirical CME MFR model. We chose two different scenarios to test our approach, with the first (the April event) consisting of three separate nearby measurements and the second (the May event) consisting of three measurements at larger longitudinal separations.

In the first case, with the nearby measurements, we do not find any major issues in our attempt to reconstruct the magnetic field measurements using our model. The seven different combinations for reconstruction show more or less the same picture, further reinforcing the validity of our description on smaller scales.

The second case illustrates the limits of our modeling approach as we are no longer able to consistently describe all three observations simultaneously. This can be attributed to large longitudinal separations of the spacecraft, which allows for a vastly different solar wind environment at the different locations. The two measurements for Wind and Bepi are strongly perturbed, which may be a hint that they only sample the CME structure at the very edge or have undergone interaction with the solar wind or heliospheric current sheet (Winslow et al. 2016). Such interactions may impose a limit to which multi-point magnetic field observations can be used for reconstructing the CME structure depending on the specific configuration of the environment.

The problem of estimating the cross-section aspect ratio and the twist has not been solved by using multi-point observations. Due to the constellation of the April event, where the model depicts all the spacecraft to lie below the flux rope axis at similar distances, it is still not entirely clear if this problem persists for other multi-point events or whether it is inherent. The specific magnetic field model that is used in our analysis may prove problematic in this aspect. Unless we are given clear observations of a crossing both below and above the flux rope axis, namely, a positive and negative unipolar B_r signature, we believe that we would not be able to resolve the degeneracy of the cross-section aspect ratio and the twist number when using our model. Further studies using different magnetic field models, specifically models that are more physical and include a proper description elliptical cross-sections (e.g., Nieves-Chinchilla et al. 2018), may resolve this problem or at least reduce the uncertainties.

In general, we believe that our approach works well as a first order approximation for multi-point observations in close proximity. Due to the different solar wind environments that can be present over longitudinal separations larger than ca. 10 deg, our model will be limited due to its reduced capability in describing the arrival times and flux rope orientation at varying positions. The next obvious step would be the construction of more complex and computationally expensive models that include deformations attributed to the solar wind environment. This will invariably require a proper 3D description of the solar wind within at least 1 AU as well as a CME model with a deformable CME structure. As our approach is valid, in principal, for any forward-based simulation model, the entire fitting pipeline can be re-used for this newer class of models. Our current model can also be used as a first step in feeding a more complicated model with initial conditions in order to limit the volume of the parameter space.

The April event represents one of the most optimal multi-point observations recorded thus far. New missions have been successfully launched within the last years, including BepiColombo, Solar Orbiter, and Parker Solar Probe. Expected launches of other missions, such as JUICE combined with the rising solar activity of cycle 25, should allow us to observe many additional multi-point events within the next decade. These observations will be important for testing our existing approach in future studies to develop models that can better describe the global ICME structure.

Acknowledgments

A.J.W, C.M., R.L.B, M.A.R., T.A., M.B. and J.H. thank the Austrian Science Fund (FWF): P31521-N27, P31659-N27, P31265-N27. We have benefited from the availability of Solar Orbiter, STEREO, Wind, and BepiColombo data, and thus would like to thank the instrument teams, and the SPDF CDAWeb and Solar Orbiter (http://soar.esac.esa.int/soar/) data archives for their distribution of data. The Solar Orbiter magnetometer was funded by the UK Space Agency (grant ST/T001062/1). This research was supported by funding from the Science and Technology Facilities Council (STFC) studentship ST/N504336/1 (E.D.) and STFC grant ST/S000361/1 (T.H.). D.H., I.R., H.-U.A. were supported by the German Ministerium für Wirtschaft und Energie and the German Zentrum für Luft- und Raumfahrt under contract 50 QW 1501. 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.

References

All Tables

All Figures

Fig. 1. In situ magnetic field measurements in local radial tangential normal (RTN) coordinates from Solar Orbiter (a/e), BepiColombo (b/f), and Wind (c/g) for the 2020 April 20 (left) and 2020 May 28 (right) event. The two bottom panels (d/h) show the in situ proton density number Np and the associated plasma beta coefficient, βp. The blue shaded area represents our estimated intervals for the unperturbed magnetic flux ropes (UMFR) which we subsequently use in our fitting approach.

In the text

Fig. 2. Spacecraft positions for Solar Orbiter, BepiColombo, Parker Solar Probe, and STEREO-A on 2020 April 20 00:00 UTC (left) and on 2020 May 29 00:00 UTC (right) in HEEQ coordinates. Furthermore shown are the planets Mercury (brown), Venus (orange), and Earth (green). Colored trails show orbital movements of the spacecraft within a time range of two months around the specifically chosen position. Two open-ended cones show the field of view (FOV) of the HI-1 (20°) and HI-2 (70°) imaging instruments aboard STEREO-A.

In the text

Fig. 3. White-light images, with the stars subtracted using a median filter, from the STEREO-A HI-1 instrument showing the April event on the 16 of April 07:29 UTC (left) and the May event on the 26 of May 15:29 UTC (right). The image of the April event shows a slightly flattened Earth-directed CME with a faint leading edge, an elliptical cavity, and a strong density pile up at the back. The image of the May event shows a structure at a higher latitude that is hard to identify in a single image (marked with the yellow arrow) and does not contain an obvious cavity.

In the text

	Fig. 4. Example of a 3DCORE simulation showing the torus-like structure at near 1 AU with a highly elliptical cross-section and two magnetic field lines at differing radial distances (twist number τ = 5).
In the text

Fig. 5. Reconstructed synthetic magnetic field measurements (dashed) with 1σ confidence interval (shaded area) using the triple point constellation for the April event and a dual-point constellation for the May event. The reconstructed magnetic field is overlayed on top of the real in situ measurements for Solar Orbiter (a/d), BepiColombo (b/e), and Wind (c/f). The blue shaded area represents the defined UMFR interval and the dotted vertical lines show the used fitting points ti. All magnetic fields are shown in local RTN coordinates.

In the text

	Fig. 6. Representative 3DCORE geometries for multiple different in situ combinations, including the inner planets and satellites using the same color codes as in Fig. 2. The initial parameters for the representative ensemble members are the same as for the representative in situ signatures in Fig. 5. The two top panels show a view from the top while the bottom panels show a view from the front.
In the text

	Fig. 7. Scatter plot showing the τ-δ degeneracy for the April event with the SolO, Bepi, and Solo+Bepi+Wind ensembles. The solid lines show the τ ∝ δ/E(δ) trend.
In the text