A 3D Hydrodynamics Study of Gravitational Instabilities in a Young Circumbinary Disc

TL;DR

This study uses 3D hydrodynamics simulations to analyze gravitational instabilities in a young circumbinary disc, revealing how a brown dwarf companion influences disc dynamics and evolution.

Contribution

It provides the first detailed comparison of GIs in circumbinary versus single-star protoplanetary discs using radiative 3D hydrodynamics simulations.

Findings

Binarity causes more violent initial GIs due to a strong one-armed density wave.

Both discs reach similar quasi-steady states after 2,500 years.

Binarity affects gravitational torques, temperature profiles, and mass transport, influencing planet formation.

Abstract

We present a 3D hydrodynamics study of gravitational instabilities (GIs) in a 0.14 M$_{\odot}$ circumbinary protoplanetary disc orbiting a 1 M$_{\odot}$ star and a 0.02 M$_{\odot}$ brown dwarf companion. We examine the thermodynamical state of the disc and determine the strengths of GI-induced density waves, nonaxisymmetric density structures, mass inflow and outflow, and gravitational torques. Results are compared with a parallel simulation of a protoplanetary disc without the brown dwarf binary companion. Simulations are performed using CHYMERA, a radiative 3D hydrodynamics code. The onset of GIs in the circumbinary disc is much more violent due to the stimulation of a strong one-armed density wave by the brown dwarf. Despite this early difference, detailed analyses show that both discs relax to a very similar quasi-steady phase by 2,500 years after the beginning of the simulations. Similarities include the thermodynamics of the quasi-steady phase, the final surface density distribution, radial mass influx, and nonaxisymmetric power and torques for spiral arm multiplicities of two or more. Effects of binarity in the disc are evident in gravitational torque profiles, temperature profiles in the inner discs, and radial mass transport. After 3,800 years, the semimajor axis of the binary decreases by about one percentage and the eccentricity roughly doubles. The mass transport in the outer circumbinary disc associated with the one-armed wave may influence planet formation.

Figures

Equations

$$Q=\frac{\kappa {\mathrm{c}}_{\mathrm{s}}}{\pi G\Sigma },$$

$$n(a)=n_o{a}^{-s}$$

$$\Sigma (r)={\Sigma }_o{r}^{-p}$$

$$A_m=\frac{({a_m}^2+{b_m}^2{)}^{\frac{1}{2}}}{\pi \int {\rho }_0\mathrm{rdrdz}},$$

$$a_m=\int \rho \cos (m\varphi )\mathrm{rdrdzd}\varphi ,$$

$$b_m=\int \rho \sin (m\varphi )\mathrm{rdrdzd}\varphi .$$

Full text

††pubyear: 2018††pagerange: A 3D Hydrodynamics Study of Gravitational Instabilities in a Young Circumbinary Disc–LABEL:lastpage

1 Introduction

A circumbinary disc is a disc around two stars orbiting the common center of mass, and a circumbinary exoplanet is a planet orbiting both stars in a binary pair. Circumbinary exoplanets presumably form in circumbinary discs. Over 3,800 exoplanets have been discovered as of November of 2018111NASA Exoplanet Archive, http://exoplanetarchive.ipac.caltech.edu, of which 21 are circumbinary exoplanets. See Haghighipour (2010) for a comprehensive review of planets in binary stars.

The question we address in this paper is how the presence of a central binary may affect the evolution of a circumbinary protoplanetary disc in the case where the disc is massive and subject to gravitational instabilities (GIs). We use a 3D radiative hydrodynamics code with realistic dust opacities to produce two simulations for the same gas disc. In one simulation, the disc orbits a single ${\operatorname{solar-type}}$ star. In the other, the ${\operatorname{solar-type}}$ star is given a relatively close brown dwarf companion. We then compare the onset and subsequent development of GIs.

GIs have long been recognized as a mechanism for producing angular momentum and mass transport in astrophysical discs (Kuiper, 1951; Toomre, 1964; Goldreich & Lynden-Bell, 1965; Pringle, 1981). GIs probably occur during the early phases of protostellar disc evolution when the disc is likely to be massive (Adams & Lin, 1993). For comprehensive reviews of GIs in protoplanetary discs, please refer to Durisen et al. (2007) and Kratter & Lodato (2016).

The Toomre $Q$ parameter

[TABLE]

where ${\rm c_{s}}$ is the sound speed, $\kappa$ is the epicyclic frequency, and $\Sigma$ is the disc surface mass density, is often used to measure a disc’s susceptibility to GIs (Toomre, 1964). For $Q$ $\lesssim$ 1.7, small density perturbations can grow exponentially due to GIs in the linear regime on a timescale comparable to the local dynamical time (Laughlin et al., 1998; Nelson et al., 1998; Pickett et al., 1998). The precise $Q$ stability limit for nonaxisymmetric modes depends in detail on the structure of the disc (Durisen et al., 2007, and references therein). In a ${\operatorname{GI-active}}$ disc, the density perturbations develop into trailing spiral waves over a wide range of radii (Papaloizou & Savonije, 1991; Pickett et al., 1998, 2000, 2003; Laughlin et al., 1998; Nelson et al., 1998; Nelson et al., 2000b). The spiral waves serve to heat the disc, thus increasing $Q$ and decreasing the disc’s susceptibility to GIs. At the same time, radiative cooling acts in the opposite direction, driving the disc towards instability. Sufficiently rapid cooling of the disc may lead to fragmentation (Boss, 2001; Gammie, 2001; Boss, 2002; Johnson & Gammie, 2003; Rice et al., 2005; Boss, 2006; Baehr et al., 2017). When cooling is not too rapid, a balance can be reached between heating and cooling, resulting in a relatively constant, but unstable, value of $Q$. When this occurs, the disc exists in a ${\operatorname{long-lived}}$, ${\operatorname{quasi-steady}}$, marginally unstable state, and the GI activity leads to sustained mass and angular momentum transport (Goldreich & Lynden-Bell, 1965; Paczynski, 1978; Lin & Pringle, 1987; Gammie, 2001; Lodato & Rice, 2004; Mejía et al., 2005; Boley et al., 2006; Durisen et al., 2007; Michael et al., 2012; Steiman-Cameron et al., 2013).

This paper investigates the extent to which the presence of an inner brown dwarf companion affects GIs, especially once the disc has settled into a ${\operatorname{quasi-steady}}$ marginally stable state. Preliminary results of these simulations were reported in Cai et al. (2013). Numerical methods employed in this study and the simulations performed are described in Section 2. Major results, discussion, and conclusion are presented in Section 3, Section 4, and Section 5, respectively.

2 Methods

2.1 3D Hydrodynamics Code

We use the ${\operatorname{grid-based}}$ 3D hydrodynamics code CHYMERA (Computational HYdrodynamics with MultiplE Radiation Algorithms) to perform the simulations. CHYMERA is an Eulerian, 2nd order accurate code both in space and time that solves the Poisson equation and the equations of hydrodynamics on a uniform cylindrical $(r,\phi,z)$ grid, including an energy equation with PdV work, net heating or cooling due to radiative flux divergence, and heating by artificial bulk viscosity. The code assumes mirror symmetry about the equatorial plane or the disc midplane. CHYMERA evolved over the years through contributions from many authors (Tohline, 1980; Yang, 1992; Pickett, 1995; Pickett et al., 2003; Mejía, 2004; Mejía et al., 2005; Cai, 2006; Boley, 2007; Michael, 2011).

In this work, the radiative cooling scheme of Boley et al. (2007b) is adopted. This uses ${\operatorname{flux-limited}}$ diffusion for optically thick regions in the $r$ and $\phi-$directions and a ${\operatorname{single-ray}}$ ${\operatorname{discrete-ordinate}}$ radiative transfer solver in the $z-$direction that treats both optically thick and thin regions. The opacity tables are obtained from D’Alessio et al. (2001). To determine the opacities, minimum and maximum grain sizes of 0.005 $\mu$m and 1 $\mu$m are used, respectively. The dust size distribution follows a

[TABLE]

${\operatorname{power-law}}$ relation, where $n_{o}$ is the normalization constant, $a$ is the grain radius, and $s$ has the value 3.5 for the current work. The ideal gas equation of state includes the internal degrees of freedom for $H_{2}$ (Boley et al., 2007a), and the ortho/para ratio of hydrogen is 3:1.

2.2 Initial Models and Simulations

The simulations follow the evolution of a 0.14 M*⊙* disc surrounding a central 1 M*⊙* star. The initial surface density profile follows the

[TABLE]

${\operatorname{power-law}}$ relation, with $p$ = $0.5$. We use the method described in Mejía et al. (2005) to produce the initial physical and thermodynamical conditions of the disc. The simulation of the protoplanetary disc around a single star will be called the “control” simulation and was already reported by Michael et al. (2011). The “circumbinary” simulation is identical except for the presence of a 0.02 M*⊙* brown dwarf companion initially orbiting the central star at a radius of 2.5 AU. The simulations reported here are computed using (512,512,64) grid elements in the ($r$,$\phi$,$z$) directions, where the $z-$axis is the rotation axis. Each ${\operatorname{grid-cell}}$ is $\approx\frac{1}{6}$ AU in both $r$ and $z$ directions. The inner boundary of the initial disc is at 5 AU, thus creating an inner hole in the disc with no mass. The outer radius of the initial disc is 40 AU. Outflow boundary conditions are enforced at the upper vertical grid boundary, the outer radial grid boundary, and an inner radial boundary at 2 AU. An initial $0.01\%$ random ${\operatorname{cell-to-cell}}$ perturbation is added to the density distribution to generate nonaxisymmetry. Both the control and the circumbinary cases are simulated for 21 ORPs $\approx$ 3,800 years, where 1 ORP = the initial outer rotation period at 33 AU ($\approx$ 179 years).

While the star remains fixed at the grid center for computational convenience, we account for the star’s acceleration through the indirect potential method (Adams et al., 1989; Nelson et al., 2000a), as discussed in Michael & Durisen (2010). In this scenario, the reference frame of the star plus grid is accelerated through inclusion of fictitious forces that account for the gravitational interactions between the star and the disc and between the star and the brown dwarf. The orbital motion of the brown dwarf companion is integrated by a Verlet method (Hut et al., 1995). Gravitational interactions of the brown dwarf with the disc and star are fully included, and the acceleration of the star by the brown dwarf is also treated by an indirect potential, as in Michael et al. (2011).

The masses of the star and brown dwarf and the orbital distance are adopted to mimic the HD 202206 system (Correia et al., 2005), where the central 1.1 M*⊙* star has two companions at 0.83 AU and 2.55 AU with lower mass limits of 17 M*♃* and 2.4 M*♃, respectively (where M♃* = Jupiter mass). Our choice of brown dwarf mass 0.02 M*⊙* is approximately 21 M*♃*, which roughly resembles that of the more massive companion of the HD 202206 system.

3 Results

Fig. 1 shows the midplane and the meridional density contours for both simulations over a range of times. Similarities and differences between the structure and physical states of the simulations, both for early and late stages, are presented in the subsections below.

3.1 Early Evolution

The early evolution, although a transient effect, is very different in the two discs. Differences can be observed for approximately the first 12 ORPs (about 2,150 years), but are dramatic over the first few ORPs. When compared with the circumbinary simulation, the control simulation has a relatively slow onset of ${\operatorname{GI-activity}}$. In the circumbinary disc, the introduction of the brown dwarf companion changes the central gravitational potential of the system, introducing a relatively strong initial $m=1$ perturbation. Therefore, the onset of GIs is hastened. The ${\operatorname{GI-activity}}$ then rapidly spreads throughout the entire circumbinary disc. The spread begins with the growth of a strong ${\operatorname{one-armed}}$ wave. Fig. 2 illustrates the outward propagation and growth of the $m=1$ wave in the circumbinary disc at early times and the appearance of other strong ${\operatorname{low-order}}$ spiral waves.

By about 5 ORPs, as shown in Figs 1 and 2, the entire circumbinary disc is ${\operatorname{GI-active}}$. Contrarily, the control disc becomes ${\operatorname{GI-active}}$ by about seven ORPs, and, as shown in Fig. 3 (bottom panel), the burst of GIs begins first in ${\operatorname{high-order}}$ disturbances, as can be seen by $m=4$ structure growing exponentially before any other.

Eventually, strong ${\operatorname{few-armed}}$ spiral structures develop through the entire disc in both cases. This initial growth is referred to as the “burst” phase (Mejía et al., 2005), a transient that is sensitive to the initial conditions. This phase lasts for a few to several ORPs. After it ends, the disc takes about 1,200 years to settle into ${\operatorname{quasi-steady}}$ behavior. The burst is a period of heating and rapid rearrangement of material in the disc. This onset for GIs in simulations of marginally gravitationally unstable discs around a central star has been described in Mejía et al. (2005), Cai et al. (2006), Boley et al. (2006), and Michael et al. (2012).

This notable difference in behavior of the control and circumbinary simulations at early times can be quantified through a Fourier decomposition of the density structure in $\sin(m\phi)$ and $\cos(m\phi)$, where $\phi$ is the azimuthal coordinate and $m$ represents the number of arms for cases where coherent spiral waves are produced. We can define a single global parameter $A_{m}$ that characterizes the relative amount of the disc’s mass involved in disturbances with $m-$armed symmetry by

[TABLE]

where $\rho_{0}$ is the axisymmetric component of the density, and

[TABLE]

[TABLE]

Fig. 3 illustrates the Fourier amplitudes, which show significant differences in the early evolution described above. In addition to the abrupt onset of GIs in the circumbinary case, the bottom panel of Fig. 3 shows that $A_{1}$ tends to be greater than or comparable to the other ${A_{m}}$’s throughout the simulation. For the control case, on the other hand, $m=2$ tends to be the dominant Fourier component throughout the nonlinear regime, while $m=1$ is relatively unimportant. At later times, the simulations have more similar distributions of amplitudes, except for $m=1$.

Fig. 4, which shows the evolution of the Toomre $Q$ parameter in the two simulations, illustrates another effect of the different ways that GIs set in. By 5$-$6 ORPs, the circumbinary disc has already expanded its radius considerably and been heated by strong spiral waves. The control disc, on the other hand, is still cooling and in fact even shrinking radially at this time, because strong GI waves have not yet erupted.

The differences between the simulations at early times suggest that initiation of ${\operatorname{GI-activity}}$ in a circumbinary disc could be relatively violent and be initially dominated by ${\operatorname{one-armed}}$ structure. It is difficult to know, however, how this might apply to real systems where the discs may transition more slowly to a marginally unstable state. Despite the early violence of the ${\operatorname{GI-onset}}$, there is no hint of fragmentation, and the circumbinary disc ultimately relaxes into a ${\operatorname{quasi-steady}}$ state similar to that seen in the control simulation.

3.2 Quasi$-$Steady State

In the earlier works of Pickett et al. (2003), Mejía et al. (2005), Boley et al. (2006), Cai et al. (2006), Michael (2011), and Steiman-Cameron et al. (2013), after about 10 to 12 ORPs, the discs asymptotically approach a ${\operatorname{quasi-steady}}$ state where heating by ${\operatorname{GI-activity}}$ is roughly balanced by radiative cooling. This state is characterized by relatively constant, marginally unstable values of $Q$, relatively steady amplitudes of the nonaxisymmetric spiral wave structure, significant torques within the disc dominated by a few ${\operatorname{low-order}}$ global spiral modes, and a resulting significant radial mass flow in the disc due to outward angular momentum transport. As shown in Figs 3 and 4, there are significant fluctuations about this average state both in time and space.

These results show that, apart from some differences in dominant spiral waves, especially due to $m=1$, the circumbinary disc achieves a ${\operatorname{quasi-steady}}$ state similar to that of the control disc, even though the onset of GIs in the first 10 ORPs is significantly different. The next few sections lay this out in quantitative detail.

3.2.1 Thermodynamic State

Fig. 4 illustrates that during the time interval of 5 to 6 ORPs, the $Q$ values are lower in the control disc from about $25-35$ AU. At these times, the control disc is still fairly axisymmetric, as shown in Figs 1 and 3. The control disc undergoes several ORPs of initial cooling before the spiral waves reach amplitudes above 10 percent. As mentioned in 3.1, the disc does contract slightly in radius, from 40 to 35 AU as it cools. The circumbinary disc, on the other hand, has already erupted into GI activity and expanded radially due to the passage of strong ${\operatorname{low-order}}$ spiral waves.

Differences persist in both discs for different radial regions of the discs when we compare later times of $10-11$ ORPs and $15-16$ ORPs. For example, in the region around $13-15$ AU, the circumbinary disc has a smooth decreasing $Q$ profile, whereas the $Q$ profile in the control disc has a local peak in $Q$ in the same region. A local peak for $10-11$ ORPs in the radial region around 30 AU is seen in both discs.

Discounting these peaks during $10-11$ ORPs, both discs appear to be approaching the ${\operatorname{quasi-steady}}$ phase by 10 ORPs given the fluctuations typical of the ${\operatorname{quasi-steady}}$ state. Although the outer control disc does have a significantly higher initial $Q$ profile, it steadily decreases within few ORPs as well.

Once the discs have settled, $Q$ values fluctuate about a relatively constant $Q$ $\approx$ 1.5 to 1.7 over the radial range $15-45$ AU, which is consistent with the discs being in a ${\operatorname{quasi-steady}}$ state of marginal instability. Differences in $Q$ at a given radius between the simulated discs are not significant starting around 14 ORPs, and the $Q$ profiles after 14 ORPs are characteristic of the asymptotic ${\operatorname{quasi-steady}}$ state for both simulations. Therefore, we consider both discs to have reached a ${\operatorname{quasi-steady}}$ state of balanced heating and cooling for the duration of 14 to 21 ORPs. It is then possible to use averages over these times to compare the asymptotic states.

Fig. 5 compares the temperature profiles in both discs at the end of simulations. Temperatures in the discs are very similar except for the inner regions, where the circumbinary disc is considerably hotter. The presence of the brown dwarf appears to produce additional heating in the inner disc. There is some material near the brown dwarf that gets rather hot, but it represents a vanishingly small amount of mass.

3.2.2 Nonaxisymmetric Structure

As discussed in Section 3.1, the Fourier $A_{m}$ values provide a global measure of the strengths of nonaxisymmetric density structures produced by GIs. Fig. 6 shows the ${\operatorname{time-averaged}}$ Fourier amplitudes $\langle A_{m}\rangle$ over the time interval $t=$ 14$-$21 ORPs, calculated over the radial region 10$-$40 AU, as a function of $m$. While $A_{m}$ values give a measure of the strengths of the nonaxisymmetric structures, these should not be taken as necessarily representing the strengths of individual, independent spiral waves with $m$ arms. For nonlinear one, two, and ${\operatorname{three-armed}}$ waves, power is expected at all values of $m$. Therefore, an individual $m$ should not necessarily be considered as representing modes (see discussion in Michael et al., 2012).

Figs 3 and 6 reveal noticeable differences in the amplitudes of structures with different $m$-symmetry, even in the ${\operatorname{quasi-steady}}$ state. The amplitude of the ${\operatorname{one-armed}}$ spiral of the circumbinary disc is higher than that of the control disc. On the other hand, after $t=7$ ORPs, the GIs in the control disc are dominated by $m=2$ over most of the simulation. As shown in Michael (2011), the $m=2$ Fourier component in the control disc is produced by a global two$-$armed spiral mode with corotation near 25 AU. The stronger $m=1$ component in the circumbinary case produces somewhat stronger $m=3$ and slightly lower $m=2$ and $4$ in the asymptotic ${\operatorname{quasi-steady}}$ phase, as can be seen in Figs 1, 3 and in the inset of Fig. 6. Nevertheless, these structural differences do not much affect the ${\operatorname{late-time}}$ $Q$ distributions shown in Fig. 4.

3.2.3 Gravitational Torques

The angular momentum and mass in this ${\operatorname{GI-active}}$ disc are transported mainly through gravitational stresses, as shown in Boley et al. (2006), Michael et al. (2012), and Steiman-Cameron et al. (2013). These papers explain in detail how we calculate the internal gravitational torques in the disc for different $m$-values by a Fourier decomposition, and they demonstrate that essentially all the mass transport in the ${\operatorname{quasi-steady}}$ phase is accounted for by these torques alone. In simulations by other researchers, local shearing disc simulations tend to yield comparable gravitational and Reynolds stresses (e.g., Gammie, 2001; Johnson & Gammie, 2003), while measured contributions from Reynolds stresses are small in other global simulations (e.g., Lodato & Rice, 2004).

Fig. 7 presents the ${\operatorname{time-averaged}}$ gravitational torques for several ${\operatorname{low-order}}$ $m$-values along with the total gravitational torques arising from the sum of all the Fourier terms. The quantity plotted represents the torque exerted by the entire disc interior to $r$ on the entire disc exterior to $r$. It is not a local torque but rather the gravitational stress integrated over a cylinder at $r$. The net torque on an annulus of the disc is the difference between the plotted values at the inner edge of the annulus and the outer edge of the annulus. Therefore, it is the slope of the torque curve that determines whether local ring annuli are losing (positive slope) or gaining (negative slope) angular momentum.

Torques shown in Fig. 7 are ${\operatorname{time-averaged}}$ over 14 to 21 ORPs (a time interval of 1,300 years) because the torques fluctuate by about ten percent on orbital time scales. This should not be surprising given the fluctuations evident in $A_{m}$ values shown in Fig. 3. Torques due to ${\operatorname{two-armed}}$ spirals are slightly larger in the control disc while torques due to ${\operatorname{three-armed}}$ spirals are slightly larger in the circumbinary disc. The circumbinary disc has larger torques for $m=1$ structures. The presence of the binary is clearly responsible for the increase. When the torques are summed over all Fourier components excepting $m=1$, it is found that total torques in both discs are very similar.

The primary result in Fig. 7 is that, averaged over the entire ${\operatorname{quasi-steady}}$ phase, the torques due to structures with $m=2$, $3$ and $4$ are remarkably similar for the two simulations, with small shifts in relative strengths of $m=2$ and $m=3$ corresponding to slight differences in $\langle A_{m}\rangle$ values in Fig. 6. The significant effect of the brown dwarf’s presence at this late phase is a strong $m=1$ torque that permeates the disc and actually introduces negative torques over $\sim$ $15-30$ AU that slightly reduce the peak torque there. At $r\gtrsim 30$ AU, the $m=1$ torque component has a broad positive peak that produces a secondary peak in the total gravitational torque. Over time, this could produce substructure in the radial mass distribution, like an annulus of enhanced surface density in the region around 35 AU. This should occur because a positive slope in Fig. 7 tends to correspond to inward transport of mass, while a negative slope tends to correspond with outward transport of mass. As a result, with the ${\operatorname{double-peak}}$ profile of the binary case, it is possible for mass to accumulate somewhere interior to the outer peak.

The $m=1$ torque is the most obvious qualitative difference caused by the binary companion in the quasi-steady phase and is fairly strong even after almost 4,000 years. Given that the $m=1$ power declines slowly but steadily over the circumbinary simulation (Fig. 3), it is unclear how long the $m=1$ torque will influence a circumbinary disc. It could be just a decaying transient due to the initial conditions or it may be a permanent feature, perhaps a forced SLING-like phenomenon (Shu et al., 1990) driven by the companion’s presence. Longer simulations could shed more light on this.

3.2.4 Radial Mass Transport

As shown in the top panel of Fig. 8, both discs settle into very similar azimuthally averaged surface density distributions. The initial ${\operatorname{power-law}}$ surface density distribution with $p=0.5$ is relaxed to $p\sim 1.5$ for the $10-40$ AU region. Even though there is continued mass transport in the ${\operatorname{quasi-steady}}$ phase, most of the density restructuring occurs during the violent ${\operatorname{burst-like}}$ onset of the GIs. It is remarkable that, despite the very different nature of the burst in these two simulations, the final surface density profiles end up quite similar.

The radial mass transport rate during the ${\operatorname{quasi-steady}}$ phase is shown in the bottom panel of Fig. 8. Average mass fluxes are obtained by differencing the total mass inside cylindrical shells at two different times and dividing the result by the time interval. The rates are plotted out to a radial distance of 50 AU, because $\sim$99$\%$ of the disc mass is within the central 50 AU for both discs.

Peak inflow occurs at $\sim 20$ AU for both discs, with mass inflow rates $\approx$ 3 $\times$ 10*-6* M*⊙* per year. The mass flux transitions from inflow to outflow around 23 AU in both discs, with inflow interior to this radius and outflow exterior to this radius. The net result is disc spreading. Spreading is a continuous process as the discs evolve. Very active mass transport is observed in the region between 10 and 40 AU, and the gravitational torque values for this region are also higher than elsewhere. About $80\%$ of the gravitational torque contribution in this region is from two, three and four$-$armed spirals in both discs (see Fig. 7), demonstrating that these structures, as can also be seen in Fig. 1, are largely responsible for the mass transport.

Overall, the mass flux profiles for the two simulations in the ${\operatorname{quasi-steady}}$ phase are quite similar. The one noticeable difference is at radii around 40 AU. A comparison of Fig. 7 and Fig. 8 shows that the difference in mass flow is associated with the secondary peak in the total gravitational torque in the circumbinary simulation arising from the $m=1$ torque component. If this difference persists for $\sim$ 104 to 105 years or more, the circumbinary disc would develop an annular region of enhanced surface density around 40 AU that will not be present in the control disc. This accumulation of mass would be due to the secondary torque peak, as discussed in Section 3.2.3.

3.3 Binary Orbit Evolution

Fig. 9 shows the ${\operatorname{time-evolution}}$ of the orbital eccentricity and the semimajor axis of the binary relative orbit. The average period for all the orbits of the binary is 0.021 ORPs. The local Keplerian period at the initial semimajor axis is 0.022 ORP. The eccentricity of the binary roughly doubles over the course of the simulation from the initial value of 0.019 to the final value of 0.040. Its average rate of increase is $\sim$ 0.001 per ORP $\approx$ 6 x 10*-6* per year while oscillating around that trend. If this growth were sustained for over 105 years, the eccentricity of the system could approach unity. The initial semimajor axis of the binary is 2.453 AU and the final semimajor axis is 2.426 AU. Over the 3,800 year duration of the simulation, the semimajor axis thus decreases by about one percent. Again this suggests significant orbital evolution through interaction with the disc if the simulation were extended to 105 years.

3.4 Inner Disc Vortex

Our analysis of midplane density structures reveals a persistent ${\operatorname{over-density}}$ in the inner regions of both discs. These ${\operatorname{banana-shaped}}$ features appear to be vortices. In fact, the velocities near the density maxima show anticyclonic motion, as expected for a ${\operatorname{long-lived}}$ vortex in a nearly Keplerian disc (Adams & Watkins, 1995). The primary focus of this work is to understand how the central binary affects the evolution of the gravitational instabilities in the bulk of the protoplanetary disc away from the inner edge. However, we want to acknowledge an indication of vortex formation in ${\operatorname{GI-active}}$ discs because it could be of interest to the community. Because of limitations in the vertical resolution of the inner disc in our simulations, we cannot rule out the possibility that these vortices are numerical artifacts (Boley et al., 2007b), but they do behave the way a vortex is expected to behave.

The red oval$-$shaped feature in Fig. 10 is the vortex in the circumbinary disc at 21 ORPs. The density structure of the vortex in the control disc is very similar. Both vortices have centers located at radii of $\sim$ 6 to 7 AU and are ${\operatorname{long-lived}}$. They appear at 4 ORPs and persist through the end of the simulations at a roughly constant density contrast with their surroundings. Their longevity and their constant density structure indicate that the vortices are stable features in the inner disc.

The vortices do vary slightly in strength throughout the simulations. The vortex observed in the control disc is not as robust as the one observed in the circumbinary disc and seems to be dissipating by 21 ORPs. The orbital periods of the vortex centers are 0.09 ORP, essentially the same as the disc rotation period at these radii. The mass of the vortex in Fig. 10 is estimated to be $1\pm 0.5$ M*♃*. The precise mass of the vortex is difficult to determine due to uncertainties inherent in defining the vortex boundaries.

If physical and not numerical in origin, the vortices in both discs are probably due to a Rossby wave instability (RWI) caused by inner edges of discs (Lovelace et al., 1999; Li et al., 2000, 2001). RWI has been shown to be responsible for the formation of vortices in protoplanetary discs (Lovelace & Romanova, 2014, and references therein). In this paper, we only present some relevant details of the vortex. Proper treatment of the vortex and a detailed analysis of the disc’s susceptibility to the RWI would require a simulation with significantly increased resolution of the inner disc, which goes beyond the scope of this paper.

4 Discussion

The reported simulations are most relevant to the early evolution of protoplanetary discs because the 0.14 M*⊙* disc is relatively massive. This work provides a comparison between GIs in a marginally unstable radiatively cooling disc around a single star and the same disc orbiting a close binary of low mass ratio contained inside the central hole of the disc.

4.1 GIs in both Discs

The initial conditions in both simulations are artificial. There are two reasons why the onset of GIs are different during the first 10 ORPs: (1) ${\operatorname{one-armed}}$ disturbances are much stronger in the circumbinary disc and (2) the initial disc is significantly out of equilibrium once the brown dwarf is included. These two effects are probably not independent. The ${\operatorname{one-armed}}$ spiral is probably induced at least in part by the initial deviation from equilibrium. In hindsight, the mass of the central star in the circumbinary disc simulation should have been reduced to 0.98 M*⊙, which would have conserved the total mass of the binary system as 1 M⊙*, equal to the mass of the primary in the control disc. Having an equal amount of mass in the inner hole in both models would have provided a somewhat better comparison. Still, the initial disc would be out of equilibrium when the central object is split into two masses even if the total mass is conserved.

Both discs are in a similar stage of sustained marginal instability in the ${\operatorname{quasi-steady}}$ phase by 14 ORPs. Strengths of the global nonaxisymmetric density features are also similar in both discs. Subtle differences exist in the nonaxisymmetric power and torques due to $m=2$, 3 and 4 armed structures. There is, however, a strong $m=1$ difference lingering even in the ${\operatorname{quasi-steady}}$ phase. Torques due to $m=1$ are much stronger and have an entirely different radial profile in the circumbinary disc than in the control disc. Our simulations show that a disc like this, with radiative cooling times typically much greater than the local dynamic time throughout the disc (Boley et al., 2006; Steiman-Cameron et al., 2013), does not fragment but seeks out a marginally unstable state of balance between GI heating and radiative cooling regardless of a significant difference in GI onset and perturbations by a central binary of low mass ratio.

The presence of the binary companion leads to excess heating in the inner disc. We do not have enough resolution in the inner disc to be certain, but it is likely that this is due to shocks produced in the inner disc through gravitational perturbations by the brown dwarf. Heating in binary protoplanetary discs can prevent formation of planetesimals (Nelson, 2000; Mayer et al., 2005). Similarly, higher temperatures in the inner region of the circumbinary disc can inhibit the planetesimal formation process.

4.2 Comparison with Other Research

There are two competing effects when the tidal forces due to binarity are strong. This has been studied mostly for discs around the individual stars in a binary. Even in this case, results remain somewhat uncertain. On the one hand, strong spiral arms induced by tidal forces produce strong shocks that heat the disc, which has the effect of raising $Q$ and suppressing disc fragmentation (Nelson, 2000). On the other hand, strong perturbations can lead to fragmentation if a disc cools quickly enough (Boss, 2006; Mayer et al., 2010).

Nelson & Marzari (2016) simulated circumbinary tori and discs meant to resemble the GG Tau system. Their simulations did produce fragmentation, but there was no comparison ${\operatorname{single-star}}$ simulation to judge the effects of binarity.

Our simulations using the rigorously tested radiative cooling scheme of Boley et al. (2006, 2007b) demonstrate that a binary of low mass ratio (0.02) placed inside a disc with a long cooling time does not make a large overall difference to the asymptotic ${\operatorname{quasi-steady}}$ behavior of the GIs. We suspect but cannot prove that this would remain true even if the mass ratio were increased to $\sim$ 0.1. On the other hand, the onset of GIs, if initiated by a burst, can be significantly different with and without a binary. We would expect this difference to be more pronounced as the binary mass ratio were increased. A ${\operatorname{one-armed}}$ spiral produced by a GI burst in a circumbinary disc should be readily visible in both continuum and molecular line emission (Ilee et al., 2011; Douglas et al., 2013; Evans et al., 2015; Evans et al., 2017).

Mutter et al. (2017a, b) performed 2D hydrodynamics simulations of the ${\operatorname{Kepler-16}}$, -34, and -35 systems to investigate the role of the disc’s ${\operatorname{self-gravity}}$. They find that GIs can have a significant impact on the formation and evolution of circumbinary planets in these systems. In their simulations, a pronounced $m=1$ global spiral wave is seen due to the disc ${\operatorname{self-gravity}}$. Their simulated orbits agree with the ${\operatorname{Kepler-16b}}$ system, but the orbits of ${\operatorname{Kepler-34}}$ and -35 systems are difficult to explain.

4.3 Mass Transport

Gravitational torques in both discs are similar in the ${\operatorname{quasi-steady}}$ phase, with peak mass inflow rates of about 3 $\times$ 10*-6* M*⊙* per year. In both cases, when the GIs become nonlinear during the burst phase there is a significant mass redistribution within a few ORPs. The only substantial difference at late times is the additional $m=1$ torque in the outer circumbinary disc, which alters the increased mass flow. If the excess $m=1$ torques persist long enough, $\sim$ 104 to 105 years, they may form an enhanced density ring in the outer disc, which could be a site for planet formation in the outer disc through the accumulation of ${\operatorname{meter-sized}}$ particles at the radial pressure maximum (Haghighipour & Boss, 2003).

4.4 The Inner Vortex and Planet Formation

Higher resolution is required for a better understanding of the inner disc vortices. Because the vortex in the circumbinary disc persists from $\sim$4 ORPs to the end of the run, it could be a favorable site for solids to accumulate and possibly even form planetesimals (Lyra et al., 2009). Barge & Sommeria (1995) suggested that anticyclonic vortices in an accretion disc can aid the concentration of dust particles, which could facilitate the formation of planetesimals. The presence of a binary companion, while not responsible for the origin of the vortex, seems to play some role in maintaining its strength.

Recent ALMA observations have suggested that a vortex in the inner disc can be an ideal site for accumulation of ${\operatorname{millimeter-sized}}$ dust particles (Pinilla et al., 2015; Raettig et al., 2015), which can ultimately grow into planetary cores. Multifrequency observations of HD 142527 have shown crescent shaped regions capable of trapping ${\operatorname{submillimeter-sized}}$ particles in the pressure maxima (Casassus et al., 2015). Thus, a ${\operatorname{GI-active}}$ disc with vortices can be a candidate site for accumulation of dust that can form planetesimals. If the inner vortex is sustained long enough, this could lead to the formation of planetary cores.

4.5 Future Work

Onset of GIs is hastened by the presence of the companion, but no fragmentation due to GIs is observed in either simulation. Because the binary mass ratio in this work is only $0.02$ and the disc is 0.14 M*⊙*, the effects of the companion on the GIs in the disc are not pronounced. More simulations involving a diverse set of the stellar and companion masses are needed to better understand GIs in circumbinary discs. With the advent of ALMA, and with its ability to better resolve spiral structures in Class 0$-$1 discs, the role of GIs in the early evolution of protoplanetary discs will hopefully be better understood. Current observations of such discs are sparse (Kratter & Lodato, 2016).

5 Conclusion

The circumbinary simulation addresses the effects of a binary companion on the evolution of a young protoplanetary disc. A ${\operatorname{low-mass}}$ companion, a brown dwarf in this case, is not able to significantly alter the ${\operatorname{long-term}}$ behavior of GIs in a protoplanetary disc once they have reached a ${\operatorname{quasi-steady}}$ state; but it affects density structures and torques associated with the $m=1$ symmetry. The increase in mass transport in some regions in the disc due to the ${\operatorname{one-armed}}$ wave could lead to the formation of an enhanced density ring, which could be a site of planet formation. The companion seems to stabilize the inner disc vortex. This could enhance the process of ${\operatorname{particle-trapping}}$ and ultimately planet formation. The orbital evolution of the binary is measured. In about 3,800 years, its semimajor axis decreases from the initial 2.453 AU to the final 2.426 AU and its eccentricity increases from the initial 0.019 to the final 0.040. If these trends persist over $\sim$ 105 years, the relative orbit would become significantly smaller and highly eccentric.

Acknowledgments: We thank the anonymous referee for providing useful suggestions that improved the paper. This research was supported in part by Lilly Endowment, Inc., through its support for the Indiana University Pervasive Technology Institute and for the Indiana MAETACyt Initiative. This work was also supported in part by the NASA Origins of Solar Systems grant NNX08AK36G. This material is based upon work supported by the National Science Foundation under Grant No. CNS-0521433 and CNS-0723054.