# Secular dipole-dipole stability of magnetic binaries

Christopher Aykroyd, Adrien Bourgoin, Christophe Le Poncin-Lafitte, Stéphane Mathis, Marie-Christine Angonin

► To cite this version:

Christopher Aykroyd, Adrien Bourgoin, Christophe Le Poncin-Lafitte, Stéphane Mathis, Marie-Christine Angonin. Secular dipole-dipole stability of magnetic binaries. Astronomy and Astrophysics - A&A, 2023, 675, pp.A32. 10.1051/0004-6361/202346171. hal-04126223

HAL Id: hal-04126223

https://hal.science/hal-04126223

Submitted on 9 Nov 2023

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Secular dipole-dipole stability of magnetic binaries

C. Aykroyd<fnref n="1" />, A. Bourgoin<fnref n="1" />,<fnref n="2" />, C. Le Poncin-Lafitte<fnref n="1" />, S. Mathis<fnref n="2" />, and M.-C. Angonin<fnref n="1" />

Received 17 February 2023 / Accepted 11 May 2023

## ABSTRACT

**Context.** The presence of strong large-scale stable magnetic fields in a significant portion of early-type stars, white dwarfs, and neutron stars is well established. Despite this, the origins of these fields remain a subject of ongoing investigation, with theories including fossil fields, mergers, and shear-driven dynamos. One potential key for understanding the formation of these fields could lie in the connection between magnetism and binarity. Indeed, magnetism can play a significant role in the long-term orbital and precessional dynamics of binary systems. In gravitational wave astronomy, the advanced sensitivity of upcoming interferometric detectors such as LISA and the *Einstein Telescope* will enable the characterisation of the orbital inspirals of compact systems, including their magnetic properties. A comprehensive understanding of the dynamics of magnetism in these systems is necessary for the interpretation of the gravitational wave signals and to avoid bias in the calibration of instruments. This knowledge can additionally be used to create new magnetic population models and provide insight into the nature and origins of their internal magnetic fields.

**Aims.** The aim of this study is to investigate the secular spin precession dynamics of binary systems under pure magnetic dipole-dipole interactions, with a focus on stars with strong, stable, and predominantly dipolar fields.

**Methods.** We employed an orbit-averaging procedure for the spin precession equations from which we derived an effective secular description. By minimising the magnetic interaction energy of the system, we obtained the configurations of spin equilibrium and their respective stabilities. Finally, we also derived a set of conditions required for the validity of our assumptions to hold.

**Results.** We show that among the four states of equilibrium, there is a single secular state that is globally stable, corresponding to the configuration where the spin and magnetic axes of one star are reversed with respect to the companions', and orthogonal to the orbital plane. Our results are compared to traditional methods of finding instantaneous states of equilibrium, in which orbital motion is generally neglected. Finally, we provide analytical solutions in the neighbourhood of the stable configuration, which can be used to derive secular orbital evolution in the context of gravitational wave astronomy.

**Key words.** magnetic fields – celestial mechanics – binaries: close – stars: neutron – white dwarfs – stars: early-type

## 1. Introduction

Close to 10% of early-type massive main sequence (MMS) stars host stable large-scale magnetic fields, ranging from $3 \times 10^2$ to $3 \times 10^4$ G (Grunhut et al. 2013; Ferrario et al. 2015b; Shultz et al. 2019). Meanwhile, it is estimated that 20–25% of the white dwarf (WD) population is magnetic (Bagnulo & Landstreet 2021), with detected fields between $10^3$ and $10^9$ G. For neutron stars (NS), surface field strengths are gauged by the order of $10^8–10^{13}$ G in classical radio pulsars (Phinney & Kulkarni 1994), and reach up to $10^{14}–10^{15}$ G in magnetars (Kaspi & Beloborodov 2017). The origins of magnetic fields are highly debated for both main-sequence stars and for compact objects; below, we present the prevailing theories for each respective class of stars.

In main-sequence late-type stars, it is considered established that external magnetic fields are driven by a dynamo action in the outer convective zone (Brun & Browning 2017; Brun et al. 2022). Conversely, this channel is less likely to be the case for hot, massive stars ($M > 1.5 M_\odot$), which maintain a radiative envelope and inner convective core: any dynamo-based explanation must resolve the challenge of transporting the magnetic field towards the outer surface faster than the stellar evolution timescale (Charbonneau & MacGregor 2001). In fact, in radiative stars, magnetic fields are believed to decay in diffusivity timescales estimated longer than their host's main-sequence

lifetime (Moss 2003). Early direct evidence via Zeeman spectropolarimetry has long since shown that chemically peculiar Ap and Bp stars, which represent around 10% of early-type A/B stars (Ferrario et al. 2015b), host strong secularly stable magnetic fields, with strengths uncorrelated with stellar rotation – as should be the case for dynamo-fed fields – and geometry largely captured by oblique dipole rotor models (see e.g. Stibbs 1950; Borra et al. 1982). These fields have therefore been suggested to have ‘fossil’ origin: fields that are remnants of a prior stellar evolution stage and are effectively frozen into the plasma (Moss 1987; Mathis et al. 2011). The exact field formation process is a topic of debate, but multiple plausible mechanisms have been proposed – ranging from accumulated magnetic flux captured from the interstellar cloud at birth, to protostar mergers and pre-main-sequence dynamos (Ferrario et al. 2009). More recent studies such as the B fields in OB stars (BOB; Morel et al. 2014) and the Magnetism in Massive Stars (MiMeS; Wade et al. 2015; Shultz et al. 2019) surveys also confirmed compatible magnetic incidence and properties on more general O/B-type massive stars. Numerical and semi-analytical magneto-hydrodynamic computations have established the existence of long-term-stable internal field configurations consistent with non-convective bodies such as radiative stars, NS, and WDs, which favours the fossil field scenario (Braithwaite & Spruit 2004; Braithwaite & Nordlund 2006; Braithwaite 2008;

<footnote n="1">¹ SYRTE, Observatoire de Paris, Université PSL, CNRS, Sorbonne Université, LNE, 61 avenue de l'Observatoire, 75014 Paris, France
e-mail: christopher.aykroyd@obspm.fr</footnote>

<footnote n="2">² Université Paris-Saclay, Université Paris-Cité, CEA, CNRS, AIM, 91191 Gif-sur-Yvette, France</footnote>
Duez et al. 2010; Fuller & Mathis 2023). These field configurations are composed of both toroidal and poloidal components that stabilise each other (Tayler 1980; Braithwaite 2009; Akgün et al. 2013); outside the star, the toroidal field is attenuated and mainly the poloidal component is visible. These results were found to reproduce the general characteristics of observations: the roughly off-centre dipolar structure, the independence from stellar spin, and finally, the strong field amplitude (Braithwaite & Spruit 2004; Duez 2011). Nevertheless, the fossil field hypothesis is not without its challenges. For example, only a small fraction of MMS stars host observable fields, with a sharp dearth of weak-field objects; the precise mechanism for field formation, stability, and evolution must explain this cutoff. It has been suggested that there are thresholds to field strength below which shear or convection instabilities develop (see e.g. Aurière et al. 2007; Gaurat et al. 2015; Jouve et al. 2020; Jermyn & Cantiello 2021), or that, in some stars, the time needed to reach an equilibrium becomes longer than the age in the main sequence, due to the Coriolis force produced by rapid rotation (Braithwaite & Cantiello 2013). An additional challenge to fossil fields is the existence of a great scarcity of magnetism in close binaries, as low as 2% incidence (Carrier et al. 2002; Alecian et al. 2014). In fact, there is a single known doubly magnetic close binary up to date, the ε-Lupi system (Pablo et al. 2019). If the fossil scenario is indeed to be the main field formation channel, it is plausible to expect a similar magnetic incidence in binaries and in single stars. However, Vidal et al. (2019) suggests that tidal instabilities in binary pairs can disrupt the magnetic fields via turbulent Joule diffusion within a few millions years, potentially explaining the scarcity of strong-field binaries. Alternatively, it has been argued that interstellar clouds with strong magnetic fields are harder to fragment (Commerçon et al. 2010), yielding selection biases towards less magnetic binary systems. Other alternatives have also been suggested to address this challenge, such as the merger scenarios (Ferrario et al. 2009, 2015b; Schneider et al. 2016, 2019). In these scenarios, coalescing main sequence stars and/or protostars would generate strong enough shear to drive dynamo action, yielding a single magnetic byproduct star. Such hypotheses are in line with the prediction that around 8% of MMS stars originate from mergers (de Mink et al. 2014), and naturally explain the lack of magnetic binaries. Nevertheless, at this stage, no channel can be completely favoured over another.

In the compact object community there is a somewhat analogous debate. On one side, classical fossil theories defend that magnetic white dwarfs (MWDs) and NS are derived from Ap/B and O-type stars respectively, and that their fields must persist from the main sequence or red giant phase (Tout et al. 2004; Ferrario & Wickramasinghe 2005; Wickramasinghe & Ferrario 2005). Another possibility lifted by Stello et al. (2016); Bagnulo & Landstreet (2021) is that internal dynamos in the convective cores of intermediate-mass and massive stars – externally invisible during the main sequence phase – might develop into strong stable fields by flux compression as the stellar core collapses into a WD. These fields would then be slowly revealed as the WD sheds its outer layers, and decay in secular Ohmic timescales. On the other side of the debate, merger theories (Tout et al. 2008; Ferrario et al. 2020) advocate an intimate link between magnetism and binarity; in such scenarios, two common-envelop stars would generate a magnetic field through differential rotation and merge to form a strong-field MWD. Closely interacting systems which failed to completely merge might instead develop into magnetic cataclysmic variables. Finally, alternative theories support the operating of some dynamo mechanism during the cooling of the WD – for example,

during the crystallisation convection of the core in a rapidly rotating WD (Isern et al. 2017; Schreiber et al. 2021). In support of the fossil hypothesis in WDs and NS is the striking similarity in magnetic flux between this group and the main sequence A/B/O stars, as well as the long field decay timescales – estimated to be of the order of tens or hundreds of billions of years (Cumming 2002; Ferrario et al. 2015a). Conversely, the progenitors of non-magnetic WDs would be low-mass stars, which are known to harbor relatively weak dynamo-driven fields. This contrasting spectral origin is consistent with the observation that MWDs are on average more massive than their non-magnetic counterparts (Liebert 1988; Bagnulo & Landstreet 2021). However, merger scenarios would also naturally explain such disparity. Against the fossil field hypothesis, it has been argued that there is an insufficient volume density of Ap/Bp stars to by itself account for the high occurrence of MWDs, as required for the classical fossil theory to hold (Kawka & Vennes 2004; Ferrario & Wickramasinghe 2005; Bagnulo & Landstreet 2021). Furthermore, many surveys have pointed out a sparsity of known detached binaries composed of a MWD plus a non-degenerate companion (Liebert et al. 2005; Ferrario et al. 2015a), whereas conversely, magnetism amongst cataclysmic variables is ubiquitous, with about a quarter of these WDs reaching the high-field range ($B \ge 1$ MG). This has propelled the suggestion that magnetism and binarity in WD systems are intrinsically connected, leading to the advent of merger hypotheses. However, as pointed out by Landstreet & Bagnulo (2020); Bagnulo & Landstreet (2021), at the present moment, at least five such detached MWD binary systems are known, and this frequency may be higher than previously thought.

Evidently, current observational data are insufficient to completely rule out one or another formation channel, and indeed, multiple of them may be at work simultaneously. Further studies of magnetic binary interactions can provide crucial insights to resolve this debate, of which binarity has shown to be a key element.

In particular, magnetism can play an important role in the dynamics of stellar, compact, and planetary systems, shaping the long-term evolution of their orbits (Bourgoin et al. 2022; Bromley & Kenyon 2022). In gravitational wave (GW) astronomy, the upcoming generation of interferometric detectors such as the Laser Interferometer Space Antenna (LISA; see Amaro-Seoane et al. 2017) and the Einstein Telescope (ET; see Maggiore et al. 2020) will provide enough sensitivity to probe the interactions of compact systems and to characterise their magnetic attributes. On one hand, this can enable the composition of new magnetic population models and bring insight to the nature and origin of internal fields. On the other hand, a careful understanding of the dynamics of magnetism in these systems is also required to avoid biases in the calibration of the instrument and in the interpretation of signals into physical parameters. Indeed, the secular (long-term) impact of magnetism on the orbits will manifest as a definitive signature on the GWs, which must be correctly accounted for (Bourgoin et al. 2022; Carvalho et al. 2022; Lira et al. 2022; Savalle et al., in prep.).

It is thus imperative to study the secular evolution of the fields themselves, their binary coupling, and the interplay with stellar orientation. In the case of stable, rigid fields, this translates to investigating the rotational dynamics of the stars, which may include their states of equilibrium. In purely tidal-driven systems, spin motion and stability has long-since been determined by Hut (1980, 1981). In the case of star-planet systems, Damiani & Lanza (2015) further explored the interaction between tides and magnetic braking – where pressure-driven
stellar winds give rise to the loss of angular momentum – and Strugarek et al. (2017) determined the relative strengths of the tidal and magnetic effects in magnetic star-planet interactions. The long-term effects of static dipole fields on stellar rotation, however, has yet to be completely explored. In this regard, the works of Pablo et al. (2019) and King et al. (1990) investigate the spin equilibrium under purely the dipolar interactions, but they neglect the interactions with orbital dynamics. In particular, they explore the cases where the obliquity between the dipole and spin axes is constrained to 0° (aligned) and 90° (perpendicular) respectively.

In this work, we direct our attention towards magnetic binary interactions – in particular, towards stars with strong, stable, and predominantly dipolar fields. We consider the magnetic moments of these stars to be aligned with the stellar spin, and investigate the secular evolution of each star’s orientation due to the mutual magnetic torques, through an effective orbit-averaged description. We provide criteria for determining whether magnetism dictates the stellar rotational motion, which may then reflect on the secular evolution of the binary’s orbits. Our study can be applied to any type of star system (MMS, WD, and NS constituents), as long as both components of the binary are magnetic and dominated by dipolar terms. In this way, it can be useful to gather further understanding of the formation processes of MMS stars as well as to ensure an efficient data processing of the LISA or ET observations in the context of compact-star binaries.

The paper is subdivided as follows. The complete physical setup and assumptions of our model are described in Sect. 2. We derive, in Sect. 3, the effective secular orientation dynamics of the spins of each star due to dipole-dipole interactions. We then provide in Sect. 4 an analysis of the equilibrium states and their respective stability, developing a simple analytical solution for spin precession which is valid for quasi-stable systems. Our results are verified numerically and applied to a system possibly satisfying our requirements (Sect. 5). Finally, we compare the contrasting results with respect to the traditional instantaneous equilibrium (Sect. 6), highlighting each method’s advantages and differences.

**Notations and conventions.** We presently introduce the notation used throughout the paper. For each vector **u** element of some vector space U, we represent its norm in light typeface $u = |u|$ and its direction by a hat $\hat{u} = u/u$. Whenever two vectors are parallel, we symbolise this relationship by $u \parallel v$; and similarly, two perpendicular vectors are denoted by $u \perp v$. We represent the vector space dual to U by a starred $U^*$. In this setting, we denote by an underscore $\underline{u} \in U^*$ the associated canonical dot-product covector, that is, the linear functional from U to $\mathbb{R}$ that satisfies $\underline{u} : v \mapsto (u \cdot v)$. For two vector spaces U and V, we represent their Cartesian product by $U \times V = \{ (u, v) ; u \in U, v \in V \}$ and their tensor product by $U \otimes V = \{ u \otimes v ; u \in U, v \in V \}$. Finally, whenever two vector subspaces are disjoint $U \cap V = \{ 0 \}$, their sum is direct and is denoted by $U \oplus V = \{ u + v \mid u \in U, v \in V \}$.

## 2. Magnetic binary model

In this section, we present the physical set up and the assumptions used throughout the paper. Then, we proceed to re-derive the instantaneous magnetically driven precession equations that govern the rotational state of the system.

We consider an isolated binary system of point-like magnetised bodies, dominated by non-relativistic motion. We assign to each body an index $\ell \in \{1, 2\}$, used throughout the paper, and we

<figure><img src="image_2.png" /><figcaption>**Fig. 1.** Binary system in the reference frame of the centre-of-mass. To simplify the drawing, the primary is placed at the origin (CM). We include the spin axis $\hat{s}_\ell$ of each star ($\ell \in \{1, 2\}$) and the basis $e_0 = (\hat{e}_x, \hat{e}_y, \hat{e}_z)$.</figcaption></figure>

refer to them as the ‘primary’ and ‘secondary’, respectively. For each star we introduce a position $\mathbf{x}_\ell$, a mass $m_\ell$, a radius $R_\ell$, a magnetic field $B_\ell$, and an intrinsic angular momentum (or spin) $\mathbf{s}_\ell$. We place ourselves in the reference frame of the centre-of-mass (CM) of the system, to which we attach a right-handed basis $e_0 = (\hat{e}_x, \hat{e}_y, \hat{e}_z)$ spanning the Euclidean tangent space $E_3 \cong \mathbb{R}^3$. The elements of $e_0$ are chosen such that $\hat{e}_x$ points towards the direction of closest approach, $\hat{e}_z$ is orthogonal to the orbital plane, and $\hat{e}_y$ completes the basis. In this frame, the system can be viewed as an effective one-body problem, parametrised by the relative separation $\mathbf{r} = \mathbf{x}_2 - \mathbf{x}_1$. In the absence of non-Keplerian perturbations, the CM frame will be inertial and the elements of $e_0$ will be static. The setup is illustrated in Fig. 1.

We consider a stable magnetic field that is rigidly frozen into each star, compatible with general observations in massive stars and compact objects. The field is assumed to be predominantly dipolar, which captures most observed topologies (see the off-centred dipole model: Borra et al. 1982; Achilleos & Wickramasinghe 1989), although quadrupoles and octupoles have been detected in some cases (Landstreet & Mathys 2000; Wickramasinghe & Ferrario 2000; Beuermann et al. 2007; Donati & Landstreet 2009; Kochukhov et al. 2011; Landstreet et al. 2017). In this scenario, we model $B_\ell$ as a centred dipole, given in function of some point $\mathbf{x}$ outside the surface of the star:

$$B_1(\mathbf{x}) = \frac{\mu_0}{4\pi} \left( \frac{3(\mu_1 \cdot (\mathbf{x} - \mathbf{x}_1))(\mathbf{x} - \mathbf{x}_1)}{|\mathbf{x} - \mathbf{x}_1|^5} - \frac{\mu_1}{|\mathbf{x} - \mathbf{x}_1|^3} \right), \quad (1a)$$

$$B_2(\mathbf{x}) = \frac{\mu_0}{4\pi} \left( \frac{3(\mu_2 \cdot (\mathbf{x} - \mathbf{x}_2))(\mathbf{x} - \mathbf{x}_2)}{|\mathbf{x} - \mathbf{x}_2|^5} - \frac{\mu_2}{|\mathbf{x} - \mathbf{x}_2|^3} \right), \quad (1b)$$

where $\mu_\ell$ is the magnetic dipole moment of each star, and $\mu_0$ is the vacuum permeability. The primary will feel the field $B_2$ of its companion at relative position $\mathbf{x} = -\mathbf{r}$, and the secondary at $\mathbf{x} = \mathbf{r}$. It is convenient to express these fields in terms of the linear map $\mathcal{B}_r : E_3 \to E_3$, which acts on the magnetic dipole moments according to:

$$B_1(\mathbf{r}) = \frac{\mu_0}{4\pi} \mathcal{B}_r(\mu_1), \qquad B_2(-\mathbf{r}) = \frac{\mu_0}{4\pi} \mathcal{B}_r(\mu_2). \quad (2)$$

We note that $\mathcal{B}_r$ can be identified with the symmetric tensor field with values in $E_3 \otimes E_3^*$,

$$\mathcal{B}_r = \frac{1}{r^3} (3 \hat{\mathbf{r}} \otimes \hat{\mathbf{r}} - \mathbf{I}), \quad (3)$$

with $\mathbf{I}$ denoting the identity in $E_3$.

As is observed in the known doubly magnetic MMS binary system $\epsilon$-Lupi (Shultz et al. 2015; Pablo et al. 2019), we examine the particular case where the fields are symmetric about the star’s
<table><caption>**Table 1.** Typical binary system physical parameters (high-field range), and the corresponding dimensionless parameters.</caption><thead><tr><th rowspan="2">Scenario</th><th colspan="6">Physical parameters</th><th colspan="4">Ratios</th></tr><tr><th>m<sub>1</sub>, m<sub>2</sub></th><th>R<sub>1</sub>, R<sub>2</sub></th><th>B<sub>1</sub><sup>p</sup>, B<sub>2</sub><sup>p</sup></th><th>P<sub>1</sub>, P<sub>2</sub></th><th>P<sub>orb</sub></th><th>γ<sup>fig</sup> = Γ<sup>fig</sup>/Γ<sup>B</sup></th><th>γ<sup>tide</sup> = Γ<sup>tide</sup>/Γ<sup>B</sup></th><th>η = F<sup>B</sup>/F<sup>N</sup></th><th>ε = P<sub>orb</sub>/τ</th></tr></thead><tbody><tr><td>MMS</td><td>10 M<sub>☉</sub></td><td>4.5 R<sub>☉</sub></td><td>10<sup>4</sup> G</td><td>5 d</td><td>5 d</td><td>2 × 10<sup>8</sup> J<sub>2</sub><sup>ℓ</sup></td><td>2 × 10<sup>6</sup> k<sub>2</sub><sup>ℓ</sup>/Q<sub>ℓ</sub></td><td>5 × 10<sup>-10</sup></td><td>2 × 10<sup>-9</sup></td></tr><tr><td>WD</td><td>1 M<sub>☉</sub></td><td>10<sup>4</sup> km</td><td>10<sup>9</sup> G</td><td>1 h</td><td>10 h</td><td>2 × 10<sup>6</sup> J<sub>2</sub><sup>ℓ</sup></td><td>7 × 10<sup>-1</sup> k<sub>2</sub><sup>ℓ</sup>/Q<sub>ℓ</sub></td><td>1 × 10<sup>-10</sup></td><td>1 × 10<sup>-7</sup></td></tr><tr><td>NS</td><td>1.4 M<sub>☉</sub></td><td>15 km</td><td>10<sup>15</sup> G</td><td>10 min</td><td>1 h</td><td>6 × 10<sup>5</sup> J<sub>2</sub><sup>ℓ</sup></td><td>7 × 10<sup>-8</sup> k<sub>2</sub><sup>ℓ</sup>/Q<sub>ℓ</sub></td><td>7 × 10<sup>-15</sup></td><td>5 × 10<sup>-7</sup></td></tr></tbody></table>

**Notes.** $\gamma^{\text{fig}} = \Gamma^{\text{fig}}/\Gamma^{\text{B}}$ the ratio between figure effects and magnetic torques (cf. Eq. (9)); $\gamma^{\text{tide}} = \Gamma^{\text{tide}}/\Gamma^{\text{B}}$ the ratio between tidal and magnetic torques (cf. Eq. (11)); $\eta = F^{\text{B}}/F^{\text{N}}$ the ratio between forces (Eq. (18)); and $\epsilon = P_{\text{orb}}/\tau$ the ratio between dynamical timescales (Eqs. (15)–(16)). We assume a binary system with two components of the same kind – that is, MMS-MMS, WD-WD or NS-NS pairs.

axis of rotation, with alignment between the magnetic dipole moment and the spin<fnref n="1" />:

$$\mu_{\ell}(t) = \mu_{\ell} \hat{s}_{\ell}(t). \qquad (4)$$

In practice, the magnetic moment amplitudes may be expressed in terms of the dipolar field evaluated at the poles $B_{\ell}^{p} := B_{\ell}(R_{\ell}\hat{s}_{\ell})$ (cf. Eq. (1)), which can be observationally estimated via spectropolarimetry (see e.g. Shultz et al. 2015). We can thus invert Eq. (1) to deduce:

$$\mu_{\ell} = \frac{2\pi}{\mu_0} B_{\ell}^{p} R_{\ell}^{3}. \qquad (5)$$

The magnetic field of each star will interact with the dipole of the companion, inducing the following torque:

$$\Gamma_{1}^{\text{B}} = \mu_{1} \times \mathbf{B}_{2}(-\mathbf{r}), \qquad (6a)$$

$$\Gamma_{2}^{\text{B}} = \mu_{2} \times \mathbf{B}_{1}(\mathbf{r}). \qquad (6b)$$

Simultaneous contributions due to gravity exist. We introduce the total torque felt by $\ell$,

$$\Gamma_{\ell} = \Gamma_{\ell}^{\text{B}} + \Gamma_{\ell}^{\text{fig}} + \Gamma_{\ell}^{\text{tide}} + (\dots), \qquad (7)$$

where $\Gamma_{\ell}^{\text{B}}$, $\Gamma_{\ell}^{\text{fig}}$ and $\Gamma_{\ell}^{\text{tide}}$ are the respective contributions from the magnetic interaction, figure effects (rigid extended-body interactions), and tides. In order to quantify the relative strength of the dipole-dipole interaction, we introduce the dimensionless parameter $\gamma$:

$$\gamma_{\ell} = \frac{|\Gamma_{\ell} - \Gamma_{\ell}^{\text{B}}|}{\Gamma_{\ell}^{\text{B}}} \le \gamma_{\ell}^{\text{fig}} + \gamma_{\ell}^{\text{tide}} + (\dots), \qquad (8)$$

where $\gamma_{\ell}^{\text{fig}} = \Gamma_{\ell}^{\text{fig}}/\Gamma_{\ell}^{\text{B}}$ and $\gamma_{\ell}^{\text{tide}} = \Gamma_{\ell}^{\text{tide}}/\Gamma_{\ell}^{\text{B}}$ are the contributions to $\gamma$ due to figure effects and due to tides, respectively. Our main interest lies in isolating the equilibrium dynamics of magnetic effects, and we shall therefore consider the regime where $\Gamma_{\ell}^{\text{B}}$ is dominant. More explicitly, we assume $\gamma_{\ell}^{\text{fig}} \ll 1$ and $\gamma_{\ell}^{\text{tide}} \ll 1$, for which we shall presently derive criteria.

The first assumption concerns the strength of figure effects. For a deformed extended body, classical gravitational torques up to quadrupole order have magnitudes roughly around $\Gamma_{\ell}^{\text{fig}} \sim (3/2)(Gm_{1}m_{2}/r) J_{2}^{\ell}(R_{\ell}/a)^{2}$, where $J_{2}^{\ell}$ is the dimensionless quadrupole moment, $a$ is the semi-major axis of the orbit, and

$G$ is the gravitational constant (see e.g. Poisson & Will 2014). The corresponding contribution to $\gamma$ is

$$\gamma_{\ell}^{\text{fig}} = \frac{\Gamma_{\ell}^{\text{fig}}}{\Gamma_{\ell}^{\text{B}}} \sim \frac{3G\mu_{0}}{2\pi} \frac{m_{1}m_{2}}{B_{1}^{p}B_{2}^{p}R_{1}^{3}R_{2}^{3}} R_{\ell}^{2} J_{2}^{\ell}. \qquad (9)$$

Equation (9) shows that the ratio $\gamma_{\ell}^{\text{fig}}$ is mainly scaled by the surface magnetic field strength, sphericity, and mean density of each stellar component. In this work, we shall be considering perfectly spherical stars with $J_{2}^{\ell} = 0$, in which case we formally have no figure torques. In practice, the cutoff to $J_{2}^{\ell}$ below which rotation is driven by magnetism is given by

$$J_{2}^{\ell} \lesssim \frac{2\pi}{3G\mu_{0}} \frac{B_{1}^{p} B_{2}^{p} R_{1}^{3} R_{2}^{3}}{m_{1} m_{2}} \frac{1}{R_{\ell}^{2}}, \qquad (10)$$

and can be used as a criteria for the domain of validity of our models.

Similarly, we turn our attention to tidal interactions, which generate torques that scale with $\Gamma_{\ell}^{\text{tide}} \sim 6(Gm_{1}^{2}R_{\ell}^{5}/a^{6})(k_{2}^{\ell}/Q_{\ell})$, where $m_{\text{m}}$ is the mass of the tide-inflicting body, $k_{2}^{\ell}$ is the gravitational Love number of $\ell$, and $Q_{\ell}$ is the tidal dissipation quality factor (see e.g. Poisson & Will 2014; Strugarek et al. 2017). Then,

$$\gamma_{\ell}^{\text{tide}} = \frac{\Gamma_{\ell}^{\text{tide}}}{\Gamma_{\ell}^{\text{B}}} \sim \frac{6G\mu_{0}}{\pi} \frac{m_{\text{m}}^{2}}{B_{1}^{p} B_{2}^{p} R_{1}^{3} R_{2}^{3}} \frac{R_{\ell}^{5} k_{2}^{\ell}}{a^{3} Q_{\ell}}, \qquad (11)$$

and the corresponding cutoff criteria for $k_{2}^{\ell}/Q_{\ell}$ is:

$$\frac{k_2^\ell}{Q_\ell} \lesssim \frac{\pi}{6G\mu_0} \frac{B_1^p B_2^p R_1^3 R_2^3}{m_m^2} \frac{a^3}{R_\ell^5}. \qquad (12)$$

Table 1 shows typical values for $\gamma_{\ell}^{\text{fig}}$ and $\gamma_{\ell}^{\text{tide}}$ in MMS, WD, and NS systems. At such distance scales ($a \sim 10^{8}-10^{9}$ km), the thresholds for $k_{2}/Q$ are well within the range to allow magnetically driven NS-NS systems, while tidal effects may have dominant contributions in MMS-MMS binaries, and WD-WDs lay somewhere in-between. However, even for the most magnetic systems (NS-NS binaries), in order for the rotational dynamics not to be dominated by figure effects, we require that the quadrupole moment $J_{2}^{\ell}$ be at most of the order of $10^{-6}$.

Regardless, under spherical, rigid star assumptions, the spins will suffer precession due to purely magnetic torques $\Gamma_{\ell}^{\text{B}}$; these torques are orthogonal to the intrinsic angular momentum, and hence the magnitude $s_{\ell}$ must be conserved. Substituting each

<footnote n="1">¹ Throughout the text, we shall assume $\mu_{\ell}$ positive, but $\mu_{\ell} < 0$ is also allowed, with appropriate sign changes in the equations.</footnote>
term and normalising by $s_\ell$ we obtain coordinate-free spin equations:

$$ \frac{d\hat{s}_1}{dt} = -\alpha_1 \mathcal{B}_r(\hat{s}_2) \times \hat{s}_1, \qquad (13a) $$

$$ \frac{d\hat{s}_2}{dt} = -\alpha_2 \mathcal{B}_r(\hat{s}_1) \times \hat{s}_2, \qquad (13b) $$

with $\alpha_\ell = \mu_0\mu_1\mu_2/(4\pi s_\ell)$. We denote Eq. (13) the ‘instantaneous’ precession system. The spin axes are constrained to the unit sphere $S_2$, yielding a total of four degrees of freedom for the coupled system of equations. Each spin precesses around a (time-dependent) axis $\hat{\omega}_\ell$ determined by the field direction $\hat{\omega}_\ell \propto \mathcal{B}_r(\hat{s}_m)$ – where $m \in \{1, 2\}$, $m \neq \ell$ is the index of the companion star – and have Larmor frequencies given by

$$ \omega_\ell = |\alpha_\ell \mathcal{B}_r(\hat{s}_m)| \sim \frac{\alpha_\ell}{r^3}. \qquad (14) $$

The orbital dynamics will induce periodic fluctuations on the separation $r$, causing the precession axis and frequency to oscillate in time with period $P_{\text{orb}}$. It is clear then that two distinct timescales will be involved in the dynamics of spin precession – (1) a timescale corresponding to the orbital period $P_{\text{orb}}$, manifesting in the ‘wobbles’ of the axis $\hat{\omega}_\ell$ and in the modulation of the frequency $\omega_\ell$; as well as (2) a timescale $\tau_\ell$ due to an average precession rate $\langle\omega_\ell\rangle$, which we define as:

$$ \tau_\ell = \frac{2\pi}{\langle\omega_\ell\rangle} \equiv \frac{2\pi b^3}{\alpha_\ell}, \qquad (15) $$

where $b = \sqrt{r_{\min}r_{\max}}$ is the geometrical mean between the separation at the pericentre and at the apocentre of the orbit. In an elliptical orbit, $b$ corresponds to the semi-minor axis.

### 3. Secular precession

We are presently interested in determining the equilibrium configurations of the precession system and their stability. However, the instantaneous equilibrium obtained by equating (13) to zero does not take into account the orbital dynamics; the strong dependence of the torques on the orbital position of each star implies that the configurations of equilibrium may largely fluctuate as the orbit evolves, which occurs in the fast timescale $P_{\text{orb}}$. Indeed, a configuration that was momentarily stable at some point in the orbit may be disrupted by the orbital motion, leading to instability. Conversely, there may exist configurations where the spins oscillate in the fast timescale $P_{\text{orb}}$, but on a longer timescale can be seen to be stable, due to an effective cancellation of the fluctuations. It is therefore in our interest to search for states of secular (long-term) equilibrium and stability. To do this we must eliminate the effects of these oscillatory terms of short period $P_{\text{orb}}$, which can be performed by employing an orbital averaging scheme to obtain the effective dynamics. Variants of this method are widely adopted for determining secular solutions for the orbital motion (Pound & Poisson 2008; Poisson & Will 2014; Gerosa et al. 2015; Will & Maitra 2017; Bourgoin et al. 2022).

For the binary systems considered in this work (magnetic MMS, WDs, and NS), the two timescales ($P_{\text{orb}}$ and $\tau_\ell$) of Eq. (13) will be distinct enough that their effects can be isolated. Indeed, if the torque strength acting on a star is small enough, then its spin axis will not be significantly affected within a single orbital revolution. Intuitively, the impact of the torques on

the spin axis is captured by the spin precession rate $\omega_\ell$. A larger precession rate (or smaller precession timescale $\tau_\ell$) implies faster changes to the axis $\hat{s}_\ell$ due to stronger torques. Thus, one may explicitly compute the characteristic time-ratio between the orbital timescale $P_{\text{orb}}$ and the precession timescale $\tau_\ell$:

$$ \epsilon_\ell = \frac{P_{\text{orb}}}{\tau_\ell} \sim \frac{5\pi}{2G\mu_0} \frac{B_1^p B_2^p R_1^3 R_2^3}{(m_1 + m_2)m_\ell} \frac{1}{R_\ell^2} \frac{P_\ell}{P_{\text{orb}}} (1 - e^2)^{-3/2}, \qquad (16) $$

where $P_\ell$ is the rotational period of body $\ell$ and $e$ the eccentricity of the orbit. To derive the above order of magnitude relation, we have considered as a rough estimate that the mean separation is given by Kepler’s third law $a^3 \sim G(m_1 + m_2)P_{\text{orb}}^2/(4\pi^2)$ and $b \sim a\sqrt{1 - e^2}$. We stress that these relationships are still valid as order-of-magnitude estimates for relativistic systems. We have also estimated the spin magnitude for each star from that of a homogeneous sphere:

$$ s_\ell \sim \frac{4\pi m_\ell R_\ell^2}{5 P_\ell}. \qquad (17) $$

In all three types of systems considered (Table 1), we obtain very low values for $\epsilon$, namely $\epsilon_{\text{MMS}} \sim 10^{-9}$, $\epsilon_{\text{WD}} \sim 10^{-7}$, and $\epsilon_{\text{NS}} \sim 10^{-6}$. We therefore place ourselves in the scenario $\epsilon \ll 1$. As discussed, in this scenario, the spin axes will suffer very little variations within the time-frame of a single orbit. We can therefore consider an effective precession dynamics which averages-out these small orbital oscillations. Conversely, for systems with $\epsilon \gtrsim 1$ the two timescales cannot be decoupled in this manner. In this way, the time-ratio parameter $\epsilon$ can be used as a criteria for the validity of the averaging procedure that follows.

For simplicity of the averaging model, we place ourselves in a classical Newtonian framework, although relativistic corrections are possible. A rough estimate of the impact of magnetic fields on the orbit can be obtained by comparing the magnetic force $F^B \sim 3\mu_0\mu_1\mu_2/(4\pi a^4)$ and the gravitational force $F^N \sim Gm_1m_2/a^2$:

$$ \eta = \frac{F^B}{F^N} \sim \frac{3\pi}{G\mu_0} \frac{B_1^p B_2^p R_1^3 R_2^3}{m_1 m_2} \frac{1}{a^2}. \qquad (18) $$

Even for the most magnetic systems considered, this ratio is limited to $F^B/F^N \lesssim 10^{-10}$ (Table 1). We can therefore consider magnetism negligible in the orbital dynamics, both for MMS stars and compact systems. In this setting, the orbital frame basis $(\hat{e}_x, \hat{e}_y, \hat{e}_z)$ is inertial, the orbits are elliptical, and the separation $\mathbf{r} = \hat{r}\hat{\mathbf{r}}$ can be parametrised as an ellipse in the centre of mass frame:

$$ \hat{r}(f) = \hat{e}_x \cos f + \hat{e}_y \sin f, \qquad r(f) = \frac{a(1 - e^2)}{1 + e \sin f}, \qquad (19) $$

where $f$ is the true anomaly, $a$ the semi-major axis, and $e$ the eccentricity of the orbit. The Keplerian solution determines the relationship between $f$ and $t$ (the time with respect to the reference pericentre passage):

$$ nt \equiv E - e \sin E \mod 2\pi, \qquad (20a) $$

$$ E = \arctan2 \left( e + \cos f, \sqrt{1 - e^2} \sin f \right), \qquad (20b) $$

where $\arctan2$ represents the 2-argument inverse tangent, $E$ is the eccentricity anomaly, and $n = \sqrt{G(m_1 + m_2)/a^3}$ is the mean angular motion.
To formalise our previous argument, we consider the binary system at some given instant $t = t_0 + h$, with a short-timescale variation $h \in [0, P_{\text{orb}}]$. In essence, we are considering that $h$ parametrises a single full orbital revolution of the binary, beginning at some time $t_0$. In this timeframe, the variations in the spin axis are bounded by

$$|\hat{s}_{\ell}(t_0 + h) - \hat{s}_{\ell}(t_0)| = \left| \int_{0}^{h} \frac{d\hat{s}_{\ell}}{dt}(t_0 + u)du \right| \quad (21a)$$

$$\le \left( \sup_{u \in [0, h]} \left| \frac{d\hat{s}_{\ell}}{dt}(t_0 + u) \right| \right) P_{\text{orb}}, \quad (21b)$$

where the supremum of the derivative of $\hat{s}_{\ell}$ can be obtained from the precession Eq. (13), with indices $\ell \in \{1, 2\}$ and $m \in \{1, 2\}$, where $m \neq \ell$:

$$\left| \frac{d\hat{s}_{\ell}}{dt} \right| = \alpha_{\ell} \left| \mathcal{B}_{\mathbf{r}}(\hat{s}_{\mathbf{m}}) \times \hat{s}_{\ell} \right| \le \frac{1}{\tau_{\ell}} \left( \frac{1+e}{1-e} \right)^{3/2}. \quad (22)$$

It is easy to see (cf. Eq. (3)) that equality can be reached when $\hat{s}_{\mathbf{m}} = \hat{\mathbf{r}}$, $\hat{s}_{\ell} \perp \hat{\mathbf{r}}$, and $r = r_{\min} = 1/(a(1-e))$.

Equation (13) may then be evaluated at time $t = t_0 + h$, and using the bounds obtained, expanded via $\hat{s}_{\ell}(t_0 + h) = \hat{s}_{\ell}(t_0) + O(\epsilon_{\ell})$, from whence:

$$\frac{d\hat{s}_{\ell}}{dt}(t_0+h) = -\alpha_{\ell} \mathcal{B}_{\mathbf{r}(t_0+h)}(\hat{s}_{\mathbf{m}}(t_0)) \times \hat{s}_{\ell}(t_0) + \frac{\alpha_{\ell}}{a^3(1-e^2)} O(\epsilon), \quad (23)$$

where we have explicitly included the temporal dependence of $\mathbf{r} = \mathbf{r}(t_0 + h)$ in the subscript of $\mathcal{B}_{\mathbf{r}}$, and defined $\epsilon = \max(\epsilon_1, \epsilon_2)$. Since the above equation is valid for any $h \in [0, P_{\text{orb}}]$, integrating over a full orbit yields the secular spin equation

$$\left\langle \frac{d\hat{s}_{\ell}}{dt} \right\rangle = -\alpha_{\ell} \langle \mathcal{B}_{\mathbf{r}} \rangle (\hat{s}_{\mathbf{m}}) \times \hat{s}_{\ell} + \frac{\alpha_{\ell}}{a^3(1-e^2)} O(\epsilon), \quad (24)$$

where the orbital averaging operator is defined for some function of time $\xi$ as

$$\langle \xi \rangle = \frac{1}{P_{\text{orb}}} \int_{0}^{P_{\text{orb}}} \xi(t_0 + h) dh = \frac{n}{2\pi} \int_{0}^{2\pi} \tilde{\xi}(f_0 + f) \left(\frac{dt}{df}\right) df, \quad (25)$$

where $\tilde{\xi}(f(t)) = \xi(t)$ is the description of $\xi$ in terms of the true anomaly, and the Jacobian factor is found via implicit differentiation of Eq. (20):

$$\frac{df}{dt} = \frac{1}{n} \frac{(1 - e^2)^{3/2}}{(1 + e \cos f)^2}. \quad (26)$$

In the Keplerian scenario the orbits are fixed, and the linear map $\mathcal{B}_{\mathbf{r}}$ – which depends purely on the radial separation $\mathbf{r}$ – is therefore $P_{\text{orb}}$-periodic. Consequently, the orbital average $\langle \mathcal{B}_{\mathbf{r}} \rangle$ will be constant, independent of the secular time $t_0$. Plugging the expressions of $\mathcal{B}_{\mathbf{r}}$ (Eq. (3)) and of the separation $\mathbf{r}$ (Eq. (19)) into Eq. (25), we obtain the effective field tensor in the orbital frame basis $(\hat{e}_{\mathbf{x}}, \hat{e}_{\mathbf{y}}, \hat{e}_{\mathbf{z}})$:

$$\langle \mathcal{B}_{\mathbf{r}} \rangle = \frac{1}{2a^3(1 - e^2)^{3/2}} (\mathbf{I} - 3 \hat{e}_{\mathbf{z}} \otimes \underline{\hat{e}_{\mathbf{z}}}). \quad (27)$$

We note that we have essentially averaged out the orbital oscillations of the magnetic field due to the short-timescale elliptical movement, obtaining a corresponding ‘average field’ operator $\langle \mathcal{B}_{\mathbf{r}} \rangle$. As a consequence, the radial component of $\mathcal{B}_{\mathbf{r}}$ has been

<figure><img src="image_3.png" /><figcaption>Fig. 2. Sample trajectories for the secular evolution of the spin axis of the primary $\hat{s}_1$, plotted against the unit sphere for different values of $\kappa = v_2/v_1$. The initial conditions are fixed, and shown for the primary as a dashed line from the origin to $\hat{s}_1(0)$. The colours are interpolated between blue and red from initial time to $t_0 = 50 v_1^{-1}$ respectively. The black axis represents the direction of $\hat{e}_z$. The behaviour for the secondary is analogous, up to a swap of initial conditions and $\kappa \mapsto \kappa^{-1} = v_1/v_2$.</figcaption></figure>

suppressed, leaving a dipolar effective field with a predominant component in the orthogonal direction $\hat{e}_z$, plus a weaker component aligned with the spin direction.

We denote the normalised term on the right-hand side of Eq. (27) by $\tilde{\mathcal{B}} = \mathbf{I} - 3 \hat{e}_z \otimes \hat{e}_z$. The final secular form of the spin precession equations is obtained by absorbing the constants together:

$$\left\langle \frac{d\hat{s}_1}{dt} \right\rangle = -\nu_1 \tilde{\mathcal{B}}(\hat{s}_2) \times \hat{s}_1, \quad (28a)$$

$$\left\langle \frac{d\hat{s}_2}{dt} \right\rangle = -\nu_2 \tilde{\mathcal{B}}(\hat{s}_1) \times \hat{s}_2, \quad (28b)$$

where we have introduced the magnetic rotational frequencies $\nu_{\ell}$, coupling magnetism, spins, and orbital parameters:

$$\nu_{\ell} = \frac{\mu_0 \mu_1 \mu_2}{4\pi} \frac{1}{s_{\ell}} \frac{1}{2a^3(1 - e^2)^{3/2}}. \quad (29)$$

The above system has four degrees of freedom, two for each spin, since the magnitudes of $\hat{s}_{\ell}$ are conserved quantities of Eq. (28). Additionally, the magnetic interaction energy is also conserved, as will be discussed in Sect. 4.1. By performing a linear transformation on the secular time variable $t_0 \mapsto v_1^{-1} t_0$ it is possible to reduce the system to a single dimensionless parameter $\kappa = v_2/v_1$, which depends only on the ratio between the two spin magnitudes. The parameter $\kappa$ will therefore completely control the dynamics of the system, producing a range of bounded trajectories such as illustrated in Fig. 2. These solutions are roughly epicyclic in nature, described by a predominant precessional motion around the axis $\hat{e}_z$ plus an important nutation component. When the respective frequencies of these two motions align as rational multiples of one another, the solutions become periodic. In the following section, we analyse the states of equilibrium of the secular dynamical system. We then proceed to analyse their stability and approximate solutions for trajectories similar to those presented in Fig. 2.

## 4. Equilibrium states and stability

A ‘secular’ equilibrium state at time $t_0$ is a pair of spins defined on an orbital period $(\hat{s}_1, \hat{s}_2) : [t_0, t_0 + P_{\text{orb}}] \to E_3 \times E_3$ such that the average change in intrinsic angular momentum is zero:

$$\left\langle \frac{d\hat{s}_1}{dt} \right\rangle(t_0) = \left\langle \frac{d\hat{s}_2}{dt} \right\rangle(t_0) = 0. \quad (30)$$
<figure><img src="image_4.png" /><figcaption>Fig. 3. Equilibrium configurations of the spin axes. These correspond to directions that are either inside the orbital plane, in some arbitrary direction (case 1, left column), or orthogonal to the orbital plane (case 2, right column). In each case, the spins may either be parallel (top) or anti-parallel (bottom). As discussed in Sect. 4.2, only the anti-parallel orthogonal configuration (bottom-right) is the only stable scenario.</figcaption></figure>

For a purely dipolar magnetic torque, it is clear from Eq. (28) that this condition can only be achieved when the terms of each cross product are either parallel or zero, which we may write more concisely in the form:

$$ \bar{\mathcal{B}}(\hat{s}_2) = \lambda_1 \hat{s}_1, \qquad \bar{\mathcal{B}}(\hat{s}_1) = \lambda_2 \hat{s}_2. \qquad (31) $$

This corresponds to a singular-value problem, which can be solved with the explicit matrix form of $\bar{\mathcal{B}}$ in the orbital-frame basis. Two classes of solutions can be determined from the singular vectors of $\bar{\mathcal{B}}$, as described below. In either class, the two spins must be parallel with each other, but they may point in the same direction or in reversed directions. Figure 3 illustrates the full set of equilibrium configurations.

### Case 1. Equilibrium configurations in the orbital plane

The first solution to the singular-value problem corresponds to an equilibrium configuration where the pair of spin axes are both contained inside the orbital plane:

$$ \hat{s}_1 = \sigma_1 \hat{p}, \qquad \hat{s}_2 = \sigma_2 \hat{p}, \qquad (32) $$

for some unit vector $\hat{p}$ in the plane of the orbit (i.e. $\hat{p} \cdot \hat{e}_z = 0$). The two unitary parameters $(\sigma_1, \sigma_2) \in \{-1, 1\} \times \{-1, 1\}$ describe the relative alignment between the two spins: they can either point in the same direction or in reversed directions. We observe that the set of all pairs $(\hat{s}_1, \hat{s}_2)$ which satisfy these conditions at a given moment in time forms a space of dimension 1.

### Case 2. Equilibrium configurations orthogonal to the orbital plane

The second solution to the singular-value problem corresponds to a configuration where both spins are orthogonal to the orbital plane – in other words, parallel to the axis $\hat{e}_z$:

$$ \hat{s}_1 = \sigma_1 \hat{e}_z, \qquad \hat{s}_2 = \sigma_2 \hat{e}_z. \qquad (33) $$

As in the previous case, the spins can be oriented in the same direction or in reverse directions according to the values of $\sigma_\ell$.

#### 4.1. Magnetic interaction energy between two dipoles

The interaction energy between two magnetic dipoles is given by the symmetric expression below (see e.g. Pablo et al. 2019):

$$ U_B(\mu_1, \mu_2) = -\frac{\mu_0}{4\pi} \mu_1 \cdot \mathcal{B}_r(\mu_2). \qquad (34) $$

We shall denote $U_B$ the ‘instantaneous magnetic energy’. Conversely, one may take the orbital average of $U_B$, normalising the resulting expression by a positive constant factor:

$$ \bar{U}_B(\hat{s}_1, \hat{s}_2) = -\hat{s}_1 \cdot \bar{\mathcal{B}}(\hat{s}_2). \qquad (35) $$

We denote $\bar{U}_B$ the ‘secular magnetic energy’. As can be straightforwardly verified, $\bar{U}_B$ is a constant of motion of the secular system – its time derivative vanishes for any pair $(\hat{s}_1, \hat{s}_2)$ of axes satisfying (28). Without considering additional forces or dissipation, the orbit-averaged precession system is conservative, and the motion is restricted to some level-curve of constant-energy $\bar{U}_B$. In astrophysical systems, we expect dissipation due to radiation, tidal forces, and internal frictions to bring the magnetic system to the lower energy states, by exchanging energy until eventually settling into a local minimum of $\bar{U}_B$. As we shall see, this local minimum corresponds to a state of stability of the physical system.

In Appendix A, we recall that any symmetric bilinear form constrained to the unit $n$-sphere $U: S_n \times S_n \to \mathbb{R}$ is bounded by its largest-magnitude eigenvalue. Applying the principle to the secular magnetic energy (cf. Eq. (35)), one obtains the bounds $-2 \le \bar{U}_B \le 2$. The spin configurations in which these bounds are actually attained correspond to $\hat{s}_\ell$ along the direction of the associated eigenvector – that is, $\hat{s}_\ell \parallel \hat{e}_z$. In fact, such configuration corresponds exactly to the equilibrium states orthogonal to the orbital plane as determined in the previous section (Eq. (33)). In particular, the energy lower bounds are reached in the anti-parallel case ($\sigma_2 = -\sigma_1$), which as we shall see, is the most stable equilibrium point. In the following section, we analyse the local stability of all the computed equilibrium states from the standpoint of the Hessian form of $\bar{U}_B$.

#### 4.2. Stability tests

The stability of each equilibrium configuration can be analysed via the local convexity of the magnetic interaction energy. We evaluate the nature of each of the extrema – minimum, maximum or saddle point – by determining the sign of the Hessian at that point. We remind the reader of its definition.

Consider a real-valued function $f$ of $n$ real variables $x = (x_1, \dots, x_n)$ with all partial second-order derivatives. The Hessian of $f$ is a matrix-valued function $\mathcal{H}: \mathbb{R}^n \to \mathcal{M}_n(\mathbb{R})$ defined as follows:

$$ \mathcal{H}(\mathbf{x}) = \begin{pmatrix} \frac{\partial^2 f}{\partial x_1^2}(\mathbf{x}) & \dots & \frac{\partial^2 f}{\partial x_1 \partial x_n}(\mathbf{x}) \\ \vdots & \ddots & \vdots \\ \frac{\partial^2 f}{\partial x_n \partial x_1}(\mathbf{x}) & \dots & \frac{\partial^2 f}{\partial x_n^2}(\mathbf{x}) \end{pmatrix}. \qquad (36) $$

Evaluating the Hessian at some point $\mathbf{x} \in \mathbb{R}^n$ provides a description of the local convexity of $f$ at that point. If $\mathbf{x}^*$ is a critical point ($\nabla f(\mathbf{x}^*) = 0$) then $f$ can be locally approximated by a quadratic function

$$ f(\mathbf{x}^* + \mathbf{u}) = f(\mathbf{x}^*) + \frac{1}{2} \mathbf{u}^T \mathcal{H}(\mathbf{x}^*) \mathbf{u} + O(\mathbf{u}^3). \qquad (37) $$

The Hessian matrix at critical point $\mathbf{x}^*$ can be decomposed into its eigenspace in order to obtain the principal directions of curvature. The sign of the corresponding eigenvalues determine whether each direction is stable or unstable.
For the problem at hand, we wish to study the critical points of the secular energy $\bar{U}_B$ as a function of spin direction. For the following, we define the relative alignment between spins $\hat{s}_1$ and $\hat{s}_2$:

$$ \sigma = \sigma_1 \sigma_2. \qquad (38) $$

### Case 1. Equilibrium configurations in the orbital plane

As we have seen, any pair of spin axes which are parallel to each other and contained within the orbital plane will correspond to a secular equilibrium state of the averaged system. Indeed, for any given $\hat{\mathbf{p}} \cdot \hat{\mathbf{e}}_z = 0$ in Eq. (32) the energy takes the same value $\bar{U}_B(\sigma_1 \hat{\mathbf{p}}, \sigma_2 \hat{\mathbf{p}}) = \sigma$, depending only on the relative orientation of the two spins. The sign of $\sigma = \sigma_1 \sigma_2$ will be positive if the spins are facing the same direction or negative if they are in opposite directions. This invariance with respect to a choice of $\hat{\mathbf{p}}$ in fact reflects the more general axis-symmetry of the system with respect to rotations around the axis $\hat{\mathbf{e}}_z$. Simultaneous rotations of both spins will leave the energy expression invariant. In the case of this particular equilibrium configuration, the act of choosing one $\hat{\mathbf{p}}$ lying in the orbital plane over another $\hat{\mathbf{p}}'$ corresponds to rotating both spins simultaneously by the angle that takes $\hat{\mathbf{p}}' \mapsto \hat{\mathbf{p}}$. The set of all unit vectors $\hat{\mathbf{p}}$ in the orbital plane corresponds to a one-dimensional level-set of constant energy (a circle). To study the stability of this level-set we can reduce the degrees of freedom of the system from four (two for each spin axis, i.e. $S_2 \times S_2$) to three (by removing the rotational degree of freedom). In this reduced three-dimensional manifold, the circle becomes a point, and the remaining three directions of the tangent space determine the convexity of the energy.

As a first step, we consider the spins parametrised in spherical coordinates,

$$ \hat{s}_{\ell} = (\cos \psi_{\ell} \sin \theta_{\ell}, \sin \psi_{\ell} \sin \theta_{\ell}, \cos \theta_{\ell}), \qquad (39) $$

with azimuth $\psi_{\ell}$ and angle $\theta_{\ell}$ measured from the north pole. In these coordinates, the secular energy takes the form

$$ \bar{U}_B = 2 \cos \theta_1 \cos \theta_2 - \cos(\psi_1 - \psi_2) \sin \theta_1 \sin \theta_2. \qquad (40) $$

As discussed above, we can see clearly that simultaneous rotations in the azimuth $\psi_1$ and $\psi_2$ do not affect the energy. We can therefore look at the difference $\psi_1 - \psi_2 := \Delta\psi$ and discard the redundant degree of freedom $\psi_1 + \psi_2$. At the equilibrium, the polar angles are $\theta_1 = \theta_2 = \pi/2$, and the azimuth is either $\Delta\psi = 0$ or $\Delta\psi = \pi$. We consider therefore the three directions of the tangent space of this reduced manifold. In this case, the Hessian matrix of the energy $\bar{U}_B$ with respect to the three variables $(\Delta\psi, \theta_1, \theta_2)$ has value:

$$ \mathcal{H}_1 = \begin{pmatrix} \sigma & 0 & 0 \\ 0 & \sigma & 2 \\ 0 & 2 & \sigma \end{pmatrix}. \qquad (41) $$

The sign of $\sigma = \pm 1$ corresponds to the two cases of the spin directions being aligned or anti-aligned. For either sign, the Hessian has both positive and negative eigenvalues $(-1, 1$ and $3\sigma)$. This implies a saddle point, which is unstable.

### Case 2. Equilibrium configurations orthogonal to the orbital plane

In the second case, we consider both spins to be orthogonal to the plane of the orbit. The spin axes $\hat{s}_{\ell}$ therefore will be lying

in the north or south pole of their corresponding unit spheres. A natural choice of coordinates for each spin in the neighbourhood of the poles is the Cartesian pair $x_{\ell}, y_{\ell}$, which with the unit-norm condition gives

$$ \hat{s}_{\ell} = (x_{\ell}, y_{\ell}, \sigma_{\ell} \sqrt{1 - x_{\ell}^2 - y_{\ell}^2}). \qquad (42) $$

In these coordinates, the secular magnetic energy can be expressed in the form:

$$ \bar{U}_B = -x_1 x_2 - y_1 y_2 + 2\sigma \sqrt{1 - x_1^2 - y_1^2} \sqrt{1 - x_2^2 - y_2^2}. \qquad (43) $$

The corresponding basis for the tangent space is the coordinate basis $dx_{\ell}$ and $dy_{\ell}$. With respect to these coordinates the Hessian of the energy takes the values:

$$ \mathcal{H}_2 = - \begin{pmatrix} 2\sigma & 0 & 1 & 0 \\ 0 & 2\sigma & 0 & 1 \\ 1 & 0 & 2\sigma & 0 \\ 0 & 1 & 0 & 2\sigma \end{pmatrix}. \qquad (44) $$

The nature of the eigenvalues of the Hessian change according to the sign of $\sigma$: when the spins are aligned ($\sigma = +1$), the Hessian is negative definite, with eigenvalues $\{-3, -3, -1, -1\}$, implying a point of maximum energy and therefore instability; when the spins are in opposing directions ($\sigma = -1$), the eigenvalues are $\{1, 1, 3, 3\}$ and the Hessian is positive definite, which implies a point of minimum energy and hence stability.

This suggests that given enough time and in the absence of stronger torques, magnetically interacting systems will naturally converge towards the stable anti-aligned orthogonal configuration. The timescale of this convergence will depend on the strength of the dissipation effects involved.

### 4.3. Solutions near equilibrium

In Appendix B, we determine solutions near the stable equilibrium state. These solutions can be plugged in to derive complete secular orbital dynamics on systems which present strong magnetic torques, such as WD or NS binaries in the context of gravitational wave emission (for further details, see e.g. Bourgoin et al. 2022; Lira et al. 2022). We present below the main results:

The spins can be decomposed as $\hat{s}_{\ell} = x_{\ell} \hat{e}_x + y_{\ell} \hat{e}_y + z_{\ell} \hat{e}_z$, where the orbital-plane components satisfy

$$ x_{\ell}(t) = \rho_{\ell}^{+} \cos(\omega_{p}^{+} t + \phi_{\ell}^{+}) + \rho_{\ell}^{-} \cos(\omega_{p}^{-} t + \phi_{\ell}^{-}), \qquad (45a) $$

$$ y_{\ell}(t) = \rho_{\ell}^{+} \sin(\omega_{p}^{+} t + \phi_{\ell}^{+}) + \rho_{\ell}^{-} \sin(\omega_{p}^{-} t + \phi_{\ell}^{-}), \qquad (45b) $$

and the orthogonal component satisfies

$$ z_1(t) = \frac{1}{2\nu_2} (\zeta + \bar{Q}_z + \rho_z \cos(\omega_z t + \phi_z)), \qquad (46a) $$

$$ z_2(t) = \frac{1}{2\nu_1} (\zeta - \bar{Q}_z - \rho_z \cos(\omega_z t + \phi_z)), \qquad (46b) $$

with real parameters $\rho_{\ell}^{\pm}, \phi_{\ell}^{\pm}, \rho_z, \phi_z$, dependent on initial conditions. The expressions of the constants $\zeta$ and $\bar{Q}_z$ and the frequencies $\omega_z, \omega_p^{+}$, and $\omega_p^{-}$ are given in the appendix.
<figure><img src="image_5.png" /><figcaption>**Fig. 4.** Evolution of the binary system in a secular timescale with dissipation introduced. Each column corresponds to a given initial condition from Table 2, chosen in a neighbourhood of an equilibrium point. On the top row: the trajectory of each spin axis $\hat{s}_l$ is plotted against the unit sphere, in blue for the primary and green for the secondary. Colours are interpolated towards red from initial time to final time of convergence. Initial conditions for the spin axes are portrayed as dashed lines from the origin to $\hat{s}_l(0)$. The orbital plane is represented in a darker shade with black contours, and the axis $\hat{e}_z$ as a vertical black arrow. On the middle row: the secular magnetic interaction energy. On the bottom row: the angular distance of each spin axis with respect to the stable equilibrium point $(\hat{e}_z, -\hat{e}_z)$.</figcaption></figure>

**5. Numerical validation**

We present in this section a numerical verification of the results that were obtained in Sect. 4. In the first part, we demonstrate the derived (in-)stability of each equilibrium configuration. In the second part, we compare the obtained analytical solutions to numerical integration, in the neighbourhood of the stable equilibrium.

**5.1. Stability of the equilibrium configurations**

We artificially introduce dissipation into the dynamics of the secular system and numerically show that it is driven towards the stable states. Such a dissipative effect must preserve the norm condition on unit vectors and manifest as a friction when orientation changes. For this we include a time-delay term into the magnetic field that is felt by each companion star. By our previous arguments in Sect. 3, it is direct to see that this term will propagate to the secular scale as follows:

$$ \left\langle \frac{d\hat{s}_1}{d\tau} \right\rangle (\tau) = \hat{s}_1(\tau) \times \tilde{B}(\hat{s}_2(\tau - \Delta\tau)), \qquad (47a) $$

$$ \left\langle \frac{d\hat{s}_2}{d\tau} \right\rangle (\tau) = \kappa \hat{s}_2(\tau) \times \tilde{B}(\hat{s}_1(\tau - \Delta\tau)), \qquad (47b) $$

where we have used the re-scaled dimensionless time $\tau \mapsto \nu_1^{-1}\tau$ and $\kappa = \nu_2/\nu_1 = s_1/s_2$ presented at the end of Sect. 3. For convenience, we consider the delay to be an order of magnitude smaller than the average torque timescale, $\Delta\tau \sim 0.1 \kappa^{-1/2}$. This provides us with dissipative effects visible on the simulation timescales.

<table><caption>**Table 2.** Initial conditions in four distinct stability simulations (a–d).</caption><thead><tr><th></th><th>θ<sub>1</sub> (°)</th><th>θ<sub>2</sub> (°)</th><th>φ<sub>1</sub> (°)</th><th>φ<sub>2</sub> (°)</th></tr></thead><tbody><tr><td>(a)</td><td>10</td><td>175</td><td>0</td><td>0</td></tr><tr><td>(b)</td><td>0</td><td>1</td><td>0</td><td>0</td></tr><tr><td>(c)</td><td>90</td><td>89</td><td>0</td><td>181</td></tr><tr><td>(d)</td><td>90</td><td>91</td><td>0</td><td>1</td></tr></tbody></table>

**Notes.** Each spin axis is parametrised in polar coordinates with azimuth $\psi_l$ and angle $\theta_l$ measured from the north pole.

In order to assess stability, we consider two metrics: the secular magnetic interaction energy $\bar{U}_B$ (see Eq. (35)); and the angular distance to the stable equilibrium point $(\hat{e}_z, -\hat{e}_z)$, which we define as:

$$ d_1 = \arg(\hat{s}_1, +\hat{e}_z), \qquad (48a) $$

$$ d_2 = \arg(\hat{s}_2, -\hat{e}_z), \qquad (48b) $$

where $\arg(\mathbf{u}, \mathbf{v}) = \arctan2(|\mathbf{u} \times \mathbf{v}|, \mathbf{u} \cdot \mathbf{v})$ is the angle between any two vectors $\mathbf{u}$ and $\mathbf{v}$.

The dimensionless system is then integrated until convergence, for four sets of initial spin conditions (Table 2), and with spin ratio $\kappa = 0.3$. Each initial condition corresponds to an unstable equilibrium state plus a small perturbation of the order of $\sim 1^\circ$. For completeness, we also include a perturbation of the stable state ($\sim 10^\circ$). The resulting trajectories are plotted in Fig. 4, together with the time evolution of the secular energy $\bar{U}_B$ and of
<figure><img src="image_6.png" /><figcaption>Fig. 5. Comparison between the analytical precession model (blue) and numerical integration (red), close to stable equilibrium. On the leftmost column: the components of the spin of the primary, in the time domain. The corresponding Discrete Fourier Transform $\mathcal{F}\{\hat{s}_1\}$ is given in absolute value (power spectrum in decibels, middle column), and in complex argument (rightmost column). The time and frequency axes are expressed in the re-scaled dimensionless units.</figcaption></figure>

the angular distances $d_1$ and $d_2$. Even for a small initial angular perturbation, the spin axes diverge from their original unstable equilibrium and converge towards the stable, anti-aligned orthogonal configuration.

## 5.2. Solutions near equilibrium

In this part, we present a brief application of the analytical solutions of the secular equations that were exposed in Sect. 4.3. The equations are expressed in dimensionless time, and we adopt a spin ratio of $\kappa = 0.3$. The following initial conditions are considered (polar coordinates): inclinations $\theta_1 = 10^\circ$ and $\theta_2 = 172.5^\circ$ from the north pole, and azimuths $\psi_1 = 0$ and $\psi_2 = 50^\circ$. Figure 5 compares the obtained analytical expression of the primary to numerical integration, decomposed in the Cartesian basis $\hat{s}_1(\tau) = x_1(\tau)\hat{e}_x + y_1(\tau)\hat{e}_y + z_1(\tau)\hat{e}_z$. Matching peaks can be observed in the power spectra of both solutions, at the obtained angular frequencies (in dimensionless units) $\omega_p^- = 0.47$ rad and $\omega_p^+ = 1.86$ rad for the orbital components $(x_1, y_1)$, and at $\omega_z = 2.34$ rad for the orthogonal component $z_1$ (cf. Appendix B).

# 6. Discussion

We have so far defined the concept of a secular equilibrium state and applied it to obtain the equilibrium dynamics of binary systems with two magnetic components. In this section, we begin by discussing the fine points between secular and instantaneous

equilibrium. Then, we apply our results to a real astrophysical scenario, the $\epsilon$-Lupi magnetic binary.

## 6.1. Comparison between instantaneous and secular equilibrium states

In Sect. 4, we defined the secular equilibrium as the spin configurations where the net torque over an orbital period is effectively zero. These states contrast with the instantaneous equilibrium, the configurations of zero torque on the instantaneous precession system (Eq. (13)). In the latter scenario, the spin dynamics are considered at a single moment in time and at a fixed orbital position. When orbital motion is introduced, an instantaneous equilibrium state may be destabilised. This can be seen in the following manner. We consider a bounded orbit parametrised by the osculating true anomaly $f = f(t)$. Expressing the spins in spherical coordinates, the instantaneous magnetic energy $U_B$ (Eq. (34)) takes the form:

$$U_B(t) = \frac{\mu_0\mu_1\mu_2}{4\pi r^3} (\hat{\mu}_1 \cdot \hat{\mu}_2 - 3 \sin\theta_1 \sin\theta_2 \cos(\psi_1 - f(t)) \cos(\psi_2 - f(t))). \quad (49)$$

Whereas the dot-product term on the right-hand side is invariant to an orbital translation of the bodies, the cosine terms oscillate with the orbital dynamics. We take two instants within the same orbital revolution, $t_1$ and $t_2$ such that the true anomaly at each
<table><caption>**Table 3.** Comparison between the states of instantaneous and of secular equilibrium.</caption><thead><tr><th colspan="5">Instantaneous equilibria</th><th colspan="5">Secular equilibria</th></tr><tr><th>Config.</th><th>Cond.</th><th>σ</th><th>U<sub>B</sub></th><th>Stable</th><th>Config.</th><th>Cond.</th><th>σ</th><th>U&#773;<sub>B</sub></th><th>Stable</th></tr></thead><tbody><tr><td>⇒</td><td>ŝ<sub>ℓ</sub> || ř</td><td>+1</td><td>-2</td><td>yes</td><td>↑↑</td><td>ŝ<sub>ℓ</sub> || ŵ<sub>z</sub></td><td>+1</td><td>+2</td><td>no</td></tr><tr><td>↑↓</td><td>ŝ<sub>ℓ</sub> ⊥ ř</td><td>-1</td><td>-1</td><td>no</td><td>⇔</td><td>ŝ<sub>ℓ</sub> ⊥ ŵ<sub>z</sub></td><td>-1</td><td>+1</td><td>no</td></tr><tr><td>↑↑</td><td>ŝ<sub>ℓ</sub> ⊥ ř</td><td>+1</td><td>+1</td><td>no</td><td>⇒</td><td>ŝ<sub>ℓ</sub> ⊥ ŵ<sub>z</sub></td><td>+1</td><td>-1</td><td>no</td></tr><tr><td>⇔</td><td>ŝ<sub>ℓ</sub> || ř</td><td>-1</td><td>+2</td><td>no</td><td>↑↓</td><td>ŝ<sub>ℓ</sub> || ŵ<sub>z</sub></td><td>-1</td><td>-2</td><td>yes</td></tr></tbody></table>

**Notes.** Left side: the instantaneous equilibria require spins either parallel or perpendicular to the orbital separation $\hat{r}$. Right side: the secular equilibria require spins either parallel or perpendicular to the axis $\hat{e}_z$. The energies $U_B$ and $\bar{U}_B$ are given in dimensionless units.

instant equals $f(t_1) \equiv (\psi_1 + \psi_2)/2$ and $f(t_2) \equiv (\psi_1 + \psi_2 + \pi)/2$. Then the difference in energy between those two instants is roughly

$$U_B(t_2) - U_B(t_1) \sim \frac{3\mu_0\mu_1\mu_2}{4\pi a^3} \sin \theta_1 \sin \theta_2, \quad (50)$$

where we have substituted $r \sim a$. We conclude that instantaneous equilibrium positions will indeed develop large energy oscillations in the timescale $t \sim P_{\text{orb}}$, particularly when the polar angle $\theta_\ell$ is large (i.e. spin axes close to orbital plane alignment). For rapidly orbiting systems, this energy fluctuation occurs very quickly and destabilises the equilibrium of the point.

There is a direct analogy between the states of secular equilibrium and those of instantaneous equilibrium. We recall the expression of the secular magnetic energy:

$$\bar{U}_B = -\hat{s}_1 \cdot \bar{\mathcal{B}}(\hat{s}_2), \quad \bar{\mathcal{B}} = \mathbb{I} - 3\hat{e}_z \otimes \hat{e}_z.$$

Such expression has been normalised by the positive constant $\mu_0\mu_1\mu_2/8\pi b^3$ as discussed in Sect. 4. We similarly normalise the expression of the instantaneous magnetic energy (cf. Eq. (34)) by the scalar $\mu_0\mu_1\mu_2/4\pi r^3$ and obtain:

$$U_B \propto + \hat{s}_1 \cdot \tilde{\mathcal{B}}_r(\hat{s}_2), \quad \tilde{\mathcal{B}}_r = \mathbb{I} - 3 \hat{r} \otimes \hat{r}.$$

We observe that the secular averaging procedure effectively produced a flip in the sign of the energy as well as a switch $\hat{r} \leftrightarrow \hat{e}_z$. Analogously to how the secular equilibrium states of the binary are given by the singular vectors of $\tilde{\mathcal{B}}$ (cf. Sect. 4), the instantaneous states are given by the singular vectors of $\tilde{\mathcal{B}}_r$. Consequently, each state of instantaneous equilibrium has a secular counterpart. These correspond to spin axes aligned with the radial direction $\hat{r}$ (resp. $\hat{e}_z$) or perpendicular to $\hat{r}$ (resp. $\hat{e}_z$). The stability of each state depends on the local convexity of the energy $U_B$ (resp. $\bar{U}_B$). We present in Table 3 a comparison between the obtained states for the two types of equilibrium.

## 6.2. $\epsilon$-Lupi

We consider as a potential application case the $\epsilon$-Lupi inner binary system. $\epsilon$-Lupi is a ternary system composed of two close-range B-type companions Aa and Ab, plus a third distant companion dubbed $\epsilon$-Lupi B. Both stars of the $\epsilon$-Lupi A inner system are magnetic, making it the first and only currently known massive binary that has two magnetic components. Furthermore, the field of each star can be captured by a dipolar model with axes roughly parallel to the spins (Shultz et al. 2015). The system also has a short orbital period of $P_{\text{orb}} \sim 4.56$ d, making it an excellent example to apply our model. Pablo et al. (2019),

Uytterhoeven et al. (2005) obtained estimates for relevant stellar and orbital parameters, from which we adopt the values for the semi-major axis $a = 29.5 R_\odot$, eccentricity of $e = 0.28$, inclination $i = 21^\circ$, stellar masses of $m_1 = 9.0 M_\odot$ for the primary and $m_2 = 7.9 M_\odot$ for the secondary, and radii $R_1 = R_2 = 4.5 R_\odot$.

Shultz et al. (2015) reports a dipolar field strength of at least $B_1^p = 600$ G and $B_2^p = 900$ G at the surface poles of the star, and projected rotational velocities at the equator $v_1 \sin i_1 = 37$ km s$^{-1}$ for the primary and $v_2 \sin i_2 = 27$ km s$^{-1}$ for the secondary. By adopting the rotational inclination for each star equal to the orbital plane inclination $i \sim 21^\circ$, we obtain rotational periods of the order of $P_1 = 2.2$ d and $P_2 = 3.0$ d.

From these physical parameters we may calculate the corresponding dimensionless ratios. The impact of figure effects can be assessed via $\gamma_\ell^{\text{fig}} = 2 \times 10^{10} J_2^\ell$, whereas for tides $\gamma_1^{\text{tide}} = 2.5 \times 10^8 (k_2^1/Q_1)$ and $\gamma_2^{\text{tide}} = 3.5 \times 10^8 (k_2^2/Q_2)$. These expressions suggest that if $\epsilon$-Lupi is both highly symmetrical with $J_2^\ell \lesssim 6 \times 10^{-11}$, and has tidal parameters $k_2^\ell/Q_\ell \lesssim 3 \times 10^{-9}$, then the system's rotation is likely driven by magnetism. In this case, the smallness of $\eta = 5 \times 10^{-12}$ and $\epsilon = 4 \times 10^{-11}$ predict that the system will be driven towards the secular stable equilibrium in a timescale proportional to energy dissipation rates. This outcome corresponds well to the expected state of $\epsilon$-Lupi based on observational data (Pablo et al. 2019).

## 7. Conclusions

This paper has presented an analysis of the secular precession dynamics of binary systems under pure magnetic dipole-dipole interactions, considering an effective description with orbit-averaged motion. In particular, we have supposed spin-dipole alignment, perfect sphericity and tidal rigidity, and then derived criteria for assessing the validity of these assumptions, as well as the relative strengths of magnetic dipole interactions, tidal torques, and figure effects. We have shown that this effective long-term description predicts a set of states of secular equilibrium which confront the traditional states of instantaneous magnetic equilibrium, where orbital dynamics are in fact neglected. Indeed, we have determined that there is a single secular state that is globally stable, corresponding to the configuration $\pm(\hat{e}_z, -\hat{e}_z)$ where the spin axes are reversed with respect to each other and orthogonal to the orbital plane. Conversely, the instantaneous state of radial alignment $\pm(\hat{r}, \hat{r})$ is in fact only momentarily stable, since the orbital motion generates energy fluctuations and destabilises the configuration. Our work can also be used to derive the long-term evolution of binary orbits providing an expected spin evolution in the absence of strong additional torques.
Our results hold for typical early-type MMS (such as the observed ϵ-Lupi system), WD, and NS binaries hosting dipolar fossil-like fields, where we expect long-term convergence towards the secularly stable state. Another interesting case of application is that of M dwarfs. For masses lower than $0.35 M_⊙$, M dwarfs are fully convective, unlike any other main-sequence stars, which renders the dynamo-driven topology unique (e.g. Dobler et al. 2006; Browning 2008). These low- and mid-mass M dwarfs are possible targets of application for our formalism since their magnetic fields often display intense dipolar components (e.g. Donati et al. 2008), although intermittent higher-order multipolar components might also be present (e.g. Kochukhov 2021, and references therein). Moreover, they have been observed forming binaries with two magnetic components (e.g. Kochukhov & Lavail 2017; Kochukhov & Shulyak 2019). An interesting example of such binary systems is that of the YY Gem system (Kochukhov & Shulyak 2019), where each star has both dipolar and multipolar components. The Zeeman-Doppler imaging analysis indeed revealed moderately complex global fields with a typical strength of 200-300 G, with dipolar components that are anti-aligned as predicted for fossil-type fields. These considerations hint that our work might be applied to any type of dipolar magnetic fields, either of fossil origin or triggered by a dynamo action, as long as its variation timescales are longer than the precession timescales. However, as mentioned in Kochukhov & Shulyak (2019), the Zeeman intensification analysis suggests that the global fields of YY Gem may only comprise a few percent of their total magnetic fields, which highlights the need to consider multipoles in future studies. In addition, we point out that in typical M dwarf binary systems, gravitational interactions may control the orientation of the spins themselves. We can in fact compute (for a typical magnetic M dwarf: $m \sim 0.35 M_⊙$, $B \sim 1$ kG) the dimensionless parameters that measure the strength of extended-body gravitational interactions with respect to magnetic torques: $\gamma_{\text{fig}} = 6 \times 10^{11} J_2$ and $\gamma_{\text{tides}} = 7 \times 10^7 k_2/Q$ (Eqs. (9), (11)). Subsequently, we obtain the (quite tight) cutoffs $J_2 \le 10^{-12}$ for the dimensionless quadrupole moment (Eq. (10)), and $k_2/Q \le 10^{-8}$ for the tidal parameters (Eq. (12)), as the necessary parameters for magnetism to control the spin evolution.

Another key assumption from this work that warrants discussion is the alignment between the spin and magnetic axes of each star, since many misaligned systems are observed in nature (e.g Shultz et al. 2019). In the case where the spin and magnetic dipole are anti-aligned ($\mu_\ell < 0$ in Eq. (5)), our general results still hold but spin directions must be flipped accordingly in the equations. More generally, this alignment constraint must be relaxed in future studies to explore the more general scenario of misaligned spin and magnetic axes. Assuming a rigid description where the magnetic axes rotate around the spins, a potentially direct case could be when the stellar rotation period $P_\ell$ is much shorter than the orbital period $P_{\text{orb}}$. In such regime, the problem may be hierarchically split into three distinct timescales $P_\ell \ll P_{\text{orb}} \ll \tau$, where $\tau$ is the ‘spin precession timescale’ as described in Sect. 2 (Eq. (15)). The dynamics could then be formulated as an effective description when seen from the longer orbital timescale $P_{\text{orb}}$, potentially reducing to the one explored in this work. Accordingly, one could expect our main results to still hold.

Further extensions of this work will also include abandoning the magnetostatic description and considering internal coupling of the fields with matter (Campbell 2018). Finally, the precession dynamics and equilibrium states of the system may be

investigated by directly taking into account not only magnetic forces but also competing figure effects and dynamical tides (see example of such combined study in Ahuir et al. 2021).

**Acknowledgements.** C.A. acknowledges the joint finantial support of Centre National d'Études Spatiales (CNES) and École Doctorale Astronomie et Astrophysique d'Ile-de-France (ED127 AAIF). This work was also supported by the Programme National GRAM, by PNPS (CNRS/INSU), by INP and IN2P3 co-funded by CNES, and by CNES LISA grants at CEA/IRFU. Authors are also grateful to S. Bouquillon for fruitful discussions. We are thankful to the anonymous reviewer for their constructive comments, which allowed us to improve our article.

## References

Achilleos, N., & Wickramasinghe, D. T. 1989, ApJ, 346, 444

Ahuir, J., Strugarek, A., Brun, A. S., & Mathis, S. 2021, A&A, 650, A126

Akgün, T., Reisenegger, A., Mastrano, A., & Marchant, P. 2013, MNRAS, 433, 2445

Alecian, E., Neiner, C., Wade, G. A., et al. 2014, Proc. Int. Astron. Union, 9, 330

Amaro-Seoane, P., Audley, H., Babak, S., et al. 2017, Laser Interferometer Space Antenna (USA: NASA)

Aurière, M., Wade, G. A., Silvester, J., et al. 2007, A&A, 475, 1053

Bagnulo, S., & Landstreet, J. D. 2021, MNRAS, 507, 5902

Beuermann, K., Euchner, F., Reinsch, K., Jordan, S., & Gänsicke, B. T. 2007, A&A, 463, 647

Borra, E. F., Landstreet, J. D., & Mestel, L. 1982, ARA&A, 20, 191

Bourgoin, A., Le Poncin-Lafitte, C., Mathis, S., & Angonin, M.-C. 2022, Phys. Rev. D, 105, 124042

Braithwaite, J. 2008, MNRAS, 386, 1947

Braithwaite, J. 2009, MNRAS, 397, 763

Braithwaite, J., & Cantiello, M. 2013, MNRAS, 428, 2789

Braithwaite, J., & Nordlund, Å. 2006, A&A, 450, 1077

Braithwaite, J., & Spruit, H. C. 2004, Nature, 431, 819

Bromley, B. C., & Kenyon, S. J. 2022, AJ, 164, 229

Browning, M. K. 2008, ApJ, 676, 1262

Brun, A. S., & Browning, M. K. 2017, Liv. Rev. Sol. Phys., 14, 4

Brun, A., Strugarek, A., Noraz, Q., et al. 2022, ApJ, 926, 21

Campbell, C. G. 2018, Magnetohydrodynamics in Binary Stars (Berlin: Springer), 456

Carrier, F., North, P., Udry, S., & Babel, J. 2002, A&A, 394, 151

Carvalho, G. A., Anjos, R. C. d., Coelho, J. G., et al. 2022, ApJ, 940, 90

Charbonneau, P., & MacGregor, K. B. 2001, ApJ, 559, 1094

Commerçon, B., Hennebelle, P., Audit, E., Chabrier, G., & Teyssier, R. 2010, A&A, 510, L3

Cumming, A. 2002, MNRAS, 333, 589

Damiani, C., & Lanza, A. F. 2015, A&A, 574, A39

de Mink, S. E., Sana, H., Langer, N., Izzard, R. G., & Schneider, F. R. N. 2014, ApJ, 782, 7

Dobler, W., Stix, M., & Brandenburg, A. 2006, ApJ, 638, 336

Donati, J. F., & Landstreet, J. D. 2009, ARA&A, 47, 333

Donati, J. F., Morin, J., Petit, P., et al. 2008, MNRAS, 390, 545

Duez, V. 2011, Astron. Nachr., 332, 983

Duez, V., Braithwaite, J., & Mathis, S. 2010, ApJ, 724, L34

Ferrario, L., & Wickramasinghe, D. T. 2005, MNRAS, 356, 615

Ferrario, L., Pringle, J., Tout, C., & Wickramasinghe, D. 2009, MNRAS, 400, L71

Ferrario, L., de Martino, D., & Gänsicke, B. T. 2015a, Space Sci. Rev., 191, 111

Ferrario, L., Melatos, A., & Zrake, J. 2015b, Space Sci. Rev., 191, 77

Ferrario, L., Wickramasinghe, D., & Kawka, A. 2020, Adv. Space Res., 66, 1025

Fuller, J., & Mathis, S. 2023, MNRAS, 520, 5573

Gaurat, M., Jouve, L., Lignières, F., & Gastine, T. 2015, A&A, 580, A103

Gerosa, D., Kesden, M., Sperhake, U., Berti, E., & O'Shaughnessy, R. 2015, Phys. Rev. D, 92, 064016

Grunhut, J., Wade, G.A., & the MiMeS Collaboration 2013, EAS Pub. Ser., 64, 67

Hut, P. 1980, A&A, 92, 167

Hut, P. 1981, A&A, 99, 126

Isern, J., García-Berro, E., Külebi, B., & Lorén-Aguilar, P. 2017, ApJ, 836, L28

Jermyn, A. S., & Cantiello, M. 2021, ApJ, 923, 104

Jouve, L., Lignières, F., & Gaurat, M. 2020, A&A, 641, A13

Kaspi, V. M., & Beloborodov, A. M. 2017, ARA&A, 55, 261

Kawka, A., & Vennes, S. 2004, IAU Symp., 224, 879

King, A., Whitehurst, R., & Frank, J. 1990, MNRAS, 244, 731

Kochukhov, O. 2021, A&ARv, 29, 1
Kochukhov, O., & Lavail, A. 2017, ApJ, 835, L4

Kochukhov, O., & Shulyak, D. 2019, ApJ, 873, 69

Kochukhov, O., Lundin, A., Romanyuk, I., & Kudryavtsev, D. 2011, ApJ, 726, 24

Landstreet, J. D., & Bagnulo, S. 2020, A&A, 634, L10

Landstreet, J. D., & Mathys, G. 2000, A&A, 359, 213

Landstreet, J. D., Bagnulo, S., Valyavin, G., & Valeev, A. F. 2017, A&A, 607, A92

Liebert, J. 1988, PASP, 100, 1302

Liebert, J., Wickramasinghe, D. T., Schmidt, G. D., et al. 2005, AJ, 129, 2376

Lira, M., Degollado, J. C., Moreno, C., & Núñez, D. 2022, Gen. Rel. Grav., 54, 146

Maggiore, M., Broeck, C. V. D., Bartolo, N., et al. 2020, J. Cosmol. Astropart. Phys., 2020, 050

Mathis, S., Duez, V., & Braithwaite, J. 2011, IAU Symp., 271, 270

Morel, T., Castro, N., Fossati, L., et al. 2014, Proc. Int. Astron. Union, 9, 342

Moss, D. 1987, MNRAS, 226, 297

Moss, D. 2003, A&A, 403, 693

Pablo, H., Shultz, M., Fuller, J., et al. 2019, MNRAS, 488, 64

Phinney, E. S., & Kulkarni, S. R. 1994, ARA&A, 32, 591

Poisson, E., & Will, C. M. 2014, Gravity: Newtonian, Post-Newtonian, Relativistic (Cambridge: Cambridge University Press)

Pound, A., & Poisson, E. 2008, Phys. Rev. D, 77, 043524

Schneider, F. R. N., Podsiadlowski, P., Langer, N., Castro, N., & Fossati, L. 2016, MNRAS, 457, 2355

Schneider, F., Ohlmann, S., Podsiadlowski, P., et al. 2019, Nature, 574, 211

Schreiber, M. R., Belloni, D., Gänsicke, B. T., Parsons, S. G., & Zorotovic, M. 2021, Nat. Astron., 5, 648

Shultz, M., Wade, G. A., Alecian, E., & Collaboration, B. 2015, MNRAS, 454, L1

Shultz, M. E., Wade, G. A., Rivinius, T., et al. 2019, MNRAS, 490, 274

Stello, D., Cantiello, M., Fuller, J., et al. 2016, Nature, 529, 364

Stibbs, D. W. N. 1950, MNRAS, 110, 395

Strugarek, A., Bolmont, E., Mathis, S., et al. 2017, ApJ, 847, L16

Tayler, R. J. 1980, MNRAS, 191, 151

Tout, C. A., Wickramasinghe, D. T., & Ferrario, L. 2004, MNRAS, 355, L13

Tout, C. A., Wickramasinghe, D. T., Liebert, J., Ferrario, L., & Pringle, J. E. 2008, MNRAS, 387, 897

Uytterhoeven, K., Harmanec, P., Telting, J. H., & Aerts, C. 2005, A&A, 440, 249

Vidal, J., Cébron, D., ud-Doula, A., & Alecian, E. 2019, A&A, 629, A142

Wade, G. A., Neiner, C., Alecian, E., et al. 2015, MNRAS, 456, 2

Wickramasinghe, D. T., & Ferrario, L. 2000, PASP, 112, 873

Wickramasinghe, D. T., & Ferrario, L. 2005, MNRAS, 356, 1576

Will, C. M., & Maitra, M. 2017, Phys. Rev. D, 95, 064003
## Appendix A: Optimisation of bilinear forms

In this appendix we recall the variational characterisation of the singular value decomposition. Consider the bilinear form $U$ taking two unit vectors $U : S_n \times S_n \to \mathbb{R}$, represented under a basis $\{\hat{e}_1, ..., \hat{e}_n\}$ via some matrix $U_{ij}$:

$$U(\hat{u}, \hat{v}) = U_{ij} u^i v^j, \qquad (A.1)$$

where $\hat{u} = u^i \hat{e}_i$, $\hat{v} = v^i \hat{e}_i$, and sum over repeated indices is presupposed.

Since $U$ is continuous on the compact domain $S_n \times S_n$, it attains both a maximum and minimum value. Indeed, consider the Lagrange function

$$g(\hat{u}, \hat{v}) = U(\hat{u}, \hat{v}) - \lambda_1(\hat{u} \cdot \hat{u} - 1) - \lambda_2(\hat{v} \cdot \hat{v} - 1), \qquad (A.2)$$

with multipliers $\lambda_1$ and $\lambda_2$. The extrema $\hat{u}_*$ and $\hat{v}_*$ necessarily satisfy the stationary condition $\nabla g(\hat{u}_*, \hat{v}_*) = 0$, which implies that for any tangent vectors $h_1, h_2$,

$$U(h_1, \hat{v}_*) = 2\lambda_1 \hat{u}_* \cdot h_1, \qquad (A.3a)$$

$$U(\hat{u}_*, h_2) = 2\lambda_2 \hat{v}_* \cdot h_2. \qquad (A.3b)$$

From the matrix representation of $U$ [Eq. (A.1)], we see that the solutions $\hat{u}_*$ and $\hat{v}_*$ correspond respectively to the left and right singular vectors of $U_{ij}$ — that is, the unique set of vectors that satisfy

$$U_{ij} u_i^* = \lambda v_i^*, \qquad U_{ij} v_j^* = \lambda u_j^*, \qquad (A.4a)$$

where $\lambda = 2\lambda_1 = 2\lambda_2$ is the corresponding singular value of $U_{ij}$.

It is direct to see that the attained value $U_* = U(\hat{u}_*, \hat{v}_*) = \lambda$. From all singular-vector pairs $(\hat{u}_*, \hat{v}_*)$, the global maximum (or minimum) of $U$ is therefore reached by the candidate pair that has the highest (or lowest) corresponding singular value $\lambda$. We note that whenever $U$ is symmetric, the singular vectors and eigenvectors of $U_{ij}$ coincide, and the singular values are equal to the magnitude of the eigenvalues.

## Appendix B: Linearised solutions

In this section, we formally derive secular solutions to the spin-spin equations,

$$\left\langle \frac{d\hat{s}_1}{dt} \right\rangle = -\nu_1 \bar{\mathcal{B}}(\hat{s}_2) \times \hat{s}_1, \qquad (B.1a)$$

$$\left\langle \frac{d\hat{s}_2}{dt} \right\rangle = -\nu_2 \bar{\mathcal{B}}(\hat{s}_1) \times \hat{s}_2, \qquad (B.1b)$$

and

$$\bar{\mathcal{B}} = \mathbb{I} - 3 \hat{e}_z \otimes \hat{e}_z,$$

in a neighbourhood of the stable equilibrium point. The solution obtained here is in fact more generally valid for any configuration close to the poles $(\hat{s}_1, \hat{s}_2) \approx (\sigma_1 \hat{e}_z, \sigma_2 \hat{e}_z)$, including the unstable aligned orthogonal case (see Sect. 4). As in the main text, we consider two indices $\ell, m \in \{1, 2\}$, with $m \neq \ell$, representing the pair of binaries permuted in some order.

In order to solve the system (B.1), we take advantage of the axial symmetry of the physical system around $\hat{e}_z$. We split the Euclidean vector space $E_3$ into an orbital plane $\Pi = \{x \hat{e}_x + y \hat{e}_y; (x, y) \in \mathbb{R}^2\}$ plus a normal line $\Lambda = \{z \hat{e}_z; z \in \mathbb{R}\}$ such that $E_3 = \Pi \oplus \Lambda$. Rotations of the plane $\Pi$ leave (B.1) invariant.

We identify the line $\Lambda$ with the reals and the orbital plane $\Pi$ with the complex plane by introducing the linear isomorphism $\Phi : \Pi \times \Lambda \to \mathbb{C} \times \mathbb{R}$, which satisfies$^2$:

$$\Phi(x \hat{e}_x + y \hat{e}_y + z \hat{e}_z) = (x + iy, z), \qquad (B.2)$$

where $i$ is the imaginary number. Eq. (B.1) can be expressed in the space $\mathbb{C} \times \mathbb{R}$ by declaring a new set of spin variables, which we define uniquely from $\Phi(\hat{s}_\ell) = (p_\ell, z_\ell)$. More explicitly, $p_\ell$ corresponds to the orbital component of $\hat{s}_\ell$, and $z_\ell$ corresponds to the projection of $\hat{s}_\ell$ on the basis element $\hat{e}_z$:

$$p_{\ell} = (\hat{s}_{\ell} \cdot \hat{e}_{x}) + i (\hat{s}_{\ell} \cdot \hat{e}_{y}), \qquad (B.3)$$

$$z_{\ell} = (\hat{s}_{\ell} \cdot \hat{e}_{z}). \qquad (B.4)$$

This identification allows us to leverage the rotational symmetries of $\Pi$ through the algebraic structures of the complex numbers. In particular, through the isomorphism, vector dot- and cross-products in $E_3$ can be described in terms of complex multiplication and conjugation$^3$.

For each body, we obtain the new and equivalent form of the spin precession equations:

$$\frac{d\mathbf{p}_\ell}{dt} = i \nu_\ell (z_\ell \mathbf{p}_m + 2z_m \mathbf{p}_\ell), \qquad (B.5)$$

$$\frac{d\mathbf{z}_\ell}{dt} = \nu_\ell \mathfrak{Im}(\mathbf{p}_\ell^* \mathbf{p}_m), \qquad (B.6)$$

where $\Re$, $\Im$ are the real and imaginary parts, and the operation $\mathbf{p} \mapsto \mathbf{p}^*$ denotes complex conjugation. These equations are coupled by the four parameters $(p_1, p_2, z_1, z_2)$ with values in $\mathbb{C}^2 \times \mathbb{R}^2$, for a total dimensionality of 6. As discussed in Sect. 3, two degrees of freedom are redundant, constrained by the unit-norm conditions

$$|\mathbf{p}_\ell|^2 + z_\ell^2 = 1, \qquad (B.7)$$

and the magnetic interaction energy defines a conserved quantity:

$$\bar{U}_B = 2z_\ell z_m - \Re(\mathbf{p}_\ell^* \mathbf{p}_m). \qquad (B.8)$$

The anti-symmetry of the equation (B.6) with respect to a swap of indices $\ell \leftrightarrow m$ allows us to determine another first integral for the problem — namely, the z-component of the total intrinsic angular momentum, $S_z = s_1 z_1 + s_2 z_2$. We introduce the quantity:

$$\zeta = \nu_2 z_1(t) + \nu_1 z_2(t) = \frac{\mu_0 \mu_1 \mu_2}{8\pi a^3 (1 - e^2)^{3/2}} \frac{1}{s_1 s_2} S_z. \qquad (B.9)$$

Since for any two complex numbers $\mathfrak{Im}(u^*v + v^*u) = 0$, it follows that the derivative of $\zeta$ vanishes. We also introduce the anti-symmetric spin:

$$Q_z(t) = \nu_2 z_1(t) - \nu_1 z_2(t), \qquad (B.10)$$

which is not a conserved quantity.

$$u \cdot v = \Re(\Phi_1^*(u)\Phi_1(v)),$$

$$\Phi_2(u \times v) = \mathfrak{Im}(\Phi_1^*(u)\Phi_1(v)).$$

with $\Phi = (\Phi_1, \Phi_2)$ and $\Phi_1^*$ the complex conjugate of $\Phi_1$.

<footnote>$^2$ The choice of real and imaginary axes is arbitrary. There is nothing special about $\hat{e}_x$ and $\hat{e}_y$. We could have equally taken any other pair of orthogonal axes in the orbital plane, and obtained equivalent results.</footnote>

<footnote>$^3$ For example, for two vectors contained in the orbital plane, $u, v \in \Pi$, the dot- and cross-product operations under the isomorphism are simply</footnote>
The spins of the two bodies can be decoupled in Eqs. (B.5) and (B.6) via differentiation and substitution of the constraints. We obtain the second-order system of equations below, purely in terms of the variables $\mathbf{p}_1, \mathbf{p}_2$, and $Q_z$:

$$ \frac{d^2 \mathbf{p}_1}{dt^2} = \alpha_+ \frac{d\mathbf{p}_1}{dt} + \beta_- \mathbf{p}_1, \qquad (B.11) $$

$$ \frac{d^2 \mathbf{p}_2}{dt^2} = \alpha_- \frac{d\mathbf{p}_2}{dt} + \beta_+ \mathbf{p}_2, \qquad (B.12) $$

$$ \frac{d^2 Q_z}{dt^2} = a - b Q_z + 3Q_z^3, \qquad (B.13) $$

with $\alpha_+, \alpha_-, \beta_+, \beta_-$ four functions defined via

$$ \alpha_{\pm}(t) = 2i\zeta \pm \frac{1}{\zeta + Q_z(t)} \frac{dQ_z}{dt}(t), \qquad (B.14) $$

$$ \beta_{\pm}(t) = \frac{3}{4} (\zeta^2 - Q_z^2(t)) \pm \frac{2i\zeta}{\zeta + Q_z(t)} \frac{dQ_z}{dt}(t), \qquad (B.15) $$

and the constants $a \in \mathbb{R}, b \in \mathbb{R}^+$ given by

$$ a = (v_2^2 - v_1^2) \zeta, \qquad b = v_1^2 + v_2^2 + 4\bar{U}_B v_1 v_2 - \frac{3}{2} \zeta^2. \qquad (B.16) $$

To solve the full system of equations, we must first solve for $Q_z$ and then plug the obtained solution in the expression of $\alpha_\pm$ and $\beta_\pm$ in order to determine $\mathbf{p}_1$ and $\mathbf{p}_2$. In fact, there is an exact analytical solution for $Q_z$ in terms of elliptic functions. The resulting expression for $Q_z$ is somewhat lengthy, and subsequently solving (B.11), (B.12) in this scenario proves to be challenging. Instead, we opt to present simpler and physically meaningful expressions for all the parameters, with domain of validity close to the poles.

We consider the expansion of Eqs. (B.11) - (B.13) around some origin $\bar{Q}_z$, with $Q_z(t) = \bar{Q}_z + \delta Q_z(t)$. For a good choice of $\bar{Q}_z$, the resulting system may be truncated at low orders in $\delta Q_z$ to yield low-order solutions $Q_z^{(0)}$, $Q_z^{(1)}$, etc. In particular, we constrain $\bar{Q}_z$ to the domain of $Q_z$ by choosing $\bar{Q}_z = Q_z(t_0)$ for a reference time $t_0$. The error incurred from this expansion will depend on the largeness of the variations $\delta Q_z$ — which are assumed small in the neighbourhood of the poles. Indeed, we consider some arbitrary spin variation from $t_0$ to $t$, namely $\delta \hat{s}_\ell(t) = \hat{s}_\ell(t) - \hat{s}_\ell(t_0)$. Due to the geometry of the unit sphere, this variation will propagate along the z component to some $|\delta z_\ell(t)| \le |\delta \hat{s}_\ell(t)| \sin \theta_\ell \le 2 \sin^2 \theta_\ell$, where $\theta_\ell$ is the maximum attained polar inclination, that is, $\sin \theta_\ell = \sup_\ell |\hat{s}_\ell(t) \cdot \hat{e}_z|$. For a solution close to the poles, the polar angle $\theta_\ell$ is assumed a small parameter. In this case, the perturbation $\delta Q_z(t) = v_2 \delta z_1(t) - v_1 \delta z_2(t)$ will also remain similarly bounded. With these considerations, we present below the solutions for low orders of $\delta Q$.

### Orthogonal component

As discussed, we begin by determining an expression for the decoupled variable $Q_z$ [Eq. (B.13)]. Until this point arbitrary, we fix the choice of $\bar{Q}_z$ to best fit (B.13) and minimise $\delta Q_z$. A good expansion parameter $\bar{Q}_z$ is in fact the zero-order constant solution $Q_z^{(0)}$, the unique real value that satisfies the third-degree algebraic equation:

$$ a - b Q_z^{(0)} + 3(Q_z^{(0)})^3 = 0, \qquad (B.17) $$

whereas the first-order solution is

$$ Q_z^{(1)}(t) = Q_z^{(0)} + \rho_z \cos(\omega_z t + \phi_z), \qquad (B.18) $$

with

$$ \omega_z = \left(b + \frac{9}{2} \bar{Q}_z^2\right)^{1/2}, \qquad (B.19a) $$

$$ \phi_z = \arctan2 \left( Q_z(0) - \bar{Q}_z, -\frac{1}{\omega_z} \frac{dQ_z}{dt}(0) \right), \qquad (B.19b) $$

$$ \rho_z = \sqrt{\left(Q_z(0) - \bar{Q}_z\right)^2 + \left(\frac{1}{\omega_z} \frac{dQ_z}{dt}(0)\right)^2}. \qquad (B.19c) $$

We see a posteriori that the choice of $\bar{Q}_z = Q_z^{(0)}$ is in fact the average of the first-order solution $Q_z^{(1)}$.

### Orbital component

For the orbital part, we consider a zero-order approximation in $\delta Q$. In this scenario, each component $\mathbf{p}_\ell^{(0)}$ obeys a constant-coefficient linear ordinary differential equation, with

$$ \alpha_+^{(0)}(t) = \alpha_-^{(0)}(t) = 2i\zeta, \qquad (B.20) $$

$$ \beta_+^{(0)}(t) = \beta_-^{(0)}(t) = \frac{3}{4} (\zeta^2 - \bar{Q}_z^2). \qquad (B.21) $$

The corresponding orbital plane solutions are given by superposed complex rotations:

$$ \mathbf{p}_1^{(0)}(t) = c_1^+ e^{i\omega_p^+ t} + c_1^- e^{i\omega_p^- t}, \qquad (B.22a) $$

$$ \mathbf{p}_2^{(0)}(t) = c_2^+ e^{i\omega_p^+ t} + c_2^- e^{i\omega_p^- t}, \qquad (B.22b) $$

with two oscillating frequencies

$$ \omega_p^{\pm} := \zeta \pm \frac{1}{2} \sqrt{\zeta^2 + 3\bar{Q}_z^2} \quad \in \mathbb{R}^+, \qquad (B.23) $$

and four constants $c_1^+, c_1^-, c_2^+, c_2^-$ in the unit complex disk $D = \{z \in \mathbb{C} : |z| \le 1\}$; they can be computed from the initial conditions as:

$$ c_{\ell}^{\pm} = \frac{\mp 1}{\omega_{\ell}^{+} - \omega_{\ell}^{-}} \left( i \frac{d\mathbf{p}_{\ell}}{dt}(0) + \omega_{\ell}^{\mp} \mathbf{p}_{\ell}(0) \right). \qquad (B.24) $$

The calculations may naturally be extended to first-order solutions $\mathbf{p}_\ell^{(1)}$, which include a higher number of harmonic frequencies.

We now return to the original Euclidean space $E_3$. The spin can be broken down into the corresponding components of the basis $e_0$ by inverting the isomorphism $\Phi$:

$$ \hat{s}_{\ell} = \Phi^{-1}(\mathbf{p}_{\ell}, z_{\ell}) = \mathrm{Re}(\mathbf{p}_{\ell}) \hat{e}_{x} + \Im(\mathbf{p}_{\ell}) \hat{e}_{y} + z_{\ell} \hat{e}_{z}. $$

From the above solutions, the component $z_\ell$ of the spins can then be retrieved from (B.18) via the relations

$$ z_1^{(1)}(t) = \frac{1}{2\nu_2} (\zeta + Q_z^{(1)}(t)), \quad z_2^{(1)}(t) = \frac{1}{2\nu_1} (\zeta - Q_z^{(1)}(t)), \qquad (B.25) $$

and the two orbital components:

$$ \mathrm{Re}(\mathbf{p}_\ell(t)) = \rho_\ell^+ \cos(\omega_\ell^+ t + \phi_\ell^+) + \rho_\ell^- \cos(\omega_\ell^- t + \phi_\ell^-), \qquad (B.26a) $$

$$ \Im(\mathbf{p}_\ell(t)) = \rho_\ell^+ \sin(\omega_\ell^+ t + \phi_\ell^+) + \rho_\ell^- \sin(\omega_\ell^- t + \phi_\ell^-), \qquad (B.26b) $$

with $\rho_\ell^+ = |c_\ell^+| \in [0, 1]$ and $\phi_\ell^+ = \arg c_\ell^+ \in [0, 2\pi]$.