A possible way to reconcile long-period precession with vortex pinning in neutron stars

TL;DR

This paper proposes a model where neutron star core vortex pinning and long-period precession coexist by assuming the core can rotate independently, addressing a key inconsistency in pulsar glitch theories.

Contribution

It introduces a three-component model allowing vortex pinning and precession to coexist, with observational constraints on interaction coefficients.

Findings

Vortex pinning can occur in the magnetically decoupled core region.

Long-period precession and vortex pinning are compatible in the same pulsar.

Constraints on interaction coefficients are derived from observational data.

Abstract

We propose a possible way to solve the problem of inconsistency between the neutron star long-period precession and superfluid vortex pinning, which is the basis of the most successful theories of pulsar glitches. We assume that the pinning takes place in the region of the neutron star core, which, being magnetically decoupled, can rotate relative to the crust. In the framework of a simple three-component model we show that these two phenomena can coexist in the same pulsar. Some constraints on the formally introduced interaction coefficients following from observation data are formulated.

Equations

$$d_t\mathbf{M}=\mathbf{K},$$

$$\mathbf{M}=I{\mathbf{\Omega }}_c+Iϵ({\mathbf{\Omega }}_c\cdot {\mathbf{e}}_c){\mathbf{e}}_c,$$

$$\mathbf{K}={\mathbf{K}}_2+{\mathbf{K}}_3,$$

$${\mathbf{K}}_2=-\delta I_m(\mathbf{\Omega }\cdot {\mathbf{e}}_m)[\mathbf{\Omega }\times {\mathbf{e}}_m],$$

$$\delta I_m\sim 10^{32}{r}_6^5{B}_{12}^2\text{mbox}{gcm}^2,$$

$${ϵ}_{\mathrm{new}}={ϵ}_{\mathrm{old}}+\delta I_m/I,$$

$$d_t\mathbf{M}=0.$$

$${\dot{\mathbf{\Omega }}}_c=ϵ({\mathbf{e}}_c\cdot {\mathbf{\Omega }}_c)[{\mathbf{e}}_c\times {\mathbf{\Omega }}_c].$$

$${\omega }_p=ϵ\cos \theta {\Omega }_c,$$

$$T_p=\frac{P}{ϵ\cos \theta }.$$

$${\mathbf{\Omega }}_c={\mathbf{\Omega }}_0=\mathrm{const},{\mathbf{\Omega }}_0\parallel {\mathbf{e}}_c.$$

$${\mathbf{\mu }}_c={\mathbf{\Omega }}_c-{\mathbf{\Omega }}_0.$$

$$V^{\parallel }=(\mathbf{V}\cdot {\mathbf{e}}_z),V^{\perp }=(\mathbf{V}\cdot {\mathbf{e}}_x)+i(\mathbf{V}\cdot {\mathbf{e}}_y),$$

$${\dot{\mu }}_c^{\parallel }=0,$$

$${\dot{\mu }}_c^{\perp }=iϵ{\Omega }_0{\mu }_c^{\perp }.$$

$$p=iϵ{\Omega }_0.$$

$$ϵ\sim \frac{P}{T_p}.$$

$$\mathbf{M}=\tilde{I}{\mathbf{\Omega }}_c+ϵ\tilde{I}({\mathbf{\Omega }}_c\cdot {\mathbf{e}}_c){\mathbf{e}}_c+{\mathbf{L}}_{\mathrm{sf}}.$$

$${\dot{\mathbf{L}}}_{\mathrm{sf}}=0.$$

$${\dot{\mathbf{\Omega }}}_c=\left(ϵ+\frac{L_{\mathrm{sf}}}{\tilde{I}{\Omega }_c\cos \theta }\right){\Omega }_c\cos \theta [{\mathbf{e}}_c\times {\mathbf{\Omega }}_c].$$

$${\dot{\mu }}_c^{\parallel }=0,$$

$${\dot{\mu }}_c^{\perp }=i\widetilde{ϵ}{\Omega }_0{\mu }_c^{\perp },$$

$$\widetilde{ϵ}=ϵ+\frac{L_{\mathrm{sf}}}{\tilde{I}{\Omega }_0}\approx I_{\mathrm{pin}}/\tilde{I}\gtrsim 10^{-2}.$$

$$p=i\widetilde{ϵ}{\Omega }_0$$

$$T_p\lesssim \frac{P}{\widetilde{ϵ}}\sim 10^{2}P.$$

$$d_t{\mathbf{M}}_c={\mathbf{N}}_{gc},$$

$$d_t{\mathbf{M}}_g={\mathbf{N}}_{cg},$$

$${\mathbf{M}}_c=I_c{\mathbf{\Omega }}_c+I_c{ϵ}_c({\mathbf{\Omega }}_c\cdot {\mathbf{e}}_c){\mathbf{e}}_c+{\mathbf{L}}_{\mathrm{sf}}.$$

$${\mathbf{M}}_g=I_g{\mathbf{\Omega }}_g.$$

$$d_{t}E=-({\mathbf{\Omega }}_c-{\mathbf{\Omega }}_g)\cdot {\mathbf{N}}_{cg}.$$

$${\dot{\mathbf{\Omega }}}_c+\left({ϵ}_c+\frac{{\omega }_{\mathrm{sf}}}{\Omega \cos \theta }\right)({\mathbf{\Omega }}_c\cdot {\mathbf{e}}_c)[\mathbf{\Omega }\times {\mathbf{e}}_c]$$

$$={\mathbf{R}}_{gc}-\frac{{ϵ}_c}{1+{ϵ}_c}({\mathbf{R}}_{gc}\cdot {\mathbf{e}}_c){\mathbf{e}}_c,$$

Full text

††pubyear: 2019††pagerange: A possible way to reconcile long-period precession with vortex pinning in neutron stars–A possible way to reconcile long-period precession with vortex pinning in neutron stars

1 Introduction

Radio pulsars are the sources of exceptionally stable pulse sequences with very slowly decreasing frequency. Some pulsars, besides smooth slow-down, demonstrate so called glitches, which are sudden increases of the frequency and its time derivative followed by smooth recovery towards pre-glitch values. The relative amplitude of the glitches $\Delta\Omega/\Omega$ lies in the range of $\sim 10^{-12}-10^{-5}$ (Espinoza et al., 2011)111see also http://www.jb.man.ac.uk/pulsar/glitches.html while the recovery time-scale is of the order of days to months (Lyne et al., 2000). The long time-scales of the glitch recovery indicate that the nucleon superfluidity supposed to be present in neutron star interiors participates in this phenomenon (Baym et al., 1969). The most successful theories of glitches are based on the assumption that the neutron superfluid vortices are pinned in some region of the neutron star interior (Haskell & Melatos, 2015). The pinned superfluid conserves angular momentum and, from time to time, releases part of it, spinning up the outer component of the star.

There are several pieces of evidence that isolated neutron stars can precess with long periods. Pulsar B1821-11 demonstrates correlated periodic variations of spin-down rate and beam shape (Stairs et al., 2000). The most favorable explanation for this is the precession of the neutron star with period $T_{p}\approx 500$ d (Ashton et al., 2016). Several pulsars show periodic variations of spin-down rate without significant correlation with beam shape (Kerr et al., 2016). The time-scales of the variations are 0.5 – 1.5 yr. Some pulsars switch between two spin down states. Despite the fact that the switching itself is a fast process ( < 1 min), there are at least six pulsars for which the probability of being in a particular state is a quasi-periodic function with a characteristic time-scale $\sim 10^{2}$ days (Lyne et al., 2010). Hence, it is natural to assume that this probability depends on the precession phase (Jones, 2012).

Besides the directly observed variations, there are observational data, which can also be interpreted as manifestations of neutron star precession but with much longer periods. Lyne et al. (2013) found a steady increase in separation between the main pulse and the interpulse of Crab pulsar at 0.62${}^{\circ}\pm 0.03^{\circ}$ per century. The authors concluded that this is a consequence of the increase of the pulsar inclination angle. The rate of increase seems to be too large for secular evolution but it can be ensured by free precession with period $\sim 10^{2}$ yr (Arzamasskiy et al., 2015). Biryukov et al. (2012) argued that anomalously large braking indices indicate that the pulsar spin-down rate can oscillate at the time-scale of $10^{3}-10^{4}$ years. This oscillation can be caused by precession with corresponding periods. The stellar magnetic field by itself should make neutron stars precess at such time-scales (Melatos, 2000).

The problem is that it is difficult for neutron vortex pinning and long-period ($\gtrsim$1 year) precession to coexist in the same neutron star. As it was first pointed out by Shaham (1977), the pinning of superfluid vortices dramatically decreases the period of precession. Thus, on the one hand, the glitch theories require about 1 percent of the total stellar moment of inertia to be contained in a pinned superfluid (Andersson et al., 2012). On the other hand, if a neutron star precess with periods $\sim$1 yr or longer, only part of $<10^{-8}$ of the total moment of inertia can be in the pinned superfluid. Shaham considered a simple model with perfectly pinned vortices. Several attempts to attack this problem with more detailed models of pinned vortex dynamics (Sedrakian et al., 1999; Link & Cutler, 2002; Alpar, 2005) or assuming that pinning is absent (Link, 2006; Kitiashvili & Gusev, 2008) have been made. However, to date the problem remains current (Jones et al., 2017).

In the present paper we suggest a way that may allow these two phenomena to be to reconciled. The paper is organized as follows. In Section 2 we introduce basic notations and formulate the problem. In Section 3 we perform a linear mode analysis for a two-component neutron star with the superfluid pinned in the crust. Section 4 contains the same analysis but for a system with the superfluid located in the core.

In Section 5 we formulate a full three-component model and perform the linear mode analysis for it. Section 6 is devoted to the derivation of quasi-stationary equations, allowing the long-term evolution of the precessing neutron stars to be studied. In Section 7 we study the ability of the model to show glitch-like events and what restrictions for formally introduced coefficients can be obtained from the glitch observation data. In Discussion we speculate on the possible physical realization of the model, discuss some observational data and draw conclusions.

2 Formulation of the problem

Let us consider a rotating neutron star. From the point of view of an observer in an inertial frame of reference the following equation should be satisfied:

[TABLE]

where $\boldsymbol{M}$ is the stellar angular momentum and $\boldsymbol{K}$ is external electromagnetic torque.

We start with the the “rigid-body” approximation assuming that the whole star rotates with a uniform angular velocity that we denote by $\boldsymbol{\Omega}_{c}$. The index “$c$” is introduced for the sake of notation unification with the subsequent sections of this paper. We assume that the star has an axisymmetric shape maintained by the crust rigidity (Cutler et al., 2003) or by the magnetic field (Wasserman, 2003). In this case, the angular momentum can be represented in the following form:

[TABLE]

where $I$ is the star’s moment of inertia, $\boldsymbol{e}_{c}$ is the unit vector directed along the star’s symmetry axis and $\epsilon$ is the oblateness parameter.

To date, there is no consensus concerning the exact form of external torque $\boldsymbol{K}$ acting on isolated neutrons stars (Beskin et al., 2015). However, it is generally accepted that it consists of two parts of different nature (Davis & Goldstein, 1970; Good & Ng, 1985; Melatos, 2000):

[TABLE]

where $\boldsymbol{K}_{2}\propto\Omega_{c}^{2}$ is the so-called anomalous torque originating from the inertia of the near-zone electromagnetic field (Good & Ng, 1985; Beskin & Zheltoukhov, 2014; Goglichidze et al., 2015) while $\boldsymbol{K}_{3}\propto\Omega_{c}^{3}$ is related to the angular momentum transfer away from the star to infinity by both particles and electromagnetic radiation (Davis & Goldstein, 1970; Jones, 1976; Beskin et al., 1983; Barsukov et al., 2009; Philippov et al., 2014).

Relation ${K}_{3}/{K}_{2}$ can be estimated as $\sim(r_{\mathrm{ns}}\Omega_{c}/c)$, i.e. it is small for the most pulsars (except the fastest millisecond ones, which are not considered in the present paper) We will ignore term $\boldsymbol{K}_{3}$ in Sections 2-5 restricting ourselves to the consideration of time-scales much smaller than $\tau_{x}\sim\Omega_{c}/K_{3}$. The effects of this term are discussed in Section 6.

In the case of an axisymmetric stellar magnetic field, anomalous torque $\boldsymbol{K}_{2}$ can be represented in the following form (Goglichidze et al., 2015):

[TABLE]

where $\boldsymbol{e}_{m}$ is the unit vector directed along the field symmetry axis,

[TABLE]

$r_{6}$ is the neutron star radius in units of $10^{6}$ cm and $B_{12}$ is the surface magnetic field in units of $10^{12}$G. In the special situation of vector $\boldsymbol{e}_{m}$ coinciding with vector $\boldsymbol{e}_{c}$, the anomalous torque can be taken into account just by redefining the oblateness parameter (Zanazzi & Lai, 2015; Goglichidze et al., 2015):

[TABLE]

where $\epsilon_{\mathrm{old}}$ characterizes the stellar matter deformarion only. In general case, the anomalous torque can be taken into account by modification of star moment of inertia tensor.

Ignoring $\boldsymbol{K}_{3}$ and including $\boldsymbol{K}_{2}$ in $\boldsymbol{M}$ (without changing the notation), we can reduce equation (1) to

[TABLE]

In Sections 2-5 we will consider only the axisymmetric case. Therefore, formula (2) is assumed to be valid, where the oblateness parameter is now calculated with formula (6). The more general case of a triaxial star is discussed in Section 6.

Equation (7) formally states the conservation of total angular momentum. It is easy to verify that the mechanical energy $E=(\boldsymbol{\Omega}\cdot\boldsymbol{M})/2$ is also conserved. It is natural because we neglect torque $\boldsymbol{K}_{3}$ describing the energy and angular momentum flows away from the star. The effects of this term are discussed in Section 6.

Substituting angular momentum expression (2) into equation (7) and changing to the co-rotating frame of reference, we obtain

[TABLE]

Here and further we will denote by a superscript dot the time derivative in the frame of reference rotating with angular velocity $\boldsymbol{\Omega}_{c}$. The solution to this equation is the uniform rotation of vector $\boldsymbol{\Omega}_{c}$ about symmetry axis $\boldsymbol{e}_{c}$ with angular frequency

[TABLE]

where $\theta$ is the angle between angular velocity $\boldsymbol{\Omega}_{c}$ and symmetry axis $\boldsymbol{e}_{c}$. This type of motion is called free precession.

Introducing pulsar rotational period $P={2\pi}/{\Omega_{c}}$ and period of pulsar precession $T_{p}={2\pi}/{\omega_{p}}$, directly from (9) we have

[TABLE]

Equation (8) has a fixed point corresponding to stationary rotation about the symmetry axis:

[TABLE]

In order to study perturbations to the stationary rotation let us introduce a departure vector:

[TABLE]

Further, it will be convenient to use the following notations, which can be applied to an arbitrary vector $\boldsymbol{V}$:

[TABLE]

where $(\boldsymbol{e}_{x},\boldsymbol{e}_{y},\boldsymbol{e}_{z})$ is an orthonormal basis fixed in co-rotating frame of reference such that $\boldsymbol{e}_{z}=\boldsymbol{e}_{c}$. Substituting (12) into (8) and neglecting the quadratic in ${\mu}_{c}^{\parallel}$ and ${\mu}_{c}^{\perp}$ terms, we obtain

[TABLE]

Free precession motion is described by the “perpendicular” equation (15). Substituting $\mu_{c}^{\perp}\propto\exp(pt)$ we find

[TABLE]

Hence, the perpendicular part of vector $\boldsymbol{\mu}$ rotates about the symmetry axis with angular frequency $\omega_{p}=\epsilon\Omega_{0}$, which in the linear approximation ($\cos\theta\approx 1$, $\Omega_{c}\approx\Omega_{0}$) coincides with expression (9). Thus, knowing pulsar period $P$ and period of pulsar precession $T_{p}$, one can estimate the oblateness parameter as

[TABLE]

As it was first pointed out by Shaham (1977), the situation dramatically changes if a neutron star contains a pinned superfluid. In this case, the stellar angular momentum expression (2) requires some modification, namely,

[TABLE]

where $\boldsymbol{L}_{\mathrm{sf}}$ is the angular momentum of the pinned superfluid and $\tilde{I}$ is the moment of inertia of the rest of the stellar matter (excluding the pinned superfluid).

The angular momentum of the superfluid is totally determined by the distribution of quantum vortices. If the vortices are pinned, their density and orientation are fixed in the frame of reference rotating with angular velocity $\boldsymbol{\Omega}_{c}$. Therefore, angular momentum $\boldsymbol{L}_{\mathrm{sf}}$ is constant in the co-rotating frame of reference:

[TABLE]

For the sake of simplicity we restrict ourselves to the consideration of a particular configuration of $\boldsymbol{L}_{\mathrm{sf}}\parallel\boldsymbol{e}_{c}$. Such a configuration is not unrealistic but it is far from the general case. The analysis of star rotation with arbitrarily oriented $\boldsymbol{L}_{\mathrm{sf}}$ was given by Shaham (1977). It shows no significantly different results.

Substituting expressions (18) and (19) into equation (7) and changing to the co-rotating frame of reference, we obtain

[TABLE]

At first sight this equation is similar to equation (8). However, the rate of vector $\boldsymbol{\Omega}_{c}$ rotation is now determined by the sum of the two terms in parentheses. The second term can be estimated as $\sim I_{\mathrm{pin}}/\tilde{I}$, where $I_{\mathrm{pin}}$ is the moment of inertia of the pinned superfluid. If the understanding of glitches as a manifestation of neutron vortex dynamics is correct, at least about 1 percent of star’s moment of inertia should be contained in the pinned superfluid ($I_{\mathrm{pin}}/\tilde{I}\gtrsim 10^{-2}$) to ensure the observed glitch activity of the pulsars (Andersson et al., 2012). As for the first term, $\epsilon$, it hardly can reach such large magnitudes. Calculations give values at least several orders of magnitude smaller (Cutler et al., 2003; Wasserman, 2003). Hence, the second term in parentheses is dominant.

After linearization we obtain

[TABLE]

where we have introduced the effective oblateness parameter

[TABLE]

Substituting $\mu_{c}^{\perp}\propto\exp(pt)$ into equation (22), we find

[TABLE]

Thus, one can see that the free precession is governed almost exclusively by the amount of pinned superfluid which does not allow the star to precess with periods substantially longer than

[TABLE]

Such short periods are many orders of magnitude smaller than the observed values.

3 Two-component system

In the previous section we assumed that the whole star except the pinned superfluid rotates as a rigid body. Let us relax this assumption. First we consider the simplest case supposing that the neutron star consists of two dynamically distinguished components: an outer c-component whose rotation is directly observed and an inner g-component (Shaham, 1977; Sedrakian et al., 1999). At this point, the components are introduced phenomenologically. The indices “c” and “g” are chosen for the sake of unification with the full model presented in Section 5. For the components we can write

[TABLE]

where $\boldsymbol{M}_{i}$ are the components angular momenta and $\boldsymbol{N}_{ij}$ are the components interaction torques, $i,j=c,g$. The c-component consists of the “normal" fraction and the pinned superfluid fraction. By the “normal" fraction we mean all the possible constituents of the c-component except the pinned neutron superfluid. It includes the ions in the crust lattice sites, electrons and possible non-superfluid neutrons. It could also include the part of stellar core matter that is rigidly connected to the crust. The “normal" fraction is assumed to be axisymmetric with moment of inertia $I_{c}$ and oblateness parameter $\epsilon_{c}$. It rotates with angular velocity $\boldsymbol{\Omega}_{c}$. The pinned superfluid contains angular momentum $\boldsymbol{L}_{\mathrm{sf}}$ directed along the “normal" fraction symmetry axis $\boldsymbol{e}_{c}$. Thus, the total c-component angular momentum has the form

[TABLE]

The g-component is the part of the core that is not rigidly connected to the crust. This component, for simplicity, is assumed to be spherically symmetric. It is characterized by moment of inertia $I_{g}$ and rotates with angular velocity $\boldsymbol{\Omega}_{g}$. Therefore, its angular momentum is equal to

[TABLE]

Since the total stellar angular momentum $\boldsymbol{M}=\boldsymbol{M}_{c}+\boldsymbol{M}_{g}$ should be conserved, $\boldsymbol{N}_{cg}=-\boldsymbol{N}_{gc}$. However, for the mechanical energy we have

[TABLE]

Therefore, the non-zero components of the interaction torque parallel to angular velocity difference $(\boldsymbol{\Omega}_{c}-\boldsymbol{\Omega}_{g})$ lead to the dissipation of mechanical energy inside the star.

Substituting expressions (28) and (29) into equations (26) and (27), changing to the frame of reference co-rotating with the c-component, and taking into account equation (19), we obtain the equations for angular velocities:

[TABLE]

where the following notations have been introduced: $\boldsymbol{R}_{ij}=\boldsymbol{N}_{ij}/I_{j}$ and $\boldsymbol{\omega}_{\mathrm{sf}}=\boldsymbol{L}_{\mathrm{sf}}/I_{c}$. Note that $\boldsymbol{\omega}_{\mathrm{sf}}$ is just a convenient notation. In the general case, it is not have the meaning of superfluid angular velocity.

This system of equations has a fixed point corresponding to the uniform rotation of both components with the same angular velocity $\boldsymbol{\Omega}_{0}$ parallel to $\boldsymbol{e}_{c}$. In order to consider the linear perturbations to the equilibrium state, it is convenient to introduce small departure vectors:

[TABLE]

Let us also assume that, with linear in $\boldsymbol{\mu_{i}}$ accuracy, the interaction torques have a simple frictional form:

[TABLE]

where coefficients $\alpha_{ij}>0$. From the angular momentum conservation it follows that

[TABLE]

Substituting expressions (33) and (34) into equations (31), (32) and neglecting the quadratically small terms, we obtain

[TABLE]

where

[TABLE]

is the c-component effective oblateness parameter. One can see that the system of equations falls into two independent parts. The “parallel” part describes the evolution of the difference in absolute values of the angular velocities while the “perpendicular” part works when the angular velocities decline from the symmetry axis. Substituting $\mu_{i}^{\parallel}\propto\exp(pt)$ into equations (36)–(37), we obtain

[TABLE]

The angular velocity difference, being excited, decays exponentially at the characteristic time-scale $\sim p_{\parallel}^{-1}$.

Since our main goal is studying the precession, the “perpendicular" modes are more interesting for us. Substituting $\mu_{i}^{\perp}\propto\exp(pt)$ into equations (38)-(39), we obtain the following characteristic equation:

[TABLE]

This equation can be solved exactly. However, it would be more informative to obtain the approximate roots corresponding to different limiting cases.

Introducing an auxiliary interaction parameter,

[TABLE]

we can rewrite equation (42) in the following form:

[TABLE]

where

[TABLE]

is the effective oblateness parameter of the star as a whole.

If the interaction between the components is weak ($\sigma\ll 1$), the roots of equation (44) can be found in the form of the expansion in $\sigma$. Keeping the linear in $\sigma$ terms, after returning to the initial interaction coefficients we obtain (Sedrakian et al., 1999):

[TABLE]

Here, $p_{1}$ is the rotational mode corresponding to the angular velocities misalignment. According to expression (46) the angular velocity difference vector , $\boldsymbol{\nu}_{cg}=\boldsymbol{\Omega}_{c}-\boldsymbol{\Omega}_{g}$, rotates with angular velocity $-\boldsymbol{\Omega}_{0}$ if it is observed from the co-rotating frame of reference. Therefore, in the inertial frame of reference, the direction of vector $\boldsymbol{\nu}_{cg}$ is constant. The absolute value of vector $\boldsymbol{\nu}_{cg}$ decreases with time. The characteristic time-scale of the mode decay is

[TABLE]

Note, that, in contrast to “parallel" mode $p_{\parallel}$, this one contains $\tilde{\epsilon}_{c}$ instead of $\epsilon_{c}$.

The second mode represents the free precession motion. The imaginary part of mode (47) is similar to mode (24). But now it has a real part and, hence, it decays with time as well as the rotational mode. The characteristic time-scale of the precession damping is

[TABLE]

The star will complete of the order of $\Omega_{0}/\alpha_{gc}\gg 1$ precession cycles before the angular velocities align with the symmetry axis but the precession, being, as in the rigid-body case, governed by the pinned superfluid, remains fast.

In the opposite limiting case of strong interaction ($\sigma\gg 1$), the approximate roots obtained as the expansion in $\sigma^{-1}$ are (Sedrakian et al., 1999)

[TABLE]

Again, we have one rotational mode, $p_{1}$, and one precession mode, $p_{2}$. The real part of mode (50) is much larger than its imaginary part. Therefore, being excited, this mode dacays before the first rotational cycle completes. Hence, the rotational mode is not oscillatory in the strong interaction limit.

The real part of (51) relates to the imaginary part approximately as $(\alpha_{ij}/\Omega_{0})^{-1}\ll 1$. Therefore, the precession mode is slowly damped as well as the weak-interaction precession mode. The imaginary part of mode (51) relates to the imaginary part of mode (47) as $I_{c}/(I_{c}+I_{g})$. Hence, the period of precession is longer in the strong-interaction limit but the precession is still governed by the pinned superfluid. Passing $(\alpha_{ij}/\Omega_{0})$ to infinity we reproduce the rigid-body precession mode (24).

In sum, the introduction of the internal dissipation in the way described above leads to the damping of the precession but the precession frequency remains high. However, the two-component approximation gives us another possible configuration.

4 Two-component system with pinning in the core

Let us shift the pinned superfluid into the inner g-component. With all the other assumptions being kept, the angular momenta of the components would have the following form:

[TABLE]

Substituting expressions (52) and (53) into equations (26) and (27), and changing to the c-component frame of reference, we obtain the equations for the angular velocities:

[TABLE]

where vector

[TABLE]

has been reintroduced.

Since the angular momentum of the pinned superfluid is now fixed in the g-component, for an observer in the c-component we have the following equation for vector $\boldsymbol{\omega}_{\mathrm{sf}}$:

[TABLE]

The rotation of the star is stationary if $\boldsymbol{\Omega}_{c}\parallel\boldsymbol{e}_{c}$, $\boldsymbol{\Omega}_{g}\parallel\boldsymbol{\omega}_{\mathrm{sf}}$ and $\boldsymbol{\Omega}_{c}=\boldsymbol{\Omega}_{g}$. Considering the linear perturbations to this equilibrium state, let us first multiply equation (57) by vector $\boldsymbol{e}_{c}$ scalarly. In terms of “parallel” and “perpendicular” components the result has the form

[TABLE]

where $\dagger$ denotes the complex conjugation. Since, in the zeroth approximation, vector $\boldsymbol{\omega}_{\mathrm{sf}}$ is parallel to $\boldsymbol{e}_{c}$, quantity ${\omega}_{\mathrm{sf}}^{\perp}$ can be treated as a small variable together with four variables $\mu_{c}^{\parallel}$, $\mu_{c}^{\perp}$, $\mu_{g}^{\parallel}$ and $\mu_{g}^{\perp}$, which were introduced earlier. Therefore, the right-hand side of equation (58) is quadratically small. Thus, operating in the linear approximation we can put

[TABLE]

Keeping this in mind, we can formulate a system of equations for the linear perturbations:

[TABLE]

Here, we have introduced the g-component effective oblateness parameter as $\tilde{\epsilon}_{g}=\omega_{\mathrm{sf}}^{\parallel}/\Omega_{0}$. Note that due to equality (59) the system of equations falls into two independent subsystems.

Since equations (60) and (61) coincide with equations (36) and (37), the “parallel" mode remains the same. Let us consider the “perpendicular” modes. Substituting $\mu_{i}^{\perp},\omega_{g}^{\perp}\propto\exp(pt)$ into equations (62)–(64) we obtain the third-order characteristic equation

[TABLE]

where we have again used the interaction parameter (43) and introduced the whole star oblateness parameter as

[TABLE]

Assuming weak interaction between the stellar components ($\sigma\ll 1$) we obtain

[TABLE]

The first mode is the rotational mode, similar to mode (46). The second and the the third ones are the precession modes, each of which represents the precession motion of a single component. Mode $p_{3}$ corresponds to the g-component fast precession caused by the superfluid pinned to it. The c-component precession mode $p_{2}$, in contrast to mode (47), contains the real (not “effective”) oblateness parameter. Hence, the c-component precession period is not constrained by the pinned superfluid and therefore it can be sufficiently long if $\epsilon_{c}$ is small enough. Strictly speaking, oblateness parameter $\epsilon_{c}$ can be both positive (oblate star) and negative (prolate star). In the second case the mode is unstable. The reason for this is as follows. Introducing internal dissipation we allow the star to convert mechanical energy into thermal energy. However, the dissipation cannot substantially affect the stellar angular momentum. Thus, the star tends to the state with minimum possible mechanical energy for a given angular momentum. This state is rotation about the major principal axis (Landau & Lifshitz, 1982). If $\epsilon_{c}>0$, this axis is $\boldsymbol{e}_{c}$. Otherwise, it is an axis perpendicular to $\boldsymbol{e}_{c}$. In the second case, one just needs to redefine the basic equilibrium state. For definiteness, further we will assume that $\epsilon_{c}>0$.

In the case of strong interaction, the approximate roots of equation (65) are the following:

[TABLE]

The first two of these are similar to the corresponding roots obtained for the c-component superfluid model (cf. (50) and (51)). It is easy to verify that combination $(\tilde{\epsilon}-{\epsilon}_{c})(\tilde{\epsilon}_{g}-\tilde{\epsilon})$ is always positive and hence the real part of mode (71) is always negative. Therefore, mode (71) is damped. The third mode arises when both $\epsilon_{c}$ and $\tilde{\epsilon}_{g}$ are non-zero. However, it is easy to see that this mode is not oscillatory.

Thus, if we assume that the pinning takes place in the internal component, the slow long-lived precession mode can exist in the case of weak interaction. It is also necessary to make sure that this mode is observable. Substituting mode (68) into equations (62) - (64), we can express

[TABLE]

where $O(\sigma)$ are small corrections due to the component interaction. Using expressions (52) and (53) it can also be shown that

[TABLE]

where $\boldsymbol{M}=\boldsymbol{M}_{c}+\boldsymbol{M}_{g}$ is the total angular momentum. We see that, if the $p_{2}$ mode is excited, vectors $\boldsymbol{\Omega}_{g}$ and $\boldsymbol{L}_{\mathrm{sf}}$ remain directed almost along the (fixed in the inertial frame of reference) stellar angular momentum $\boldsymbol{M}$ while vectors $\boldsymbol{e}_{c}$ and $\boldsymbol{\Omega}_{c}$ are declined from it. In the case of an infinitely small parameter $\sigma$, this mode corresponds to the independent precession of the external c-component.

5 Three-component system

Having considered the basic idea, we turn to the full model. Let us assume that the neutron star consists of three dynamically distinguished components that we will denote as c-, g-, and r-component.

The c-component is the outer component. The c-component as before is assumed to have axisymmetric shape. Its angular momentum can be represented as

[TABLE]

The g-component is an inner component. It consists of a “normal" fraction and a pinned superfluid fraction. The “normal" fraction, for more generality, is assumed to have axisymmetric shape. It is characterized by moment of inertia $I_{g}$ and oblateness parameter $\epsilon_{g}$. The “normal" fraction rotates with angular velocity $\boldsymbol{\Omega}_{g}$. The pinned superfluid contains angular momentum $\boldsymbol{L}_{\mathrm{sf}}$. Thus, the total angular momentum of the g-component can be represented as

[TABLE]

where $\boldsymbol{e}_{g}$ is unit vector directed along the g-component symmetry axis. We will assume that vector $\boldsymbol{L}_{\mathrm{sf}}$ is directed along the g-component symmetry axis, i.e.

[TABLE]

where $\boldsymbol{\omega}_{\mathrm{sf}}$ is defined by expression (56).

The r-component is the other inner component. It rotates with angular velocity $\boldsymbol{\Omega}_{r}$. Its possible asymmetries would not lead to substantially new effects but would complicate the expressions. Hence, we will assume this component to be spherically symmetric. Therefore, the angular momentum of the r-component is assumed to have the form

[TABLE]

The angular momentum conservation law can be represented in the form of three equations:

[TABLE]

where $i,\ j=c,\ g,\ r$ and $\boldsymbol{N}_{ij}$ is the torque acting on the $i$th component from the $j$th component. In the general case, all the components are supposed to interact with each other.

The further analysis is similar to that performed for the two-component system. We will describe it briefly mentioning the main steps and focusing only on new features. Substituting expressions (76)-(79) into equations (80) and changing to the frame of reference co-rotating with the c-component we obtain the following equations for the angular velocities:

[TABLE]

From the angular momentum conservation it is follows that

[TABLE]

In order to close this system, equation (57) should be added.

It is assumed that if

[TABLE]

and $\boldsymbol{\Omega}_{c}\parallel\boldsymbol{\omega}_{\mathrm{sf}}\parallel\boldsymbol{e}_{c}$, the system is in stable equilibrium. The reasoning leading to approximate equality (59) remains valid. Hence, if the star rotates near the fixed point, there are seven small variables, namely $\mu_{c}^{\parallel}$, $\mu_{c}^{\perp}$, $\mu_{g}^{\parallel}$, $\mu_{g}^{\perp}$, $\mu_{r}^{\parallel}$, $\mu_{r}^{\perp}$ and $\omega_{\mathrm{sf}}^{\perp}$.

We assume that there are no preferred directions for the components interaction except the angular velocities. However, since the exact interaction mechanisms are not specified at this point, we consider a more general form of vectors $\boldsymbol{R}_{ij}$. Namely,

[TABLE]

where $\xi_{ij}=\beta_{ij}+i\gamma_{ij}$ and $\alpha_{ij}$, $\beta_{ij}$, $\gamma_{ij}$ are the phenomenological interaction constants. The particular case of simple linear friction corresponds to $\beta_{ij}=\alpha_{ij},\gamma_{ij}=0$. If we consider, for instance, the mutual friction interaction between the $i$- and $j$-component, one of which is superfluid, the coefficients can be represented in the following form (Sedrakian et al., 1999; Barsukov et al., 2013):

[TABLE]

where $x$ is the coupling parameter. From relations (84) it follows that

[TABLE]

The linearized system of equations has the form:

[TABLE]

where the g-component effective oblateness parameter is equal to

[TABLE]

As before, the system falls into “parallel” and “perpendicular” parts.

Let us first consider “parallel” modes. Substituting $\mu_{i}^{\parallel}\propto\exp(pt)$ into equations (90)–(92) we obtain the following second-order characteristic equation:

[TABLE]

Here, we have introduced three interaction parameters

[TABLE]

it is easy to see that $\sigma_{ij}=\sigma_{ji}$. The exact roots of equation (98) can be given. However, instead of that we focus on the one particular case of

[TABLE]

which appears to be the most suitable for the glitch model considered in Section 7. With relations (89) the roots can be represented in the form

[TABLE]

In this particular case, $p_{+}\gg p_{-}$.

Let us turn to the “perpendicular" part. We do not give here the fourth-order characteristic equation because of its awkwardness. By analogy with the “parallel" case we can introduce three interaction parameters

[TABLE]

In the weak-interaction limit ($|Z_{ij}|\ll 1$) the approximate expressions for the modes can be represented in the form:

[TABLE]

With the introduction of the third component the rotational mode becomes degenerate. It splits into two modes when the dissipation corrections $\delta p_{1}$ and $\delta p_{2}$ are taken into account. They can be found as the roots of the following equation:

[TABLE]

To make some estimations by analogy with "parallel" modes we assume that $|Z_{cg}|\gg|Z_{cr}|,|Z_{gr}|$. This allows us to treat relations $Z_{cr}/Z_{cg}$ and $Z_{gr}/Z_{cg}$ as small parameters and represent the roots of equation (108) in the form of expansion in them. The resulting expressions formally coincide with expressions (101) and (102) with corresponding replacements of $\alpha_{ij}$ by $\xi_{ij}$ and $\epsilon_{g}$ by $\tilde{\epsilon}_{g}$. Since our goal is estimating the damping time-scales, for the sake of brevity, we neglect all the oblateness parameters, assuming their smallness. To the first order in the small interaction parameters we have

[TABLE]

Another possible case is $\tilde{\epsilon}_{g}\gg$1, for which we have

[TABLE]

Thus, assuming that the cg-interaction is much stronger than the other two, we obtain the hierarchy of damping time-scales. Namely, the real part of $p_{+}$ and $p_{1}$ is determined by the strength of the strongest interaction (i.e. cg-intraction) while $p_{-}$ and $p_{2}$ modes are damped rather due to the cr- and gr-interactions.

Mode $p_{3}$ corresponds to the c-component slow precession. It is almost similar to mode (68) obtained in the framework of the two-component model. However, let us note that the precession period, in the general case, would be equal to

[TABLE]

Mode $p_{4}$ corresponds to the fast precession of the g-component. We do not take into account the possible non-sphericity of the r-component. However, one can see that in the weak-interaction limit each component precesses almost independently. Hence, the precession mode for the r-component can easily be obtained from $p_{3}$ or $p_{4}$ by the corresponding interchange of the indices.

We restrict ourselves to considering only the weak-interaction limit between all pairs of components. We have seen in Section 4 that the strong-interaction limit between the c- and g-component does not allow the observed c-component to precess with long period.

Hence, the gc-interaction should be weak. On the other hand, we will see in sec. 7 that the cg-interaction should be the strongest among all pair interactions. Thus, the only possible case is the weak interaction between all the components. Recall that by weak interaction we mean that

[TABLE]

6 Quasi-stationary evolution

We have seen that the long-period precession mode can arise if we assume that the region of superfluid pinning can rotate relative to the crust. This result was obtained for a freely rotating neutron star perturbed from stable equilibrium. Such a formulation of the problem, however, could be not quite realistic. As it was mentioned at the beginning of the paper, strictly speaking, the rotation of neutron stars is not free. Even isolated neutron stars rotate under the action of electromagnetic torque $\boldsymbol{K}$ caused by the rotation of the strong magnetic field anchored to the star (Davis & Goldstein, 1970; Beskin et al., 1983; Melatos, 2000). The part of this torque that we have denoted by $\boldsymbol{K}_{2}$ can be easily taken into account by redefinition of the stellar moment of inertia tensor (see Section 2). As for the second part, $\boldsymbol{K}_{3}$, we have not taken it into account because its effects are negligibly small at precession time-scales. Now we are going to formulate the equations that allow us to consider the evolution of pulsar rotation at the pulsar life time-scales.

We will consider the three-component model introduced in the previous section. Equations (82), (83) and (57) remain the same. Equation (81) should be replaced by a slightly modified one that is more convenient for subsequent consideration. First, let us re-establish the third-order external torque term $\boldsymbol{K}_{3}$, which we put on the right-hand side of equation (81). After simple rearrangement the equation takes the following form:

[TABLE]

where we have introduced vector

[TABLE]

Multiplying equation (115) by $\boldsymbol{e}_{c}$ and substituting the obtained expression for $(\boldsymbol{e}_{c}\cdot\boldsymbol{\Omega}_{c})$ back into equation (115) we get

[TABLE]

where by $\boldsymbol{\mathrm{RHS}}$ we mean the right-hand side of equation (115). Recall that the c-component oblateness parameter is supposed to be of the order of $10^{-8}$ or smaller. Therefore, the second term on the right-hand side of equation (117) is negligibly small. Thus, the c-component equation becomes

[TABLE]

The c-component angular velocity vector time derivative can be represented as a sum of two terms:

[TABLE]

where $\boldsymbol{e}_{\Omega}=\boldsymbol{\Omega}_{c}/\Omega_{c}$. Here, the first term represents the change of the angular velocity absolute value while the second term arises due to change of its orientation. Again it will be convenient to introduce the “parallel” and “perpendicular” parts of the vectors. However, in contrast to previous sections, we define the parts relative to the direction of vector $\boldsymbol{\Omega}_{c}$ instead of $\boldsymbol{e}_{c}$. Let $\boldsymbol{e}_{x},\ \boldsymbol{e}_{y},\ \boldsymbol{e}_{z}$ be an orthonormal basis oriented such that $\boldsymbol{e}_{z}\parallel\boldsymbol{e}_{\Omega}$ and $\boldsymbol{e}_{x}\parallel\dot{\boldsymbol{e}}_{\Omega}$. Having fixed the basis, we can introduce the following notations:

[TABLE]

where $\boldsymbol{V}$ is an arbitrary vector. Taking into account expression (119) and (120), we can represent equation (118) as two equations:

[TABLE]

The next step is linearization. Let us introduce departure vectors

[TABLE]

We will assume that all the components rotate with almost the same angular velocity and the g-component rotates almost about its symmetry axis $\boldsymbol{e}_{g}$. However, the angle between $\boldsymbol{\Omega}_{c}$ and $\boldsymbol{e}_{c}$ is now not necessarily small. Multiplying equation (57) by vector $\boldsymbol{e}_{\Omega}$ we obtain

[TABLE]

Since, in the zeroth approximation, vector $\boldsymbol{\omega}_{\mathrm{sf}}$ is parallel to $\boldsymbol{e}_{\Omega}$, quantity ${\omega}_{\mathrm{sf}}^{(\perp)}$ can be treated as a small variable together with the components of the angular velocity departure vectors. Hence, the first term on the right-hand side of equation (124) is straightforwardly quadratically small. According to equation (122) time derivative $\dot{{e}}_{\Omega}$ is linearly small. Therefore, the second term on the right-hand side of equation (124) is quadratically small as well. Thus, the parallel part of vector $\boldsymbol{\omega}_{\mathrm{sf}}$ with linear accuracy can be treated as a constant.

The interaction between the components is described by vectors $\boldsymbol{R}_{ij}$. It is assumed that in the small departure approximation (${\nu}_{ij}/\Omega_{0}\ll 1$) they can be represented in the following form:

[TABLE]

where $\xi_{ij}=\beta_{ij}+i\gamma_{ij}$ and $\alpha_{ij}$, $\beta_{ij}$, $\gamma_{ij}$ are the interaction constants. If $\boldsymbol{e}_{\Omega}\approx\boldsymbol{e}_{c}$, this representation, up to linear terms, reproduces expressions (86) and (87). Hence, the interaction constants $\alpha_{ij}$, $\beta_{ij}$, $\gamma_{ij}$ coincide with the coresponding interaction constants introduced in Section 5.

Substituting definition (123) into equations (82), (83) and (57) and neglecting quadratically small terms, we obtain

[TABLE]

where

[TABLE]

The equation for ${\nu}_{gr}^{(\parallel)}={\nu}_{gc}^{(\parallel)}-{\nu}_{rc}^{(\parallel)}$ and ${\nu}_{gr}^{(\perp)}={\nu}_{gc}^{(\perp)}-{\nu}_{rc}^{(\perp)}$ can be obtained from equations (127) – (130). Equations (127)-(131) together with equations (121) and (122) form a closed system describing the rotation evolution of our three-component neutron star.

We have linearized the system of equations with respect to $\nu_{gc}^{(\parallel)}$, $\nu_{rc}^{(\parallel)}$, ${\nu}_{gc}^{(\perp)}$, ${\nu}_{rc}^{(\perp)}$ and ${\omega}_{\mathrm{sf}}^{(\perp)}$. However, equations (121), (122) and (127) - (131) contain different combinations with vector $\boldsymbol{\Omega}_{c}$ which is not small. In its present form this system of equations is quite difficult to solve. Fortunately, the analysis performed in Section 5 gives us a way to simplify this system.

Let us highlight the following:

The observed quasi-periodic processes that are interpreted as a manifestation of the long-period precession mode have the time-scale $T_{p}\sim 10^{2}$ d or longer and, hence, $\epsilon_{c}\sim 10^{-8}$ or smaller. 2. 2.

If at least 1 percent of the total moment of inertia is contained in the g-component pinned superfluid, according to expression (97) the g-component effective oblateness parameter $\tilde{\epsilon}_{g}$ cannot be smaller that $\sim 10^{-2}$. The estimations that will be done in Section 7 provide $\tilde{\epsilon}_{g}\sim 10$. 3. 3.

The interactions between the components are assumed to be weak but not too weak. It is required that $\alpha_{ij}/\Omega_{c},|\xi_{ij}|/\Omega_{c}\ll 1$ for all pairs of the components. Let us also constrain the coefficient values from the other side assuming that

[TABLE]

From points 1 – 3 it is follows that there is a gap between the long-period precession time-scale ($\sim T_{p}$) and the time-scales of all the other modes except the long-period precession one. In other words, if we introduce the internal relaxation time-scale $\tau_{r}$ as a smallest real part of all the other modes except the long-period precession mode, condition $\tau_{r}\ll T_{p}$ is assumed to be satisfied. Such a hierarchy of time-scales allows us to consider the behaviour of the system on the different time-scales separately.

Let us first consider the evolution of the system at the internal relaxation time-scale. At this time-scale, the absolute value of angular velosity $\boldsymbol{\Omega}_{c}$ and its orientation relative to vector $\boldsymbol{e}_{c}$ can change only slightly. Therefore, we can treat vector $\boldsymbol{\Omega}_{c}$ as well as vector $\boldsymbol{S}$ as constants. In this case, equations (127) – (131) become a closed system of linear equations. The system is inhomogeneous. Therefore, the general solution consists of a particular solution to system (127) – (131) and the general solution to this system with zeroth $\boldsymbol{S}$. The first can be found by putting all the time derivatives in equations (127) – (131) equal to zero. We will call this solution “quasi-stationary" where by “quasi” we mean that it is an approximate solution, valid as long as vector $\boldsymbol{\Omega}_{c}$ can be treated as a constant. The general solution of homogeneous equations is a sum of modes proportional to $\exp(pt)$. It is not difficult to see that they are the same modes as were obtained in Section 5 (except the long-period precession mode). Indeed, making the following formal change of variables:

[TABLE]

and renaming $\Omega_{0}$ to $\Omega_{c}$, we can obtain system of equation (127) – (131) with zeroth $\boldsymbol{S}$ from system of equations (91) – (96) with zeroth $\epsilon_{c}$. Since all the real parts of modes (101), (102) and (104) – (107) are negative, small perturbations decay to the quasi-stationary solution.

The quasi-stationary values of $\nu_{gc}^{(\parallel)}$, $\nu_{rc}^{(\parallel)}$, ${\nu}_{gc}^{(\perp)}$, ${\nu}_{rc}^{(\perp)}$ and ${\omega}_{\mathrm{sf}}^{(\perp)}$ can be obtained from equations (127) – (131) with simple algebra. The expressions, however, are quite cumbersome. We do not give exact expressions, restricting ourselves to estimations. Let us first look at the “parallel" part of the system of equations. Since the time derivatives are negligibly small, the non-zeroth value of $S^{(\parallel)}$ can be counterbalanced only by interaction torques $R_{ij}^{(\parallel)}$, i.e. $R_{ij}^{(\parallel)}\sim S^{(||)}$ and $\nu_{ij}^{(\parallel)}\sim S^{(||)}/\mathrm{min}(\alpha_{jk})$. The precession term does not contribute to $S^{(\parallel)}$. Therefore, $S^{(\parallel)}\sim K_{3}/I_{c}\sim(I/I_{c})(\Omega_{c}/\tau_{c})$ where $\tau_{c}$ is the characteristic age of the pulsar. For the “perpendicular" part we have ${\nu}_{ij}^{(\perp)},\ {\omega}_{\mathrm{sf}}^{(\perp)}\sim{S}^{(\perp)}/\Omega_{c}$. The leading term in ${S}^{(\perp)}$ is the precession term. Hence, ${\nu}_{ij}^{(\perp)},\ {\omega}_{\mathrm{sf}}^{(\perp)}\sim\epsilon_{c}\Omega_{c}$. Thus, we see that the quasi-stationary values of these quantities are small as it was assumed.

Let us turn to the long-period precession time-scale. If the small perturbations are not excited during stellar evolution, they should all be damped on this time-scale. However, we obviously can no longer ignore the evolution of vector $\boldsymbol{\Omega}_{c}$. It produces corrections to the quasi-stationary expressions arising to compensate the time derivatives in equations (127) – (131). The time derivatives, in turn, arise due to changing the components of vector $\boldsymbol{\Omega}_{c}$ contained in the quasi-stationary expressions. Hence, $\dot{\nu}_{ij}\sim\nu_{ij}(\dot{\Omega}_{c}/\Omega_{c})$. According to equations (121) and (122) we have $\dot{\Omega}_{c}\sim\alpha_{ij}\nu_{ij}$ and ${\dot{e}}_{\Omega}\sim\nu_{ij}$. Hence, $\dot{\nu}_{ij}\sim\nu_{ij}^{2}$, i.e. the corrections are quadratically small. Thus, we can continue to use the quasi-stationary expressions for small departures at the long-period time-scales if the three following conditions are satisfied. First, the departures from the solid-body rotation should be small. Secondly, the internal relaxation time-scale should be much smaller than the slow precession period $T_{p}$ (ensured by conditions (133)). Thirdly, the time-dependent perturbations are not excited during the star’s life-time. The last condition is obviously violated if the neutron star is glitching. This case is the subject of the next section.

The quasi-stationary approximation allows us to exclude the rotation of the internal components from consideration. Let us first look at the “parallel" part. The quasi-stationary expressions for $\nu_{gc}^{(\parallel)}$ and $\nu_{rc}^{(\parallel)}$, being expressed from equations (127) and (128), can be substituted into the right-hand side of equation (121). However, a more illustrative way to obtain the same result is the following. Taking a look at equations (127) and (128), we notice that the three last terms on their right-hand sides (cf. equation (122)) are exactly equal to $\dot{\Omega}_{c}$. After transposing $\dot{\Omega}_{c}$ to the left-hand sides of the equations one can compose the combination $I_{c}$(121)$+(1+\epsilon_{g})I_{g}$(127)$+I_{r}$(128). Neglecting quadratically small time derivatives $\dot{\nu}_{ij}^{(\parallel)}$ and taking into account relations (84), we finally obtain

[TABLE]

where we have denoted $\tilde{I}=I_{c}+(1+\epsilon_{g})I_{g}+I_{r}$. According to this equation the neutron star is braked by the external torque as if it is a rigid body with the moment of inertia equal to $\tilde{I}$. It does not contain the moment of inertia of the pinned superfluid because the superfluid does not slow down as long as it remains pinned. The rigid-body braking is a general feature of the quasi-stationary approximation (Barsukov et al., 2014, section 2.3).

Let us turn to the “perpendicular" part. According to equation (122), the right-hand side of equation (131) is equal to $-\dot{{e}}_{\Omega}$. Therefore, from equation (131) in the quasi-stationary approximation ($\dot{{\omega}}_{\mathrm{sf}}^{(\perp)}\approx 0$) we have

[TABLE]

Substituting equation (138) into equations (129) and (130) (with zeroth time derivatives), we obtain the system of linear equations for $\nu_{rc}^{(\perp)}$ and $\omega_{\mathrm{sf}}^{(\perp)}$, allowing us to express these quantities as functions of $\dot{{e}}_{\Omega}$. The solution is the following:

[TABLE]

We see that all three quantities (138) – (140) contain the same term $i\dot{{e}}_{\Omega}$. The last two contain some corrections to it that, however, are small in the weak-interaction limit. Substituting equations (138) and (139) into equation (122) and making some rearrangement, we finally obtain the equation for $\dot{e}_{\Omega}$, which after returning to the vector form can be represented as

[TABLE]

where coefficients $B$ and $\Gamma$ are determined by equality

[TABLE]

Here we have used notation

[TABLE]

However, all steps of the derivation of equation (141) remain valid if we assume a more general form for the c-component angular momentum

[TABLE]

instead of form (76), where $\hat{\boldsymbol{\epsilon}}$ is an arbitrary symmetric tensor describing small deviations of the c-component from spherical symmetry.

Since the weak-interaction limit is assumed ($|\xi_{ij}|/\Omega_{c}\ll 1$), the coefficients in equation (141) approximately are equal to

[TABLE]

Replacing expression (142) by expressions (145), we ignore the difference between the rotation of the two internal components (the second term in equation (139) is neglected).

If the internal components do not interact with the c-component at all, both coefficients $B$ and $\Gamma$ are equal to zero. In this case, equation (141) takes the form

[TABLE]

With $\boldsymbol{K}_{3}=0$ and equality (143) it reproduces the rigid-body free precession equation (8). The influence of $\boldsymbol{K}_{3}$ at the precession evolution is beyond the scope of the present paper. This question has already been well studied (Goldreich, 1970; Melatos, 2000; Arzamasskiy et al., 2015; Goglichidze et al., 2015). The electromagnetic torque can force the precession amplitude to both increase or decrease at the time-scales of pulsar characteristic age $\tau_{c}$. It is difficult to formulate any more specific statement because the result depends strongly on the mutual orientation of the stellar magnetic axis and its principal axes as well as on the exact form of the external electromagnetic torque $\boldsymbol{K}_{3}$.

Next, let us switch off the external torque but take into account the interaction with the internal components. If $\boldsymbol{K}_{3}=0$, according to equation (137) we have $\dot{\boldsymbol{\Omega}}_{c}=\Omega_{c}\dot{\boldsymbol{e}}_{\Omega}$. Since the interaction is assumed to be weak, we keep only the linear in $B$ and $\Gamma$ terms. The equation in the axisymmetric case would be the following:

[TABLE]

The interaction with the internal components modifies the rigid-body precession equation (8) in two ways. Coefficient $\Gamma$ produced by non-dissipative interaction renormalizes the precession period (cf. with expression (10)):

[TABLE]

It is easy to see that in the small-angle approximation this coincides with value (113) obtained in the framework of the linear analysis. Considering the influence of coefficient $B$, let us multiply equation (147) by vector $\boldsymbol{e}_{c}$:

[TABLE]

Here, $(\boldsymbol{e}_{c}\cdot\boldsymbol{\Omega}_{c})>0$ by choosing the direction of $\boldsymbol{e}_{c}$; coefficient $B$ calculated with formula (145) is positive because interaction (125) should reduce the difference between the rotation the components. Therefore, angular velocity $\boldsymbol{\Omega}_{c}$ tends to align with the symmetry axis $\boldsymbol{e}_{c}$. The characteristic time-scale of the alignment is

[TABLE]

This is the same time-scale that we obtained for $p_{3}$ mode in Section 5.

The quasi-stationary evolution formalism confirms the results obtained in the framework of the linear mode analysis for the more general case of an arbitrary angle between the angular velocity and the c-component symmetry axis. However, the main advantage of this formalism is that it allows us to obtain equations (137) and (141), which can be used for studying the stellar rotation evolution at the time-scales comparable to the pulsar life time.

7 Glitch-like event

We have formulated a model of a rotating neutron star with pinned superfluid that allows the neutron star to precess with long period. The model should also be examined for the ability to demonstrate glitch-like behaviour.

We assume that the glitch relaxation (the longest glitch stage) is governed by internal relaxation processes. As it was discussed in the previous section, at internal relaxation time-scale $\tau_{r}$, vectors $\boldsymbol{\Omega}_{c}$ and $\boldsymbol{S}$ can be treated as constants in equations (127) – (131). In this case, the full solution can be represented as a sum of the quasi-stationary solution (considered in the previous section) and a time-dependent solution to the homogeneous system of equations (with zeroth $S^{(\parallel)}$ and $S^{(\perp)}$). The last, in turn, is a sum of the linear modes studied in Section 5 (except the $p_{3}$ mode, which has already been taken into account in the quasi-stationary solution; see Section 6). Generally speaking, the glitch can excite both “parallel" and “perpendicular" modes. However, in the present paper, we will consider only the “parallel" perturbations, which allow us to reproduce the main glitch manifestation – the pulsar rotation frequency jump. Thus, the system of equations is reduced to

[TABLE]

The exact glitch-triggering mechanism is not known at present (Haskell & Melatos, 2015). The mechanism is likely based on some non-linear process that cannot be described by simple linear equations. We will model it in the following way. It is assumed that at $t=0$ the superfluid pinned in the g-component instantly releases a small part of the stored angular momentum $\Delta\boldsymbol{L}_{\mathrm{sf}}$. Since we are interested in studying the “parallel" modes, released angular momentum $\Delta\boldsymbol{L}_{\mathrm{sf}}$ is supposed to be parallel to $\boldsymbol{\Omega_{c}}$. The whole released angular momentum is assumed to be instantly injected into the “normal" fraction of the g-component such that the total g-component angular momentum is conserved. Thus, at $t=0$ the value of $\omega_{\mathrm{sf}}^{(\parallel)}$ decreases by $\Delta L_{g}/I_{g}$ while angular velocity $\Omega_{g}$ increases by $\Delta L_{g}/I_{g}$. The value of $\omega_{\mathrm{sf}}^{(\parallel)}$ is not contained in equations (151) – (153). Thus, right after the triggering we have: $\Omega_{c}=\Omega_{0}$, $\nu_{rc}^{(\parallel)}=0$, $\nu_{gc}^{(\parallel)}=\Delta L_{g}/I_{g}$. We can consider these values as initial conditions and study how the system evolves.

To find the c-component rotation response to the glitch in the g-component we can substitute the solution to equations (152) – (153) into equation (151) and then integrate it. The result can be represented in the following form:

[TABLE]

where $\Delta\Omega={\Delta\Omega_{\infty}}/{(1-Q)}$, $\Delta\Omega_{\infty}={\Delta L_{g}}/{\tilde{I}}$,

[TABLE]

and coefficients $p_{+}$ and $p_{-}$ are the two “parallel" modes. If one wants to relate the solution to observed pulsar glitches, function (154) should have the form shown in Fig. 1.

Values $1/p_{+}$ and $1/p_{-}$ should be interpreted as glitch rise and glitch relaxation time-scales respectively. Hence, it is required that $p_{+}\gg p_{-}$. This condition is satisfied if one of the pair interactions is much stronger than the other two (see Section 5). Let us first assume that assume that the strongest one is the cg-interaction, i.e. inequality (100) is satisfied. In this case, coefficients $p_{+}$ and $p_{-}$ can be calculated with expressions (101) and (102). If we want to reproduce the angular velocity behaviour sketched in Fig. 1, we should also ensure that (Lyne et al., 2000)

[TABLE]

Substitutting expressions (101) and (102) into expression (155) we obtain

[TABLE]

If we choose another interaction to be the strongest one assuming, for instance, that $\sigma_{cr}\gg\sigma_{cg},\sigma_{gr}$, the estimations for $p_{+}$ and $p_{-}$ can be obtained from expressions (101) and (102) by the corresponding interchange of indices. However, condition (156), in this case, requires fine tuning of the interaction parameters, which is unlikely to be able to be maintained for long time since the parameters evolve with the neutron star’s internal temperature. Thus, the case of inequality (100) is most plausible.

From the physical point of view the picture is as follows. The glitch suddenly increases the g-component “normal" fraction angular velocity making angular velocity lag $\nu_{gc}^{(\parallel)}$ exceed its quasi-stationary value. Right after that the c-component spins up due to the cg-interaction at the (observationally unresolved) $p_{+}$-time-scale. The spinning up is followed by spinning down due to the reaction of the r-component. Thus, in the framework of the proposed model, the glitch relaxation is provided by the presence of the r-component. Indeed, putting $I_{r}=0$ in expression (157), we obtain $Q=0$ and, hence, $\Delta\Omega_{\infty}=\Delta\Omega$.

Observations could give us some constraints for $p_{+}$ and $p_{-}$ and, hence, for the interaction coefficient parameters. From expression (101) it is follows that

[TABLE]

where $\tau_{\mathrm{glitch}}$ is the characteristic glitch rise time. Supposing that $\beta_{ij}\sim\alpha_{ij}$ and taking into account condition (100), we can relate $\sigma_{cg}$ to $B$:

[TABLE]

The relaxation time-scale $\tau_{\mathrm{g.relax.}}$ allows us to estimate the other interaction coefficients. From expression (102), we have

[TABLE]

8 Discussion

We have formulated a simple model of a rotating neutron star supposing the existence of three abstract dynamically distinguished components.

Let us now speculate on the possible physical counterparts of these components.

The outer c-component can be represented by the neutron star’s crust and part of its core strongly coupled with the crust. Therefore, the moment of inertia of the c-component can be estimated as $I_{c}\sim(10^{-2}-10^{-1})I$.

The role of the g-component can be performed by tangles of closed flux tubes that could be formed from chaotic small-scale magnetic field after protons became superconductive. Alternatively it could be a torus composed of closed flux tubes (Gügercinoğlu & Alpar, 2014). If, for instance, the characteristic cross-section $S_{\mathrm{tor}}$ of the region occupied by the toroidal field is of the order of 1 km2, then $I_{g}\sim\rho_{p}S_{\mathrm{tor}}r_{\mathrm{ns}}^{3}\sim 10^{-3}I$, where $\rho_{p}$ is the proton mass density. Some of the superfluid neutron vortices located in the core can be pinned to the closed flux tubes. On the one hand, this interaction prevents the tangles or torus collapsing. On the other hand, when the critical rotational lag is reached, the vortices unpin, triggering the glitch. These pinned vortices carry angular momentum $\boldsymbol{L}_{\mathrm{sf}}$. If we assume that 1 percent of the total stellar moment of inertia is contained in the pinned superfluid, then ${L}_{\mathrm{sf}}/I_{g}\Omega_{0}\approx 10^{-2}(I/I_{g})$. Assuming that a thin ring is a good approximation for the g-component mass distribution, the real oblateness parameter can be estimated as $\epsilon_{g}\approx 2$. Hence, the second term in expression (97) is dominant and $\tilde{\epsilon}_{g}\sim 10^{-2}(I/I_{g})\sim 10$.

Since the flux tubes are closed in the core, the g-component, being magnetically decoupled, can rotate with an angular velocity different from $\boldsymbol{\Omega}_{c}$ (Glampedakis & Lasky, 2015). As it has been argued, the cg-interaction should be quite strong to ensure the rapid c-component spin up but not too strong to allow the g-component to rotate relative to the crust. We guess that this interaction can be implemented by the viscosity of crust-core interface region with the possible inclusion of the Ekman pumping mechanism. Unfortunately, there are not many observational data for glitch rise times at present. The Crab pulsar is the only pulsar for which partially resolved spin-up data have been obtained (Shaw et al., 2018). The delayed frequency increase lasted for 0.5–1.7 days. However, this phenomenon was observed only in three large Crab glitches. Therefore, it could be more correct to apply the upper limit for the fast unresolved part of the glitch spin-up for estimating the cg-interaction coefficients. For the 2017 glitch Shaw et al. (2018) placed the upper limit $\tau_{\mathrm{glitch}}<6$ h. Hence, for the Crab pulsar with $P=3.3\times 10^{-2}$ s (Manchester et al., 2005)222http://www.atnf.csiro.au/research/pulsar/psrcat/, using estimation (158), we obtain $\sigma_{cg}>2.4\times 10^{-7}$. For the Vela pulsar much stiffer restrictions were obtained. For 2004 glitch, Dodson et al. (2007) obtained that $\tau_{\mathrm{glitch}}<30$s and, hence, $\sigma_{cg}>4.7\times 10^{-4}$, where $P=8.9\times 10^{-2}$ s for the Vela pulsar22footnotemark: 2. The value of $\sigma_{cg}$ can vary substantially from one pulsar to another due to different internal temperatures and magnetic flux tube organization. Moreover, we do not state that the presented model describes all glitching pulsars. In some of them the g- and the c-component could be tightly coupled. This part, however, cannot precess with long periods, as it was shown earlier.

From the observed glitch relaxation time-scales (Lyne et al., 2000), using expression (160), we obtain the following estimation for the corresponding interaction parameters:

[TABLE]

The weakly interacting r-component is likely composed of superfluid matter. We suppose this component consists of some of the neutrons superfluid that is not pinned and possibly some of the normal matter coexisting with it and weakly coupled with c-component. Hence, in order of magnitude, we can estimate that $I_{r}\sim I$. Strictly speaking, it would be more consistent to treat the superfluid and normal fractions of the r-component as two separate components because of their weak interaction. It would complicate the calculations and would probably lead to a more complex glitch relaxation but qualitatively the model remains the same.

It is generally accepted that the neutrons in the inner crust should be superfluid and, hence, the superfluid vortices can pin to the crust lattice. According to theoretical calculations, about 1 percent of the total moment of inertia is contained in the crust superfluid. If all that superfluid is pinned, the crust’s effective oblateness parameter would be of the order of $10^{-2}$ and, hence, long-period precession becomes impossible. However, as shown by Link & Cutler (2002), the Magnus force, acting on the crust vortices could be enough to unpin them. Thus, we assume that the crust superfluid is unpinned in precessing stars. Strictly speaking, in this case, it should be included in the model as an additional component. Again, it can modify coefficients $B$ and $\Gamma$ and complicate the glitch relaxation. The exact effects of mentioned additional components are the subject of future study.

Ashton et al. (2017) found that the modulations in spin-down rate and beam shape of pulsar B1828-11 become faster. If the modulations are induced by precession, the rate of precession period decrease is $\dot{T}_{p}\approx-10^{-2}\mbox{s s}^{-1}$. According to expression (10) this may indicate the gradual increase of the star’s oblateness, which would be quite counter-intuitive. Alternatively, the variation of precession period could be caused by a decrease in angle $\theta$. In the case of a rigidly rotating star, the angle evolution can be caused only by electromagnetic torque $\boldsymbol{K}$; this process would occur on the same time-scale as the pulsar spin-down and, hence, it cannot be the source of such fast variation. However, the pulsar could align due to internal dissipation. Taking the time derivative of precession period (148) and using equation (149), we can obtain the estimation for the corresponding coefficient:

[TABLE]

Since the spin-down rate is not affected by internal dissipation (see equation (137)) in the quasi-stationary approximation, the rotation frequency can be treated as a constant at the considered time-scale. Due to the symmetry of the problem there are two configurations with $\theta=5^{\circ}$ and $\theta=89^{\circ}$, which fit the observational data equally well (Arzamasskiy et al., 2015). For the small angle configuration with estimations (162) and (159) we formally obtain $B\approx 0.2$ and $\sigma_{cg}\sim 1-10$. This case is not self-consistent because the weak-interaction condition is not satisfied. In the case of the large angle configuration we have $B\approx 3\times 10^{-5}$ and $\sigma_{cg}\sim 10^{-4}-10^{-3}$. This, according to estimation (158), leads for the pulsar with $P=0.4$ s to $\tau_{\mathrm{glitch}}\sim 10^{2}-10^{3}$ s. These values seem to be quite plausible. The configuration with large angle $\theta$ assumes a small angle between the symmetry axis and the magnetic dipole axis. Therefore, speculatively we can guess that, some time in the past, a starquake occurred, which led to the rearrangement of the crust almost along the magnetic axis. This occasion started the precession. The phenomenon of decreasing $T_{p}$, however, could have an alternative explanation. For example, the crust superfluid, which is unpinned due to precession, can gradually re-pin, which caused the increase of the crust’s effective oblateness.

In the present paper we proposed a way in which the vortex-pinning-based pulsar glitches can be reconciled with the long-period precession. Basically it is assumed that the superfluid is pinned in the region located in the stellar interior, which has sufficient freedom to rotate relative to the crust. This allows the pinned superfluid’s angular momentum to stay aligned almost along the stellar rotation axis and, hence, to not affect the precession period. At the same time, this region should be sufficiently effectively coupled with the crust to ensure the sufficiently rapid spin-up of the crust after a glitch there has happened. Since the precession and the glitches are governed by different regions of the star, it is natural that the glitch has almost no effect on the precession characteristics (Jones et al., 2017). The exact glitch-triggering mechanism, however, is beyond the scope of the model in its present form. It is rather brought into it by hand. Therefore, the model by itself does not allow us to predict the sizes and the waiting times of the glitches but these questions could be addressed in the existing core-located glitch models. The main new point that we assume is the possible weakness of cg-interaction. This assumption does not affect the “normal" fraction spin-down relative to the pinned superfluid, which is usually supposed to be the key process for glitch triggering.

Unfortunately, our model in its present form does not allow us to obtain the whole range of observed recovery fraction values. According to estimation (157) it is of the order of unity. It is not so bad for young pulsars but the older ones demonstrate a wide range of recovery parameters (Yu et al., 2013). However, we believe that this discrepancy arises due to oversimplification of the model. In particular, the injection of the angular momentum from the pinned superfluid into the “normal" fraction requires more self-consistent consideration. Thus, further research is required.

Acknowledgments

The authors are grateful to A I Tsygan, K Y Kraav, M V Vorontsov, and A N Yurkin for fruitful discussions. The authors are also grateful to Dr. Ashton for many useful comments.