Mapping of steady-state electric fields and convective drifts in geomagnetic fields – Part 2 : The IGRF

A method of mapping electric fields along geomagnetic field lines is applied to the IGRF (International Geomagnetic Reference Field) model. The method involves integrating additional sets of first order differential equations simultaneously with those for tracing a magnetic field line. These provide a measure of the rate of change of the separation of two magnetic field lines separated by an infinitesimal amount. From the results of the integration Faraday’s law is used to compute the electric field as a function of position along the field line. Examples of computations from a software package developed to implement the method are presented. This is expected to be of use in conjugate studies of magnetospheric phenomena such as SuperDARN (Super Dual Auroral Radar) observations of convection in conjugate hemispheres, or comparison of satellite electric field observations with fields measured in the ionosphere.


Introduction
In an accompanying paper (Walker and Sofko, 2016, hereafter Paper 1) a new method for mapping electric fields in geomagnetic field models has been described.The method uses an additional set of first order differential equations for w ⊥ , the perpendicular separation of two closely spaced field lines.These differential equations use the spatial gradient tensor of the geomagnetic field to compute how the separation of the field lines changes with distance measured along the field line.For magnetostatic fields the field lines are equipotentials.This means that Faraday's law can be used to deduce how a component of the electric field changes with distance along the field line.If two values of w ⊥ are calculated then the total electric field can be found.
Electric field mapping is important for conjugate studies of ionospheric and magnetospheric convection.In this paper we discuss its application to the IGRF (International Geomagnetic Reference Field) (Finlay et al, 2010;Finlay and Thébault, 2015), presenting expressions for the gradient tensor, and computing examples of electric field mapping.
Because we have limited this paper only to the Earth's internal field, the results will only be suitable for use at fairly low altitude or at lower latitudes, where the external fields are small.Future work will incorporate the external fields in the form of the Tsyganenko (1987Tsyganenko ( , 1995Tsyganenko ( , 1996) ) models.

Differential equations
The method of calculating the electric field mapping is described in detail in Paper 1.It involves integrating a set of nine simultaneous first order differential equations: Here x 1 , x 2 , x 3 are the Cartesian coordinates of a point on a field line, and s is the distance measured along the field line.
The separation vectors w (1) i and w (2) i represent the perpendicular separation of two closely spaced field lines, normalized with respect to their initial values.The geomagnetic field Published by Copernicus Publications on behalf of the European Geosciences Union.
B i is assumed to be modelled by a well-behaved function of x i whose spatial derivatives are known.The unit vector in the magnetic field direction is μi ≡ B i /B.The gradient tensor of μi is (4)

Coordinate systems
The IGRF is fixed in a frame rotating with the Earth.The coordinate system in which it is expressed is the XY Z system of geomagnetism, which is a spherical polar system in which X is the horizontal northward component of B, Y its eastward component, and Z its component inwards along r.
Computationally, the most straightforward way to integrate Eqs. ( 1), (2), and ( 3) is in a suitable Cartesian system.The obvious coordinate system to use is the geographic system GEO (Russell, 1971) in which z is aligned with the spin axis, x is in the equatorial plane through the Greenwich meridian, and y is eastward.These are the only two coordinate systems used in this paper.Only when the external magnetic field is introduced will other coordinate systems be introduced.In this section we provide some results necessary for the use of these coordinate systems.
In the XY Z system the gradient operator is (5) If we define dx = −r dθ (6) dy = r sin θ dφ (7) we can write and The transformation of B i and T ij from the orthogonal curvilinear XY Z system to the GEO system is the only required coordinate transformation that involves non-Cartesian tensors.We can avoid the need to invoke the formalism of general tensor notation by calculating these transformations explicitly.
The transformation of the operator ∂/∂x i (and hence the magnetic field) from spherical polar to geographic coordinates is where where Note that l ij is not a tensor but a 3 × 3 transformation matrix.
For transformations between Cartesian systems l ij is constant: for the curvilinear XY Z system it is not.Thus and the 3 × 3 × 3 transformation matrix S ij k ≡ ∂l ik /∂x j can be represented by its i = 1, 2, 3 components

The electric field
In Eqs.(1), (2), and (3) the initial magnitudes of w (1) i and w (2) i , are unity.So long as they are not collinear, they can be in any direction normal to the magnetic field at the initial point on the magnetic field line.For convenience we can take the initial values as those lying respectively in the magnetic meridian and in the longitudinal direction.Thus they and the magnetic field are initially mutually perpendicular.
During the integration they remain perpendicular to the magnetic field, but not necessarily to each other.
The electric field components are found using Faraday's law as described in paper 1.The more elementary discussion there means that if E (1) i and E (2) i are the covariant components of the electric field along w (1) i and w (2) i then the scalar products, E (1) i w (1) i and E (2) i w (2) i are constant along a field line.The covariant components of E can be found from this and converted to contravariant components using where θ is the angle between w (1) i and w (2) i .The total electric field can then be found from the parallelogram rule.
3 The IGRF and its spatial derivatives

Definition of the IGRF
The internal magnetic field of the earth (Finlay et al, 2010) can be described by a magnetic scalar potential, expressed as a spherical harmonic expansion: where where a is the radius of the Earth, r the radial distance from its centre, θ the co-latitude, φ the longitude measured eastwards, ξ = cos θ , and the P m n are the Schmidt quasinormalized associated Legendre polynomials of degree n and order m with normalization constants The International Geomagnetic Reference Field (IGRF) is characterized by the coefficients g m n and h m n that have been computed by determining the best fit to data from the extensive world wide network of geomagnetic observatories.For dates earlier than the year 2000 N = 10, and for dates thereafter N = 13.The coefficients g and h have units of nanotesla.We choose to measure all lengths in units of Earth radii, so that, in what follows, a = 1.Thus, in our calculations, the units of V are nanotesla × Earth-radii, the units of B are nanotesla and the units of the magnetic gradient tensor are nanotesla per Earth-radius.

Expressions for the magnetic field and gradient tensor
The magnetic field is given by and thus the components of B are The elements of the tensor ∂B i /∂x j in the geomagnetic system of coordinates are then We require B i and the tensor ∂B i /∂x j in the geographic system.These can be found using the transformations (Eq.11) and (Eq.15).The coding of these is easy in Python, using the Scipy Einstein summation function to sum on repeated suffices.

Examples of electric field mapping
The method of calculation is essentially the same as that described in Paper 1, except that the magnetic field and the gradient tensor are found from the IGRF rather than the Harris model.The integration is carried out in the GEO system.The initial values for the integration are the coordinates in the GEO system of the starting point on the field line, and the components of the unit vectors perpendicular to the magnetic field in the azimuthal and meridional directions.The normalization of the field line separations is in terms of the initial infinitesimal separation element, hence the initial values have magnitude unity.At each step of the integration the magnetic field and gradient tensor are computed in the XY Z system and transformed to the GEO system using the transformations of Sect.2.2.The associated Legendre polynomials and their first derivatives are obtained from the standard Scientific Python (scipy) package.Their second derivatives are found from the second order differential equation for the polynomials.After each step the field line coordinates and the normalized separation vector are recorded in the GEO system.At the conclusion of the integration procedure any initial electric field can be mapped to each point along the field line, as described in Sect.2.3.
We present here some examples of mapping in the IGRF.It should be borne in mind that, because the external field is not included, the techniques are only valid for mid to low invariant latitudes, or for low altitudes.

Mapping to the ionosphere
In order to map to the ionosphere, either from the conjugate ionosphere in the opposite hemisphere, or from a spacecraft position, after each step in the integration we evaluate the difference between the calculated height and the height at which to terminate the integration: When f (s) changes sign we exit from the integration process and use the last two values of f (s) as the initial values of a root finding process to find the value of the step that makes f (s) zero.The root finding process is regula falsi, the rule of false position (Press et al., 1989, Sect. 15.1).Although a first order process it generally converges after only a few steps.
To illustrate the procedure we select a starting point at Sanae Antarctica (71 • , 40 22 S, 2 • 50 26 W ).Although this site is at a relatively high geographic latitude, it is a plasmapause station with corrected geomagnetic latitude 60.69 • , based on IGRF 2005 (Cafarella et al., 2008).At times of low magnetic activity, the field lines are quite well represented by the IGRF.
Figure 1 shows the normalized field line separations plotted as functions of position along the field line.The "az- imuthal" direction means that the initial separation is measured in the direction w az = μ × r; the "meridional" direction means that the initial direction of the separation is in the direction w az × μ.It must be emphasized that these separations do not remain in these orientations as the integration progresses as would be the case with a dipole field.Nor do the two directions remain normal to one another.While the integration procedure ensures that they each remain normal to B, by the time the opposite hemisphere is reached the angle between the two directions is 80.4 • .The electric field can be found at points along the field line by defining its value at the starting point, by using Eqs.( 52) and (53) of Paper 1.As an example, let its component in the meridian direction be 40 mV m −1 and in the azimuthal direction 30 mV m −1 .Figure 2 shows a Mayavi2 (Ramachandran and Varaquaux, 2011) visualization of this mapping from Sanae Antarctica to the opposite hemisphere.The visualization software allows the image to be rotated as desired.Four different views are shown as the image is rotated.It can be seen how the field line reaches a conjugate point in the iono- sphere of the southern coast of Greenland.The field is large at the ionospheres and becomes much smaller in the equatorial plane.Visualizations such as this do not greatly improve our understanding when the field is approximately dipolar: it will be much more useful when the external fields have been incorporated in the model.The components and total electric field as a function of distance measured along the magnetic field line are shown in Fig. 3.This allows us to gain a more quantitative idea of the variation along the magnetic field line.Note, in particular, how much stronger the electric field is off Greenland compared with that at Sanae.This a consequence of the fact the Sanae is sufficiently close to the South Atlantic magnetic anomaly for the magnetic field to be significantly smaller than the dipole field for the magnetic latitude.
It should be noted that, while the electric field components are inversely proportional to the field line separation, the magnetic field is inversely proportional to the cross sectional area of the flux tube, which is of the order of the square of the field line separation.The convection velocity, with magnitude E/B therefore scales as cross-sectional area divided by field line separation: the convection velocity E/B increases with decreasing magnetic field.
Another representation is shown in Fig. 4.Here the projections of the electric field on the three coordinate planes of the GEO system are shown.This gives a good idea of how the components of the electric field vary along the field line.The actual values of the electric fields in the GEO system at the beginning and end points for this particular case are

Mapping to the equatorial surface
We define the intersection with the magnetic equatorial surface to be the point on the field line where r reaches its maximum value, i.e.where dr/ds = 0. Then and dx/dr, dy/dr, dz/dr are given by Eq. ( 1).We thus have the conditions which we could have stated directly.The left hand side of this equation is the function of s whose root must be found in the convergence process.The result of this for the same point of origin as the previous example is that, at the equator, This field is compared with the initial field in the right hand panel of Fig. 5.We need not show any graphical information for this case since it is the same as the initial half of the trace for the previous case.

Discussion and conclusions
The purpose of this paper has been to implement the method of electric field mapping described by Walker and Sofko (2016) for the IGRF.Expressions for the gradient tensor of the geomagnetic field for this case have been derived.From these a set of first order differential equations for the perpendicular distance between two magnetic field lines separated by an infinitesimal amount can be found.This distance is normalized in terms of the value at the starting point to give a separation vector.These equations can be evaluated numerically simultaneously with the equations for tracing a magnetic field line.The component of the electric field in the direction of the separation vector can be mapped along the field line by applying Faraday's Law, noting that the electric field is always normal to the magnetic field.If this is done for two different initial separation vectors then the total electric field can be found.
The results of calculations show that the method is viable and generally useful.Currently the limitations are that it can only be used where external magnetic fields are small enough to be neglected, that is not between points on field lines that extend into a region where the external fields are important.This limitation will be removed by applying the method to the Tsyganenko field model.
Of course, the method only applies to electrostatic fields.In practice this means that scale of time variation of the electric field must be long compared with the Alfvén transit time along the field line.So long as we are only considering conditions where the IGRF is a good representation of the total magnetic field this is reasonable.Extending the method to greater radii and more active conditions will require implementation using the Tsyganenko models.In such cases it will be necessary to proceed with caution, ensuring that the electric fields induced by the changing magnetic fields are small compared with the electrostatic field.
A software package is being developed for general use and will be released under an open source licence.

Figure 2 .
Figure 2. Different aspects of the geomagnetic field line originating from Sanae Antarctica (Geographic coordinates −71.67 • , −2.84 • ) and electric field vectors (shown in red).Initial electric field has meridional component 40 mV m −1 and azimuthal component −30 mV m −1 .Equator and Greenwich meridian shown in red, 180 • meridian shown in green.Lines of latitude and longitude are separated by 15 • .

Figure 3 .
Figure 3. Electric field as a function of distance measured along the field line.

E 0 =Figure 4 .Figure 5 .
Figure 4. Projection of electric field on the coordinate planes for the case shown in Fig. 2.