Open Access
Volume 675, July 2023
Article Number A32
Number of page(s) 15
Section Celestial mechanics and astrometry
Published online 30 June 2023

© The Authors 2023

Licence Creative CommonsOpen Access article, published by EDP Sciences, under the terms of the Creative Commons Attribution License (, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

This article is published in open access under the Subscribe to Open model. Subscribe to A&A to support open access publication.

1 Introduction

Close to 10% of early-type massive main sequence (MMS) stars host stable large-scale magnetic fields, ranging from 3 × 102 to 3 × 104 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 103 and 109 G. For neutron stars (NS), surface field strengths are gauged of the order of 108–1013 G in classical radio pulsars (Phinney & Kulkarni 1994), and reach up to 1014–1015 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), 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; 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, Ferrario et al. 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 occurance 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 ≥ 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 û = u/u. Whenever two vectors are parallel, we symbolise this relationship by u || υ; and similarly, two perpendicular vectors are denoted by uυ. We represent the vector space dual to U by a starred U*. In this setting, we denote by an underscore uU* the associated canonical dot-product covector, that is, the linear functional from U to ℝ that satisfies u : υ → (u · υ). For two vector spaces U and V, we represent their Cartesian product by U × V = {(u, υ); uU, υV} and their tensor product by UV = {uυ; uU, υV}. Finally, whenever two vector subspaces are disjoint UV = {0}, their sum is direct and is denoted by UV = {u + υ | uU, υ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 ∊ {1, 2}, used throughout the paper, and we refer to them as the ‘primary’ and ‘secondary’, respectively. For each star we introduce a position x, a mass m, a radius R, a magnetic field B, and an intrinsic angular momentum (or spin) s. We place ourselves in the reference frame of the centre-of-mass (CM) of the system, to which we attach a right-handed basis e0 = (êx, êy, êz) spanning the Euclidean tangent space E3 ≅ ℝ3. The elements of e0 are chosen such that êx points towards the direction of closest approach, êz is orthogonal to the orbital plane, and êy completes the basis. In this frame, the system can be viewed as an effective one-body problem, parametrised by the relative separation r = x2x1. In the absence of non-Keplerian perturbations, the CM frame will be inertial and the elements of e0 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 as a centred dipole, given in function of some point x outside the surface of the star: (1a) (1b)

where µ is the magnetic dipole moment of each star, and µ0 is the vacuum permeability. The primary will feel the field B2 of its companion at relative position x = −r, and the secondary at x = r. It is convenient to express these fields in terms of the linear map ℬr : E3E3, which acts on the magnetic dipole moments according to: (2)

We note that ℬr can be identified with the symmetric tensor field with values in E3, (3)

with I denoting the identity in E3.

As is observed in the known doubly magnetic MMS binary system e-Lupi (Shultz et al. 2015; Pablo et al. 2019), we examine the particular case where the fields are symmetric about the star’s axis of rotation, with alignment between the magnetic dipole moment and the spin1:


In practice, the magnetic moment amplitudes may be expressed in terms of the dipolar field evaluated at the poles (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: (5)

The magnetic field of each star will interact with the dipole of the companion, inducing the following torque: (6a) (6b)

Simultaneous contributions due to gravity exist. We introduce the total torque felt by , (7)

where , and 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 γ: (8)

where and are the contributions to γ 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 is dominant. More explicitly, we assume and , 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 , where 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 γ is (9)

Equation (9) shows that the ratio 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 , in which case we formally have no figure torques. In practice, the cutoff to below which rotation is driven by magnetism is given by (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 , where mm is the mass of the tide-inflicting body, is the gravitational Love number of , and Q is the tidal dissipation quality factor (see e.g. Poisson & Will 2014; Strugarek et al. 2017). Then, (11)

and the corresponding cutoff criteria for is: (12)

Table 1 shows typical values for and in MMS, WD, and NS systems. At such distance scales (a ~ 108–109 km), the thresholds for k2/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 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 ; these torques are orthogonal to the intrinsic angular momentum, and hence the magnitude s must be conserved. Substituting each term and normalising by s we obtain coordinate-free spin equations: (13a) (13b)

with α = µ0µ1µ2/(4πs). We denote Eq. (13) the ‘instantaneous’ precession system. The spin axes are constrained to the unit sphere 𝒮2, yielding a total of four degrees of freedom for the coupled system of equations. Each spin precesses around a (time-dependent) axis determined by the field direction – where 𝔪 ∊ {1, 2}, m ≠ is the index of the companion star – and have Larmor frequencies given by (14)

The orbital dynamics will induce periodic fluctuations on the separation r, causing the precession axis and frequency to oscillate in time with period Porb. 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 Porb, manifesting in the ‘wobbles’ of the axis and in the modulation of the frequency ω; as well as (2) a timescale τ due to an average precession rate 〈ω〉, which we define as: (15)

where 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.

thumbnail 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 ŝ of each star ( ∊ {1, 2}) and the basis e0 = (êx, êy, êz).

Table 1

Typical binary system physical parameters (high-field range), and the corresponding dimensionless parameters.

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 Porb. 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 Porb, 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 Porb, 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 (Porb and τ) 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 ω. A larger precession rate (or smaller precession timescale τ) implies faster changes to the axis ŝ due to stronger torques. Thus, one may explicitly compute the characteristic time-ratio between the orbital timescale Porb and the precession timescale τ: (16)

where P is the rotational period of body 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 and . 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: (17)

In all three types of systems considered (Table 1), we obtain very low values for ϵ, namely ϵMMS ~ 10−9, ϵWD ~ 10−7, and ϵNS ~ 10−6. We therefore place ourselves in the scenario ϵ ≪ 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 ϵ ≳ 1 the two timescales cannot be decoupled in this manner. In this way, the time-ratio parameter ϵ 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 FB ~ 3µ0µ1µ2/(4πa4) and the gravitational force FN ~ Gm1 m2/a2: (18)

Even for the most magnetic systems considered, this ratio is limited to FB/FN ≲ 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 (êx, êy, êz) is inertial, the orbits are elliptical, and the separation can be parametrised as an ellipse in the centre of mass frame: (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): (20a) (20b)

where arctan2 represents the 2-argument inverse tangent, E is the eccentric anomaly, and is the mean angular motion.

To formalise our previous argument, we consider the binary system at some given instant t = t0 + h, with a short-timescale variation h ϵ [0, Porb]. In essence, we are considering that h parametrises a single full orbital revolution of the binary, beginning at some time t0. In this timeframe, the variations in the spin axis are bounded by (21a) (21b)

where the supremum of the derivative of ŝ can be obtained from the precession Eq. (13), with indices ∊ {1, 2} and 𝔪 ∊ {1, 2}, where 𝔪 ≠ : (22)

It is easy to see (cf. Eq (3)) that equality can be reached when , and r = rmin = 1/(a(1 − e)).

Equation (13) may then be evaluated at time t = t0 + h, and using the bounds obtained, expanded via ŝ (t0 + h) = ŝ (t0) + O(ϵ), from whence: (23)

where we have explicitly included the temporal dependence of r = r(t0 + h) in the subscript of ℬr, and defined ϵ = max (ϵ1, ϵ2). Since the above equation is valid for any h ϵ [0, Porb], integrating over a full orbit yields the secular spin equation (24)

where the orbital averaging operator is defined for some function of time ξ as (25)

where is the description of ξ in terms of the true anomaly, and the Jacobian factor is found via implicit differentiation of Eq. (20): (26)

In the Keplerian scenario the orbits are fixed, and the linear map ℬr – which depends purely on the radial separation r – is therefore Porb-periodic. Consequently, the orbital average (ℬr) will be constant, independent of the secular time t0. Plugging the expressions of ℬr (Eq. (3)) and of the separation r (Eq. (19)) into Eq. (25), we obtain the effective field tensor in the orbital frame basis (êx, êy, êz): (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 (ℬr). As a consequence, the radial component of ℬr has been suppressed, leaving a dipolar effective field with a predominant component in the orthogonal direction êz, plus a weaker component aligned with the spin direction.

We denote the normalised term on the right-hand side of Eq. (27) by . The final secular form of the spin precession equations is obtained by absorbing the constants together: (28a) (28b)

where we have introduced the magnetic rotational frequencies ν, coupling magnetism, spins, and orbital parameters: (29)

The above system has four degrees of freedom, two for each spin, since the magnitudes of ŝ 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 t0 it is possible to reduce the system to a single dimensionless parameter κ = ν2 /ν1, which depends only on the ratio between the two spin magnitudes. The parameter κ 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 ê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.

thumbnail Fig. 2

Sample trajectories for the secular evolution of the spin axis of the primary ŝ1, plotted against the unit sphere for different values of κ = ν2 /ν1. The initial conditions are fixed, and shown for the primary as a dashed line from the origin to ŝ1(0). The colours are interpolated between blue and red from initial time to respectively. The black axis represents the direction of êz. The behaviour for the secondary is analogous, up to a swap of initial conditions and κ → κ−1 = ν1 /ν2.

4 Equilibrium states and stability

A ‘secular’ equilibrium state at time t0 is a pair of spins defined on an orbital period (ŝı, ŝ2) : [t0, t0 + Porb] → E3 × E3 such that the average change in intrinsic angular momentum is zero: (30)

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: (31)

This corresponds to a singular-value problem, which can be solved with the explicit matrix form of in the orbital-frame basis. Two classes of solutions can be determined from the singular vectors of , 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.

thumbnail 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.

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: (32)

for some unit vector in the plane of the orbit (i.e. ). The two unitary parameters (σ1, σ2) ∊ {−1, 1} × {−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 (ŝ1, ŝ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 êz: (33)

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

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): (34)

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

We denote ŪB the ‘secular magnetic energy’. As can be straightforwardly verified, ŪB is a constant of motion of the secular system – its time derivative vanishes for any pair (ŝ1, ŝ2) ofaxes 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 Ū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 Ū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 : 𝒮n ×𝒮n → ℝ is bounded by its largest-magnitude eigenvalue. Applying the principle to the secular magnetic energy (cf. Eq. (35)), one obtains the bounds −2 ≤ ŪB ≤ 2. The spin configurations in which these bounds are actually attained correspond to ŝ along the direction of the associated eigenvector – that is, ŝ || ê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 antiparallel case (σ2 =−σ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 Ū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 = (x1,…, xn) with all partial second-order derivatives. The Hessian of f is a matrix-valued function ℋ : ℝn → ℳn(ℝ) defined as follows: (36)

Evaluating the Hessian at some point x ∊ ℝn provides a description of the local convexity of f at that point. If x* is a critical point (∇f (x*) = 0) then f can be locally approximated by a quadratic function (37)

The Hessian matrix at critical point 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 ŪB as a function of spin direction. For the following, we define the relative alignment between spins ŝ1 and ŝ2: (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 in Eq. (32) the energy takes the same value , depending only on the relative orientation of the two spins. The sign of σ = σ1σ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 in fact reflects the more general axis-symmetry of the system with respect to rotations around the axis ê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 lying in the orbital plane over another corresponds to rotating both spins simultaneously by the angle that takes . The set of all unit vectors 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. 𝒮2 × 𝒮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, (39)

with azimuth ψ and angle θ measured from the north pole. In these coordinates, the secular energy takes the form (40)

As discussed above, we can see clearly that simultaneous rotations in the azimuth ψ1 and ψ2 do not affect the energy. We can therefore look at the difference ψ1ψ2 := ∆ψ and discard the redundant degree of freedom ψ1 + ψ2. At the equilibrium, the polar angles are θ1 = θ2 = π/2, and the azimuth is either ∆ψ = 0 or ∆ψ = π. We consider therefore the three directions of the tangent space of this reduced manifold. In this case, the Hessian matrix of the energy ŪB with respect to the three variables (∆ψ, θ1, θ2) has value: (41)

The sign of σ = ±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σ). 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 ŝ 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, y, which with the unit-norm condition gives (42)

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

The corresponding basis for the tangent space is the coordinate basis dx and dy. With respect to these coordinates the Hessian of the energy takes the values: (44)

The nature of the eigenvalues of the Hessian change according to the sign of σ: when the spins are aligned (σ = +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 (σ = −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 ŝ = x êx + y êy + z êz, where the orbital-plane components satisfy (45a) (45b)

and the orthogonal component satisfies (46a) (46b)

with real parameters , , ρz, ϕz, dependent on initial conditions. The expressions of the constants ζ and and the frequencies ωz, , and are given in the appendix.

thumbnail 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 ŝ 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 ŝ (0). The orbital plane is represented in a darker shade with black contours, and the axis ê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 (êz, − êz).

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: (47a) (47b)

where we have used the re-scaled dimensionless time and κ = ν2/ν1 = s1/s2 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, Δτ ~ 0.1 κ−1/2. This provides us with dissipative effects visible on the simulation timescales.

In order to assess stability, we consider two metrics: the secular magnetic interaction energy ŪB (see Eq. (35)); and the angular distance to the stable equilibrium point (êz, − êz), which we define as: (48a) (48b)

where arg(u, v) = arctan2 (|u ×v|, u · v) is the angle between any two vectors u and v.

The dimensionless system is then integrated until convergence, for four sets of initial spin conditions (Table 2), and with spin ratio κ = 0.3. Each initial condition corresponds to an unstable equilibrium state plus a small perturbation of the order of ~1º. For completeness, we also include a perturbation of the stable state (~10º). The resulting trajectories are plotted in Fig. 4, together with the time evolution of the secular energy ŪB and of the angular distances d1 and d2. 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.

Table 2

Initial conditions in four distinct stability simulations (a–d).

thumbnail 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 ℱ{ŝ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.

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 κ = 0.3. The following initial conditions are considered (polar coordinates): inclinations θ1 = 10° and θ2 = 172.5° from the north pole, and azimuths ψ1 = 0 and ψ2 = 50°. Figure 5 compares the obtained analytical expression of the primary to numerical integration, decomposed in the Cartesian basis ŝ1 (τ) = x1 (τ) êx + y1 (τ) êy + z1 (τ) êz. Matching peaks can be observed in the power spectra of both solutions, at the obtained angular frequencies (in dimensionless units) rad and rad for the orbital components (x1, y1), and at ωz = 2.34 rad for the orthogonal component z1 (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 ϵ-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 UB (Eq. (34)) takes the form: (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, t1 and t2 such that the true anomaly at each instant equals f (t1) = (ψ1 + ψ2)/2 and f (t2) = (ψ1 + ψ2 + π)/2. Then the difference in energy between those two instants is roughly (50)

where we have substituted r ~ a. We conclude that instantaneous equilibrium positions will indeed develop large energy oscillations in the timescale t ~ Porb, particularly when the polar angle θ 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:

Such expression has been normalised by the positive constant µ0µ1µ2/8πb3 as discussed in Sect. 4. We similarly normalise the expression of the instantaneous magnetic energy (cf. Eq. (34)) by the scalar µ0µ1µ2/4πr3 and obtain:

We observe that the secular averaging procedure effectively produced a flip in the sign of the energy as well as a switch . Analogously to how the secular equilibrium states of the binary are given by the singular vectors of (cf. Sect. 4), the instantaneous states are given by the singular vectors of . Consequently, each state of instantaneous equilibrium has a secular counterpart. These correspond to spin axes aligned with the radial direction (resp. êz) or perpendicular to (resp. êz). The stability of each state depends on the local convexity of the energy Ub (resp. Ub). We present in Table 3 a comparison between the obtained states for the two types of equilibrium.

Table 3

Comparison between the states of instantaneous and of secular equilibrium.

6.2 ϵ-Lupi

We consider as a potential application case the ϵ-Lupi inner binary system. ϵ-L·upi is a ternary system composed of two close-range B-type companions Aa and Ab, plus a third distant companion dubbed ϵ-L·upi B. Both stars of the ϵ-L·upi 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 Porb ~ 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, eccentricity of e = 0.28, inclination ι = 21°, stellar masses of m1 = 9.0 M for the primary and m2 = 7.9 M for the secondary, and radii R1 = R2 = 4.5 R.

Shultz et al. (2015) reports a dipolar field strength of at least G and G at the surface poles of the star, and projected rotational velocities at the equator v1 sin i1 = 37 km s−1 for the primary and v2 sin i2 = 27 km s−1 for the secondary. By adopting the rotational inclination for each star equal to the orbital plane inclination ι ~ 21°, we obtain rotational periods of the order of P1 = 2.2 d and P2 = 3.0 d.

From these physical parameters we may calculate the corresponding dimensionless ratios. The impact of figure effects can be assessed via , whereas for tides and . These expressions suggest that if ϵ-·upi is both highly symmetrical with , and has tidal parameters , then the system's rotation is likely driven by magnetism. In this case, the smallness of η = 5 × 10−12 and ϵ = 4 × 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 ϵ-·upi 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 ±(êz, −ê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 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.35M, 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 ~ 0.35 M, B ~ 1 kG) the dimensionless parameters that measure the strength of extended-body gravitational interactions with respect to magnetic torques: γfig = 6 × 1011 J2 and γtides = 7 × 107k2/Q (Eqs. (9), (11)). Subsequently, we obtain the (quite tight) cutoffs J2 ≤ 10−12 for the dimensionless quadrupole moment (Eq. (10)), and k2/Q ≤ 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 (µ < 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 is much shorter than the orbital period Porb. In such regime, the problem may be hierarchically split into three distinct timescales PPorbτ, where τ 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 Porb, 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).


C.A. acknowledges the joint finantial support of Centre National d’Études Spatiales (CNES) and École Doctorale Astronomie et Astrophysique d’Île-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.

Appendix A Optimisation of bilinear forms

In this appendbix we recall the variational characterisation of the singular value decomposition. Consider the bilinear form U taking two unit vectors U : 𝒮n × 𝒮n → ℝ, represented under a basis {ê1,…, ên} via some matrix Uij: (A.1)

where û = uiêi, , and sum over repeated indices is presupposed.

Since U is continuous on the compact domain 𝒮n × 𝒮n, it attains both a maximum and minimum value. Indeed, consider the Lagrange function (A.2)

with multipliers λ1 and λ2. The extrema û* and , necessarily satisfy the stationary condition , which implies that for any tangent vectors h1, h2, (A.3a) (A.3b)

From the matrix representation of U [Eq. (A.1)], we see that the solutions û* and correspond respectively to the left and right singular vectors of Uij — that is, the unique set of vectors that satisfy (A.4a)

where λ = 2λ1 = 2λ2 is the corresponding singular value of Uij.

It is direct to see that the attained value . From all singular-vector pairs , the global maximum (or minimum) of U is therefore reached by the candidate pair that has the highest (or lowest) corresponding singular value λ. We note that whenever U is symmetric, the singular vectors and eigenvectors of Uij 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, (B.1a) (B.1b)


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 (ŝ1, ŝ2) ≈ (σ1êz, σ2êz), including the unstable aligned orthogonal case (see Sect. 4). As in the main text, we consider two indices , m ∊ {1, 2}, with m ≠ , 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 êz. We split the Euclidean vector space E3 into an orbital plane Π = {x êx + y êy; (x,y) ∊ ℝ2} plus a normal line Λ = {zêz; z ∊ ℝ} such that E3 = Π ⊕ Λ. Rotations of the plane Π leave (B.1) invariant.

We identify the line Λ with the reals and the orbital plane Π with the complex plane by introducing the linear isomorphism Φ : Π × Λ → ℂ × ℝ, which satisfies2: (B.2)

where i is the imaginary number. Eq. (B.1) can be expressed in the space ℂ × ℝ by declaring a new set of spin variables, which we define uniquely from Φ(ŝ) = (p, z). More explicitly, p corresponds to the orbital component of ŝ, and z corresponds to the projection of ŝ on the basis element êz: (B.3) (B.4)

This identification allows us to leverage the rotational symmetries of Π through the algebraic structures of the complex numbers. In particular, through the isomorphism, vector dot-and cross-products in E3 can be described in terms of complex multiplication and conjugation3.

For each body, we obtain the new and equivalent form of the spin precession equations: (B.5) (B.6)

where 𝔑𝔢, ℑ𝔪 are the real and imaginary parts, and the operation pp* denotes complex conjugation. These equations are coupled by the four parameters (p1, p2, z1, z2) with values in ℂ2 × ℝ2, for a total dimensionality of 6. As discussed in Sect. 3, two degrees of freedom are redundant, constrained by the unit-norm conditions (B.7)

and the magnetic interaction energy defines a conserved quantity: (B.8)

The anti-symmetry of the equation (B.6) with respect to a swap of indices ↔ m allows us to determine another first integral for the problem — namely, the z-component of the total intrinsic angular momentum, Sz = s1z1 + s2z2. We introduce the quantity: (B.9)

Since for any two complex numbers ℑ𝔪(u*v + v*u) = 0, it follows that the derivative of ζ vanishes. We also introduce the anti-symmetric spin: (B.10)

which is not a conserved quantity.

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 p1, p2, and Qz: (B.11) (B.12) (B.13)

with α+, α_,β+,β_ four functions defined via (B.14) (B.15)

and the constants a ∊ ℝ, b ∊ ℝ+ given by (B.16)

To solve the full system of equations, we must first solve for Qz and then plug the obtained solution in the expression of α± and β± in order to determine p1 and p2. In fact, there is an exact analytical solution for Qz in terms of elliptic functions. The resulting expression for Qz 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 , with . For a good choice of , the resulting system may be truncated at low orders in δQz to yield low-order solutions , , etc. In particular, we constrain to the domain of Qz by choosing for a reference time t0. The error incurred from this expansion will depend on the largeness of the variations δQz — which are assumed small in the neighbourhood of the poles. Indeed, we consider some arbitrary spin variation from t0 to t, namely δŝ(t) = ŝ(t) − ŝ(t0). Due to the geometry of the unit sphere, this variation will propagate along the z component to some |δz(t)| ≤ |δŝ(t)| sin θ ≤ 2 sin θ, where θ is the maximum attained polar inclination, that is, sin θ = supt |ŝ(t) · êz|. For a solution close to the poles, the polar angle θ is assumed a small parameter. In this case, the perturbation δQz(t) = v2δz1(t) − v1δz2(t) will also remain similarly bounded. With these considerations, we present below the solutions for low orders of δQ.

Orthogonal component

As discussed, we begin by determining an expression for the decoupled variable Qz [Eq. (B.13)]. Until this point arbitrary, we fix the choice of to best fit (B.13) and minimise δQz. A good expansion parameter is in fact the zero-order constant solution , the unique real value that satisfies the third-degree algebraic equation: (B.17)

whereas the first-order solution is (B.18)

with (B.19a) (B.19b) (B.19c)

We see a posteriori that the choice of is in fact the average of the first-order solution .

Orbital component

For the orbital part, we consider a zero-order approximation in δQ. In this scenario, each component obeys a constant-coefficient linear ordinary differential equation, with (B.20) (B.21)

The corresponding orbital plane solutions are given by superposed complex rotations: (B.22a) (B.22b)

with two oscillating frequencies (B.23)

and four constants , , , in the unit complex disk D = {z ∊ ℂ : |z| ≤ 1}; they can be computed from the initial conditions as: (B.24)

The calculations may naturally be extended to first-order solutions , which include a higher number of harmonic frequencies.

We now return to the original Euclidean space E3. The spin can be broken down into the corresponding components of the basis e0 by inverting the isomorphism Φ:

From the above solutions, the component z of the spins can then be retrieved from (B.18) via the relations (B.25)

and the two orbital components: (B.26a) (B.26b)

with and .


  1. Achilleos, N., & Wickramasinghe, D. T. 1989, ApJ, 346, 444 [NASA ADS] [CrossRef] [Google Scholar]
  2. Ahuir, J., Strugarek, A., Brun, A. S., & Mathis, S. 2021, A&A, 650, A126 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  3. Akgün, T., Reisenegger, A., Mastrano, A., & Marchant, P. 2013, MNRAS, 433, 2445 [Google Scholar]
  4. Alecian, E., Neiner, C., Wade, G. A., et al. 2014, Proc. Int. Astron. Union, 9, 330 [CrossRef] [Google Scholar]
  5. Amaro-Seoane, P., Audley, H., Babak, S., et al. 2017, Laser Interferometer Space Antenna (USA: NASA) [Google Scholar]
  6. Aurière, M., Wade, G. A., Silvester, J., et al. 2007, A&A, 475, 1053 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  7. Bagnulo, S., & Landstreet, J. D. 2021, MNRAS, 507, 5902 [NASA ADS] [CrossRef] [Google Scholar]
  8. Beuermann, K., Euchner, F., Reinsch, K., Jordan, S., & Gänsicke, B. T. 2007, A&A, 463, 647 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  9. Borra, E. F., Landstreet, J. D., & Mestel, L. 1982, ARA&A, 20, 191 [NASA ADS] [CrossRef] [Google Scholar]
  10. Bourgoin, A., Le Poncin-Lafitte, C., Mathis, S., & Angonin, M.-C. 2022, Phys. Rev. D, 105, 124042 [NASA ADS] [CrossRef] [Google Scholar]
  11. Braithwaite, J. 2008, MNRAS, 386, 1947 [Google Scholar]
  12. Braithwaite, J. 2009, MNRAS, 397, 763 [NASA ADS] [CrossRef] [Google Scholar]
  13. Braithwaite, J., & Cantiello, M. 2013, MNRAS, 428, 2789 [NASA ADS] [CrossRef] [Google Scholar]
  14. Braithwaite, J., & Nordlund, Å. 2006, A&A, 450, 1077 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  15. Braithwaite, J., & Spruit, H. C. 2004, Nature, 431, 819 [Google Scholar]
  16. Bromley, B. C., & Kenyon, S. J. 2022, AJ, 164, 229 [NASA ADS] [CrossRef] [Google Scholar]
  17. Browning, M. K. 2008, ApJ, 676, 1262 [Google Scholar]
  18. Brun, A. S., & Browning, M. K. 2017, Liv. Rev. Sol. Phys., 14, 4 [Google Scholar]
  19. Brun, A., Strugarek, A., Noraz, Q., et al. 2022, ApJ, 926, 21 [NASA ADS] [CrossRef] [Google Scholar]
  20. Campbell, C. G. 2018, Magnetohydrodynamics in Binary Stars (Berlin: Springer), 456 [Google Scholar]
  21. Carrier, F., North, P., Udry, S., & Babel, J. 2002, A&A, 394, 151 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  22. Carvalho, G. A., Anjos, R. C. D., Coelho, J. G., et al. 2022, ApJ, 940, 90 [NASA ADS] [CrossRef] [Google Scholar]
  23. Charbonneau, P., & MacGregor, K. B. 2001, ApJ, 559, 1094 [NASA ADS] [CrossRef] [Google Scholar]
  24. Commerçon, B., Hennebelle, P., Audit, E., Chabrier, G., & Teyssier, R. 2010, A&A, 510, L3 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  25. Cumming, A. 2002, MNRAS, 333, 589 [NASA ADS] [CrossRef] [Google Scholar]
  26. Damiani, C., & Lanza, A. F. 2015, A&A, 574, A39 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  27. de Mink, S. E., Sana, H., Langer, N., Izzard, R. G., & Schneider, F. R. N. 2014, ApJ, 782, 7 [Google Scholar]
  28. Dobler, W., Stix, M., & Brandenburg, A. 2006, ApJ, 638, 336 [NASA ADS] [CrossRef] [Google Scholar]
  29. Donati, J. F., & Landstreet, J. D. 2009, ARA&A, 47, 333 [Google Scholar]
  30. Donati, J. F., Morin, J., Petit, P., et al. 2008, MNRAS, 390, 545 [Google Scholar]
  31. Duez, V. 2011, Astron. Nachr., 332, 983 [Google Scholar]
  32. Duez, V., Braithwaite, J., & Mathis, S. 2010, ApJ, 724, L34 [Google Scholar]
  33. Ferrario, L., & Wickramasinghe, D. T. 2005, MNRAS, 356, 615 [NASA ADS] [CrossRef] [Google Scholar]
  34. Ferrario, L., Pringle, J., Tout, C., & Wickramasinghe, D. 2009, MNRAS, 400, L71 [NASA ADS] [CrossRef] [Google Scholar]
  35. Ferrario, L., de Martino, D., & Gänsicke, B. T. 2015a, Space Sci. Rev., 191, 111 [Google Scholar]
  36. Ferrario, L., Melatos, A., & Zrake, J. 2015b, Space Sci. Rev., 191, 77 [CrossRef] [Google Scholar]
  37. Ferrario, L., Wickramasinghe, D., & Kawka, A. 2020, Adv. Space Res., 66, 1025 [NASA ADS] [CrossRef] [Google Scholar]
  38. Fuller, J., & Mathis, S. 2023, MNRAS, 520, 5573 [CrossRef] [Google Scholar]
  39. Gaurat, M., Jouve, L., Lignières, F., & Gastine, T. 2015, A&A, 580, A103 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  40. Gerosa, D., Kesden, M., Sperhake, U., Berti, E., & O’Shaughnessy, R. 2015, Phys. Rev. D, 92, 064016 [CrossRef] [Google Scholar]
  41. Grunhut, J., Wade, G.A., & the MiMeS Collaboration 2013, EAS Pub. Ser., 64, 67 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  42. Hut, P. 1980, A&A, 92, 167 [NASA ADS] [Google Scholar]
  43. Hut, P. 1981, A&A, 99, 126 [NASA ADS] [Google Scholar]
  44. Isern, J., García-Berro, E., Külebi, B., & Lorén-Aguilar, P. 2017, ApJ, 836, L28 [NASA ADS] [CrossRef] [Google Scholar]
  45. Jermyn, A. S., & Cantiello, M. 2021, ApJ, 923, 104 [NASA ADS] [CrossRef] [Google Scholar]
  46. Jouve, L., Lignières, F., & Gaurat, M. 2020, A&A, 641, A13 [EDP Sciences] [Google Scholar]
  47. Kaspi, V. M., & Beloborodov, A. M. 2017, ARA&A, 55, 261 [Google Scholar]
  48. Kawka, A., & Vennes, S. 2004, IAU Symp., 224, 879 [CrossRef] [Google Scholar]
  49. King, A., Whitehurst, R., & Frank, J. 1990, MNRAS, 244, 731 [NASA ADS] [Google Scholar]
  50. Kochukhov, O. 2021, A&ARv, 29, 1 [Google Scholar]
  51. Kochukhov, O., & Lavail, A. 2017, ApJ, 835, L4 [Google Scholar]
  52. Kochukhov, O., & Shulyak, D. 2019, ApJ, 873, 69 [Google Scholar]
  53. Kochukhov, O., Lundin, A., Romanyuk, I., & Kudryavtsev, D. 2011, ApJ, 726, 24 [NASA ADS] [CrossRef] [Google Scholar]
  54. Landstreet, J. D., & Bagnulo, S. 2020, A&A, 634, L10 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  55. Landstreet, J. D., & Mathys, G. 2000, A&A, 359, 213 [NASA ADS] [Google Scholar]
  56. Landstreet, J. D., Bagnulo, S., Valyavin, G., & Valeev, A. F. 2017, A&A, 607, A92 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  57. Liebert, J. 1988, PASP, 100, 1302 [NASA ADS] [CrossRef] [Google Scholar]
  58. Liebert, J., Wickramsinghe, D. T., Schmidt, G. D., et al. 2005, AJ, 129, 2376 [NASA ADS] [CrossRef] [Google Scholar]
  59. Lira, M., Degollado, J. C., Moreno, C., & Núñez, D. 2022, Gen. Rel. Grav., 54, 146 [NASA ADS] [CrossRef] [Google Scholar]
  60. Maggiore, M., Broeck, C. V. D., Bartolo, N., et al. 2020, J. Cosmol. Astropart. Phys., 2020, 050 [CrossRef] [Google Scholar]
  61. Mathis, S., Duez, V., & Braithwaite, J. 2011, IAU Symp., 271, 270 [NASA ADS] [Google Scholar]
  62. Morel, T., Castro, N., Fossati, L., et al. 2014, Proc. Int. Astron. Union, 9, 342 [CrossRef] [Google Scholar]
  63. Moss, D. 1987, MNRAS, 226, 297 [CrossRef] [Google Scholar]
  64. Moss, D. 2003, A&A, 403, 693 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  65. Pablo, H., Shultz, M., Fuller, J., et al. 2019, MNRAS, 488, 64 [Google Scholar]
  66. Phinney, E. S., & Kulkarni, S. R. 1994, ARA&A, 32, 591 [NASA ADS] [CrossRef] [Google Scholar]
  67. Poisson, E., & Will, C. M. 2014, Gravity: Newtonian, Post-Newtonian, Relativistic (Cambridge: Cambridge University Press) [CrossRef] [Google Scholar]
  68. Pound, A., & Poisson, E. 2008, Phys. Rev. D, 77, 043524 [NASA ADS] [CrossRef] [Google Scholar]
  69. Schneider, F. R. N., Podsiadlowski, P., Langer, N., Castro, N., & Fossati, L. 2016, MNRAS, 457, 2355 [Google Scholar]
  70. Schneider, F., Ohlmann, S., Podsiadlowski, P., et al. 2019, Nature, 574, 211 [CrossRef] [Google Scholar]
  71. Schreiber, M. R., Belloni, D., Gänsicke, B. T., Parsons, S. G., & Zorotovic, M. 2021, Nat. Astron., 5, 648 [NASA ADS] [CrossRef] [Google Scholar]
  72. Shultz, M., Wade, G. A., Alecian, E., & Collaboration, B. 2015, MNRAS, 454, L1 [NASA ADS] [CrossRef] [Google Scholar]
  73. Shultz, M. E., Wade, G. A., Rivinius, T., et al. 2019, MNRAS, 490, 274 [NASA ADS] [CrossRef] [Google Scholar]
  74. Stello, D., Cantiello, M., Fuller, J., et al. 2016, Nature, 529, 364 [Google Scholar]
  75. Stibbs, D. W. N. 1950, MNRAS, 110, 395 [Google Scholar]
  76. Strugarek, A., Bolmont, E., Mathis, S., et al. 2017, ApJ, 847, L16 [Google Scholar]
  77. Tayler, R. J. 1980, MNRAS, 191, 151 [NASA ADS] [Google Scholar]
  78. Tout, C. A., Wickramasinghe, D. T., & Ferrario, L. 2004, MNRAS, 355, L13 [CrossRef] [Google Scholar]
  79. Tout, C. A., Wickramasinghe, D. T., Liebert, J., Ferrario, L., & Pringle, J. E. 2008, MNRAS, 387, 897 [NASA ADS] [CrossRef] [Google Scholar]
  80. Uytterhoeven, K., Harmanec, P., Telting, J. H., & Aerts, C. 2005, A&A, 440, 249 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  81. Vidal, J., Cébron, D., Ud-Doula, A., & Alecian, E. 2019, A&A, 629, A142 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  82. Wade, G. A., Neiner, C., Alecian, E., et al. 2015, MNRAS, 456, 2 [Google Scholar]
  83. Wickramasinghe, D. T., & Ferrario, L. 2000, PASP, 112, 873 [NASA ADS] [CrossRef] [Google Scholar]
  84. Wickramasinghe, D. T., & Ferrario, L. 2005, MNRAS, 356, 1576 [NASA ADS] [CrossRef] [Google Scholar]
  85. Will, C. M., & Maitra, M. 2017, Phys. Rev. D, 95, 064003 [CrossRef] [Google Scholar]


Throughout the text, we shall assume µ positive, but µ < 0 is also allowed, with appropriate sign changes in the equations.


The choice of real and imaginary axes is arbitrary. There is nothing special about êx and êy. We could have equally taken any other pair of orthogonal axes in the orbital plane, and obtained equivalent results.


For example, for two vectors contained in the orbital plane, u, v ∊ Π, the dot- and cross-product operations under the isomorphism are simply

with Φ = (Φ1, Φ2) and the complex conjugate of Φ1.

All Tables

Table 1

Typical binary system physical parameters (high-field range), and the corresponding dimensionless parameters.

Table 2

Initial conditions in four distinct stability simulations (a–d).

Table 3

Comparison between the states of instantaneous and of secular equilibrium.

All Figures

thumbnail 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 ŝ of each star ( ∊ {1, 2}) and the basis e0 = (êx, êy, êz).

In the text
thumbnail Fig. 2

Sample trajectories for the secular evolution of the spin axis of the primary ŝ1, plotted against the unit sphere for different values of κ = ν2 /ν1. The initial conditions are fixed, and shown for the primary as a dashed line from the origin to ŝ1(0). The colours are interpolated between blue and red from initial time to respectively. The black axis represents the direction of êz. The behaviour for the secondary is analogous, up to a swap of initial conditions and κ → κ−1 = ν1 /ν2.

In the text
thumbnail 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.

In the text
thumbnail 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 ŝ 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 ŝ (0). The orbital plane is represented in a darker shade with black contours, and the axis ê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 (êz, − êz).

In the text
thumbnail 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 ℱ{ŝ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.

In the text

Current usage metrics show cumulative count of Article Views (full-text article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.

Data correspond to usage on the plateform after 2015. The current usage metrics is available 48-96 hours after online publication and is updated daily on week days.

Initial download of the metrics may take a while.