Exploring the biases of a new method based on minimum variance for interplanetary magnetic clouds

Magnetic clouds (MCs) are twisted magnetic structures ejected from the Sun and probed by in situ instruments. They are typically modeled as flux ropes (FRs). The determination of the FR global characteristics requires the estimation of the FR axis orientation. Among the developed methods, the minimum variance (MV) is the most flexible, and features only a few assumptions. However, as other methods, MV has biases. We aim to investigate the limits of the method and extend it to a less biased method. We first identified the origin of the biases by testing the MV method on cylindrical and elliptical models with a temporal expansion comparable to the one observed in MCs. Then, we developed an improved MV method to reduce these biases. In contrast with many previous publications we find that the ratio of the MV eigenvalues is not a reliable indicator of the precision of the derived FR axis direction. Next, we emphasize the importance of the FR boundaries selected since they strongly affect the deduced axis orientation. We have improved the MV method by imposing that the same amount of azimuthal flux should be present before and after the time of closest approach to the FR axis. We emphasize the importance of finding simultaneously the FR axis direction and the location of the boundaries corresponding to a balanced magnetic flux, so as to minimize the bias on the deduced FR axis orientation. This method can also define an inner flux-balanced sub-FR. We show that the MV results are much less biased when a compromise in size of this sub-FR is achieved. For weakly asymmetric field temporal profiles, the improved MV provides a very good determination of the FR axis orientation. The main remaining bias is moderate (lower than 6 degrees) and is present mostly on the angle between the flux rope axis and the plane perpendicular to the Sun-Earth direction.


Introduction
The Sun release mass and magnetic field in a permanent solar wind, and also as transients called coronal mass ejections (CMEs). Coronal remote white-light observations using coronagraphs have shown that the distribution of mass in some CMEs present images consistent with twisted structures (e.g., Krall 2007;Gopalswamy et al. 2013;Vourlidas et al. 2013;Wood et al. 2017). When CMEs are observed in the interplanetary medium they are called interplanetary CMEs (ICMEs). The association between CMEs and ICMEs was well established since several decades (e.g., Sheeley et al. 1985). In fact, twisted flux tubes, or flux ropes (FRs), are present also in several other systems in the heliosphere, such as the Sun atmosphere, the solar wind, and different locations of planetary magnetospheres (e.g., Fan 2009;Imber et al. 2011;Smith et al. 2017;Pevtsov et al. 2014;Kilpua et al. 2017). Flux ropes can store and transport magnetic energy and, because their magnetic field lines can be strongly twisted, Send offprint requests to: P. Démoulin FRs can also contain and transport important amounts of magnetic helicity (e.g., Lynch et al. 2005;Dasso 2009;Sung et al. 2009;Démoulin et al. 2016).
When FRs are observed in situ by a crossing spacecraft, they present a large and coherent rotation of the magnetic field vector. A particular set of events that present these FR characteristics corresponds to magnetic clouds (MCs). They are also characterized by a stronger magnetic field and a lower proton temperature than the typical solar wind (e.g., Burlaga et al. 1981;Gosling 1990). Thus, they have a low plasma beta, so that magnetic forces are expected to be dominant.
Interplanetary MCs have been systematically observed from the 80's, and several models have been proposed to describe their magnetic structure. The simplest one and generally used to model the field in MCs is an axially symmetric cylindrical magneto-static FR solution, with a relaxed linear force-free field, the Lundquist's model (Lundquist 1950;Goldstein 1983). This model describes relatively well the field distribution for a significant number of observed MCs (e.g., Burlaga & Behannon 1982; Article number, page 1 of 15 arXiv:1809.00522v1 [astro-ph.SR] 3 Sep 2018 A&A proofs: manuscript no. demoulin_MV_arXiv Lepping et al. 1990;Burlaga 1995;Burlaga et al. 1998;Lynch et al. 2003;Dasso et al. 2005;Lynch et al. 2005;Dasso et al. 2006). Still, the Lundquist solution is known to have difficulties in fitting the magnetic field strength, in particular it was found that it frequently overestimates the axial component of the field near the FR axis (e.g., Gulisano et al. 2005).
Several other models have been developed. These models include more general FR properties, such as non-linear forcefree field (Farrugia et al. 1999), non-force free magnetic configurations (e.g., Mulligan et al. 1999;Hidalgo et al. 2000Hidalgo et al. , 2002Cid et al. 2002;Nieves-Chinchilla et al. 2016), and models with non-cylindrical cross section (e.g., Vandas & Romashets 2003;Nieves-Chinchilla et al. 2009, 2018a. In particular, the elliptical model of Vandas & Romashets (2003) provides a better fit to some observed MCs having a field strength more uniform than in the Lundquist solution. This indicates the existence of some flat FRs (Vandas et al. 2005). On the other hand, from a superposed epoch analysis, Masías-Meza et al. (2016) showed that while slow MCs present a symmetric profile of the magnetic field, fast ones present a magnetic intensity profile with a significantly stronger B near the front than at the rear.
The models presented above were typically compared and fitted to MC in situ observations, which are limited to the 1D cut provided by the spacecraft trajectory inside the FR. This allows a local reconstruction of the FR cross section, and then it is possible to make estimations of global magnetohydrodynamic (MHD) quantities, such as magnetic fluxes, twist, helicity, and energy (e.g., Dasso et al. 2003;Leamon et al. 2004;Dasso et al. 2005;Mandrini et al. 2005;Qiu et al. 2007;Démoulin et al. 2016;Wang et al. 2016). However, all these estimations depend directly on the FR orientation, and thus on the quality of the method used to get it. Moreover, the orientation itself is an important property of a given FR in order to compare it with its solar origin (e.g., Nakwacki et al. 2011;Isavnin et al. 2013;Palmerio et al. 2017).
Different methods to estimate the FR orientation from spacecraft observations have been developed and applied to MCs. Some authors have used the so-called Grad-Shafranov method to get the FR orientation, which consists in applying the Grad-Shafranov formalism, valid for describing general MHD magnetostatic equilibria invariant in one direction (Sonnerup & Guo 1996;Hau & Sonnerup 1999;Sonnerup et al. 2006;Isavnin et al. 2011). This method use both the magnetostatic constraints and the folding of the same magnetic field line when it is crossed by the spacecraft during its in-bound and out-bound travel across the FR. However, the method cannot recover the FR axis for the simplest magnetic configurations, such as symmetric FRs or when the spacecraft crossed the FR close to its axis. Hu & Sonnerup (2002) argued that observed FRs are typically asymmetric and that this asymmetry removes the above difficulty. In fact a detailed study, Démoulin et al., in preparation, show that for FRs typical of MCs the Grad-Shafranov method provides a biased orientation, which depends on the symmetry properties of the magnetic field components.
Finally, the method that has less hypotheses, mainly because it does not introduce the details of a given model, is the Minimum Variance method (MV). The MV method just requires a well ordered large scale variance of B in the three spatial directions. MV has been extensively used to find the orientation of structures in the interplanetary medium (see e.g., Sonnerup & Cahill 1967;Burlaga & Behannon 1982;Hausman et al. 2004;Siu-Tapia et al. 2015). Several authors have shown that the MV method estimates quite well the orientation of the FR axis, when the distance between the axis and the spacecraft trajectory in the MC is small with respect to the FR radius (see e.g., Klein & Burlaga 1982;Bothmer & Schwenn 1998;Farrugia et al. 1999;Xiao et al. 2004;Gulisano et al. 2005Gulisano et al. , 2007Ruffenach et al. 2012Ruffenach et al. , 2015. Other authors use the MV method to get a first order approximation for the MC orientation. Then, they use this estimation as a seed in a non-linear least squares fit of a magnetic model to the data. This approach is expected to improve the cloud orientation (e.g., Lepping et al. 1990Lepping et al. , 2003Lepping et al. , 2006Dasso et al. 2005).
The most crucial information to introduce in any of the methods to get the FR orientation is the start and end time of the FR. This is a major problem because different authors generally define the FR boundaries at different places/times (e.g., Riley et al. 2004;Russell & Shinde 2005;Dasso et al. 2006;Al-Haddad et al. 2013), with the consequent difference on the determination of both the axis orientation and the parameters for the proposed models. One of the reasons at the origin of these different boundaries is the partial erosion of the FR due to magnetic reconnection. If it happens in the FR front, this creates after the closed FR a 'back' which presents mixed signatures of a FR and stationary solar wind (Dasso et al. 2006). This erosion process also can occur at the FR rear, producing a mixed region before the beginning of the real closed FR (Ruffenach et al. 2015). This lack of exact information on the FR boundary can have an important influence on the proper determination of the FR axis direction. Indeed, comparing different methods show a large dispersion of the estimated FR orientation (Al-Haddad et al. 2013;Janvier et al. 2015).
In Sect. 2 we present a geometrical description of FRs in the solar wind, proposing new angles to determine its orientation, and review the MV method and its application to get the orientation of MCs. In Sect. 3 we generate synthetic linear force-free FRs with symmetry of translation along their axis, including the possibility of expansion and cylindrical/elliptical cross-sections, and emulate the observed time-series observed by spacecraft. Then, in Sect. 4 the application of several variants of the MV technique are applied to the synthetic clouds, finding biases produced by different MCs properties. In Sect. 5 a deeper analysis of the bias introduced due to a boundary selection is done, and we proposed a new MV method to minimize this bias. Our conclusions are given in Sect. 7, and in the Appendix we show the expected coupling between field components when MV is applied.

Flux rope axis and frame
The coordinates generally used for the analysis of data provided by spacecraft located in the vicinity of Earth are defined in the Geocentric Solar Ecliptic (GSE) system of reference. It is defined with an orthogonal base of unit vectors (x GSE ,ŷ GSE ,ẑ GSE ).
x GSE points from the Earth toward the Sun,ŷ GSE is also in the ecliptic plane and in the direction opposite to the Earth rotation motion around the Sun, andẑ GSE points to the north pole of the heliosphere (Fig. 1a). A similar coordinate system is the heliographic radial tangential normal (RTN) system of reference (e.g., Fränz & Harper 2002) where the following could also be transcripted.
The FR axis orientation is classically defined with respect to the GSE system in spherical coordinates with the polar axis chosen as +ẑ GSE , by two angles: the longitude (ϕ axis ) and the latitude (θ axis ). However these angles are not the natural ones to describe Definitions of the angles of the FR axis. (a) Schema showing the observation, here GSE, and FR frames. The FR frame is defined by the vectorsx FR ,ŷ FR ,ẑ FR whereẑ FR is along the FR axis and in the direction of the axial magnetic field,x FR is orthogonal toẑ FR and in the plane defined byx GSE andẑ FR , andŷ FR completes the orthonormal direct basis. Two rotations are needed to pass from the GSE frame to the FR frame. They are defined by the angles i and λ which are respectively the inclination and the location angle (spherical coordinates with the polar axisx GSE ). (b) Schema showing the meaning of the angle λ in the plane of the FR axis shown in panel (a) in light blue.ẑ pr,FR , in green, is the projection ofẑ FR in the plane orthogonal tox GSE , see panel (a).
the FR axis orientation with respect to its translation direction (≈ −x GSE ). Indeed, the geometrical configuration of the spacecraft crossing is defined both by the closest approach distance and the angle between the spacecraft trajectory and the FR axis orientation. A rotation of the FR around the spacecraft trajectory does not change the geometry of the spacecraft crossing, while it changes both ϕ axis and θ axis . This implies that the same crossing geometry is present along curves defined in the {ϕ axis , θ axis } space. Even worse, this choice of reference system has the disadvantage to set the polar axis (θ = ±90 • ) alongẑ GSE which is both a possible and an un-particular axis direction for MCs while the polar axis is singular as it corresponds to any values of ϕ axis . Therefore, the coordinates {ϕ axis , θ axis } are not appropriate to study the orientation of the FR.
The motion of MCs is mainly radial away from the Sun, especially away from the corona (where limited deflection could occur). Then, the radial direction is a particular direction. We then set a new spherical coordinate system having its polar axis alongx GSE (Fig. 1a). This direction could correspond to the spacecraft crossing a FR leg. However, such crossing does not allow to detect the rotation of the magnetic field as the FR is typically only partially crossed on one side. It implies that FRs with an axis direction almost parallel tox GSE are mostly not present in MC data sets. Then, the directionx GSE can be used as the polar axis of the new spherical coordinate system.
The local axis direction of the FR is calledẑ FR (with B z,FR > 0 in the central region of the FR). We define the inclination on the ecliptic (i) and the location (λ) angles as shown in Fig. 1. More precisely, let us define the unit vectorẑ pr,FR , in green in Fig. 1a, along the projection of the FR axis on the plane defined by (ŷ GSE ,ẑ GSE ), then i is the angle fromŷ GSE toẑ pr,FR . Similarly, λ is the angle fromẑ pr,FR toẑ FR .
If the FR axis is located in a plane, i is the inclination of this plane (in light blue) on the ecliptic as shown in Fig. 1a. And if we further suppose that the distance to the Sun is monotonously decreasing from the FR apex to any of the FR leg when following the FR axis, the angle λ is evolving monotonously along the FR from −90 • in one of the FR leg, to λ = 0 at the apex, then toward λ = 90 • in the other leg. It implies that λ can be used implicitly as a proxy of the location where the spacecraft intercepts the FR (as shown in Fig.1c of Janvier et al. 2013). Finally, the above conclusions extend to other coordinate systems such as RTN.
We next define the FR frame.ẑ FR corresponds to the FR axis, as defined before. Since the speed of MCs is much higher than the spacecraft speed, we assume a rectilinear spacecraft trajectory defined by the unit vectord. We defineŷ FR in the direction z FR ×d andx FR completes the right-handed orthonormal base (x FR ,ŷ FR ,ẑ FR ) which defines the FR frame (Fig. 1a).

Minimum variance technique applied to flux ropes
The MV method finds the directions in which the projection of a series of N vectors has an extremum mean quadratic deviation (e.g., Sonnerup & Scheible 1998). This method can be applied to the time series of the magnetic field B measured in situ across MCs, and it provides an estimation of the FR axis orientation as follows.
The mean quadratic deviation, or variance, of the magnetic field B in a given direction defined by the unit vectorn is: where the summation is done on the N data points of the time series and B is the mean magnetic field. The MV method finds the directionn where σ 2 n is extremum. The constraint |n| = 1 is incorporated with the Lagrange multiplier variational method. The resulting set of equations in matrix form is where λ L is the Lagrange multiplier and The matrix m i j is symmetric with real positive eigenvalues: 0 ev min ev int ev max . Their corresponding orthogonal eigenvectors (x MV ,ẑ MV ,ŷ MV ) are the directions of minimum, intermediate and maximum variation of the magnetic field, respectively.
The magnetic configuration of MCs is assumed to be a FR with some generic properties which have direct implications on the eigenvalues and eigenvectors found by the MV applied to MC data as follows. First, the axial field, B z,FR , is stronger on the FR axis and decline to a low value near the boundary. Second, the azimuthal field vanishes on the axis (to have no singular electric current) and typically grows to a magnitude comparable to the axial field strength at the FR periphery.
For trajectories passing close enough from the FR axis, the azimuthal field is mostly in theŷ FR direction. Then, the variance of B y,FR is expected to be about twice the one of B z,FR . The variance of B x,FR is the lowest one. Then the eigenvectors (x MV ,ŷ MV ,ẑ MV ) provide an estimation of (x FR ,ŷ FR ,ẑ FR ), taking into account that the direction of the eigenvectors provided by MV could need a change of sign in order to satisfy the convention defined before for the FR frame. Next, sinceŷ FR is perpendicular tod (defining the spacecraft trajectory), then the angle Article number, page 3 of 15 A&A proofs: manuscript no. demoulin_MV_arXiv betweenŷ MV andd allows to test how close is the MV frame from the FR frame. Finally, as the spacecraft trajectory is further from the FR axis, this estimation of the FR frame is expected to be worse. This is quantitatively tested in Sects. 4-6 with the models presented in Sect. 3.

Test models
We describe in this section the geometrical aspects of the spacecraft trajectory, and the simulation of measured magnetic field along the trajectory. The magnetic models include key ingredients observed or consistent with observations of MC: typical magnetic profiles, expansion, flatness of the FR cross section and asymmetry due to magnetic reconnection. These models are described in the FR frame, with the origin of the reference system set at the FR axis.

Geometry of the flux rope crossing
During the spacecraft crossing of the FR, its small acceleration has a negligible effect on the measured quantities (i.e., on the time where they are measured, Démoulin et al. 2008). Then, we consider a motion with a constant velocity V c during the spacecraft crossing. Supposing that the FR moves along −x GSE , the spacecraft position, r S /C (t), in the FR frame is: where t is the time with t = 0 set at the time when the spacecraft has the closest distance, |y p |, to the FR axis. Thus, t < 0 corresponds to the in-bound trajectory (i.e., when the spacecraft is entering the FR going toward its axis) and t > 0 to the out-bound one (i.e., when the spacecraft is going away from the FR axis).
In order to provide a generic description of a FR crossing we suppose here that the 3D magnetic field evolution, B FR (x, y, z, t), is known in the FR frame (with an analytical or a numerical simulation). This magnetic field is transformed to the GSE frame by two rotations (with angles i FR and λ FR ). Then, inserting the spacecraft trajectory r S /C (t), Eq. (4), in B GSE (r), limits the description of B to the one observed along the spacecraft trajectory, noted B obs (t), which is only a function of time. After sampling the time uniformly, this time series models the synthetic magnetic field observed by a spacecraft crossing the magnetic field configuration comparable to in situ data obtained in MCs.
The direct MV method, Sect. 2.2, and later on including some additional procedures, Sect. 5, are applied below to the above synthetic field. The estimated angles, i MV and λ MV , are next compared to i FR and λ FR allowing to test the performance of the MV method on an ensemble of models by scanning the parameter space.

Expanding magnetic field
We consider a general magnetic field configuration, B 0 (x 0 , y 0 , z 0 ), associated to an element of fluid located at (x 0 , y 0 , z 0 ) and defined at time t 0 . Below, to be coherent with Eq. (4), we simply set t 0 = 0 at the time when the FR axis distance to the spacecraft is minimum. The self-similar isotropic expansion by the factor e(t) implies x = e(t) x 0 , y = e(t) y 0 , z = e(t) z 0 where x, y, z are the new coordinates of the same element of fluid, but at time t, with e(0) = 1. Only an isotropic expansion is considered here since it preserves the force balance of force-free configurations while keeping the magnetic field structure, while a non-isotropic expansion introduces forces which deform the FR so that it would require a numerical simulation to determine the corresponding expanded configuration.
The expanded magnetic field, B, is function of space and time: The term e 2 (t) at the denominator is included to preserve the magnetic flux. Equation (5) describes the temporal evolution in the 3D space. Introducing Eq. (4) into Eq. (5), this limits the description of B to the observed field along the spacecraft trajectory with B obs (t) only function of time. The expansion is driven by the decrease of total pressure in the surrounding of the FR as it is moving away from the Sun (Démoulin & Dasso 2009a). This total pressure can be approximated by a power law of the distance from the Sun to the FR center (D(t)). The force balance at the FR boundary implies that e(t) is also a power law of D(t): e(t) = (D(t)/D 0 ) ζ where ζ is the normalized expansion factor (e.g., Démoulin et al. 2008) and D 0 is a reference distance taken here as the spacecraft distance from the Sun. The velocity of the FR core, V c , is approximately constant during the spacecraft crossing, as in Eq. (4), then D(t) = V c t + D 0 with the spacecraft being located at D = D 0 at t = 0. In this frame work e(t) is: with which corresponds to the time difference between the observed time and the time when the axis reaches D 0 , normalized by the travel time from the Sun to the distance D 0 at velocity V c . The linear approximation on the right side of Eq. (6) provides a good approximation for most MCs because ζ ≈ 1 and τ max << 1 with τ max being half of the FR crossing time (Démoulin et al. 2008). It implies that the expansion typically introduces only a small modification to the observed field profile compared to a case without expansion.

Cylindrical flux rope model
We consider in this subsection a cylindrically symmetric magnetic field configuration for the FR, with its axis along the z direction in the FR frame. Then, at a time t = 0, B 0 (r) is only function of the distance to the axis x 2 0 + y 2 0 , or equivalently of the relative position of an element of fluid from the FR axis, ρ, that is the distance normalized by the FR radius, b 0 . Considering only expansion (e.g., no magnetic diffusion) ρ is preserved during the self-similar isotropic expansion: where b 0 and b(t) are the radius of the FR at t = 0 and t respectively. As other spatial variables (Sect. 3.2), b(t) is simply related to b 0 as b(t) = e(t) b 0 . ρ is a Lagrangian marker of the relative distance to the FR axis. It follows the magnetic field as  the expansion occurs. In the FR frame, the magnetic field can be written simply as The spacecraft rectilinear trajectory satisfies y p = constant. We define the impact parameter as p = y p /b 0 . Using Eq. (4), Eq. (8) is rewritten along the spacecraft trajectory as: where the subscript "obs" has been added to specify that ρ is determined along the simulated spacecraft trajectory so ρ obs is function of time (and it corresponds to the different elements of fluid observed by the spacecraft), while ρ, defined in Eq. (8), is a general Lagrangian marker, so independent of time. Since the simulated magnetic field along the spacecraft trajectory is B 0 (ρ obs (t)) /e 2 (t), the magnitudes of |τ| and ζ determine how strongly the expansion affects the modeled magnetic field (compared to a case without expansion).
In the FR frame, B = (−B θ sin θ, B θ cos θ, B z ) where B θ (ρ) and B z (ρ) are the azimuthal and axial components respectively.
Its components along the spacecraft trajectory are: For a linear force-free field (Lundquist's FR): where B 0 is the axial field strength of the FR axis. We define α as the first zero of the Bessel function J 0 , α ∼ 2.4, then ρ obs = 1 defines the reversal of the axial field component. Observed MCs have typically a weak axial field component at their boundary so ρ obs ≈ 1, while some cases with ρ obs 1 have been observed (Vandas & Geranios 2001).
Examples of simulated magnetic field profiles are shown in Fig. 2a,b. Observations provide ζ ≈ 0.8 ± 0.2 at 1 AU and ζ ≈ 0.9 ± 0.2 and ζ ≈ 1.05 ± 0.3 in the inner and outer heliosphere, respectively, for unperturbed MCs (not overtaken by a fast stream or by another MC, Démoulin et al. 2008;Gulisano et al. 2010Gulisano et al. , 2012. Then, we show the cases ζ = 0, 1, 2, which are respectively the cases without, with typical and with twice as large of an expansion rate than observed. As ζ or/and the FR radius increase, the simulated observed field is more asymmetric. Finally, the case ζ = 1 introduces only a moderate asymmetry in all B components for not extremely large FR radius (with a radius 0.1 AU at 1 AU).
The B x,FR and B z,FR temporal profiles are nearly symmetric in time while B y,FR is nearly antisymmetric. A non-zero impact parameter p introduces a component B x,FR , whose magnitude increases with |p|, see Eq. (11). This property can be used in observations to estimate |p| if a B θ (ρ) profile is assumed (see Sect. 4.3 in Démoulin & Dasso 2009b). Based on different models and the analysis of synthetic FRs, an empirical method to estimate |p| from the observed B x,FR profile was developed in Gulisano et al. (2007).

Elliptical flux rope model
Vandas & Romashets (2003) generalized the linear-force free field of Lundquist to a FR with an elliptical cross-section. Along and across the trajectory, more precisely in the x and y directions of the FR frame, the maximum FR extension is 2 a 0 and 2 b 0 , respectively. The Lundquist's solution is recovered for b 0 = a 0 . Since the encountered or overtaking solar wind tends to compress the FR in the radial direction away from the Sun, then the relative sizes are expected to satisfy a 0 ≤ b 0 .
With invariance along the FR axis, the linear force-free field equations are reduced to the Helmholtz equation Vandas & Romashets (2003) solved this equation with elliptic cylindrical coordinates, one of the few coordinate systems where the Helmholtz equation has separable solutions. They set the magnetic field to be tangential to the elliptical boundary with a vanishing axial component. For all b 0 /a 0 values, they found an analytical solution expressed with the even Mathieu function of zero order.
The above elliptical model can be set in self similar expansion as done above for the Lundquist's field. The sizes become a = a 0 e(t) and b = b 0 e(t). The cartesian components of the elliptical model solution (B x , B y , B z ) are first expressed in function of x, y. Then, the use of Eqs. (4) and (5) provide the simulated crossing of the FR. The impact parameter p is generalized to p = y p /b 0 , with y p being the minimum signed distance to the FR axis. Finally, Eq. (8) is generalized to with a 0 , b 0 defined at t = 0 and a, b at t. The selected boundary condition implies B z,FR (ρ = 1) = 0. Finally, along the spacecraft trajectory ρ obs (t) is still expressed by Eq. (10). An example of field component profiles of this elliptical structure is shown in Fig. 2c with b/a = 2. The main difference with the circular case, Fig. 2b, is a flatter B(t) profile. Indeed, since the magnetic tension is less important due to the FR elongation in the y direction, the outward gradient of the magnetic pressure is lower then B is more uniform.

Defining flux rope boundaries
In the absence of magnetic diffusion and reconnection a given value of ρ (Eq. (8)) traces a selected cylindrical shell within the FR. Along the spacecraft trajectory the front and rear boundaries are located at ρ front and ρ rear , respectively. Below, we consider eroded FRs, with 0 < ρ front ≤ 1 and 0 < ρ rear ≤ 1, in order to include possible events where the peeling of the FR via magnetic reconnection can be present (see, e.g., Dasso et al. 2006;Ruffenach et al. 2012Ruffenach et al. , 2015. The temporal series of B is computed in the time interval [t F , t R ] from Eq. (11), or its equivalent for the elliptical case. Because the closest approach is set at t = 0, Eq. (4), one has the ordering: t F < 0 < t R . t F and t R are found by solving Eq. (10), with ρ obs (t) = ρ front and ρ rear , respectively. We suppose below that p < ρ front and p < ρ rear , so that the FR is crossed by the spacecraft.
The case of the front boundary and ζ ≥ 0 (expanding FR) is simple. ρ obs (t) in the in-bound branch is a decreasing function of time (for t < 0 in Eq. (10), the numerator is a decreasing function of time while the denominator is always increasing with time). This implies that ρ obs (t F ) = ρ front has a unique solution when p < ρ front .
The case of the rear boundary is more case-dependent. This is because both the numerator and the denominator of Eq. (10) are both increasing with time for ζ > 0 and t > 0. For ζ < 1, the numerator is growing faster than the denominator so ρ obs (t R ) = ρ rear has a unique solution (when p < ρ rear ). For larger ζ values there are two t R solutions or even none for ζ 4. The physical solution in the first case is the one closer to the axis while in the second case the expansion rate is so large that the simulated spacecraft never crosses the rear boundary. This last case is not observed as the maximum ζ values observed within unperturbed MCs is 1.5 in the inner heliosphere, 1.1 at 1 AU, and 1.7 in the outer heliosphere (Démoulin et al. 2008;Gulisano et al. 2010Gulisano et al. , 2012.

Main procedure
A classical test of the quality of the MV method to distinguish directions is to analyze the separation between the variance (eigenvalues) associated with the three eigenvalues ev min , ev int and ev max of the matrix m i j defined by Eq. (3) (see references in Sect. 1). Indeed, if two directions have similar variances, so similar eigenvalues, say for example, x and z directions, a rotation of the frame around the third eigenvector (y) is expected to provide also similar variances in the two rotated directions (x and z ). This implies that the x, z directions are not well defined.
However, a large separation of the three eigenvalues is a necessary but not sufficient condition to have a precise determination of the FR axis direction. A systematic bias could be present due to the mixing of the field components to achieve extremum variances in the eigenvector directions. To understand this possible bias, let us analyze the Lundquist's FR with a non-negligible impact parameter p as an example. |B x,FR | has a shape comparable to B z,FR (Fig. 2b). Then, the MV combines B x,FR and B z,FR by a rotation of frame to produce the flattest possible B x,MV (so with the lowest variance). This produces a bias in the MV frame orientation which increases with |p|, confirming the results of Gulisano et al. (2007). Indeed, in general, the MV method will mix components which are correlated, see Appendix A.
More generally, this bias and other ones are present for the axis direction obtained when the MV is applied to expanding FRs. Below, we first identify, then correct, these biases as much as possible. They are measured by the positive and acute angles θ x , θ y , θ z between the x, y, z directions of the MV and FR frames, respectively (these angles can be defined as signed, but the absolute value is sufficient for our purpose). Next, we compute ∆i = i MV − i FR and ∆λ = λ MV − λ FR which quantify the  difference of axis orientations found by the MV and the correct one simulated with the FR model. Finally, we quantify the bias on the axial field by the ratio B z,max,MV /B z,max,FR .
From the tests made, we conclude that a large separation of the three eigenvalues is needed but this is far from being sufficient to have a precise axis determination. For example, for |p| 0.6, we show in Fig. 3b a case where the two lowest eigenvalues are more separated when |p| increases, up to |p| ≤ 0.6, while the bias in the axis orientation increases (θ x , θ z , ∆λ in the left panels).

Minimum variance applied to B orB
The MV technique, described in Sect. 2.2, can be applied to the time series of B orB = B/|B| defined within the modeled FR (and later on with in situ data). These approaches are called MV0 and MV1 and were used, for example, by Siscoe & Suey (1972) and Gulisano et al. (2007), respectively.
MV0, applied to the Lundquist's model has an orientation bias increasing linearly with p (Fig. 3a, see θ x , θ z and ∆λ variations). This bias is reduced with MV1 (Fig. 3b). Next, θ y = ∆i = 0 for this simple test because B x,FR , B z,FR are symmetric and B y,FR is antisymmetric with time so that they cannot be mixed by MV (see Appendix A). These tests agree with the results of Gulisano et al. (2007) and precise the origin of the bias since the Article number, page 7 of 15 A&A proofs: manuscript no. demoulin_MV_arXiv   use of i and λ is more adapted to separate the relevant magnetic components than the longitude and latitude of the FR axis. This further justifies the choice of i and λ angles because, when using them, the bias of the orientation remains mainly on λ.
Moreover, the origin of the biases is illustrated for ζ = 0 with the B components of the Lundquist's FR drawn in the FR frame (black lines) and MV0 frame (orange lines) in Fig. 4a. B x,FR and B z,FR are combined to provide a flat B x,MV (Fig. 4a, left panel). More generally, for a cylindrically symmetric FR, B x = B θ (ρ obs ) p/ρ obs in the FR frame from Eq. (11), then B x = constant if the azimuthal component B θ is linear with the radius ρ obs . In this case, the variance of B x vanishes and the MV0 method would exactly associate this minimum variance direction tox FR . In brief, there would be no orientation biases. However, in general this linear dependence is only approximately true close to the FR center for cylindrically symmetric FR (using a Taylor expansion of B θ ).
The above bias decreases when normalizing B to unity since |B| decreases away from the FR axis, compensating partly the decrease of B θ . It implies that B θ /|B| is more linear with ρ obs than B θ . This effect is directly seen on the B y component for p = 0 since B θ = |B y | and it is still well present for p = −0.5 when comparing the central panels of Fig. 4. Consequently, B x /|B| is more uniform than B x and the MV1 frame is closer to the FR frame than the MV0 one (Fig. 3). Then, the B components in the MV1 frame (green lines in Fig. 4) are closer to the ones in the FR frame (black lines) than when MV0 is applied to B of the Lundquist's field (orange lines).  . The results are shown for two cases of flux unbalance between the FR front and rear. The front limit is located where B z,FR = 0 (ρ front = 1) and the rear limit is set at a fraction ρ rear of the FR radius. For comparison the case ρ rear = 1 is shown in Fig. 3b. Left column: θ x , θ y , θ z are the angles in degrees between the x, y, z directions, respectively, of the FR and MV frames. Middle left column: ∆i and ∆λ quantify the difference of axis orientations between the FR and the one found by MV, so the orientation biases. Middle right column: ratios of the eigenvalues ev min , ev int and ev max . Right column: ratio of the maximum values of B z,MV to B z,FR (amplified by a factor of 100 for the first ratio). The horizontal dotted lines and the yellow regions are landmarks.

Sign and magnitude of biases
The MV bias identified above is due to a mix between B x,FR and B z,FR . By definition of the orientation of the z axis, B z,FR ≥ 0 in the FR core. From Eq. (11), B x,FR has the sign of −p B θ . Then, the bias ∆λ has also the sign of −p B θ . For example, Fig. 3 shows a case with positive helicity, B θ B z ≥ 0, and p ≥ 0 so that ∆λ ≤ 0. Other cases are simply obtained by changing the sign of ∆λ accordingly to the signs of p and magnetic helicity. The effect of expansion for a typical ζ value, ≈ 1, introduces only a small extra bias on i, |∆i| < 4 • , as seen by comparing Fig. 5a to Fig. 3b. The expansion induces an asymmetry in the field components (Fig. 2), then the MV method combines a part of B y,FR with B z,FR , Fig. 6a, introducing a bias for θ y and ∆i values. Other quantities (θ x , θ z , ∆λ, the ratio of eigenvalues and B z,max,MV ) are not significantly affected by the expansion. This small effect of expansion is a generic result for all the models tested.
Increasing the aspect ratio b/a enhances more significantly the bias (Fig. 5b, we note the change of vertical scale on the left panels). In fact, since the B(t) profile gets flatter with increasing b/a value (Fig. 2c, right panel) MV1 is converging to MV0 as b/a is larger, so the bias increases (as shown above with Fig. 3 comparing MV0 to MV1). In contrast, as b/a increases, the ratio ev min /ev int decreases (Fig. 5, compare middle right panels), so this ratio is again not a reliable indicator of the quality of the axis direction determination. Furthermore, as b/a increases, B x,FR (t) behaves closer to B z,FR (t) then both MV0 and MV1 rotate the frame combining B z,FR and B x,FR , trying to get B x,MV as flat as possible (compare the two left panels of Fig. 6). This further increases |∆λ| so the bias in the determined axis orientation (Fig. 5b).
The bias on the axis orientation introduced by the MV implies a systematic overestimation of the axial component B z , as shown in Fig. 4, which increases with |p| (Fig. 3, right panels).
This effect is stronger for MV0 than MV1 as expected from the orientation bias results. This bias on B z also increases slightly with b/a value (Fig. 5, right panels). This is illustrated with an example shown in the right panels of Fig. 6. The MV method rotates the MV frame so that the bump shape of B x,FR is removed in B x,MV . This strengthens B z,MV independently of the sign of p and of magnetic helicity.

Defining the flux rope boundaries
The boundaries of the time interval selected in the in situ data to find the FR orientation are typically set where abrupt changes of the magnetic field and plasma parameters are detected. However, there are frequently some ambiguities in setting them so that the selected boundaries depend on the author's specific criteria. For example, the MC observed on October 18-20 1995 by the Wind spacecraft was studied by several groups of authors which set different MC boundaries (Lepping et al. 1997;Larson et al. 1997;Janoo et al. 1998;Collier et al. 2001;Hidalgo et al. 2002;Dasso et al. 2006).
The main origin of these discrepancies is that the measured magnetic field components and the plasma data, such as proton temperature, plasma composition and properties of high energy particles, do not always agree on the extension of the MC, and more generally of the magnetic ejecta (see, e.g., Russell & Shinde 2005). A part of this ambiguity is due to the partial magnetic reconnection of the ejected solar FR with the magnetic field of the encountered solar wind or of an overtaking stream. This reconnection does not affect the field and plasma on the same time scale. For example a change of magnetic connectivity is affecting on short time scales (on the order of few tens of minutes) the propagation of high energy particles while the mix up of MC and solar wind protons from the distribution core takes days to be mixed after reconnection occurred. The above processes imply that different boundaries are often defined in different studies of the same MC (e.g., Riley et al. 2004;Al-Haddad et al. 2013).
The selection of boundaries have large implications on the derived FR orientation (e.g., Dasso et al. 2006). For example, the axis orientation derived from MV and a least square fit to a Lundquist's model can be significantly different and without coherence along the axis of a MC observed by four spacecraft (ACE, STEREO A and B, Wind) for a set of boundaries (Farrugia et al. 2011), while consistent results are obtained between both methods and along the FR axis when refined time intervals are used . The important effect of the boundaries is further shown in Janvier et al. (2015) by the difference in orientations found with the same set of MCs analyzed by Lynch et al. (2005) and Lepping & Wu (2010), while the FR axes were determined in both cases by the same method (fit of B data to Lundquist's model).

Testing the effect of boundary location
The effect of the boundary selection on the determined axis is illustrated in Fig. 7 using MV1 on a Lundquist's FR without expansion (in order to focus the analysis only on the boundary effects). We set the front at ρ front = 1 (where B z,FR = 0) and allows ρ rear < 1. The FR is crossed only for |p| < ρ rear so the large values of ∆i, ∆λ and B z,max,MV /B z,max,FR obtained for large |p| values in Fig. 3b are absent in figures with ρ rear < 1 and the comparison of biases should be done for the same p value.
Already with ρ rear = 0.9 (Fig. 7a), so only with 10% in radius removed at the FR rear, there is a clear increase of MV1 bias (compare to Fig. 3b where ρ rear = 1 is the only different parameter). θ y , which was null in Fig. 3b, is comparable to θ x and θ z (left panel of Fig. 7a). This introduces mainly a rotation of the estimated FR axis by changing its inclination on the ecliptic plane (∆i > 0). There is a mix up of B y,FR with the other components in order to increase the variance of B y,MV . This is possible as the components are no longer antisymmetric or symmetric with respect to the closest approach time, then they are partly correlated (see Appendix A).
As ρ rear is further decreased the bias on i increases, and dominates the bias on λ for all p values. For example, with ρ rear = 0.5, so with half of the FR rear removed, this introduces the main bias on i with a rotation of the determined FR axis on the order of 24 − 35 • compared to the known axis direction (Fig. 7b, middle left panel). Even worse, this bias on i is also present for low impact parameter cases (bias ≈ 24 • ). This increasing bias is associated with an increasing separation of the two lowest eigenvalues (Fig. 7b, middle right panel). It shows once more that the eigenvalue separation is not a good criteria for estimating the precision of the axis direction defined by MV.
The location of the FR boundaries also introduces biases on other parameters of the FR such as its radius (by changing the geometry of the crossing) or its axial field strength (by mixing the field components). Still, it has only a moderate effect on B z,max,MV /B z,max,FR since MV1 mixes B y,FR with B z,FR and B y,FR ≈ 0 where B z,FR is maximum (Fig. 6), so that B z,max,MV /B z,max,FR is comparable to the case ρ rear = 1 for the same p value (compare the left panel of Fig. 3 to those of Fig. 7).

Boundaries selected with flux conservation
Since ∇.B = 0, the same amount of azimuthal magnetic flux should be present in the FR front and rear (e.g., Dasso et al. 2006). This implies that in the FR frame the same amount of flux crosses the y-z plane before and after the time of the closest spacecraft approach to the FR axis, then for a locally straight FR axis: However, the FR axis of observed MC is expected to be bend toward the Sun. Therefore, following field lines from the front to the rear part, the same amount of flux is located in a slightly longer region along the axis in the front region than at the rear. Nevertheless, the curvature radius of the axis, R c , is large, typically several AUs (e.g., Janvier et al. 2015), compared to the distance r to the FR axis (r 0.1 AU), then the correction in r/R c , is dominant near the border of the FR, but it is still small. Moreover, the magnetic torque balance is expected to distribute equally the twist along the FR, at least locally. Then, the hypothesis of local invariance by translation along the FR axis is expected to be a good approximation, at least away from the FR legs (see Owens et al. 2012).
Within the above approximations, the accumulated magnetic flux between t F and t R , in the FR frame and per unit length along the axial direction, is Since dF y /dz includes the observed velocity component V x this flux computation includes automatically the details of the local FR expansion. It is also valid for any FR cross-section shape. Next, for a given boundary time, either t F or t R , dF y (t F , t R )/dz = 0 defines the other FR boundary by imposing the flux balance (as it should be the case for a FR). In case of multiple solutions, this implies that more than one FR, or a FR and another structure, are present. Then, the FR, selected by the initial choice of boundary, is defined by minimizing |t F − t R | so selecting only one FR. Another equivalent formulation is to remark that ∇.B = 0 implies for a configuration invariant in z: where A(x, y) is the z component of vector potential. Integrating any field line implies that it is located on a surface A(x, y) = constant. This property relate any location in the FR front to its corresponding rear position. Along the spacecraft trajectory, so y = y p (defined in Eq. (4)), A(x, y p ) can be derived by integration of B y . Next, taking the origin of A at the front boundary this implies A(t, y p ) = −dF y (t F , t)/dz where the x coordinate has been replaced by time for marking the position. Then, Eq. (15) can also be expressed with the vector potential component A, as done, for example, in the Grad-Shrafranov equation for magnetostatic field. The above property of flux conservation was used to define the back region of MCs (the region which was part of the FR before reconnection at the front occurred, Dasso et al. 2006Dasso et al. , 2007 and, more generally, the amount of magnetic flux reconnected with the background field (Ruffenach et al. 2015). Below, we further propose that Eq. (15) can also be used in a broader context to define coherent sub-FRs within the core of the physical FR. The aim is to test the stability and to limit the bias of the FR axis determined by the MV variance.

Defining the flux rope core
We have shown in Sect. 4.2, for FRs with cylindrical symmetry, that the MV bias decreases as the azimuthal field component B θ or B θ /|B| (MV0 or MV1) has a more linear dependance with the distance to the FR axis. For non singular distributions of the electric current and cylindrically symmetric FRs, B θ vanishes on the axis. A Taylor expansion of B θ (r) implies that B θ (r) is linear for sufficiently small r values, so the bias of MV is expected to be lower when MV is applied to the core of the FR. However, restricting the core has the disadvantage of considering less data points, and for real FRs it implies including relatively more small structures and noise/fluctuations that can compete with the FR structure and produce an additional undesired bias. Then, we use Eq. (15) to define sub-regions corresponding to inner FRs in the simulated data. For the models considered, Eq. (15) is simply solved by setting ρ front = ρ rear = ρ b . When real observations are analyzed, getting the boundaries for this flux-balanced inner sub-FR requires to know its orientation (i.e., B y in the FR frame is needed to compute F y ). Thus, it is needed to use a method that feedbacks on the best estimations of the orientation while checking that the flux balance criterion is met. Next, we describe how to impose the flux balance with data. The following method was tested with the above models where the flux balance is rigorously defined with ρ front = ρ rear . The method has two main steps: fixing either the front or rear boundary and scanning the other boundary. More precisely, the method first tests if there is more azimuthal flux at the rear than at the front of the MC, then, if it is the case it defines the FR rear time. If no flux balance is achieved in the previous step, the end time is then used as the rear boundary and the method defines the FR front time when the azimuthal flux is conserved.
We next describe this method with more details. Let us suppose first that the front boundary is fixed (e.g., it can be defined by a sharp jump in the magnetic and plasma data, or it can be also located inside the MC to study the FR core). The MV method is recurrently applied to the data limited in between t F and the scanned values of t R , then the flux balance is computed from Eq. (15), in the MV frame (function of t R ). The lowest value of t R which satisfies the flux balance is selected as the FR rear boundary. If no flux balance is achieved, the same method is applied by fixing t R at a sharp jump in the magnetic and plasma data (at the latest time compatible with MC properties), or earlier on, again to study only the FR core. Then, t F is scanned. Finally, the largest value of t F which satisfies the flux balance is retained.
In summary, present method defines the FR, or its core, which is present at the spacecraft crossing by imposing the azimuthal flux balance. Since the MV eigenvectors are computed with flux balance, the bias on the deduced FR axis orientation is minimized. Finally, this method of scanning boundaries to impose flux balance can be implemented in other methods, such as fitting a flux rope model. The limitation is the computation time, as a non-linear least-square fit is required for each tested boundary. Still, the computation time can be limited by including in the method a more efficient search of the zero of a function than a simple scan (the simplest being bisection or Newton's root finding).

Results with selecting the flux rope core
As expected, the results of Fig. 8a show that the MV1 bias is significantly reduced as a smaller central part of the Lundquist's FR is considered up to the limit where the FR is not crossed (ρ b = p). This contrast with the eigenvalues being closer as ρ b decreases, further showing another example where the ratios of eigenvalues are not good indicators of the quality of the deduced FR axis. Finally, as a natural consequence of MV1 and FR frames becoming closer with a lower ρ b value, the axial field is also better recovered when ρ b is decreased (right panel of Fig. 8a).
The above trends with ρ b are also obtained with MV0 but still with a larger bias than with MV1. Next, similar results are obtained when the expansion is included with a magnitude as typically observed (Fig. 8b). However, this improved axis determination is less effective as the FR is flatter, so as b/a is larger (Fig. 8c). Of course in application to observations, perturbations of the FR field and the limitation to too few data points will constrain the decrease of ρ b so close to the impact parameter p as in Fig. 8.

Biases due to fluctuations
The effects of fluctuations can separate two groups: -a change of position for the boundary determined by the flux balance, -a change of the variances in the x,y,z directions (with fixed boundaries).
Let us call b y the fluctuation in the y direction. The flux balance writes: which is simpler to write in function of x: Let us fix t F . t R is changing by ∆t, corresponding to the size ∆x, to satisfy the flux balance with fluctuations. The change of flux due to the main field B y is about B y,a ∆x where B y,a is the amplitude of the field near the FR boundary. The change of flux due to a sinusoidal fluctuation of wavelength l is at most b y,a 2/π, where b y,a is the amplitude of the fluctuation, which corresponds to the contribution of half wavelength (the others cancel). Inserting these orders of magnitude in Eq. (18), and diving by the FR radius b to normalize, the magnitude of the boundary change ∆x satisfies: Both the last two fractions are typically << 1, say on the order of 0.1, so ∆x/b < 0.01, which will give a very small bias on the axis (Fig. 7a is for ∆x/b = 0.1, so the effect would be much smaller, likely a factor of ten smaller). In conclusion, fluctuations would change significantly the computed boundary, so the axis determined, unless they are both large in amplitude and with a large wavelength. As an example, to get the effect shown in Fig. 7a we need ∆x/b = 0.1, then with b y,a /B y,a ≈ l/b both ratio needs to be ≈ 0.4, and this is the most favorable case where half the wavelength perturbation is not compensated (odd number of half wavelength across the FR).
For the second effect, the variance of Eq. (1) is rewritten with fluctuations as: where the linear terms in b n are neglected because they contribute only by at most half wavelength (same as in previous paragraph). Then, the contributions of the fluctuations is mostly to increase the variances of all B components. If isotropic fluctuations are introduced they have no effect in the determination of the variance extremums so on the eigenvectors. This is not the case for fluctuations orthogonal to B (Alfvenic fluctuations are expected in a low β plasma). Still, could the eigenvectors be significantly rotated? Let suppose a rotation of x, y directions to x , y in the FR frame. Reducing b 2 n in the x direction would need a significant correlation of the perturbation in the x, y directions (see Appendix A). Moreover, this rotation would significantly increase (B n − B n ) 2 by mixing B y with B x counteracting the above hypothetical decrease of variance due to fluctuations. The same consideration can be done with a rotation from x, z directions to x , z , with the difference that it is less costly to mix x and z component for the variance as they have more similar behavior than compared with the y component.
In conclusion, fluctuations cannot change significantly the i angle. The effect is expected to be larger on λ angle but still weak unless the fluctuations have a long wavelength, comparable to the FR radius and are of large amplitude. This is very difficult to achieve as we already know that the effect of the expansion, which also can be thought as a kind of large-scale perturbation, is small.

Conclusion
An accurate determination of the flux rope (FR) axis direction from the interplanetary data of B is important for statistical studies determining the global axis shape, for the determination of the FR physical properties (e.g., the internal distribution of the magnetic twist, the amount of magnetic flux and helicity) as well as for comparing the MC global properties and its orientation itself to the ones determined from its solar source region. MV is one method to estimate the direction of the FR axis. However, as any present method, MV is known to have biases. Our main aim was to identify them and to correct them as much as possible.
The orientation of the axis is better defined by the location angle λ and the inclination i, as defined in Fig. 1, rather than the traditional latitude and longitude of the axis, usingẑ GSE as the polar angle. For example, the biases of the MV method are of different origins and magnitudes for λ and i while the bias origins are mixed for latitude and longitude. The angles λ and i are then important to identify, then to correct, the MV biases. For that, we test MV results with several methods and models.
The eigenvalues ratio provided from applying MV to the FR observations quantifies how different/similar are the variances in the three different directions. The eigenvalues should be sufficiently separated to identify accurately the eigenvectors found by MV. However, this eigenvalues ratio is not a proxy for quantifying the accuracy of the estimation of the real FR axis direction, in contrast to the claim of several previous publications.
MV was first applied to the simulated B time series computed with a Lundquist's model and its generalisation to an elliptical cross-section. This introduces an increasing bias on the axis determination as the simulated spacecraft crosses the FR further away from its axis. This bias is reduced when the MV is applied to B/B time series, confirming the previous results of Gulisano et al. (2007). By analyzing cylindrical models we identify the origin of this bias as due to the departure of a linear dependence of the azimuthal field components with the distance to the FR axis. Including expansion in the modeled FR, with an expansion rate as it is typically observed, does not have a significant effect on the orientation obtained from MV.
Another bias is due to the selection of the time interval used to apply MV. In MC observations there are frequently several sets of plausible boundaries depending on which magnetic and plasma data are used to define abrupt change of temporal behavior. We verify on modeled synthetic FR that the selection of the time interval boundaries is crucial to minimize the bias in the determined FR orientation. We next solve this difficulty by imposing that the same amount of azimuthal magnetic flux should be measured before and after the closest approach of the FR axis. The estimation of the flux balance is computed in the MV frame scanning alternately the rear and the front boundary. We find, as expected, that the MV bias for the estimation of the FR axis direction is minimized when the azimuthal flux balance is satisfied. This flux balance technique can be applied to other methods, such as FR model fitting. We also anticipate a decrease in the axis orientation bias since the FR models have an intrinsic flux balance, so that the selected part of the data should also satisfy this constraint.
The above constraint of flux conservation also allows us to explore a variable part of the FR core as input to MV. Indeed, the test made with cylindrical models shows that the bias in the axis orientation is lower when an inner part of the FR is selected. This effect is weaker for FRs with flatter cross sections. In applications to observations, the limitation to the FR core is further justified by the expected presence of larger perturbations at the FR periphery, due to interactions with the external ambient solar wind.
However in observations, the presence of B fluctuations implies that the analyzed time interval should be sufficiently large and does not affect the axis determination too much. Then, a compromise should be taken between too small and too large time intervals. This compromise is MC case dependent since it is linked to the intensity and distribution of B perturbations. This can be realized in studied MCs by analyzing the evolution of i and λ in function of the position of one boundary (the other being computed by flux balance).
In the present study we limit the application of the MV to MCs with approximately symmetric B(t) observed profiles. This corresponds to about 70% the MCs observed at 1 AU with an asymmetry parameter below 15% (Figure 3h of Lepping et al. 2018) or with a distortion parameter (DiP) of 0.5 ± 0.07 (Table 4 of Nieves-Chinchilla et al. 2018b). Next, on top of the expansion, which introduces only a weak asymmetry, an intrinsic asymmetry exist especially for the faster MCs ). This spatial intrinsic asymmetry introduces another bias in the axis determined by the MV. Its correction needs both the development of the MV method, and also the development of asym-metric models to test the amount of remaining bias. This will be the object of another study.
Still, the present results are important for applications to mostly slow MCs. For FR close to cylindrical symmetry and in expansion, MV applied to B/B time series combined with the flux conservation gives results where the remaining biases are mostly present on the location angle λ with a magnitude typically proportional to the impact parameter. It implies that, with a known axis shape, there is a bias on how far the spacecraft crossing occurred from the axis apex, see Janvier et al. (2015). The sign of the bias depends on the sign of the product between the impact parameter and of the magnetic helicity.
At the opposite, the inclination i is well recovered (with a bias less than 4 • for an impact parameter below 0.7). This implies that the direction of the FR axis projected orthogonally to the radial direction (from the Sun) is well determined. Then, the application of the MV method to B/B time series, called MV1, is recommended for studying the amount of global rotation of the FR during its transit from the Sun to the in situ observation position. where the summation is done on the N data points of the time series and B is the mean magnetic field.
Let us apply Eq. (A.1) to the FR magnetic components observed (or simulated) by a spacecraft, with B x (t), B y (t), B z (t) the components in the FR frame (defined in Sect. 2.1). This defines a series of N data points B k i (i = x, y, z). We call the variances of these field components Var i . We define also the correlations between two specific components, as Appendix A.1: Is the FR y-direction well recovered by MV?
Let us suppose we are in the FR frame, and let us consider a 2D rotation by an angle δ around a fixed x axis, so supposingx MV = x FR , in order to analyze if MV will mix the components B y with B z when this rotation is permitted. Thus, the new component B y in the rotated frame (y ) will be B y = B y cos δ − B z sin δ (A.3) Then, the variance in the new rotated y axis, Var y (δ), can be written as Var y (δ) = Var y (cos δ) 2 + Var z (sin δ) 2 − 2 sin δ cos δ Cor y,z = Var y − Var z 2 cos 2δ − Cor y,z sin 2δ + Var y + Var z 2 This is of the form a cos 2(δ−δ b )+b with a, b two constants with a > 0 if Var y > Var z . δ b is the bias angle of the MV: the rotation around the axisx MV =x FR needed to maximize Var y (δ). More precisely: a = Var y − Var z 2 2 + Cor 2 y,z tan 2δ b = − 2 Cor y,z Var y − Var z (A.4) A typical order of magnitude of δ b can be estimated with the field B y = B 0 sin t, B z = B 0 cos t with t in the interval [−π/2, π/2] (this keeps the global shape of the B y and B z profiles, in particular their symmetry, and it is sufficient for an order of magnitude estimation). This implies Var y ≈ B 2 0 /2, Var z ≈ B 2 0 /10 and Cor y,z = 0. Indeed, in practice we have often in MCs the ordering: |Cor y,z | << Var z << Var y (A.5) It implies δ b ≈ Cor y,z /Var y << 1 (A.6) A similar analysis can be done with x replacing z, so considering a rotation by an angle δ around the z axis, so supposingẑ MV =ẑ FR , and thus mixing B x with B y . As above, we have often in MCs the ordering: |Cor x,y | << Var x << Var y (A.7) This implies δ b ≈ −Cor x,y /Var y << 1 (A.8) which is even better satisfied than Eq. (A.6) in MCs when B x is lower than B z , so for a low impact parameter. We conclude that the y-direction, both orthogonal to the FR axis and spacecraft trajectory, is expected to be well recovered by MV when flux balance is achieved. This is summarized by a low expected bias for the i angle.
Appendix A.2: Is the FR axis direction well recovered by MV?
The same analysis done in the previous sub-section can be applied to the mix between B x with B z . Here it is clearer to look at the bias of the x axis, as a minimum of variance is expected nearby while a saddle point is expected around the z axis (variance decreases when changing in some directions, but increases in other directions). This is not a problem as the variance matrix is symmetric, then it has orthogonal eigenvectors, so when the bias on both x-and y-axis are known, the bias on the z-axis can be deduced.
Let us consider a rotation by an angle δ y around the y axis, so supposingŷ MV =ŷ FR . B x = B x cos δ y − B z sin δ y (A.9) Then, Var x (δ y ) is written as Var x (δ y ) = Var x (cos δ y ) 2 + Var z (sin δ y ) 2 − 2 sin δ y cos δ y Cor x,z = − Var z − Var x 2 cos 2δ y − Cor x,z sin 2δ y + Var x + Var z 2 just as in previous subsection with y changed to x and extracting a negative sign in front of cos 2δ y to write Var x (δ y ) as −a cos 2(δ y − δ y,b ) + b so that a minimum of Var x (δ y ) is present close to the x-axis for small δ y,b values (with a, b two constants with a > 0 if Var z > Var x ). δ y,b is the bias angle of the MV for this rotation: the rotation around the axisŷ MV =ŷ FR needs to minimize Var x (δ y ). More precisely: is typically much less satisfied in MCs, so the bias of the rotation δ y,b around the y-axis, as provided by the minimum of variance, could be important, and it increases with Cor x,z so with the impact parameter. Sincex MV andẑ MV are orthogonal this conclusion applies also to the axis direction estimated by MV. We conclude that the FR axis, or z, direction is expected to have more bias than the y direction, except when the impact parameter is small (which implies a small B x component, then a small Cor x,z ). This is summarized by a larger bias for the λ angle than for i except for small p values.