Interacting Young M-dwarfs in Triple System -- Par 1802 Binary System Case Study

TL;DR

This study investigates how a third companion in the Par 1802 triple system can induce tidal heating and mass transfer, explaining the observed temperature differences between the binary stars through detailed dynamical and evolutionary modeling.

Contribution

The paper presents a novel dynamical and evolutionary model showing that mass transfer and tidal heating can account for temperature discrepancies in young binary stars influenced by a tertiary companion.

Findings

Mass transfer and tidal heating explain temperature differences.

Predicted orbital configuration of the tertiary companion.

Dynamical evolution impacts pre-main-sequence stellar properties.

Abstract

The binary star Par 1802 in the Orion Nebula presents an interesting puzzle in the field of stellar dynamics and evolution. Binary systems such as Par 1802 are thought to form from the same natal material and thus the stellar members are expected to have very similar physical attributes. However, Par 1802's stars have significantly different temperatures despite their identical (within $3\%$) masses of about $0.39$ solar mass. The leading proof-of-concept idea is that a third companion gravitationally induced the two stars to orbit closer than their Roche-limit, which facilitated heating through tidal effects. Here we expand on this idea and study the three-body dynamical evolution of such a system, including tidal and pre-main-sequence evolution. We also include tidal heating and mass transfer at the onset of Roche-limit Crossing. We show, as a proof-of-concept, that mass transfer combined with tidal heating can naturally explain the observed temperature discrepancy. We also predict the orbital configuration of the possible tertiary companion. Finally we suggest that the dynamical evolution of such a system has pervasive consequences. We expect an abundance of systems to undergo mass transfer during their pre-main-sequence time, which can cause temperature differences.

Equations

$$R_{\mathrm{Roche},1}=2.7\times R_1{\left(\frac{m_1+m_2}{m_1}\right)}^{-1/3}.$$

$$a_i(1-e_i^2)\approx a_f,$$

$$a_f\sim 2a_i(1-e_i)\sim 2R_{\mathrm{Roche}},$$

$$R_1(t)\sim \frac{a_f}{2\times 2.7}{\left(\frac{m_1+m_2}{m_1}\right)}^{-1/3}\sim 1.57R_{\odot }.$$

$$t_k\sim \frac{16}{15}\frac{\sqrt{m_1+m_2}}{\sqrt{G}{m}_3}\frac{a_2^3}{a_1^{3/2}}(1-e_2^2{)}^{3/2},$$

$$ϵ=\frac{a_1}{a_2}\frac{e_2}{1-e_2^2}<0.1.$$

$$\frac{a_2}{a_1}>2.8{\left(1+\frac{m_2}{m_1+m_2}\right)}^{2/5}\frac{(1+e_2{)}^{2/5}}{(1-e_2{)}^{6/5}}\left(1-\frac{0.3i}{180^{\circ }}\right).$$

$$a_1(1-e_1)>R_{\mathrm{Roche},1}.$$

$$m_1+m_2=0.776M_{\odot }=\text{constant},$$

$$\dot{M}=\frac{\Delta M}{\Delta t},$$

$$Q_j=\frac{4}{3}\frac{k_1}{(1+2k_1{)}^2}\frac{G{m}_j}{R_j^3}\frac{t_V}{n_1},$$

$$Q_j{ϵ}_{0,j}$$

$$Q_j{ϵ}_{1,j}$$

$$Q_j{ϵ}_{2,j}$$

$$Q_j{ϵ}_{5,j}$$

$$Q_j{ϵ}_{8,j}$$

$$Q_j{ϵ}_{9,j}$$

$$\Sigma (x)=+1\vee -1,$$

$${\dot{E}}_{\text{orbit},j}$$

$${\dot{E}}_{\text{rotation},j}$$

$${\dot{E}}_j=-({\dot{E}}_{\text{orbit},j}+{\dot{E}}_{\text{rotation},j}).$$

$$d{T}_j={\left(\frac{{\dot{E}}_j/2}{4\pi R_j^2{\sigma }_{\text{SB}}}+T_{\text{eff},j}^4\right)}^{1/4}-T_{\text{eff},j},$$

$$\Delta T=(d{T}_1+T_{\text{eff},1})-(d{T}_2+T_{\text{eff},2}).$$

Full text

††pubyear: 2019††pagerange: Interacting Young M-dwarfs in Triple System - Par 1802 Binary System Case Study–References

1 Introduction

Basic physical properties of a star such as mass and temperature are strongly determined by its initial mass and chemical composition. Equal-mass components of binary stars are typically considered to be coeval and identical ‘twins’ formed from the same protostellar core, and are expected to have identical physical attributes. However, a small number of identical twin systems seem to be inconsistent with this narrative, featuring a large temperature difference between its component stars. These systems are rare, with $\sim 5$ such systems having been observed and studied (e.g., Cargile et al., 2008; Gómez Maqueo Chew et al., 2009; Gómez Maqueo Chew et al., 2019; Gillen et al., 2017; Wang et al., 2009). We suggest that this type of system is a natural result of dynamical evolution in three-body systems.

As a case study, we focus on the Par 1802 binary system (Stassun et al., 2008). This system is a relatively young (few-million-year-old, see Cargile et al., 2008; Gómez Maqueo Chew et al., 2012; Stassun et al., 2014) pre-main-sequence ‘twin’ composed of two M-dwarf stars, each with identical (to within $3\%$) masses of $0.391$ M⊙ and $0.385$ M⊙ (Gómez Maqueo Chew et al., 2012; Stassun et al., 2014). However, the temperatures of the binary stars are $3675$ K and $3365$ K, which surprisingly differ by $\sim 300$ K (or $9.2\%$) (Gómez Maqueo Chew et al., 2012; Stassun et al., 2014) despite their otherwise identical physical attributes.

The majority ($\sim 75\%$ Salpeter, 1955; Henry et al., 2006) of all the stars in our Galaxy are low mass (about $10\%-50\%$ the mass of the Sun) and cold (about $30\%-65\%$ the temperature of the Sun) stars that belong to the M-dwarf (red dwarf) class. Approximately $20-40\%$ of all M-dwarfs have a binary companion (Janson et al., 2012; Bergfors et al., 2010; Fischer & Marcy, 1992; Raghavan et al., 2010), and in a survey of $151$ contact binaries using data from the Hipparcos satellite mission, Pribulla & Rucinski (2006) found that $42\%\pm 5\%$ are at least in triple systems. Therefore, many binary stars in our galaxy are in fact in triple configurations (Tokovinin, 1997; Tokovinin et al., 2006; Eggleton et al., 2007; Griffin, 2012; Rappaport et al., 2013).

The prevalence of M-dwarf binaries in triple systems motivated Gómez Maqueo Chew et al. (2012) and Stassun et al. (2014) to study Par 1802 in the presence of a third star. In particular, they proposed that the temperature difference between the components of Par 1802 is likely due to different spin evolution and uneven mass accretion, perhaps caused by pre-main-sequence tidal evolution between the two stars. However, their work did not include the dynamical evolution of the three body system, which can dramatically alter the binary’s orbital parameters (e.g. Naoz, 2016). Here we follow the dynamical evolution of a Par 1802-like system and show that a distant tertiary can naturally explain all of the observed properties. Moreover we constrain the orbital configuration of the third companion.

Dynamical stability of such a system requires a hierarchical configuration where a third body on a far wider (outer) orbit orbits around the tighter (inner) binary of Par 1802. In this configuration, the triple body orbital parameters change on secular timescales. In particular, gravitational perturbations from the tertiary star can cause eccentricity excitations of the inner orbit, i.e., Kozai-Lidov cycles (Kozai, 1962; Lidov, 1962; Naoz, 2016). The two non-resonant orbits can be described by two “wires” where the line-densities are inversely proportional to orbital velocity. Gravitational potential can then be expanded in the semi-major axis ratio, as the hierarchical configuration drives the ratio between the inner and outer body’s semi-major axes to a small parameter. Here, we employ the hierarchical three-body secular approximation up to the octupole-level of approximation, known as the Eccentric Kozai Lidov (EKL) mechanism (e.g. Naoz, 2016). It was recently shown that the EKL mechanism plays a pivotal role in the evolution of triple stellar systems (Thompson, 2011; Naoz et al., 2013a; Naoz & Fabrycky, 2014; Naoz, 2016; Stephan et al., 2016; Bataille et al., 2018; Moe & Kratter, 2018). Gravitational perturbations from a third companion excite high eccentricities in the inner binary, thus the stars spend longer times near each other, yielding to stronger tides that can shrink the binary separation and circularising it. This process is very effective in forming short period inner binaries such as Par 1802 (e.g. Naoz & Fabrycky, 2014).

This paper provides a theoretical background (Section 2), describes the setup for the dynamical simulations (Section 3), presents and analyses results of the simulations (Section 4) with particular emphasis on predictions of the third companion’s orbital parameters. We investigate the interaction at Roche-limit and explain the observed temperature difference between the stars of the inner binary (Section 5) through a proof-of-concept mass transfer analysis and a tidal heating model. A discussion of our results is offered in Section 6 before drawing our conclusions.

2 Theory

2.1 Pre-main-sequence Evolution

Tidal evolution of such a compact binary system is very sensitive not only to orbital separation, but also to stellar radii. As we are modelling the tidal evolution during the stars’ pre-main-sequence lifetime, it is therefore very important to also consider the dramatic evolution in stellar radii during this time. To do this, we consider the calculations of Chabrier & Baraffe (1997) for the radius evolution of small ($<0.8\ M_{\odot}$) pre-main-sequence stars, wherein they consider stellar evolution as a function of mass and metallicity using grainless non-grey atmosphere models. We interpolate their results, assuming solar metallicity, for pre-main-sequence stars to find radius for each star as a function of its mass and of time while modelling the tidal evolution. The radii of these stars dramatically contracts during their pre-main-sequence lifetime, with initial inflated radii of approximately $11$ times the radius of the Sun, to final radii of approximately $1.5$ times the radius of the Sun (see Figure 1). This in turn has dramatic effects on the tidal evolution, with the tidal forces weakening greatly as the star contracts.

2.2 Analytical Expectations for the Tertiary companion

Par 1802 is a very tight binary with semi-major axis $a\sim 0.049$ au and on a circularised orbit of eccentricity $e\sim 0.02$ (Gómez Maqueo Chew et al., 2012). The temperature difference between the stars implies that mass transfer may have taken place during Roche-limit crossing phase of the evolution. The Roche-limit, $R_{\rm Roche,1}$ of a body of mass $m_{1}$ and radius $R_{1}$ with respect to an orbiting body of mass $m_{2}$ is defined as:

[TABLE]

The numerical pre-factor is highly uncertain and here we adopt $2.7$ following numerical simulations from Guillochon et al. (2011) and Liu et al. (2013). As the M-dwarf goes through pre-main-sequence contraction the Roche-limit contracts. Thus, a system may have a short phase of mass transfer during the evolution before settling down on a stable tight configuration.

The two stars probably began their life on a wider orbit than currently observed, avoiding merging at birth. The tertiary induces large values of eccentricity that eventually, with tidal evolution, drive the binaries to a tight configuration (e.g., Eggleton et al., 1998; Kiseleva et al., 1998; Eggleton & Kiseleva-Eggleton, 2001; Fabrycky & Tremaine, 2007; Perets & Fabrycky, 2009; Thompson, 2011; Shappee & Thompson, 2013; Naoz & Fabrycky, 2014). The final semi-major axis from this process is most likely the result of tidal shrinking that conserved angular momentum (e.g., Ford et al., 2000). In other-words, we can write:

[TABLE]

where subscript “$i$” refers to initial configuration (before tidal effects) and subscript “$f$” refers to final configuration, $a$ is the semi-major axis, and $e$ is the eccentricity. High eccentricity migration implies that $e\to 1$ during the evolution, thus,

[TABLE]

assuming that the closest approach before disruption can be the Roche-limit. In our case, where mass transfer may have taken place, $R_{\rm Roche}$ represents the maximum approach for interaction.

Adopting the observed semi-major axis of Par 1802 $a_{f}=0.049$ au, and using Equation (3), we find that the pericentre crosses the Roche-limit at a radius:

[TABLE]

We adopt pre-main-sequence evolution of M-dwarfs following Chabrier & Baraffe (1997), and can estimate the time the star reached this value (the vertical line in Figure 1). We thus estimate that the limit on the Roche-limit crossing time was at $1.4\times 10^{6}~{}\text{yrs}$, in other words, about $0.6$ Myrs ago.

3 Three-Body Dynamical Evolution for pre-main sequence stars

Here we discuss the secular three body evolution in the presence of pre-main sequence contraction of the stellar radius and General Relativity (GR) precession and tides. We begin by describing the EKL evolution and our Monte Carlo implementation. In Section 5.1 we discuss the possibility of mass transfer during Roche-limit crossing and calculate the resultant temperature of the stars.

3.1 Three-Body Dynamical Evolution

We study the hierarchical three-body dynamics, tidal, general relativity, and pre-main-sequence evolution of the binary system Par 1802 in the presence of a tertiary companion. We model the dynamical evolution of these systems by performing Monte-Carlo simulations of $7500$ different initial systems for $2$ Myrs, the approximate age of the Par 1802 system (Gómez Maqueo Chew et al., 2012; Stassun et al., 2014). Our Monte-Carlo implementation first involves randomly selecting $7500$ different sets of initial conditions (that satisfy the stability criteria and constraints described in Section 3.2). We then solve the three-body secular equations up to the octupole level of approximation (following Naoz et al., 2013a) through numerical integration with the Runge-Kutta 4th order method. We note that a variable step-size was used in our numerical integration due to the presence of highly eccentric orbits. Physical processes such as tidal effects and GR precession for the inner and outer orbits (e.g., Naoz et al., 2013b) are included in the numerical integration. Additionally, the effect of shrinking stellar radius during the pre-main-sequence evolution of an M-dwarf, which is of particular interest to the case of the Par 1802 binary, is included in the numerical integration.

The quadrupole level of approximation timescale for such a system is

[TABLE]

where $G$ is the gravitational constant, $a_{1}$ ($a_{2}$) are the inner (outer) semi-major axis, and $e_{1}$ ($e_{2}$) are the inner (outer) eccentricity. The tidal evolution of the inner binary is determined following Eggleton & Kiseleva-Eggleton (2001) through differential equations that account for distortion due to both tides and rotation. Tidal dissipation follows the equilibrium tidal model in Eggleton et al. (1998), with viscous time scale $t_{V}=0.5$ for all simulations (see Naoz, 2016, for full set of equations). This is a simplified treatment of tidal effects as $t_{V}$ is expected to vary with stellar radius and mass. However, the dependence of $t_{V}$ on radius and mass is unknown. The apsidal motion constants (which is twice the Love parameter) for this polytrope is $0.014$. Hence, this tidal evolution approach is adopted in our integration since it is self consistent with the secular framework, accounts for precession due to tidal torques and oblateness, and enables us to qualitatively understand the physical effect of tides in the inner binary.

While we adopted equilibrium tides, dynamical tides were shown to have significant effect on the orbital evolution (see for review Mathis, 2019). For example, in a fully convective M-type star such as Par 1802, dynamical tides in the convection zone results in a frequency-dependent tidal dissipation (e.g., Greenberg, 2009; Vick & Lai, 2018) that, at resonances, may be significantly larger than equilibrium tide dissipation (Braviner & Ogilvie, 2015). Dynamical tides have numerous potential consequences on dynamical evolution (Witte & Savonije, 2002; Auclair-Desrotour et al., 2015), such as causing the binary orbit to rapidly circularise, abrupt shrinking of $a$, and changes to the orbital period. However, we stress that the equilibrium tide model may capture the qualitative behaviour of the system (e.g., Vick & Lai, 2018)

3.2 Initial conditions

The numerical setup incorporates data from observations of Par 1802, selects initial distributions for unknown orbital quantities, and finally applies constraints and stability criteria. The masses of the constituent stars of Par 1802 ($m_{1}$, $m_{2}$) are taken to be $0.391$ $M_{\odot}$ and $0.385$ $M_{\odot}$ from Gómez Maqueo Chew et al. (2012), and the radii ($R_{1}$, $R_{2}$) taken to be $11.71$ $R_{\odot}$ and $11.64$ $R_{\odot}$, respectively from Baraffe et al. (1998). The stellar spins for both stars in the inner binary are taken to be $2$ days, and the spin-orbit angle was randomly sampled from an uniform distribution. The mass of the third outer star $m_{3}$ was randomly sampled from an uniform distribution between $0.1$ $M_{\odot}$ and $0.8$ $M_{\odot}$.

The eccentricity of the inner (outer) orbit $e_{1}$ ($e_{2}$) were sampled from a thermal distribution between [math] and $1$. The inclination $i$ between the inner and outer orbit’s angular momenta was assumed to be isotropic (i.e., uniform in cosine). The period of the inner and outer orbit was assumed to be the log-normal distribution of Duquennoy & Mayor (1991). Note that this period distribution represents the final periods of binaries population, rather than the initial one. However, if we extrapolate for a recent study for more massive stars, the final period distribution, at wide orbits, is a signature of the birth distribution (Rose et al., 2019).

We also adopt the following stability criteria. The first represents the relative importance of the octupole and quadrupole terms in the system’s Hamiltonian (Naoz, 2016):

[TABLE]

The second stability criteria follows Mardling & Aarseth (2001):

[TABLE]

Equation 6 is numerically similar to this criteria (Naoz et al., 2013a). It is noted that the dynamics of systems in the upper limit of Equation 6 may be dominated by effects other than the EKL mechanism, especially for a strong perturber (Antonini et al., 2014; Bode & Wegg, 2014). Our work does not take these cases into account because Equation 6 and 7 are numerically very close, and therefore it is expected that the vast majority of stable systems can be faithfully described by the EKL mechanism. We also note that these stability criteria deem a system unstable if at any point in time it violates the stability criterion, but do not provide a timescale for growth of the instability (e.g., Mylläri et al., 2018). In practice we expect more systems to be stable for $2$ Myr of our integration timescale than considered here.

Additionally, we require the pericentre distance to be greater than initial Roche-limit, which ensures that the inner binary does not merge before the tertiary companion can affect the system:

[TABLE]

Together, these constraints and stability criteria allow us to model systems in which the EKL mechanism may facilitate tidal heating during the pre-main-sequence evolution of the system. We depict the initial conditions of the system in Figure 2, after these criteria have been applied.

4 Three-body evolution Results

4.1 Inner Orbit Dynamical Outcomes

We identified two general dynamical outcomes. As expected, the outcome is dependent on the presence and strength of eccentricity excitations. EKL in the presence of tides tends to circularise and tighten systems. Larger eccentricity excitations, due to the EKL mechanism, results in a more efficient production of such systems. The different dynamical outcomes are depicted in Figure 3, and are summarised in Figure 4. In the latter Figure we consider the initial distribution of $a_{1}$ (top) and $e_{1}$ (bottom) and over-plot the $15\%$ of systems that crossed their Roche-limit at any point throughout the $2$ Myr evolution (blue line). The green vertical lines in Figure 4 represents the observed parameters for Par 1802. We then use the observed $a_{1}$, up to a factor of $2$, to constrain the orbital parameters of the third companion.

4.1.1 Roche-limit crossing systems

About $15\%$ of all systems crossed their Roche-limit. We identify two subtypes of Roche-limit crossing systems:

• •

Single eccentricity excitation. Depicted in top left panel in Figure 3. These systems are in the most efficient part of the parameter space of EKL. Here, the eccentricity undergos a single dip throughout the evolution, with pericentre dipping below the Roche-limit of the inner binary for a time $\Delta{t}$ and then exiting the Roche-limit. Mass transfer may occur within the Roche-limit (see Section 5.1). These systems are of interest since Par 1802 is observed to be circularised ($e_{1}=0.0166$) with small semi-major axis ($a_{1}=0.0496$ au) (Gómez Maqueo Chew et al., 2012).

• •

Long term eccentricity oscillations. Depicted in top right panel in Figure 3. In these systems EKL oscillations take place, but the eccentricity excitation is not high enough to cause the system to quickly circularise. Tides may cause an increase of $a_{1}$. In general, for these systems the inner binary does not circularise and tighten, and therefore these systems are not of interest as they do not match observations.

4.1.2 Non-Roche-limit crossing systems

In many cases, the EKL eccentricity excitation are either taking place on longer timescales than the age of the system (e.g., bottom right panel in Figure  3), or had not yet reached high eccentricity values (bottom left panel in Figure 3). These systems are not of interest since the inner binary does not circularise or tighten, and its final configuration does not match observations.

4.2 Tertiary predictions

We focus on Roche-limit crossing systems since their behavior matches the observations of Gómez Maqueo Chew et al. (2012) in which the inner binary is circularised and tightened. We show in Figure 4 that the vast majority of Roche-limit crossing systems (those with a single eccentricity excitation) have a final inner binary configuration that is circularised and tightened. Qualitatively, we note, that these systems should be a generic result of the EKL evolution. About $10\%$ of all systems that cross their Roche-limit are within in a factor of $2$ of Par 1802 observed semi-major axis of $0.0496\pm 0.0008$ au and have a circularised orbit. We use this subset to constrain the parameter space for the Par 1802 hypothetical companion. We choose to constrain our predictions using a factor of $2$ due to the uncertainties of our model, such as in the prefactor of the Roche-limit in Equation 1 and the choice of our viscous timescale in Section 3.1.

We predict the median semi-major axis for Par 1802 outer binary to likely be $350$ au.

As shown in the top panel of Figure 5, Roche-limit crossing systems (blue) have $a_{2}$ that ranges between $2-2500$ au with a median of $55$ au. This is distinctly different from the distribution for $a_{2}$ for all simulated systems (black), which has a median of $140$ au. Due to the young age of Par 1802 (and thus the short time of simulation), the final distribution for $a_{2}$ for all systems (black) closely resembles the initial conditions (e.g., top right panel in Figure 2). This result is consistent with the result for massive stars found by Rose et al. (2019).

The Par 1802-like systems (the $10\%$ of systems within factor of $2$ of the observed $a_{1}$) have a bi-modal distribution in $a_{2}$. The low value $a_{2}$ peak (median of $\sim 5$ au) corresponds to low $e_{2}$ systems. The higher $a_{2}$ values represent the majority of systems and has a median of $\sim 375$ au.

As expected the final distribution of $e_{2}$ resembles the initial distribution. The eccentricity of Par 1802-like system’s outer orbit corresponds to the outer orbit semi-major axis, with the rule of thumb that lower eccentricity requires lower (closer in) semi-major axis values. The median eccentricity in this case is $0.5$ (bottom panel in Figure 5).

The inclination of the system (top panel in Figure 6) also behaves as expected, and we find a bi-model distribution (consistent with Fabrycky et al., 2007; Naoz & Fabrycky, 2014). However, due to the short timescales of the integration, the bi-modal distribution is more closely located near $90^{\circ}$ rather than near the Kozai-angles, being centered at about $70^{\circ}$ degrees and $110^{\circ}$ degrees. The inclination of the system is correlated with $\epsilon$ (bottom panel of Figure 6), a parameter that describes the strength of the EKL mechanism. As expected, the EKL mechanism is efficient for systems that cross the Roche-limit and particularly efficient for Roche-limit-crossing systems that are within a factor of 2 of the observed $a_{1}$.

5 Interaction at Roche-Limit

5.1 Mass transfer

At the Roche-limit, mass transfer occurs. As a proof-of-concept we assume a conservative mass inversion where the binary mass

[TABLE]

and that mass transfer only occurs while the inner binary stars are within the Roche-limit. We note that Eggleton (1983) adopted the Roche Lobe radius (where a surface mass element of a star becomes unbound) as the criterion for mass transfer, with mass transfer occurring slightly earlier. Here, we instead adopt the Roche-limit as the criterion for mass transfer, which is more conservative (smaller $\Delta{t}$) than the Roche Lobe radius criterion by a small factor of $\sim 1.2$. This factor is well within the uncertainties in the model, and therefore will not lead to any significant qualitative changes to the result.

Thus, the mass transfer rate is loosely defined as:

[TABLE]

where $\Delta M=m_{1}-m_{2}$. In our proof-of-concept model, the mass transfer rate varies for different lengths of time $\Delta{t}$ spent within the Roche-limit. Since we assume mass inversion, the fraction of binary mass transferred is constant between different systems. Realistically we may expect a constant mass transfer rate rather than constant fraction of binary mass transferred.

The time that the systems spend in each-other Roche-limit is about $1000$ yr which results in a constant mass transfer rate of $\sim 6\times 10^{-6}$ M⊙ yr*-1* (see Figure 7). The Par 1802-like systems spend about a factor of $3$ longer time inside the Roche-limit, which results in a mass transfer rate of $\sim 3\times 10^{-6}$ M⊙ yr*-1*. These mass transfer rates are not unreasonable to assume and we thus proceed with this proof-of-concept calculation. We note that alternative models exist (Hamers & Dosopoulou, 2019; Dosopoulou & Kalogera, 2016; Sepinsky et al., 2007) that incorporate orbital evolution, and will be included in future work.

5.2 Tidal Heating

As the pre-main sequence stars undergo mass transfer they also feel the tidal forces that shrink and circularised the orbit (see for example Figure 3). The tidal model adopted throughout the three body dynamical evolution followed the equilibrium tidal model from Hut (1980) that assumes a constant viscous time (e.g., Eggleton et al., 1998; Eggleton & Kiseleva-Eggleton, 2001; Naoz, 2016, see the latter for a complete set of equations). However we are particularly interested in the response of the pre-main sequence stars to the energy deposited within them due to tides during the Roche-limit crossing period.

Tides on these polytropes that take orbital energy are transferred into heat in the star. This may cause the inflation of the star. Here we adopt a constant phase-lag tidal heating model (Goldreich & Soter, 1966; Wisdom, 2008; Ferraz-Mello et al., 2008), often used for brown dwarfs (e.g., Heller et al., 2010). Since the stars in questions did not begin their thermonuclear fusion reaction, this tidal model is adequate.111Note that dynamical tides have a significant effect on tidal dissipation, and may allow more efficient tidal heating (see for review Mathis, 2019). The tidal heating model used here is based on the static equilibrium tidal model and is sufficient for our proof-of-concept analysis. We note that the constant phase-lag model adopts a constant quality factor for each of the stars $Q_{j}$, with $j=1,2$. However, since $a_{1}$ and $R_{j}$ of the stars ($j=1,2$) vary in time, the quality factor at the onset of Roche-lobe crossing is dramatically different than the quality factor at the existing point. The evolution of $Q$ strongly impacts the tidally driven orbital evolution of systems (Heller, 2018).

For the following proof-of-concept calculation we adopt two quality factors to calculate the excess temperature due to tidal effects. First, we consider the relation between the quality factor and the viscous time:

[TABLE]

from Naoz et al. (2016), with constant viscous timescale $(t_{V})_{j}=0.5$ typical for M-dwarfs, classical apsidal motion constant $k_{1}=0.014$, and where the mean motion is $n_{1}={2\pi}/{P_{1}}$. Since $R_{j}$ and $a_{i}$ varies in time, $Q_{j}$ varies in time.

Tidal heating takes place throughout the tidal evolution. However, it is most significant when the stars are within each-other Roche-Lobe regimes (see upper panel of Figure 7 for time scale). We thus limit ourselves to two relevant points in time of the evolution: the onset of the Roche-limit, and exit from this regime. We find that for the majority of cases $Q_{j,\rm enter}>>Q_{j,\rm exit}$, where the subscript “enter” (“exit”) refers to entering (exiting) the Roche-limit. While tidal heating takes place for both stars, the $3\%$ mass difference results in slight changes in the amount of heating per unit mass, with the less massive star being heated more.

Under the tidal heating, constant phase-lag model, the potential of the less massive body is expressed in terms of periodic contributions of tidal frequencies at various phase-lags (e.g., Heller et al., 2010), which can be expanded in terms of these phase lags. The phase lags $\epsilon_{\eta,j}$ for the $j$th body with $\eta=0,1,2,5,8,9$ are:

[TABLE]

with $\Sigma(x)$ defined as the algebraic sign of $x$:

[TABLE]

and $\Omega_{\text{s},j}$ is the rotational period associated with the spin of the stars.

The change in orbital energy of star $1$ due to star $2$ (or vice versa) is then:

[TABLE]

where $j,l=1,2$, and the change in rotational energy is

[TABLE]

where $k_{L}=0.5$ is the dynamical Love parameter and $\psi_{j}$ is the obliquity of the $j$th body. All of the relevant variables are obtained through the dynamical evolution results described in Section 4.1, and the masses of the stars invert between Roche-limit entry and exit (see Section 5.1). Finally the change in total energy due to tides is

[TABLE]

The temperature change for each body due to tidal heating is then evaluated from the Stefan-Boltzmann equation

[TABLE]

where $\sigma_{\text{SB}}$ is the Stefan-Boltzmann constant, and $T_{\text{eff},j}$ is the temperature at Roche-crossing for the $j$th body of mass $m_{j}$ adopting Baraffe et al. (1998) model. We then estimate the above equation at two times: when the system enters and exits the Roche-Lobe regime222In reality, heating takes place throughout the time the stars spend in the Roche-limit. However, since we find that the maximum heating takes place when the stars exit their Roche-limit, we evaluate heating only at these two points in time.. In other words, $dT_{j}=dT_{j,\rm{enter}}+dT_{j,\rm exit}$, where subscript enter (exit) refers to entering (exiting) the Roche-Lobe regime. We note that $T_{j,\rm eff}$ at the entering and exiting points of the Roche-limit are roughly the same because the short amount of time the stars spends in each-other Roche-lobe regime (see Figure 7). The temperature difference is then estimated by:

[TABLE]

In Figure 8 we show a histogram of the final excess temperatures as a result from the heuristic calculation of tidal heating. We find that Roche-limit crossing systems have a wide range of $\Delta{T}$, with a maximum $\Delta{T}_{\text{max}}=910$ K (Figure 8). We predict a maximum $\Delta{T}_{\text{predict}}\approx 115$ K.

The luminosity of the two stars is estimated by summing the total energy change due to tides for each star using Equation 16 at entry and exit from the Roche-limit. The luminosity difference between Par 1802-like stars (with $a_{1}$ within a factor of 2 from observation) has a median of $\sim 3\times 10^{25}$ W, which is comparable to the observed luminosity difference of $\sim 7\times 10^{25}$ W (Gómez Maqueo Chew et al., 2012).

Since the constant phase-lag tidal heating model typically adopts a constant $Q$, we complete a second set of calculations of $\Delta{T}$ for an assumed constant $Q=500$. We predict a maximum $\Delta{T}\approx 180$ K for constant $Q=500$. Both predictions ($\sim 115$ K for $Q$ following Equation 11, and $\sim 180$ K for constant $Q=500$) are comparable to the observation of $\Delta{T}_{\text{observed}}\approx 300$ K (Stassun et al., 2014). With $Q=500$, the luminosity difference between Par 1802-like stars has a median of $\sim 5\times 10^{26}$ W, comparable to observation.

6 Discussion

The Par 1802 M-dwarf binary system is an interesting quandary for stellar evolution since one of the component stars has unexpectedly higher temperature compared to the other despite their almost identical (within $3\%$) masses. Previous work by Stassun et al. (2014) and Gómez Maqueo Chew et al. (2012) have proposed that the temperature difference was caused by the presence of a third star. Here we focused on the dynamical evolution of such a system and have shown that observed quantities can be naturally reproduced. Furthermore, we predict the orbital configuration of such a tertiary companion.

We studied the secular evolution of the binary M-dwarf system Par 1802 in the presence of a third, far away star. We numerically solved the three-body hierarchical secular equation to the octupole level of approximation (i.e., the EKL mechanism Naoz, 2016). Tidal effects, GR precession, and pre-main-sequence evolution were self-consistently included in our simulations.

During the evolution, the EKL mechanism could cause large eccentricity excitations in the inner orbit, which prompted the inner binary to tighten and circularise (see top left plot in Figure 3). In some of these systems, eccentricity excitations caused the inner binary stars to cross their Roche-limits (see Figure 3). During that time, mass transfer could take place, resulting in the aforementioned temperature difference. As the stars’ radii contracted their Roche-radii shrunk, effectively halting mass transfer, and the stars exit the Roche-lobe regime.

We found that $\sim 15\%$ of all simulations were Roche-limit crossing systems with the majority ($\sim 80\%$) of these having circularised and their separation shrunk due to tides (see Figure 4). We focused on systems that evolved to within a factor of $2$ of the observed Par 1802 separation (of $\sim 0.05$ au). These systems represented about $10\%$ of all Roche-limit crossing systems, and about $1.5\%$ of all simulated systems. Since our initial conditions for semi-major axis and eccentricity featured broad distributions (with $a_{1}$ varying between $\sim 0.2-380$ au and $e_{1}$ varying between $\sim 0.1-0.9$), this $1.5\%$ fraction is statistically significant and suggests that systems similar to Par 1802 can exist in the presence of a tertiary companion.

For these systems similar to Par 1802, we constrain the orbital configuration for the Par 1802 tertiary star. In particular, we predict that the outer orbit is most likely to have a separation of $350$ au (top panel Figure 5). We find that the predicted median eccentricity of the tertiary is $e_{2}=0.5$ (bottom panel Figure 5). The inclination, (as expected, e.g., Naoz & Fabrycky, 2014), has a bi-modal distribution with peaks centered around $70^{\circ}$ and $110^{\circ}$.

During the time the system crossed its Roche-limit, mass transfer took place and tidal heating was probably most efficient. The two stars of Par 1802 have an observed temperature difference of $\sim 300$ K (Gómez Maqueo Chew et al., 2012). We conducted a proof-of-concept calculation, assuming conservative mass transfer during the Roche-limit crossing time. For simplicity the two inner masses were set to be the observed masses, and thus, to match observations the mass ratio inverted over the course of the evolution. These assumptions determined the mass transfer rate (see Figure 7). Following the Heller et al. (2010) tidal heating model we find about $0.2\%$ of all systems have temperature difference within a factor of $2$ of the observed value. These systems represent a small fraction ($5\%$) of the Par 1802-like systems (i.e., systems with final inner binary separation within a factor of 2 of the observed separation). Since our mass exchange and tidal heating analysis are proof-of-concept and only accounts for heating at entry and exit from the Roche-limit, we realistically expect a higher fraction of Par 1802-like systems to exhibit the observed temperature difference.

Finally, our results suggest that secular interaction in hierarchical triple body systems, during their pre-main-sequence contractions, may lead to a generic feature: tight binaries with temperature difference. Par 1802 is a young system ($\sim 2$ Myr), which results in limited time for the EKL mechanism to operate. We thus expect similar, more pervasive behaviour for general triple systems. In particular, the presence of a tertiary star may explain the $4$ other twin binaries with significant temperature differences discovered so far (Gómez Maqueo Chew et al., 2009; Gómez Maqueo Chew et al., 2019; Gillen et al., 2017; Wang et al., 2009). This dynamically-induced temperature difference may affect conclusions of the coevality of binary systems (e.g., Kraus & Hillenbrand, 2009).

Acknowledgements

S.J.C. would like to acknowledge the support of the Lau Endowment as part of the Undergraduate Research Scholars Program of the UCLA Undergraduate Research Center - Sciences, the Undergraduate Research Fellows Program, and the Maggie Gilbert Research Award through UCLA College Honors. S.N. thanks Howard and Astrid Preston for their generous support. We would like to thank Professor Stassun, Professor Gómez Maqueo Chew, Professor Heller, and Dr Hamers for their thoughtful comments and suggestions.