# A phenomenological wobbling model for isolated pulsars and the braking index

E. C. A. Araujo, V. A. de Lorenci, P. Peter, L. S. Ruiz

► To cite this version:

E. C. A. Araujo, V. A. de Lorenci, P. Peter, L. S. Ruiz. A phenomenological wobbling model for isolated pulsars and the braking index. Monthly Notices of the Royal Astronomical Society, 2024, 527, pp.7956-7964. 10.1093/mnras/stad3531. insu-04851596

HAL Id: insu-04851596

https://insu.hal.science/insu-04851596v1

Submitted on 20 Dec 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# A phenomenological wobbling model for isolated pulsars and the braking index

E. C. A. Araujo<sup>1</sup>,<sup>2</sup> V. A. De Lorenci<sup>2,3★</sup> P. Peter<sup>3</sup> and L. S. Ruiz<sup>4,5</sup>

Accepted 2023 November 14. Received 2023 November 12; in original form 2023 August 2

## ABSTRACT

An isolated pulsar is a rotating neutron star possessing a very high magnetic dipole moment, thus providing a powerful radiating mechanism. These stars lose rotational energy *E* through various processes, including a plasma wind originating from a highly magnetized magnetosphere and the emission of magnetic dipole radiation (MDR). Such phenomena produce a time decreasing angular velocity $\Omega(t)$ of the pulsar that is usually quantified in terms of its braking index. Although these mechanisms are widely acknowledged as the primary drivers of the spin evolution of isolated pulsars, it is plausible that other contributing factors influencing this effect have yet to be comprehensively investigated. Most of young isolated pulsars present a braking index different from that given by the MDR and plasma wind processes. Working in the weak field (Newtonian) limit, we take in this work a step forward in describing the evolution of such a system by allowing the star's shape to wobble around an ellipsoidal configuration as a backreaction effect produced by its rotational deceleration. It is assumed that an internal damping of the oscillations occurs, thus introducing another form of energy loss in the system, and this phenomenon may be related to the deviation of the braking index from the models based on $\dot{E} \sim -\Omega^4$ predictions. Numerical calculations suggest that the average braking index for typical isolated pulsars can be thus explained by a simple phenomenological model.

**Key words:** methods: analytical—methods: statistical—stars: neutron—stars: oscillations (including pulsations)—pulsars: general.

## 1 INTRODUCTION

The identification of isolated pulsars (Hewish et al. 1968) with rotating neutron stars presenting a high surface magnetic field was suggested long ago (Gold 1968), where predictions about their spin evolution were also anticipated. Shortly thereafter, the pulsar detected in the Crab nebula was measured (Richards 1968) to slow down. Commonly accepted models describing the loss of energy leading to the observed spin evolution of pulsars are based on the emission of magnetic dipole radiation (MDR) (Pacini 1967, 1968; Gunn & Ostriker 1969; Hamil et al. 2015; Lyne et al. 2015; Shaw et al. 2022) and plasma wind (Goldreich & Julian 1969; Ruderman & Sutherland 1975; Beskin, Gurevich & Istomin 1984; Xu & Qiao 2001; Spitkovsky 2006; Yue et al. 2007). As these processes carry away angular momentum, the star's angular velocity decreases with time, leading to a slowing-down behaviour that has been detected in virtually all known pulsars and in some cases monitored for more than 50 yr (Hobbs, Lyne & Kramer 2010; Lyne et al. 2015; Namkham, Jaroenjittichai & Johnston 2019; Parthasarathy et al. 2019, 2020; Shaw et al. 2022).

When only these processes are considered, the rate at which the rotational energy is radiated away from the pulsar is described by the power-law decay (Shapiro & Teukolsky 1983)

$$ \dot{E} = -\beta I \Omega^4, \quad (1) $$

where $\beta$ is a parameter whose physical meaning depends on the adopted model, *I* is the moment of inertia of the star, and $\Omega = 2\pi v$ is its angular velocity, with $v$ its rotation frequency. In equation (1) and in what follows, a dot over a physical quantity represents its time derivative. It should be noticed that depending on the choice for the acceleration potential in wind models (Ruderman & Sutherland 1975; Yue et al. 2007; Li et al. 2014), different power rates for the loss of rotational energy are possible (Zhang et al. 2022). Here, the analysis is applicable to the models obeying the power law described by equation (1).

The relationship between $\dot{E}$, $\Omega$, and *I* is given by the torque equation $\dot{\Omega} = (I\Omega)^{-1}\dot{E}$, leading to

$$ \dot{\Omega} = -\beta \Omega^3. \quad (2) $$

The parameter $\beta$ can be evaluated for the popular MDR and wind mechanisms, with similar orders of magnitude. For instance, in the case of the MDR model (in units where $c = 1 = \mu_0$, with $\mu_0$ the magnetic vacuum permeability) it follows that $\beta = (2\mu^2\sin^2\alpha)/(3I)$,

© 2023 The Author(s).

<footnote n="1"><sup>1</sup>Institute of Myology, 47 Bd de l'Hôpital, F-75013 Paris, France</footnote>

<footnote n="2"><sup>2</sup>Instituto de Física e Química, Universidade Federal de Itajubá, Itajubá, Minas Gerais 37500-903, Brazil</footnote>

<footnote n="3"><sup>3</sup>GR&ECO-Institut d'Astrophysique de Paris, CNRS and Sorbonne Université, UMR 7095 98 bis Boulevard Arago, F-75014 Paris, France</footnote>

<footnote n="4"><sup>4</sup>Instituto de Matemática e Computação, Universidade Federal de Itajubá, Itajubá, Minas Gerais 37500-903, Brazil</footnote>

<footnote n="5"><sup>5</sup>CFisUC, Departamento de Física, Universidade de Coimbra, P-3004-531 Coimbra, Portugal</footnote>

<footnote>*E-mail: delorenci@unifei.edu.br</footnote>
<figure><img src="image_2.png" /><figcaption>**Figure 1.** The pulsar geometric configuration: the neutron star is rotating with angular velocity $\Omega(t) = \dot{\phi}$ around the z axis while its magnetic axis is inclined by an angle $\alpha$ leading to precession of the field lines. The rotation induces a deviation from sphericity so the original sphere of radius $R$ (dashed grid) turns into an ellipsoid whose semi-axes depend on time through the function $b(t)$ related to the ellipticity.</figcaption></figure>

where $\mu$ is the magnitude of the magnetic dipole moment $\mu$ of the star, making an angle $\alpha$ with its rotation axes, determined by $\Omega$, as illustrated in Fig. 1. In our system of units, $\mu$ is related with the magnetic field $B$ through $\mu = \frac{1}{2}BR^3$. It is worth noting that, apart from an overall numerical factor, the wind model yields the same behaviour $\beta \sim B^2R^6/I$, so that the phenomenological model below does not depend on the dominant energy loss mechanism.

More generally, the different causes behind the slowing-down phenomenon can be encapsulated by the power-law formula (Goldwire & Michel 1969) $\dot{\Omega} = -K\Omega^n$, where $n$ is the so-called braking index of the pulsar. It is usually defined as if $K$ were constant, namely

$$n \doteq \frac{\Omega\dot{\Omega}}{\Omega^2}, \qquad (3)$$

and is actually constant only for constant $K$.

In the context of the models mentioned above, assuming $K = \beta$ constant, one naturally gets $n = 3$. Other models have been suggested, in which for instance the energy loss is dominated by gravitational waves radiation; these imply different values for the braking index, e.g. $n = 5$ in the case of Ref. Ostriker & Gunn (1969).

Precise measurements of the braking index of several isolated pulsars have revealed values that are consistently different than that predicted by assuming a constant $K$ (see for instance Refs. Gourgouliatos & Cumming 2015; Lyne et al. 2015; Parthasarathy et al. 2020; Lower et al. 2021). Consequences of assuming $K$ as time-dependent function (Blandford & Romani 1988) have been examined (Magalhaes, Miranda & Frajuca 2012; Lyne et al. 2015), and depending on the way $K$ evolves in time, the braking index can be $n \neq 3$. Furthermore, assuming (Hamil et al. 2015) $I = I[\Omega(t)]$ only is not enough to explain the present observational data for the known isolated pulsar. It was recently shown (Hamil, Stone & Stone 2016) that a MDR-based model with a time-dependent inclination angle $\alpha = \alpha(t)$ would be a possible way to explain the phenomenon.

All the models described above are based on the simplifying assumption that the kinetic energy $E_K$ of the body is only due to its rigid rotation, i.e. $E_K = \frac{1}{2}I\Omega^2$. However, a neutron star cannot be strictly considered as rigid and even though the rotation is very slow from the point of view of relativistic effects, with typical surface velocities of the order of$^1 v \simeq (10^{-4} - 10^{-2})c$, the shape should be allowed to depend on time as this rotation may naturally induce a flattening of the poles. In such a scenario, the kinetic energy acquires new contributions that need to be taken into account.

It is the purpose of this work to introduce a more complete description taking into account the variation of the internal potential energy of a self-gravitating body (Ragazzo & Ruiz 2015). In particular, we explore the consequences of allowing the star shape to evolve in time, under its coupling with the rotation of the body. It is assumed that energy can be lost in this process, a phenomenon that could be relevant in the explanation of the measured slowing-down of isolated pulsars. Embedding our model in a general relativistic context goes beyond the scope of this work, and as we consider the quasi-rigid rotation of the star, we restrict attention to the Newtonian weak field limit; given the slow velocities involved, we expect this approximation to be meaningful.

In the next section, the basic assumptions of the model are described and the coupled system of non-linear differential equations governing the evolution of the star are derived. The results of our numerical analysis are presented in Section 3, where some suggestive solutions are studied. In particular, we show that a pulsar evolution with $n < 3$ can easily be reproduced. A comparison between our results and the available data describing the behaviour of the Crab pulsar is given in Section 4, before a few final and concluding remarks in Section 5.

## 2 THE MODEL

Suppose the star is described as a mass distribution that is slightly deformed compared to a spherical body and is rotating around the z-axis with time-dependent angular velocity $\Omega(t)$. We assume the volume of the star to be that of the non-rotating sphere $\frac{4}{3}\pi R^3$, its actual shape being ellipsoidal with two equal semi-axes in the $(x, y)$-plane slightly larger than the sphere radius, i.e. $(1+b)R$ (see Fig. 1). Even for the large velocities involved in a rotating neutron star, we do not expect large deviations from sphericity and thus demand that $b \ll 1$ (Baym & Pines 1971; Cutler, Ushomirsky & Link 2003; Rencoret, Aguilera-Gómez & Reisenegger 2021). As a result, the volume of the ellipsoid matches that of the sphere to second order in $b$ provided the semi-axis in the direction of rotation is $(1 - 2b)R$. The location of an arbitrary element of mass $dm$ in the body will thus be described by the set $R' = \{X', Y', Z'\}$ such that

$$\left(\frac{X'}{1+b}\right)^2 + \left(\frac{Y'}{1+b}\right)^2 + \left(\frac{Z'}{1-2b}\right)^2 \le R^2.$$

As discussed in the introduction, we expect the quantity $b$ may not necessarily be constant, so we anticipate that $b = b(t)$. It is convenient to use a coordinate system related to that of the embedding sphere, i.e. the set $\{x, y, z\}$ such that $x^2 + y^2 + z^2 = r^2$, with $r \le R$. That is, we implement the coordinate transformation $X' = (1+b)x$, $Y' = (1+b)y$, and $Z' = (1-2b)z$. The volume element, as expected, is $dV = dX'dY'dZ' = dxdydz + O(b^2)$.

<footnote n="1">¹The Crab is among the fastest rotating isolated pulsars with known braking index, with velocity $v_{crab} \sim 1.2 \times 10^{-2}c$.</footnote>
Let us implement a rotation $R_{\hat{z}}(\varphi)$ to the body about the z-axis by an angle $\varphi$. The position $\mathbf{r}'$ of the mass element $dm$ is then given by the rotation applied to its location $\mathbf{R}'$, i.e. $\mathbf{r}'(\varphi, b) = R_{\hat{z}}(\varphi) \cdot \mathbf{R}'$. One therefore gets

$$ \mathbf{r}' = \underbrace{\begin{pmatrix} \cos \varphi & \sin \varphi & 0 \\ -\sin \varphi & \cos \varphi & 0 \\ 0 & 0 & 1 \end{pmatrix}}_{R_{\hat{z}}(\varphi)} \underbrace{\begin{pmatrix} 1+b & 0 & 0 \\ 0 & 1+b & 0 \\ 0 & 0 & 1-2b \end{pmatrix}}_{\mathbf{R}'} \begin{pmatrix} x \\ y \\ z \end{pmatrix}, \quad (4) $$

whose modulus we denote by $r'$.

Since we consider a pulsar, i.e. a rotating neutron star, one now needs to assume from that point on that both the rotation angle $\varphi$ and the flattening of the poles depend on time. The absolute velocity, $v = \|\mathbf{dr}'/dt\|$, of each mass element in such system is given by

$$ v = \sqrt{(x^2 + y^2 + 4z^2)\dot{b}^2 + (x^2 + y^2)\Omega^2(1+b)^2}, \quad (5) $$

where $\Omega = \dot{\varphi}$ is the angular velocity of $dm$ around the z-axis. The kinetic energy $E_K$ of the star is obtained by integrating $\frac{1}{2}v^2dm$ over the whole body volume, leading to

$$ E_K = \frac{1}{2} \dot{b}^2 \int_V (x^2 + y^2 + 4z^2) \rho(r) dV + \frac{1}{2} \Omega^2 (1+b)^2 \int_V (x^2 + y^2) \rho(r) dV, \quad (6) $$

where $\rho(r)$ is the mass-density; it should be noted at this point that this function depends here only on the radial distance and not on time as we assumed a constant star volume. In that case, the fluid making the star satisfies the continuity equation (mass conservation in a fixed volume), namely $\partial_t \rho + \nabla(\rho v) = 0$, with the velocity being orthoradial but depending only on $r$, so that $\nabla(\rho v) = 0$, and therefore $\partial_t \rho = 0$.

The moment of inertia $I$ of the neutron star, seen as an idealized spherical mass distribution rotating about the z-axis, is defined by

$$ I \doteq \int (x^2 + y^2) \rho(r) dV, $$

which, because of the spherical symmetry, is also expressible as

$$ I = 2 \int \rho x^2 dV = 2 \int \rho y^2 dV = 2 \int \rho z^2 dV. \quad (7) $$

For the spherical unperturbed configuration, we assume a constant density $\rho(r) = \langle\rho\rangle$ in order to approximate the moment of inertia by $I \sim \frac{2}{5}MR^2$, so that $\beta \sim B^2R^4/M$, with a numerical prefactor depending on the model.

The kinetic energy now reads

$$ E_K = \frac{3}{2} I \dot{b}^2 + \frac{1}{2} I (1 + 2b) \Omega^2 + \mathcal{O}(b^2). \quad (8) $$

As the system evolves, the body will be allowed to oscillate. Its potential energy $E_p$ can be expanded about $b=0$, as

$$ E_p \approx E_p(0) + \frac{1}{2} \kappa b^2, \quad (9) $$

where we have used that the potential energy is minimized for the spherical configuration, so that $(\partial E_p/\partial b)_{b=0} = 0$. In equation (9), we noted the elastic constant as

$$ \kappa = \left( \frac{\partial^2 E_p}{\partial b^2} \right)_{b=0} = 3I\gamma, \quad (10) $$

thereby defining the coefficient $\gamma$. The seemingly arbitrary relationship between $\gamma$ and the elastic constant $\kappa$ in (10) is justified a posteriori to lead to a simplification of the corresponding equation of motion (16a) for $b(t)$. The leading contribution to the elastic constant

can be obtained by assuming the spherical approximation, which leads to $\kappa \approx 24GM^2/(5R)$, such that $\gamma \approx 4GM/R^3$.

Neglecting higher order terms, the Lagrangian of the system reads, up to a constant,

$$ \mathcal{L} = \frac{1}{2} I [3\dot{b}^2 + (1 + 2b)\Omega^2 - 3\gamma b^2]. \quad (11) $$

It should be noticed at this point that the $\mathcal{O}(b^2)$ term in the kinetic contribution, namely $\frac{1}{2}Ib^2\Omega^2$, has been neglected because it is assumed to be very small when compared to $\gamma Ib^2$ coming from the potential energy contribution. This corresponds to assuming $\Omega^2 \ll GM/R^3$, a condition that is related to the slow rotation Newtonian hypothesis, satisfied for the physical system under consideration.

The Euler-Lagrange equations stemming from (11) must be supplemented by dissipation terms (Baillieul & Levi 1987; Ragazzo & Ruiz 2015, 2017; Caressa & Bersani 2020), in order to account for the radiation. For a system with degrees of freedom $q_i$ and Lagrangian $\mathcal{L}(q_i, \dot{q}_i)$, one introduces a so-called dissipation function $D(q_i, \dot{q}_i)$ such that $\dot{q}_i D_{,q_i}$ represents the rate of energy loss per unit time (power dissipated by the 'force' $D_{,q_i} = \partial D/\partial \dot{q}_i$). The Euler-Lagrange equations are modified to include this new force as a (negative) source term, namely

$$ \frac{d}{dt} \left( \frac{\partial \mathcal{L}}{\partial \dot{q}_i} \right) - \frac{\partial \mathcal{L}}{\partial q_i} = - \frac{\partial D}{\partial \dot{q}_i}. $$

This interpretation stems from the fact that defining the energy in the usual way through $H = \dot{q}_i p^i - \mathcal{L}$, with the momenta given by $p^i = \partial \mathcal{L}/\partial \dot{q}_i$, and using the above equation, one finds the rate of energy loss as $dH/dt = -\dot{q}_i D_{,q_i}$, as announced.

In the case at hand, introducing the dissipation function $D(\dot{\varphi}, \dot{b})$, the modified equations of motion read

$$ 3I\ddot{b} + 3I\gamma b - I\Omega^2 = -\frac{\partial D}{\partial \dot{b}}, \quad (12) $$

$$ \frac{d}{dt} [I(1 + 2b)\Omega] = -\frac{\partial D}{\partial \dot{\varphi}}, \quad (13) $$

with $D$ taking into account the radiation emission losses and the damping of the body oscillations.

Our simplified model relies on internal dissipation processes associated with the quadrupole moment tensor, and we demand that the oscillations have a small amplitude such that they should remain linear in their time derivative. These requirements can be achieved with the following prescription for the dissipation function (Caressa & Bersani 2020)

$$ D = \frac{1}{4} \beta I \Omega^4 + \frac{3}{2} \sigma I \dot{b}^2, \quad (14) $$

thereby defining our final phenomenological parameter $\sigma$.

Now, defining the total energy $E = E_K + E_p$, and using the above results, it is straightforward to evaluate the energy losses, namely

$$ \dot{E} = -\dot{\varphi} \frac{\partial D}{\partial \dot{\varphi}} - \dot{b} \frac{\partial D}{\partial \dot{b}} = -\beta I \Omega^4 - 3\sigma I \dot{b}^2, \quad (15) $$

which is the equation that governs the energy balance of the system.

In a scenario where MDR and plasma wind are the only processes behind the loss of energy of a pulsar, with relative strengths to be defined by the explicit model under consideration, equation (1) would hold and the external torque $\tau_{\text{ext}} = -\beta I \Omega^3$ would be the only responsible for the star slowdown. However, in the more complete scenario under investigation in this work, the evolution of the system is governed by the set of coupled equations of motion given by equations (12) and (13) which, after inserting equation (14), can be presented in the more compact form as
$$ \dot{b} + \sigma \dot{b} + \gamma b = \frac{1}{3} \Omega^2, \qquad (16a) $$

$$ \dot{\Omega} = -\frac{2\Omega\dot{b}}{(1+2b)} - \frac{\beta}{(1+2b)}\Omega^3. \qquad (16b) $$

From the point of view of physical units, the parameters $\beta$, $\gamma$, and $\sigma$ are expressed respectively in $s$ ($[\beta] = T$), $s^{-2}$ ($[\gamma] = T^{-2}$), and $s^{-1}$ ($[\sigma] = T^{-1}$).

The dissipation contribution in the equation of motion (16a) for $b(t)$ naturally imposes a time-evolving equilibrium point for this function, which is given by $\Omega^2/(3\gamma)$. It is worth mentioning that had the $O(b^2)$ kinetic term been kept in equation (11), this equilibrium point would have been corrected to $b(t) \approx [\Omega^2/(3\gamma)][1 + \Omega^2/(3\gamma)]$, a much higher order correction: our initial approximation to neglect $O(b^2)$ contribution to the kinetic energy is justified for the physical systems discussed in this work.

Before closing this section, a few words about angular momentum conservation are in order. First, in our model, the quantity $\mathcal{I} \doteq I(1+b)^2$ is identified as the time-dependent effective moment of inertia of the body, thus making equation (16b) the equation of motion relating the total angular momentum $L = \mathcal{I}\Omega$ with the external torque produced by the radiation emission. Naturally, in the absence of external torque, the angular momentum is a conserved quantity, i.e. when $\tau_{\text{ext}} = -\partial D/\partial \dot{\phi} = 0$. As expected, internal processes, as those described by the second term in the rhs of equation (14), do not interfere with the angular momentum conservation law. Finally, it should be emphasized that when the moment of inertia is allowed to vary with time, $\frac{1}{2}\mathcal{I}\Omega^2$ will not be the only contribution to the kinetic energy of the body, as clearly emphasized by equation (8). The time evolution of the angular momentum is governed by $\Omega(t)$, and also by $\mathcal{I}(t)$ through $b(t)$. These functions are solutions of the coupled differential equations of motion (16) that naturally follow from the Lagrangian method.

### 3 MODELLING A PULSAR SLOWDOWN

Quantities like the mass of the pulsar, its radius, or the strength of the field at its magnetic pole are not known with great precision, and these values can also be model dependent. For instance, the mass of a pulsar, like Crab and others, is usually taken to be approximately $1.4 M_{\odot}$, with $M_{\odot}$ the solar mass (Shapiro & Teukolsky 1983). The goal of this work is to test if our theoretical model is able to produce acceptable solutions to the problem of pulsars slowdown, i.e. if a braking less than 3 is possible when the oscillations described by $b(t)$ are taken into account.

Using the results obtained in the last section, the parameters $\beta$ and $\gamma$ can be conveniently expressed in terms of $M_{\odot}$ and the typical values for the radius and the magnetic dipole field of a certain class of known pulsars, namely

$$ \beta \approx 1.485 \times 10^{-18} \Upsilon_{\text{m}} \left(\frac{M_{\odot}}{M}\right) \left(\frac{B}{10^{8}\text{T}}\right)^{2} \left(\frac{R}{10\text{km}}\right)^{4} \text{s}, \qquad (17) $$

$$ \gamma \approx 5.307 \times 10^8 \left( \frac{M}{M_{\odot}} \right) \left( \frac{10 \text{km}}{R} \right)^3 s^{-2}, \qquad (18) $$

where we have used the expressions for $I$ and $\mu$ assuming a spherical star, and $\Upsilon_{\text{m}}$ is a numerical value that depends on the specific choice of MDR or plasma wind models. For instance, MDR model leads to $\Upsilon_{\text{MDR}} = (5/12)\sin^2\alpha$, while a plasma-wind model (Xu & Qiao 2001; Li et al. 2014) leads to $\Upsilon_{\text{w}} = (5/4)(\Delta\phi/\Delta\Phi)\cos^2\alpha$, where $\Delta\phi$ denotes the acceleration electric potential whose maximum value

is $\Delta\Phi$. It is interesting to notice that assuming $\Delta\phi \approx \Delta\Phi/3$ results in $\Upsilon_{\text{w}} = (5/12)\cos^2\alpha$, which coincides with $\Upsilon_{\text{MDR}}$ for $\alpha = \pi/4$.

In the subsequent numerical calculations and for definiteness, we assume specific values for the neutron star model. To begin with, we assumed the mass $M$ to be $M \rightarrow 1.4 M_{\odot}$. Now, conveniently setting $\gamma \approx 1.583 \times 10^8 s^{-2}$ in the simulations and using equation (18), it follows that the star radius in this model is $R = 1.674 \times 10^4$ m. Similarly, by setting $\beta \approx 3.542 \times 10^{-16}$ s and plugging in the above values of $M$ and $R$ into equation (17), using the MDR model with a misalignment angle $\alpha = (\pi/4)$ rad to estimate the parameter $\Upsilon_{\text{m}}$, namely $\Upsilon_{\text{MDR}} = 5/24$, we obtain a magnetic dipole generating a field amplitude $B = 1.428 \times 10^9$ T, which denotes the magnitude of the field at the pole of the star (Shapiro & Teukolsky 1983).

As already discussed, the value of $\Upsilon_{\text{m}}$ is the only one which explicitly depends on the specific model describing the radiative processes involved in the loss of rotational energy of the pulsar. Given the degeneracy in the choice of the parameters, a different value for the model-dependent parameter $\Upsilon_{\text{m}}$ can be associated with the same value of $\beta$. The remaining parameter $\sigma$ is associated to the dissipation processes during the oscillations of the quadrupole moment of the body, and can be adjusted in the numerical calculations in order to obtain the braking index of the pulsar.

Following the model described in the previous section, the evolution of the system is governed by the coupled non-linear differential equations given by equations (16). The initial angular velocity is set to be $\Omega(0) = 188.5 \text{ rad} \cdot \text{s}^{-1}$. The initial deformation of the body, $b(0)$, is assumed to be the equilibrium value of $b(t)$ in equation (16b), i.e. $b(0) = \Omega(0)^2/(3\gamma) \approx 7.482 \times 10^{-5}$, for which $\dot{b}(0)$ was set to zero.

In a fashion similar to that present in the analysis in the existing literature for the calculation of the braking index of the Crab pulsar (Lyne et al. 2015), we consider here the following method, that can be applied to both simulated or measured data:

(i) Let $(t, \Omega) = (t_k, \Omega_k)$, for $k = 1, \dots, N$, denote the complete time-series for the angular velocity of the pulsar;

(ii) For each $j$ such that $1 < j \le N$, the cumulative braking index at time $t_j$ is computed by fitting the points $(t_k, \Omega_k)$, where $k = 1, \dots, j$, with the 3rd-degree polynomial

$$ P_{c,j}(t) = \sum_{k=0}^{3} a_{k,j} (t - t_{c,j})^k, $$

where $t_{c,j} = \frac{1}{2}(t_j - t_0)$ is the half time of the interval $[t_0, t_j]$, and the corresponding braking index at $t_j$ is given by equation (3), reading here

$$ n_c(t_j) = \frac{\ddot{P}_{c,j}(t_{c,j}) P_{c,j}(t_{c,j})}{\dot{P}_{c,j}^2(t_{c,j})} = 2 \frac{a_{2,j} a_{0,j}}{a_{1,j}^2}; \qquad (19) $$

(iii) The local braking index at time $t_j$ must be determined over each data subset with a fixed size $\Delta N \in \mathbb{N}$ as illustrated below:

<figure><img src="image_3.png" /></figure>

In this case, for each $\Delta N < j \le N$, the local braking index at time $t_j > t_{\Delta N}$ is computed by fitting the points $(t_k, \Omega_k)$, where $k = j - \Delta N, \dots, j$, with the 3rd-degree polynomial

$$ P_{L,j}(t) = \sum_{k=0}^{3} b_{k,j} (t - t_{L,j})^k, $$

where $t_{L,j} = \frac{1}{2}(t_j - \Delta N + t_j)$ is the centre of the interval $[t_j - \Delta N, t_j]$, and the corresponding local braking index is again given by equation
<figure><img src="image_4.png" /><figcaption>**Figure 2.** Cumulative [upper panel, equation (19)] and local [lower panel, equation (20)] braking indices calculated from the simulation data for some representative values of the dissipation parameter $\sigma$. Note that there is no direct relationship between the magnitude of the dissipation process and the order of the braking indices when small variations of $\sigma$ are considered. However, on average, a more intense dissipation process (larger $\sigma$ values) leads to smaller values for the braking index.</figcaption></figure>

(3), namely

$$n_{\mathrm{L}}(t_{\mathrm{j}})=\frac{\ddot{P}_{\mathrm{L}, \mathrm{j}}(t_{\mathrm{L}, \mathrm{j}}) P_{\mathrm{L}, \mathrm{j}}(t_{\mathrm{L}, \mathrm{j}})}{\dot{P}_{\mathrm{L}, \mathrm{j}}^{2}(t_{\mathrm{L}, \mathrm{j}})}=2 \frac{b_{2, \mathrm{j}} b_{0, \mathrm{j}}}{b_{1, \mathrm{j}}^{2}} \quad (20)$$

An extrapolation-algorithm, based on the explicit midpoint rule, with step size control and order selection [see Section II.9 from Ref. (Hairer, Nørsett & Wanner 1993)] was used to numerically integrate the coupled system described by equations (16a) and (16b), leading to the results depicted in Figs 2 and 3. The integration spans a time window of about 5 yr, which was enough to obtain solutions with stable braking indices. In fact, after a short time of instability, $n(t)$ eventually behaves as a slowly evolving function of time, as it can be confirmed by direct inspection of Fig. 2, where some solutions presenting positive braking indices were selected.

The magnitude of the dissipation process associated with the quadrupole oscillations is dominant in determining the behaviour of the braking index of the system. Processes for which $\sigma$ is of the order of 0.01 Hz lead to braking indices around $n = 3$, which is the expected result when the pulsar's rotational energy is taken away only by means of magnetic dipole radiation. However, for higher values of $\sigma$, richer scenarios appear, as shown in Fig. 2. In particular, when $\sigma \approx 2$ Hz, the solutions exhibit braking indices around 2.5. Small variations of $\sigma$ lead to different solutions for $n(t)$. On the other hand, this function does not seem to be very sensitive to small variations of the other parameters. The local behaviour of the braking index, depicted in the lower panel of Fig. 2, was obtained using a moving average ($\Delta N$) of 400 d, which explains why it starts after the cumulative index (upper panel). When a sufficiently high dissipation process is taken into account, the simulations suggest that even negative braking indices are possible solutions. This is an aspect that deserves further examination.

<figure><img src="image_5.png" /><figcaption>**Figure 3.** Time evolution of $b(t)$ as the pulsar slows down. The upper inset depicts the behaviour of the normalized $b(t)$, obtained when its time-evolving equilibrium point $\Omega(t)^2/(3\gamma)$ is subtracted, thus showing that it is a highly oscillating function in time but on a small amplitude scale. The lower inset presents a close-up of the normalized curve (about day 2000) in a scale of milliseconds, showing that the oscillations occur with a well defined frequency, approximately the natural frequency $\sqrt{\gamma}/(2\pi)$ of the harmonic oscillator described by equation (12).</figcaption></figure>

The behaviour of the angular velocity is very similar for all solutions examined in Fig. 2. If the curves corresponding to the angular velocities for these models were included in a same plot, almost no visual difference would be seen. In fact, it can be shown that for any instant of time in the simulations, the difference between the angular velocities of any of these solutions is smaller than $10^{-5}$ Hz.

The behaviour of $b(t)$ for the model with $\sigma = 1.995 \text{ s}^{-1}$ is shown in Fig. 3. In the plot scale it looks like a slowly decreasing monotonic function of time. However, a more detailed examination shows that $b(t)$ is a highly oscillatory function around the time-dependent equilibrium point $\Omega^2(t)/(3\gamma)$, as highlighted in the inserts. Indeed, due to their mutual coupling, both $b(t)$ and $\Omega(t)$ decompose into a slow monotonically decreasing component and a fast (and tiny) oscillatory component; the slow component can be extracted out by calculating the difference $b(t) - \Omega^2(t)/(3\gamma)$. Thus, as the system loses energy by means of radiation emission and oscillation damping, as effectively described by equation (15), it will slow-down its rotation frequency and also the amplitude of the oscillations.

The evolution of the rotation frequency corresponding to the model with $\sigma = 1.995 \text{ s}^{-1}$ is depicted in Fig. 4(a), which is a solution presenting a braking index of approximately 2.5. The other panels in the figure show the residuals of the first, second, and third order, which were obtained following the standard procedure [see for instance the analysis for the Crab pulsar (Lyne et al. 2015)]: the data set is fitted by means of a k-degree polynomial, which can be written as $v(t) = \sum_{i=0}^{k} c_i (t - t_0)^i + \delta v^{(k)}$, where $t_0$ is chosen, for instance, to be the medium time of the data set, the coefficients $c_i$ are obtained by the fitting procedure, and the time-dependent function $\delta v^{(k)}$ is the k-th order residual obtained when the k-th order fitting polynomial is subtracted from the data. For instance, Fig. 4(b) depicts the first-order residual $\delta v^{(1)} = [v(t)]_{\text{data}} - [c_0 + c_1(t - t_0)]$, where $t_0 = 9.299 \times 10^7 \text{ s} \approx 1076 \text{ d}$, $c_0 \approx 29.97 \text{ Hz}$, and $c_1 \approx -3.757 \times 10^{-10} \text{ s}^{-2}$.

As it is well known, the gradual deceleration of the rotational frequency of an isolated pulsar experiences sudden and occasional changes, followed by a period of partial recovery to the previous
<figure><img src="image_6.png" /><figcaption>**Figure 4.** Rotational frequency (a) and residuals of the first (b), second (c), and third-order (d) shown as functions of time (in days) for the simulation data corresponding to a dissipation parameter $\sigma = 1.995 \text{ s}^{-1}$. The cumulative braking index corresponding to this simulation data, after about 6 yr of integration time, achieves a value of approximately 2.5, as can be inferred by direct inspection of Fig. 2. The initial frequency in this simulation was 30.00 Hz.</figcaption></figure>

ous regular rotational rate. This phenomenon is referred to as a glitch and is attributed to complex processes taking place within the star's internal fluid (Lyne, Pritchard & Graham Smith 1993; Antonopoulou, Haskell & Espinoza 2022). If we examine the residuals at the second and third order in Fig. 4, we can observe irregularities that cannot be attributed to numerical errors. They can be interpreted as very short moments in time during which the angular velocity is suddenly changed before the star returns to its original state. That could be interpreted as micro-glitches: a more realistic model would describe for instance various rotating shells, all of which would be subject to equations similar to those presented here and somehow interacting. Could such a more elaborate model enhance this phenomenon to the level of the observed glitches?

<figure><img src="image_7.png" /><figcaption>**Figure 5.** Rotation frequency and residuals of the Crab pulsar, according to data collected from 1982 February to 2022 August (http://www.jb.man.ac.uk/pulsar/crab.html). The starting time in this figure corresponds to MJD 45015, for which the measured frequency was $\nu = 30.0592241133 \text{ Hz}$. The residuals of second and third orders, depicted in panels (c) and (d), respectively, clearly show a rich glitch activity of the pulsar in this period.</figcaption></figure>

# 4 A NOTE ABOUT GLITCHES: THE BEHAVIOUR OF THE CRAB PULSAR

Having described our simple model, one wants to compare with the existing data relevant to the dynamical range under investigation. The best example one can think of is provided by the enormous amount of data available concerning the Crab pulsar.

The Crab pulsar (PSR B0531+21) is an isolated neutron star whose angular velocity deceleration has been measured since the 1970s (Gullahorn et al. 1977). Monthly spaced pulsar timing measurements have been taken by Jodrell Bank Observatory since 1982 (Lyne et al. 1993). In Fig. 5(a), the rotation frequency measured for the Crab pulsar is shown as a function of time, from MJD 45 015 (1982 February 15) to MJD 59 806 (2022 August 15) (the data are available at
<figure><img src="image_8.png" /><figcaption>**Figure 6.** Braking index as function of time for the Crab pulsar. Here the index is calculated using the cumulative and local procedures discussed in the last section, as depicted in upper and lower panels, respectively. A moving average with $\Delta N = 30$ was here used to obtain the local index.</figcaption></figure>

http://www.jb.man.ac.uk/pulsar/crab.html). The residuals are shown from top to bottom [panels (b) to (d)], where the coefficients of the third-order polynomial fitting are $t_0 = 6.389712 \times 10^8$ s = 7395.5 d, $c_0 \approx 31.36$ Hz, $c_1 \approx -4.205 \times 10^{-10}$ s$^{-2}$, $c_2 \approx 6.728 \times 10^{-21}$ s$^{-3}$, and $c_3 \approx -1.260 \times 10^{-31}$ s$^{-4}$.

This system has occasional glitches, in which the star is spinnéd up for a short period of time and returns to its former rotation frequency within an interval of about 20 d. The glitch activities can be clearly seen in Fig. 5(c) and 5(d). They contribute massively to the cumulative braking index. The cumulative and local indices agree with the constant value $n \approx 2.5$ if and only if they are computed in an interval not containing a glitch. However, the cumulative braking index decreases monotonically in time as multiple intervals are included in the data set so that its final value reaches $n \approx 2.3$ (Lyne et al. 2015).

Cumulative and local behaviours of the braking index with time are shown in Fig. 6, upper and lower panels, respectively. In particular, the local index as a function of time (in months) is shown in lower panel of Fig. 6, where the moving average was calculated using a window containing 30 successive measurements. It can be seen that after each glitch, the braking index returns approximately to the value it had before the glitch, and this happens in less than a month. However, its influence in the local braking index calculation goes way longer, an aspect that is dependent on the choice of $\Delta N$.

It is noteworthy that in the simulations leading to Fig. 4, the chosen values for the physical parameters of the system are comparable to those typically employed to describe the Crab pulsar. Particularly, if we use the data from Crab timing measurements (Lyne et al. 1993) and adjust a cubic polynomial to it, as given by equation (2), we obtain $\beta = 3.573 \times 10^{-16}$ s, which is approximately the value used in the simulations described in Section 3. Additionally, the obtained local braking index approximately aligns with the measured value for the Crab pulsar. This parallel suggests that the simulated system shares similarities with the Crab pulsar, potentially indicating comparable underlying dissipative mechanisms at work, which would lead us

to suggest that $\sigma \approx 2$s is expected in systems like Crab. However, in order to establish a conclusive connection between the simulated system and the Crab pulsar, further simulations and comparisons would be necessary.

## 5 FINAL REMARKS

In this work, we explored the idea that as radiation is emitted by an isolated pulsar, its energy is continuously driven away, causing a slow-down of its spin, and a possible modification of the shape of the star. Although the radiation emission (MDR or plasma wind) is largely considered in the literature, adding a perturbation in its ellipsoidal shape by means of small oscillations has never been considered. As the star cannot be strictly rigid, oscillations are expected: they are produced almost in a stationary regime as it is linked to the spin slow-down process. These oscillations must be dissipated by internal phenomena leading to a secondary form of energy loss by the star. The possible mechanisms behind such dissipation of energy were not considered in details in this work. Instead, it was assumed that the effect is described by an effective damping process that is dependent of the velocity square of the quadrupole moment oscillations, leading to a forced (by means of radiation emission) and damped linear differential equation governing the evolution of the body oscillations. In planetary tide theory terminology, this equation describes a Kelvin–Voigt damping of the quadrupole moment oscillations, endowed with a deformation inertia term (Correia, Ragazzo & Ruiz 2018). Additionally, the equation of motion for the angular velocity couples to the quadrupole oscillations, as described by equation (16), thus extending beyond previous treatments (Pacini 1968; Gunn & Ostriker 1969; Magalhaes et al. 2012; Hamil et al. 2015). As a consequence of this description, solutions presenting braking indices below the predicted value for a simple radiative model ($n = 3$) were found by means of numerical calculations. In particular, we found that there exist choices for the phenomenological parameters for which the solutions exhibit values similar to those measured in isolated young pulsars (Gourgouliatos & Cumming 2015; Lyne et al. 2015).

A potential direction for further research would involve establishing a connection between the damping effect, which is effectively described by the second term in the dissipation function $D$, as defined in equation (14), and physical processes that could be associated with energy loss driven by quadrupole oscillations. Among the possibilities, thermal and gravitational radiation emissions are the most anticipated causes [see e.g. Chau (1967) for a calculation of the continuous emission of gravitational radiation from pulsars]. Obtaining an actual value for the parameter $\sigma$ for a given star from first principles requires a complete model which, as far as we know, is not yet available.

It should be noticed that the braking index calculated by means of the cumulative method is highly dependent on the initial conditions of the system. Furthermore, if glitches occur during the evolution, as is the case in most of the isolated young pulsars, they can significantly contribute to the value of this index. This aspect can be appreciated, for instance, in the case of the Crab pulsar (Lyne et al. 2015), where the value of $n$ calculated by means of the local method results in $n_L = 2.51$, while the cumulative method leads to $n_C = 2.34$. If the data set is restricted to the period from 1982 onwards, the cumulative method would result in a different value, while the local index would not be significantly affected, as discussed in the previous section. This suggests that the local method provides a more robust index to describe the slow-down of isolated pulsars.
The exact reason behind the occurrence of glitches in an isolated pulsar is still a matter of investigation (Espinoza et al. 2014; Antonopoulou et al. 2022). Most likely, they are associated with redistribution of mass in short time intervals activated by resonance phenomena throughout the evolution of the system. After a glitch, the system approximately returns smoothly to its former state. In the idealized model we investigated here, it is assumed that the shape of the star evolves in time, as governed by the oscillating function $b(t)$. Thus, after the initial transient, the ellipsoid describing the star's surface will oscillate with an almost constant amplitude. In this scenario, localized sub-micro glitches are expected to occur all the time, as suggested by the zigzags in the second and third-order residuals appearing in Fig. 4. A more elaborate model could shed more light on this important issue. Consider a multilayer model for instance. In such a model, resonance effects between the different layer oscillations could lead to significant angular momentum redistribution (Anderson & Itoh 1975; Haskell & Melatos 2015), and possibly to macroscopic glitches, that would then have to be compared to those observed in isolated pulsars. This is an issue that deserves investigation.

A note about the meaning of the $b(t)$ function is in order. First of all, $b(t)$ gives the magnitude of the ellipticity of the star at a given time $t$. It is expected that the equilibrium form of a young pulsar [also denoted as non-spherical relaxed configuration (Gittins, Andersson & Jones 2021)] is reached when its crust solidifies at a given time $t'$ during its evolution. At this instant, according to the model here investigated, the pulsar's ellipticity is determined by the local equilibrium expression $b(t') = \Omega(t')^2/(3\gamma)$. For times $t > t'$, as the angular velocity of the star decreases, e.g. due to the loss of its rotational energy, stresses will arise and the crust shape will evolve from its initial equilibrium state. Numerical solutions for the coupled system discussed in the previous sections show that $b(t)$ is a decreasing oscillatory function of time, whose magnitude deviates from its initial value $b(0) \approx 10^{-4}$ by about $10^{-9}$ per year, while its amplitude of oscillation is of the order of $10^{-18}$. Thus, $10^{-9}$ should be understood as the ellipticity relative to its value at the initial equilibrium configuration. It is interesting to notice that recent publications [see for instance Gittins et al. (2021) and references therein] have argued that a maximum deviation ('mountains') from the equilibrium shape of a neutron star may be something related to a relative ellipticity of about $10^{-8}$ (the modulus of the difference between the ellipticity of the relaxed and strained star). Thus, it is expected that as the pulsar evolves in time the star crust will break when stresses overcome a certain critical value. This is a possible scenario for the occurrence of macro glitches (Baym & Pines 1971).

Among the possible applications of this work, it should be mentioned that the experimental knowledge of the rotation frequency curve of a given pulsar could be used as a starting point to find the best set of physical parameters behind its behaviour. It should be noted however that the damping effects over the oscillations due to emission of thermal radiation and quadrupole gravitational radiation for instance, are not yet fully understood for such systems and also deserve further investigation. Models assuming different forms for the dissipation function and its implications in the possible values of $n$ could be of great value in such investigations.

## ACKNOWLEDGEMENTS

VADL is supported in part by the Brazilian research agency CNPq (Conselho Nacional de Desenvolvimento Científico e Tecnológico) under Grant No. 302492/2022-4. LSR is supported in part by CFIS UC projects (UIDB/04564/2020 and UIDP/04564/2020), and ENGAGE SKA (POCI-01-0145-FEDER-022217), funded by COMPETE 2020 and FCT, Portugal, and also FAPEMIG (Fundação de Amparo à Pesquisa do Estado de Minas Gerais) under Grant No. RED-00133-21.

## DATA AVAILABILITY

The data underlying Section 4 are available in Jodrell Bank Centre for Astrophysics, at http://www.jb.man.ac.uk/pulsar/crab.html (Lyne et al. 1993).

## REFERENCES

* Anderson P. W., Itoh N., 1975, Nature, 256, 25

* Antonopoulou D., Haskell B., Espinoza C. M., 2022, Rep. Prog. Phys., 85, 126901

* Baillieul J., Levi M., 1987, Phys. D: Nonlin. Phenomena, 27, 43

* Baym G., Pines D., 1971, Ann. Phys., 66, 816

* Beskin V. S., Gurevich A. V., Istomin I. N., 1984, Ap&SS, 102, 301

* Blandford R. D., Romani R. W., 1988, MNRAS, 234, 57P

* Caressa P., Bersani A., 2020, Math. Mech. Solids, 26, 785

* Chau W.-Y., 1967, ApJ, 147, 664

* Correia A. C. M., Ragazzo C., Ruiz L. S., 2018, Celest. Mech. Dyn. Astron., 130, 51

* Cutler C., Ushomirsky G., Link B., 2003, ApJ, 588, 975

* Espinoza C. M., Antonopoulou D., Stappers B. W., Watts A., Lyne A. G., 2014, MNRAS, 440, 2755

* Gittins F., Andersson N., Jones D. I., 2021, MNRAS, 500, 5570

* Gold T., 1968, Nature, 218, 731

* Goldreich P., Julian W. H., 1969, ApJ, 157, 869

* Goldwire H. C. J., Michel F. C., 1969, ApJ, 156, L111

* Gourgouliatos K. N., Cumming A., 2015, MNRAS, 446, 1121

* Gullahorn G. E., Isaacman R., Rankin J. M., Payne R. R., 1977, AJ, 82, 309

* Gunn J. E., Ostriker J. P., 1969, Nature, 221, 454

* Hairer E., Nørsett S., Wanner G., 1993, Solving Ordinary Differential Equations: Non stiff problems. Solving Ordinary Differential Equations II: Stiff and Differential-algebraic Problems, Springer-Verlag, Berlin

* Hamil O., Stone J. R., Urbanec M., Urbancová G., 2015, Phys. Rev. D, 91, 063007

* Hamil O., Stone N. J., Stone J. R., 2016, Phys. Rev. D, 94, 063012

* Haskell B., Melatos A., 2015, Int. J. Mod. Phys. D, 24, 1530008

* Hewish A., Bell S. J., Pilkington J. D. H., Scott P. F., Collins R. A., 1968, Nature, 217, 709

* Hobbs G., Lyne A. G., Kramer M., 2010, MNRAS, 402, 1027

* Li L., Tong H., Yan W. M., Yuan J. P., Xu R. X., Wang N., 2014, ApJ, 788, 16

* Lower M. E. et al., 2021, MNRAS, 508, 3251

* Lyne A. G., Pritchard R. S., Graham Smith F., 1993, MNRAS, 265, 1003

* Lyne A. G., Jordan C. A., Graham-Smith F., Espinoza C. M., Stappers B. W., Weltevrede P., 2015, MNRAS, 446, 857

* Magalhaes N. S., Miranda T. A., Frajuca C., 2012, ApJ, 755, 54

* Namkham N., Jaroenjittichai P., Johnston S., 2019, MNRAS, 487, 5854

* Ostriker J. P., Gunn J. E., 1969, ApJ, 157, 1395

* Pacini F., 1967, Nature, 216, 567

* Pacini F., 1968, Nature, 219, 145

* Parthasarathy A. et al., 2019, MNRAS, 489, 3810

* Parthasarathy A. et al., 2020, MNRAS, 494, 2012

* Ragazzo C., Ruiz L. S., 2015, Celest. Mech. Dyn. Astron., 122, 303

* Ragazzo C., Ruiz L. S., 2017, Celest. Mech. Dyn. Astron., 128, 19

* Rencoret J. A., Aguilera-Gómez C., Reisenegger A., 2021, A&A, 654, A47

* Richards D., 1968, IAU Astron. Telegram Circ. No. 2114

* Ruderman M. A., Sutherland P. G., 1975, ApJ, 196, 51
Shapiro S. L., Teukolsky S. A., 1983, Black Holes, White Dwarfs and Neutron
Stars. The Physics of Compact Objects. Wiley-VCH, Weinheim

Shaw B. et al., 2022, MNRAS, 513, 5861

Spitkovsky A., 2006, ApJ, 648, L51

Xu R. X., Qiao G. J., 2001, ApJ, 561, L85

Yue Y. L., Xu R. X., Zhu W. W., 2007, Adv. Space Res., 40, 1491

Zhang C.-M. et al., 2022, Universe, 8, 628

This paper has been typeset from a LaTeX/LATEX file prepared by the author.