# Magnetic angle evolution in accreting neutron stars

Anton Biryukov1,2 and Pavel Abolmasov3
1Sternberg Astronomical Institute, Moscow State University, 13 Universitetsky
pr., Moscow, 119234, Russia
2Kazan Federal University, 18 Kremlyovskaya str., Kazan, 420008, Russia
3Department of Physics and Astronomy, FI-20014 University of Turku, Finland
Mailto<EMAIL_ADDRESS>
(Accepted XXX. Received YYY; in original form ZZZ)

###### Abstract

The rotation of a magnetised accreting neutron star (NS) in a binary system is
described by its spin period and two angles: spin inclination $\alpha$ with
respect to the orbital momentum and magnetic angle $\chi$ between the spin and
the magnetic moment. Magnetospheric accretion spins the NS up and adjusts its
rotation axis, decreasing $\alpha$ to nearly perfect alignment. Its effect
upon the magnetic angle is more subtle and relatively unstudied. In this work,
we model the magnetic angle evolution of a rigid spherical accreting NS. We
find that the torque spinning the NS up may affect the magnetic angle while
both $\alpha$ and $\chi$ significantly deviate from zero, and the spin-up
torque varies with the phase of the spin period. As the rotation axis of the
NS is being aligned with the spin-up torque, the magnetic axis becomes
misaligned with the rotation axis. Under favourable conditions, magnetic angle
may increase by $\Delta\chi\sim 15\degr-20\degr$. This orthogonalisation may
be an important factor in the evolution of millisecond pulsars, as it
partially compensates the $\chi$ decrease potentially caused by pulsar
torques. If the direction of the spin-up torque changes randomly with time, as
in wind-fed high-mass X-ray binaries, both the rotation axis of the NS and its
magnetic axis become involved in a non-linear random-walk evolution. The
ultimate attractor of this process is a bimodal distribution in $\chi$ peaking
at $\chi=0$ and $\chi=90\degr$.

###### keywords:

stars: neutron – accretion – X-rays: binaries

††pubyear: 2021††pagerange: Magnetic angle evolution in accreting neutron
stars–A

##  1 Introduction

Magnetic dipole is usually a good approximation for the magnetic field of a
magnetised star, at least in a certain range of radii, much larger than
stellar radius and much smaller than light cylinder. Magnetic dipole is
described by a single vector quantity (magnetic moment), that sets the
direction of the magnetic axis. Magnetic angle $\chi$ is the angle between
spin and magnetic axes of a neutron star (NS, see Figure 1). Though often
ignored, magnetic angle is an important quantity affecting the observational
properties of a NS as well as its rotational evolution. Non-zero value of
magnetic angle is critical for the observations of the pulsating radiation
from any type of radiating NS, accretion- or rotation-powered.

Magnetic angles, their distribution and evolution were comprehensively
investigated for rotation-powered pulsars. Their $\chi$ were found to be
nearly isotropically distributed in the broad interval from almost zero to
nearly 90 degrees (Lyne & Manchester, 1988; Rankin, 1993; Tauris & Manchester,
1998; Malov & Nikitina, 2011). At the same time, their evolution is still a
subject of discussion. Pulsar magnetic angle is thought to be either
decreasing on the spin-down time scale (Philippov et al., 2014), or,
alternatively, increasing toward 90 degrees (Beskin et al., 1993; Novoselov et
al., 2020). In any case, for an isolated pulsar, magnetic angle as a function
of time is completely determined by a single process, the pulsar losses
themselves.

On the other hand, there are several populations of NSs where the rotation of
the star is strongly affected by the effects of binary interaction. These
objects include millisecond pulsars (MSP, Manchester (2017)) and X-ray
binaries. The latter are normally classified according to the mass of the
donor star as low- and high-mass X-ray binaries (LMXBs and HMXBs
correspondingly, Bhattacharya & van den Heuvel 1991). HMXB containing
magnetised NSs are usually observed as X-ray pulsars (Caballero & Wilms,
2012). In turn, most LXMBs have weakly magnetised accretors that does not
allow to study their magnetic angles observationally (their magnetic fields
are either buried by accreting matter or decay due to Joule losses in the
crust, see later Section 2.3). However, accreting millisecond pulsars (Patruno
& Watts, 2021), apparently having relatively large magnetospheres, also fall
into this category. As the amount of observational data on NSs in binary
systems increases, there is a growing demand for a theory describing the
physics of NS rotation affected by mass accretion and additional external
torques.

Unfortunately, magnetic angle distributions in MSP and X-ray binaries are
known much worse than for classical radiopulsars. Moreover, they are also
affected by multiple processes, most importantly by the high-angular-momentum
accretion from the binary companion.

For MSP, Chen et al. (1998) find an intriguing bi-modal distribution in
magnetic angle. Considering a sample of 11 objects, they found that some
pulsars are close to magnetic alignment, while the others are close to the
orthogonal rotation. If the torques acting on the NS do not distinguish
between the two magnetic poles and depend smoothly on magnetic angle, these
are the two expected equilibrium points. Using a different line of reasoning,
Lipunov & Shakura (1980) come to the same conclusion. However, it is hard to
explain why a significant number of objects should evolve towards orthogonal
rotation.

More recently, Johnson et al. (2014) considered a larger sample of objects
using different data (including gamma-ray data from Fermi/LAT for 40 objects)
and techniques. Apparently, the pulsars considered show a broad distribution
in magnetic angle with no significant bi-modality, though a large number of
objects still show very large (above 70°) magnetic angles. This agrees with
the estimates made by Bogdanov et al. (2008) that argue for large magnetic
angles ($\chi\gtrsim 40^{\circ}$) in some MSP.

MSPs are commonly understood as “recycled”, i. e. spun-up by accretion after
the end of their initial period of pulsar activity (Tauris & van den Heuvel,
2006). This interpretation is supported by the existence of a population of
accreting MSPs (Hessels, 2008; Patruno & Watts, 2021). There is also at least
one MSP, PSR J1023+0038, showing pulsations both in X-ray and in radio,
switching between radio-bright/X-ray-faint and radio-faint/X-ray-bright states
(Archibald et al., 2009). It is tempting to interpret these two states as
accretor and ejector regimes, but the recent study by Veledina et al. (2019)
shows that most likely both states in this object are rotationally powered but
require the presence of an accretion disc.

As for X-ray pulsars, the constraints upon their magnetic field geometry are
still highly uncertain and model-dependent. Leahy (1990) estimated the
geometry of the radiating areas of 15 X-ray pulsars using a simple non-
relativistic hot spot model. Magnetic angles found in this study cluster
around $\chi\sim 24\degr$. Later, Annala & Poutanen (2010) considering a
larger sample of X-ray pulsars with relativistic models, came to a mean value
of $\chi_{\rm mean}=40\pm 4\degr$.

In this paper, we present an attempt to trace the NS’s magnetic angle
evolution during magnetospheric accretion. In Section 2, we derive the
equations of rotational evolution including different external torques and
build a phenomenological model. We show that the crucial factor for the
evolution of magnetic angle is the modulation of the spin-up torque with the
spin phase. In Section 3, we calculate the rotational evolution for four
different cases aimed to represent an MSP progenitor or an X-ray pulsar. We
discuss the results and consider certain extensions and implications of the
model in Section 4. Conclusions are made in Section 5.

##  2 Rotational evolution of an accreting neutron star

###  2.1 General model

Let us consider the rotation of a solid spherical NS with variable moment of
inertia $I(t)$. We will use a coordinate frame $xyz$ associated with the star
as follows: axis $z$ is aligned with the magnetic moment of the star, $x$ and
$y$ are perpendicular to $z$ and form a right-handed set of vectors. The
orientation of $yz$ plane is defined by the spin $\mn@boldsymbol{\it\Omega}$
of the NS: we assume that the spin axis lies within the plane (see Figure 1).
The considered reference frame $xyz$ co-rotates with the star.

Figure 1: Main axes and angles of an accreting neutron star. Unit vectors
along the magnetic, disc, and spin axes are $\mn@boldsymbol{\it m}$,
$\mn@boldsymbol{\it d}$, and $\mn@boldsymbol{\it s}$, respectively. The angles
shown are $\chi$ (magnetic angle), $\alpha$ (rotation axis inclination with
respect to the axis of the disc), and $\theta$ (instantaneous angle between
$\mn@boldsymbol{\it d}$ and $\mn@boldsymbol{\it m}$).

In this frame, the rotational evolution of the NS is described in vector form
by Euler’s equations (Landau et al., 1976)

$\frac{\mathrm{d}\mn@boldsymbol{\it
L}(t)}{\mathrm{d}t}+\mn@boldsymbol{\it\Omega}(t)\mn@boldsymbol{\it\times}\mn@boldsymbol{\it
L}(t)=\mn@boldsymbol{\it N}(t),$ (1)

where $\mn@boldsymbol{\it\Omega}$ is the instantaneous angular velocity of the
star,

$\mn@boldsymbol{\it L}(t)=I(t)\mn@boldsymbol{\it\Omega}(t)$ (2)

is its full angular momentum, and $\mn@boldsymbol{\it N}$ is the external
torque.

We ignore the possible oblateness of the star for the following reasons.

Expected values of magnetic and other non-hydrostatic ellipticities of NSs are
very low $\epsilon\sim 10^{-8}-10^{-10}$. Their effect on the rotation of the
star is restricted to very slow precession on the time scales or $\sim
10^{5}-10^{7}$ years (Goldreich, 1970; Melatos, 2000; Haskell et al., 2008),
that is much longer than the spin-disc alignment time scale (see Equation 16
below).

Hydrostatic deformation related to rotation may reach few per cent if the spin
period is of the order several milliseconds, but the deformation always adds
an angular momentum component aligned with $\mn@boldsymbol{\it\Omega}$ and
hence may be considered as a correction to the variable scalar moment of
inertia. In Appendix A below we consider this issue in more details.

Substituting Equation (2) into Equation (1), we get an equation for the
rotational evolution of a spherical NS in the most general form

$I\dfrac{\mathrm{d}\mn@boldsymbol{\it\Omega}}{\mathrm{d}\it
t}=\mn@boldsymbol{\it N}-\dfrac{\mathrm{d}\it I}{\mathrm{d}\it
<EMAIL_ADDRESS>(3)

The full torque $\mn@boldsymbol{\it N}$ acting on an accreting star consists
of three components

$\mn@boldsymbol{\it N}=\mn@boldsymbol{\it N}_{\mathrm{acc}}+\mn@boldsymbol{\it
N}_{\mathrm{mag}}+\mn@boldsymbol{\it N}_{\mathrm{psr}},$ (4)

where $\mn@boldsymbol{\it N}_{\mathrm{acc}}$ describes the NS spin-up due to
accretion, $\mn@boldsymbol{\it N}_{\mathrm{mag}}$ is magnetic braking caused
by disc-magnetosphere interaction, and $\mn@boldsymbol{\it N}_{\mathrm{psr}}$
corresponds to pulsar losses. In principle, interaction between the disc and
the magnetic field also creates additional precession and warping torques (see
Lipunov & Shakura 1980; Wang 1981a; Lai et al. 2011; Lai 2014). It is easy to
show, however, that if the interaction region between the magnetosphere and
the disc is small, the warping torque is also small and thus we omit it in our
basic simulations. Nevertheless, we will still consider the effect of the
warping torque below in Section 4.5.

In the co-rotating $xyz$ frame we use, $\mn@boldsymbol{\it N}_{\mathrm{acc}}$
rotates around $\mn@boldsymbol{\it\Omega}$ with NS spin frequency $-\Omega$,
while the other two components of the torque change relatively slowly
retaining their direction in the frame of the star on a spin time scale.

At the same time, we expect that, if the angle $\theta$ between the NS
magnetic moment and accretion disc axis changes with spin phase, both the
instantaneous mass accretion rate and, therefore, spin-up torque acting on the
NS are modulated. Thus, we will consider the spin-up torque
$\mn@boldsymbol{\it N}_{\mathrm{acc}}$ in the form

$\mn@boldsymbol{\it N}_{\mathrm{acc}}=N_{0}f(\cos\theta)\mn@boldsymbol{\it
d},$ (5)

where

$N_{0}=\dot{M}\sqrt{GM_{*}r_{\mathrm{m}}},$ (6)

and $\mn@boldsymbol{\it d}$ is the unit vector in the direction of the disc
axis. Here, $\dot{M}$ is mass accretion rate, $M_{*}$ is the mass of the NS,
and $r_{\rm m}$ is the radius of the magnetosphere. Vector $\mn@boldsymbol{\it
d}$ rotates around the spin axis $\mn@boldsymbol{\it\Omega}$ in the direction
opposite to the rotation of the star so that

$\dot{\mn@boldsymbol{\it
d}}=-\mn@boldsymbol{\it\Omega}\times\mn@boldsymbol{\it d}.$ (7)

Geometrically

$\mn@boldsymbol{\it d}=\left(\begin{array}[]{c}\sin\alpha\sin\psi\\\
\sin\chi\cos\alpha-\cos\chi\sin\alpha\cos\psi\\\
\cos\chi\cos\alpha+\sin\chi\sin\alpha\cos\psi\\\ \end{array}\right),$ (8)

where $\alpha$ is the angle between $\mn@boldsymbol{\it d}$ and
$\mn@boldsymbol{\it\Omega}$, or spin inclination with respect to the disc.
Here, $\psi=-\int\Omega\mathrm{d}t$ is the rotational phase equal, up to the
constant of integration, to the dihedral angle between the
$\mn@boldsymbol{\it\Omega}-\mn@boldsymbol{\it m}$ and
$\mn@boldsymbol{\it\Omega}-\mn@boldsymbol{\it d}$ planes, and
$\mn@boldsymbol{\it m}=\mn@boldsymbol{\it\mu}/\mu$ is the $z$-axis basis unit
vector directed along the magnetic moment $\mn@boldsymbol{\it\mu}$, so that

$\mn@boldsymbol{\it m}=\left(\begin{array}[]{c}0\\\ 0\\\ 1\\\
\end{array}\right).$ (9)

We assume that at the time $t=0$ all the three vectors ($\mn@boldsymbol{\it
d}$, $\mn@boldsymbol{\it m}$, and $\mn@boldsymbol{\it\Omega}$) lie in the same
plane $yz$.

Radius $r_{\rm m}$ in Equation (6) is the size of the magnetosphere calculated
according to the classical formulae (see for instance section 6.3 of Frank et
al. 2002)

$r_{\mathrm{m}}=\xi\left(\dfrac{\mu^{4}}{2GM_{*}\dot{M}^{2}}\right)^{1/7},$
(10)

where $\xi\approx 0.5$ is a dimensionless constant dependent on the unknown
details of the disc-magnetosphere interaction. Dimensionless function
$f(\cos\theta)$ in Equation (5) describes the modulation of the accretion
torque as a function of the angle $\theta$ between the magnetic axis and disc
axis $\mn@boldsymbol{\it d}$. A phenomenological representation for this
function will be introduced below in Section 2.2.

We take the accretion term, Equation (5), into account only when the inner
edge of the accretion disc is inside the co-rotation radius
$r_{\mathrm{co}}=(GM_{*}/\Omega^{2})^{1/3}$. If this condition is violated,
the NS is likely to enter the so-called propeller regime when accretion is
suppressed by centrifugal forces (Illarionov & Sunyaev, 1975).

As the processes responsible for the spin-down of accreting magnetised NSs are
poorly constrained and probably related to the open magnetic field lines co-
rotating with the star (Lovelace et al., 1995), we adopt the following
expression for the magnetospheric spin-down torque

$\mn@boldsymbol{\it
N}_{\mathrm{mag}}=-\dfrac{\mu^{2}}{3r^{3}_{\mathrm{co}}}\mn@boldsymbol{\it
s},$ (11)

where $\mn@boldsymbol{\it s}=\mn@boldsymbol{\it\Omega}/\Omega$ is a unit
vector directed along the spin axis:

$\mn@boldsymbol{\it s}=\left(\begin{array}[]{c}0\\\ \sin\chi\\\ \cos\chi\\\
\end{array}\right).$ (12)

Though Equation (11) was initially obtained in the framework of disc-
magnetosphere interaction model (Ghosh et al., 1977), its general form also
appears in the models involving magnetospheric outflows (Illarionov &
Kompaneets, 1990; Rappaport et al., 2004). The importance of magnetospheric
outflows as a spin-down mechanism for accreting NS was confirmed by MHD
simulations (Romanova et al., 2009; Zanni & Ferreira, 2013; Parfrey et al.,
2017). Note the difference between the two spin-down mechanisms: pulsar torque
(considered below), possibly enhanced by the increase of the number of open
magnetic field lines (Parfrey et al., 2016), and magnetospheric ejections that
we associate with $\mn@boldsymbol{\it N}_{\rm mag}$.

The direction of this spin-down torque has not been studied for inclined
rotators, and it is difficult to predict from the basic principles. As the
outflows are formed by centrifugal forces, we will assume that
$\mn@boldsymbol{\it N}_{\mathrm{mag}}$ is always directed along the spin axis
of the neutron star. This torque depends on the size of the NS magnetopshere
and can be included into the calculations only if
$r_{\mathrm{m}}<r_{\mathrm{LC}}=c/\Omega$, where $r_{\mathrm{LC}}$ is the
light cylinder radius. Otherwise, the magnetosphere is in the ejector regime
(Shvartsman, 1970; Illarionov & Sunyaev, 1975) and is not affected by the
possible presence of the accretion disc.

Finally, for the pulsar spin-down term, we use the results obtained by
Spitkovsky (2006) and Philippov et al. (2014)

$\displaystyle{\mn@boldsymbol{\it
N}}_{\mathrm{psr}}=K_{\mathrm{psr}}\left\\{\left(1+\sin^{2}\chi\right)\mn@boldsymbol{\it
s}+((\mn@boldsymbol{\it m\times s})\mn@boldsymbol{\it\times
s})\cos\chi\right\\},$ (13)

where

$K_{\mathrm{psr}}=-\dfrac{\mu^{2}}{r_{\mathrm{LC}}^{3}}$ (14)

is the characteristic pulsar braking torque. The first term in Equation (13)
describes the torque spinning the NS down while the second one is related to
the evolutionary secular decrease of obliquity. Another possible component of
the pulsar loss is the so-called “anomalous torque” (e.g. Beskin &
Zheltoukhov, 2014). This torque is proportional to
$\mn@boldsymbol{\it\Omega}\times\mn@boldsymbol{\it m}$ and causes a forced
precession of the NS spin axis around its magnetic moment. However, as we will
show in Section 2.2, the contribution of this precession term to all the
evolutionary equations is exactly zero. Therefore, we ignore the anomalous
torque in our calculations.

The more important issue, however, is the validity of Equation (13) during
accretion. It seems clear that accreting plasma changes the structure of the
magnetosphere near the polar caps and affects the distribution of
magnetospheric currents and, therefore, the character of the torques acting on
the star. However, pulsar spin-down time scale,

$\tau_{\rm psr}\sim\frac{c^{3}I}{\mu^{2}\Omega^{2}}\simeq 10^{8}I_{45}P_{\rm
s}^{2}\mu^{-2}_{30}\hskip 5.69046pt{\rm\,yr},$ (15)

is much longer than that of the spin-up from the accretion

$\tau_{\mathrm{acc}}\sim I\Omega/N_{0}\simeq
600I_{45}\dot{M}_{1}^{-6/7}\left(\dfrac{M_{*}}{\rm
M_{\sun}}\right)^{-9/7}\mu^{-2/7}_{30}{\rm\,yr}.$ (16)

Here $I_{45}=I/10^{45}{\rm g\,cm^{2}}$, $P_{\rm s}=P/1$s,
$\mu_{30}=\mu/10^{30}{\rm\,G\,cm^{3}}$, and mass accretion rate
$\dot{M}_{1}=\dot{M}/\dot{M}_{\rm Edd,{\sun}}$ is normalised by the Eddington
accretion rate for one Solar mass, namely $\dot{M}_{\mathrm{Edd,{\sun}}}=4\pi
G\rm M_{\sun}/\varkappa_{\rm T}c\approx 1.6\times
10^{17}{\rm\,g\,s^{-1}}\approx 2.5\times 10^{-9}\ \rm M_{\sun}\,yr^{-1}$. Even
for MSP $\tau_{\mathrm{psr}}\gtrsim 10^{8}{\rm\,yr}$ (Manchester, 2017)
because of smaller magnetic fields $\mu_{30}\sim 10^{-4}-10^{-2}$.

Hence, pulsar torque in the form (13) somehow affects the NS obliquity
evolution only during the ejector stage. So, we always keep the pulsar losses
on during our simulation for the sake of simplicity.

It is also possible that the presence of an accretion disc increases the
pulsar spin-down torques (Parfrey et al., 2016). We discuss this possibility
in Section 4.1.

###  2.2 Averaged torque

In the reference frame we use, accretion torque (5) rotates rapidly around the
spin axis of the NS. Its projection upon the spin axis may vary with the
rotational phase $\psi$, that is of primary importance for the evolution of
obliquity $\chi$. As the angle $\theta$ sets the relative orientation of the
disc and the star magnetosphere we propose that accretion torque depends on it
as $N_{\mathrm{acc}}\propto f(\cos\theta)$, where, according to the spherical
cosine rule,

$\cos\theta=\cos\alpha\cos\chi+\sin\alpha\sin\chi\cos\psi.$ (17)

The function $f(\cos\theta)$ is unknown and depends on the details of disc-
magnetosphere interaction. As the north and south magnetic poles are
indistinguishable from the point of view of magnetohydrodynamics, the unknown
function $f(\cos\theta)$ should be an even function of its argument.

Torque modulations may be related to the changes in the configuration of the
magnetosphere and its interface with the disc on dynamical time scales. The
mean value of $N_{\rm acc}$ should not be affected by these variations, as the
mean mass accretion rate is set by the boundary conditions (mass transfer rate
from the donor star), and the net angular momentum depends, apart from
$\dot{M}$, on the size of the magnetosphere, that is itself a weak function of
inclination. Though in general accreting NSs should be inclined, and accretion
upon inclined dipoles was studied both analytically (Wang, 1981a; Bozzo et
al., 2018) and numerically (Romanova et al., 2009), the problem of mass
accretion upon a misaligned oblique rotating dipole (where both $\alpha$ and
$\chi$ are non-zero) is, to our knowledge, relatively unstudied. Notable
exceptions are Lai et al. (2011) who considered the aligning torques
analytically, and Romanova et al. (2020), where magnetospheric accretion upon
a misaligned oblique magnetic dipole is considered by means of 3D
magnetohydrodynamic simulations. One of the interesting findings of Romanova
et al. (2020) is that the mass accretion rate upon the simulated star varies
with time at about spin period. As we show below, the only case when magnetic
angle changes due to accretion torque is this misaligned oblique case
($\alpha\neq 0$ and $\chi\neq 0$), and the accretion torque variations within
the spin period plays crucial role in the evolution of $\chi$.

Figure 2: Accretion torque modulation factor $f(\eta,\psi)$ as a function of
spin phase $\psi$ for various $\eta$. It was calculated by combining Equations
(18) and (17). Constant $\eta$ describing the modulation amplitude was set to
0.1, 0.3, 0.5, 0.75 and 0.99. Magnetic angle $\chi$ and spin-disc angle
$\alpha$ are equal to 60∘ and 45∘, respectively.

For the modulation function $f$, we propose the simplest possible non-trivial
even-order polynomial form

$f(\cos\theta)=A(\eta,\alpha,\chi)(1-\eta\cos^{2}\theta),$ (18)

where $0<\eta<1$111We presume that only positive values of $\eta$ less then 1
are physically reasonable. This corresponds to the case when the instantaneous
accretion torque is stronger for larger misalignment (smaller $\cos\theta$).
Indeed, it seems much easier for the matter to permeate the magnetosphere if
the magnetic axis of the star lies closer to the disc plane. Moreover,
$\eta<1$ means that, independently of the spin phase, accretion torque can not
be directed opposite to the angular momentum in the disc. is a constant
describing the accretion torque modulation within the spin period, and
$A(\eta,\alpha,\chi)$ is a normalization factor. The latter can be calculated
by averaging $f(\cos\theta)$ over the NS spin period assuming $\alpha$ and
$\chi$ vary negligibly within one revolution of the star,

$A(\eta,\alpha,\chi)=\left(\dfrac{1}{2\uppi}\int_{0}^{2\uppi}(1-\eta\cos^{2}\theta)d\psi\right)^{-1}.$
(19)

This multiplier ensures that the mean accretion torque value remains
unchanged. Taking the integral in (19), we get an analytic expression for $A$:

$A(\eta,\alpha,\chi)=\left[1-\dfrac{\eta}{2}(\sin^{2}\chi\sin^{2}\alpha+2\cos^{2}\chi\cos^{2}\alpha)\right]^{-1}.$
(20)

The modulation factor (18) is shown as a function of spin phase $\psi$ in
Figure 2. For all the calculations below we set $\eta=0.99$ in order to
estimate the maximal effect from the proposed modulation.

Equation describing the evolution of the rotation velocity of the star may be
then obtained by projecting (3) upon the instantaneous direction of
$\mn@boldsymbol{\it s}$ with subsequent averaging over the spin period

$I\langle\dot{\Omega}\rangle=\langle\mn@boldsymbol{\it N\cdot
s}\rangle-\dot{I}\Omega,$ (21)

Hereafter $\langle x\rangle$ means

$\langle x\rangle=\dfrac{1}{2\pi}\int_{0}^{2\pi}xd\psi.$ (22)

As a result, we get, for a spherical star,

$I\left\langle\dfrac{\mathrm{d}\Omega}{\mathrm{d}t}\right\rangle=N_{0}\cos\alpha+N_{\mathrm{mag}}+K_{\mathrm{psr}}(1+\sin^{2}\chi)-\dfrac{\mathrm{d}\it
I}{\mathrm{d}\it t}\Omega.$ (23)

Similar equation for the evolution of the star-disc inclination angle $\alpha$
can be obtained in the same way by projecting (3) onto the unit vector
$\mn@boldsymbol{\it e}$ which always lies within the $\mn@boldsymbol{\it
s}-\mn@boldsymbol{\it d}$ plane, orthogonal to the angular velocity vector
$\mn@boldsymbol{\it e}\perp\mn@boldsymbol{\it s}$, and makes an acute angle
with $\mn@boldsymbol{\it d}$:

$\mn@boldsymbol{\it e}=\left(\begin{array}[]{c}\sin\psi\\\
-\cos\chi\cos\psi\\\ \sin\chi\cos\psi\\\ \end{array}\right).$ (24)

As $\mn@boldsymbol{\it\dot{\Omega}\cdot e}=-\Omega\dot{\alpha}$, projection
along $\mn@boldsymbol{\it e}$ results in

$I\Omega\left\langle\dfrac{d\alpha}{dt}\right\rangle=-\langle\mn@boldsymbol{\it
N\cdot e}\rangle=-N_{0}\sin\alpha.$ (25)

At last, evolution of the magnetic angle is calculated similarly by projecting
(3) onto $\mn@boldsymbol{\it m}$. After the averaging and using the
relationship $d\Omega_{z}/dt=\dot{\Omega}\cos\chi-\Omega\sin\chi\dot{\chi}$,
we finally obtain

$\begin{array}[]{l}I\Omega\left\langle\dfrac{d\chi}{dt}\right\rangle=\eta
A(\eta,\alpha,\chi)N_{0}\sin^{2}\alpha\cos\alpha\sin\chi\cos\chi\\\
\qquad\qquad\qquad+K_{\mathrm{psr}}\sin\chi\cos\chi\\\ \end{array}$ (26)

This equation is the main theoretical result of the paper. The evolution of
the magnetic angle, unlike that of $\Omega$ and $\alpha$, is affected by the
variations of the spin-up torque within the spin period of the NS.

Note that the pulsar anomalous torque contributions to all the three equations
vanish after the period-averaging procedure. For the equations for $\Omega$
and $\chi$, this is trivial, as pulsar torque is always orthogonal to
$\mn@boldsymbol{\it s}$ and $\mn@boldsymbol{\it m}$. Multiplying the anomalous
torque by $\mn@boldsymbol{\it e}$ produces a triple product
$\mn@boldsymbol{\it e}\cdot(\mn@boldsymbol{\it s}\times\mn@boldsymbol{\it
m})=\sin\chi\sin\psi$ that becomes zero after averaging over $\psi$. Therefore
this torque does not affect the evolution of $\alpha$ as well.

Equations (23), (25), and (26), together with (5), (11), and (13) completely
describe the rotational evolution of a spherical accreting neutron star.

###  2.3 Evolution of the physical parameters

It has been shown that millisecond pulsars are on average $\sim 0.2\rm
M_{\sun}$ heavier than normal isolated neutron stars, that likely reflects the
amount of matter accreted by the NS (Zhang et al., 2011; Cheng et al., 2014).
This additional mass affects the inertia of the star.

Here, we will use an approximate phenomenological relationship between the
moment of inertia of the NS and its mass, valid for a broad range of equations
of state and masses (Biryukov et al., 2017a)

$I\simeq 10^{45}\ \left(\dfrac{M_{*}}{\rm M_{\sun}}\right){\rm\,g\,cm^{2}}.$
(27)

The last term in (23) contributes to total spin-down as

$-\dfrac{1}{I}\dfrac{\mathrm{d}\it I}{\mathrm{d}\it t}\Omega\approx-3\times
10^{-14}\mbox{ }\dot{M}_{1}M_{*,1.5}^{-1}P^{-1}_{0.01}\ \mathrm{rad\,s^{-2}},$
(28)

where the mass of the star $M_{*}$ is normalised by $1.5\rm M_{\sun}$, and its
spin period as $P_{0.01}=P/0.01$ s.

Another parameter affected by accretion is magnetic moment. We assume that the
accreted matter changes the magnetic moment of the star as

$\mu=\mu_{0}\left(1+\dfrac{\Delta M_{*}}{\cal M}\right)^{-14/11}$ (29)

where $\mu_{0}$ is the initial magnetic moment, $\Delta M_{*}$ is the amount
of the matter accreted by the star and

${\cal M}=1.1\times
10^{-5}\dot{M}_{1}^{1/7}\mu_{0,30}^{3/14}R_{\mathrm{12.5}}^{3}\,\rm M_{\sun}.$
(30)

Here, initial magnetic moment $\mu_{0,30}$ is expressed in the units of
$10^{30}$ G cm3, and the radius is normalised as
$R_{\mathrm{12.5}}=R_{*}/12.5$ km. The law (29) approximately reproduces the
magnetic field burial scaling predicted by Melatos & Phinney (2001) for
Eddington mass accretion rate and pulsar-scale magnetic field (the mass scale
depends on mass accretion rate, as well as on magnetic field, very weakly).

At the same time, one may expect stellar magnetic field decay due to
dissipative processes within it crust (Reisenegger, 2003). However, we have
not taken this effect into account, as the field burial time scale

$\tau_{\mathrm{bur}}=-\dfrac{\mu}{\dot{\mu}}\approx 3.4\times
10^{3}\dot{M}_{1}^{-6/7}\mu_{0,30}^{3/14}R_{\mathrm{12.5}}^{3}\ {\rm\,yr}$
(31)

is significantly shorter than the decay time scale $\tau_{\mathrm{dec}}\gtrsim
10^{6}$yr (Gullón et al., 2014; Igoshev & Popov, 2015; Biryukov et al.,
2017b).

There are other mechanisms possibly important for the rotational evolution of
NSs, that we do not consider in this paper. As we already mentioned in Section
2.1, non-sphericity of the star is probably a minor issue, though deformation
in the direction non-coplanar with $\mathbf{s}$ and $\mathbf{m}$ may lead to
free precession and a more complicated evolutionary path. Presence of higher
multipoles, especially misaligned with the dipole, is another important factor
that requires a separate study. The importance of higher multipoles was
stressed, for example, by Pétri (2019). Probably, the most relevant
contribution of higher multipoles would be the deformation of the NS not
aligned with the axis of the magnetic dipole. In our calculations, we consider
NS as a rigid body, that may be incorrect if the core of the star is weakly
coupled to the crust (Casini & Montemayor, 1998). If the crust itself is non-
rigid (as in the model by Ruderman 1991) or the magnetic field lines are
allowed to slip through the crust, the rotational evolution of the NS should
also become much more complicated. All these effects, however, are not usually
included in the calculations of NS evolution due to the large uncertainties in
their parameters.

##  3 Results

Table 1: Parameters of the NS evolution models considered in the paper. Initial spin period $P_{0}=2\pi/\Omega_{0}$. Modulation constant $\eta=0.99$ are the same for all models. Initial magnetic moment $\mu_{0,30}$ is normalised on $10^{30}$ G cm3, while accretion rate $\dot{M}_{1}$ is taken in Eddington rates $\dot{M}_{\rm Edd,{\sun}}$ for one Solar mass. Model |  $P_{0}$, sec | $\mu_{0,30}$ | $\dot{M}_{1}$
---|---|---|---
A | 0.1 | 5.0 | 0.01
B | 0.1 | 5.0 | 1.0
C | 5.0 | 1.0 | 1.0
D | 100.0 | 1.0 | 1.0

Bringing together all the physics described in the previous section allows to
reconstruct the rotational evolution of an accreting NS, including magnetic
angle (Equation 26). In our model, there are two processes affecting $\chi$:
the modulation term and pulsar losses. While the last is well known in the
context of the evolution of radiopulsars, it is usually not important during
the accretion stage. The modulation term, however, is proportional to
$\sin^{2}\alpha$ and thus is only important during the rotation axis alignment
stage. The resulting change in $\chi$ may be under favourable circumstances
comparable to the initial misalignment $\alpha$ and important for the
resulting distribution in magnetic angles.

In this Section, we consider four different scenarios of NS evolution in
binary systems undergoing mass transfer at a fixed rate. The models differ in
the initial conditions and in the mass transfer rates. Their parameters are
given in Table 1. The initial mass of the NS in all the models is $1.5\rm
M_{\sun}$. Modulation factor $\eta$ for the accretion torque was everywhere
fixed to 0.99. We use the radius of the NS equal to 12.5km, but this does not
affect the results much as the radius of the magnetosphere is always larger.
The moment of inertia is given by Equation (27), and the magnetic moment of
the star decreases with time according to the burial law (29).

We track the evolution for $10^{7}-10^{8}$ yr depending on the model until the
mass of the NS changes considerably (up to $\sim 0.1\rm M_{\sun}$). The four
models differ in the magnetic moment, mass accretion rate, and initial NS
spin, and may be loosely associated with the different scenarios of NS
rotational evolution. Lower mass accretion rates are expected in binary
systems with lower-mass donor stars (Patruno & Watts, 2021).

Note that we do not pretend here on capturing all the details of the binary
evolution. In combination with binary evolution codes, out results will
hopefully allow to reproduce a realistic picture of the evolving magnetic
angle distributions in real X-ray binary populations.

First model (A) corresponds to the case of a putative MSP progenitor with a
low-mass donor. Indeed, accreting millisecond pulsars usually have very low
mass donor stars and relatively low mass accretion rates $\sim 0.01-0.1$
$\dot{M}_{\rm Edd,\,{\sun}}$ (Patruno & Watts, 2021; Tauris et al., 2012). At
the same time, accretion at a mean rate smaller than $\sim 0.01M_{\rm
Edd,\,{\sun}}$ is insufficient to spin up an average MSP progenitor to a
typical observed period of $\sim 10$ms (Pan et al., 2013).

In this model initial parameters of the NS are similar to those of a newly
born radiopulsar: spin period $\sim 0.1$s and initial surface magnetic field
$5\times 10^{12}$G. The adopted mass accretion rate is
$0.01\dot{M}_{\mathrm{Edd},{\sun}}$, that allows to track the rotational
evolution for more that tens of millions of years.

The second model (B) describes a similar situation but with a larger mass
accretion rate of about Eddington. This is closer to the evolution of an X-ray
pulsar, though real X-ray pulsars also show strong modulations (orbital and
stochastic) of mass accretion rate we did not include in the model. We
consider chosen $\dot{M}$ as an average accretion rate within the system
evolutionary timescale.

The remaining two models (C and D) explore the case of an initially slowly
rotating NS spun up at the Eddington rate. Thus, model D may represent the
case of the rapid spin-up observed for NGC 300 ULX1 (Carpano et al., 2018;
Vasilopoulos et al., 2018).

For each of the models, we have considered different initial values of
$\alpha$ and $\chi$, arranged in a two-dimensional grid with a step of one
degree along each coordinate.

Figure 3: Rotational evolution of accreting NSs, according to the models A-D
listed in Table 1. Solid lines represent the cases with initial
$\alpha_{0}=45^{\circ}$, dashed are for $\alpha_{0}=135^{\circ}$. Thin dotted
lines mark the initial values of the quantities. Dipole magnetic moment of the
star $\mu_{30}$ is normalised to $10^{30}{\rm\,G\,cm}^{3}$.

In Figure 3, selected evolutionary trajectories are shown. For each model, we
show the cases with initial magnetic angle $\chi_{0}=60\degr$ and two
inclination angles, $\alpha_{0}=45\degr$ and $\alpha_{0}=135\degr$. For all
the models, a short-term ($\sim 10^{4}$ years for model A, and about $10^{3}$
yr for all the others) magnetic angle counter-alignment episode is clearly
seen at the beginning of the accretion stage.

The evolution of $\chi$ is related either to this spin-disc alignment stage or
to the pulsar torque. For the model A, the NS starts as an ejector, and the
rotational evolution for the first $\sim 25$ kyr is indistinguishable from
that of a classic radiopulsar. Model B starts with a propeller stage several
hundred years long. On average, $\chi$ increases by several degrees (up to
$\Delta\chi\sim 15\degr$) during the alignment stage and settles at a constant
level when $\alpha$ stops evolving.

Indeed, the second term in (26) vanishes as $\alpha\to 0$, while the pulsar
torque is inefficient due to the small spin period and buried magnetic field.
The contribution of the pulsar torque may however become important if the
torque is enhanced by the opening of magnetic field lines considered by
Parfrey et al. (2016). We discuss this possibility in Section 4.1.

If the NS initially rotates backwards with respect to the angular momentum of
the disc ($\alpha_{0}>90\degr$), the coupled evolution of $\chi$ starts with
magnetic alignment, that is then (as $\alpha$ becomes smaller) compensated by
rapid orthogonalisation, resulting in identical final magnetic angle values
for $\alpha_{0}$ and $90\degr-\alpha_{0}$. As we will show in Section 4.2,
this is a general result valid as long as we neglect pulsar losses.

In Figure 4, we show the overall change in magnetic angle
$\Delta\chi=\chi-\chi_{0}$ after $10^{5}$ years of evolution for different
$\chi_{0}$ and $\alpha_{0}$. It is clearly seen that for all the models, there
is systematic magnetic orthogonalisation. Under our assumption of relatively
strong accretion torque modulation, magnetic angle $\chi$ changes for up to
$15$ degrees during the spin-disc alignment stage and then remains constant
unless pulsar losses become important.

The mirror symmetry of the plots relative to the $\chi_{0}=90\degr$ line is
related to the symmetry of the magnetic poles. Both positive values of
$\Delta\chi$ at $\chi_{0}<90\degr$ and negative at $\chi_{0}>90\degr$
correspond to the increasing misalignment of the magnetic axis with the axis
of rotation.

Figure 4: Difference in magnetic angles $\Delta\chi=\chi-\chi_{0}$ after
$10^{5}$ years of evolution for different models as function of $\chi_{0}$ and
$\alpha_{0}$. Top left plot is for model A, top right – B, bottom left – C and
bottom right – D respectively (see Table 1 for details.)

In addition, in Figure 5 we show how the difference between final and initial
magnetic angles depends upon the adopted accretion rate $\dot{M}_{1}$. Namely,
we calculated $\Delta\chi$ with the initial parameters of model A, but for
different values of $\dot{M}_{1}$ = 0.01, 0.03, 0.1, 0.5, and 1.0. First of
these values corresponds to the original model A as it is defined in Table 1,
while the last one is identical to model B. The maps shown in Figure 5 outline
the result of the first $10^{5}$ years of evolution. The map changes smoothly
with the increasing mass accretion rate, mostly due to decreasing contribution
from pulsar losses, that in our model tend to compensate the orthogonalisation
during the spin alignment stage. The larger the mass accretion rate, the
smaller is the overall contribution of the pulsar losses. The last two panels
of the figure are practically unaffected by pulsar losses, and show the effect
of the coupled evolution during the alignment stage.

Figure 5: The same as Figure 4 for models with the mass accretion rates
intermediate between models A and B. All the other parameters are the same as
in A/B.

##  4 Discussion

###  4.1 Pulsar losses during accretion

Should pulsar torques remain the same during the accretion stage? On one hand,
accretion is often claimed to interfere with pulsar emission mechanisms (e.g.
Shvartsman, 1971; Melatos & Mastrano, 2016; Papitto & de Martino, 2020). On
the other hand, independently of the physical conditions in the magnetosphere,
open magnetic field lines should still provide a torque determined by the
magnetic stresses at the light cylinder (Contopoulos & Spitkovsky, 2006).

As it was shown by Parfrey et al. (2016) for the case of an aligned rotator,
the presence of a conducting disc opens up some of the closed magnetic field
lines and thus increases the spin-down by a factor proportional to the squared
magnetic flux through the open field lines. For an accreting NS with $r_{\rm
m}\sim r_{\rm co}\ll r_{\rm LC}$, this results in $N_{\rm Parf}\sim N_{\rm
psr}\left(r_{\rm LC}/r_{\rm m}\right)^{2}\sim N_{\rm psr}\left(r_{\rm
LC}/r_{\rm co}\right)^{2}\left(r_{\rm co}/r_{\rm m}\right)^{2}\sim N_{\rm
mag}\left(r_{\rm co}/r_{\rm LC}\right)\left(r_{\rm co}/r_{\rm m}\right)^{2}$.
The contribution of this enhanced spin-down torque is normally smaller than
the estimated magnetospheric losses unless the NS is far from equilibrium
rotation and $r_{\rm m}\ll r_{\rm co}$. But, if the aligning pulsar torque
(the second term in Equation 13) is scaled similarly, magnetic angle very
rapidly evolves to zero during the accreting millisecond pulsar stage.

Figure 6: Magnetic angle evolution when magnetic alignment torque enhancement
described in (Parfrey et al., 2016) is taken into account. For each considered
model from the Table 1, the case of $(\chi_{0},\alpha_{0})$ =
$(60^{\circ},135^{\circ})$ is shown. The evolution without enhancement is
shown by grey dashed lines, while enhancement scenarios are represented by
solid black lines.

Indeed, as pulsar obliquity evolution is driven physically by the same torque
as the pulsar spin-down, one would expect the magnetic alignment time scale
$\tau_{\mathrm{psr}}$ (15) to decrease by the same factor as the spin-down
time scale.

So, torque enhancement decreases it by the factor

$\left(\dfrac{r_{\mathrm{m}}}{r_{\mathrm{LC}}}\right)^{2}\approx 3.5\times
10^{-3}\mu_{30}^{8/7}{\dot{M}_{1}}^{-4/7}\left(\dfrac{M_{*}}{\rm
M_{\sun}}\right)^{-2/7}P_{\rm s}^{-2}.$ (32)

Taking this effect into account leads to faster evolution of $\chi$ in an
accreting system. In Figure 6 we show the evolution of $\chi$ for all the
models from Table 1. Initial conditions are $\chi_{0}=60^{\circ}$ and
$\alpha_{0}=135^{\circ}$ that allows comparison with Figure 3. Quite
expectedly, the impact of the enhancement factor strongly depends on the
initial rotation period and on the magnetic field of the star. If the initial
spin period of the NS is relatively large ($P_{0}\gtrsim 1\,$s, models C and
D), then the star is additionally aligned by $\sim 20\degr$ during the
accretion stage and then stops due to magnetic field burial. On the other
hand, for the relatively small initial period $P_{0}=0.1$ s (models A and B),
magnetic angle decreases effectively to zero before the expected magnetic
field burial time.

We conclude that magnetic alignment torque enhancement, if present, is a very
strong effect. If present on the spin-up stage in MSP progenitors, it should
probably make all the magnetic angles in these sources too small for the MSPs
to be observed as pulsars. However, more detailed analysis of this issue is
required.

###  4.2 Analytic solution for the coupled orthogonalisation process

In our approach, there are two main processes affecting magnetic angle:
magnetic alignment by pulsar torques and $\chi-\alpha$ coupling during the
alignment of the rotation axis with the disc. First effect is well studied and
plays important role in the evolution of normal radiopulsars on their spin-
down time scales (typically $\sim 10^{6}-10^{8}$ years). For accreting
systems, the most important stage of $\chi$ evolution is the alignment of the
rotation axis. The overall change of $\chi$ may be comparable to $\alpha$ but
strongly depends on the modulation within the spin phase. If we neglect the
pulsar torque and exclude time from equations (25) and (26), magnetic angle
becomes effectively a function of $\alpha$ only

$\frac{d\chi}{d\alpha}=-\eta
A(\eta,\alpha,\chi)\sin\alpha\cos\alpha\sin\chi\cos\chi,$ (33)

where $A(\eta,\alpha,\chi)$ is given by Equation (20). Note that, as $A$ does
not change with the changing sign of $\cos\alpha$, and the differential of
$\mathrm{d}\alpha$ appears in combination
$\cos\alpha\sin\alpha\mathrm{d}\alpha$, any solution of Equation (33) will
have the symmetry between $\alpha$ and $\uppi-\alpha$. After the end of the
alignment stage, the solution $\chi(\alpha=0)$ depends only on $\chi_{0}$ and
$\sin\alpha_{0}$, that justifies the independence of the final results on the
sign of $\cos\alpha_{0}$ in Figure 3.

For $\eta\ll 1$, Equation (33) may be integrated analytically, providing an
expression for the evolution of $\chi$

$\displaystyle\tan\chi\simeq\tan\chi_{0}\exp\left\\{\frac{\eta}{2}\left(\cos^{2}\alpha-\cos^{2}\alpha_{0}\right)\right\\}.$
(34)

Or, in terms of $\Delta\chi=\chi-\chi_{0}$

$\Delta\chi\simeq\dfrac{\eta}{2}\cos^{2}\chi_{0}\left(\cos^{2}\alpha-\cos^{2}\alpha_{0}\right).$
(35)

The changes in magnetic angle depend only on initial inclination, initial
magnetic angle, and spin-phase modulation. As cosine is a decreasing function,
$\chi$ increases with alignment if and only if $\alpha<90\degr$.

###  4.3 Spin-disc alignment from the observational point of view

The evolution of a NS in a binary system may consist of one or more $\chi$
alignment episodes driven by pulsar torque, and one or more magnetic angle
leaps due to modulation coupling with $\alpha$. The time scales for these
stages are profoundly different. Pulsar magnetic alignment time is
approximately the same as the pulsar spin-down time (15), while rotation axis
is aligned with the disc on the spin-up time scale (16).

Normally, spin-disc alignment of the NS is a very short stage compared to the
lifetime of the source. Besides, the overall change in magnetic angle is
limited by a quantity of about $\Delta\chi_{\rm
max}\lesssim\frac{\eta}{2}\Delta\alpha$. In favourable conditions, if the
modulation factor $\eta\simeq 1$ and initial $\alpha_{0}\sim 20\div 40$ or
$140\div 160\degr$, we expect the magnetic angle to change by a considerable
amount of 10$\div$15 degrees (see Figure 4).

Apparently, such alignment episodes sometimes happen. An important example is
NGC300 X-1 (Carpano et al., 2018), known as a peculiar Be/X-ray binary with a
NS spun up from about 100 to 20s in several years in an episode of rapid
accretion (Vasilopoulos et al., 2018). Such a strong change in observed period
allows for an equally dramatic change in the direction of the axis of
rotation. In fact, the state in which the system existed before the outburst
of 2010 suggests that the NS was fed by a low-density wind carrying little
angular momentum.

###  4.4 Wind accretion

During wind accretion, the direction of the angular momentum captured by the
NS is not as stable as in the case of accretion through the inner Lagrangian
point. As it was proposed by Wang (1981b) and qualitatively confirmed by
numerical simulations (Fryxell & Taam, 1988; Ruffert, 1997, 1999; MacLeod &
Ramirez-Ruiz, 2015), inhomogeneities in density and velocity of the wind make
the direction of the angular momentum in the disc or in the quasi-spherical
accretion flow (as in Shakura et al. 2012) a random quantity changing on the
time scales much smaller than the alignment time scale. In this case, the
equilibrium rotation period of the NS should be longer than for disc-fed
accretion at equal rate, and there is no real alignment of the rotation axis.
Instead, each separate episode of stable spin-up is associated with evolution
(alignment) of $\alpha$ and related modulation-mediated evolution of $\chi$.

On the time scales smaller than the time scale of the variations in the
inflowing angular momentum, Equation (33) is still valid. On longer time
scales, the evolution of magnetic angle is described by a stochastic equation
following from (26) in the absence of pulsar losses

$\frac{d}{dt}\left(\ln\tan\chi\right)=\eta
A\frac{N_{0}}{I\Omega}\sin^{2}\alpha\cos\alpha$ (36)

where $\alpha$ is now a random function of time. We can think of the evolution
of $\ln\tan\chi$ as a random walk process caused by a white noise in the
right-hand side of Equation (36). Let us consider that the direction of the
external angular momentum $\mn@boldsymbol{\it N}$ is stable on some time scale
$\tau_{\alpha}\ll\tau_{\rm align}$, and completely random and isotropic on
longer times.

If the orientation of the inflowing angular momentum is random and isotropic,
the mean value of the right-hand side of Equation (36) is zero, and there is
no average trend in magnetic angle evolution. However, dispersion of the
derivative of $\ln\tan\chi$ is non-zero and independent of $\chi$, that means
dispersion of the quantity itself grows with time approximately linearly, as

$\begin{array}[]{lcc}\displaystyle\left\langle\left(\ln\tan\chi\right)^{2}\right\rangle&\simeq&\displaystyle\eta^{2}A^{2}\left(\frac{N_{0}}{I\Omega}\right)^{2}\left\langle\left(\cos^{2}\alpha\sin\alpha\right)^{2}\right\rangle\tau_{\alpha}t\\\
\displaystyle&=&\displaystyle\frac{8}{105}\eta^{2}A^{2}\left(\frac{N_{0}}{I\Omega}\right)^{2}\tau_{\alpha}t,\\\
\end{array}$ (37)

where $\tau_{\alpha}$ is the characteristic correlation time scale of
$\alpha$. Magnetic angle is involved in a random walk process that is linear
in $\ln\tan\chi$, and thus tends to form a uniform distribution in this
quantity. Uniform distribution in $\ln\tan\chi$ corresponds to

$\displaystyle\frac{dN}{d\chi}\propto\frac{1}{\left|\sin\chi\cos\chi\right|}.$
(38)

We expect such a stochastic regime to be realised in HMXB systems where wind-
fed accretion is the most probable scenario at low mass accretion rates (see
for instance Chaty 2011). Sufficiently long accretion in this regime will
leave the NS with a random orientation of the rotation axis, but the
distribution in magnetic angle will have much more objects with $\chi$ close
to either zero (or 180°, as the magnetic poles are indistinguishable) or 90°.
Thus, we expect many wind-accreting HMXBs to have nearly-aligned or nearly-
orthogonal NSs.

###  4.5 Magnetic orthogonalisation due to magnetosphere-disc interaction

Lai et al. (2011) have considered interaction of a protoplanetary disc with a
magnetised star. Such a system is physically similar to that considered in
this work. Unlike most of the studies in the field, Lai et al. (2011) assumed
that both the magnetic moment and the angular momentum in the disc are
inclined with respect to the spin axis of the star.

Though the initial assumptions of this paper are based on a simple resistive
model similar to Ghosh et al. (1977), their results give the impression of
what may be the effect of disc interaction with the magnetosphere. The real
penetration depth is likely much smaller than in resistive models, implying
that the effect should also be not so profound.

According to Lai et al. (2011), the spin of the central star is affected by a
back-reaction torque owing to the interaction of the surface currents in the
disc with the poloidal field of the star. The torque may be split into two
components, warping and precession. Warping torque equals

$\mn@boldsymbol{\it
N}_{\mathrm{warp}}=\dfrac{1}{6}\xi^{-7/2}N_{0}\cos^{2}\chi\cos\alpha\
(\mn@boldsymbol{\it d}\times(\mn@boldsymbol{\it s}\times\mn@boldsymbol{\it
d})).$ (39)

After averaging over the spin period, it introduces additional terms to Eqs
(23) and (25),

$I\left\langle\dfrac{d\Omega}{dt}\right\rangle_{\mathrm{warp}}=\dfrac{N_{0}}{6}\xi^{-7/2}\cos^{2}\chi\sin^{2}\alpha\cos\alpha$
(40)

and

$I\Omega\left\langle\dfrac{d\alpha}{dt}\right\rangle_{\mathrm{warp}}=\dfrac{N_{0}}{6}\xi^{-7/2}\cos^{2}\chi\cos^{2}\alpha\sin\alpha,$
(41)

but does not contribute to Equation (26) describing the evolution of magnetic
angle. The other, precessional, back-reaction torque vanishes after averaging
for the same reasons as the other precession torques (see Section 2.2).

The aligning warping torque given by Equation (41) is positive for
$\alpha<\uppi/2$ and thus slows down the alignment of the spin axis.

We have included the warping terms (40) and (41) to our simulations and
additionally run all the four models with them. The resultant NS spin
evolution is shown in Figure 7. It is clearly seen that adding a warping
torque allows $\alpha$ to remain relatively large for a longer time. As the
magnetic angle evolves due to the modulation term at a rate proportional to
$\sin^{2}\alpha\cos\alpha$ (see Equation 26), magnetic angle orthogonalisation
becomes more profound. Depending on the particular initial configuration,
$\Delta\chi$ may reach values several times larger than during conventional
evolution considered in our paper. Maximal values of $|\Delta\chi|\simeq
50\degr$ are observed for initially small magnetic angles $\chi_{0}\lesssim
5-10^{\circ}$ for a broad range of $\alpha_{0}$.

However, the contribution of the warping torque should be treated with
caution, as it relies on the rather poorly known details of magnetic field
interaction with the disc. In particular, it treats the disc as planar, while
a real accretion disc around an inclined rotator, being subject to warping
torques, should become warped. A more elaborate version should take into
account the shape of the disc as well as for a more realistic dynamics of
disc-magnetosphere interaction, involving instabilities and magnetic field
reconnection.

Figure 7: Magnetic angle evolution taking into account back-reaction warp
torque described in (Lai et al., 2011). We show the same set of models as in
Figure 3. Alignment history without warp torque is shown by thick grey lines
in each panel, while warp scenarios are represented by black lines. Dashed
lines are used for $\alpha_{0}=135\degr$, solid for $\alpha_{0}=45\degr$.

##  5 Conclusions

The main issue we addressed in this paper was whether it is possible to change
the magnetic angle of a NS by accretion. We find this possible under quite
specific circumstances: when the rotation axis is misaligned with the spin-up
torque, and the spin-up torque itself is a function of the spin phase. Under
favourable conditions, spin-disc alignment leads to misalignment of the
magnetic axis, but always at a lower rate.

If the direction of the angular momentum of the accreting matter is stable (as
in the case of Roche-lobe overflow), this alignment is only a short episode in
the beginning of accretion. Outside this short alignment stage, evolution of
the magnetic angle is affected only by pulsar aligning torque or its
modifications. It is unclear if the aligning torque should be enhanced due to
the opening of magnetic field lines by the accretion disc.

Accretion from a strongly inhomogeneous stellar wind results in a different
picture. Strong variations of velocity and density make the net angular
momentum of the accreting matter a random quantity. We show that, in this
case, magnetic angle experiences a random walk, related to the random-walk
behaviour of the rotation axis. The rate of this random walk depends strongly
on the magnetic angle itself, leading to a higher probability to find an
object in either aligned or orthogonal state.

In all the scenarios, the rate of magnetic angle alignment is proportional to
$\sin\chi\cos\chi$, suggesting that any distribution in $\chi$ sufficiently
affected by accretion should have an excess of aligned and orthogonal
rotators. However, the full picture is more complex and requires taking into
account the initial distribution in magnetic angle, accretion histories, and
the correct form of the spin-down torque during the accretion stage.

## Acknowledgements

AB acknowledges the Russian Government Program of Competitive Growth of Kazan
Federal University supporting the work.

## Data availability

The data underlying this article will be shared on reasonable request to the
corresponding author.

## References

  * Annala & Poutanen (2010) Annala M., Poutanen J., 2010, A&A, 520, A76
  * Archibald et al. (2009) Archibald A. M., et al., 2009, Science, 324, 1411
  * Beskin & Zheltoukhov (2014) Beskin V. S., Zheltoukhov A. A., 2014, Physics Uspekhi, 57, 799
  * Beskin et al. (1993) Beskin V. S., Gurevich A. V., Istomin Y. N., 1993, Physics of the pulsar magnetosphere. Cambridge, New York: Cambridge University Press
  * Bhattacharya & van den Heuvel (1991) Bhattacharya D., van den Heuvel E. P. J., 1991, Phys. Rep., 203, 1
  * Biryukov et al. (2017a) Biryukov A., Astashenok A., Beskin G., 2017a, MNRAS, 466, 4320
  * Biryukov et al. (2017b) Biryukov A., Astashenok A., Karpov S., Beskin G., 2017b, in Journal of Physics Conference Series. p. 012044 (arXiv:1711.05081), doi:10.1088/1742-6596/932/1/012044
  * Bogdanov et al. (2008) Bogdanov S., Grindlay J. E., Rybicki G. B., 2008, ApJ, 689, 407
  * Bozzo et al. (2018) Bozzo E., Ascenzi S., Ducci L., Papitto A., Burderi L., Stella L., 2018, A&A, 617, A126
  * Caballero & Wilms (2012) Caballero I., Wilms J., 2012, Mem. Soc. Astron. Italiana, 83, 230
  * Carpano et al. (2018) Carpano S., Haberl F., Maitra C., Vasilopoulos G., 2018, MNRAS, 476, L45
  * Casini & Montemayor (1998) Casini H., Montemayor R., 1998, ApJ, 503, 374
  * Chaty (2011) Chaty S., 2011, in Schmidtobreick L., Schreiber M. R., Tappert C., eds, Astronomical Society of the Pacific Conference Series Vol. 447, Evolution of Compact Binaries. p. 29 (arXiv:1107.0231)
  * Chen et al. (1998) Chen K., Ruderman M., Zhu T., 1998, ApJ, 493, 397
  * Cheng et al. (2014) Cheng Z., Zhang C.-M., Zhao Y.-H., Wang D.-H., Pan Y.-Y., Lei Y.-J., 2014, Chinese Astron. Astrophys., 38, 294
  * Contopoulos & Spitkovsky (2006) Contopoulos I., Spitkovsky A., 2006, ApJ, 643, 1139
  * Frank et al. (2002) Frank J., King A., Raine D. J., 2002, Accretion Power in Astrophysics: Third Edition
  * Fryxell & Taam (1988) Fryxell B. A., Taam R. E., 1988, ApJ, 335, 862
  * Ghosh et al. (1977) Ghosh P., Lamb F. K., Pethick C. J., 1977, ApJ, 217, 578
  * Goldreich (1970) Goldreich P., 1970, ApJ, 160, L11
  * Gullón et al. (2014) Gullón M., Miralles J. A., Viganò D., Pons J. A., 2014, MNRAS, 443, 1891
  * Haskell et al. (2008) Haskell B., Samuelsson L., Glampedakis K., Andersson N., 2008, MNRAS, 385, 531
  * Hessels (2008) Hessels J. W. T., 2008, in Wijnands R., Altamirano D., Soleri P., Degenaar N., Rea N., Casella P., Patruno A., Linares M., eds, American Institute of Physics Conference Series Vol. 1068, American Institute of Physics Conference Series. pp 130–134 (arXiv:0903.0493), doi:10.1063/1.3031183
  * Igoshev & Popov (2015) Igoshev A. P., Popov S. B., 2015, Astronomische Nachrichten, 336, 831
  * Illarionov & Kompaneets (1990) Illarionov A. F., Kompaneets D. A., 1990, MNRAS, 247, 219
  * Illarionov & Sunyaev (1975) Illarionov A. F., Sunyaev R. A., 1975, A&A, 39, 185
  * Johnson et al. (2014) Johnson T. J., et al., 2014, ApJS, 213, 6
  * Lai (2014) Lai D., 2014, in European Physical Journal Web of Conferences. p. 01001 (arXiv:1402.1903), doi:10.1051/epjconf/20136401001
  * Lai et al. (2011) Lai D., Foucart F., Lin D. N. C., 2011, MNRAS, 412, 2790
  * Landau et al. (1976) Landau L., Lifshitz E., Sykes J., Bell J., 1976, Mechanics. Butterworth-Heinemann, Elsevier Science, https://books.google.fi/books?id=e-xASAehg1sC
  * Leahy (1990) Leahy D. A., 1990, MNRAS, 242, 188
  * Lipunov & Shakura (1980) Lipunov V. M., Shakura N. I., 1980, Soviet Astronomy Letters, 6, 14
  * Lovelace et al. (1995) Lovelace R. V. E., Romanova M. M., Bisnovatyi-Kogan G. S., 1995, MNRAS, 275, 244
  * Lyne & Manchester (1988) Lyne A. G., Manchester R. N., 1988, MNRAS, 234, 477
  * MacLeod & Ramirez-Ruiz (2015) MacLeod M., Ramirez-Ruiz E., 2015, ApJ, 803, 41
  * Malov & Nikitina (2011) Malov I. F., Nikitina E. B., 2011, Astronomy Reports, 55, 878
  * Manchester (2017) Manchester R. N., 2017, Journal of Astrophysics and Astronomy, 38, 42
  * Melatos (2000) Melatos A., 2000, MNRAS, 313, 217
  * Melatos & Mastrano (2016) Melatos A., Mastrano A., 2016, ApJ, 818, 49
  * Melatos & Phinney (2001) Melatos A., Phinney E. S., 2001, Publications of the Astronomical Society of Australia, 18, 421
  * Novoselov et al. (2020) Novoselov E. M., Beskin V. S., Galishnikova A. K., Rashkovetskyi M. M., Biryukov A. V., 2020, MNRAS, 494, 3899
  * Pan et al. (2013) Pan Y., Zhang C., Wang N., 2013, in Zhang C. M., Belloni T., Méndez M., Zhang S. N., eds, Proceedings of the International Astronomical Union Vol. 290, Feeding Compact Objects: Accretion on All Scales. pp 291–292 (arXiv:1304.0073), doi:10.1017/S1743921312020066
  * Papitto & de Martino (2020) Papitto A., de Martino D., 2020, arXiv e-prints, p. arXiv:2010.09060
  * Parfrey et al. (2016) Parfrey K., Spitkovsky A., Beloborodov A. M., 2016, ApJ, 822, 33
  * Parfrey et al. (2017) Parfrey K., Spitkovsky A., Beloborodov A. M., 2017, MNRAS, 469, 3656
  * Patruno & Watts (2021) Patruno A., Watts A. L., 2021, Astrophysics and Space Science Library, 461, 143
  * Pétri (2019) Pétri J., 2019, MNRAS, 485, 4573
  * Philippov et al. (2014) Philippov A., Tchekhovskoy A., Li J. G., 2014, MNRAS, 441, 1879
  * Rankin (1993) Rankin J. M., 1993, Astrophys. J., 405, 285
  * Rappaport et al. (2004) Rappaport S. A., Fregeau J. M., Spruit H., 2004, ApJ, 606, 436
  * Reisenegger (2003) Reisenegger A., 2003, ArXiv Astrophysics e-prints,
  * Romanova et al. (2009) Romanova M. M., Ustyugova G. V., Koldoba A. V., Lovelace R. V. E., 2009, MNRAS, 399, 1802
  * Romanova et al. (2020) Romanova M. M., Koldoba A. V., Ustyugova G. V., Blinova A. A., Lai D., Lovelace R. V. E., 2020, arXiv e-prints, p. arXiv:2012.10826
  * Ruderman (1991) Ruderman R., 1991, ApJ, 382, 576
  * Ruffert (1997) Ruffert M., 1997, A&A, 317, 793
  * Ruffert (1999) Ruffert M., 1999, A&A, 346, 861
  * Shakura et al. (2012) Shakura N., Postnov K., Kochetkova A., Hjalmarsdotter L., 2012, MNRAS, 420, 216
  * Shvartsman (1970) Shvartsman V. F., 1970, Izvestiia Vysshaia Uchebn. Zaved., Radiofizika, 13, 1852
  * Shvartsman (1971) Shvartsman V. F., 1971, Soviet Ast., 15, 342
  * Spitkovsky (2006) Spitkovsky A., 2006, Astrophys. J. Lett., 648, L51
  * Tauris & Manchester (1998) Tauris T. M., Manchester R. N., 1998, MNRAS, 298, 625
  * Tauris & van den Heuvel (2006) Tauris T. M., van den Heuvel E. P. J., 2006, Formation and evolution of compact stellar X-ray sources. pp 623–665
  * Tauris et al. (2012) Tauris T. M., Langer N., Kramer M., 2012, MNRAS, 425, 1601
  * Vasilopoulos et al. (2018) Vasilopoulos G., Haberl F., Carpano S., Maitra C., 2018, A&A, 620, L12
  * Veledina et al. (2019) Veledina A., Nättilä J., Beloborodov A. M., 2019, ApJ, 884, 144
  * Wang (1981a) Wang Y. M., 1981a, Space Sci. Rev., 30, 341
  * Wang (1981b) Wang Y.-M., 1981b, A&A, 102, 36
  * Zanni & Ferreira (2013) Zanni C., Ferreira J., 2013, A&A, 550, A99
  * Zhang et al. (2011) Zhang C. M., et al., 2011, A&A, 527, A83

##  Appendix A Spheroidal neutron star

Let us consider a spheroidal NS having variable ellipticity $\epsilon$
(positive or negative) along an arbitrary axis $\mn@boldsymbol{\it l}$
($|\mn@boldsymbol{\it l}|=1$). Then the equation for the rotational evolution
of such a NS may be written as

$I\dfrac{\mathrm{d}\mn@boldsymbol{\it\Omega}}{\mathrm{d}\it
t}=\mn@boldsymbol{\it N}-\dfrac{\mathrm{d}\it I}{\mathrm{d}\it
t}\mn@boldsymbol{\it\Omega}-\epsilon
I\Omega_{l}\left[\left(\dfrac{\dot{I}}{I}+\dfrac{\dot{\Omega}_{l}}{\Omega_{l}}+\dfrac{\dot{\epsilon}}{\epsilon}\right)\mn@boldsymbol{\it
l}+\mn@boldsymbol{\it\Omega}\times\mn@boldsymbol{\it l}\right],$ (42)

where $\Omega_{l}$ is the scalar product
$\mn@boldsymbol{\it\Omega}\cdot\mn@boldsymbol{\it l}$.

The right-hand side of Equation (42) consists of three terms: external torque
$\mn@boldsymbol{\it N}$, the term proportional to $\mathrm{d}{\it
I}/\mathrm{d}{\it t}$ and related to the variable moment of inertia of the
star, and the term arising from non-sphericity. The latter can be also
rewritten as

${\mn@boldsymbol{\it N}}_{\epsilon}=-(\epsilon\Omega_{l}\dot{I}+\epsilon
I\dot{\Omega}_{l}+\dot{\epsilon}I\Omega_{l})\mn@boldsymbol{\it l}+\epsilon
I\Omega_{l}(\mn@boldsymbol{\it\Omega}\times\mn@boldsymbol{\it l}).$ (43)

Since the ellipticity of a neutron star is expected to be very small (i.e.
$|\epsilon|\ll 1$, see for instance Haskell et al. 2008) and changes slowly,

$|\epsilon\Omega_{l}\dot{I}|\ll|\Omega\dot{I}|$ (44)

and

$|\epsilon I\dot{\Omega}_{l}|\ll|I\dot{\Omega}|\simeq
I\Omega\tau^{-1}_{\mathrm{align}},$ (45)

as the alignment time $\tau_{\mathrm{align}}$ and the spin-up time scale are
determined by the same external torque $\mn@boldsymbol{\it N}$.

Therefore, the first term in $\mn@boldsymbol{\it N}_{\epsilon}$, aligned with
the axis $\mn@boldsymbol{\it l}$, may be omitted while considering the
evolution of an accreting neutron star during spin-disc alignment, as the
contribution of the alignment itself is always larger.

If the deformation is aligned with the magnetic moment (or, equivalently,
$\mn@boldsymbol{\it l}=\mn@boldsymbol{\it m}$), the second term in Equation
(43), proportional to $\mn@boldsymbol{\it\Omega}\times\mn@boldsymbol{\it l}$,
does not contribute to the evolutionary equations for the reasons discussed in
Section 2.2 for the anomalous pulsar torque and valid for any torque
orthogonal to the magnetic and spin axes at the same time.

The only possible case when the precession term may become important for the
long-term evolution of the rotation parameters is the case when
$\mn@boldsymbol{\it l}$ is not co-planar with the spin and magnetic axes.

But this situation requires a triaxial rather than a spheroidal star (Melatos,
2000), which is beyond of the scope of current paper.