SPIRAL PATTERNS IN PLANETESIMAL CIRCUMBINARY DISKS

Planet formation scenarios and the observed planetary dynamics in binaries pose a number of theoretical challenges, especially concerning circumbinary planetary systems. The orbital motion of the central binary induces the formation of global structures in the gas–dust protoplanetary disk. In particular, a cavity is formed in the central part of the circumbinary disk; though, the cavity can be penetrated by streams of matter, maintaining the accretion activity of the binary components. Numerical simulations of the gas accretion in such systems were performed using the smoothed particle hydrodynamics (SPH) method (Artymowicz & Lubow 1996; Bate & Bonnell 1997; Sotnikova & Grinin 2007) and finite-difference schemes (Gunther & Kley 2002; Ochi et al. 2005; Hanawa et al. 2010; Kaigorodov et al. 2010).

In a general theoretical context, various spiral wave modes, which may emerge due to tidal interaction with planets in gaseous circumstellar disks, were revealed in Lubow & Ogilvie (1998). An empirical five-parameter fitting relation describing the spiral pattern excited by a planet in a gaseous disk was proposed in Muto et al. (2012). In numerical simulations of gaseous protoplanetary disks in binary systems, as well as in simulations of disks with embedded planets, a variety of spiral patterns were found to form (Nelson 2000; Kley & Nelson 2008; Paardekooper et al. 2008; Marzari et al. 2009, 2012; Picogna & Marzari 2013). Recent ALMA observations provide remarkable examples of spiral structures in gaseous disks excited by embedded planets (Christiaens et al. 2014) or by central stellar binaries (Takakuwa et al. 2014).

Spiral density waves may also emerge in the inner zones of cool accretion disks around components of semidetached binary systems. The possibility of observations of such patterns, obtained in three-dimensional gas-dynamical simulations, was discussed in Bisikalo et al. (2004). The pattern formation in such disks is due to the differential orbital precession of the accreting particles, caused by the perturbations from the donor star.

It is well known that a spiral structure can form if a planet is present in a gas-free planetesimal disk (Wyatt 2005; Mustill & Wyatt 2009; Matthews et al. 2014). Wyatt (2005) showed that a tightly wound spiral starts to propagate through the disk if a planet is introduced inside it. As in the theory of Bisikalo et al. (2004), the spiral is generated by the differential precession of the neighboring particle orbits. To estimate the perturbation effect of a planet on the orbits of the planetesimals, Wyatt (2005) used the theory of secular variation of the complex eccentricity coupled with the variation of the longitude of pericenter using the Lagrange planetary equations. He applied this secular theory to explain a spiral pattern resolved by Clampin et al. (2003) in the image of the HD 141569 disk and to interpret the results of a numerical simulation of the planetesimal disk evolution in the presence of a giant planet. An important effect considered by Mustill & Wyatt (2009) is the dynamical heating of the disk as a result of the spiral wave propagation caused by a planet being introduced into the disk: the planetesimal velocity dispersion increases. This phenomenon, hindering the planetesimal coalescence, may affect the further process of planetary formation. Of course, for the circumbinary case considered below, a similar effect may be of even greater importance.

The protoplanetary disks strongly radiate in the infrared and submillimeter spectral ranges because they contain small-sized dust particles. The fine dust fraction decreases with time due to the particles' coagulation; on the other hand, numerical simulations show that the collisional destruction of planetesimals is also an effective process, causing the dust to persist over millions of years (Birnstiel et al. 2009). Many images of gas-rich, transitional, and debris disks have been obtained (see, e.g., Padgett et al. 1999; Grady et al. 2005; Kalas et al. 2005; Mathews et al. 2012; Mayama et al. 2012; Kennedy et al. 2012). In these images, global structures were detected, namely spiral arms (Grady et al. 2001; Hashimoto et al. 2011; Christiaens et al. 2014; Takakuwa et al. 2014), ring-like gaps (Weinberger et al. 1999), bright rings (Kalas et al. 2005), inner holes (Mathews et al. 2012; Mayama et al. 2012), warps (Heap et al. 2000), and density clamps (Greaves et al. 1998).

The formation of spiral waves in the gaseous disk depends on the viscous gas properties. When the gaseous component disappears, the accretion stops and any further structures are then formed in the planetesimal ensemble. In this paper, we explore the dynamical stirring of a planetesimal circumbinary disk in the epoch when the gas component disappears. We show that the orbital precession of the particles forming the disk generates prominent spiral patterns, which we describe analytically.

In our analysis, we assume that the timescale for the secular wave propagation is greater than the lifetime of the gas-rich disk (and by far greater than the characteristic timescales of gas depletion and planetesimal formation). Whether this assumption is realistic for any particular system depends on the values of a number of parameters, which we identify. We demonstrate that the secular wave propagation timescale can, in principle, be long enough for our assumption to be valid.

We emphasize that in this paper we mostly consider circumbinary disks (where the planetesimals are moving around both stars that are situated at the disk center; thus, it is an outer disk, with respect to the binary). Circumstellar disks, where the planetesimals are moving around one of the binary components (thus, it is an inner disk), are only discussed briefly in Section 5.

The paper is structured as follows. In Section 2, we consider a secular theory for the dynamics of planetesimals in circumbinary disks and derive the necessary formulas for the secular evolution of the planetesimal eccentricity and longitude of pericenter. An analytical expression describing the geometry of the spiral density wave is given here. In Section 3, we provide a formula for the secular timescale of the wave propagation and discuss whether this timescale can exceed the lifetime of a gas-rich disk in realistic systems. In Section 4, the radial (with respect to the central binary) secular "oscillations" of the planetesimal eccentricity are investigated analytically and compared to the known results from numerical simulations; we also consider the time evolution of individual orbits and, again, make comparisons with the known results from numerical simulations. In Section 5, a secular theory for the circumstellar case is briefly discussed. In Section 6, we explore correspondence between the analytical spiral pattern and the numerical–experimental density waves in the circumbinary case. Our conclusions are summarized in Section 7.

In this section, we develop a secular theory for the dynamics of planetesimals in circumbinary disks following theoretical approaches by Heppenheimer (1978) and Moriwaki & Nakagawa (2004). A secular perturbation theory, providing analytical formulas for a particle's eccentricity e and the longitude of periastron ϖ, was derived for the circumstellar (circumprimary or circumsecondary) case by Heppenheimer (1978; see also Whitmire et al. 1998; Thébault et al. 2006). It was used to analytically describe how the circumstellar disk of a young star is stirred by a companion star. On the other hand, the circumbinary case was considered by Moriwaki & Nakagawa (2004). They presented approximate differential equations for the secular evolution of the eccentricity vector; see the Appendix in Moriwaki & Nakagawa (2004). We combine the approaches of Heppenheimer (1978) and Moriwaki & Nakagawa (2004) to derive the necessary explicit analytical formulas for the secular evolution of the planetesimal eccentricity and longitude of pericenter.

Hereafter we adopt a barycentric frame; m₁ and m₂ are the masses of the binary components (we set ${{m}_{1}}\geqslant {{m}_{2}}$), a_b is the binary semimajor axis, e_b is the binary eccentricity, and a is the semimajor axis of the planetesimal orbit. The masses are measured in solar units, distances in astronomical units (AU), and time in years. Thus, the gravitational constant G is equal to $4{{\pi }^{2}}$.

The hierarchical three-body problem, the "binary-planetesimal," is considered: the distance of the planetesimal from the binary barycenter is assumed to be much greater than the size of the binary. In fact, it is superfluous to consider a non-hierarchical case because a large central chaotic circumbinary zone exists at all eccentricities of the planetesimal if $\mu \gtrsim 0.05$ (Shevchenko 2015), where $\mu ={{m}_{2}}/({{m}_{1}}+{{m}_{2}})$ is the mass parameter of the binary. This relative mass threshold has an important physical meaning (Shevchenko 2015): above the threshold, the tertiary (planetesimal) can diffuse, even starting from small eccentricities, following the sequence of the overlapping p:1 mean motion resonances between the tertiary and the binary, up to ejection from the system; close encounters with other bodies are not required for the escape. (Note that, on the other hand, in the case of inner circumstellar orbits, stable orbits always exist inside the Hill spheres of the binary components.) In what follows, we assume that $\mu \gtrsim 0.05$.

We use an averaged perturbing function expansion presented in Moriwaki & Nakagawa (2004) for the circumbinary case. This is a power-law expansion in the ratio of the binary and planetesimal semimajor axes and in the eccentricities. The expansion is up to the third order in the ratio of the semimajor axes and up to the fourth order in the eccentricities, inclusively. Note that the ratio of the binary and planetesimal semimajor axes is assumed to be small, since the hierarchical problem is considered. We integrate analytically the corresponding equations of motion (see Equation (A7) in Moriwaki & Nakagawa 2004), and thus straightforwardly deduce the formulas for the secular evolution of the eccentricity e and the longitude of periastron ϖ of the circumbinary planetesimal. They turn out to be very similar to those in the circumstellar case, presented in Heppenheimer (1978), and read

$$\begin{eqnarray}&&e={{e}_{{\rm max} }}|\;{\rm sin} \frac{ut}{2}|,\end{eqnarray} \tag{ 1 }$$

$$\begin{eqnarray}&&{\rm tan} \varpi =-\frac{{\rm sin} \;ut}{1-{\rm cos} \;ut},\end{eqnarray} \tag{ 2 }$$

where t is time,

$$\begin{eqnarray}&&u=\frac{3\pi }{2}\frac{{{m}_{1}}{{m}_{2}}}{{{({{m}_{1}}+{{m}_{2}})}^{3/2}}}\frac{a_{b}^{2}}{{{a}^{7/2}}}\left( 1+\frac{3}{2}e_{b}^{2} \right),\end{eqnarray} \tag{ 3 }$$

and ${{e}_{{\rm max} }}=2{{e}_{f}}$, where e_f is the so-called forced eccentricity

$$\begin{eqnarray}&&{{e}_{f}}=\frac{5}{4}\frac{({{m}_{1}}-{{m}_{2}})}{({{m}_{1}}+{{m}_{2}})}\frac{{{a}_{b}}}{a}{{e}_{b}}\frac{\left( 1+\frac{3}{4}e_{b}^{2} \right)}{\left( 1+\frac{3}{2}e_{b}^{2} \right)}.\end{eqnarray} \tag{ 4 }$$

Using a new variable, namely, $y=\frac{ut}{2}$, we rewrite Equation (2) in the form

$$\begin{eqnarray}\begin{array}{ccccccccccccccc} {\rm if}\ y\geqslant -\pi \ {\rm and}\ y\leqslant -\frac{\pi }{2},{\rm then}\ \varpi & = & y+5\frac{\pi }{2}; \\ {\rm if}\ y\geqslant -\frac{\pi }{2}\ {\rm and}\ y\leqslant 0,{\rm then}\ \varpi & = & y+\frac{\pi }{2}; \\ {\rm if}\ y\geqslant 0\ {\rm and}\ y\leqslant \frac{\pi }{2},{\rm then}\ \varpi & = & y+3\frac{\pi }{2}; \\ {\rm if}\ y\geqslant \frac{\pi }{2}\ {\rm and}\ y\leqslant \pi ,{\rm then}\ \varpi & = & y-\frac{\pi }{2}. \\ \end{array}\end{eqnarray} \tag{ 5 }$$

Thus, u can be regarded as a "precession rate" (though in a modified fashion) of an individual orbit.

Based on the analytical dependences of the eccentricity and the longitude of periastron on the semimajor axis (Equations (1) and (5)), let us proceed to describe the circumbinary disk structure emerging due to the secular perturbations from the central binary. Namely, let us derive an analytical formula for the spiral pattern in the gas-free case. This can be done by tracing the radial distance r of the apocenter at the maximum eccentricity of the particle's orbit  (in the course of its secular evolution) as a function of the polar angle θ. The apocenters are attained at $\varpi =\pi$ mod $2\pi$, i.e., at the polar angle $\theta =\pi$ mod $2\pi$. At these points the analytical solution is known (Equations (1) and (5)). Taking into account that at $\theta =0$ mod $2\pi$ the spiral corresponds to particle's minimum eccentricity (zero), and interpolating between the consecutive points $\theta =\pi n$ ($n=1,2,3,...$), one arrives at

$$\begin{eqnarray}&&r(\theta ,t)={{\left( \frac{At}{\theta } \right)}^{2/7}}+B(1-{\rm cos} \;\theta ),\end{eqnarray} \tag{ 6 }$$

where the constants

$$\begin{eqnarray}\begin{array}{ccccccccccccccc} A & = & \frac{3\pi }{2}\frac{{{m}_{1}}{{m}_{2}}}{{{({{m}_{1}}+{{m}_{2}})}^{3/2}}}a_{b}^{2}\left( 1+\frac{3}{2}e_{b}^{2} \right), \\ B & = & \frac{5}{4}\frac{{{m}_{1}}-{{m}_{2}}}{{{m}_{1}}+{{m}_{2}}}{{a}_{b}}{{e}_{b}}\frac{\left( 1+\frac{3}{4}e_{b}^{2} \right)}{\left( 1+\frac{3}{2}e_{b}^{2} \right)}. \\ \end{array}\end{eqnarray} \tag{ 7 }$$

Thus, the pattern represents a shifted "power-law spiral." The main term, proportional to ${{\theta }^{-2/7}}$, represents a generalized lituus (a crook). (Note that the classical lituus is a little bit different, being given by $r\propto {{\theta }^{-1/2}}$.) The shifting term $\propto (1-{\rm cos} \;\theta )$ represents the cardioid.

The observability of the derived secular spiral pattern depends on an interplay of several dynamical and physical timescales. As determined from observations of stars in young clusters and associations, the protoplanetary disks start their evolution in a gas-rich phase, which lasts some $\sim {{10}^{6}}$–10⁷ yr (Haisch et al. 2001; Williams & Cieza 2011). Then, in a relatively short-lived ($\sim {{10}^{5}}$ yr) "transitional disk" phase (Cieza et al. 2007; Williams & Cieza 2011), the gas is lost and the pure debris disk is left, consisting of planetesimals, possibly with planets or planetary embryos already embedded. The planetesimals form (by dust accretion) already in the first gas-rich phase on a short timescale of $\sim {{10}^{4}}$ yr (Weidenschilling 1997); the planets or planetary embryos may also form in the first (gas-rich) phase (see, e.g., Kokubo & Ida 1996; Ikoma et al. 2000).

In our analysis, we assume that the timescale of the secular wave propagation across the disk is greater than the lifetime of the gas-rich disk, and is by far greater than the characteristic timescales of gas depletion and planetesimal formation.

From Equations (6) and (7) one can readily find the secular timescale T_s, namely, the time needed for the density wave to propagate across the disk. For this purpose we set $\theta =\pi$ (the opposite value of the singular value $\theta =0$) and define the time $t={{T}_{s}}$ as a function of $r={{r}_{{\rm disk}}}$ (the disk's outer edge radius) as

$$\begin{eqnarray}&&{{T}_{s}}=\frac{\pi }{A}{{({{r}_{{\rm disk}}}-2B)}^{7/2}}\approx \frac{\pi }{A}r_{{\rm disk}}^{7/2}.\end{eqnarray} \tag{ 8 }$$

The approximate equality is valid because the problem is hierarchical.

Thus, the secular timescale depends on a number of parameters and rather strongly on the disk radius, increasing sharply with r_disk and decreasing (rather moderately) with a_b, μ, and m₁. The binary eccentricity e_b is relatively unimportant.

For Kepler-16 (m₁ = 0.69 ${{M}_{\odot }}$, m₂ = 0.2029 ${{M}_{\odot }}$, ${{a}_{b}}=0.2243$ AU, ${{e}_{b}}=0.159$), if one sets ${{r}_{{\rm disk}}}=30$ AU, one has ${{T}_{s}}=1.1\cdot {{10}^{7}}$ yr. For an "M-dwarf–brown dwarf" system (m₁ = 0.5 ${{M}_{\odot }}$, m₂ = 0.025 ${{M}_{\odot }}$, ${{a}_{b}}=1$ AU, ${{e}_{b}}=0.1$), if one sets ${{r}_{{\rm disk}}}=100$ AU, one has ${{T}_{s}}=2.0\cdot {{10}^{8}}$ yr. Therefore, the secular timescale can vary significantly from system to system.

For the secular spiral to be observable, one must assume that the timescale for the secular wave propagation is greater than the lifetime of the gas-rich disk ($\sim {{10}^{6}}$–10⁷ yr; see the Introduction). Whether this assumption is realistic for any particular system depends on the values of the problem parameters: r_disk, a_b, μ, and m₁. As we have just seen, the secular wave propagation timescale can, in principle, be long enough; it may well exceed 10⁷ yr, depending on the problem parameters.

Numerical simulations of the dynamical stirring of planetesimal disks on secular timescales were performed in Moriwaki & Nakagawa (2004), Meschiari (2012a, 2012b), and Paardekooper et al. (2012) in various problem settings. In particular, graphs of the eccentricity and longitude of periastron of a circumbinary particle as a function of its semimajor axis were constructed numerically. Let us show that such graphs can also be constructed analytically using the theory described above.

We consider Kepler-16 as an actual binary example.  This binary is remarkable, hosting the first ever discovered circumbinary planet in a double main sequence star system. The binary's parameters are ${{m}_{1}}=0.69\;{{M}_{\odot }}$, ${{m}_{2}}=0.2026\;{{M}_{\odot }}$, ${{e}_{b}}\,=0.159$, and  ${{a}_{b}}=0.2243$ AU (Doyle et al. 2011). Analytically constructed diagrams, as given by Equations (1) and (5) at $t={{10}^{5}}$ yr, are shown in Figure 1. Radial "waves" of eccentricity, well known from numerical diagrams in Moriwaki & Nakagawa (2004), Meschiari (2012a, 2012b), and Paardekooper et al. (2012), are clearly recognizable.

Zoom In Zoom Out Reset image size

Our diagrams exactly match those constructed in Moriwaki & Nakagawa (2004) and Meschiari (2012b) for models with the same parameters when the gas component is absent. A direct comparison of our diagrams for models with the gas component with those constructed in Meschiari (2012a) and Paardekooper et al. (2012) shows that the presence of gas increases the timescale of eccentricity pumping.

An analogous approach can be used to describe secular oscillations of a circumbinary orbit of a particle (a planetesimal or a planet) with a dependence on time at a fixed semimajor axis. In this respect, one can compare the theory described above with the results of numerical simulations performed in Dunhill & Alexander (2013), Leung & Lee (2013), and Meschiari (2014) for the planet Kepler-16b. We set the binary parameters to be the same as cited above in this section and a = 0.7016 AU is fixed for the planet. The initial conditions for the planetary motion are taken from Doyle et al. (2011).

In Figure 2, the time dependences for the eccentricity and the longitude of periastron, given by Equations (1) and (5), are constructed. As follows from a direct comparison with the corresponding numerical–experimental graphs in Leung & Lee (2013) and Meschiari (2014), our secular theory formulas match these graphs perfectly.

Zoom In Zoom Out Reset image size

As opposed to the circumbinary case, circumstellar disks are defined as those for which planetesimals move around one of the binary components. Thus, the disk is an inner structure, with respect to the binary itself. In this section we briefly consider the circumstellar case.

Planetesimals orbiting a component of a binary star are subject to perturbations due to the stellar companion. The perturbations pump the eccentricities of planetesimals and inhibit their accumulation process (Marzari & Scholl 2000; Thébault et al. 2004). In Thébault et al. (2006), the timescales for the inward propagation of the orbital crossing "wave" are calculated in gas-free models. In the more realistic models, the density wave inward propagation is damped by the gas drag and collisions.

In circumstellar disks, strong spiral arms can be present, subject to retrograde precession (Kley & Nelson 2008; Paardekooper et al. 2008). Oscillations of planetesimal orbital parameters (the eccentricity and the longitude of pericenter) as functions of time at a fixed semimajor axis are presented graphically in Kley & Nelson (2008), Marzari et al. (2009), Beaugé et al. (2010), Marzari et al. (2012), Müller & Kley (2012), and Leiva et al. (2013) for various values of disk parameters (mass, viscosity, radiative transfer properties). Possible mechanisms accounting for the spiral structure in a circumstellar gaseous disk perturbed by an outer body (i.e., planet, star, brown dwarf) were studied in Quillen et al. (2005) using two-dimensional hydrodynamic simulations. Numerically simulated diagrams of the eccentricity and the longitude of periastron of planetesimals as a function of the semimajor axis at a fixed time are presented in Marzari & Scholl (2000), Paardekooper et al. (2008), Paardekooper & Leinhardt (2010), and Leiva et al. (2013). The choice of parameters for the simulations usually corresponds to the binary systems γ Cep and α Cen, possessing circumstellar planets.

What are the basic expressions describing the secular evolution in the circumstellar case? The expressions for the eccentricity and longitude of periastron are virtually the same as in the circumbinary case and are given by Equations (1), (2), and (5), but the formulas for the u and e_f parameters, entering these expressions, are different. We quote them as given in Heppenheimer (1978), Whitmire et al. (1998), and Thébault et al. (2006), but arrange in our notations as

$$\begin{eqnarray}&&u=\frac{3\pi }{2}\frac{{{m}_{2}}}{m_{1}^{1/2}}\frac{{{a}^{3/2}}}{a_{b}^{3}}{{\left( 1-e_{b}^{2} \right)}^{-3/2}},\end{eqnarray} \tag{ 9 }$$

$$\begin{eqnarray}&&{{e}_{f}}=\frac{5}{4}\frac{a}{{{a}_{b}}}\frac{{{e}_{b}}}{\left( 1-e_{b}^{2} \right)}.\end{eqnarray} \tag{ 10 }$$

Here m₁ is the primary mass (around which the particle orbits) and m₂ the perturbing mass; ${{m}_{1}}\gt {{m}_{2}}$. Note that in the treatment by Heppenheimer (1978), the perturbing function is not expanded in the eccentricities, but solely in the ratio of the planetesimal and binary semimajor axes (up to the third order inclusive). The ratio of the planetesimal and binary semimajor axes is assumed to be small, since the hierarchical problem is considered.

Thus, the secular evolution in the circumstellar case is described analytically by Equations (1), (2), (5), (9), and (10).

Let us consider in more detail the formation of patterns in circumbinary disks. As is clear from Figure 1, the formation of secular structures in the disk is inevitable because the maximum eccentricities correspond to the particles' groupings at apocenters due to the low velocities of the particles there. Therefore, the disk radial densities are non-uniform at any fixed time.

Let us provide an illustrative example. Fixing the time of evolution to $t={{10}^{4}}$ yr, we calculate analytically the evolved planetesimal orbits in the semimajor axis a interval from 3${{a}_{b}}$ to $16{{a}_{b}}$ with a distance step ${\Delta }a=0.01{{a}_{b}}$. (The chosen inner radius corresponds to roughly the typical size of the circumbinary chaotic zone; see Shevchenko 2015.) The planetesimal orbits at each a are shown by dots, with the time step equal to 1% of the planetesimal orbital period. The binary parameters are set to be ${{m}_{1}}={{M}_{\odot }}$, ${{m}_{2}}=0.2\;{{M}_{\odot }}$, ${{e}_{b}}=0.4$, and ${{a}_{b}}=1$ AU; this choice roughly corresponds to the "most eccentric" model in Moriwaki & Nakagawa (2004). The resulting diagram is presented in Figure 3(a).

Zoom In Zoom Out Reset image size

We have also performed an SPH code numerical simulation corresponding to this model (i.e., the parameters and initial conditions are the same). The SPH code realized in Sotnikova (1996) has been used. The number of particles was set to 30,000. The kernel for smoothing the hydrodynamic quantities was chosen in the form of a spline (Monaghan & Gingold 1983), and the smoothing length was set to be constant. On starting the simulation, the SPH particles are placed on circular orbits; the velocities are Keplerian. The particles have equal sizes and masses and are distributed in accordance with the surface density profile ${\Sigma }\sim {{r}^{-1}}$. The total mass of the planetesimal disk is assumed to be negligible compared to that of the binary; therefore, the self gravity of the disk is neglected. In the gas-free case, the system of equations in the applied SPH code reduces to the equations of motion of the particles in the gravitational potential of the central binary. The results of the simulations in the gas-free case are presented in Figure 3(b).

In both panels of Figure 3, spiral patterns are evident. They are tightly bound in the direction of the companion motion. Note a close resemblance of the patterns in the panels, though panel (a) is constructed using an analytical theory, whereas panel (b) is a pure simulation.

In a separate figure (Figure 4), we illustrate graphically how well our analytical spiral, given by Equation (6), matches the density wave emerging in the numerical simulations in the gas-free case. Namely, we superimpose the analytical spiral on the numerical–experimental plot taken from Figure 3(a). In Figure 4, the analytical curve is depicted by a thick solid line. A good agreement between the analytical spiral pattern and the numerical density wave is evident.

Zoom In Zoom Out Reset image size

Consider now a planetesimal disk containing gas. Gas drag acting on a planetesimal depends on the gas and planetesimal parameters (Whipple 1973; Weidenschilling 1977). We consider large particles and assume that the free path of gas molecules is less than the planetesimal radius. The planetesimals are subject to gas drag; to calculate the effect of the gas drag, we complement the system of equations in the SPH code with the continuity equation for gas and the ideal gas equation of state in the isothermal case; namely, gas pressure $P={{c}^{2}}{{\rho }_{g}}$, where c is the sound speed and ${{\rho }_{g}}$ is the gas density. To simulate the gas presence in the planetesimal disk, we add an extra acceleration term

$$\begin{eqnarray}&&\frac{d{{{\boldsymbol{v}} }_{p}}}{dt}=-\frac{3}{8}{{C}_{D}}\frac{{{\rho }_{g}}}{{{\rho }_{p}}s}|{{v}_{p}}-{{v}_{g}}|({{{\boldsymbol{v}} }_{p}}-{{{\boldsymbol{v}} }_{g}})\end{eqnarray} \tag{ 11 }$$

in the equations of motion of planetesimals, as was done in Whipple (1973) and Weidenschilling (1977). Here ${{\rho }_{g}}$ and ${{\rho }_{p}}$ are gas and planetesimal densities, s is the radius of a planetesimal, and v_p and v_g are planetesimal and gas velocities, respectively. The dimensionless coefficient ${{C}_{D}}=0.44$, as in Whipple (1973), for a regime with a  Reynolds number $\operatorname{Re}\gt 800$.

In our model the sound speed is 1.2 km s⁻¹, corresponding to a temperature $\approx 100$ K; this gives the disk's effective semithickness $h=H/r=0.06$ at r = 1 AU, where H is the disk's semithickness and r is the distance from the center of mass. (The disk's semithickness is controlled by the gas temperature, namely, $H={{2}^{1/2}}c/$Ω, where Ω is the Keplerian angular velocity.) We describe the gaseous disk's effective viscosity in the α-parameterization (Shakura & Sunyaev 1973; Pringle 1981): $\nu =\alpha cH$, where α is the Shakura–Sunyaev viscosity parameter; $\alpha \approx 0.03$, like in Artymowicz & Lubow (1996). The model's other parameters are as follows: the planetesimal radius s = 1 km, the planetesimal density ${{\rho }_{p}}=1$ g cm⁻³, and the gas total mass in the disk ${{M}_{g}}={{10}^{-2}}$ ${{M}_{\odot }}$.

Our SPH simulations of the gas-containing disk show that the presence of residual gas slows the wave propagation, however, with the wave pattern remaining self-similar. If the presence of the gas is great enough, the spiral pattern does not emerge altogether. Note that, since the work of Marzari & Scholl (2000), it is well known that the local gas damping, if high enough, causes an alignment of the apsidal lines of circumbinary orbits of equal-sized planetesimals. Therefore, increasing the gas amount must suppress the secular spiral structure formation in the planetesimal disk.

We limit our analysis of gas-containing disks with these short considerations. The problem of pattern formation in the two-component (gas plus planetesimals) media deserves a detailed separate study, which we leave for future work. Let us give just one illustrative example (Figure 5). The initial data for the planetesimal disk are the same as were used for the construction of Figure 3(b), but a small amount of gas is introduced (the model parameters are as described above). As follows from the simulations, it requires 10⁴ yr to form the same spiral pattern as emerging in the gas-free case in $8\cdot {{10}^{3}}$ yr. In Figure 5, it is evident that gas amount is insufficient in the given model to suppress the spiral pattern completely. On the other hand, one can see that an elliptical structure is also formed (at the inner edge of the disk), with the major axis perpendicular to the apsidal line of the central binary; this is in accord with Marzari & Scholl (2000) on the alignment of planetesimal orbits due to gas drag.

Zoom In Zoom Out Reset image size

Note that the disk self gravity is ignored in our simulations. The self-gravity effect depends on the masses of planetesimals and on their spatial concentration; thus, we assume that the masses and/or concentration are small enough. However, one may speculate that the self-gravity effect, when non-negligible, generally suppresses the spiral pattern, because, as follows from simulations (Marzari et al. 2009), it leads to low eccentricities of the disk particles.

Concerning the circumstellar case, the "secular" spiral pattern is basically analogous to that derived earlier in this section for the circumbinary case. A detailed analysis for the circumstellar case, along with comparisons with the numerous simulation results cited in Section 5, will be performed elsewhere.

Thus, we have combined the approaches of Heppenheimer (1978) and Moriwaki & Nakagawa (2004) to derive explicit analytical formulas for the secular evolution of the eccentricities and longitudes of pericenters of planetesimals in circumbinary disks when the gas component disappears. We have shown that, if the binary is eccentric and its components have unequal masses, spiral density waves are generated. They engulf the disk on a secular timescale, which may well exceed 10⁷ yr, depending on the problem parameters. Being transient (on the secular timescale), their observed presence may betray a system's young age.

Our analytical results, given by Equations (1)–(5) and Equations (9) and (10), explain qualitatively and quantitatively many numerical–experimental diagrams, such as those presented in Marzari & Scholl (2000), Moriwaki & Nakagawa (2004), Thébault et al. (2006), Paardekooper et al. (2008), Paardekooper & Leinhardt (2010), Meschiari (2012a, 2012b), Paardekooper et al. (2012), and Leiva et al. (2013) for the eccentricity and the longitude of periastron as a function of the semimajor axis at a fixed time, as well those as presented in Marzari et al. (2009), Beaugé et al. (2010), Müller & Kley (2012), Leung & Lee (2013), Leiva et al. (2013), and Meschiari (2014) for the eccentricity and the longitude of periastron as a function of time at a fixed semimajor axis.

For the gas-free case, we have derived an analytical formula for the spiral pattern; this formula describes a modified "lituus" (a shifted power-law spiral). Its form is in perfect agreement with our results from modeling the precessing planetesimal orbits, as well as our results of direct SPH simulations performed for the same disk parameters and initial conditions. Using the SPH scheme, we have observed that the effect of residual gas is to slow down the spiral wave propagation.

In our analysis, we have assumed that the timescale for the secular wave propagation is greater than the lifetime of the gas-rich disk ($\sim {{10}^{6}}$–10⁷ yr). Whether this assumption is realistic for any particular system depends on the values of a number of parameters, which we have identified. We have demonstrated that the secular wave propagation timescale can be, in principle, long enough: it may exceed 10⁷ yr, depending on the problem parameters.

One should underline that the circumbinary spiral pattern emerges if the binary is eccentric and its components have unequal masses. If any of these two conditions are violated, concentric circular waves ("secular circles") are formed instead of the spiral wave, as is evident from the results of numerical simulations presented in Figure 1 in Lines et al. (2014; for the case of an equal-mass binary) and Figure 5 in Rollin et al. (2014; for the case of a zero-eccentricity binary). Transitions (possibly involving thresholds in the values of the problem parameters) from these two limiting "circular" cases to the "spiral" one deserve a separate study.

Another important issue that has not been addressed here and deserves a study is the propagation of the secular spiral wave in the presence of planets (which themselves can generate patterns in the disk). Indeed, modern theories favor scenarios in which planets form in gas-rich environments (see, e.g., the review by Zhou et al. 2012), i.e., the planets should already be present when the secular spiral wave starts to propagate.

Finally, note that the derived one-armed spiral pattern is a secular one; i.e., it does not concern any orbital resonances. On the other hand, the presence of multi-armed spiral patterns in astrophysical disks is usually associated with Lindblad resonances; see Binney & Tremaine (2008) for the case of galactic dynamics and Murray & Dermott (1999) for the case of planetary ring dynamics. Is it possible to observe multi-armed structures in the case of planetesimal circumbinary disks? We think that this is not possible, because all low-order resonances are deep inside the central chaotic zone around the binary (if $\mu \gtrsim 0.05$; see Shevchenko 2015). Conversely, observable multi-armed structures can be generated by resonances in the circumstellar case, i.e., in the inner disks, because the stellar companions are inside their individual stability regions.

We are grateful to the referee for useful remarks. We wish to thank Alice Quillen for useful remarks and comments and Natalia Sotnikova for providing the SPH code software. This work was supported in part by the Russian Foundation for Basic Research (project No. 14-02-00319 and 14-02-00464) and the Programme of Fundamental Research of the Russian Academy of Sciences "Fundamental Problems of the Solar System Study and Exploration."