Precession and accretion in circumbinary discs: the case of HD 104237

Abstract

1 INTRODUCTION

It is widely recognized that most Sun-like stars form in binary systems (e.g. Duquennoy & Mayor 1991; Raghavan et al. 2010; Janson et al. 2012; De Rosa et al. 2014) and the binary fraction increases towards earlier spectral type. It is also known that binarity increases towards younger ages in pre-main-sequence objects (e.g. Reipurth et al. 2014). It should therefore be unsurprising that binaries are common in Herbig Ae/Be systems (68 ± 11 per cent of the full Herbig Ae/Be population, with 35 per cent being spectroscopic binaries; Corporon & Lagrange 1999; Baines et al. 2006), where intermediate-mass (2 ≲ M_⋆ ≲ 8 M_⊙) stars at young ages are accompanied by gaseous circumstellar (or as it may be, circumbinary) discs.

It is also well known that circumbinary discs can affect the dynamical evolution of the binary through resonant interactions, acting as an angular momentum reservoir and changing the orbital elements of the binary (Artymowicz et al. 1991; Artymowicz & Lubow 1994). The direction of influence is not one-way however, and the potential of the binary dictates the form of the inner gap or cavity formed in the disc (Artymowicz & Lubow 1996) and how gas is able to accrete through it (Artymowicz & Lubow 1994).

The importance of how binaries interact with the inner part of the disc is beginning to gain widespread recognition. The accretion process inside the central cavity is vital for understanding how both stars and their discs evolve and there are an increasing number of observations seeking to characterize this process (at first from SED modelling and now from resolved mm/sub-mm images; e.g. Jensen & Mathieu 1997; Andrews et al. 2011; Harris et al. 2012). Very recently, Dutrey et al. (2014) found evidence of ongoing planet formation in the disc around the young binary GG Tau A, identifying accretion streams flowing from the outer disc and feeding a subdisc around one star. Understanding how gas flows from the outer circumbinary disc into the central cavity is therefore crucial to understanding planet formation processes in these systems. In particular, close binaries (common in Herbig systems) offer the opportunity to study this process as a function of the binary's orbital phase as the dynamical times in these systems are short.

Due to the non-Keplerian potential it creates, even a circular or equal-mass binary can create an eccentric shape to the disc cavity (e.g. MacFadyen & Milosavljević 2008), and this is indeed observed in widely studied objects such as AK Sco (Gómez de Castro et al. 2013). An unequal-mass binary can also cause the cavity to become strongly decentred from the binary centre of mass as it takes on this eccentric shape. In addition, the non-Keplerian potential of the binary can force the disc's orbital elements to change in the same way that a circumbinary planet's elements osculate (Leung & Lee 2013).

For solar-type and intermediate-mass binary stars on the main sequence, the mass ratio q_b is roughly flat for most stellar types and companion masses (e.g. Raghavan et al. 2010; Janson et al. 2012), but with a strong bias towards equal masses at small binary separations (De Rosa et al. 2014). However, Kouwenhoven et al. (2005) found that for a sample of young A-stars, the mass ratio is biased towards unequal masses, indicating that accretion at young ages plays an important role in changing the primordial mass ratio into the one observed in older systems. This is supported by simulations that find accretion occurs primarily on to the secondary (e.g. Artymowicz 1983; Bate 2000; de Val-Borro et al. 2011; Bate 2014), because its Roche lobe is further from the barycentre of the system and the relative velocity between the accreting gas and the secondary is lower. This allows the secondary to accrete more gas, pushing the mass ratio towards unity (Cuadra et al. 2009; Roedig et al. 2011; Clarke 2012; Young, Baird & Clarke 2015). However, moving towards higher eccentricity can change this, and Roedig et al. (2011) found that for high eccentricity, the binary components can accrete at comparable rates with a mass ratio of q_b = 1/3. This should affect the mass ratio of the final binary, and Halbwachs et al. (2003) indeed found that for stellar ‘twins’, there is a small preference towards lower eccentricities, indicating that more eccentric binaries have less equal-mass ratios. However, this has not been studied in detail to date.

Simulations studying accretion in eccentric binary systems have shed a little light on this process. Günther & Kley (2002) found that the accretion rate was dependent upon orbital phase for a number of eccentric systems. More recently, Shi et al. (2012) found that an eccentric binary drives a disc precession in high-resolution magneto-hydrodynamic (MHD) simulations, but did not run the simulation for more than a fraction of the precession time-scale. Gómez de Castro et al. (2013) also found an eccentric and decentred disc cavity in simulations of GG Tau A, but again the time-scale of their simulations was short.

The modelling work presented in this paper has been motivated by the Herbig Ae system HD104237 (DX Cha). HD104237 is a nearby (d = 116 ± 7 pc; Perryman et al. 1997) spectroscopic binary (SB2) with a semi-major axis a_b = 0.22 au, eccentricity e = 0.6 and orbital period T_b = 20 d (Böhm et al. 2004; Garcia et al. 2013). The total stellar mass is 3.6 M_⊙ and the mass ratio inferred from the spectral type of the primary is q_b = 0.64. The system is inclined by i ≲ 20° to the plane of the sky (Böhm et al. 2004; Grady et al. 2004; Garcia et al. 2013).

Because of its near face-on geometry, proximity and brightness, this system is currently a unique and prime target for detailed study of the accretion process in a close binary system with long baseline interferometry. AMBER/VLTI observations show that most of the K-band continuum flux arises from a ring of radius R = 0.5 au (Tatulli et al. 2007), in agreement with the location of the dust sublimation radius in the circumbinary disc (Garcia et al. 2013). However significant unresolved emission is present within this radius, in excess of the expected contributions from the photospheres. The stars in this system probably do not hold individual circumstellar discs due to outer tidal truncation at closest binary approach.

This unresolved flux likely arises in compact structures (e.g. accretion streams) inside the tidally disrupted circumbinary disc (Garcia et al. 2013). The AMBER observations also show that most (90 per cent) of the HI Brγ emission is unresolved and comes from the immediate vicinity of the stars. An increase in Brγ emission is observed on each binary component at periastron passage, likely related to an accretion burst (Garcia et al. 2013), similarly to what has been observed previously. We have recently observed this system with PIONIER/VLTI in the H band at different orbital phases, with the aim of reconstructing the morphology of this tidally disrupted circumbinary disc and study the evolution of the accretion process on to the central binary with orbital phase. The description and analysis of these data, in particular their detailed comparison to the results of the modelling work presented here, will be presented in a forthcoming paper (Dougados et al., in preparation).

To build on and extend the findings of work by previous authors (e.g. Shi et al. 2012; Gómez de Castro et al. 2013), we perform smoothed particle hydrodymanics (SPH) simulations of this system, with an emphasis on exploring the long-term evolution of accretion patterns in the system (over more than 1000 orbital periods). The structure of the paper is as follows. In Section 2 we describe the numerical parameters of our simulations and the method used to carry them out, and in Section 3 we describe the results. In Section 4 we apply these findings to the general case of a binary of any given eccentricity and compare the precession time-scale of our SPH simulation with the predictions of Leung & Lee (2013). We discuss the implications for the study of discs around stellar binaries, in particular the role of accretion in eccentric systems, and the limitations of our methods and model in Section 5, before presenting our conclusions in Section 6.

2 SIMULATIONS

We have simulated the HD 104237 system using the hybrid SPH/N-body code Gadget-2 (Springel 2005), modified to closely follow the dynamics of the binary and to include a Navier–Stokes viscosity (Dunhill, Alexander & Armitage 2013). We model the binary components as N-body particles with a mass ratio q_b = 0.64 (so that M_A = 2.2 M_⊙ and M_B = 1.4 M_⊙), semi-major axis a_b = 0.22 au and eccentricity e_b = 0.6.

We model the disc using SPH particles, with N = 2 million particles for our main simulation. We parametrize the disc viscosity using the Shakura & Sunyaev (1973) α-disc prescription, with α = 10⁻². The disc initially extends between 1.5 < R < 20 a_b, with a surface density profile that approximately follows Σ ∝ R^−1.7, but is not a strict power law.¹ We impose an isothermal equation of state on the disc, so that the temperature (and thus the disc scale-height H) is a function only of radius from the barycentre of the binary. We choose the temperature profile so that it results in a flaring disc where H/R ∝ R^1/4, normalized so that H/R = 0.05 at R = 1 au. We model the stellar accretion using a simple sink particle method, where each star is modelled as a point mass with an associated sink radius, R_sink = 0.03 a_b. Any SPH particle found within this radius has its mass and momentum added to the star, and the SPH particle is removed from the simulations. We discuss the physicality of this implementation and its implications in Section 5.1. For calculations of the accretion rate discussed in Section 3 we refer to the rate at which SPH particles are ‘accreted’ by the stars in this way. The binary components are modelled gravitationally as point masses and their orbit is free to evolve through the simulation. As the disc has a low mass this only occurs at a level below one per cent of the initial orbital parameters on the time-scales of our simulation, and we therefore do not report these changes.

Our disc has a mass M_d = 5 × 10⁻³M_b = 0.018 M_⊙. Hales et al. (2014) estimated the total disc mass to be M_d = 0.04 M_⊙ from their SED modelling. Their model suggests a characteristic radius of 90 au, and our simulated disc is much smaller than this, so our disc is likely to be overmassive. The smaller radial extent of our disc is due to numerical reasons, but our neglect of all but the inner regions is justified as the outer part of the disc should not play any dynamical role in what happens in the centre. However, it does affect the mass in the inner regions as the disc spreads viscously, and we therefore normalize the accretion rates discussed in Section 3 to account for the fact that while the surface density in our disc decreases with time, in the real disc around HD 104237 this process occurs only on time-scales much longer than we simulate. Since we are well below the self-gravitating limit the mass effectively becomes a scaling factor and does not affect the hydrodynamics. SED models of HD 104237 are consistent with a non-flaring disc (Fang et al. 2013; Hales et al. 2014), while our equation of state enforces a flared disc. The flaring is not strong, and H/R only goes from 0.03 at the inner edge to 0.07 at the outer edge of our disc, so we do not expect this to affect the results of our simulations. Our model also does not include stellar winds or the jet associated with the system.

We use a set of code units such that the distance unit R₀ = a_b, the mass unit M₀ = M_A + M_B and the time unit T₀ = T_b/2π where T_b is the orbital period of the binary. This choice of units sets the gravitational constant G = 1 internally to the code.

After the initial conditions settle, the action of the binary drives an eccentric inner cavity in the disc. Material drawn in as a tidal tail by the primary is then kicked across the cavity, and pushes the far cavity wall away from the binary. A snapshot showing this process is shown in Fig. 1, and we show how this evolves as a function of orbital phase in Fig. 2. The eccentric cavity shape is then key to the evolution of the system in our simulations.

Figure 1.

Snapshot of our reference simulation with N = 2 million particles after 159 binary orbits, with the binary at pericentre. The clearly eccentric cavity shape is formed by the unequal-mass, eccentric binary ‘slingshotting’ material across the cavity and pushing the far wall of the cavity further from the binary.

Open in new tabDownload slide

Figure 2.

Series of snapshots showing the accretion process as a function of the binary's orbital phase, with time T in units of the binary period T_b. The binary is at pericentre in the uppermost left panel (the same snapshot as in Fig. 1), and at apocentre in the panel third from left in the middle row. The primary (larger circle) can be seen to drag disc material in from the nearest cavity edge at apocentre passage, which is then thrown across the cavity due to the binary's eccentric orbit. This stream is then partially dragged back towards the binary by the secondary when it next reaches apocentre, and this material is accreted by the binary components near pericentre. As the cavity precesses with respect to the binary orbit, which binary component plays the role of creating the accretion stream and which drags it back to the binary will switch.

Open in new tabDownload slide

The inner region of the disc then begins to precess with respect to the orbit of the binary, as shown in Fig. 3. This precession is due to the highly non-Keplerian nature of the potential close to an unequal-mass eccentric binary such as HD 104237, and is a natural consequence of the fact that orbits in a non-Keplerian potential are not closed ellipses. To date, this effect has only been studied for discs around circular binaries in the context of supermassive black hole (SMBH) binaries (MacFadyen & Milosavljević 2008; Shi et al. 2012; D'Orazio, Haiman & MacFadyen 2013; Farris et al. 2014). We find that the precession is uniform between 4 ≲ R ≲ 6a_b. We ran the simulation for 1500 binary orbits before arbitrarily halting it.

Figure 3.

Series of snapshots showing the precession of the cavity. Times are in units of the binary orbit and the binary is at pericentre in each snapshot. At these early stages in the simulation the precession rate is changing so the time interval between the panels is not uniform (see Fig. 5).

Open in new tabDownload slide

As our simulation is only of moderate resolution (especially inside the cavity), we wanted to conclusively rule out the possibility that this precession rate was affected by any anomalous resolution or viscosity issues in the SPH. We therefore took a snapshot after approximately T = 450 T_b and resampled it, replacing each SPH particle with four new particles at positions randomly sampled from the smoothing kernel and the parent particle's smoothing length. We then ran this higher resolution version of the simulation for a further 180 binary orbits. Although the up-sampling introduces noise into the new simulation at the very start, this is soon dissipated and the simulation follows the same course as the original. A comparison between the original and up-sampled simulations at the time we halted the higher resolution run can be seen in Fig. 4. The only significant difference between the two simulations is that in the higher-resolution version the accretion rate is lower as we are better able to resolve the individual accretion streams.

Figure 4.

Comparison between the original (top panel) and higher-resolution run (bottom panel) at the end of the latter run. The accretion streams within the cavity are noticeably better-resolved in the higher-resolution version, and the cavity shape is slightly more eccentric. We compare cavity precession rates and the measured accretion rates between the two runs in Figs 5 and 7.

Open in new tabDownload slide

3 RESULTS

Figure 5.

Evolution of the disc eccentricity (e_d, upper panel), semi-major axis (a_d, middle panel) and argument of periapse (ϖ_d, lower panel). Thin solid lines are for our main simulation and bold lines are for the shorter higher-resolution version. These values are averaged between 4.5 < R < 5.5 a_d, where the disc precesses uniformly. In the lower panel, we also plot a best fit to the precession, where |$\dot{\varpi } = 0{^{\circ}_{.}} 48\,\,T_{\mathrm{b}}^{-1}$|⁠.

Open in new tabDownload slide

How the accretion proceeds inside the cavity is strongly dependent upon the position angle of the cavity. This results in a periodicity in the accretion signature in our simulation that is wholly tied to the precession of the disc. In Fig. 6 we show the accretion as a function of orbital phase at the same four points in the simulation as the panels of Fig. 3, showing that the peak accretion rate per orbit, and which component it is on to, changes with the position angle of the cavity. To generate this figure we average over 10 orbits either side of the snapshots shown in Fig. 3 in order to reduce noise and highlight this short-term trend.

Figure 6.

Accretion rates as a function of orbital phase at four different times during the simulation. Times 0 and 1 correspond to the binary at pericentre and 0.5 corresponds to the binary at apocentre. The accretion rates are averaged over 10 binary orbits, and smoothed using a moving-average algorithm. It is noteworthy that the binary component with the highest accretion rate changes as the simulation evolves, as this is dictated by the precession of the disc cavity. Also of note is that the highest accretion peak occurs just before pericentre, and the lower accretion peak occurs at apocentre. This is due to how the accretion proceeds within the cavity (see Fig. 2), but is not necessarily correct due to our neglect of proper thermodynamic treatment (see Section 5).

Open in new tabDownload slide

Fig. 7 shows how these peak accretion rates evolve through the simulation. As we neglect the outer regions in our simulation, the surface density of the disc decreases over time due to viscous spreading. In the real HD 104237 disc, the outer disc acts as a mass reservoir so this does not occur. To account for this, we normalize the accretion rate to the remaining disc mass within R ≤ 5a_b. We include a comparison with the higher-resolution resample as an inset panel, showing how increasing the resolution changes the absolute accretion rate but that the trend remains unchanged.

Figure 7.

Measured accretion rates in our main simulation, with the accretion rates for the up-sampled higher resolution run inset. The periodic change in the accretion rate is dictated by the changing position angle of the cavity, which can be easily seen by comparing this with the lower panel of Fig. 5. The accretion rates here are a factor of 10³−10⁴ higher than that measured for the real system (⁠|$\dot{M} \sim 10^{-8}$| M_⊙ yr⁻¹; Grady et al. 2004). This is primarily due to the low resolution of our runs (the inset panel clearly shows that the increased resolution run has accretion rates that are lower by an order of magnitude) and the higher disc mass in our simulation. The higher resolution runs, in which the accretion is better (but not fully) resolved, show good agreement with the general pattern of the primary run, with both primary and secondary rates turning over at T ∼ 550 T_b as the cavity precesses. Grey bars indicate the times when the binary and disc cavity position angles coincide (i.e. ϖ_b = ϖ_d).

Open in new tabDownload slide

We see from these three figures that there are two strong periodicities obvious in the accretion rates we find in the simulation. First, on a per-orbit time-scale as shown in Fig. 6, the binary component which sees the highest accretion rate peaks just before the binary reaches pericentre, with the lesser accretor peaking at pericentre. This is due to the fact that the binary accretes by pulling material back in from the accretion stream that has been thrown out across the cavity, and this reaches the binary just before pericentre and is accreted, with a little remaining gas then falling on to the other component just as the binary reaches pericentre (see Fig. 2).

On longer time-scales, and in a pattern dictated entirely by the non-Keplerian potential of the binary, we see accretion variability as the cavity precesses and the process described above is modulated by the changing amounts of gas available in the accretion streams and switching over of the binary component that creates the streams (see Fig. 3). The time-scale for this process to change is half the precession period of the disc, corresponding to approximately 20 years for the real HD 104237 disc.²

Taking the precession-period averaged accretion rates |$\langle \dot{M}_{\mathrm{A}} \rangle / \langle \dot{M}_{\mathrm{B}} \rangle$| where |$\dot{M}_{\mathrm{A}}$| and |$\dot{M}_{\mathrm{B}}$| are the accretion rates on to the primary and secondary respectively, we find values between 0.9 and 1.3. This means that this precession mechanism allows the primary to accrete at the same rate as, and at times at an even higher rate than, the secondary. In purely physical terms then, this provides a viable mechanism for allowing an eccentric binary to accrete matter while maintaining an unequal-mass ratio, something that has proven problematic for simulations of binary formation (e.g. Clarke 2012; Young et al. 2015).

We believe that the first variability, on single-orbit time-scales, is indeed observed in the real system. As reported by Garcia et al. (2013), the equivalent width of the Brγ line increases around pericentre passage, and this is interpreted as evidence of enhanced accretion at this phase of the orbit. The observability of the second periodicity depends upon how the accretion luminosity and physical accretion rate are related, which is far beyond the physics included in our SPH simulation. However, it is instructive to see at what level we might expect to see the accretion to vary on these 20 year time-scales, so we use the most basic possible model for accretion luminosity – gravitational energy lost is converted into luminosity – to make very rough estimates. Although this model is known to be incorrect in many ways (e.g. Da Rio et al. 2014), we do not necessarily trust all aspects of the accretion process that we do model in the simulations, so it would be over-complicating matters to apply a model too sophisticated.

Although it is possible that the stars have small sub-discs around them which would provide a mechanism to delay accretion on to the stars themselves, we think this unlikely and do not include it in our model. This is because the truncation and stellar radii are very similar, so any disc would have to be very small indeed (Pichardo, Sparke & Aguilar 2005; Garcia et al. 2013). Instead we model the accretion as falling directly from the corotation radius of each star (Da Rio et al. 2014). As the two components of the binary cannot be distinguished at pericentre, we sum their luminosities and plot the resultant accretion luminosity as a function of time in Fig. 8.

Figure 8.

Estimated luminosity evolution of the HD 104237 binary using equation (2). We use values of |$\dot{M}$| taken from Fig. 7, scaled down by a factor of 10³ to match the observed accretion rate of the system (Grady et al. 2004). Values for the stellar masses were taken from Garcia et al. (2013) and stellar radii from Fumel & Böhm (2012). As the components cannot be distinguished at pericentre (where accretion peaks), we sum the luminosities of the components. The variability we find using this very simplified model is less than the uncertainty in the luminosity of the primary (e.g. Garcia et al. 2013), so we would not expect the signal to be observable for the case of HD 104237.

Open in new tabDownload slide

We find that the variability expected even for this most optimistic model (where we assume all lost gravitational energy is radiated as observable signature) is lower than the uncertainty in the luminosity of the primary star of the system (⁠|$L_{\mathrm{A}} = 30^{+5}_{-4}\,\mathrm{L}_{\odot }$|⁠; Garcia et al. 2013). We therefore do not expect the signature of precession-driven accretion variability to be observable over decade-long time-scales if the accuracy of these measurements does not improve dramatically. However, it is possible that this signature could be observed for eccentric short-period binaries with higher accretion rates, as the signal would be stronger in such cases. Furthermore, if the uncertainties in the observed stellar luminosities are systematic rather than intrinsic, and apply equally to measurements taken years apart, it is possible that this ΔL_acc ∼ 2 L_⊙ variability will indeed be observable.

4 APPLICATION

Our primary result is that precession of the circumbinary disc can cause long-term periodic accretion variability. In order to make this more widely applicable, we ran a large number of simple three-body integrations of a test particle orbiting binaries of different mass ratios and semi-major axes. We varied the initial eccentricity and semi-major axis of the test particle orbit, and for each configuration we ran 100 integrations with the test particle's initial orbital phase drawn at random from a uniform distribution between 0 and 2π. We explored test particle semi-major axes in the range 2.5 ≤ a_t ≤ 8.5, and eccentricities 0.1 ≤ e_t ≤ 0.5. The test particle eccentricity made no difference to the precession rates of the particles in the range we explored, so we average the results over the eccentricities. We tested six different binaries, with mass ratios q_b = 1/3 and 1/2, each with eccentricities e_b = 0, 0.3 and 0.6. In each case we integrated the system for T = 10 000 T_b, but stopped if scattering from a close interaction caused e_t ≥ 0.8. For each of the six binary configurations we ran 5000 integrations, for a total of 30 000 runs.

We then measured the change in the test particle's argument of periapse, and binned these values to create a histogram of precession rates for each configuration of binary and test particle orbit. Fitting a Gaussian to this gives us a simple estimate of how the precession rate of the test particle varies with its own and the binary's orbital parameters. We plot these precession rates in Fig. 9, with the standard deviation of the fitted Gaussian as an error bar to show the width of the precession rate distribution for each set of parameters.

Figure 9.

Apsidal precession rates for test particles around binaries with mass ratio q_b = 1/3 (left panel) and q_b = 1/2 (right panel), and binary eccentricities e_b = 0 (red points), 0.3 (blue points) and 0.6 (green points). Red and green points are, respectively, offset slightly negatively and positively on the x-axis to increase visibility. Error bars are the standard deviation of normal distributions fitted to histograms of precession rates measured from the test particle runs. Dashed lines in each case give the Leung & Lee (2013) prescription for the precession rates given in equation (3). We find excellent agreement, except for the case of circular binaries at large semi-major axis, and at small semi-major axis.

Open in new tabDownload slide

We find that this prescription gives a very good approximation to our test particle precession rates, even for high eccentricity (e_b = 0.6) binaries. In fact, we only find significant disagreement for two cases: at low a_t, and for large a_t around circular binaries. The former is easily understood, as Leung & Lee (2013) derive equation (3) by taking only the lowest power of a_b/a_t in an expansion. The latter disagreement is more puzzling, but as we are interested in using their theory to predict precession time-scales for discs around high eccentricity binaries, finding the source of this discrepancy where e_b = 0 is beyond the scope of this paper and we do not pursue it further.

We also repeated the test particle experiment for a binary that matches that in our SPH simulation, performing a further 8000 integrations of test particles with 4 ≤ a_t ≤ 8.5 and 0.1 ≤ e_t ≤ 0.8, and plot the results in Fig. 10 along with the values found for our disc shown in Fig. 5 and the Leung & Lee predictions using equation (3). Again, we find good agreement with the prescription of Leung & Lee (2013), both for the test-particle integrations and for the SPH disc. We conclude that while resonances with the eccentric binary that act to carve the gap in the first place may have some effect on the stability of gas orbits at the cavity edge, their strength is far exceeded by the non-Keplerian potential of the binary (and similarly for the effect of pressure gradients in the disc). As the only stable orbits in this potential are precessing, the gas orbits do so too and their motion is well approximated by that expected for a massless test particle.

Figure 10.

As Fig. 9 but for a binary matching the binary in our SPH simulations. The dashed line give the Leung & Lee (2013) prescription for the precession rates given in equation (3), and the red circle is the measured rate from our simulations (values taken from Fig. 5). We find good agreement between our SPH precession rate and that predicted by Leung & Lee (2013), consistent with the uncertainty from the test particle runs.

Open in new tabDownload slide

5 DISCUSSION

5.1 Numerical limitations and omissions from the model

There are a number of considerations that must be made in interpreting our results, given the simple disc model and the moderate (and in places poor) resolution of our primary simulation. In the region of uniform precession (R ∼ 5 a_b), we resolve the disc vertical structure into four mid-plane smoothing lengths, narrowly satisfying the requirements of Nelson (2006) to avoid under-estimating the mid-plane density. At larger and smaller radii (most especially inside the cavity), we fall below this and so our vertical disc structure in these regions is likely incorrect. However, vertical resolution is not the primary problem within the cavity as the gas there no longer behaves as if it were in a disc. The poor vertical resolution at large radius equally does not affect our results, as the high-resolution re-run is resolved here and shows no differences.

Of greatest concern is the poor resolution inside the cavity. In our main simulation with N = 2 × 10⁶ particles, each accretion stream typically consists of only a few thousand SPH particles – as our code uses 50 nearest neighbours within the kernel, this is poor resolution indeed. In particular, in most of the stream the width is at best the same as the smoothing lengths of the particles in it, and wider than that towards the end of the stream being accreted.

The situation is less dire in the up-scaled run with higher resolution, where an accretion stream is typically ≳10 000 particles and is nearly always more than one smoothing length across. Given that Figs 5 and 7 show that both the disc and accretion follow the same pattern in the low- and high-resolution runs, the poor resolution within the cavity does not seem to affect the results significantly. Indeed, the only difference is that the accretion rate is lower by approximately an order of magnitude, and this is a result of the better resolved accretion process and lower individual particle mass in the high-resolution case.

Although the accretion rate is lowered by increasing resolution, the pattern followed by the relative accretion rates on to the components (i.e. |$\dot{M}_{\mathrm{A}}/\dot{M}_{\mathrm{B}}$|⁠) remains unchanged. In calculating the potential change in accretion luminosity in Section 3, we scale the raw accretion rates from our simulation to match the observed accretion rate seen in HD 104237. Although these raw numbers are unconverged between our simulations, it is instead this pattern in the relative accretion rates that dictates the signal shown in Fig. 8. Therefore, it is unlikely that this signal is negatively affected by the unresolved accretion in our primary simulation.

Beyond these resolution issues, our assumption of a locally isothermal equation of state is problematic. This assumption may be valid far from the central binary, where the propagation time for photons scattering through the optically thick disc is longer than the dynamical time of the binary, but closer to it creates issues. In particular, Marzari et al. (2012) showed that the inclusion of a radiative treatment reduces the disc eccentricity, and has something of an effect on the disc's argument of periapse, though they did not run their simulations for long enough to investigate the precession effect we find. This reduced eccentricity could prevent the accretion variability pattern seen in our simulations from occurring, as the process relies on the pericentre of the cavity being significantly closer than its apocentre.

We also neglect proper consideration of how the accretion process itself proceeds. Correct thermodynamic treatment should consider where the gas shocks and how efficiently it is able to radiate this away, and as a result this can affect which binary component the gas is accreted on to (e.g. Clarke 2012; Young et al. 2015). We also make no checks for boundness before swallowing gas particles, and both sink particles have the same accretion radius. Our neglect of good thermodynamic treatment means that particles being accreted do not have correct energies in any case, so testing if they are bound to the sink particle would be meaningless.

The size of any subdiscs around the stars is expected to be very small indeed, due to tidal interactions. A disc around the primary would be truncated at R_trunc = 4 R_⊙, and the secondary disc at R_trunc = 3 R_⊙ (Pichardo et al. 2005; Garcia et al. 2013). Comparing this with the stellar radii themselves (respectively 3.3 and 2.5 R_⊙; Böhm et al. 2004; Fumel & Böhm 2012) clearly shows that any discs will be small indeed. Our neglect of them may well be justified then. If the stars have similar corotation radii, as suggested by Garcia et al. (2013), then gas can be accreted directly on to the stars from this radius and so similar sink radii would be justified – however it is unknown if this is the case. It is also not known what level of magnetospheric interaction occurs between the stars, and this could strongly influence how the accretion proceeds.

5.2 Interpretation and observability

We have found that circumbinary discs precessing around a strongly asymmetric binary will undergo changing accretion patterns on time-scales tied to the precession. In the case of HD 104237, the precession period is approximately 40 years, and so we expect that if it is observable, these changes will be seen on a time-scale about half of this, around 20 years. In Fig. 8 we made a very simple model of how the accretion luminosity at pericentre will evolve due to this precession effect.

The disc precession effect will be testable within a few years, given the high rate found in our simulations. In 2 years the disc will have precessed by approximately 18°, and by 36° in 4 years, which should be easily detectable with VLTI interferometry. The strong decentring of the cavity with respect to the binary should also be detectable. It may also be possible to determine the current state of the system by monitoring the Brγ emission of each component over one binary period to determine which component is currently accreting at a higher rate – however, this would require consistent monitoring over decades to give a high level of certainty about where in the precession-accretion cycle HD 104237 currently sits. A more detailed comparison of predictions of these simulations in the context of the HD104237 system will be conducted in a forthcoming paper (Dougados et al., in preparation).

It is possible for this mechanism to play a key role in setting the mass ratio for eccentric binaries. Instead of primarily accreting on to the secondary, gas is able to accrete equally on to both components and maintain the mass ratio status-quo rather than driving the mass ratio towards unity, as is often seen in simulations of binary accretion (e.g. Artymowicz 1983; Bate 2000; de Val-Borro et al. 2011; Bate 2014). Indeed, in our simulation we find that the ratio of precession-period averaged accretion rates between the primary and secondary has values between 0.9 and 1.3, showing that in fact the primary can in fact accrete a little more gas using this mechanism than the secondary can.

Bate (2014) found that with the inclusion of opacity effects in simulations of molecular cloud fragmentation it is possible to reproduce the observed mass ratio distribution well. Importantly, in these simulations there is a bias towards equal-mass ratios in low eccentricity systems, and higher eccentricities correlate with more unequal masses. This is attributed in that paper to dissipation caused by circumbinary discs damping the eccentricity of equal-mass systems, and given the resolution constraints of such large simulations this seems likely to be the case.

More noteworthy for our result is the work of Halbwachs et al. (2003), who found in a survey of solar type (F − K) spectroscopic binaries that high-mass ratio systems (q_b ≥ 0.8) have predominantly lower eccentricities than binaries of unequal masses. This is tentative evidence that the components of eccentric binaries accrete at comparable rates over long time-scales, as suggested by our simulations, allowing them to preserve their non-unity mass ratios. Further work is required in quantifying both how widely this effect operates in terms of binary mass ratio and eccentricity, and how significant the eccentricity differences are between observed equal- and unequal-mass binary systems.

This phenomenon has applications for SMBH binaries in galactic centres, and indeed to date such circumbinary precession has only been explored in this context (e.g. MacFadyen & Milosavljević 2008; Shi et al. 2012; D'Orazio et al. 2013; Farris et al. 2014). However, given that the time-scales in these systems are of the order of tens to hundreds of Myr, it clearly is not going to be directly observable. It is still possible though that this mechanism helps maintain unequal masses for such binaries in the same manner as for stellar binaries. This is especially relevant as we expect these binaries to grow large eccentricities driven by disc interactions (e ∼ 0.6; Roedig et al. 2011). Therefore, it is worth considering this effect when modelling observable signatures of candidate SMBH binaries rather than assuming that the majority of the accretion is on to the secondary as is standard (e.g. Tanaka 2013).

6 CONCLUSIONS

We have presented the results of SPH simulations of the circumbinary accretion disc of the Herbig Ae binary star HD 104237 (DX Cha). We have found that the eccentric, unequal-mass binary causes the disc to become significantly eccentric at the inner edge and precesses uniformly with a rate |$\dot{\varpi } = 0{^{\circ}_{.}} 48\,T_{\mathrm{b}}^{-1}$|⁠, giving a precession period of approximately 40 years. This eccentric shape to the cavity is consistent with observations that the cavity is not centred on the binary centre of mass, and the accretion pattern on single-orbit time-scales is also consistent with the observational evidence (Garcia et al. 2013). We find that this accretion pattern changes with the binary precession, as the position angle of the cavity with respect to the binary orbit changes which binary component accretes at the highest rate.

We model the changing accretion signature on to the binary as an accretion luminosity and find that, if the signal is observable, it should be seen on time-scales of 20 years (half the precession period). Although the change in luminosity is less than the uncertainty in current measurements of the stellar luminosities, if these uncertainties are systematic rather than intrinsic then the signal should still be observable.

We predict that this process will occur whenever a binary of sufficient eccentricity and extreme mass ratio is accompanied by an accretion disc. We find excellent agreement between the precession rates in our simulations and those predicted by Leung & Lee (2013), showing that their prescription can be used to predict the long-term accretion variability time-scale for any such system. While we are unable to provide limits on what values for these parameters are necessary for the process to occur, we suggest that the observed trend towards increased eccentricity for unequal-mass binaries is evidence of this process at work (Halbwachs et al. 2003).

We have used splash (Price 2007) for the SPH visualization in Figs 1 to 4. We would like to thank the anonymous referee for very helpful comments that helped us to improve the paper. We acknowledge support from CONICYT-Chile through ALMA-CONICYT (311200007), FONDECYT (1141175), Basal (PFB0609) and Anillo (ACT1101) grants and from the Millennium Science Initiative (Chilean Ministry of Economy), through grant ‘Nucleus RC130007’. ACD also acknowledges support from the Science & Technology Facilities Council (STFC) in the form of a PhD studentship.

The calculations were run on the DiRAC Complexity system, operated by the University of Leicester IT Services, which forms part of the STFC DiRAC HPC Facility (www.dirac.ac.uk). This equipment is funded by BIS National E-Infrastructure capital grant ST/K000373/1 and STFC DiRAC Operations grant ST/K0003259/1. DiRAC is part of the UK National E-Infrastructure.

REFERENCES