Differential rotation in neutron stars with open and closed magnetic topologies

ABSTRACT

1 INTRODUCTION

Differential rotation between the rigid crust and the multiple fluid components of a neutron star has been proposed as one possible driver of rotational glitches (Warszawski & Melatos 2011; Andersson et al. 2012; Haskell & Melatos 2015). Among other factors, the amount of differential rotation is limited by the magnetic field threading the stellar interior. In general, hydromagnetic forces couple the charged components of the multi-fluid interior to the crust, which is decelerated by electromagnetic braking. Eventually, if the coupling is rapid, one expects the crust and the charged components to approach corotation, although they may rotate differentially with respect to the superfluid neutrons. The global effectiveness of the coupling depends on the magnetic topology. It has been argued that open topologies cause the charged fluid components to corotate with the decelerating crust on the Alfvén time-scale, so that all the magnetically coupled parts of the system spin-down together (Easson 1979; Melatos 2012; Glampedakis & Lasky 2015). In contrast, magnetic flux tubes that close inside the star decouple from the regions threaded by open field lines, allowing differential rotation to persist up to long, viscous time-scales in parts of the star.

Numerical and experimental studies of a magnetized, conducting fluid rotating differentially in a spherical Couette configuration have an extensive history. Quantitative results have been obtained concerning the role of the driving shear and shell thickness (Dormy, Cardin & Jault 1998; Nakabayashi, Zheng & Tsuchida 2002b; Nakabayashi & Tsuchida 2005), topology of the applied field (Dormy et al. 1998; Hollerbach, Daix Canet & Fournier 2007; Schmitt et al. 2008), conducting boundary conditions (Dormy et al. 1998; Hollerbach et al. 2007; Soward & Dormy 2010), superrotating shear layers (Dormy, Jault & Soward 2002; Nataf et al. 2006; Soward & Dormy 2010), non-axisymmetric instabilities (Hollerbach, Junk & Egbers 2006; Hollerbach 2009; Gissinger, Ji & Goodman 2011), and magnetorotational instabilities (Gissinger et al. 2011). In all the above studies, the magnetic field induced by the shear flow is small compared to the applied magnetostatic field. This approximation is suitable in various geophysical and laboratory applications.

Previous studies of how the magnetic topology affects differential rotation in the neutron-star context are based on analytic calculations (Easson 1979; Melatos 2012; Goglichidze & Barsukov 2019). In this paper, we approach the problem numerically for the first time without limiting the analysis to small induced magnetic fields. Using the finite-difference, Godunov-type, ideal-magnetohydrodynamics (ideal-MHD) solver pluto (Mignone et al. 2007), we perform a structured sequence of numerical experiments to measure and compare the angular velocity shear sustained by open and closed magnetic topologies over the short and long term in a differentially rotating spherical shell. The system represents crudely the rigid crust and multi-fluid outer core of a neutron star. We emphasize, however, that the model is highly idealized: its parameters are unrealistic astrophysically due to computational limitations, the boundary conditions do not capture the full complexity of the interaction between different stellar layers, and the single-fluid MHD equations of motion approximate the full, multi-component physics (Glampedakis, Andersson & Samuelsson 2010). A similar, idealized model, but without a magnetic field, has been used in the past successfully to study neutron star turbulence (Peralta et al. 2005, 2006; Peralta & Melatos 2009), superfluid spherical Couette flow (Peralta et al. 2009), pulsar glitch statistics (Melatos & Peralta 2007), and pulsar glitch recovery (Howitt, Haskell & Melatos 2016).

The paper is organized as follows. In Section 2, we formulate the problem in an MHD context and follow previous authors in emphasizing the important role played by magnetic-field-aligned coordinates and field line integrals when interpreting the numerical results to follow (Easson 1979; Melatos 2012; Glampedakis & Lasky 2015). Section 3 outlines the numerical method, the initial and boundary conditions, and the dimensionless parameters in the problem. The latter are ordered the same as in a neutron star, although their dynamic range is much smaller due to computational limitations; for example, the magnetic coupling time-scale is kept shorter than the viscous time-scale, as expected in neutron stars. Section 4 validates the code in the unmagnetized case and sets a baseline against which to compare the magnetized flow. Sections 5 and 6 calculate the magnetized flow in several closed and open geometries as a function of the magnetic field strength and the driving shear. The degree of differential rotation is quantified in each scenario. In Section 7, we relate the degree of differential rotation to field-line-averaged tension and use the same diagnostic to study the geometric evolution of the magnetized flow.

2 EQUATIONS OF MOTION

We model the stellar interior in an idealized fashion as a spherical Couette system, i.e. a differentially rotating, conducting spherical shell containing an incompressible viscous fluid obeying the equations of ideal MHD. The spherical Couette geometry and boundary conditions are transplanted from previous studies of the global flow pattern in a neutron star (Peralta et al. 2005, 2006; Howitt et al. 2016). The simulated system is related but not identical to that considered by Easson (1979), Melatos (2012) and Glampedakis & Lasky (2015). The main difference is that we neglect the presence of a superfluid, neutral component and consider a dynamical magnetic field that evolves non-linearly, in contrast to previous perturbative treatments (Easson 1979; Glampedakis & Lasky 2015).

2.1 Magnetic coordinates

Certain key invariants associated with the presence or absence of differential rotation involve integrals of MHD variables along open or closed magnetic field lines, as specified in Section 2.2. The line integrals are easier to calculate and interpret if we switch to magnetic coordinates. The use of magnetic coordinates is widespread in the analysis of complicated magnetic topologies in tokamaks (Lifshits 1989).

2.2 Field line integrals

3 NUMERICAL METHOD

3.1 PLUTO

The pluto solver (Mignone et al. 2007) integrates a closed set of MHD conservation laws using a finite-volume formalism. Volume-averaged conserved quantities are converted into primitive variables. Flux differences at cell interfaces are computed by solving the associated Riemann problem. We adopt a static grid stretched radially close to the inner and outer boundaries (to resolve viscous boundary layers). The typical grid resolution ranges between 252 and 300 points in radius, 206 points in latitude, and up to 20 points in azimuth ϕ. (The initial and boundary conditions enforce axisymmetry. The reader is referred to Section 7.3 for some preliminary tests of non-axisymmetric configurations.) We select the pluto option of a Lax–Friederichs–Riemann solver to compute the flux differences and a second-order Runge–Kutta algorithm to advance the solution in time. The time-step is controlled by the Courant–Friederichs–Lewy condition. The typical step is of order O(10⁻³) in dimensionless units, spanning the range 1.4 × 10⁻³ ≲ Δt ≲ 3.2 × 10⁻³. The maximum adaptive time-step adjustment is fixed to 1.1. We verify that the flow is resolved through a sample of higher-resolution tests.

3.2 Initial and boundary conditions

We start the simulations with |$\boldsymbol{v} = 0$|⁠, such that incompressibility |$(\nabla \cdot \boldsymbol{v} = 0)$| is satisfied at the outset and henceforth. The two spheres at r = R_i (inner sphere) and R_o (outer sphere) rotate with angular velocities Ω_i > Ω_o and Ω_o, respectively. We define the rotational shear ΔΩ = (Ω_i − Ω_o)/Ω_i, and neglect the back-reaction of the magnetic and viscous torques, so that ΔΩ is fixed.

At the boundaries, we impose no-penetration (v_r = 0, v_θ = 0) and no-slip (v_ϕ = Ω_{i, o}R_{i, o}sin θ) conditions for all t. We model the fluid and the inner and outer boundaries as perfect conductors. The magnetic field is anchored to the boundaries and frozen into the fluid. We denote with B₀ the characteristic magnetic field strength at t = 0.

Taken at face value, neutron star interiors have high Reynolds numbers (Re ≈ 10^11;Mastrano & Melatos 2005; Melatos & Peralta 2007). It is possible to argue to the contrary that Re is low in the frame of the stellar crust, given a shear of 10⁻¹⁰ ≤ ΔΩ/Ω ≤ 10⁻⁵ as inferred from measurements of rotational glitches (Espinoza et al. 2011). However, this assumption is debatable, as glitches may relax only a small fraction of the shear, and the poorly known effective viscosity of the neutron superfluid component may be lower than expected. Either way, computational limitations restrict us to Re = 500 in order to avoid unresolvable turbulent eddies and non-axisymmetric instabilities (Hollerbach & Skinner 2001; Nakabayashi, Tsuchida & Zheng 2002a; Gissinger et al. 2011). For the same reasons, the shell cannot be too thick (Nakabayashi et al. 2002b; Nakabayashi & Tsuchida 2005). We take δ = (R_o − R_i)/R_i = 1.86 in this paper (Peralta et al. 2005). A brief discussion of how higher Reynolds numbers (or a tilted magnetic axis) lead to hard-to-resolve non-axisymmetric and even turbulent flows is presented in Section 7.3.

3.3 Magnetic and viscous time-scales

4 VALIDATION: UNMAGNETIZED FLOW

We test the numerical code by comparing its output with unmagnetized pseudospectral simulations obeying the same initial and boundary conditions (Peralta et al. 2005, 2009; Howitt et al. 2016). As well as validating pluto for our problem against an independent solver, this also provides a control experiment that sets a baseline against which we compare the magnetized flow.

In spherical Couette flow, in the regime of fast rotation and low viscosity, the Taylor–Proudman theorem (Proudman 1956) states that the flow is approximately azimuthal and columnar (i.e. gradients are small parallel to the rotation axis). The imposed differential rotation drives a secondary flow via Ekman pumping, with v_r/v_ϕ∝Re^−1/2. Fig. 1 confirms this behaviour for ΔΩ = 0.1. The right half of the meridional plane displays angular velocity contours, which are approximately columnar. A cylindrical Stewartson layer touches the inner sphere at the equator. Fluid rotates with v_ϕ ≈ Ω_iR_i inside the Stewartson layer and v_ϕ ≈ Ω_oR_o outside the Stewartson layer. The left half of the meridional plane displays streamlines of the meridional recirculation. The fluid is pumped from the inner, faster sphere down towards the equator, then outwards to the outer, slower sphere, and thence to the poles. Close to the equator, the streamlines trace two circulation cells per hemisphere, where the fluid is sucked out of the equatorial plane, driven to higher latitudes, and recycled back to the equator.

Figure 1.

Primary rotation and secondary flow driven by Ekman pumping for an unmagnetized, viscous fluid in spherical Couette geometry with ΔΩ = 0.1 at t = 1.4t_μ. (Right half) Contours of the angular velocity Ω (colour bar in units of Ω_i). The red curves mark the contour levels in steps of 0.01Ω_i. (Left half) In-plane streamlines. The fluid satisfies the Taylor–Proudman theorem; the angular velocity is almost independent of the latitude in much of the shell. The figure looks identical for B₀ = 0 and 1.5 × 10⁻⁵ ≤ B₀ ≤ 1.5 × 10⁻², if |$\boldsymbol{B}$| is purely toroidal at t = 0.

Open in new tabDownload slide

Next we verify that the flow in a purely toroidal magnetic field evolves, as if the magnetic field is absent. To do so, we produce a version of Fig. 1 with ΔΩ = 0.1, 1.5 × 10⁻⁵ ≤ B₀ ≤ 1.5 × 10⁻², and |$\boldsymbol{B}_p = 0$| at t = 0. The result (not shown) is completely indistinguishable from Fig. 1. This is because (i) the magnetic tension term in equation (5) vanishes locally for |$\boldsymbol{B} = B_{\phi }\hat{\phi }$|⁠, independent of B₀, and (ii) the form of the induction equation guarantees |$\boldsymbol{B}_{\rm p} = 0$| at t > 0, given |$\boldsymbol{B}_{\rm p} = 0$| at t = 0.

5 MAGNETIZED FLOW: COROTATION WITH THE CRUST

We now turn on the magnetic field and explore numerically the effect of its topology on the flow.

Figure 2.

Magnetic topology of the initial dipole configuration given by (12). Near the poles, open field lines connect the differentially rotating inner and outer spheres. At lower latitudes, open field lines connect the northern and southern hemispheres of the outer sphere without touching the inner sphere. Some small, closed poloidal loops also exist at r ≈ R_o, θ ≈ π/2.

Open in new tabDownload slide

Magnetic tension enforces corotation of the fluid with the outer sphere, along those field lines that do not touch the inner sphere. Fig. 3 displays the angular velocity Ω, and linear velocity components v_r and v_θ on the equatorial plane, as functions of r. When the magnetic field is weak, i.e. B₀ ≤ 1.5 × 10⁻³, the solution cannot be distinguished from the unmagnetized one (blue curve) or else resembles it closely (green curve). However, for B₀ ≥ 5.9 × 10⁻³, there are clear differences. Fig. 3(a), which plots Ω(r), shows that |$\boldsymbol{B}$| enforces corotation with the outer sphere over a wider range of radii (0.65 ≲ r ≲ 1.54) than in the unmagnetized flow (1.0 ≲ r ≲ 1.54).

Figure 3.

Approximate corotation with the outer sphere for an initial dipole magnetic field (initial strength 1.5 × 10⁻⁵ ≤ B₀ ≤ 1.5 × 10⁻²). (Panel a) Angular velocity as a function of radius r at the equator (θ ≈ π/2). As B₀ increases, more fluid corotates with the outer shell. (Panel b) Radial velocity component as a function of r at θ ≈ π/2. (Panel c) Latitudinal velocity component as a function of r at θ ≈ π/2. (Panel d) Angular velocity contours for B₀ = 1.5 × 10⁻², with poloidal magnetic field lines overplotted (red contours). The closed poloidal loops at θ ≈ π/2 are pushed closer to the outer boundary by the open field lines with endpoints at θ ≈ 5π/12 and ≈ 7π/12 on the outer sphere. All panels are snapshots at t = 1.4t_μ for ΔΩ = 0.1.

Open in new tabDownload slide

Fig. 3(b) shows the corresponding behaviour for v_r. In the weak field regime B₀ ≤ 1.5 × 10⁻³ (blue and green curves), v_r peaks due to Ekman pumping at the Stewartson layer near the inner sphere (r ≈ 0.75) as in the unmagnetized flow in Fig. 1. For B₀ ≥ 5.9 × 10⁻³, however, v_r is suppressed near the inner sphere and peaks near the outer sphere (r ≈ 1.3). The peak shifts because the secondary circulation cells move away from the inner sphere towards the outer sphere. The maximum |v_r| in this radial jet increases with B₀. The v_θ component in Fig. 3 (c) is between one and two orders of magnitude smaller than v_r. The equatorial symmetry of the system imposes v_θ = 0 at the equator, separating the circulation in the two hemispheres. For B₀ ≤ 1.5 × 10⁻³, v_θ closely follows the unmagnetized solution. The peak in v_θ does not move monotonically with B₀. Increasing B₀ from 5.9 × 10⁻³ up to 1.2 × 10⁻² produces local extrema of v_θ at r ≈ 1.1 and 1.4, respectively, which are suppressed when B₀ increases further to 1.5 × 10⁻².

Fig. 3(d) demonstrates how Ω(r, θ) (colour scale) relates to the poloidal magnetic field lines (red contours). It generalizes Fig. 3(a) to show Ω(r, θ) throughout the spherical shell for B₀ = 1.5 × 10⁻² initially and ΔΩ = 0.1, i.e. when the magnetic field is dynamically important. In Appendix  A, for completeness, we present analogous plots for lower magnetic fields in the range 1.5 × 10⁻⁵ ≤ B₀ ≤ 1.5 × 10⁻². The dark blue region, in which the fluid corotates with the outer boundary, is threaded by magnetic field lines that connect the northern and southern outer hemispheres without touching the inner shell. Magnetic tension enforces approximate corotation with the outer shell. Open field lines that touch the outer shell at latitudes θ ≈ 5π/12 and 7π/12 are deformed close to the equator, squeezing the poloidal loop on the outer boundary. Accordingly, Ω rises slightly at r ≈ 1.3 as seen in Fig. 3(a). The angular velocity exceeds Ω_i > Ω_o by up to |$20{{\ \rm per\ cent}}$| close to the poles, where magnetic field lines connect the inner and the outer sphere.

The reader may wonder why v_θ does not vanish at θ = π/2 in Fig. 3(c), as expected from the north–south symmetry. Do we introduce a symmetry-breaking perturbation into the system (beyond those that arise unavoidably from numerical errors)? The answer is no: v_θ is not exactly zero because it is sampled slightly away from the equator, in a region where the gradient ∂v_θ/∂θ is high. The curves in Fig. 3(c) are for θ = 1.5733 rad ≠ π/2, which is as close as one gets to the equator with the grid placement and resolution chosen (206 points in θ). By interpolating the solution across the equator, we obtain v_θ ≈ 0 to a good approximation at θ = 1.5708 = π/2 exactly, with |v_θ(θ = π/2)| ≤ 0.01|v_θ(θ = 1.5733)|. Furthermore, the v_θ pattern is north–south symmetric but switches in sign about the equator, as expected for equatorial symmetry. Secondary flows circulate fluid from the equator towards the north pole in the northern hemisphere, and from the equator to the south pole in the southern hemisphere, as in unmagnetized spherical Couette flow.

6 MAGNETIZED FLOW: DIFFERENTIALLY ROTATING, MAGNETICALLY COUPLED CRUST AND CORE

Figure 4.

Magnetic topology of the initial dipole configuration given by (13). Near the poles, open field lines connect the inner and outer spheres. At lower latitudes, the magnetic field lines connect the inner sphere to itself.

Open in new tabDownload slide

Fig. 5 displays Ω, v_r, and v_θ as functions of radius in the equatorial plane for 1.5 × 10⁻⁵ ≤ B₀ ≤ 1.5 × 10⁻² and ΔΩ = 0.1. The physics is similar to Fig. 3, but the results are different because the magnetic topology is different. For B₀ ≤ 1.5 × 10⁻³, when the magnetic field is unimportant dynamically, the flow closely resembles the unmagnetized case. For 5.9 × 10⁻³ ≤ B₀ ≤ 1.5 × 10⁻², the closed equatorial zone strives to corotate with the inner sphere, in marked contrast to Fig. 3. In doing so, it overshoots, and Ω peaks at Ω ≈ 1.05 in the interval 0.7 ≲ r ≲ 0.8. This occurs because currents are free to flow through the rigid, conducting walls of the spherical shell (Dormy et al. 2002; Soward & Dormy 2010). In order to compensate, the fluid generates additional currents via rotation, allowing the fluid to deviate from Ferraro’s law of isorotation (Ferraro 1937) and to attain max (Ω) > Ω_i (Nataf et al. 2006). The peak of Ω moves from r ≈ 0.7 at B₀ = 5.9 × 10⁻³ to r ≈ 0.8 at B₀ = 1.5 × 10⁻², as the magnetic coupling to the inner sphere strengthens. As in Fig. 3(b), the peak of v_r ≈ 0.1rΩ at B₀ = 5.9 × 10⁻³ shifts outwards as B₀ increases because the principal circulation cells of the meridional flow move closer to the outer boundary. The v_θ component in Fig. 5 (c) is small, with |v_θ| ≈ 10⁻³rΩ, again because of equatorial symmetry.

Figure 5.

(Panel a) Angular velocity as a function of r at the equator (θ ≈ π/2). As the magnetic field strengthens, one finds max(Ω) >Ω_i. (Panel b) Radial velocity of the fluid as a function of the radius at θ ≈ π/2. (Panel c) Latitudinal velocity v_θ(r) for θ ≈ π/2. B₀ varies in the range [1.5 × 10⁻⁵, 1.5 × 10⁻²]. The plots are snapshots taken at t = 1.4t_μ for ΔΩ = 0.1.

Open in new tabDownload slide

We now double the shear to ΔΩ = 0.2 in Fig. 6. The greater the shear, the more magnetic field lines joining the differentially rotating shells are stretched, and the harder it is for the magnetic tension to enforce corotation. As a result, we find Ω > Ω_i for a narrower range of r (0.6 ≲ r ≲ 0.8; cf. Fig. 5a) and higher values of B₀ (B₀ > 7.7 × 10⁻³; cf. Fig. 5a). Interestingly, the peak of Ω is not only narrower but taller, reaching max (Ω) ≈ 1.2 [cf. max (Ω) ≈ 1.05 in Fig. 5a]. Something similar happens to the radial jet in Fig. 6 (b). Its profile exhibits a broad plateau at 0.9 ≲ r ≲ 1.2 for B₀ = 1.3 × 10⁻², with max (v_r) ≈ 0.25, which is broader and faster than the radial jet in Fig. 5(b). The latitudinal velocity v_θ remains like in Fig. 5(b) due to equatorial symmetry. Contours Ω(r, θ) for different B₀ with ΔΩ = 0.1 and 0.2 are presented in Appendix  A for completeness.

Figure 6.

Same as in Fig. 5 but with higher rotational shear (ΔΩ = 0.2) and t = 1.6t_μ.

Open in new tabDownload slide

The astrophysical interpretation of the above results is as follows. If the inner sphere corresponds to a distinct, rigid, stellar component, such as a crystalline colour superconducting core (Alford, Bowers & Rajagopal 2001; Mannarelli, Rajagopal & Sharma 2007; Alford et al. 2008), then Figs 5 and 6 imply that there is a region of differential rotation at r ≈ R_i. The differential rotation persists because the spin-down torque on the outer sphere maintains Ω_i ≠ Ω_o. The magnitude of the lag, Ω_i − Ω_o, depends on the moment of inertia of the inner sphere, e.g. the size and composition of the crystalline superconducting core. On the other hand, if the inner core is not rigid, r = R_i corresponds to a mathematical surface introduced to promote numerical stability in the simulation (Nakabayashi et al. 2002b; Nakabayashi & Tsuchida 2005; Peralta et al. 2005), and the region r ≤ R_i in the real star contains the same fluid as R_i ≤ r ≤ R_o. In the latter scenario, the differential rotation is expected to be weaker than in Figs 5 and 6, as the field lines connecting the outer and inner spheres enforce approximate corotation in the stellar volumes 0 ≤ r ≤ R_i and R_i ≤ r ≤ R_o. However, exact corotation is never achieved for two reasons. First, the hydromagnetic torques do not act instantaneously; in the time t_A that Alfvén waves take to propagate into the core, the outer sphere spins down electromagnetically, and one has |$\Omega _{\rm i} - \Omega _{\rm o} \approx t_{\rm A}\dot{\Omega }_{\rm o}$|⁠. Secondly, taking the magnetic topology in Fig. 4 as an example, there is a region near the equator that is magnetically disconnected from the outer sphere. This is the region threaded by magnetic field lines that start and end at r = R_i in the (artificial) simulation and close within the core without touching r = R_o in the real star. Corotation in a magnetically disconnected volume takes longer to achieve, as discussed in Section 7.

7 EVOLUTION OF THE MAGNETIC GEOMETRY

7.1 Approximately closed equatorial flux tubes

Does an initially open magnetic geometry develop closed (or approximately closed) flux tubes under differential rotation? This question is important because once a flux tube closes, the fluid inside can rotate differentially on the long, viscous time-scale t_μ, even if the rest of the star corotates (Melatos 2012; Glampedakis & Lasky 2015). Indeed, we show below that the fluid inside a closed flux tube evolves, as if one has |$|\boldsymbol{B}| \approx 0$| in the tube; it adjusts to the boundary conditions on the flux tube surface (set by the rest of the flow) on the time-scale t_μ. The qualifier ‘approximately’ is important because, strictly speaking, a change in topology is impossible in ideal MHD, where field lines cannot break and reconnect. In astrophysical reality, dissipative processes permit reconnection.

Figs 7 and 8 show the meridional field lines evolving with time for two choices of ΔΩ and B₀. At t = 0.56t_μ (Fig. 7) and t = 0.32t_μ (Fig. 8), two closed meridional loops appear just above and below the equator. They persist until t = 1.4t_μ and 1.6t_μ (end of the simulations). Inside the volume of revolution whose meridional cross-section is enclosed by these loops, the field lines have |B_ϕ| ≫ |B_r|, |B_θ| and are approximately ‘closed’ in the sense described below. Meridional loops, which resemble peninsulas (i.e. which do not pinch off completely), also develop at high latitudes, to the left of the vertical Stewartson layer which is tangent to the inner sphere at the equator (Peralta & Melatos 2009). Note that a field line is not closed in three dimensions in general, just because its projection on to the meridional plane is closed. For example, the two small, meridional loops near the equator in Figs 7 and 8 are projections of predominantly toroidal field lines that perform multiple revolutions without closing, as described below. Moreover, the meridional loops themselves appear to change their connectedness when the density of contours in the plot increases or the grid resolution in the simulations increases. These small-scale numerical effects are accompanied by small-scale, non-ideal-MHD diffusive effects due to resistivity and viscosity in realistic physical flows. Viscosity and boundary conditions on the surface of these regions control the velocity gradients inside, as opposed to magnetic tension elsewhere, as discussed in Section 7.2.

Figure 7.

Evolution of the magnetic geometry for ΔΩ = 0.1 and B₀ = 1.5 × 10⁻².

Open in new tabDownload slide

Figure 8.

Evolution of the magnetic geometry for ΔΩ = 0.2 and B₀ = 1.3 × 10⁻².

Open in new tabDownload slide

In what sense are the equatorial flux tubes enclosed by the loops to be regarded as approximately closed? Inside the closed flux tube, the toroidal magnetic field component typically dominates the poloidal one. We parametrize the magnetic field lines by [r(χ), θ(χ), ϕ(χ)], where χ is the arc length in magnetic coordinates. For 0 ≤ ϕ(χ) ≤ 2π, the functions r(χ) and θ(χ) oscillate periodically in χ, with the oscillation period χ_max. In order to establish whether the magnetic field lines in the flux tube are closed or not, we compare the values of r(χ), θ(χ), and ϕ(χ) at χ = χ_min = 0 and χ = χ_max. We consider the magnetic field line to be approximately closed if |$|r(0) - r(\chi _{\max })|/r(0) \le 0.5{{\ \rm per\ cent}}$|⁠, |$|\theta (0) - \theta (\chi _{\max })|/\theta (0) \le 0.1{{\ \rm per\ cent}}$| and the field line completes at least one orbit [ϕ(χ_max) ≥ 2π].

Table 1 presents examples of three categories of field lines. The field lines termed ‘open’ connect the inner sphere to itself at different latitudes. Field lines termed ‘periodic’ are open, predominatly toroidal, and orbit around the inner sphere, with r(χ), θ(χ) oscillating periodically as |ϕ(χ)| increases. In this case, χ_max corresponds to the oscillation period of r(χ) and θ(χ). ‘Closed’ field lines return approximately to where they started with |$|r(0) - r(\chi _{\max })|/r(0) \le 0.5{{\ \rm per\ cent}}$|⁠, |$|\theta (0) - \theta (\chi _{\max })|/\theta (0) \le 0.1{{\ \rm per\ cent}}$|⁠, and ϕ(χ_max) ≥ 2π.

In Figs 9 and 10, we plot the coordinates of intersection at ϕ = 0 of the approximately closed field lines (green dots), open field lines that are predominantly toroidal and touch neither r = R_i nor R_o (blue dots), and open poloidal field lines (red dots) as defined in the previous paragraph. In both figures, we find that at least one closed toroidal field line forms from the initial open topology (13). Around the closed line, there is a bundle of predominantly toroidal (|B_r|, |B_θ| ≪ |B_ϕ|) but open field lines, for which r(χ) and θ(χ) oscillate periodically. The oscillation centre stays within the black dot–dashed circle over many periods. The blue dots that fall within the solid green circle containing the green dot deserve a special mention. They also indicate open and periodic field lines, whose barycentre tends towards the closed field line [over many oscillation periods (≈10)]. The solid green circle and dot–dashed black circles are smaller for ΔΩ = 0.1 than for ΔΩ = 0.2; the higher shear produces a fatter toroidal flux tube. A detailed study of the equatorial flux tubes is presented in Appendix B.

Figure 9.

Geometry of the magnetic field at t = 1.4t_μ when the rotational shear is ΔΩ = 0.1 and the initial configuration is a central point dipole of the form (13) with B₀ = 1.5 × 10⁻². The red dots are the footpoints of open field lines connecting the inner and the outer spheres for θ ≤ 0.7 rad. For θ ≥ 1.0 rad, the magnetic field lines start at the inner sphere and do not touch the outer sphere. The blue dots represent intersection points with the half-plane ϕ = 0 for open field lines that do not touch the boundaries and that are almost toroidal. The green dot indicates a closed toroidal field line. The blue dots falling within the green circle drift towards the green dot within ten oscillation periods.

Open in new tabDownload slide

Figure 10.

As in Fig. 9 but for t = 1.6t_μ, ΔΩ = 0.2, and B₀ = 1.3 × 10⁻².

Open in new tabDownload slide

7.2 Field-line-averaged tension

The approximately closed equatorial flux tubes identified in Section 7.1 are not special in the sense that they are out of the ordinary. But they are important because, inside them, the differential rotation is determined differently qualitatively to the rest of the flow. That is, viscosity and boundary conditions on the surface of these regions control the velocity gradients inside, as opposed to magnetic tension elsewhere.

Figure 11.

Absolute value of I_B/I_T along open field lines with footpoints at r = R_i and 0.1 ≤ θ ≤ 1.3 rad, versus B₀ at t = 1.4t_μ with ΔΩ = 0.1. The initial configuration is given by (13). When the ratio reaches unity, magnetic forces are comparable to viscous and inertial forces in a field-line-averaged sense.

Open in new tabDownload slide

Figure 12.

As in Fig. 11 with ΔΩ = 0.2 at t = 1.6t_μ.

Open in new tabDownload slide

We now evaluate the ratio I_B/I_T for field lines in the closed flux tube. Our choice of the criteria |$|r(0) - r(\chi _{\max })|/r(0) \le 0.5{{\ \rm per\ cent}}$|⁠, |$|\theta (0) - \theta (\chi _{\max })|/\theta (0) \le 0.1{{\ \rm per\ cent}}$| to distinguish between purely toroidal, closed field lines and periodic, predominantly toroidal field lines is consistent with I_B/I_T. For ΔΩ = 0.1 (Fig. 9), we obtain I_B/I_T = O(10⁻²) for periodic field lines (blue dots) and I_B/I_T = O(10⁻⁴) for the closed field line (green dot). By doubling the shear (Fig. 10), the ratio I_B/I_T reaches O(10⁻³) for periodic field lines and I_B/I_T = O(10⁻⁶) for the closed field line.

Field line integrals reveal the internal structure of magnetic flux tubes. By inspecting the order of magnitude of the ratio I_B/I_T within the flux tube, it is possible to (i) explore the internal structure of the flux tube and (ii) verify that the field-line-averaged magnetic tension inside the tube is small. Under condition (ii), the angular velocity of the fluid enclosed in the flux tube is determined by the boundary conditions on the surface of the flux tube, which are set by the surrounding fluid. The matching to the flux tube surface occurs via viscous forces on a time-scale longer than t_A.

7.3 Non-axisymmetry

Non-axisymmetry arises naturally even in unmagnetized spherical Couette flow (e.g. see the phase plane in fig. 1 in Nakabayashi et al. 2002a, among many other results), and in superfluid spherical Couette flow (Melatos & Peralta 2007).

In our paper, we restrict the simulations to relatively low Reynolds number (Re < 10³) and to parallel rotation and magnetic axes. We do this partly for simplicity in response to computational constraints and partly because our main application (neutron stars) is arguably a low-Reynolds-number problem with 10⁻¹⁰ ≤ ΔΩ/Ω ≤ 10⁻⁵ from measurements of rotational glitches (Espinoza et al. 2011). However, the low-Re assumption is certainly debatable. One can also argue that the problem is high-Re because the poorly known effective viscosity of the neutron superfluid may be lower than expected, and glitches may only relax a small fraction of the underlying shear. In that case, the flow may be fully turbulent (Melatos & Peralta 2007), and a different style of study would be needed.

Assuming low Re and parallel axes, the solution is axisymmetric, and the problem is effectively bidimensional. To test this assertion, we have performed three-dimensional runs, and verified that the solution is independent of the azimuthal coordinate ϕ. We have performed test runs with 180 points in the radial direction, 150 points in the θ direction, and 20 points in the ϕ direction for the configuration ΔΩ = 0.1 and B₀ = 1.5 × 10⁻². Both the radial velocity and the azimuthal velocity (not plotted here) are independent of the ϕ variable at r = 1.04 and θ = 0.61, 0.81, and 1.21 rad. So are all the other dependent variables, e.g. magnetic field components. The system starts axisymmetric and remains so as it evolves.

There are two astrophysically realistic ways to drive non-axisymmetry in the flow pattern and magnetic geometry. One way would be to tilt the magnetic axis with respect to the rotational axis, as implied by the radio polarization swings observed in pulsars. The numerical complexity of the problem depends strongly on B₀, Re, the tilt angle α, and the grid resolution in θ and ϕ required to resolve turbulent meridional flows and non-axisymmetries. We perform some three-dimensional test runs with resolution 202 × 176 × 50 (in r, θ, and ϕ, respectively) to investigate the behaviour informally; a full study lies outside the scope of the paper. We report in Fig. 13 the radial velocity at θ = 0.71, 0.90, and 1.08 rad for Re = 500, α = 0.3 rad, ΔΩ = 0.1, B₀ = 1.5 × 10⁻², r = 0.99, and t = 7.2t_E, where t_E is the Ekman time-scale. There is a clear dependence on ϕ at all three latitudes as well as hints of turbulence. In order to resolve the turbulent flow properly, a fine grid is required in r, θ, and ϕ, beyond the capacity of the computational resources at our disposal.

Figure 13.

Radial velocity as a function of ϕ for α = 0.3 rad, ΔΩ = 0.1, B₀ = 1.5 × 10⁻², r = 0.99, and t = 7.2t_E, at three latitudes (see the legend).

Open in new tabDownload slide

A second astrophysically plausible way to obtain non-axisymmetric flows is to increase the Reynolds number (Hollerbach 2009; Gissinger et al. 2011). Already for Re ≈2000, non-axisymmetric features emerge even in unmagnetized spherical Couette flow (Nakabayashi & Tsuchida 2005) and superfluid spherical Couette flow (Peralta et al. 2009). Again, a finer grid resolution is required to resolve the flow. Once one reaches Re ≈ 10⁵, the flow is fully turbulent (Nakabayashi et al. 2002a; Peralta et al. 2009). In an ordinary neutron star, where hydrodynamic forces dominate the Lorentz force typically, the magnetic field alsobecomes tangled too, once the flow is fully turbulent. The turbulent fate of the flow is clear qualitatively in this scenario, but several interesting points of detail arise, whose study we postpone for future work, when we gain access to more computational resources.

8 CONCLUSION

An important open question in neutron star astrophysics is how the topology of the star’s internal magnetic field affects the degree of differential rotation in the interior. This paper explores numerically an idealized model of a neutron star consisting of two differentially rotating, concentric spheres containing a magnetized, incompressible fluid described by ideal MHD using the solver pluto (Mignone et al. 2007). Previous analytic studies suggest that the fluid is forced into corotation with the rigid crust, if the field is open, but corotation is harder to maintain everywhere, if the field is closed (Easson 1979; Melatos 2012; Glampedakis & Lasky 2015). Differential rotation between the multiple fluid components and the rigid crust is a possible trigger for observed rotational irregularities like glitches (Warszawski & Melatos 2011; Andersson et al. 2012; Haskell & Melatos 2015). Numerical investigations of magnetized, spherical Couette flow have been carried out previously in the limit of small magnetic Reynolds number, where the magnetic perturbations induced by the flow are small compared to the externally imposed magnetic field (Dormy et al. 1998, 2002; Hollerbach et al. 2007; M. Soward & Dormy 2010; Gissinger et al. 2011). In this paper, we consider the opposite regime.

When the topology of the field is open and most of the internal fluid is threaded by field lines that connect the northern and southern hemispheres of the outer sphere (i.e. rigid crust), the fluid corotates with the outer sphere after a few Alfvén time-scales t_A. When open magnetic field lines connect the differentially rotating, inner and outer spheres, some of the fluid at intermediate radii rotates faster than both boundaries. If the inner sphere corresponds astrophysically to a solid inner core (Alford et al. 2001, 2008; Mannarelli et al. 2007), fluid along magnetic field lines connecting the inner and outer spheres rotates differentially, as electromagnetic and viscous torques maintain the crust–core rotational lag. If neutron stars do not have solid cores, regions in the inner sphere that are connected magnetically to the crust approach corotation up to corrections of order |$\dot{\Omega }_{\rm o}t_{\rm A}$|⁠.

It is shown that a central point dipole evolves to generate approximately closed toroidal flux tubes just above and below the equator. The equatorial flux tubes divide into a core (containing at least one approximately closed field line and open, predominantly toroidal and periodic field lines that asymptote to the closed field line) and an annular sheath surrounding it, filled with periodic field lines. For lower rotational shear, the periodic field lines stay within the sheath and out of the core. For higher rotational shear, some field lines leave the sheath, enter the core, and asymptote to the closed field line. The approximately closed equatorial flux tubes develop on roughly the viscous time-scale t_μ and coincide roughly with closed loops in the magnetic field projected into the meridional plane. Inside the approximately closed flux tubes, the magnetic tension integrated along magnetic field lines averages approximately to zero, and the field-line-averaged MHD momentum equation reduces to its hydrodynamical counterpart. Consequently, fluid within the closed flux tube corotates with the surrounding fluid on the viscous time-scale t_μ ≫ t_A. The magnetic tension does not average to zero along open field lines over a range of ΔΩ and B₀ values.

We close with two caveats. First, the results of the simulations depend on the initial configurations; different configurations are certain to evolve differently. The fields (12) and (13) are used widely in the general literature (Dormy et al. 1998, 2002; Nataf et al. 2006; Hollerbach et al. 2007; Schmitt et al. 2008; Soward & Dormy 2010) and neutron-star models (Bocquet et al. 1995; Braithwaite & Spruit 2004). They lead to two phenomena, which are potentially important astrophysically. (i) In some regions, the fluid rotates faster than both the inner and outer spheres; and (ii) in some magnetically disconnected regions the angular velocity is determined predominantly by viscous forces instead of magnetic tension, as predicted by others analytically (Easson 1979; Melatos 2012; Glampedakis & Lasky 2015). It is too early to say whether or not these phenomena are generic to most initial configurations, but their existence for one plausible configuration is interesting for future neutron star modeling. Alternative scenarios, e.g. fully developed turbulence, require a different sort of study.

The second caveat is that the field evolution presented in this paper has been observed before in previous studies of magnetized spherical Couette flow, albeit usually in the regime where the field component induced by the flow is small compared to the applied field. Nevertheless the results are an instructive addition to the neutron star literature. For a long time now, it has been customary for neutron star models to assume (for simplicity) a simple internal magnetic geometry, e.g. uniform magnetization or star-centred dipole. There is nothing wrong with this as a theoretical device of course, but Section 7 reminds readers that the reality is likely to be more complicated. In particular it illustrates that even relatively small amounts of crust–core differential rotation produce complicated magnetic geometries, a point that is fairly straightforward but nevertheless has not received much attention via numerical simulations in the literature. In these more complicated regions, the physics governing differential rotation is qualitatively different to elsewhere, with viscosity dominating magnetic tension as discussed above. The results in Section 7 complement other studies on complicated internal magnetic fields, e.g. Braithwaite & Spruit (2004) and Braithwaite & Spruit (2006), who did not consider differential rotation but did consider very interesting topologies; Akgün et al. (2018), Pons & Viganò (2019), and Carrasco et al. (2019), who focused on the crust and simulate the core; Drummond & Melatos (2017a,b), who looked at the related problem of spontaneous emergence of complicated neutron vortex structures in a type II superconductor; Ruderman, Zhu & Chen (1998), who studied surface multipoles and their tectonic evolution; and Sur, Haskell & Kuhn (2020), who studied the turbulent equilibrium of the magnetized stellar interior.

ACKNOWLEDGEMENTS

FA acknowledges the support of the University of Melbourne through a Melbourne Research Scholarship. The authors thank the anonymous referee for pointing out several relevant references on magnetized spherical Couette flow and for helpful suggestions on how to restructure and sharpen parts of the presentation.

REFERENCES

APPENDIX A:

ANGULAR VELOCITY PROFILES

In this appendix, we show for completeness Ω(r, θ) for different B₀, ΔΩ, and for the initial magnetic field topologies given by (12) and (13). In Figs A1–A3, we plot the contours of the angular velocity (colour bar in units of Ω_i). The red curves are the contour levels of Ω.

Figure A1.

Contour plots of the angular velocity of the fluid at t = 1.4t_μ for ΔΩ = 0.1 and an initial dipole given by (12). Contours of Ω are represented as red solid curves. As the initial magnetic field strength rises, we observe max (Ω) > Ω_i near the poles. For 0.4 ≲ θ ≲ 5.9 rad, the fluid corotates approximately with the outer sphere.

Open in new tabDownload slide

Figure A2.

Contours of the angular velocity of the fluid (red curves) for ΔΩ = 0.1 at t = 1.4t_μ. As the initial magnetic field strengthens, Ω > Ω_i for B₀ ≥ 5.9 × 10⁻³. A small north–south asymmetry is observed at the poles for B₀ ≥ 1.2 × 10⁻³.

Open in new tabDownload slide

Figure A3.

As in Fig. A2 but for ΔΩ = 0.2 at t = 1.6t_μ.

Open in new tabDownload slide

Fig. A1 corresponds to the initial magnetic field configuration given by (12). As expected, a weak magnetic field does not lead to appreciable differences between the magnetized and the unmagnetized cases. For B₀ ≥ 1.5 × 10⁻³, differences are clearly observed. The fluid spins faster than the inner sphere at the poles for B₀ = 5.9 × 10⁻³, and Ω keeps growing up to B₀ = 8.8 × 10⁻³. For 0.4 ≲ θ ≲ 5.9 rad, the fluid is approximately in solid body rotation with the outer sphere, independently of B₀. The border separating regions with max (Ω) > Ω_i from corotating regions is delimited, for increasing θ, by the first magnetic field line that does not touch the inner sphere. The presence of poloidal closed field lines in the initial configuration does not prevent corotation with the outer sphere.

Figs A2 and A3 show Ω(r, θ) across the spherical shell for ΔΩ = 0.1 and 0.2 and 1.3 × 10⁻⁵ ≤ B₀ ≤ 1.5 × 10⁻². For ΔΩ = 0.1 and B₀ ≈ 10⁻⁵, the solution is indistinguishable from the unmagnetized case in Section 4. For B₀ ≥ 5.9 × 10⁻³, we obtain Ω > Ω_i in certain regions. For B₀ ≥ 1.2 × 10⁻², a small north–south asymmetry is observed, which seems numerical in origin.

In Fig. A3, for ΔΩ = 0.2, we find max (Ω) > Ω_i for B₀ ≥ 1.0 × 10⁻². Unlike ΔΩ = 0.1, the locations with max (Ω) > Ω_i lie near the equator rather than at the poles (Soward & Dormy 2010). As for (12), the fluid rotates differentially along open field lines connecting the inner and outer spheres.

APPENDIX B:

DETAILED STRUCTURE OF THE APPROXIMATELY CLOSED EQUATORIAL FLUX TUBES

In Section 7.1, we find that field lines that lie within the green circles in Figs 9 and 10 stay within the green circles and indeed asymptote to the nearly closed field lines marked by a green dot. We now ask what happens to field lines that lie outside the green circle but inside the dot–dashed black circle. The question is studied in Fig. B1 for ΔΩ = 0.1. We choose a new colour for each blue dot outside the solid green circle in Fig. 9 and follow 12 field line revolutions for each.

Figure B1.

Successive intersection points with the ϕ = 0 half-plane of the blue dots outside the solid green circle shown in Fig. 9 for ΔΩ = 0.1 and B₀ = 1.5 × 10⁻². Every intersection point corresponds to a 2π-revolution in ϕ. There are 12 revolutions shown per field line. Each colour marks a separate field line. The black dot–dashed circle is densely filled with points, but no field line strays into the green circle.

Open in new tabDownload slide

We see in Fig. B1 that the successive intersection points with the half-plane ϕ = 0 stay within the black dot–dashed circle and swirl around the green circle without entering the green circle. In other words, the geometry near the green dot separates into two parts: a core, where field lines asymptote to the green dot, and a separate annular sheath around the core. Similarly, in Fig. B2, we track 12 revolutions for each of the blue dots outside the green circle in Fig. 10 for ΔΩ = 0.2. The magnetic geometry is richer than for ΔΩ = 0.1. The magnetic field lines spread across a larger volume, and the clustering of intersection points is less dense. Some points (red hexagons) enter the green circle after starting outside it unlike in Fig. B1, and their barycentre tends subsequently towards the closed field line. In summary, a magnetic field given by (13) evolves to generate a toroidal equatorial flux tube. Inside the flux tube, we find that one purely toroidal, closed field line is surrounded by a bundle of predominantly toroidal but open field lines. In order to qualify as closed, the field line meets the conditions |$|r(0) - r(\chi _{\max })|/r(0) \le 0.5{{\ \rm per\ cent}}$| and |$|\theta (0) - \theta (\chi _{\max })|/\theta (0) \le 0.1{{\ \rm per\ cent}}$|⁠. In Section 7.2, we show that these thresholds can be related to field line integrals; if r(χ) and θ(χ) satisfy the above conditions, the field-line-averaged magnetic tension vanishes.

Figure B2.

As in Fig. B1 but for ΔΩ = 0.2 and B₀ = 1.3 × 10⁻². Unlike in Fig. B1, some of the intersection points enter the green solid circle.

Open in new tabDownload slide