The Grad-Shafranov Reconstruction of Toroidal Magnetic Flux Ropes: Method Development and Benchmark Studies

TL;DR

This paper introduces a new method for reconstructing toroidal magnetic flux ropes in space plasmas using the Grad-Shafranov equation, validated through benchmark tests with known solutions and noise, with publicly available code.

Contribution

The paper develops a novel two-step numerical approach for toroidal GS reconstruction, including an optimal axis orientation and major radius determination, with validation and benchmarking.

Findings

Errors in axis orientation are about 9 degrees.

Major radius estimation error is approximately 22%.

Numerical solutions have less than 10% relative error.

Abstract

We develop an approach of Grad-Shafranov (GS) reconstruction for toroidal structures in space plasmas, based on in-situ spacecraft measurements. The underlying theory is the GS equation that describes two-dimensional magnetohydrostatic equilibrium as widely applied in fusion plasmas. The geometry is such that the arbitrary cross section of the torus has rotational symmetry about the rotation axis $Z$, with a major radius $r_0$. The magnetic field configuration is thus determined by a scalar flux function $\Psi$ and a functional $F$ that is a single-variable function of $\Psi$. The algorithm is implemented through a two-step approach: i) a trial-and-error process by minimizing the residue of the functional $F(\Psi)$ to determine an optimal $Z$ axis orientation, and ii) for the chosen $Z$, a $\chi^2$ minimization process resulting in the range of $r_0$. Benchmark studies of known analytic solutions to the toroidal GS equation with noise additions are presented to illustrate the two-step procedures and to demonstrate the performance of the numerical GS solver, separately. For the cases presented, the errors in $Z$ and $r_0$ are 9$^\circ$ and 22\%, respectively, and the relative percent error in the numerical GS solutions is less than 10\%. We also make public the computer codes for these implementations and benchmark studies.

Tables

Table 1. Table 1. : Comparison of the major geometrical parameters for the benchmark case. \ilabel tbl:para

Benchmark	$Z$, $[r_{sc},t,n]$	$r_0$ (AU)
Exact	[0.05076, 0.2538, 0.9659]	1.02
GS	[-0.09366, 0.3134, 0.9450]	0.80
Error	9∘	22%

Table 2. Table 2. : Comparison of the Outputs of the Numerical GS Solver with the Exact Solution ( 𝙽𝙻 = 0.0 ) {\tt{NL}}=0.0) for θ 0 = 0 subscript 𝜃 0 0 \theta_{0}=0 (first 3 rows), and θ 0 = 10 ∘ subscript 𝜃 0 superscript 10 \theta_{0}=10^{\circ} (last row). \ilabel tbl:solver

$\mathrm{𝙽𝙻}$	$R_f$	$B_{\varphi ,max}$ (nT)	$⟨E⟩$	${\Phi }_p$ (${10}^{12}$Wb/radian)
0.0	-	11.3	-	22.5
0.01	0.10	11.3	5.5%	22.4
0.1	0.28	10.5	9.5%	23.4
0.01	0.10	11.3	5.2%	23.1

Equations

$$R\frac{\partial }{\partial R}\left(\frac{1}{R}\frac{\partial \Psi }{\partial R}\right)+\frac{{\partial }^2\Psi }{\partial Z^2}=-{\mu }_0{R}^2\frac{dp}{d\Psi }-F\frac{dF}{d\Psi }.\text{ilabel}eq:GSt$$

$$\mathbf{B}=\frac{1}{R}\nabla \Psi \times {\hat{\mathbf{e}}}_{\varphi }+\frac{F(\Psi )}{R}{\hat{\mathbf{e}}}_{\varphi },\text{ilabel}eq:B$$

$$\frac{1}{r}\frac{\partial }{\partial r}\left(r\frac{\partial \Psi }{\partial r}\right)+\frac{1}{r^2}\frac{{\partial }^2\Psi }{\partial {\theta }^2}-\frac{1}{R}\left(\cos \theta \frac{\partial \Psi }{\partial r}-\frac{\sin \theta }{r}\frac{\partial \Psi }{\partial \theta }\right)=-{\mu }_0{R}^2\frac{dp}{d\Psi }-F\frac{dF}{d\Psi }.\text{ilabel}eq:GSrth$$

$${\Psi }_{sc}(r,\theta ={\theta }_0)={\int }_{r(1)}^{r}R{B}_{\theta }dr,\text{ilabel}eq:Psi0$$

$$Res=\frac{[{\sum }_i(F_i^{1st}-F_i^{2nd}{)}^2{]}^{\frac{1}{2}}}{|\Delta F|},\text{ilabel}eq:Rf$$

$${\chi }^2=\sum _{\nu =X,Y,Z}\sum _{i=1}^N\frac{(b_{\nu i}-B_{\nu i}{)}^2}{{\sigma }_{\nu i}^2}.\text{ilabel}eq:chi2$$

$$Q=1-chi2cdf({\chi }^2,\mathtt{dof}),\text{ilabel}\mathtt{eq}:\mathtt{Q}$$

$${\Psi }_{sc}^{\prime }={\int }_{r(1)}^r{R}^{\prime }{B}_{\theta }dr={\int }_{r(1)}^{r}R{B}_Z\cos {\theta }_{0}dr.$$

$$\Psi =\frac{C\gamma }{8}[(R^2-R_a^2{)}^2-R_b^4]+\frac{C}{2}[(1-\gamma )R^2]Z^2-\frac{1}{2}A{Z}^2,\text{ilabel}eq:PsiRZ$$

$$\Psi ={\Psi }_0\left[x^2-1+\frac{1-\gamma }{\gamma }\frac{1+{ϵ}^2}{{ϵ}^2}(1+\frac{2ϵ}{1+{ϵ}^2}x)y^2-\frac{1}{2}\frac{A}{{\Psi }_0/r_0^2}{y}^2\right].\text{ilabel}eq:Psixy$$

$$\tilde{\mathbf{B}}=\mathbf{B}+\mathtt{randn}()\ast \mathtt{NL}\ast ⟨|\mathbf{B}|⟩,\text{ilabel}\mathtt{eq}:\mathtt{dB}$$

$$E=\frac{|\psi -\Psi |}{⟨|\Psi |⟩}\times 100\%,\text{ilabel}eq:E$$

$$R=|{\mathbf{r}}_{sc}-\mathbf{O}|.\text{ilabel}eq:R$$

$$\frac{r^{\prime }-r_o}{Z_r}=\frac{t^{\prime }-t_o}{Z_t}=\frac{n^{\prime }-n_o}{Z_n}.$$

$$r_o=r^{\prime }-\frac{Z_r}{Z_t}(t^{\prime }-t_o)\text{ilabel}eq:ro$$

$$n_o=n^{\prime }-\frac{Z_n}{Z_t}(t^{\prime }-t_o).\text{ilabel}eq:no$$

$$t_o=\frac{({\mathbf{r}}_{sc}-{\mathbf{r}}_{op})\cdot \hat{Z}}{|Z{|}^2/Z_t}+t^{\prime },\text{ilabel}eq:to$$

$$u_i^{j+1}$$

$$v_i^{j+1}$$

$$a_i^j=rh{s}_i^j-{\left(\frac{{\partial }^{2}u}{\partial r^2}\right)}_i^j+\sin {\theta }^j{v}_i^j-\left(\frac{1}{r_i}-\frac{\cos {\theta }^j}{R_i}\right){\left(\frac{\partial u}{\partial r}\right)}_i^j.$$

$${\tilde{u}}_i^j=\frac{1}{3}[k_1{u}_{i+1}^j+k_2{u}_i^j+k_3{u}_{i-1}^j],$$

$$f_y=\min \left\{0.7,\frac{{\theta }^j-{\theta }_0}{{\theta }^{n_y}-{\theta }_0}\right\}.$$

Full text

\setlastpage\inarticletrue

{opening}

1 Introduction

\ilabel

s:intro Magnetic flux rope modeling based on in-situ spacecraft measurements plays a critical role in characterizing this type of magnetic and plasma structures. Simply put, it provides the most direct, definitive and quantitative evidence for the existence of such structures and their characteristic configuration of that of a magnetic flux rope with winding magnetic field lines embedded in space plasmas of largely magnetohydrostatic equilibrium. Such analysis dated back to early times of the space age, especially with the discovery of Magnetic Clouds (MCs) from in-situ solar wind data (see, e.g., Burlaga, 1995, and references therein). Among these modeling methods, employing in-situ magnetic field and plasma time-series data across such structures, the so-called Grad-Shafranov (GS) reconstruction method stands out as one (and the only one) truly two-dimensional (2D) method that derives the cross section of a flux rope in complete 2D configuration, or more precisely, $2\frac{1}{2}$D, with two transverse magnetic field components lying on the cross-sectional plane and the non-vanishing axial component perpendicular to the plane.

The conventional GS method applies to a flux rope configuration of translation symmetry, i.e., that of a straight cylinder with a fixed axis, but of an arbitrary 2D cross section perpendicular to it. Therefore the field lines are winding along such a central axis lying on distinct and nested flux surfaces defined by an usual flux function in 2D geometry. The GS method for a straight-cylinder geometry was first proposed by Sonnerup and Guo (1996), later further developed to its present form by Hau and Sonnerup (1999) and applied to magnetopause current sheet crossings (see also, Hu and Sonnerup, 2000, 2003). It was first applied to the flux rope structures in the solar wind by Hu and Sonnerup (2001), at first to the small-scale ones of durations $\sim$30 minutes, then to the large-scale MCs with detailed descriptions of the procedures tailored toward this type of GS reconstruction in Hu and Sonnerup (2002). Since then, the GS reconstruction method has been applied to the solar wind in-situ measurements of MCs by a number of research groups (e.g., Hu et al., 2016; Vemareddy et al., 2016; Wang et al., 2016; Hara et al., 2016; Sharma et al., 2013; Möstl et al., 2012; Kilpua et al., 2009; Möstl et al., 2009a, b, 2008; Liu et al., 2008; Du, Wang, and Hu, 2007). For a detailed review of the works related to GS reconstruction of magnetic flux rope structures, see Hu (2017).

Challenges facing the in-situ flux rope modeling including GS reconstruction stem from the variabilities in the configuration, properties and origins of magnetic flux ropes, concerning the MCs. For example, Kahler, Krucker, and Szabo (2011) examined a number of MC events at 1 AU, interpreted as magnetic flux ropes using relatively simple models of axi-symmetric cylindrical configuration. By comparing directly the modeled field-line lengths with the ones measured by traversing energetic electrons from the Sun to 1 AU (Larson et al., 1997), they concluded that the MC flux rope configuration, interpreted by the commonly known linear force-free model (Lundquist, 1950), is not consistent with such measurements. On the other hand, we showed in Hu, Qiu, and Krucker (2015) that for the same set of measurements, the field-line length estimates from the GS reconstruction results agree better with such measured path lengths from electron burst onset analysis. In addition to these unique measurements for the purpose of validating flux rope models, we also attempted indirect means by relating the in-situ GS flux rope model outputs with the corresponding solar source region properties. In an early work (Qiu et al., 2007), we established certain correlation between the magnetic flux contents and the corresponding flare reconnection flux on the Sun. Following that work, Hu et al. (2014) further extended the analysis to derive magnetic field-line twist distributions inside MCs based on GS reconstruction results, and hinted at the formation mechanism of flux ropes, at least partially, due to morphology in flares or magnetic reconnection sequence, thus leading to the variability in the twist distributions as observed from in-situ data. Capitalizing on these findings based on both in-situ flux rope modeling and observational analysis on the Sun, theoretical investigations (Priest, Longcope, and Janvier, 2016; Priest and Longcope, 2016) were also attempted very recently to probe the formation of flux ropes due to magnetic reconnection, as manifested by solar flares. Therefore, it is imperative to further develop the existing approaches of flux rope modeling to account for such variabilities in order to shed light on the important question regarding the origination and formation of magnetic flux ropes from the Sun.

In the present study, we intend to address the variability concerning the configuration of a magnetic flux rope, by extending the applicability of the GS reconstruction method to the geometry of a torus. We acknowledge that such an extension is not meant to be a replacement of the cylindrical flux-rope model, but an addition or an alternative to the toolset of flux rope modeling. The advantage of such a configuration over a straight-cylinder has to be assessed on a case-by-case basis. Sometimes it offers a useful and complementary alternative to the straight-cylinder model, especially when the latter model fails (see, e.g., Section \irefsubsec:solver).

A word of caution is that we only use a section of the torus to approximate the local structure of the flux rope in the vicinity of the spacecraft path across the toroidal section. Otherwise it would have implied that the flux ropes as detected in-situ possess a closed configuration with complete detachment from the Sun which has generally been refuted (e.g., Burlaga, 1995). However, a number of numerical simulations have utilized a closed magnetic configuration similar to that of a typical tokamak (or spheromak) to initiate CMEs close to the Sun (e.g., Shiota and Kataoka, 2016). In fusion sciences, the confined plasma experiments always have a closed geometry, e.g., a tokamak of axi-symmetric toroidal configuration (Freidberg, 1987). In this study, we try to tap into the wealth of knowledge in fusion plasma science describing 2D configurations in ideal magnetohydrodynamic (MHD) equilibria under such a geometry.

Somewhat as done before (see, e.g., Sonnerup and Teh, 2008, 2009), we adopt the practice of presenting basic theoretical consideration, analysis procedures and benchmark studies first in this presentation, but leave some more comprehensive benchmark studies and application to real events to a follow-up publication. This serves the purpose of not overwhelming the reader and ourselves, but guaranteeing a relatively short and focused report of the new development of this technique to benefit the user community.

The article is organized as follows. The GS equation in the toroidal geometry and the basic setup of the reconstruction frame are described in Section \irefsec:GSeq. Then a recipe in terms of a two-step reconstruction procedure is described in detail in Section \irefsec:proc. Benchmark studies of the basic procedures and the performance of the numerical GS solver are given in Section \irefsec:bench. We conclude in the last section, followed by several appendices laying out additional details and a special situation to be considered. We emphasize that the focus of this article is to allow interested readers to perform their own case studies and to devise their own computer codes if they choose to, facilitated by the detailed descriptions and the auxiliary material including the complete set of computer codes implemented in Matlab.

2 Grad-Shafranov Equation in Toroidal

Geometry

\ilabel

sec:GSeq

Equivalent to the GS equation in a Cartesian geometry on which the traditional GS reconstruction method is based, there is a GS equation in the so-called toroidal geometry of rotational symmetry, given in a usual cylindrical coordinate $(R,\phi,Z)$:

[TABLE]

As illustrated in Figure \ireffig:torcoord, the above GS equation describes the space plasma structure in quasi-static equilibrium of rotational symmetry, i.e., that of a torus. The configuration is fully characterized by a cross section of such a torus rotating around the rotation axis, $Z$, thus yielding invariance in the azimuthal $\phi$ direction, i.e., $\partial/\partial\phi\approx 0$. Under this geometry, the magnetic field vector is

[TABLE]

where the (poloidal) flux function $\Psi$ characterizes the transverse field components and has the unit of Wb/radian. The plasma pressure $p$ and the composite function $F=RB_{\phi}$, appearing in the right-hand side of equation (\irefeq:GSt), become functions of $\Psi$ only. Therefore similar to the straight-cylinder case, the 2D magnetic field components plus the out-of-plane one ($B_{\phi}$) are derived from the spacecraft measurements along its path across (along $-{\mathbf{r}}_{sc}$ in Figure \ireffig:torcoord) by solving the toroidal GS equation (\irefeq:GSt) over certain cross-sectional domain. In practice, the numerical GS solver is implemented for the GS equation written in the alternative $(r,\theta)$ coordinate (Freidberg, 1987):

[TABLE]

In this geometry, there are two main geometrical parameters to be determined, the orientation of the rotation axis $Z$, and the major radius $r_{0}$, whereas in the straight-cylinder case, only one parameter, namely the axis orientation of the cylinder, is to be determined. We note that the major radius can be either defined as the radial distance between the rotation axis and the geometrical center of the cross section of the torus or the distance to the location where the poloidal (transverse) magnetic field vanishes. We adopt the former in this study since it is the convention for plasma confinement studies (Freidberg, 1987). Inevitably, the parameter space is much enlarged in the present case and the reconstruction procedures are more evolved as to be described in the following section.

3 Procedures of Toroidal GS Reconstruction

\ilabel

sec:proc The procedures are presented for the most general cases of arbitrary orientation of the $Z$ axis and a relatively wide range of major radii of the torus. The analysis is primarily performed in the spacecraft or Sun centered $r_{sc}tn$ coordinate system (to distinguish from the local spherical coordinate $r,\theta,-\phi$; see Figure \ireffig:RSZ), where the radial direction is always along the Sun-spacecraft line, assuming a radially propagating solar wind carrying the structure.

We present a two-step recipe that is based on an extensive benchmark study of known analytic solutions. We stress that this is the best approach we have found so far, based on our experience and largely empirical studies. It is our intention to present what we have devised, deemed an optimal approach, and to deliver the reconstruction code to the user community for a timely release, for the purpose of much enhanced and collective effort in further validation and application of the toroidal GS reconstruction beyond the limitations of a solo effort. This is also the reason for our concise presentation of a “recipe” accompanied by the computer codes to enable others to either repeat the results or to generate their own.

3.1 The Most General Case

\ilabel

subsec:general

As illustrated in Figure \ireffig:RSZ, left panel, the most general case corresponds to a torus of arbitrary major radius and $Z$ axis orientation, whose central rotation axis intersects the $r_{sc}t$ plane at point $O^{\prime}$. Then relatively speaking, the spacecraft is moving along $-\mathbf{r}_{sc}$ across the torus, viewed in the frame of reference moving with the structure, usually the deHoffmann-Teller (HT) frame (taking the radial component only) that is well determined from the solar wind measurements (Hu and Sonnerup, 2002). The setup of such a local reconstruction frame $R\phi Z$ or $RSZ$ in Cartesian is shown in Figure \ireffig:RSZ, right panel, where the latter $R$ in $RSZ$ is fixed corresponding to the radial distance from $Z$ axis at the point of exit of the spacecraft from the torus. The spacecraft path along $-\mathbf{r}_{sc}$ with spatially distributed data points (via the usual transformation of a constant HT frame speed, $V_{HT}$) is rotated onto the light-shaded $RZ$ plane, where the spacecraft path is projected approximately onto the dimension $r$ of $\theta\approx\theta_{0}=Const$ in the alternative $(r,\theta,-\phi)$ coordinate. A cross section is obtained by solving the GS equation (\irefeq:GSrth) on the light-shaded plane, utilizing spacecraft measurements along $r$ at $\theta\approx\theta_{0}$, as spatial initial values, similar to the straight-cylinder case. However the distinction here is that due to the toroidal geometry, the projection onto the cross-sectional plane is not as straightforward as before. A rotation, rather than a simple direct projection, has to be performed. For brevity and completeness, we describe the details of determining both the origins $O$ and the corresponding radial distance array $R$ along the spacecraft path in the Appendix \irefapp:R, and also describe the details of the numerical GS solver in the local spherical polar coordinate $(r,\theta)$ in the Appendix \irefapp:solver.

In what follows, we describe, in details, the two-step procedures for determining the $Z$ axis orientation and its location in terms of its intersection with the $r_{sc}t$ plane, $O^{\prime}$, which, in turn, yields the size of the major radius of the torus. As before, this is implemented in a trial-and-error process with the location of $O^{\prime}$ distributed over a finite-size grid on the $r_{sc}t$ plane, each denoted by the pair $(\rho,\Theta)$, as shown in Figure \ireffig:RSZ, left panel. Then all possible $Z$ axis orientations are enumerated at each $O^{\prime}$ location. The current implementation is such that $\rho\in[0,1)$ AU of a uniform grid size 0.05 AU, and $\Theta\in[0,2\pi)$ of a uniform grid size $\pi/20$ for a spacecraft located at a radial distance 1 AU from the Sun, but excluding $\Theta=0$ and $\pi$ (see Section \irefsubsec:de). At each location $O^{\prime}$, a trial $Z$ axis of a unit vector is varied with its arrow tip running over a hemisphere of unit radius. Associated with each arrow tip, a residue (see equation \irefeq:Rf) is calculated based on the theoretical consideration of finding a functional $F$ that best satisfies the requirement of being a single-valued function of $\Psi$, based on the GS equation (\irefeq:GSt) (omitting plasma pressure for the time being; in other words, considering low $\beta$ plasma configuration only).

1. i)

The first step is to determine the $Z$ axis orientation via a minimization procedure of the residue defined in equation (\irefeq:Rf). This is done by a trial-and-error process as before, but over the finite-size grid on the $r_{sc}t$ plane. As shown in Figure \ireffig:RSZ (left panel), at each grid point $(\rho,\Theta)$, the trial unit $Z$ axis is varying over a hemisphere of unit radius. For each trial $Z$ axis, the local reconstruction frame is set up as shown in Figure \ireffig:RSZ, right panel, then the usual transformation from time-series data to spatially distributed data along the spacecraft path is performed, together with proper projection (rotation in the present case) to obtain data along the “projected” spacecraft path at $\theta\approx\theta_{0}$ on the light-shaded cross-sectional plane. Then the flux function along $r$ at $\theta=\theta_{0}$ is calculated

[TABLE]

implying $\Psi(r(1),\theta_{0})=0$. Conforming to the straight-cylinder case, a residue $Res$ is calculated following exactly the same definition as given in Hu et al. (2004) to quantitatively assess the satisfaction of the requirement that the functional $F$ be single-valued across the toroidal flux rope, i.e., quantifying the deviation between the $F$ values measured along the overlapping inbound (denoted “1st”) and outbound (“2nd”) branch along the spacecraft path, as represented by circles and stars in Figure \ireffig:RSZ (right panel; see Figure \ireffig:pta003 for an example), respectively:

[TABLE]

where the index $i$ runs through an abscissa spanning the range of $\Psi$ value of overlapping branches, and the normalization factor $\Delta F$ represents the corresponding range of the functional value $F$ over the two branches. Then the optimal $Z$ axis orientation is chosen as the direction of minimum $R_{f}$ within certain error bound, among the set of locations of $O^{\prime}$. 2. ii)

The second step is to re-run Step I with the chosen $Z$ axis orientation and a proper evaluation of $\chi^{2}$ with measurement uncertainties over the $(\rho,\Theta)$ grid. The quantity $\chi^{2}$ is defined according to Press et al. (2007) to evaluate the goodness-of-fit between the measured magnetic field $\mathbf{B}$ and the GS model output $\mathbf{b}$ along the spacecraft path, with given uncertainties (e.g., those available from NASA CDAWeb for Wind spacecraft measurements) $\sigma$:

[TABLE]

Often a reduced $\chi^{2}$ value is obtained by dividing the above by the degree-of-freedom ($\tt{dof}$) of the system. Since in producing $\mathbf{b}$, a polynomial fit of order $m$ (usually 2 or 3) is performed for $F(r,\theta=\theta_{0})$ versus $\Psi_{sc}(r,\theta=\theta_{0})$, it follows ${\tt{dof}}=3N-m-1$. Then the usage of this step is to yield a unique pair $(\rho_{min},\Theta_{min})$ at which the corresponding reduced $\chi^{2}$ value reaches minimum, $\chi^{2}_{min}$, for the $Z$ axis orientation determined in Step I. In addition, a quantity $Q$ indicating the probability of a value greater than the specific $\chi^{2}$ value is also obtained

[TABLE]

where the function $\tt{chi2cdf}$ is the cumulative distribution function of $\chi^{2}$ as implemented, for example, in Matlab. The associated uncertainty bounds can be assessed for various output based on the standard $\chi^{2}$ statistics (Press et al., 2007).

These are the two essential steps we develop to carry out the GS reconstruction in a general toroidal geometry that have been implemented in Matlab (the code is included in the auxiliary material accompanying this article). The additional details, such as the construction of the reference frame $RSZ$ illustrated in Figure \ireffig:RSZ, and the final step of computing the numerical solution of $\Psi$ over an annular region on the cross section of the torus, utilizing equation (\irefeq:GSrth), are given in the Appendices. In short, the coordinate system $RSZ$ as illustrated in Figure \ireffig:RSZ, right panel, is used to obtain the projection of $\mathbf{r}_{sc}$ onto $r$ at $\theta\approx\theta_{0}$. Afterwards, the working coordinate system is switched to $(r,\theta)$ in which both Steps I and II are carried out.

We also caution that the toroidal GS reconstruction we present here applies to the situation of a spacecraft path exiting into the “hole” of the torus, but not to a situation of a spacecraft path submerged within the torus, i.e., not crossing through into the “hole”. This particular case would yield a “projected” spacecraft path departing significantly from a single coordinate line $\theta=\theta_{0}$ which renders a numerical solution to the GS equation impossible. We discuss in Appendix \irefapp:hodo in more detail what the indications are in terms of the magnetic hodograms from in-situ spacecraft measurements for such paths.

3.2 A Degenerated Case

\ilabel

subsec:deBefore we proceed to present benchmark studies of GS reconstruction of general toroidal configurations, following the aforementioned steps, we single out one special case that needs special treatment. This is the case of the rotation axis $Z$ being along $\mathbf{r}_{sc}$, i.e., for $\Theta=0,\pi$, in Figure \ireffig:RSZ. In this case, a degeneracy occurs such that the residue remains the same for all trial axis lying on the plane spanned by $\mathbf{r}_{sc}$ and the true $Z$ axis.

Such degeneration can be understood as follows. As illustrated in Figure \ireffig:Zd, left panel, all calculations are simply carried out in the plane spanned by $\mathbf{r}_{sc}$ and the true rotation axis. Then consider two cases: one with a $Z$ being perpendicular to $\mathbf{r}_{sc}$ and the other with $Z^{\prime}$ arbitrarily chosen as shown. For the former, the composite functional is $F=RB_{\phi}$, and the flux function along the spacecraft path is, according to equation (\irefeq:Psi0), $\Psi_{sc}=\int_{r(1)}^{r}RB_{Z}dr$ ($r\equiv R$). For the other case, correspondingly, $F^{\prime}=R^{\prime}B_{\phi}=RB_{\phi}\cos\theta_{0}$, and

[TABLE]

Therefore, it results $F/F^{\prime}=\Psi_{sc}/\Psi^{\prime}_{sc}$, given that the field components $B_{\phi}$ and $B_{\theta}=B_{Z}$ remain the same, and the above integrals are always evaluated along $\mathbf{r}_{sc}=r\hat{\mathbf{r}}$, in both cases. Since both functional values $F$ and $\Psi_{sc}$ change with the $Z$ axis orientation in the same proportion, the residue of $F(\Psi)$ remains unchanged for any trial $Z$ axis in the plane. An example of such a residue map is shown in Figure \ireffig:Zd, right panel, where the residue remains the same for any $Z$ axis that is lying on the plane spanned by the true $Z$ axis (along $n$) and $\mathbf{r}_{sc}$. Note that this behavior does not change with added noise since the derivation shown above still applies no matter whether or not noise is added.

In practice, such a degenerated case presented above in Section \irefsubsec:de can either be run separately or simply excluded, considering that this case may be encompassed by the uncertainty regions of the most general cases discussed in Section \irefsubsec:general, as to be illustrated below in benchmark studies. Alternatively, since such degeneracy only affects Step I the most, one may still include these grid points along $\Theta=0$ and $\pi$ in Step II, once an optimal $Z$ axis has been determined.

4 Benchmark Studies

\ilabel

sec:bench The benchmark studies of the reconstruction procedures are carried out against a set of analytic solutions to the GS equation (\irefeq:GSt) that has been well studied in fusion plasmas. In particular, such solutions were given by Freidberg (1987) for 2D toroidal configurations (for additional details and variations, see Cerfon and Freidberg (2010)). We provide below such analytic formulas in terms of the flux function as a function of space in the $(R,\phi,Z)$ coordinates and associated parameters defining the overall geometry that forms the basis of analysis in this Section.

From Freidberg (1987) (Chapter 6, pp. 162-167), an exact solution to GS equation (\irefeq:GSt) exists for a special known functional form of the right-hand side, i.e., $FF^{\prime}=A=Const$ and $-\mu_{0}p^{\prime}=C=Const$, and can be written

[TABLE]

with $R_{a}^{2}=r_{0}^{2}(1+\epsilon^{2})$ and $R_{b}^{2}=2r_{0}^{2}\epsilon$ where the ratio between the minor and major radii of the torus is $\epsilon=a/r_{0}$. The geometry of the cross section of the torus is completely determined by the parameters $r_{0}$ and $a$, which define the center $R=r_{0}$ and the boundary $R=r_{0}\pm a$ of the cross section at $Z=0$. The other constant $\gamma=\frac{\kappa^{2}}{1+\kappa^{2}}$ is related to the plasma “elongation”, $\kappa$, in confinement devices, defined as the ratio between the area of the plasma cross section and $\pi a^{2}$ (Freidberg, 1987). Now we start to deviate from Freidberg (1987) referenced above since our purpose is to utilize the solution provided by equation (\irefeq:PsiRZ), but not to follow the subsequent analysis of the properties of such a solution.

By normalizing both spatial dimensions by $r_{0}$, i.e., $R=xr_{0}$ and $Z=yr_{0}$, we obtain

[TABLE]

We choose the following parameter values to obtain solutions that yield reasonable geometric dimensions and magnetic field magnitude consistent with in-situ MC observations at 1 AU: $\epsilon=0.1$, $\gamma=0.8$, and $\Psi_{0}/r_{0}^{2}=1$ nT, $A=-40$ or -10 nT. Then the transverse field components $B_{R}$ and $B_{Z}$ (equivalently, $B_{r}$ and $B_{\theta}$) are obtained from equation (\irefeq:B). The axial field $B_{\phi}$ is determined from $F^{2}=2A\Psi+B_{0}^{2}$, where the integration constant $B_{0}$ is arbitrarily chosen.

The time-series data for analysis are obtained by flying a virtue spacecraft through such a torus along a pre-set path, in the direction opposite to $\mathbf{r}_{sc}$ (for different perspectives, see Figures \ireffig:RSZ and \ireffig:psi332). Then the magnetic field vectors $\mathbf{B}$ extracted from the analytic solution described above along this path in $r_{sc}tn$ coordinate are further modified by adding normally distributed noise component-wise up to certain level characterized by the quantity $\tt{NL}$:

[TABLE]

where the random number generator $\tt{randn()}$ yields numbers from a normal distribution of zero mean and unit standard deviation. Therefore each magnetic field component in the time series for the following analysis carries a constant standard deviation in each case $\sigma=\tt{NL}*\langle|\mathbf{B}|\rangle$.

Note that in the following benchmark studies, we omit the pressure gradient in the right-hand side of the GS equation completely, although the exact solution we test against does include a finite pressure distribution ($C\neq 0$; otherwise the solution is trivial). This is based on the consideration that in real applications to mostly low $\beta$ flux rope structures in the solar wind, the plasma pressure is usually less important and carries relatively larger measurement uncertainties. So the current GS model outputs for the toroidal geometry, i.e., the determination of $Z$ and $r_{0}$, are primarily based on the magnetic field measurements. The measurements of plasma pressure, of course, will be included in applications to real events.

4.1 Determination of $Z$ and $r_{0}$

\ilabel

subsec:Z

In this section, we present one example of benchmark studies to show the results of determining the orientation and location of the rotation axis $Z$, i.e., in turn, the major radius $r_{0}$, following the steps outlined in Section \irefsec:proc. Although a number of additional benchmark studies was carried out, based on the analytic solution of equation (\irefeq:Psixy) with different configurations, i.e., different virtue spacecraft paths across, and different noise levels, it is not possible to study and present all cases in an exhaustive manner. Therefore we choose to present one case and are providing the computer codes in Matlab to encourage the interested users to repeat or generate new results, and to follow up with their own studies.

Figure \ireffig:psi332 shows the overall configuration of this benchmark case in the $RZ$ plane, on which the exact solution is shown within the rectangular domain. The “projected” spacecraft path is along the red line of an approximately constant $\theta\approx\theta_{0}=3.0^{\circ}$ formed with the horizontal line intersecting the cross section at $R=R_{0}$ (see also Figure \ireffig:RSZ). The exact $Z$ axis orientation and major radius of the torus are noted in the title of the figure. The synthetic time-series data for analysis are obtained along $-\mathbf{r}_{sc}$ from the analytic solution shown with additional noise according to equation (\irefeq:dB) for $\tt{NL}=0.025$ in this case. The resulting time series are shown in Figure \ireffig:Zr0, right panel, together with the GS model output to be further discussed.

We carried out the analysis following the two steps delineated in Section \irefsec:proc for the most general case, i.e., $Z$ not along $\mathbf{r}_{sc}$ and not parallel to $n$ either, in this case. From Step I, we calculated the residue at each $(\rho,\Theta)$ grid point based on equation (\irefeq:Rf) and found the minimum value, $\min(Res)=0.37$. The corresponding residue map at this particular location where the minimum value was obtained is shown in Figure \ireffig:Zr0, left panel. The distribution of residues on this residue map exhibits multiple local minima, in the form of a string of “islands”, each of value $\min(Res)+1$. Sometimes, they often merge and form one elongated shape enclosing a number of grid points. The general rule-of-thumb based on our experiments and experience is that the optimal $Z$ axis orientation should always be chosen near the middle of either one large single contour or one single “island” located near the middle of the group, as in the present case. Such an axis is chosen, usually through an interactive, manual process, as marked by the cross symbol, which is $[-0.09366,0.3134,0.9450]$ in $r_{sc}tn$ coordinate.

With this chosen $Z$ axis, we subsequently carried out Step II. The results of the reduced $\chi^{2}$ distribution and the corresponding $Q$ values are shown in Figure \ireffig:chi2a and b, respectively. Equation (\irefeq:chi2) can be used to evaluate the reduced $\chi^{2}$ values by replacing the variables $\mathbf{B}$ and $\mathbf{b}$ by the ones normalized by $\sqrt{\tt{dof}}$. As stated in Press et al. (2007), such defined reduced $\chi^{2}$ values tend to a distribution of mean 1 and standard deviation $\sqrt{2/\tt{dof}}$ (maximum $\sqrt{2}$). A value $\sim 1$ indicates a “moderately good” fit. Correspondingly, the probability of such a “good” fit, $Q$, has to be significant, e.g., $>0.1$. Therefore, in Figure \ireffig:chi2, two contours of levels 1 and $1+\sqrt{2}$ are shown for the $\chi^{2}$ distribution, and a number of contours are shown for $\log_{10}Q$, with the innermost one of value $Q=0.9$. Combined, the contours of values $\chi^{2}=1$ and $Q=0.9$ indicate the extent of uncertainty in the location of $Z$, i.e., the uncertainty in major radius. Both the exact location and one selected location of $Z$ where $\chi^{2}$ reaches minimum are enclosed by the innermost contours. The corresponding major radii for these two locations are 1.02 AU and 0.80 AU, respectively. The resulting GS model output $\mathbf{b}$ components (together with $\mathbf{B}$) along the spacecraft path for the chosen $Z$ axis orientation and location of minimum $\chi^{2}_{\min}=0.711$ ($Q=1$), are shown in Figure \ireffig:Zr0, right panel.

As summarized in Table \ireftbl:para, the two major geometrical parameters, namely, the rotation axis $Z$ and major radius $r_{0}$, were determined through the above procedures and are compared with the exact values of this benchmark case. The absolute error in the $Z$ axis orientation is 9∘ and that in $r_{0}$ is about 22%. The latter can be regarded as an uncertainty estimate in $r_{0}$, since the separation between the exact and selected $Z$ axis locations spans approximately the half-width of the maximum extent of the innermost contours in Figure \ireffig:chi2.

4.2 Accuracy of the GS Solver

\ilabel

subsec:solver

We present, separately in this section, the benchmark studies on the accuracy of the numerical GS solver with details given in the Appendix \irefapp:solver. The purpose is to test the implementation of the solver in the code, and to assess its performance in terms of error estimates under idealized condition of an exact set of $Z$ axis and $r_{0}$, independently from Section \irefsubsec:Z.

Two cases of two different $\tt{NL}$ values are considered, for a geometry of the spacecraft path parallel to $R$, i.e., $\theta_{0}=0$, so that a direct point-by-point comparison between the exact and numerical GS solutions can be made with minimal interpolation effect. Such an exact solution is shown in Figure \ireffig:psisc (left panel) where the solution is given on the grid in $RZ$ coordinate, while the right panel shows the corresponding numerically calculated flux function values along the spacecraft path for the two cases indicated by the legend. The time series for the two cases of different levels of noise added to the exact solution are shown in Figure \ireffig:B003 for (a) $\tt{NL}=0.01$, and (b) $\tt{NL}=0.1$, respectively. Case (b) is used as an extreme example to illustrate the effect of noise (see additional results below). We observe that a real event in terms of derived quantities is close to case (a) or somewhere in-between case (a) and (b).

This is demonstrated by the corresponding $F(\Psi)$ plots and the fitting residues $R_{f}$ (Hu et al., 2004) in Figure \ireffig:pta003 along the spacecraft path. Case (a) resembles what one gets from real data with a typical and relatively small fitting residue that is considered acceptable (usually when $R_{f}<0.20$), indicating reasonable satisfaction of the requirement that the functional $F(\Psi)$ be single-valued. On the other hand, in case (b), the data scattering is large and the fitting residue exceeds 0.20, indicating that the satisfaction of $F(\Psi)$ being single-valued is questionable. The fitting polynomials are of 2nd order in these cases, while the 1st-order polynomials yield similar results. In practice, such reconstruction results for case (b) with this metric value would have been rejected.

The numerical GS reconstruction results for the two cases are shown in Figure \ireffig:map003 (a) and (b), respectively, in the usual format. Compared with the exact solution in Figure \ireffig:psisc, there are clear distortions due to noise and numerical errors. The deviations seemingly increase with increasing noise levels. The maximum axial field is 11.3 nT and 10.5 nT, respectively, the location of which is also different from that of the exact solution. **The areas of the strongest $B_{\phi}$ seem to be distorted or shrunk compared with Figure \ireffig:psisc (left panel), due to the errors which directly affect the evaluation of $F(\Psi)$ in obtaining $B_{\phi}$. ** To further assess, quantitatively, the numerical errors, Figure \ireffig:psi003 shows the contour plots of the flux function, with both the exact solution $\Psi$ and the numerical solution $\psi$, overplotted on the same set of contour levels, for both cases. It becomes clearer that case (a) solution agrees better with the exact solution than case (b). The range of the $\psi$ values, representing the amount of poloidal flux $\Phi_{p}$, for both cases, is well recovered, as indicated by the colorbar. This agrees with Figure \ireffig:psisc, right panel, where the calculated flux functions along the spacecraft path for both cases, although case (b) exhibits slightly larger errors, agree with the exact values well. This indicates the effectiveness of low-pass filtering we carry out at the beginning of the analysis in processing the time-series data.

We also quantify the error by calculating the relative percent error between the exact and numerical solutions, defined as:

[TABLE]

after interpolating the numerical solution $\psi$ (obtained on a set of $(r,\theta)$ grid) onto the set of $RZ$ grid on which the exact solution is defined. The corresponding results are shown in Figure \ireffig:err003 (a) and (b), respectively, in terms of contour plots of $E$ at certain levels between 1% and 30%. The overall pattern is that surrounding the initial line, i.e., the spacecraft path at $Z=0$ in these cases, the errors are generally small, especially for case (a), mostly $<5\%$, to greater vertical extent. The errors increase with increasing distance away from the initial line and toward corners of the computational domain. In case (b), the performance of the solver in the lower half domain ($Z<0$) is comparable to that in case (a), although that in the upper half domain is much worse.

**We also supply the time-series data from case (a) to the standard straight-cylinder GS solver to check the effect of the toroidal geometry and the specific magnetic field profile in this case. The axial orientation is determined as $z=[-0.1710,0.9838,0.05440]$, in the $r_{sc}tn$ coordinate, primarily along $t$ (or $\phi$) direction, in this case. The corresponding field-line invariant $P_{t}=p+B_{z}^{2}/2\mu_{0}$ versus the flux function $A$ and the functional fitting is shown in Figure \ireffig:GS0a, yielding a fitting residue $R_{f}=0.12$ of acceptable quality. The reconstruction result, however, fails to yield a flux rope solution, as shown in Figure \ireffig:GS0b. It shows an X-line type geometry, rather than an O-line type, i.e., that of a two and a half dimensional magnetic flux rope (or island). This is due to the peculiar magnetic field profile in this case (see Figure \ireffig:B003a), where the magnetic field magnitude decreases significantly toward the center, down by about a half, resulting in such a configuration of an X-line with much weaker field strength in the middle. **

In summary, the various quantities derived from the toroidal GS solutions are given in Table \ireftbl:solver, whereas the straight-cylinder GS solver fails to yield the flux rope solution. As discussed above, case (b) generally exhibits more significant errors than case (a), not surprisingly, due to its higher level of noise, while case (a) yields fairly accurate results in this limited set of outputs. Overall the errors in these quantities are limited within 10%, with the case (b) outputs approaching the limit, which likely represents an extreme-case scenario.

**In addition, we also examine a case of $\theta_{0}=10^{\circ}$ for ${\tt{NL}}=0.01$, as one example of nonzero $\theta_{0}$, such that the spacecraft is crossing along a slanted path. Figure \ireffig:psi010 shows the comparison of exact and numerical solutions, and the corresponding error evaluation by the quantity $E$. The results are similar to the case of $\theta_{0}=0$ of the same noise level. Because the underlying numerical scheme is exactly the same as laid out in the Appendix \irefapp:solver, the computation is still limited within an annular region. The corresponding set of outputs is also listed in Table \ireftbl:solver (last row), for which the exact value of $\Phi_{p}$ is 23.2 TWb/radian due to a slightly different boundary. **

5 Conclusions and Discussion

\ilabel

sec:summ In conclusions, we have developed a practical approach for Grad-Shafranov (GS) reconstruction of magnetic flux ropes in toroidal geometry, i.e., that of ring-shaped structures of rotational symmetry. We devised a recipe to derive the unknown geometrical parameters, i.e., the orientation of the rotation axis $Z$ and the major radius of the torus $r_{0}$, from in-situ spacecraft data and the toroidal GS equation. The algorithm utilizes uncertainty estimates associated with the spacecraft measurements to carry out proper $\chi^{2}$ minimization of the deviation between the measured magnetic field components and GS model outputs. Benchmark studies with analytic solutions to the GS equation and added noise of known variances were carried out and are presented to illustrate the procedures and to show the performance of the numerical GS solver in the toroidal geometry. Although shown separately and still limited, the results indicate an absolute error of 9∘ in $Z$ axis orientation, and a relative error of about 22% for the major radius in one case, while the relative percent errors in numerical GS solutions are generally less than 10%. The straight-cylinder GS solver failed to yield the flux-rope solution for this particular case.

We also make the computer codes written in Matlab publicly available, accompanying this publication, which can also be downloaded from the shared Dropbox folder111https://www.dropbox.com/sh/wd5btkbldu5xvga/AABHQjCRRUH1NpEprmnKsccOa?dl=0. The codes can generate most of the results presented in the main text, and are also ready for applications to real events. **The included Readme file outlines the command-line execution of the codes in Matlab to generate the results presented here with little need to modify the codes. ** We encourage the potential users to run the codes and to communicate with the author on any issues that may arise.

We will present additional and more comprehensive benchmark studies in a follow-up presentation, together with examples of applications to real events (Hu, 2016). The limitation of the current study is somewhat idealized conditions including adding the artificial noise of normal distributions. The best approach to overcome this might be to perform a more complete benchmark study by utilizing the numerical simulation data, for example, that of Riley et al. (2004), where a toroidal flux rope was propagated to 1 AU with synthetic data taken along two separate spacecraft paths across the structure. Those data were utilized in assessing the cylindrical flux rope models, and will be re-examined by the current toroidal GS model. A more comprehensive benchmark study combining Sections \irefsubsec:Z and \irefsubsec:solver will be presented.

The current implementation relies on the availability of reliable estimate of measurement uncertainties, for example, associated with magnetic field, which were usually derived from the corresponding higher resolution data. The utilization of such uncertainty estimates in real events will be further investigated, especially by using multiple time-series data from multiple spacecraft across the same structure. As demonstrated in the benchmark studies here, the contour of reduced $\chi^{2}\approx 1$ outlines the extent of uncertainties in GS model output. A more complete assessment of such uncertainties associated with various output parameters of the GS reconstruction results will be carried out in the forthcoming study.

Appendix A Calculation of $R$ for a Given $Z$ at $O^{\prime}$

\ilabel

app:R We present one approach here the calculation of the array $R$ for each point denoted by a vector $\mathbf{r}_{sc}$ along the spacecraft path across the torus, for a given $Z$ axis of components $(Z_{r},Z_{t},Z_{n})$ at location $O^{\prime}$, as illustrated in Figure \ireffig:RSZ. This is the distance between the origin $O$, given by the vector $\mathbf{O}$ and $\mathbf{r}_{sc}$ (note all vectors are given in the ${r}_{sc}tn$ coordinate):

[TABLE]

Then the key step is to derive $\mathbf{O}$ for each $\mathbf{r}_{sc}$, realizing that it is changing along $Z$ except for $Z$ being perpendicular to $\mathbf{r}_{sc}$. It is trivial for the special case when all $O$s coincide with one point along $Z$ (becoming $O^{\prime}$ when $Z$ is perpendicular to the $r_{sc}t$ plane). So the following is for a general case and for $Z_{t}\neq 0$.

From the known fact that both $O$ and $O^{\prime}$, denoted by vector components $(r_{o},t_{o},n_{o})$ and $(r^{\prime},t^{\prime},n^{\prime})$, respectively, are along $Z$, it follows

[TABLE]

For $Z_{t}\neq 0$, we obtain

[TABLE]

and

[TABLE]

By substituting them into $(\mathbf{r}_{sc}-\mathbf{O})\cdot\hat{Z}=0$ and rearranging the terms, we obtain

[TABLE]

where quantities on the right-hand side are all known with $\mathbf{r}_{op}=(r^{\prime},t^{\prime},n^{\prime})$. Then the vector $\mathbf{O}$ is fully determined from equations (\irefeq:ro) and (\irefeq:no) above. So is the array of $R$ from equation (\irefeq:R) along the spacecraft path.

Similar set of formulas can be obtained for the cases of $Z_{r}\neq 0$ or $Z_{n}\neq 0$.

Appendix B The Numerical GS Solver

\ilabel

app:solver The numerical GS solver for the toroidal GS reconstruction is in direct analogy to the straight-cylinder case (see, e.g., Hau and Sonnerup, 1999), i.e., the approach by the Taylor expansion, utilizing the GS equation (\irefeq:GSrth) for evaluating the 2nd-order derivative in $\theta$.

To lay out the implementation of the numerical scheme in the code, we denote $u_{i}^{j}=\Psi$ and $v_{i}^{j}=B_{r}$, where the indices $i$ and $j$ represent uniform grids along dimensions $r$ and $\theta$, with grid sizes $h$ and $\Delta\theta$, respectively. It is set $\Delta\theta=0.01h$, and $\theta^{j}=(j-j_{0})\Delta\theta+\theta_{0}$ ($j=1:n_{y}$), where the index of the grid at $\theta=\theta_{0}$, i.e., along the projected spacecraft path, is denoted $j_{0}$. Changing $j_{0}$ will allow the spacecraft path where the initial data are derived to shift away from the center line of the computational domain. Then the solutions to the GS equation can be obtained through usual Taylor expansions in $\theta$ (truncated at the 2nd-order term with respect to $\Psi$), both upward and downward from the initial line ($\theta=\theta_{0}$). For example, for the upper half annular region $j\geq j_{0}$, noting the relations $\frac{\partial\Psi}{\partial\theta}=rRB_{r}$, $\frac{\partial\Psi}{\partial r}=RB_{\theta}$, and $R=R_{0}+r\cos\theta$, we obtain (further denoting $rhs=-FF^{\prime}$, as a known function of $u$ via the functional fitting $F(\Psi)$, e.g., see Figure \ireffig:pta003):

[TABLE]

where the term $a_{i}^{j}$ involves the 2nd-order derivative in $\theta$ and is evaluated via the GS equation,

[TABLE]

As usual, the partial derivatives in $r$ are evaluated by 2nd-order centered finite difference for inner grid points and one-sided finite difference for boundary points.

Also similar to the usual straight-cylinder case, smoothing of the solution at each step is necessary to suppress the growth of numerical error. The same scheme is applied as follows to inner grid points only (Hu, 2001; Hu and Sonnerup, 2002) and for the upper half domain ($j\geq j_{0}$):

[TABLE]

where the coefficients are $k_{1}=k_{3}=f_{y}$, and $k_{2}=3-2f_{y}$, with

[TABLE]

The same applies to $v$, and similarly to the lower half domain.

Appendix C The Hodograms for the Cases of Submerged Spacecraft

Paths

\ilabel

app:hodo

These are the cases that cannot be dealt with by the toroidal GS reconstruction technique developed here. These had been traditionally analyzed by a fitting method to fit the spacecraft measurements along its embedded path to a theoretical toroidal flux rope model (see., e.g., Marubashi et al., 2015). As we discussed earlier and demonstrate further below, the “projected” spacecraft path takes a peculiar shape and the measured magnetic field components possess certain features as indicated by the associated hodogram pairs obtained from the usual minimum variance analysis (Sonnerup and Scheible, 1998).

We again demonstrate these cases by utilizing the analytic solutions presented in Section \irefsec:bench. However here the spacecraft path is specially taken, not to exit into the “hole” of the torus, but to be along the green line in Figure \ireffig:torcoord. Two such cases are presented in Figure \ireffig:embedded: (a) the spacecraft path is perpendicular to $Z$ so that the “projected” path is double-folded onto itself, resulting in a situation where the spacecraft is entering and exiting the cross section along the same path but is only half-way through, and (b) the spacecraft path is traversing along a slanted path, resulting in a warped non-overlapping path across about half of the cross section. For both cases, the magnetic field components change in time and show clear features of symmetry or anti-symmetry, and possess significant radial components, persistently $\sim 10$ nT throughout the intervals. This is because that the spacecraft is nearly encountering the same set of field lines during its inbound and outbound passages, and of the up-down symmetry in these cases. These features are clearly demonstrated by the corresponding hodogram pairs shown in Figure \ireffig:hodos. Especially in Case (a), the $B_{1}$ versus $B_{2}$ hodogram exhibits a nearly closed loop while the other one is double-folded, due to completely folded path. Case (b) also displays significant rotation in $B_{1}$, about 180 degrees. It is worth noting that this type of pattern in Case (a) is rarely reported in in-situ magnetic field measurements, except for the case of Romashets and Vandas (2003) where a nearly 360 degree rotation in the magnetic field was seen in the MC interval. In other words, we caution that for this type of configuration of a glancing pass by a spacecraft through a torus, the magnetic field signatures as demonstrated here need to be considered for proper modeling of these configurations.

The current implementation of the numerical GS solver cannot solve for a solution over a significant portion of the cross section because the “projected” spacecraft path is no longer along a single constant coordinate dimension, i.e., that of $\theta\approx\theta_{0}=const$, across the whole cross-sectional domain. A word of caution is that when interpreting the measured time series in the $r_{sc}tn$ coordinate, they have to be taken along the actual spacecraft path $\mathbf{r}_{sc}$ shown in Figure \ireffig:torcoord, not the “projected” ones on the $RZ$ plane shown in Figure \ireffig:embedded. Another important observation from these preliminary analysis is that the field rotation is actually more significant as indicated by the hodogram pairs in these cases of “glancing” passage of the spacecraft, contrary to general perceptions one may have. Although this provides proof of merits of flux rope model fitting to in-situ spacecraft data under the toroidal geometry, we urge that such fitting better be done in the way of equation (\irefeq:chi2) with the mathematical rigor of proper uncertainty estimates for quantitative and more objective assessment of the goodness-of-fit.

Acknowledgments

QH acknowledges partial support from NASA grants NNX14AF41G, NNX12AH50G, and NRL contract N00173-14-1-G006 (funded by NASA LWS under ROSES NNH13ZDA001N). The author benefits greatly from decade-long collaboration with Prof. Jiong Qiu. The author also acknowledges illuminating discussions with the LWS FST team members on flux ropes, in particular, Drs. M. Linton, T. Nieves-Chinchilla, B. Wood, and the PSI group. The author is also grateful for a few site visits to NRL hosted by Dr. M. Linton.