EDP Sciences
Highlight
Free Access
Issue
A&A
Volume 604, August 2017
Article Number A38
Number of page(s) 16
Section Interstellar and circumstellar matter
DOI https://doi.org/10.1051/0004-6361/201630084
Published online 31 July 2017

© ESO, 2017

1. Introduction

The development of circumstellar disks around single stars and their formation during the collapse of molecular clouds were studied extensively in the past three decades (McKee & Ostriker 2007; Andrews et al. 2010). However, at the same time it was found that a substantial fraction of protoplanetary disks form around binaries (Duquennoy & Mayor 1991; Kraus & Hillenbrand 2009).

The assumed mechanism by which a interstellar cloud evolves into a binary is either by multiple fragmentation or stellar capture (Bate 2000; Bate et al. 2002; Wolf et al. 2001). During the collapse both components can develop accretion disks in the same way as a single star system (Regály et al. 2011; Müller & Kley 2012). In addition to that, the dust and gas from outside of the binary orbit form a third disk, the so-called cirumbinary disk (Rodriguez et al. 2010; Romero et al. 2012; Lines et al. 2015), which revolves around the centre of mass of the binary. Furthermore, it is possible for a newly developed single star system to capture a transiting star to form a binary. In this case the primary could retain a part of its original disk and the captured star would accrete a disk of its own. The remaining matter would form the circumbinary disk.

Regardless of the exact formation mechanism the major difference from single star systems is that the orbiting binary generates a time-dependent non-axisymmetric gravitational potential, inducing strong tidal forces on the surrounding gas and dust (Artymowicz & Lubow 1994). This leads to a change in basic physical properties of the circumbinary disk, which should in turn result in unique structures such as a cavity in the disk centre, accretion arms, density waves, and spiral arms (Günther & Kley 2002; Günther et al. 2004). Previously it was shown that some of these structures, such as the cavity and parts of the accretion arms, are observable with the current generation of instruments (Ruge et al. 2015). There are also examples of more complex systems, such as GG Tau, which is a triple star system (Piétu et al. 2011; Di Folco et al. 2014). Some of the disk structures, for example spiral arms potentially associated with embedded planets, have also been observed (Muto et al. 2012; Garufi et al. 2013; Stolker et al. 2016).

The presence of a planet further complicates the case (Pierens & Nelson 2013). A binary can also feature two types of planetary orbits: planets that revolve around a single binary component and planets that orbit around both stars (Welsh et al. 2012; Orosz et al. 2012; Schwamb et al. 2013). Furthermore, the growth of planetesimals from smaller objects might be hindered (Meschiari 2012; Marzari et al. 2013). The migration behaviour of planets in a binary system can also be problematic since the planet can be ejected if it enters a 4:1 resonance with the binary (Kley & Haghighipour 2014; Pierens & Nelson 2008).

As mentioned above, the binary exerts a tidal force on the surrounding disk, which should result in characteristic structures. However, specific questions regarding the exact nature of these structures are still open: how do the characteristic spatial scales and timescales depend on parameters such as the mass and binary separation? How does the overall density and temperature structure compare to single star systems? Can structures induced by the binary be observed with the current or the next generation of instruments/observatories?

To answer the above questions, we simulate the density distribution of a circumbinary disk, using the 2D hydrodynamic code Fosite (Illenseer & Duschl 2009). Subsequently, radiative transfer simulations are performed to obtain the corresponding temperature structure of these disks with the code Mol3D (Ober et al. 2015). Furthermore, scattered light and re-emission maps of these systems are generated. In each step the unique structures and quantities generated by the asymmetrical gravitational potential of the binary are discussed. Finally, a feasibility study to observe these quantities with the Atacama Large Millimeter/submillimetre Array (ALMA) and the future European Extremely Large Telescope (E-ELT) operating at optical/infrared wavelengths is performed.

The paper is divided into three parts. First, we provide a general overview of the underlying 2D hydrodynamic simulations and 3D radiative transfer calculations (Sect. 2). Subsequently, the structures caused by the binary-disk interaction are analysed (Sect. 3.1). In the third part (Sects. 3.2 and 3.3) we discuss the feasibility of observing distinctive structures resulting from the binary-disk interaction.

2. Hydrodynamic and radiative transfer simulations

In this chapter the applied hydrodynamic and radiative transfer simulation software is introduced. Furthermore, the general simulation set-up is presented and discussed.

2.1. Two-dimensional hydrodynamic simulation

We simulate the density distribution of the gas and dust around a binary system for a fixed parameter set. We consider the heating by the star by calculating the temperature for an initial density distribution after the system has reached a quasi-stationary stage. This density distribution is calculated using a time- and spatial-constant temperature profile. With the new time-independent temperature profile a second density distribution, which is used in our study, is calculated. With this iteration a more realistic hydrodynamic simulation can be performed compared to a simulation with a constant temperature distribution. Similar work was performed by Günther & Kley (2002) and Günther et al. (2004). The major difference between these studies and our approach is the calculation of the temperature distribution. In case of Günther & Kley (2002) and Günther et al. (2004), the heating via the viscosity and the stellar irradiation were considered in each time step of the hydrodynamic simulation. By avoiding the calculation of the temperature in each time step of the hydrodynamic simulation we can decrease the required simulation time and achieve a higher spatial resolution. The viscous heating of the disk can be neglected, since the stellar irradiation is the dominant form of heating in circumbinary disks as was shown by Günther et al. (2004).

The hydrodynamic simulations are performed with Fosite (Illenseer & Duschl 2009). The code solves the 2D continuity and Navier-Stokes equations for the vertically integrated surface density Σ, pressure Π, and gas velocity v. The quantity T is the viscous stress tensor and Φ the gravitational potential of the binary, since the self-gravitation of the disk is neglected. The hydrostatic scale height H of the disk depends on the midplane temperature Tc and the gravitational potential and can be written as follows (see Günther & Kley 2002): (3)where cs is the speed of sound and Mi the masses of the components.

Table 1

Hydrodynamic simulation parameters and ranges.

Binary system:

The masses of the primary and secondary star are chosen to be equal to reduce the parameter space and inner radius rin of the computational domain (see below). A small inner radius is chosen to trace the disk structures up to the binary orbit. Furthermore, with equally massive components the system obtains a point symmetry, making it easy to spot possible non-physical behaviour in the results. As was shown by Bate et al. (2002), the mass tends to be accreted by the lower mass component, as it sweeps a larger area of the disk. This mechanism could lead to binary systems with mass ratio q ≈ 1. Long period systems with mass ratios ~1 have also been observed (Raghavan et al. 2010). The mass ratio of the primary and secondary is q = Msec/Mprim and we denote from now on Mprim = Msec = MB.

Three different semi-major axes a are applied to study their influence on the structures generated by the binaries, such as the inner cavity, accretion arms, and radial inhomogeneities in the density distribution (a = 10 AU, 20 AU, and 30 AU). Furthermore, for a separation of 20 AU three different stellar masses MB are applied to examine the influence of the stellar mass on the disk structures. The full set of parameters is summarised in Table 1.

The eccentricity of the orbit is set to zero, which reduces the parameter space and simplifies the system. Simultaneously, it allows one to reduce the inner radius rin (see below), since the grid centre is located at the centre of mass of the binary. At this point it is important to note that although there are some examples of binary systems with eccentricity ε < 0.1, the majority of them have a short orbital period P < 1 yr. Systems with larger values of P tend to have smaller eccentricities (Raghavan et al. 2010).

Temperature:

All calculations are performed in a locally isothermal mode based on a temperature profile Tc. We derive the temperature profile by performing a low-resolution hydrodynamical simulation for each disk set-up with a constant temperature profile as a starting condition. Subsequently, the dust temperature distribution is calculated with the 3D radiative transfer code Mol3D (Ober et al. 2015, see Sect. 2.2). The resulting azimutally averaged temperature profile is then used as default temperature Tc for the high resolution hydrodynamic calculation.

Initial and boundary conditions:

The initial disk mass is concentrated in a Gaussian surface density profile with its maximum at 150 AU. The mass is fixed at Mdisk = 0.02 M for all hydrodynamic simulations. The gas is set in a near Keplerian motion that accounts for the potential asymmetry with a higher velocity to counter additional gravitational attraction.

In order to reduce the impact of boundary conditions on the density distribution at the relevant region (inner 300 AU), the outer radius is set to 1000 AU for all simulations. The density as well as the gas velocities are sufficiently small at this boundary so as not to influence the relevant areas of the disk. Furthermore, the choice of the inner radius has an impact on the disk structure (Pierens & Nelson 2013). It cannot be chosen to be zero because the sources of gravitational potential have to be outside of the simulated region. An inner radius of rin = a·0.6 AU is used, where a is the semi-major axis of the binary. This ensures that no singularities occur inside the grid. It also permits us to trace the accretion flow nearly up to the binary orbits. This is important, as the biggest differences to a circumstellar disk are expected to be in the inner disk region. The inner boundary condition allows for a mass flux from inside the computational domain through the boundary, but prevents any flow in the opposite direction. This has a distinct disadvantage of preventing the circumprimary or circumsecondary disks from forming. The lack of matter in the direct vicinity of both stars should result in short wavelength radiation deficiency in the radiative transfer simulations. The accretion rate at the inner boundary amounts to ~ 10-8M/ yr and hence does not significantly change the mass of the binary or the disk for the duration of the simulation.

Grid.

In order to achieve a sufficient resolution for the subsequent radiative transfer simulation, especially in the inner regions, a polar grid (r,φ) with logarithmic spacing for the r coordinate and a grid resolution of (Nr = 508,Nφ = 508) is chosen.

Simulation time.

The evolution of the system is calculated over a period of tsim = 5 × 104 yr for all simulations. For the subsequent temperature calculations, density distributions with time steps between 2.2 × 104 yr and 2.7 × 104 yr are used. This decision is based on two constraints. First, the gas has to be sufficiently mixed to ensure that the initial state does not influence the resulting distribution. The timescale related to the advective mixing processes can be estimated using the dynamical timescale, (4)Here, vφ is the φ component of the velocity and Ω is the angular velocity. By choosing a simulation time step corresponding to over 100 binary periods this condition is met. The second constraint is to avoid losing too much of the initial mass. The timescale for this process is the viscous timescale, which can be estimated using the radial drift velocity vr, (5)This timescale is only reached in the simulation in the innermost regions of the disk. The viscous timescale at 200 AU, which is the outer radius of the disk in radiative transfer simulations, is of the order of 106 yr. The averaged accretion rate is of the order of 5 × 10-10M/ yr, which is consistent with observations (Bary & Petersen 2014).

The result of the hydrodynamic simulations are the circumbinary disk surface density Σ, scale height H, and the positions of the binary components. Based on those, radiative transfer simulation are performed.

2.2. Radiative transport simulation

In the following we describe the radiative transfer simulation. The density distribution from the previous section is used to derive the corresponding temperature distribution and to calculate synthetic scattered light and thermal re-emission maps. These provide the basis for the assessment of whether the generated structures can be observed.

The radiative transfer simulations are performed with the 3D continuum and line radiative transfer code Mol3D (Ober et al. 2015). First, a temperature is calculated solving the radiative transfer equation using a Monte Carlo algorithm and taking the optical properties of the dust into account. The temperature is then used to simulate emission maps. In addition, scattered light maps are calculated with the Monte Carlo algorithm as well.

Grid.

A spherical grid is used with a logarithmic scaling of the r coordinate to ensure that small structures near the binaries, in particular the temperature gradient at the inner disk rim and accretion arms, are sufficiently resolved.

While the inner radius rin is identical to that of the hydrodynamic simulations, the outer radius rout is set to 200 AU because most of the gas and all significant density structures are located inside this radius. All simulations are performed for a face-on inclination of i = 0°.

Gas and dust distribution.

Mol3D employs a perfectly mixed gas and dust with a gas-to-dust mass ratio of 100:1. Based on the isothermal hydrodynamically calculated 2D surface density distribution, the 3D density distribution ϱ is constructed from (6)where z is the z coordinate (Lynden-Bell 1969).

The normalised density distribution calculated with Fosite is scaled accordingly, allowing us to consider three different disk masses Mdisk = 10-1M, 10-2M, 10-3M. This provides the basis for studying the effect of optical depth on the emission maps.

Dust.

The dust is considered to be homogeneous spheres of radius a and with the following size distribution (Dohnanyi 1969): (7)Here, a minimum radius of amin = 5 nm and maximum radius amax = 250 nm (Mathis et al. 1977) are assumed. A composition of 62.5% astronomical silicate and 37.5% graphite is used with optical data from Weingartner & Draine (2001). Applying Mie theory, the optical properties are calculated at 100 logarithmically distributed wavelengths λsim in the interval between 0.05 and 2000μm (Wolf & Voshchinnikov 2004).

Binary luminosity.

To conduct the radiative transfer simulations of the dust density distribution constructed beforehand, the luminosities of the binary components are required. The stars are considered to be main sequence stars with a same chemical composition as the sun. The stellar luminosities corresponding to the masses are derived from the main sequence tracks from Siess et al. (2000).

The full parameter set for the radiative transfer simulation is shown in Table 2.

Table 2

Radiative transfer simulation parameters and ranges.

3. Results

In this chapter, we present the results of the simulations described in Sect. 2. We discuss the most notable features of binary systems and their influence on the radiation emitted by the surrounding disk, i.e. observable quantities.

3.1. Characteristic structures in the density distribution of circumbinary disks

First we present the results of the hydrodynamic simulations and discuss the structures caused by the disk-binary interaction and the mechanism creating them.

Figure 1 shows an image of the surface density distribution for a semi-major axis a = 30 AU. The most notable feature of the distribution is the cavity in the inner region of the disk. The cavity appears if the binary can exert enough torque on the surrounding disk to counter the torque generated by the viscous stresses. As can be gathered from Fig. 2, the radius of the cavity is approximately 2−2.5 times the binary separation. This is consistent with the theoretical predictions (Artymowicz & Lubow 1994). As mentioned in Sect. 2.1 the choice of equally massive stars and an eccentricity of ε = 0 leads to an almost circular cavity centred on the barycentre of the binary. Since the accretion arms have to stretch from the inner edge of the disk to the binary orbit, larger cavities result in longer accretion arms. The accretion arms also show a distinct point symmetry that remains undisturbed for at least 100 binary orbits. This finding is in agreement with the results of Günther & Kley (2002).

thumbnail Fig. 1

Normalised surface density for a = 30 AU, MB = 1 M. The crosses indicate the positions of the binary components.

Open with DEXTER

The hydrostatic scale height H only depends on the temperature and gravitational potential in our simulations. These are either provided as a fixed initial condition or only changes slightly for radii ra during the simulation. Consequently, the hydrostatic scale height profile is smooth and monotonically is increasing. This means that in the current model only the shadowing effects that result from the increase of the surface density, in contrast to variation of H, can be studied.

Another striking feature of the investigated binary systems are the density waves that extend radially from the inner edge of the disk. In Fig. 2 the surface density Σ, azimuthally averaged and normalised by the total mass is shown for three different values of the binary semi-major axis a. We identify two important quantities: the wavelength of the surface density oscillation λa, measured as the distance between two maxima, and its magnitude a, measured as the difference of the maximum and its subsequent minimum. One can clearly see that both quantities depend on the semi-major axis a. To quantify this dependence, the wavelength λa was calculated by first averaging over the distance of the maxima in each time step of the simulation and averaging the result over 100 time steps in a time interval of 104 years. In Table 3 the resulting wavelength λa and magnitude a are shown. One can see a trend of increasing wavelength λa with increasing a. The same applies to the magnitude a of the density oscillation, i.e. an increase of the semi-major axis leads to an increase of the magnitude of the density oscillation.

thumbnail Fig. 2

Normalised by the total mass and azimuthally averaged surface density for semi-major axes a = 10 AU, 20 AU, and 30 AU; MB = 1 M.

Open with DEXTER

thumbnail Fig. 3

Azimuthally averaged normalised surface density for a = 20 AU, and MB = 0.5 M,1 M, and 1.5 M.

Open with DEXTER

Table 3

Wavelength and magnitudes for different semi-major axes a and binary masses MB.

In Fig. 3 the impact of binaries with different masses on the circumbinary disk is illustrated. We find that the radius of the cavity does not depend on the binary mass. By applying the same procedure as for the calculation of wavelength λa and the magnitude a, the values of λMB and MB can be determined. Here, λMB is the oscillation wavelength of the surface density and MB magnitude with regard to binary mass MB. We find that with increasing binary mass MB the wavelength of the oscillation λMB decreases. At the same time the trend for the mass-dependent magnitude MB appears to be more complex, as the values for MB = 0.5 M and MB = 1.5 M are smaller than for MB = 1 M. This suggests that there is a value of MB for which the magnitude reaches a maximum. Similar density waves were shown in Fig. 2 of the work from Günther et al. (2004).

At this point it appears necessary to discuss the theoretical basis of these density waves. Starting with the 2D continuity and Navier-Stokes equations we calculate (8)which is a 1D wave equation in polar coordinates (full derivation can be found in Appendix A). The interpretation of this result is the following: while the gas is accreted inward with velocity v, it provides a moving frame for density waves. These density waves are driven by the binary in the centre and propagate outward with sound speed cs. This interpretation also explains the trend of wavelength (λa, λMB) with regard to the semi-major axis a and binary mass MB. The system behaves as a driven harmonic oscillator that swings with the frequency of the instigator f. This frequency is coupled with the binary period P. Neglecting the dispersion, one would get the formula for the wavelength λ = cs/fcs·P. By increasing the semi-major axis a or decreasing the binary mass MB, we increase the orbital period and with that the wavelength (λa, λMB). The magnitude of the oscillation (a, MB) should be proportional to the amount of torque that the binary can exert on the disk. This seems to be the case with increasing values of the semi-major axis a. However, one would assume that a binary with higher mass MB would be able to generate a stronger torque. Instead, we find that the magnitude MB first increases with increasing binary mass, but then decreases for value of MB = 1.5 M. The exact nature of this phenomenon should be investigated in a future study.

3.2. Observability of characteristic structures: analysis of ideal observations

The goal of this study is to investigate whether the characteristic structures in a circumbinary disk can be observed with currently operating and future instruments/observatories. For this purpose we simulate various observable quantities, such as scattered light and re-emission images, as well as the spectral energy distribution (SED) in the wavelength range between 0.05 μm and 2000 μm of the circumbinary disks discussed in Sect. 3.1. A distance of 140 parsec to the object is assumed for all simulations.

In the hydrodynamic simulations it was shown that the most striking differences from an undisturbed circumstellar disk are the large gas and dust depleted cavities in the centre, accretion arms in the cavity, and density waves at the inner rim of the disk. Of course, a central cavity is also characteristic for circumstellar transitional disk. However, in the case of a transitional disk this cavity results from disk evolution. With the calculated surface density and scale height it is now possible to derive a 3D density distribution ϱ(r,θ,φ) according to Eq. (6). This density is post-processed with Mol3D to derive observable quantities (see Sect. 2.2 for details).

Figure 4 shows an exemplary density distribution for a = 20 AU and MB = 1 M. Although the hydrostatic scale height H is fairly smooth as a function of the radial distance to the central star, a variation of the disk profile in the plane perpendicular to the midplane can be observed. From Eq. (6) it is obvious that the variation of the surface density causes a shadow that screens parts of the inner disk edge (Fig. 5). Accretion arms caused by the presence of a binary, although not very dense, have enough mass to reach optical depths τ ≳ 1 (see Fig. 5).

thumbnail Fig. 4

Normalised dust density distribution for a = 20 AU, MB = 1 M. The crosses indicate the positions of the binary components.

Open with DEXTER

thumbnail Fig. 5

Left: optical depth τ at λ = 526 nm calculated in the midplane. Right: resulting midplane temperature of the same system. Black line denotes where the optical depth reaches values of τ = 1 and τ = 104, respectively.

Open with DEXTER

In the following we discuss the influence of the binary on the observational properties of the disk.

3.2.1. Spectral energy distribution

The cavity in the centre of the disk is clearly visible in all calculated density distributions and emission maps and its observability is the topic of Sect. 3.3. Because of the cavity, circumbinary disks tend to have less dust, which can be heated to a high temperature. This results in much lower flux at near- and mid-infrared spectrum in comparison to circumstellar disks (Fig. 6). Because the cavity radius scales with the semi-major axis a, the flux at short wavelengths is in indirect proportion to a (Fig. 7). The flux of the corresponding circumstellar disk (i.e. same disk mass and net stellar mass) is higher for all wavelengths (Fig. 6). It was simulated using the standard set-up presented in Sect. 2.1 and a star with the mass of Mstar = 2 M. The luminosity of this star is higher than the combined luminosity of two stars with individual masses Mstar = 1 M. This leads to an overall higher temperature and thus to higher fluxes. This opens a potential way for distinguishing between disks around binaries and single stars: the outer edges of the circumbinary disk behave nearly Keplerian. Through the line broadening one can determine the mass of the central star. For the same mass of the central star (or binary), the decreased flux due to the cavity and the lower disk temperature (due to the lower net stellar luminosity) result in a lower infrared emission in the case of the circumbinary disk.

thumbnail Fig. 6

Spectral energy distribution for two simulated disks with Mdisk = 10-1M, stellar mass of 2 M. Blue: single star; red: binary with two components of MB = 1 M. Inner cavity radius for the binary ≈ 40 AU.

Open with DEXTER

thumbnail Fig. 7

Dependence of the SED on the semi-major axes a = 10,20,30 AU (MB = 1 M).

Open with DEXTER

thumbnail Fig. 8

Dependence of the SED on the binary disk mass. a = 20 AU, MB = 1 M, and Mdisk = 10-1−10-3M.

Open with DEXTER

At longer far-infrared to millimetre wavelengths the differences of the fluxes arising from circumbinary disks with different cavity radii (i.e. different semi-major axes of the central binary) become negligible (Fig. 7). This is expected since, for large distances from the central binary, the disks structures are similar.

The dependence on the disk mass can be seen in Fig. 8. Both the short wavelength and long wavelength radiation increases with disk mass, although the increase is caused by different mechanisms. The presence of more hot dust in the vicinity of the binary leads to higher levels of short wavelength radiation. For long wavelength radiation, for which the disk is optically thin in the vertical direction, higher dust mass leads to a higher column density of radiating particles.

3.2.2. Surface brightness distribution

So far, we only considered the synthetic SEDs of the simulated disks. However, as both the contributions of the scattered and re-emitted radiation depend solely on the radial density distribution, an interpretation of the SED to reveal the characteristic disk structures discussed in Sect. 3.1 is rather limited. We now investigate the impact of these characteristic structures on spatially resolved observations in various wavelength regimes.

Disk mass.

The effect the disk mass has on observable features is best seen in the (sub)millimetre regime, since the observed column density is much higher. Figure 9 shows the two surface brightness distribution for two disk masses at λ = 324μm. The cavity in the centre of the disk is not affected by the mass increase. As the disks are optically thin in this wavelength range, the re-emission flux scales approximately linearly with the disk mass.

thumbnail Fig. 9

Surface brightness distribution at λ = 324μm for disk masses Mdisk = 10-1M (left) and Mdisk = 10-3M (right).

Open with DEXTER

Flux maximum.

Figure 10 shows the spatial flux density distribution at different wavelengths. The near- and mid-infrared radiation reaches its maximum along the accretion arms, which are the innermost structures, directly exposed to the stellar radiation. The thermal re-emission from the inner disk edge, where the density waves are strongest reaches its maximum at submillimetre wavelengths.

thumbnail Fig. 10

Surface brightness distribution at λ = 10.5μm (left) and 324μm (right) (a = 20 AU, MB = 1 M).

Open with DEXTER

Density waves.

In Sect. 3.1 we presented a quantitative description of the dependence of the density waves generated by the binary on the semi-major axis a and the binary mass MB. We now investigate the wave structure in the resulting surface brightness distributions employing a similar algorithm as in Sect. 3.1. At first, we average the brightness distribution in azimuthal direction. Subsequently, we determine the locations of maxima and minima and calculate the flux difference between the maxima and their surrounding minima. From both differences, the smaller one is chosen, which ensures that in case of sufficient sensitivity at least one flux maximum is detected. After performing this procedure for all maxima, the largest difference for each parameter configuration, defined by the binary mass MB, disk mass Mdisk, semi-major axis a, and wavelength λ, is chosen. This allows us to derive the lower boundary for the sensitivity required to detect at least one flux maximum.

Figure 11 shows the flux differences as a function of the wavelength of the emitted radiation. Each figure represents a different semi-major axis a and the colours denote different disk masses for each binary configuration. Figure 12 depicts the same quantities for different binary masses MB.

We start at the long wavelength end of these graphs to find an explanation for the trend observed here. At wavelengths λ> 1 mm the disks are optically thin (τ ≈ 1 in the most dense regions for Mdisk = 10-1M). The linear increase of the flux towards shorter wavelengths is due to the higher fluxes at these wavelengths (see e.g. Fig. 8). Additionally, the amplitude of the density waves is the highest at the inner disk rim (Figs. 2 and 3). As the flux maximum moves inward for shorter wavelengths, the flux at the location of the largest magnitudes increases, resulting in the further increase of the flux differences. At about 100 μm the flux difference decreases abruptly. This is the wavelength at which the disk becomes optically thick. At even shorter wavelengths, one no longer traces the full column density, but only the wavelengths-dependent photosphere of the disk. This interpretation also explains why this decrease happens at shorter wavelengths for disks with less mass.

The remaining feature to be discussed is the bump at wavelengths between the increase of the optical depth and about 10 μm. In this wavelength range the flux difference shows a more complex dependence on the parameters than in the cases discussed previously. It increases with disk mass, but not linearly (in a log-log diagram) as before. Here, the temperature at the inner disk rim influences the maximum differences. These are larger for smaller semi-major axis a values and higher binary masses MB. In both cases the temperature increases owing to the proximity of the binary orbit to the disk inner rim or owing to higher stellar luminosity. At these short wavelengths the regions of the disk contributing most to the flux no longer have azimuthal symmetry that is necessary for the applied algorithm. For this reason these parts of the graphs are not discussed further.

thumbnail Fig. 11

Flux difference for MB = 1 M and a = 10, 20, 30 AU (see Sect. 3.1 for details).

Open with DEXTER

thumbnail Fig. 12

Flux difference for MB = 0.5, 1.5 M and a = 20 AU (see Sect. 3.1 for details).

Open with DEXTER

thumbnail Fig. 13

Surface brightness distribution at λ = 4.5μm (left) and λ = 20μm (right) for MB = 1 M and a = 10 AU.

Open with DEXTER

thumbnail Fig. 14

Surface brightness distribution at λ = 4.5μm (left) and λ = 20μm (right) for MB = 1 M and a = 20 AU.

Open with DEXTER

thumbnail Fig. 15

Surface brightness distribution at λ = 4.5μm (left) and λ = 20μm (right) for MB = 1 M and a = 30 AU.

Open with DEXTER

thumbnail Fig. 16

Surface brightness distribution at λ = 4.5μm (left) and λ = 20μm (right) for MB = 0.5 M and a = 20 AU.

Open with DEXTER

thumbnail Fig. 17

Surface brightness distribution at λ = 4.5μm (left) and λ = 20μm (right) for MB = 1.5 M and a = 20 AU.

Open with DEXTER

Accretion arms.

As discussed earlier, the appearance of the disk at infrared wavelengths is dominated by the thermal re-emission and scattered light of the accretion arms (see Fig. 10). In this part of the disk, the gravitational potential differs significantly from a standard Keplerian potential. This leads to a high temporal variability and to the lack of the azimuthal symmetry. For those reasons, we abstain from attempting to quantify our results as we did with density waves. Instead, we only outline general trends with regard to the various parameter spaces.

Figures 1317 show the surface brightness distribution for semi-major axes a = 10, 20, 30 AU, MB = 1 M, and masses MB = 0.5, 1.5 M, a = 20 AU at two different wavelength, λ = 4.5μm and λ = 20μm.

Here we can note a major trend that was partially mentioned before, namely that greater values of the semi-major axis a lead to greater distances of disk edge to the binary orbit and therefore lower fluxes. Similarly, a higher binary mass leads to a higher stellar luminosity and consequently to a higher thermal re-emission flux. A further trend is the increase of flux with wavelength (see SEDs in Figs. 68). In Table 4 the fluxes originating from inside a circle with radius of 2 × a are compiled. From these we can deduce that the binary mass MB and thus its luminosity has the biggest impact on the fluxes.

Table 4

Integrated flux of the scattered and re-emitted radiation in units of mJy originating from inside a circle with the radius of 2 × a.

Taking the trends discussed above into account, one can derive predictions for the observability of circumbinary disk features (at this stage we disregard the influence of a realistic observing instruments; this is the topic of the following sections).

  • The accretion arms can easiest be observed in the infrared, wherewe find the highest fluxes for them. Systems with smallersemi-major axis values a and higher binary masses MB are preferable for such observations.

  • The flux of a disk with a higher dust mass is significantly larger, which leads to the conclusion that more massive disks are more suitable for observations in infrared and in millimetre radiation.

  • For observations at submillimetre/millimetre wavelengths we find that higher values of the semi-major axis a result in larger flux differences for the density waves. The binary mass has no significant impact on the flux difference, but one should expect the overall flux to be higher for systems with higher binary masses.

3.3. Simulated observations: specific instrument studies

We now study synthetic observations based on the analysis of ideal surface brightness distributions. First, we consider simulated observations at submillimetre/millimetre wavelengths. Subsequently, simulated observations with the future E-ELT are presented, which will operate at optical to infrared wavelengths.

3.3.1. ALMA

In this subsection we evaluate the observability of the characteristic structures discussed in Sect. 3.2 with the (sub)-millimetre interferometer ALMA. In particular these are the inner cavity and density waves. At first, we specify basic technical data assumed in this observational feasibility study.

ALMA set-up.

We perform simulations within the frame of capabilities in Observational Cycle 4, which employs 40 of the 12 m antennas in nine different configurations with maximum baselines between 155 m and 12 644 m (seven wavelength bands between 0.32 mm and 3.6 mm). Despite the fixed number of antennas and configurations, the results presented here are of general nature because they only depend on the simulated scattered light and re-emission maps and the chosen angular resolution and sensitivity. The synthetic observations are performed with the ALMA simulation tool kit CASA (McMullin et al. 2007). The on-source time for all simulated observations is set to 6 h with the thermal noise option enabled and recommended water vapours values applied. A declination representative for the nearby star-forming region in Taurus is used (δ ≈ 22°).

Density waves.

The goal of this study is to determine whether the density waves can be observed with ALMA and, if so, which configuration and band are the most suitable to perform these observations with regard to the considered binary system. We begin with an example of a synthetic ALMA observation (see Fig. 18). For this, we choose the shortest possible wavelength ALMA is capable of observing (λ = 320μm, band 10) and the configuration with the maximum baseline in that band (D = 3697 m), which results in an angular resolution of θres = 0.024″. A quick look reveals that, at least for the semi-major axis a = 30 AU and disk mass of Mdisk = 0.1 M, the sensitivity and angular resolution are sufficient for the detection of the density waves (see Fig. 18). As a next step we quantify the flux difference between the wave maxima and minima in units of sensitivity σ, employing an algorithm similar to that applied in Sect. 3.2.2. The results are shown in Figs. 1923.

We begin the discussion of the results with the highest values of the flux difference, which are measured in units of sensitivity σ. For all configurations the maximum is reached at λ = 849μm. Similarly, the higher values of the flux difference are concentrated at longer wavelengths and configurations with shorter maximum baselines. However, the resolution at those wavelengths and configurations is not sufficient to resolve the density waves. What we detect is the flux increase with the radius, which corresponds to the inner disk cavity. The combination of wavelength and configuration for which we can detect at least one density maximum are labeled with the appropriate value of flux difference. The detectability of the density waves increases with the semi-major axis a. This is not surprising taking into account the results of Table 3. The wavelength and the magnitude of the density wave increase with the semi-major axis a. However, it is surprising, that the detectability is higher for systems with lower binary mass. From Table 3 we can infer that the wavelength is increasing for decreasing binary mass, which results in more pronounced maxima. However, the magnitude for a binary mass MB = 1 M is almost a factor of two greater than for MB = 0.5 M. From this we conclude that the wavelength of the density wave has a greater influence on the detectability than the wave magnitude.

Finally, we want to compare our results to a similar study performed by Ruge et al. (2015). One of the major differences of our work compared to the results of Ruge et al. (2015) is the lack of density waves. There are three major reasons for this. The first reason is the difference in the simulations. Ruge et al. (2015) conducted a SPH simulation of the inner disk regions. The hydrodynamic simulation performed in this study achieves a better resolution, which allows us to track smaller structures. The second reason is the parameters of the systems. A semi-major axis of a = 2 AU, considered by Ruge et al. (2015), should result in much smaller density wave wavelength than studied here. The last reason are the considered observing wavelengths. The shortest wavelength employed by Ruge et al. (2015) was 750 μm. As can be seen in Fig. 11, the flux difference between the maximum and minimum decreases with wavelength. Thus, a longer wavelength should decrease the feasibility of detecting the density waves. The variability in flux distribution on the inner disk rim and spiral arm seen in Fig. 7 in the study of Ruge et al. (2015) is caused by the eccentric orbit of the binary (ε = 0.3).

3.3.2. European Extremely Large Telescope

In Sect. 3.2.2 it became clear that the accretion arms are best observed at infrared wavelengths. To simulate an observation with the E-ELT instrument METIS, we convolve the ideal maps (Figs. 1317) with the wavelength-dependent point spread function of an ideal 39 m telescope with a circular aperture. The goal here is to determine whether it is possible to observe the highly variable structures, such as the inner disk rim and accretion arms.

The brightness ratio between the individual binary component and the brightest structures amounts to ~ 102− ~ 104 in the considered cases. For this reason, the binaries are removed from the maps before the convolution (see Figs. 2428).

All images at λ = 4.5μm are dominated by the flux from the binary. However, in all systems, with the exception of the case MB = 1 M and a = 10 AU, it is still possible to spatially resolve the accretion arms. Similarly, at λ = 20μm one would only detect a bright ring.

thumbnail Fig. 18

Synthetic ALMA observation at λ = 320μm with configuration # 7; MB = 1 M and a = 30 AU.

Open with DEXTER

thumbnail Fig. 19

Flux difference between density wave maximum and minimum in units of sensitivity for a = 10 AU, MB = 1 M, and Mdisk = 10-1M.

Open with DEXTER

Table 5

Integrated flux of the scattered and re-emitted radiation in units of mJy originating from inside a circle with the radius of 2 × a.

thumbnail Fig. 20

Flux difference between density wave maximum and minimum in units of sensitivity for a = 20 AU, MB = 1 M, and Mdisk = 10-1M.

Open with DEXTER

thumbnail Fig. 21

Flux difference between density wave maximum and minimum in units of sensitivity for a = 30 AU, MB = 1 M, and Mdisk = 10-1M.

Open with DEXTER

thumbnail Fig. 22

Flux difference between density wave maximum and minimum in units of sensitivity for a = 20 AU, MB = 0.5 M, and Mdisk = 10-1M.

Open with DEXTER

thumbnail Fig. 23

Flux difference between density wave maximum and minimum in units of sensitivity for a = 20 AU, MB = 1.5 M, and Mdisk = 10-1M.

Open with DEXTER

In all but one system the accretion arms are resolved in both wavelength regimes. In the case of MB = 1 M and a = 10 AU, we only detect a bright ring that corresponds to the inner disk rim. For λ = 4.5μm we find two small flux maxima that are caused by hot dust in the direct vicinity of the binary at the edge of the computational domain. In Table 5 the fluxes originating from inside a circle with radius of 2 × a are compiled. We find no major differences from the original maps caused by the convolution process (see Table 4).

4. Summary

The goal of this work was to study the observability of characteristic structures in circumbinary disks, generated by binary-disk interaction. To achieve this goal, a 2D hydrodynamic simulation was performed with the code Fosite to calculate a surface density distribution. The results were applied to a 3D grid-based Monte Carlo code Mol3D and a temperature profile; in addition, resulting scattered light and re-emission maps were simulated. These images were used to simulate imaging observations in the near to mid-infrared wavelength range (E-ELT) and at submillimetre/millimetre wavelengths (ALMA).

We found that the torque exerted on the circumbinary disk by the binary generates unique features that can be used to distinguish them from protoplanetary disks around single stars. In particular these are the inner cavity, accretion flows from disk inner rim to the binary orbit, and the density waves on the disk inner edge. We quantified the dependence of those features on binary mass MB, semi-major axis a, and disk mass Mdisk. Furthermore, we derived a wave equation governing the behaviour of the density waves, which can be transformed into a standard 1D wave solution in two dimensions.

We have shown that those features alter the thermal and radiation balance of the disk to a sufficiently high degree, which can be seen in the re-emission and scattered light maps. Here we once again studied the dependence of the observable features on the parameters and derived predictions with regard to the successful observation of those features.

thumbnail Fig. 24

Flux after convolution at λ = 4.5μm (left) and λ = 20μm (right) for MB = 1 M and a = 10 AU.

Open with DEXTER

thumbnail Fig. 25

Flux after convolution at λ = 4.5μm (left) and λ = 20μm (right) for MB = 1 M and a = 20 AU.

Open with DEXTER

thumbnail Fig. 26

Flux after convolution at λ = 4.5μm (left) and λ = 20μm (right) for MB = 1 M and a = 30 AU.

Open with DEXTER

thumbnail Fig. 27

Flux after convolution at λ = 4.5μm (left) and λ = 20μm (right) for MB = 0.5 M and a = 20 AU.

Open with DEXTER

thumbnail Fig. 28

Flux after convolution at λ = 4.5μm (left) and λ = 20μm (right) for MB = 1.5 M and a = 20 AU.

Open with DEXTER

Subsequently synthetic observations of the simulated circumbinary disks with ALMA and E-ELT were generated. We could show that ALMA configurations with the highest angular resolution are capable of observing the density waves generated by the binary with a = 20 AU and a = 30 AU. Furthermore, we made a parameter study that showed that even a configuration with a far smaller maximum baseline is sufficient to observe the inner cavity of those disks. Simulated observations with E-ELT have shown that it will be possible to resolve the accretion arms at the considered wavelengths of 4.5μm and 20μm.

Further studies of this kind would need to include systems with binary mass ratio q = Msec/Mprim ≠ 1 and an eccentricity value greater than 0, as they comprise the majority of the known systems (Raghavan et al. 2010) and those parameters have a large impact on the disk structure (Ruge et al. 2015). Furthermore, this study was focused only on dust emissions. However, the gas of the circumbinary disk is influenced in the same manner. If realised, the observation of binary-induced structures in molecular lines would allow one to observe velocity fields in the disk and in the accretion flows in the disk centre.

Acknowledgments

We thank all the members of the Astrophysics Department Kiel for helpful discussions and remarks and for their language corrections. This study was funded by the German Science Foundation (DFG), grant: WO 857/12-1.

References

Appendix A: Derivation of the wave equation

In this section we present a detailed derivation of the wave equation mentioned in Sect. 3.1.

Since the wave propagates in r direction we neglect the φ component and terms containing the derivative with regard to φ, which leaves us with the 2D continuity and the r componentof the Navier-Stokes equations in polar coordinates, The parameter Σ is the vertically integrated surface density; vr is the radial drift velocity, which is in this case generated by the disks differential rotation combined with the sheer viscosity; and cs is the sound speed of the gas. By applying the time derivative on Eq. (A.1) we get (A.3)The time derivation of the surface density can be substituted using Eq. (A.1) once again and Eq. (A.2) to substitute Σtvr, (A.4)which can be further simplified to (A.5)Assuming the drift velocity vr and its derivation rvr are much smaller than the sound speed cs and its derivative rcs, Eq. (A.5) can be reduced to Eq. (8) A full summary of the solution for the simplified Eq. (8) can be found in Polyanin (2002), pp. 296, 297.

All Tables

Table 1

Hydrodynamic simulation parameters and ranges.

Table 2

Radiative transfer simulation parameters and ranges.

Table 3

Wavelength and magnitudes for different semi-major axes a and binary masses MB.

Table 4

Integrated flux of the scattered and re-emitted radiation in units of mJy originating from inside a circle with the radius of 2 × a.

Table 5

Integrated flux of the scattered and re-emitted radiation in units of mJy originating from inside a circle with the radius of 2 × a.

All Figures

thumbnail Fig. 1

Normalised surface density for a = 30 AU, MB = 1 M. The crosses indicate the positions of the binary components.

Open with DEXTER
In the text
thumbnail Fig. 2

Normalised by the total mass and azimuthally averaged surface density for semi-major axes a = 10 AU, 20 AU, and 30 AU; MB = 1 M.

Open with DEXTER
In the text
thumbnail Fig. 3

Azimuthally averaged normalised surface density for a = 20 AU, and MB = 0.5 M,1 M, and 1.5 M.

Open with DEXTER
In the text
thumbnail Fig. 4

Normalised dust density distribution for a = 20 AU, MB = 1 M. The crosses indicate the positions of the binary components.

Open with DEXTER
In the text
thumbnail Fig. 5

Left: optical depth τ at λ = 526 nm calculated in the midplane. Right: resulting midplane temperature of the same system. Black line denotes where the optical depth reaches values of τ = 1 and τ = 104, respectively.

Open with DEXTER
In the text
thumbnail Fig. 6

Spectral energy distribution for two simulated disks with Mdisk = 10-1M, stellar mass of 2 M. Blue: single star; red: binary with two components of MB = 1 M. Inner cavity radius for the binary ≈ 40 AU.

Open with DEXTER
In the text
thumbnail Fig. 7

Dependence of the SED on the semi-major axes a = 10,20,30 AU (MB = 1 M).

Open with DEXTER
In the text
thumbnail Fig. 8

Dependence of the SED on the binary disk mass. a = 20 AU, MB = 1 M, and Mdisk = 10-1−10-3M.

Open with DEXTER
In the text
thumbnail Fig. 9

Surface brightness distribution at λ = 324μm for disk masses Mdisk = 10-1M (left) and Mdisk = 10-3M (right).

Open with DEXTER
In the text
thumbnail Fig. 10

Surface brightness distribution at λ = 10.5μm (left) and 324μm (right) (a = 20 AU, MB = 1 M).

Open with DEXTER
In the text
thumbnail Fig. 11

Flux difference for MB = 1 M and a = 10, 20, 30 AU (see Sect. 3.1 for details).

Open with DEXTER
In the text
thumbnail Fig. 12

Flux difference for MB = 0.5, 1.5 M and a = 20 AU (see Sect. 3.1 for details).

Open with DEXTER
In the text
thumbnail Fig. 13

Surface brightness distribution at λ = 4.5μm (left) and λ = 20μm (right) for MB = 1 M and a = 10 AU.

Open with DEXTER
In the text
thumbnail Fig. 14

Surface brightness distribution at λ = 4.5μm (left) and λ = 20μm (right) for MB = 1 M and a = 20 AU.

Open with DEXTER
In the text
thumbnail Fig. 15

Surface brightness distribution at λ = 4.5μm (left) and λ = 20μm (right) for MB = 1 M and a = 30 AU.

Open with DEXTER
In the text
thumbnail Fig. 16

Surface brightness distribution at λ = 4.5μm (left) and λ = 20μm (right) for MB = 0.5 M and a = 20 AU.

Open with DEXTER
In the text
thumbnail Fig. 17

Surface brightness distribution at λ = 4.5μm (left) and λ = 20μm (right) for MB = 1.5 M and a = 20 AU.

Open with DEXTER
In the text
thumbnail Fig. 18

Synthetic ALMA observation at λ = 320μm with configuration # 7; MB = 1 M and a = 30 AU.

Open with DEXTER
In the text
thumbnail Fig. 19

Flux difference between density wave maximum and minimum in units of sensitivity for a = 10 AU, MB = 1 M, and Mdisk = 10-1M.

Open with DEXTER
In the text
thumbnail Fig. 20

Flux difference between density wave maximum and minimum in units of sensitivity for a = 20 AU, MB = 1 M, and Mdisk = 10-1M.

Open with DEXTER
In the text
thumbnail Fig. 21

Flux difference between density wave maximum and minimum in units of sensitivity for a = 30 AU, MB = 1 M, and Mdisk = 10-1M.

Open with DEXTER
In the text
thumbnail Fig. 22

Flux difference between density wave maximum and minimum in units of sensitivity for a = 20 AU, MB = 0.5 M, and Mdisk = 10-1M.

Open with DEXTER
In the text
thumbnail Fig. 23

Flux difference between density wave maximum and minimum in units of sensitivity for a = 20 AU, MB = 1.5 M, and Mdisk = 10-1M.

Open with DEXTER
In the text
thumbnail Fig. 24

Flux after convolution at λ = 4.5μm (left) and λ = 20μm (right) for MB = 1 M and a = 10 AU.

Open with DEXTER
In the text
thumbnail Fig. 25

Flux after convolution at λ = 4.5μm (left) and λ = 20μm (right) for MB = 1 M and a = 20 AU.

Open with DEXTER
In the text
thumbnail Fig. 26

Flux after convolution at λ = 4.5μm (left) and λ = 20μm (right) for MB = 1 M and a = 30 AU.

Open with DEXTER
In the text
thumbnail Fig. 27

Flux after convolution at λ = 4.5μm (left) and λ = 20μm (right) for MB = 0.5 M and a = 20 AU.

Open with DEXTER
In the text
thumbnail Fig. 28

Flux after convolution at λ = 4.5μm (left) and λ = 20μm (right) for MB = 1.5 M and a = 20 AU.

Open with DEXTER
In the text

Current usage metrics show cumulative count of Article Views (full-text article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.

Data correspond to usage on the plateform after 2015. The current usage metrics is available 48-96 hours after online publication and is updated daily on week days.

Initial download of the metrics may take a while.