Physics of Eclipsing Binaries. III. Spin–Orbit Misalignment

At first glance it is very tempting to think that stars in stellar systems would have their rotational axes aligned with the orbital axis: the total angular momentum during a protostellar cloud contraction is conserved, so we expect a high degree of retained symmetry. Yet this is not what we observe in nature.

Looking at our solar system alone, we see that misalignment abounds. The angle between the axes of rotation and revolution (or, conversely, between the equatorial and orbital planes) is called axial tilt or obliquity. If obliquity is 0, then the two axes are aligned. Starting with the Sun, its equator is tilted to the ecliptic by 725. Earth's equator is tilted on average 2344 and it precesses at a rate of ∼50''/yr, so obliquity changes as a function of time. Solar system planets all orbit very close to the ecliptic plane (within ∼7°), but their rotational axes are nowhere near aligned, ranging from 02 to 82° (Astronomical Almanac 2017).

Transiting exoplanet host stars exhibit a wide range in their obliquities, from completely aligned (i.e., HD 189733; Winn et al. 2006), to moderately aligned (i.e., XO-3 at 37°; Hirano et al. 2011), to perpendicular (WASP-7; Albrecht et al. 2012), and even retrograde (WASP-17; Anderson et al. 2010). Giant extrasolar planets on very close, eccentric orbits (the so-called hot and warm Jupiters) also demonstrate a large ranges and oscillatory behavior in their obliquities (Dawson & Chiang 2014). Two good reviews on the methodology and results of obliquity measurements in exoplanet systems are done Winn & Fabrycky (2015) and Triaud (2017).

Lastly, there are two shining examples for misalignment among EB stars: DI Her (Albrecht et al. 2009; Philippov & Rafikov 2013) and CV Vel (Albrecht et al. 2014). DI Her is the current record holder, with sky-projected spin–orbit angles of β_p = 72° ± 4° and β_s = −84° ± 8° for the primary and secondary stars, respectively, while CV Vel features a misaligned primary at β_p = −52° ± 6° and a (sky-projected) aligned secondary at β_s = 3° ± 7°.

The most commonly used method to measure misalignment is to acquire spectroscopic observations of a binary star during the eclipse. The measured radial velocity is a weighted average of individual radial velocities across the visible elements of the star. In a misaligned system, the transiting star no longer passes along the parallel of the eclipsed star, resulting in an asymmetry in the Rossiter–McLaughlin effect (RME; McLaughlin 1924; Rossiter 1924). Furthermore, as the obliquity affects the distortion of the eclipsed star, it will also affect the intensities and weighting of the individual eclipsed elements (most notably due to gravity brightening). Together, these deviations from the aligned case are called the anomalous RME.

The analytical formalism for binary systems with misaligned rotational and orbital angular velocity vectors has been discussed previously by Limber (1963), Kruszewski (1966), and Avni & Schiller (1982). It was not until now, however, that this effect has been built into an eclipsing binary (EB) modeling code. A few recent reports (Winn et al. 2006; Albrecht et al. 2012; Harding et al. 2013; Triaud et al. 2013) use the anomalous RME to measure obliquity, but none of these provide any technical insight into their treatment of tidally and gravitationally distorted binary systems. Our goal is to provide that here, along with a publicly available tool to model misaligned cases. We discuss different aspects of modeling binary systems with spin–orbit misalignment. We focus on properly defining the model and its parameters in Section 2 and discuss equilibrium (generalized Lagrange) points of the misaligned potential in Section 3. In Section 4, we demonstrate the misalignment treatment with a toy model as implemented in the open source software package PHOEBE (Prša et al. 2016), and on two observed systems: DI Her and Kepler-13. In the Appendix, we provide further technical details about the model, along with mathematical tools to compute poles, areas, and volumes that are used for a robust synthesis of observables (light curves, radial velocity curves, and spectral line profiles).

In this section we introduce the framework for misaligned spin axes in binaries. We use the following nomenclature: vectors are denoted with boldface (i.e., ${\boldsymbol{r}}$); vectors with unit magnitude (unit vectors) are denoted with a  $\hat{{\boldsymbol{}}}$ symbol (i.e., $\hat{{\boldsymbol{r}}}$); vector magnitudes (norms) are slanted (i.e., $r\equiv \parallel {\boldsymbol{r}}\parallel$); and fractional, unitless values are denoted with greek letters (i.e., ρ = r/a).

Let two stars (labeled A and B) rotate in the xy plane of the inertial coordinate system about the common center of the mass (labeled C) with an angular frequency ω_L. The angular momentum ${\boldsymbol{L}}$ of the binary system points in the direction of the z-axis of the inertial coordinate frame. Let star A rotate uniformly with angular velocity ${{\boldsymbol{\omega }}}_{{\rm{S}}}={\omega }_{{\rm{S}}}\hat{{\boldsymbol{S}}}$ about its center of mass, where $\hat{{\boldsymbol{S}}}$ is the unit spin vector. We introduce a canonical coordinate system with the origin at the center of star A, the x-axis pointing toward the center of star B, and the z-axis aligned with orbital axis (i.e., the z-axis of the inertial system). We denote the orthogonal vector basis of the canonical coordinate system by ($\hat{{\boldsymbol{\imath }}},\hat{{\boldsymbol{\jmath}}},\hat{{\boldsymbol{k}}}$). A schematic of the binary system and the canonical coordinate system is depicted in Figure 1.

Zoom In Zoom Out Reset image size

The lobes of the stars are defined as the surfaces of constant pressure and density. These lobes are approximated by isosurfaces of the potential V, which is written in the canonical coordinate system as:

$$\begin{eqnarray}V({\boldsymbol{r}}) & = & -\displaystyle \frac{{{GM}}_{{\rm{A}}}}{r}-\displaystyle \frac{{{GM}}_{{\rm{B}}}}{\parallel {\boldsymbol{r}}-d\hat{{\boldsymbol{\imath }}}\parallel }\\ & & +\ \displaystyle \frac{{{GM}}_{{\rm{B}}}}{{d}^{2}}(\hat{{\boldsymbol{\imath }}}\cdot {\boldsymbol{r}})-\displaystyle \frac{1}{2}{\omega }_{{\rm{S}}}^{2}\parallel (\hat{{\boldsymbol{S}}}\cdot {\boldsymbol{r}})\hat{{\boldsymbol{S}}}-{\boldsymbol{r}}{\parallel }^{2},\end{eqnarray} \tag{ 1 }$$

with ${\boldsymbol{r}}$ denoting the position and d denoting the distance between the stars. The masses of stars A and B are labeled by M_A and M_B, respectively. A detailed discussion and derivation of the potential V for a circular orbit can be found in Limber (1963) and its generalization to the non-circular case, as used here, is presented in Avni & Schiller (1982) and is summarized in Appendix A. Equation (1) can be further simplified by introducing a dimensionless potential ${\rm{\Omega }}({\boldsymbol{\rho }})$, where ${\boldsymbol{\rho }}={\boldsymbol{r}}/a$:

$$\begin{eqnarray}&&V(a{\boldsymbol{\rho }})=-\displaystyle \frac{{{GM}}_{{\rm{A}}}}{a}{\rm{\Omega }}({\boldsymbol{\rho }};q,F,\delta ,\hat{{\boldsymbol{S}}}),\end{eqnarray} \tag{ 2 }$$

where:

$$\begin{eqnarray}&&{\rm{\Omega }}({\boldsymbol{\rho }};q,F,\delta ,\hat{{\boldsymbol{S}}})=\displaystyle \frac{1}{\rho }+q\left(\displaystyle \frac{1}{\parallel {\boldsymbol{\rho }}-\delta \,\hat{{\boldsymbol{\imath }}}\parallel }-\displaystyle \frac{\hat{{\boldsymbol{\imath }}}\cdot {\boldsymbol{\rho }}}{{\delta }^{2}}\right)\\ &&\qquad \qquad \qquad \qquad +\,\displaystyle \frac{1}{2}(1+q){F}^{2}\parallel (\hat{{\boldsymbol{S}}}\cdot {\boldsymbol{\rho }})\hat{{\boldsymbol{S}}}-{\boldsymbol{\rho }}{\parallel }^{2}.\end{eqnarray} \tag{ 3 }$$

This form (for the aligned, circular, synchronous case) was first proposed by Kopal (1978) and generalized to eccentric, asynchronous orbits by Wilson (1979). Here, $q={M}_{{\rm{B}}}/{M}_{{\rm{A}}}$ is the mass ratio, F = ω_S/ω_L is the synchronicity parameter, and δ ≡ d/a is fractional instantaneous separation. We use Kepler's third law to replace ${\omega }_{{\rm{L}}}^{2}=({{GM}}_{{\rm{A}}}+{{GM}}_{{\rm{B}}})/{a}^{3}$. Note that the Kopal potential is invariant to the sign of the vector $\hat{{\boldsymbol{S}}}$: ${\rm{\Omega }}{| {}_{-\hat{{\boldsymbol{S}}}}={\rm{\Omega }}| }_{\hat{{\boldsymbol{S}}}}$.

Next, we introduce a rotated coordinate system about the x-axis w.r.t. the canonical coordinate system so that vector $\hat{{\boldsymbol{S}}}$ lies in the new xz plane. The vector basis ($\hat{{\boldsymbol{\imath }}}^{\prime} ,\hat{{\boldsymbol{\jmath}}}^{\prime} ,\hat{{\boldsymbol{k}}}^{\prime}$) of the rotated coordinate system is related to the canonical vector basis ($\hat{{\boldsymbol{\imath }}},\hat{{\boldsymbol{\jmath}}},\hat{{\boldsymbol{k}}}$) by the following relations:

$$\begin{eqnarray}\hat{{\boldsymbol{\imath }}}^{\prime} & = & \hat{{\boldsymbol{\imath }}},\\ \hat{{\boldsymbol{\jmath}}}^{\prime} & = & \hat{{\boldsymbol{\jmath}}}\cos \alpha +\hat{{\boldsymbol{k}}}\sin \alpha ,\\ \hat{{\boldsymbol{k}}}^{\prime} & = & -\hat{{\boldsymbol{\jmath}}}\sin \alpha +\hat{{\boldsymbol{k}}}\cos \alpha ,\end{eqnarray} \tag{ 4 }$$

where $\alpha =-\arctan (\hat{{\boldsymbol{S}}}\cdot \hat{{\boldsymbol{\jmath}}}/\hat{{\boldsymbol{S}}}\cdot \hat{{\boldsymbol{k}}})$. The positions are denoted by ${\boldsymbol{r}}^{\prime} =x^{\prime} \hat{{\boldsymbol{\imath }}}^{\prime} +y^{\prime} \hat{{\boldsymbol{\jmath}}}^{\prime} +z^{\prime} \hat{{\boldsymbol{k}}}^{\prime}$. Consequently, the vector $\hat{{\boldsymbol{S}}}$ can be given as:

$$\begin{eqnarray}&&\hat{{\boldsymbol{S}}}=\sin \beta \,\hat{{\boldsymbol{\imath }}}^{\prime} +\cos \beta \,\hat{{\boldsymbol{k}}}^{\prime} ,\end{eqnarray} \tag{ 5 }$$

with the angle $\beta =\arcsin (\hat{{\boldsymbol{S}}}\cdot \hat{{\boldsymbol{\imath }}}^{\prime} )\in [-\pi /2,\pi /2]$. We can concentrate on β ≥ 0 without any loss of generality because negative β corresponds to the mirroring across the xy plane ($z\mapsto -z$), cf. Appendix B.

We present our analysis on the rescaled Kopal potential $\widetilde{{\rm{\Omega }}}$ where $\rho ^{\prime} \mapsto \rho ^{\prime} \delta$ and ${\rm{\Omega }}\mapsto {\rm{\Omega }}\delta$, defined in the rotated coordinate system as

$$\begin{eqnarray}\widetilde{{\rm{\Omega }}}({\boldsymbol{\rho }}^{\prime} ;q,b,\beta ) & = & {\rm{\Omega }}({\boldsymbol{\rho }}^{\prime} \delta ;q,F,\delta ,\hat{{\boldsymbol{S}}})\delta \\ & = & \displaystyle \frac{1}{\rho ^{\prime} }+q\left(\displaystyle \frac{1}{\parallel {\boldsymbol{\rho }}^{\prime} -\hat{{\boldsymbol{\imath }}}^{\prime} \parallel }-\hat{{\boldsymbol{\imath }}}^{\prime} \cdot {\boldsymbol{\rho }}^{\prime} \right)\\ & & +\ \displaystyle \frac{1}{2}b[{(x^{\prime} \cos \beta -z^{\prime} \sin \beta )}^{2}+y{{\prime} }^{2}],\end{eqnarray} \tag{ 6 }$$

where b = (1 + q)F²δ³ is an auxiliary parameter. The shape of the star is fully determined by the value of $\widetilde{{\rm{\Omega }}}$: it corresponds to the isosurface of the potential. We refer to this shape as the lobe and denote it with $\widetilde{{ \mathcal L }}$. Figure 2 depicts the $\widetilde{{\rm{\Omega }}}$ contours that correspond to q = 1, b = 2 (F = 1, δ = 1) and β = 0.3π (left) and the lobe that corresponds to $\widetilde{{\rm{\Omega }}}=3.6$ (right).

Zoom In Zoom Out Reset image size

Equilibrium points of the Kopal potential Ω (Equation (3)) are defined as the points where the gradient of the potential equals 0. The simplified case of the aligned (β = 0), synchronous (F = 1) equipotential has been extensively studied. The corresponding equilibrium points are called Lagrange points. There are 5 such points: L₁, L₂ and L₃ lie on the x-axis, while L₄ and L₅ lie in the xy plane, forming an equilateral triangle with the two massive bodies. L₄ and L₅ can thus be computed analytically, while L₁, L₂, and L₃ are computed⁶ numerically. For a recent study of the analytic properties of Lagrange points, see Seidov (2004).

In the aligned and non-synchronous ($F\ne 1$) case there are still five equilibrium points ${{\boldsymbol{G}}}_{i}(q,b)$, i = 1, ..., 5 (Kallrath & Milone 2009), which are generalizations of the Lagrange points. In the F → 1 limit, these points are identical to Lagrange points. For the purposes of binary star physics, the first three generalized Lagrange points are our focus. Obtaining their values for arbitrary q and b using generic algorithms for solving nonlinear equations can be time-consuming and in some cases unstable. Because of this, we developed a specialized numerical solver to stabilize and speed up the process, which we have implemented into PHOEBE 2.0 (Prša et al. 2016). The solver is based on novel analytical approximations of generalized Lagrange points in difference regimes of parameters (q, b) that are further polished via the Newton–Raphson scheme or with the Laguerre method (Širca & Horvat 2012).

For the reduced potential $\widetilde{{\rm{\Omega }}}$ (Equation (6)), the condition ${\rm{\nabla }}\widetilde{{\rm{\Omega }}}({{\boldsymbol{K}}}_{i})=0$ that determines the equilibrium points ${{\boldsymbol{K}}}_{i}\,=(x^{\prime} ,y^{\prime} ,z^{\prime} )$ can be written as

$$\begin{eqnarray}\left(b\left[\begin{array}{ccc}{\cos }^{2}\beta & 0 & -\sin \beta \cos \beta \\ 0 & 1 & 0\\ -\sin \beta \cos \beta & 0 & {\sin }^{2}\beta \end{array}\right]\right.\\ &&\left.-\,\left(\displaystyle \frac{1}{{r}_{1}^{3}}+\displaystyle \frac{q}{{r}_{2}^{3}}\right)\mathrm{id}\right){{\boldsymbol{K}}}_{i}+q\left(\displaystyle \frac{1}{{r}_{2}^{3}}-1\right)\hat{{\boldsymbol{\imath }}}=0,\end{eqnarray} \tag{ 7 }$$

with b = (1 + q)F²δ³, ${r}_{1}\equiv \rho ^{\prime} =\parallel {{\boldsymbol{K}}}_{i}\parallel$ and ${r}_{2}=\parallel {{\boldsymbol{K}}}_{i}-\hat{{\boldsymbol{\imath }}}\parallel$. The equilibrium points are crucial for understanding the global behavior of the potential and for determining the necessary condition for lobe existence. The solutions of Equation (7) can be generally divided into two groups:

• (i)  
  The points outside the xz plane: For $b-q\,\geqslant {(\sqrt{2}| \cos \beta | )}^{-3}$, we find that there are two equilibrium points of the reduced Kopal potential $\widetilde{{\rm{\Omega }}}$:

  $$\begin{eqnarray}&&\left(\xi ,\pm \sqrt{2\xi -\displaystyle \frac{{\xi }^{2}}{{\cos }^{2}\beta }},-\xi \tan \beta \right),\end{eqnarray} \tag{ 8 }$$

  with $\xi =\tfrac{1}{2}{[b-q]}^{-\tfrac{2}{3}}$, which is a further generalization of the Lagrange points ${{\boldsymbol{G}}}_{\mathrm{4,5}}$. If, on the other hand, $b-q\lt {(\sqrt{2}| \cos \beta | )}^{-3}$, there are no real equilibrium points.
• (ii)  
  The points in the xz plane. In the aligned case (β = 0), there are exactly three equilibrium points ${{\boldsymbol{G}}}_{i}(q,b)$ for i = 1, 2, 3 and all are saddle points on the x-axis. In the misaligned case ($\beta \ne 0$), however, there can be more than three points and they can be of a different type (minimum, maximum, or a saddle point).

The points in the latter group determine the smallest value of the potential for which the lobe of the primary star exists. We discuss this group of equilibrium points next.

3.1. Phenomenology

Finding the equilibrium points of the reduced potential $\widetilde{{\rm{\Omega }}}$ (Equation (6)) in the xz plane is a non-trivial computational task. Figure 3 depicts the equilibrium points as a function of the misalignment parameter β, β ∈ [−π/2, π/2] and several values of q and b. The color of the points corresponds to the value of β. We see that the position of the equilibrium points varies continuously with β over a large range of values, but there are some discontinuities represented by the lack of points. Note the symmetry in the positions of points across the x-axis.

Zoom In Zoom Out Reset image size

We distinguish three types of equilibrium points in the xz plane: saddle points, local minima, and local maxima. Let $\{{h}_{k}({\boldsymbol{r}}):k\,=\,1,2\}$ be the eigenvalues of the Hessian matrix $({{\rm{\nabla }}}_{{xz}}\otimes {{\rm{\nabla }}}_{{xz}})\widetilde{{\rm{\Omega }}}({\boldsymbol{r}})$. These are proportional to the local principal curvatures, and the type of an equilibrium point is determined by the sign of the principal curvatures, $\mathrm{sign}({h}_{k}({{\boldsymbol{K}}}_{i}))$ (do Carmo 2016). The types of equilibrium points are important for drawing qualitative conclusions about the shape of the nearby isosurfaces. Figure 4 depicts the same equilibrium points as Figure 3, but here the colors denote their type. The saddle points are our focus, as they determine the separatrix, i.e. they yield the limiting value of the potential for which the lobes exist. If we consider equilibrium points of any given type as "branches," we see that, by perturbing a certain parameter, branches can cross, meaning that parts of the branches change their type.

Zoom In Zoom Out Reset image size

3.2. A Method for Finding Equilibrium Points

We developed and present here an efficient method to obtain a subset of equilibrium points in the xz plane. The location of the equilibrium points that correspond to the misaligned potential is found by tracing the variation in the x-location of the equilibrium points that correspond to the aligned potential as the misalignment parameter β is varied. This procedure is significantly faster than the general nonlinear root-finding algorithm employed in constructing Figures 3 and 4. The method is applicable for the values of β = 0 up to the value β_C,i for the ith equilibrium point at which the Hessian becomes singular.

Let us denote with ${ \mathcal I }$ an interval around the aligned value β = 0 at which the location of the equilibrium point ${{\boldsymbol{K}}}_{i}$ varies smoothly with β at any constant q and b. We can then write

$$\begin{eqnarray}&&{\rm{\nabla }}\widetilde{{\rm{\Omega }}}({{\boldsymbol{K}}}_{i}(\beta );q,b,\beta )=0\qquad \forall \beta \in { \mathcal I },\end{eqnarray} \tag{ 9 }$$

where ${{\boldsymbol{K}}}_{i}$ is the ith generalized Lagrange point at β = 0:

$$\begin{eqnarray*}&&{{\boldsymbol{K}}}_{i}(\beta =0)={{\boldsymbol{G}}}_{i}(q,b).\end{eqnarray*}$$

By differentiating Equation (9) w.r.t. β, we obtain a differential equation that determines the equilibrium point manifold as a function of β:

$$\begin{eqnarray}\displaystyle \frac{d}{d\beta }{{\boldsymbol{K}}}_{i}(\beta ) & = & -{[({\rm{\nabla }}\otimes {\rm{\nabla }})\widetilde{{\rm{\Omega }}}({{\boldsymbol{K}}}_{i}(\beta );q,b,\beta )]}^{-1}\\ & & \cdot \ {\rm{\nabla }}\displaystyle \frac{\partial }{\partial \beta }\widetilde{{\rm{\Omega }}}({{\boldsymbol{K}}}_{i}(\beta );q,b,\beta ).\end{eqnarray} \tag{ 10 }$$

To obtain the ith equilibrium point for a given (q, b, β₀), we integrate Equation (10) over the range [0, β₀] with the initial condition ${{\boldsymbol{K}}}_{i}(0)={{\boldsymbol{G}}}_{i}(q,b)$. This can be done numerically by an ordinary differential equation integrator, e.g., a 4th-order Runge–Kutta (Širca & Horvat 2012), as long as the Hessian (${\rm{\nabla }}\otimes {\rm{\nabla }})\widetilde{{\rm{\Omega }}}$ is non-singular. The Hessian is for all three Lagrange points non-singular on the entire β₀∈ [−π/2, π/2] range only for q = F = δ = 1; for other parameter values, β₀∈ [−β_C,i, β_C,i] for the ith equilibrium point. Typically, the smallest is β_C,3. Thus, a singular Hessian pinpoints the transition between the types of equilibrium points, as depicted in Figure 4.

3.3. Minimal Value of the Kopal Potential for the Existence of Lobes

With ${{\boldsymbol{K}}}_{i}(q,b;\beta )$ being the misaligned generalizations of their aligned counterparts G_i(q, b), the existence of a detached primary star lobe is determined by the values of ${{\boldsymbol{K}}}_{1}(q,b;\beta )$ and ${{\boldsymbol{K}}}_{2}(q,b;\beta )$. If the potential is smaller than either value, the equipotential will not delimit a closed surface. The minimal value of the reduced Kopal potential ${\widetilde{{\rm{\Omega }}}}_{\min }$ for which a detached primary lobe exists is thus the maximal value of the potentials at equilibrium points ${{\boldsymbol{K}}}_{1}$ and ${{\boldsymbol{K}}}_{2}$:

$$\begin{eqnarray}&&{\widetilde{{\rm{\Omega }}}}_{\min }(q,b,\beta )=\max \{\widetilde{{\rm{\Omega }}}({{\boldsymbol{K}}}_{1}(q,b;\beta )),\widetilde{{\rm{\Omega }}}({{\boldsymbol{K}}}_{2}(q,b;\beta ))\}.\end{eqnarray} \tag{ 11 }$$

The lobe at ${\widetilde{{\rm{\Omega }}}}_{\min }$(q, b, β) thus represents a generalized Roche lobe. Figure 5 demonstrates how the value of the potential changes with β for ${{\boldsymbol{K}}}_{1}(\beta ;q,b)$ and ${{\boldsymbol{K}}}_{2}(\beta ;q,b)$. As already pointed out by Avni & Schiller (1982), the curves of the potential values associated with both equilibrium points can intersect, meaning that the roles of the equilibrium points in constraining the lobes can change as β changes. The angle of intersection, referred to as the critical angle, was closely analyzed by Avni & Schiller. The minimum value of the potential for which a lobe exists at a certain angle β is the maximum value of both of these curves.

Zoom In Zoom Out Reset image size

The orbit of a binary system is described in the canonical coordinate system, where the orbital plane coincides with the xy plane. The plane-of-sky is the plane perpendicular to the line of sight. The corresponding coordinate system, spun by unit vectors $\hat{{\boldsymbol{u}}},\hat{{\boldsymbol{v}}}$, and $\hat{{\boldsymbol{w}}}$, is oriented so that $\hat{{\boldsymbol{u}}}$ and $\hat{{\boldsymbol{v}}}$ lie in the plane-of-sky and point toward east and north, respectively, and $\hat{{\boldsymbol{w}}}$ points toward the observer. To place the orbit in space w.r.t. the observer, we use three angles: longitude of the ascending node Ω, argument of periastron ω, and inclination i. The transformation from the orbital plane to the plane-of-sky is given by the following rotation:

$$\begin{eqnarray*}&&{{\boldsymbol{R}}}_{{uw}}(i,{\rm{\Omega }})={{\boldsymbol{R}}}_{w}({\rm{\Omega }}){{\boldsymbol{R}}}_{u}(-i),\end{eqnarray*}$$

where ${{\boldsymbol{R}}}_{w}$ and ${{\boldsymbol{R}}}_{u}$ are rotation matrices about the w and u axes, respectively. The definitions of rotation matrices here in use are given in Appendix E. The direction of the angular momentum $\hat{{\boldsymbol{L}}}$ is then given by

$$\begin{eqnarray*}&&\hat{{\boldsymbol{L}}}={{\boldsymbol{R}}}_{{uw}}(i,{\rm{\Omega }})\hat{{\boldsymbol{w}}}={[-\sin i\sin {\rm{\Omega }},\sin i\cos {\rm{\Omega }},\cos i]}^{T}.\end{eqnarray*}$$

Likewise, the spin vector $\hat{{\boldsymbol{S}}}$ can be written as

$$\begin{eqnarray*}&&\hat{{\boldsymbol{S}}}={{\boldsymbol{R}}}_{{uw}}(i^{\prime} ,{\rm{\Omega }}^{\prime} )\hat{{\bf{w}}}={[-\sin i^{\prime} \sin {\rm{\Omega }}^{\prime} ,\sin i^{\prime} \cos {\rm{\Omega }}^{\prime} ,\cos i^{\prime} ]}^{T},\end{eqnarray*}$$

where i' and Ω' denote the inclination and longitude of the ascending node w.r.t. the rotated coordinate system, and are related to the orbital inclination i and longitude of the ascending node Ω by

$$\begin{eqnarray*}&&{\rm{\Omega }}^{\prime} ={\rm{\Omega }}+{\rm{\Delta }}{\rm{\Omega }},\qquad i^{\prime} =i+{\rm{\Delta }}i.\end{eqnarray*}$$

The angle differences ΔΩ and Δi uniquely describe spin misalignment. The spin vector $\hat{{\boldsymbol{S}}}$ in the canonical coordinate system can then be written as

$$\begin{eqnarray*}&&\hat{{\boldsymbol{S}}}={{\boldsymbol{R}}}_{w}(-\varpi ){{\boldsymbol{R}}}_{u}(i){{\boldsymbol{R}}}_{w}({\rm{\Delta }}{\rm{\Omega }}){{\boldsymbol{R}}}_{u}(-i-{\rm{\Delta }}i)\hat{{\boldsymbol{w}}},\end{eqnarray*}$$

where ϖ = ω + ω₀ + ν is the argument o latitude of the considered star and is a sum of the true anomaly ν, the argument of periastron ω, and its positional offset ω₀, which depends on the star: ω₀ = π for the primary star and ω₀ = 0 for the secondary star.

4.1. Implementation in PHOEBE

PHOEBE is an open source modeling suite developed for the analysis of single, binary, and multiple stellar systems. Its initial version, released in 2005 and described in Prša & Zwitter (2005), was built on top of the widely used Wilson & Devinney (1971) code and it was specifically designed for the modeling of EB stars (hence the name, PHysics Of Eclipsing BinariEs). The updated version, PHOEBE 2.0, was released in 2016 and described by Prša et al. (2016). It constitutes a complete rewrite that generalizes the algorithms to single and multiple stellar systems. This work further expands the functionality of PHOEBE 2 by the implementation of misalignment described in the sections above and in the Appendices, and is accompanied by the PHOEBE 2.1 release. It is available in its entirety at http://phoebe-project.org.

PHOEBE 2.1 introduces spectral line profiles as a new type of data set. The line profiles are computed from fiducial spectral lines (i.e. a Gaussian or a Lorentzian profile) at the user-provided rest wavelength, Doppler-shifted at each local surface element and weighted by the passband brightness distribution across the visible surfaces. Line profiles are provided in normalized flux units and do not include any slopes due to continuum or passband effects.

Using PHOEBE, we demonstrate the effect of misalignment on astrophysical observables (light curves, radial velocity curves, and spectral line profiles) for a toy-model binary system with a misaligned secondary star. The parameters of the toy model are given in Table 1. Figure 6 showcases the comparison between the misaligned system and the aligned system with the matching equivalent radius (the radius of the sphere that has the same volume as the bounding equipotential) of each component. Spin misalignment clearly has a significant effect on all observables. A telltale sign of misalignment is an asymmetry in light curves, although asymmetries can arise from other physical effects as well, such as ellipsoidal variation and reflection in eccentric systems, spots, etc. Radial velocity curves are similar, with the telltale difference obvious in the eclipses (the RME). Figure 7 depicts line profiles for several phases in the aligned system (dashed line), the misaligned system (solid line), and the purely spherical system (dashed–dotted line) with matching equivalent radii. As expected, the differences are the largest during eclipses, as the main driver for the line profile is the sum of the local intensities weighted by the projected surface element area. Even outside the eclipses, though, the widths of the line profiles can be substantially different (top left panel) because of the modified surface brightness distribution across the disk of the secondary star. Other effects, such as relativistic gravitational redshift (Takeda & Ueno 2012), convective blueshift (Shporer & Brown 2011), microturbulence, and macroturbulence (Steffen et al. 2013), also affect the line profiles; while those can in principle be included in the computation within the PHOEBE framework, we did not include them in the simulation in order to quantify the influence of misalignment by itself.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Table 1.  Principal Parameters of the Toy Model of a Misaligned Binary Star

Parameter:		System:
	Primary Star:		Secondary Star:
Semimajor Axis a[R⊙]		3.98
Period P[days]		0.65
Mass Ratio q		0.7
Eccentricity		0
Inclination i [◦]		80
Long. of Ascending Node Ω[◦]		0
Systemic Velocity γ[km s−1]		0
Atmosphere	Blackbody		Blackbody
Equivalent Radius R[R⊙]	1.3		0.8
Effective Temperature Teff[K]	6500		5500
Synchronicity Parameter F	1		3.61
LD Model	Logarithmic		Logarithmic
LD Coefficient xLD	0.5		0.5
LD Coefficient yLD	0.5		0.5
Gravity Darkening βgrav	0.32		0.32
ΔΩ[◦]	0		45
Δi[◦]	0		45

Download table as:  ASCIITypeset image

Due to the selected misalignment parameters of the toy model, the corresponding light curve has a brighter out-of-eclipse region than the aligned model light curve. This is because the hotter polar regions are tilted toward the observer, thus contributing excess flux w.r.t. the aligned case. The generally asymmetric excess flux is more prominent in distorted (i.e., close) and/or rapidly rotating systems, where gravity darkening causes a significant variation of surface brightness across the stellar disk(s).

Radial velocity curves exhibit a well-known RME (McLaughlin 1924; Rossiter 1924), which describes a deviation from the dynamical (i.e., center of mass) radial velocity curve due to eclipses that block certain parts of the star and thus induce a bias in the photometrically weighted mean radial velocity curve for each component. The RME is symmetric for the aligned case, but is generally⁷ asymmetric in the misaligned case. This is depicted in Figure 6: the effect is symmetrical near timestamp 0, where the aligned primary star is eclipsed, and asymmetrical near timestamp 0.32, where the misaligned secondary star is eclipsed.

In consequence, light and radial velocity curves in conjunction allow us to solve for both misalignment parameters, Δi and ΔΩ. If spectral line profiles are also available, further improvement in the accuracy of these two parameters can generally be attained (Albrecht et al. 2007; Philippov & Rafikov 2013).

The treatment of misalignment in PHOEBE 2.1 is warranted whenever the tidal and rotational distortion of misaligned stars are non-negligible. Depending on the precision of acquired data points and the degree of misalignment, this detailed treatment may or may not be warranted and spherical models might be adequate in terms of precision and superior in terms of computation time. We are not aware of any other public codes that deal with misalignment in deformed stars. The computational time cost is only marginally impacted by the addition of misalignment.

4.1.1. The DI Herculis System

Albrecht et al. (2009) reported that DI Her is strongly misaligned, with the spin axes nearly perpendicular to the orbital axis. We use this example to further test and demonstrate the implementation of misaligned binary systems in PHOEBE. Figure 8 depicts the RV curves synthesized using the DI Her parameters summarized in Table 2 plotted over the observed radial velocities from Albrecht et al. (2009). The misalignment parameters are taken from Philippov & Rafikov (2013). We did not refit the data, as that is beyond the scope of the current paper; we only report qualitative agreement with the published results.

Zoom In Zoom Out Reset image size

Table 2.  Principal Parameters of the DI Herculis System

Parameter:		System:
	Primary Star:		Secondary Star:
Semimajor Axis a[R⊙]		42.8731
Period P[days]		10.550164
Mass Ratio q		0.815
Eccentricity		0.489
Inclination i[◦]		89.3
long. of ascending node Ω[◦]		330.2
Systemic Velocity γ[km s−1]		9.1
Time of Sup. Conjunction [days]		2442233.3481
Atmosphere	Blackbody		Blackbody
Equivalent Radius R[R⊙]	2.68		2.48
Effective Temperature Teff[K]	17300		15400
Synchronicity Parameter F	8.4819		9.8487
mass $m[{M}_{\odot }]$	5.1		4.4
LD Model	Logarithmic		Logarithmic
LD Coefficient xLD	0.5		0.5
LD Coefficient yLD	0.5		0.5
Gravity Darkening βgrav	1		1
${\rm{\Delta }}{\rm{\Omega }}{[}^{^\circ }]$	72		−84
${\rm{\Delta }}i{[}^{^\circ }]$	62		100

Download table as:  ASCIITypeset image

4.1.2. The Kepler-13Ab System

Kepler-13Ab is a transiting hot Jupiter system with an A-type host star. It was first discovered as a misaligned system by Szabó et al. (2011), and is the only such system ever found without the accompanying Rossiter–McLaughlin effect. Since then the system has been widely studied, yet the models feature inconsistent parameter values. For this paper, we take a representative sample of these values, given in Table 3, and create a model light curve of this system.

Table 3.  Principal Parameters of the Kepler-13 System

Parameter:		System:
	Kepler 13A :		Kepler 13b:
Semimajor Axis a[R⊙]		7.36a
Period $P[\mathrm{days}]$		1.763586522a
Mass Ratio q		0.0037
Eccentricity		0.0
Inclination i [◦]		85.82a
Atmosphere	Interpolated		Blackbody
Equivalent Radius $R[{R}_{\odot }]$	1.69a		0.144b
Effective Temperature Teff[K]	7650b		2750b
Synchronicity Parameter F	0.591c		1.0
Mass m[M⊙]	1.72b		0.0063b
LD Model	Interpolated		Logarithmic
LD Coefficient xLD	⋯		0.5
LD Coefficient yLD	⋯		0.5
Gravity Darkening βgrav	0.32		0.32
${\rm{\Delta }}{\rm{\Omega }}{[}^{^\circ }]$	4.52d		⋯
${\rm{\Delta }}i{[}^{^\circ }]$	58.6e		⋯

Notes.

^aMüller et al. (2013). ^bShporer et al. (2014). ^cSzabó et al. (2014). ^dMasuda (2015). ^eJohnson et al. (2014).

Download table as:  ASCIITypeset image

For comparison purposes we computed the light curves of the Kepler-13 system using different geometrical models: Roche, rotating star, and spherical star. We kept the volume and misalignment of each star fixed between these different models. The light curves are normalized w.r.t. their corresponding integrals in order to make them more comparable to each other. The results are depicted in Figure 9. The differences between light curves are the largest in the ingress and egress of the primary eclipse. In the middle panel we compare light curves computed by different distortion models with the light curve obtained by the spherical model while keeping the volume and the degree of misalignment constant. We see that the largest discrepancy is in the Roche model, of the order of ∼3 × 10⁻⁴, and it changes the signs depending on alignment. The bottom panel depicts the comparison of light curves of the misaligned model with the light curve of the aligned model, where we again see that the largest differences are in the Roche model, approximately equal to ∼6 × 10⁻⁴. Thus, the effect in Kepler-13 is under 1 mmag, however that is well within Kepler's precision reach of ∼20–30 ppm. Note that computing these differences accurately requires a sufficiently precise eclipsing algorithm, which is provided by PHOEBE.

Zoom In Zoom Out Reset image size

This paper summarizes the mathematical formalism of binary systems with the misaligned spin and orbital axes and introduces a new version of the modeling suite PHOEBE that implements this formalism. The topic has been studied in the past, i.e. by Avni & Schiller (1982), but to the best of our knowledge this is the first public implementation for the Roche-based geometry.

Beyond the anticipated systematic treatment of misalignment in EB and extrasolar planet systems, a thorough study of the parameter space can yield some very interesting and readily testable predictions. For example, by varying the misalignment parameter β, we can find a local minimum near the lobes as depicted in Figure 10. Such islands of stability could harbor Trojan objects that are synchronized with the rotation of the primary star and reminiscent of the features seen in Tabby's star, KIC 8462852 (Boyajian et al. 2016). Such hypotheses clearly merit further investigation beyond the scope of this introductory paper. The complexities of the nonlinear space spun by the parameters of misaligned objects predict many unexpected and intuition-challenging scenarios to exist in nature.

Zoom In Zoom Out Reset image size

We would like to thank Simon Albrecht for his attentive review of the manuscript in which he pointed out several deficiencies and suggested the discussion of photometric effects of misalignment in the case of Kepler-13. This work was supported by the NSF AAG grant #1517474, which we gratefully acknowledge. A.P. and K.H. also acknowledge partial funding from the Slovenian Research Agency Grant P1-0188. K.H. acknowledges the NASA ADAP grant 16-ADAP16-0201. K.C. is supported under NASA NESSF Fellowship #NNX15AR87H.

A binary system consists of two stars, labeled A and B. Their positions in the inertial (center of mass) coordinate system are denoted by ${{\boldsymbol{r}}}_{{\rm{A}}}$ and ${{\boldsymbol{r}}}_{{\rm{B}}}$. We assume that the center of mass of the binary system is at rest or moving with a constant velocity. Star A rotates as a rigid body about a misaligned axis $\hat{{\boldsymbol{S}}}$ with the angular velocity ${{\boldsymbol{\omega }}}_{S}$. The rigid body assumption asserts that every point on the primary star lobe co-rotates with the star. The equation of motion that describes the dynamics of the particle at position ${\boldsymbol{r}}$ is given by

$$\begin{eqnarray*}&&\ddot{{\boldsymbol{r}}}=-{\rm{\nabla }}U-\displaystyle \frac{1}{\rho }{\rm{\nabla }}p,\end{eqnarray*}$$

where p is the pressure, ρ is the particle density, and U is the gravitational potential of both stars:

$$\begin{eqnarray*}&&U=\displaystyle \frac{{{GM}}_{{\rm{A}}}}{\parallel {\boldsymbol{r}}-{{\boldsymbol{r}}}_{{\rm{A}}}\parallel }+\displaystyle \frac{{{GM}}_{{\rm{B}}}}{\parallel {\boldsymbol{r}}-{{\boldsymbol{r}}}_{{\rm{B}}}\parallel }.\end{eqnarray*}$$

We now introduce a canonical coordinate system that is centered in star A: its x-axis points toward star B, its z-axis is aligned with the revolution axis ${\boldsymbol{L}}$, and it co-rotates with the center of star A in orbit about the common center of mass. We express vector ${\boldsymbol{r}}$ as the sum of the vector to the center of star A (${{\boldsymbol{r}}}_{{\rm{A}}}$) and the vector relative to the center of star A (${\boldsymbol{r}}^{\prime}$):

$$\begin{eqnarray*}&&{\boldsymbol{r}}={{\boldsymbol{r}}}_{{\rm{A}}}+{\boldsymbol{r}}^{\prime} ,\qquad \ddot{{\boldsymbol{r}}}={\ddot{{\boldsymbol{r}}}}_{{\rm{A}}}+{\ddot{{\boldsymbol{r}}}}^{{\prime} }.\end{eqnarray*}$$

The term ${\ddot{{\boldsymbol{r}}}}_{{\rm{A}}}$ describes the acceleration of star A caused by gravity:

$$\begin{eqnarray*}&&{\ddot{{\boldsymbol{r}}}}_{{\rm{A}}}={{GM}}_{{\rm{B}}}\displaystyle \frac{{{\boldsymbol{r}}}_{{\rm{B}}}-{{\boldsymbol{r}}}_{{\rm{A}}}}{{\parallel {{\boldsymbol{r}}}_{{\rm{B}}}-{{\boldsymbol{r}}}_{{\rm{A}}}\parallel }^{3}}.\end{eqnarray*}$$

The term $\ddot{{\boldsymbol{r}}}^{\prime}$ corresponds to the acceleration relative to the center of the primary star. To express it, we introduce a third coordinate system S that co-rotates with the primary star itself about the rotation axis $\hat{{\boldsymbol{S}}}$. The relative vector $\ddot{{\boldsymbol{r}}}^{\prime}$ is then

$$\begin{eqnarray}&&\ddot{{\boldsymbol{r}}}^{\prime} ={(\ddot{{\boldsymbol{r}}}^{\prime} )}_{S}+2{{\boldsymbol{\omega }}}_{S}\times {(\dot{{\boldsymbol{r}}}^{\prime} )}_{S}+{\dot{{\boldsymbol{\omega }}}}_{S}\times {\boldsymbol{r}}^{\prime} +{{\boldsymbol{\omega }}}_{S}\times {{\boldsymbol{\omega }}}_{S}\times {\boldsymbol{r}}^{\prime} .\end{eqnarray} \tag{ 12 }$$

We assume that the angular velocity is constant in time, so ${\dot{{\boldsymbol{\omega }}}}_{S}=0$, and both the velocity and the acceleration in the co-rotating frame S are zero: ${(\ddot{{\boldsymbol{r}}}^{\prime} )}_{S}=0$, ${(\dot{{\boldsymbol{r}}}^{\prime} )}_{S}=0$. The equation of motion thus takes the following form:

$$\begin{eqnarray}&&{{GM}}_{{\rm{B}}}\displaystyle \frac{{{\boldsymbol{r}}}_{{\rm{B}}}-{{\boldsymbol{r}}}_{{\rm{A}}}}{\parallel {{\boldsymbol{r}}}_{{\rm{B}}}-{{\boldsymbol{r}}}_{{\rm{A}}}{\parallel }^{3}}+{{\boldsymbol{\omega }}}_{S}\times {{\boldsymbol{\omega }}}_{S}\times {\boldsymbol{r}}^{\prime} =-{\rm{\nabla }}U-\displaystyle \frac{1}{\rho }{\rm{\nabla }}p.\end{eqnarray} \tag{ 13 }$$

The first term in Equation (13) can be written as the radial gradient ${{\rm{\nabla }}}_{{\boldsymbol{r}}^{\prime} }$ of the potential, where ${{\rm{\nabla }}}_{{\boldsymbol{r}}^{\prime} }\equiv d/d{\boldsymbol{r}}^{\prime}$ operates in the canonical coordinate system spun by the basis vectors $(\hat{{\boldsymbol{\imath }}},\hat{{\boldsymbol{\jmath}}},\hat{{\boldsymbol{k}}})$:

$$\begin{eqnarray*}&&{{GM}}_{{\rm{B}}}\displaystyle \frac{{{\boldsymbol{r}}}_{{\rm{B}}}-{{\boldsymbol{r}}}_{{\rm{A}}}}{\parallel {{\boldsymbol{r}}}_{{\rm{B}}}-{{\boldsymbol{r}}}_{{\rm{A}}}{\parallel }^{3}}\equiv {{GM}}_{B}\displaystyle \frac{\hat{{\boldsymbol{\imath }}}}{{d}^{2}}={{\rm{\nabla }}}_{{\boldsymbol{r}}^{\prime} }\left\{{{GM}}_{{\rm{B}}}\displaystyle \frac{\hat{{\boldsymbol{\imath }}}\cdot {\boldsymbol{r}}^{\prime} }{{d}^{2}}\right\},\end{eqnarray*}$$

with d being the distance between the stars, and similarly, using the triple product rule ${\boldsymbol{a}}\times {\boldsymbol{b}}\times {\boldsymbol{c}}={\boldsymbol{b}}({\boldsymbol{a}}\cdot {\boldsymbol{c}})-{\boldsymbol{c}}({\boldsymbol{a}}\cdot {\boldsymbol{b}})$, the last term in Equation (12) can be written as:

$$\begin{eqnarray*}{{\boldsymbol{\omega }}}_{S}\times {{\boldsymbol{\omega }}}_{S}\times {\boldsymbol{r}}^{\prime} & \equiv & {\omega }_{S}^{2}[(\hat{{\boldsymbol{S}}}\cdot {\boldsymbol{r}}^{\prime} )\hat{{\boldsymbol{S}}}-{\boldsymbol{r}}^{\prime} ]\\ & = & {{\rm{\nabla }}}_{{\boldsymbol{r}}^{\prime} }\left\{\displaystyle \frac{1}{2}{\omega }_{S}^{2}{[(\hat{{\boldsymbol{S}}}\cdot {\boldsymbol{r}}^{\prime} )\hat{{\boldsymbol{S}}}-{\boldsymbol{r}}^{\prime} ]}^{2}\right\}.\end{eqnarray*}$$

Using these two expressions in the equation of motion yields

$$\begin{eqnarray*}&&{{\rm{\nabla }}}_{{\boldsymbol{r}}^{\prime} }\left\{U+{{GM}}_{{\rm{B}}}\displaystyle \frac{\hat{{\boldsymbol{\imath }}}\cdot {\boldsymbol{r}}^{\prime} }{{d}^{2}}+\displaystyle \frac{1}{2}{[(\hat{{\boldsymbol{S}}}\cdot {\boldsymbol{r}}^{\prime} )\hat{{\boldsymbol{S}}}-{\boldsymbol{r}}^{\prime} ]}^{2}\right\}\\ &&\quad =\,-\displaystyle \frac{1}{\rho }{{\rm{\nabla }}}_{{\boldsymbol{r}}^{\prime} }p.\end{eqnarray*}$$

We can now readily recognize the expression within the curly braces as the (negative) potential V of the misaligned binary system, as given in Equation (1).

The reduced potential has several useful symmetries that have been implicitly used in the paper:

$$\begin{eqnarray*}\widetilde{{\rm{\Omega }}}(x^{\prime} ,y^{\prime} ,z^{\prime} ;q,b,\beta ) & = & \widetilde{{\rm{\Omega }}}(x^{\prime} ,-y^{\prime} ,z^{\prime} ;q,b,\beta ),\\ \widetilde{{\rm{\Omega }}}(x^{\prime} ,y^{\prime} ,z^{\prime} ;q,b,\beta ) & = & \widetilde{{\rm{\Omega }}}(x^{\prime} ,y^{\prime} ,z^{\prime} ;q,b,\beta +\pi ),\\ \widetilde{{\rm{\Omega }}}(x^{\prime} ,y^{\prime} ,z^{\prime} ;q,b,\beta ) & = & \widetilde{{\rm{\Omega }}}(x^{\prime} ,y^{\prime} ,-z^{\prime} ;q,b,-\beta ).\end{eqnarray*}$$

The poles are defined as the radii of the lobes $\widetilde{{ \mathcal L }}$ along the positive and negative direction of the spin vector $\hat{{\boldsymbol{S}}}$. In order to investigate this in more detail, we further rotate a coordinate system with the vector basis $(\hat{{\boldsymbol{\imath }}}^{\prime} ,\hat{{\boldsymbol{\jmath}}}^{\prime} ,\hat{{\boldsymbol{k}}}^{\prime} )$, given by Equation (4), so that the new z-axis is aligned with the spin vector. The vector basis of this new coordinate system is:

$$\begin{eqnarray}\hat{{\boldsymbol{\imath }}}^{\prime\prime} & = & \hat{{\boldsymbol{\imath }}}^{\prime} \cos \beta -\hat{{\boldsymbol{k}}}^{\prime} \sin \beta ,\\ \hat{{\boldsymbol{\jmath}}}^{\prime\prime} & = & \hat{{\boldsymbol{\jmath}}}^{\prime} ,\\ \hat{{\boldsymbol{k}}}^{\prime\prime} & = & \hat{{\boldsymbol{\imath }}}^{\prime} \sin \beta +\hat{{\boldsymbol{k}}}^{\prime} \cos \beta ,\end{eqnarray} \tag{ 14 }$$

and position denoted by ${\boldsymbol{\rho }}^{\prime\prime} =r(\sin \theta \cos \phi \,\hat{{\boldsymbol{\imath }}}^{\prime\prime} +\sin \theta \sin \phi \,\hat{{\boldsymbol{\jmath}}}^{\prime\prime} \,+\cos \theta \,\hat{{\boldsymbol{k}}}^{\prime\prime} )$. In this coordinate system, the reduced Kopal potential can be rewritten as

$$\begin{eqnarray}&&W(r,\theta ,\phi ;q,b,\beta )=\widetilde{{\rm{\Omega }}}({\boldsymbol{\rho }}^{\prime\prime} ;q,b,\beta ),\\ &&=\ \displaystyle \frac{1}{r}+\displaystyle \frac{1}{2}{{br}}^{2}{\sin }^{2}\theta -{qr}(\sin \beta \cos \theta +\cos \beta \sin \theta \cos \phi )\\ &&\quad +\ \displaystyle \frac{q}{\sqrt{1-2r(\sin \beta \cos \theta +\cos \beta \sin \theta \cos \phi )+{r}^{2}}}.\end{eqnarray} \tag{ 15 }$$

We see that W does not have a quadratic term for the distance from the origin, which is associated with the centrifugal contribution to the potential. In general, the lobe is not symmetric across z = 0 for the new coordinate frame, so the poles in the positive and negative directions of the rotating axis are not equal.

For a given set of parameters (q, b, β) and the reference potential value ${\widetilde{{\rm{\Omega }}}}_{0}$, the pole in the positive direction of the spin, r₊, and in the negative direction, r₋, are defined as

$$\begin{eqnarray}&&W({r}_{+},0,0;q,b,\beta )=W({r}_{-},\pi ,0;q,b,\beta )={\widetilde{{\rm{\Omega }}}}_{0}.\end{eqnarray} \tag{ 16 }$$

This yields the following equation for the poles:

$$\begin{eqnarray}&&\displaystyle \frac{1}{{r}_{\pm }}+q\left(\displaystyle \frac{1}{\sqrt{1\mp 2\sin \beta {r}_{\pm }+{r}_{\pm }^{2}}}\mp \sin \beta {r}_{\pm }\right)={\widetilde{{\rm{\Omega }}}}_{0}.\end{eqnarray} \tag{ 17 }$$

In the case of lobes with spin–orbit misalignment, it is more meaningful to discuss the diameter r₊ + r₋ along the rotation axis as a measure of the characteristic size of the object. These equations are solved in PHOEBE by employing the standard Newton–Raphson method.

In the limit of large potential reference values, ${\widetilde{{\rm{\Omega }}}}_{0}\gg 1$, the poles can be approximated by a power series in $s={\widetilde{{\rm{\Omega }}}}_{0}^{-1}$, obtained by the inverse series method. The expansion of poles ${r}_{\pm }$ in the positive (+) and negative (−) direction are identical up to the 4th degree in s:

$$\begin{eqnarray*}{r}_{\pm } & = & s+{{qs}}^{2}+{q}^{2}{s}^{3}\\ & & +\ \displaystyle \frac{1}{2}{{qs}}^{4}(2{q}^{2}+3{\sin }^{2}\beta -1)+O({s}^{5}),\end{eqnarray*}$$

with the difference between the poles found only in terms of degree 5 and higher:

$$\begin{eqnarray*}&&{r}_{+}-{r}_{-}=q(-3+5{\sin }^{2}\beta ){s}^{5}+O({s}^{6}).\end{eqnarray*}$$

We see that, in the limit of a large potential, the lobe size depends only weakly on the misalignment parameter.

In general, the pole ${r}_{\pm }$ can be obtain by integrating the differential of the pole w.r.t. the reciprocal potential s, given by

$$\begin{eqnarray*}&&\displaystyle \frac{{{dr}}_{\pm }}{{ds}}=\displaystyle \frac{g{(g+{{qr}}_{\pm }\mp {qg}\sin \beta {r}_{\pm }^{2})}^{2}}{{g}^{3}\pm q({g}^{3}-1)\sin \beta {r}_{\pm }^{2}+{{qr}}_{\pm }^{3}},\end{eqnarray*}$$

for $s\in [0,{\widetilde{{\rm{\Omega }}}}_{0}^{-1}]$ and the initial condition r_±  (s = 0) = 1 and $g=\sqrt{1\mp 2\sin \beta {r}_{\pm }+{r}_{\pm }^{2}}$.

We now turn our attention to the surface area ($\widetilde{A}$) and volume ($\widetilde{V}$) of the lobes $\widetilde{{ \mathcal L }}$ and the derivative of the volume w.r.t. the value of the potential (${\widetilde{V}}_{,{\widetilde{{\rm{\Omega }}}}_{0}})$. We present a numerical method to compute these quantities using spherical coordinates (r, θ, ϕ) and the reduced Kopal potential W (Equation (15)). We write a partial derivative of a function f w.r.t. variable x as f_,x = ∂f/∂x.

If r(θ, ϕ) is known, the quantities are given by the following integrals:

$$\begin{eqnarray}\widetilde{A} & = & {\displaystyle \int }_{0}^{\pi }\dot{\widetilde{A}}(\theta )\,d\theta \\ & & \dot{\widetilde{A}}(\theta )=\displaystyle \frac{2}{3}{\displaystyle \int }_{0}^{\pi }r(\theta ,\phi )\sqrt{{r}_{,\phi }^{2}+{\sin }^{2}\theta ({r}^{2}+{r}_{,\theta }^{2})}\,d\phi ,\end{eqnarray} \tag{ 18 }$$

$$\begin{eqnarray}\widetilde{V} & = & {\displaystyle \int }_{0}^{\pi }\dot{\widetilde{V}}(\theta )\,d\theta \\ & & \dot{\widetilde{V}}(\theta )=\displaystyle \frac{2}{3}\sin \theta {\displaystyle \int }_{0}^{\pi }{r}^{3}(\theta ,\phi )\,d\phi ,\end{eqnarray} \tag{ 19 }$$

$$\begin{eqnarray}{\widetilde{V}}_{,{\widetilde{{\rm{\Omega }}}}_{0}} & = & {\displaystyle \int }_{0}^{\pi }{\dot{\widetilde{V}}}_{,{\widetilde{{\rm{\Omega }}}}_{0}}(\theta )\,d\theta \\ & & {\dot{\widetilde{V}}}_{,{\widetilde{{\rm{\Omega }}}}_{0}}(\theta )=2{\displaystyle \int }_{0}^{\pi }\displaystyle \frac{{r}^{2}(\theta ,\phi )}{{P}_{,r}(r(\theta ,\phi ),\theta ,\phi )}\,d\phi ,\end{eqnarray} \tag{ 20 }$$

where we took into account the symmetry over the xz plane. The derivatives of the radius r w.r.t. spherical angles (θ, ϕ) are given by:

$$\begin{eqnarray}{r}_{,\theta }(r,\theta ,\phi ) & = & -\displaystyle \frac{{W}_{,\theta }(r,\theta ,\phi )}{{P}_{,r}(r,\theta ,\phi )},\\ {r}_{,\phi }(r,\theta ,\phi ) & = & -\displaystyle \frac{{W}_{,\phi }(r,\theta ,\phi )}{{P}_{,r}(r,\theta ,\phi )}.\end{eqnarray} \tag{ 21 }$$

The derivative ${\widetilde{V}}_{,{\widetilde{{\rm{\Omega }}}}_{0}}$ is needed in the volume conservation process, whereby we find the value of the potential ${\widetilde{{\rm{\Omega }}}}_{0}$ corresponding to a certain volume ${\widetilde{V}}_{0}$ as other parameters are fixed. This is analogous to calculating the inverse of ${\widetilde{{\rm{\Omega }}}}_{0}={\widetilde{V}}^{-1}({\widetilde{V}}_{0})$ by the Newton–Raphson method:

$$\begin{eqnarray}&&{\widetilde{{\rm{\Omega }}}}_{0,k+1}={\widetilde{{\rm{\Omega }}}}_{0,k}-\displaystyle \frac{\widetilde{V}({\widetilde{{\rm{\Omega }}}}_{0,k})-{\widetilde{V}}_{0}}{{V}_{,{\widetilde{{\rm{\Omega }}}}_{0}}({\widetilde{{\rm{\Omega }}}}_{0,k})}.\end{eqnarray} \tag{ 22 }$$

We perform the calculation of $\widetilde{A}$, $\widetilde{V}$, and ${\widetilde{V}}_{,{\widetilde{{\rm{\Omega }}}}_{0}}$ using two techniques: the integration across the surface and the asymptotic approximation in the limit of small lobes (large values of the potential). We explain both methods below.

D.1. Integration Over the Surface

The quantities $\widetilde{A}$, $\widetilde{V}$, and ${\widetilde{V}}_{,{\widetilde{{\rm{\Omega }}}}_{0}}$ are written as definite integrals of their derivatives $\dot{\widetilde{A}}$, $\dot{\widetilde{V}}$, and ${\dot{\widetilde{V}}}_{,{{\rm{\Omega }}}_{0}}$ over the azimutal angle θ ∈ [0, π] per Equations (18)–(20). The derivatives are given as integrals over the polar angles ϕ ∈ [0, π]. We start the calculation by first approximating the derivatives and then integrate them across the total range of azimutal angle.

The integrals defining derivatives are calculated by discretizing the polar angle domain and using the Legendre–Gauss quadrature (Širca & Horvat 2012), whereby an integral of a function g over the interval [0, π] is approximated by:

$$\begin{eqnarray}&&{\int }_{0}^{\pi }g(\phi )\,d\phi =\displaystyle \sum _{k=1}^{n}{u}_{i}g({\phi }_{i})+\displaystyle \frac{{\pi }^{2n+1}{(n!)}^{4}}{(2n+1){[(2n)!]}^{3}}{g}^{(2n)}(\zeta ),\end{eqnarray} \tag{ 23 }$$

where u_i and ϕ_i are appropriately chosen weights and nodes, respectively, and ζ ∈ [0, π]. The weights and nodes are given by

$$\begin{eqnarray*}&&{u}_{i}=\displaystyle \frac{\pi }{2}{w}_{i},\quad {\phi }_{i}=\displaystyle \frac{\pi }{2}(1+{x}_{i}),\end{eqnarray*}$$

where w_i and x_i are standard Legendre–Gauss weights and nodes, respectively, determined for functions integrated over the range [−1, 1].

The radius r(θ, ϕ) of the lobe at arbitrary angles θ and ϕ can be obtained by integrating

$$\begin{eqnarray}&&\displaystyle \frac{{dr}}{d\theta }\equiv {r}_{,\theta }(r,\theta ,{\phi }_{i}),\end{eqnarray} \tag{ 24 }$$

with the initial condition r(θ = 0, ϕ) = r₊, where the derivative r_{, θ} is given by Equation (21). For each polar angle ϕ_i we introduce a radius r_i(θ) = r(θ, ϕ_i) along the azimuthal angle θ. Using Legendre–Gauss quadrature (Equation (23)) we approximate the integrals defining derivatives $\dot{\widetilde{A}}$, $\dot{\widetilde{V}}$, and ${\dot{\widetilde{V}}}_{,{\widetilde{{\rm{\Omega }}}}_{0}}$ as sums over the set of functions ${\{{r}_{i}\}}_{i=1}^{n}$. Then, by taking into account Equation (24), we rewrite $\widetilde{A}$, $\widetilde{V}$, and ${\widetilde{V}}_{,{\widetilde{{\rm{\Omega }}}}_{0}}$ as a solution to n + 3 ordinary differential equations:

$$\begin{eqnarray}&&\displaystyle \frac{{{dr}}_{i}}{d\theta }={r}_{,\theta }({r}_{i}(\theta ),\theta ,{\phi }_{i})\quad i=1,\,\ldots ,\,n,\end{eqnarray} \tag{ 25 }$$

$$\begin{eqnarray}&&\displaystyle \frac{d\widetilde{A}}{d\theta }=2\displaystyle \sum _{i=1}^{n}{u}_{i}{r}_{i}(\theta )\\ &&\times \,\sqrt{{r}_{,\phi }^{2}({r}_{i}(\theta ),\theta ,{\phi }_{i})+{\sin }^{2}\theta ({r}_{i}^{2}(\theta )+{r}_{,\theta }^{2}({r}_{i}(\theta ),\theta ,{\phi }_{i}))},\end{eqnarray} \tag{ 26 }$$

$$\begin{eqnarray}&&\displaystyle \frac{d\widetilde{V}}{d\theta }=\displaystyle \frac{2}{3}\sin \theta \displaystyle \sum _{i=1}^{n}{u}_{i}{r}_{i}^{3}(\theta ),\end{eqnarray} \tag{ 27 }$$

$$\begin{eqnarray}&&\displaystyle \frac{d{\widetilde{V}}_{,{\widetilde{{\rm{\Omega }}}}_{0}}}{d\theta }=2\displaystyle \sum _{i=1}^{n}{u}_{i}\displaystyle \frac{{r}_{i}^{2}(\theta )}{{W}_{,r}({r}_{i}(\theta ),\theta ,{\phi }_{i})},\end{eqnarray} \tag{ 28 }$$

which are integrated for θ ∈ [0, π] with the initial conditions

$$\begin{eqnarray}&&{r}_{i}(\theta =0)={r}_{+}\quad i=1,\,\ldots ,\,n,\end{eqnarray} \tag{ 29 }$$

$$\begin{eqnarray}&&\widetilde{A}(\theta =0)=\widetilde{V}(\theta =0)={\widetilde{V}}_{,{\widetilde{{\rm{\Omega }}}}_{0}}(\theta =0)=0.\end{eqnarray} \tag{ 30 }$$

The quantities in question are then obtained at θ = π:

$$\begin{eqnarray*}\widetilde{A} & = & \widetilde{A}(\theta =\pi ),\\ \widetilde{V} & = & \widetilde{V}(\theta =\pi ),\\ {\widetilde{V}}_{,{\widetilde{{\rm{\Omega }}}}_{0}} & = & {\widetilde{V}}_{,{\widetilde{{\rm{\Omega }}}}_{0}}(\theta =\pi ).\end{eqnarray*}$$

Note that r_i(θ = π) = r₋, which can be used as a numerical check of integration.

D.2. The Limit of Small Lobes

In the limit of large ${\widetilde{{\rm{\Omega }}}}_{0}$, the radius can be expressed as a power series of $s={{\widetilde{{\rm{\Omega }}}}_{0}}^{-1}$ using the inverse series technique, which we can symbolically write as

$$\begin{eqnarray}&&r(\theta ,\phi ;q,b,s,\beta )=\displaystyle \sum _{k=1}^{\infty }{a}_{k}(\theta ,\phi ;q,b,\beta ){s}^{k}.\end{eqnarray} \tag{ 31 }$$

By plugging Equation (31) into the formulæ for area $\widetilde{A}$ (Equation (18)) and volume $\widetilde{V}$ (Equation (19)), we obtain their own expansions in s and write them as

$$\begin{eqnarray}&&\widetilde{A}=\displaystyle \frac{4\pi }{{{\widetilde{{\rm{\Omega }}}}_{0}}^{2}}{ \mathcal A }(s;q,b,\beta )\qquad \widetilde{V}=\displaystyle \frac{4\pi }{3{\widetilde{{\rm{\Omega }}}}_{0}^{3}}{ \mathcal V }(s;q,b,\beta ),\end{eqnarray} \tag{ 32 }$$

with auxiliary expressions

$$\begin{eqnarray}&&{ \mathcal A }(s;q,b,\beta )=1+2{qs}+3{q}^{2}{s}^{2}+\ {s}^{3}\left(\displaystyle \frac{2b}{3}+4{q}^{3}\right)\\ &&\,\,\,\,+\ 5{s}^{4}\left(\displaystyle \frac{2{bq}}{3}+{q}^{4}\right)+{s}^{5}(10{{bq}}^{2}+6{q}^{5})\\ &&\,\,\,\,+\ {s}^{6}\left({b}^{2}+\displaystyle \frac{70{{bq}}^{3}-{bq}}{3}+{bq}\cos (2\beta )\right.\\ &&\,\,\,\,+\,7{q}^{6}+2{q}^{2})\,+\ldots ,\end{eqnarray} \tag{ 33 }$$

$$\begin{eqnarray}{ \mathcal V }(s;q,b,\beta ) & = & 1+3{qs}+6{q}^{2}{s}^{2}+{s}^{3}(b+10{q}^{3})\\ & & +{s}^{4}(6{bq}+15{q}^{4})+21{s}^{5}({{bq}}^{2}+{q}^{5})\\ & & +\ \displaystyle \frac{2}{5}{s}^{6}(4{b}^{2}+140{{bq}}^{3}+3{bq}\cos (2\beta )\\ & & -\,{bq}+70{q}^{6}+6{q}^{2})+\ldots .\end{eqnarray} \tag{ 34 }$$

The influence of misalignment is thus very weak in the limit of small lobes and it affects only the terms of the 6th degree and higher in the series expansion.

The expression for ${\widetilde{V}}_{,{\widetilde{{\rm{\Omega }}}}_{0}}$, needed in the volume conservation procedure defined by Equation (22), is obtained by taking the derivative of the volume $\widetilde{V}$ (Equation (32)) w.r.t. the reference potential value ${\widetilde{{\rm{\Omega }}}}_{0}$:

$$\begin{eqnarray}&&{\widetilde{V}}_{,{\widetilde{{\rm{\Omega }}}}_{0}}=-\displaystyle \frac{4\pi }{3{{\widetilde{{\rm{\Omega }}}}_{0}}^{4}}\left(3{ \mathcal V }(s;q,b,\beta )+s\displaystyle \frac{d{ \mathcal V }(s;q,b,\beta )}{{ds}}\right).\end{eqnarray} \tag{ 35 }$$

In the paper we use the following convention for the rotation matrices in three dimensions with rotation angle ϕ:

$$\begin{eqnarray*}{{\boldsymbol{R}}}_{u}(\phi )=\left[\begin{array}{ccc}1 & 0 & 0\\ 0 & \cos \phi & -\sin \phi \\ 0 & \sin \phi & \cos \phi \end{array}\right]\\ {{\boldsymbol{R}}}_{v}(\phi )=\left[\begin{array}{ccc}\cos \phi & 0 & \sin \phi \\ 0 & 1 & 0\\ -\sin \phi & 0 & \cos \phi \end{array}\right]\\ {{\boldsymbol{R}}}_{w}(\phi )=\left[\begin{array}{ccc}\cos \phi & -\sin \phi & 0\\ \sin \phi & \cos \phi & 0\\ 0 & 0 & 1\end{array}\right].\end{eqnarray*}$$