# On the energy dissipation rate at the inner edge of circumbinary discs

Caroline Terquem, John C. B. Papaloizou

<figure><img src="image_1.png" /></figure>

## To cite this version:

Caroline Terquem, John C. B. Papaloizou. On the energy dissipation rate at the inner edge of circumbinary discs. Monthly Notices of the Royal Astronomical Society, 2017, 464 (2), pp.2429-2440. 10.1093/mnras/stw2485. hal-01554675

HAL Id: hal-01554675

https://hal.science/hal-01554675v1

Submitted on 27 Aug 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_3.png" /></figure>

# On the energy dissipation rate at the inner edge of circumbinary discs

Caroline Terquem<sup>1,2*</sup> and John C. B. Papaloizou<sup>3*</sup>

Accepted 2016 September 28. Received 2016 September 27; in original form 2016 August 5

## ABSTRACT

We study, by means of numerical simulations and analysis, the details of the accretion process from a disc on to a binary system. We show that energy is dissipated at the edge of a circumbinary disc and this is associated with the tidal torque that maintains the cavity: angular momentum is transferred from the binary to the disc through the action of compressional shocks and viscous friction. These shocks can be viewed as being produced by fluid elements that drift into the cavity and, before being accreted, are accelerated on to trajectories that send them back to impact the disc. The rate of energy dissipation is approximately equal to the product of potential energy per unit mass at the disc's inner edge and the accretion rate, estimated from the disc parameters just beyond the cavity edge, that would occur without the binary. For very thin discs, the actual accretion rate on to the binary may be significantly less. We calculate the energy emitted by a circumbinary disc taking into account energy dissipation at the inner edge and also irradiation arising there from reprocessing of light from the stars. We find that, for tight PMS binaries, the SED is dominated by emission from the inner edge at wavelengths between 1–4 and 10 µm. This may apply to systems like CoRoT 223992193 and V1481 Ori.

**Key words:** accretion, accretion discs – hydrodynamics – binaries: general – stars: pre-main-sequence.

## 1 INTRODUCTION

About 50 per cent of pre-main sequence (PMS) stars in T associations have been found to be in binary systems (Duchêne & Kraus 2013). Therefore, a large number of PMS stars evolve through accretion from circumbinary discs, and those discs are also the birthplace of planets.

Discs around PMS binary systems were first detected at the same time as discs around single stars in the late 1980s and 1990s through the excess emission in the infrared and at submillimetric wavelengths (Bertout, Basri & Bouvier 1988). The first circumbinary disc to be imaged, using interferometry and then adaptive optics, was that around GG Tau (Dutrey, Guilloteau & Simon 1994; Roddier et al. 1996). Jensen & Mathieu (1997) reported a deficit of near-infrared emission from a few short-period PMS binaries, confirming the theoretical prediction that binaries clear up a cavity in their circumbinary disc (Lin & Papaloizou 1979). Recently, several circumbinary discs have been imaged with the Atacama Large Millimeter Array (ALMA) (Czekala et al. 2015, 2016), and it is believed that ALMA may be able to detect structures in circumbinary

discs produced by the tidal potential of the binary (Ruge et al. 2015).

Numerical simulations have shown that the cavity cleared up by the binary does not prevent accretion on to the stars (Artymowicz & Lubow 1996), as gas spirals in along streams. Those streams connect to circumstellar discs that ultimately are accreted on to the stars. Such circumstellar discs in binary systems have been imaged using adaptive optics (Mayama et al. 2010) and gas in the cavity of the disc around GG Tau has been imaged with ALMA (Dutrey et al. 2014).

Because of the presence of streams and circumstellar discs in the cavity being fed by the circumbinary disc, PMS binaries are complex systems. Emission from such systems is indeed observed to have contribution from various sources and to often be variable, but it is not always understood what produces the emission and the variability (Jensen et al. 2007; Boden et al. 2009; Bary & Petersen 2014; Gillen et al. 2014; Ardila et al. 2015; Messina et al. 2016). Here, we investigate a new source of (variable) emission, namely energy emission from the inner edge of the cavity.

Although the clearing of a cavity by the binary has been extensively studied, the physical conditions in the region of the cavity edge have yet to be elucidated. Here, we show that the action of the gravitational torque from the binary, which truncates the disc and produces the cavity, is accompanied by energy dissipation at the disc's inner edge. This energy dissipation, for the most part,

<footnote n="1"><sup>1</sup>Physics Department, University of Oxford, Keble Road, Oxford OX1 3RH, UK</footnote>

<footnote n="2"><sup>2</sup>Institut d'Astrophysique de Paris, UPMC Univ. Paris 06, CNRS, UMR7095, 98 bis bd Arago, F-75014 Paris, France</footnote>

<footnote n="3"><sup>3</sup>DAMTP, University of Cambridge, Wilberforce Road, Cambridge CB3 0WA, UK</footnote>

<footnote>* E-mail: caroline.terquem@physics.ox.ac.uk (CT); jcbp2@damtp.cam.ac.uk (JCBP)</footnote>
comes from shocks produced by fluid elements which are ultimately accreted by the binary, and that suffer some oscillations in their distance to the centre of mass of the binary resulting in them impacting the disc's inner edge on their way in. In addition, the edge of the cavity may behave as a hard wall which absorbs the radiation from the stars and re-emits it at longer wavelength. These two sources of energy produce an excess of emission which is in the mid-infrared for spectroscopic binaries.

The plan of the paper is as follows. In Section 2, we present hydrodynamic simulations of a disc orbiting around a binary system. We calculate the rate of energy dissipation associated with shocks which are found to be highly localized at the inner edge of the disc. As already noted by previous authors, we find that the rate of accretion on to the binary tends to be suppressed, compared to what would be expected from the value of the surface density just beyond the cavity edge in the absence of the binary, when the disc's aspect ratio decreases below 0.05 or so. However, the energy dissipation at the disc's inner edge is not affected by this suppression. In Section 3, we develop an analysis to elucidate the origin of this energy dissipation. We show that it is associated with the tidal torque which maintains the cavity: angular momentum is transferred from the binary to the disc through the action of compressional shocks and viscous friction. The rate of energy dissipation is given by the binary angular velocity times the torque exerted by the binary, and does not directly depend on the actual accretion rate on to the binary. In Section 4, we present particle simulations to illustrate the dynamics near the disc's inner edge in more detail and the accretion flow on to the binary. These simulations show that, before being accreted, particles are in general accelerated on to a trajectory that sends them back towards the circumbinary disc. Shocks that result at the disc's inner edge together with the action of viscous friction circularize the orbits of the particles and dissipate energy. We find that the energy released during one collision of a unit mass particle with the disc's inner edge is of the order of the potential energy per unit mass there, leading to agreement with the results of the hydrodynamic simulations. When accretion on to the binary is reduced, the particles which end up being accreted undergo more collisions, such that the resulting total energy dissipation rate stays the same. In Section 5, we calculate the spectral energy distribution (SED) of a circumbinary disc taking into account the energy dissipated at the inner edge and also irradiation of the inner edge by the stars. We show that, for tight binaries having a separation ~10 R⊙, emission from the edge completely dominates the SED in the mid-infrared (from 1–4 to 10 µm). In Section 6, we summarize and discuss our results. We conclude that the processes presented in this paper may help to explain the excess of emission in the mid-infrared observed in some PMS binary systems, like CoRoT 223992193 (Gillen et al. 2014) and V1481 Ori (Messina et al. 2016).

## 2 HYDRODYNAMIC SIMULATIONS

We have performed 2D hydrodynamic simulations of a disc orbiting exterior to a binary system in circular orbit using NIRVANA (see, e.g. Ziegler & Yorke 1997). This code has been frequently used to simulate discs interaction with orbiting bodies (e.g. Nelson et al. 2000).

### 2.1 Governing equations

The governing equations solved numerically express the conservation of mass and momentum for a razor thin disc. In a non-rotating

<table><caption>**Table 1.** Parameters of models for which results are described in the text. The first column gives the model label, the second the adopted disc aspect ratio and the third gives the value of $\alpha$. Models A, B and C have the same effective kinematic viscosity while for model D this is a factor of two smaller.</caption><thead><tr><td>Model</td><td>H/r</td><td>α</td></tr></thead><tbody><tr><td>A</td><td>0.1</td><td>0.01</td></tr><tr><td>B</td><td>0.07</td><td>0.02</td></tr><tr><td>C</td><td>0.05</td><td>0.04</td></tr><tr><td>D</td><td>0.05</td><td>0.02</td></tr></tbody></table>

frame with origin at the location of the centre of mass of the binary, these take the form:

$$ \frac{\partial \Sigma}{\partial t} + \nabla \cdot (\Sigma \mathbf{v}) = 0, \qquad (1) $$

$$ \Sigma \left[ \frac{\partial \mathbf{v}}{\partial t} + (\mathbf{v} \cdot \nabla) \mathbf{v} \right] = -\Sigma \nabla \Phi + \nabla \cdot \mathbf{T} + \Sigma \mathbf{F}, \qquad (2) $$

where $\Sigma$ is the surface density, $\mathbf{v}$ is the velocity, $\mathbf{T}$ is the stress tensor which has contributions from the vertically integrated pressure, $\Pi$, and the standard Navier–Stokes viscous stress tensor. The gravitational potential is $\Phi$, and a drag force per unit mass $\mathbf{F}$ is also included. Note that the latter is set to zero in our grid-based simulations but may be retained in our analytic discussion given in Section 3 below. In the simulations presented below, we will use an $\alpha$ prescription for the kinematic viscosity (Shakura & Sunyaev 1973).

### 2.2 Computational setup and initial conditions

We adopt a system of dimensionless units for which the total mass of the binary system $M_1 + M_2$ is unity. In these units, the mass of the primary is $M_1 = 7/12$ and the mass of the secondary is $M_2 = 5/12$. We choose the unit of distance such that the separation, $D$, of the binary components, taken to be in circular orbit, is unity and the unit of time to be such that the binary angular frequency $\omega = \sqrt{G(M_1 + M_2)/D^3}$ is unity. Thus, in these units, the binary orbital period is $2\pi$. In a system of dimensionless cylindrical coordinates $(r, \varphi)$, the radial computational domain is taken to be $[r_i, r_o]$ with $r_i = 0.7$ and $r_o = 10.5$. Thus, the binary orbits interior to the domain. The $\varphi$ domain is $[0, 2\pi]$. The boundary condition is taken to be inflow at $r = r_i$ with the boundary at $r = r_o$ taken to be rigid. We have performed simulations on an equally spaced grid with $N_r$ grid points in the radial direction and $N_\varphi$ grid points in the azimuthal direction. For the simulations presented here, we have adopted $N_r = 384$ and $N_\varphi = 512$ and checked that the same behaviour is obtained adopting twice the resolution in both $r$ and $\varphi$.

The discs are assumed to be locally isothermal with a constant aspect ratio $H/r$. Thus $\Pi = \Sigma H^2 \Omega^2$, where $\Omega$ is the angular velocity. For kinematic viscosity, $\nu$, we adopt an $\alpha$ prescription (Shakura & Sunyaev 1973) with $\alpha$ taken to be constant, i.e. $\nu = \alpha H^2 \Omega$. The value of these parameters for the models considered here are given in Table 1. In order to handle shocks, an artificial viscosity was adopted (see e.g. Stone & Norman 1992).

For initial surface density profile, we adopt $\Sigma = \Sigma_0(r/2.5)^2$ for $r < 2.5$ and $\Sigma = \Sigma_0(2.5/r)$ for $r > 2.5$, where $\Sigma_0$ is a constant. This is taken to be $9 \times 10^{-4}$ in dimensionless units. This form exhibits a cavity, of the type eventually maintained in the simulations, with radius $r \sim 2.5$. Note that, as self-gravity is neglected and the binary is
<figure><img src="image_4.png" /><figcaption>**Figure 1.** The left-hand panel shows a contour plot of $\log_{10}(10^4\Sigma)$, with $\Sigma$ expressed in dimensionless units, for model A, after 208 binary orbits. The right-hand panel represents the rate of dissipation of energy we associate with shocks (see text). In both panels, the yellow circles near the centre indicate the centres of the stellar components and the red crosses indicate the three collinear Lagrangian points. The primary star is on the right with coordinates (0.41, 0), whereas the secondary star is on the left with coordinates (-0.58, 0). The Lagrangian points are $L_3$, $L_1$ and $L_2$ from right to left.</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>**Figure 2.** As in Fig. 1 but for model B after 210 binary orbits.</figcaption></figure>

taken to be in a fixed circular orbit, this value of $\Sigma_0$ is arbitrary. Thus, an arbitrary constant scaling may be applied to the surface density, disc mass and quantities such as the rate of energy dissipation in the disc presented below.

## 2.3 Numerical results

Each of models A–D was run for several hundred orbits. After about 50 orbits a cavity was formed in which the binary resided. This cavity has a radius between 1.5 and 3 and is in general not strictly circular. Inside this cavity, the surface density was dramatically reduced. However, accretion into the binary system occurs through streams of material passing through the exterior Lagrangian points $L_2$ and $L_3$. These features are illustrated in surface density contour plots for models A–D presented in Figs 1–4. These are shown after 208 binary orbits for models A, C and D and after 210 orbits for model B. But note that the pattern looks similar at all times. These figures also show the distribution of the rate of dissipation of energy associated with shocks. This is here traced by evaluating the rate of

energy dissipation associated with artificial viscosity. We see that this is highly localized in radius in regions near the boundary of the cavity, with $1.5 < r < 3$, as would be expected for dissipation associated with tidal torques that maintain the cavity.

In Fig. 5, we show the mass of the disc as a function of time in code units for models A–D. Also shown is the evolution of the mass of a disc with the same kinematic viscosity prescription as models A–C but with the binary absent. The slopes of these curves give the accretion rate into the binary system. The accretion rate averaged over time for $600 < t < 1600$, ($\dot{M}$), for models A–D, is given in Table 2 in units of ($\dot{M}_0$), the latter being the average accretion rate for the disc without the binary system.

For models A and B, this quantity is, respectively, 2.0 and 1.43, indicating that accretion is if anything modestly enhanced by the presence of the binary. This is in line with results given by Shi & Krolik (2015). For model C that has the same kinematic viscosity prescription but the smallest $H/r$, it is reduced to 0.44, indicating that accretion is somewhat suppressed by the presence of the binary. For model D, which has a same value of $H/r$ as model C but a value
<figure><img src="image_6.png" /><figcaption>**Figure 3.** As in Fig. 1 but for model C after 208 binary orbits.</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>**Figure 4.** As in Fig. 1 but for model D after 208 binary orbits.</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>**Figure 5.** The left-hand panel shows the mass of the disc as a function of time in code units for model A (solid line), and model B (dotted line). The dashed line, is for a disc with the same kinematic viscosity as these models but with the binary absent. The right-hand panel shows the mass of the disc as a function of time for model C (solid line) and model D (dotted line). The dashed line is for a disc with the same kinematic viscosity as model C but with the binary absent.</figcaption></figure>

of $\alpha$ that has been reduced by a factor of 2, the averaged accretion rate ratio is reduced to 0.087. This is a factor of ~5 smaller than the value for model C, as compared to an expected reduction of a factor of 2. This situation arises because at early times, measured in units of the characteristic viscous time of the model, the accretion rates are increasing, taken together with the fact that the characteristic

viscous time is a factor of 2 longer for model D as compared to the other models. If we compare the average accretion rate in model C over the time interval $300 < t < 800$ that contracts the time limits by a factor of 2, with the above rate for model D evaluated as an average for $600 < t < 1600$, it is smaller by a factor ~2, as expected. We remark that $\langle \dot{M}_0 \rangle = 1.63 \dot{M}_f$, where the fiducial accretion rate
<table><caption>**Table 2.** Mean accretion rates for $600 < t < 1600$, for the models discussed in the text, expressed in terms of the mean accretion rate evaluated over the same time interval for a disc in which the central binary was absent. Note that in units of the fiducial value $\dot{M}_f = 3\pi\nu\Sigma_0$, with $\nu = \alpha(H/r)^2 D^2\omega$ evaluated for model A (see Table 1), $\langle \dot{M}_0 \rangle = 1.63\dot{M}_f$.</caption><thead><tr><th>Model</th><th>⟨M⟩/⟨M<sub>0</sub>⟩</th></tr></thead><tbody><tr><td>A</td><td>2.0</td></tr><tr><td>B</td><td>1.43</td></tr><tr><td>C</td><td>0.44</td></tr><tr><td>D</td><td>0.087</td></tr></tbody></table>

$\dot{M}_f = 3\pi\nu\Sigma_0$ with $\nu = \alpha(H/r)^2 D^2\omega$ evaluated for model A. In addition, we note that a reduction in the ratio $\langle \dot{M} \rangle/\langle \dot{M}_0 \rangle$, similar to what we find here as $H/r$ is decreased, has been reported in smoothed particle hydrodynamics simulations by Ragusa, Lodato & Price (2016). This corresponds to a suppression of the accretion rate into the binary system as a result of the action of tidal torques that is not seen in the simulations with larger $H/r$.

However, it is important to note that in spite of this suppression, in a steady state we expect that the mass accretion rate at large distances should be the same as that on to the central binary. Although it is impracticable to perform simulations for long enough with an adequate dynamic range in such cases, we expect that when suppression is large, the disc surface density just outside the cavity edge will significantly exceed that at large distances and that in the case of constant kinematic viscosity (see Section 3.3 below) the ratio of these will be approximately equal to $\langle \dot{M}_f \rangle/\langle \dot{M} \rangle$. When considering the rate of energy dissipation due to tidal torques, as in this paper, the value of the surface density just outside the cavity edge, where these torques are applied, is the relevant one. This leads to the fiducial rate of energy dissipation $\langle \dot{M}_f \rangle D^2\omega^2$, rather than $\langle \dot{M} \rangle D^2\omega^2$ (see below). This can be estimated without connecting conditions near the cavity edge to conditions at large distance.

As indicated above, angular momentum is transferred from the binary system to the disc through that action of compressional shocks and viscous friction. This transport maintains the cavity inside which the binary resides. The associated dissipation occurs independently of whether accretion into the binary system is suppressed or not.

The upper panel of Fig. 6 shows the total rate of shock dissipation as a function of time for $r > 1.6$, in code units, for models A and B. The lower panel gives the corresponding quantity for models C and D. We note that the fiducial rate of dissipation associated with tides, indicated here through the rate of energy dissipation associated with artificial viscosity, is $\epsilon_f = \dot{M}_f D^2\omega^2 = 3\pi\nu\Sigma_0 D^2\omega^2$, with $\nu = \alpha(H/r)^2 D^2\omega$ evaluated for model A, and is equal to $8.5 \times 10^{-7}$. For models A–C, the average dissipation rate is $\sim 2\epsilon_f$ although most of this is concentrated in localized spikes. In model D, for which the effective viscosity was a factor of 2 smaller, the average rate of dissipation $\sim 0.5\epsilon_f$. These results indicate that, unlike the accretion rate, the tidal dissipation is not greatly inhibited for small $H/r$. This occurs at approximately the same rate in all cases and this rate is approximately $\langle \dot{M}_f \rangle D^2\omega^2$. Note that, as indicated above, the reduction by a factor of 4 rather than a factor of 2 in the case of model D as compared to model C can be accounted for by the increase in the rate of evolution that occurs as a function of time for models A–C.

<figure><img src="image_9.png" /><figcaption>**Figure 6.** The upper panel shows the total rate of shock dissipation as a function of time for $r > 1.6$, in code units, for model A (solid line) and model B (dashed line). The lower panel gives the corresponding plot for models C (solid line) and D (dashed line). Note that, in code units, the fiducial rate of dissipation given by $\dot{M}_f D^2\omega^2 = 3\pi\nu\Sigma_0 D^2\omega^2$, with $\nu = \alpha(H/r)^2 D^2\omega$ evaluated for model A (see Table 2), is $8.5 \times 10^{-7}$.</figcaption></figure>

Fig. 6 indicates that the rate of shock dissipation is variable on a time-scale on the order of a few $\omega^{-1}$.

# 3 THE ORIGIN OF TIDAL DISSIPATION

We noted above that the rate of energy dissipation associated with the tidal interaction is approximately given by $\langle \dot{M}_f \rangle D^2\omega^2$, independently of whether accretion on to the central binary is significantly suppressed. We now explore a simple picture of this dissipation.

## 3.1 Rate of change of the Jacobi invariant

We begin by considering an infinitesimal fluid element of mass $m$ occupying a varying area $A$ with surface density $\Sigma$, going on to derive an equation governing the evolution of the Jacobi invariant associated with it. The constancy of $m$ is consistently implied by the continuity equation (1). From equation (2) we obtain

$$m \frac{Dv}{Dt} = \frac{D}{Dt} \int_A \Sigma v \, dA = -m \nabla \Phi + F_T + mF. \quad (3)$$

Here

$$F_T = \int_A (\nabla \cdot T) \, dA = \int_{\Gamma} T \cdot \hat{n} \, dl, \quad (4)$$

where the $i$th component of $T \cdot \hat{n}$ is $T_{ij} \hat{n}_j$, and the second integral is obtained using Green's theorem. It is taken over the boundary $\Gamma$ of the fluid element, with $\hat{n}$ being the unit vector normal to the
boundary and pointing outward, and so it is derived from the stresses due to exterior fluid. Thus, it could act like a drag force per unit mass. For the sake of brevity, we set $f = F + F_T/m$.

Using the fact that for our setup the gravitational potential $\Phi$ is a function of $\varphi$ and $t$ through the combination $\varphi - \omega t$, we obtain an equation for the rate of change of the Jacobi invariant $E - \omega J$, where $E$ is the kinetic plus potential energy and $J$ is the angular momentum of the fluid element. We begin by noting that we can obtain an equation for the rate of change of the kinetic energy of a fluid element from equation (2) under the form:

$$ \frac{1}{2} \frac{D}{Dt} \int_A \Sigma v^2 \, dA = -mv \cdot \nabla \Phi + mv \cdot f - \epsilon, \quad (5) $$

where

$$ \epsilon = \int_A (\nabla v : T) \, dA, \quad (6) $$

with $\nabla v : T = T_{ij}(\partial v_i / \partial x_j)$, represents the rate of doing work by internal stresses. This contains the rate of dissipation due to viscosity when that is present. We may also write equation (5) as

$$ \frac{DE}{Dt} = \frac{D}{Dt} \int_A \Sigma \left( \frac{1}{2} v^2 + \Phi \right) dA = m \frac{\partial \Phi}{\partial t} + m v \cdot f - \epsilon. \quad (7) $$

From the $\varphi$-component of equation (3), we obtain:

$$ \frac{DJ}{Dt} = m \frac{D(r v_{\varphi})}{Dt} = -m \frac{\partial \Phi}{\partial \varphi} + m r f \cdot \dot{\varphi}, \quad (8) $$

where $\dot{\varphi}$ is the unit vector in the azimuthal direction. Multiplying equation (8) by $\omega$ and subtracting the resulting equation from equation (7) we obtain

$$ \frac{D(E - \omega J)}{Dt} = m (v - r \omega \dot{\varphi}) \cdot f - \epsilon. \quad (9) $$

Note that this form applies in the inertial frame.

## 3.2 Relation between energy dissipation and binary torque for an oscillating fluid element

Suppose a fluid element near some inner disc edge undergoes cyclic behaviour such that, starting from a circular orbit, it moves inwards and as a result of interacting with the binary and surrounding fluid it then returns to the same edge where it is circularized. Then, the total change in Jacobi constant is zero and from equation (9) we have

$$ \int (-mv \cdot f + \epsilon) \, dt = \int -mr\omega\dot{\varphi} \cdot f \, dt = \omega \int \langle T \rangle dt. \quad (10) $$

Here, the time integrals are over the duration of the process. Equation (10) states that the mean rate of energy dissipation through the action of the force $f$, as measured in the inertial frame, together with the mean rate of energy dissipation by internal stresses, is the product of the binary orbital frequency and time average of the torque exerted by the binary, $\langle T \rangle$, taken over the duration of the process. That is the time-averaged torque exerted by the binary follows from equation (8). As by assumption, the fluid elements angular momentum is not changed through this process, the torque is communicated to material in contact.

As we know accretion into the binary system is ongoing, the fluid element will subsequently move away from the edge into the binary. A similar discussion could be applied to such a final stage if one knows the change in $E - \omega J$, the Jacobi invariant, during this process. As it is dynamical, we shall assume this is small so that equation (10) applies throughout over extended time periods.

## 3.3 A standard accretion disc model

We now connect the above discussion to standard accretion disc modelling that performs azimuthally averaging. Within this framework, we consider a standard steady accretion disc with an inner boundary at which the above torque is applied by the binary and through which there is a constant accretion rate $\dot{M}$. Angular momentum conservation gives

$$ \frac{d}{dr} \left( |\dot{M}| r^2 \Omega + 2\pi\nu \Sigma r^3 \frac{d\Omega}{dr} \right) = -\langle T \rangle_1, \quad (11) $$

where $\Omega$ is the local angular velocity and $\langle T \rangle_1$ is the applied time-averaged torque per unit radius. If we suppose that the torque is applied in a thin layer at an inner edge, integrating through this layer we get on the exterior + side where $r = r_+$:

$$ 2\pi\nu \Sigma r^3 \left. \frac{d\Omega}{dr} \right|_+ = -\langle T \rangle, \quad (12) $$

where $\langle T \rangle$ is again the time-averaged total applied torque and $|_+$ indicates evaluation at $r = r_+$. However, integrating equation (11) in the main body of the disc we have

$$ |\dot{M}| r^2 \Omega + 2\pi\nu \Sigma r^3 \frac{d\Omega}{dr} = \text{constant} = |\dot{M}| r^2 \Omega|_+ - \langle T \rangle. \quad (13) $$

This equation determines the $\Sigma$ profile and rate of viscous energy dissipation as in a standard disc model. Note that, as compared to the standard solution with a zero torque applied at the boundary, there is the additional $\langle T \rangle$ term. This is responsible for an additional viscous dissipation rate in the disc amounting to $\Omega|_+ \langle T \rangle$. This is comparable to, but significantly less than, the tidal dissipation rate estimated above as $\omega \langle T \rangle$. This is because the dissipation envisaged there cannot all be associated with a simple shear viscosity. The difference $(\omega - \Omega|_+) \langle T \rangle$ must represent the dissipation rate associated with compressional shocks (Papaloizou & Pringle 1977).

In addition, equation (12) implies that the rate of dissipation due to tidal torques acting on a Keplerian disc is $\omega \langle T \rangle = 3\pi\nu \Sigma r^2 \omega \Omega|_+$. This is similar to the fiducial rate seen in the simulations and given by

$$ \epsilon_f = \dot{M}_f D^2 \omega^2 = 3\pi\nu \Sigma_0 D^2 \omega^2, \quad (14) $$

where $\Sigma_0$ should be thought of as being the value of $\Sigma$ just outside the cavity edge and $\nu$ is taken to be $\alpha(H/r)^2 D^2 \omega$.

However, when tidal torques are effective at restricting accretion, as noted above, in a steady state the surface density just exterior to the gap edge will significantly exceed the surface density at large distances. For the simple example where $\nu$ is constant, equation (13) implies that at larger distances $\Sigma \to \Sigma_\infty$, where

$$ 3\pi\nu \Sigma_\infty = |\dot{M}|. \quad (15) $$

Then

$$ \frac{\Sigma_+}{\Sigma_\infty} = \frac{\langle T \rangle}{|\dot{M}| r^2 \Omega|_+}. \quad (16) $$

Thus, in this case, the mean energy dissipation rate due to tidal torques will exceed $|\dot{M}| \omega r^2 \Omega|_+$ by a factor approximately equal to the steady state accretion rate estimated adopting the value of the surface density just exterior to the cavity edge to the actual steady state accretion rate.

# 4 PARTICLE SIMULATIONS

In order to gain more insight into the nature and magnitude of the tidal dissipation occurring near the cavity boundary, we use a
particle code to simulate the accretion flow around the binary. This is based on the fact that, apart from where shocks occur, the pressure force is small compared to the gravitational force (or, equivalently, the sound speed is small compared to the Keplerian velocity), which implies, to a good approximation, that the motion should be ballistic (e.g. Shi et al. 2012).

## 4.1 Computational setup

The particles are set down at some distance $r$ from the binary centre of mass $O$ and are initiated the local Keplerian angular velocity $\Omega_K = [G(M_1 + M_2)/r^3]^{1/2}$ in the inertial frame centred on $O$. The particles move under the gravitational potential due to the stars. To force them to move towards the cavity, we add a drag force per unit mass $F$:

$$F = -\eta |a| \frac{v}{|v|}, \quad (17)$$

where $\eta < 1$ is a free parameter, and $a$ is the acceleration in the inertial frame centred on $O$. This form of the drag force ensures that it acts opposite to the velocity and, as the acceleration results mainly from gravitational forces, is $\eta$ times smaller than the gravitational force per unit mass. In practice, we only apply $F$ to particles beyond the inner radius of the cavity.

This drag force mimics the action of disc viscosity, modelled using the $\alpha$ prescription above, but more realistically resulting from turbulent processes arising from the magnetorotational instability. In a realistic situation, this would enable particles to lose angular momentum and get accreted on to the binary.

The parameter $\eta$ defined through equation (17) can be approximated by $\eta \sim |v_r|/v_\omega$. In a standard accretion disc, $|v_r| \sim v/r$ where $v = \alpha c_s H$ is the kinematic viscosity, with $c_s$ being the sound speed. Using $c_s = (H/r)v_\omega$, we obtain $\eta \sim \alpha(H/r)^2$, which gives $\eta \simeq 10^{-4}$ for models A–C and $\eta = 5 \times 10^{-5}$ for model D. However, at the disc’s inner edge, the length-scale is given by $H$ rather than $r$, and the magnitude of the radial velocity in low surface density regions may be better approximated by $|v_r| \sim v/H$. This would then lead to $\eta \sim \alpha(H/r)$, i.e. $\eta \simeq 10^{-3}$ for all the models A–D. Given that there is some uncertainty in the value of $\eta$ that should be used, we have varied $\eta$ between $10^{-5}$ and a few times $10^{-3}$.

Before being accreted by one of the stars, a particle may be accelerated on to a trajectory that sends it back to the circumbinary disc. When that happens, we assume that shocks recircularize the orbit of the particle and accordingly we reset its velocity to be the Keplerian value.

## 4.2 Numerical results

The accretion flow is shown in Fig. 7. As for the hydrodynamic simulations, we have $M_1 = 7/12$, $M_2 = 5/12$ and distances are given in units of the binary separation ($D=1$). The particles were initially set up on circular orbits with $r$ between 2 and 4.4. The small black crosses represent the position of the particles and the large red crosses the position of the two stars, with the primary on the right-hand side, as in Figs 1–4. The red triangles on the figure indicate the position of the Lagrangian points $L_3$, $L_1$ and $L_2$ from right to left. The value $\eta = 10^{-3}$ was adopted for the calculation shown in Fig. 7, but the accretion flow is very similar for smaller values of $\eta$. Although it takes longer for particles to drift inwards when $\eta$ is reduced, the overall structure of the flow is not affected.

The rotation of the binary is anticlockwise. Because the binary is rotating faster than the particles, in the corotating frame the motion

<figure><img src="image_10.png" /><figcaption>**Figure 7.** Accretion flow around a binary system in the corotating frame simulated with a particle code. The mass of the primary (right red cross) is 7/12 and that of the secondary (left red cross) is 5/12. The separation is $D=1$. The red circles indicate the outer edge of the circumstellar discs, whose radius is about $0.4D$. The red triangles indicate the position of the Lagrangian points $L_3$, $L_1$ and $L_2$ from right to left.</figcaption></figure>

of the particles is clockwise. Particles that get close enough to one of the stars are captured on to orbits around the star, and circumstellar discs form. Tidal truncation of these discs by the other star limits their outer radii, which are shown as the red circles on the figure. We find here that the radius $r_d$ of these discs is 0.4 times the binary separation, close to the value of one-third predicted by theoretical studies (Paczynski 1977; Papaloizou & Pringle 1977). Particles within these discs have been removed from the figure.

As can be seen from the plot, the mean radius of the cavity cleared out by the binary is about 1.9, close to the value of 2 expected from theory (see Lubow & Artymowicz 2000 and references therein). The details of the flow beyond this radius cannot be captured by the particle code, as hydrodynamical effects become important there.

The flow shown in Fig. 7 is very similar to that obtained from our hydrodynamical simulations (see, e.g. Fig. 4) and also those of, e.g. Hanawa, Ochi & Ando (2010), who considered a binary with a mass ratio similar to that adopted here.

## 4.3 Collisions at the edge of the cavity

As has been noted in previous studies of circumbinary discs (McFadyen & Milosavljevic 2008; Shi & Krolik 2015), material entering the inner cavity does not in general fall directly on to one of the circumstellar discs. A particle originating from the edge of the cavity may be accelerated by the gravitational torque on a trajectory that propels it back towards the circumbinary disc. Shocks may then occur that result in circularization of the trajectory, the excess kinetic energy being released and radiated away.

We model this process using the particle code in the following way. After a particle has entered the cavity, which is assumed to have a radius $r_{in}$, if it is flung back out to the disc, i.e. is found with $r \simeq r_{in}$ again, its velocity is reset to the Keplerian value at this radius.
<figure><img src="image_11.png" /><figcaption>**Figure 8.** The radius *r* (in units of the binary separation *D*) of a particle moving inside the cavity as a function of time (in units of *ω*⁻¹) for *η* = 4 × 10⁻³ (upper left panel), 10⁻³ (upper right panel), 10⁻⁴ (lower left panel) and 10⁻⁵ (lower right panel). In this calculation, the velocity of the particle is reset to the Keplerian value if it is flung back to the disc and reaches *r* = 2.4 after having gone down to *r* = 1.6 (these values of *r* are indicated by the dotted lines). The trajectories that are shown are typical but note that trajectories may differ significantly from one particle to another. As *η* is reduced, the particle undergoes more collisions before crossing the binary orbit.</figcaption></figure>

The particle then resumes drifting towards the cavity. We note $\delta E_K$ the kinetic energy which is released, and which is the difference between the kinetic energy *E* of the particle when it reaches *r*<sub>in</sub> and the kinetic energy *E*<sub>in</sub> of the same particle undergoing Keplerian motion there. To evaluate $\delta E_K$, we need in principle to know the value of *r*<sub>in</sub>. However, if we vary this radius between 1.5 and 2.5 for example, *E*<sub>in</sub> is found to vary only by a factor of 0.6. Since *E* itself varies only by a factor of order unity when the particle goes from *r* = 1.5 to *r* = 2.5, $\delta E_K$ does not depend very sensitively on *r*<sub>in</sub>. In addition, as can be seen from Fig. 7, the cavity is not actually circular. Therefore, it is not essential for the argument to know *r*<sub>in</sub> precisely.

If the distance between the particle and one of the stars becomes smaller than the radius $r_d$ of the circumstellar discs, it is then supposed to be accreted by the star. We fix $r_d = 0.4$.

The radius of a particle moving inside the cavity as a function of time is shown in Fig. 8. It moves slowly inwards in an orbit with

steadily increasing eccentricity. It is then flung outwards ultimately being reset on a circular orbit at larger radii and the process then repeats until the particle is accreted by the binary. As can be seen, the particle undergoes both small and large excursions, indicating that there should be a corresponding range of shock strengths present as we found to be the case in our hydrodynamic simulations (see Fig. 6). In the calculations presented in Fig. 8, the velocity of the particle is reset to the Keplerian value if it is flung back to the disc and reaches $r = 2.4$ after having gone down to $r = 1.6$. Therefore, circularization proceeds only after a large excursion occurs. Because we have chosen $r = 2.4$ rather arbitrarily as the radius where circularization takes place, we cannot calculate the real time between two collisions with this model. However, we can calculate the number of times a particle collides with the disc’s inner edge before being accreted. For a fixed value of $\eta$, this number is obtained by averaging the results corresponding to many trajectories calculated for different initial conditions (for a given $\eta$, the trajectories
for two different particles may be significantly different, which is why an averaging is required).

We find that, for $\eta = 4 \times 10^{-3}$, on average particles get flung back out to the disc having undergone one or two large radial excursions. For higher values of $\eta$, the drift force is rather large and particles get accreted straight away without colliding with the disc. When collisions occur, particles hit the inner edge of the circumbinary disc at $r \sim r_{in}$ with a velocity which is on the order of the Keplerian velocity at this radius but with a different orientation. Therefore, $\delta E_K \sim E_{in}$. If the mass of particles drifting through the circumbinary disc's inner edge per unit time is $\dot{M}$, then the total kinetic energy released per unit time is $\Delta E_K \sim G(M_1 + M_2)\dot{M}/r_{in}$. This is essentially the same as given by equation (14) if we identify $\dot{M} = \dot{M}_f$, the fiducial accretion rate which in steady state is $3\pi\nu\Sigma_0$ (see Section 2.2 and Table 2), and it corresponds to what is found in our hydrodynamic simulations in the cases of models A and B.

If $\eta$ is reduced to $10^{-4}$ or $10^{-5}$, particles move more slowly inwards and undergo more large radial excursions leading to more energy dissipated per particle as it is accreted. Suppose this is $\mathcal{N}\delta E_K$ with $\mathcal{N} > 1$. In our particle simulations, all the particles are forced inwards by the drag force given by equation (17), and therefore they all end up being accreted by one of the stars. Therefore, the total kinetic energy released per unit time will be $\Delta E_K \sim \mathcal{N}G(M_1 + M_2)\dot{M}/r_{in}$.

However, we know from the more realistic hydrodynamic simulations presented in Section 2 and from the results by Ragusa et al. (2016) that the accretion rate on to the binary is actually reduced below the fiducial value $\dot{M}_f$ as $H/r$ and $\alpha$ are reduced, which here corresponds to reducing $\eta$. On the other hand, as seen in the hydrodynamic simulations and confirmed by the analysis done in Section 3, the rate of energy dissipation at the disc's inner edge is determined by the tidal torque which maintains the cavity and is expected to be given by $\Delta E_K \sim G(M_1 + M_2)\dot{M}_f/r_{in}$. This leads us to identify $\mathcal{N} \sim \dot{M}_f/\dot{M}$. In other words, the total energy which is dissipated, being tapped from the binary, depends only on the size of the cavity and the surface density at the inner edge, not on the actual accretion rate on to the binary.

The above discussion implies that when $H/r$ and $\alpha$ (or equivalently $\eta$) are decreased, as particles are accreted more slowly, each of them has to dissipate more energy, and therefore undergo more collisions, which is exactly what is seen in the results presented in Fig. 8. Thus, in the case $\eta = 10^{-4}$ or $10^{-5}$, the tidal torques and related dissipation rate required to maintain the cavity (see equation 14) can be maintained with an accretion rate on to the binary that is significantly less than the fiducial rate $\dot{M}_f$. This would correspond to models C and D considered above where the accretion rate on to the binary is significantly less than $\dot{M}_f$.

Fig. 9 shows a sketch of the trajectories of the particles in the frame corotating with the binary near the edge of the circumbinary disc. Particles that make it through one of the Lagrangian points enter into orbital motion around one of the stars, whereas other particles are propelled back to the disc after having drifted inside the cavity. As can be seen from this figure, the shocks tend to be localized at some specific locations rather than be uniformly distributed in azimuth. Similar features were observed in the hydrodynamic simulations presented in Section 2, although the location of the shocks changes with time.

## 5 DISC'S SPECTRAL ENERGY DISTRIBUTION

In this section, we calculate the contribution from the collisions described above to the disc's SED. We model the circumbinary

<figure><img src="image_12.png" /><figcaption>**Figure 9.** Sketch of the trajectories of the particles in the frame corotating with the binary near the edge of the circumbinary disc. The black arrow indicates the (anticlockwise) direction of the binary's rotation. The green lines show trajectories of particles that drift inwards but are propelled back to the disc and collide at the inner edge. The red lines show trajectories of particles that go through $L_2$ or $L_3$. They arrive at either of these points with a small speed and are subsequently accelerated on an anticlockwise trajectory around one of the stars.</figcaption></figure>

<figure><img src="image_13.png" /><figcaption>**Figure 10.** Schematic view of the circumbinary disc's inner edge and cavity.</figcaption></figure>

disc's inner edge as a 'wall' with radius $r_{in}$ and vertical extension $2H_{in}$, as illustrated in Fig. 10.

At radii $r > r_{in}$, the disc is heated by irradiation from the central stars and by dissipation of gravitational energy due to accretion. These sources contribute the energies per unit time per unit surface area $\sigma T_{irr}^4$ and $\sigma T_{acc}^4$, respectively, where $\sigma$ is the Stefan-Boltzmann constant and (Lynden-Bell & Pringle 1974; Friedjung 1985)

$$T_{irr}^4 = \frac{2T_*^4}{3\pi} \left(\frac{R_*}{r}\right)^3, \qquad (18)$$

$$T_{acc}^4 = \frac{3GM_*\dot{M}_f}{8\pi\sigma r^3}. \qquad (19)$$

In equation (19), $\dot{M}_f$ is what has been defined in Section 2 as the fiducial accretion rate, which in steady state is $3\pi\nu\Sigma$. Here, we have assumed $r > R_*$ and replaced the two stars with luminosities, radii and masses $L_{1,2}$, $R_{1,2}$ and $M_{1,2}$, respectively, by one 'equivalent' star located at the centre of mass with temperature $T_*$ given by $\sigma T_*^4 = L_1/(4\pi R_1^2) + L_2/(4\pi R_2^2)$, radius $R_*$ such that
$\sigma T_*^4 = (L_1 + L_2)/(4\pi R_*^2)$ and mass $M_* = M_1 + M_2$. Therefore, the surface temperature $T_{\text{disc}}$ of the disc at $r > r_{\text{in}}$ can be calculated from

$$T_{\text{disc}} = (T_{\text{irr}}^4 + T_{\text{acc}}^4)^{1/4}. \qquad (20)$$

Note that equation (18) assumes that the disc is of constant aspect ratio, i.e. $H/r$ is constant. If the disc were flared, the temperature due to irradiation would be larger at larger radii, but since we are only concerned here with the SED near the disc's inner edge, we will ignore such complications. We have also assumed in using equation (18) that the entire star could be seen from any point at the surface of the disc. This would not be the case if $H_{\text{in}}$ were comparable to or larger than $R_*$, in which case the temperature due to irradiation would be significantly smaller than that given by equation (18). As this would only reinforce the conclusions presented below, we will also ignore this effect.

At $r = r_{\text{in}}$, the disc's inner wall is heated by irradiation from the central stars and by the energy released through the action of tidal torques. That release occurs through shock dissipation and viscous friction in the hydrodynamic model and collisions that result in circularization just beyond the disc edge in the particle model. These sources, respectively, contribute the energies per unit time per unit surface area $\sigma T_{\text{in,irr}}^4$ and $\sigma T_{\text{in,coll}}^4$, with

$$T_{\text{in,irr}}^4 = \frac{L_1 + L_2}{4\pi\sigma r_{\text{in}}^2}, \qquad (21)$$

$$T_{\text{in,coll}}^4 = \frac{\Delta E_K}{\sigma S_{\text{in}}}, \qquad (22)$$

where $S_{\text{in}} = 4\pi f r_{\text{in}} H_{\text{in}}$ is the surface of the wall over which collisions occur. The factor $f < 1$ accounts for the fact that shocks are not distributed uniformly in azimuth. From the point of view of the particle model, we have assumed that the particles that are propelled back to the disc collide and are thermalized at $r = r_{\text{in}}$. This is consistent with the fact that the hydrodynamic simulations presented in Section 2 (see also Shi & Krolik 2015) show that the edge of the cavity is very sharp, with the mass density dropping by several orders of magnitude over a few scaleheights. Therefore, the particles that drift in the cavity and are propelled back to the disc do not penetrate significantly through the disc.

At locations where shocks occur, i.e. over a surface of area $S_{\text{in}}$, the temperature $T_{\text{in}}$ at the disc's inner edge can be calculated from

$$T_{\text{in}} = (T_{\text{in,irr}}^4 + T_{\text{in,coll}}^4)^{1/4}, \qquad (23)$$

whereas at locations where there are no shocks, i.e. over a surface area $(1-f)S_{\text{in}}/f$, it is given by $T_{\text{in,irr}}$.

Hereafter we will take

$$\Delta E_K = \frac{G(M_1 + M_2)\dot{M}_f}{r_{\text{in}}}, \qquad (24)$$

as this gives the correct order of magnitude for the energy released by the collisions (see previous sections and equation 14). Here again, $\dot{M}_f$ is the fiducial accretion rate, as this is what sets the level of energy dissipation at the disc inner edge, even when the accretion rate on to the binary is reduced.

The flux received from the disc at a wavelength $\lambda$ is then $\lambda F_{\lambda} = (\cos i/L^2)E_{\lambda}$, with

$$E_{\lambda} = \int_{r_{\text{in}}}^{r_{\text{disc}}} \lambda B_{\lambda} [T_{\text{disc}}(r)] 2\pi r dr + 4\pi r_{\text{in}} H_{\text{in}} \lambda [f B_{\lambda}(T_{\text{in}}) + (1-f)B_{\lambda}(T_{\text{in,irr}})], \qquad (25)$$

where $B_{\lambda}$ is the Planck function, $L$ is the distance to the binary system and $i$ is the angle between the light of sight and the perpendicular to the orbital plane. In the calculations below, we will adopt $r_{\text{in}} = 2D$, $H_{\text{in}} = 0.05 r_{\text{in}}$ and $r_{\text{disc}} = 10^3 r_{\text{in}}$ (the value of $r_{\text{disc}}$ is not important as long as it is large enough to capture all the contribution to the flux at the wavelengths we are interested in).

Fig. 11 shows the flux $\lambda F_{\lambda}$ as a function of $\lambda$ for the parameters derived for the eclipsing binary CoRoT 223992193 (Gillen et al. 2014): $L_1 = 0.275 L_{\odot}$, $L_2 = 0.182 L_{\odot}$, $R_1 = 1.3 R_{\odot}$, $R_2 = 1.1 R_{\odot}$, $M_1 = 0.7 M_{\odot}$, $M_2 = 0.5 M_{\odot}$, $D = 0.05$ au, $i = 85^{\circ}$ and $L = 756$ pc. Curves have been calculated for $\dot{M}_f = 10^{-7}$, $10^{-8}$ and $10^{-9} M_{\odot} \text{yr}^{-1}$ and for $f=1$ and $f=0.1$. As pointed out in Section 4, $f$ is smaller than 1, but we have done a calculation with $f=1$ to show that the results do not depend very sensitively on the value of $f$. This is because, although $T_{\text{in,coll}}$ given by equation (22) increases slightly when $f$ decreases (as $f^{-1/4}$), the surface that re-radiates the energy dissipated through shocks is reduced. The net effect is therefore not very significant.

Fig. 11 shows that at wavelengths between 1–4 and 10 µm, the SED is completely dominated by emission from the disc's inner edge. This emission is itself dominated by the shocks produced by fluid particles being propelled back to the disc for large fiducial accretion rates exceeding about $10^{-8} M_{\odot} \text{yr}^{-1}$. The maximum of the excess moves towards smaller wavelengths when the surface over which shocks occur is decreased. For smaller accretion rates, irradiation by the stars is the main source of emission from the disc's inner edge. For $\dot{M}_f < 10^{-9} M_{\odot} \text{yr}^{-1}$, the SED is indistinguishable from that corresponding to $\dot{M}_f = 10^{-9} M_{\odot} \text{yr}^{-1}$, as the contribution to the SED from accretion becomes negligible.

## 6 SUMMARY AND DISCUSSION

In this paper, we have shown that energy has to be dissipated in the circumbinary disc for the tidal torque to maintain the cavity. Energy dissipation happens at the inner edge of the disc through compressional shocks produced by fluid that is propelled back to the disc after drifting inwards into the cavity. The rate of energy dissipation does not depend on the actual accretion rate on to the binary, only on what we call the fiducial accretion rate through the disc, $\dot{M}_f$, which is determined by the state variables in the disc just beyond the cavity edge. This would be the accretion rate on to the central object if it was a single star rather than a binary and the material near the cavity edge was part of a putative exterior steady state disc. Energy dissipation occurs at the fiducial rate even when the accretion rate on to the binary is reduced below $\dot{M}_f$ (which happens when the aspect ratio and viscosity of the disc is reduced). The rate of energy dissipation is typically $G(M_1 + M_2)\dot{M}_f/r_{\text{in}}$, where $r_{\text{in}}$ is the inner radius of the circumbinary disc.

Hydrodynamic simulations [see also the magnetohydrodynamic (MHD) simulations performed by Shi & Krolik (2015)] show that the edge of the circumbinary disc is very sharp. Therefore, emission of the energy released by the shocks is localized at the edge. In addition, this edge is irradiated by the central stars and this also contributes to the SED of the system.

In a tight PMS binary with a separation of ~10 $R_{\odot}$, emission from the disc's inner edge completely dominates the SED in the mid-infrared, i.e. for $\lambda \sim 1-4$ to 10 µm. For large values of the accretion rate $\dot{M}_f$ exceeding about $10^{-8} M_{\odot} \text{yr}^{-1}$, the emission predominantly comes from the shocks derived from fluid elements being propelled back towards the disc. For lower accretion rates, irradiation from the central star becomes the main source of emission from the inner edge.
<figure><img src="image_14.png" /><figcaption>**Figure 11.** Flux $\lambda F_{\lambda}$ (in units $10^{-13}$ erg cm$^{-2}$ s$^{-1}$) as a function of $\lambda$ (in µm). The dotted lines represent the disc without inner wall, i.e. only the integral is taken into account in calculating $E_{\lambda}$ in equation (25). The dashed lines show the contribution of the inner wall when only irradiation is present, i.e. equation (25) without the integral and $T_{\text{in}} = T_{\text{in, irr}}$. The dot-dashed lines show the contribution of the inner wall when only tidal torques are present, i.e. equation (25) with only the $fB_{\lambda}(T_{\text{in}})$ term and $T_{\text{in}} = T_{\text{in, coll}}$. The solid lines represent the total flux, i.e. disc and inner wall with both irradiation and tidal torques. The plots on the left- and right-hand sides correspond to $f=1$ and $f=0.1$, respectively. The upper, middle and lower plots correspond to $\dot{M}_{f} = 10^{-7}$, $10^{-8}$ and $10^{-9}$ M$_{\odot}$ yr$^{-1}$, respectively. The curves correspond to the parameters of the eclipsing binary CoRoT 223992193 (Gillen et al. 2014): $L_{1} = 0.275$ L$_{\odot}$, $L_{2} = 0.182$ L$_{\odot}$, $R_{1} = 1.3$ R$_{\odot}$, $R_{2} = 1.1$ R$_{\odot}$, $M_{1} = 0.7$ M$_{\odot}$, $M_{2} = 0.5$ M$_{\odot}$, $D = 0.05$ au, $i = 85^{\circ}$ and $L = 756$ pc.</figcaption></figure>

When shocks dominate, we expect the SED to display some variability, as the shock dissipation rate varies with time (see Fig. 6) on a time-scale of a few $\omega^{-1}$, where $\omega$ is the angular velocity of the binary, i.e. comparable to the binary period. Even when irradiation dominates, we would expect some variability on a time-scale comparable to the binary period as the distance between the edge of the disc and the stars varies as the stars move around their centre of mass.

## Application to CoRoT 223992193

CoRoT 223992193 is a double-lined PMS eclipsing binary which was discovered by Gillen et al. (2014). Its age has been estimated

to be between 3.5 and 6 Myr. The SED of this system presents a mid-infrared excess between 4 and 10 µm. Gillen et al. (2014) found that this excess could not be produced by a circumbinary disc unless it had an accretion rate of $10^{-7}$ M$_{\odot}$ yr$^{-1}$, four orders of magnitude higher than the accretion rate derived from the $H\alpha$ emission from the stars. It was then proposed that the excess may be due to thermal emission from dust present in the vicinity of the two stars. As little as $10^{-13}$ M$_{\odot}$ of dust would be needed to produce the observed flux. This dust may also be responsible for the out of eclipse variability observed in this system (Terquem, Sørensen-Clark & Bouvier 2015).

Here, we comment that the excess emission observed for CoRoT 223992193 is very similar to what would be expected from the inner edge of a circumbinary disc. In this system, emission would not be
dominated by the shocks as the accretion rate on to the stars is only
$10^{-11}$ M$_{\odot}$ yr$^{-1}$ and the fiducial accretion rate through the disc is
probably not that much larger (as $H/r$ is not expected to be smaller
than 0.05 or so in a disc around a PMS binary). Emission would
then be produced by irradiation from the stars. With $H_{\text{in}} = 0.05$, the
disc's inner edge would produce a flux with a maximum of about
$3.5 \times 10^{-13}$ erg cm$^{-2}$ s$^{-1}$ at a wavelength between 3 and 4 µm,
and a flux about 3.5 times smaller at 10 µm (see the lower plots
on Fig. 11). The maximum value of the flux would be increased up
to $10^{-12}$ erg cm$^{-2}$ s$^{-1}$ if we adopted $H_{\text{in}} = 0.15$ instead. Such an
enhanced value of the disc's aspect ratio near the edge of the cavity
is consistent with 3D MHD simulations (Shi & Krolik 2015) and
may be produced by magnetic fields that get relatively stronger in
the cavity. The peak of the excess measured for CoRoT 223992193
is about $3 \times 10^{-12}$ erg cm$^{-2}$ s$^{-1}$ at 4 µm and decreases by the same
factor as in our model at 10 µm (see fig. 13 of Gillen et al. 2014).
Given that we have calculated irradiation of the disc's inner edge in a
very crude way, our results are close enough to the observations that
emission from the disc's inner edge appears as a serious candidate
for producing the observed excess. We note that, as the system is
seen with an inclination $i = 85^{\circ}$ (Gillen et al. 2014) which, although
large, is such that $90^{\circ} - i$ in radians exceeds $H/r$, accordingly the
inner rim of the disc would not be obscured by the disc itself.

Application to V1481 Ori

V1481 Ori is a spectroscopic binary with a period of about 4.4 d in the Orion nebula cluster. It is associated with a strong infrared excess which suggests the presence of an accreting circumbinary disc (Messina et al. 2016 and references therein). Recent spectroscopic and photometric observations performed by Messina et al. (2016) in the *V* band at ~0.7 µm and in the *I* band have shown that the luminosity ratio of the stellar components varies in a way that can be explained by the presence of a hotspot on the secondary, which is speculated to be produced by accretion from a circumstellar disc. These observations also show that the minimum value of the luminosity ratio is 30 per cent larger than what would be expected in the case of no accretion. Messina et al. (2016) conclude that there has to be an additional source of emission that is visible at all phases.

Here, we suggest that this excess of energy may be produced by
the inner edge of the circumbinary disc. Given the short wavelength
at which this excess occurs, it would have to be produced by shocks
localized over a small fraction of the surface at the inner edge,
as illustrated on the upper right plot of Fig. 11. The peak of the
excess on this plot is at about 1 µm, so that there is significant
emission in the bands in which Messina et al. (2016) have made
the observations. This model requires a rather large accretion rate
in the circumbinary disc which, as pointed out above, is suggested
by the observations.

ACKNOWLEDGEMENTS

It is a pleasure to thank Julian Krolik for stimulating discussions
about accretion on to binary systems.

REFERENCES

Ardila D. R., Jonhs-Krull C., Herczeg G. J., Mathieu R. D., Quijano-Vodniza
A., 2015, ApJ, 811, 131

Artymowicz P., Lubow S. H., 1996, ApJ, 467, 77

Bary J. S., Petersen M. S., 2014, ApJ, 792, 64

Bertout C., Basri G., Bouvier J., 1988, ApJ, 330, 350

Boden A. F., Akeson R. L., Sargent A. I., Carpenter J. M., Ciardi D. R.,
Bary J. S., Skrutskie M. F., 2009, ApJ, 696, L111

Czekala I., Andrews S. M., Jensen E. L. N., Stassun K. G., Torres G., Wilner
D. J., 2015, ApJ, 806, 154

Czekala I., Andrews S. M., Torres G., Jensen E. L. N., Stassun K. G., Wilner
D. J., Latham D. W., 2016, ApJ, 818, 156

Duchêne G., Kraus A., 2013, ARA&A, 51, 269

Dutrey A., Guilloteau S., Simon M., 1994, A&A, 286, 149

Dutrey A. et al., 2014, Nature, 514, 600

Friedjung M., 1985, A&A, 146, 366

Gillen E. et al., 2014, A&A, 62, A50

Hanawa T., Ochi Y., Ando K., 2010, ApJ, 708, 485

Jensen E. L. N., Mathieu R. D., 1997, AJ, 114, 301

Jensen E. L. N., Dhital S., Stassun K. G., Patience J., Herbst W., Walter
F. M., Simon M., Basri G., 2007, AJ, 134, 241

Lin D. N. C., Papaloizou J., 1979, MNRAS, 188, 191

Lubow S. H., Artymowicz P., 2000, in Mannings V., Artymowicz A. P.,
Russell S. S., eds, Protostars and Planets IV. Univ. Arizona Press,
Tucson, AZ, p. 731

Lynden-Bell D., Pringle J. E., 1974, MNRAS, 168, 603

Mayama S. et al., 2010, Science, 327, 306

McFadyen A. I., Milosavljevic M., 2008, ApJ, 672, 83

Messina S., Parihar P., Biazzo K., Lanza A. F., Distefano E., Melo C. H. F.,
Bradstreet D. H., Herbst W., 2016, MNRAS, 457, 3372

Nelson R. P., Papaloizou J. C. B., Masset F., Kley W., 2000, MNRAS,
318, 18

Paczynski B., 1977, ApJ, 216, 822

Papaloizou J., Pringle J. E., 1977, MNRAS, 181, 441

Ragusa E., Lodato G., Price D. J., 2016, MNRAS, 460, 1243

Roddier C., Roddier F., Northcott M. J., Graves J. E., Jim K., 1996, ApJ,
463, 326

Ruge J. P., Wolf S., Demidova T., Grinin V., 2015, A&A, 579, 110

Shakura N. I., Sunyaev R. A., 1973, A&A, 24, 337

Shi J.-M., Krolik J. H., 2015, ApJ, 807, 131

Shi J.-M., Krolik J. H., Lubow S. H., Hawley J. F., 2012, ApJ, 749, 118

Stone J. M., Norman M. L., 1992, ApJS, 80, 753

Terquem C., Sørensen-Clark P. M., Bouvier J., 2015, MNRAS, 454, 3472

Ziegler U., Yorke H. W., 1997, Comput. Phys. Commun., 101, 54

This paper has been typeset from a TeX/LATeX file prepared by the author.