Analytical model of precessing binaries using post-Newtonian theory
in the extreme mass-ratio limit I: General Formalism

Consider the PN spin precession equations, which to second order in the spins of the individual bodies include the spin-orbit, spin-spin, and quadrupole-monopole couplings Khan et al. (2019); Kesden et al. (2015); Gerosa et al. (2015); Chatziioannou et al. (2017). To the leading PN order, the equations take the form

	$\displaystyle\dot{\vec{S}}_{1}$	$\displaystyle=\vec{\Omega}_{1}\times\vec{S}_{1}\,,$		(1)
	$\displaystyle\dot{\vec{S}}_{2}$	$\displaystyle=\vec{\Omega}_{2}\times\vec{S}_{2}\,,$		(2)
	$\displaystyle\dot{\hat{L}}$	$\displaystyle=-L^{-1}\left(\dot{\vec{S}}_{1}+\dot{\vec{S}}_{2}\right)\,,$		(3)

where

$$\vec{\Omega}_{1}=\vec{\Omega}^{\rm SO}_{1}+\vec{\Omega}^{\rm SS}_{1}+\vec{\Omega}^{\rm QM}_{1}\,,$$

(4)

with

	$\displaystyle\vec{\Omega}^{\rm SO}_{1}$	$\displaystyle=\frac{1}{r^{3}}\left(2+\frac{3}{2}\frac{m_{2}}{m_{1}}\right)L\hat{L}\,,$		(5)
	$\displaystyle\vec{\Omega}^{\rm SS}_{1}$	$\displaystyle=\frac{1}{2r^{3}}\left[S_{2}-3\left(\vec{S}_{2}\cdot\hat{L}\right)\hat{L}\right]\,,$		(6)
	$\displaystyle\vec{\Omega}^{\rm QM}_{1}$	$\displaystyle=-\frac{3}{2r^{3}}\frac{m_{2}}{m_{1}}\left(\vec{S}_{1}\cdot\hat{L}\right)\hat{L}\,.$		(7)

In the above expressions, $L=\mu\sqrt{Mr}$, the masses of the two BHs are $(m_{1},m_{2})$, $\mu=m_{1}m_{2}/M$ is the reduced mass, and $M=m_{1}+m_{2}$ is the total mass. The expression for $\vec{\Omega}_{2}$ is obtained by taking $1\leftrightarrow 2$ in Eqs. (4)-(7). The spins can be written in terms of dimensionless quantities as $\vec{S}_{A}=m_{A}^{2}\vec{\chi}_{A}$, with $A=(1,2)$. Eq. (3) is due to the fact that the total angular momentum, $\vec{J}=L\hat{L}+\vec{S}$ (where $\vec{S}=\vec{S}_{1}+\vec{S}_{2}$ is the total spin vector) is conserved on the precession timescale.

Now, consider taking the EMRI limit of these equations. Let the larger object have mass and spin $(m_{1},\vec{S}_{1})$, while the smaller object has mass and spin $(m_{2},\vec{S}_{2})$. The limit is given by $m_{1}\gg m_{2}$, with $q=m_{2}/m_{1}\ll 1$. From this,

	$\displaystyle M$	$\displaystyle=(1+q)m_{1}\,,\qquad\mu=\frac{qm_{1}}{1+q}\,,$		(8)
	$\displaystyle L$	$\displaystyle=\frac{qm_{1}^{2}}{v}\,,\qquad\vec{S}_{2}=q^{2}m_{1}^{2}\vec{\chi}_{2}=q^{2}\vec{\sigma}_{2}\,,$		(9)
	$\displaystyle\vec{J}$	$\displaystyle=\vec{S}_{1}+\frac{qm_{1}^{2}}{v}\hat{L}+q^{2}\vec{\sigma}_{2}\,,$		(10)

where $v=(M/r)^{1/2}$ is the orbital velocity, and $\vec{\sigma}_{2}=\vec{\chi}_{2}m_{1}^{2}$ is a “renormalized” spin vector. The inclusion of quadrupole-monopole interactions introduces an additional constant of motion on the precession timescale, specifically

	$\displaystyle M^{2}\chi_{\rm eff}$	$\displaystyle=\left(1+\frac{m_{2}}{m_{1}}\right)\left(\vec{S}_{1}\cdot\hat{L}\right)+\left(1+\frac{m_{1}}{m_{2}}\right)\left(\vec{S}_{2}\cdot\hat{L}\right)\,,$		(11)
		$\displaystyle=\left(\vec{S}_{1}\cdot\hat{L}\right)+q\left[\left(\vec{S}_{1}+\vec{\sigma}_{2}\right)\cdot\hat{L}\right]+q^{2}\left(\vec{\sigma}_{2}\cdot\hat{L}\right)\,,$		(12)

where we have recast all expressions to highlight all $q$ factors. Inserting these into Eqs. (1)-(3), and expanding in $q\ll 1$ gives

	$\displaystyle\frac{d\vec{\chi}_{1}}{d\tau}$	$\displaystyle=\left[q\vec{\Omega}_{1}^{(1)}+q^{2}\vec{\Omega}_{1}^{(2)}+{\cal{O}}(q^{3})\right]\times\vec{\chi}_{1}\,,$		(13)
	$\displaystyle\frac{d\vec{\chi}_{2}}{d\tau}$	$\displaystyle=\left[\vec{\Omega}_{2}^{(0)}+{\cal{O}}(q)\right]\times\vec{\chi}_{2}\,,$		(14)
	$\displaystyle\frac{d\hat{L}}{d\tau}$	$\displaystyle=\left[\vec{\Omega}_{L}^{(0)}+q\vec{\Omega}_{L}^{(1)}+{\cal{O}}(q^{2})\right]\times\hat{L}$		(15)

with $\tau=t/m_{1}$, and

	$\displaystyle\vec{\Omega}_{1}^{(1)}$	$\displaystyle=v^{5}\left(2-\frac{3}{2}v\chi_{\rm eff}\right)\hat{L}\,,$		(16)
	$\displaystyle\vec{\Omega}_{1}^{(2)}$	$\displaystyle=-v^{5}\left(\frac{9}{2}-3v\chi_{\rm eff}\right)\hat{L}+\frac{v^{6}}{2}\vec{\chi}_{2}$		(17)
	$\displaystyle\vec{\Omega}_{2}^{(0)}$	$\displaystyle=\frac{3}{2}v^{5}\left(1-v\chi_{\rm eff}\right)\hat{L}+\frac{v^{6}}{2}\vec{\chi}_{1}$		(18)
	$\displaystyle\vec{\Omega}_{L}^{(0)}$	$\displaystyle=\frac{v^{6}}{2}\left(4-3v\chi_{\rm eff}\right)\vec{\chi}_{1}\,,$		(19)
	$\displaystyle\vec{\Omega}_{L}^{(1)}$	$\displaystyle=\frac{3}{2}v^{6}\left(-3+2v\chi_{\rm eff}\right)\vec{\chi}_{1}+\frac{3}{2}v^{6}\left(1-v\chi_{\rm eff}\right)\vec{\chi}_{2}\,,$		(20)

We have used the definition of $\chi_{\rm eff}$ in Eq. (11) to replace all instances of $\vec{S}_{1}\cdot\hat{L}$. It is worth noting that Eqs. (13)-(15) are expanded to different orders in $q$. The reason for this is two-fold: 1) it is a necessity to employ MSA when we study solutions to the precession equations in the EMRI limit in Sec. III, and 2) it is necessary to ensure that $\vec{J}$ is conserved on the precession timescale, due to the different scalings of the spin vectors with $q$ in Eq. (10). Indeed, it is straightforward to check that $d\vec{J}/dt={\cal{O}}(q^{3})$, and thus $\vec{J}$ is conserved at this order in the mass ratio (neglecting radiation reaction).

At this stage, there is a useful simplification that can be performed, by replacing $\vec{\chi}_{1}$ with $\vec{J}/m_{1}^{2}$. The expansion of $\vec{J}$ in the EMRI limit is given by Eq. (10). Using this expansion to replace $\vec{\chi}_{1}$ in Eqs. (19)-(20) will produce a remainder ${\cal{O}}(q^{2})$, since the linear term in the mass ratio is proportional to $\hat{L}$ and $\hat{L}\times\hat{L}=0$. Such ${\cal{O}}(q^{2})$ remainder can be neglected, being higher order than what we are presently considering. Likewise, performing the same procedure with Eq. (18) yields a remainder ${\cal{O}}(q)$, which can also be neglected. Thus, to the order we are working in $q$, we may replace $\vec{\chi}_{1}$ with $\vec{J}/m_{1}^{2}$ without loss of accuracy. Doing so completely eliminates the precessing $\vec{\chi}_{1}$ from Eqs. (14)-(15), and replaces it with the fixed $\vec{J}$. As a result, the evolution of $\chi_{1}$ decouples from the other angular momenta. The latter, namely $\vec{\chi}_{2}$ and $\hat{L}$, are the only quantities we need to take into account when solving for the precession equations in the EMRI Limit.

Gravitational waveforms from spin precessing binaries were first studied in the PN approximation in Apostolatos et al. (1994), with the first frequency domain templates being computed by Lang & Hughes Lang and Hughes (2006) using the stationary phase approximation (SPA) Bender and Orszag (1999). Standard SPA breaks down for spin precessing waveforms Klein et al. (2014), and is not suitable for actual searches. However, the Lang & Hughes waveform provides the necessary qualitative features to capture the effects of spin precession of the frequency domain waveform, and results in conservative estimates on parameter estimation Lang et al. (2011); Yagi and Tanaka (2010). As a matter of simplicity, we limit ourselves to this waveform model. We remark however that the formalism and the analysis of the precession equations, contained in Sec. III & IV.5, is general enough to be used in any precessing waveform.

In the frequency domain, the Lang & Hughes waveform for spin precessing quasi-circular binaries takes the form

$$\tilde{h}_{I}(f)=\sqrt{\frac{5}{96}}\frac{{\cal{M}}^{5/6}}{\pi^{2/3}D_{L}}A_{I}(f)f^{-7/6}e^{i\Phi_{I}(f)}$$

(21)

with frequency $f$, chirp mass ${\cal{M}}=\mu^{3/5}M^{2/5}$, luminosity distance $D_{L}$, and the subscript $I$ labelling the two possible data streams from the three arms of the LISA detector²²2While our analysis is independent of the specific GW detector, for concreteness we will focus on LISA, for which EMRIs are a primary target.. The amplitude function $A_{I}(f)$ is a slowly varying function of frequency through the precession orbital angular momentum of the binary and the motion of the LISA constellation. For the present analysis, the phase is more important than the amplitude, but explicit expression for the latter can be found in Eq. (2.29) of Lang and Hughes (2006). The waveform phase $\Phi(f)$ is a sum of multiple contributions, specifically

$$\Phi(f)=\Psi(f)-\varphi_{{\rm pol},I}(f)-\varphi_{D}(f)-\delta\Phi(f)\,,$$

(22)

with $\Psi(f)$ the Fourier phase arising from the evolution of orbital quantities under radiation reaction, $\varphi_{\rm pol}(f)$ the polarization phase due to the time-(frequency-)varying polarization basis of the GWs, $\varphi_{D}$ the Doppler phase due to the motion of the LISA constellation, and $\delta\Phi$ the waveform precession phase due to the precessing orbital angular momentum. The first and last of these contributions are most important to the present analysis, and are the focus of the remainder of this paper. Expressions for the polarization and Doppler phases can be found in Eqs. (2.30) & (2.31) in Lang and Hughes (2006), respectively.

In the time domain, the Fourier phase $\Psi$ is given by

$$\Psi(t)=2\pi ft-2\phi(t)$$

(23)

with $\phi(t)$ the orbital phase, and for any given value of $f$. The factor of two multiplying the orbital phase comes from the leading PN approximation of GWs, specifically the quadrupole approximation. Coordinate time $t$ and orbital phase $\phi$ may be computed in TaylorF2 approximates in terms of $v=(2\pi MF)^{1/3}$, with $F$ the orbital frequency. Application of the stationary phase approximation to Eq. (23) enforces $f=2F$. Thus $v=(\pi Mf)^{1/3}$, and the Fourier phase takes the form

$$\Psi(f)=2\pi ft_{c}-2\phi_{c}-\frac{\pi}{4}+\frac{3}{128\eta v^{5}}\left[1+\sum_{n}\psi_{n}v^{n}\right]\,.$$

(24)

The PN coefficients of the phase are well known for generic mass binaries, see for example Damour et al. (2001). We will provide expressions in the EMRI limit in Sec. IV.3.

The evolution of the waveform precession phase is directly related to the precession equation of the orbital angular momentum through Eq. (28) of Apostolatos et al. (1994), specifically

$$\frac{d\delta\Phi}{dt}=\frac{(\hat{L}\cdot\vec{N})}{1-(\hat{L}\cdot\vec{N})^{2}}\left(\hat{L}\times\vec{N}\right)\cdot\frac{d\hat{L}}{dt}\,,$$

(25)

where $\vec{N}$ is the line of sight from the detector to the source. Since the LISA constellation is not fixed relative to the Solar System’s barycenter, $\vec{N}$ is generally a function of time. Further, because of the dependence on $\hat{L}$, the solution for $\delta\Phi$ is intricately tied to the solutions to the precession equations. Solving for these equations is the main goal of this work, and the detailed procedure will be explained in the next sections. Once one obtains the solution within the time domain, the waveform precession phase can be computed in the frequency domain by application of the stationary phase condition.

To leading order in MSA, the evolution of the orbital angular momentum is given by Eq. (29). The presence of the cross product therein couples the $x$- and $y$-components of $\hat{L}^{(0)}$ together. However, these can be decoupled by re-casting the equations in second order form, specifically

$\displaystyle\frac{\partial^{2}\hat{L}_{x}^{(0)}}{\partial\tau^{2}}$

$\displaystyle=-\omega_{L}^{2}\hat{L}_{x}\,,\qquad\frac{\partial^{2}\hat{L}_{y}^{(0)}}{\partial\tau^{2}}=-\omega_{L}^{2}\hat{L}_{y}\,,$

(33)

The general solutions are

	$\displaystyle\hat{L}_{x}^{(0)}(\tau,\tilde{\tau})$	$\displaystyle=X(\tilde{\tau})\cos[\alpha(\tau)]+Y(\tilde{\tau})\sin[\alpha(\tau)]\,,$		(34)
	$\displaystyle\hat{L}_{y}^{(0)}(\tau,\tilde{\tau})$	$\displaystyle=X(\tilde{\tau})\sin[\alpha(\tau)]-Y(\tilde{\tau})\cos[\alpha(\tau)]\,,$		(35)

where $(X,Y)$ are purely functions of the long timescale $\tilde{\tau}$, and the precession angle $\alpha$ is defined by

$$\frac{d\alpha}{d\tau}=\omega_{L}\,,$$

(36)

which reduces to $\alpha(\tau)=\omega_{L}\tau$ since $\omega_{L}$ is a constant in the absence of radiation reaction. However, once radiation reaction is included, one needs to directly integrate Eq. (36) to obtain the full evolution of $\alpha(\tau)$, i.e.

$$\alpha(v)=\alpha_{c}+\int dv\frac{\omega_{L}(v)}{dv/d\tau}\,.$$

(37)

The $z$-component can trivially be solved to obtain

$$\hat{L}_{z}^{(0)}(\tilde{\tau})=Z(\tilde{\tau})$$

(38)

which does not depend on the short timescale $\tau$. The functions $(X,Y,Z)$ cannot be determined at this order in the MSA, but due to the normalization of the orbital angular momentum obey

$$[X(\tilde{\tau})]^{2}+[Y(\tilde{\tau})]^{2}+[Z(\tilde{\tau})]^{2}=1\,.$$

(39)

To leading order, the precession equation for $\vec{\chi}_{2}^{(0)}(\tau,\tilde{\tau})$ is given in Eq. (30), with the frequencies $\omega_{SL}$ and $\omega_{SJ}$ given in Eq. (III). The original precession equations possess a sufficient number of constants of motion to define a co-precessing reference frame, and the same holds true at ${\cal{O}}(q^{0})$. We can exploit this fact to define a frame where the unit orbital angular momentum is fixed to be

$\displaystyle\hat{L}^{(0)}_{\rm co-p}=\left[X(\tilde{\tau}),Y(\tilde{\tau}),Z(\tilde{\tau})\right]\,,$

(40)

while $\hat{J}$ still defines the z-axis of the frame. To define a suitable basis in the co-precessing frame, allow two basis vectors to be $\hat{J}$ and

	$\displaystyle\hat{P}$	$\displaystyle=\frac{\hat{L}^{(0)}\times\hat{J}}{|\hat{L}^{(0)}\times\hat{J}|}$
		$\displaystyle=\frac{\left[X\sin\alpha-Y\cos\alpha,-X\cos\alpha-Y\sin\alpha,0\right]}{\sqrt{1-Z^{2}}}\,,$		(41)

the latter of these is orthogonal to the $LJ$-plane. A third can then be defined by

	$\displaystyle\hat{Q}$	$\displaystyle=\frac{\hat{J}\times\hat{P}}{|\hat{J}\times\hat{P}|}$
		$\displaystyle=\frac{\left[X\cos\alpha+Y\sin\alpha,X\sin\alpha-Y\cos\alpha,0\right]}{\sqrt{1-Z^{2}}}\,,$		(42)

which, along with $\hat{L}$, spans the $LJ$-plane. The secondary’s spin vector can then be written as

$\displaystyle\vec{\chi}_{2}^{(0)}$

$\displaystyle=\chi_{2,P}(\tau)\hat{P}+\chi_{2,Q}(\tau)\hat{Q}+\chi_{2,J}(\tau)\hat{J}\,,$

(43)

with the components satsifying

	$\displaystyle\frac{\partial\chi_{2,P}}{\partial\tau}$	$\displaystyle=\omega_{R}(\tilde{\tau})\chi_{2,J}-\omega_{Q}(\tilde{\tau})\chi_{2,Q}\,,$		(44)
	$\displaystyle\frac{\partial\chi_{2,Q}}{\partial\tau}$	$\displaystyle=\omega_{Q}(\tilde{\tau})\chi_{2,P}\,,$		(45)
	$\displaystyle\frac{\partial\chi_{2,J}}{\partial\tau}$	$\displaystyle=-\omega_{R}(\tilde{\tau})\chi_{2,P}\,.$		(46)

where

	$\displaystyle\omega_{R}(\tilde{\tau})$	$\displaystyle=\omega_{SL}\sqrt{[X(\tilde{\tau})]^{2}+[Y(\tilde{\tau})]^{2}}\,,$		(47)
	$\displaystyle\omega_{Q}(\tilde{\tau})$	$\displaystyle=\omega_{SJ}-\omega_{L}+\omega_{SL}Z(\tilde{\tau})\,.$		(48)

This system can be decoupled to obtain a single equation for $\chi_{2,P}$, specifically

$$\frac{d^{2}\chi_{2,P}}{d\tau^{2}}+\omega_{P}^{2}(\tilde{\tau})\chi_{2,P}=0\,,$$

(49)

where

$\displaystyle\omega_{P}(\tilde{\tau})$

$\displaystyle=\sqrt{[\omega_{Q}(\tilde{\tau})]^{2}+[\omega_{R}(\tilde{\tau})]^{2}}\,,$

(50)

The solution is then

$$\chi_{2,P}=\chi_{P,0}\cos[\gamma_{P}(\tau)]+\dot{\chi}_{P,0}\sin[\gamma_{P}(\tau)]\,,$$

(51)

with $[\chi_{P,0},\dot{\chi}_{P,0}]$ set by initial conditions, and the spin angle $\gamma_{P}$ is defined through

$$\frac{\partial\gamma_{P}}{\partial\tau}=\omega_{P}(\tilde{\tau})\,.$$

(52)

Much like the precession angle $\alpha$, $\gamma_{P}$ must be integrated under the effect of radiation reaction. In general, $[\chi_{P,0},\dot{\chi}_{P,0}]$ are functions of the long time variable $\tilde{\tau}$. However, at the order in $q$ we are working, this dependence is irrelevant and they may be taken to be constants.

The remaining two components satisfy equations that depend on $v$, even after they are transformed to $\gamma_{P}$ being the dependent variable, i.e.

	$\displaystyle\frac{d\chi_{2,Q}}{d\gamma_{P}}$	$\displaystyle=-\frac{\omega_{Q}(v)}{\omega_{P}(v)}\chi_{2,P}(\gamma_{P})\,,$		(53)
	$\displaystyle\frac{d\chi_{2,J}}{d\gamma_{P}}$	$\displaystyle=-\frac{\omega_{R}(v)}{\omega_{P}(v)}\chi_{2,P}(\gamma_{P})\,.$		(54)

To solve these exactly, one would have to solve integrals of the form

$\displaystyle\chi_{2,J}\sim\int d\tau\;\omega_{R}[v(\tau)]e^{i\gamma_{P}(\tau)}$

(55)

which are non-trivial. These can be evaluated by realizing that $\omega_{P}$ in Eq. (52) is always positive definite, and thus, there are no critical points of the phase $\gamma_{P}$. Then, the integral can be approximated by Laplace’s method via repeated integration by parts, the first application of which produces a correction of order

$$\frac{1}{\omega_{P}}\frac{dv}{d\tau}\frac{\partial\omega_{SL}}{\partial v}\sim{\cal{O}}(v^{8})$$

(56)

while the leading order term is of order $\omega_{R}/\omega_{P}\sim 1+{\cal{O}}(v)$. Thus, $\chi_{2,J}$ can be solved to high precision by assuming $\nu_{R}=\omega_{R}/\omega_{P}$ does not evolve appreciably on the radiation reaction timescale, with the result being

$\displaystyle\chi_{2,J}$

$\displaystyle=\chi_{J,0}+\nu_{R}\left[-\chi_{P,0}\sin\gamma_{P}+\dot{\chi}_{P,0}\cos\gamma_{P}\right]+{\cal{O}}(v^{8})\,,$

(57)

where $\chi_{J,0}$ is an integration constant. Likewise,

$$\chi_{2,Q}=\chi_{Q,0}+\nu_{Q}\left[\chi_{P,0}\sin\gamma_{P}-\dot{\chi}_{P,0}\cos\gamma_{P}\right]+{\cal{O}}(v^{7})\,,$$

(58)

with $\nu_{Q}=\omega_{Q}/\omega_{P}$, and the components of the secondary’s spin in the co-precessing frame are now solved.

As a final point before proceeding, the solutions in Eqs. (51) & (57)-(58) possess four integration constants, but $\vec{\chi}_{2}$ has only three components. One of these is redundant, which can be discovered by inserting these solutions and Eq. (43) into the original precession equations in Eq. (30). This is only satisfied if $\chi_{J,0}=\omega_{Q}\chi_{Q,0}/\omega_{R}$, thus eliminating the anomalous degree of freedom.

We are now left with solving Eq. (III), which must simultaneously give us $\hat{L}^{(1)}(\tau)$ and $[X(\tilde{\tau}),Y(\tilde{\tau}),Z(\tilde{\tau})]$. Rearranging terms, we may write

	$\displaystyle\frac{\partial\hat{L}^{(1)}}{\partial\tau}-\omega_{L}\hat{e}_{z}\times\hat{L}^{(1)}$	$\displaystyle=-\frac{\partial\hat{L}^{(0)}}{\partial\tilde{\tau}}-\omega_{L}^{(1)}\hat{e}_{z}\times\hat{L}^{(0)},$
		$\displaystyle+v\omega_{SL}\vec{\chi}^{(0)}_{2}\times\hat{L}^{(0)}\,.$		(59)

The left hand side above indicates that $\hat{L}^{(1)}$ will oscillate with phase variable $\alpha$ on the short timescale. The artificial resonance arises due to the right hand side of the above equation containing two types of terms, specifically those that oscillate solely with $\alpha$, and those that couple both $\alpha$ and $\gamma_{P}$. The former of these are what cause the artificial resonance, and can be seen if one re-writes the source term above, which only depends on $\hat{L}^{(0)}$ and $\vec{\chi}_{2}^{(0)}$, into a harmonic decomposition. More specifically

$$\frac{\partial\hat{L}^{(0)}}{\partial\tilde{\tau}}=\hat{C}_{1}(\tilde{\tau})\cos\alpha+\hat{S}_{1}(\tilde{\tau})\sin\alpha+Z^{\prime}(\tilde{\tau})\hat{e}_{z}\,,$$

(60)

with

	$\displaystyle\hat{C}_{1}(\tilde{\tau})$	$\displaystyle=[X^{\prime}(\tilde{\tau}),-Y^{\prime}(\tilde{\tau}),0]\,,$
	$\displaystyle\hat{S}_{1}(\tilde{\tau})$	$\displaystyle=[Y^{\prime}(\tilde{\tau}),X^{\prime}(\tilde{\tau}),0]$		(61)

where ^′ corresponds to differentiation with respect to $\tilde{\tau}$. Further,

	$\displaystyle\hat{e}_{z}\times\hat{L}^{(0)}$	$\displaystyle=\hat{C}_{2}(\tilde{\tau})\cos\alpha+\hat{S}_{2}(\tilde{\tau})\sin\alpha\,,$		(62)
	$\displaystyle\hat{C}_{2}(\tilde{\tau})$	$\displaystyle=[Y(\tilde{\tau}),X(\tilde{\tau}),0]\,,$
	$\displaystyle\hat{S}_{2}(\tilde{\tau})$	$\displaystyle=[-X(\tilde{\tau}),Y(\tilde{\tau}),0]\,.$		(63)

Lastly, we may write

	$\displaystyle\vec{\chi}_{2}^{(0)}\times\hat{L}^{(0)}$	$\displaystyle=\sum_{j,k=-1}^{1}\hat{E}_{jk}(\tilde{\tau})e^{i(j\alpha+k\gamma_{P})}$
		$\displaystyle=\hat{C}_{3}(\tilde{\tau})\cos\alpha+\hat{S}_{3}(\tilde{\tau})\sin\alpha$
		$\displaystyle+\sum_{j=-1}^{1}\sum_{k\neq 0}\hat{E}_{jk}(\tilde{\tau})e^{i(j\alpha+k\gamma_{P})}$		(64)

where

	$\displaystyle\hat{C}_{3}(\tilde{\tau})$	$\displaystyle=\chi_{Q,0}{\cal{F}}(\tilde{\tau})\left[Y(\tilde{\tau}),X(\tilde{\tau}),0\right]\,,$		(65)
	$\displaystyle\hat{S}_{3}(\tilde{\tau})$	$\displaystyle=\chi_{Q,0}{\cal{F}}(\tilde{\tau})\left[-X(\tilde{\tau}),Y(\tilde{\tau}),0\right]\,,$		(66)
	$\displaystyle{\cal{F}}(\tilde{\tau})$	$\displaystyle=\frac{\omega_{Q}(\tilde{\tau})}{\omega_{R}(\tilde{\tau})}-\frac{Z(\tilde{\tau})}{\sqrt{[X(\tilde{\tau})]^{2}+[Y(\tilde{\tau})]^{2}}}\,.$		(67)

The $[\hat{C}_{j}(\tilde{\tau}),\hat{S}_{j}(\tilde{\tau}),Z^{\prime}(\tilde{\tau})]$ produce the artificial resonance and, thus, the sum total of them must vanish in order to remove this. Thus, the equations that $[X(\tilde{\tau}),Y(\tilde{\tau}),Z(\tilde{\tau})]$ must satisfy are

	$\displaystyle Z^{\prime}(\tilde{\tau})$	$\displaystyle=0\,,$		(68)
	$\displaystyle X^{\prime}(\tilde{\tau})$	$\displaystyle=\omega_{XY}(\tilde{\tau})Y(\tilde{\tau})\,,$		(69)
	$\displaystyle Y^{\prime}(\tilde{\tau})$	$\displaystyle=-\omega_{XY}(\tilde{\tau})X(\tilde{\tau})\,,$		(70)

with

$$\omega_{XY}(\tilde{\tau})=-\omega_{L}^{(1)}+v\chi_{Q,0}\frac{\omega_{SL}}{\omega_{R}(\tilde{\tau})}(\omega_{SJ}-\omega_{L})\,.$$

(71)

Eq. (68) implies that $Z$ is a constant, and thus, $[\omega_{Q},\omega_{R},\omega_{XY}]$ are also constant as a result of the normalization condition in Eq. (39). As such, we will drop the dependence on $\tilde{\tau}$ from the quantities. Choosing $Z=\hat{L}_{z,0}$, the solutions for $X(\tilde{\tau})$ and $Y(\tilde{\tau})$ are

	$\displaystyle X(\tilde{\tau})$	$\displaystyle=\hat{L}_{x,0}\cos[\lambda(\tilde{\tau})]-\hat{L}_{y,0}\sin[\lambda(\tilde{\tau})]\,,$		(72)
	$\displaystyle Y(\tilde{\tau})$	$\displaystyle=-\hat{L}_{y,0}\cos[\lambda(\tilde{\tau})]-\hat{L}_{x,0}\sin[\lambda(\tilde{\tau})]\,.$		(73)

where $\lambda(\tilde{\tau})$, which we call the renormalization angle, is defined by

$$\frac{d\lambda}{d\tilde{\tau}}=\omega_{XY}\,,$$

(74)

and $[\hat{L}_{x,0},\hat{L}_{y,0},\hat{L}_{z,0}]$ are constants set by initial data. The function $\hat{L}^{(0)}(\tau,\tilde{\tau})$ is now completely determined and we have removed the artificial resonance in Eq. (III).

We are now left with the task of determining $\hat{L}^{(1)}$, or more specifically its dependence on the short timescale $\tau$. The only remaining source on the right hand side of Eq. (III.3) is the double sum in Eq. (III.3), i.e.

$$\frac{\partial\hat{L}^{(1)}}{\partial\tau}-\omega_{L}\hat{e}_{z}\times\hat{L}^{(1)}=v\omega_{SL}\hat{\cal{S}}(\tau,\tilde{\tau})$$

(75)

with

$$\hat{\cal{S}}(\tau,\tilde{\tau})=\sum_{j}\sum_{k\neq 0}\hat{E}_{jk}(\tilde{\tau})e^{i[j\alpha(\tau)+k\gamma_{p}(\tau)]}\,.$$

(76)

It is simplest to consider each component of this separately. The z-component of the above equation decouples, and after integrating, produces

$$\hat{L}^{(1)}_{z}=\frac{v\omega_{SL}}{\omega_{P}}\sqrt{1-\hat{L}_{z,0}^{2}}\left[\chi_{P,0}\sin\gamma_{P}+\dot{\chi}_{P,0}\left(1-\cos\gamma_{P}\right)\right]\,.$$

(77)

The $x-$ and $y-$ components are coupled in the same way as in Sec. III.1. Decoupling them produces

$\displaystyle\frac{\partial^{2}\hat{L}_{x,y}^{(1)}}{\partial\tau^{2}}+\omega_{L}^{2}\hat{L}_{x,y}^{(1)}$

$\displaystyle=v\omega_{SL}{\cal{T}}_{x,y}(\tau,\tilde{\tau})$

(78)

where

	$\displaystyle{\cal{T}}_{x}(\tau,\tilde{\tau})$	$\displaystyle=\frac{\partial}{\partial\tau}\hat{\cal{S}}_{x}(\tau,\tilde{\tau})-\omega_{L}\hat{\cal{S}}_{y}(\tau,\tilde{\tau})$		(79)
	$\displaystyle{\cal{T}}_{y}(\tau,\tilde{\tau})$	$\displaystyle=\frac{\partial}{\partial\tau}\hat{\cal{S}}_{y}(\tau,\tilde{\tau})+\omega_{L}\hat{\cal{S}}_{x}(\tau,\tilde{\tau})\,,$		(80)

The sources ${\cal{T}}_{x,y}$ can be written in harmonics of $\alpha$ and $\gamma_{P}$, specifically

	$\displaystyle{\cal{T}}_{x,y}$	$\displaystyle=\Omega_{+}\left\{C_{x,y}^{(+)}(\tilde{\tau})\cos[\delta_{+}(\tau)]+S_{x,y}^{(+)}(\tilde{\tau})\sin[\delta_{+}(\tau)]\right\}$
		$\displaystyle+\Omega_{-}\left\{C_{x,y}^{(-)}(\tilde{\tau})\cos[\delta_{-}(\tau)]+S_{x,y}^{(-)}(\tilde{\tau})\sin[\delta_{-}(\tau)]\right\}\,,$		(81)

with $\Omega_{\pm}=2\omega_{L}\pm\omega_{P}$, $\delta_{\pm}=\alpha\pm\gamma_{P}$, and the functions

	$\displaystyle C_{x}^{(\pm)}(\tilde{\tau})=S_{y}^{(\pm)}(\tilde{\tau})$	$\displaystyle=A_{\pm}\left[-\chi_{P,0}Y(\tilde{\tau})\mp\dot{\chi}_{P,0}X(\tilde{\tau})\right]\,,$		(82)
	$\displaystyle S_{x}^{(\pm)}(\tilde{\tau})=-C_{y}^{(\pm)}(\tilde{\tau})$	$\displaystyle=A_{\pm}\left[\pm\dot{\chi}_{P,0}Y(\tilde{\tau})-\chi_{P,0}X(\tilde{\tau})\right]\,,$		(83)

where

$$A_{\pm}=\frac{L_{z,0}(\omega_{P}\pm\omega_{Q})\pm\omega_{SL}\mp L_{z,0}^{2}\omega_{SL}}{2\omega_{P}\sqrt{1-L_{z,0}^{2}}}$$

(84)

are constants on the precession timescale. It is straightforward to show that the solutions to Eq. (78) are

	$\displaystyle\hat{L}_{x,y}^{(1)}(\tau,\tilde{\tau})$	$\displaystyle=\ell_{x,y}\cos[\alpha(\tau)]+\dot{\ell}_{x,y}\sin[\alpha(\tau)]$
		$\displaystyle+\frac{v\omega_{SL}}{\omega_{P}}\left\{C_{x,y}^{(-)}(\tilde{\tau})\cos[\delta_{-}(\tau)]+S_{x,y}^{(-)}(\tilde{\tau})\sin[\delta_{-}(\tau)]\right\}$
		$\displaystyle-\frac{v\omega_{SL}}{\omega_{P}}\left\{C_{x,y}^{(+)}(\tilde{\tau})\cos[\delta_{+}(\tau)]+S_{x,y}^{(+)}(\tilde{\tau})\sin[\delta_{+}(\tau)]\right\}$		(85)

where $[\ell_{x,y},\dot{\ell}_{x,y}]$ are constants set by initial conditions. In general, these will be functions of the long timescale $\tilde{\tau}$, which would then be determined by proceeding to next order in our MSA. Since we are truncating our analysis at linear order $q$, we take these to be constants instead.