Periodic Fast Radio Bursts with Neutron Star Free Precession

Consider an NS with mass M, radius R, dipolar magnetic field of strength B_p, dipole moment $p=\tfrac{1}{2}{B}_{{\rm{p}}}{R}^{3}$ and axis $\hat{{\boldsymbol{p}}}$, spin period P, and spin frequency $\omega =2\pi /P$. The NS could also have complex quadrupole field and internal fields (see below). For simplicity, we assume the NS is homogeneous, with a constant density $\rho =3M/(4\pi {R}^{3})$ and moment of inertia $I=\tfrac{2}{5}{{MR}}^{2}$. We define ${\bar{R}}_{6}=R/({10}^{6}\,\mathrm{cm})$ and ${\bar{M}}_{1.4}=M/(1.4\,{M}_{\odot })$.

In the frame rotating with the NS, the equations of motion describing the evolution of the NS spin vector ${\boldsymbol{\omega }}=\omega \hat{{\boldsymbol{\omega }}}$ is (Zanazzi & Lai 2015)

$$\begin{eqnarray}&&\displaystyle \frac{{\rm{d}}{{\boldsymbol{L}}}_{\mathrm{eff}}}{{\rm{d}}t}+{\boldsymbol{\omega }}\times {{\boldsymbol{L}}}_{\mathrm{eff}}=0,\end{eqnarray} \tag{ 1 }$$

where ${{\boldsymbol{L}}}_{\mathrm{eff}}={{\boldsymbol{I}}}_{\mathrm{eff}}\cdot {\boldsymbol{\omega }}$ is the effective angular momentum of the NS. The effective moment of inertia tensor ${{\boldsymbol{I}}}_{\mathrm{eff}}$ takes account of the nonsphericity of the NS due to rotation and internal magnetic fields, as well as the inertia from the near-zone fields corotating with the NS (Davis & Goldstein 1970; Goldreich 1970; Zanazzi & Lai 2015). Take I_i to be the eigenvalues of ${{\boldsymbol{I}}}_{\mathrm{eff}}$ (effective principal moments of inertia), with ${\hat{{\boldsymbol{I}}}}_{i}$ their associated unit eigenvectors (effective principal axes). For simplicity, we assume I₁ = I₂, but ${\epsilon }_{\mathrm{eff}}=({I}_{3}-{I}_{1})/{I}_{1}\ne 0$ (biaxial NS; we assume $| {\epsilon }_{\mathrm{eff}}| \ll 1$ throughout). Then Equation (1) has the solution (Landau & Lifshitz 1969; Goldreich 1970)

$$\begin{eqnarray}&&\hat{{\boldsymbol{\omega }}}=\sin \theta \cos {\varphi }_{\omega }{\hat{{\boldsymbol{I}}}}_{1}+\sin \theta \sin {\varphi }_{\omega }{\hat{{\boldsymbol{I}}}}_{2}+\cos \theta {\hat{{\boldsymbol{I}}}}_{3},\end{eqnarray} \tag{ 2 }$$

where

$$\begin{eqnarray}&&{\varphi }_{\omega }(t)=({\epsilon }_{\mathrm{eff}}\omega \cos \theta )t+{\varphi }_{\omega 0}\end{eqnarray} \tag{ 3 }$$

is the precession phase of $\hat{{\boldsymbol{\omega }}}$ around ${\hat{{\boldsymbol{I}}}}_{3}$, with ${\varphi }_{\omega 0}={\varphi }_{\omega }(0)$, while θ is the angle between $\hat{{\boldsymbol{\omega }}}$ and ${\hat{{\boldsymbol{I}}}}_{3}$ ($\cos \theta =\hat{{\boldsymbol{\omega }}}\cdot {\hat{{\boldsymbol{I}}}}_{3}$). Notice $\omega$ and θ are constants of motion for Equation (1). The NS precession period ${P}_{\mathrm{prec}}$ is then

$$\begin{eqnarray}&&{P}_{\mathrm{prec}}=\displaystyle \frac{P}{{\epsilon }_{\mathrm{eff}}\cos \theta }.\end{eqnarray} \tag{ 4 }$$

When ${I}_{1}\ne {I}_{2}$, Equation (1) can be solved with qualitatively similar dynamics, except the magnitude of $\omega$ oscillates and θ nutates with time (Landau & Lifshitz 1969; Melatos 1999; Zanazzi & Lai 2015; see also Levin et al. 2020 for why triaxiality does not make the precession rate vary for a nearly spherical rotator).

We postulate that the observed 16.35 day period seen in FRB 180916 is the precession period ${P}_{\mathrm{prec}}$. This constrains ${\epsilon }_{\mathrm{eff}}$ to be

$$\begin{eqnarray}&&{\epsilon }_{\mathrm{eff}}=\displaystyle \frac{7.1\times {10}^{-7}}{\cos \theta }\left(\displaystyle \frac{P}{1\,{\rm{s}}}\right)\left(\displaystyle \frac{16.35\,\mathrm{days}}{{P}_{\mathrm{prec}}}\right).\end{eqnarray} \tag{ 5 }$$

There are several contributions to the nonsphericity parameter ${\epsilon }_{\mathrm{eff}}$. Two primary sources are intrinsic to the NS. The first arises from the internal magnetic field of strength ${B}_{\star }$, leading to a deformation of order

$$\begin{eqnarray}&&{\epsilon }_{\mathrm{mag}}=\beta \displaystyle \frac{{R}^{4}{B}_{\star }^{2}}{{{GM}}^{2}}=1.9\times {10}^{-6}\beta {\left(\displaystyle \frac{{B}_{\star }}{{10}^{15}{\rm{G}}}\right)}^{2}\displaystyle \frac{{\bar{R}}_{6}^{4}}{{\bar{M}}_{1.4}^{2}},\end{eqnarray} \tag{ 6 }$$

where $\beta$ is a constant satisfying $| \beta | \ll 1$ (either $\beta \gt 0$ or $\beta \lt 0$), with a value that depends on the magnetic field's topology (Mastrano et al. 2013); a complex internal field can yield $| \beta | \ll 1$. The second deformation source is an elastic crust that has a rotational bulge with principal axis ${\hat{{\boldsymbol{I}}}}_{3}$ misaligned with $\hat{{\boldsymbol{\omega }}}$, formed when the crust crystallized at a higher rotational frequency (e.g., Goldreich 1970; Cutler et al. 2003). Assuming the NS has a uniform shear modulus μ, the deformation from elasticity is of order (assuming $19\mu \ll 2\rho {gR}$, where $g={GM}/{R}^{2};$ Munk & MacDonald 1975)

$$\begin{eqnarray}\begin{array}{rcl}{\epsilon }_{\mathrm{elast}} & \simeq & \left(\displaystyle \frac{19\mu }{2\rho {gR}}\right)\left(\displaystyle \frac{15{\omega }^{2}}{16\pi G\rho }\right)\\ & = & 2.0\times {10}^{-11}\left(\displaystyle \frac{\mu }{{10}^{30}\,\mathrm{dynes}\,{\mathrm{cm}}^{-2}}\right){\left(\displaystyle \frac{1{\rm{s}}}{P}\right)}^{2}\displaystyle \frac{{\bar{R}}_{6}^{7}}{{\bar{M}}_{1.4}^{3}}.\end{array}\end{eqnarray} \tag{ 7 }$$

Although the frozen-in rotational deformation could be larger when the NS was born rotating fast, the NS's spin-down causes the crust to experience stress and break before the NS has slowed to near its present rotation rate (e.g., Baym & Pines 1971; Cutler et al. 2003). Since the NS may have experienced many crust-quakes over its lifetime, it is reasonable to assume that the frozen-in rotational deformation is of order the present rotational deformation.

In addition to the intrinsic deformations ${\epsilon }_{\mathrm{mag}}$ and ${\epsilon }_{\mathrm{elast}}$, the near-zone fields corotating with the NS induce a precessional torque (Goldreich 1970), and this effect can be incorporated into the effective deformation parameter (Zanazzi & Lai 2015). The dipole field gives (Melatos 1999, 2000; Zanazzi & Lai 2015)

$$\begin{eqnarray}&&{\epsilon }_{{\rm{p}}}=\displaystyle \frac{3{B}_{{\rm{p}}}^{2}{R}^{5}}{20{{Ic}}^{2}}=1.5\times {10}^{-7}{\left(\displaystyle \frac{{B}_{{\rm{p}}}}{{10}^{15}{\rm{G}}}\right)}^{2}\displaystyle \frac{{\bar{R}}_{6}^{3}}{{\bar{M}}_{1.4}}.\end{eqnarray} \tag{ 8 }$$

Equation (8) takes into account the inertia of the field exterior to the NS in vacuum; including the inertia of the field inside the NS (Beskin & Zheltoukhov 2014), or the effect of magnetosphere plasma (Arzamasskiy et al. 2015), modifies Equation (8) by factors of order unity. Similarly, a quadrupolar magnetic field with strengths specified by B_∥ and ${B}_{\delta }$ leads to deformations of order

$$\begin{eqnarray}&&{\epsilon }_{\parallel }=\displaystyle \frac{4}{105}{\left(\displaystyle \frac{{B}_{\parallel }}{{B}_{{\rm{p}}}}\right)}^{2}{\epsilon }_{{\rm{p}}},\hspace{5mm}{\epsilon }_{\delta }=\displaystyle \frac{16}{945}{\left(\displaystyle \frac{{B}_{\delta }}{{B}_{{\rm{p}}}}\right)}^{2}{\epsilon }_{{\rm{p}}};\end{eqnarray} \tag{ 9 }$$

see Zanazzi & Lai (2015) for details and definitions of B_∥ and ${B}_{\delta }$. For magnetic field strengths ($B\sim {10}^{15}\,{\rm{G}}$) and spin periods ($P\sim 1\,{\rm{s}}$) typical of magnetars, we see ${\epsilon }_{\mathrm{mag}}$, ${\epsilon }_{{\rm{p}}}$, ${\epsilon }_{\parallel }$, and ${\epsilon }_{\delta }$ are all feasible ways to effectively deform the NS to give a spin precession period ${P}_{\mathrm{prec}}=16.35\ \mathrm{days}$, but elastic deformation ${\epsilon }_{\mathrm{elast}}$ requires $P\sim 1\,\mathrm{ms}$ to get ${\epsilon }_{\mathrm{elast}}\sim {\epsilon }_{\mathrm{eff}}$. Since this is much shorter than a typical magnetar P value, we will not consider ${\epsilon }_{\mathrm{elast}}$ for the remainder of this work.

In Equation (1), we have neglected the radiative torque, which works to spin down the NS and secularly align $\hat{{\boldsymbol{\omega }}}$ with $\hat{{\boldsymbol{p}}}$. This is valid as long as the shift in the NS precession phase due to spin-down over the course of the observation (${\rm{\Delta }}{\varphi }_{\omega }\sim 2\pi {t}_{\mathrm{obs}}^{2}/[{P}_{\mathrm{prec}}{t}_{\mathrm{sd}}]$) is less than unity, where

$$\begin{eqnarray}&&{t}_{\mathrm{sd}}=\displaystyle \frac{3{c}^{3}I}{2{p}^{2}{\omega }^{2}}=145{\left(\displaystyle \frac{P}{1{\rm{s}}}\right)}^{2}{\left(\displaystyle \frac{{10}^{15}{\rm{G}}}{{B}_{{\rm{p}}}}\right)}^{2}\displaystyle \frac{{\bar{M}}_{1.4}}{{\bar{R}}_{6}^{4}}\,\mathrm{yr}\end{eqnarray} \tag{ 10 }$$

is the spin-down time for the NS, and ${t}_{\mathrm{obs}}=410$ days is the length of time that FRB 180916 was observed. Requiring ${\rm{\Delta }}{\varphi }_{\omega }\lesssim 1$ gives the constraint ${t}_{\mathrm{sd}}\gtrsim 2\pi {t}_{\mathrm{obs}}^{2}/{P}_{\mathrm{prec}}\equiv {t}_{\mathrm{sd},\min }$ for our precession model.⁴

Figure 1 depicts the constraints on the NS spin period P and the strengths of various magnetic field components (internal, dipole, and quadrupole) in order for magnetic deformations (both intrinsic and effective) to produce ${P}_{\mathrm{prec}}=16.35\ \mathrm{days}$. For spin periods $P\sim 0.1-10\ {\rm{s}}$, a range of magnetic field values (B ∼ 10¹⁴ − 10¹⁷) are required, depending on which deformation mechanism dominates ${\epsilon }_{\mathrm{eff}}$. Figure 1 also shows the condition ${t}_{\mathrm{sd}}={t}_{\mathrm{sd},\min }$. For the NS to stably precess over the observed duration of FRB 180916 (${t}_{\mathrm{obs}}$) with a given poloidal field ${B}_{{\rm{p}}}$, the P value must lie somewhat above the red line.

Zoom In Zoom Out Reset image size

Over timescales comparable to ${t}_{\mathrm{sd}}$, the NS spin frequency $\omega$ and precession angle θ evolve due to the radiative torque. When ${t}_{\mathrm{sd}}\gg {P}_{\mathrm{prec}}$, the evolutionary equations for $\omega$ and θ are (Goldreich 1970; Zanazzi & Lai 2015)

$$\begin{eqnarray}&&\displaystyle \frac{{\rm{d}}\omega }{{\rm{d}}t}=-\displaystyle \frac{\omega }{{t}_{\mathrm{sd}}}\left[{\sin }^{2}\psi +{\sin }^{2}\theta \left(1-\displaystyle \frac{3}{2}{\sin }^{2}\psi \right)\right],\end{eqnarray} \tag{ 11 }$$

$$\begin{eqnarray}&&\displaystyle \frac{{\rm{d}}\theta }{{\rm{d}}t}=-\displaystyle \frac{1}{{t}_{\mathrm{sd}}}\cos \theta \sin \theta \left(1-\displaystyle \frac{3}{2}{\sin }^{2}\psi \right).\end{eqnarray} \tag{ 12 }$$

Figure 2 depicts an example of the evolution of P, θ, ${P}_{\mathrm{prec}}$, and χ (the magnetic inclination angle, or angle between $\hat{{\boldsymbol{\omega }}}$ and the dipole axis $\hat{{\boldsymbol{p}}}$) over timescales comparable to ${t}_{\mathrm{sd}}$. For the example given, P and ${P}_{\mathrm{prec}}$ increase, while θ and χ decrease, with time.

Zoom In Zoom Out Reset image size

The central engine of FRBs and the radiation mechanism are uncertain. Given the millisecond timescale of the radio bursts, it is natural that most models associate FRB emissions to magnetized neutron stars (e.g., Lyubarsky 2014; Cordes & Wasserman 2016; Katz 2016; Beloborodov 2017; Lu & Kumar 2018; Margalit et al. 2019). Here we consider a simple geometric model to illustrate how NS precession affects the arrival times of radio bursts from NSs. We use the terminology "arrival times" for FRB emission in our model, but note there are no propagation effects.

Figure 3 presents the setup for our emission model, in the frame corotating with the NS (body frame), with (effective) principal axis ${\hat{{\boldsymbol{I}}}}_{i}$ defining an orthogonal coordinate system. The NS spin axis $\hat{{\boldsymbol{\omega }}}$ is inclined to ${\hat{{\boldsymbol{I}}}}_{3}$ by an angle θ ($\cos \theta =\hat{{\boldsymbol{\omega }}}\cdot {\hat{{\boldsymbol{I}}}}_{3}$), and precesses about ${\hat{{\boldsymbol{I}}}}_{3}$ at the period ${P}_{\mathrm{prec}}$ (Equation (4); see Equation (2)). An observer views the NS in a direction $\hat{{\boldsymbol{n}}}$ constant in the inertial frame, but rotating about $\hat{{\boldsymbol{\omega }}}$ in the body frame with inclination ν ($\cos \nu =\hat{{\boldsymbol{n}}}\cdot \hat{{\boldsymbol{\omega }}}$) and spin period P. Note that in the body frame of the NS, $\hat{{\boldsymbol{n}}}$ satisfies the equation ${\rm{d}}\hat{{\boldsymbol{n}}}/{\rm{d}}t+{\boldsymbol{\omega }}\times \hat{{\boldsymbol{n}}}=0$. Since $\hat{{\boldsymbol{\omega }}}$ evolves over a timescale much longer than $\hat{{\boldsymbol{n}}}$ ($| {\rm{d}}\hat{{\boldsymbol{\omega }}}/{\rm{d}}t| /| {\rm{d}}\hat{{\boldsymbol{n}}}/{\rm{d}}t| \sim {\epsilon }_{\mathrm{eff}}\ll 1$), we can treat $\hat{{\boldsymbol{\omega }}}$ as approximately constant to obtain

$$\begin{eqnarray}\begin{array}{rcl}\hat{{\boldsymbol{n}}}(t) & =\ & \displaystyle \frac{\sin \nu \cos {\varphi }_{n}}{\sin \theta }({\hat{{\boldsymbol{I}}}}_{3}\times \hat{{\boldsymbol{\omega }}})\times \hat{{\boldsymbol{\omega }}}\\ & & -\,\displaystyle \frac{\sin \nu \sin {\varphi }_{n}}{\sin \theta }({\hat{{\boldsymbol{I}}}}_{3}\times \hat{{\boldsymbol{\omega }}})+(\cos \nu )\hat{{\boldsymbol{\omega }}},\end{array}\end{eqnarray} \tag{ 13 }$$

where

$$\begin{eqnarray}&&{\varphi }_{n}(t)=\omega t+{\varphi }_{n0}\end{eqnarray} \tag{ 14 }$$

is the rotation phase of $\hat{{\boldsymbol{n}}}$ around $\hat{{\boldsymbol{\omega }}}$, with ${\varphi }_{n0}={\varphi }_{n}(0)$. The NS's dipole axis $\hat{{\boldsymbol{p}}}$ is fixed in the body frame, inclined to ${\hat{{\boldsymbol{I}}}}_{3}$ by an angle ψ ($\cos \psi =\hat{{\boldsymbol{p}}}\cdot {\hat{{\boldsymbol{I}}}}_{3}$). For concreteness, we take $\hat{{\boldsymbol{p}}}$ to lie in the plane spanned by ${\hat{{\boldsymbol{I}}}}_{1}$ and ${\hat{{\boldsymbol{I}}}}_{3}$. The inclination between $\hat{{\boldsymbol{n}}}$ and $\hat{{\boldsymbol{p}}}$ is specified by the angle ξ ($\cos \xi =\hat{{\boldsymbol{n}}}\cdot \hat{{\boldsymbol{p}}}$). We consider a phenomenological emission model as an example, where we assume the radiation is emitted from a cone centered at $\hat{{\boldsymbol{p}}}$ with opening angle $\alpha$, with the emission intensity I tapering off as $\hat{{\boldsymbol{n}}}$ becomes more misaligned with $\hat{{\boldsymbol{p}}}$:

$$\begin{eqnarray}&&I={I}_{0}\exp \left(-\displaystyle \frac{{\xi }^{2}}{2{\alpha }^{2}}\right).\end{eqnarray} \tag{ 15 }$$

Zoom In Zoom Out Reset image size

Figure 4 shows two examples of the FRB emission pattern produced in our model. Although reproducing the periodicity of FRB 180916 requires ${P}_{\mathrm{prec}}=16.35\ \mathrm{days}$, the NS spin period P and effective ellipticity ${\epsilon }_{\mathrm{eff}}$ remain unconstrained (but related; see Equation (4)). To leave P unconstrained, and to add stochasticity to our simple emission model, we evaluate $\hat{{\boldsymbol{\omega }}}(t)$ (Equation (2)) at N times (N = 400 for the example), which we denote by t_i, spread linearly between ${t}_{1}=0$ to ${t}_{N}={t}_{\mathrm{obs}}$. We then pick ${\varphi }_{n}({t}_{i})={\varphi }_{i}$ (Equation (14)) randomly from a uniform distribution over the interval $[0,2\pi ]$. The observer's orientation ${\hat{{\boldsymbol{n}}}}_{i}=\hat{{\boldsymbol{n}}}({t}_{i})$ is then evaluated with Equation (13), and the FRB emission I at time $t={t}_{i}$ is computed with Equation (15).

Zoom In Zoom Out Reset image size

We see from Figure 4 that, despite the simplicity of our model, it does well in reproducing the spacing of the periodic bursts, as well as the clustering of bursts over the precession phase. The left panel of Figure 4 shows some of the ±2.6 day intervals around multiples of ${P}_{\mathrm{prec}}$ (light cyan bands) have no detectable bursts, while other epochs have multiple detectable bursts. More specifically, when $t/{P}_{\mathrm{prec}}\approx \mathrm{integer}$, some draws at times t_i get no instances of ${\hat{{\boldsymbol{n}}}}_{i}\approx \hat{{\boldsymbol{p}}}$, while other draws get multiple instances of ${\hat{{\boldsymbol{n}}}}_{i}\approx \hat{{\boldsymbol{p}}}$, due to the changing phase of $\hat{{\boldsymbol{n}}}$ around $\hat{{\boldsymbol{\omega }}}$. Notice that no bursts are detected at intervals away from integer multiples of ${P}_{\mathrm{prec}}$. This occurs because over most of the spin precession phase ${\varphi }_{\omega }$, $\hat{{\boldsymbol{\omega }}}$ lies far from $\hat{{\boldsymbol{p}}}$, and $\hat{{\boldsymbol{n}}}$ closely follows $\hat{{\boldsymbol{\omega }}}$. The right panel of Figure 4 shows that the burst intensities I are clustered around the phase 0.5, with a spread that depends on the model parameters (angles ψ, θ, and ν; see Figure 3). Over this spread, the burst intensities vary by two orders of magnitude, with little dependence on ${\varphi }_{\omega }$. All these features were seen in FRB 180916 (The CHIME/FRB Collaboration et al. 2020).

Figure 5 is the same as Figure 4, except instead of drawing a single value of the rotational phase ${\varphi }_{n}$ at $t={t}_{i}$, we pick N_n linearly spaced values spanning the interval $[0,2\pi ]$. With many more points sampled for ${\varphi }_{n}$, we see the FRB emission is confined to the light cyan epochs. The intensity profile shape with the FRB phase and the amount of clustering around phase 0.5 depend on the model parameters (ψ, θ, ν).

Zoom In Zoom Out Reset image size