Multi point analysis of coronal mass ejection flux ropes using combined data from Solar Orbiter, BepiColombo and Wind

The recent launch of Solar Orbiter and BepiColombo opened a brief window in which these two spacecraft were positioned in a constellation that allows for the detailed sampling of any Earth-directed CMEs. Fortunately, two such events occurred with in situ detections of an ICME by Solar Orbiter on the 19th of April and the 28th of May 2020. These two events were subsequently also observed in situ by BepiColombo and Wind around a day later. We attempt to reconstruct the observed in situ magnetic field measurements for all three spacecraft simultaneously using an empirical magnetic flux rope model. This allows us to test the validity of our flux rope model on a larger and more global scale and allows for cross-validation of the analysis with different spacecraft combinations. Finally, we can also compare the results from the in situ modeling to remote observations obtained from the STEREO-A heliospheric imagers. We make use of the 3D coronal rope ejection model in order to simulate the ICME evolution. We adapt a previously developed ABC-SMC fitting algorithm for the application to multi point scenarios. We show that we are able to generally reconstruct the flux ropes signatures at three different spacecraft positions simultaneously using our model in combination with the flux rope fitting algorithm. For the well-behaved 19th of April ICME our approach works very well. The 28th of May ICME, on the other hand, shows the limitations of our approach. Unfortunately, the usage of multi-point observations for these events does not appear to solve inherent issues, such as the estimation of the magnetic field twist or flux rope aspect-ratios due to the specific constellation of the spacecraft positions. As our general approach can be used for any fast forward simulation-based model we give a blueprint for future studies using more advanced ICME models.


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, a 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 a strongly south orientated magnetic field over the duration of multiple hours. As such, the understanding of the structure of these magnetic obstacles is key to link-Article number, page 1 of 11 arXiv:2103.16187v2 [astro-ph.SR] 17 May 2021 ing 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 there has been a focus on the study of multi point events, i.e. 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 on board 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 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), that aim to mimic the general properties that ICMEs are thought to have. These models serve as a bridge in 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 semiempirical 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 de-scribes 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 will show if it is even possible to describe multiple observations using the global assumptions within our model and, if not, give a hint on 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 which can further the understanding on the underlying event and the used model.
Section 2 introduces the two events for which we will test our approach. The first event, the April 19 ICME, represents three very close by 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 Section 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 Section 6 we discuss the results and the conclusions we can draw.

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 will make use of is the 2020 April 19 ICME which we henceforth refer to as 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 Bepi-Colombo on the April 20 at 01:34 UTC and 03:09 UTC respectively. For an extensive description of this event we refer the reader to Davies et al. (2021) as we will only focus on the details that are of interest for our analysis.
The second event that we have chosen is the 2020 May 28 ICME which we henceforth refer to as the May event. Similar 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 in- In situ magnetic field measurements in local radial tangential normal (RTN) coordinates from Solar Orbiter (1a/e), BepiColombo (1b/f) and Wind (1c/g) for the 2020 April 20 (left) and 2020 May 28 (right) event. The two bottom panels (1d/h) show the in situ proton density number N p and the associated plasma beta coefficient β p . The blue shaded area represents our estimated intervals for the unperturbed magnetic flux ropes (UMFR) which we will be subsequently using for our fitting approach.
strument (Ogilvie et al. 1995) for both events. The magnetic field measurements are shown in their respective local radialtangential-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 will be the intervals that 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. 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 Section 2.1. Unfortunately, PSP was located on the back 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. 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, the transversal component B T unipolar and the radial component B R negative near zero. According to the standard classification schemes this corresponds to a low inclination lefthanded South-East-North flux rope (Bothmer & Schwenn 1998;Davies et al. 2021). Due to the textbook like structure of this flux rope it represents a unique opportunity to test the validity of our model with multiple close-by 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 one 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 will partially alleviate these issues by applying 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 hours 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 edgeon 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 Article number, page 3 of 11 A&A proofs: manuscript no. main  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 few hours. Both the April 1 and May 2 event are part of the HELCATS catalog. Despite only being observed from a single vantage point one can extract rough estimates for the propagation speed and direction by tracking the CME evolution. Under the assumption of

Model
In this paper we use an updated version of the 3DCORE model to analyze our two selected events. While we will give a brief general overview we will focus on the changes to the implementation that we have made in order to facilitate analyzing multi point events. The two primary changes that are worth mentioning are the "twist reparametrzitation" where we substituted the magnetic field twist parameter using a different related parameter which has beneficial numerical properties and the introduction of a new noise model. For the interested reader we refer to W21 for a detailed description of the original implementation.

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, that is described in Riley & Crooker (2004), is approximated by stretching the cross-section of the torus into an elliptical shape. 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 in 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 the usage Monte-Carlo based methods which require in the order of thousands or millions of simulation runs.
In total the model contains 12 controllable parameters: the initialization time t 0 , propagation direction and orientation lat, lon, inc, cross-section width d 1au , cross-section aspect ratio δ, initialization distance r 0 , initial speed v 0 , magnetic field strength b 1au 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 1au r(t) 1.14 and b(t) = b 1au r(t) −1.68 , where r(t) is the distance of the CME front from the Sun.

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 exists as we are not directly capable of measuring the twist of the magnetic field, and instead 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 nonlinear 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 postprocessing step.

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 which is easy to implement but fails to fully capture the small scale behaviour of the magnetic field that is seen in the observations. We will therefore 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 the type of fluctuations within our model 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) ∼ N(0, 1) and its corresponding Fourier transformμ(k) = F (µ(t)) one can generate a random time series: Article number, page 5 of 11 A&A proofs: manuscript no. main so thatδ(t) ≈ 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 that is 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 little 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 will act 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 δ ∼ N(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 of 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 as the the point-to-point fluctuations are more accurately modelled to decrease 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 are used in our noise model. But as our summary statistic that 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.

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 which 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 one 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 Section 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 if 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.

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: 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 negative trade-off of another and the end result clearly shows if a certain observation cannot be properly fitted.
As a consequence, a simulation output for N observations x = (x 1 , . . . , 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 k 0 = ρ(0, y k ) and e k j+1 is set as a quantile fraction of e k j .

Adaptive Time Markers
An important aspect of the fitting algorithm is to determine if a simulation run can be considered as valid in the first place. In a large ensemble simulation a majority of runs will not generate any valid results due to missing one or many of the defined observers or producing 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 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 1 and t e just after t K one is able to closely control 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 which produce signatures of just the right length. As one adds additional constraints by adding observations the efficiency of the fitting algorithm can suffer significantly. For this reason we implemented an adaptive time marker scheme, in which the the t s and t e markers are set at a large temporal distance to the 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 considerably helps the sampling process 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.

Bias
An inherent property of our 3DCORE model is that it is an empirical model and does not contain any physical calculations. As we will 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 their 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 not hitting any of the predefined observers or runs that generate magnetic fields that are too strongly time shifted are automatically rejected before computing any error metrics.
While the size of the shapes is limited by using bounds for the model parameters this nonetheless leads to an extremely strong bias towards the largest allowed shapes as they do not fall victim to this selection bias. The two parameters in our model that control the size of our structure are the d 1AU (torus width) and δ (torus height) parameters. The width parameter d 1AU 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 mention previously we know that we cannot constrain the δ parameter using a single in situ measurement. This generally leads to linear bias p(δ) ∝ δ as a structure that has a larger aspect ratio δ has a larger chance of hitting the observer.
Under certain circumstances this bias does not occur, such as a CME that directly propagates towards the observer, which makes this issue hard to handle. In general this must be done at a case by case basis by looking at the model simulations that are accepted and rejected due to hitting or missing the observers. The conclusion is that our analysis will normally generate extremely large estimates for the δ parameter while the correct interpretation is that nothing about this parameter is actually known. One hope is that this bias would disappear when adding multiple observations which would largely resolve this issue.

Results
In this section we present the results of our paper in which 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 this then leads to a total of seven different possible data combinations which can be used. This allows us to study in detail how the results using our method differs when using either the full data set or smaller sub sets. 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.
The model initialization time t 0 is set to a specific date for both events, the initialization distance r 0 is fixed to 35 solar radii, and the initial speed v 0 and background solar wind speed v s w are bounded within a smaller range. All other parameters are largely unbounded within a sensible range. For the April event we set t 0 = 2020-04-15T23:00 UTC, v 0 ∈ [350,750] and v sw ∈ [275,375]. For the May event we use t 0 = 2020-05-26T14:00 UTC, v 0 ∈ [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 largely 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 the Figure 5 we used an extremely low number of points for these two measurements which are located towards the endpoints of the flux rope signature.
The results for the April event show 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 are 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 can be seen alone from the spread. Since we only took the SolO and Wind measurements into account 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 Figure  the same as in Figure 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 for the May event we show the individual SolO and Wind fit including the combined SolO + Wind fit. 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 0 , 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 showing 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 Section 4.3 this can be explained due to 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. Taking a closer look at the spacecraft positions in the April event in Figure 6 we are able to see that all three spacecraft are at a similar distance below the inferred flux rope axis. This furthermore reduces the amount of information that can be extracted on the aspect ratio within our magnetic field model.
Lastly we can also compare some of the results to the remote observations shown in Section 2.1. 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 28 deg and −3 to −10 deg where the estimated value is −2/ − 6 from the images. For the May event we infer a HEEQ longitude and latitude of −23 to 20 deg and 1 to 18 deg where the estimated value is +36/ + 11.

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 close by measurements and the second (the May event) consisting of three measurements at larger longitudinal separations.
In the first case, with the close by 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 will persist for other multi point events or 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, i.e. a positive and negative unipolar B r signature we believe that we will 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 crosssections (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 not properly being able to describe 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 due to the solar wind environment. This will invariably require a proper 3D description of the solar wind within at least 1 AU and a CME model with a deformable CME structure. As our approach is, in principal, valid 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 best multi point observations so far recorded. 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 and future studies in order to develop models that can better describe the global ICME structure.