The early life of millisecond magnetars

TL;DR

This paper investigates the early evolution of millisecond magnetars, focusing on how viscous dissipation and electromagnetic torques influence the alignment of magnetic and spin axes, which impacts gravitational wave and electromagnetic emissions.

Contribution

It introduces a combined model considering both viscous dissipation and electromagnetic torques affecting magnetic axis evolution in young millisecond magnetars.

Findings

Viscous dissipation tends to orthogonalize magnetic and spin axes.

Electromagnetic torques promote alignment of the axes.

An interplay exists between these effects influencing the star's evolution.

Abstract

Some neutron stars may be born spinning fast and with strong magnetic fields---the so-called \emph{millisecond magnetars}. It is important to understand how a star's magnetic axis moves with respect to the spin axis in the star's early life, as this effects both electromagnetic and gravitational wave emission. Previous studies have highlighted the importance of viscous dissipation within the star in this process. We advance this program by additionally considering the effect of the electromagnetic torque. We find an interesting interplay between the viscous dissipation, which makes the magnetic axis orthogonalise with respect to the spin, verses magnetic torques that tend to make the magnetic axis align with the spin axis. We present some results, and highlight areas where our model needs to be made more realistic.

Equations

$$\rho (r,\theta ,\varphi )={\rho }_0(r)+\delta {\rho }_B(r,\theta )+\delta {\rho }_{\alpha }(r,\theta ,\varphi ).$$

$$\omega =\alpha \cos \chi \frac{I_{zz}-I_{xx}}{I_0}\equiv \alpha {ϵ}_B\cos \chi ,$$

$${\mathbf{V}}_{\mathrm{inertial}}=\alpha {\mathbf{e}}_z^{(\alpha )}\times \mathbf{r}+\omega {\mathbf{e}}_z^{(B)}\times \mathbf{r}+\mathbf{\xi }$$

$$\frac{\partial \mathbf{v}}{\partial t}+(\mathbf{v}\cdot \nabla )\mathbf{v}+2\mathbf{\Omega }\times \mathbf{v}+\frac{\mathrm{d}\mathbf{\Omega }}{\mathrm{d}t}\times \mathbf{r}+\mathbf{\Omega }\times (\mathbf{\Omega }\times \mathbf{r})=-\nabla H-\nabla \Phi +\frac{1}{4\pi \rho }(\nabla \times \mathbf{B})\times \mathbf{B},$$

$$H(\mathbf{r})\equiv {\int }_{\hat{P}=0}^{\hat{P}=P(\mathbf{r})}\frac{\mathrm{d}\hat{P}}{\rho (\hat{P})},$$

$$\frac{\partial \rho }{\partial t}=-\nabla \cdot (\rho \mathbf{v}),$$

$${\nabla }^2\Phi =4\pi G\rho ,$$

$$\frac{\partial \mathbf{B}}{\partial t}=\nabla \times (\mathbf{v}\times \mathbf{B}),$$

$$H=H(\rho ),$$

$$\nabla \cdot \mathbf{B}=0.$$

$${ϵ}_{\alpha }\sim \frac{{\Omega }^2{R}_{\ast }^3}{GM}\sim 0.21\frac{R_6^3{f}_{\mathrm{kHz}}^2}{M_{1.4}},$$

$${ϵ}_B\sim \frac{B^2{R}_{\ast }^4}{G{M}^2}\sim 1.9\times 10^{-6}\frac{B_{15}^2{R}_6^4}{M_{1.4}^2},$$

$$\rho ={\rho }_0+\delta {\rho }_{\alpha }+\delta {\rho }_B+\delta {\rho }_{\alpha \alpha }+\delta {\rho }_{\alpha B}+\delta {\rho }_{BB}+\dots$$

$$\mathbf{B}={\mathbf{B}}_0+\delta {\mathbf{B}}_{\alpha }+\delta {\mathbf{B}}_B+\dots$$

$$0=-\nabla H_0-\nabla {\Phi }_0,$$

$${\nabla }^2{\Phi }_0=4\pi G{\rho }_0,$$

$$H_0=H_0({\rho }_0).$$

$${\rho }_0={\rho }_0(r)={\rho }_c\frac{R_{\ast }\sin \left(\frac{\pi r}{R_{\ast }}\right)}{\pi r}.$$

$$0=-\nabla \delta H_B-\nabla \delta {\Phi }_B+\frac{1}{4\pi {\rho }_0}(\nabla \times {\mathbf{B}}_0)\times {\mathbf{B}}_0,$$

$${\nabla }^2\delta {\Phi }_B=4\pi G\delta {\rho }_B,$$

$$\delta H_B=\frac{\mathrm{d}H}{\mathrm{d}\rho }\delta {\rho }_B,$$

$$\nabla \cdot {\mathbf{B}}_0=0.$$

$${\mathbf{B}}_0=B_{\varphi }{\mathbf{e}}_{\varphi }=\Lambda {\rho }_{0}r\sin \theta {\mathbf{e}}_{\varphi },$$

$${ϵ}_B=k_B\frac{{\Lambda }^2}{G}\approx -0.019\frac{{\Lambda }^2}{G},$$

$${\mathbf{\Omega }}_{\alpha }\times ({\mathbf{\Omega }}_{\alpha }\times \mathbf{r})=-\nabla \delta H_{\alpha }-\nabla \delta {\Phi }_{\alpha },$$

$${\nabla }^2\delta {\Phi }_{\alpha }=4\pi G\delta {\rho }_{\alpha },$$

$$\delta H_{\alpha }=\frac{\mathrm{d}H}{\mathrm{d}\rho }\delta {\rho }_{\alpha }.$$

$$\delta {\rho }_{\alpha }=-\frac{5\pi {\alpha }^2}{16G}{j}_2\left(\frac{\pi r}{R_{\ast }}\right)\left(\sin (2\chi )\sin (2\theta )\cos \lambda +{\sin }^2\chi {\sin }^2\theta \cos (2\lambda )\right),$$

$$\frac{\mathrm{d}{\mathbf{\Omega }}_{\alpha }}{\mathrm{d}t}\times \mathbf{r}+\omega {\mathbf{e}}_z\times ({\mathbf{\Omega }}_{\alpha }\times \mathbf{r})+{\mathbf{\Omega }}_{\alpha }\times (\omega {\mathbf{e}}_z\times \mathbf{r})=-\nabla \delta H_{\alpha B}-\nabla \delta {\Phi }_{\alpha B}-\frac{\delta {\rho }_{\alpha }}{4\pi {\rho }_0^2}(\nabla \times {\mathbf{B}}_0)\times {\mathbf{B}}_0$$

$$+\frac{1}{4\pi {\rho }_0}[(\nabla \times \delta {\mathbf{B}}_{\alpha })\times {\mathbf{B}}_0+(\nabla \times {\mathbf{B}}_0)\times \delta {\mathbf{B}}_{\alpha }],$$

Full text

I Introduction

In this paper we describe recent efforts to describe how the angle between the spin axis and the magnetic axis evolves, early in the life of a millisecond magnetar. It is based upon two recent papers Lander and Jones (2017, 2018), and also more recent Lander and Jones work.

This problem is of interest for several reasons. Firstly, millisecond magnetars have been considered as possible sources of gravitational waves, but only if the magnetic axis is misaligned with the spin axis (Dall’Osso, Shore, and Stella, 2009; Dall’Osso, Stella, and Palomba, 2018; Cutler, 2002). Secondly, there have been claims that the injection of energy from a rapidly spinning-down magnetar may partially power some gamma-ray bursts (Metzger et al., 2011). Finally, the results of our analysis could be compared against the observed distribution of inclination angles in the observed Galactic magnetar population. We do not discuss these interesting issues in what follows, but rather describe our attempts thus far to build a model of a spinning highly magnetised star that contains many (but not yet all) of the important ingredients in this problem.

II Precession of a magnetised star

We model our star as consisting of a spherically symmetric background configuration, with density $\rho_{0}(r)$, on which is superimposed two perturbations. One is due to the magnetic field, and produces a density perturbation $\delta\rho_{B}(r,\theta)$, while the other is due to rotation, and produces a density perturbation $\delta\rho_{\alpha}(r,\theta,\phi)$:

[TABLE]

Note that in writing down the density field in this way, we are using spherical polar coordinates orientated such that the $z$-axis coincides with the assumed axis of axisymmetry of the magnetic field.

The idea is that such a star will precess, if the magnetic and rotational and magnetic axes are misaligned, as described in some detail in Mestel and Takhar (1972). One can use axes that rotate with the star, such that the $z$-axis remains always along the magnetic axis, while the rotation axis is observed to move in a cone about the $z$-axis, at the (slow) precession frequency.

The effective ellipticity of such a star comes entirely from the magnetic deformation, such that this slow precession has frequency

[TABLE]

where $\epsilon_{B}$ is the magnetically-induced ellipticity, and $\chi$ is the angle between the (invariant) angular momentum axis and the magnetic axis. This angle is sometimes refereed to as the inclination angle, or, in this particular context, the wobble angle. Note that we only consider very young stars, before the onset of crust solidification, so we do not allow for any elastic component in the moment of inertia tensor.

It follows that the centrifugal bulge itself moves in a cone around the magnetic axis at the slow precession frequency. This is represented schematically in Figure 1. This means that the total velocity field of a fluid element is the sum of the motion one would expect for a classical rigid body, plus an extra piece due to this effective density wave propagating around the star. We can therefore write the velocity of a fluid element as:

[TABLE]

The classical rigid body motion is the sum of the rotation, at a rate $\alpha$, about the fixed angular momentum axis, plus the slow precessional motion at rate $\omega$ about the (instantaneous) location of this axis. The extra ‘non-rigid’ motions are encoded in the so-called $\xi$-motions, named after the symbol originally used to denote them in Mestel and Takhar (1972). It is these motions that we wish to compute, as they must play a key role in damping the precessional motion.

Mestel and collaborators noted that one can without too much difficult compute the density perturbation associated with the motion of the centrifugal bulge. Then one can write down the equation of mass conservation (equation (6) below) to try to compute the velocity field associated with these motions. This however only provides one equation for three unknowns (the three components of the displacement/velocity field). This was noted as a problem, and a suggestion made that beyond-linear-order perturbation theory might be relevant, but little progress was subsequently made in solving it.

The resolution to this problem is provided in (Lander and Jones, 2017). The key was to first assemble the full set of equations describing this magnetic precession problem. These consist of the Euler equation:

[TABLE]

the definition of the enthalpy:

[TABLE]

the equation of mass conservation:

[TABLE]

Poisson’s equation for the gravitational field:

[TABLE]

the time evolution equation for a magnetic field in a perfectly conducting fluid:

[TABLE]

an equation of state:

[TABLE]

and of course the solenoidal constraint on the magnetic field:

[TABLE]

In these equations, $\bf v$ is the fluid velocity, $\mathbf{\Omega}$ the rotation rate, $\bf r$ the position vector of the fluid element under consideration, $H$ the enthalpy, $\Phi$ the gravitational potential, $\bf B$ the magnetic field strength, and $\rho$ the density.

We would like a solution to these equations that represents our precessing magnetised star. This is clearly a formidable problem to solve. The key lies in recognising that there are two small parameters that we can exploit when using perturbation theory. One measures the rotational deformation of the star:

[TABLE]

where $f_{\rm kHz}=f/$kHz. The other small parameter is related to the size of the magnetic deformation:

[TABLE]

where $B_{15}=B/10^{15}$ G, $M_{1.4}=M/1.4M_{\odot}$, and $R_{6}=R/10^{6}$ cm. Clearly, these will both be small for realistic rotation rates and magnetic field strengths.

This immediately suggests a $2$-parameter perturbative scheme. In an obvious notation we can expand the various quantities that appear in our equations. For instance, for the density we can write:

[TABLE]

For the magnetic field we can write:

[TABLE]

We can then assemble coupled sets of equations, order by order.

To zeroth order we have:

[TABLE]

[TABLE]

[TABLE]

This simply represents a static spherical star. If we specialise to an $n=1$ polytrope for the equation of state, we have the well known density profile:

[TABLE]

To first order in $\epsilon_{B}$ we have:

[TABLE]

[TABLE]

[TABLE]

[TABLE]

This represents a static star deformed by a static magnetic field. Using techniques described in Lander and Jones (2009) one can find a solution:

[TABLE]

where the $\Lambda$ is a coefficient that can be determined for a given star and a given magnetic field configuration. We consider toroidal magnetic fields here, so as to give prolate deformations. For the field employed here we have

[TABLE]

where the coefficient that appears comes from a full numerical solution of the equations Lander and Jones (2009).

To first order in $\epsilon_{\alpha}$ we have:

[TABLE]

[TABLE]

[TABLE]

This represents a star deformed by rotation, with the complication that the rotation axis moves slowly relative to the star. At this level of approximation, the solution is simply a ‘snap-shot’, i.e. given by the density field appropriate to the instantaneous location of the spin axis relative to the star. Explicitly:

[TABLE]

where $\lambda=\phi+\omega t$.

If we work to order $\epsilon_{\alpha}\epsilon_{B}$ we have:

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

This constitutes a set of nine equations in the nine unknowns $(\delta\rho_{\alpha B},\delta H_{\alpha B},\delta\Phi_{\alpha B},{\boldmath\dot{\xi}},\delta{\bf B}_{\alpha}$). The set represents the combined effect of rotation and magnetic fields, and, crucially, allows computation of the all-important $\xi$-motions.

A method of solution is set out in detail in Lander and Jones (2017), where the solution to the problem is written in terms of spherical harmonics. A solution retaining harmonics up to and including $l=4$-terms is given in their section 7.1. A graphical representation is reproduced in Figure 2. It is this displacement field that can now be used to investigate dissipation of the precessional motion.

III Formalism for evolving the spin vector

We now need to write down the equations that will allow us to see how the precessional motion evolves in time, when electromagnetic torques and internal dissipation are added.

The magnetic axis spins around the angular momentum axis at a rate

[TABLE]

which can be thought of as the magnitude of the spin rate.

We know that the star’s energy will decrease, due to both electromagnetic torque, and viscous dissipation:

[TABLE]

while the angular momentum decreases only by virtue of the electromagnetic torques:

[TABLE]

The star’s energy will be a function of both the angular momentum and the wobble angle, so we can write:

[TABLE]

This can be rearranged to give us an evolution equation for $\chi$:

[TABLE]

If we consider spin-down by a vacuum dipole, one can show (see discussion in Cutler and Jones (2001) for the analogous gravitational wave torque case):

[TABLE]

We then obtain

[TABLE]

and

[TABLE]

which to our order of working may be rewritten

[TABLE]

We now have two (coupled) evolutions equations for the wobble angle $\chi$ and the spin rate $\Omega$, equations (40) and (42), which can then be solved if we specify the nature of the torques and the damping mechanisms.

IV Bulk viscosity

To account for viscous dissipation one needs to add terms involving shear viscosity (with coefficient $\eta$) and bulk viscosity (with coefficient $\zeta$) to the Euler equation, to give:

[TABLE]

where $\sigma_{ab}$ measures the amount of shear in the motion:

[TABLE]

see e.g. Lindblom and Owen (2002).

Back-of-the-envelope estimates readily show that in the regime of interest to us, bulk viscosity is much more important than shear, so we will neglect shear dissipation; see Lander and Jones (2018). Also, the damping timescales are long compared with the free precession period, so we can use the dissipation-free eigenfunctions computed above to estimate the rate of energy dissipation.

The bulk viscosity coefficient itself takes the form

[TABLE]

In this equation, $\omega$ is the (angular) frequency of the free precession, while $\tau$ is the timescale of the microphysical relaxation process that supplies the bulk viscosity. Note that at fixed $\omega$, this coefficient has a maximum when $\omega=1/\tau$, i.e. there is a ‘resonance’-like feature in the bulk viscosity magnitude.

We will consider dissipation due to modified Urca reactions which introduces a steep temperature dependence in $\tau$, as described in Reisenegger and Goldreich (1992):

[TABLE]

This means that the stellar temperature will be important when attempting to evolve the precessional motion. We will use a simple formula for cooling, also due to modified Urca reactions, as described in Page, Geppert, and Weber (2006):

[TABLE]

V The electromagnetic torque

We will consider the overly-simple case of vacuum dipole spin-down, for which:

[TABLE]

with $\lambda(\chi)$ a function that encodes how the torque varies with $\chi$, with $\lambda(\chi)=\sin^{2}\chi$ for the vacuum dipole considered here. Realistic torques will include a magnetospheric component, which will give non-zero spin-down even for aligned $(\chi=0)$ stars. We will discuss this point in the Section ‘Future Work’ below.

VI Time evolutions

We have now assembled all the ingredients needed to evolve the quantities $\Omega,\chi$, with the temperature evolution being given by equation (46). Before presenting the results of such evolutions, we can gain a little insight as follows.

The electromagnetic torque tends to make the star align. For a prolate star, the internal dissipation tends to make it orthogonalise. It follows that these two effects will balance when their respective timescales of operation are equal:

[TABLE]

where we have defined a timescale of the form

[TABLE]

By combing equations given above, it can be shown that this occurs when the spin frequency takes a critical value:

[TABLE]

where for convenience we have defined

[TABLE]

where $g(\chi)$ is a dimensionless function that encodes how the rate of viscous dissipation varies with $\chi$. For small $\chi$, it can be shown that $g(\chi)\sim\sin^{2}\chi\sim\chi^{2}$; see Lander and Jones (2018) for details. Stars with spin frequencies above this tend to orthogonalise, those with frequencies below this tend to align. This equation therefore defines a critical curve for orthogonalisation in the spin frequency–temperature plane. A family of such curves, for three different magnetic field strengths, is given in Figure 3, where the $\chi$-dependent factors have been set equal to unity.

One can use such curves to anticipate how a star may evolve.

Consider the possible trajectories sketched, again schematically, in Figure 4. Stars will be born hot and spinning fast, somewhere in the top right of the plot, marked by position 1. The star will them cool rapidly, moving left in the plot, to position 2. If the star spins down significantly it may fail to enter the orthogonalisation region, and so never othogonalise. This is the case for the trajectory that ends at the point 3a, relevant for a high field strength orthogonalisation curve. However, for most weakly magnetised stars, the larger dotted orthogonalisation curve applies. In this case the star may enter and go through a period of orthogonalisation, denoted by 3b in the plot. It will eventually cool through the orthogonalisation region, exiting to evolve to 3c, where the lower temperatures and slower rotation rate will make subsequent evolution slow. Such a star may be driven to orthogonality.

The results of actual time evolutions are shown in Figure 5. This shows the phase space of birth spin frequency and magnetic field strength. The colour scale is used to indicate the nature of the stellar configuration after various times since birth. The blue colour scale shows the time taken (in log seconds) for states to align, for those stars that do indeed align. Orange is used to show the time taken for stars to orthogonalise. Note that we evolve to timescales of order a million years. However, we only trust our model prior to the formation of the crust, which occurs on a timescale of a hundred years or less. We show such long term evolutions simply to make clearer the nature of the solutions.

We can summarise the main features as follows. Below a field strength of about $10^{14}$ G, all stars orthogonalise, regardless of birth spin frequency. This corresponds to the orthogonalisation region of Figure 4 being large, so that all stars enter and undergo a significant period of orthogonalisation. Above $10^{14}$ G, there is a (field-dependent) spin frequency below which stars align, and above which star orthogonalise. This can be understood again in terms of Figure 4: only those stars that are born spinning fast enough can enter the orthogonalisation region, and undergo a significant period of orthogonalisation.

VII Future work

The above results illustrate the interplay of internal dissipation and electromagnetic torques that play key roles in the early life of a millisecond magnetar. However, the model is highly idealised, and a number of issues will need to be addressed before firm conclusions can be drawn. We will now discuss the most important outstanding issues.

The form of the spin-down torque. The model and results above assume a vacuum dipole model for the electromagnetic spin-down torque. However, it has been known for a long time that a significant fraction of the star’s angular momentum is likely to be carried away by particles accelerated within the pulsar’s magnetosphere (Goldreich and Julian, 1969). In fact, in the paper Lander and Jones (2018) we presented results for such a torque, with angular dependence $\lambda(\chi)=1+\sin^{2}\chi$. Such a torque is qualitatively different from the vacuum dipole one, as in the limit $\chi\rightarrow 0$ the vacuum torque goes to zero, while the magnetospheric one does not. However, there was an inconsistency in our previous treatment, as we employed equation (39),which is valid for the vacuum dipole torque, but not for the magnetospheric one. We aim to correct this in a future publication.

Neglect of buoyancy. In perturbing the Euler equation above, we have implicitly assumed that the perturbations obey the same equation of state as the background (unperturbed) star. In general this will not be the case. In particular, when the reactions responsible for the bulk viscosity are slow compared with the mode period, we might expect the eigenfunctions themselves to be different from those calculated in Lander and Jones (2017), with a smaller divergence $\nabla\cdot\xi$, as described in Lasky and Glampedakis (2016). This could be modelled by shutting off the effects of bulk viscosity when the reaction timescales becomes sufficiently long, as was done in Lasky and Glampedakis (2016).

Opacity to neutrinos. When the star is sufficiently hot, it will be opaque to neutrinos. This could be accounted for in future work using the prescription of Lasky and Glampedakis (2016), who switched bulk viscosity off above a critical temperature.

Evolution of the temperature. Our prescription for evolving the temperature was particularly simple. In reality, additional cooling effects may be important in the early life of magnetars, as described in Metzger et al. (2011).

We intend addressing all of these issues in future work (Lander & Jones, in preparation). In the meantime, we feel the scheme here, with its competition between electromagnetic alignment and dissipative orthogonalisation, has taken an important step forward in modelling the early life of a millisecond magnetar.

VIII Acknowledgement

DIJ acknowledges funding from STFC through grant number ST/M000931/1.