\section{Introduction}

The Kerr geometry, which describes the spacetime outside an isolated rotating black hole according to general relativity, possesses a  \emph{hidden symmetry}, not corresponding to any spacetime isometry, responsible for the integrability of geodesic motion and for the separability of various field equations. The spacetime possesses only two independent Killing vector fields, $\xi^\mu$ and $\eta^\mu$, generating the isometries of time-translation and rotation about the spin axis, leading to the energy $E=-\xi^\mu p_\mu$ and the axial component of angular momentum $L_z=\eta^\mu p_\mu$, respectively, as conserved quantities for geodesic motion   with momentum $p_\mu$, \emph{i.e.} $p^\nu\nabla_\nu p_\mu=0$.   It was thus rather unexpected when Carter \cite{PhysRev.174.1559} found a further nontrivial constant of the motion which is quadratic in the momentum, of the form $Q=K^{\mu\nu}p_\mu p_\nu$, where $K^{\mu\nu}=K^{\nu\mu}$ is a Killing tensor satisfying $\nabla_{(\mu}K_{\nu)\rho}=0$, a generalization of the Killing equation $\nabla_{(\mu}\xi_{\nu)}=0$.   Adding the Carter constant $Q$ to $E$, $L_z$ and $m^2=-g^{\mu\nu}p_\mu p_\nu$ (conserved for geodesics in any background), there are four Poisson-commuting constants of motion, ensuring that geodesic motion in the Kerr spacetime is a fully integrable dynamical system.   This implies the absence of chaos in the motion, and it reduces the problem of solving for geodesic trajectories to the task of performing one-dimension integrals.

While the geodesic equation describes the motion of a structureless monopolar test body in a fixed background spacetime, an important generalization is to allow the test body (while still having negligible mass and thus negligible influence on the gravitational field) to have a finite size and nontrivial structure.   In the case where such an ``extended test body'' has a size (length scale) $l$ which is small compared to the radius of curvature $R$ of the background, $l\ll R$, it is usefully characterized by a centroid worldline $x=z(\tau)$, with tangent $v^\mu\triangleq \dd z^\mu/\dd\tau$, and a tower of gravitational multipole moments defined along the worldline.   These begin with the momentum $p_\mu$ as the monopole and the spin (relativistic angular momentum) tensor $S_{\mu\nu}=-S_{\nu\mu}$ as the dipole.   Using only the fact that the body is described by a test stress-energy tensor $T_{\mu\nu}$ which is conserved within the background, $\nabla_\mu T^{\mu\nu}=0$, and certain definitions of the multipole moments as spatial integrals of $T^{\mu\nu}$ (which reduce locally to standard definitions in special relativity), one can show that the monopole $p_\mu$ and dipole $S_{\mu\nu}$ must evolve along the worldline according to the Mathisson-Papapetrou-Dixon (MPD) equations \cite{Mathisson:1937zz,1951RSPSA.209..248P,1970RSPSA.314..499D},
\begin{subequations}\label{MPD_intro}
\begin{align}
      \frac{\text{D}p^\mu}{\dd\tau}&=-\frac{1}{2}R\tud{\mu}{\nu\alpha\beta}v^\nu S^{\alpha\beta}+\ldots,\\
      \frac{\text{D}S^{\mu\nu}}{\dd\tau}&=2 p^{[\mu}v^{\nu]}+\ldots,
\end{align}
\end{subequations}
reflecting local conservation of momentum and angular momentum, where the dots represent corrections due to the quadrupole and higher multipole moments.   These are to be supplemented by a condition of the form $w^\mu S_{\mu\nu}=0$ for some timelike vector field $w^\mu$ (setting to zero the mass-dipole moment in the frame of $w^\mu$), which fixes a choice of centroid worldline; a convenient choice is the Tulczyjew-Dixon condition $p^\mu S_{\mu\nu}=0$ \cite{Tulczyjew:1957aa,Tulczyjew:1959aa}.   Along with such a condition, the pole-dipole MPD equations (Eqs.~(\ref{MPD_intro}) with the dots dropped, neglecting the quadrupole and higher corrections) provide a closed set of evolution equations for the worldline $z(\tau)$ and the momentum $p^\mu(\tau)$ and spin $S^{\mu\nu}(\tau)$ along it, describing the motion of spinning test body a background curved spacetime.

As shown by Dixon \cite{1970RSPSA.314..499D,1977GReGr...8..197E}, and as was central to his construction of the multipole moments, if the background has a Killing vector $\xi^\mu$, then the quantity
\begin{equation}
\mathcal P_\xi=p_\mu\xi^\mu+\frac{1}{2}S^{\mu\nu}\nabla_\mu\xi_\nu
\end{equation}
is exactly conserved along any worldline when $p^\mu$ and $S^{\mu\nu}$ are evolved by the MPD equations (\ref{MPD_intro}), to all orders in the multipole expansion, for arbitrary quadrupole and higher moments.  (See also, e.g., the earlier derivation by Souriau for the pole-dipole system \cite{Souriau:1974}, and the insightful exposition by Harte \cite{Harte:2014wya}.) Considering a background Kerr spacetime, one is then naturally lead to wonder whether the hidden symmetry leads to conserved quantities for the MPD dynamics, including a generalization of the Carter constant to the case of spinning extended test bodies.

This question was answered for the case of the pole-dipole MPD equations by R\"udiger \cite{doi:10.1098/rspa.1981.0046,doi:10.1098/rspa.1983.0012}. First he showed that the quantity
$
Q_Y={}^*Y^{\mu\nu}S_{\mu\nu}
$
built from the Killing-Yano tensor $Y_{\mu\nu}$ of Kerr is conserved, up to remainders quadratic in the spin tensor and quadrupolar corrections. He further showed that there is indeed a generalized Carter constant of the form
\begin{equation}\label{Q_dipole_intro}
Q^{(2)}=K^{\mu\nu}p_\mu p_\nu+L^{\mu\nu\rho}S_{\mu\nu}p_\rho +\ldots,
\end{equation}
which is conserved along (\ref{MPD_intro}) with $p^\mu S_{\mu\nu}=0$, up to remainders quadratic in the spin tensor and quadrupolar corrections that were not determined. A generalization of this result to the Kerr-Newman (charged spinning black hole) spacetime was independently discovered by Gibbons \emph{et al.}~\cite{Gibbons:1993ap} using a supersymmetric description of spinning particle dynamics.    The existence of these constants of motion has been shown by Witzany to allow the separation of a Hamilton-Jacobi equation for the pole-dipole system in Kerr, leading to analytic expressions for the fundamental frequencies of the motion \cite{Witzany:2019nml} using a Hamiltonian formalism for spinning particles  \cite{Witzany:2018ahb}.

Our purpose in this paper is to explore whether such ``hidden constants'' exist for test bodies with spin-induced quadrupole moments moving in a Kerr background.  Dixon's generalizations of the equations of motion (\ref{MPD_intro}) to the quadrupolar order in the test body's multipole expansion read
\begin{subequations}\label{MPD_quad_intro}
\begin{align}
      \frac{\text{D}p^\mu}{\dd\tau}&=-\frac{1}{2}R\tud{\mu}{\nu\alpha\beta}v^\nu S^{\alpha\beta}-\frac{1}{6}J^{\alpha\beta\gamma\delta}\nabla^\mu R_{\alpha\beta\gamma\delta}+\ldots,
      \\
      \frac{\text{D}S^{\mu\nu}}{\dd\tau}&=2 p^{[\mu}v^{\nu]}+\frac{4}{3}R\tud{[\mu}{\alpha\beta\gamma}J^{\nu]\alpha\beta\gamma}+\ldots,
\end{align}
\end{subequations}
where $J^{\mu\nu\rho\sigma}(\tau)$ is the quadrupole tensor, having the same symmetries as the Riemann tensor, and the ellipses here represent octupolar and higher corrections.  As has been developed and applied in a number of works (see e.g.\ \cite{Porto:2005ac,Porto:2008jj,Steinhoff:2011sya,Steinhoff:2014kwa,Marsat_2015,Bohe:2015ana,Levi:2015msa,Bini:2015zya,Vines:2016unv}), the form of $J$ appropriate to describe a spin-induced quadrupole moment is given by
\begin{equation}\label{quad_intro}
J^{\mu\nu\rho\sigma}=\kappa\frac{3 p\cdot v}{(p^2)^2}p^{[\mu}S^{\nu]\lambda}S^{[\rho}{}_{\lambda}p^{\sigma]},
\end{equation}
where $\kappa$ is a response coefficient controlling the magnitude of the quadrupolar deformation, proportional to the square of the spin.  Typical values for $\kappa$ for a neutron star are in the range 4 to 8 \cite{Laarakkers:1997hb}, while for a black hole $\kappa_\text{BH}=1$.

As the central results of this paper, we establish that two quantities, $Q_Y$ and $Q^{(2)}$, are conserved up to cubic-in-spin or octupolar corrections,
\begin{equation}
\frac{\dd Q_Y}{\dd\tau}=\mathcal O(\mathcal S^3),
\qquad
\frac{\dd Q^{(2)}}{\dd\tau}=\mathcal O(\mathcal S^3),
\end{equation}
along the motion of a ``quadrupolar test black hole'', governed by (\ref{MPD_quad_intro})--(\ref{quad_intro}) with $\kappa=1$ and with the Tulczyjew condition $p_\mu S^{\mu\nu}=0$, in a background Kerr spacetime, for arbitrary orbital and spin orientations.  The first quantity $Q_Y$ is R\"udiger's linear-in-spin constant, unmodified,
\begin{equation}
Q_Y={}^*Y^{\mu\nu}S_{\mu\nu},
\end{equation}
where ${}^*Y^{\mu\nu}=\tfrac{1}{2}\epsilon^{\mu\nu\alpha\beta}Y_{\alpha\beta}$ is the dual of the Kerr spacetime's nontrivial Killing-Yano tensor $Y^{\mu\nu}$.  The second quantity $Q^{(2)}$ is quadratic in $p$ and $S$ and generalizes R\"udiger's constant (\ref{Q_dipole_intro}) to the quadrupolar order for a test black hole; it is given explicitly by
\begin{alignat}{3}
Q^{(2)}&=
Y_{{\mu}\rho}Y^\rho_{\;\; \nu}p^{\mu} p^{\nu}+4 \xi^\lambda\epsilon_{\lambda \mu\sigma[\rho}Y_{\nu]}{}^\sigma S^{\mu\nu}p^\rho \nonumber
\\
&\quad-\bigg[
g_{\mu\rho}\Big(\xi_\nu \xi_\sigma-\frac{1}{2}g_{\nu\sigma}\xi^2\Big)-\frac{1}{2}Y_\mu{}^\lambda\Big(Y_\rho{}^\kappa R_{\lambda\nu\kappa\sigma}+\frac{1}{2}Y_\lambda{}^\kappa R_{\kappa\nu\rho\sigma}\Big)
\bigg]S^{\mu\nu}S^{\rho\sigma},
\end{alignat}
where $\xi^\mu$ is the timelike Killing vector, and it reduces to the Carter constant $K_{\mu\nu}p^\mu p^\nu=Y_{{\mu}\rho}Y^\rho_{\;\; \nu}p^{\mu} p^{\nu}$ for geodesic motion when the test body's spin is set to zero.

After reviewing details of the motion of quadrupolar test bodies in curved spacetime in Sec.~\ref{sec:QTB}, we develop the constraint equations for the existence of the conserved quantities in tensorial form in Sec.~\ref{sec:ConstrainEqs}.
In Sec.~\ref{sec:CBB}, we discuss covariant algebraic and differential relations amongst basic fields, ``covariant building blocks,'' characterizing the Kerr geometry, which play a central role in our solutions to the constraint equations.  We use these to reduce the tensorial constraint equations to a system of scalar equations in Sec.~\ref{sec:Reduction}, and we derive our solutions for the special black-hole case $\kappa=1$ in Sec.~\ref{sec:BH}.  In Sec.~\ref{sec:NS}, we investigate the case $\kappa\ne 1$, for non-black-hole bodies such as neutron stars with spin-induced quadrupoles, concluding that there is no solution to the constraints for $\kappa\ne 1$.  Finally we summarize our findings, some aspects of their broader context, and future directions in Sec.~\ref{sec:outlooks}.

We use the same conventions as adopted previously in \cite{Compere:2021kjz}. The main results of this paper have been crosschecked using \textit{Mathematica}, which was also used to performed the numerical evaluations of Section \ref{sec:NS}. The various notebooks can be found on the GitHub repository \url{https://github.com/addruart/generalizedCarterConstant}.

\section{Quadrupolar test bodies in curved spacetime}\label{sec:QTB}

In this section we review the motion of test-bodies in curved spacetime and we discuss the problem of finding conserved quantities associated to the corresponding dynamical system.

\subsection{Motion of test bodies in curved spacetime}

\paragraph{Evolution equations.} In General Relativity, the motion of an extended test body over a curved background is described by the Mathisson-Papapetrou-Dixon (MPD) equations \cite{Mathisson:1937zz,1951RSPSA.209..248P,1970RSPSA.314..499D}:
\begin{subequations}\label{MPD}
\begin{align}
    \frac{\text{D}p^\mu}{\dd\tau}&=-\frac{1}{2}R\tud{\mu}{\nu\alpha\beta}v^\nu S^{\alpha\beta}+\mathcal F^\mu,\label{MPD_1}\\
    \frac{\text{D}S^{\mu\nu}}{\dd\tau}&=2 p^{[\mu}v^{\nu]}+\mathcal L^{\mu\nu}.\label{MPD_2}
\end{align}
\end{subequations}
Here, $v^\mu$ denotes the four-velocity of the test body, $p^\mu$ its four-impulsion and $S^{\mu\nu}$ its spin dipole tensor. $\mathcal F^\mu$ and $\mathcal L^{\mu\nu}$ are respectively the force and torque terms that include corrections to the equations of motion due to quadrupole and higher moments. From now on, we assume that $\tau$ is the proper time, yielding $v_\mu v^\mu=-1$.

Taken alone, the MPD equations do not form a closed set of equations. Roughly speaking, this comes from the fact that the MPD equations arise from the skeletonization of a compact body stress-energy tensor above an arbitrary worldline belonging to the body's worldtube \cite{Mathisson:1937zz}. They shall be supplemented by a so-called \textit{spin supplementary condition} (SSC), which will play the role of fixing this worldline \cite{dixon1979}. In this paper, we choose to work with Tulczyjew SSC  \cite{Tulczyjew:1957aa,Tulczyjew:1959aa}
\begin{align}
    p_\mu S^{\mu\nu}=0.\label{SSC}
\end{align}
Among other consequences extensively described in \cite{Compere:2021kjz}, enforcing this SSC allows in particular to express the dipole spin tensor solely in terms of a spin vector $S^\mu$ which is orthogonal to the impulsion, $S_\mu p^\mu=0$.

\paragraph{Spin-induced quadrupole approximation.}
We will consider only spin-induced multipole moments and work in the quadrupole approximation, \textit{i.e.} neglecting octupole and higher moments. This is the relevant approximation for addressing spin-squared interactions: considering only spin-induced multipole terms, the $2^n$-pole scales as $\mathcal O(\mathcal S^n)$, with $\mathcal S^2\triangleq\frac{1}{2}S_{\mu\nu}S^{\mu\nu}$.

At the level of the equations of motion, this corresponds to choose the force and torque given by
\begin{align}
    \mathcal F^\mu&=-\frac{1}{6}J^{\alpha\beta\gamma\delta}\nabla^\mu R_{\alpha\beta\gamma\delta},\qquad
    \mathcal L^{\mu\nu}=\frac{4}{3}R\tud{[\mu}{\alpha\beta\gamma}J^{\nu]\alpha\beta\gamma}.
\end{align}
The quadrupole tensor $J^{\mu\nu\rho\sigma}$ possesses the same algebraic symmetries as the Riemann tensor.

We further particularize our study by considering only a quadrupole moment that is induced by the spin of the body, discarding the possible presence of some intrinsic quadrupole moment. This \textit{spin-induced quadrupole} can be shown to take the form \eqref{quad_intro}  \cite{Porto:2005ac,Porto:2008jj,Steinhoff:2011sya,Steinhoff:2014kwa,Marsat_2015,Bohe:2015ana,Levi:2015msa,Bini:2015zya,Vines:2016unv}. Here specialized to the case $v_\alpha v^\alpha=-1$ and at leading order in the spin expansion using Eq. \eqref{p_of_v}, its form reduces to
\begin{align}
    J^{\mu\nu\rho\sigma}=\frac{3\kappa}{\mu}v^{[\mu}S^{\nu]\lambda}S\tdu{\lambda}{[\rho}v^{\sigma]} = -\frac{3\kappa}{\mu} v^{[\mu}\Theta^{\nu][\rho}v^{\sigma]},\qquad\text{where }\Theta^{\alpha\beta}\triangleq S^{\alpha\lambda}S\tud{\beta}{\lambda}.\label{spin_induced_Q}
\end{align}
Here $\kappa$ is a free coupling parameter that equals 1 for a Kerr black hole and takes another value if the test-body is another compact object, \textit{e.g.} a neutron star.

Under the Tulczyjew SSC, the spin tensor can be solely expressed in terms of a spin vector defined as $S^\alpha\triangleq\frac{1}{2}\epsilon^{\alpha\beta\gamma\delta}\hat p_\beta S_{\gamma\delta}$. This relation can be inverted as
\begin{equation}
    S^{\alpha\beta}=2S^{[\alpha}\hat p^{\beta]*},
\end{equation}
were we have introduced the Hodge duality $A^*_{\mu\nu}\triangleq\frac{1}{2}\epsilon_{\mu\nu\rho\sigma}A^{\rho\sigma}$ (which is here specialized to the outer product of vectors, $l^{[\mu}m^{\nu]*}\triangleq\frac{1}{2}\epsilon^{\mu\nu\rho\sigma}l_\rho m_\sigma$).
One can show that this implies the following decomposition for $\Theta^{\alpha\beta}$:
\begin{align}
    \Theta^{\alpha\beta}&=\Pi^{\alpha\beta}\mathcal S^2-S^\alpha S^\beta
\end{align}
with $\Pi^\alpha_\beta\triangleq\delta^\alpha_\beta+\hat p^\alpha\hat p_\beta$ the projector on the hypersurface orthogonal to the timelike unit vector $\hat p^\alpha\triangleq\frac{p^\alpha}{\sqrt{-p_\alpha p^\alpha}}$ and $\mathcal S^2=\frac{1}{2}S_{\alpha\beta}S^{\alpha\beta}=S_\alpha S^\alpha$. Moreover, one has the identities
\begin{subequations}
\begin{align}
    \mathcal F^\mu&=\frac{\kappa}{2\mu}\hat p^\alpha\Theta^{\beta\gamma}\hat p^\delta\nabla^\mu R_{\alpha\beta\gamma\delta}\label{F}+\mathcal O(\mathcal S^3) ,\\
    \mathcal L^{\mu\nu}&=\frac{2\kappa}{\mu}R\tud{\nu}{\alpha\beta\gamma}v^{[\mu}\Theta^{\alpha]\beta}v^\gamma-\qty(\mu\leftrightarrow\nu),\label{L}\\
    \mathcal L^{\mu\nu}v_\nu&=\frac{\kappa}{\mu}\qty(\hat p^\mu \hat p^\nu R_{\nu\alpha\beta\gamma}+R\tud{\mu}{\alpha\beta\gamma})\Theta^{\alpha\beta} \hat p^\gamma+\mathcal O(\mathcal S^4)
= \frac{\kappa}{\mu} \Pi^{\mu\nu} R_{\nu \alpha\beta\gamma}\Theta^{\alpha\beta}\hat p^\gamma +\mathcal O(\mathcal S^4).\label{Lv}
\end{align}
\end{subequations}

\paragraph{Conservation of the spin, mass; relation between four-velocity and impulsion.} We define the invariant and kinematic masses, respectively given by
\begin{align}
    \mu^2\triangleq-p_\alpha p^\alpha,\qquad \mathfrak m\triangleq -v_\alpha p^\alpha.
\end{align}
Differentiating the SSC \eqref{SSC} yields
\begin{align}
    \mu^2 v^\mu-\mathfrak m p^\mu=\frac{1}{2}S^{\mu\nu}R_{\nu\lambda\rho\sigma}v^\lambda S^{\rho\sigma}-\mathcal L^{\mu\nu}p_\nu- S^{\mu\nu}\mathcal F_\nu.
\end{align}
Contracting this equation with $v_\mu$ provides us with
\begin{align}
    \mu^2=\mathfrak m^2+\mathcal O(\mathcal S^3).\label{masses}
\end{align}
The expression of the 4-impulsion in terms of the 4-velocity reads as \begin{align}
    p^\mu=\mu v^\mu-\frac{1}{2\mu}S^{\mu\nu}R_{\nu\lambda\rho\sigma}v^\lambda S^{\rho\sigma}+\mathcal L^{\mu\nu}v_\nu+\mathcal O(\mathcal S^3).
\label{p_of_v}
\end{align}
In the quadrupole approximation, $\mu$ is no longer conserved at $\mathcal O(\mathcal S^2)$, since
\begin{align}
    \dv{\mu}{\tau}=-v_\mu\mathcal F^\mu+\mathcal O(\mathcal S^3).
\end{align}
However, notice that, provided we assume\footnote{This condition is automatically satisfied for the spin-induced quadrupole.}
\begin{align}
    \frac{\text{D}}{\dd\tau}J^{\alpha\beta\gamma\delta}=\mathcal O(\mathcal S^3),
\end{align}
one can still define a mass-like quantity, given by
\begin{align}
    \tilde\mu\triangleq \mu-\frac{1}{6}J^{\alpha\beta\gamma\delta}R_{\alpha\beta\gamma\delta},
\end{align}
which is quasi-conserved, namely
\begin{align}
    \dv{\tilde\mu}{\tau}=\mathcal O(\mathcal S^3).
\end{align}
Moreover, one can perturbatively invert \eqref{p_of_v} to obtain an expression of the four-velocity in terms of the impulsion and the spin:
\begin{align}
    v^\mu = \hat p^\mu +(D^{\mu\nu}-\frac{1}{\mu}\mathcal L^{\mu\nu})\hat p_\nu+\mathcal O(\mathcal S^3),
    \label{v_of_p}
\end{align}
with
\begin{align}
    \hat p^\alpha\triangleq\frac{p^\alpha}{\mu}=\frac{p^\alpha}{\tilde\mu}+\mathcal O(\mathcal S^2),\qquad D\tud{\mu}{\nu}&\triangleq\frac{1}{2\mu^2}S^{\mu\lambda}R_{\lambda\nu\rho\sigma}S^{\rho\sigma}.
\end{align}
Eq. \eqref{v_of_p} will play a central role when we will work out the conservation equations in the following sections.

The spin magnitude $\mathcal S^2$ is exactly conserved \cite{Bailey:1975fe}
\begin{equation}
\frac{d}{d\tau}(\mathcal S^2) = 0.
\end{equation}

\subsection{Searching for conserved quantities: Rüdiger's procedure}

In two papers published in the early 80s, Rüdiger described a scheme for constructing quantities conserved along the motion driven by the MPD equations \cite{doi:10.1098/rspa.1981.0046,doi:10.1098/rspa.1983.0012}. The basic guideline followed in his scheme was to enforce directly the conservation equation on a generic Ansatz for the conserved quantity, and to subsequently solve the constraints obtained. This procedure was extensively reviewed and discussed in \cite{Compere:2021kjz}. We provide here a short summary of its main steps, which would allow the reader to get familiar with our terminology and notations.

\begin{itemize}[label=$\diamond$]
    \item \textbf{Step 1: postulate an Ansatz for the conserved quantity.} The conserved quantity should be a function of the dynamical variables $p^\mu$ and $S^{\mu}$. It is therefore a function $Q\qty(x^\mu,S^{\mu},p^\mu)$. Assuming its analyticity, it can be expanded as
\begin{align}
    Q\qty(x^\mu,S^\alpha,p^\mu)=\sum_{\substack{s,p\geq0\\s+p>0}}Q^{[s,p]}\qty(x^\mu,S^\alpha,p^\mu)\label{ansatz_all_orders}
\end{align}
with
\begin{align}
    Q^{[s,p]}\qty(x^\mu,S^\alpha,p^\mu)
    \triangleq Q^{[s,p]}_{\alpha_1\ldots\alpha_{s}\mu_1\ldots\mu_p}(x^\mu)S^{\alpha_1}\ldots S^{\alpha_{s}}p^{\mu_1}\ldots p^{\mu_p}.\label{ansatz_term}
\end{align}
Expressions like this one -- that is, tensorial quantities fully contracted with occurrences of the impulsion and the spin -- will often appear in the following computations. It is useful to enable a distinction between them by introducing a grading allowing the counting of the number of occurrences of both the spin vector $S^{\mu}$ and the impulsion $p^\mu$, which is provided by the notation $[s,p]$. More generally, we define:

\begin{definition}
A fully-contracted expression of the type
\begin{align}
    T_{\alpha_1\ldots\alpha_s\mu_1\ldots\mu_p}\ell_s^{\alpha_1}\ldots\ell_s^{\alpha_s}\ell_p^{\mu_1}\ldots\ell_p^{\mu_p}\label{fully_contracted}
\end{align}
where $\ell_s^\alpha=S^\alpha,s^\alpha$ (the relaxed spin vector $s^\mu$ will be defined below) and $\ell_p^\mu=p^\mu,\hat p^\mu$ is said to be of grading $[s,p]$. Equivalently, $s$ (resp. $p$) will be referred to as the spin (resp. momentum) grading of this expression.
\end{definition}

Since we have only included the quadrupole term in the equations of motion but neglected all the $\mathcal{O}\qty(\mathcal S^3)$ terms, it is not self-consistent to look at quantities which are conserved beyond second order in the spin magnitude. We therefore restrict our analysis to Ansätze that contain terms of of spin grading at most equal to two. Historically, Rüdiger didn't consider the full set of possible Ansätze originating from this discussion, but only the two restricted cases
\begin{align}
    Q^{(1)}&\triangleq\sum_{p=1}Q^{[s,p]}\triangleq X_\mu p^\mu+W_{\mu\nu}S^{\mu\nu},\label{linear_ansatz}\\
    Q^{(2)}&\triangleq\sum_{p=2}Q^{[s,p]}\triangleq K_{\mu\nu}p^\mu p^\nu+L_{\mu\nu\rho}S^{\mu\nu}p^\rho+M_{\mu\nu\rho\sigma}S^{\mu\nu}S^{\rho\sigma}. \label{quadratic_ansatz}
\end{align}
We will refer to them as respectively the \textit{linear} and the \textit{quadratic} invariants in $p^\mu$. They are homogeneous in the number of occurrences of $p^\mu$ and $S^{\mu\nu}$ they contain. As long as we consider the MPD equations at linear order in the spin magnitude or at quadratic order with the quadrupole coupling of the test body being the one of a black hole ($\kappa=1$), it turns out that considering only these two types of ansatzes will be enough to derive a complete set of conserved quantities. However, a more general ansatz will be necessary to consider arbitrary quadrupole couplings ($\kappa\neq 1$), as discussed in Section \ref{sec:NS}.

\item \textbf{Step 2: write down the conservation equation.}
We only require our quantity to be conserved up to second order in the spin magnitude,
\begin{align}
    \dot Q\triangleq v^\lambda\nabla_\lambda Q\stackrel{!}{=}\mathcal O\qty(\mathcal S^3).
\end{align}

\item \textbf{Step 3: expand the conservation equation using the equations of motion.}
The next step is to plug the explicit form of the Ansatz chosen in the conservation equation, and to use the MPD equations \eqref{MPD} to replace the covariant derivatives of the impulsion and of the spin tensor. The occurrences of the four-velocity are replaced by the means of Eq. \eqref{v_of_p}.

\item \textbf{Step 4: express the conservation equation in terms of independent variables.}
The presence of a SSC make the variables $p^\mu$, $S^{\alpha\beta}$ not independent among themselves. We turn to an independent set of variables in two steps: (i) we use the relation $S^{\alpha\beta}=2S^{[\alpha}\hat p^{\beta]*}$ to replace all the spin tensors $S^{\mu\nu}$ by the spin vectors $S^{\mu}$ and (ii) we replace the occurrences of the spin vector by the \textit{relaxed spin vector} $s^\alpha$ defined through
\begin{align}
    S^\alpha=\Pi^\alpha_\beta s^\beta.
\end{align}
It allows to relax the residual constraint $S_\mu p^\mu=0$ by considering a spin vector possessing a non-vanishing component along the direction of the impulsion. Physical quantities will be independent of this component. It is introduced in order to decouple the conservation equation. For convenience, we scale the unphysical component of the relaxed spin vector such that $s_\alpha s^\alpha\sim S_\alpha S^\alpha=\mathcal S^2$. Notice that we have the useful identity
\begin{align}
    S^{[\alpha}p^{\beta]}=s^{[\alpha}p^{\beta]}\qquad\Rightarrow\qquad S^{\alpha\beta}=2s^{[\alpha}\hat p^{\beta]*}.
\end{align}

\item \textbf{Step 5: infer the independent constraints.} The conservation equation takes now the form of a sum of fully-contracted expressions of the type \eqref{fully_contracted}, involving only the \textit{independent} dynamical variables $p^\mu$ and $s^\alpha$:
\begin{align}
    \dot Q=\sum_{\substack{s,p\geq 0\\s+p>0}}T^{[s,p]}_{\alpha_1\ldots\alpha_s\mu_1\ldots\mu_p}s^{\alpha_1}\ldots s^{\alpha_s} p^{\mu_1}\ldots p^{\mu_p}\stackrel{!}{=}\mathcal O\qty(\mathcal S^3).
\end{align}
The conservation equation is then equivalent to the requirement that all the terms of different gradings $[s,p]$ vanish independently:
\begin{align}
    T^{[s,p]}_{\alpha_1\ldots\alpha_s\mu_1\ldots\mu_p}s^{\alpha_1}\ldots s^{\alpha_s} p^{\mu_1}\ldots p^{\mu_p}\stackrel{!}{=}\mathcal O\qty(\mathcal S^3).
\end{align}
$s^\alpha$ and $p^\mu$ being arbitrary, this is equivalent to the \emph{constraint equations}
\begin{align}
    T^{[s,p]}_{(\alpha_1\ldots\alpha_s)(\mu_1\ldots\mu_p)}\stackrel{!}{=}\mathcal O\qty(\mathcal S^3).\label{generic_tensor_cst}
\end{align}

\item \textbf{Step 6: find a solution and prove uniqueness.} This final step is non-systematic. For the simplest cases (linear invariant with black hole quadrupole coupling, quadratic invariant at first order in the spin magnitude), it will be sufficient to work only with the tensorial constraints \eqref{generic_tensor_cst}. However, for more involved cases (linear invariant with arbitrary quadrupole coupling, quadratic invariant at second order in the spin), the tensorial relations will become so cumbersome that turning to another formulation of the problem will appear to be fruitful. This will be the purpose of the covariant building blocks for Kerr introduced in Section \ref{sec:CBB}.
\end{itemize}

\section{Constraint equations: tensorial formulation}\label{sec:ConstrainEqs}
In this section, we will apply the aforementioned procedure to derive the tensorial constraint equations that should be obeyed for ensuring the conservation of the two quantities \eqref{linear_ansatz}, \eqref{quadratic_ansatz}.

\subsection{Linear constraint}
Following Rüdiger, we start from the Ansatz \eqref{linear_ansatz} for the linear invariant:
\begin{align}
    Q^{(1)}\triangleq X_\mu p^\mu+W_{\mu\nu}S^{\mu\nu}.
\end{align}
Notice that $W_{\mu\nu}$ should be  a skew-symmetric tensor. The time variation of \eqref{linear_ansatz} is given by
\begin{align}
    \dot Q^{(1)}&=v^\lambda\qty(\nabla_\lambda X_\mu p^\mu+X_\mu\nabla_\lambda p^\mu+\nabla_\lambda W_{\mu\nu}S^{\mu\nu}+W_{\mu\nu}\nabla_\lambda S^{\mu\nu}).\label{Q1}
\end{align}
Applying Rüdiger's procedure, the conservation equation $\dot Q^{(1)}=\mathcal O\qty(\mathcal S^3)$ reduces to the following set of equations:
\begin{subequations}\label{linear_cst}
\begin{align}
    [0,2]:&\qquad\nabla_\mu X_\nu\hat p^\mu\hat p^\nu=\mathcal O\qty(\mathcal S^3),\label{[0,2]}\\
    [1,2]:&\qquad\nabla_\mu Y_{\alpha\nu}s^\alpha\hat p^\mu\hat p^\nu-\frac{1}{2}X^\lambda R^*_{\lambda\nu\beta\rho}s^\beta\hat p^\nu\hat p^\rho=\mathcal O\qty(\mathcal S^3),\label{[1,2]}\\
    [2,2]:&\qquad\frac{\kappa}{2\mu}X^\lambda\nabla_\lambda R_{\nu\alpha\beta\rho}s^\alpha s^\beta\hat p^\nu\hat p^\rho+Y_{\mu\nu}\mathcal L^{*\mu\nu}=\mathcal O\qty(\mathcal S^3),\label{[2,2]}\\
    [2,4]:&\qquad\qty(\nabla_\lambda X_\mu-2 W_{\lambda\mu})\qty(\mu D\tud{\lambda}{\nu}-\mathcal L\tud{\lambda}{\nu})\hat p^\mu\hat p^\nu=\mathcal O\qty(\mathcal S^3).\label{[2,4]}
\end{align}
\end{subequations}
Here, we have introduced the notation $Y_{\mu\nu}\triangleq W^*_{\mu\nu}$.
We therefore have the following proposition:
\begin{prop}
For any pair $(X_\mu,W_{\mu\nu})$ satisfying the constraint equations \eqref{linear_cst} and assuming the MPD equations \eqref{MPD} are obeyed, the quantity $Q^{(1)}$ \eqref{Q1} will be conserved up to second order in the spin parameter, \textit{i.e.} $\dot Q^{(1)}=\mathcal O\qty(\mathcal S^3)$.
\end{prop}
Two independent classes of solutions to these constraint equations can be constructed:
\begin{itemize}[label=$\diamond$]
    \item For $X_\mu\neq 0$, the $[0,2]$ equation \eqref{[0,2]} requires that $X_\mu$ should be a Killing vector field,
    \begin{align}
        \nabla_{(\mu}X_{\nu)}=0.
    \end{align}
    In this case, making use of the Kostant formula
    \begin{align}
        \nabla_\alpha\nabla_\beta X_\mu=R_{\mu\nu\alpha\beta}X^\nu\, ,
    \end{align}
    which holds for any Killing vector $X_\mu$, the $[1,2]$ constraint \eqref{[1,2]} reduces to
    \begin{align}
        \nabla_\mu\qty(W_{\alpha\beta}-2\nabla_\alpha X_\beta)\hat p^\mu S^{\alpha\beta}=\mathcal O\qty(\mathcal S^3).
    \end{align}
    It is clear that this constraint as well as the $[2,4]$ constraint \eqref{[2,4]} are solved by
    \begin{align}
        W_{\alpha\beta}=\frac{1}{2}\nabla_\alpha X_\beta.
    \end{align}
    A little more work is necessary to show that the remaining constraint Eq. \eqref{[2,2]} also holds for this value of $Y_{\alpha\beta}$. At the end of the day, we have recovered the well-known conservation of the quantity
    \begin{align}
        Q^{(1)}=X_\mu p^\mu+\frac{1}{2}\nabla_\mu X_\nu S^{\mu\nu}.
    \end{align}
    The conservation can be shown to be exact and to hold at any order of the multipolar expansion \cite{1977GReGr...8..197E}.
    \item A second, independent solution may be obtained by considering $X_\mu=0$. In this case, the constraint equations \eqref{linear_cst} reduce to
        \begin{subequations}\label{linear_cst2}
        \begin{align}
            [1,2]:&\qquad\nabla_\mu Y_{\alpha\nu}s^\alpha\hat p^\mu\hat p^\nu=\mathcal O\qty(\mathcal S^3),\label{[1,2]Q}\\
            [2,2]:&\qquad Y_{\mu\nu}\mathcal L^{*\mu\nu}=\mathcal O\qty(\mathcal S^3),\\
            [2,4]:&\qquad W_{\lambda\mu}\qty(\mu D\tud{\lambda}{\nu}-\mathcal L\tud{\lambda}{\nu})\hat p^\mu\hat p^\nu=\mathcal O\qty(\mathcal S^3).
        \end{align}
        \end{subequations}
       Eq. \eqref{[1,2]Q} enforces $Y_{\mu\nu}$ to be a Killing-Yano tensor, up to second order corrections in the spin parameter:
        \begin{align}
            \nabla_{(\mu}Y_{\nu)\alpha}=\mathcal O\qty(\mathcal S^2).
        \end{align}
        We consequently recover Rüdiger's linear invariant \cite{doi:10.1098/rspa.1981.0046},
        \begin{align}\label{QY}
            Q_Y=Y_{\alpha\beta}^* S^{\alpha\beta},
        \end{align}
        which is well-known to be conserved at linear order in the spin magnitude. The conservation at second order assuming induced quadrupoles will be discussed in Section \ref{sec:linCBB}.
\end{itemize}

\subsection{Quadratic constraint}

It is useful to decompose the quadratic quantity \eqref{quadratic_ansatz} as
\begin{align}
    Q^{(2)}=Q_R+ Q^\text{quad}\label{quasi_inv}
\end{align}
where
\begin{align}
    Q_R\triangleq K_{\mu\nu}p^\mu p^\nu+L_{\mu\nu\rho}S^{\mu\nu}p^\rho,\qquad
    Q^\text{quad}\triangleq M_{\alpha\beta\gamma\delta}S^{\alpha\beta}S^{\gamma\delta}.
\end{align}
Here, the tensors satisfy the following identities
\begin{align}
    K_{\mu\nu}=K_{(\mu\nu)},\qquad
    L_{\mu\nu\rho}=L_{[\mu\nu]\rho},\qquad
    M_{\alpha\beta\gamma\delta}=M_{[\alpha\beta]\gamma\delta}=M_{\alpha\beta[\gamma\delta]}=M_{\gamma\delta\alpha\beta}.
\end{align}

The linear-in-spin quantity $Q^\text{lin}$ has been extensively studied in \cite{Compere:2021kjz}. At linear order in the spin magnitude, the variation of $Q^{(2)}$ and the variation of $Q^\text{lin}$ coincide. It was shown to be given by
\begin{equation}
\dot{Q}^{(2)}=  \mu^3  U_{\mu\nu\rho}\hat p^{\mu}\hat p^{\nu}\hat p^\rho + 2 \mu^2 {}^\star V_{\alpha\mu\nu\rho}s^\alpha \hat p^\mu \hat p^\nu \hat p^\rho + \mathcal O(\mathcal S^2). \label{vars1}
\end{equation}
The explicit expressions of the tensors $U_{\mu\nu\rho}$ and $V_{\alpha\mu\nu\rho\sigma}$ can be found in \cite{Compere:2021kjz}.
The conservation conditions at zeroth and first order $U_{(\mu\nu\rho)}=0$, ${}^\star V\tud{\alpha}{(\mu\nu\rho)}=0$ are unchanged by the presence of quadrupolar terms in the MPD equations \eqref{MPD}. In \cite{Compere:2021kjz}, we showed that, at first order in the spin parameter, the only non-trivial stationary and axisymmetric solution to these equations above a Kerr background was Rüdiger's quadratic quasi-invariant $Q_R$. This solution corresponds to
\begin{align}
    K_{\mu\nu}&=Y_{\mu\lambda}Y\tud{\lambda}{\nu},\quad
    L_{\alpha\beta\gamma}=\frac{2}{3}\nabla_{[\alpha}K_{\beta]\gamma}+\frac{4}{3}\epsilon_{\alpha\beta\gamma\delta}\nabla^\delta\mathcal Z,\label{R2}
\end{align}
where $Y_{\alpha\beta}$ is Kerr's Killing-Yano tensor and where we have defined the scalar $\mathcal Z\triangleq \frac{1}{4}Y^*_{\alpha\beta}Y^{\alpha\beta}$. We can compactly write $L_{\mu\nu\rho}S^{\mu\nu}p^\rho=\qty(\epsilon_{\mu\nu\rho\sigma}\xi_\lambda Y^{\lambda\sigma}+{}^\star Y_{\mu\nu}\xi_\rho)S^{\mu\nu}p^\rho=4 \xi^\lambda\epsilon_{\lambda \mu\sigma[\rho}Y_{\nu]}{}^\sigma S^{\mu\nu}p^\rho$.

From this point, we will assume that the zeroth and first orders in the spin magnitude are solved by Rüdiger's solution \eqref{R2}. This completely cancels the $\mathcal O(\mathcal S^0)$ and $\mathcal O(\mathcal S^1)$ terms in the constraint equations. The presence of quadrupole terms in the MPD equations \eqref{MPD} will only appear at quadratic order in $\mathcal S$. Hence, we are left with only one constraint, which is of grading $[2,3]$. The derivation of this quadratic constraint is too long to be provided in the main text and can be found instead in Appendix \ref{app:quadratic_constraint}. We have now demonstrated the following proposition:
\begin{prop}
Any tensor $N_{\alpha\beta\gamma\delta}$ possessing the same algebraic symmetries as the Riemann tensor and satisfying the constraint equation
\begin{align}
\bigg[4\nabla_\mu &N_{\alpha\nu\beta\rho}+2\kappa\nabla_{[\alpha}\mathcal M^{(1)}_{\vert\mu\vert \nu]\beta\rho} +\kappa\qty(g_{\alpha\mu}Y_{\lambda\nu}-g_{\mu\nu}Y_{\lambda\alpha})\xi_\kappa\,^*\!R\tud{\lambda\kappa}{\beta\rho}\nonumber\\
    &+\qty(2\kappa Y_{\alpha\mu}\xi_\lambda+\qty(2-\kappa)\qty(Y_{\lambda\mu}\xi_\alpha+Y_{\alpha\lambda}\xi_\mu)+{3\kappa}g_{\alpha\mu}\nabla_\lambda\mathcal Z)\,^*\!R\tud{\lambda}{\nu\beta\rho}-3\kappa g_{\mu\nu}\nabla_\lambda\mathcal Z\,^*\!R\tud{\lambda}{\alpha\beta\rho}\nonumber\\
    &+(3\kappa-2)\nabla_\mu\mathcal Z R^*_{\nu\alpha\beta\rho}\bigg]s^\alpha s^\beta \hat p^\mu \hat p^\nu \hat p^\rho
    \stackrel{!}{=}\mathcal O( \mathcal S^3),\label{developped_cst}
\end{align}
where\footnote{This notation $\mathcal M^{(1)}_{\alpha\beta\gamma\delta}$ will become clearer later on.}
\begin{align}
    \mathcal M^{(1)}_{\alpha\beta\gamma\delta}&\triangleq K_{\alpha\lambda}R\tud{\lambda}{\beta\gamma\delta}\label{M1}
\end{align}
gives rise to a quantity
\begin{align}
    \mathcal Q^{(2)}=Q_R+M_{\alpha\beta\gamma\delta}S^{\alpha\beta}S^{\gamma\delta},\qquad M_{\alpha\beta\gamma\delta}\triangleq\, ^*\!N^*_{\alpha\beta\gamma\delta}.\label{explicit_quantity}
\end{align}
which is conserved up to second order in the spin parameter for the  MPD equations with spin-induced quadrupole \eqref{MPD}, \textit{i.e.} $\dot{\mathcal Q}^{(2)}=\mathcal O\qty(\mathcal S^3)$.
\end{prop}

Our next goal with be to find a way to disentangle the $\kappa=1$ and the $\kappa\neq 1$ problems. Without loss of generality, we set
\begin{align}
N_{\alpha\beta\gamma\delta}={\triangleq N^\text{BH}_{\alpha\beta\gamma\delta}}+(\kappa-1) N^{\text{NS}}_{\alpha\beta\gamma\delta},
\end{align}
Because $\kappa$ is \textit{a priori} arbitrary, the constraint \eqref{developped_cst} turns out to be equivalent to the two independent equations
\begin{align}
    &\bigg[4\nabla_\mu N^{\text{BH}}_{\alpha\nu\beta\rho}+2\nabla_{[\alpha}\mathcal M^{(1)}_{\vert\mu\vert \nu]\beta\rho} +\qty(g_{\alpha\mu}Y_{\lambda\nu}-g_{\mu\nu}Y_{\lambda\alpha})\xi_\kappa\,^*\!R\tud{\lambda\kappa}{\beta\rho}+\big(2 Y_{\alpha\mu}\xi_\lambda+\qty(Y_{\lambda\mu}\xi_\alpha+Y_{\alpha\lambda}\xi_\mu)\nonumber\\
    &\quad+{3}g_{\alpha\mu}\nabla_\lambda\mathcal Z\big)\,^*\!R\tud{\lambda}{\nu\beta\rho}-3 g_{\mu\nu}\nabla_\lambda\mathcal Z\,^*\!R\tud{\lambda}{\alpha\beta\rho}+\nabla_\mu\mathcal Z R^*_{\nu\alpha\beta\rho}\bigg]s^\alpha s^\beta \hat p^\mu \hat p^\nu \hat p^\rho
    =\mathcal O( \mathcal S^3)\label{developped_cst_BH}
\end{align}
and
\begin{align}
    &\bigg[4\nabla_\mu N^{\text{NS}}_{\alpha\nu\beta\rho}+2\nabla_{[\alpha}\mathcal M^{(1)}_{\vert\mu\vert \nu]\beta\rho} +\qty(g_{\alpha\mu}Y_{\lambda\nu}-g_{\mu\nu}Y_{\lambda\alpha})\xi_\kappa\,^*\!R\tud{\lambda\kappa}{\beta\rho}+\big(2 Y_{\alpha\mu}\xi_\lambda-\qty(Y_{\lambda\mu}\xi_\alpha+Y_{\alpha\lambda}\xi_\mu)\nonumber\\
    &\quad+{3}g_{\alpha\mu}\nabla_\lambda\mathcal Z\big)\,^*\!R\tud{\lambda}{\nu\beta\rho} -3 g_{\mu\nu}\nabla_\lambda\mathcal Z\,^*\!R\tud{\lambda}{\alpha\beta\rho}+3\nabla_\mu\mathcal Z R^*_{\nu\alpha\beta\rho}\bigg]s^\alpha s^\beta \hat p^\mu \hat p^\nu \hat p^\rho=\mathcal O( \mathcal S^3).\label{developped_cst_NS}
\end{align}
In the continuation, we will refer to these two problems are respectively the ``black hole problem'' ($\kappa=1$) and the ``neutron star problem'' ($\kappa\neq 1$). Their resolutions are independent and will be addressed separately. Notice that the overall quasi-conserved quantity is given by
\begin{align}
Q^{(2)}=Q_\text{R}+\ Q_{\text{BH}}+(\kappa-1)\, Q_\text{NS}. \label{conserved_splitting}
\end{align}
The contributions $Q_{\text{BH}}$ and $ Q_\text{NS}$ can be directly computed from the corresponding $N_{\alpha\beta\gamma\delta}$ tensor through Eq. \eqref{explicit_quantity}.

\section{Kerr covariant formalism: generalities}\label{sec:CBB}

In this Section, we will show that the very structure of Kerr spacetime allows us to reduce the  differential constraint equations \eqref{linear_cst}-\eqref{developped_cst_BH}-\eqref{developped_cst_NS} to \textit{purely algebraic} relations. It is then possible to find a unique non-trivial solution to the black hole constraint \eqref{developped_cst_BH}, as will be demonstrated in Section \ref{sec:BH}. It also enables to provide an algebraic way for solving the  $\kappa\neq 1$ linear and quadratic problems (\textit{i.e.} Eq. \eqref{linear_cst} and \eqref{developped_cst_NS}, respectively), as will be discussed in Section \ref{sec:NS}.

\subsection{Covariant building blocks for Kerr}
In Kerr spacetime, the constraint equations \eqref{linear_cst}-\eqref{developped_cst_BH}-\eqref{developped_cst_NS} can be \textit{fully} expressed in terms of the basic tensors that live on the manifold (that is the metric $g_{\mu\nu}$, the Levi-Civita tensor $\epsilon_{\mu\nu\rho\sigma}$ and the Kronecker symbol $\delta^\mu_\nu$) and of three additional tensorial structures: the timelike Killing vector field $\xi^\mu$, the complex scalar
\begin{align}\label{defR}
    \mathcal R\triangleq r+ia\cos\theta
\end{align}
and the 2-form
\begin{align}
    N_{\alpha\beta}\triangleq -i G_{\alpha\beta\mu\nu}l^\mu n^\nu.
\end{align}
We use the convention $\varepsilon_{tr\theta\phi}=-1$ as in \cite{Compere:2021kjz}. Here, $l^\mu$ and $n^\nu$ are the two principal null directions of Kerr:
\begin{align}
    l^\mu\triangleq\frac{1}{\Delta}\mqty(r^2+a^2,&\Delta,&0,&a),\qquad n^\mu\triangleq\frac{1}{2\Sigma}\mqty(r^2+a^2,&-\Delta,&0,&a)
\end{align}
and $G\tdu{\alpha\beta}{\gamma\delta}$ is (four times) the projector
\begin{align}
    G\tdu{\alpha\beta}{\gamma\delta}\triangleq 2\delta^{[\gamma}_\alpha\delta^{\delta]}_\beta-i\epsilon\tdu{\alpha\beta}{\gamma\delta}.
\end{align}
Notice that we have the property
\begin{equation}
N_{\alpha\beta}=\frac{2}{\xi^2} ( \nabla_{[\alpha}\mathcal R\xi_{\beta]^*}+i\nabla_{[\alpha}\mathcal R\xi_{\beta]}).
\end{equation}
The Killing-Yano and Riemann tensors can be written algebraically in terms of these objects:
\begin{align}
    Y_{\alpha\beta}=-\frac{1}{2}\mathcal R N_{\alpha\beta}+c.c.,\qquad R_{\alpha\beta\gamma\delta}=M\Re\qty(\frac{3N_{\alpha\beta}N_{\gamma\delta}-G_{\alpha\beta\gamma\delta}}{\mathcal R^3}).
\end{align}
Moreover, they obey  the following closed differential relations,
\begin{align}
    i\nabla_\alpha\mathcal R= N_{\alpha\beta}\xi^\beta,\quad i\nabla_\gamma\qty(\mathcal R N_{\alpha\beta})=G_{\alpha\beta\gamma\delta}\xi^\delta,\quad i\nabla_\alpha\xi_\beta=-\frac{M}{2}\qty(\frac{N_{\alpha\beta}}{\mathcal R^2}-\frac{\bar N_{\alpha\beta}}{\bar{\mathcal{R}}^2}).\label{diff_to_alg}
\end{align}
All the derivatives appearing in the constraints can consequently be expressed in terms of purely algebraic relations between the covariant building blocks.

\subsection{Some identities}

Let us first derive some useful identities. Many of them can be found in \cite{floyd_phdthesis,Yasui:2011pr}. We have the algebraic identities
\begin{align}
    N_{\alpha\beta}N\tud{\beta}{\gamma}=-g_{\alpha\gamma},\quad N_{\alpha\beta}\bar N^{\alpha\beta}=0.
\end{align}
Notice that this first relation yields
\begin{align}
    N_{\alpha\beta}N^{\alpha\beta}=4.
\end{align}
Both $N_{\alpha\beta}$ and $G\tdu{\alpha\beta}{\gamma\delta}$ are self-dual tensors:
\begin{align}
    N^*_{\alpha\beta}=i N_{\alpha\beta},\qquad ^*G\tdu{\alpha\beta}{\gamma\delta}={G^*}\tdu{\alpha\beta}{\gamma\delta}=i G\tdu{\alpha\beta}{\gamma\delta}.
\end{align}
This leads to the relations
\begin{align}
    ^*R_{\alpha\beta\gamma\delta}=R^*_{\alpha\beta\gamma\delta}=-M\Im\qty(\frac{3N_{\alpha\beta}N_{\gamma\delta}-G_{\alpha\beta\gamma\delta}}{\mathcal R^3}),\qquad \bar N^*_{\alpha\beta}=-i\bar N_{\alpha\beta}.
\end{align}
Given the identities just derived, the only non-trivial contraction of the 2-form that can be written is
\begin{align}
    h_{\mu\nu}\triangleq N\tdu{\mu}{\alpha}\bar N_{\nu\alpha}.
\end{align}
It is a real, symmetric and traceless tensor:
\begin{align}
    h_{\mu\nu}=h_{(\mu\nu)}=\bar h_{\mu\nu},\qquad h^\mu_\mu=0.
\end{align}
Using the previous identities, one shows that
\begin{align}
    \mathcal Z=-\frac{1}{2}\Im\qty(\mathcal R^2).
\end{align}
This yields
\begin{align}
    \nabla_\alpha\mathcal Z=-\Re\qty(\mathcal R \xi^\lambda N_{\lambda\alpha}).
\end{align}
The Killing tensor can be written as
\begin{align}
    K_{\mu\nu}=-\frac{1}{2}\qty(\Re\qty(\mathcal R^2)g_{\mu\nu}+\abs{\mathcal R^2}h_{\mu\nu}).
\end{align}
Its trace is simply
\begin{align}
    K&=-2\Re\qty(\mathcal R^2).
\end{align}
Other useful identities include
\begin{align}
   &N_{\lambda\kappa}G\tud{\lambda\kappa}{\beta\rho}=4N_{\beta\rho},\qquad N_{\lambda\kappa}\bar G\tud{\lambda\kappa}{\beta\rho}=\bar N_{\lambda\kappa}G\tud{\lambda\kappa}{\beta\rho}=0,\nonumber\\
    &\bar N_{\lambda\kappa}\bar G\tud{\lambda\kappa}{\beta\rho}=4\bar N_{\beta\rho},\qquad h_{\alpha\lambda}N\tud{\lambda}{\beta}=\bar N_{\alpha\beta}.\label{lastId}
\end{align}

\subsection{Basis of contractions}

Our goal is now to rewrite Eqs. \eqref{linear_cst}-\eqref{developped_cst_BH}-\eqref{developped_cst_NS} as \textit{scalar} (that is, fully-contracted) equations involving only contractions between the Kerr covariant building blocks and the dynamical variables $s^\alpha$ and $\hat p^\alpha$. We define
\begin{align}
    \mathcal S^2&\triangleq s_\alpha s^\alpha,\qquad
    \mathcal P^2\triangleq -\hat p_\alpha \hat p^\alpha,\qquad
    \mathcal A\triangleq s_\alpha\hat p^\alpha.
\end{align}
We will naturally set $\mathcal P^2=1$ at the end of the computation, but we find useful to keep this quantity explicit in the intermediate algebra. Since $\mathcal A$ is unphysical, it will disappear from any physical expression but it will appear in intermediate computations.

We further define the following quantities at least linear in either $\hat p^\mu$ or $s^\mu$,
\begin{subequations}\label{ABCD}
\begin{align}
    A&\triangleq N_{\lambda\mu}\xi^\lambda\hat p^\mu,\qquad
    B\triangleq N_{\alpha\mu}s^\alpha \hat p^\mu,\qquad
    C\triangleq N_{\lambda\alpha}\xi^\lambda s^\alpha,\qquad
    D\triangleq h_{\lambda\alpha}\xi^\lambda s^\alpha,\qquad
    E\triangleq -\xi_\alpha \hat p^\alpha,\\
    E_s&\triangleq-\xi_\alpha s^\alpha,\qquad
    F\triangleq h_{\lambda\mu}\xi^\lambda \hat p^\mu,\qquad
    G\triangleq h_{\alpha\mu}s^\alpha\hat p^\mu,\qquad
    H\triangleq h_{\mu\nu}\hat p^\mu\hat p^\nu,\qquad
    I\triangleq h_{\alpha\beta} s^\alpha s^\beta .
\end{align}
\end{subequations}
Because of the algebraic identities derived above, these scalars form a spanning set of scalars built from contractions among the Kerr covariant building blocks. Any higher order contraction between building blocks will reduce to a product of the ones provided in the above list with coefficients that may depend upon $M$, $a$ and $\mathcal R$. The quantities $A,B,C$ are complex while the others are real. Notice that we don't have to include $\xi^2$ in our basis of building blocks, since it is a function of $\mathcal R$ and $M$,
\begin{align}
    \xi^2=-1+2M\Re\qty(\mathcal R^{-1}).
\end{align}
It can consequently be written in terms of the other quantities.

We further define the following quantity independent from $M$:
\begin{equation}
J \triangleq \xi^\alpha (h_{\alpha\beta}+g_{\alpha\beta})\xi^\beta =  \frac{2a^2\sin^2\theta}{r^2+a^2 \cos^2\theta}
\end{equation}
where the last expression is evaluated in Boyer-Lindquist coordinates. We can now use $J$ as a Kerr covariant substitute for $a$: we will consider in what follows quantities built from \eqref{ABCD}, the complex scalar $\mathcal R$, the mass $M$ and $J$.

\subsection{A $\mathbb Z_2$ grading}

We now define a $\mathbb Z_2$ grading $\{,\}$ as follows. We note that the determining equations for the covariant building blocks for Kerr from Eq. \eqref{defR} to Eq. \eqref{lastId} are invariant under the following $\mathbb Z_2$ grading: $\{g_{\alpha\beta}\}=\{M\}=\{x^\mu\}=\{\nabla_\mu\}=\{\mathcal R\}=\{\mathcal Z\}=\{ G_{\alpha\beta\gamma\delta}\}=\{h_{\mu\nu}\}=\{K_{\mu\nu}\}=+1$ and $\{N_{\alpha\beta}\}=\{\xi^\alpha\}=\{Y_{\alpha\beta}\}=-1$.

Further assigning $\{s^\mu\}=\{p^\mu\}=+1$, we deduce that
\begin{equation}
\{A\}=\{C\}=\{G\}=\{H\}=\{I\}=+1,\qquad \{B\}=\{D\}=\{E\}=\{E_s\}=\{F\}=-1.
\end{equation}
Since the constraints \eqref{developped_cst_BH}, \eqref{developped_cst_NS} have grading $+1$, the odd quantities will have to be combined in pairs in order to build a solution to the constraint.

We define the $(s,p)^\pm$ grading of an expression as the $s$ numbers of $s^\alpha$ and $p$ number of $\hat p^\alpha$ factors in the expression with the sign $\pm$ indicating the $\mathbb Z_2$ grading. The complete list of the lowest $s+p=1$ and $s+p=2$ grading spanning elements is given in Table \ref{tab:spbasis}. The list of spanning elements of grading $(s,p)^\pm$ for $s+p \geq 3$ is obtained iteratively by direct product of the lower order basis elements. For example, the independent real terms of grading $(2,1)^+$ are obtained from $(2,0)^+ \times (0,1)^+$, $(2,0)^- \times (0,1)^-$, $(1,1)^+ \times (1,0)^+$ and $(1,1)^- \times (1,0)^-$ with duplicated elements suppressed.

\begin{table}[!htb]
    \centering
    \begin{tabular}{|c|l|c|}\hline
     $(s,p)^\pm$ & Spanning set & Real dimension of the spanning set \\ \hline
        $(1,0)^+$ &  $C$ & 2 \\
       $(1,0)^-$ & $D$,\; $E_s$ &2 \\
$(0,1)^+$& $A$ & 2\\
$(0,1)^-$& $E$,\; $F$ &2 \\
$(2,0)^+$ &$I$,\; $\mathcal S$,\; $(1,0)^+ \times (1,0)^+$,\; $(1,0)^- \times (1,0)^-$ & 8\\
$(2,0)^-$& $(1,0)^+ \times (1,0)^-$& 4\\
$(1,1)^+$& $G$, \;$\mathcal A$, \; $(1,0)^+ \times (0,1)^+$ , \; $(1,0)^- \times (0,1)^-$  & 10\\
$(1,1)^-$& $B$, \; $(1,0)^+ \times (0,1)^-$,\; $(1,0)^- \times (0,1)^+$   & 10\\
$(0,2)^+$& $H$,\; $\mathcal P$,\;  $(0,1)^+ \times (0,1)^+$,\; $(0,1)^- \times (0,1)^-$  &8\\
$(0,2)^-$&  $(0,1)^+ \times (0,1)^-$ & 4\\
\hline    \end{tabular}
    \caption{Spanning set of elements with $(s,p)^\pm$ grading with $s+p=1$ and $s+p=2$.}
    \label{tab:spbasis}
\end{table}

Of prime importance for solving the linear and quadratic constraint equations will be the elements of gradings $(2,2)^+$ and $(2,3)^+$. Their respective spanning sets contain $118$ and $284$ elements, which are explicitly listed in the appended \textit{Mathematica} notebooks.

\subsection{The $\alpha$-$\omega$ basis}

We note that the covariant building blocks all depend on $\mathcal R$ through real and imaginary parts of expressions containing fractions of $\mathcal R$ and $\bar{\mathcal R}$. We find therefore natural to define the objects ($n,p \in\mathbb Z$ and $K=1,A,B,C,\ldots,J$ or any combination of these objects):
\begin{align}
    \alpha_K^{(n,p)}\triangleq\Re\qty(\frac{K\bar{\mathcal R}^n}{\mathcal R^p}),\qquad\omega_K^{(n,p)}\triangleq\Im\qty(\frac{K\bar{\mathcal R}^n}{\mathcal R^p}).
\end{align}
They satisfy the following properties:
\begin{align}
    \alpha_{iK}^{(n,p)}&=-\omega_K^{(n,p)},\qquad
    \omega_{iK}^{(n,p)}=\alpha_K^{(n,p)},\qquad
    \alpha_{\bar K}^{(n,p)}=\alpha_K^{(-p,-n)},\qquad
    \omega_{\bar K}^{(n,p)}=-\omega_K^{(-p,-n)}.\label{ao_prop}
\end{align}
Moreover, one has
\begin{align}
    &\abs{\mathcal R}^2\alpha_K^{(n,p)}=\alpha_K^{(n+1,p-1)},\qquad
    \Re\qty(\mathcal R^2)\alpha_K^{(n,p)}=\frac{1}{2}\qty[\alpha_K^{(n,p-2)}+\alpha_K^{(n+2,p)}],\nonumber\\
    &
    \alpha_{\mathcal R^k K}^{(n,p)}=\alpha_K^{(n,p-k)},\qquad
    \alpha_{\bar{\mathcal R}^kK}^{(n,p)}=\alpha_K^{(n+k,p)}.
\end{align}
The same properties hold with $\omega$ instead of $\alpha$. Finally, $\alpha$ and $\omega$ are linear in their subscript argument with respect to real-valued functions.
Let us denote $\ell^\mu=\hat p^\mu$ or $s^\mu$. Then, for any $T\triangleq T_{\mu_1\ldots\mu_k}\ell^{\mu_1}\ldots\ell^{\mu_k}$, we define the operator $\hat\nabla$ as
\begin{align}
    \hat \nabla T\triangleq \hat p^\lambda\nabla_\lambda\qty(T_{\mu_1\ldots\mu_k})\ell^{\mu_1}\ldots\ell^{\mu_k}.
\end{align}
Making use of the identities
\begin{align}
    \hat\nabla \mathcal R^n=in\mathcal R^{n-1}A,\qquad\hat\nabla\bar{\mathcal R}=-in\bar{\mathcal R}^{n-1}\bar A,
\end{align}
we get the following relations:
\begin{align}
    \hat\nabla\alpha_K^{(n,p)}&=\alpha_{\hat\nabla K}^{(n,p)}+n\omega_{K\bar A}^{(n-1,p)}+p\omega_{KA}^{(n,p+1)},\qquad
    \hat\nabla\omega_K^{(n,p)}=\omega_{\hat\nabla K}^{(n,p)}-n\alpha_{K\bar A}^{(n-1,p)}-p\alpha_{KA}^{(n,p+1)}.\label{deriv_ao}
\end{align}

We use dimensions such that $G=c=1$. Given the large amount of definitions, we find useful to summarize the mass dimensions $[,]$ of all quantities in order to keep track of the powers of the mass $M$ that can arise. We have the following mass dimensions  $[\nabla_\mu]=-1$, $[g_{\mu\nu}]=[\xi^\alpha]=[N_{\alpha\beta}]=[G_{\alpha\beta\mu\nu}]=0$, $[x^\mu]=[M]=[\mathcal R]=[Y_{\alpha\beta}]=1$ and $[K_{\alpha\beta}]=[K]=2$. We deduce
\begin{eqnarray}
[X]=0,\qquad [\alpha_X^{(n,p)}]=[\omega_X^{(n,p)}]=n-p,
\end{eqnarray}
where $X$ is any function of the set $A,B,C,D,E,E_s,F,G,H,I$ defined in Eqs. \eqref{ABCD}.

\section{Kerr covariant formalism: reduction of the constraints}\label{sec:Reduction}

\subsection{Linear constraint}
\label{sec:linCBB}
In this section, we will reduce the linear constraint equations in the case where $Y_{\mu\nu}$ is Kerr Killing-Yano tensor. Using the explicit form of $\mathcal L_{\mu\nu}$ as defined in \eqref{L} and expressing the quantity $   \mathcal L^*_{\mu\nu} Y^{\mu\nu}$ in terms of covariant building blocks we find after evaluation
\begin{align}
    \mathcal L^*_{\mu\nu} Y^{\mu\nu}=0,
\end{align}
and therefore the $[2,2]$ constraint is automatically fulfilled. The $[2,4]$ constraint can be rewritten
\begin{align}
    -\mu Y_{\alpha\beta}p^{[\alpha}D^{\beta]*\lambda}\hat p_\lambda+Y_{\alpha\beta}p^{[\alpha}\mathcal L^{\beta]*\lambda}\hat p_\lambda=\mathcal O\qty(\mathcal S^3).
\end{align}
A direct computation shows that
\begin{align}
    \mu Y_{\alpha\beta}p^{[\alpha}D^{\beta]*\lambda}\hat p_\lambda&=-3M\qty(\mathcal A H+\mathcal P^2 G)\omega_B^{(1,3)},\\
    Y_{\alpha\beta}p^{[\alpha}\mathcal L^{\beta]*\lambda}\hat p_\lambda&=3\kappa M\qty(\mathcal A H+\mathcal P^2 G)\omega_B^{(1,3)}.
\end{align}
Using these identities and defining $\kappa\triangleq 1+\delta\kappa$, the $[2,4]$ constraint takes the very simple form
\begin{align}
    -3M\delta\kappa\qty(\mathcal AH+\mathcal P^2G)\omega_B^{(1,3)}=\mathcal O\qty(\mathcal S^3).
\end{align}
It is automatically fulfilled for the test body being a black hole, because $\delta\kappa=0$ in this case. However, if the test body is a neutron star, $\delta\kappa\neq 0$ and the $[2,4]$ constraint is not obeyed anymore. Therefore, $Q_Y$ is not anymore a constant of the motion at second order in the spin in the NS case.

A way to enable the $[2,4]$ constraint to be solvable in the neutron star case is to supplement the Ansatz for the conserved quantity with a term
\begin{align}
    Q^{(1)}_\text{NS}=\delta\kappa M_{\alpha\beta\mu\gamma\delta}S^{\alpha\beta}S^{\gamma\delta}p^\mu.
\end{align}
The conservation equation will then acquire a correction given by
\begin{align}
    \dot Q^{(1)}_\text{NS}&=\delta\kappa v^\lambda\nabla_\lambda\qty(M_{\alpha\beta\mu\gamma\delta}S^{\alpha\beta}S^{\gamma\delta}p^\mu)\\
    &=\delta\kappa\hat p^\lambda\nabla_\lambda M_{\alpha\beta\mu\gamma\delta}S^{\alpha\beta}S^{\gamma\delta}p^\mu+\mathcal O\qty(\mathcal S^3)\\
    &=4\delta\kappa\hat p^\lambda\nabla_\lambda N_{\alpha\beta\mu\gamma\delta}s^\alpha\hat p^\beta s^\gamma\hat p^\delta p^\mu+\mathcal O\qty(\mathcal S^3)
\end{align}
where $N_{\alpha\beta\mu\gamma\delta}={}^\star M_{\alpha\beta\mu\gamma\delta}^\star$. In our scalar notation, it corresponds to supplement the $[2,4]$ constraint with a term $4\delta\kappa\hat \nabla N$, with $N$ being of grading $[2,3]$. The constraint to be solved then takes the simple form
\begin{align}
    \boxed{
    \hat\nabla N=\frac{3M}{4}\qty(\mathcal A H+\mathcal P^2 G)\omega_B^{(1,3)}.
    }
\end{align}

It is useful to summarize the discussion by the two following statements:
\begin{main_result}
Rüdiger's linear invariant $Q_Y=Y_{\alpha\beta}^* S^{\alpha\beta}$ is still conserved for the MPD equations at second order in the spin magnitude for spin-induced quadrupoles, \textit{i.e.} $\dot Q_Y=\mathcal O\qty(\mathcal S^3)$ provided that $\delta\kappa=0$, \textit{i.e.} if the test body possesses the multipole structure of a black hole.
\end{main_result}

\begin{preliminary_result}
Any tensor $N_{\alpha\beta\mu\gamma\delta}$ possessing the algebraic symmetries $N_{\alpha\beta\mu\gamma\delta}=N_{\gamma\delta\mu\alpha\beta}=N_{[\alpha\beta]\mu\gamma\delta}=N_{\alpha\beta\mu[\gamma\delta]}$ and satisfying the constraint equation
\begin{align}
\hat\nabla N=\frac{3M}{4}\qty(\mathcal A H+\mathcal P^2 G)\omega_B^{(1,3)}\label{lin_cst}
\end{align}
will give rise to a quantity
\begin{align}
    \mathcal Q^{(1)}=Q_Y+\delta\kappa M_{\alpha\beta\mu\gamma\delta}S^{\alpha\beta}S^{\gamma\delta}p^\mu,\qquad M_{\alpha\beta\mu\gamma\delta}=\,^*\!N^*_{\alpha\beta\mu\gamma\delta}
\end{align}
which is conserved up to second order in the spin parameter for the spin-induced quadrupole MPD equations \eqref{MPD}, \textit{i.e.} $\dot Q^{(1)}=\mathcal O\qty(\mathcal S^3)$, regardless to the value taken by $\delta\kappa$.
\end{preliminary_result}

\subsection{Quadratic constraint}

\subsubsection{Some identities}
Before going further on, it is useful to notice that all the covariant building blocks combinations that will appear in our equations will not be linearly independent. Actually, a direct computation shows that
\begin{subequations}\label{lin_dep_terms}
\begin{align}
    2\abs{B}^2+2\mathcal AG+\mathcal P^2I-\mathcal S^2H&=0,\\
    \qty(\mathcal A H+\mathcal P^2G)\omega^{(1,3)}_C&=\qty(\mathcal AG+\mathcal P^2I)\omega^{(1,3)}_A-\qty(\mathcal A F+\mathcal P^2 D)\omega^{(1,3)}_B\nonumber\\
    &\quad+\qty(\mathcal A^2+\mathcal P^2\mathcal S^2)\omega^{(1,3)}_{\bar A}+\qty(\mathcal AE+\mathcal P^2E_s)\omega^{(1,3)}_{\bar B},\\
    \omega^{(1,3)}_{\bar AB^2}&=-\abs{B}^2\omega^{(1,3)}_A-\qty(\mathcal AF-EG+E_sH+\mathcal P^2D)\omega_B^{(1,3)}.
\end{align}
\end{subequations}
Moreover, let us mention that the identities
\begin{align}
    &\omega^{(0,k)}_K\alpha_L^{(n,p)}=\frac{1}{2}\qty[\omega_{KL}^{(n,p+k)}-\omega_{\bar KL}^{(n-k,p)}],\\
    &\Re\qty(\mathcal R^2)\Im\qty(\frac{K}{\mathcal R^4})=\frac{1}{2}\qty(\omega^{(0,2)}_K+\omega^{(2,4)}_K),\nonumber\\
    &\abs{\mathcal R}^2\Im\qty(\frac{K}{\mathcal R^4})=\omega^{(1,3)}_K
    \end{align}
will be useful in the following computations.

\subsubsection{Reducing the $\mathcal M^{(1)}$ contribution}
Our goal is here to compute the contribution
\begin{align}
    \text{DM}\triangleq 2\nabla_{[\alpha|}\mathcal M^{(1)}_{\mu|\nu]\beta\rho}s^\alpha s^\beta\hat p^\mu\hat p^\nu\hat p^\rho
\end{align}
in some details, as a proof of principle of the computations to follow, which will not be developed in full details. Noticing the identity
\begin{align}
    \nabla_\mu K_{\alpha\beta}=2\epsilon_{\lambda\rho\mu(\alpha}Y\tud{\lambda}{\beta)}\xi^\rho,\label{DK}
\end{align}
we get
\begin{align}
    \nabla_\mu \mathcal M^{(1)}_{\alpha\nu\beta\rho} &=K_{\alpha\lambda}\nabla_\mu R\tud{\lambda}{\nu\beta\rho}+\nabla_\nu\mathcal Z R^*_{\mu\alpha\beta\rho}-\qty(Y_{\mu\nu}\xi_\lambda+g_{\mu\nu}\nabla_\lambda\mathcal{Z}+Y_{\lambda\mu}\xi_\nu)R\tud{*\lambda}{\alpha\beta\rho}\nonumber\\
    &\quad+\qty(2\xi_\nu Y_{\lambda\alpha}-2\xi_\lambda Y_{\nu\alpha}+g_{\alpha\nu}\nabla_\lambda\mathcal Z)R\tud{*\lambda}{\mu\beta\rho}+\qty(2g_{\mu\nu}Y_{\kappa\alpha}-g_{\alpha\nu} Y_{\kappa\mu})\xi_\lambda R\tud{*\lambda\kappa}{\beta\rho}.\label{DM}
\end{align}
Making use of Eq. \eqref{DM}, one can show that
\begin{align}
    \text{DM}&=\bigg[2K_{\mu\lambda}\nabla_{[\alpha}R\tud{\lambda}{\nu]\beta\rho}+\nabla_\nu\mathcal Z R^*_{\alpha\mu\beta\rho}+\qty(2\xi_\nu Y_{\lambda\mu}+g_{\mu\nu}\nabla_\lambda\mathcal Z)R\tud{*\lambda}{\alpha\beta\rho}\nonumber\\
    &\quad-\qty(Y_{\lambda\alpha}\xi_\mu+\xi_\alpha Y_{\lambda\mu}+g_{\mu\alpha}\nabla_\lambda\mathcal Z)R\tud{*\lambda}{\nu\beta\rho}+\qty(g_{\mu\alpha}Y_{\kappa\nu}-g_{\mu\nu}Y_{\kappa\alpha})\xi_\lambda R\tud{*\lambda\kappa}{\beta\rho}\bigg]s^\alpha s^\beta\hat p^\mu\hat p^\nu\hat p^\rho.
\end{align}
Using the various identities derived above, the relations of Appendix \ref{app:CBB_identities} and performing some simple algebra, one can express this contribution in terms of linearly independent quantities as
\begin{align}
    \text{DM}&=-\frac{M}{4}\qty(\mathcal A^2+\mathcal P^2\mathcal S^2)\qty(5\omega^{(0,2)}_A+4\omega^{(1,3)}_{\bar A}+3 \omega^{(2,4)}_A)-\frac{M}{2}\qty(\mathcal AE+\mathcal P^2E_s)\qty(\omega^{(0,2)}_B-\omega^{(1,3)}_{\bar B}-3\omega^{(2,4)}_B)\nonumber\\
    &\quad+\frac{3M}{2}\qty(2\mathcal A F-EG+E_sH+2\mathcal P^2D)\omega^{(1,3)}_B-\frac{9M}{4}\omega^{(0,2)}_{AB^2}-\frac{15M}{4}\omega^{(2,4)}_{AB^2}.
\end{align}

\subsubsection{The black hole constraint equation}
Making use of the notations introduced above, the constraint equation \eqref{developped_cst_BH} can be written as
\begin{align}
    4\hat\nabla N^\text{BH}+\text{DM}+\Upsilon=\mathcal O\qty(\mathcal S^3),\label{scalar_like_constraint}
\end{align}
where
\begin{align}
    \Upsilon&\triangleq\bigg[\qty(g_{\alpha\mu}Y_{\lambda\nu}-g_{\mu\nu}Y_{\lambda\alpha})\xi_\kappa\,^*\!R\tud{\lambda\kappa}{\beta\rho}+\qty(2 Y_{\alpha\mu}\xi_\lambda+\qty(Y_{\lambda\mu}\xi_\alpha+Y_{\alpha\lambda}\xi_\mu)+3g_{\alpha\mu}\nabla_\lambda\mathcal Z)\,^*\!R\tud{\lambda}{\nu\beta\rho}\nonumber\\
    &\quad-3 g_{\mu\nu}\nabla_\lambda\mathcal Z\,^*\!R\tud{\lambda}{\alpha\beta\rho}+\nabla_\mu\mathcal Z R^*_{\nu\alpha\beta\rho}\bigg]s^\alpha s^\beta \hat p^\mu \hat p^\nu \hat p^\rho.
\end{align}
Using the scalar basis introduced above and the identities \eqref{lin_dep_terms} yields
\begin{align}
    \Upsilon&=\frac{M}{2}\qty(\mathcal A^2+\mathcal P^2\mathcal S^2)\qty[3\omega^{(0,2)}_A{+}\omega^{(1,3)}_{\bar A}]+\frac{M}{2}\qty(\mathcal AE+\mathcal P^2E_s)\qty[8\omega^{(0,2)}_B-3\omega^{(1,3)}_{\bar B}]\nonumber\\
    &\quad+\frac{3M}{2}\omega^{(0,2)}_{AB^2}+\frac{9M}{2}\abs{B}^2\omega^{(1,3)}_A-\frac{3M}{2}\qty(\mathcal AF+\mathcal P^2D)\omega^{(1,3)}_B.
\end{align}
In summary, Eq. \eqref{scalar_like_constraint} can be written as
\begin{align}
    \hat \nabla N^\text{BH}=\Upsilon_\text{BH},\label{BH_scalar_constraint}
\end{align}
with the source term
\begin{align}
\Upsilon_\text{BH}=-\frac{\text{DM}+\Upsilon}{4}.
\end{align}

\subsection{The neutron star constraint equation}
Repeating the very same procedure, the neutron star constraint \eqref{developped_cst_NS} reduces to the scalar-like equation
\begin{align}
  \hat\nabla N^\text{NS}=\Upsilon_\text{NS}
    \label{source_BH}
\end{align}
with source
\begin{align}
\Upsilon_\text{NS}&=\frac{3M}{16}\qty(\mathcal A^2+\mathcal P^2\mathcal S^2)\qty(\omega_A^{(0,2)}+\omega_A^{(2,4)}+2\omega_{\bar A}^{(1,3)})-\frac{3M}{8}\qty(\mathcal AE+\mathcal P^2E_s)\qty(\omega_B^{(0,2)}+\omega_B^{(2,4)})\nonumber\\
&\quad+\frac{15M}{16}\qty(\omega_{AB^2}^{(0,2)}+\omega_{AB^2}^{(2,4)})+\frac{15M}{8}\omega_{\bar AB^2}^{(1,3)}+\frac{3}{4}M\qty(\mathcal AF+\mathcal P^2D)\omega_B^{(1,3)}.\label{source_NS}
\end{align}

\section{Solution for the quadratic invariant in the black hole case}
\label{sec:BH}

We will now try to find a quadratic conserved quantity for the $\delta\kappa=0$ case. This corresponds to find a solution to the black hole constraint equation \eqref{BH_scalar_constraint}. In order to reach this goal, we will postulate an Ansatz for the fully-contracted quantity $N$ appearing in the left-hand side of Eq. \eqref{BH_scalar_constraint} and then use the covariant building blocks formulation to constrain the Ansatz coefficients.

\subsection{The Ansatz}
Let us consider the following Ansatz
\begin{align}
    N^\text{BH}_{\alpha\beta\gamma\delta}\triangleq\sum_{A=1}^4 \Lambda_A N_{\alpha\beta\gamma\delta}^{(A)}\label{ansatz}
\end{align}
where $\Lambda_A$ are arbitrary coefficients and where
\begin{align}
    N_{\alpha\beta\gamma\delta}^{(A)}\triangleq \,^*\!\mathcal M^{*(A)}_{\alpha\beta\gamma\delta}.
\end{align}
The quantity $\mathcal M^{(1)}_{\alpha\beta\gamma\delta}$ has been defined in Eq. \eqref{M1}, and we introduce
\begin{align}
    \mathcal M^{(2)}_{\alpha\beta\gamma\delta}&\triangleq Y\tud{\lambda}{\alpha}Y\tud{\sigma}{\gamma}R_{\lambda\beta\sigma\delta},\qquad
        \mathcal M^{(3)}_{\alpha\beta\gamma\delta}\triangleq g_{\alpha\gamma}\xi_\beta\xi_\delta,\qquad
    \mathcal M^{(4)}_{\alpha\beta\gamma\delta}\triangleq g_{\alpha\gamma}g_{\beta\delta}\xi^2.
\end{align}

Using the identities derived in Appendix \ref{app:CBB_identities}, one can show that the directional derivatives of the $N^{(A)}$ are given by
\begin{subequations}
\begin{align}
    \hat\nabla N^{(1)}&=\frac{M}{4}\qty(\mathcal A^2+\mathcal P^2\mathcal S^2)\qty(\omega^{(0,2)}_A+2\omega^{(1,3)}_{\bar A}+3\omega^{(2,4)}_A)\nonumber\\
    &\quad-\frac{M}{2}\qty(\mathcal AE+\mathcal P^2E_s)\qty(5\omega^{(0,2)}_B-2\omega^{(1,3)}_{\bar B}+3\omega^{(2,4)}_B)\quad-\frac{9M}{2}\abs{B}^2\omega^{(1,3)}_A\nonumber\\
    & -\frac{3M}{2}\qty(\mathcal AF-EG+E_sH+\mathcal P^2D)\omega^{(1,3)}_B+\frac{9M}{4}\omega^{(0,2)}_{AB^2}+\frac{15M}{4}\omega^{(2,4)}_{AB^2},\label{DN1_Kerr}\\
    \hat\nabla N^{(2)}&=-\frac{M}{4}\qty(\mathcal A^2+\mathcal P^2\mathcal S^2)\omega^{(0,2)}_A+\frac{M}{2}\qty(\mathcal A E+\mathcal P^2 E_s)\omega^{(0,2)}_B-\frac{3M}{4}\omega^{(0,2)}_{AB^2},\\
    \hat\nabla N^{(3)}&=\frac{M}{2}\qty[\qty(\mathcal S^2\mathcal P^2+\mathcal A^2)\omega^{(0,2)}_A+\qty(\mathcal A E+\mathcal P^2 E_s)\omega^{(0,2)}_B],\label{DN3_Kerr}\\
    \hat\nabla N^{(4)}&=M\qty(\mathcal S^2\mathcal P^2+\mathcal A^2)\omega^{(0,2)}_A.\label{DN4_Kerr}
\end{align}
\end{subequations}

\subsection{Solution to the constraint}

We will now look for a solution to the black hole constraint equation \eqref{BH_scalar_constraint} using the Ansatz \eqref{ansatz}, \textit{i.e.} we are seeking for specific values of the parameters $\Lambda_A$ such that Eq. \eqref{scalar_like_constraint} is fulfilled. More explicitly, one therefore requires
\begin{align}
    \sum_{A=1}^4\Lambda_A \text{DN}^{(A)}-\Upsilon_\text{BH}\stackrel{!}{=}0.
\end{align}
The left-hand side of this equation takes the form of a first order polynomial, homogeneous in the ten \textit{linearly independent} elements (as it can be shown through a direct computation)
\begin{align}
    \omega^{(0,2)}_A,\quad\omega^{(0,2)}_B,\quad\omega^{(0,2)}_{AB^2},\quad\omega^{(1,3)}_A,\quad\omega^{(1,3)}_B,\quad\omega^{(1,3)}_{\bar A},\quad\omega^{(1,3)}_{\bar B},\quad\omega^{(2,4)}_A,\quad\omega^{(2,4)}_B,\quad\omega^{(2,4)}_{AB^2}.
\end{align}
Because all the combinations of these elements implied in the constraint equation are linearly independent, all the coefficients appearing in front of these expressions   should vanish independently.

\begin{table}[b]
    \centering
    \begin{tabular}{c|c|c|c|c|c|c|c|c|c|c}
 & $\omega^{(0,2)}_A$ & $\omega^{(0,2)}_B$& $\omega^{(0,2)}_{AB^2}$ & $\omega^{(1,3)}_A$  & $\omega^{(1,3)}_B$ & $\omega^{(1,3)}_{\bar A}$ & $\omega^{(1,3)}_{\bar B}$ & $\omega^{(2,4)}_A$ & $\omega^{(2,4)}_B$ & $\omega^{(2,4)}_{AB^2}$\\
 \hline
 $\text{DN}^{(1)}$& \checkmark& \checkmark&\checkmark &\checkmark &\checkmark &\checkmark &\checkmark &\checkmark &\checkmark &\checkmark \\\hline
 $\text{DN}^{(2)}$&\checkmark &\checkmark & \checkmark & & & & & & & \\\hline
 $\text{DN}^{(3)}$& \checkmark & \checkmark & & & & & & & & \\\hline
 $\text{DN}^{(4)}$& \checkmark & & & & & & & & & \\\hline
 $\Upsilon_\text{BH}$& \checkmark &\checkmark  &\checkmark  &\checkmark  &\checkmark  &\checkmark  &\checkmark  & \checkmark & \checkmark & \checkmark
\end{tabular}
    \caption{Structure of the distribution of the different types of contractions in the various contribution to the black hole constraint equation.}
    \label{tab:terms}
\end{table}

All the terms do not appear in all the contributions, as depicted in Table \ref{tab:terms}. In order to fix the values of the Ansatz coefficients, let us proceed along the following sequence:
\begin{itemize}[label=$\diamond$]
    \item \underline{\textbf{$\omega^{(1,3)}_A$ term}:} this contribution reads
    \begin{align}
        3M\qty(-6\Lambda_1+\frac{3}{2})\abs{B}^2\omega^{(1,3)}_A.
    \end{align}
    One therefore requires
    \begin{align}
      \Lambda_1=\frac{1}{4}\label{L1_kappa}\, .
    \end{align}
    \item \underline{\textbf{$\omega^{(1,3)}_{\bar A}$, $\omega^{(1,3)}_{\bar B}$, $\omega^{(1,3)}_B$, $\omega^{(2,4)}_A$, $\omega^{(2,4)}_B$ and $\omega^{(2,4)}_{AB^2}$ terms}:} their coefficients consistently vanish when \eqref{L1_kappa} is fulfilled.
    \item \underline{\textbf{$\omega^{(0,2)}_{AB^2}$ term}:} this contribution reads
    \begin{align}
        3M\qty(3\Lambda_1-\Lambda_2\frac{1}{4})\omega^{(0,2)}_{AB^2}.
    \end{align}
    Using \eqref{L1_kappa}, this yields
    \begin{align}
       \Lambda_2=\frac{1}{2}.\label{L2}
    \end{align}
    \item \underline{\textbf{$\omega^{(1,3)}_{\bar A}$, $\omega^{(1,3)}_{\bar B}$ and $\omega^{(1,3)}_B$ terms:}} their coefficients consistently vanish when Eqs. \eqref{L1_kappa} and \eqref{L2} are fulfilled.
    \item \underline{\textbf{$\omega^{(0,2)}_{B}$ term}:} this contribution reads
    \begin{align}
        M\qty(\mathcal AE+\mathcal P^2 E_s)\qty(-10\Lambda_1+2\Lambda_2+2\Lambda_3+\frac{7}{2})\omega^{(0,2)}_B.
    \end{align}
    Using \eqref{L1_kappa} and \eqref{L2}, this yields
    \begin{align}
       \Lambda_3=-1.\label{L3}
    \end{align}
    \item \underline{\textbf{$\omega^{(0,2)}_{A}$ term}:} this contribution reads
    \begin{align}
        M\qty(\mathcal A^2+\mathcal P^2 \mathcal S^2)\qty(\Lambda_1-\Lambda_2+2\Lambda_3+4\Lambda_4+\frac{1}{4})\omega^{(0,2)}_A.
    \end{align}
    Using Eqs. \eqref{L1_kappa}, \eqref{L2} and \eqref{L3}, this finally yields
    \begin{align}
      \Lambda_4=\frac{1}{2}.\label{L4}
    \end{align}
\end{itemize}
In conclusion, the Ansatz \eqref{ansatz} gives a coherent solution to the constraint equation \eqref{BH_scalar_constraint} only if
\begin{align}
    \Lambda_1=\frac{1}{4},\quad\Lambda_2=\frac{1}{2},\quad\Lambda_3=-1,\quad\Lambda_4=\frac{1}{2}.\label{coeff_values}
\end{align}
More explicitly, it corresponds to set
\begin{align}
    \boxed{
    M_{\alpha\beta\gamma\delta}=-g_{\alpha\gamma}\qty(\xi_\beta \xi_\delta-\frac{1}{2}g_{\beta\delta}\xi^2)+\frac{1}{2}Y\tdu{\alpha}{\lambda}\qty(Y\tdu{\gamma}{\kappa}R_{\lambda\beta\kappa\delta}+\frac{1}{2}Y\tdu{\lambda}{\kappa}R_{\kappa\beta\gamma\delta}).}\label{sol_M}
\end{align}

\subsection{Uniqueness of the solution}

We now address the uniqueness to the non-trivial solution \eqref{sol_M} to the constraint \eqref{BH_scalar_constraint} derived above.
If one adds an additional piece to our Ansatz, it will satisfy an homogeneous equation since all source terms have been cancelled by the Ansatz. Demonstrating uniqueness of the non-trivial solution \eqref{sol_M} therefore amounts to prove that
\begin{align}
    \nabla_{(\mu}N_{\nu|(\alpha\beta)|\rho)}=0\label{homogeneous_eqn}
\end{align}
does not admit any non-trivial solution in Kerr spacetime. We call such a tensor field a Young tableau $(2,2)$ Killing tensor. A trivial Killing tensor is defined as a Killing tensor which is given by a cross-product. Such a trivial Killing tensor would add to the quadratic conserved quantity a product of conserved quantities that are already defined. There is only trivial Killing tensor of symmetry type $(2,2)$ namely $N_{\mu\nu\alpha\beta}=Y_{\mu\nu}Y_{\alpha\beta}$ which correspond to add the product $(Q_Y)^2$ defined in \eqref{QY} to the quadratic conserved quantity.
We checked explicitly by solving the partial differential equations analytically using a \textit{Mathematica} notebook that no non-trivial such tensor exists in a perturbative series expansion in $a$ around $a=0$ assuming that it only depends upon $r$ and $\theta$.

\subsection{Summary of the results}
Let us summarize the results we have obtained about the quadratic invariants. Our discussion can be compactified in the two following propositions:
\begin{main_result}
The quadratic invariant
\begin{align}
Q^{(2)}_\text{BH}=Q_R+\qty[-g_{\alpha\gamma}\qty(\xi_\beta \xi_\delta-\frac{1}{2}g_{\beta\delta}\xi^2)+\frac{1}{2}Y\tdu{\alpha}{\lambda}\qty(Y\tdu{\gamma}{\kappa}R_{\lambda\beta\kappa\delta}+\frac{1}{2}Y\tdu{\lambda}{\kappa}R_{\kappa\beta\gamma\delta})]S^{\alpha\beta}S^{\gamma\delta}\label{quadratic_rudiger}
\end{align}
is conserved for the MPD equations at second order in the spin magnitude for spin-induced quadrupole, \textit{i.e.} $\dot Q^{(2)}_\text{BH}=\mathcal O\qty(\mathcal S^3)$ provided that $\delta\kappa=0$, \textit{i.e.} if the test body possesses the multipole structure of a black hole. Here, $Q_R=  K_{\mu\nu}p^\mu p^\nu+L_{\mu\nu\rho}S^{\mu\nu}p^\rho$ with $L_{\mu\nu\rho}$ given in Eq. \eqref{R2} is Rüdiger's quadratic invariant \cite{doi:10.1098/rspa.1983.0012,Compere:2021kjz}.
\end{main_result}

\begin{preliminary_result}
Any tensor $N_{\alpha\beta\gamma\delta}$ possessing the same algebraic symmetries than the Riemann tensor and satisfying the constraint equation
\begin{align}
\hat\nabla N=\Upsilon_\text{NS},
\end{align}
where the source term $\Upsilon_\text{NS}$ is given in Eq. \eqref{source_NS}
will give rise to a quantity
\begin{align}
    Q^{(2)}=Q^{(2)}_\text{BH}+\delta\kappa M_{\alpha\beta\gamma\delta}S^{\alpha\beta}S^{\gamma\delta},\qquad M_{\alpha\beta\gamma\delta}=\,^*\!N^*_{\alpha\beta\gamma\delta}
\end{align}
which is conserved up to second order in the spin parameter for the spin-induced quadrupole MPD equations \eqref{MPD}, \textit{i.e.} $\dot Q^{(2)}=\mathcal O\qty(\mathcal S^3)$, regardless to the value taken by $\delta\kappa$.
\end{preliminary_result}

We notice that $\Upsilon_\text{NS}\eval_{a=0}=0$ in the Schwarzschild case by explicit evaluation of \eqref{source_NS}. A direct consequence is that the deformation of Rüdiger's quadratic invariant constructed in the black hole case is still quasi-conserved for arbitrary $\kappa$:
\begin{main_result}
In Schwarschild spacetime ($a=0$), the deformation of Rüdiger's quadratic invariant $Q_\text{BH}^{(2)}$
given in Eq. \eqref{quadratic_rudiger}
is still conserved for the MPD equations up to $O(\mathcal S^3)$ corrections for arbitrary ($\kappa \in \mathbb R$)  spin-induced quadrupole.
\end{main_result}
Notice that the conservation does not hold for Rüdiger's linear invariant $Q_Y$. The Kerr case $a \neq 0$ will be further discussed in Section \ref{sec:NS}.

\section{Neutron star case around Kerr: a no-go result}
\label{sec:NS}

We summarize in Table \ref{tab:sources} the three constraint equations discussed previously. They all take the form
\begin{align}
    \hat\nabla N=\Upsilon,
\end{align}
with $\Upsilon$ being of grading $[s,p]^+$. It implies that $N$ should be of grading $[s,p-1]^+$.

\begin{table}[!hbt]
\begin{tabular}{c|c|c|c|c}
\textsc{problem} & \textsc{source term $\Upsilon$} & \textsc{equation} & $\Upsilon\eval_{a=0}$ & \textsc{source grading}\\\hline
NS linear invariant & $\Upsilon_\text{lin}$ & \eqref{lin_cst} & $\neq 0$ & $[2,4]^+$ \\\hline
BH quadratic invariant & $\Upsilon_\text{BH}$ & \eqref{source_BH} & $\neq 0$& $[2,3]^+$ \\\hline
NS quadratic invariant & $\Upsilon_\text{NS}$ & \eqref{source_NS} & $0$ & $[2,3]^+$
\end{tabular}
\caption{Source terms for the various constraint equations $\hat\nabla N=\Upsilon$ studied in the paper.}
\label{tab:sources}
\end{table}

Let $\qty{K_\mathfrak{a}}$ be a basis of linearly independent and dimensionless functions build from the (manifestly real) functions $\mathcal A, \mathcal P^2,\mathcal S^2,\Re A,\Im A,\ldots,\Re C,\Im C,D,\ldots,I$. Given that $N$ is dimensionless and given the structure of the source terms $\Upsilon$, we propose the following Ansatz,
\begin{align}
    N=\sum_\mathfrak{a}\sum_{(k,l)\in\mathbb Z^2}K_\mathfrak{a}M^{l}\qty(C_\mathfrak{a}^{(k,l)}f_\mathfrak{a}^{(k,l)}(J)\alpha_1^{(k,k+l)}+D_\mathfrak{a}^{(k,l)}g^{(k,l)}_\mathfrak{a}(J)\omega_1^{(k,k+l)}).\label{generic_ansatz}
\end{align}
Here, $C_\mathfrak{a}^{k,l}$ and $D_\mathfrak{a}^{k,l}$ are numerical coefficients and $f_\mathfrak{a}^{k,l}(J)$ and $g_\mathfrak{a}^{k,l}(J)$ are smooth functions of $J$.

We can work with dimensionless quantities by first introducing the dimensionless variables
\begin{align}
    \tilde r=\frac{r}{M},\qquad\tilde a=\frac{a}{M}.
\end{align}
We notice that the $K_\mathfrak{a}$'s are left unchanged and do not depend anymore on $M$, whereas $\mathcal R\triangleq M\tilde {\mathcal R}$, with $\tilde {\mathcal R}\triangleq \tilde r+i\tilde a$. This yields
\begin{align}
    \alpha^{(n,p)}_1=M^{n-p}\Re\qty(\frac{\bar{\tilde{\mathcal{R}}}}{\tilde{\mathcal{R}}})\triangleq M^{n-p}\tilde{\alpha}^{(n,p)},\qquad
    \omega^{(n,p)}_1=M^{n-p}\Im\qty(\frac{\bar{\tilde{\mathcal{R}}}}{\tilde{\mathcal{R}}})\triangleq M^{n-p}\tilde{\omega}^{(n,p)}.
\end{align}
Each derivative of the term present in the Ansatz scales as $M^{-1}$ times a manifestly dimensionless quantity. All the source terms appearing earlier can be written as $\Upsilon=M^{-1}\tilde\Upsilon$, with $\tilde{\Upsilon}$ being an dimensionless quantity.
This implies that Eq. \eqref{generic_ansatz} reduces to
\begin{align}
   N=\sum_\mathfrak{a}\sum_{(k,l)\in\mathbb Z^2} K_\mathfrak{a}\qty(C_\mathfrak{a}^{k,l}f_\mathfrak{a}^{k,l}(J)\tilde\alpha^{(k,k+l)}+D_\mathfrak{a}^{k,l}g^{k,l}_\mathfrak{a}(J)\tilde\omega^{(k,k+l)}),\label{generic_ansatz_dimless}
\end{align}
which contain only terms that are explicitly independent of $M$. We can further define $\tilde \nabla = M \hat \nabla$ the dimensionless derivative operator and the constraints take the dimensionless form $\tilde\nabla N = \tilde \Upsilon$.

\subsection{Perturbative expansion in $a$ of the constraint equations}

Instead of addressing the non-linear problem in $a$ we will perform a perturbative series in $a$. For any smooth function $f$ of $a$, we define
\begin{align}
    \qty(f)_n\triangleq\dv[n]{f}{a}\eval_{a=0}.
\end{align}
The constraint equation then becomes an infinite hierarchy of equations
\begin{align}
\qty(\tilde{\nabla} N)_n=\tilde\Upsilon_n,\qquad\forall n\geq 0.
\end{align}

Let us describe the $n=0$ and $n=1$ equations.
Since $J\propto a^2$, the functions $f_\mathfrak{a}^{k,l}(J)$ and $g_\mathfrak{a}^{k,l}(J)$ do not contribute  and can be set to one without loss of generality.

\paragraph{$n=0$ equation.} Noticing the identities
\begin{align}
    &\tilde\omega^{(k,k+l)}\eval_{a=0}=\tilde\alpha^{(k,k+l+1)}_A\eval_{a=0}=\tilde\alpha^{(k-1,k+l)}_{\bar A}\eval_{a=0}=0,
    \\
    &\tilde\alpha^{(k,k+l)}\eval_{a=0}=\tilde r^{-l},\quad\tilde\omega^{(k,k+l+1)}_A\eval_{a=0}=-p_r \tilde r^{-(l+1)},
    \quad\tilde\omega^{(k-1,k+l)}_{\bar A}\eval_{a=0}=p_r \tilde r^{-(l+1)}
\end{align}
and making use of Eq. \eqref{deriv_ao}, the $n=0$ constraint becomes
\begin{align}
    \sum_\mathfrak{a}\sum_{(k,l)\in\mathbb Z^2} C_\mathfrak{a}^{(k,l)}\qty[\qty({\tilde\nabla K_\mathfrak{a}})_0\tilde{r}^{-l}-l \qty(K_\mathfrak{a})_0 p_r \tilde r^{-(l+1)}]=\qty(\tilde\Upsilon)_0.
\end{align}
It does not depend on the terms involving $\omega$'s contributions. Moreover, denoting
\begin{align}
C^{(l)}_\mathfrak{a}\triangleq\sum_{k\in\mathbb Z}C^{(k,l)}_\mathfrak{a},
\end{align}
this equation can be further simplified to
\begin{align}
\sum_\mathfrak{a}\sum_{l\in\mathbb Z} C^{(l)}_\mathfrak{a}\qty[\qty({\tilde\nabla K_\mathfrak{a}})_0-l \qty(K_\mathfrak{a})_0 p_r \tilde r^{-1}]\tilde r^{-l}=\qty(\tilde\Upsilon)_0.
\label{eq_LO}
\end{align}

\paragraph{$n=1$ equation.} Following an identical procedure and denoting
\begin{align}
    D_\mathfrak{a}^{(l)}\triangleq\sum_{k\in\mathbb Z}\qty(2k+l)D_\mathfrak{a}^{(k,l)},
\end{align}
the $n=1$ constraint equation can be shown to take the form
\begin{align}
    &\sum_\mathfrak{a}\sum_{l\in\mathbb Z}\bigg\lbrace C_\mathfrak{a}^{(l)}\qty[\qty({\tilde\nabla K_\mathfrak{a}})_1-\qty(K_\mathfrak{a})_1lp_r\tilde r^{-1}]\tilde r^{-l}\nonumber\\
    &-D_\mathfrak{a}^{(l)}\qty[\qty({\tilde\nabla K_\mathfrak{a}})_0x+\qty(K_\mathfrak{a})_0\qty(p_\theta-\qty(l+1)p_rx\tilde r^{-1})]\tilde r^{-(l+1)}\bigg\rbrace=(\tilde\Upsilon)_1.\label{eq_NLO}
\end{align}

\paragraph{Numerical evaluation.}
Eqs. \eqref{eq_LO} and \eqref{eq_NLO} have been numerically evaluated using \textit{Mathematica} in order to try to fix the values of the coefficients $C_\mathfrak{a}^{(l)}$ and $D_\mathfrak{a}^{(l)}$ that would enable a possible solution to the neutron star cases.  We have only looked for ``polynomial'' solutions to these equations, \textit{i.e.} solutions for which the coefficients of the ansatz are non-vanishing only over a finite interval $[l_\text{min},l_\text{max}]$. Given the size of the expressions involved, the only computationally reasonable solving method available to us was the following: let us denote $N$ the number of terms present in the left-hand side of \eqref{eq_LO} (resp. \eqref{eq_NLO}) for a given $[l_\text{min},l_\text{max}]$. Eq. \eqref{eq_LO} (resp. \eqref{eq_NLO}) was then evaluated $N+1$ times at different random values of its variables and parameters, resulting into a linear system of $N+1$ algebraic equations in $N$ variables (the coefficients $C_\mathfrak{a}^{(l)}$ and $D_\mathfrak{a}^{(l)}$) that was then solved using the built-in numerical equation solver of \textit{Mathematica}.

This procedure has been proof-tested by reproducing the coefficients corresponding to the black hole quadratic invariant from the source term $\Upsilon_\text{BH}$. It has then been used to attempt to find a solution to both the linear and the quadratic neutron star problems, with $[l_\text{min},l_\text{max}]=[-10,100]$. No solution has been found, discarding \textit{a priori} the existence of polynomial-type solutions.

In the \textit{Mathematica} notebooks appended as supplementary material, the interval of $l$ is reduced to $[l_\text{min},l_\text{max}]=[0,5]$ in order to reduce the computational time for the interested reader. The four notebooks related to this section are:
\begin{itemize}[label=$\diamond$]
    \item \texttt{\lstinline{Quadratic_BH_LO_final.nb}}: check of the numerical evaluation of the $n=0$ equation: reproduction of the black hole quadratic invariant from $\Upsilon_\text{BH}$;
    \item \texttt{\lstinline{Linear_NS_LO_final.nb}}: attempt of finding a polynomial solution to the $n=0$ equation for the neutron star linear problem (source term $\Upsilon_\text{lin}$);
    \item \texttt{\lstinline{Quadratic_BH_NLO_final.nb}}: check of the numerical evaluation of the $n=1$ equation: reproduction of the black hole quadratic invariant from $\Upsilon_\text{BH}$;
    \item \texttt{\lstinline{Quadratic_NS_NLO_final.nb}}: attempt of finding a polynomial solution to the $n=1$ equation for the neutron star quadratic problem (source term $\Upsilon_\text{NS}$).
\end{itemize}

\section{Discussion and outlook}
\label{sec:outlooks}

The main results obtained in this paper are as follows. At second order in the spin magnitude and for spin-induced quadrupole with black hole-type coupling ($\kappa=1$): (i) the linear Rüdiger invariant $Q_Y$ is still quasi-conserved and (ii) a deformation of Rüdiger's quadratic invariant $Q_\text{BH}^{(2)}$ exists such that the deformed Rüdiger's quadratic invariant is also quasi-conserved. Finally, (iii) the quasi-conservation of the deformed quadratic invariant can only be extended to arbitrary coupling ($\kappa\neq 1$) around the  Schwarszchild spacetime ($a=0$).

All our attempts to find solutions to the constraint equations in the case of an arbitrary spin-induced coupling in generic Kerr spacetime have failed. Let us notice that, even if someone would succeed in solving them, the quasi-invariants so obtained would not be of direct astrophysical  interest. This arises from the fact that, except in the special case where the test body is itself a black hole, the spin-induced term is not the only contribution to the quadrupole. Tidal-type contributions will also arise which will break the quasi-conservation obtained for spin-induced quadrupoles only.

Various extensions of this work would be interesting to explore. One could investigate whether deformations of the quasi-invariants that have studied in this paper still exist at higher orders in the multipolar expansion for black hole-type couplings. Another possible direction would be to understand the link between the existence of these quasi-conserved quantities and the separability of the associated Hamilton-Jacobi equation at second order in the spin magnitude, thus pushing the analysis of Witzany \cite{Witzany_2019} to the next order in the multipole expansion. This would also enable to compute the corresponding shifts in the fundamental frequencies of the action-angle variables description of the finite size particle, which are of direct relevance for modeling EMRIs involving spinning secondaries.  It would also be enlightening to explore the relationship between the new constants for test black holes in Kerr found in this paper and those for arbitrary-mass-ratio binary black holes at second-post-Newtonian order (including the spin-induced quadrupole effects) found by Tanay \emph{et al.}\ in Ref.~\cite{Tanay:2020gfb}, ensuring the integrablity of the system at that order.  Finally, if the structure of the 4-dimensional Kerr covariant building blocks can be generalized to higher dimensions and more generic spacetimes, it could bring new insights on hidden symmetries of such spacetimes \cite{Frolov:2017kze}.

\section*{Acknowledgments}
A. D. is a Research Fellow and G. C. is Senior Research Associate of the F.R.S.-FNRS. G. C. acknowledges support from the FNRS research credit No. J003620F and the IISN convention No. 4.4503.15.  JV is grateful to Lars Andersson, \'Eanna Flanagan, Abraham Harte and Leo Stein for helpful discussions, and especially to Steffen Aksteiner for confirmation via spinorial methods of the conservation of the generalized Carter constant for a quadrupolar test black hole in Kerr.