Global MHD simulations of Neptune's magnetosphere

A global magnetohydrodynamic (MHD) simulation has been performed in order to investigate the outer boundaries of Neptune's magnetosphere at the time of Voyager 2's flyby in 1989 and to better understand the dynamics of magnetospheres formed by highly inclined planetary dipoles. Using the MHD code Gorgon, we have implemented a precessing dipole to mimic Neptune's tilted magnetic field and rotation axes. By using the solar wind parameters measured by Voyager 2, the simulation is verified by finding good agreement with Voyager 2 magnetometer observations. Overall, there is a large‐scale reconfiguration of magnetic topology and plasma distribution. During the “pole‐on” magnetospheric configuration, there only exists one tail current sheet, contained between a rarefied lobe region which extends outward from the dayside cusp, and a lobe region attached to the nightside cusp. It is found that the tail current always closes to the magnetopause current system, rather than closing in on itself, as suggested by other models. The bow shock position and shape is found to be dependent on Neptune's daily rotation, with maximum standoff being during the pole‐on case. Reconnection is found on the magnetopause but is highly modulated by the interplanetary magnetic field (IMF) and time of day, turning “off” and “on” when the magnetic shear between the IMF and planetary fields is large enough. The simulation shows that the most likely location for reconnection to occur during Voyager 2's flyby was far from the spacecraft trajectory, which may explain the relative lack of associated signatures in the observations.


Introduction
Within the solar system, the planets Neptune and Uranus are examples of "ice giants," where an important component of their internal composition is volatile ices such as water, ammonia, and methane [e.g., Guillot, 2005;Arridge et al., 2012;. Both planets are both very important for a whole host of solar system, planetary, and heliophysics science goals and test our theoretical understanding of planetary formation [e.g., Tsiganis et al., 2005], dynamo physics [e.g., Soderlund et al., 2013], and magnetospheric dynamics [Bagenal, 1992]. Furthermore, the study of both Neptune and Uranus may provide insights into the behavior of exoplanets more generally, where a significant number are thought to be Neptune-like [Batalha et al., 2013].
The primary source of experimental observations concerning the ice giants comes from Voyager 2, which is the only spacecraft thus far to visit Uranus andNeptune in 1986 and1989, respectively Miner, 1986, 1989]. These flybys provided the bulk of the evidence that the ice giants are fundamentally different from Jupiter and Saturn Connerney, 1993;Mauk and Fox, 2010]. In addition to their composition, both exhibit magnetic moments that are significantly tilted when compared to the planetary rotation axis (∼60 ∘ at Uranus and ∼47 ∘ at Neptune). The rotation axes of the planets are also tilted relative to the normal of their orbital plane (∼97.5 ∘ at Uranus and ∼23 ∘ at Neptune) [Bagenal, 1992;Holme and Bloxham, 1996]. Both planets' magnetic fields also exhibit considerable higher-order moments compared to the other giant planets [Russell and Dougherty, 2010]. Such large magnetic moment tilts are thought to be generated by a complex dynamo within each planet which is not fully understood. It is thought both Uranus and Neptune are composed of a mixture of ices and rocks [Connerney, 1993;Soderlund et al., 2013], which causes the formation of complex, high order, nonaxisymmetric magnetic fields [Ness et al., 1989;Connerney et al., 1991;Holme and Bloxham, 1996].
Neptune's specific combination of a relatively large tilt in its rotation axis and a similarly tilted magnetic moment leads to the formation of a very dynamic magnetosphere which represents a unique challenge to our understanding of planetary magnetospheres. The Voyager 2 observations showed that Neptune rotates On approach to Neptune, Voyager 2's trajectory took it first through the bow shock, then exiting the magnetosheath into the magnetosphere through the dayside cusp region near the subsolar point [Ness et al., 1989;Szabo et al., 1991]. At this time the magnetic moment was serendipitously found to be approximately parallel to the solar wind flow. Voyager 2 remained inside the magnetosphere for more than 38 h and thus more than 2 Neptunian days. Compared to the other outer planets, Neptune's magnetosphere was found to be relatively empty [Belcher et al., 1989] although Triton, with its atmosphere [Broadfoot et al., 1989], may act as a source Zhang et al., 1991]. However, the role of Triton as a plasma source in Neptune's magnetosphere is ultimately not yet known .
Due to its distance ∼30 AU from the Sun, on approach to Neptune Voyager 2 found both the solar wind density and magnetic field strength to be low, with average values of 6.78 × 10 −3 cm −3 and 0.18 nT, respectively [Arridge, 2015]. The solar wind speed was largely similar to values observed in the inner heliosphere, v avg ∼ 461 km s −1 [Arridge, 2015]. This, in combination with Neptune's relatively large magnetic dipole moment (2.2 × 10 17 Tm −3 , 25 times greater than that of Earth), means that the magnetosphere is relatively large and that the rotation period cannot be ignored in the context of the flyby. Its rotation leads to two distinct configurations of the magnetosphere, which alternate every half rotation (∼8 h) [Belcher et al., 1989;Selesnick, 1990;Bagenal, 1992].
On the basis of the measurements made by Voyager 2, it was concluded that the magnetosphere goes back and forth from an Earth-like to a pole-on configuration, as shown in Figure 1. With the magnetic axis almost perpendicular to the solar wind flow, the Earth-like configuration is thought to resemble the Earth's magnetosphere, with the standard dayside plasma region, cusps, and tail containing a plasma sheet. In the pole-on configuration, the magnetic axis is almost parallel to the solar wind direction, causing a rarefied cusp region on the dayside, and northern and southern current sheets separated by a rarefied cusp region in the tail [Ness et al., 1989;Belcher et al., 1989;Selesnick, 1990;Lepping, 1994;Schulz et al., 1995]. Although they appear to be two separate current sheets in the 2-D slice shown in Figure 1, it is thought they may be connected, forming a cylindrical structure [Voigt and Ness, 1990;Schulz and McNab, 1996].
Following the Voyager 2 observations, several attempts have since been made to model Neptune's magnetosphere analytically and to predict aspects of its dynamics. For example, by using a simple model for the convection electric field [Selesnick and Richardson, 1986] and calculating the convection velocity inside the corotating frame, Selesnick [1990] found that the rotation could cause a particle accelerator effect, which may transport particles either toward or away from the planet.
Regarding the energy input into Neptune's magnetosphere, magnetopause reconnection is thought be significant [Desch et al., 1991], but highly dependent on season, time of Neptune day, and interplanetary magnetic field (IMF) orientation [Masters, 2015], and is expected to be "bursty" and infrequent [Huddleston et al., 1997]. Subsequent processes within the magnetosphere, such as tail reconnection, may suffer complications due to winding of magnetic field lines in a similar mechanism as proposed by Cowley [2013] for Uranus, who argued a seasonal dependence on the ability to close open flux in the Uranian system. Previous attempts at modeling Neptune's magnetosphere have been performed by Voigt and Ness [1990] and Schulz et al. [1995]. Voigt and Ness [1990] used a two-dimensional magnetohydrodynamic (MHD) equilibrium magnetosphere to explore the energy balance within the magnetosphere. They used an empirical model of the magnetosphere [Voigt, 1981] to verify the accuracy of their MHD model. By finding linear solutions to the Grad-Shavranov equation, they suggest the free energy within Neptune's magnetosphere remains constant throughout Neptune's daily rotation. Schulz et al. [1995] modeled Neptune's magnetic topography using the source surface model [Schulz and McNab, 1996] for different angles of attack (Ψ), i.e., the largest angle between the solar wind velocity and dipole moment. Both models [Voigt, 1981;Schulz et al., 1995] predict for large Ψ; the magnetotail undergoes a distinct change in topology: the tail current sheet no longer connects to the magnetopause current sheet but closes in on itself, forming a cylinder. However, neither model is intrinsically time dependent. This is shown more specifically in Schulz et al. [1995], who note that it is not fully self consistent as their model does not achieve thermodynamic equilibrium. Although there have been science cases put forward for new missions to visit both Uranus [Arridge et al., 2014] and Neptune [Agnor et al., 2009;Hansen et al., 2009;Christophe et al., 2012;, in situ observational data for the foreseeable future remain limited to the Voyager 2 data set. Computer simulation offers an alternative approach to experimental measurements and can be useful in providing more insight into Neptune's magnetosphere and complement existing models and observations. To better understand the general physics of magnetospheres formed by highly inclined, rapidly precessing dipoles and the magnetosphere of Neptune in particular, we have implemented a precessing dipole into the global MHD code Gorgon [Ciardi et al., 2007] and use this code to simulate the dynamics of Neptune's magnetosphere as observed during the Voyager 2 flyby. Our aim in this initial investigation is to understand the gross morphology of the magnetosphere of a highly inclined precessing dipole and specifically the change in shape and location of the outer magnetospheric boundaries on the dayside, the possible location and variation of magnetic reconnection on the magnetopause, and the changes in morphology of the magnetotail as a function of planetary rotation.
Since both the properties of Neptune's ionosphere and the impact of Triton on the Neptune's magnetosphere are very poorly understood and highly uncertain, and since the focus of our investigation is not on the inner magnetosphere, we have not included a Triton-like plasma source into the simulation, and the ionosphere is treated as a conducting shell. This therefore represents the initial step in modeling Neptune's magnetosphere but reveals important insight into its physical behavior. The simulation is verified by comparing the results with the observations taken by Voyager 2. We find that the Neptunian magnetosphere is highly dynamic, with many processes being inherently time dependent and regulated by the daily rotation.
The remainder of the paper is structured as follows. First there is a short overview of Gorgon and the parameters used in this simulation. Then the results are analyzed in the context of the magnetospheric configuration, dayside bow shock, magnetopause reconnection, and tail dynamics. Finally, the simulation is compared to magnetic field and plasma observations acquired during Voyager 2's flyby, followed by a discussion of the findings and concluding remarks.

Computational Approach
So-called global MHD codes are a common tool used to simulate the interaction of the solar wind with planetary magnetospheres. Most effort has been directed at reproducing the dynamics and behavior of the Earth's 10.1002/2015JA022272 magnetosphere [e.g., Lyon et al., 2004;Tóth et al., 2006;Raeder et al., 2008;Janhunen et al., 2012], but global MHD codes have also been used to study other planets, for example, the magnetospheres of Mercury [e.g., Kabin, 2000;Ip, 2002;Jia et al., 2015], Jupiter [e.g., Ogino et al., 1998;Moriguchi et al., 2008;Chanė et al., 2013], and Saturn [e.g., Fukazawa et al., 2007;Kidder et al., 2009;Jia et al., 2012]. The significant nonalignment of Neptune's and Uranus' dipole and rotational axes presents a different challenge to modeling their magnetospheres. In this regard, Uranus is more simple because the rotational axis is approximately contained within the planet's orbital plane. At the time of the Voyager 2 flyby of Uranus, the rotational axis was approximately antiparallel to the solar wind flow. The consequence of this is that if one changes into a frame rotating with the dipole, the solar wind flow is unaffected, but the IMF rotates around the x axis. This symmetry was exploited in the global simulations of Uranus' magnetosphere performed by Tóth et al. [2004]. In this work, they accounted for the rotation of the magnetic axis by utilizing the rotational symmetry of Uranus' magnetosphere at the time of Voyager 2's flyby, recasting the MHD equations to include rotational effects, effectively solving the MHD equations in a corotating frame. In the case of Neptune, this approach cannot be used because as discussed in section 1, the rotational axis does not lie in the Neptune orbital plane, being offset by ∼28 ∘ . The simplest approach to modeling Neptune's interaction with the solar wind requires a precessing dipole (i.e., whose axial orientation changes with time) to be included specifically in the code.

The Gorgon Code
To model Neptune's magnetosphere, we use Gorgon, a 3-D resistive MHD code originally developed to simulate laboratory plasmas. It has been used to simulate wire array Z-pinches [Chittenden et al., 2004;Jennings, 2006;Jennings et al., 2010], to model the physics of magnetic tower jets produced in astrophysical laboratory experiments [Ciardi et al., 2007] and laser-plasma interactions [Smith et al., 2007]. Gorgon solves the MHD equations (1) to (6) using a finite volume scheme on a 3-D uniform Eulerian Cartesian grid; here we model a simple fully ionized quasi-neutral H + plasma.
Equations (1) and (2) describe the mass continuity and momentum conservation equations, respectively, with the mass density = ( m p + m e ) n ≃ m p n, ⃗ v the fluid velocity, P p,e the proton/electron pressure, ⃗ J the current density, and ⃗ B the magnetic field.
The momentum equation (2) has two pressure terms (P p and P e ), for the proton and electron pressures. This is because the proton (equation (3)) and electron (equation (4)) energy equations are solved separately. The pressure is given by the ideal gas law, P p,e = nk B T p,e = 2 3 p,e , where represents internal energy density.
The next two terms in equation (4) are the Ohmic heating and the proton-electron energy equili- . The equilibration rate and the resistivity are based on the electron proton collision frequency calculated from the Spitzer model assuming an isotropic magnetization. Since the resistivity of the plasma is small (for T = 20 eV, ∼ 10 −5 Ω), the Ohmic heating term is negligible but has not explicitly been disabled. Reconnection arises through numerical resistivity The final term (Λ) is the optically thin radiation loss term due to bremsstrahlung. This is also a negligibly small component of the energy equation due to the conditions of a space plasma. A Von Neumann artificial viscosity is introduced in order to correctly capture the shock jump conditions and the deBar correction terms for Eulerian codes are introduced to improve the overall energy conservation [Benson, 1992].
The magnetic field evolution (equation (5)) is solved by adopting a vector potential representation. Using a staggered grid, with vertex centered ⃗ A field components and face centered ⃗ B field components on cubic cells, enables a differencing scheme which is conservative in ⃗ ∇ ⋅ ⃗ B. Therefore, no divergence cleaning is required to satisfy the magnetic divergence condition in this formulation.
Regions below a cutoff density (referred to as floor regions) of 6 × 10 −6 protons per cm 3 or 10 −26 kg/m 3 are treated numerically as vacuum such that the pressure, velocity, and current density are treated as zero and the vector potential solution resorts to the vacuum wave solution This expression can be solved explicitly using a Courant-Friedrichs-Lewy (CFL) condition time step of Δx 3c . The speed of light can however be relaxed to one tenth of its physical value without influencing the dynamical behavior of the plasma. Above the cutoff density, the hydrodynamics and ⃗ v × ⃗ B advection equations are solved using a second order van Leer scheme, with a variable time step set to 50% of the CFL condition for the largest magnetosonic speed. To increase the time step, a limit of 2 × 10 6 m/s is placed on the maximum speed of Alfvén waves that need to be resolved, with wave speeds above this damped using the method of Boris [1970].

Neptune Simulation Setup
We simulate a region −110 < x NSO < 40, −60 < y NSO < 60, −60 < z NSO < 60 (distances in Neptune radii, R N ). The simulation axes are equivalent to the Neptune Solar Orbital (NSO) coordinate system, with the x direction sunward, y opposite to the direction of Neptune's orbital motion, and z completes the right-handed set. To clarify the terminology used in this paper, we define north to be in the +z direction (i.e., northern regions are where z > 0) and south to be in the −z direction. Dusk is defined to be in the +y direction and dawn in the −y direction. Simulations were performed with a resolution of 0.3 R N and 0.5 R N . Both exhibited the same physical results and so only the high-resolution run is presented here.
The outer boundaries of simulation contain free-flow (Dirichlet) conditions, with exception of the +x boundary the solar wind, where a steady solar wind and IMF is set. The size of the tailward x boundary was also found to be sufficient to ensure it had no influence on the observed dynamics. Solar wind is input from the +x edge and allowed to propagate through the simulation. The solar wind conditions are kept steady throughout the simulation. They are based on the observations made by Voyager 2 just prior to crossing the bow shock: n P = 4.6 × 10 −3 cm −3 , v sw = 400 kms −1 , T = 0.5 eV, B x = 0 nT, B y = 0.07 nT, and B z = 0.11 nT [Szabo and Lepping, 1995]. It is worth noting that these solar wind conditions are approximately 75% lower (with the exception of the speed) than the average expected value [Slavin and Holzer, 1981] due to the variable nature of the solar wind, and hence, the size and shape of the magnetosphere is expected to differ from typical Neptunian magnetosphere. Furthermore, running the simulation with a steady solar wind does not capture such variability. In order to verify with Voyager 2 observations, the lower than average solar wind parameters are used.
Though Neptune's magnetic field is distinctly nondipolar near Neptune, the field becomes roughly dipolar at distances greater than 4 R N [Ness et al., 1989;Connerney et al., 1991]. Since here we are concerned only with the global structure of the magnetosphere and the behavior of the outer magnetospheric boundaries, we therefore approximate the field as a planet centered dipole, with strength 2.2 × 10 17 T m 3 (0.14 G R N 3 ) [Connerney et al., 1991]. This follows the approach of previous theoretical studies [Schulz and McNab, 1996;Voigt and Ness, 1990]. Although the Voyager 2 observations indicate that this dipole is offset by 0.5 R N from the center of the planet, centering the dipole simplifies the model but does not impact on the basic features of the simulated magnetosphere that are the focus of this study.
The inner boundary is spherical with a radius of 3.2 planetary radii, centered on the planet (since the grid is Cartesian, it is actually a "ragged" sphere). This region is designated a floor region, meaning all plasma density and velocity is numerically removed and so the surface acts as a mass sink. Inside this region, the vector potential inside a nested sphere of radius 1.86 planetary radii is forced to behave as a precessing magnetic dipole. The surface of this region of forced vector potential is therefore equivalent to a perfectly conducting sphere. In the floor region outside the inner boundary, the vector potential evolution is solved using the vacuum wave solution. The changing field of the dipole thus propagates rapidly via the wave solution before reaching the plasma, in which the fields are then updated through ⃗ v × ⃗ B and resistive diffusion. The simulation is initialized as an empty box (numerically set to vac ). The precessing dipole source begins at northern midnight: angled at colatitude (from +z axis) of 18.7 ∘ and azimuth (from +x axis) of 16.1 ∘ . It rotates about an axis angled at a colatitude of 28.3 ∘ and azimuth of 16.1 ∘ , which corresponds to the orientation of Neptune's rotation axis at the season Voyager 2 arrived. The solar wind is input from the start, filling the empty box and interacting with the precessing dipole to first form the magnetosheath and magnetopause. It runs for a total of 307,500 s (5.3 Neptune days). Figure 2 shows mass density slices in the y = 0 plane and z = 0 plane every half rotation from the 2.5 days into the simulation. By day 4 to day 5, the simulation has reached a quasiperiodic state and in the results presented below, we analyze output generated from this interval to avoid any phenomena associated with the initialization of the simulation.

Simulation Results
In this section we discuss the basic features and morphology of the magnetosphere. We discuss the overall configuration, the outer magnetospheric boundaries, magnetopause reconnection, and tail dynamics.

Magnetospheric Reconfiguration
We first discuss the overall configuration of the magnetosphere at different times of day. Figure 3a shows a cut of the plasma density in the y = 0 plane. Overlaid are magnetic field lines, whose color indicates the location in the y direction (blue y < 0, red y > 0). This corresponds to time step t = 290, 000 s, which is at the start of the fifth day. At this time, the angle between the magnetic dipole and the solar wind flow is the nearest to being perpendicular. The solar wind flows from left to right, and the solar wind magnetic field points in the [+y, +z] direction. The bow shock forms ahead of the magnetosphere, and the magnetosheath corresponds to the dark red region. Closed field lines on the dayside present an obstacle to the solar wind and two well-developed cusp regions form. These cusps are tilted asymmetrically, due to the tilt of the dipole. The plasma density in the closed dayside field region is lower than in the magnetosheath. A boundary layer of plasma persists antisunward of the cusps. Inside the magnetosphere, relatively empty plasma lobes are connected to the polar regions and map into the magnetotail as expected. The magnetotail lobes sandwich the plasma sheet, which largely corresponds to closed magnetic field regions on the nightside. Figure 3b shows the equivalent cut in the z = 0 plane at the same time. Again, the overlaid lines are magnetic field lines, but their color corresponds to their position in the z direction. The overall configuration of the magnetosphere therefore resembles that of the Earth. Figures 3c and 3d show the configuration of the magnetosphere half a planetary rotation earlier (time step 261,000 s) in the y = 0 and z = 0 planes, respectively. In this case, the dipole is almost parallel to the solar wind velocity. The bow shock and magnetosheath still form, though due to the dipole tilt, there is a cusp near the subsolar magnetopause, corresponding to Neptune's southern magnetic pole. North of this region, within the magnetosphere, there is a low density plasma region (compared to the sheath), which exists on closed magnetic field lines. South of this region, the magnetosphere is filled with magnetic field lines emanating from the southern magnetic pole which extend into the tail, acting as a barrier between the magnetosheath and the tail plasma sheet (the light red plasma region on closed field lines which extends into the tail). Finally, in the northern tail region, there is also a rarefied lobe plasma region corresponding to the magnetic field emanating from Neptune's northern magnetic pole, which also acts to contain the northern closed field plasma region. The dayside closed field region is limited to the vicinity of the planet by the magnetic field emanating from the northern magnetic pole and does not extend deep into the magnetotail. Figures 3c and 3d therefore show that the structure of the magnetotail for the pole-on configuration is somewhat different from that shown in Figure 1. Before we discuss this, we examine the boundaries of the magnetosphere. Szabo and Lepping [1995] used the inbound crossing to categorize the shock as "a low , high Mach number, and strong quasi-perpendicular shock," which was moving away from the planet at the time of the crossing. and Ness et al. [1989] bow shock models. The pole-on bow shock agrees well with the empirical models. It is much more blunt than the Earth-like bow shock, which follows the shape of the magnetopause.

Dayside Bow Shock and Magnetopause Profile
Using remote observation of the shock, Cairns et al. [1991] deduced the position and flaring of the shock is controlled by the rotation of Neptune's magnetic dipole. The crossing of the bow shock by Voyager 2 was measured both by its plasma experiment [Belcher et al., 1989] and its magnetic field experiment [Ness et al., 1989]. Two cylindrically symmetric parabolic models of the bow shock were obtained, which largely agree but have minor differences in standoff distance and flaring. They are derived taking the two crossings of the bow shock as observed by Voyager 2 and hence do not account for the rotation of the dipole nor the solar wind conditions, hence, only give a rough guide for Neptune's bow shock.
The location of the bow shock in the simulation can be obtained by locating the point where the plasma is compressed, indicating an increase in mass density. By scanning for this increase along the x direction for each value of y or z, 2-D slices of the bow shock position are obtained, as shown in Figure 4. The magnetopause location is identified using a method similar to Palmroth [2003], using momentum streamlines to identify the surface. Figure 4 shows the cuts of the shock on the dayside in the y = 0 and z = 0 planes for the Earth-like and pole-on configurations, as well as the two empirical models by Belcher et al. [1989] and Ness et al. [1989]. A good agreement is shown in Figure 4 between the empirical models and the pole-on configuration with regard to the location and shape, although the simulation predicts a blunter bow shock than the models. According to Richardson et al. [1994], Voyager 2 was in the magnetosheath for approximately 3.5 h; the observed bow shock location is therefore closer to the standoff distance of the pole-on bow shock.
Furthermore, the simulation reveals the bow shock is not cylindrically symmetric in either the pole-on or Earth-like configuration, being deflected in both the y = 0 and z = 0 planes. This asymmetric flaring could be due to the flaring of the magnetopause due to the IMF orientation, which is known to alter the shape in Earth's magnetopause [Lin et al., 2010]. The Earth-like configuration in Figure 4a has a cusp like dip at approximately z = 20 R N , which reflects the more complex magnetopause surface shape.
The Earth-like standoff distance is roughly 4 R N closer to Neptune than in the pole-on case and exhibits a similar amount of flaring. This result confirms that the shock is controlled by Neptune's precessing dipole, as explored by Cairns et al. [1991], but we find an opposite behavior with the standoff distance larger for the pole-on configuration because the magnetopause is more blunt.

Reconnection
Reconnection at Neptune has been studied before by Desch et al. [1991], Selesnick [1990], Huddleston et al. [1997], and Masters [2015], using either models or remote observations. Early studies focus on purely antiparallel reconnection [e.g., Selesnick, 1990] as the only form of reconnection: now it is well understood that reconnection can still occur for smaller magnetic shear angles [e.g., Paschmann et al., 2013;Eastwood et al., 2013]. The simulation shows evidence of reconnection occurring on the magnetopause and that it is modulated by the daily rotation (see below). Though reconnection is fundamentally a kinetic process, whose entire dynamics cannot be captured by MHD, it is generally thought that global MHD simulations with numerical resistivity accurately predict the location of reconnection sites [Komar et al., 2015].
To illustrate the reconnection process, Figure 5 shows snapshots in time of the magnetopause. Three-dimensional magnetic field lines and a slice of the velocity magnitude of the plasma are shown. The slice contains the IMF orientation (which is constant), and the solar wind flow direction, and centered through Neptune's origin (x = y = z = 0). It also shows the dipole moment direction. Reconnection is most likely to occur for antiparallel magnetic field configurations, which geometrically will be located in this plane if they arise. In fact, reconnection is found to be confined to the Northern Hemisphere as we now describe.
The eight panels show the progression from just before reconnection occurs in Figure 5a, to when it stops in Figure 5h, with 5000 s (= 0.086 Neptune days) between each figure. Reconnection "turns on" as the dipole moment rotates from the Earth-like to the pole-on configuration, approximately 0.2 Neptune day ( N ) before pole-on, shown in Figure 5c. There are high velocity jets shown in light red visible at the northern magnetopause surface, with kinked magnetic field lines. After that, in Figures 5d to 5f, reconnection continues to occur with the same high velocity jets and highly kinked magnetic field lines forming and moving downstream. These kinked field lines loop around themselves, indicating a flux rope type structure. In Figure 5g and later, approximately 0.2 N after pole-on, the field lines are smooth and no jets are visible, suggesting no reconnection is occurring in this plane. At this time, the local shear between the magnetospheric field and the IMF is 120 ∘ on the northern dayside. Hence, the reconnection appears to be controlled by the shear angle which, for the IMF condition observed, is largest for the configuration midway between the pole-on and Earth-like configurations. This gives a window of approximately 0.5 N (8 h) during Neptune's rotation where reconnection occurs. This reconnection site is in some ways similar to the case at Earth with lobe reconnection under northward IMF. However, the effect of daily rotation heavily regulating the reconnection may suggest that there is little or no global circulation.

Magnetotail
Neptune's magnetic topology has been modeled with the source surface model [Schulz and McNab, 1996] and empirical "stretch" model [Voigt and Ness, 1990]. It is predicted that near the pole-on case, for angles such that the solar wind velocity and dipole moments are parallel, the tail neutral sheet current is expected to close in on itself, rather than connect to the magnetopause current sheet.
By looking at the B x component in the tail, aspects of the magnetotail shape and configuration can be inferred. Figure 6 shows this at a distance of 40 R N downstream of Neptune, for the Earth-like and pole-on cases (see Figure 3). The location of the magnetopause in this plane, calculated with the method explained above, is also shown. In both configurations, there are two hemispheres of oppositely directed B x , separated by a current sheet, which is represented as B x = 0. In the top left and bottom right corners of the figures, the draped B x of the IMF can be seen. Overall, the tail is not cylindrically symmetric. In the Earth-like case, it is elongated in the −y,+z direction, suggesting that the IMF acts to skew the tail in its magnetic field direction.
The pole-on case shows a more symmetric tail magnetopause, which has been slightly shifted downward in the z direction.
We do not find a current system which closes in on itself in our simulation, as can be seen in Figure 6b. In both the Earth-like ( Figure 6a) and pole-on cases, the tail current sheet extends to each side of the magnetopause, connecting to the magnetopause current sheet. Figure 6b gives the appearance that the blue region is enclosed by the red region (top left), possibly indicating a lobe like region which could suggest the tail current sheet is almost separated from the magnetopause current sheet. However, this is a region where reconnection occurs, containing the highly kinked field lines which connect the IMF to the planetary field at the dayside cusp. This time step corresponds to Figure 5d, where the kinked field lines can be observed at x = −40 R N . These field lines do not connect to the dayside magnetic pole. In fact, they are connected to the nightside pole and are of the same topology as the blue region in Figure 6b. Furthermore, this is supported by the position of the magnetopause, which finds that the red region at top left is outside the magnetopause.

Comparison With Voyager 2 Flyby
Voyager 2 arrived at Neptune from roughly the ecliptic plane, passing near the subsolar bow shock and magnetopause. Neptune's rotation was such that Voyager 2 entered the magnetosphere through the cusp, in the near pole-on configuration. Neptune's gravity assist deflected the trajectory southward, where Voyager 2 crossed a plasma sheet, the magnetopause on the southern flank, and the bow shock again. As shown in Figure 7b by the first Voyager 2 position (point 1), the Voyager 2 spacecraft is very close to the magnetopause during the pole-on configuration, in keeping with the observation that Voyager 2 entered Neptune's magnetosphere through a cusp. As Voyager 2 arrives at its closest approach to Neptune, the dipole is now in its Earth-like configuration, as shown in Figure 7a at point 2. In the next Voyager 2 position, shown in Figure 7b, point 3, the dipole is back to pole-on, and Voyager 2 is positioned in the magnetotail. In subsequent positions (4-6), Figure 7 shows that Voyager 2 passes through the blue southern lobe region, indicating a very low mass density. The outbound (or southern) magnetopause crossing is not contained within the simulation domain. This can be seen in the current density plots in Figure 8, which shows the magnetopause as a peak in current density, indicating the magnetopause current sheet.
We have recreated Voyager 2's flyby by interpolating the magnetic field and plasma parameters along its trajectory, in accordance with its trajectory. Figures 9 and 10 show the inbound and outbound trajectories, respectively. For the outbound trajectory ( Figure 10), we only compare the magnetic field data because for large parts of Voyager 2's trajectory inside the magnetosphere, the ion flux was too low in order to calculate plasma parameters , except for in some isolated cases nearby the planet. Plasma observations are shown for the inbound trajectory, since there is complete data in the solar wind and magnetosheath, but no measurements were made after crossing the dayside magnetopause, except for an interval deep in the magnetosphere near to closest approach.
For the inbound trajectory in Figure 9 the simulated number density (n), bulk velocity (v), and the proton and electron temperatures (T p and T e ) are in good agreement with the Voyager 2 flyby data. No Voyager 2 ion measurements are available shortly after crossing the magnetopause, but cold dense plasma was detected closer to the planet. The simulation does not capture this, due to the fact it is initialized as an empty box, with no sources of mass inside the magnetosphere. The inbound magnetometer data agree well with observations too, capturing the bow shock, magnetopause, and the time varying magnetic field within the outer magnetosphere. Gorgon predicts a bow shock approximately 2 R N closer than that measured by Voyager 2, also seen in the bow shock profiles in Figure 4. The magnetopause is found to be in the same place, as seen in all components, but most noticeably the change in sign of the B y component.
There is also a good agreement in all components of the magnetic field in the outer magnetosphere on the outbound flyby, as shown in Figure 10. However, the B x is slightly smaller in this region than in the observations. The peak in B x and B z at around 25 August 17:00 is captured by the simulation and is due to the rotation of the magnetic dipole. Nearer the outbound magnetopause crossing 26 August 02:00, there is a difference in sign in the B y component between Gorgon and the Voyager 2 observation.
Gorgon predicts the flank magnetopause to be further out than Voyager 2, to the extent that it is not contained within the simulation domain. The change in sign in B y on 26 August 09:00 could have indicated a magnetopause. However, since no decrease in B x is observed, and the magnetic field lines drawn from this region are still connected to Neptune, it is still in the magnetosphere. This is also supported by our method of finding the magnetopause using momentum streamlines. For the Voyager 2 comparison, it is worth noting that the simulation is run with a steady solar wind, and hence, there is no observation on how the solar wind has varied once Voyager 2 has entered the magnetosphere, and this could account for discrepancies with the simulation and observed data. Figure 9. Voyager 2 observations compared to Gorgon. Starting in the solar wind, the Voyager 2 probe observes the bow shock at approximately 14:30 on 24 August, followed by the magnetopause at 18:00. Shortly after crossing the magnetopause, the density was too low in order for Voyager 2 to measure the plasma data. Gorgon data agree well, capturing the shock ratio well in the number density (n), velocity (v), ion temperature (T i ), and magnetic field (B) and predicting the position of the magnetopause. The location of the shock is closer to Neptune in the simulation, and Gorgon does not capture the high number density close to the planet.

Discussion
Overall, the simulation performs well and shows good agreement both the shock models and the Voyager 2 flyby for the large scale features. It shows that Neptune undergoes global reconfiguration during its daily rotation, with effects that are intrinsically time dependent in nature. These effects would be difficult to capture in a steady state or equilibrium model. For this reason, we find a difference in magnetotail magnetic topology compared to that previously thought and shown in Figure 1. The simulation shows that the magnetosphere is more asymmetric in the pole-on case, with the closed field line region north of Neptune remaining near to the planet, rather than stretching downstream to form a tail plasma sheet. Furthermore, the sunward magnetic cusp (corresponding to the south magnetic pole) creates a region of rarefied plasma between the southern magnetopause and the current sheet. The tail current sheet does not detach from the magnetopause to close in on itself during near pole-on configurations, as was expected from models of the magnetic topography [Schulz et al., 1995;Voigt and Ness, 1990]. Instead, the tail current sheet remains attached to the magnetopause current sheet. The difference can be due to the fact that neither of these models are time dependent and that the closed tail current sheet does not form due to the rotation of Neptune's magnetic field.
As predicted in previous studies, the window for reconnection to occur is found to be small, occurring intermittently on the northern dawnside region for the given IMF and season. The position of the reconnection site could suggest a reason for the lack of dynamics measured by Voyager 2 on its outbound trajectory, southward of the planet. The fact that reconnection occurs on the northern dawnside magnetopause could also be a factor in influencing the magnetotail by preventing the dayside closed field region extending deep into the nightside magnetosphere. Since we have only simulated the solar wind conditions as measured by Voyager 2, further studies could investigate whether how the magnetotail is affected by different solar wind conditions and different planetary season.
Though a good agreement is found with the Voyager 2 flyby on the inbound trajectory with respect to the outer boundaries, the simulation does not capture the plasma close to the planet. The simulation was initialized empty of plasma with no mass sources other than the solar wind; hence, we expect the plasma to be low in this region. Additionally, the outbound magnetopause is not captured. However, the solar wind is kept steady at the last known parameters that were observed. The arrival of a completely new and unmeasured solar wind front may explain any discrepancy in the location of the outbound magnetopause crossing.
As stated in section 1, the goal of this initial study is to understand the nature of the outer magnetospheric boundaries (including reconnection), the magnetotail morphology and the large-scale behavior of magnetospheres created by highly inclined precessing dipoles. Future work to explore the inner magnetosphere may benefit from a planetary magnetic field model containing higher-order moments and a separate investigation of Neptune's ionosphere and the possible role of Triton. Furthermore, simulations run under more typical and ideal solar wind and seasonal conditions could investigate the questions put forward by previous attempts at modeling Neptune's magnetosphere. However, this may be frustrated by the lack of experimental data.

Conclusions
We have used a global MHD simulation containing a precessing dipole to simulate the interaction between the solar wind and Neptune's magnetosphere for the first time. We used the solar wind conditions observed by Voyager 2 as a steady solar wind for the duration of the simulation. The simulation was verified by comparing the magnetic field components between simulation and Voyager 2, where a good agreement was found, leading us to conclude that within the scope of the Voyager 2 measurements, the simulation accurately reproduces the outer boundaries and magnetic environment.
We find that instead of a relatively symmetrical magnetosphere, there is a highly skewed distribution of magnetic field and plasma. During the pole-on configuration, the dayside plasma region trapped on the closed field lines remains near the planet, instead of being dragged out into a tail plasma sheet. The high field strength emanating from the pole creates a large rarefied plasma region, similar to the lobes in an Earth-like configuration, deflecting the plasma sheet and current sheet in the tail southward.
The sunward cusp extends its high field strength rarefied plasma region into the southern dayside, before extending into the tail, creating a highly skewed magnetopause and bow shock. By extracting the bow shock from the simulation, we find that the pole-on bow shock agrees well with the empirical models consistent with Voyager 2 observations. The bow shock has a larger standoff distance when the dipole is pole-on, than when it is Earth-like, differing by approximately 4 R N . Reconnection was found to occur in the northern region of Neptune's magnetosphere due to the season and IMF orientation. It is heavily modulated by Neptune's daily rotation and only emerges during the transition from pole-on to Earth-like.
Finally, previous models predict that Neptune's tail current sheet should close in on itself during the pole-on configuration, rather than connect to the magnetopause current system. We find this not to be the case for the solar wind measured during the Voyager 2 flyby, as the tail current sheet remains connected to the magnetopause current sheet. However, future work examining the structure of Neptune's magnetosphere for different seasons is required to fully understand the full range of possible configurations.