text:
\section{Introduction}

Current astrophysical observations suggest that General Relativity (GR) may have to be modified on large scales. Dark energy, dark matter and even inflation have been suggested as natural consequences of certain modified gravity theories (see e.g. Ref.~\cite{fR}). Perhaps similar modifications will be necessary in the non-linear, dynamical regime, the so-called {\emph{strong-field}}, such as when black holes (BHs) merge and neutron stars (NSs) collide. In this regime, physical phenomena is not well-described by a leading-order weak-field and slow-motion expansion of the Einstein equations. Instead, one must either retain high-orders in such perturbative expansion, or solve the full set of Einstein equations numerically. The only way to determine whether modifications to GR in the strong field are necessary is to make observations and test GR in this regime.

Modified gravity theories have been tested accurately in the Solar System and with binary pulsar observations~\cite{TEGP,will-living,stairs}. While the former allow for tests in the weak-field regime only, the latter has led to tests when gravitational fields are moderately strong. Binary pulsar observations, however, are still not capable of probing the non-linear regime of GR. For example, the orbital velocity of the double binary pulsar J0737-3039~\cite{burgay,lyne,kramer-double-pulsar} is roughly $10^{-3}$ smaller than the speed of light, and thus, its orbital behavior can be well-approximated by a weak-field, slow-velocity expansion of the Einstein equations, retaining only the leading and first subleading terms. On the other hand, GR will soon be tested during compact binary late coalescence through gravitational wave (GW) observations~\cite{TEGP,will-living,gair-living}, which will allow for a probe of the fully non-linear and dynamical, strong-field regime.

Recently, there has been a suggestion that binary pulsar observations may constrain deviations from GR to such a degree that GW tests will not be necessary in the future~\cite{damour-GW-vs-pulsar}. Such a suggestion emerged from studies of certain scalar-tensor theories~\cite{fujii}, the constraints on which indeed cannot be improved with second-generation GW detectors, such as Advanced LIGO~\cite{eardley,will1977,zaglauer,will1994,arun-dipole}. But of course, this suggestion depends strongly on the particular modified gravity model considered. For example, the same type of GW detectors may be able to place constraints on massive graviton propagation that are three orders of magnitude stronger than current binary pulsar constraints~\cite{will1998,arunwill,keppel,delpozzo,cornishsampson,mirshekari,sutton}.

By studying NSs in dynamical Chern-Simons (CS) gravity~\cite{CSreview}, we here find another counterexample to this suggestion. Currently, the only constraints on this theory come from Solar System
observations~\cite{alihaimoud-chen}, by comparing
the CS correction to gravitomagnetic precession to observations with
Gravity Probe B, and from table top~\cite{kent-CSBH} experiments,
by requiring that no CS corrections are present above the smallest gravitational
length scales sampled experimentally on Earth. These tests lead to
comparable constraints, namely $\sqrt{\alpha} < 10^8 \; {\rm{km}}$,
where $\alpha$ is the dimensional CS coupling constant. Such an
incredibly weak
constraint is perhaps not surprising, given the weakness of the
gravitational field in the Solar System. Recent work has shown that future GW observations of BH binary inspirals could improve upon Solar System constraints by as much as seven orders of magnitude~\cite{sopuerta-yunes-DCS-EMRI,pani-DCS-EMRI,canizares,kent-CSGW}. In this paper, we study whether current isolated and binary pulsar observations can already constrain dynamical CS gravity. We will see that indeed this is not possible and that only future GW observations of compact binary coalescence can place stringent constraints on this theory.

Dynamical CS modified gravity is a parity-violating, quadratic-curvature theory that is defined by modifying the Einstein-Hilbert action with a term that is the product of a dynamical scalar field and the Pontryagin density (contraction of the Riemann tensor and its dual). The scalar field has dynamics through the addition of a standard kinetic term and a potential to the action. The Pontryagin term in the action is involved in anomaly cancellation~\cite{polchinski2}. Such a term also appears naturally in heterotic superstring theory~\cite{polchinski1,polchinski2,alexandergates} and loop quantum gravity~\cite{taveras,calcagni,Mercuri:2009zt}. Dynamical CS gravity also arises naturally in effective field theories of inflation~\cite{weinberg-CS}. Historically, CS gravity was introduced without scalar field dynamics, assuming the field was given {\emph{a priori}}~\cite{jackiw}. Such a formulation was plagued with problems, such as metric instability~\cite{Alexander:2007kv} and overconstrained field equations. The modern (dynamical) incarnation of the theory restores the dynamics of the scalar field through a standard kinetic term in the action~\cite{Smith:2007jm}, while treating the model as an effective theory, thus avoiding instabilities~\cite{yunespretorius}.

Dynamical CS gravity has only recently begun to be studied in detail. Non-rotating BHs are described by the Schwarzschild metric, but rotating ones are not. Analytic slowly-rotating BH solutions have been constructed to linear order in spin~\cite{yunespretorius,konnoBH} and to quadratic order in spin~\cite{kent-CSBH} within the small-coupling approximation, i.e. linearizing all expressions in the CS coupling. Slowly-rotating NS solutions were first constructed in~\cite{yunes-CSNS} to linear order in spin and within the small-coupling approximation. This work was extended in~\cite{alihaimoud-chen}, where NS solutions were obtained still to linear order in spin, but without imposing the small-coupling approximation.

The purpose of this paper is to investigate whether any meaningful constraints on dynamical CS gravity can be obtained from isolated~\cite{1.97NS} and binary pulsar observations~\cite{stairs}. In order to achieve this goal, one must first study isolated NSs solutions to quadratic-order in spin, since the latter will introduce modifications to the NS mass as measured by an observer at spatial infinity. One must then place such NS solutions in binary systems and study the conservative and dissipative corrections to the evolution of the binary. The former come from the CS deformation of the quadrupole moment of each individual NS, as well as from the scalar dipole-dipole interaction of the binary components. The latter are caused by the modification to the emitted gravitational and scalar radiation. These corrections can be calculated once one finds the scalar dipole charge induced by each individual NS and its associated CS quadrupole moment deformation. The charge can be extracted by studying the asymptotic behavior of the scalar field at spatial infinity to linear order in spin. The quadrupole moment deformation appears at quadratic-order in spin. All throughout this paper, we treat dynamical CS gravity as an
effective field theory. Such a treatment is definitely valid for the
systems considered here, as explained in detail in~\cite{kent-CSBH}. In turn, this implies that we can use the weak-coupling
approximation to linearize all expression in the CS coupling constant. Such a treatment guarantees that the field
equations are second order, and the theory is ghost-free~\cite{kent-CSBH}.

\subsection{Executive Summary}

The first half of this paper focuses on finding isolated NS solutions to quadratic order in spin in the small-coupling approximation. We follow Hartle's approach~\cite{hartle1967}, in which one treats rotation perturbatively, i.e.~one expands in powers of the product of the NS mass and the NS rotational angular frequency. As expected, we find that CS corrections at quadratic order in spin modify the NS mass monopole and quadrupole moments, where the former can be absorbed by a redefinition of mass. At linear order in spin, the CS corrections appear at much higher multipole order than quadrupole, since the correction to the current dipole moment can be absorbed by a redefinition of spin. Therefore, CS corrections at quadratic order in spin can be larger than those at linear order if NSs are spinning moderately fast, yet sufficiently slowly for the small-rotation expansion to hold.

\begin{figure}[b]
\begin{center}
\includegraphics[width=8.5cm,clip=true]{Q-M-APR-Dup.eps}
\caption{\label{fig:mu-M}
(Color online)
Numerical (symbols) and fitted (curve) results as functions of NS mass with various EoSs.
We plot the dimensionless scalar dipole susceptibility $\bar{\mu}$ (top left), the
quadrupole correction $Q$ (top right), the mass shift $\delta M$ (bottom
left), and the angular momentum shift $\delta J$ (bottom right), which are
defined in Eqs.~\eqref{eq:mu-mubar-relation},
\eqref{eq:QA}, and \eqref{eq:asymptotic-deltaM-deltaJ} respectively.
The y-axis in the top-right panel is normalized by $J^2/M_*$, the absolute value of the NS quadrupole moment
in the point particle limit in GR, while the y-axis of the bottom-left panel is normalized by $\Omega_\mathrm{1ms}$,
the angular frequency of a NS with a period of 1ms defined by ($\Omega_\mathrm{1ms} \equiv 2 \pi / \mathrm{1ms}$).
Moreover, $Q$, $\delta M$ and $\delta J$ are also re-scaled to $\zeta = 0.1$. Observe how the mass shift is always
negative, while the quadrupole moment deformation is always positive.}
\end{center}
\end{figure}
Figure~\ref{fig:mu-M} shows the scalar dipole charge $\bar{\mu}$ (top left) and the CS corrections to the mass $\delta M$ (bottom left), spin angular momentum $\delta J$ (bottom right) and quadrupole moment $Q$ (top right) as a function of the GR mass parameter $M_{*}$ in solar masses $M_{\odot}$. In this figure, $\zeta$ is the dimensionless coupling constant of dynamical CS gravity\footnote{This quantity is proportional to the fourth power of the Chern-Simons natural length scale and it is such that as $\zeta \to 0$ one recovers GR. See Eq.~\eqref{zeta-def} for a more precise definition.}, $\Omega_*$ is the NS angular velocity, $\Omega_\mathrm{1ms}$ is the angular frequency of a NS with a period of 1ms, $J$ is the NS spin angular momentum in GR and we used the APR~\cite{APR}, SLy~\cite{SLy,shibata-fitting}, Lattimer-Swesty (LS220)~\cite{LS, ott-EOS} and Shen~\cite{Shen1,Shen2,ott-EOS} equations of state (EoSs). Observe that the CS corrections \emph{reduce} the observed mass, but they \emph{increase} the observed angular momentum and quadrupole moment with increasing $M_{*}$.

Such corrected NS observables can then be contrasted with NS observations to constrain dynamical CS gravity. Observations of the massive millisecond pulsar J1614-2230~\cite{1.97NS} require that one be able to construct NSs with observed masses larger than 1.93$M_\odot$. Unfortunately, but perhaps not surprisingly, the allowed maximum NS mass in dynamical CS gravity depends not only on the EoS but also on the magnitude of the dimensionless CS coupling parameter $\zeta \propto \xi_{{\mbox{\tiny CS}}} M_{*}^{2}/\mathcal{R}_*^{6}$. Figure~\ref{fig:MR-CS} shows the mass-radius relation for different EoSs with different choices of $\zeta$.
\begin{figure}[tb]
\begin{center}
\includegraphics[width=8.5cm,clip=true]{deltaM-R.eps}
\caption{\label{fig:MR-CS}
(Color online)
Mass-radius relations for GR (solid curves) and CS with different
coupling strengths ($\zeta=0.1,0.5,1$ curves are respectively dotted,
dashed, and dot-dashed). Each panel uses a different EoS:
APR (top left), SLy (top right),  LS220 (bottom left) and Shen
(bottom right). The horizontal dashed line corresponds to the lower bound on the NS mass,
$1.93M_\odot$, provided by observations of PSR J1614-2230. The
mass-radius relations depend on the spin of the NS, which we have set
to match that for PSR J1614-2230. Since APR,  SLy and Shen EoSs can each
produce NSs above the mass bound for all $\zeta < 1$, there is a clear degeneracy between
EoS and the CS correction to the mass-radius relation. Thus, one cannot constrain the theory
from observations of PSR J1614-2230.}
\end{center}
\end{figure}
Therefore, such a test of dynamical CS gravity cannot be performed due to degeneracies with the EoS.

The CS corrections to the NS quadrupole moment and the moment of inertia induce modifications to the orbital evolution of binary pulsars, which then can be contrasted with binary pulsar observations to constrain the modified theory. Observations of the double binary pulsar J0737-3039~\cite{burgay,lyne,kramer-double-pulsar} require that the evolution of certain Keplerian elements~\cite{damour-taylor,TEGP,stairs} agree with GR to within a certain observational uncertainty. The CS modification to the quadrupole moment leads to corrections to these evolution equations, which we explicitly calculate in this paper. We find that the rate of periastron advance $\dot{\omega}$ is the best observable to constrain CS gravity, because its CS modification enters at first PN order relative to the GR behavior, i.e. this modification is suppressed by a single factor of ${\cal{O}}(m/a)$, where $m$ is the total mass of the binary and $a$ is the semi-major axis of the binary. However, for J0737-3039, this still implies a suppression of the CS correction of $\mathcal{O}(10^{-8})$ relative to the measured GR $\dot{\omega}$. Moreover, this CS effect is even smaller than the GR spin-orbit correction to $\dot{\omega}$, which enters at 0.5PN order and depends on the NS's moments of inertia. Therefore, current binary pulsar observations are not accurate enough to allow for tests of dynamical CS gravity.

Unlike for scalar-tensor theories~\cite{damour-GW-vs-pulsar}, our results suggest that dynamical CS gravity cannot be constrained well with binary pulsar observations using current data. Instead, GW observations may be the only way to place stringent constraints. We have therefore found a concrete example of a modified gravity theory that cannot be stringently constrained with current electromagnetic observations.

\subsection{Organization}

The organization of the paper is as follows. In Sec.~\ref{sec:DCS}, we introduce the basics of dynamical CS gravity and explain the approximations that we use throughout the paper. In Sec.~\ref{sec:decomposition}, we explain the coordinate systems that we use and impose an ansatz on the metric and the stress-energy momentum tensor of the matter field. In Sec.~\ref{sec:field-equations}, we derive both GR and CS field equations at each order in spin and derive the exterior solutions, modulo integration constants, which are discussed in Sec.~\ref{sec:exterior}. In Sec.~\ref{sec:interior}, we explain how to construct the interior solution. This amounts to finding the matching conditions at the NS surface and the initial conditions at the center. Then, we explain the numerical procedure that we use to solve these sets of differential equations. In Sec.~\ref{sec:solutions}, we present the numerical results obtained by solving the interior equations and matching the solutions to the exterior ones. We derive the evolution of the NS binary in Sec.~\ref{sec:binary-evolution}, and in Sec.~\ref{sec:applications} we apply the results to the massive millisecond pulsar J1614-2230 and to the double binary pulsar PSR J0737-3039. We conclude in Sec.~\ref{sec:conclusions} with a discussion of possible avenues for future work.

All throughout the paper, we mostly follow the conventions of Misner, Thorne and Wheeler~\cite{MTW}. We use the Greek letters $(\alpha, \beta, \cdots)$ to denote spacetime indices. The metric is denoted by $g_{\mu \nu}$ and it has signature $(-,+,+,+)$. We use geometric units, with $G=c=1$.

\section{Dynamical Chern-Simons Gravity}
\label{sec:DCS}

\subsection{Basics}

The action in dynamical CS gravity is given by~\cite{CSreview}
\ba
S &\equiv & \int d^4x \sqrt{-g} \Big[ \kappa_g R + \frac{\alpha}{4} \vartheta R_{\nu\mu \rho \sigma} {}^* R^{\mu\nu\rho\sigma}  \nonumber \\
& &  - \frac{\beta}{2} \nabla_\mu \vartheta \nabla^{\mu} \vartheta + \mathcal{L}_{{\mbox{\tiny mat}}} \Big]\,.
\label{action}
\ea
Here, $\kappa_g \equiv (16\pi)^{-1}$, $g$ denotes the determinant of the metric $g_{\mu\nu}$ and $R_{\mu\nu \rho \sigma}$ is the Riemann tensor. ${}^* R^{\mu\nu\rho\sigma}$ is the dual of the Riemann tensor which is defined by~\cite{CSreview}
\be
{}^* R^{\mu\nu\rho\sigma} \equiv \frac{1}{2} \varepsilon^{\rho \sigma \alpha \beta} R^{\mu\nu}{}_{\alpha \beta}\,,
\ee
where $\varepsilon^{\rho \sigma \alpha \beta}$ is the Levi-Civita tensor. $\vartheta$ is a scalar field while $\alpha$ and $\beta$ are coupling constants. $\mathcal{L}_{{\mbox{\tiny mat}}}$ denotes the matter Lagrangian density. Following Refs.~\cite{yunespretorius,quadratic,kent-CSBH}, we have neglected the potential of the scalar field for simplicity. We take $\vartheta$ and $\beta$ to be dimensionless while we set $\alpha$ to have a dimension of (length)$^2$~\cite{yunespretorius}.

The field equations in dynamical CS gravity are given by~\cite{CSreview}
\be
G_{\mu\nu} + \frac{\alpha}{\kappa_g} C_{\mu\nu} =\frac{1}{2\kappa_g} (T_{\mu\nu}^\mathrm{mat} + T_{\mu\nu}^\vartheta)\,,
\label{field-eq}
\ee
where $G_{\mu\nu}$ is the Einstein tensor and $T_{\mu\nu}^\mathrm{mat}$ is the matter stress-energy tensor. The C-tensor and the stress-energy tensor for the scalar field are defined by
\begin{align}
C^{\mu\nu} & \equiv  (\nabla_\sigma \vartheta) \epsilon^{\sigma\delta\alpha(\mu} \nabla_\alpha R^{\nu)}{}_\delta + (\nabla_\sigma \nabla_\delta \vartheta) {}^* R^{\delta (\mu\nu) \sigma}\,, \\
\label{eq:Tab-theta}
T_{\mu\nu}^\vartheta & \equiv  \beta (\nabla_\mu \vartheta) (\nabla_\nu \vartheta) -\frac{\beta}{2} g_{\mu\nu} \nabla_\delta \vartheta \nabla^\delta \vartheta\,.
\end{align}

The evolution equation of the scalar field is given by
\be
\square \vartheta = -\frac{\alpha}{4 \beta} R_{\nu\mu \rho \sigma} {}^*R^{\mu\nu\rho\sigma}\,.
\label{scalar-wave-eq}
\ee
By using this equation, one can show that
\ba
\nabla_\nu C^{\mu \nu} &=& -\frac{1}{8} (\nabla^\mu \vartheta) R_{\alpha \beta \rho \sigma} {}^*R^{\beta \alpha\rho\sigma} \nonumber \\
 &=& \frac{1}{2\kappa} \nabla_{\nu} T_\vartheta^{\mu\nu}\,.
\ea
Together with the Bianchi identity, if we take the divergence of Eq.~\eqref{field-eq}, we end up with
\be
\nabla_{\nu} T_{\mbox{\tiny mat}}^{\mu\nu}=0\,.
\label{mat-cons}
\ee
Thus, the equation of motion for a test particle is not modified.

The evolution equation for the scalar field, Eq.~\eqref{scalar-wave-eq}, admits a flat background solution of the form $\vartheta = C_{\mu} x^{\mu}$, for $C_{\mu} = {\rm{const}}$, i.e.~a solution to the homogeneous equation in a flat background. The constants $C_{\mu}$ control the strength of the CS modification, and thus, one would expect them to be proportional to $\alpha/\beta$. Although such a solution is in principle allowed, $C_{\mu}$ would have to be prescribed {\emph{a priori}}, as some form of cosmological term. Moreover, such a term would ruins the $\vartheta$ shift-invariance of the equations of motion, which is a desirable property that allows us to treat $\vartheta$ as massless. Also, such a scalar field would have an infinite energy density, when its stress-energy tensor is integrated over the entire manifold. Finally, dynamical CS gravity with such a homogeneous scalar field is very similar to the non-dynamical theory, and thus, it is already severely constrained by cosmological observations~\cite{Dyda:2012rj}. For the rest of this paper, we will ignore the homogeneous solution for $\vartheta$, and instead, concentrate on dynamically generated scalar fields.

\subsection{Small-Coupling, Slow-Rotation Approximations}
\label{sec:approx}

We here work in the small-coupling and slow-rotation approximations. The former implies that we treat the action given in Eq.~\eqref{action} as defining an effective theory, which requires the CS quadratic term (the second term) to be much smaller than the Einstein-Hilbert term. For isolated NSs, the small-coupling approximation is valid if the following inequality holds:
\be
\label{zeta-def}
\zeta \equiv \frac{\xi_{\mbox{\tiny CS}}  M_*^2}{\mathcal{R}_*^6} \ll 1\,, \quad \xi_{\mbox{\tiny CS}} \equiv \frac{\alpha^2}{\beta \kappa_g}\,,
\ee
where $M_*$ and $\mathcal{R}_*$ are the mass and radius of the NS respectively. $\xi^{1/4}$ corresponds to the characteristic length scale of the theory while $\sqrt{\mathcal{R}_*^3/M_*}$ corresponds to the curvature length scale of the NS. This definition of the dimensionless small coupling constant $\zeta$ is the same as $\zeta^2$ defined in~\cite{alihaimoud-chen} modulo an $\mathcal{O}(1)$ numerical factor. Figure~7 of~\cite{alihaimoud-chen} shows that the small coupling approximation is valid when $\zeta \ll 1$. Solar System experiments require $\xi_{\mbox{\tiny CS}}^{1/4} \leq \mathcal{O}(10^8)$km~\cite{alihaimoud-chen}. In this paper, we work to linear order in $\zeta$.

Dynamical CS gravity must thus have cut-off scale beyond which one cannot treat it as an effective model. Ref.~\cite{kent-CSBH} estimated this scale by calculating the critical length scale at which loop corrections to the second term in Eq.~\eqref{action} due to $n$-point interactions become of order unity. Requiring that this length scale be smaller than the smallest scales probed by table-top experiments led to the requirement that $\sqrt{\alpha} < {\cal{O}}(10^{8} \; {\rm{km}})$. This requirement is of the same order as the current constraints from Solar System experiments. For the NS systems we will consider here, both this requirement and the small coupling condition $\zeta \ll 1$ are satisfied.

We further assume that NSs are slowly rotating (i.e. $|\chi | \ll 1$ where $\chi =J/M^{2}$ is the dimensionless spin parameter, with $J$ the magnitude of the NS spin angular momentum and $M$ is its mass) and consider terms up to quadratic order.

All physical quantities, $A$, can be expanded bivariately as
\be
A = \sum_{m,n} \chi'^m \alpha'^n A_{(m,n)}\,,
\ee
where $\chi'$ and $\alpha'$ are book keeping parameters to label the order of the slow-rotation and the small-coupling approximations, respectively. Notice that $A_{(m,n)} \propto \chi^m \alpha^n$. For the spherically symmetric case, $R_{\nu\mu \rho \sigma} {}^* R^{\mu\nu\rho\sigma} = 0$, and thus, there is no source to the inhomogeneous equation of motion for $\vartheta$. Therefore, for the scalar field, $\vartheta_{(0,n)}=0$ and
\be
\vartheta = \alpha' \chi' \vartheta_{(1,1)} + \mathcal{O}(\alpha' \chi'^3)\,.
\ee
Notice that $\vartheta_{(2m,1)}=0$ due to parity.

\section{Spacetime and Matter Decomposition}
\label{sec:decomposition}

\subsection{Metric}

Following Hartle~\cite{hartle1967}, we start with the metric ansatz given by
\begin{align}
ds^2 &= -e^{\bar{\nu}(r)} \left(1+2 \bar{h}(r,\theta) \right) dt^2 \nonumber \\
&  + e^{\bar{\lambda}(r)} \left( 1+\frac{2\bar{m}(r,\theta)}{r-2\bar{M}(r)} \right) dr^2 \\
&  + r^2 \left( 1+2\bar{k}(r,\theta) \right) \left[ d\theta^2 + \sin^2 \theta \left( d\phi - \bar{\omega}(r,\theta) dt \right)^2 \right]\,, \nonumber
\label{metric-ansatz-rth}
\end{align}
where $\bar{\nu}$ and $\bar{\lambda}$ are $\mathcal{O}(\alpha'^0 \chi'^0)$ quantities that only depend on $r$ while $\bar{h}$, $\bar{k}$, $\bar{m}$ and $\bar{\omega}$ refer to perturbations that depend\footnote{These quantities do not depend on $t$ and $\phi$ because we are searching for stationary and axisymmetric NS solutions.} on both $r$ and $\theta$. $\bar{M}$ is related to $\bar{\lambda}$ by
\be
\bar{M}(r) \equiv \frac{\left( 1-e^{-\bar{\lambda}(r)} \right)r}{2}\,.
\ee
The coordinates $(t,r,\theta,\phi)$ are Hartle-Thorne coordinates. In particular, $(r,\theta)$ denote ordinary polar coordinates.

Since we are treating rotation perturbatively, one must be careful with the choice of polar coordinates~\cite{hartle1967}. A perturbative approach is valid only when the fractional change in a quantity between rotating and non-rotating cases is small. If one were to use $(r,\theta)$ coordinates, this condition could not be met for the density $\rho$ and the pressure $p$ near the surface. This is because the rotation changes the shape of a star, and hence there are points on the NS surface where these quantities vanish in the non-rotating case while they acquire finite values in the rotating case, leading to infinite fractional changes. Instead, we adopt the coordinates $(R,\Theta)$ as proposed by Hartle~\cite{hartle1967}, which are defined via
\be
\rho \left[ r(R,\Theta ), \Theta \right] = \rho (R) = \rho_{(0,0)}(R), \quad \Theta = \theta\,.
\ee
In other words, the new radial coordinate $R$ is chosen such that the density at $\rho \left[ r(R,\Theta ), \Theta \right]$ in the rotating configuration is the same as $\rho_{(0,0)}(R)$ in the non-rotating configuration. By construction, the density and the pressure in the new coordinates contain the non-rotating part only:
\be
\rho(R) = \rho_{(0,0)}(R), \quad p(R)=p_{(0,0)}(R)\,.
\ee
We expand $r(R,\Theta)$ as
\be
r(R,\Theta) = R + \xi(R,\Theta)\,,
\ee
where
\be
\xi(R,\Theta) = \alpha'^2 \chi'^2 \xi_{(2,2)}(R,\Theta) + \mathcal{O}(\alpha'^0 \chi'^2, \alpha'^2 \chi'^4)\,.
\label{Eq:xi}
\ee
Notice that we have neglected $\mathcal{O}(\alpha'^0 \chi'^2)$
quantities (pure GR, quadratic in spin effects) in
Eq.~\eqref{Eq:xi}. Since in this work we are interested in the
detectability of the CS corrections, and the mentioned terms would
simply add linearly but not modify the CS corrections at
$\mathcal{O}(\alpha'^2 \chi'^2)$, we ignore the $\mathcal{O}(\alpha'^0
\chi'^2)$ term henceforth.

After the coordinate transformation, the new metric is given by
\ba
ds^2 &=& - \left[\left( 1+2h+\xi \frac{d \nu}{dR} \right)e^\nu - R^2 \omega^2 \sin^2 \Theta \right]dt^2 \nonumber \\
& & -2 R^2 \omega \sin^2 \Theta dt d\phi + \left[ R^2 (1+2k) + 2 R \xi \right] \sin^2\Theta d\phi^2 \nonumber \\
& & + e^\lambda \left(1 +\frac{2m}{R-2M} + \xi \frac{d\lambda}{dR} + 2 \frac{\partial \xi}{\partial R} \right) dR^2 \nonumber \\
& & + 2e^\lambda \frac{\partial \xi}{\partial \Theta} dRd\Theta + \left[ R^2 (1+2k) + 2R\xi \right]d\Theta^2 \nonumber \\
& & + \mathcal{O}(\alpha'^0 \chi'^2, \alpha'^2 \chi'^4)\,,
\label{metric-ansatz-RTh}
\ea
where
\be
M(R) \equiv \frac{\left( 1-e^{-\lambda (R)} \right)R}{2}\,.
\label{M-lambda}
\ee
Each quantity in the above equation is bivariately expanded as
\begin{align}
\nu (R) &= \nu_{(0,0)} (R)\,, \nonumber \\
\lambda (R) &= \lambda_{(0,0)} (R)\,, \nonumber \\
\omega (R,\Theta) &= \chi' \omega_{(1,0)} (R,\Theta)+\alpha'^2 \chi' \omega_{(1,2)} (R,\Theta) \nonumber \\
& \quad + \mathcal{O}(\alpha'^0 \chi'^3, \alpha'^2 \chi'^3) \,, \nonumber \\
h (R,\Theta) &= \alpha'^2 \chi'^2 h_{(2,2)} (R,\Theta)+ \mathcal{O}(\alpha'^0 \chi'^2, \alpha'^2 \chi'^4)\,, \nonumber \\
m (R,\Theta) &= \alpha'^2 \chi'^2 m_{(2,2)} (R,\Theta)+ \mathcal{O}(\alpha'^0 \chi'^2, \alpha'^2 \chi'^4)\,, \nonumber \\
k (R,\Theta) &= \alpha'^2 \chi'^2 k_{(2,2)} (R,\Theta)+ \mathcal{O}(\alpha'^0 \chi'^2, \alpha'^2 \chi'^4)\,.
\end{align}
Due to parity, the $(t,t)$, $(R,R)$, $(\Theta,\Theta)$, $(\phi,\phi)$ and $(R,\Theta)$ metric components only have terms proportional to even powers of $\chi'$ while $(t,\phi)$ component only contains terms that are odd powers in $\chi'$. Also, $\mathcal{O}(\alpha')$ terms do not appear in the metric because of the structure of the field equation and the fact that $\vartheta$ is proportional to $\alpha$ in the small-coupling approximation.

Notice that the coordinate deformation $\xi(R,\Theta)$ is well defined only inside the star. We take it to be constant outside the star. This means that the exterior metric in $(t,r,\theta,\phi)$ coordinates can be obtained simply by replacing $R \to r$ and $\Theta \to \theta$ in the exterior metric in $(t,R,\Theta,\phi)$ coordinates.

\subsection{Neutron Star Stress-Energy Tensor}

In this paper, we assume that the matter field inside a NS is a perfect fluid and that the NS is rotating uniformly. The stress-energy tensor for the matter field $T_{\mu\nu}^{\mbox{\tiny mat}}$ is given by
\be
T_{\mu\nu}^{\mbox{\tiny mat}} = (\rho + p ) u_\mu u_\nu + p g_{\mu\nu}\,,
\ee
where the four-velocity $u^\mu$ is given by
\be
u^\mu = (u^0, 0,0,\Omega_* u^0)\,.
\ee
$\Omega_*$ is the constant angular velocity of the NS. By using the normalization condition $u_\mu u^\mu = -1$, we obtain the time component of the four-velocity $u^0$ as
\ba
e^{\nu/2} u^0 &=& 1 + \alpha'^2 \chi'^2 \left\{ e^{-\nu} R^2 \omega_{(1,2)} [ \omega_{(1,0)}-\Omega_* ] \sin^2(\Theta) \right. \nonumber \\
& & \left. - h_{(2,2)} - \frac{1}{2} \xi_{(2,2)} \frac{d\nu}{dR} \right\} + \mathcal{O}(\alpha'^0 \chi'^2, \alpha'^2 \chi'^4)\,. \nonumber \\
\ea

As mentioned in Eq.~\eqref{mat-cons}, the stress-energy tensor of this
minimally-coupled fluid is divergence free, which contributes another
equation to the system. In order to close the system, one needs a
relationship between density $\rho$ and pressure $P$, an EoS.
In this paper we only consider one-parameter EoSs,
i.e.~$P=P(\rho)$ without any entropy dependence. The EoSs we
use are APR~\cite{APR}, SLy~\cite{SLy,shibata-fitting},
Lattimer-Swesty with nuclear incompressibility of 220MeV
(LS220)~\cite{LS, ott-EOS} and Shen~\cite{Shen1,Shen2,ott-EOS}.
For the latter two, we use a temperature of 0.01MeV and an electron
fraction of 30$\%$. The
EoS and conservation of matter stress-energy tensor close the
system. We discuss the matter equations of motion order-by-order
below.

\section{Modified Field Equations}
\label{sec:field-equations}

In this section, we derive the modified field equations  expanded in $\chi'$ and $\alpha'$.

\subsection{Zeroth Order in Spin}

First, we consider equations at $\mathcal{O}(\chi'^0)$. As mentioned previously, in the spherically symmetric case there is no CS correction, and hence the field equations reduce to the Einstein equations for a non-rotating star. The only non-vanishing components are the $(t,t)$, $(R,R)$, $(\Theta,\Theta)$ and $(\phi,\phi)$ ones, but  the last two components are linearly dependent. The first two components yield
\ba
\label{tt-zeroth}
\frac{d M}{dR} &=& 4 \pi R^2 \rho\,, \\
\label{RR-zeroth}
\frac{d \nu}{dR} &=& 2\frac{4 \pi R^3 p + M}{R(R-2M)}\,.
\ea
Together with the EoS, we need one more equation to close the system of differential equations. Instead of using the $(\Theta,\Theta)$ or $(\phi,\phi)$ component of the Einstein equations, one can use the $R$ component of the conservation equations of the matter stress-energy tensor [Eq.~\eqref{mat-cons}] which is given by
\be
\frac{d\nu}{dR} = -\frac{2}{\rho + p} \frac{dp}{dR}\,.
\label{R-cons-zeroth}
\ee
By combining Eqs.~\eqref{RR-zeroth} and~\eqref{R-cons-zeroth}, we obtain the Tolman-Oppenheimer-Volkoff (TOV) equation
\be
\frac{dp}{dR} = -\frac{(4\pi R^3 p + M) (\rho + p)}{R(R-2M)}\,.
\label{TOV-zeroth}
\ee

\subsection{First Order in Spin}

\subsubsection{GR}

At $\mathcal{O}(\alpha'^0 \chi')$, the only non-vanishing component of the Einstein equations is the $(t,\phi)$ one. By using Eqs.~\eqref{tt-zeroth}--\eqref{TOV-zeroth}, the associated equation can be simplified to
\ba
&& \frac{\partial^2 \tilde{\omega}_{(1,0)}}{\partial R^2} +4 \frac{1 - \pi R^2 (\rho+p) e^{\lambda}}{R}  \frac{\partial \tilde{\omega}_{(1,0)}}{\partial R} \nonumber \\
& & + \frac{e^\lambda}{R^2} \left( \frac{\partial^2 \tilde{\omega}_{(1,0)}}{\partial \Theta^2} +3 \cot\Theta \frac{\partial \tilde{\omega}_{(1,0)}}{\partial \Theta} \right) \nonumber \\
& &  -16 \pi (\rho+p) e^\lambda \tilde{\omega}_{(1,0)}=0\,,
\label{1st-GR}
\ea
where
\begin{equation}
\label{eq:tildeomega-def}
\tilde{\omega}_{(1,0)} \equiv \Omega_* - \omega_{(1,0)}\,.
\end{equation}
We simplify the equation further by performing a Legendre decomposition:
\begin{equation}
\tilde{\omega}_{(1,0)} (R,\Theta) = \sum_\ell \tilde{\omega}_\ell (R) \left( - \frac{1}{\sin\Theta} \frac{\partial}{\partial \Theta} P_\ell (\cos\Theta) \right)\,,
\label{Legendre-omega}
\end{equation}
where $P_\ell$ is the $\ell$-th Legendre polynomial. Then, Eq.~\eqref{1st-GR} becomes
\ba
&& \frac{d^2 \tilde{\omega}_{\ell}}{d R^2} +4 \frac{1 - \pi R^2 (\rho+p) e^\lambda }{R}  \frac{d \tilde{\omega}_{\ell}}{d R} \nonumber \\
& &- \left[ \frac{(\ell+2) (\ell -1)}{R^2} +16 \pi (\rho+p) \right] e^\lambda \tilde{\omega}_{\ell}=0\,.
\ea
By imposing asymptotic flatness and regularity at the center, one can show that $\tilde{\omega}_{\ell}$ must vanish for all $\ell$ except $\ell = 1$~\cite{hartle1967}. Thus, the above equation reduces to
\ba
\frac{d^2 \tilde{\omega}_1}{d R^2} + 4\frac{1- \pi R^2 (\rho +p) e^{\lambda}}{R}\frac{d \tilde{\omega}_1}{d R}  -16 \pi (\rho +p) e^{\lambda} \tilde{\omega}_1=0\,. \nonumber \\
\label{omega1RR}
\ea

\subsubsection{CS: Scalar Evolution Equation}

Next, we look at the scalar evolution equation at $\mathcal{O}(\alpha' \chi')$. By using Eqs.~\eqref{tt-zeroth}--\eqref{TOV-zeroth}, the equation can be written as
\ba
& &  \frac{\partial^2 \vartheta_{(1,1)}}{\partial R^2} + \frac{1 + e^\lambda \left[ 1 -4\pi R^2 (\rho-p)  \right]}{R} \frac{\partial \vartheta_{(1,1)}}{\partial R} \nonumber \\
& & + \frac{e^\lambda}{R^2} \left( \frac{\partial^2 \vartheta_{(1,1)}}{\partial \Theta^2} + \cot \Theta \frac{\partial \vartheta_{(1,1)}}{\partial \Theta} \right) \nonumber \\
&=&  8\pi \frac{\alpha}{\beta}  \delta \; \; e^{(\lambda - \nu)/2} \left( \sin\Theta \frac{\partial^2 \tilde{\omega}_{(1,0)}}{\partial R \partial \Theta} + 2 \cos\Theta \frac{\partial \tilde{\omega}_{(1,0)}}{\partial R} \right)\,, \nonumber \\
\label{scalar-evol-eq}
\ea
where
\be
\delta \equiv \rho- \frac{M}{(4/3) \pi R^3}
\ee
denotes the density shift from the average. Since we are only interested in stationary and axisymmetric solutions, we consider a $\vartheta$ that only depends on the coordinates $R$ and $\theta$. As in Eq.~\eqref{Legendre-omega}, we perform a Legendre decomposition of $\vartheta$:
\be
\vartheta_{(1,1)} (R,\Theta) = \sum_\ell \vartheta_\ell (R) P_\ell (\cos \Theta)\,,
\ee
and find that the $\ell = 1$ term of Eq.~\eqref{scalar-evol-eq} is the only one with a source term.
The only homogeneous solution that satisfies both asymptotic flatness
and regularity at the center is the trivial $\vartheta_\ell = 0$ for $\ell \neq 1$.
For $\ell =1$, Eq.~\eqref{scalar-evol-eq} becomes
\ba
&& \frac{d^2 \vartheta_1}{d R^2} + \frac{1 + e^\lambda \left[ 1 -4\pi R^2 (\rho-p)  \right]}{R} \frac{d \vartheta_1}{d R}  - 2\frac{e^\lambda}{R^2} \vartheta_1 \nonumber \\
&= & 16\pi \frac{\alpha}{\beta}  \delta e^{(\lambda - \nu)/2} \frac{d \tilde{\omega}_1}{d R}\,.
\label{vartheta1RR}
\ea

\subsubsection{CS: Field Equation}

Next, we consider the modified Einstein equation at $\mathcal{O}(\alpha'^2 \chi')$. As in GR, the only non-vanishing component of the field equations is the $(t,\phi)$ one. Notice that since $\vartheta$ is of $\mathcal{O}(\alpha' \chi')$, there is no $\mathcal{O}(\alpha'^2 \chi')$ contribution from $T_{\mu\nu}^\vartheta$. The $(t,\phi)$ component of the field equations is given by
\begin{equation}
\begin{aligned}
& \frac{\partial^2 \omega_{(1,2)}}{\partial R^2} + 4 \frac{1-\pi R^2 (\rho +p) e^\lambda}{R} \frac{\partial \omega_{(1,2)}}{\partial R}\\
& + \frac{e^\lambda}{R^2} \left( \frac{\partial^2 \omega_{(1,2)}}{\partial \Theta^2} + 3 \cot \Theta \frac{\partial \omega_{(1,2)}}{\partial \Theta} \right) -16 \pi (\rho +p) e^\lambda \omega_{(1,2)} \\
{}={}& \frac{128 \pi^2 \alpha e^{(\nu+\lambda)/2}}{R^3 \sin\Theta} \left[ \delta R \frac{\partial^2 \vartheta_{(1,1)}}{\partial R \partial \Theta} + \left( R \frac{d\rho}{dR} - \delta \right) \frac{\partial \vartheta_{(1,1)}}{\partial \Theta} \right]\,.
\end{aligned}
\label{1st-CS}
\end{equation}
As for $\tilde{\omega}_{(1,0)}$, we decompose $\omega_{(1,2)}$ via
\begin{equation}
\label{eq:w-def}
\omega_{(1,2)} (R,\Theta) = \sum_\ell w_\ell (R) \left( - \frac{1}{\sin\Theta} \frac{\partial}{\partial \Theta} P_\ell (\cos\Theta) \right)\,.
\end{equation}
Again, we find that Eq.~\eqref{1st-CS} with $\ell \neq 1$ becomes a homogeneous equation, leading to $w_\ell = 0$ for $\ell \neq 1$ after imposing asymptotic flatness and regularity at the center. For $\ell =1$, the equation reduces to
\ba
&& \frac{d^2 w_1}{d R^2} + 4 \frac{1-\pi R^2 (\rho +p) e^\lambda}{R} \frac{d w_1}{d R}  -16 \pi (\rho +p) e^\lambda w_1 \nonumber \\
&=& - \frac{128 \pi^2 \alpha e^{(\nu+\lambda)/2}}{R^3} \left[ \delta \; R \frac{d \vartheta_1}{d R} + \left( R \frac{d\rho}{dR} - \delta \right) \vartheta_1 \right]\,. \nonumber \\
\label{w1RR}
\ea
Notice that Eq.~\eqref{omega1RR} corresponds to the homogeneous version of Eq.~\eqref{w1RR}.

\subsection{Second Order in Spin}

In this subsection, we derive the stress-energy conservation and the field equations at quadratic order in spin. As mentioned previously, we do not consider the GR part since $\mathcal{O}(\alpha'^0 \chi'^2)$ terms do not appear in any equations at $\mathcal{O}(\alpha'^2 \chi'^2)$. Following the previous subsection, we perform the Legendre decompositions:
\begin{subequations}
\begin{align}
\label{eq:h-decomp}
h_{(2,2)} (R,\Theta) &= \textstyle\sum_\ell h_\ell (R) P_\ell (\cos \Theta)\,, \\
\label{eq:m-decomp}
m_{(2,2)} (R,\Theta) &= \textstyle\sum_\ell m_\ell (R) P_\ell (\cos \Theta)\,, \\
\label{eq:k-decomp}
k_{(2,2)} (R,\Theta) &= \textstyle\sum_\ell k_\ell (R) P_\ell (\cos \Theta)\,, \\
\label{eq:xi-decomp}
\xi_{(2,2)} (R,\Theta) &= \textstyle\sum_\ell \xi_\ell (R) P_\ell (\cos \Theta)\,.
\end{align}
\end{subequations}
We use the gauge freedom of the theory to set $k_0(R)=0$~\cite{hartle1967}.

\subsubsection{CS: Stress-Energy Conservation}
\label{sec:CS-SET-cons}

The non-vanishing components of the stress-energy conservation equations $\nabla^\mu T_{\mu\nu}^{\mbox{\tiny mat}}=0$ are $\nu = R$ and $\Theta$. First, let us look at the $\Theta$ component. By performing a Legendre decomposition, one finds that $\ell = 2$ is the only non-vanishing mode. By using the previously obtained results, one arrives at the algebraic condition
\be
\xi_2 = -\frac{R(R-2M)}{3 (4\pi R^3 p + M)} (3 h_2 - 2 e^{-\nu} R^2 \tilde{\omega}_1 w_1)
\label{xi2}
\ee
in the interior of the star.

Next, let us look at the $R$ component. The only non-vanishing modes are the $\ell =0$ and $\ell =2$ ones. For the latter, we obtain Eq.~\eqref{xi2} after we integrate the differential equation once with respect to $R$. For the former, again by integrating once, we obtain the algebraic relation
\be
h_0 = -\frac{ 4\pi R^3 p + M}{R(R-2 M)} \xi_0 - \frac{2}{3} e^{-\nu} R^2 \tilde{\omega}_1 w_1 + h_{0c}\,,
\label{xi0}
\ee
valid in the interior of the star, where $h_{0c}$ is an integration constant that corresponds to $h(R=0)$.

\vspace{10mm}

\subsubsection{CS: Field Equations}

Let us now derive the modified field equations at $\mathcal{O}(\alpha'^2 \chi'^2)$. The only non-vanishing modes are the $\ell=0$ and $\ell=2$ ones. Let us first focus on the $\ell =0$ case, where our dependent variables are $h_0$, $m_0$ and $\xi_0$. Eq.~\eqref{xi0} already determines $h_0$, and thus, we are left with 2 degrees of freedom ($m_0$ and $\xi_0$). Following Ref.~\cite{hartle1967}, we use the $(t,t)$ and $(R,R)$ components of the field equations, and Eq.~\eqref{xi0}, to find
\bw
\allowdisplaybreaks[1]
\begin{align}
\label{m0R}
\frac{d m_0}{d R} &= -4\pi R^2 \frac{d\rho}{d R} \xi_0 + \frac{1}{6}R^2 e^{-(\nu + \lambda)} \left( 16 \pi \xi_{\mbox{\tiny CS}} \delta\frac{d \tilde{\omega}_1}{d R} -R^2 \frac{d w_1}{d R} \right) \frac{d \tilde{\omega}_1}{d R} \nonumber \\*
& -\frac{16\pi}{3} e^{-\nu} R^4  (\rho +p) \tilde{\omega}_1 w_1 + \frac{2 \pi}{3}  \beta e^{-\lambda} R^2 \left( \frac{d \vartheta_1}{d R} \right)^2  + \frac{4\pi}{3} \beta \vartheta_1^2 \nonumber \\*
& -\frac{4\pi}{3} \alpha e^{-(\nu + \lambda)/2} \left\{ R\left( 9e^{-\lambda} +1+8 \pi R^2 p \right) \frac{d \vartheta_1}{d R} \frac{d \tilde{\omega}_1}{d R} -2\left[ 2e^{-\lambda} + 3 + 16 \pi R^2 (\rho + p) \right] \vartheta_1 \frac{d \tilde{\omega}_1}{d R} \right. \nonumber \\*
& -\left. 64\pi R^2 (\rho+p) \tilde{\omega}_1 \frac{d \vartheta_1}{d R}  -32 \pi R \left[ R \frac{d\rho}{dR} + (\rho + p) \left[ 2 - (1 + 8 \pi R^2 p) e^{\lambda} \right] \right] \vartheta_1 \tilde{\omega}_1  \right\} \,, \\
\label{h0R}
\frac{d \xi_0}{d R} &= \frac{1}{3[1-(1+8\pi R^2 p)e^\lambda ]} \left\{ \frac{6(1+8 \pi R^2 p)}{R} e^{2\lambda} m_0 + 3 \frac{e^{-\lambda} + 16\pi R^2 p - (1-8 \pi R^2 \rho) (1+8\pi R^2 p)e^{\lambda}}{R}e^{\lambda} \xi_0 \right. \nonumber \\*
& + \left. R^4 e^{-\nu} \frac{d \tilde{\omega}_1}{d R} \frac{d w_1}{d R} +4R^3 e^{-\nu} \left( \frac{d \tilde{\omega}_1}{dR} w_{1} + \tilde{\omega}_1 \frac{dw_1}{dR} \right)
+ 4\left( 3 e^{-\lambda} -1-8\pi R^2 p \right) R^2 e^{-(\nu -\lambda)} \tilde{\omega}_1 w_1 +4\pi \beta R^2 \left( \frac{d \vartheta_1}{d R} \right)^2 \right. \nonumber \\*
& - \left. 8\pi \beta e^{\lambda} \vartheta_1^2  + 8\pi \alpha e^{-(\nu-\lambda)/2} \left[ \left( 3 e^{-\lambda} -1-8 \pi R^2 p  \right) R \left( \frac{d \vartheta_1}{d R} \frac{d \tilde{\omega}_1}{d R} +16 \pi e^{\lambda} (\rho +p) \vartheta_1 \tilde{\omega}_1 \right) \right. \right. \nonumber \\*
& + \left. \left. 2 [8 \pi R^2 (\rho + p)-1] \vartheta_1 \frac{d \tilde{\omega}_1}{d R}   \right] \right\}\,.
\end{align}
\allowdisplaybreaks[0]
Using Eqs.~\eqref{xi0}, \eqref{m0R} and~\eqref{h0R}, one can also obtain
\begin{align}
\label{h0R2}
\frac{d h_0}{d R} &= \frac{8 \pi R^2 p + 1}{R^2} e^{2\lambda} m_0 - \frac{4\pi R^2 (\rho +p)}{R-2M}h_0 + \frac{1}{6} R^3 e^{-\nu} \frac{d \tilde{\omega}_1}{d R} \frac{d w_1}{d R} - \frac{8\pi}{3} \frac{e^{-\nu} R^4 (\rho +p)}{R(R-2M)} \tilde{\omega}_1 w_1 -\frac{4\pi}{3} \beta \frac{e^{\lambda}}{R} \vartheta_1^2 + \frac{2\pi}{3} \beta R \left( \frac{d \vartheta_1}{d R} \right)^2 \nonumber \\
& +\frac{4\pi}{3} \alpha e^{(\lambda - \nu)/2} \left[ \left( 3 e^{-\lambda} -1-8 \pi R^2 p  \right) \left( \frac{d \vartheta_1}{d R} \frac{d \tilde{\omega}_1}{d R} +16 \pi e^{\lambda} (\rho +p) \vartheta_1 \tilde{\omega}_1 \right) +2 \frac{8 \pi R^2 (\rho + p)-1}{R} \vartheta_1 \frac{d \tilde{\omega}_1}{d R}  \right]\,.
\end{align}

We now focus on the $\ell = 2$ mode, where we have 4 unknown functions: $h_2$, $m_2$, $k_2$ and $\xi_2$. Eq.~\eqref{xi2} already gives us $\xi_2$, and thus, we need 3 more equations to close the system. We choose one of these to be the $(\Theta,\Theta)$ component of the field equations minus the $(\phi,\phi)$ component. This yields
\begin{align}
m_2 &= -e^{-\lambda} R h_2 -\frac{8\pi}{3} \beta e^{-\lambda} R \vartheta_1^2 - \frac{e^{-(\nu + \lambda)} R^5}{3} \left( e^{-\lambda} \frac{d w_1}{d R} \frac{d \tilde{\omega}_1}{d R} + 16 \pi (\rho + p) w_1 \tilde{\omega}_1 \right)  - \frac{16\pi}{3} \alpha e^{-(\nu + 3 \lambda)/2} R \left[ e^{-\lambda} R \frac{d \vartheta_1}{d R} \frac{d \tilde{\omega}_1}{d R} \right. \nonumber \\
& - \left. \left[ e^{-\lambda} + 4 \pi R^2 (\rho + p) \right] \vartheta_1 \frac{d \tilde{\omega}_1}{d R} - 8\pi R^2 (\rho + p) \tilde{\omega}_1 \frac{d \vartheta_1}{d R} - 8 \pi \left( R^2 \frac{d\rho}{dR} - e^{\lambda} \left( 4 \pi R^3 p + M \right) (\rho+p) \right) \vartheta_1 \tilde{\omega}_1 \right]\,.
\label{m2}
\end{align}
For the remaining 2 equations, we use the $(R,\Theta)$ and $(R,R)$ components of the field equations. By using Eq.~\eqref{xi2}, they yield
\begin{align}
 \frac{d k_2}{d R} &= -\frac{d h_2}{d R} + \frac{\left( 3 e^{-\lambda}-1-8 \pi R^2 p\right) e^{\lambda}}{2R}h_2 + \frac{\left( e^{-\lambda} + 1 + 8\pi R^2 p \right) e^{2\lambda}}{2R^2}m_2 - \frac{8\pi}{3} \beta \vartheta_1 \frac{d \vartheta_1}{d R} \nonumber \\
& + \frac{16\pi}{3} \alpha e^{-(\nu + \lambda)/2} \left( \frac{d \vartheta_1}{d R} \frac{d \tilde{\omega}_1}{d R} - \frac{\vartheta_1}{R} \frac{d \tilde{\omega}_1}{d R} +8 \pi e^{\lambda} (\rho + p) \vartheta_1 \tilde{\omega}_1 \right)
\label{h2k2}
\intertext{and}
\frac{d h_2}{d R} &= \frac{[3-4 \pi R^2 (\rho +p)]e^{\lambda}}{R} h_2 - \frac{ 1 + (1 + 8 \pi R^2 p) e^{\lambda}}{2} \frac{d k_2}{dR} + 2 \frac{e^\lambda}{R} k_2 + \frac{e^{2 \lambda}(1+8\pi R^2 p)}{R^2} m_2 +\frac{4\pi}{3} \beta R\left( \frac{d \vartheta_1}{dR} \right)^2 \nonumber \\
& + \frac{4\pi}{3} \beta \frac{e^\lambda}{R} \vartheta_1^2  - \frac{e^{-\nu}}{6}  R^3 \frac{d \tilde{\omega}_1}{d R} \frac{d w_1}{d R}  + \frac{8\pi}{3} e^{\lambda-\nu} R^3 (\rho +p) \tilde{\omega}_1 w_1 + \frac{8 \pi}{3} \alpha  e^{(\lambda -\nu )/2 } \left( 3 e^{-\lambda} - 1 - 8 \pi R^2 p \right) \frac{d \vartheta_1}{dR} \frac{d \tilde{\omega}_1}{dR} \nonumber \\
& - \frac{16\pi}{3} \alpha  e^{(\lambda -\nu )/2} \frac{1+4\pi R^2 (\rho +p)}{R} \vartheta_1 \frac{d \tilde{\omega}_1}{dR} - \frac{64\pi^2}{3} \alpha e^{(3 \lambda - \nu )/2} \left( 3 e^{-\lambda} - 1 - 8 \pi R^2 p \right) (\rho + p) \vartheta_1 \tilde{\omega}_1\,.
\label{RR2}
\end{align}
These equations are technically valid only inside the NS, because they
were developed from the matter conservation equations
[Eqs.~\eqref{xi2}-\eqref{xi0}]. However, we have checked that the
field equations outside the star are identical to those above after setting $\rho=0=p$.
\ew

\section{Isolated Neutron Star Solution: \\* Exterior Fields}
\label{sec:exterior}

In this section, we analytically solve the equations derived in the previous section outside the NS, where $\rho = 0 = p$. We use the superscript ``ext'' to refer to exterior quantities. We impose asymptotic flatness at spatial infinity as a boundary condition.

\subsection{Exterior Solutions}

\subsubsection{GR}

At $\mathcal{O}(\alpha'^0 \chi'^0)$, we solve Eqs.~\eqref{tt-zeroth} and~\eqref{RR-zeroth} to obtain
\ba
\left( e^\nu \right)^\mathrm{ext} = \left( e^{-\lambda} \right)^\mathrm{ext} &=& 1-\frac{2M_*}{R} \nonumber \\
&\equiv & f(R)\,,
\label{nu-lambda-ext}
\ea
where $M_* = \mathrm{const.}$ is the mass of the non-rotating NS, which is to be determined by matching $\nu^\mathrm{ext}$ or $\lambda^\mathrm{ext}$ with $\nu^\mathrm{int}$ or $\lambda^\mathrm{int}$ at the surface. At $\mathcal{O}(\alpha'^0 \chi'^1)$, we substitute the above equation in Eq.~\eqref{omega1RR} and obtain
\be
\tilde{\omega}_1^\mathrm{ext} = \Omega_* - \frac{2J}{R^3}\,,
\label{omega-ext}
\ee
where $J$ is an integration constant that corresponds to the spin angular momentum of the NS in GR. We define the mass and the spin angular momentum by expanding the metric about $R = \infty$ and extracting the appropriate coefficients~\cite{alihaimoud-chen}.

\subsubsection{CS: First Order in Spin}

Let us substitute Eqs.~\eqref{nu-lambda-ext} and~\eqref{omega-ext} in
the equation for $\vartheta_{1}$ [Eq.~\eqref{vartheta1RR}]. The solution to this equation is
\begin{align}
\vartheta_1^\mathrm{ext} & = \frac{5}{8} \frac{\alpha}{\beta} \frac{J}{M_*^2} \frac{1}{R^2} \left\{ 1+2\frac{M_*}{R} + \frac{18}{5} \frac{M_*^2}{R^2} \right. \nonumber \\
& + \left. C_\vartheta \frac{R^2}{M_*^2} \left[ 1 + \frac{R}{2M_*} \left( 1- \frac{M_*}{R} \right) \ln f(R) \right] \right\}\,,
\label{vartheta1-ext}
\end{align}
where $C_\vartheta$ is an integration constant that is to be
determined by matching this solution with the interior one at the NS
surface under continuous and differentiable boundary
conditions. Notice that an integration constant has been set to
0 by asymptotic flatness.

The modified field equation at $\mathcal{O}(\alpha'^2 \chi')$ can be solved by substituting $\vartheta_1^\mathrm{ext}$ in Eq.~\eqref{w1RR} to obtain
\begin{align}
w_1^\mathrm{ext} &= \frac{2 J_{\mbox{\tiny CS}}}{R^3} - \frac{5}{8} \frac{\xi_{\mbox{\tiny CS}} J}{M_* R^6} \left\{ 1+ \frac{12}{7} \frac{M_*}{R} + \frac{27}{10}\frac{M_*^2}{R^2} \right. \nonumber \\
& - \left. \frac{5}{32} C_\vartheta \frac{R^5}{M_*^5} \left[ 1+ \frac{M_*}{R} - \frac{54}{5} \frac{M_*^3}{R^3} \right. \right. \nonumber \\
& + \left. \left. \frac{R}{2M_*} \left( 1 - \frac{64}{5} \frac{M_*^3}{R^3} + \frac{48}{5} \frac{M_*^4}{R^4} \right) \ln f(R) \right] \right\}\,.
\label{w1-ext}
\end{align}
Here, $J_{\mbox{\tiny CS}}$ is an integration constant that corresponds to a CS
correction to the spin angular momentum (there is also a correction to
the spin angular momentum from the log
term). Eqs.~\eqref{vartheta1-ext} and~\eqref{w1-ext} agree with those
found in~\cite{alihaimoud-chen}.

\subsubsection{CS: Second Order in Spin}

Next, we find the exterior solutions at $\mathcal{O}(\alpha'^2
\chi'^2)$. For the $\ell =0$ mode, we can solve the exterior version
of Eqs.~\eqref{m0R} and~\eqref{h0R2} to find
\bw
\ba
\label{h0_ext}
h_0^\mathrm{ext} &=&  \frac{5}{768} \frac{\xi_{\mbox{\tiny CS}} J^2}{M_*^3 R^5 f(R)} \left\{ 1 + 100 \frac{M_*}{R} +66 \frac{M_*^2}{R^2} + \frac{684}{7} \frac{M_*^3}{R^3} - 648 \frac{M_*^4}{R^4}    -30 C_\vartheta \frac{R^4}{M_*^4} \left[ 1-2 \frac{M_*}{R} \right. \right. \nonumber\\
& & \left. \left. {} - \frac{4}{3} \frac{M_*^2}{R^2} - \frac{2}{3} \frac{M_*^3}{R^3} + \frac{24}{5} \frac{M_*^4}{R^4}  + \frac{R}{2 M_*} \left( 1-3 \frac{M_*}{R} + \frac{1}{3} \frac{M_*^2}{R^2} + \frac{1}{3} \frac{M_*^3}{R^3} + \frac{34}{5} \frac{M_*^4}{R^4} - \frac{24}{5} \frac{M_*^5}{R^5} \right) \ln f(R) \right] \right. \nonumber \\
& & \left. {}- \frac{15}{2} C_\vartheta^2 \frac{R^4}{M_*^4} \left[1 +\frac{R}{6M_*} \left( 1- \frac{M_*}{R} \right)\ln f(R) - \frac{R^2}{12 M_*^2} f(R) [\ln f(R)]^2  \right] \right\} + \frac{2J J_{\mbox{\tiny CS}}}{R^4 f(R)} -  \frac{M_{\mbox{\tiny CS}}}{R f(R)} \,, \\
\label{m0_ext}
m_0^\mathrm{ext} &=&  -\frac{25}{768} \frac{\xi_{\mbox{\tiny CS}} J^2}{M_*^4 R^3} \left\{ 1 + 3\frac{M_*}{R} + \frac{322}{5} \frac{M_*^2}{R^2} + \frac{198}{5} \frac{M_*^3}{R^3} + \frac{6276}{175} \frac{M_*^4}{R^4} - \frac{17496}{25} \frac{M_*^5}{R^5} \right.  \nonumber \\
& & \left. {}+ 2 C_\vartheta \frac{R^2}{M_*^2} \left[ 1+\frac{M_*}{R} - \frac{4}{5} \frac{M_*^2}{R^2} - \frac{72}{5} \frac{M_*^3}{R^3} + \frac{R}{2M_*} \left( 1 + \frac{1}{5} \frac{M_*^2}{R^2} - \frac{78}{5} \frac{M_*^3}{R^3} +24 \frac{M_*^4}{R^4}  \right) \ln f(R) \right] \right. \nonumber \\
& & \left. {}-\frac{1}{2} C_\vartheta^2  \frac{R^4}{M_*^4} \left[1 + \frac{R}{M_*} \left( 1 - 2 \frac{M_*}{R} + \frac{1}{2} \frac{M_*^2}{R^2} \right) \ln f(R) + \frac{R^2}{4 M_*^2}  \left( 1-\frac{M_*}{R} \right) f(R) [\ln f(R)]^2  \right] \right\} \nonumber \\
& & {}-\frac{2J J_{\mbox{\tiny CS}}}{R^3} +  M_{\mbox{\tiny CS}}\,,
\ea
\ew
where $M_{\mbox{\tiny CS}}$ is an integration constant that corresponds to a part of the correction to the mass if one expands the above solutions about $R = \infty$. For the $\ell =2$ mode, we obtain similar solutions, shown in Appendix~\ref{l2}, but with another integration constant $C_Q$. If one is to deal with a BH, one needs to impose regularity at the horizon. By setting $J_{\mbox{\tiny CS}} =0$, $C_\vartheta=0$, $M_{\mbox{\tiny CS}}=0$ and $C_Q=(3015/14336) \xi_{\mbox{\tiny CS}} (J^2/M_*^8)$, the exterior metric in $(t,r,\theta,\phi)$ coordinates reduces to the previously found BH solution~\cite{kent-CSBH}.

\subsection{Asymptotic Behavior at Spatial Infinity:\\* Scalar Dipole Charge and mass\\* Quadrupole Moment Deformation}

The asymptotic behaviors of $h_0^\mathrm{ext}$ and $w_1^\mathrm{ext}$ about spatial infinity is
\begin{equation}
\label{eq:asymptotic-deltaM-deltaJ}
h_0^\mathrm{ext} = -\frac{\delta M}{R} + \mathcal{O}\left( \frac{1}{R^{2}} \right), \quad w_1^\mathrm{ext} = \frac{2 \delta J}{R^3} + \mathcal{O}\left( \frac{1}{R^4} \right)
\end{equation}
with
\ba
\delta M &=& M_{\mbox{\tiny CS}} + \frac{25}{1536} \xi_{\mbox{\tiny CS}} \frac{C_\vartheta^2 J^2}{M_*^7}\,,  \\
\delta J &=& J_{\mbox{\tiny CS}} - \frac{25}{384} \xi_{\mbox{\tiny CS}} \frac{C_\vartheta J}{M_*^4}\,.
\ea
Physically, $\delta M$ and $\delta J$ correspond to the CS corrections in the mass and the angular momentum respectively.

Let us now define the observable mass $\tilde{M}$ and angular momentum $\tilde{J}$ as measured by an observer at spatial infinity via
\be
\tilde{M} \equiv M_* + \delta M, \quad \tilde{J} \equiv J + \delta J\,,
\ee
such that
\be
g_{tt}^{\mathrm{ext}} \sim -\left( 1- \frac{2 \tilde{M_*}}{R} \right)\,,
\qquad
g_{t\phi}^{\mathrm{ext}} \sim - \frac{2 \tilde{J}}{R} \sin^2 \theta\,,
\ee
near spatial infinity. The exterior metric and the scalar field obtained in the previous subsection can be rewritten in terms of $\tilde{M}$ and $\tilde{J}$ by replacing $M_* \to \tilde{M}- \delta M$ and $J \to \tilde{J}-\delta J$.

\if

With these replacements, the metric perturbations now have the asymptotic behaviors as
\ba
w_1^\mathrm{ext} &=& \frac{5}{24} \xi_{\mbox{\tiny CS}} \frac{(C_\vartheta-3) \tilde{J}}{\tilde{M} R^6} + \mathcal{O}\left( \frac{1}{R^7} \right)\,,  \\
h_0^\mathrm{ext} &=& \frac{5}{6912} \xi_{\mbox{\tiny CS}} \frac{(C_\vartheta-3)^2 \tilde{J}^2}{\tilde{M}^3 R^5} + \mathcal{O}\left( \frac{1}{R^6} \right)\,, \\
m_0^\mathrm{ext} &=& -\frac{25}{6912} \xi_{\mbox{\tiny CS}} \frac{(C_\vartheta-3)^2 \tilde{J}^2}{\tilde{M}^4 R^3} + \mathcal{O}\left( \frac{1}{R^4} \right)\,, \\
h_2^\mathrm{ext} &=& \frac{5}{6912} \xi_{\mbox{\tiny CS}} \frac{(C_\vartheta-3)^2 \tilde{J}^2}{\tilde{M}^3 R^5} + \mathcal{O}\left( \frac{1}{R^6} \right)\,, \\
k_2^\mathrm{ext} &=& \frac{1}{3840} \frac{25 \xi_{\mbox{\tiny CS}} C_\vartheta (C_\vartheta-8) \tilde{J}^2 + 2048 C_Q \tilde{M}^8}{\tilde{M}^5 R^3} \nonumber \\
& & + \mathcal{O}\left( \frac{1}{R^4} \right)\,, \\
m_2^\mathrm{ext} &=& \frac{1}{3840} \frac{25 \xi_{\mbox{\tiny CS}} C_\vartheta (C_\vartheta-8) \tilde{J}^2 + 2048 C_Q \tilde{M}^8}{\tilde{M}^5 R^2} \nonumber \\
& & + \mathcal{O}\left( \frac{1}{R^3} \right)\,.
\ea
From Eq.~\eqref{metric-ansatz-rth},

\f

The full exterior metric can then be written as $g_{\mu \nu}^{\mathrm{ext}} = g_{\mu \nu}^{\mathrm{ext},{\mbox{\tiny GR}}} + g_{\mu \nu}^{\mathrm{ext},{\mbox{\tiny CS}}}$, where $g_{\mu \nu}^{\mathrm{ext},{\mbox{\tiny GR}}}$ is the Hartle-Thorne metric with mass $\tilde{M}$ and angular momentum $\tilde{J}$, while $g_{\mu \nu}^{\mathrm{ext},{\mbox{\tiny CS}}}$ is a CS correction. The latter, in $(t,r,\theta,\phi)$ coordinates, has the following asymptotic behavior about spatial infinity
\ba
\label{eq:gtt-CS-asymp}
g_{tt}^{\mathrm{ext}, {\mbox{\tiny CS}}} &=&  \frac{1}{3840} \frac{25 \xi_{\mbox{\tiny CS}} C_\vartheta (C_\vartheta-8) \tilde{J}^2 + 2048 C_Q \tilde{M}^8}{\tilde{M}^5 r^3} \nonumber \\
& & {}\times (3 \cos^2 \theta -1) + \mathcal{O}\left( \frac{1}{r^4} \right) \\
g_{rr}^{\mathrm{ext}, {\mbox{\tiny CS}}} &=&  \frac{1}{3840} \frac{25 \xi_{\mbox{\tiny CS}} C_\vartheta (C_\vartheta-8) \tilde{J}^2 + 2048 C_Q \tilde{M}^8}{\tilde{M}^5 r^3} \nonumber \\
& & {}\times (3 \cos^2 \theta -1) + \mathcal{O}\left( \frac{1}{r^4} \right) \\
g_{\theta \theta}^{\mathrm{ext}, {\mbox{\tiny CS}}} &=&  \frac{1}{3840} \frac{25 \xi_{\mbox{\tiny CS}} C_\vartheta (C_\vartheta-8) \tilde{J}^2 + 2048 C_Q \tilde{M}^8}{\tilde{M}^5 r} \nonumber \\
& & {}\times (3 \cos^2 \theta -1) + \mathcal{O}\left( \frac{1}{r^2} \right) \\
g_{\phi \phi}^{\mathrm{ext}, {\mbox{\tiny CS}}} &=&  \sin^2 \theta g_{\theta\theta}^{\mathrm{ext}, {\mbox{\tiny CS}}} \\
g_{t \phi}^{\mathrm{ext}, {\mbox{\tiny CS}}} &=&  \frac{5}{24} \xi_{\mbox{\tiny CS}} \frac{(3- C_\vartheta) \tilde{J}}{\tilde{M} r^4} \sin^2 \theta + \mathcal{O}\left( \frac{1}{r^5} \right)\,.
\ea
One can read off the correction to the gravitational potential per unit mass from Eq.~\eqref{eq:gtt-CS-asymp} as
\be
\label{eq:UA_CS}
\delta U^{{\mbox{\tiny CS}}} = -\frac{3}{r^{3}} Q \hat{S}^{i}\hat{S}^{j} n^{<ij>}
\ee
with
\begin{equation}
Q \equiv \frac{1}{7680} \frac{25 \xi_{\mbox{\tiny CS}} C_\vartheta (C_\vartheta-8) \tilde{J}^2 + 2048 C_Q \tilde{M}^8}{\tilde{M}^5} \,.
\label{eq:QA}
\end{equation}
Here, $\hat{S}^{i}$ is the unit spin angular momentum vector of the NS and $n^i$ is the unit vector to the field point. Notice that this definition of $\hat{S}^{i}$ \emph{differs} from that used in~\cite{quadratic}, because
in the latter $\hat{S}^{i}$ was not a proper unit vector. $Q$ corresponds to the correction to the quadrupole moment\footnote{This definition of quadrupole moment is different from the usual one given in e.g.~\cite{poisson-quadrupole,laarakkers} by a factor of 2. Although somewhat unconventional, we choose here to continue using the definitions of Ref.~\cite{kent-CSGW}.} and in the non-spinning BH case, Eq.~\eqref{eq:QA} reduces to $Q_\mathrm{BH} = (201/3584) (\xi_{\mbox{\tiny CS}}/\tilde{M}^4) \tilde{M}^3 \tilde{\chi}^2$, which agrees with~\cite{kent-CSGW}.

Finally, the asymptotic behavior of the scalar field is
\be
\vartheta^\mathrm{ext} =  \frac{5(3- C_\vartheta)}{24} \frac{\alpha}{\beta} \tilde{\chi} \frac{ \cos\theta}{ r^2} + \mathcal{O}\left( \frac{1}{r^3} \right)\,,
\ee
where $\tilde{\chi} \equiv \tilde{J}/\tilde{M}^2$. Comparing this to $\vartheta = \mu^i n_i/r^2 + \mathcal{O}(r^{-3})$ (see Eq.~(57) of Ref.~\cite{quadratic}) where $\mu^i$ is the scalar dipole charge, we can extract $\mu^i$:
\begin{equation}
\label{eq:mu-J-relation}
\mu^i = \frac{5(3-C_\vartheta)}{24} \frac{\alpha}{\beta} \tilde{\chi} \hat{S}^i\,.
\end{equation}
 For later convenience, we introduce the dimensionless scalar charge $\bar{\mu}$ in terms of $C_\vartheta$ and the compactness of the NS $C$:
\begin{equation}
\label{eq:mu-mubar-relation}
\bar{\mu} \equiv \frac{3-C_\vartheta}{3} \frac{1}{C^3}, \quad C \equiv \frac{\tilde{M}}{\mathcal{R}_*}\,.
\end{equation}
By using this $\bar{\mu}$, Eq.~\eqref{eq:mu-J-relation} can be rewritten as
\be
\label{eq:mu-S-relation}
\mu^{i} = \frac{5}{8}\frac{\alpha}{\beta}C^{3}\tilde{\chi} \hat{S}^{i}\bar{\mu}\,.
\ee
From Eq.~(59) of Ref.~\cite{quadratic}, we expect $\mu^i$ to be proportional to $C^3$, and hence, we have factored out $C^3$ from the definition of $\bar{\mu}$ above. In the non-spinning BH case, $\bar{\mu}$ reduces to $\bar{\mu}_\mathrm{BH} = 8$, because then $C^{3} = 1/8$, so that $\bar{\mu}_{\mathrm{BH}} C^{3} = 1$.

Notice that the difference between $\tilde{M}$ and $M_*$ (and also between $\tilde{J}$ and $J$) in Eqs.~\eqref{eq:gtt-CS-asymp}--\eqref{eq:mu-S-relation} is irrelevant to the order of approximation considered here. Hence, we can freely set $\tilde{M} = M_*$ (and $\tilde{J} = J$) in these equations.

\section{Isolated Neutron Star Solution:\\* Interior Fields}
\label{sec:interior}

In order to determine $\bar{\mu}$, $Q$, $\delta M$ and $\delta J$, we need to calculate the integration constants that appear in the definitions of these quantities. This can be achieved by matching the interior solutions to the exterior ones at the NS surface. In this section, we first explain the matching conditions at the surface. We then obtain the asymptotic behavior of the metric perturbations at the NS center to obtain the initial conditions for the numerical integration of the interior solution. We conclude with an explanation of the numerical procedure adopted in this paper. We use the superscript ``int'' to refer to interior quantities.

\subsection{Boundary Conditions at the Surface}

In the previous subsection, we imposed asymptotic flatness at spatial infinity to find the exterior solutions but 4 integration constants remain, which are to be determined by the boundary conditions at the NS surface. At $\mathcal{O}(\alpha'^0 \chi'^0)$, we solve Eqs.~\eqref{tt-zeroth},~\eqref{RR-zeroth} and~\eqref{TOV-zeroth}, together with the EoS, for $\nu$, $\lambda$, $\rho$ and $p$. At the NS surface, we impose the continuity condition
\be
p(\mathcal{R}_*) =0, \quad e^{\nu (\mathcal{R}_*)} = e^{-\lambda (\mathcal{R}_*)} = 1-\frac{2 M_*}{\mathcal{R}_*}\,.
\label{BC-zeroth}
\ee
The first equation determines $\mathcal{R}_*$ while the second equation determines $M_*$.

At $\mathcal{O}(\alpha'^0 \chi')$, we solve Eq.~\eqref{omega1RR} for $\tilde\omega_{1}$. Since this is a second-order differential equation, we need to impose two boundary conditions at the surface. We choose one to be the continuity condition for $\tilde{\omega}_1$, i.e.
\be
\tilde{\omega}_1 (\mathcal{R}_*) = \Omega_* - \frac{2 J}{\mathcal{R}_*^3}\,.
\label{BC-omega1-1}
\ee
The other condition can be obtained by integrating Eq.~\eqref{omega1RR} from $\mathcal{R}_*-\epsilon$ to $\mathcal{R}_* + \epsilon$ and taking the limit as $\epsilon \rightarrow 0$:
\ba
\left[  \frac{d \tilde{\omega}_1}{dR} (\mathcal{R}_*) \right] &=& \lim_{\epsilon \rightarrow 0} \left\{ \frac{d \tilde{\omega}_1^\mathrm{ext}{}}{dR}  (\mathcal{R}_* + \epsilon) - \frac{ d \tilde{\omega}_1^\mathrm{int}{}}{d R}  (\mathcal{R}_* - \epsilon) \right\}  \nonumber \\
&=& \lim_{\epsilon \rightarrow 0} \int^{\mathcal{R}_* + \epsilon}_{\mathcal{R}_* - \epsilon} \frac{d^2 \tilde{\omega}_1}{d R^2} dR\,, \nonumber \\
&=& \lim_{\epsilon \rightarrow 0} \int^{\mathcal{R}_* + \epsilon}_{\mathcal{R}_* - \epsilon} \left\{ - 4\frac{1- \pi R^2 (\rho +p) e^{\lambda}}{R}\frac{d \tilde{\omega}_1}{d R} \right. \nonumber \\
& & \left.  + 16 \pi (\rho +p) e^{\lambda} \tilde{\omega}_1\right\} dR\,,
\label{omega-junc}
\ea
where we have introduced the notation
\be
[ A(\mathcal{R}_*) ] \equiv \lim_{\epsilon \to 0} \left[ A^\mathrm{ext} (\mathcal{R}_* + \epsilon) - A^\mathrm{int} (\mathcal{R}_*- \epsilon) \right]
\ee
for any quantity $A(R)$. Since the integrand is bounded, the right hand side of Eq.~\eqref{omega-junc} vanishes, leading to
\be
\left[ \frac{d \tilde{\omega}_1}{d R} (\mathcal{R}_*) \right] = 0\,.
\label{BC-omega1-2}
\ee

Let us introduce the moment of inertia $I$, defined by~\cite{hartle1967}
\be
\label{eq:I-definition}
I \equiv \frac{J}{\Omega_*}\,.
\ee
Equation~\eqref{BC-omega1-1} can then be rewritten as
\be
\tilde{\omega}_1 (\mathcal{R}_*) = \Omega_* \left( 1 - \frac{2 I}{\mathcal{R}_*^3} \right)\,.
\label{BC-omega1-3}
\ee
By Eqs.~\eqref{tt-zeroth},~\eqref{RR-zeroth},~\eqref{omega1RR} and~\eqref{BC-omega1-2}, $I$ can be expressed in integral form as~\cite{hartle1967,kalogera-psaltis}
\be
I = \frac{8\pi}{3} \frac{1}{\Omega_*} \int_0^{\mathcal{R}_*} \frac{e^{-(\nu + \lambda )/2} R^5 (\rho + p) \tilde{\omega}_1}{R-2M} dR\,.
\label{I}
\ee
This expression is valid both in GR, as well as in dynamical CS gravity~\cite{alihaimoud-chen}.

As in the $\mathcal{O}(\alpha'^0 \chi')$ case, the scalar evolution equation [Eq.~\eqref{vartheta1RR}] is a second-order differential equation. Since all terms in this equation are bounded at the surface, we impose the continuity and differentiability conditions
\be
[ \vartheta_1 (\mathcal{R}_*) ] = 0 = \left[ \frac{d \vartheta_1}{d R} (\mathcal{R}_*) \right]\,.
\label{BC-vartheta}
\ee
At $\mathcal{O}(\alpha'^2 \chi')$, the non-vanishing component of the field equations is given by Eq.~\eqref{w1RR}, which is also a second-order differential equation, but here the situation is slightly different. First, we impose the continuity condition at the surface
\be
[ w_1 (\mathcal{R}_*) ] = 0\,.
\label{BC-w-1}
\ee
Next, we integrate Eq.~\eqref{w1RR} from $\mathcal{R}_*-\epsilon$ to $\mathcal{R}_* + \epsilon$ and take the limit $\epsilon \to 0$ as
\ba
& & \left[  \frac{d w_1}{dR} (\mathcal{R}_*) \right] \nonumber \\
&=& \lim_{\epsilon \rightarrow 0} \int^{\mathcal{R}_* + \epsilon}_{\mathcal{R}_* - \epsilon} \frac{d^2 w_1}{d R^2} dR\,, \nonumber \\
&=& \lim_{\epsilon \rightarrow 0} \int^{\mathcal{R}_* + \epsilon}_{\mathcal{R}_* - \epsilon} \left\{ - 4 \frac{1-\pi R^2 (\rho +p) e^\lambda}{R} \frac{d w_1}{d R} \right. \nonumber \\
& & \left. +16 \pi (\rho +p) e^\lambda w_1 \right. \nonumber \\
& & \left.  - \frac{128 \pi^2 \alpha e^{(\nu+\lambda)/2}}{R^3} \left[ \delta R \frac{d \vartheta_1}{d R} + \left( R \frac{d\rho}{dR} - \delta \right) \vartheta_1 \right] \right\} dR\,, \nonumber \\
&=& - 128 \pi^2 \alpha \lim_{\epsilon \rightarrow 0} \int^{\mathcal{R}_* + \epsilon}_{\mathcal{R}_* - \epsilon} \frac{e^{(\nu+\lambda)/2} \vartheta_1}{R^2} \frac{d\rho}{dR} dR\,.
\label{w-junc}
\ea
where in the last equality, we only kept the term in the integrand that is not bounded\footnote{The factor $d \rho /dR$ is not bounded since $\rho$ is discontinuous at the surface.}. This equation can be further simplified by integrating by parts:
\ba
\left[  \frac{d w_1}{dR} (\mathcal{R}_*) \right] &=& - 128 \pi^2 \alpha \lim_{\epsilon \rightarrow 0} \int^{\mathcal{R}_* + \epsilon}_{\mathcal{R}_* - \epsilon}  \left\{ \frac{d}{dR} \left( \frac{e^{(\nu+\lambda)/2} \vartheta_1}{R^2} \rho \right) \right. \nonumber \\
& & \left. - \rho \frac{d}{dR} \left( \frac{e^{(\nu+\lambda)/2} \vartheta_1}{R^2} \right)  \right\} dR\,.
\ea
Since the integrand in the second term is bounded, this term vanishes, so that one finally finds the jump condition
\ba
\left[  \frac{d w_1}{dR} (\mathcal{R}_*) \right] &=& 128 \pi^2 \alpha  \frac{\vartheta_1 (\mathcal{R}_*)}{\mathcal{R}_*^{2}} \rho (\mathcal{R}_*)\,.
\ea
This agrees with the condition found in~\cite{alihaimoud-chen} in the limit $\epsilon \to 0$.
Since the density of the NS at the surface is typically about 7 orders of magnitude lower than the mean density, the above jump condition shows that $d w_1/dR$ is almost continuous at the surface~\cite{alihaimoud-chen}. Therefore, in the numerical calculation below, we adopt the differentiability condition~\cite{alihaimoud-chen}
\be
\left[  \frac{d w_1}{dR} (\mathcal{R}_*) \right] = 0\,.
\label{BC-w-2}
\ee
At $\mathcal{O}(\alpha'^2 \chi'^2)$, the evolution equations are first-order. We thus impose the continuity conditions
\begin{equation}
\begin{aligned}
\left[h_0 (\mathcal{R}_*) \right] &= 0 = [m_0 (\mathcal{R}_*) ]\,,\\
[h_2 (\mathcal{R}_*) ] &= 0 = [k_2 (\mathcal{R}_*) ]\,.
\end{aligned}
\label{BC-2nd-order}
\end{equation}

\subsection{Asymptotic Behavior at the Stellar Center}
\label{sec:asymptotia-center}

Before we numerically integrate the interior equations, we must understand the asymptotic behavior of the pressure, density, metric perturbations, and the scalar field at the NS center. First, let us focus on the pure GR case. In the non-spinning sector, we can use Eqs.~\eqref{tt-zeroth},~\eqref{RR-zeroth} and~\eqref{TOV-zeroth} to  asymptotically expand $\rho$, $p$, $M$ and $\nu$ about $R \ll M$ via
\begin{align}
\label{rho0}
\rho (R) &= \rho_c + \rho_2 R^2 + \mathcal{O}(R^3)\,, \\
\label{p-center}
p (R) &= p_c - \frac{2\pi}{3} (\rho_c + p_c) (\rho_c + 3 p_c) R^2 + \mathcal{O}(R^3)\,, \\
\label{center-M}
M(R) &= \frac{4\pi}{3} \rho_c R^3 + \frac{4\pi}{5} \rho_2 R^5 + \mathcal{O}(R^6)\,, \\
\label{nu0}
\nu (R) &= \nu_c + \frac{4\pi}{3} (\rho_c + 3 p_c) R^2 + \mathcal{O}(R^3) \quad (R \to 0^{+})\,,
\end{align}
where $\rho_c$, $p_c=p(\rho_c)$ and $\nu_c$ are the density, pressure and $\nu$ at the NS center, respectively. $\rho_2$ can be expressed in terms of $\rho_c$ and $p_c$ through Eq.~\eqref{p-center}, the TOV equation and the EoS. $\nu_c$ is determined by matching the interior solution with the exterior one at the surface. To first-order in spin, the solution to Eq.~\eqref{omega1RR} that is regular at the center is
\be
\tilde{\omega}_1 (R) = \tilde{\omega}_c + \frac{8\pi}{5} (\rho_c + p_c) \tilde{\omega}_c R^2 + \mathcal{O}(R^3) \quad (R \to 0^{+})\,,
\label{omega0}
\ee
where $\tilde{\omega}_c$ is a constant that is to be determined by the matching condition at the surface.

Let us now focus on the scalar field and the CS corrections to the metric perturbation. For the former, the solution to Eq.~\eqref{vartheta1RR} that is regular at the center has the asymptotic form
\be
\vartheta_1 (R) = \vartheta'_c R + \frac{2\pi}{15} (5\rho_c - 3p_c) \vartheta'_c R^3 + \mathcal{O}(R^4) \quad (R \to 0^{+})\,.
\label{vartheta0}
\ee
$\vartheta'_c$ is a constant that corresponds to $d\vartheta/dR$ at the center. From Eq.~\eqref{w1RR}, the metric perturbation at first order in spin has the asymptotic behavior
\ba
w_1 (R) &=& w_c + \frac{8\pi}{5} \left[ (\rho_c + p_c) w_c -16\pi \alpha e^{\nu_c/2} \rho_2 \vartheta'_c \right] R^2  \nonumber \\
& &{}+ \mathcal{O}(R^3) \quad (R \to 0^{+})\,.
\label{w0}
\ea
Again, $\vartheta'_c$ and $w_c$ are to be determined by matching at the surface.
For the $\ell =0$ mode corrections at quadratic order in spin, Eqs.~\eqref{xi0}--\eqref{h0R} yield
\ba
\label{center-h0}
h_0 (R) &=& h_{0c} + \frac{128\pi^2}{3} \alpha (\rho_c + p_c) (\rho_c + 3 p_c) e^{-\nu_c/2} \tilde{\omega}_c \vartheta'_c R^2 \nonumber \\
& & {}+ \mathcal{O}(R^3), \\
\label{center-m0}
m_0 (R) &=& \frac{2\pi}{3} \vartheta'_c \left[ 64 \pi \alpha (\rho_c + p_c) e^{-\nu_c/2} \tilde{\omega}_c + \beta \vartheta'_c \right] R^3 \nonumber \\
& & {}+ \mathcal{O}(R^4), \\
\label{center-xi0}
\xi_0 (R) &=& -\frac{\tilde{\omega}_c e^{-\nu_c/2} }{3\pi(\rho_c + 3 p_c)}  \left[ 64 \pi^2 \alpha (\rho_c + p_c)  \vartheta'_c + e^{-\nu_c/2} w_c \right]  R \nonumber \\
& & {}+ \mathcal{O}(R^2)  \quad (R \to 0^{+})\,.
\label{center0}
\ea
Similarly, Eqs.~\eqref{xi2} and~\eqref{m2}--\eqref{RR2} yield the asymptotic behaviours for the $\ell =2$ mode:
\ba
\label{center-h2}
h_2 (R) &=& h_{2c}'' R^2 + \mathcal{O}(R^3)\,, \\
\label{center-k2}
k_2 (R) &=& -\frac{1}{3} \left[ 3 h_{2c}'' -128 \pi^2 \alpha (\rho_c + p_c) e^{-\nu_c/2} \vartheta'_c \tilde{\omega}_c \right. \nonumber \\
& & \left. + 8 \pi \beta \vartheta'_c{}^2 \right] R^2 + \mathcal{O}(R^3)\,, \\
\label{center-m2}
m_2 (R) &=& -\frac{1}{3} \left[ 3 h_{2c}'' -128 \pi^2 \alpha (\rho_c + p_c) e^{-\nu_c/2} \vartheta'_c \tilde{\omega}_c \right. \nonumber \\
& & \left.  + 8 \pi \beta \vartheta'_c{}^2 \right] R^3 + \mathcal{O}(R^4)\,, \\
\label{center-xi2}
\xi_2 (R) &=& - \frac{3 h_{2c}'' - 2 e^{-\nu_c} \tilde{\omega}_c w_c}{4\pi (\rho_c + 3 p_c)} R + \mathcal{O}(R^2)  \quad (R \to 0^{+})\,. \nonumber \\
\ea
Here, $h_{2c}''$ is a constant that needs to be determined by matching at the NS surface.

\subsection{Numerical Method}

Let us now explain how we solve the interior equations numerically to obtain the corresponding interior solutions.

\subsubsection{Solution to $\mathcal{O}(\chi'{}^0)$}

First, we must obtain the GR solutions at zeroth order in spin. There are 4 unknown functions that need to be determined, $\nu$, $\lambda$ [or equivalently $M$ through Eq.~\eqref{M-lambda}], $\rho$ and $p$. The 4 equations that we need to solve are Eqs.~\eqref{tt-zeroth},~\eqref{RR-zeroth},~\eqref{TOV-zeroth}, given the EoS. Notice that Eqs.~\eqref{tt-zeroth},~\eqref{TOV-zeroth} and the EoS form a closed system for $M$, $\rho$ and $p$. We solve these equations as an initial value problem using an adaptive step-size, fourth-order Runge-Kutta method (from the GSL library~\cite{gsl}) with an accuracy of $10^{-3}$ on each step, starting at $R=R_\epsilon$ toward $R=\mathcal{R}_*$, with the initial conditions given in Eqs.~\eqref{rho0} and~\eqref{center-M}, where $R_\epsilon$ is the core radius $R_\epsilon/\mathcal{R}_* \ll 1$. The radius of the NS can be obtained by finding the radius $\mathcal{R}_*$ where $p(\mathcal{R}_*)=0$ and the mass of the NS is given by $M_* = M(\mathcal{R}_*)$. We repeat the calculation for various $\rho_c$ to obtain a mass-radius relation.

The remaining equation, Eq.~\eqref{RR-zeroth}, can be solved for $\nu$, again using an adaptive step-size, fourth-order Runge-Kutta method with an accuracy of $10^{-3}$ per step, from $R=R_\epsilon$ to $\mathcal{R}_*$ and with the initial condition given in Eq.~\eqref{nu0}. However, one must be careful that $\nu(\mathcal{R}_*)$ satisfies the boundary condition given in Eq.~\eqref{BC-zeroth}. Suppose we obtained the trial solution $\nu_\mathrm{tr}(R)$ by using the initial condition $\nu(R_\epsilon) =(4\pi/3) (\rho_c + 3 p_c) R_\epsilon^2$ (i.e. $\nu_c = 0$). Since Eq.~\eqref{RR-zeroth} is shift-invariant, $\nu_\mathrm{tr}(R)$ plus a constant $C_\nu$ is also a solution. The new solution $\nu_\mathrm{tr}(R)+C_\nu$ will satisfy the correct boundary condition provided
\be
e^{\nu_\mathrm{tr}(\mathcal{R}_*)+C_\nu} = f(\mathcal{R}_*)\,,
\ee
which then yields~\cite{yunespsaltis}
\be
C_\nu = \ln f(\mathcal{R}_*) - \nu_\mathrm{tr} (\mathcal{R}_*)\,.
\ee

\subsubsection{Solution to $\mathcal{O}(\chi')$}
\label{GR-1st}

Next, we solve Eq.~\eqref{omega1RR}, a second-order differential equation, for $\tilde{\omega}_1$.  We take advantage of the fact that Eq.~\eqref{omega1RR} is linear and homogeneous~\cite{kalogera-psaltis} and solve this equation with the initial condition given in Eq.~\eqref{omega0}. We choose a specific value for $\tilde{\omega}_c$ and solve the equation from $R=R_\epsilon$ to $\mathcal{R}_*$ using a fourth-order Runge-Kutta method to obtain the trial solution $\tilde{\omega}_\mathrm{tr}(R)$ and the trial moment of inertia $I_\mathrm{tr}$ calculated from Eq.~\eqref{I}. The solution we seek is one that satisfies the boundary condition of Eq.~\eqref{BC-omega1-3}. It is improbable that $\tilde{\omega}_\mathrm{tr}$ will satisfy this boundary condition. However, due to the homogeneity of the differential equation, the product of $\tilde{\omega}_\mathrm{tr}$ and a constant is also a solution. Hence, we can construct a new solution via $C_{\tilde{\omega}} \tilde{\omega}_\mathrm{tr}$, where $C_{\tilde{\omega}}$ is a constant, which will satisfy the boundary condition, provided that
\be
C_{\tilde{\omega}} \tilde{\omega}_\mathrm{tr} (\mathcal{R}_*) = \Omega_* \left( 1 - \frac{2 C_{\tilde{\omega}} I_\mathrm{tr}}{\mathcal{R}_*^3} \right)\,,
\ee
which then yields
\be
C_{\tilde{\omega}} = \frac{\Omega_* \mathcal{R}_*^3}{\tilde{\omega}_\mathrm{tr}(\mathcal{R}_*) \mathcal{R}_*^3 + 2 I_\mathrm{tr} \Omega_*}\,.
\ee
Clearly, $\Omega_{*}$ (like $\rho_{c}$ or $p_{c}$) is a quantity that must be specified {\emph{a priori}} and controls
how rapidly the NS is rotating.

The scalar field satisfies Eq.~\eqref{vartheta1RR}, which is not homogeneous or shift-invariant. Since the differential equation is linear, we take the following approach~\cite{hartle1967}. First, we solve Eq.~\eqref{vartheta1RR} with arbitrary values for $\vartheta_1 (R_\epsilon)$ and $\vartheta_1' (R_\epsilon)$ such that the solution satisfies regularity at the NS center [i.e. with an arbitrary value of $\vartheta_1'$ in Eq.~\eqref{vartheta0}] to obtain a particular solution $\vartheta_p^\mathrm{int} (R)$ in the interior region. Next, we look for a homogeneous solution $\vartheta_h^\mathrm{int} (R)$ in the interior region by solving Eq.~\eqref{vartheta1RR} with a vanishing source term on the right-hand-side with arbitrary $\vartheta_1 (R_\epsilon)$ and $\vartheta_1' (R_\epsilon)$ such that $\vartheta_h^\mathrm{int} (R)$ satisfies regularity at the NS center. The asymptotic behavior of $\vartheta_h^\mathrm{int} (R)$ that is regular at the NS center is also given by Eq.~\eqref{vartheta0}. With $\vartheta_p^\mathrm{int} (R)$ and $\vartheta_h^\mathrm{int} (R)$ at hand, one can construct a generic solution to Eq.~\eqref{vartheta1RR} in the interior region that satisfies regularity at the NS center as
\be
\vartheta_1^\mathrm{int} (R) = \vartheta_p^\mathrm{int} (R) + C_\vartheta^h \vartheta_h^\mathrm{int} (R)\,,
\label{theta-sol-generic}
\ee
where $C_\vartheta^h$ is a constant. The constants $C_\vartheta$ in Eq.~\eqref{vartheta1-ext} and $C_\vartheta^h$ in Eq.~\eqref{theta-sol-generic} are determined by the matching condition at the NS surface given in Eq.~\eqref{BC-vartheta}. One can solve Eq.~\eqref{BC-vartheta} for  $C_\vartheta$ and $C_\vartheta^h$ algebraically in terms of $\vartheta_p^\mathrm{int} (\mathcal{R}_*)$ and $\vartheta_h^\mathrm{int} (\mathcal{R}_*)$.

The CS correction to the gravitomagnetic sector of the metric is controlled by Eq.~\eqref{w1RR}, which is a second order differential equation for $w_1$. Again, this equation is not shift invariant or scale invariant, and thus, we must solve it as we did the scalar field above. The matching conditions at the NS surface are given in Eqs.~\eqref{BC-w-1} and~\eqref{BC-w-2}, while the asymptotic behavior of the solution that is regular at the NS center is given in Eq.~\eqref{w0}. The asymptotic behavior of the homogeneous solution is given by setting $\alpha = 0$ in Eq.~\eqref{w0}, which is the same as Eq.~\eqref{omega0} with $\tilde{\omega}_c$ being replaced by $w_c$.

We have that the solutions for $\vartheta_1$ and $w_1$ obtained in this way are identical to solving Eqs.~\eqref{vartheta1RR} and~\eqref{w1RR} via a purely numerical shooting method.

\subsubsection{Solution to $\mathcal{O}(\chi'{}^2)$}

For the $\ell=0$ mode at $\mathcal{O}(\alpha'^2 \chi'^2)$, we need to solve Eqs.~\eqref{m0R} and~\eqref{h0R} to obtain $h_{0c}$ in Eq.~\eqref{xi0}. Since these are two coupled first-order differential equations, we can solve them as an initial value problem with the initial conditions given in Eqs.~\eqref{center-m0} and~\eqref{center-xi0}. Then, we impose the continuity condition at the surface to obtain $M_{\mbox{\tiny CS}}$ in Eq.~\eqref{m0_ext}. The interior solution for $h_0$ can be obtained through Eq.~\eqref{xi0}. By imposing the continuity of $h_0$ at the surface, one can determine $h_{0c}$.

For the $\ell = 2$ mode, we need to solve Eqs.~\eqref{h2k2} and~\eqref{RR2} for $h_2$ and $k_2$. We take the same steps explained above for solving Eqs.~\eqref{vartheta1RR} and~\eqref{w1RR}. We first obtain the particular solutions $h_p^\mathrm{int}$ and $k_p^\mathrm{int}$ with an arbitrary choice of $h_c''$ in Eqs.~\eqref{center-h2} and~\eqref{center-k2}. Next, we look for homogeneous solutions $h_p^\mathrm{int}$ and $k_p^\mathrm{int}$. The initial conditions at the NS center is given by setting $\alpha=0$ and $\beta=0$ in Eqs.~\eqref{center-h2} and~\eqref{center-k2}. As in Eq.~\eqref{theta-sol-generic}, one can construct generic solutions for $h_2$ and $k_2$ that are regular at the NS center via
\ba
h_2^\mathrm{int} (R) &=& h_p^\mathrm{int} (R) + C_Q^h h_h^\mathrm{int} (R)\,, \\
k_2^\mathrm{int} (R) &=& k_p^\mathrm{int} (R) + C_Q^h k_h^\mathrm{int} (R)\,,
\ea
where $C_Q^h$ is a constant that is to be determined, together with another constant $C_Q$ in Eqs.~\eqref{h2-ext} and~\eqref{k2-ext}, with the matching condition at the NS surface given in Eq.~\eqref{BC-2nd-order}. One can solve Eq.~\eqref{BC-2nd-order} for $C_Q$ and $C_Q^h$ algebraically in terms of $h_2(\mathcal{R}_*)$ and $k_2(\mathcal{R}_*)$. We have checked that the solutions obtained in this fashion are identical to numerically solving Eqs.~\eqref{h2k2} and~\eqref{RR2} via a Riccati method~\cite{dieci1,dieci2,takata}\footnote{We had difficulty in stably carrying out the shooting method due to the nearly scale-invariant structure of Eqs.~\eqref{h2k2} and~\eqref{RR2} for most of the region in parameter space.}.

\section{Numerical Solutions}
\label{sec:solutions}

In this section, we show the numerically obtained NS solutions using 4 different EoSs: APR~\cite{APR}, SLy~\cite{SLy,shibata-fitting}, LS220~\cite{LS, ott-EOS} and Shen~\cite{Shen1,Shen2,ott-EOS}.

We begin by showing that the results obtained here reproduce previously obtained results, but for a wider range of EoSs. The left panel of Fig.~\ref{fig:MR} shows results at zeroth order in spin, i.e.~the NS mass-radius relation. Observe that the SLy curve agrees with~\cite{shibata-fitting}. Notice also that the 4 EoSs we adopt in this paper are consistent with the existence of the recently found $(1.97\pm0.04 ) M_\odot$ millisecond pulsar J1614-2230~\cite{1.97NS}. The right panel of Fig.~\ref{fig:MR} presents results at linear order in spin, i.e.~the moment of inertia as a function of mass. The APR curve agrees exactly with the results of~\cite{pani-NS-EDGB}). The left panel of Fig.~\ref{fig:theta} shows $\vartheta_1 /(\Omega_{*} \alpha/\beta)$ as a function of $R/\mathcal{R}_*$ for a NS of mass $M_*=1.4M_\odot$. Observe that $\vartheta_1$ has its maximum near the NS surface. The behavior of these curves is consistent with Fig.~5 in~\cite{alihaimoud-chen}.  The top right panel of Fig.~\ref{fig:theta} shows $\omega_1 /\Omega_{*}$ versus $R/\mathcal{R}_*$, with
\begin{equation}
\label{eq:omega-1-def}
\omega_1 \equiv \Omega_* - \tilde{\omega}_1\,.
\end{equation}
One can see that frame-dragging is strongest near the NS center.  These curves are consistent with Fig.~6 of~\cite{alihaimoud-chen}. In the bottom right panel of the same figure, we show $w_1/(\xi_{\mbox{\tiny CS}} \Omega_*)$, which is also peaked in the inner region. This behavior is also consistent with Fig.~6 of~\cite{alihaimoud-chen}.

\begin{figure*}[ht]
\begin{center}
\begin{tabular}{l}
\includegraphics[width=8.5cm,clip=true]{MR.eps}  \qquad
\includegraphics[width=8.5cm,clip=true]{IM.eps}
\end{tabular}
\caption{\label{fig:MR}
(Color online)
The mass-radius relation (left) and the moment-of-inertia
[Eq.~\eqref{eq:I-definition}] vs. mass relation (right) in GR for NSs
with various EoSs: APR (solid), SLy (dotted), LS220 (dashed), and Shen
(dot-dashed).}
\end{center}
\end{figure*}

\begin{figure*}[tp]
\includegraphics[width=9.0cm,clip=true]{Theta.eps}
\includegraphics[width=8.0cm,clip=true]{CS-Omega-Dup.eps}
\caption{\label{fig:theta}
(Color online)
Left: Scalar field profile as a function of NS radius with $M_* =
1.4M_\odot$ and for various EoSs. The linear behaviour of $\vartheta$ near
$R\to 0$ is smooth and differentiable, after multiplying by
$P_{1}(\cos\Theta)$.
Right: Correction to the metric rotation $\omega$ as a function of radius, at
linear order in spin, for various EoSs with $M_* = 1.4M_\odot$. The top
plot shows the $\mathcal{O}(\alpha'^{0})$ (GR) contribution $\omega_1/\Omega_*$
[Eq.~\eqref{eq:omega-1-def}], while the bottom plot shows the $\mathcal{O}(\alpha'^{2})$ (CS)
contribution $w_1/(\xi_{\mbox{\tiny CS}}\Omega_*)$ [Eq.~\eqref{eq:w-def}].
In all cases, only the $\ell=1$ multipole moment
survives at linear order in spin in the small coupling
limit. See Sec.~\ref{sec:asymptotia-center} for details on asymptotia.}
\end{figure*}

{\renewcommand{\arraystretch}{1.2}
\begin{table*}[ptb]
\begin{centering}
\begin{tabular}{cr@{\hspace{1em}}r@{.}l r@{.}l r@{.}l r@{.}l r@{.}l r@{.}l}
\hline
\hline
\noalign{\smallskip}
& EoS &  \multicolumn{2}{c}{$a$} &  \multicolumn{2}{c}{$b$}
&  \multicolumn{2}{c}{$c$} &  \multicolumn{2}{c}{$d$} &  \multicolumn{2}{c}{$e$} &  \multicolumn{2}{c}{Error (\%)} \\
\hline
\multirow{4}{*}{$\bar{\mu}$} &APR & 2&90 (0.022)  & -2&30 (0.067) & 1&54 (0.074) & -0&546 (0.036) & 0&0748 (0.0062) & 0&11\\
& SLy & 2&82 (0.027)  & -2&31 (0.087) & 1&62 (0.10) & -0&581 (0.052) & 0&0777 (0.0095) & 0&12\\
& LS220 & 5&87 (0.071) & -8&07 (0.25) & 6&60 (0.31) & -2&63 (0.17) & 0&400 (0.033) & 0&49\\
& Shen & 5&35 (0.049) & -6&86 (0.16) & 5&39 (0.20) & -2&07 (0.099) & 0&307 (0.018) & 0&36\\
\hline
\multirow{4}{*}{$[Q/(J^2/M_*)] (0.1/\zeta)$}& APR & -0&802 (0.034)  & -0&675 (0.10) & 0&578 (0.11) & -0&192 (0.053) & 0&0150 (0.0091) & 0&16\\
& SLy & -1&79 (0.055)  & 1&51 (0.18) & -1&37 (0.20) & 0&675 (0.10) & -0&138 (0.018) & 0&23\\
& LS220 & 6&72 (0.21) & -16&8 (0.72) & 15&6 (0.91) & -6&66 (0.49) & 1&06 (0.097) & 1&2\\
& Shen & 5&01 (0.13) & -13&0 (0.42) & 11&8 (0.49) & -4&94 (0.25) & 0&778 (0.045) & 0&67\\
\hline
\multirow{4}{*}{$-(\delta M/M_\odot) (\Omega_\mathrm{1ms}/\Omega_*)^2 (0.1/\zeta)$ }& APR & -2&42 (0.046) & -0&503 (0.14) & 0&752 (0.15) & -0&266 (0.069) & 0&0234 (0.012) & 0&18\\
& SLy & -3&48 (0.15)  & 2&18 (0.48) & -2&29 (0.55) & 1&33 (0.27) & -0&297 (0.0487) & 0&84\\
& LS220 & 4&53 (0.24) & -14&0 (0.85) & 12&3 (1.1) & -4&81 (0.58) & 0&673 (0.11) & 2&0\\
& Shen & 3&00 (0.16) & -10&1 (0.51) & 8&60 (0.61) & -3&24 (0.30) & 0&445 (0.055) & 1&5\\
\hline
\multirow{4}{*}{$(\delta J/J) (0.1/\zeta)$ }& APR & -1&43 (0.030)  & -1&73 (0.088) & 1&92 (0.095) & -0&815 (0.044) & 0&141 (0.0074) & 0&13\\
& SLy & -0&242 (0.26)  & -6&30 (0.80) & 7&80 (0.91) & -3&97 (0.44) & 0&766 (0.078) & 1&0\\
& LS220 & 5&22 (0.11) & -16&1 (0.37) & 15&9 (0.46) & -7&27 (0.24) & 1&30 (0.047) & 0&56\\
& Shen & 3&94 (0.097) & -12&9 (0.32) & 12&1 (0.37) & -5&18 (0.18) & 0&869 (0.032) & 0&98\\
\noalign{\smallskip}
\hline
\hline
\end{tabular}
\end{centering}
\caption{
\label{table:C1}
Numerical coefficients for the fitting formula of $\bar{\mu}$, $Q$,
$\delta M$, and $\delta J$ as functions of mass using the functional
form in Eq.~\eqref{C1-C2-fit}. Standard errors on each coefficient are
in parentheses. For definitions of dimensionless scalar dipole susceptibility $\bar{\mu}$,
quadrupole correction $Q$, mass shift $\delta M$, and angular momentum
shift $\delta J$, see Eqs.~\eqref{eq:mu-mubar-relation},
\eqref{eq:QA}, and \eqref{eq:asymptotic-deltaM-deltaJ}. The last column shows the maximum fractional error between values that are obtained numerically and from the fitting formula.}
\end{table*}
}

\begin{figure*}[tp]
\begin{center}
\includegraphics[width=8.5cm,clip=true]{CS-h0-Dup.eps} \qquad
\includegraphics[width=8.5cm,clip=true]{CS-h2-Dup.eps}
\caption{\label{fig:CS-h0}
(Color online)
$\ell=0$ (left) and $\ell=2$ (right) metric perturbations at second order in spin as functions of
radius, for various EoSs with $M_* = 1.4M_\odot$:
$h_0 /(\xi_{\mbox{\tiny CS}} \Omega_*^2)$ (top left) and $m_0 /(\xi_{\mbox{\tiny CS}}\Omega_*^2)$
(bottom left) as defined in Eqs.~\eqref{eq:h-decomp} and
\eqref{eq:m-decomp}; $h_2 /(\xi_{\mbox{\tiny CS}} \Omega_*^2)$ (top right) and
$k_2 /(\xi_{\mbox{\tiny CS}} \Omega_*^2)$ (bottom right)
as defined in Eqs.~\eqref{eq:h-decomp} and \eqref{eq:k-decomp}.}
\end{center}
\end{figure*}

\begin{figure*}[htp]
\begin{center}
\begin{tabular}{l}
\includegraphics[width=8.5cm,clip=true]{CS-m2-int.eps}
\includegraphics[width=8.5cm,clip=true]{CS-m2-ext.eps}
\end{tabular}
\caption{\label{fig:CS-m2}
$\ell=2$ metric perturbation at second order in spin as a function of
radius, for various EoSs with $M_* = 1.4M_\odot$:
$m_2 /(\xi_{\mbox{\tiny CS}} \Omega_*^2)$ [Eq.~\eqref{eq:m-decomp}] for the
interior (left) and the exterior (right). Notice the difference in
scales. The features in the interior correspond to the locations of
nuclear phase transitions (see Fig.~\ref{fig:rho-r}).}
\end{center}
\end{figure*}

\begin{figure}[htp]
\begin{center}
\includegraphics[width=8.5cm,clip=true]{rho-r-Dup.eps}
\caption{\label{fig:rho-r}
Mass density $\rho$ (top) and radial derivative $d\rho/dR$ (bottom) in
the NS interior, for various EoSs with $M_* = 1.4M_\odot$. Recall
that $R$ is defined by surfaces of constant $\rho$, so there is no
angular dependence to $\rho$. Nuclear phase transitions are visible as
sudden changes in slope. These give rise to the features in the metric
functions seen in Fig.~\ref{fig:CS-m2}.}
\end{center}
\end{figure}

Let us now present new results, valid to quadratic order in spin. Figures~\ref{fig:CS-h0} and~\ref{fig:CS-m2} show the metric perturbations at quadratic order in spin for the $\ell = 0$ and $\ell=2$ modes respectively. The interior behavior of $h_0$, $h_{2}$ and $k_{2}$ is similar to that of $\vartheta_1$. The strange behavior of $m_2$ in the interior reflects the strange behavior of $d\rho(R)/dR$, plotted in Fig.~\ref{fig:rho-r}. The oscillatory structure in these quantities disappears for a polytropic EoS. The oscillations in the EoS, which then propagate to $m_2$, are due to nuclear phase transitions~\cite{APR,SLy,LS,Shen1,Shen2}. We recall that we used an adaptive step-size, fourth-order Runge-Kutta method with an accuracy of $10^{-3}$ per step to solve the differential equations; therefore, the strange behavior alluded to above is not a numerical artifact and it is well-resolved.

Figure~\ref{fig:mu-M} already presented $\bar{\mu}$, $[Q/(J^2/M_*)] (0.1/\zeta)$, $\delta M (\Omega_\mathrm{1ms}/\Omega_*)^2 (0.1/\zeta)$ and $(\delta J/J) (0.1/\zeta)$ as a function of $M_*/M_\odot$ for various EoSs. This behavior can be nicely fit by
\ba
\label{C1-C2-fit}
A_i &=& \exp\left( a_i + b_i x  +c_i x^2  + d_i x^3 + e_i x^4 \right)\,, \\
A_i &\equiv & \left( \bar{\mu}, \ \frac{Q}{J^2/M_*} \frac{0.1}{\zeta}, \ - \frac{\delta M}{M_\odot} \frac{\Omega_\mathrm{1ms}^2}{\Omega_*^2} \frac{0.1}{\zeta}, \ \frac{\delta J}{J} \frac{0.1}{\zeta} \right)\,, \nonumber \\
 \ea
with $x \equiv M_*/M_\odot$. The estimated numerical coefficients are shown in Table~\ref{table:C1}, together with standard errors. Figure~\ref{fig:mu-M} also shows these fits as a function of $M_*/M_\odot$. Notice that the CS correction increases the spin angular momentum which is consistent with the result reported in~\cite{alihaimoud-chen}. Moreover, the CS correction increases the quadrupole moment while it decreases the observed mass.

\section{Neutron Star Binary Evolution}
\label{sec:binary-evolution}

Let us now study the evolution of a NS binary with component masses $m_A$, radii $\mathcal{R}_{* A}$, spin angular momenta $\chi_A m_A^2 \hat{S}_A^i$, dimensionless scalar charges $\bar{\mu}_A$ and quadrupole moment shifts $Q_A$. We use the subscript $A =(1,2)$ to denote the $A^\mathrm{th}$ binary component. We also use the subscript ``12'' to denote relative differences. The binary orbit can be described by 5 intrinsic and extrinsic parameters: semimajor axis $a$, eccentricity $e$, the angle of periastron $\omega$, inclination $\iota$, and the angle of the ascending node $\Omega$ (not to be confused with the gravitomagnetic metric perturbation $\omega$ or the NS angular velocity $\Omega_*$).

A binary's orbital evolution is affected by dissipative effects and conservative effects. The former can be obtained by calculating the energy flux and the angular momentum flux radiated via gravitational and scalar radiation, which we explain in detail in Appendix~\ref{sec:E-flux}. In this section, we focus on the latter. Here, we are only interested in CS corrections on the conservative effect relative to the leading GR contribution, and hence, we consider CS corrections to Newtonian dynamics. We follow Gauss' perturbation
method~\cite{smart,robertson-noonan,TEGP} with the conventions and notation
of~\cite{TEGP}. One must then calculate the perturbing accelerations
$\delta a_{A}^{i}$ and from them the relative perturbing acceleration
$\delta a_{12}^{i}$ (in any convenient coordinate system). This vector is then
decomposed by projecting onto a \emph{time-varying}
orthonormal triad with $e_{1}^{i}=n_{12}^{i}$ and
$e_{2}^{i}=\hat{L}^{i}$ (and $e_{3}=e_{2}\times e_{1}$ so as to complete the
orthonormal triad). In this triad, the components of
$\delta a_{12}^{i}$ are defined as~\cite{TEGP}
\begin{subequations}
\label{eq:accel-decomp}
\begin{align}
\mathscr{R} &\equiv \delta a_{12}^{i} e_{1,i}\,, \\
\mathscr{W} &\equiv \delta  a_{12}^{i} e_{2,i}\,, \\
\mathscr{S} &\equiv \delta a_{12}^{i} e_{3,i}\,,
\end{align}
\end{subequations}
where inner products are taken with a flat Euclidean metric.
With this decomposition, the osculating orbital elements evolve
secularly as
\begin{subequations}
\label{eq:oscul-evol}
\begin{align}
\dot\omega &= - \frac{p \mathscr{R}}{he}\cos\phi +
\frac{(p+r) \mathscr{S} }{h e}\sin\phi -\dot\Omega \cos\iota\,, \\
\dot{e} &= \frac{1-e^{2}}{h}\left[ a \mathscr{R}\sin\phi +
  \frac{\mathscr{S}}{e} \left( \frac{ap}{r}-r \right) \right]\,, \\
\dot{a} &= \frac{2a^{2}}{h} \left( \frac{\mathscr{S} p}{r}+\mathscr{R} e\sin\phi \right)\,, \\
{\textstyle\frac{d}{dt}\iota} &= \frac{\mathscr{W} r}{h}\cos(\omega+\phi)\,, \\
\dot\Omega &= \frac{\mathscr{W} r}{h}\sin(\omega+\phi) \csc\iota\,.
\end{align}
\end{subequations}
Here,
\be
p \equiv a(1-e^2)
\ee
is the semi-latus rectum, $r$ and $\phi$ are the quantities related to the instantaneous orbital elements, given by
\ba
r &\equiv & \frac{p}{1+e \cos\phi}\,, \\
r^2 \frac{d \phi}{dt} &\equiv & h \equiv \sqrt{m p}\,,
\ea
with $m \equiv m_1 + m_2$ and $h$ is the orbital angular momentum per unit mass.
In all of Eqs.~\eqref{eq:oscul-evol}, the right-hand sides are to be
orbit averaged as defined in Eq.~\eqref{eq:orbit-ave-def} (this is
appropriate when the time derivatives of the osculating elements are
much smaller than the orbital timescale and there are no
resonances). Of course, all of the $\phi$ dependence in $\mathscr{R},
\mathscr{W}, \mathscr{S}$ and $r$ must be included in the orbit
averaging.

In Sec.~\ref{sec:bin-evol-metric}, we calculate the secular
evolution from the near-zone metric deformation due to the CS
correction to the quadrupole. In order to calculate the secular
evolution due to the scalar field, we derive the effective
dipole interaction in Sec.~\ref{sec:dipole-interaction}. With the
acceleration due to the scalar interaction, we calculate the
scalar's correction to the binary's secular evolution.

\subsection{Metric Quadrupole Correction}
\label{sec:bin-evol-metric}

Since both the mass monopole and current dipole of each body
is reabsorbed into physically measured quantities, the leading
order CS correction is the quadrupole deformation to the metric given by Eq.~\eqref{eq:QA}. From Eq.~\eqref{eq:UA_CS}, we can
find the acceleration on body 1 due to body 2,
\begin{equation}
\label{eq:aU_1i}
\delta a^{U}_{1i} = -\partial_{i}\delta U_{2}^{{\mbox{\tiny CS}}} = -15 Q_{2}
\frac{n_{12}^{<ijk>}}{r_{12}^{4}} \hat{S}_{2}^{j} \hat{S}_{2}^{k}\,,
\end{equation}
where $r_{12}$ is the separation of the binary components and $n_{12}^i$ is the unit vector from body 2 to body 1.
The acceleration on body 2 due to body 1 is simply the above
expression with $1\leftrightarrow 2$. This gives the relative
acceleration
\begin{equation}
  \delta a^{U}_{12i} = -15 \frac{n_{12}^{<ijk>}}{r_{12}^{4}}
\left(Q_{1}\hat{S}_{1}^{j}\hat{S}_{1}^{k} +
  Q_{2}\hat{S}_{2}^{j}\hat{S}_{2}^{k} \right)\,.
\end{equation}

The orbit averages are performed on eccentric orbits lying in the
$x$---$y$ plane with the major axis along the $\hat{x}$-axis, parametrized
as~\cite{PetersMathews}
\begin{align}
x_{1}^{i} &= d_{1}( \cos \phi, \sin \phi, 0)\,, \\
x_{2}^{i} &= -d_{2}( \cos \phi, \sin \phi, 0)\,,
\end{align}
\begin{align}
d_{1} &= \frac{m_{2}}{m} d\,, &
d_{2} &= \frac{m_{1}}{m} d\,, &
d &= \frac{a(1-e^{2})}{1+e\cos\phi}\,,
\end{align}
\begin{equation}
\dot\phi = \frac{\left[ma(1-e^{2})\right]^{1/2}}{d^{2}}\,,
\end{equation}
from which $v$ and all other derived quantities can be calculated,
e.g.~$n_{12}^{i}= (x_{1}^{i} - x_{2}^{i})/d = (\cos\phi,\sin\phi,0)$.

Taking the above expression, decomposing as in
Eqs.~\eqref{eq:accel-decomp}, and performing orbit averaging with the period
$P=2\pi a^{3/2}/M^{1/2}$, we find
\begin{widetext}
\begin{subequations}
\begin{align}
\ave{ \dot\omega}_h &=
\frac{3}{a^{7/2}\sqrt{m}(1-e^{2})^2}
Q_{1} \left[ -1 + \frac{3}{2}\left(
    \hat{S}_{1,x}^{2}+\hat{S}_{1,y}^{2} \right)
  -\hat{S}_{1,z}\cot\iota \left( \hat{S}_{1,x}\sin\omega +
    \hat{S}_{1,y}\cos\omega \right)
 \right] +(1\leftrightarrow 2)\,,
\\
\ave{ {\textstyle \frac{d}{dt}}\iota }_h &=
\frac{3}{a^{7/2}\sqrt{m}(1-e^{2})^2}
Q_{1} \hat{S}_{1,z} \left( \hat{S}_{1,x}\cos\omega -
    \hat{S}_{1,y}\sin\omega  \right) +(1\leftrightarrow 2)\,,
\\
\label{eq:ave-dotOmega}
\ave{ \dot\Omega }_h &=
\frac{3}{a^{7/2}\sqrt{m}(1-e^{2})^2}
Q_{1}
\hat{S}_{1,z}\csc\iota \left( \hat{S}_{1,x}\sin\omega +
    \hat{S}_{1,y}\cos\omega \right) +(1\leftrightarrow 2)\,,
\end{align}
\end{subequations}
and $\langle\dot{e}\rangle_h = 0 = \langle \dot{a}\rangle_h$, and where
the $i=x,y,z$ components of $\hat{S}_{A,i}$ are taken in the
(non-inertial) coordinate system where the binary's orbit is in the
$x$---$y$ plane, with pericenter along the $+\hat{x}$
direction.
\end{widetext}

\subsection{Dipole Interaction and Scalar Force Correction}
\label{sec:dipole-interaction}

To derive the effective dipole interaction, we start by finding an
effective interaction Lagrangian between one of the body's scalar dipole moments
and the scalar field. This comes from the cross-interaction part of the kinetic term of the scalar field in the action, i.e., by decomposing the kinetic term of the scalar field in the Lagrangian density as
\be
\mathcal{L}_\mathrm{kin} = -\frac{\beta}{2} (\partial_\mu \vartheta)(\partial^\mu \vartheta) = \mathcal{L}_\mathrm{kin,1} + \mathcal{L}_\mathrm{kin,2} + \mathcal{L}_\mathrm{int}\,,
\ee
where
\ba
\mathcal{L}_\mathrm{kin,A} & \equiv & -\frac{\beta}{2} (\partial_\mu \vartheta_A)(\partial^\mu \vartheta_A)\,, \\
\label{eq:L_int}
\mathcal{L}_\mathrm{kin,int} & \equiv & -\beta (\partial_\mu \vartheta_1)(\partial^\mu \vartheta_2)\,,
\ea
and $\mathcal{L}_\mathrm{kin,int}$ corresponds to the dipole interaction Lagrangian density. By substituting $\vartheta_A = \mu_A^i n_{A,i}/r_A^2 = - \mu_A^i \partial_i (1/r_A)$ into Eq.~\eqref{eq:L_int}, we obtain
\ba
\mathcal{L}_\mathrm{kin,int} &=& -\beta \mu_1^j \mu_2^k \partial_{ij} \left( \frac{1}{r_1} \right) \partial_{ik} \left( \frac{1}{r_2} \right) \nonumber \\
&=& -\beta \mu_1^j \mu_2^k \partial^{(1)}_{ij} \left( \frac{1}{r_1} \right) \partial^{(2)}_{ik} \left( \frac{1}{r_2} \right)\,,
\ea
where $\partial^{(A)}_{i}$ denotes the partial derivative with respect to $x_A^i$.

The dipole interaction Lagrangian can be obtained by performing the volume integral of $\mathcal{L}_\mathrm{kin,int}$ as
\begin{align}
L_\mathrm{int} &= -\beta \mu_1^j \mu_2^k \int \partial^{(1)}_{ij} \left( \frac{1}{r_1} \right) \partial^{(2)}_{ik} \left( \frac{1}{r_2} \right) d^3x \nonumber \\
&= -\beta \mu_1^j \mu_2^k \partial^{(1)}_{ij} \partial^{(2)}_{ik} \int \frac{1}{r_1} \frac{1}{r_2} d^3x\,.
\end{align}
By applying Hadamard regularization, as explained e.g.~in Appendix B of~\cite{quadratic}, the above integration can be performed to yield
\begin{align}
L_\mathrm{int} &= 2 \pi \beta \mu_1^j \mu_2^k \partial^{(1)}_{ij} \partial^{(2)}_{ik} r_{12} = 12 \pi \beta \mu_1^j \mu_2^k \frac{n_{12\left\langle jk \right\rangle}}{r_{12}^3} \nonumber \\
&= 4\pi\beta \frac{1}{r_{12}^{3}} \left[ 3(\mu_{1}\cdot n_{12})
  (\mu_{2}\cdot n_{12}) - (\mu_{1}\cdot\mu_{2})  \right] \,.
\end{align}
As there are no derivatives on particle locations, this gives an
effective pairwise interaction potential
$$U_{\mathrm{int}} = -L_{\mathrm{int}}.$$
This expression agrees with Eq.~(5) of~\cite{kent-CSGW}. We can find the force with (minus) the particle derivative of the effective pairwise interaction potential as
\begin{equation}
F^{(\vartheta)}_{A,i} = -\partial^{(A)}_{i} U_{\mathrm{int}} \,.
\end{equation}

From the above we can compute the relative dipole-dipole force
(with $\mu$ the reduced mass)
\begin{align}
a^{(\vartheta)}_{12,i} &= \frac{75}{128} \frac{1}{\mu}\frac{\xi_{\mbox{\tiny CS}}}{m^4} \chi_{1} \chi_{2}
C_1^3 C_2^3 \bar{\mu}_1 \bar{\mu}_2 \left( \frac{m}{r_{12}} \right)^{4}
 \left\{\hat{S}_{1}^{i}(n_{12}\cdot\hat{S}_{2})
 \right.
 \nonumber \\
 &+ \left. \frac{1}{2}
n_{12}^{i} \left[   (\hat{S}_{1}\cdot \hat{S}_{2}) - 5 (n_{12}\cdot
  \hat{S}_{1})(n_{12}\cdot \hat{S}_{2}) \right]
\right\} - (1 \leftrightarrow 2)\,.
\end{align}
Decomposing $a_{12}^{(\vartheta)}$, inserting this into the Gauss
equations, and averaging, we find
\begin{widetext}
\begin{subequations}
\begin{align}
\ave{\dot{\omega}}_\vartheta &= \frac{75}{256} \frac{1}{\mu}\frac{\xi_{\mbox{\tiny CS}}}{m^4} \frac{\chi_{1}\chi_{2}}{(1-e^{2})^{2}} C_1^3 C_2^3 \bar{\mu}_1 \bar{\mu}_2 \left( \frac{m}{a} \right)^{7/2}
\left\{
\frac{1}{2}\left(\hat{S}_{1,x}\hat{S}_{2,x}
 +\hat{S}_{1,y}\hat{S}_{2,y}\right)
-\hat{S}_{1,z}\hat{S}_{2,z} \right.\nonumber\\
&\qquad\left.{}-\cot\iota \left[
\hat{S}_{1,z}\left(
\hat{S}_{2,x}\sin\omega + \hat{S}_{2,y}\cos\omega
\right)
\right]
\right\}+(1\leftrightarrow 2)\,, \\
\ave{{\textstyle \frac{d}{dt}}\iota}_\vartheta &= \frac{75}{256} \frac{1}{\mu}\frac{\xi_{\mbox{\tiny CS}}}{m^4} \frac{\chi_{1}\chi_{2}}{(1-e^{2})^{2}} C_1^3 C_2^3 \bar{\mu}_1 \bar{\mu}_2 \left( \frac{m}{a} \right)^{7/2}
\left[ \hat{S}_{1,z} \left( \hat{S}_{2,x}\cos\omega -
    \hat{S}_{2,y}\sin\omega \right) \right] +(1\leftrightarrow 2)\,, \\
\ave{\dot{\Omega}}_\vartheta &= \frac{75}{256} \frac{1}{\mu}\frac{\xi_{\mbox{\tiny CS}}}{m^4} \frac{\chi_{1}\chi_{2}}{(1-e^{2})^{2}} C_1^3 C_2^3 \bar{\mu}_1 \bar{\mu}_2 \left( \frac{m}{a} \right)^{7/2}
\csc\iota
\left[ \hat{S}_{1,z} \left( \hat{S}_{2,x}\sin\omega +
    \hat{S}_{2,y}\cos\omega \right) \right] +(1\leftrightarrow 2) \,,
\end{align}
\end{subequations}
and $\langle \dot{e} \rangle_\vartheta = 0 = \langle \dot{a} \rangle_\vartheta$.
\end{widetext}

{\renewcommand{\arraystretch}{1.2}
\begin{table*}
\begin{centering}
\begin{tabular}{rcc@{\hspace{1.5em}}c}
\hline
\hline
\noalign{\smallskip}
Timing parameter &&   Pulsar A & Pulsar B  \\
\hline
Spin frequency &$f_\nu$ (Hz) & 44.054069392744(2) & 0.36056035506(1) \\
Orbital period &$P_{b}$ (day) & 0.10225156248(5) & --- \\
Eccentricity &$e$ & 0.0877775(9) & --- \\
Projected semimajor axis &$(a/c) \sin \iota$ (s) & 1.415032(1) &1.5161(16) \\
Longitude of periastron &$\omega$ ($^\circ$) & 87.0331(8) & $87.0331+180.0$ \\
\hline
Advance rate of periastron & $\langle \dot{\omega} \rangle$ ($^\circ$ /yr) & 16.89947(68) & [16.96(5)]\footnote{An independent parameter fit of $\langle \dot{\omega} \rangle$ for pulsar B is consistent with the more precise result of pulsar A~\cite{kramer-double-pulsar}.} \\
Orbital decay rate &$\dot{P}_{b}$ & $-1.252(17) \times 10^{-12}$ & --- \\
Gravitational redshift &$\gamma$ (ms) & 0.3856(26) & --- \\
Shapiro delay range &$r$ ($\mu$s) & 6.21(33) & --- \\
Shapiro delay shape &$s$ & 0.99974(-39,+16) & --- \\
\hline
Inclination &$\iota$ ($^\circ$) &\multicolumn{2}{c}{88.69(-76,+50)} \\
Mass function &($M_\odot$) & 0.29096571(87) & 0.3579(11) \\
Mass ratio && \multicolumn{2}{c}{1.0714(11)} \\
Mass ($M_\odot$) && 1.3381(7) & 1.2489(7) \\
\noalign{\smallskip}
\hline
\hline
\end{tabular}
\end{centering}
\caption{Timing parameters of the double binary pulsar PSR J0737-3039~\cite{kramer-double-pulsar}. Measurement uncertainties on last digits are shown in parentheses.}
\label{table:double-pulsar}
\end{table*}
}

\begin{table}
\begin{centering}
\begin{tabular}{cc}
\hline
\hline
\noalign{\smallskip}
EoS & $\langle \dot{\omega} \rangle_h/\langle \dot{\omega} \rangle_{\mbox{\tiny GR}}$  \\
\noalign{\smallskip}
\hline
\noalign{\smallskip}
APR & -$3.7 \times10^{-9} $ \\
SLy & -$5.1 \times10^{-8} $ \\
LS & -$3.2 \times10^{-9} $ \\
Shen & -$3.0 \times10^{-8} $ \\
\noalign{\smallskip}
\hline
\hline
\end{tabular}
\end{centering}
\caption{\label{table:ratio}
$\langle \dot{\omega} \rangle_h/\langle \dot{\omega} \rangle_{\mbox{\tiny GR}}$ for PSR J0737-3039 with various NS EoSs, where we set $\zeta = 1$ (this ratio is linearly proportional to $\zeta$).}
\end{table}

\section{Applications to NS Observations}
\label{sec:applications}

In this section, we apply the results derived in the previous sections to observed NS systems. We consider the recently found massive millisecond pulsar J1614-2230~\cite{1.97NS} and the double binary pulsar system PSR J0737-3039~\cite{burgay,lyne,kramer-double-pulsar}. In this section, we consider the maximum CS corrections allowed within the weak coupling approximation, i.e.~$\zeta =1$. Clearly, such a value for $\zeta$ violates the small coupling approximation. But as we shall see, even with such a large $\zeta$ value, which leads to a strong GR deviation, NS observations will still not be accurate enough to allow for a bound on the theory.

\subsection{Massive Millisecond Pulsar J1614-2230}

From measurements of the Shapiro time delay, the mass of the pulsar J1614-2230 has been determined to be $(1.97 \pm 0.04)M_\odot$~\cite{1.97NS}. As mentioned previously, the effect of a CS modification is to \emph{reduce} the magnitude of the observed mass relative to the GR expectation. In Fig.~\ref{fig:MR-CS}, we plotted the mass-radius relation in both GR and dynamical CS gravity. We have set the spin period of the NS to be the one observed for J1614-2230, $P_\mathrm{spin} = 3.1508076534271(6)$ms. When making these plots, we did not include $\mathcal{O}(\alpha'^0 \chi'^2)$ contribution, but we have checked that such contributions only affect the results to 0.2$\%$ at most. As expected, the maximum mass decreases as we increase the CS dimensionless coupling constant $\zeta$. Therefore, we can try to place bounds on the theory by requiring that the maximum observed mass be greater than 1.93$M_\odot$. Of course, such a method to test modified gravity theories is not new; it has been considered in e.g.~\cite{eling-AE-NS} for Einstein-Aether theory and~\cite{pani-EDGB-NS} for Einstein-Dilaton-Gauss-Bonnet theory.

The horizontal dashed line in Fig.~\ref{fig:MR-CS} is the lower bound on the mass of PSR J1614-2230. For the LS220 EoS, the CS maximum NS mass with $\zeta=1$ becomes less than 1.93$M_\odot$. This means that we can place a meaningful constraint on the theory \emph{assuming} that this is the correct EoS. However, the CS maximum NS mass with $\zeta=1$ with the APR, SLy and Shen EoSs are all above 1.93$M_\odot$. Therefore, whether a meaningful constraint of dynamical CS gravity can be placed depends strongly on the EoS. Since the EoS has not been observationally measured, this degeneracy prevents us from constraining dynamical CS gravity with mass-radius relations.

\subsection{Double Pulsar Binary PSR J0737-3039}

The observed timing parameters of PSR J0737-3039~\cite{burgay,lyne,kramer-double-pulsar} are summarized in Table~\ref{table:double-pulsar}. One can test GR from measurements of the post-Keplerian (PK) parameters~\cite{damour-taylor}: the (orbital averaged) advance rate of the periastron $\langle \dot{\omega} \rangle$, the orbital decay rate $\dot{P}$, gravitational redshift parameter $\gamma$, and the range $r$ and the shape $s$ of the Shapiro time delay.

Let us first look at the CS effect on $\langle \dot{\omega} \rangle$. Since the spin of the secondary pulsar is 100 times smaller than the primary one, we only consider CS corrections that depend on $\chi_1$, meaning we only keep $\langle \dot{\omega} \rangle_h$ and neglect $\langle \dot{\omega} \rangle_\vartheta$. By taking the ratio of the former to the GR expression for $\langle \dot{\omega} \rangle_{\mbox{\tiny GR}} = 3 (2 \pi/P)^{5/3} m^{2/3} (1-e^2)^{-1}$~\cite{stairs,TEGP}, we obtain
\begin{align}
\label{eq:dotomega-ratio}
\frac{\langle \dot{\omega} \rangle_h}{\langle \dot{\omega} \rangle_{\mbox{\tiny GR}}} &= \frac{Q_1}{m^3} \frac{1}{1-e^2} \frac{m}{a} \left[ -1 + \frac{3}{2}\left(
    \hat{S}_{1,x}^{2}+\hat{S}_{1,y}^{2} \right)\right. \nonumber \\
& \quad\left.  -\hat{S}_{1,z}\cot\iota \left( \hat{S}_{1,x}\sin\omega +
    \hat{S}_{1,y}\cos\omega \right)
 \right]\,.
\end{align}
Since this is proportional to $m/a$, the CS correction to $\langle \dot{\omega} \rangle$ is of 1PN order relative to GR.

$\dot{P}$ is proportional to $\dot{E}$ which is given in Eqs.~\eqref{eq:dotEtheta} and~\eqref{eq:dotEh} and both are proportional to $(m/a)^7$. By comparing this with the GR expectation, which is proportional to $(m/a)^5$, the CS correction to $\dot{P}$ is of 2PN relative order. Similarly, one can show that the CS correction to $\gamma$, $r$, and $s$ are also of relative 2PN order. Therefore, the dominant CS correction to the binary evolution is in $\langle \dot{\omega} \rangle$ and we can neglect CS corrections to all other PK parameters. We can also neglect the CS corrections to the spin-precession equations, which appear at 2PN relative order (see Appendix~\ref{sec:spin-precession}).

In order to place constraints on dynamical CS gravity, one needs to know the orientation of the spin of the primary pulsar since Eq.~\eqref{eq:dotomega-ratio} depends on $\hat{S}_1^i$. Unfortunately, this quantity is currently unconstrained. In 2004, Ref.~\cite{jenet-ransom} attempted to measure this quantity by modeling the intensity variation of pulsar B caused by the emission from pulsar~A~\cite{jenet-ransom}. Doing so, one obtained $\theta_{S_1}=167^\circ \pm 10^\circ$ and  $\phi_{S_1} = 246^\circ \pm 5^\circ$, where $\theta_{S_1}$ is the angle between $\hat{S}_1^i$ and the unit orbital angular momentum vector $\hat{L}^i$, and $\phi_{S_1}$ is the angle between $\hat{S}_1^i$ projected onto the orbital plane and the direction of the ascending node. However, since then, Ref.~\cite{manchester} found no observable major profile changes in the emission of pulsar A, and hence the model of Ref.~\cite{jenet-ransom} is ruled out.  Nevertheless, following Ref.~\cite{lattimer-schutz}, we adopt the above values for $\hat{S}_1^i$ to give an order of magnitude estimate of the possible magnitude of the constraint, if this quantity were measured in the future.

By using Kepler's Law $a^3 = m (P/2 \pi)^2$ and the CS quadrupole moment deformation found in Sec.~\ref{sec:solutions}, we can evaluate Eq.~\eqref{eq:dotomega-ratio} for NSs with various EoSs. The results for $\zeta =1$ are summarized in Table~\ref{table:ratio}. One sees that in the small coupling approximation, the ratio of $|\langle \dot{\omega} \rangle_h / \langle \dot{\omega} \rangle_{\mbox{\tiny GR}}|$ is of $\mathcal{O}(10^{-8})$ at most. For other choice of the unit vector $\hat{S}_1^i$, $|\langle \dot{\omega} \rangle_h / \langle \dot{\omega} \rangle_{\mbox{\tiny GR}}|$ increases by at most 0.075\%. This means that in order to place meaningful constraints on dynamical CS gravity from binary pulsar observations, we need to observe $\langle \dot{\omega} \rangle$ \emph{and} at least 2 other PK parameters (or the mass functions or the mass ratio) within $10^{-8}$ accuracy (the latter are needed to determine the masses of the pulsars). With the current observational data, we can determine the masses of the pulsars using the mass ratio and $s$, and test dynamical CS gravity with $\langle \dot{\omega} \rangle$. However, since the former only have accuracies of $\mathcal{O}(10^{-3})$, we cannot place any constraint on this theory from the double binary pulsar system. This conclusion also holds for other binary pulsar systems.

Notice that $|\langle \dot{\omega} \rangle_h/\langle \dot{\omega} \rangle_{\mbox{\tiny GR}}|$ as shown in Table~\ref{table:ratio} is smaller than $|\langle \dot{\omega} \rangle_\mathrm{SO}/\langle \dot{\omega} \rangle_{\mbox{\tiny GR}}| \sim 10^{-5}$ where $\langle \dot{\omega} \rangle_\mathrm{SO}$ is the advance rate of the periastron due to GR spin-orbit coupling and it is given by~\cite{damour-schaefer,lattimer-schutz}
\ba
\frac{\langle \dot{\omega} \rangle_\mathrm{SO}}{\langle \dot{\omega} \rangle_{\mbox{\tiny GR}}} &=& -\frac{1}{6 \sqrt{1-e^2}} \frac{4 m_1 + 3 m_2}{m_1} \frac{I_1 f_{\nu_1}}{m^2} \left( \frac{m}{a} \right)^{1/2} \nonumber \\
& & \times (2 \cos \theta_{S_1} + \cot \iota \sin \theta_{S_1} \sin \phi_{S_1})\,,
\ea
where $f_{\nu_1}$ is the spin frequency of pulsar A. Since this ratio depends on the moment of inertia $I$, which encodes information of the internal structure of the NS, it would be difficult to test dynamical CS gravity even if the measurement accuracies of the PK parameters improved, unless $I$ were determined to high accuracy.

The results obtained here suggests that pulsar binaries cannot currently be used to constrain dynamical CS gravity. If so, GW observations are the only way to test dynamical CS gravity in the dynamical, strong-field regime, as demonstrated in~\cite{kent-CSGW}.

\section{Conclusions and Discussions}
\label{sec:conclusions}

In this paper, we extended the previous work in~\cite{yunespsaltis,alihaimoud-chen} to construct slowly-rotating NSs in dynamical CS gravity in the small-coupling and slowly-rotating approximations to quadratic order in spin. At this order, we found a negative CS mass correction and a positive CS quadrupole moment deformation. We applied the former correction to test the theory by requiring that the observed maximum mass be greater than the lower bound on the mass of the massive millisecond pulsar J1614-2230. Unfortunately, we could not obtain any meaningful constraint due to degeneracies with the unknown EoS of nuclear matter. Next, we used the quadrupole deformation to derive the correction to the evolution of the NS binary. Among all PK parameters, we found that the dominant CS correction appears in the advance rate of the periastron $\langle \dot{\omega} \rangle$. We applied our results to the double binary pulsar PSR J0737-3039~\cite{burgay,lyne,kramer-double-pulsar} but found that the CS correction is too small to be constrained, i.e.~it is of relative 1PN order. This correction is even smaller than the one arising from GR spin-orbit coupling, which in turn depends on the moment of inertia $I$. Hence, we would not be able to constrain the theory unless we knew $I$ to high precision. The results obtained here indicate that GW observations are the only way to test the theory in the strong and dynamical field regime~\cite{kent-CSGW}.

One possible way to extend our work is to construct gravitational waveforms for NS binaries. This can be done by following the analysis of Ref.~\cite{kent-CSGW}. As in the BH binary case discussed in the reference, the corrections to GR waveforms will enter at 2PN order. However, one needs to be careful about constructing gravitational waveforms for spinning NS binaries in GR, since the spin of the NS induces a quadrupole deformation. This deformation contains information about the internal structure of the NS, which also enters at 2PN order~\cite{poisson-quadrupole}. In order to obtain 2PN waveforms, one must thus first construct slowly-rotating NS solutions in GR to quadratic-order in spin. This spin-induced 2PN effect will be strongly degenerate with the 2PN CS correction. However, this degeneracy might be broken as follows. The \emph{tidally} induced quadrupole moment enters at 5PN order~\cite{flanagan-hinderer}, but is enhanced by a factor of $\mathcal{O}(\mathcal{R}_*/M_*)^5$. The tidal effect may be detected using future ground-based GW interferometers~\cite{flanagan-hinderer,read-markakis-shibata,hinderer-lackey-lang-read,lackey-kyutoku,damour-nagar-villain}. Since the internal structure can be determined to some extent from tides, the degeneracy between the spin-induced term and the CS correction might be broken. There is also a GR spin-spin interaction term at 2PN order in the phase of the gravitational waveform. Again, the degeneracies between spins and the CS correction term might be broken if the binaries are precessing. However, we expect that BH binaries are more useful in constraining the theory than NS binaries because (i) the degeneracy with the spin-induced term would not be broken completely, (ii) the CS deformation to the quadrupole moment for a NS is smaller than that of a stellar-mass BH due to cancellation between the first and the second terms in Eq.~\eqref{eq:QA}, and (iii) the spin of the NS binaries just before coalescence is expected to be rather small~\cite{bildsten-cutler,hinderer-lackey-lang-read,damour-nagar-villain}, making the CS correction even smaller. Also, the radius of a NS is larger than that of a stellar-mass BH, and hence, the latter allows tests at smaller length scales, which means that BH binaries have a stronger potential to constrain the theory.

Another possible avenue of future work is to study oscillations of NSs in dynamical CS gravity. By considering perturbations of the gravitational, scalar and the matter fields, one can investigate e.g. f-, p-, r- and w-modes~\cite{kokkotas-living}. The first two modes can be studied under the Cowling approximation, where one neglects the perturbations of the gravitational and the scalar fields (f- and p-modes in alternative theories of gravity have been studied in e.g.~\cite{sotani-fp-ST} for scalar-tensor theory and~\cite{sotani-fp-TeVeS} for tensor-vector-scalar (TeVeS) theory). This amounts to solving the perturbed matter equation $\nabla^{\mu} \delta T_{\mu\nu}^\mathrm{mat} =0$, where $\delta T_{\mu\nu}^\mathrm{mat}$ is the perturbed matter stress-energy tensor. At zeroth order in spin, since the background is exactly the same as GR, we expect that f- and p-modes to be identical to those in GR. To study w-modes, one needs to include gravitational and scalar field perturbations. W-modes in alternative theories of gravity have been investigated in e.g.~\cite{sotani-w-ST} for scalar-tensor theory and~\cite{sotani-w-TeVeS} for TeVeS. In these theories, one finds that axial metric perturbations decouple from scalar and polar metric ones. Hence, w-modes can be studied by solving the axial perturbation equation as an eigenvalue problem. In dynamical CS gravity, we expect that polar perturbations will decouple from the rest, as in the case for BHs~\cite{cardoso-gualtieri,garfinkle,molina}. However, again, since the background is the same as in GR, the polar perturbation equations should be identical to those in GR. At linear order in spin, one can look at r-modes by studying the toroidal oscillation of the NS in the Cowling approximation. As in GR~\cite{kojima,beyer} and in TeVeS~\cite{sotani-toroidal}, the spectrum for toroidal oscillations of slowly-rotating stars should also be continuous in dynamical CS gravity. However, this might be an artifact of the slow-rotation limit, since in GR, the spectrum becomes discrete for fast-rotating NSs~\cite{ruoff,gaertig}. In principle, if the maximum and/or the minimum of the toroidal oscillation frequencies and the moment of inertia are independently determined, one can break the degeneracy between the gravitational theory and the EoS, and hence test GR~\cite{sotani-toroidal}. This looks challenging from an observational point of view.  At quadratic order in spin, there can be corrections to f-, p-, and w-modes, but obviously, they are suppressed by $\mathcal{O}(\zeta \chi^2)$.

\acknowledgments
We would like to thank Will Farr for pointing out the Riccati method.
NY acknowledges support from NSF grant PHY-1114374, as well as support provided by the National Aeronautics and Space Administration from grant NNX11AI49G, under sub-award 00001944. TT is supported by the Grant-in-Aid for Scientific Research (Nos. 21244033, 21111006, 24103006 and 24103001). LCS acknowledges that support for this work was provided by the National Aeronautics and Space Administration through Einstein Postdoctoral Fellowship Award Number PF2-130101 issued by the Chandra X-ray Observatory Center, which is operated by the Smithsonian Astrophysical Observatory for and on behalf of the National Aeronautics Space Administration under contract NAS8-03060. Some calculations used the computer algebra-systems \textsc{Maple}, in combination with the \textsc{GRTensorII} package~\cite{grtensor}. Other calculations were carried out with the \textsc{xTensor} package for \textsc{Mathematica}~\cite{2008CoPhC.179..586M,2009GReGr..41.2415B}.

redpajama_set_name: RedPajamaArXiv