\section{Introduction}

Lorentz invariance, which claims \Reply{the equivalence between any two inertial reference frames in
formulating the laws of physics}, is a fundamental principle in both general
relativity (GR) and the Standard Model of particle physics. However, in
pursuing a unified theory of gravity and quantum particles, violation of
Lorentz invariance was suggested \cite{Kostelecky:1988zi,
Kostelecky:1989jp, Kostelecky:1991ak, Will:1972zz, Will:2018bme}, and has
been treated as a possible suppressed effect emerging from the Planck scale
where the unified theory lives. Therefore, searching for Lorentz-invariance
violation in high-precision terrestrial experiments and astrophysical
observations~\cite{Kostelecky:2008ts, Will:2018bme, Shao:2016ezh} not only
serves as a necessary test of our current best knowledge of the laws of
physics, but also provides us the chance to obtain information on the
underlying theory of quantum gravity that is otherwise unattainable
directly in the experiments and observations with limited energy scales.

A practical framework to study Lorentz violation without digging into the
extensive intricacies of the underlying theory is the Standard-Model
Extension (SME) developed since 1990s \cite{Colladay:1996iz,
Colladay:1998fq, Kostelecky:2003fs}. Aiming at guiding experimental
searches of Lorentz violation, the SME framework is constructed at the
level of effective field theory and treats all possible Lorentz-violating
operators as perturbations on top of GR and the Standard Model. In this
work we are studying effects of Lorentz violation on \Reply{the rotation} of stars.
This belongs to the gravitational sector of the SME framework, and for
simplicity, we only consider the modification to gravity generated by the
so-called minimal gravitational SME \cite{Bailey:2006fd}, given by the
action
\begin{eqnarray}
S = \frac{1}{16\pi} \int \sqrt{-g} \, d^4x \, \left( R + k^{\alpha\beta\gamma\delta} R_{\alpha\beta\gamma\delta} \right) + S_k + S_m,
\label{action}
\end{eqnarray}
where the Ricci scalar $R$ represents the usual Einstein-Hilbert term for
GR, $R_{\alpha\beta\gamma\delta}$ is the Riemann tensor, and $k^{\alpha\beta\gamma\delta}$ is the
tensor field that breaks local Lorentz invariance when it acquires a
nonzero vacuum expectation value in an underlying theory. The action $S_k$
describes the dynamics of the Lorentz-violating field $k^{\alpha\beta\gamma\delta}$ at
the level of effective field theory. The symmetry-breaking mechanism is
important but generally unspecified~\cite{Kostelecky:2003fs, Bluhm:2004ep}.
One of the central tasks in the gravitational SME is obtaining approximate
but general expressions for the contribution of $S_k$ in the field
equations using geometric properties of the spacetime manifold like the
Bianchi identity and the diffeomorphism invariance \cite{Bailey:2006fd,
Bailey:2014bta}. The last term $S_m$ is the action of conventional matter
and in this work we consider it Lorentz invariant. Note that the
geometrized unit system where $G=c=1$ is used. We will stick to this unit
system except when units are specified explicitly. Another explanation
about the notations in the remaining of the paper is that repeated indices
are summed even when they both are subscripts or superscripts. Also note
that the Greek letters run over spacetime indices while the Latin letters
are restricted to spatial indices only.

The weak field solution of the metric to the field equations obtained from
action \rf{action} is calculated by \citet{Bailey:2006fd}. Especially, as
the starting point of our study, the anisotropic modification to the
Newtonian potential energy between two point particles, $A$ and $B$, is
\begin{eqnarray}
\delta U = -\frac{1}{2} \bar s^{ij} \frac{ m_A m_B}{|{\boldsymbol{x}}_A - {\boldsymbol{x}}_B|^3} (x_A^i-x_B^i) (x_A^j - x_B^j) ,
\label{pecorrect}
\end{eqnarray}
where $x_A^i$ and $x_B^i$ are the components of the position vectors
${\boldsymbol{x}}_A$ and ${\boldsymbol{x}}_B$. The quantities $\bar s^{ij}$
are defined as the spatial components of
\begin{eqnarray}
\bar s^{\alpha\beta} \equiv 2 \left( \bar k^{\alpha\mu\beta\nu} g_{{\mu\nu}} - \frac{1}{4} g^{\alpha\beta} \bar k^{\lambda\mu\kappa\nu} g_{{\mu\nu}} g_{\lambda\kappa} \right),
\end{eqnarray}
with $\bar k^{\alpha\beta\gamma\delta}$ being the vacuum expectation value of the
tensor field $k^{\alpha\beta\gamma\delta}$~\cite{Kostelecky:2003fs}. In practical
calculations aiming at testing the theoretical predictions against
experimental results, the vacuum expectation value $\bar k^{\alpha\beta\gamma\delta}$
is taken to be constant around the experimental setting in approximate
inertial frames, and its components, as well as the combinations $\bar
s^{\alpha\beta}$, are called coefficients for Lorentz violation.

Stringent constraints have been set on the SME coefficients for Lorentz
violation using an extensive number of laboratory experiments and
astrophysical observations \cite{Kostelecky:2008ts}. In fact, they include
results from considering effects of Eq.~\rf{pecorrect} on \Reply{the rotation} of the
Sun \cite{Bailey:2006fd} and two isolated millisecond pulsars
\cite{Shao:2014oha}. In this work, \Reply{we are going to consider the two isolated millisecond pulsars in Ref.~\cite{Shao:2014oha} and set new constraints on the coefficients $\bar s^{ij}$.} Our work
complements the results in the literature as we perform a rigorous study on
the rotation of a spheroidal star under the anisotropic gravitational
self-energy caused by Eq.~\rf{pecorrect} to serve as the theoretical basis.
In addition, our calculation clarifies
the fact that the effect of Lorentz violation adds to stationary spinning
stars as well as free-precessing stars. While the former case is the
Lorentz-violating precession considered in the literature and has been used
to set constraints on the coefficients for Lorentz violation, the latter
case has not been investigated to our knowledge. When Lorentz violation
modifies the rotation of an otherwise free-precessing neutron star (NS),
the pulsar signal and associated continuous gravitational waves (GWs)
emitted by the star~\cite{Gao:2020zcd} will change accordingly. A
significant part of our work is devoted to this new topic.

Before ending the introduction, we are obligated to point out that the
gravitational sector of the SME is not the only framework to study Lorentz
violation in gravity. The celebrated parameterized post-Newtonian (PPN)
formalism also includes coefficients describing possible preferred-frame
effects in metric gravitational theories \cite{Will:2018bme, Will:2014kxa}.
The first study of Lorentz-violating effects on the spins of the Sun and of
millisecond pulsars was done by \citet{Nordtvedt:1987} considering the specific PPN
modification\footnote{Notice that we follow the usage of $\alpha_2$ by
\citet{Will:2018bme}, related to the one used by \citet{Nordtvedt:1987}
via $\alpha_2 = 2 \alpha^{\rm Nordtvedt}_2$. }
\begin{eqnarray}
\delta U_{\rm PPN} = \frac{\alpha_2}{2} \frac{ m_A m_B}{|{\boldsymbol{x}}_A - {\boldsymbol{x}}_B|^3} w^{i} w^j (x_A^i-x_B^i) (x_A^j - x_B^j),
\label{ppndeu}
\end{eqnarray}
to the Newtonian potential. The coefficient $ \alpha_2$ controls the size of
Lorentz violation while the ``absolute'' velocity $\boldsymbol{w}$ with
repect to the preferred inertial frame picks up a special direction and
breaks Lorentz invariance~\cite{Will:1972zz}. Our study is originally
inspired by Nordtvedt's work~\cite{Nordtvedt:1987}, and our result recovers
his solution when the replacement~\cite{Bailey:2006fd,Shao:2014oha}
\begin{eqnarray}
\bar s^{ij} \rightarrow -\alpha_2 w^i w^j ,
\label{replace}
\end{eqnarray}
is applied and proper approximations are made. The SME framework is more
generic than the PPN formalism in terms of Lorentz
violation~\cite{Bailey:2006fd}. Detailed comparisons will be presented in
relevant paragraphs of the paper.

The organization of the paper is as follows. We start with calculating the
anisotropic gravitational self-energy due to Eq.~\rf{pecorrect} for a
spheroidal star in Sec.~\ref{sec2a}. Then, in Sec.~\ref{sec2b}, we
investigate the solutions to the rotational equations of motion thoroughly,
where both perturbation method and numerical calculation are employed. The
observational consequences are discussed in Sec.~\ref{sec3} with regard to NSs. Section \ref{sec3a} considers stationary spinning stars
affected by Lorentz violation and obtains constraints on the coefficients
for Lorentz violation from observations of two solitary millisecond pulsars. Sections \ref{sec3b} and
\ref{sec3c} consider free-precessing stars affected by Lorentz violation
and provide preliminary signal templates of pulsar pulses and continuous
GWs for fitting observational data in the future. Finally, conclusions are
summarized in Sec.~\ref{sec4}. Appendix~\ref{app1} displays some
expressions for uniform spheroids, which are useful for estimating
numerical values.

\section{rotation of a spheroid under Lorentz-violating gravity}
\label{sec2}

\subsection{Anisotropic gravitational self-energy}
\label{sec2a}

\begin{figure}
    \centering
    \includegraphics[width=6cm]{fig_eulerangles1}
    \caption{Euler angles $\alpha, \, \beta$ and $\gamma$ connecting the
    $X$-$Y$-$Z$ frame and the $x$-$y$-$z$ frame. The line $MN$ is the
    intersection of the $X$-$Y$ plane and the $x$-$y$ plane. }
    \label{euler1}
\end{figure}

For a nonspherical star, the integral of the potential energy correction in
Eq.~\rf{pecorrect} depends on the orientation of the star, causing a torque
during its rotation. Specifically speaking, we calculate $\delta U$ in the
body frame of the star at any instant,
\begin{eqnarray}
\delta U = -\frac{1}{4} \bar s^{ij} \int d^3x \, d^3x' \rho(\boldsymbol{x}) \rho(\boldsymbol{x'}) \frac{ \left( x^i - x^{\prime\,i} \right) \left( x^j - x^{\prime\,j} \right) }{|\boldsymbol{x} - \boldsymbol{x'}|^3} ,
\label{pecorrect1}
\end{eqnarray}
where $\rho$ is the density of the star, and is assumed to be independent of
time in the body frame. The orientation dependence goes into $\delta U$ as the
star rotates. If we set up an inertial frame $X$-$Y$-$Z$, then $\bar
s^{ij}$ $(i,j=x,y,z)$ in the body frame is related to $\bar s^{IJ}$
$(I,J=X,Y,Z)$ in the inertial frame by a rotation transformation
\begin{eqnarray}
\bar s^{ij} = R^{iI} R^{jJ} \bar s^{IJ} .
\end{eqnarray}
Noticing that $\bar s^{IJ}$ are the constant coefficients for Lorentz
violation, $\delta U$ therefore depends on the orientation of the body through
the rotation matrix $R^{iI}$. The orientation of the body is conveniently
described by the Euler angles (see Fig. \ref{euler1}), which are the
kinematic quantities to be solved from the equations of motion. For later
use, we write the elements of the rotation matrix as
\begin{eqnarray}
R^{iI} = \hat {\boldsymbol{e}}_i \cdot \hat {\boldsymbol{e}}_I ,
\label{rmat}
\end{eqnarray}
where $\{ \hat {\boldsymbol{e}}_i \}$ $(i=x,y,z)$ is the basis of the body
frame $x$-$y$-$z$, and $ \{ \hat {\boldsymbol{e}}_I \}$ $(I=X,Y,Z)$ is the
basis of the inertial frame $X$-$Y$-$Z$. The inner products can be computed
using the relations~\cite{landau1960course}
\begin{eqnarray}
\hat {\boldsymbol{e}}_x &=& \left(\cos\alpha\cos\gamma-\sin\alpha\cos\beta\sin\gamma \right) \hat {\boldsymbol{e}}_X
\nonumber \\
 && + \left(\sin\alpha\cos\gamma+\cos\alpha\cos\beta\sin\gamma \right) \hat {\boldsymbol{e}}_Y + \sin\beta\sin\gamma \,\hat {\boldsymbol{e}}_Z ,
\nonumber \\
\hat {\boldsymbol{e}}_y &=& \left(-\cos\alpha\sin\gamma-\sin\alpha\cos\beta\cos\gamma \right) \hat {\boldsymbol{e}}_X
\nonumber \\
 && + \left(-\sin\alpha\sin\gamma+\cos\alpha\cos\beta\cos\gamma \right) \hat {\boldsymbol{e}}_Y + \sin\beta\cos\gamma \,\hat {\boldsymbol{e}}_Z ,
\nonumber \\
\hat {\boldsymbol{e}}_z &=& \sin\alpha\sin\beta \,\hat {\boldsymbol{e}}_X - \cos\alpha\sin\beta \,\hat {\boldsymbol{e}}_Y + \cos\beta \,\hat {\boldsymbol{e}}_Z .
\label{trans}
\end{eqnarray}

Back to the potential energy correction in Eq.~\rf{pecorrect1}, to proceed
the calculation analytically as much as possible, we focus on the simple
case where the star is a spheroid described by
\begin{eqnarray}
\frac{x^{2} + y^{2}}{a_1^2} + \frac{z^{2}}{a_3^2}  = 1 ,
\label{spheroid}
\end{eqnarray}
in the body frame and its density $\rho$ is axisymmetric about the $z$-axis (for an ellipsoid, see e.g., Ref.~\cite{Xu:2021dcw}). In such case, one can show that $\bar
s^{xy}, \, \bar s^{xz}$ and $\bar s^{yz}$ do not appear in $\delta U$, and
that $\bar s^{xx}$ and $\bar s^{yy}$ appear as the combination $\bar s^{xx}
+ \bar s^{yy}$. Because the trace $\bar s^{xx} + \bar s^{yy} + \bar s^{zz}$
is rotationally invariant and hence does not contribute to the anisotropic
correction, the contribution from $\bar s^{xx} + \bar s^{yy}$ can be
represented by $\bar s^{zz}$ so that the true anisotropic correction of the
potential energy is simply
\begin{eqnarray}
\delta U = C \, \bar s^{zz},
\label{pecorrect2}
\end{eqnarray}
with the constant $C$ being
\begin{eqnarray}
C = \frac{1}{4} \int d^3x \, d^3x' \rho(\boldsymbol{x}) \rho(\boldsymbol{x'}) \frac{ \left( x - x' \right)^2 - \left( z - z' \right)^2 }{|\boldsymbol{x} - \boldsymbol{x'}|^3} .
\end{eqnarray}

The result \rf{pecorrect2} deserves a comparison with that of
\citet{Nordtvedt:1987}. Treating the potential correction \rf{ppndeu} as a
perturbation on top of the Newtonian potential, \citet{Nordtvedt:1987} used
the tensor virial relation to obtain the anisotropic gravitational
self-energy for a spheroid star as
\begin{eqnarray}
\delta U_{\rm PPN} = -\frac{1}{2} \alpha_2 T_{\rm rot} ( \boldsymbol{w} \cdot \hat{\boldsymbol{\Omega}} )^2 ,
\label{ppndeu1}
\end{eqnarray}
where $T_{\rm rot}$ is the rotational kinetic energy and
$\hat{\boldsymbol{\Omega}}$ is the unit vector along the angular velocity of
the star.

First of all, we point out that with the replacement \rf{replace},
Eq.~\rf{pecorrect2} recovers the form of Eq.~\rf{ppndeu1} when the star
spins stationarily along the $z$-axis. As Nordtvedt obtained, for
nonrelativistic stars, the constant $C$ is equal to $T_{\rm rot}/2$ by
virtue of the tensor virial relation
\begin{eqnarray}
K^{ij} + \delta^{ij} P - U^{ij} = 0,
\end{eqnarray}
where $\delta^{ij}$ is the Kronecker delta, and
\begin{eqnarray}
K^{ij} &\equiv& \int \rho v^i v^j d^3x,
\nonumber \\
P &\equiv& \int p d^3x ,
\nonumber \\
U^{ij} &\equiv& \int \rho x^i \partial_j \Phi d^3x,
\end{eqnarray}
with ${\boldsymbol{v}}$ being the velocity field inside the star, $p$ being
the pressure inside the star, and $\Phi$ being the usual Newtonian
potential
\begin{eqnarray}
\Phi ({\boldsymbol{x}}) = -\int d^3x' \frac{ \rho(\boldsymbol{x'}) }{|\boldsymbol{x} - \boldsymbol{x'}|} .
\label{newp}
\end{eqnarray}
The result $C = T_{\rm rot}/2$ is straightforward to prove once we notice
the integral identity
\begin{eqnarray}
U^{ij} = \frac{1}{2} \int d^3x \, d^3x' \rho(\boldsymbol{x}) \rho(\boldsymbol{x'}) \frac{ \left( x^i - x^{\prime\,i} \right) \left( x^j - x^{\prime\,j} \right) }{|\boldsymbol{x} - \boldsymbol{x'}|^3},
\end{eqnarray}
and for a star spinning along the $z$-direction,
\begin{eqnarray}
K^{zz} = 0, \quad T_{\rm rot} = \frac{1}{2} \left( K^{xx} + K^{yy} \right) .
\end{eqnarray}
In Appendix~\ref{app1}, the constant $C$ is exhibited for uniform spheroids
by utilizing the Newtonian potential inside an ellipsoid of uniform density
(see e.g., Refs.~\cite{Poisson:2014, 1962ApJ...136.1037C}).

Secondly, we notice that Eq.~\rf{pecorrect2} and Eq.~\rf{ppndeu1} are
different for a generally rotating spheroidal star whose angular velocity
is not aligned with its symmetric axis. The origin of the difference is
that Eq.~\rf{ppndeu1} requires the star to be a stationary fluid in
equilibrium while Eq.~\rf{pecorrect2} applies as long as the star has an
axisymmetric mass distribution. Whether the star is a fluid or a rigid body
does not affect the validity of Eq.~\rf{pecorrect2}. Now comes the
question: How could the state of the star influence the alignment between
its angular velocity and symmetric axis? In fact, if the star is a rigid
body, then the direction of its angular velocity can be different from
its symmetric axis. But when it is a fluid, the deformation is caused by
rotation, and therefore its symmetric axis must be aligned with its angular
velocity as long as the system is in equilibrium.

To better understand the context and the limitation for Eq.~\rf{ppndeu1},
let us consider a freely rotating star, namely when there is no torque. If
modeled as a rigid body, depending on whether the angular velocity is
aligned with the symmetric axis, two kinds of solution exist: the
stationary spinning solution where the star spins around the fixed
symmetric axis, and the free-precessing solution where the star spins
around its symmetric axis while that axis rotates around the fixed
direction of the conserved angular momentum~\cite{landau1960course}. If
modeled as a fluid, the free-precessing solution cannot exist because even
if the angular velocity is not aligned with the symmetric axis initially,
the deformation caused by rotation will eventually change the symmetric
axis to the direction of angular velocity, leaving the star in the state of
stationary spin.

When the torque caused by Lorentz violation is taken into consideration,
the two types of solution for a freely rotating star, namely the stationary
spinning solution and the free-precessing solution, might be treated as the
zeroth-order solution for applying perturbation method. The details are
discussed in Sec.~\ref{sec2b}. Here we just point out that the effect of
the Lorentz-violating torque is forcing the star to precess around one of
the principal axes of the $\bar s^{ij}$ tensor. When this effect is added
to a stationary spinning star, we call the motion {\it{single
precession}}. When the forced precession due to Lorentz violation is added
to a free-precessing star, we call the motion {\it{twofold precession}}.

In the single-precession solution, the precession is caused by Lorentz
violation, which is assumed to be small. Therefore, the total angular
velocity can be approximated as aligned with the symmetric axis of the
star. This enables us to put Eq.~\rf{pecorrect2} in the form of
Eq.~\rf{ppndeu1} at the leading order of Lorentz violation via the
replacement \rf{replace}. So the single-precession solution is the one
found by \citet{Nordtvedt:1987}. However, for our result,
there is no need to restrict it to fluid stars as Eq.~\rf{pecorrect2}
applies to fluids as well as rigid bodies. The other solution that we
defined, the twofold-precession solution, as its name suggests, is a
superposition of the free precession where the symmetric axis of the star
precesses around the angular momentum, and the forced precession where the
angular momentum precesses around one of the principal axes of the $\bar
s^{ij}$ tensor. Though only applicable to stars if they are modeled as rigid
bodies, the twofold-precession solution describes a new Lorentz-violating
effect in the rotation of stars that has not been studied in the literature
to our knowledge.

Before moving to study the solutions in detail, \Reply{a brief clarification on
the model of NS deformation might be useful.} We will adopt the
model described by \citet{2001MNRAS.324..811J} where the deformation of the
star has two contributions: the centrifugal deformation and the Coulomb
deformation. The centrifugal deformation is the fluid deformation that
scales with the square of the angular velocity, while the Coulomb
deformation describes the rigid deviation from the spherical shape
sustained by the electrostatic force. The fact that the electrostatic
interaction is much weaker than the gravitational interaction in NSs
implies that the rigid Coulomb deformation is much smaller than the fluid
centrifugal deformation. Using the oblateness defined as
\begin{eqnarray}
\epsilon = \frac{ I^{zz} - I^{xx} }{ I^{xx} } ,
\label{obl}
\end{eqnarray}
where $I^{xx}$ and $I^{zz}$ are the eigenvalues of the moment of inertia tensor,
to characterize the deformations of spheroidal stars, an estimation for the centrifugal
deformation is
\begin{eqnarray}
\hspace{-0.5cm}
\epsilon_{f} \approx \frac{|{\boldsymbol{\Omega}}|^2 R^3}{M} \approx 2.1 \times 10^{-3} \left( \frac{{|{\boldsymbol{\Omega}}|}/{2\pi} }{100\,{\rm Hz} } \right)^2 \left( \frac{R}{10 \, {\rm km}} \right)^3 \left( \frac{1.4M_{\odot}}{M} \right) ,
\label{epf}
\end{eqnarray}
where ${\boldsymbol{\Omega}}$ is the angular velocity, $R$ is the radius, and
$M$ is the mass of the NS, while the oblateness $\epsilon_r$ caused by the
Coulomb deformation is about five orders of magnitude smaller
\cite{2001MNRAS.324..811J}. Therefore, when applying the single-precession
solution to NSs, the rigid deformation $\epsilon_r$ can be neglected and the
deformation is described by the fluid oblateness $\epsilon_f$. However, the fact
that NSs can have rigid deformation described by $\epsilon_r$ is vital when the
twofold-precession solution is to apply. Similar to the free-precessing NSs
considered in Refs.~\cite{Zimmermann:1979ip, 2001MNRAS.324..811J,
Gao:2020zcd, Gao:2020dwy}, NSs with Lorentz-violating gravity under the
twofold precession produce modulated pulsar signals and in the meantime are
sources of continuous GWs, providing potential new tests of Lorentz
violation. This will be the central topic in Sec.~\ref{sec3}.

Finally, we point out that Lorentz violation itself also deforms
stars~\cite{Altschul:2006uu, Xu:2019gua}. But then the deformation due to Lorentz violation
contributes at the second order in terms of the coefficients for Lorentz
violation to the anisotropic gravitational self-energy. To keep the
analysis clear and tractable, we neglect any Lorentz-violating correction
to the structure of stars in this work.

\subsection{Rotation of a spheroidal star}
\label{sec2b}

We are ready to write down the equations of motion and solve the rotation
of a spheroidal star. First, the coefficients for Lorentz violation
naturally fix a convenient inertial frame, namely the one that diagonalizes
the $\bar s^{IJ}$ matrix. We will use it as the inertial frame $X$-$Y$-$Z$. Please note that the inertial frame widely used in the literature is the canonical Sun-centered
celestial-equatorial frame~\cite{Kostelecky:2002hh}. It is generally different from the frame used here as the off-diagonal coefficients for Lorentz violation unlikely happen to vanish in the canonical Sun-centered celestial-equatorial frame.

With the rotation matrix~\rf{rmat}, the anisotropic
gravitational self-energy \rf{pecorrect2} is
\begin{eqnarray}
\delta U = C \left( (\bar s^{XX} \sin^2\alpha + \bar s^{YY} \cos^2\alpha) \sin^2\beta + \bar s^{ZZ} \cos^2\beta \right) ,
\label{anisopot}
\end{eqnarray}
where $\bar s^{XX}, \, \bar
s^{YY}$, and $\bar s^{ZZ}$ are the eigenvalues of the $\bar s^{IJ}$ matrix.
Then, to express the rotational kinetic energy in terms of the Euler angles
and their derivatives, we employ the kinematic relation between the
velocity components in the body frame and the Euler
angles~\cite{landau1960course},
\begin{eqnarray}
\Omega^x &=& \dot \alpha \sin\beta \sin\gamma + \dot \beta \cos\gamma ,
\nonumber \\
\Omega^y &=& \dot \alpha \sin\beta \cos\gamma - \dot \beta \sin\gamma ,
\nonumber \\
\Omega^z &=& \dot \alpha \cos\beta + \dot \gamma,
\label{angvel1}
\end{eqnarray}
where dot denotes time derivative. With $I^{xx} = I^{yy}$, the rotational
kinetic energy simplifies to
\begin{eqnarray}
T_{\rm rot} = \frac{1}{2} I^{xx} \left( \dot \alpha^2 \sin^2\beta + \dot \beta^2 \right) + \frac{1}{2} I^{zz} \left( \dot \alpha\cos\beta + \dot \gamma \right)^2 .
\end{eqnarray}
The Euler-Lagrange equations from $L = T_{\rm rot} - \delta U$ can be derived,
and there are two first integrals,
\begin{eqnarray}
E &=& T_{\rm rot} + \delta U  ,
\nonumber \\
\Omega^z &=& \dot \alpha \cos\beta + \dot \gamma .
\label{firstint}
\end{eqnarray}
The total energy $E$ is not much useful in helping simplify the other
equations. The $z$-component of the angular velocity, $\Omega^z$, can be used
to eliminate $\dot \gamma$ from the other two Euler-Lagrange equations so that
they become
\begin{eqnarray}
\sin^2\beta \, \ddot \alpha + \sin{2\beta} \, \dot \alpha \dot \beta - \frac{I^{zz} \Omega^z }{I^{xx}} \sin\beta \, \dot \beta &=& - \frac{1}{I^{xx}} \frac{\partial \delta U}{\partial \alpha},
\nonumber \\
\ddot \beta - \frac{1}{2} \sin{2\beta} \, \dot \alpha^2 + \frac{I^{zz} \Omega^z }{I^{xx}} \sin\beta \, \dot \alpha &=& - \frac{1}{I^{xx}} \frac{\partial\delta U}{\partial \beta} .
\label{eleq}
\end{eqnarray}

First, we consider an illustrative case where only $\bar s^{ZZ}$ is
nonzero. This is the case where the tensor $\bar s^{ij}$ degenerates
to a vector as in the correspondence \rf{replace}. With $\bar s^{XX} = \bar
s^{YY} = 0$, $\delta U$ is independent of $\alpha$ so that solutions with $\ddot
\alpha = 0, \, \dot \beta = 0$ are consistent with Eqs.~\rf{eleq} given
\begin{eqnarray}
\dot \alpha = \frac{I^{zz} \Omega^z \pm \sqrt{\left( I^{zz} \Omega^z \right)^2 - 8 \bar s^{ZZ} C I^{xx} \cos^2\beta } }{2 I^{xx} \cos\beta } .
\label{sol1}
\end{eqnarray}
When $\bar s^{ZZ} = 0$, the ``$+$'' sign recovers the
free-precessing solution and the ``$-$'' sign gives the stationary spinning
solution. Correspondingly, for $\bar s^{ZZ} \ne 0$, the ``$+$'' sign gives
what we call {\it{twofold precession}}, and the ``$-$'' sign gives what we
call {\it{single precession}}. The single-precession solution, once
approximated at the leading order of $\bar s^{ZZ} \leftrightarrow - \alpha_2
|\boldsymbol{w}|^2$, is the one obtained by \citet{Nordtvedt:1987}.

The solutions \rf{sol1} hold when the matrix $\bar s^{IJ}$ has only one
nonvanishing eigenvalue whose eigenvector is chosen to be the $Z$-axis.
Similar solutions exist when the matrix $\bar s^{IJ}$ has two same
eigenvalues whose eigenvectors are chosen to be the $X$-axis and the
$Y$-axis, because in such case $\delta U$ is also independent of $\alpha$. In
fact, only one nonvanishing eigenvalue is a special case of two same
eigenvalues where the same eigenvalues are zero. In conclusion, when the
matrix $\bar s^{IJ}$ has two same eigenvalues, we can set up the
$X$-$Y$-$Z$ frame properly so that solutions similar to Eq.~\rf{sol1} with
a constant $\beta$ exist.

Now we can discuss the solutions in general when the matrix $\bar s^{IJ}$
has three different eigenvalues. In such case, $\delta U$ depends on both
$\alpha$ and $\beta$. Solutions with a constant $\beta$ do not exist, but as
Lorentz violation is supposed to be small, we use the ansatz
\begin{eqnarray}
\alpha &=& \alpha_0 + a t + \alpha^{(1)} ,
\nonumber \\
\beta &=& \beta_0 + \beta^{(1)},
\label{ansatz}
\end{eqnarray}
to find perturbative solutions at the leading order of Lorentz violation.
In the ansatz \rf{ansatz}, $\alpha_0$ and $\beta_0$ are constants, while
$\alpha^{(1)}$ and $\beta^{(1)}$ are functions of time assumed to be much
smaller than 1. The constant $a$, representing the precessing angular
velocity for \Reply{free rotation}, takes the values
\begin{eqnarray}
a = \begin{cases}
\frac{\omega}{\cos\beta_0}, & {\rm free \ precession} , \\
0 ,  & {\rm stationary \ spin} ,
\end{cases}
\end{eqnarray}
with $\omega \equiv I^{zz} \Omega^z /I^{xx}$ according to Eq.~\rf{sol1} when
$\bar s^{ZZ} = 0$. Substituting derivatives of the ansatz \rf{ansatz} into
Eqs.~\rf{eleq}, and approximating $\alpha$ and $\beta$ in the equations with
$\alpha_0+at$ and $\beta_0$, we find
\begin{eqnarray}
\ddot \alpha^{(1)} \pm \frac{\omega}{\sin\beta_0} \, \dot \beta^{(1)} &=& - \frac{1}{I^{xx} \sin^2\beta_0} \frac{\partial \delta U}{\partial \alpha} \Big|_{\alpha=\alpha_0+at,\, \beta=\beta_0} ,
\nonumber \\
\ddot \beta^{(1)} \mp \omega\sin\beta_0 \, \dot \alpha^{(1)} &=& - \frac{1}{I^{xx} } \frac{\partial \delta U}{\partial \beta} \Big|_{\alpha=\alpha_0+at,\, \beta=\beta_0} .
\label{eleq1}
\end{eqnarray}
Note that the upper signs are for the perturbation to a free-precessing
star, generating twofold-precession solutions, while the lower signs are
for the perturbation to a stationary spinning star, generating
single-precession solutions. We will stick to this sign convention when we
write expressions containing upper and lower signs. Equations \rf{eleq1}
are two coupled oscillation equations for $\dot \alpha^{(1)}$ and $\dot
\beta^{(1)}$ with driven forces. Using Eq. \rf{anisopot} to write out the
right-hand sides of Eqs.~\rf{eleq1}, the solutions can be found as
\begin{eqnarray}
\dot \alpha^{(1)} &=& \mp a^{(1)} \mp \frac{A}{\sin\beta_0} \sin{(\omega t + \varphi)} + \eta_\alpha \, \tilde C \cos{2(at+\alpha_0)} ,
\nonumber \\
\dot \beta^{(1)} &=&  A \cos{(\omega t + \varphi)} + \eta_\beta \, \tilde C \sin{2(at + \alpha_0)} ,
\label{sol2}
\end{eqnarray}
with $A$ and $\varphi$ being two integral constants of the homogeneous
solutions. The constants $a^{(1)}, \, \tilde C, \, \eta_\alpha$, and $\eta_\beta$
are
\begin{eqnarray}
a^{(1)} &=& \frac{C \left( 2 \bar s^{ZZ} - \bar s^{XX} - \bar s^{YY} \right) \cos\beta_0}{I^{zz} \Omega^z} ,
\nonumber \\
\tilde C &=& \frac{C \left( \bar s^{XX} - \bar s^{YY} \right) }{I^{zz} \Omega^z} ,
\nonumber \\
\eta_\alpha &=& \begin{cases}
\frac{\cos\beta_0 \left(2+\cos^2\beta_0 \right)}{4-\cos^2\beta_0} , & {\rm twofold \ precession}, \\
\cos\beta_0 , & {\rm single \ precession},
\end{cases}
\nonumber \\
\eta_\beta &=& \begin{cases}
\frac{3\sin\beta_0 \cos^2\beta_0}{4-\cos^2\beta_0} , & ~~~~{\rm twofold \ precession}, \\
\sin\beta_0 , & ~~~~{\rm single \ precession} .
\end{cases}
\label{consdef}
\end{eqnarray}
We point out that when $\bar s^{XX} = \bar s^{YY} = 0$, the inhomogeneous
solutions in Eqs.~\rf{sol2} recover the solutions \rf{sol1} at the leading
order of Lorentz violation.

The constant $a^{(1)}$, generated by the constant term on the right-hand
side of the second equation in Eqs.~\rf{eleq1}, represents the forced precession
due to Lorentz violation. It can be absorbed into $a$ by redefining $a$ to
be
\begin{eqnarray}
a = \begin{cases}
\frac{\omega}{\cos\beta_0} - a^{(1)}, & {\rm twofold \ precession} , \\
a^{(1)}, & {\rm single \ precession}  .
\end{cases}
\label{precangvel}
\end{eqnarray}
Note that the forced precession acts oppositely on a free-precessing star
and a stationary spinning star.

Now we can use this new $a$ in solutions \rf{sol2}, and the benefit is that
the second term in $\dot \beta^{(1)}$ is no longer constant for the
single-precession solution. As $\beta$ at most changes from $0$ to $\pi$ by
definition, its rate of change really should not contain any constant term.
Keeping in mind that the definition \rf{precangvel} is used, then the
ansatz \rf{ansatz}, together with the solutions for $\alpha^{(1)}$ and
$\beta^{(1)}$,
\begin{eqnarray}
\alpha^{(1)} &=& \pm \frac{A}{\omega \sin\beta_0} \cos{(\omega t + \varphi)} + \eta_\alpha \frac{\tilde C}{2a} \sin{2(at+\alpha_0)} ,
\nonumber \\
\beta^{(1)} &=& \frac{A}{\omega}\sin{(\omega t + \varphi)} - \eta_\beta \frac{\tilde C}{2a} \cos{2(at + \alpha_0)} ,
\label{sol3}
\end{eqnarray}
gives the perturbation solutions with integral constants $A , \, \varphi, \,
\alpha_0$, and $\beta_0$.

\begin{figure*}
    \centering
    \includegraphics[width=17cm]{fig_fig2.pdf}
    \caption{Example I of twofold precession (denoted as tp-I hereafter).
    The left plots are the solutions of the Euler angles $\alpha$ and $\beta$ as
    functions of time. The right plot shows the trajectories of the heads
    of $\hat {\boldsymbol{e}}_x$ (green) and $\hat {\boldsymbol{e}}_z$
    (blue) in the $X$-$Y$-$Z$ frame while their tails are fixed at the
    origin. The arrows mark the two vectors at $t=0$. The initial values,
    $\alpha|_{t=0}=\frac{\pi}{2}$, $\beta|_{t=0} \approx 0.701$, $\dot
    \alpha|_{t=0} \approx 1.509$, and $\dot \beta|_{t=0} = 0$ are adopted by
    setting $\alpha_0=\frac{\pi}{2}$, $\beta_0 = 0.7$, $A=0$, and
    $\varphi=\frac{\pi}{2}$ in Eqs.~\rf{integraltoinitial}. Other parameters
    used are $\gamma|_{t=0} = 0$, $\Omega^z = 1$, $I^{zz}/I^{xx} = 1.1$, and
    $\left\{ \bar s^{XX}, \bar s^{YY}, \bar s^{ZZ} \right\} = \{ 0.02,
    0.01, -0.04\}$. Time and time derivatives are dimensionless under the
    time unit $t_c$ defined in Eq.~\rf{unitt}. }
    \label{fig2}
\end{figure*}
\begin{figure*}
    \centering
    \includegraphics[width=17cm]{fig_fig3.pdf}
    \caption{Example II of twofold precession (denoted as tp-II
    hereafter). See the caption of Fig.~\ref{fig2} for the meaning of
    illustration. The initial values, $\alpha|_{t=0}=\frac{\pi}{2}$,
    $\beta|_{t=0} \approx 1.065$, $\dot \alpha|_{t=0} \approx 0.889$, and $\dot
    \beta|_{t=0} = 0$ are adopted by setting $\alpha_0=\frac{\pi}{2}$, $\beta_0 =
    0.7$, $A=0.4$, and $\varphi=\frac{\pi}{2}$ in Eqs.~\rf{integraltoinitial}.
    Other parameters are the same as those in Fig.~\ref{fig2}.}
    \label{fig3}
\end{figure*}
\begin{figure*}
    \centering
    \includegraphics[width=17cm]{fig_fig4.pdf}
    \caption{Example I of single precession. See the caption of
    Fig.~\ref{fig2} for the meaning of illustration. The initial values,
    $\alpha|_{t=0}=\frac{\pi}{2}$, $\beta|_{t=0} \approx 0.671$, $\dot
    \alpha|_{t=0} \approx -0.068$, and $\dot \beta|_{t=0} = 0$ are adopted by
    setting $\alpha_0=\frac{\pi}{2}$, $\beta_0 = 0.7$, $A=0.01$, and
    $\varphi=\frac{\pi}{2}$ in Eqs.~\rf{integraltoinitial}. Other parameters
    are the same as those in Fig.~\ref{fig2}. }
    \label{fig4}
\end{figure*}
\begin{figure*}
    \centering
    \includegraphics[width=17cm]{fig_fig5.pdf}
    \caption{Example II of single precession. See the caption of
    Fig.~\ref{fig2} for the meaning of illustration. The initial values,
    $\alpha|_{t=0}=\frac{\pi}{2}$, $\beta|_{t=0} \approx 1.146$, $\dot
    \alpha|_{t=0} \approx -0.008$, and $\dot \beta|_{t=0} = 0$ are adopted by
    setting $\alpha_0=\frac{\pi}{2}$, $\beta_0 = 1.4$, $A=0.01$, and
    $\varphi=\frac{\pi}{2}$ in Eqs.~\rf{integraltoinitial}. Other parameters
    are the same as those in Fig.~\ref{fig2}. }
    \label{fig5}
\end{figure*}

The perturbation solutions are useful in analytically illustrating how
Lorentz violation described by the matrix $\bar s^{IJ}$ with three
different eigenvalues affects a freely rotating spheroidal star. But a
complete discussion of the rotational motion has to involve numerical
solutions, because there are solutions to Eqs.~\rf{eleq} unable to be
described by the perturbation solutions even at the leading order of
Lorentz violation. To see this, let us consider solving Eqs.~\rf{eleq} with
initial values given by
\begin{eqnarray}
\alpha|_{t=0} &=& \alpha_0 \pm \frac{A}{\omega \sin\beta_0} \cos{ \varphi} + \eta_\alpha \frac{\tilde C}{2a} \sin{2\alpha_0} ,
\nonumber \\
\beta|_{t=0} &=& \beta_0 + \frac{A}{\omega}\sin{\varphi} - \eta_\beta \frac{\tilde C}{2a} \cos{2\alpha_0} ,
\nonumber \\
\dot \alpha|_{t=0} &=& a \mp \frac{A}{\sin\beta_0} \sin{\varphi} + \eta_\alpha \tilde C \cos{2\alpha_0} ,
\nonumber \\
\dot \beta|_{t=0} &=& A \cos\varphi + \eta_\beta \tilde C \sin{2\alpha_0} .
\label{integraltoinitial}
\end{eqnarray}
The fact that the initial values $\alpha|_{t=0}$, $\beta|_{t=0}$, $\dot
\alpha|_{t=0}$, and $\dot \beta|_{t=0}$ might be arbitrarily assigned indicates
that the corresponding integral constants $A, \, \varphi, \, \alpha_0$, and $\beta_0$
should also be able to take any values. However, the perturbative approach
restricts $A$ to be small, losing solutions where $\dot \alpha|_{t=0}$ is at
the order of $\omega/\cos\beta_0$ but significantly different from it. These
solutions belong to the case of twofold precession, but because the
directions of angular momentum around which \Reply{free precession} happens
deviate from the $Z$-axis too much, \Reply{free precession} around the
$Z$-axis is no longer able to serve as the zeroth-order solution for the
perturbation method. More seriously, the perturbative approach assumes
$\alpha^{(1)}$ and $\beta^{(1)}$ to be much smaller than 1. But for the
single-precession solution in Eqs.~\rf{sol3}, this assumption is very much
questionable as a small quantity $a = a^{(1)}$ appears in the denominators
of the second terms of $\alpha^{(1)}$ and $\beta^{(1)}$. Somehow the
perturbative \Reply{single-precession solution} does mimic numerical solutions in
one period of the forced precession when $\beta_0$ is small ($\beta_0 \mathrel{\rlap{\lower4pt\hbox{\hskip1pt$\sim$}
1$). But when $\beta_0$ gets close to $\pi/2$, the approximation fails
completely.

Figures \ref{fig2}--\ref{fig5} display examples to compare the perturbation
solutions with numerical solutions, and also illustrate the cases where
perturbation method is invalid. In the figures the time and angular
velocities have dimensionless values by employing a time unit $t_c$ defined
as
\begin{eqnarray}
t_c = \sqrt{\frac{I^{xx}}{C}} .
\label{unitt}
\end{eqnarray}
Then, for demonstrating purpose, the integral constant $\Omega^z$ is taken to
be 1, the ratio $I^{zz}/I^{xx}$ is taken to be 1.1, and the eigenvalues of
$\bar s^{ij}$ are taken to be $\bar s^{XX} = 0.02, \, \bar s^{YY}=0.01, \,
\bar s^{ZZ} =-0.04$.

Before turn to observational consequences, let us address the fact that our
discussion above assumes that the forced precession due to Lorentz
violation is around the $Z$-axis. Depending on the initial values, the
forced precession can be around the other two principal axes of the matrix
$\bar s^{IJ}$. An interesting result from our study of the numerical
solutions is that without loss of generality, if the three different
eigenvalues of $\bar s^{ij}$ satisfy $\bar s^{ZZ} < \bar s^{YY} < \bar
s^{XX}$, then \Reply{the forced precession} around the $Y$-axis is unstable. This
is very similar to the Dzhanibekov effect (also called the tennis racket
theorem) in free rotation~\cite{Louis:1851}. Once this is stated, we point
out that the above discussion equally applies to \Reply{the forced precession} around
the $X$-axis by changing the indices
\begin{eqnarray}
Z \rightarrow X, \quad X \rightarrow Y, \quad Y \rightarrow Z ,
\end{eqnarray}
in Eqs.~\rf{consdef} and keeping in mind that $(\alpha, \beta, \gamma)$ now refer
to the Euler angles between the $x$-$y$-$z$ frame and the $Y$-$Z$-$X$ frame
(see Fig.~\ref{euler2}).

\section{observational consequences}
\label{sec3}

\begin{figure}
    \centering
    \includegraphics[width=5.5cm]{fig_eulerangles2}
    \caption{Euler angles $\alpha, \, \beta$, and $\gamma$ connecting the
    $Y$-$Z$-$X$ frame and the $x$-$y$-$z$ frame. The line $PQ$ is the
    intersection of the $Y$-$Z$ plane and the $x$-$y$ plane. }
    \label{euler2}
\end{figure}

NSs observed as pulsars provide tests against the above
predicted single-precession motion due to Lorentz-violating
gravity~\cite{Nordtvedt:1987, Shao:2013wga, Shao:2014oha,
Shao:2019nso}. In Sec.~\ref{sec3a} we will discuss constraints set on the
coefficients for Lorentz violation by applying the single-precession
solution to two observed solitary
millisecond pulsars~\cite{Shao:2013wga}. On the other hand, the
twofold-precession solution, developed from free precession once
Lorentz-violating gravity presents, has not found its tests as evidences of
free-precessing NSs are yet preliminary (see e.g.
Ref.~\cite{Stairs:2000zz}). \Reply{Based on the hypothesis that NSs might possess tiny rigid deformation,} the modulations on pulsar signals and the
emission of continuous GWs due to free precession have been studied in the
literature to predict observational signatures for searching such
NSs~\cite{Zimmermann:1979ip, 2001MNRAS.324..811J, Gao:2020zcd,Gao:2020dwy}.
Following similar considerations, we investigate in Secs.~\ref{sec3b}
and~\ref{sec3c} observational signatures in pulsar pulses and continuous GW
signals from NSs undergoing the twofold-precession motion.

\subsection{Single-precession motion and constraints on Lorentz violation}
\label{sec3a}

We start with considering the change of the angle between the spin axis
$\hat {\boldsymbol{e}}_z$ and a fixed direction in the $X$-$Y$-$Z$ frame when a star
is in the state of single precession. Calling that angle $\lambda$ and
describing the fixed direction using the spherical angular coordinates
$(\theta_o, \, \phi_o)$ in the $X$-$Y$-$Z$ frame, we have the relation
\begin{eqnarray}
\cos\lambda &=& \sin\theta_o \sin\beta \sin{(\alpha - \phi_o)} + \cos\theta_o \cos\beta ,
\end{eqnarray}
where $\alpha$ and $\beta$ are Euler angles used in Sec.~\ref{sec2}. Assume that
the changes in $\alpha$ and $\beta$ are much smaller than 1 during a certain
time interval, then the change of $\cos\lambda$ can be approximated as
\begin{eqnarray}
\Delta \cos\lambda &\approx& \left( \sin\theta_o \cos\beta \sin{(\alpha - \phi_o)} - \cos\theta_o \sin\beta \right) \Delta \beta
\nonumber \\
&& + \sin\theta_o \sin\beta \cos{(\alpha - \phi_o)} \Delta \alpha.
\end{eqnarray}
For estimations, we use the perturbation solution \rf{ansatz} to
approximate the changes in $\alpha$ and $\beta$ as
\begin{eqnarray}
\Delta \alpha &\approx&\dot \alpha \Delta t = \left( 1 + \frac{A}{a^{(1)} \sin\beta_0} \sin{\varphi} + \frac{\cos\beta_0 \tilde C}{a^{(1)}} \cos{2\alpha_0} \right) a^{(1)} \Delta t,
\nonumber \\
\Delta \beta &\approx&\dot \beta \Delta t = \left( \frac{A}{a^{(1)}} \cos{\varphi} + \frac{\sin\beta_0 \tilde C}{a^{(1)}} \sin{2\alpha_0} \right) a^{(1)} \Delta t ,
\end{eqnarray}
where the start of the time interval $\Delta t$ has been set at $t=0$.
Because $A,\, \tilde C$, and $a^{(1)}$ are at the same order, $\Delta \alpha, \,
\Delta\beta$, and hence $\Delta \cos\lambda$ are proportional to $a^{(1)} \Delta t$ with
factors of order unity. Here we neglect any situation where $(\theta_o,\,
\phi_o)$ and $\{A,\, \varphi, \, \alpha_0, \, \beta_0\}$ are fine-tuned to vanish the
factor for $\Delta \cos\lambda$. Once an observation puts a bound on the
change of $\lambda$ during a certain time interval, corresponding constraints
on the combination $2\bar s^{ZZ} - \bar s^{XX} - \bar s^{YY}$ is
\begin{eqnarray}
\left| 2\bar s^{ZZ} - \bar s^{XX} - \bar s^{YY} \right| \lesssim  \left| \frac{I^{zz} \Omega^{z} }{C \Delta t \cos\beta_0} \, \Delta\cos\lambda \right| .
\label{est}
\end{eqnarray}

The estimation \rf{est} depends on $\beta_0$, which roughly characterizes the
angle between the symmetric axis of the star and the Lorentz-violating
principal axis around which the forced precession happens. Checked with
numerical results, we find that when $\beta_0 \lesssim 1.1 $, the estimation
\rf{est} gives correct orders of magnitude for $\left| 2\bar s^{ZZ} - \bar
s^{XX} - \bar s^{YY} \right|$. But as we mentioned in Sec.~\ref{sec2b},
when $\beta_0$ approaches $\pi/2$, the estimation \rf{est} fails because the
perturbation method breaks down. A semianalytical relation based on
Eq.~\rf{est} and numerical results is
\begin{eqnarray}
\left| 2\bar s^{ZZ} - \bar s^{XX} - \bar s^{YY} \right| \lesssim \eta \left| \frac{I^{zz} \Omega^{z} }{C \Delta t} \, \Delta\cos\lambda \right| ,
\label{est1}
\end{eqnarray}
where $\eta$ can be approximated by $1/\cos\beta_0$ for $\beta_0 \lesssim 1.1 $
but then only increases to $10$ when $\beta_0$ approaches $\pi/2$.

\begin{table*}
    \caption{Observational quantities from PSRs~B1937+21 and
    J1744$-$1134~\cite{Shao:2013wga} to constrain Lorentz violation. The
    pulse profile of PSR~B1937+21 consists of a main pulse and an
    interpulse, both of which put bounds on the change of $\lambda$ and hence
    constrain Lorentz violation. The pulse width is taken at $50\%$ intense
    level in practice. The last row, conservative bound on $\left| 2\bar s^{ZZ} - \bar s^{XX} - \bar s^{YY} \right|$, is obtained by setting $\eta = 10$.  \label{tab2}}
    \def1.3{1.3}
    \begin{tabularx}{\textwidth}{p{5.7cm}p{4.5cm}p{4.5cm}p{3cm}}
      \hline   \hline
        & Main pulse of PSR B1937+21 & Interpulse of PSR B1937+21 & PSR J1744$-$1134 \\
      \hline
      Spin period (ms)  & $1.6$ & $1.6$ & $4.1$ \\
      $W$ (deg)         & $8.3$ & $10$ & $13$ \\
      Bound on $\Delta W$ ($10^{-3}$ deg)     & $-48$ & $53$ & $20$ \\
      $\lambda$ (deg)  & $100$ & $100$ & $95$ \\
      $\chi$ (deg) & 75 & $105$ & 51 \\
      Bound on $\Delta \lambda$ ($10^{-3}$ deg)   & $-19$ & $-5.2$ & 1.6 \\
      Bound on $\left| 2\bar s^{ZZ} - \bar s^{XX} - \bar s^{YY} \right|$  & $10^{-14}\eta$ & $10^{-15}\eta$ & $10^{-16}\eta$ \\
      Conservative bound on $\left| 2\bar s^{ZZ} - \bar s^{XX} - \bar s^{YY} \right|$  & $10^{-13}$ & $10^{-14}$ & $10^{-15}$ \\
      \hline
  \end{tabularx}
\end{table*}

Now we are ready to apply Eq.~\rf{est1} to two solitary pulsars: PSRs~B1937+21 and
J1744$-$1134. Their observed pulse profiles over 15 years were thoroughly
analyzed by \citet{Shao:2013wga}. Here we simply use the conclusion that
the change in the angle $\lambda$, now being the angle between the spin axis
and the line of sight, is bounded by
\begin{eqnarray}
\Delta \lambda = \frac{1}{2} \frac{\sin{\frac{W}{2}}}{\cot\lambda\cos{\frac{W}{2}}+\cot\chi} \, \Delta W ,
\label{dela}
\end{eqnarray}
where $W$ is the pulse width, $\Delta W$ is the change in $W$ during the 15
years, and $\chi$ is the angle between the spin axis of the NS and the
symmetric axis of the pulsar beam which has been assumed to take the shape
of a narrow cone~\cite{Lorimer:2005misc}. The values of those quantities
for each pulsar are listed in Table \ref{tab2}. The derived bounds on $\Delta
\lambda$ and hence on the combination $\left|
2\bar s^{ZZ} - \bar s^{XX} - \bar s^{YY} \right|$ are shown in the last three
rows of the table (the last row is the same as the last second row except that $\eta$ is set to $10$, which is suggested as the upper limit of $\eta$ by our numerical solutions). As an estimation, a uniform density of $\rho_{\rm NS} \sim 10^{15} {\rm g}/{\rm cm}^3$ and a fluid deformation of $\epsilon_f \sim 10^{-3}$ have been used to get
\begin{eqnarray}
\frac{I^{zz} }{C} \approx \frac{15}{4\pi} \frac{1}{ \rho_{\rm NS} \epsilon_f} \sim 10^{-5} \, {\rm s}^2,
\end{eqnarray}
for the NSs.

Under the correspondence \rf{replace}, the constraints in Table~\ref{tab2}
are consistent with those in Ref.~\cite{Shao:2013wga} for the PPN $\alpha_2$
parameter. Then we do notice that they are 3 to 5 orders of magnitude
better than the {\it global} results in Ref.~\cite{Shao:2014oha} where the
same two solitary pulsars were used but the test was done together with
another 11 binary pulsars in order to obtain global constraints on the
coefficients for Lorentz violation instead of the ``maximal-reach'' ones as
is done here~\cite{Tasson:2019kuw}. This shows that, (i) observations of
solitary pulsars are more sensitive to Lorentz violation than those of
binary pulsars, and (ii) the correlations between different coefficients
for Lorentz violation can severely degrade the
constraints~\cite{Shao:2014oha}.

\subsection{Twofold-precession motion and pulsar signal}
\label{sec3b}

To construct pulsar signals from a NS under the twofold-precession motion,
we adopt the cone model to describe the radiation beam
\cite{1984A&A...132..312G, Lorimer:2005misc}. Figure \ref{figcone}
illustrates a half radiation cone in the $X$-$Y$-$Z$ frame. In principle, the radiation comes from two opposite sides of a NS and is in the shape of a doule cone. But for clarity, we will only track the half cone as shown in Fig.~\ref{figcone} and analyze the signal from it. The signal from the opposite half cone can be obtained simply by reversing the unit vector along the axis of the cone in the following analysis.

In the cone model, signals are observed when the line of sight
is inside the radiation cone. Mathematically, it means
\begin{eqnarray}
\cos\rho < \hat {\boldsymbol{m}} \cdot \hat {\boldsymbol{n}} ,
\label{ineq}
\end{eqnarray}
where $\rho$ is the semi-open angle of the cone, $\hat {\boldsymbol{m}}$ is
the unit vector along the axis of the cone, and $\hat {\boldsymbol{n}}$ is
the unit vector along the line of sight. The axis of the cone, generally
believed to be aligned with the magnetic dipole moment of the star, is
fixed in the body frame and therefore can be described by two constant
angular coordinates $(\chi, \delta)$ as
\begin{eqnarray}
\hat {\boldsymbol{m}} = \sin\chi \cos\delta \, \hat {\boldsymbol{e}}_{x} + \sin\chi \sin\delta \, \hat {\boldsymbol{e}}_{y} + \cos\chi \, \hat {\boldsymbol{e}}_{z} .
\end{eqnarray}
The direction of the line of sight, if we neglect the proper motion of the
star relative to the observer, is fixed in the $X$-$Y$-$Z$ frame and will be
described by two constant angular coordinates $(\theta_o, \phi_o)$ as
\begin{eqnarray}
\hat {\boldsymbol{n}} = \sin\theta_o \cos\phi_o \, \hat {\boldsymbol{e}}_{X} + \sin\theta_o \sin\phi_o \, \hat {\boldsymbol{e}}_{Y} + \cos\theta_o \, \hat {\boldsymbol{e}}_{Z} .
\end{eqnarray}
Using Eqs.~\rf{trans}, the product of $ \hat {\boldsymbol{m}}$ and $ \hat
{\boldsymbol{n}}$ can be written as a function of time due to the
time-dependent Euler angles,
\begin{eqnarray}
\hat {\boldsymbol{m}} \cdot \hat {\boldsymbol{n}} &=&  \sin\theta_o \sin\chi \cos{(\alpha-\phi_o)} \cos{(\gamma+\delta)}
\nonumber \\
&& - \sin\theta_o \sin\chi \sin{(\alpha-\phi_o)} \cos\beta\sin{(\gamma+\delta)}
\nonumber \\
&& + \cos\theta_o \cos\chi \cos\beta
\nonumber \\
&& + \sin\theta_o \cos\chi \sin{(\alpha-\phi_o)} \sin\beta
\nonumber \\
&&  + \cos\theta_o \sin\chi \sin\beta\sin{(\gamma+\delta)} .
\label{pulsemn}
\end{eqnarray}
The angles $\phi_o$ and $\delta$ can be absorbed into the initial values of
$\alpha$ and $\gamma$ in the above expression. For given $\theta_o, \, \chi$, and
$\rho$, once the rotation of the star is known, the inequality \rf{ineq} can
be solved to predict the time intervals during which signals are observed.

\begin{figure}
    \centering
    \includegraphics[width=8.5cm]{fig_radcone}
    \caption{Cone model for pulsar radiation~\cite{1984A&A...132..312G,
    Lorimer:2005misc}. The unit vector along the magnetic dipole moment is
    denoted as $\hat {\boldsymbol{m}}$. It is the symmetric axis of the
    cone. The semi-open angle of the cone is denoted as $\rho$. While the
    star rotates, $\hat {\boldsymbol{m}}$ changes its orientation, and at
    some time, the cone encloses a unit vector $\hat {\boldsymbol{n}}$, so
    the radiation is received by the observer in the direction of $\hat
    {\boldsymbol{n}}$. The colatitude and the azimuth for $\hat
    {\boldsymbol{m}}$ and $\hat {\boldsymbol{n}}$ in the $X$-$Y$-$Z$ frame
    are denoted as $(\Theta, \, \Phi)$ and $(\theta_o, \, \phi_o)$ respectively. }
    \label{figcone}
\end{figure}

A quick review of the observational signatures of free-precessing pulsars
is heuristic to explore the twofold-precession case; more details can be
found in e.g. Ref.~\cite{Gao:2020zcd}. For an axisymmetric NS, the
free-precessing solution has
\begin{eqnarray}
\dot \alpha &=& \frac{I^{zz} \Omega^z}{I^{xx} \cos\beta},
\nonumber \\
\dot \gamma &=& \Omega^z - \dot\alpha \cos\beta = - \epsilon \Omega^z ,
\label{angfreqfre}
\end{eqnarray}
with the angle $\beta$ being constant. Therefore, the expression $\hat
{\boldsymbol{m}} \cdot \hat {\boldsymbol{n}}$ as shown in Eq.~\rf{pulsemn}
is simply a sum of four sinusoids of angular frequencies $\dot\alpha-\dot\gamma,
\, \dot\alpha, \, \dot\alpha+\dot\gamma$, and $\dot\gamma$. The analysis is further
simplified by the fact that the deformation $\epsilon$ is extremely small for
NSs so that the frequency components $\dot\alpha\pm\dot\gamma$ are basically the
same as $\dot\alpha$, leaving the pulsar period to be $P \approx 2\pi/ \dot
\alpha$ and modulated with a period of $P_{P} \approx 2\pi/|\dot\gamma|$. Here we
assume that the directions of the $z$-axis and the $Z$-axis are properly
chosen so that $\Omega^z$ and $\cos\beta$ are positive.

An exact analytical expression of the pulsar period $P$ in terms of the
Euler angles was derived in Ref.~\cite{2001MNRAS.324..811J} by calculating
the rate of the azimuthal angle of $\hat {\boldsymbol{m}}$ in the
$X$-$Y$-$Z$ frame. Denote the angular coordinates of $\hat
{\boldsymbol{m}}$ in the $X$-$Y$-$Z$ frame as $(\Theta, \Phi)$, then their
relation with $(\chi, \delta)$ can be found by transforming the Cartesian
components of $\hat {\boldsymbol{m}}$ under the rotation \rf{rmat}. The
result is
\begin{eqnarray}
\cos\Theta &=& \cos\beta \cos\chi + \sin\beta \sin\chi \sin{(\gamma+\delta)} ,
\nonumber \\
\tan\Phi &=& \tan{(\alpha + \Delta\Phi)} ,
\label{tranang}
\end{eqnarray}
where
\begin{eqnarray}
\tan\Delta\Phi = \frac{\cos\beta \sin{(\gamma+\delta)} - \sin\beta \cot\chi}{\cos{(\gamma+\delta)}} .
\end{eqnarray}
The pulsar period, defined as the time intervals between $\hat
{\boldsymbol{m}}$ consecutively passing through the plane formed by the
$Z$-axis and the observer, is then
\begin{eqnarray}
P \equiv \frac{2\pi}{\dot \Phi} = \frac{2\pi}{\dot\alpha + \Delta\dot\Phi} \approx \frac{2\pi}{\dot\alpha} \left( 1 - \frac{\Delta\dot\Phi}{\dot\alpha} \right) ,
\label{pulsep}
\end{eqnarray}
where
\begin{eqnarray}
\Delta\dot\Phi = - \dot\gamma \frac{\sin\beta\cot\chi\sin{(\gamma+\delta)} - \cos\beta}{\cos^2{(\gamma+\delta)} + \left( \sin\beta \cot\chi - \cos\beta\sin{(\gamma+\delta)} \right)^2} \,.
\end{eqnarray}
Equation \rf{pulsep} not only confirms that the pulsar period $P$ is
approximately generated by the angular frequency $\dot\alpha$ and that the
modulation on $P$ involves sinusoidal functions of $\gamma$, but also
describes the time evolution of $P$ quantitatively, ready to fit
observational data, say, to extract the parameters $\dot \alpha$, $\dot \gamma$,
$\beta$, and $\chi$.

Besides the pulsar period, the widths of pulse signals also encode
information on the motion of the star. With the cone model of the radiation
beam, it is defined as the change of the azimuthal angle of $\hat
{\boldsymbol{m}}$ in the $X$-$Y$-$Z$ frame during the time when $\hat
{\boldsymbol{n}}$ is inside the cone, namely
\begin{eqnarray}
W \equiv \Phi_2 - \Phi_1 ,
\end{eqnarray}
where $\Phi_1$ and $\Phi_2$ are the azimuthal angle $\Phi$ solved from the
equation of the inequality \rf{ineq} with
\begin{eqnarray}
\hat {\boldsymbol{m}} \cdot \hat {\boldsymbol{n}} = \sin\Theta\sin\theta_0\cos{(\Phi-\phi_o)} + \cos\Theta\cos\theta_o .
\label{mninertial}
\end{eqnarray}
Fortunately, because $|\dot\gamma|$ is much smaller than $\dot\Phi \approx
\dot\alpha$, the angle $\Theta$, given by the first equation in Eqs.~\rf{tranang}, can
be treated as unchanged during $\Phi$ changing from $\Phi_1$ to $\Phi_2$ so
that the roots $\Phi_1$ and $\Phi_2$ are symmetric about the azimuthal angle
$\phi_o$ of the observer, namely,
\begin{eqnarray}
\Phi_1 \approx \phi_o - \frac{W}{2}, \quad \Phi_2 \approx \phi_o + \frac{W}{2},
\end{eqnarray}
and hence, as given in Refs.~\cite{1984A&A...132..312G, Lorimer:2005misc},
setting the expression~\rf{mninertial} equal to $\cos\rho$ leads to an
analytical expression for $W$,
\begin{eqnarray}
\cos{\frac{W}{2}} \approx \frac{\cos\rho - \cos\Theta \cos\theta_o}{\sin\Theta \sin\theta_o} .
\label{pulsewid}
\end{eqnarray}
When the time dependence of $\Theta$ is put back in Eq.~\rf{pulsewid}, it
describes the modulated pulse width at the period $P_P \approx 2\pi/|\dot
\gamma|$. By fitting this equation to observational data, the parameters
$\dot\gamma, \, \rho$, and $\theta_o$ can be extracted.

\begin{figure*}
    \centering
    \includegraphics[width=16cm]{fig_mnfig1}
    \caption{Time evolutions of $\hat {\boldsymbol{m}} \cdot \hat
    {\boldsymbol{n}}$ using the numerical solutions in Figs.~\ref{fig2} and \ref{fig3}. The constant angles $\theta_o$ and $\chi$ are taken to be $0.8\,
    {\rm rad}$ and $0.5\, {\rm rad}$ in the plots. The red dashed lines
    represent the cosine value of $\rho = 0.6\, {\rm rad}$, illustrating the
    time segments when $\hat {\boldsymbol{m}} \cdot \hat {\boldsymbol{n}}$
    is above $\cos\rho$ so that the signal is observed under the
    corresponding setup. \Reply{With our choice of $\epsilon_r = 0.1$ in the numerical solutions, the time unit is $t_c \sim 10^{-3}\,$s if the density of the NS is $\rho_{\rm NS} \sim 10^{15} \, {\rm g}/{\rm cm}^3$.} }
    \label{mnfig1}
\end{figure*}

Let us now turn to the Lorentz-violating twofold-precession case. Similar
characteristics of pulse signals exist. In fact, Eqs.~\rf{pulsep} and
\rf{pulsewid} apply straightforwardly as long as $|\dot\gamma| \ll \dot\alpha$.
This is exactly the regime where the perturbation solution of the twofold
precession works. Neglecting the oscillating terms, $\alpha^{(1)}$ and
$\beta^{(1)}$, the perturbation solution has
\begin{eqnarray}
\dot \alpha &\approx& a = \frac{I^{zz} \Omega^z}{I^{xx}\cos\beta_0} - a^{(1)} ,
\nonumber \\
\dot\gamma &=& \Omega^z - \dot\alpha \cos\beta \approx -\epsilon \Omega^z + \cos\beta_0 \, a^{(1)} ,
\label{tpconstant}
\end{eqnarray}
where the free-precessing angular frequencies, $\dot \alpha$ and $\dot \gamma$,
are shifted due to Lorentz violation. By fitting Eqs.~\rf{pulsep} and
\rf{pulsewid} to observational data, we will be able to extract the two shifted
frequencies, $\dot\alpha$ and $\dot\gamma$, and the angular parameters $\beta_0$,
$\chi$, $\rho$, and $\theta_o$.

\begin{figure}
    \centering
    \includegraphics[width=8cm]{fig_mnfft.pdf}
    \caption{Fourier transformations for $\hat {\boldsymbol{m}} \cdot \hat
    {\boldsymbol{n}}$ in Fig.~\ref{mnfig1}. The dashed red vertical lines mark $\dot \alpha \approx 1.438$
    and $|\dot \gamma| = 0.1$, which are the values from the free-precessing
    equations \rf{angfreqfre} with $I^{zz}/I^{xx} = 1.1$, $\Omega^z = 1$, and
    $\beta =0.7$. \Reply{As we have taken $\epsilon_r = 0.1$, the frequency unit is $1/t_c \sim 1000 \, {\rm rad/s}$ if the density of the NS is $\rho_{\rm NS} \sim 10^{15} \, {\rm g}/{\rm cm}^3$.} }
    \label{mnfig2}
\end{figure}

\begin{figure*}
    \centering
    \includegraphics[width=15cm]{fig_mnfig2}
    \caption{Discrete time series of pulsar period and time segment
    associated to pulse width, solved from $\hat {\boldsymbol{m}} \cdot
    \hat {\boldsymbol{n}} = \cos\rho$ with $\hat {\boldsymbol{m}} \cdot \hat
    {\boldsymbol{n}}$ from Fig.~\ref{mnfig1} and $\rho = 0.6\, {\rm rad}$.
    The left plots are for the solution tp-I where the lines are generated
    by Eqs.~\rf{pulsep} and \rf{pulsewid}. The right plots are for the
    solution tp-II where Eqs.~\rf{pulsep} and \rf{pulsewid} fail to
    approximate the discrete points and the lines are spline interpolations
    of the discrete points.}
    \label{mnfig3}
\end{figure*}

A crucial defect here is that once $\dot\alpha$ and $\dot\gamma$ are extracted,
we are unable to deduce $\Omega^z$, $\epsilon$, and $a^{(1)}$ simultaneously. One
more piece of information on the three quantities is required from the
observational data. It involves details of the pattern of $\hat
{\boldsymbol{m}} \cdot \hat {\boldsymbol{n}}$ that are beyond the
perturbation approximation of $\dot\alpha$ and $\dot\gamma$ in
Eqs.~\rf{tpconstant}, and inevitably invokes numerical solutions of the
twofold-precession motion. Numerical calculations for $\hat {\boldsymbol{m}}
\cdot \hat {\boldsymbol{n}}$ are also necessary for another important
reason. When $\dot\gamma$ is not much smaller than $\dot\alpha$, the perturbation
solution breaks down, and to make matters worse, Eqs.~\rf{pulsep} and
\rf{pulsewid} are also not valid any more, so there are no analytical
templates to fit the pulsar period and the pulse width. In this case, the
inequality~\rf{ineq} needs to be solved numerically to deduce discrete time
series of pulsar period and pulse width to fit observational data. The
rotation of the star, the parameters $\epsilon$, $\chi$, $\rho$, and $\theta_o$, and
the coefficients for Lorentz violation $\bar s^{XX}$, $\bar s^{YY}$, and $\bar
s^{ZZ}$, might be determined via elaborate numerical calculations and
careful treatments. Of course, other measurements, like the polarization
properties of the radiation, might provide extra information and can be
combined to derive a full solution. This lays beyond the scope of the current
work.

In Fig.~\ref{mnfig1}, the time evolutions of $\hat {\boldsymbol{m}}
\cdot \hat {\boldsymbol{n}}$ using the two numerical solutions in
Fig.~\ref{fig2} (solution tp-I) and Fig.~\ref{fig3} (solution tp-II) are displayed, while
Fig.~\ref{mnfig2} depicts the corresponding angular frequency spectra. A
decisive characteristic shows up in the case of the solution that
invalidates the perturbation approach (namely, solution tp-II). That is the
extra smallest frequency component appearing in the lower plot of
Fig.~\ref{mnfig2}. This frequency corresponds to the average angular
velocity at which the angular momentum of the star precesses about the
$Z$-axis. Though no simple approximate expression relates it to the
coefficients for Lorentz violation, numerical calculations indicate that
this angular frequency has the same order of magnitude as the coefficients
for Lorentz violation under the dimensionless parametrization using the
time unit defined in Eq.~\rf{unitt}, \Reply{which for realistic NSs has an estimation}
\begin{eqnarray}
t_c \approx \sqrt{\frac{15}{4\pi} \frac{1}{\rho_{\rm NS} \epsilon_r}} \sim 1 {\rm s},
\end{eqnarray}
with density $\rho_{\rm NS} \sim 10^{15}\, {\rm g}/{\rm cm}^3$ and rigid
deformation $\epsilon_r \sim 10^{-8}$ assumed.\footnote{\Reply{We must point out that because we use $I^{zz}/I^{xx} = 1.1$ for illustrative purposes in our numerical examples, meaning $\epsilon_r = 0.1$, the time unit $t_c$ in the plots is really at the order of $10^{-3}\,$s if the density is $\rho_{\rm NS} \sim 10^{15} \,{\rm g}/{\rm cm}^3$.}} The extra frequency component does not
show up in the case of solution tp-I. Because in that case, the angular
momentum of the star is almost aligned with the $Z$-axis, so the precession
of the angular momentum around the $Z$-axis caused by Lorentz violation
degenerates with the free precession, generating a total precession rate
$\dot\alpha$ approximately given by the first equation in Eqs.~\rf{tpconstant}
where the Lorentz-violating contribution $a^{(1)}$ is unable to be decoupled.

In Fig.~\ref{mnfig3}, discrete time series of pulsar period and of time
segment corresponding to pulse width calculated by solving the inequality
\rf{ineq} using the numerical templates of $\hat {\boldsymbol{m}} \cdot
\hat {\boldsymbol{n}}$ from Fig.~\ref{mnfig1} are plotted. For the case of
the solution tp-I (left panels of Fig.~\ref{mnfig3}), the smooth curves
that approximate the discrete points are generated using Eqs.~\rf{pulsep}
and \rf{pulsewid}. They fit very well. For the case of the solution tp-II
(right panels of Fig.~\ref{mnfig3}), Eqs.~\rf{pulsep} and \rf{pulsewid}
fail to approximate the discrete series, so we have to interpolate the
points to reveal the trends of them. Note that in both cases there is a period about $40\, t_c$ originated from the angular frequency between $0.1\, t_c^{-1}$
and $0.2\,t_c^{-1}$ in Fig.~\ref{mnfig2}, while in the right panels a
period about $180\, t_c$ shows up---more perceivable in the pulsar-period plot than in the pulse-width plot---reflecting the characteristic angular
frequency between $0.03 \, t_c^{-1}$ and $0.04 \, t_c^{-1}$ in the lower
panel of Fig.~\ref{mnfig2}.

\subsection{Twofold-precession motion and continuous GWs}
\label{sec3c}

The quadrupole GW radiated by a freely precessing rigid body was calculated
in Refs.~\cite{Zimmermann:1979ip, Zimmermann:1980ba}. To investigate the
quadrupole radiation by a rigid body undergoing the twofold-precession
motion due to Lorentz violation, we generalize Zimmermann's
calculation~\cite{Zimmermann:1980ba} to any rotating rigid body with
torques.

\begin{figure*}
    \centering
    \includegraphics[width=17cm]{fig_gw.pdf}
    \caption{Continuous GWs from NSs subjected to the solutions
    tp-I and tp-II. The angular coordinates of the
    observer are taken to be $\theta_o = 0.8\, {\rm rad}$ and $\phi_o = 0$ in
    the plots. Time and distance are parametrized using the time unit $t_c$
    defined in Eq.~\rf{unitt}, \Reply{which is at the order of $10^{-3}\,$s with our choice of $\epsilon_r = 0.1$ in the numerical solutions and $\rho_{\rm NS} \sim 10^{15} \, {\rm g}/{\rm cm}^3$ assumed for the density of the NS.} }
    \label{gwfig1}
\end{figure*}

The metric perturbation in the $X$-$Y$-$Z$ frame is given by the quadrupole
formula
\begin{eqnarray}
h^{IJ} = - \frac{2}{r} \ddot I^{IJ} ,
\end{eqnarray}
where $r$ is the distance from the star to the observer, and the double
dots denote the second time derivative. The time dependence of the
inertial-frame components of the moment of inertia tensor originates in the
rotation matrix via
\begin{eqnarray}
I^{IJ} = R^{iI} R^{jJ} I^{ij},
\end{eqnarray}
while the body-frame components of the moment of inertia tensor are
constant. The derivatives of $R^{iI}$ can be derived by noticing that the
body-frame axes undergo an infinitesimal rotation
\begin{eqnarray}
\hat {\boldsymbol{e}}_i \rightarrow \hat {\boldsymbol{e}}_i + dt \, \boldsymbol{\Omega} \times \hat {\boldsymbol{e}}_i ,
\end{eqnarray}
at the instant when the body rotates with an angular velocity
$\boldsymbol{\Omega}$. The change of $R^{iI}$ is therefore $dt \, \left(
\boldsymbol{\Omega} \times \hat {\boldsymbol{e}}_i \right) \cdot \hat
{\boldsymbol{e}}_I $, giving
\begin{eqnarray}
\frac{ dR^{iI} }{dt} = \left( \boldsymbol{\Omega} \times \hat {\boldsymbol{e}}_i \right) \cdot \hat {\boldsymbol{e}}_I = \epsilon^{IJK} \Omega^J R^{iK} ,
\label{1stdrmat}
\end{eqnarray}
where $\epsilon^{IJK}$ is the Levi-Civita symbol. Equation \rf{1stdrmat} can be
used to calculate $\dot I^{IJ}$. When combined with the equation of motion
\begin{eqnarray}
\dot I^{IJ} \Omega^J + I^{IJ} \dot\Omega^J = \Gamma^I ,
\end{eqnarray}
$\boldsymbol{\dot\Omega}$ can be solved in terms of $\boldsymbol{\Omega}$ and the
torque $\boldsymbol{\Gamma}$. Then, the demanded second derivative $\ddot
I^{IJ}$ can be found to take the form
\begin{eqnarray}
\ddot I^{IJ} = R^{iI} R^{jJ} A^{ij} ,
\end{eqnarray}
with the help of the derivative of Eq.~\rf{1stdrmat}. The body-frame tensor
$A^{ij}$ is lengthy in general, but much simplified when the body frame
diagonalizes the moment of inertia tensor, which is the case in our setup. The components then read
\begin{eqnarray}
A^{xx} &=& 2\left( \Delta_2 \left(\Omega^y\right)^2 - \Delta_3 \left(\Omega^z\right)^2 \right) ,
\nonumber \\
A^{yy} &=& 2\left( \Delta_3 \left(\Omega^z\right)^2 - \Delta_1 \left(\Omega^x\right)^2 \right) ,
\nonumber \\
A^{zz} &=& 2\left( \Delta_1 \left(\Omega^x\right)^2 - \Delta_2 \left(\Omega^y\right)^2 \right) ,
\nonumber \\
A^{xy} &=& \left( \frac{ \left(\Delta_3 \right)^2}{I^{zz}} + \Delta_1 - \Delta_2 \right) \Omega^x \Omega^y + \frac{\Delta_3}{I^{zz}} \Gamma^z,
\nonumber \\
A^{xz} &=& \left( \frac{ \left(\Delta_2 \right)^2}{I^{yy}} + \Delta_3 - \Delta_1 \right) \Omega^x \Omega^z + \frac{\Delta_2}{I^{yy}} \Gamma^y ,
\nonumber \\
A^{yz} &=& \left( \frac{ \left(\Delta_1 \right)^2}{I^{xx}} + \Delta_2 - \Delta_3 \right) \Omega^y \Omega^z + \frac{\Delta_1}{I^{xx}} \Gamma^x ,
\label{gwa}
\end{eqnarray}
where $( \Omega^x, \, \Omega^y, \, \Omega^z)$ are the body-frame velocity components
given in terms of the Euler angles and their derivatives in
Eqs.~\rf{angvel1}, and the body-frame torque components $(\Gamma^x, \, \Gamma^y,
\, \Gamma^z)$ can be calculated from the anisotropic potential $\delta U$ via
\begin{eqnarray}
\Gamma^x &=& -\frac{\sin\gamma}{\sin\beta} \, \partial_\alpha \delta U - \cos\gamma \, \partial_\beta \delta U + \cot\beta\sin\gamma \, \partial_\gamma \delta U,
\nonumber \\
\Gamma^y &=& -\frac{\cos\gamma}{\sin\beta} \, \partial_\alpha \delta U + \sin\gamma \, \partial_\beta \delta U + \cot\beta\cos\gamma \, \partial_\gamma \delta U,
\nonumber \\
\Gamma^z &=& - \partial_\gamma \delta U .
\end{eqnarray}
The symbols $\Delta_1,\, \Delta_2$, and $\Delta_3$ follow the definition of
\citet{Zimmermann:1980ba},
\begin{eqnarray}
\Delta_1 = I^{yy}-I^{zz}, \quad \Delta_2 = I^{zz}-I^{xx}, \quad \Delta_3 = I^{xx}-I^{yy}.
\end{eqnarray}
It is plain that Eqs.~\rf{gwa} go back to the result of
\citet{Zimmermann:1980ba} when $\boldsymbol{\Gamma} = 0$.

\begin{figure}
    \centering
    \includegraphics[width=8cm]{fig_gwfft.pdf}
    \caption{Fourier transformations for the continuous GWs shown in
    Fig.~\ref{gwfig1}. The dashed red vertical lines mark $\dot \alpha \approx
    1.438$ and $2\dot \alpha \approx 2.876 $, which are the angular
    frequencies of the continuous GWs from a free-precessing NS with
    $I^{zz}/I^{xx} = 1.1, \, \Omega^z = 1$, and $\beta =0.7$. \Reply{As we have taken $\epsilon_r = 0.1$, the frequency unit is $1/t_c \sim 1000 \, {\rm rad/s}$ if the density of the NS is $\rho_{\rm NS} \sim 10^{15} \, {\rm g}/{\rm cm}^3$.} }
    \label{gwfig2}
\end{figure}

As we restrict ourselves to spheroids, the tensor $A^{ij}$ is further
simplified by $I^{xx} = I^{yy}$ and $\Gamma^z = 0$. In the absence of Lorentz
violation, the star precesses freely, and the tensor $\ddot I^{IJ}$ has a
simple form,
\begin{eqnarray}
\ddot I^{IJ} &=& \frac{ 1}{2}\Delta_2 \, \dot\alpha^2 \sin{2\beta}
\nonumber \\
&& \times \begin{pmatrix}
2\cos{2\alpha} \tan\beta & 2\sin{2\alpha} \tan\beta & -\sin\alpha \\
2\sin{2\alpha} \tan\beta & -2\cos{2\alpha} \tan\beta & \cos\alpha \\
-\sin\alpha & \cos\alpha & 0
\end{pmatrix} ,
\label{gwwf}
\end{eqnarray}
which tells that $h^{IJ}$ is a spherical wave with exactly two
frequency components, namely $\dot\alpha$ and $2\dot\alpha$. When Lorentz
violation presents and the star is subjected to the twofold-precession
motion, the tensor $\ddot I^{IJ}$ becomes complicated as $\dot\beta$ no
longer vanishes. But if we only consider motions given by the perturbation
solution, we expect Eq.~\rf{gwwf} to be a fair approximation and that the
continuous GW still has two frequency components, $\dot\alpha$ and $2\dot\alpha$,
with $\dot\alpha$ being the Lorentz-violating shifted angular frequency given
by the first equation in Eqs.~\rf{tpconstant}.

For an observer with colatitude $\theta_o$ and azimuth $\phi_o$ in the
$X$-$Y$-$Z$ frame, the
basis to decompose the two physical degrees of freedom for the continuous
GW can be taken as
\begin{eqnarray}
\boldsymbol{e}_{+} \equiv \hat {\boldsymbol{\theta}}_o \otimes \hat {\boldsymbol{\theta}}_o  - \hat {\boldsymbol{\phi}}_o \otimes \hat {\boldsymbol{\phi}}_o , \quad \boldsymbol{e}_{\times} \equiv \hat {\boldsymbol{\theta}}_o \otimes \hat {\boldsymbol{\phi}}_o  + \hat {\boldsymbol{\phi}}_o \otimes \hat {\boldsymbol{\theta}}_o ,
\end{eqnarray}
where
\begin{eqnarray}
\hat {\boldsymbol{\theta}}_o &=& \cos\theta_o \cos\phi_o \, \hat{\boldsymbol{e}}_X + \cos\theta_o \sin\phi_o \, \hat{\boldsymbol{e}}_Y - \sin\theta_o \hat{\boldsymbol{e}}_Z,
\nonumber \\
\hat {\boldsymbol{\phi}}_o &=& -\sin\phi_o \, \hat{\boldsymbol{e}}_X + \cos\phi_o \, \hat{\boldsymbol{e}}_Y  .
\label{thphinertial}
\end{eqnarray}
Then the ``+'' and the ``$\times$'' components of the continuous GW
are
\begin{eqnarray}
h_{+} &=& \frac{1}{2} \left( \hat {\theta}_o^I \, \hat \theta_o^J  - \hat {\phi}_o^I \, \hat {\phi}_o^J \right) h^{IJ} = - \frac{1}{r} \left( \hat {\theta}_o^I \, \hat \theta_o^J  - \hat {\phi}_o^I \, \hat {\phi}_o^J \right) \ddot I^{IJ},
\nonumber \\
h_{\times} &=& \frac{1}{2} \left( \hat {\theta}_o^I \, \hat \phi_o^J  + \hat {\phi}_o^I \, \hat {\theta}_o^J \right) h^{IJ} = - \frac{1}{r} \left( \hat {\theta}_o^I \, \hat \phi_o^J  + \hat {\phi}_o^I \, \hat {\theta}_o^J \right) \ddot I^{IJ},
\label{hpluscross}
\end{eqnarray}
where $\hat {\theta}_o^I$ and $\hat {\phi}_o^I$ are the components of $\hat
{\boldsymbol{\theta}}_o$ and $\hat {\boldsymbol{\phi}}_o$ in the $X$-$Y$-$Z$
frame as shown in Eqs.~\rf{thphinertial}.

When continuous GWs from NSs are detected, Eqs.~\rf{hpluscross} can be used
as template waveforms to match observational data. Similar to the analysis
of pulsar signal, to extract the coefficients for Lorentz violation, $\bar
s^{XX}$, $\bar s^{YY}$, $\bar s^{ZZ}$, and the NS parameters, $\epsilon$,
$\Omega^z$, $\theta_o$, $r$, from observational data, elaborate numerical
calculations are required. Using the two numerical solutions in
Figs.~\ref{fig2} and \ref{fig3} to calculate $\ddot I^{IJ}$, we plot the
corresponding waveforms in Fig.~\ref{gwfig1} and their Fourier
transformations in Fig.~\ref{gwfig2}. Tiny high frequency components in the
spectra are found. They are featured in the subplots of Fig.~\ref{gwfig2}
and turn out to be exactly three times of the fundamental frequencies,
which are represented by the peaks near $1.438 \,t_c^{-1}$ and can be
regarded as the average values of $\dot \alpha$ for the twofold-precession
motions. Frequencies higher than the third harmonic exist, but their
amplitudes are even much smaller. We point out that unlike the pulsar
signal of solution tp-II, there is no low frequency at the order of the
coefficients for Lorentz violation in the continuous GW spectra. However, the
advantage of continuous GW signal is that both solution tp-I and solution
tp-II exhibit the third harmonic, characterizing the deviation of the waveform
from that of a free-precessing NS, though in the case of solution tp-I, the
amplitude of the third harmonic is further suppressed. These Lorentz-violating
features are worth looking for and may reveal new physics beyond the
current understanding.

\section{conclusion}
\label{sec4}

Lorentz violation modifies the Newtonian potential by an anisotropic
correction, generating a torque on spheroidal stars that forces the
otherwise conserved angular momentum of the star to precess around a
preferred direction defined by the coefficients for Lorentz violation. To solve
the motion rigorously, we first calculate the anisotropic gravitational
self-energy of the star caused by Lorentz violation in the minimal
gravitational SME in Sec.~\ref{sec2a}. The result is proved to be
equivalent to that of \citet{Nordtvedt:1987} when the star is treated as a
stationary spinning fluid star in equilibrium so the tensor virial relation
holds. Discrepancy occurs when the star possesses a rigid deformation which
invalidates the tensor virial relation. Besides stationary spin, free
precession is also a solution of motion for such stars in the absence of
torques.

Then in Sec.~\ref{sec2b}, \Reply{the forced precession} caused by the
Lorentz-violating torque on stationary spinning stars and on
free-precessing stars is explicitly calculated using perturbation method
to solve the rotational equations of motion. Interestingly, we find that
the direction of the forced precession on stationary spinning stars is
opposite to the direction of the forced precession on free-precessing stars
as shown in Eq.~\rf{precangvel}. Numerical solutions are explored to check
the validity of the perturbation approach. Initial values for which the
perturbation approach fails are identified. The study of the solutions is
finished by clarifying that the preferred directions around which \Reply{the
forced precession} happens are exactly the eigenvectors of the matrix $\bar
s^{ij}$. An interesting result is that when the matrix $\bar s^{ij}$
acquires 3 different eigenvalues, the forced precession is unstable if it
is around the eigenvector corresponding to the middle eigenvalue. This is
similar to the well-known Dzhanibekov effect (or the tennis racket
theorem).

After the solutions of motion are studied thoroughly, we apply them to
explore observational consequences for NSs in Sec.~\ref{sec3}.
\Reply{Section \ref{sec3a} treats the two solitary pulsars in Ref.~\cite{Shao:2013wga} as stationary
spinning spheroids at zeroth order, and sets bounds on the coefficients for Lorentz violation by attributing any possible tiny alteration of the NSs'
orientations hidden in observational uncertainties to the forced precession
due to Lorentz violation.} With the connection
between the SME coefficients $\bar s^{ij}$ and the PPN coefficient $\alpha_2$
shown in Eq.~\rf{replace}, the constraints obtained here are consistent
with the ones for the $\alpha_2$ coefficient in Ref.~\cite{Shao:2013wga}. However, our ``maximal-reach'' constraints obtained from the
two solitary pulsars are 3 to 5 orders of magnitude better than those in
Ref.~\cite{Shao:2014oha}. The exactly same two solitary pulsars are used
there but together with another 11 binary pulsars for a global analysis.
Notice that here we are using a different coordinate frame from that of
Ref.~\cite{Shao:2014oha} so a plain comparison is only heuristic.
Nevertheless, this suggests that the orbital motions of the binary pulsars
are less sensitive to Lorentz violation than the rotational motions of the
solitary pulsars. Therefore, we urge observers to analyze the stability of
pulsar spins---possibly via the pulse profile stability as it was done by
\citet{Shao:2013wga}---with more suitable systems, and put tighter bounds
on the coefficients for Lorentz violation.

In Sec.~\ref{sec3b} and Sec.~\ref{sec3c}, pulsar signals and continuous GWs
from Lorentz-violating-affected free-precessing NSs are investigated. When
the angular momentum of the star is close to the preferred direction around
which the forced precession happens, the spectra of the signals are very
much like those from a free-precessing NS. The forced precession does
shift frequencies in the spectra, but the contribution from it is
practically unable to be decoupled from the free-precessing frequency
components. When the angular momentum of the star makes a relatively large
angle to the preferred direction around which the forced precession
happens, decisive characteristics that are worth
searching for show up in the spectra. For pulsar signals, the signature is an extremely low frequency
component as shown in the lower plot of Fig.~\ref{mnfig2}. The observables,
namely the pulsar period and the pulse width, are not only modulated by the
average rate of the Euler angle $\gamma$, but also modulated by this frequency
(right plots in Fig.~\ref{mnfig3}). The polarization characteristics of
pulsar pulses are interesting to investigate in future studies, and they could provide more information on the rotation of the NS and its radiation properties. For GWs,
the signature is the third harmonic shown in the subplots of
Fig.~\ref{gwfig2}. Free-precessing spheroidal NSs emit continuous GWs only
at the first and the second harmonics of the rate of the Euler angle $\alpha$.
But once affected by Lorentz violation, the continuous GW emitted by the
star is no longer a simple sum of sinusoidal functions. The waveform
generally involves all harmonics of the average rate of the Euler angle
$\alpha$, with the amplitudes decrease rapidly after the first two of them.

The characteristic frequencies in the spectra of the signals are
qualitative supports for Lorentz violation if observed. To extract
quantitative values of the coefficients for Lorentz violation, as well as
the physical parameters of the NS, numerical calculations are necessary to
fit observational data. Tentative candidates of free-precessing NSs are
proposed in pulsar observations \cite{Stairs:2000zz, Shabanova:2001ud},
while searches of continuous GWs have not yet confidently identified any
positive detections~\cite{Pisarski:2019vxw, Covas:2020nwy,
Dergachev:2020fli, Papa:2020vfz, Steltner:2020hfd, Zhang:2020rph}. Our study supplies necessary templates for potential new tests of Lorentz violation that take
advantage of the two state-of-the-art observation channels in the era of
multimessenger astronomy. Once free-precessing NSs are unambiguously
identified, the new tests using their modulated pulsar signals and
continuous GWs are bound to enrich our fundamental knowledge on Lorentz
symmetry of our spacetime.

\begin{acknowledgments}
We thank Quentin G.\ Bailey and V.\ Alan Kosteleck\'y for discussions, \Reply{and anonymous referees for comments}. This work was supported
by \Reply{the National SKA Program of China (2020SKA0120300),} the National Natural Science Foundation of China (11975027, 11991053,
11721303), the Young Elite Scientists Sponsorship Program by the China
Association for Science and Technology (2018QNRC001), the Max Planck
Partner Group Program funded by the Max Planck Society, and the
High-Performance Computing Platform of Peking University. It was partially
supported by the Strategic Priority Research Program of the Chinese Academy
of Sciences through the Grant No. XDB23010200. R.X. is supported by the
Boya Postdoctoral Fellowship at Peking University.
\end{acknowledgments}