Current closure through the neutron star crust

TL;DR

This paper investigates how electric currents flow within neutron star crusts, revealing differences between millisecond and slower pulsars, and assessing the potential for crustal decoupling under strong magnetic fields.

Contribution

It models the internal electric current paths in neutron stars, highlighting how they depend on pulsar spin periods and magnetic field strengths, and evaluates associated crustal stresses.

Findings

Current reaches the crust base in millisecond pulsars.

In slower pulsars, currents do not penetrate beyond 100 m.

Maxwell stresses are below crustal failure threshold in typical pulsars, but may cause decoupling in magnetars.

Abstract

Force-free pulsar magnetospheres develop a large scale poloidal electric current circuit that flows along open magnetic field lines from the neutron star to the termination shock. The electric current closes through the interior of the neutron star where it provides the torque that spins-down the star. In the present work, we study the internal electric current in an axisymmetric rotator. We evaluate the path of the electric current by requiring the minimization of internal Ohmic losses. We find that, in millisecond pulsars, the current reaches the base of the crust, while in pulsars with periods of a few seconds, the bulk of the electric current does not penetrate deeper than about $100$ m. The region of maximum spin-down torque in millisecond pulsars is the base of the crust, while in slowly spinning ones it is the outer crust. We evaluate the corresponding Maxwell stresses and find that, in typical rotation-powered radio pulsars, they are well below the critical stress that can be sustained by the crust. For magnetar-level fields, the Maxwell stresses near the surface are comparable to the critical stress and may lead to the decoupling of the crust from the rest of the stellar rotation.

Tables

Table 1. Table 1: Summary of the models studied. The first column is the name of the model. P 𝑃 P and P ˙ ˙ 𝑃 \dot{P} are the period and period derivative, B 𝐵 B is the spin-down dipole magnetic field of an orthogonal rotator in vacuum corresponding to the period and period derivative mentioned before, θ pc subscript 𝜃 pc \theta_{\rm pc} is the semi-opening angle of the polar cap, the depth is the lowest value of r 𝑟 r for I = 1 2 ​ I max 𝐼 1 2 subscript 𝐼 max I=\frac{1}{2}I_{\rm max} , N tot , T subscript 𝑁 tot T N_{\rm tot,~{}T} is the torque obtained through the integration of equation 14 using the boundary condition from Timokhin ( 2006 ) and N tot , G subscript 𝑁 tot G N_{\rm tot,~{}G} is the torque using the polynomial approximation of Gralla et al. ( 2016 ) , N align subscript 𝑁 align N_{\rm align} is the torque evaluated through equation 18 , P Ohm subscript 𝑃 Ohm P_{\rm Ohm} is the total Ohmic losses using the polynomial approximation of Gralla et al. ( 2016 ) . The last column indicates whether the conductivity is set equal to a constant (c) or depends on radius according to equation 11 (v).

Model	$P$	$\dot{P}$	$B$	${\theta }_{\mathrm{pc}}$	depth	$N_{\mathrm{tot},T}$	$N_{\mathrm{tot},G}$	$N_{\mathrm{align}}$	$P_{\mathrm{Ohm}}$	$\sigma$
	(s)		(G)	(deg)	($\times 100$ cm)	(dyn cm)	(dyn cm)	(dyn cm)	(erg/s)
A1	0.01	9.77e-18	${10}^{10}$	9.24	740	8.7e32	8.5e32	8.7e32	1.2e19	v
A2	0.01	9.77e-18	${10}^{10}$	9.24	280	8.8e32	8.6e32	8.7e32	2.5e16	c
B1	0.1	9.77e-15	${10}^{12}$	2.91	310	8.8e32	8.6e33	8.7e33	4.9e19	v
B2	0.1	9.77e-15	${10}^{12}$	2.91	100	8.8e32	8.6e33	8.7e33	6.6e16	c
C1	1	9.77e-16	${10}^{12}$	0.92	100	8.8e32	8.7e30	8.7e30	2.1e16	v
C2	1	9.77e-16	${10}^{12}$	0.92	40	8.7e32	8.5e30	8.7e30	1.7e13	c
D1	5	1.95e-12	${10}^{14}$	0.41	50	7.1e32	7.1e32	7.0e32	9.4e17	v
D2	5	1.95e-12	${10}^{14}$	0.41	20	6.8e32	6.6e32	7.0e32	4.7e14	c
SGR 1806$-$20	7.5	4.95e-12	$2\times {10}^{15}$	0.33	40	8.4e32	8.0e34	8.2e34	6.2e19	v

Equations

$$\mathbf{j}=\frac{1}{2\pi }\nabla I\times \nabla \varphi +j_{\varphi }\hat{\varphi },$$

$$I=\frac{c}{2}r\sin \theta B_{\varphi },$$

$$\mathbf{E}=-r\sin \theta \Omega \hat{\mathbf{\varphi }}\times \mathbf{B}/c+\frac{\mathbf{j}}{\sigma }.$$

$$\mathbf{j}\cdot \mathbf{E}=-r\sin \theta \Omega \left(\hat{\mathbf{\varphi }}\times \mathbf{B}/c\right)\cdot \mathbf{j}+\frac{{\mathbf{j}}^2}{\sigma }=\frac{1}{c}\left(\mathbf{j}\times \mathbf{B}\right)\cdot \mathbf{v}+\frac{j^2}{\sigma },$$

$$P_{\mathrm{Ohm}}\equiv {\int }_V\frac{j^2}{\sigma }\mathrm{d}V={\left(\frac{c}{4\pi }\right)}^2{\int }_V\frac{{\left(\nabla \times \mathbf{B}\right)}^2}{\sigma }\mathrm{d}V.$$

$$\frac{c}{4\pi }\nabla \times \left(\frac{\nabla \times \mathbf{B}}{\sigma }\right)\equiv \frac{c}{4\pi }\nabla \times \left(\frac{\mathbf{j}}{\sigma }\right)=0$$

$$\nabla \times \left(\frac{\mathbf{j}}{\sigma }\right)=0,$$

$$\frac{{\partial }^{2}I}{\partial r^2}-\frac{1}{\sigma }\frac{\partial \sigma }{\partial r}\frac{\partial I}{\partial r}-\frac{\cos \theta }{r^2\sin \theta }\frac{\partial I}{\partial \theta }+\frac{1}{r^2}\frac{{\partial }^{2}I}{\partial {\theta }^2}=0$$

$$I(\Psi )=\frac{\Psi }{{\Psi }_0}\left[2-\frac{\Psi }{{\Psi }_0}-\frac{1}{5}{\left(\frac{\Psi }{{\Psi }_0}\right)}^3\right],$$

$$\rho =\left(1+{\left(\frac{r_{\mathrm{ns}}-r}{r_{\mathrm{ns}}-r_{\mathrm{in}}}\right)}^4\frac{{\rho }_{\mathrm{in}}}{{\rho }_{\mathrm{out}}}\right){\rho }_{\mathrm{out}}.$$

$$\sigma (r)={\sigma }_{\mathrm{out}}{\left(\frac{\rho (r)}{\rho (r_{\mathrm{ns}})}\right)}^{9/10},$$

$${\mathbf{B}}_p(r,\theta )=B{r}_{\mathrm{ns}}^3\left(\frac{\cos \theta }{r^3}\hat{\mathbf{r}}+\frac{\sin \theta }{2r^3}\hat{\mathbf{\theta }}\right).$$

$${\mathbf{F}}_L(r,\theta )$$

$$\mathbf{N}(r,\theta )$$

$$N_{\mathrm{tot}}={\int }_V(N_r\cos \theta -N_{\theta }\sin \theta )dV.$$

$$M_{r\varphi }=\frac{B_r(r,\theta )B_{\varphi }(r,\theta )}{4\pi }\text{mbox}and{M}_{\theta \varphi }=\frac{B_{\theta }(r,\theta )B_{\varphi }(r,\theta )}{4\pi }$$

$${\tau }_{br}=\left(0.0195-\frac{1.27}{\Gamma -71}\right)n_i\frac{Z^2{e}^2}{\alpha }$$

$$N_{\mathrm{align}}=0.94\times \frac{2}{3}\frac{\Omega }{c}{\left(\frac{1.23r_{\mathrm{ns}}^{3}B}{r_{\mathrm{lc}}}\right)}^2.$$

$$E(\mathbf{B})={\int }_V\frac{{\left(\nabla \times \mathbf{B}\right)}^2}{\sigma }dV$$

$$V(\mathbf{B},\mathbf{h})$$

$$\mathbf{A}=\frac{\nabla \times \mathbf{B}}{\sigma },$$

$$\nabla \cdot \left(\mathbf{A}\times \mathbf{h}\right)=-\mathbf{A}\cdot \nabla \times \mathbf{h}+\mathbf{h}\cdot \nabla \times \mathbf{A}.$$

$${\int }_V\mathbf{A}\cdot \nabla \times \mathbf{h}dV=$$

$${\int }_V\mathbf{h}\cdot \left(\nabla \times \mathbf{A}\right)dV-{\int }_V\nabla \cdot \left(\mathbf{A}\times \mathbf{h}\right)dV=$$

$${\int }_V\mathbf{h}\cdot \left(\nabla \times \mathbf{A}\right)dV+{\int }_{\partial V}\left(\mathbf{A}\times \mathbf{h}\right)\cdot d\mathbf{S}$$

$$\nabla \times \left(\frac{\nabla \times \mathbf{B}}{\sigma }\right)=0,$$

Full text

††pagerange: Current closure through the neutron star crust–A††pubyear: -

1 Introduction

A rotating magnetized neutron star is surrounded by a plasma-filled electrically conducting force-free magnetosphere in which the magnetic field is energetically dominant and governs its overall dynamics (Goldreich & Julian, 1969). Time-dependent electrodynamic and magnetohydrodynamic numerical simulations relax to a steady-state ideal force-free magnetic field configuration that corotates with the neutron star (Spitkovsky, 2006; Komissarov, 2006; Tchekhovskoy et al., 2013). This magnetospheric solution was first obtained by Contopoulos et al. (1999) in the case of an axisymmetric rotator.

An important characteristic of the steady-state configuration is that magnetic field lines that cross the light-cylinder contain a certain distribution of poloidal electric current that forms a large scale electric circuit. This is associated with a toroidal magnetic field component that reflects the fact that magnetic field lines are swept backwards with respect to the stellar rotation. Without it, plasma ‘frozen into’ these field lines beyond the light cylinder would move faster than the speed of light. This electric current distribution is the only one that guarantees smooth crossing of the light cylinder by the magnetic field, and in that sense, it is an ‘eigenfunction’ of the problem.

The ‘generator’ (or ‘battery’) of the magnetospheric electric circuit is the neutron star rotation, the ‘wires’ are the magnetic flux surfaces, and the ‘loads’ are finite dissipation regions at large distances (near and beyond the light cylinder, and the termination shock at very large distances). The electric current closes through the stellar interior. It penetrates deep inside the crust where it generates the torques necessary to spin down the neutron star (force-free conditions must be abandoned there).

The crust comprises an exceptionally strong ion lattice, nevertheless, it can only sustain finite stresses (Strohmayer, 1991; Chamel & Haensel, 2008; Horowitz et al., 2015). If the spin-down torque is exerted on a very thin volume, the Maxwell stresses could in principle exceed the yield limit of the crust, and the crust would yield. While the magnetospheric solution is obtained by assuming an ideal plasma with infinite conductivity, the crust has a high but finite conductivity $\sigma$ ranging between $10^{20}$ and $10^{27}$ s*-1* (Potekhin et al., 2015). We note that even if there is a finite resistivity in the magnetosphere (Li et al., 2012; Kalapotharakos et al., 2012), the big picture does not change qualitatively.

In this paper we obtain the flow of the magnetospheric electric current inside the neutron star crust and calculate the transfer of magnetospheric spindown torque into the stellar interior. The plan of the paper is as follows. In section 2, we derive the governing equations that describe the flow of electric current in the stellar interior. In section 3, we solve these equations numerically and present results for several pulsar models. We discuss their implications in section 4, and present our conclusions in section 5.

2 Problem setup

2.1 The equation for the electric current in the crust

Let us consider an axisymmetric stationary configuration. In what follows, we will work in spherical coordinates $(r,\theta,\phi)$ centered on the neutron star and aligned with the axis of symmetry (which coincides with the axis of rotation and the magnetic axis). The electric current density may in general be expressed as

[TABLE]

where $I=I(r,\theta)$ is the electric current that passes through a ring perpendicular to and concentric with the axis of symmetry passing through position $(r,\theta)$. Notice that $I$ is related to the toroidal component of the magnetic field ${\bf B}$ as

[TABLE]

and we further assume that the poloidal component of the magnetic field is a dipole. The azimuthal component of the electric current $j_{\phi}$ in eq. (1) is due to the corrotation of the internal space-charge, namely $j_{\phi}=r~{}\sin\theta~{}\Omega~{}\nabla\cdot{\bf E}$, and does not enter in our calculations below. $\Omega$ is the stellar angular velocity and the electric field is given by Ohm’s law

[TABLE]

Here $c$ is the speed of light and $\sigma$ the electric conductivity of the crust. The power per unit volume in the crust is given by:

[TABLE]

where ${\bf v}=r~{}\sin\theta~{}\Omega~{}{\bf{\bm{\hat{\phi}}}}$ is the velocity of the crust at $(r,\theta)$ for an observer in the lab frame. The term $\left({\bf j}\times{\bf B}\right)\cdot{\bf v}/c$ in eq. (4) expresses the work per unit volume and time done by the Lorentz force that spins-down the pulsar. The $j^{2}/\sigma$ term is the Ohmic thermal losses per unit volume and time, due to the finite conductivity of the crust.

Our goal is to obtain the distribution $I(r,\theta)$ in the stellar interior. We will approach this question by applying a Fermat-type principle. We propose that the current inside the crust of the neutron star will distribute itself so that it minimises the total Ohmic thermal losses. This allows us to formulate a minimisation equation.

[TABLE]

where we have used that ${\bf j}=(c/4\pi)\nabla\times{\bf B}$. By demanding that $P_{\rm Ohm}$ is minimized, we obtain the condition

[TABLE]

(see Appendix A). The same result is obtained if we start from Ohm’s law in the stellar interior (eq. 3) and realise that, for a stationary configuration, $\nabla\times{\bf E}={\bf 0}$ (Faraday’s law), this yields

[TABLE]

which is identical to eq. (6). With the help of eq. (1), eq. (6) then takes the form

[TABLE]

where we assume that the electric conductivity of the crust is a function of radius only $\sigma=\sigma(r)$.

The two approaches are interconnected. Equation 6 essentially describes an Ohmic eigenmode (Chanmugam & Gabriel, 1972) corresponding to the zero eigenvalue, or equivalently infinite decay time, subject to given boundary conditions. The infinite decay time is imposed here by setting $\nabla\times{\bf E}=0$. The Ohmic thermal power minimisation derivation, starts from Ohm’s law as well. Here we assume that among all possible electric current configurations that are compatible with the boundary conditions the one that will survive the longest is the one that has the lowest Ohmic thermal losses. We note here that while magnetic field energy is converted into heat as described by the term $j^{2}/\sigma$, this power is replenished by the fact that we enforce time-independent boundary conditions.

We note that the magnetospheric current closing through the crust is not related to the Hall current (Goldreich & Reisenegger, 1992). The latter is associated with the structure of the crustal magnetic field, and it can drive magnetic field evolution, especially for magnetic fields above $10^{14}$ G. Here we assume for simplicity that the magnetic field in the crust is a pure dipole, and that the source of the field (i.e. the associated azimuthal electric current) lies interior to the crust.

2.2 Boundary conditions

We will integrate eq. (8) in a computational domain $r_{\rm in}\leq r\leq r_{\rm out}$ and $0\leq\theta\leq\theta_{\rm out}$ of the crust. We first need to specify $I(r,\theta)$ at the boundaries of the domain. We set the outer radial boundary at the radius of the star, namely $r_{\rm out}=r_{\rm ns}$, and the inner radial boundary at the inner crust radius $r_{\rm in}=0.9r_{\rm ns}$. Along the axis, $I(r,0)=0$. $I(r_{\rm in},\theta)=0$ as we assume that the current is contained within the crust, (this is a reasonable assumption due to the high resistivity that prevents the current from penetrating below the crust). $I(r,\theta>\theta_{\rm out})=0$. This is justified by the fact that the magnetspheric poloidal current flows only along open field lines, and no current flows in the dead zone. Inside the crust, the majority of the current stays below the polar cap region, and does not spread to much lower latitudes. Actually, as we shall see, part of the current spreads beyond the polar region, but as long as the latitudinal boundary $\theta_{\rm out}$ is taken to be sufficiently large, its particular value is not important. For computational convenience we take $\theta_{\rm out}=2\theta_{\rm pc}$, where $\theta_{\rm pc}\approx(1.23\ r_{\rm ns}/r_{\rm lc})^{1/2}$ is the latitude of the footpoint of the last open field line on the star111We have also integrated eq. (8) with $\theta_{\rm out}=3\theta_{\rm pc}$ and the difference in the solution was minimal.. Here, $r_{\rm lc}=c/\Omega$ is the radius of the light cylinder. The prefactor $1.23$ in the above expression is based on the most detailed numerical solution of the axisymmetric problem to date, namely that of Timokhin (2006).

The distribution of electric current $I(r_{\rm out},\theta)$ along the surface of the neutron star is provided by the magnetospheric solution. This has been thoroughly investigated by several authors (Contopoulos et al., 1999; Gruzinov, 2005; Timokhin, 2006). Here we use the model with $x_{0}=0.992$ from Fig. 3 of Timokhin (2006). In that paper, $I$ is given as a function of $\Psi$, the poloidal magnetic flux, which on the surface of the neutron star is defined as $\Psi(\theta)=Br_{\rm ns}^{2}\sin^{2}\theta/2\pi$. Here and below, $B$ refers to the value of the magnetic field at the poles of the star. This allows us to obtain the distribution $I=I(r_{\rm out},\theta)$. Based on Timokhin’s solution, the maximum value of $I$ is $I_{\rm max}=0.87I_{\rm sm}$, where $I_{\rm sm}=1.23\times 2\pi Br_{\rm ns}^{3}r_{\rm lc}^{-2}$ corresponds to a split-monopole solution with the same amount of open magnetic flux (Michel, 1973). Therefore, the maximum value used in the present work is $I_{\rm max}=1.07\times 2\pi Br_{\rm ns}^{3}r_{\rm lc}^{-2}$. We implemented that by sampling the corresponding curve ανδ constructing a table with $\Psi~{}-~{}I$ pairs. Gralla et al. (2016) have shown that this solution is approximated to high accuracy by the polynomial expression

[TABLE]

where $\Psi_{0}=1.23$ to agree with the normalisation adopted above. In our approach we have solved the equation using expressions for the boundary condition, noting a deviation between the solutions of $\sim 1\%$.

Outside the polar cap, $I(r_{\rm out},\theta)$ drops to zero as a step function. We have smoothened this abrupt drop within a layer of width $0.01\theta_{\rm pc}$. This smooths out the current density near the surface of the star, but leaves the flow of current deeper in the crust mostly unaffected.

2.3 Neutron star parameters

We consider a neutron star radius $r_{\rm ns}=10$ km and we adopt a typical ground state structure for the crust as described in Chamel & Haensel (2008) (§ 3 see Figure 4). Atoms are fully ionized at mass densities higher than about $\rho\sim 10^{4}$ g cm*-3*. The so-called ocean extends up to densities of $10^{6}$ g cm*-3*. This forms a layer from about a few meters up to $100$ m, depending on the temperature of the neutron star (Potekhin et al., 2015). Below the ocean, the outer crust consists of a body-centered iron 56Fe cubic lattice with the composition of the nuclei becoming more neutron rich as a result of electron capture. The inner crust region extends from $\rho_{\rm nd}\sim 4\times 10^{11}$ to about $10^{14}$ g cm*-3*. At the bottom of the crust, some calculations predict various “pasta” phases of non-spherical nuclei, such as slabs or cylinders (Horowitz et al., 2015). Such pasta layer are believed to be highly resistive due to the anisotropic structure of the lattice and a low electron fraction (Pons et al., 2013). In the present work, we consider only the part of the crust ranging from the base of the ocean down to the crust-core boundary. This corresponds to densities ranging from $\rho_{\rm out}=1.3\times 10^{6}$ g cm*-3* to $\rho_{\rm in}=1.3\times 10^{14}$ g cm*-3* at the outer and inner crust boundary, respectively. Finally, we express the density of the crust as a function of the depth from the neutron star surface (Chamel & Haensel, 2008) with the following analytical expression

[TABLE]

The expression for the electric conductivity $\sigma(\rho)$ is taken from analytical fits of the numerical solutions obtained using the codes developed by Potekhin et al. (2015)222The codes are available at http://www.ioffe.ru/astro/conduct/index.html. We have chosen the following set of parameters: ion charge number (atomic number) Z$=26$, mass number A$=56$, impurity parameter Z${}_{\rm imp}=0.1$, range of densities $10^{6}\leq\rho\leq 10^{14}$ g cm*-3*, and temperature $T=10^{7}$ K. We have experimented with two magnetic field values, $B=10^{10}$ G and $10^{12}$ G. The differences between the two conductivity estimates are minimal (see the green and blue curves in Fig. 1), and we approximate them by the following analytical power-law expression

[TABLE]

with $\sigma_{\rm out}=2.5\times 10^{20}$ s*-1* (red line in Fig. 1). We note the small deviation between the expressions derived in Potekhin et al. (2015) and the analytical expression used here. We have verified that they have minimal impact on the electric current flow inside the crust (less than $1\%$ deviation). To assess the importance of the conductivity profile, we also integrated eq. (8) for a constant conductivity ($\sigma=10^{24}$ s*-1*) and we report the differences below.

3 Results

We solve eq. (8) with the Gauss-Seidel numerical method. The algorithm for this elliptic solver is provided in Numerical Recipes (Vetterling et al., 1988). We initialize the scheme with a trial distribution $I(r,\theta)$ and we repeat the iterative procedure until convergence is achieved. We implemented an $r-\theta$ numerical grid with a uniform resolution of $160\times 400$, and we find that the solution converges after $10^{6}$ iterations.

Based on the solutions that we have derived, we can evaluate the torque exerted on the pulsar, the corresponding stresses, and the Ohmic heating. In order to calculate these quantities, we assume that the internal magnetic field is a dipole, namely

[TABLE]

First, we calculate the Lorentz force per unit volume using eq. (1),

[TABLE]

Here we have to note that because of ${\bf E}=\left(r_{ns}/r_{lc}\right){\bf B}_{p}<<{\bf B}_{p}$, the electrostatic term $\rho_{e}{\bf E}$, which exists in eq. (13) is $\left(r_{ns}/r_{lc}\right)^{2}$ times smaller than the calculated one, so in limit of our numerical error is negligible.

Then, we calculate the torque per unit volume

[TABLE]

Because of axial symmetry only the torque component parallel to the axis of symmetry is non-zero. Thus, the total torque is given by the integral

[TABLE]

Finally, we calculate the components of the Maxwell stresses. The diagonal components correspond to pressure terms, while the off-diagonal components lead to the deformation of the crust due to shear-stresses. $M_{r\theta}$ involves the $B_{r}$ and $B_{\theta}$ components which are not due to the pulsar spin-down current.

[TABLE]

are due to the spin-down current which is associated with a toroidal magnetic field $B_{\phi}$ (see eq. 2). The breaking stress limit of the crust is

[TABLE]

(Chugunov & Horowitz, 2010), where $\Gamma=Z^{2}e^{2}/\left(\alpha k_{B}T\right)$ is the Coulomb coupling parameter, $\alpha=\left[3/\left(4\pi n_{\mathrm{i}}\right)\right]^{1/3}$ is the ion sphere radius, $n_{\mathrm{i}}$ is the ion number density, $k_{\mathrm{B}}$ is the Boltzmann constant and $e$ is the electron charge. A Maxwell stress comparable to $\tau_{br}$ may lead to crust yielding and deformation. We evaluate the breaking stress for densities in the range $\rho_{\rm nd}<\rho\leq\rho_{\rm in}$ using the table from Douchin & Haensel (2001) and for densities $\rho_{\rm out}\leq\rho\leq\rho_{\rm nd}$ using the results of Haensel & Pichon (1994). The breaking stress at the base of the crust is $\tau_{br}(\rho=10^{14}$g cm${}^{-3})=2\times 10^{29}$ dyn cm*-2*, at the neutron drip point $\tau_{br}(\rho=4\times 10^{11}$g cm${}^{-3})=1.4\times 10^{27}$ dyn cm*-2* and at the base of the ocean $\tau_{br}(\rho=10^{6}$ g cm${}^{-3})=10^{20}$ dyn cm*-2*. These results at the base of the crust and at the neutron drip point are consistent with the estimates of Cumming et al. (2004); Lander & Gourgouliatos (2019). The value near the upper boundary of our integration domain (the base of the ocean) depends on temperature through the Coulomb coupling parameter $\Gamma$.

The most important parameter of the problem is the pulsar period. This determines the size of the polar cap, and thus the boundary conditions on the surface. We have integrated eq. (8) for 5 different choices of the period $P=0.01,~{}0.1,~{}1,~{}5,~{}7.5$ s, thus exploring configurations that range from a rapidly rotating millisecond pulsar to a slowly rotating magnetar. The flow pattern of the electric current in the crust is independent of the strength of the magnetic field. Nevertheless, the physical quantities that we evaluate below depend on it. For this reason, we have assigned realistic values to the magnetic field ranging from $10^{10}$ to $10^{15}$ G to allow a direct comparison. The combinations employed are shown in detail in Table 1. For each magnetic model we integrate eq. (8) both for a constant and a variable conductivity given by eq. (11). We also report the deepest point reached by the electric current flow line that corresponds to $I(r,\theta)=I_{\rm max}/2$. This yields an estimate of the electric current penetration ’half-depth’, namely how deep $50\%$ of the current reaches inside the star. We evaluate the Ohmic power using eq. (5).

As a consistency test, we calculate the torque exerted on the star by integrating eq. (15), and comparing it with the spin-down torque of an aligned rotator in the force-free approximation appropriately corrected (Contopoulos & Spitkovsky, 2006)

[TABLE]

The correction factor of $0.94$ is due to the integral of the electromagnetic luminosity (Gralla et al., 2016). Once this is taken into account the volume integral of the spin-down torque $N_{tot}$ and $N_{\rm align}$ expression are for most models within a $2\%$ difference from each other. The results obtained using the boundary condition from the numerical solution of Timokhin (2006) and the polynomial fit of Gralla et al. (2016) give the same results within $2\%$ of each other, as expected since the two solutions agree at this level of accuracy.

4 Discussion

4.1 Electric current flowlines and Joule heating

In all models, there is a significant penetration of the electric current in the crust (Fig. 2). This is rather prominent in model A (millisecond pulsar) where the electric current practically reaches the base of the crust, with $50\%$ of the current reaching depths greater than $740$ m. Here the solution is affected by the boundary condition enforced at the base of the crust that does not allow the electric current to proceed any deeper. If this constraint were to be relaxed assuming the rest of the star had a similar conductivity, the current would formally reach into the core. In models B, C, D and SGR 1806$-$20, the current travels to a much smaller depth which scales with the radius of the polar cap. In these models the polar cap radii are smaller than the crust radius and the boundary condition at the base of the crust does not play any significant role. We remark further that a constant conductivity calculation yields a depth attained by the current approximately equal to $0.4$ times that of a realistic conductivity calculation (Fig. 3).

The paths of the electric current illustrate how the minimisation of Ohmic losses is achieved. Ohmic losses are larger for higher electric current densities, yet for lower ones, the same total current imposed on the boundary has to travel a longer distance inside the crust that eventually leads to a larger integration volume. Thus, if the conductivity is kept constant, the current will follow a path compromising these two effects. Once the conductivity varies with depth, the current will travel even deeper as this will allow it to cross a region of lower resistivity and thus suffer less Ohmic losses, despite the total path being longer. Joule heating is higher near the surface and decreases towards the base of the crust. This variation is more pronounced for the realistic conductivity profile, as there the conductivity increases by several orders of magnitude as one approaches the base of the crust. The maximum Joule heating occurs at the rim of the polar cap. This is because the bulk of the current enters the crust through this location leading to formally infinite current density. The total Joule heating scales with the magnetic field as $\propto B^{2}$, and with period as $\propto P^{-2.5}$.

We note that the total Ohmic losses occuring in the crust are negligible compared to the total radiated spin-down power (typically 10 orders of magnitude smaller). This implies that the coupling between the crust and the electric current is strong. Furthermore, while the coupling with the crust is essential for the pulsar spin-down, its effect on the global pulsar electric circuit is minimal.

4.2 Torque

As we saw in the previous section, the global torque approximates quite accurately the spin-down torque calculated through the magnetosphere. An interesting point here, is that there is a significant amount of localized spin-up torque that is of course overwhelmed by the spin-down torque. The reversal occurs along the surface $\partial I/\partial r=0$ where the current flow-lines become radial. The spin-up torque is mostly exerted on the part of the star where $\theta>\theta_{\rm pc}$, due to the spreading of the electric current at latitudes smaller than that of the polar cap (see Fig. 4).

The maximum torque per unit volume occurs in the region below the rim of the polar cap, and is due to the high electric current density there. In the rapidly spinning model (A1), the bulk of the spin-down torque is exerted close to the base of the crust. For slower rotating models, the bulk of the torque is exerted closer to the surface.

4.3 Maxwell Stresses

The distribution of the magnetospheric current of the force-free solution is such so that the current enters the star through the central part and most of the area of the polar cap and leaves the star through a narrow ring and mostly through a current sheet flowing on the separatrix between the open and closed magnetic field lines (here we have assumed that the magnetic moment and the angular momentum are parallel). Because of the singularity in the density of the current sheet, the Lorentz force and torque per unit volume become formally infinite at the edge of the polar cap. This is illustrated in Figures 2-6 by the convergence of the current flow lines at the edge of the polar cap. We note however that these are integrable singularities and the physical quantities associated to them (i.e. torque, net force) remain finite once we integrate over the corresponding volume. The possibility of crust yielding does not depend on the local value of the force density, but rather by comparing the Maxwell stress to the the breaking stress, i.e. eqs. (16) and (17). Indeed, Maxwell stresses remain finite and are a few orders of magnitude below $\tau_{br}$ for models A1 and C1. The crust does not yield either in model B1, but the ratio becomes $M_{r\phi}/\tau_{br}=0.2$ at the outermost layer of the integration domain. On the contrary, in model D1 the maximum shear stress becomes $M^{\rm max}_{r\phi}=5.1\times 10^{20}$ erg cm*-3* close to surface, which is high enough to cause crust yielding for a few meters ($\sim 10$ m) below the surface, see Fig. 5. Quite remarkably, the stress does not peak below the rim of the polar cap, where the current sheet enters the neutron star and the electric current density is the highest, but at some intermediate angle $\theta\approx 0.8\theta_{pc}$. This is because the Maxwell stress is proportional to $B_{\phi}=2I/(cr\sin\theta)$, which becomes maximum at some intermediate angle. The Maxwell stress is zero on the axis ($\theta=0$). This is because the $I(r_{rmns},\theta)$ becomes proportional to $\sin^{2}\theta$ as $\theta\to 0$ (Timokhin, 2006). Thus, if there is a part of the crust more likely to yield, this will be a ring of semi-opening angle $\approx 0.8\theta_{pc}$, rather than the region where the torque reverses from spin-down to spin-up, or even the axis.

The possibility of crust yielding and the maximum depth where this could occur depend strongly on the detailed physics of the outer crust. As the maximum stresses appear near the conventional surface of the neutron star ($\rho\sim 10^{6}$ g cm*-3*), this essentially lies at the interface between the ocean and the ion lattice. A hotter neutron star could have a deeper ocean. In practice, this implies that these stresses act on the fluid part of the crust, where eq. (17) is no longer applicable.

Pushing the question of crust yielding to the extreme, we have also considered the magnetar with the highest known magnetic field SGR 1806$-$20 (Woods et al., 2007), which has a long period $P=7.54$ s and thus a very small polar cap $\theta_{\rm pc}=0.33^{\rm o}$. Its inferred dipole magnetic field is $B=2\times 10^{15}$ G, and is the most prominent candidate for crust yielding. Integrating eq. (8) we find that the current reaches a depth of only 40 m beneath the surface. The spin-down torque $N_{\rm tot,~{}G}=8.0\times 10^{34}$ erg is in agreement within $2\%$ with $N_{\rm align}=8.2\times 10^{34}$ erg. We find that $M^{\rm max}_{r\phi}=1.1\times 10^{23}$ erg cm*-3* near the surface, which implies that the magnetospheric current will be extremely high to cause crust yielding to about 30 m below the surface (see Fig. 6). Such an event will be energetically unimportant compared to the energy that could potentially be released by magnetar activity. Nevertheless, it may impact the coupling between the magnetosphere and the crust, and therefore, the spin-down efficiency. This could be related to the higher timing irregularities that are observed in strongly magnetised neutron stars and magnetars (Hobbs et al., 2010). The timing noise in neutron stars with polar magnetic fields below $10^{12}$ G is independent of the magnetic field strength, whereas, in neutron stars with magnetic fields above this value it tends to increase and scale strongly with the magnetic field strength (Tsang & Gourgouliatos, 2013). It is conceivable that such behavior is related with crust yielding near the surface. Stronger magnetic fields lead to deeper crust failure. Given the episodic nature of crust failure (Thompson et al., 2017), the loss and recovery of the spin-down current coupling with the neutron star could manifest itself as erratic variations of the spin-down, i.e. practically as timing noise.

A possible consequence of the shallow penetration of the electric current and consequently the inefficient coupling between the spin-down current and the crust, can be the lack of isolated neutron stars with long periods, with a cut-off period in the range of $24$s across the entire pulsar population (Tan et al., 2018). This effect has been previously attributed to magnetic field decay (Pons et al., 2013), alignment between the magnetic axis (Johnston & Karastergiou, 2017) and observational selection effects (Faucher-Giguère & Kaspi, 2006). In the current picture, we note that a pulsar with a rotation period of $15$s will have a polar cap opening angle of $0.2^{\rm o}$ and the electric current will penetrate to a depth of $30$m which could lead to poor coupling and inefficient spin-down, therefore, it would be even harder for these pulsars to move to lower periods.

4.4 Twisted magnetospheres

In rotation powered radio pulsars, electric currents are associated with their spin-down. On the contrary, strongly magnetised neutron stars may have electric currents that support twisted magnetospheric structures (Beloborodov, 2009). According to this paradigm, electric current bundles form near the surface of the star and accelerate particles that bombard the surface of magnetars, thus generating X-ray emission. To determine the crustal current that supports such structures, one can follow the internal magnetic field evolution and solve self-consistently for the crust and the magnetosphere (Akgün et al., 2018), or alternatively consider an MHD equilibrium state taking into account the internal and external field (Glampedakis et al., 2014). Under the approach presented in this work, one can use the minimisation technique proposed to determine the minimum crustal electric current required to generate such a bundle. Using order-of-magnitude estimates, we find that the current supporting the bundle will be very much higher than the spin-down current, scaling approximately by a factor $(r_{\rm lc}/l_{\rm b})^{2}$, where $l_{\rm b}$ is the size of the bundle. Assuming that the bundle is comparable to the thickness of the crust, and considering a slowly spinning magnetar, this factor could be on the order of $10^{12}$. This would thus bring the magnetic energy dissipation rate in the range of $10^{32}$ erg s*-1*. This implies that a non-neglible fraction of the bundle energy may be dissipated inside the crust.

5 Conclusions

In this study we have explored the closure of the magnetospheric electric current through the neutron star crust in the simplest case of axisymmetry and steady-state. We have treated the crust and the magnetosphere as a global electric circuit, where the stellar rotation generates a poloidal electric current along the “infinitely conducting magnetic field wires” in the magnetosphere. This is very different from previous studies which proposed that the current penetrates only within a thin surface layer in a manner similar to the interaction of electromagnetic waves with the surface charges of a perfect conductor (Michel, 1991; Beskin et al., 1993; Beskin & Nokhrina, 2007). Skin-depth penetration refers to the interaction of an externally generated electromagnetic wave with a conductor. In our case, the magnetic field already penetrates deep into the stellar interior because the latter is the source of the magnetic field, and the stellar rotation is the generator (battery) of the large scale poloidal electric current. We must acknowledge, though, that, while the axisymmetric case offers important insight on the overall properties of this current, a more interesting and complicated situation arises once the three-dimensional magnetosphere is considererd. Indeed, if the magnetic and rotation axes are not aligned, calculating the electric current distribution in the stellar interior becomes highly non-trivial.

We also remark, that the adopted dipolar form for the magnetic field, while being the norm in models of pulsar magnetospheres, it could be a simplified picture of the realistic magnetic field structure. More complex magnetic fields are likely to be present in magnetars and even older neutron stars (Gourgouliatos & Hollerbach, 2018). If this is the case, the solution for the magnetosphere and consequently the crustal electric current would become more complicated (Gralla et al., 2017).

We have found that the magnetospheric current, responsible for the pulsar spin-down, enters the crust and reaches its base only if we consider rapidly rotating millisecond pulsars. In the case of slower pulsars, with periods longer than $1$ s, the bulk of the electric current reaches depths less than $100$ m. Even in the case of shallow current, the crust remains within its elastic limit without yielding, provided the magnetic field of the star is $10^{12}$ G or less. Ohmic losses in the crust are found to be orders of magnitude below the spin-down power.

While slower spinning neutron stars (old radio pulsars close to the death line, young magnetars) have rather small polar caps and the whole magnetospheric current closes through a narrow and shallow region of the crust, we find no stresses that exceed the elastic limit, except possibly in the outer few meters of the crust (as e.g. in SGR 1806$-$20). Energetically, such an effect may be insignificant compared to the X-ray power radiated by strongly magnetised neutron stars. Nevertheless, it affects the efficiency of the coupling and generates timing noise due to torque variations.

Acknowledgements

VK acknowledges COST ACTION PHAROS (CA16214) for STSM Grant 41713 that funded a visit to Durham University, and the Department of Mathematical Sciences of Durham University for hospitality and funding. KNG thanks Andrew Cumming and Anthony Yeates for useful discussions on the formulation of the problem and Jose Pons for suggesting the potential application to twisted magnetospheres of magnetars. We thank the referee of this paper for pointing us the correction factor of $0.94$ in eq. (18), and removing a systematic discrepancy of about $5\%$ between the torque evaluated through the integral of eq. (15) and analytically through eq. (18).

Appendix A Minimisation of Ohmic losses

Here we present the derivation equation 6 through a minimisation principle. Let

[TABLE]

and consider a variation ${\bf h}$ that vanishes at the boundary of $V$ so that ${\bf h}|_{\partial V}={\bf 0}$, as the magnetic field ${\bf B}$ is given on the boundaries of the domain. We then define

[TABLE]

$E({\bf B})$ will have a minimum, as $E>0$ if $V({\bf B},{\bf h})=0$. Let us further define

[TABLE]

and use from vector calculus the identity:

[TABLE]

Then equation 20 becomes:

[TABLE]

where we have used the divergence theorem. The second integral in the last equation is zero as ${\bf h}$ vanishes on the boundary of the domain. The first integral needs to be zero for any choice of ${\bf h}$. This is possible only if $\nabla\times{\bf A}={\bf 0}$ thus

[TABLE]

which is eq. (6).