Eruption and propagation of twisted flux ropes from the base of the solar corona to 1 au

Interplanetary Coronal Mass Ejections (ICMEs) originate from the eruption of complex magnetic structures occurring in our star's atmosphere. Determining the general properties of ICMEs and the physical processes at the heart of their interactions with the solar wind is a hard task, in particular using only unidimensional in situ profiles. Thus, these phenomena are still not well understood. In this study we simulate the propagation of a set of flux ropes in order to understand some of the physical processes occurring during the propagation of an ICME such as their growth or their rotation. We present simulations of the propagation of a set of flux ropes in a simplified solar wind. We consider different magnetic field strengths and sizes at the initiation of the eruption, and characterize their influence on the properties of the flux ropes during their propagation. We use the 3D MHD module of the PLUTO code on an Adaptive Mesh Refinement grid. The evolution of the magnetic field of the flux rope during the propagation matches evolution law deduced from in situ observations. We also simulate in situ profiles that spacecraft would have measured at the Earth, and we compare with the results of statistical studies. We find a good match between simulated in situ profiles and typical profiles obtained in these studies. During their propagation, flux ropes interact with the magnetic field of the wind but still show realistic signatures of ICMEs when analyzed with synthetic satellite crossings. We also show that flux ropes with different shapes and orientations can lead to similar unidimensional crossings. This warrants some care when extracting magnetic topology of ICMEs using unidimensional crossings.


Introduction
An interplanetary coronal mass ejection (ICME) corresponds to the magnetically dominated plasma that is the interplanetary counterpart of eruptions occurring in the solar corona.When these magnetic structures, called magnetic ejecta (ME), are expelled from the solar corona, they propagate in the interplanetary medium.During the propagation, ICMEs and their different substructures can be probed by in situ instruments on board spacecraft such as ACE (Chiu et al. 1998) and Wind (Harten & Clark 1995), situated at the first Lagrange point (L1); the Parker Solar Probe (Fox et al. 2016) and Solar Orbiter (Müller 2020), which span a wide range of heliospheric distances; or STEREO A and B (Kaiser & Adams 2007), which sample other heliographic longitudes than the observatories at L1.
If ICMEs are fast enough (i.e., faster than the local solar wind), they may accumulate solar wind material ahead of them and form a sheath.In that case the solar wind encountered during the propagation of the ICME does not have enough time to flow around the ICME.It thus creates a region of heated and compressed solar material.Moreover, the expansion of the ME during its propagation enhances the accumulation of the solar wind material (Siscoe & Odstrcil 2008).Recent studies (Moissard et al. 2019;Kilpua et al. 2020) have shown that this region is more turbulent than the local solar wind by computing the power spectrum of the magnetic field fluctuations and deducing its anisotropy in the sheath and in the solar wind.They found that the orientation of the magnetic field fluctuations and their magnitudes are more erratic in the sheath than in the solar wind.However, ICMEs do not always have an observable sheath; in fact, no-sheath ICMEs are comparable to the slow ICMEs with a sheath (Regnault et al. 2020).
If the speed of the ICME is faster than the local fast mode speed, then they can develop a fast-forward shock ahead of the sheath.Not all ICMEs have a shock in front of their sheaths.Salman et al. (2020) studied the difference between ICME sheaths that have a front shock with those that do not.They found that the magnetic and plasma properties of sheaths with a front shock are higher in magnitude than for those without it.Moreover, all shocks detected in the in situ data are not necessarily ICME-driven, but they can be related to a stream interacting region (SIR; Kilpua et al. 2015), which is formed by the interaction of a fast solar wind stream with the slow solar wind stream.
Since it is a magnetically dominated plasma, the ME is characterized by a low plasma β (thermal pressure over magnetic pressure), typically ≈0.1.We also observe a decrease in the proton temperature and density compared to the sheath area, due to the expansion of the ME along its propagation.This expansion is mainly due to the decrease in the solar wind pressure with distance (Démoulin 2009).The radial expansion of the ME produces a monotonic decrease in the speed profile (Klein & Burlaga 1982;Farrugia et al. 1993), although this decrease is not always present (e.g., in the case of a fast stream or another ICME overtaking the primary structure; see more details in Gulisano et al. 2010).The ratio of the fluctuations in the magnetic field to the magnitude is lower (Masías-Meza et al. 2016) because MEs are more coherent magnetic structures than the solar wind, so the magnetic field is smoother.
We sometimes observe a smooth rotation of one of the components of the magnetic field inside the ME.If this is observed along with a low proton temperature and a high magnetic field intensity, then the ME is called a magnetic cloud (MC, Burlaga et al. 1981).This smooth rotation of one component of the magnetic field is interpreted as a trace of the flux rope (FR) structure of the MC.Such a shape is often used to describe the magnetic structure of ICMEs.For example, the Lundquist (1951) model assumes a FR structure with a circular cross section and an invariant axis.
One-third of MEs are detected as MCs (Richardson & Cane 2004;Wu & Lepping 2011), although this low proportion may be directly linked to the crossing of the structure.Kilpua et al. (2011) have shown that the observation of the FR signature is highly impacted by the trajectory of the spacecraft through the ICME.Using multi-spacecraft observations of ICMEs, Jian et al. (2006), Démoulin et al. (2013), and Zhang et al. (2013) also found that the trajectory of the spacecraft through the ICME is important for the observation of MC signatures.
We see here one limitation of in situ data analysis for the study of the property of the propagating ICME.The 1D crossing of the 3D complex structure of the ICME only provides limited information (Riley et al. 2004;Al-Haddad et al. 2019).Some reconstruction techniques are able to get the global properties of the ICME from the 1D crossings, but it comes with the price of significant approximation of the magnetic structure (e.g., force-freeness, cylindrical symmetry) following for instance the Lundquist (1951) or Gold & Hoyle (1960) models.According to Riley & Crooker (2004), a force-free flux rope should not keep a circular cross section during its propagation, which challenges then the applicability of the cylindrical symmetry assumption.
Moreover, studying the evolution of the properties of ICMEs is a hard task due to the low number of spacecraft in the interplanetary medium.Ideally, one would need to be able to measure the properties of the same ICME at different distances from the Sun.These ICMEs are rare and are mostly studied in single case studies (see, e.g., Kilpua et al. 2011 for a review of multi-spacecraft observations with the STEREO mission).Janvier et al. (2019) managed to perform a statistical study of ICMEs detected close to Mercury with the Messenger probe, close to Venus with the Venus Express probe, and close to the Earth with the ACE probe.However, Messenger and Venus Express do not have any in situ measurements of the properties of the interplanetary plasma (e.g., temperature, speed, and proton density), which limits the conclusions drawn by the paper to the magnetic properties of ICMEs.
Simulations provide an interesting complementary approach to the study of the propagation of a flux rope in the interplanetary medium.In a simulation, we can perform a 3D analysis of the ICME structure.This allows a more detailed analysis of the structure of the ICME, but also of the physical processes happening during the propagation of a flux rope magnetic structure in a simulated solar wind.We can then compare the properties of the simulated ICMEs with properties of actual ICMEs, for example by doing synthetic crossings through the ICME and comparing them with observations.Simulations of ICMEs have already been performed.However, most of them focus either on the trigger part of the ICME (see, e.g., Aulanier 2010;Kliem et al. 2012) or on the propagation of the ICME up to the Earth.The EUHFORIA simulations insert a spheromak magnetic structure at 0.1 au and study its propagation up to the Earth (Verbeke et al. 2019).Shen et al. (2014) performed a Sun-to-Earth simulation of the propagation of a magnetic structure that is initiated closer to the Sun but not anchored to it.Inserting the magnetic structure at a certain distance from the Sun cannot properly take into account the different physical mechanisms linked to the eruption of the flux rope (e.g., its impact on the rotation).Török et al. (2018) took into account the effect of the eruption mechanism on the properties of the propagating flux rope when performing a Sun-to-Earth simulation of an eruptive prominence event.They did a relaxation of a flux rope structure (anchored at the solar surface), triggered its eruption with converging flows in the photosphere, and then studied its propagation up to the Earth.
Similarly, the aim of this paper is mainly to investigate the propagation of a flux rope that erupted within the same simulation domain.We run a parametric study in order to highlight the physical mechanisms that happen during the propagation of a flux rope in the solar wind.
In Sect. 2 we describe the numerical setup and the flux rope model used in this study.In Sect. 3 we investigate the evolution of the magnetic structure during its propagation.In Sect. 4 we compare the results of the simulations with results obtained with in situ measurements.The results are discussed in Sect. 5.

The PLUTO code
In this study we used the PLUTO code to solve the 3D magnetohydrodynamics (MHD) equations in an adaptive mesh refinement (AMR) spherical grid (Mignone et al. 2012).PLUTO is a multi-physics and multi-solver code using a finite-volume Godunov-type scheme.Such a scheme requires computing the flux through the interfaces of each elementary volume paving the simulation domain.To do so, we used the Harteen, Lax, van Leer solver (Toro 2009).
The MHD equations are written in their conservative form: where B is the magnetic field intensity, ρ the mass density, u the speed, g the gravitational acceleration, and where p tot = p th +B 2 /2 corresponds to the total pressure (thermal + magnetic).
Finally, E is the total energy density which is given by To close the system of equations, we used an ideal equation of state for : ).The orange areas correspond to level 4 of the AMR grid, purple to level 3, and red to level 2. Level = 1 does not appear here.
The ratio of specific heats, γ, is set to 1.05 to maintain a heated corona.This value is chosen to sufficiently heat the low corona so that it can drive a solar wind (see Sect. 2.2 for more details).
Equations are solved in a spherical geometry (r,θ,ϕ) in a grid that ranges from 1 to 420 R for r, from 0 to π for θ, and from 0 to 2π for ϕ.The initial grid is stretched in the r direction so that the spatial resolution decreases while going away from the Sun following a logarithm.Other simulations (Török et al. 2018;Poedts et al. 2020) use a non-uniform grid in r with a higher spatial resolution closer to the Sun.Moreover, the AMR grid adapts its refinement in function of the gradient of the total energy density.This allows us to keep a finer grid in the region of the propagating ICME while having a coarser grid far from it.We used four levels of AMR with a refinement ratio set to 2. At maximum spatial resolution (level 4), the step in θ and ϕ is 0.012 rad and ≈0.012 R in the r direction close to the Sun and ≈2.5 R at 210 R .We did not prescribe any resistivity in this model.However, these ideal MHD equations can still lead to a change in connectivity of magnetic field lines due to numerical resistivity.
Figure 1 shows two equatorial slices of the speed magnitude at different times during the propagation of a simulated ICME.Boxes with different colors show the level of refinement in the equatorial slice.We see here that the high resolution (in orange) closely follows the ICME, shown in yellow, during its propagation in the heliosphere.

Solar wind model
We used the solar wind model of Réville et al. (2015), initialized using the Parker model (Parker 1958) with a polytropic equation of state with P the pressure, ρ the density, and γ = 1.05 the ratio of specific heat.This assumption produces a polytropic heating of the full corona and makes it almost isothermal at ≈10 6 K.The surface magnetic field is set up to be a dipole with a maximum field strength at the surface of 2.9 G.This configuration is close to the magnetic configuration of the Sun during the solar minimum (Réville & Brun 2017).The polytropic assumption is used only during the initialization, the model then uses an ideal equation of state (see Sect. 2.1).This model drives a solar wind that has the same speed in all direction; we thus calibrate the model so that we have a speed of ≈ 440 km s −1 at 1 au.We also paid attention to the mass loss rate and the angular momentum loss of the simulated solar wind.The model is calibrated so that the mass loss rate and the angular momentum (averaged from 3 to 250 R ) correspond to the values measured in the solar system (see, e.g., Finley et al. 2018 for the angular momentum and Schwadron & McComas 2008 for the mass loss rate).For more details about the solar wind model, see Réville et al. (2015).
Once the simulation evolves in time, the hydrodynamic properties (speed, density) and the magnetic properties will interact with each other and produce a self-consistent MHD solar wind.

Flux rope model
In this study we used the analytical formulation of a flux rope described in Titov et al. (2014).The modified Titov-Démoulin (hereafter TDm) flux rope provides an analytical formulation of a ring current that has a magnetic flux rope structure.The model is built so that the magnetic structure is force-free if the ambient magnetic field is perpendicular to the torus plane and distributed uniformly along its toroidal segment.We call this magnetic field B ⊥ .A horizontal magnetic field can meet these requirements, but a bipolar magnetic configuration allows us to have a magnetic field that is approximately uniform along the toroidal segment.
Figure 2 shows a schematic of the TDm flux rope with the different parameters for the TDm model: R and a are the major and minor radius of the torus, respectively; d is the depth at which the torus is buried below the solar surface, represented by the (x,y) plane in this figure.(Titov et al. 2014).The magnetic field of the twisted flux rope appears in green.B ⊥ correspond to the ambient magnetic field.
In order to be force-free, the intensity of the ring current should be equal to the Shafranov current I s (Shafranov 1966) which is given by with µ the magnetic permeability and l i is the internal selfinductance per unit length of the rope that is ≈1 in this case.
In this work, we study the effect of the initial magnetic field on the properties of the propagating flux rope.To do so we define the ζ parameter such that with I the intensity of the ring current of the flux rope.
The ζ parameter allows us to control the initial magnetic field of the flux rope.It thus means that a flux rope initialized with ζ 1 is not in a force-free state.The TDm flux ropes in this study are initialized out of equilibrium.This has the advantage of allowing us to skip the costly relaxation step and to immediately simulate the eruptive phase.A more thorough treatment of the initial relaxation phase is deferred to future works.
The original TDm model is proposed by Titov et al. (2014) in two flavors.The difference between the two cases lies in the current density distribution inside the torus.Case 1 has a "hollowcore" current distribution, meaning that the twist of the magnetic field, thus the current, is concentrated at the outer edge of the torus.This current distribution has been observed in MCs by Lanabere et al. (2020).They studied the distribution of twist in FR observed at 1 au (assuming a Lundquist model) and they found a constant twist inside the core of the FR and a twist a higher by a factor 2 at the boundary.Case 2 has a parabolic distribution of the current, meaning that the magnetic twist is distributed over the whole cross section of the flux rope.In this study we only use the case 1 flux ropes; we will leave the propagation of the case 2 flux ropes for another study.The code developed to insert a TDm flux rope in the PLUTO data files has been released and is available online1 .Figure 3 shows an example of a TDm initialization in the global wind model.On the right side of the figure, a spherical slice shows B r (blue corresponds to negative B r ).The colored lines correspond to magnetic field lines chosen in the negative polarity, and anchored in a 0.9a width circle at the center of the flux rope.The shape of the ring current that appears in yellow and white is typical of case 1, as mentioned earlier.On the left side, a zoomed-out view of the global solar wind magnetic field is shown with white streamlines.One can see here that the TDm flux rope is initialized within the equatorial streamer region where the density is higher and the magnetic field is in a closed configuration.

Boundary conditions
The simulation domain is surrounded by two layers of cells (called ghost cells), whose values set the behavior at the boundary.The outer boundary conditions (r = 420 R ) are the Neumann boundary conditions, meaning a constant gradient.At the poles (θ = 0 and π) the conditions are axisymmetric, which means that the components parallel to the θ direction change sign and the others remain unchanged.At the ϕ = 0 and 2π boundaries, the condition is also of Neumann type which is enough for the asymmetric solar wind.The ICMEs simulated in this study propagate in the direction radially opposite to this boundary to avoid any spurious effect due to its presence.The inner boundary condition (r = R ) is set differently close to the flux rope and far away from it.We defined the ratio with P SW mag and P TDm mag the magnetic pressure of the solar wind (before the insertion of the TDm) and of the TDm flux rope only, respectively.
Cells with χ BC < 0.5 are far from the flux rope, and the boundary conditions are those of the classic solar wind model.In that case the poloidal (r and θ) speed is set to be parallel to the poloidal magnetic field to limit the production of current and to ensure that the Sun is a perfect conductor.
Cells with χ BC > 0.5 are close to the flux rope and the boundary conditions are set such that there is a fixed gradient for all the components of the magnetic field to extend the magnetic field lines in the two layers of ghosts cells.We applied a reflective condition on the r component on the speed, simulating a very dense solar surface.The θ component of the velocity had a Neumann-type boundary condition at this boundary.
The thick flux ropes (see Sect. 2.5) had a Neumann-type inner boundary condition for the pressure, while the thin ones had a condition on the pressure that ensures hydrostatic equilibrium in the ghost cells.In the end, we find that this condition has a low impact on the property of the flux rope.

Flux rope properties for the parametric study
To study the effects of the flux rope properties on its propagation, we investigate in the following the influence of the initial magnetic field strength and thickness of the flux rope.We initialized six different TDm flux ropes as follows: we chose two sets of TDm flux ropes, with one set corresponding to a thickness a of 0.1 R and the other to a thickness of 0.05 R .
Each of these sets then has three different initial magnetic fields, for which the parameters are summarized in Table 1.All are initialized at the same location on the Sun's surface, at (θ 0 ,ϕ 0 ) = ( π 2 ,π) aligned along the solar equator.The ∆ parameter (last column of Table 1) controls the width of the ring current (see the yellow ring in the left panel of Fig. 3); it is set to be a 10 so that the condition δ 1 is met (assumption used for the analytic development of the TDm; see Titov et al. 2014).All six flux ropes have a major radius R = 0.3 R and have d = 0.15 R .
The different initial magnetic fields range from 16 G to 84 G, which typically corresponds to the magnetic field observed in prominences (Casini et al. 2003).The magnetic flux of the flux ropes ranges from 1.3 to 4.4 × 10 21 Mx, which is close to the flux of large sunspots (Zwaan 1987).Therefore, the six flux ropes have realistic (i.e., near solar) values of magnetic field intensity and flux, although it can be argued that the largest prominences are generally quiescent sun regions rather than large sunspots regions.

Eruption
In this study the flux ropes are initialized in an unstable state and erupt quickly after their insertion.Figure 4 shows 3D visualizations of the TDm flux rope 1 min and 8 min (physical time) after the insertion of the flux rope structure for the thick2 case.We recall that the speed gained by the FR comes from Lorentz forces (no ad hoc speed is added to the simulation) and that the initial FR has magnetic field and flux similar to that of the magnetic  structure found in the solar atmosphere.However, it evolves in the medium that did not have time to relax due to the insertion of the FR which is unstable from the beginning.The colored streamlines correspond to the magnetic field lines of the flux rope starting from a circle of radius 0.9a at the negative polarity where we find the footpoints of the flux rope.In the right panel, 8 min after the insertion, the thick2 FR has a speed of ≈900 km s −1 at 1.5 R .For comparison, the thin3 FR has a speed of ≈4400 km s −1 and reaches 1.5 R 1.5 min after its insertion.This case is more magnetized (see Table 1), and thus produces the largest initial Lorentz forces propelling the eruption.According to Chen (2011), plane-of-sky CME speeds can occasionally reach 3500 km s −1 .The thin3 case therefore corresponds to an extreme case in terms of speed during the eruption.
During the eruption some field lines change their connectivity, with some becoming the low-lying loops shown in the red rectangle.These loops are reminiscent of the post-flare loops, and are formed by reconnection happening during the eruption of the flux rope in the trailing region indicated with the orange rectangle.This is expected according to the standard model of solar eruptive flares in 2D (Carmichael 1964;Sturrock 1966;Hirayama 1974;Kopp & Pneuman 1976 type) or in 3D A14, page 5 of 13 (see Janvier 2017 and references therein).This kind of structure resembles the loops that are observed in ultraviolet images after an eruption (Aulanier et al. 2012 and references therein; e.g., Su et al. 2006;Warren et al. 2011).We note however that the reconnection happening in these simulations is numerically driven since we do not have any explicit resistivity in the simulation.
During the early propagation of the thick flux ropes, we observe two high speed flows starting from the inner boundary of the simulation.These high speed streams cover a small area in the wake of the ICME and reach typical speeds of the bulk of the ICME itself.Because of the high speed, a region of low pressure (10 orders of magnitude smaller than the solar wind pressure) is formed.A threshold beyond which the pressure value is set to 10 −12 (in PLUTO units) for the cells that encountered this situation.Interestingly, thin flux ropes have almost no high speed streams in their wake.

Propagation
In this section we discuss the 3D magnetic structure of TDm FRs and their evolution during their propagation of up to 210 R .Panel a) of Fig. 5 shows a radial cut along the θ = π 2 , ϕ = π line for the B ϕ in red and B θ in green.In panel b), the left side shows the magnetic structure of the thick1 FR and the right side shows the thin3 FR.From the top to the bottom in each panel, the flux ropes are at 15, 30, 50, and 210 R .The magnetic field lines start close to the front of the flux rope (purple) and go close to its end (green).The colored dots in panel a) are the seeds of the streamlines of the corresponding color in panel b).We can see some magnetic field lines that have reconnected with the open flux of the solar wind.
We observe that the ICMEs propagate radially.This is expected because the B of the solar wind has an axisymmetry and also because the simulated ICMEs do not encounter any other ICMEs.The encounter of one ICME with another can lead to their deflection as shown in some events (e.g., Lugaz et al. 2012), but it is not the case here.
In Fig. 5 the two sides (see white arrows in Fig. 5) of the thin3 case (right side of panel b)) are included in the equato-rial plane, while the right (left) side of thick1 is above (below) the equatorial plane.This suggests a rotation of the thick case compared to the thin case happening before 15 R .Nevertheless, we can see that the thick1 case seems to continue to rotate slowly after 15 R .A similar orientation is found for the FRs of the same set.It suggests that the initial magnetic field intensity (and thus the speed of the ICME) has no significant effect on the rotation of the FR during its propagation in our simulation.Lynch et al. (2009) studied the rotation of a flux rope formed by reconnection in numerical simulations, and they found a rotation for their left-handed flux rope (which goes in the same direction as the thick1 case) that reaches 50 • at 3.5 R , which is a similar angle to the one we observe, but at 15 R (see the left side of panel b) in Fig. 5).
The observed rotation of the TDm flux rope is also in agreement with Kay & Opher (2015) who showed, with CMEs propagating in an extrapolated magnetic field with the PFSS method (potential field source surface, Schatten et al. 1969;Altschuler & Newkirk 1969), that the rotation of the flux rope mainly occurs close to the Sun (<10 R ).When they do not consider the torque applied on CMEs during their propagation between 10 and 210 R , they only found a 10% error on the orientation of the CME at 210 R , implying that most of the rotation is due to the torque applied before 10 R .Moreover, Isavnin et al. (2014) studied the evolution of 14 CMEs in the heliosphere with data-driven simulations, and showed that ≈57% of the rotation happens below 30 R .We essentially find the same result in the thick1 case, which rotates quickly before 15 R and then rotates slowly up to 210 R .
The observed rotation of the FR could be caused by a kink instability, as discussed by Manchester et al. (2017).However, the thinnest FRs, the thin ones, do not seem to rotate as much as the thick and the kink instability is known to be predominant for thinner flux ropes (Titov et al. 2014).
On the other hand, Shiota et al. ( 2010) (along with Cohen et al. 2010) suggested that the reconnection with the ambient magnetic field (in particular inside the helmet streamer) is causing a rotation of the flux rope.This may be the process at heart here, since we can see in Fig. 4 that the flux rope partly reconnects with the magnetic field within the streamer.
A14, page 6 of 13 This could explain why we do not see a great deal of rotation happening after 15 R since the tip of the streamer is roughly at 10 R .
We see in the left panel of Fig. 5 that the B ϕ profile for both thick1 and thin3 FRs does not change its shape from 15 R to 210 R .Conversely, we observe that B θ goes progressively to lower values (compared to the minimum value of B ϕ ) during the propagation.This is clearly visible for the thin3 FR while this effect is less pronounced for the thick1 one.Such accumulation of B θ could be the trace of magnetic field lines at the rear of the ME that are reconnecting with each other.This reconnection process seems to be more predominant with the thin FRs than for the thick ones.
The magnetic structure of the FR from the same set (thin or thick) is very similar, but it seems that the two different sets show different magnetic structure.If we compare the magnetic field lines between the thick1 case and the thin3 FR in the right panel of Fig. 5, they do not seem very similar.However, if we focus on the part of the streamlines close to the (θ = π 2 , ϕ = π) line, we observe that the streamlines close to the leading edge of the FR (in purple) are roughly vertical, while those close to the center of the FR (in white, and light orange and green) are more horizontal.This matches the model described in Lundquist (1951) for a FR whose axis is aligned with the ecliptic plane.
This can be seen more easily at 210 R , but it is present at each distance presented in Fig. 5. Nevertheless, streamlines show significant differences between the two FRs when considering the sides of the FR (even if we rotate the FR).

Evolution of plasma and magnetic properties of the flux rope with distance
To study the flux rope properties, we need to identify the position of the flux rope.To do so, we developed an algorithm that fits a Gaussian on B ϕ close to the nose of the simulated ICME at each time step.According to the typical Lundquist (1951) model, we expect the toroidal component (B ϕ ) to be higher than the poloidal components (B r and B θ ) at the center of the flux rope.Even though the shape of B ϕ , is not necessarily expected to be Gaussian, it is still a good approximation of the location of the flux rope.To enhance the automatic detection, we do a two-step Gaussian fitting.After the first fit, we find the position of the local minimum of B ϕ around the position of the fitted Gaussian.We then use this position as an initial guess for the second Gaussian fit.Finally, the position of the FR corresponds to the position of the second fitted Gaussian.This fit is illustrated in Fig. 5 in the left panel, third row.In order to help the fitting algorithm find the position of the FR at the next time step, the speed at the center of the first A14, page 7 of 13 10 0 Heliospheric distances (au) Fig. 6.Evolution of the magnetic field (in nT) at the position determined by the Gaussian fitting of the ME as a function of the distance from the Sun (in au).The solid colored lines correspond to the simulated flux ropes.The black dashed line corresponds to a power law deduced from the observations in Winslow et al. (2015).Finally, the maximum magnetic field measured during an ICME encounter at Messenger (MES), Venus Express (VEX), Parker Solar Probe (PSP) and Solar Orbiter (SolO) (Winslow et al. 2015;Good & Forsyth 2016;Möstl et al. 2017Möstl et al. , 2020) ) are also shown in this plot.
fitted Gaussian is used to predict the next position of the ICME, and thus facilitates the automated detection.
Figure 6 shows the evolution of the magnetic field amplitude measured at the position determined by the Gaussian fitting for each time step as a function of the propagation distance in loglog scale.The values for the six simulated flux ropes are represented by the solid colored lines.We see that these intensities decrease as a power law along the flux rope propagation.
Similar laws have been derived to determine the magnetic field of the flux rope at different distances from the Sun, and are found in the literature (e.g., Bothmer & Schwenn 1997;Liu et al. 2005;Wang et al. 2005;Leitner et al. 2007;Winslow et al. 2015, and references therein).In Fig. 6 the black dashed line corresponds to a power-law fit of the maximum magnetic field B max measured during an ICME encounter performed by Winslow et al. (2015).The authors found a power law with an exponent of α B = −1.89± 0.14 for the evolution of B max as a function of the distance with the 95% confidence interval.We see that the six simulated flux ropes follow roughly a similar trend.Finally, the symbols correspond to the maximum magnetic field measured by Messenger, Venus Express, Parker Solar Probe, and Solar Orbiter during an ICME encounter using the Winslow et al. (2015), Good & Forsyth (2016) and the HELIO4CAST ICMECAT v11 23 catalog (Möstl et al. 2017(Möstl et al. , 2020)).We point out that the thick1 FR rotates according to the 3D visualization displayed in Fig. 5.This might introduce a bias in the automated detection of the position of the flux rope (and thus the measurement of B).Nevertheless, this FR is found to rotate slowly after 15 R , meaning that the magnitude of B might Notes.We recall that Winslow et al. (2015) found α = −1.89± 0.14 based on maximum B measurements during ICME encounters with in situ data.
Table 3. Power-law exponent for the evolution of the magnetic field, the speed, the density, the temperature, and the β (when applicable) with distance from the Sun, and their errors averaged over the six TDm flux ropes, from the Scolini et al. ( 2021) study (third row) and from the Liu et al. (2005) study and the associated error (fourth and last row).
Source In order to compare the evolution of the TDm flux rope with the observations more quantitatively, we perform a power-law fit of the evolution of the simulated B as a function of distance.The corresponding slope exponents (α) are reported in Table 2.
We see that the slopes of all six TDm flux ropes fit within the 95% confidence interval of the Winslow et al. (2015) power law.There is thus a good agreement between the power law deduced from the observations and the magnetic field evolution of TDm flux ropes.We also find that the greatest magnetic field strengths in our sample, the thick3 (purple line), lies on the upper range observed for ICMEs.This means that the magnetic field of the simulated ICMEs have a realistic, if somewhat high, magnetic field strength at different distances from the Sun between Mercury and the Earth.To summarize, we find that the eruption of an initially out of equilibrium flux rope and its propagation in a quasi-isothermal solar wind produces magnetic field amplitudes that agree with the in situ observations of ICMEs at the different distances from the Sun.
Proceeding similarly, we compared the evolution of the magnetic field, the speed, the density, the temperature, and the β parameter with distance to the in situ measurements made during ICME encounters from Liu et al. (2005), but also to the result of Scolini et al. (2021), which used EUHFORIA simulations with a spheromak magnetic structure (Verbeke et al. 2019).Table 3 shows the fit parameter for the ICMEs we simulate (averaged over all the FRs), and gives the associated error on the row just below.The next row shows the fitting parameters from Scolini et al. (2021).Finally, the last two rows correspond to the fitting parameters and the associated error from Liu et al. (2005).Liu et al. (2005) studied ICMEs seen by different probes and performed a power-law fit on their magnetic and plasma parameters.We find in Table 3 a good agreement between the results of Liu et al. (2005) and of these simulations for the density and the temperature.The speed of the simulated ICMEs decreases faster than in the observations, but the values are still close.This A14, page 8 of 13 Fig. 7. Synthetic crossings at 1 au of the three thick flux ropes.Shown are (from the top to bottom) the magnetic field and its components (see legend for the color-coding) in nT, the speed in km s −1 , the density in cm −3 , and the β parameter as a function of time.The green and the blue areas correspond respectively to the sheath and the ME.
difference can be explained for instance by the interaction of the FR with the solar wind.In particular, in our model the wind is a bit denser than the solar wind, which makes the propagation of the flux rope more difficult.
There is also a good match between the evolution of B, V, and n with the distance in Scolini et al. (2021) and in the study presented in this paper (see third row of Table 3).Despite the differences between the TDm + PLUTO and the Linear Force Free Spheromak model (LFSS) + EUHFORIA model (see Introduction), the power laws are found to be similar in the two models.This shows the robustness of the result.
However, if we compare the evolution of T and β we find some differences.We find that α T and α β are −0.37 and 1.04, respectively, in the PLUTO simulations, while their values are −1.19 and 0.11 in the EUHFORIA simulations.These differences are probably due to the solar wind model, which is not the same in the PLUTO and EUHFORIA simulations.As described in Sect.2.2, the solar wind model used in this study is quasiisothermal, which means that the temperature is almost constant in the whole simulation domain of PLUTO.This explains why the temperature evolution is different compared to the EUHFO-RIA simulations, which have a more realistic description of the thermal expansion of the solar wind.

Synthetic crossing
We now compare the properties of the propagation of a TDm flux rope in a quasi-isothermal solar wind with in situ profiles made by spacecraft during an ICME encounter.To do so, synthetic crossings are performed through the flux rope during its propagation, and are compared with the observations performed at 1 au by ACE in Regnault et al. (2020).This study used the superposed epoch method on more than 300 ICMEs and computed the most probable profile for different physical parameters, which can then be used as a typical crossing profile to assess the realism of our modeled ICMEs.
Figures 7 and 8 show synthetic crossings at 210 R at (ϕ = π, θ = π 2 ).The green and blue areas show the sheath and magnetic ejecta (ME) areas respectively, as they would have been defined in actual in situ data according to the typical signatures detailed in Sect. 1. From top to bottom, we show the magnetic field and its components in nT, the speed in km s −1 , the density in cm −3 and the plasma β.We observe in these two figures the typical expected ICME signatures, which we now describe in detail.
The beginning of the sheath, a region of compressed and heated plasma, is highlighted by a sudden increase in the speed and the density.This increase is also present for the magnetic field, but it is smoother.We note that the increase in the speed and density are not as sudden as we can observe in the in situ measurements where there is a discontinuity in the data at the beginning of the sheath (see, e.g., ICME examples in Masías-Meza et al. 2016).The smooth increase in V and n might come from the HLL Riemann solver that we chose, which is quite diffusive.Moreover, it might also be caused by the spatial resolution, which can cause numerical diffusion, especially that far from the Sun where the grid is coarser.As mentioned earlier, the finer spatial resolution at 1 au is ≈2.5 R .
More specifically, the density in the sheath is increased compared to the value of the solar wind by a factor of ≈3 for thick1, 4 for thick2, 5.3 for thick3, 4 for thin1, 5 for thin2, and 6.7 for thin3.Janvier et al. (2014) studied the change in properties of the solar wind after shocks detected in the in situ data.In particular, A14, page 9 of 13 Fig. 8. Synthetic crossings at 1 au of the three thin flux ropes.The physical parameters and the color-coding are the same as in Fig. 7. they found that ICME-driven shocks produce an increase in density that lies between 1.4 and 4.6.Hence, cases thick1, thick2, and thin1 agree well with these observations.We also observe that the faster the ICME propagates, the denser the sheath.The solar wind accumulates more in front of the flux rope since it has less time to evacuate on the sides.
We observe that the β parameter decreases after the beginning of the sheath and stays lower than the solar wind value up to the end of the magnetic ejecta.However, in Regnault et al. (2020) the authors found that the plasma β remains almost constant in the sheath and starts to decrease only in the ME.Moreover, we find a value for β that is very high in the simulation compared to that from the in situ data.This is due to the effect of our quasi-isothermal approximation, which leads to unrealistically high thermal pressures and increases the value of β.At 1 au, the pressure of the modeled ICMEs are about four orders of magnitude too high compared to Masías-Meza et al. ( 2016), Nieves-Chinchilla et al. ( 2018) and Regnault et al. (2020).A β parameter that remains >1 inside the ME for all the FR (instead of <1 as observed in the in situ profile) means that in the simulation the plasma forces are dominant (compared to the magnetic forces).For instance, this extra plasma pressure could cause an overexpansion (compared with expansion caused by the magnetic pressure only).These issues are the main weaknesses of this model.However, they could be alleviated using a more realistic solar wind heating process (e.g., based on Alfvén-wave heating, see Réville et al. 2020), which is beyond the scope of the present paper and will be explored in future studies.
The synthetic B profile of thick FRs shows a double bump structure in the ME, one right after the beginning of the ME and a second one close to the end.Conversely, thin cases show only one bump.However, we observe the smooth rotation of the θ component of the magnetic field in each case, showing the sig-nature of a magnetic flux rope.The beginning of the ME corresponds to the start of the rotation of B θ .The fact that the rotation is smooth and almost symmetric (e.g., from ≈20 nT to −20 nT for the thick2 case), and the fact that B ϕ is at its maximum (in absolute value) almost when B θ = 0, matches very well with the idea of a flux rope according to Lundquist (1951).This agrees with the similarities discussed in Sect.3.2 between the part of the streamlines that are close to the FR for the thin and thick FRs.However, this local agreement does not imply that the two magnetic structures are actually similar, as we can see in Fig. 5.
We also observe that B θ crosses 0 almost when B ϕ is at its maximum at the very beginning of the thin MEs.However, we see that the B θ rotation for the thin FRs is clearly asymmetric, as pointed out in Sect.3.2.
The beginning of the ME closely matches the location where the density goes back to a value similar to that of the solar wind for the thick1 and thick2 cases; however, this is not the case for the thin MEs.This could be explained by reconnection happening between the sheath and the ME magnetic field allowing plasma to enter the FR.
We can see that the magnitude of B of thin1 and thick1 is similar, ≈20 nT at the front part of the ME, while the initial B presented in Table 1 shows that the thick1 case is roughly two times more magnetized than the thin1.The same behavior is observed for the thick2 and thin2 FRs.This can be explained by a greater expansion of the magnetic structure of thin FRs during their propagation and thus a greater decrease in the inner B.
Regnault et al. (2020) used the relative speed of the ICME compared to the solar wind to build different groups of ICMEs.They found that the "relatively fast" ICMEs (compared to the local solar wind) have an asymmetric magnetic field profile with a higher magnetic field at the front of the ME, as shown in their Fig. 4.However, in our simulations, where all six of the A14, page 10 of 13 simulated FRs would be considered relatively fast, we do not find such an asymmetry.
We observe that the speed follows a monotonically decreasing profile.This is due to the radial expansion of the flux rope, as often observed in the in situ data (see, e.g., ICME examples in Farrugia et al. 1993;Masías-Meza et al. 2016;Regnault et al. 2020).Moreover, the density goes back to pre-ICME solar wind values for most of the magnetic ejecta, as is observed in Regnault et al. (2020).The part of the ME that shows the smooth rotation closely matches the region of the lowest density for the thin1 and thick1 MEs.This low density probably highlights the part of the magnetic ejecta that has the shape of a FR, as suggested by the smooth rotation of B θ .The remainder of the magnetic ejecta corresponds to the eroded part of the flux rope due to reconnection with the solar wind, as described by Dasso et al. (2006) and Ruffenach et al. (2012) in their Figs. 1 and 6, respectively.As shown in Fig. 4, the initial B has an impact on the velocity of the FR during its eruption.This relation between the initial B and the speed is still clear at 1 au in one set of the simulated FR.
In Regnault et al. (2020) the authors found that the mean, the median, and the most probable value of the relative duration of the ME compared with the sheath range from 1.5 to 3.5.In our simulation, the duration of the sheath and the ME are shown in days at the top of each column in Figs.7 and 8.By computing the ratio of the ME to sheath duration, we find 2.6, 3.2, and 5.0 for the thick1, thick2, and thick3 events, respectively, while we have 3.8, 6.55, and 5.2 for the thin1, thin2, and thin3 events.We thus find a good match between the simulation and the observation, except for the thick3, thin2, and thin3 FRs which have a higher ratio than the average.Finally, we also observe that the sheath is shorter (compared to the ME) for faster FRs for the thick set.We interpret that as an enhanced compression of the sheath due to the higher speed.This behavior is less clear for the thin set if we take the duration of the sheath, using the B θ rotation as a start for the ME.However, we can see that the width at half height of the density peak is smaller for faster thin FRs, suggesting a higher compression for faster thin FRs than for the thick FRs.

Effect of rotation on the identification of the magnetic ejecta
As illustrated in Fig. 5, the thick1 FR rotates during its propagation, while the thin FR axis remains aligned with the equator.The relative orientation of thick FRs with respect to thin FRs is of importance if we want to compare the components of B.
Figure 9 shows the effect of the rotation on the B components at 1 au.From left to the right, B and its components are shown for the thin2 FR, for the thin2 with a 60 • rotation on the components and finally for the thick2 case for comparison.The 60 • degree angle is found simply finding the angle that maximizes the symmetry (estimated by eye) of the B θ rotation looking at angles between 0 and 90 • with steps of 10 • degrees.We note that ±10 • degrees does not significantly change the results, which we describe below.
The rotation of the magnetic field components of the thin FRs (illustrated here only for thin2) allows us to obtain synthetic crossings that are very similar between the thin and the thick FRs.In both cases (middle and right panels) the B θ rotation is symmetric and B ϕ is at its maximum when B θ is close to 0.
We also note that B θ starts its rotation later.Instead of starting at t = 1.92 days, it starts at t = 2.13 days.Since we use the B θ smooth rotation as a parameter to define the ME boundaries, the rotation of the magnetic field component has an impact on their definitions.In the new coordinate system, the rotated ME has a lower density compared to its non-rotated version.This lower density is in better agreement with the typical sheath-ME transition, as observed in the typical profiles (Regnault et al. 2020).Finally, it also produces a lower ME to sheath duration ratio, which allows the thin FRs to match the Janvier et al. ( 2014) study more closely.
Figures 5 and 9 show that, in our simulations, the flux ropes that are tilted with respect to the equator (the thick cases, according to the 3D visualization) show FR signatures that are clearer than the FR that is aligned with the equator while crossed at the nose.If we apply a rotation on the FRs that are aligned with the equator, we recover clearer FR signatures.
We see here that the synthetic crossing of the B components at the nose of the FR shows similar features between the thin3 and the thick1 FRs if we apply a rotation to one of the FRs.Nevertheless, the magnetic structure presented in Fig. 5 shows significant differences in the sides of the FR.This highlights that the resulting profiles are very dependent on the crossing location within the structure.To lift this ambiguity, multiple crossings of the ICME are necessary.

Summary and discussion
We presented a novel numerical model that simulates the propagation of a TDm flux rope in a quasi-isothermal solar wind from the low corona to 2 au.We applied this model to a series of TDm considering for now a solar wind with a dipolar magnetic field configuration resembling solar minimum phases.The flux ropes are initialized out of equilibrium such that they erupt quickly after their insertion.
During their propagation, we compare the evolution of physical properties of the ICMEs (such as their magnetic field, speed, and density) with evolution laws of these parameters deduced from in situ observations.The evolution of the magnetic field of the simulated ICMEs follows the same trend as the evolution of B max measured by space probes during an ICME encounter at different distances (Winslow et al. 2015).In addition, if we directly compare the value of the magnetic field in the simulations with B max measured by different probes (such as Messenger, Venus Express, Solar Orbiter, and Parker Solar probe) we find that the magnetic field of the simulated ICMEs is comparable to actual measurements in ICMEs, with the thick3 flux rope that lies on the upper edge of the range of observations.We also find a good agreement on the evolution of the density and the temperature with the study of Liu et al. (2005).
We find a good match between the evolution of B, V, n with the distance between the PLUTO + TDm simulation and the EUHFORIA + spheromak simulation.This shows that even if the initial magnetic structures are different between the two simulations and are initialized at different distances, the evolution of B in the heliosphere is similar.However, we do find differences in the temperature and the β parameter which are due to the quasi-isothermal assumption done in this study and not in the simulations of Scolini et al. (2021).
The change in the initial thickness of the TDm flux rope has a significant impact on the properties of the flux rope while it propagates.The first main effect is that flux ropes with different initial sizes are subject to different rotations.Thicker flux ropes rotate more in this study, which suggests that the kink instability is not the main process at the heart of the observed rotation.We interpret the rotation as a consequence of reconnection happening between the magnetic structure and the ambient solar wind.The main part of the rotation takes place in the low corona (<15 R ).
A14, page 11 of 13 .Synthetic crossing at 1 au of the thin2 and the thick2 From the left to the right, B and its components (same color code as in Fig. 7, as for the colored area) are shown for the thin2 FR, for the thin2 with a 60 • of rotation on the components and finally for the thick2 case for comparison.
We also compare synthetic crossings done close to the nose of the simulated ICMEs with the superposed-epoch analysis of Regnault et al. (2020).Overall, the typically expected ICME signatures are observed in the simulation.We observe an increase in the magnetic field, the speed, and the density in the sheath area (with a smooth increase for B) due to the compression of the solar wind plasma that accumulates in front of the propagating flux rope.However, we find that the increase in B at the beginning of the sheath is smoother than for n and V.In the magnetic ejecta we find a lower density compared to the sheath area and a monotonic decrease in the speed, which are interpreted as signs of the expansion of the simulated flux rope during its propagation.We observe a smooth rotation of B θ , which is the trace of the TDm flux rope that is still observable at 1 au.The β parameter starts to decrease at the beginning of the sheath, while the β value of the sheath of the superposed epoch profile in Regnault et al. (2020) remains the same as the solar wind.We also see that its value is higher than that measured at 1 au.This is due to the quasi-isothermal assumption that produces a strong thermal pressure over the whole simulation domain, and thus increases the β parameter.This highlights the need for a better solar wind model, for example the WindPredict-AW model (Réville et al. 2020) that has a more realistic thermal treatment of the solar corona.
The two sets of flux rope thin and thick do not show any clear difference in their 3D magnetic structure close to their noses, except for their rotation during their propagation.Similarly, the synthetic crossings do not show any clear difference between the two sets.Even the components of the magnetic field (once the rotation is corrected) are fairly similar.There is thus no clear effect of the initial size of the flux rope on the observational properties at 1 au close to the nose of the FR.The duration of the modeled events is similar at 1 au disregarding their initial size.However, differences are expected between the two sets of FRs when far from the nose, highlighting the effect of the crossing on ICME properties.
The initial magnetic field does not significantly impact the rotation of the flux rope during its propagation, as FR sharing the same size exhibit the same rotation pattern.Moreover, the change in initial magnetic field have an impact on the speed of the simulated ICMEs due to greater Lorentz forces happening during the eruption of the flux rope.It also leads to a higher compression of the solar wind plasma, and thus a higher density in the sheath area.
To conclude, the eruption and propagation of an initially out of equilibrium TDm flux rope manages to produce a magnetic structure that has the typical signatures seen in the in situ data.The results obtained during this study highlight the difficulty in determining the 3D magnetic structure of an ICME because of the degeneracy of unique in situ profiles for each ICME.It thus highlights the need for multi-spacecraft observations of ICMEs in order to determine their global properties.Thanks to the launch of the Parker Solar Probe (Fox et al. 2016), Solar Orbiter (Müller 2020), andBepi Colombo (Benkhoff et al. 2021) spacecraft, we have now more chances to characterize ICME propagation with multi-point in situ observations, as described in (Hadid et al. 2021).One example of such a study using Solar Orbiter, Bepi Colombo, and Wind is carried out in Davies et al. (2021).
As a follow-up to this study, more numerical work will be done, first to insert the TDm FR in a more realistic solar wind and second to perform the relaxation of the flux rope before its eruption is triggered.This could allow us to study in more detail the physical mechanism that leads to the eruption of such a flux rope in a full 3D MHD model, as well as the initial rotation in the low corona when the erupting structure interacts with the ambient magnetic field and how its topology impacts the property of the propagating FR.

Fig. 1 .
Fig.1.Snapshots of the equatorial slice of the PLUTO simulation with the propagating TDm flux rope.The color gradient shows the velocity magnitude in PLUTO units (1 ≈ 438 km s −1 ).The orange areas correspond to level 4 of the AMR grid, purple to level 3, and red to level 2. Level = 1 does not appear here.

Fig. 2 .
Fig.2.Two-dimensional schematic of the TDm setup.The yellow line corresponds to the inner boundary of the simulation domain and the blue line corresponds to the edge of the flux rope of minor radius a and major radius R. The d parameter controls how much the flux rope is buried in the solar surface(Titov et al. 2014).The magnetic field of the twisted flux rope appears in green.B ⊥ correspond to the ambient magnetic field.

Fig. 3 .
Fig. 3. Two visualizations of the initiation of the TDm flux rope in the solar wind model.The left panel shows the global magnetic field of the solar wind with white streamlines.The sphere corresponds to the inner limit of the simulation domain r = R and the colors correspond to B r .Blue means B r < 0 and red B r > 0. The right panel shows a zoomed-in view of the TDm flux rope.The colored streamlines correspond to the magnetic field lines of the flux rope starting from a circle of radius 0.9a chosen at the footpoints of the flux rope in the negative polarity.The transparent red-yellow slice shows log( J B ) (with J = ∇ × B the current density) in the plane perpendicular to the torus plane.The yellow ring shows the annular distribution of the current specific to case 1 of the TDm.

Notes.
From left to right: Name of the case, magnetic field intensity at the core of the initial flux rope in G, magnetic flux of the flux rope in units of 10 21 Mx, minor radius in R of the torus, and value of the ∆ parameter.

Fig. 4 .
Fig. 4. Visualization of the eruption of the thick2 flux rope.The spherical slice and the streamline colors are the same as in Fig. 3.Here the white streamlines start from the θ = π 2 , ϕ = π line along the radial direction.Panel a) shows the flux rope 1 min after (physical time) its eruption and panel b) 8 min after.

Fig. 5 .
Fig. 5. Evolution of the FR magnetic structure during its propagation.Panel a): 1D cuts of B θ (in green) and B ϕ (in red) at the nose of the thick1 (left) and thin3 (right) flux rope at 15, 30, 50, and 210 R (from top to bottom).The black dashed line for each FR at 50 R is a Gaussian fit of B ϕ (see Sect. 4.1).The dots show the position of the seed streamline displayed in panel b) (of the corresponding color) projected on the B ϕ profile.The x-axis decreases from left to right to facilitate comparison with Figs. 7 and 8. Panel b): Magnetic structure of the thick1 (left) and thin3 (right) cases.All streamlines originate from seeds that are on the (θ = π 2 , ϕ = π) line.The radial position of each seed is shown in panel a).From purple to dark green, the streamlines start from the front of the flux rope and move to the back.
Fig.9.Synthetic crossing at 1 au of the thin2 and the thick2 From the left to the right, B and its components (same color code as in Fig.7, as for the colored area) are shown for the thin2 FR, for the thin2 with a 60 • of rotation on the components and finally for the thick2 case for comparison.

Table 1 .
Flux rope properties for the six different initial TDm.

Table 2 .
Slopes of the linear regression and associated error bars of the evolution law of the magnetic field for the six cases in this study.