Gravitational torque in circumbinary discs: global radial oscillations

Circumbinary discs (CBDs) arise in many astrophysical settings, including young stellar binaries and supermassive black hole binaries. Their structure is mediated by gravitational torques exerted on the disc by the central binary. The spatial distribution of the binary torque density (so-called excitation torque density) in CBDs is known to feature global large-amplitude, quasi-periodic oscillations, which are often interpreted in terms of the local resonant Lindblad torques. Here we investigate the nature of these torque oscillations using 2D, inviscid hydrodynamic simulations and theoretical calculations. We show that torque oscillations arise due to the gravitational coupling of the binary potential to the density waves launched near the inner cavity and freely propagating out in the disc. We provide analytical predictions for the radial periodicity of the torque density oscillations and verify them with simulations, showing that disc sound speed and the multiplicity of the density wave spiral arms are the key factors setting the radial structure of the oscillations. Resonant Lindblad torques play no direct role in determining the radial structure and periodicity of the torque oscillations and manifest themselves only by driving the density waves in the disc. We also find that vortices forming at the inner edge of the disc can provide a non-trivial contribution to the angular momentum transport in the CBD. Our results can be applied to understanding torque behaviour in other settings, e.g. discs in cataclysmic variables and X-ray binaries.

A characteristic feature of CBDs is the cavity that forms at the disc centre where the binary orbits (Artymowicz & Lubow 1996;Armitage & Natarajan 2002).The radius of the cavity  cav is typically found to lie in the interval (2 − 5) b , where  b is the binary semimajor axis, depending on the binary mass ratio, eccentricity, as well as the disc aspect ratio and viscosity (Pelupessy & Portegies Zwart 2013;Ragusa et al. 2020;Hirsh et al. 2020;Dittmann & Ryan 2022).The cavity arises due to the gravitational (tidal) torques exerted by the binary on the disc (Artymowicz & Lubow 1994), similar to the origin of gaps carved out by massive planets in protoplanetary discs (Papaloizou & Lin 1984).The value of  cav is set by the competition between the injection of angular momentum into the disc by the binary torque and the internal (viscous) disc stresses.
It is important to recognize that tidal binary-disc coupling oper-★ E-mail: rrr@damtp.cam.ac.uk (RRR) ates as a nonlocal, two-stage process (Lunine & Stevenson 1982;Greenberg 1983;Goldreich & Nicholson 1989;Rafikov & Petrovich 2012;Petrovich & Rafikov 2012): first, the torque due to the nonaxisymmetric component of the binary potential excites a density wave in the disc, adding angular momentum to the wave at the rate (per unit radius ) given by the excitation torque density d ex /d, defined explicitly in Section 2.3.After propagating over some distance, the wave dissipates (via linear or nonlinear processes) ultimately transferring angular momentum to the disc fluid at the rate (per unit ) given by the deposition torque density d dep /d, which is different from d ex /d.It is set the balance of the deposition torque density and the internal stresses that sets the value of  cav .However, it is clear that the excitation torque density is still a crucial part of the disc dynamics: after all, d ex /d is one of the factors determining the spatial distribution of d dep /d (see Section 3.2).For this reason, understanding the behaviour (radial profile) of d ex /d is very important for clarifying the properties of CBDs.
A number of past studies explored d ex /d as a function of  in simulations of CBDs for a variety of binary and disc parameters and physical assumptions (MacFadyen & Milosavljević 2008;Farris et al. 2011Farris et al. , 2014;;Shi et al. 2012;Noble et al. 2012Noble et al. , 2021;;Roedig et al. 2012;D'Orazio et al. 2013;Zilhão et al. 2015;Miranda et al. 2017;Dempsey et al. 2020;Tiede et al. 2022;Wang et al. 2022a,b;Mahesh et al. 2023), generally agreeing on a set of common features.First, d ex /d is quite small inside the cavity, which is easily explained by the low surface density of gas in that region.Second, starting from  ∼  cav and extending far in the disc, d ex /d features prominent quasi-periodic oscillations in , attaining both positive and negative Radial profile of the torque density d ex /d extracted from a simulation (blue) of a mass ratio  b = 10 −3 binary with a circumbinary disc, see Section 2 for the setup and details.A globally isothermal equation of state with the aspect ratio ℎ 0 = 0.1 (at  =  b , see Eq. ( 2)) is adopted here.Note the oscillatory, quasi-periodic behaviour of d ex /d.The salmon curve shows the resonant torque d 2 /d (depending on  via the Airy function) given by Eq. ( 20), with the amplitude chosen arbitrarily to match the first peak of d ex /d; see Section 3.2 for details.
values.The amplitude of these oscillations decays with the distance but they can often be traced out to  ≳ (5 − 10) b .Third, the radial integral of d ex /d converges to a positive value as  → ∞, indicating that the gravitational torque due to the binary drives an outward flow of angular momentum through the CBD.These features are illustrated in Fig. 1, where we show d ex /d (blue curve) extracted from one of our simulations (described in detail in Section 3) of a CBD around a binary with a mass ratio  b = 10 −3 .Understanding the spatial structure of d ex /d and, more specifically, the nature of its global radial oscillations, is the primary goal of our present work.Starting with the seminal work of MacFadyen & Milosavljević (2008), the standard explanation for the excitation torque density oscillations has been to associate d ex /d with the azimuthal  = 2 component of the local torque density d 2 /d excited in the vicinity of 3 : 2 mean motion resonance with the binary.The use of the  = 2 component of the torque density is naturally motivated by the ubiquitous presence of two prominent spiral arms in many CBD simulations (MacFadyen & Milosavljević 2008).The radial dependence enters the resonant torque eigenfunction via the argument of an Airy function (Meyer-Vernet & Sicardy 1987), which has oscillatory behaviour; see Fig. 1 and Eq. ( 20) for the explicit representation of d 2 /d.Because of that, it was suggestive to relate d 2 /d to the radial oscillations of d ex /d.This explanation of d ex /d oscillations has also been invoked in Farris et al. (2011) and Shi et al. (2012).
However, it has been noted already in MacFadyen & Milosavljević (2008) that the radial structure of d 2 /d provides a rather poor fit to the spatial distribution of d ex /d on global scales.Indeed, because of the Airy function dependence, the radial scale of d 2 /d oscillation decreases with  for  ≳  b , see salmon curve in Fig. 1 and Section 5.1.This behaviour is not supported by simulations, which generally show a different evolution of the radial period of d ex /d, as evident from Fig. 1.
In this study we provide an alternative explanation for d ex /d oscillations, which successfully reproduces their spatial variation.Our approach is closely related to the recent study (Cimerman & Rafikov 2023a) of the so-called torque wiggles in disc-planet interaction -small amplitude, quasi-periodic oscillations of the excitation torque density arising far from the planet.That work, focusing on the regimes of linear coupling (i.e. b → 0) and a smooth underlying disc around the planet (i.e.no gap), has shown that the global oscillations of d ex /d result from the gravitational coupling of the binary to the freely-propagating planet-driven density wave.Although in the binary case the underlying disc structure is very different due to the presence of a central cavity, and  b can be as high as unity, we show that the nature of d ex /d oscillations in CBDs is pretty much the same as in Cimerman & Rafikov (2023a).
In the course of this study we also found that for  b approaching unity the spatial structure of the surface density perturbation and the angular momentum budget in the inner disc exhibit curious anomalies, which are caused by the vortex-driven modes (Coleman et al. 2022a) excited by the vortices forming at the cavity edge.These vortices may be intimately related to lumps -non-axisymmetric overdensities that are frequently found to form at the cavity edge in long-term simulations of CBDs (Shi et al. 2012;Noble et al. 2012).
Our work is organized as follows.After describing our physical and numerical setup and methodology in Section 2, we present our results on the torque density behaviour in Section 3, including theoretical understanding of the origin of d ex /d oscillations in Sections 3.3&3.4.In Section 4 we describe the emergence of vortices in our simulations, their contribution to the angular momentum transport in the disc and their possible role in the lump formation at the cavity edge.We provide a discussion of our results in Section 5, including the role of resonant torques (Section 5.1) and comparison with the planetary case (Section 5.2), and summarize our findings in Section 6.

PROBLEM SETUP AND METHODOLOGY
To elucidate the nature of the torque density oscillations we design a particular CBD setup that makes it easy to reveal their nature using numerical simulations.
First, since it takes only a short amount of time for the d ex /d oscillations to develop in the CBD -basically, the sound crossing time across the simulation domain -we do not run most of our simulations for too long, see Section 2.2.This allows us to avoid the complications that arise in long-term CBD simulations: development of the eccentricity in the inner disc (MacFadyen & Milosavljević 2008;Farris et al. 2014), emergence of the lump at the cavity edge (Shi et al. 2012;Noble et al. 2012), etc.Second, since we are interested only in the gravitational torques in the CBD outside the cavity edge, we do not need to resolve the physics of what is happening inside the cavity particularly accurately.For that reason we use a numerical setup with an excised central region (Section 2.2).
Third, to reduce the number of physical processes that may complicate the interpretation of our results, in this study we mainly focus on inviscid CBDs (although we also describe some simulations with non-zero viscosity in Section 4.2), which did not receive much attention so far.Also, to lower the number of relevant parameters we use simple thermodynamic assumptions, as outlined in Section 2.1.

Physical setup
We consider a binary of two point-masses  1 ≥  2 , with total mass  b =  1 +  2 and introduce the mass-ratio  b =  2 / 1 ≤ 1.The binary orbit is circular and fixed in time, i.e.  b = const and binary eccentricity  b = const = 0. We work in the inertial (i.e.non-rotating) frame, with the origin of the polar coordinates (, ) coinciding with the centre of mass of the binary (note the difference with Cimerman & Rafikov (2023a) who worked in a frame centred on and co-moving with one of the components of the binary).The binary orbital frequency is , such that the coordinates of the (circular) binary components are given by R The gravitational potential of the binary is: where the individual contributions are We do not soften the potential, since the binary components never enter the simulation domain (see below).
The binary orbit lies within the cavity of a two-dimensional, coplanar circumbinary disc.We neglect self-gravity of the disc gas and do not account for the back-reaction it may have on the binary.
Our treatment of disc thermodynamics generally assumes a locally isothermal equation of state (EoS)  =  2 s ()Σ, where , Σ and  s are the (vertically integrated) pressure, surface density and sound speed of disc gas.The sound speed is set using the powerlaw temperature profile  () ∝  − .For  ≫  b , where the local Ω K () ≈ √︁   b / 3 , the disc scale-height is  =  s /Ω K and the disc aspect ratio ℎ = / can be expressed as where ℎ 0 is the nominal value of ℎ at  =  b .However, in most of our calculations (except in Section 5), we use the globally isothermal EoS with radially independent  and  s (i.e. = 0).This is done partly for simplicity, since in this case the constant sound speed  s (or, equivalently, the value of ℎ 0 ) becomes the only physical parameter characterizing disc thermodynamics.Another reason for choosing this EoS is that it preserves the angular momentum flux carried by the density waves in the disc.As was shown recently by Miranda & Rafikov (2019, 2020a), wave angular momentum conservation breaks down for other, more general thermodynamic setups, e.g. a locally isothermal EoS with radially varying  s .Avoiding this outcome somewhat simplifies the interpretation of our results.
Our choice of the EoS with  = 0 has a slight disadvantage that the disc aspect-ratio ℎ = / ∝  1/2 , see equation (2), i.e. the disc is flared.However, in our analysis we use the data only for  ≤ 8 b , for which the maximum aspect ratio is ℎ( = 8 b ) ≃ 0.28; the disc flaring beyond that point is not worrisome.
As mentioned earlier, most of our runs do not include explicit viscosity.However, we do consider viscosity characterized via a dimensionless -parameter (Shakura & Sunyaev 1973) in some illustrative simulations in Section 4.2.

Numerical setup
Our numerical setup is similar to the one used in Cimerman & Rafikov (2023a), with the important difference that now we work in the inertial polar frame centred on the system barycentre.We perform our hydrodynamic simulations using Athena++ (Stone et al. 2020).The code uses a Godunov scheme to solve the hydrodynamic equations in a conservative form fully accounting for the binary potential Φ b .We use Roe's approximate Riemann solver, second-order Runge-Kutta time-stepping and second-order spatial reconstruction.We run most of our simulations for 100 binary orbits, which is a relatively short time.We find this duration to be sufficient for the d ex /d pattern to fully develop (and for vortices to appear in high- b runs), so that we can perform a meaningful averaging of the fluid variables.Apart from lowering the numerical cost of the runs, such a short duration does not allow gas streamlines near the cavity edge to develop substantial eccentricity (which often happens in long-term CBD simulations) simplifying our analysis.In Section 5.3 we discuss how our results could be modified in longer term simulations.

Initial conditions
Our initial setup for the CBD includes a well-developed cavity, i.e. we do not track its formation.We follow the prescription of Muñoz et al. (2019) for the radial profile of a cavity in the CBD and adopt the following axisymmetric surface density distribution at  = 0, illustrated in Fig. 2: where  cav = 2.45 b ; the maximum of Σ(,  = 0) is about 0.41Σ 0 and is at  ≈ 3.26 b .A density floor of Σ fl = 10 −9 Σ 0 is used to avoid numerical instability in the cavity region.The prescription (3) for the cavity shape is based on the steady-state profiles of Σ() attained in viscous simulations of  b = 1 binaries (Miranda et al. 2017).Since here we mainly focus of the inviscid case and often on lower  b < 1, we do not expect Σ(,  = 0) to be in perfect equilibrium and Σ might evolve somewhat during our runs.However, this is not going to affect our results for d ex /d and their interpretation, see Section 5.
Since we consider inviscid discs, we initialize radial velocity of the flow as   (,  = 0) = 0.The azimuthal velocity   (,  = 0) = Ω() in the globally isothrmal setup is initially given by the radial force balance with where the factor  =  b (1 +  b ) −2 /4 accounts for the quadrupolar correction to the axisymmetric component of the binary gravity.The use of Eqs. ( 3) and (4) minimizes the evolution of the background disc properties after the start of the simulation.

Grid and boundary conditions
We use a numerical grid that extends from  =  b to  = 70 b in radius with logarithmic spacing and covers the full azimuthal angle 0 ≤  < 2 uniformly.The number of grid cells is   = 1400 and   = 2200 cells in radial and azimuthal directions, respectively.The large outer radius allows us to minimize its impact on the dynamics in the inner disc; in our analysis we use only the data for  ≤ 8 b .
The inner boundary condition (BC) is a 'diode' or 'outflow only' BC (e.g.Miranda et al. 2017): the values of the fluid variables in the ghost cells are always set to achieve zero radial gradient for Σ and   .The radial velocity in the ghost cells is set depending on its sign in the active cell: for   < 0, the same value is set in the ghost cell enabling free inflow; for   > 0 the ghost cell is set to −  < 0 to give zero net-flux of matter across the boundary and avoid mass injection into the simulation domain.
The values of the fluid variables in the outer ghost cells are set constant and equal to the initial data.However our results are not sensitive to this choice, since the outer radius is set artificially large and density waves do not reach the boundary during the duration of the most runs we consider (so wave reflection is not an issue).

Torque and angular momentum flux
The excitation torque density due to the gravity of the binary acting on the disc is This formula is used to compute d ex /d in the barycentric frame (cf.Cimerman & Rafikov 2023a) with Σ(, ) provided by our simulations.
We also define the angular momentum flux of the perturbation in the disc as where   =   − ũ  and ũ  =   Σ  /⟨Σ⟩  .In other words, we consider perturbations to the (generally time-dependent) densityweighted azimuthally-averaged   as a background.
We extract azimuthal averages of surface density and momentum, as well as d ex /d and   during run-time and save 100 snapshots of these processed data per orbit.This allows us to adequately average over potential time-variability of these quantities.

Units
We use one binary orbital period  b = 2Ω −1 b as a unit of time.Lengths are expressed in units of binary separation  b and surface density Σ is normalized by Σ 0 .
For the unit of   we use a generalization of the (integrated) onesided torque injected in the disc by a planet in linear disc-planet interaction (Goldreich & Tremaine 1980): This expression was also used in Cimerman & Rafikov (2023a), but now it features the reference surface density Σ 0 instead of Σ at the location of the secondary (as in the planetary case) -in the binary case the secondary is orbiting inside the cavity where Σ → 0. Use of   ,b provides us with a dimensionally suitable unit of   that would reduce to the planetary case in the linear limit  b → 0 and without a cavity.Finally, the torque density d ex /d is measured in units of   ,b / b .

RESULTS
Now we present the results of our globally isothermal CBD simulations with ℎ 0 = 0.1 for several values of the binary mass ratio  b : 10 −3 , 10 −2 , 0.1, 0.3, 1.The lowest value of  b corresponds to the so-called thermal mass (Goodman & Rafikov 2001), for which   = (/) 3 and below which the disc-perturber coupling is in the linear regime.This allows us to make meaningful comparison with the low-mass planetary case, see Section 5.3 for details.At the other end of spectrum, for  b approaching unity, the disc-binary interaction can be highly non-linear.
After describing the general morphology of the binary-induced perturbation in the disc (Section 3.1), we present the results for the d ex /d behaviour in Section 3.2, followed by the discussion of the nature of torque oscillations in Sections 3.3 & 3.4.

Perturbation to the disc structure
The top row of Fig. 3 shows maps of the surface density of the CBDs at  = 50 b as they are being perturbed by binaries with four different values of  b .One can see the density waves propagating through the disc from the edge of the central cavity, which become more visible as  b increases.The cavity itself maintains a roughly axisymmetric shape for  b = 10 −2 , 0.1, but starts to develop some asymmetries for  b = 0.3 at the time shown.In the equal mass case (panel (d)), the cavity is considerably elongated by what appears to be an  = 2 perturbation -this shape is different from an eccentric  = 1 shape often found in CBD simulations at later times.
The structure of the surface density perturbation Σ induced by the binary can be more easily seen in the bottom panels of Fig. 3, where we subtract the axisymmetric component of the surface density ⟨Σ⟩  from Σ(, ) and plot Σ/⟨Σ⟩  −1.This clearly reveals the nonaxisymmetric structure of Σ for which the binary is responsible.To facilitate the interpretation of these maps, in Fig. 4 we also plot as a function of  the strength of several low- Fourier harmonics 1 of Σ measured in our simulations at a fixed time  = 50 b .
For low  b ≤ 10 −2 the perturbation has a very regular shape in the form of two tightly-wrapped spiral arms propagating away from the centre.This is the classical density wave with  = 2 azimuthal periodicity that is launched in the CBD by the binary (see also Section 3.3).Indeed, the number of arms is consistent with Fig. 4b, which shows that  = 2 Fourier component Σ 2 dominates in the CBD out to  ≲ 4.5 b (Fig. 4a shows that Σ 2 is even more prominent MNRAS 000, 1-20 ( 2022) for  b = 10 −3 ).Moreover, this pattern is stationary in the frame co-rotating with the binary.Assuming a linear perturbation in the form , stationary in the co-rotating frame, the radial wavenumber   is given in the WKB limit by where Ω P is the (azimuthal) pattern frequency of the wave mode, which is equal to Ω b for the binary-driven waves.The phase of the perturbation is preserved along the characteristics given by d/d = −  ()/, i.e. along the curves where  ,0 is a constant and we used the relevant root of Eq. ( 8) for   .This expression is valid outside the Outer Lindblad Resonance (OLR) at which   = 0; the corresponding radius  OLR is such that Ω( OLR ) = Ω P /( + 1).We illustrate the WKB prediction (10) for the location of the wavecrest of the perturbation with the black dotted curve in Fig. 3e for  = 2, Ω P = Ω b and by choosing  ,0 to roughly align with the wave peak in the inner disc.One can see that it traces the shape of the spiral arm very well over large distances in the disc, confirming its identification with the density wave.
According to Fig. 4b, for  b = 10 −2 , the  = 1 component starts to dominate Σ outside 4.5 b , this can clearly be seen in Fig. 3e right at the outer edge (  ≃ 8 b ).Also, note that at small  each winding of the black dotted WKB curve encompasses two spiral arms, while in the outer disc it encompasses only one.We believe that this transition from  = 2 to  = 1 dominance occurs as a result of the spiral arm merger driven by the nonlinear effects: one of the two initial arms has a higher amplitude and propagates somewhat faster (Goodman & Rafikov 2001) than the lower-amplitude one, eventually overtaking it and merging with it -this process is better visualized in Fig. 9b.
To support this interpretation we note that in the lower  b = 10 −3 simulation such an arm merger occurs much further out in the disc.This is because density wave amplitude is smaller in lower  b systems (which can be deduced from the vertical scale of the different panels in Fig. 4) so that the nonlinear wave evolution is slower as well.
For higher  b = 0.1, the two tightly-wound spiral arms are still clearly visible in the inner disc and merge into a single arm in the outer disc, see Fig. 3f.However, some perturbations to their shape start to emerge.Note that in the near-cavity part of the disc the perturbation map reveals a three-fold symmetry, which does not appear as a strong  = 3 component in Fig. 4c.This is simply because in Fig. 3f the perturbation is normalized by ⟨Σ⟩  , which becomes very small inside the cavity.The origin of this  = 3 component will become clear in Section 4.
The trend of increasing irregularity continues at higher  b = 0.3, with the tightly-wrapped two-arm spiral arms losing their coherence across the whole domain and a one-sided perturbation emerging at  ≲ 4 b .Both effects are visible in Fig. 3g, and Fig. 4d also reveals a strong  = 1 component of Σ in this radial range.
In the equal-mass case ( b = 1) the perturbation pattern looks quite different from the low- b runs.Figure 3h shows that while Σ 2 is still the dominant component of the perturbation (see Fig. 4e), it is no longer tightly-wrapped as in e.g. the  b = 10 −2 case.What we see instead are the strong, open  = 2 spiral arms that are also clearly perturbed by a more-tightly wound pattern (sharp 'feathers' emanating from the prominent open spiral arms).
An explanation of these trends in the behaviour of Σ(, ) as a function of  b -a gradual transition from tightly-wrapped to more open spiral arms and the general increase of irregularity of Σ with growing  b -will be provided in Section 4.

Torque Density Profiles
Next we discuss our results on the torques characterizing disc-binary coupling.In the left column of Fig. 5, we show the radial profiles of d ex /d for five different values of  b that we consider, assuming the CBD setup with the initial Σ() distribution in the form (3). By examining these profiles, we can make the following observations.First, d ex /d exhibits oscillatory behaviour with  for all  b that we examined.Moreover, oscillations are clearly quasi-periodic, i.e. the radial spacing of the locations where d ex /d = 0 follows a particular pattern, at least approximately.Second, the amplitude of d ex /d oscillations, when normalized by   ,b , decreases with  b .However, in the lowest- b runs,  b = 10 −3 and 10 −2 , the (scaled) amplitude is almost the same, and there is an overall similarity of the scaled d ex /d profiles.Third, the amplitude of d ex /d oscillations decays with  either from the start (for high  b ) or beyond  ≈ 2.5 b (for  b ≲ 0.1).This decay with  is faster for higher  b .Fourth, the radial period of d ex /d oscillations is somewhat larger for higher  b closer to the cavity (even though there are fewer clear oscillation −0.001 0.000 0.001 (a) q b = 0.001 dT ex /dR 0.0000 0.0001 (c) q b = 0.01 0.0000 0.0001 0.0002 b given by Eq. ( 7).See Section 3.2 for details.
cycles in these runs).These trends in the behaviour of d ex /d with  and  b will be explained in Sections 3.3 & 3.4.
In the right column of Fig. 5 we show other important angular momentum metrics: the integrated gravitational torque  ex () = ∫  0 (d ex ( ′ )/d) d ′ (blue dashed curve), and angular momentum flux   (orange curve) defined by Eq. ( 6).The integrated torque  ex () is a measure of the angular momentum injected (up to radius ) by the non-axisymmetric component of the binary gravity into the density waves propagating through the disc.Angular momentum flux   is the amount of angular momentum that these density waves carry across radius  per unit of time.Thus, the two metrics are related (see below).
One can see that  ex () saturates at a finite and positive value for all  b .This implies that gravitational torque due to the binary endows density waves with a positive angular momentum, as expected (Meyer-Vernet & Sicardy 1987).At the same time,   decays as  increases.This happens since the density waves not only accumulate angular momentum due to the binary torque but also lose it as they dissipate due to the their non-linearity (Goodman & Rafikov 2001;Rafikov 2002a).As a result, while initially   tends to replicate the oscillatory features of  ex () (especially at low  b ), eventually the two curves diverge from each other and   falls below  ex ().This divergence is faster for higher  b because the density waves become more nonlinear as  b increases, as we noted in Section 3.1.
If the binary torque were the only driver of the perturbation in the disc, the difference between  ex () and   would constitute the deposition torque density d dep /d due to the binary (defined in Section 1), which characterizes the amount of angular momentum passed by the density wave from the binary to the disc fluid (Rafikov 2002b;Dong et al. 2011;Cimerman & Rafikov 2021;Dempsey et al. 2020).In the CBD context, d dep /d should be positive as the gravitational torque of the binary drives positive angular momentum flux through the disc.However, this is not the case in Fig. 5j for  b = 1, which shows that   () exceeds  ex () over a radial range (2 − 3.3) b ; in fact, a slight increase of   () above  ex () near the first peak can be seen also for  b = 0.3 in panel (i) of that figure.Interestingly, for these values of  b one also finds   not decaying to zero rapidly, as for lower  b ; instead, it falls off with  rather slowly.These intriguing anomalies will be explained in Section 4.

Origin of the Torque Density Oscillations
Figure 6 provides a graphic explanation for the origin of the torque oscillations, using the results of the  b = 10 −3 simulation as a basis.In panel (a) we plot the relative non-axisymmetric perturbation amplitude Σ/⟨Σ⟩  −1, similar to Fig. 3e, but now in Cartesian (, ) coordinates.The black dotted curve again shows the linear WKB prediction ( 9)-(10) for shape of the wake, which matches the shape of the two spiral arms very well.We have chosen a low value of  b to illustrate our point simply because we want to minimize the effects of wave nonlinearity, which would otherwise complicate our interpretation.
In panel (b) we reproduce d ex /d profile (as in Figs. 1, 5a) on the same radial scale as in panel (a).Vertical violet and green dashed lines running through both panels mark the radial locations where the peak of each spiral arm in panel (a) crosses the line2  =  2 running from  1 to  2 .One can see that these lines cross d ex /d curve near its every second node, matching the radial periodicity of d ex /d remarkably well.Thus, this figure 3 clearly shows that the periodicity of d ex /d in  is set entirely by the shape of the two binary-driven spiral arms as they wrap around the origin.In fact, this connection can be gleaned already from MacFadyen & Milosavljević (2008) by noting that the locations of d ex /d peaks in their Fig. 4 exhibit the same radial periodicity as the locations, at which the two-armed spiral wake crosses the line  = 0 in their Fig. 5.  pattern in the form Σ ∝ Σ(,  = 0) exp{i [ −  2 () − φ ()]}, i.e. a single -th azimuthal harmonic of constant amplitude relative to Σ(,  = 0) given by Eq. ( 3), which forms spiral arm(s) with the shape given by the linear WKB equation (10) in the frame co-rotating with the binary, very similar to the arms in Fig. 6a.This calculation is done for  = 1, 2 and several values of  b , with the true d ex /d from a  b = 10 −3 simulation plotted for illustration as a black dashed line One can see that  = 2 calculation in panel (b) reproduces the periodicity of the true d ex /d extremely well, independent of  bthe value of  b affects only the amplitude of torque oscillations (in this exercise the amplitude of the imposed perturbation is independent of  b ).At the same time, for  = 1 in panel (a) the periodicity has twice the radial length scale of the d ex /d from simulations.This is consistent with the idea that it is the two-armed spiral in our simulation with the shape closely following Eq.( 10) that sets the periodicity of d ex /d.This toy model also makes it obvious that the radial periodicity of d ex /d is directly related to the azimuthal wavenumber of the dominant Fourier harmonic of the density wave perturbation.These findings will be explained further in Section 3.4 In conclusion we also point out that in Fig. 6, the green markers, associated with the stronger spiral arm, steadily move closer to the violet ones (for the weaker arm) as  increases.This is yet another illustration of the nonlinear wave propagation (the front of the stronger shock travels faster catching up with the weaker one) which will ultimately result in the merger of these arms into one.

Theoretical analysis
Next, we provide a more quantitative analysis of the d ex /d oscillations in CBDs.In Appendix A, we demonstrate that the definition (5) for d ex /d can be re-cast in the following form: where  2 () is the azimuthal angle towards  2 ,   =   /,  = 1, 2, and () 1/2 () are the Laplace coefficients defined by Eq. (A3).The functions   () absorb all information about the spatial structure of Σ in the disc.Since for all  ≥ 1 the axisymmetric component of Σ(, ) does not contribute to   , we can replace Σ(, ) in Eq. ( 12) with Σ(, ) -the perturbation of the surface density in the CBD.
Note that according to Eq. ( 12), the time-averaged torque density ⟨d ex /d⟩  (which is what we display in Fig. 5) is determined only by the component of Σ, which is stationary in the frame co-rotating with the binary.In other words, only Σ(, ) in the form Σ(, − 2 ()) can give rise to non-zero   , while all other components of Σ(, ) varying in the co-rotating frame would result in no net ⟨d ex /d⟩  , given a long enough averaging period.This point will be discussed further in Section 4 (see also Fig. 10).
Motivated by this observation, we now consider the following WKB-inspired azimuthal Fourier decomposition of Σ(, ), which is stationary in the frame co-rotating with the binary: where φ () is given by Eq. ( 10), and |Σ  ()| is the amplitude of the -th harmonic propagating according to its WKB dispersion relation.Note that the Fourier components shown in Fig. 4 are based on instantaneous Σ and are somewhat different (especially for high  b ) from Σ  computed from Σ time-averaged in the binary frame, see Section 4.
As illustrated by Fig. 4, in CBDs the perturbation Σ is typically dominated by one or two azimuthal harmonics Σ  .The most common situation is the dominance of  = 2 pattern (e.g. for  b = 10 −3 , 10 −2 or 1; see also Fig. 10), although  = 1 or 3 harmonics may also play a role, at least over a range of radii (e.g. for 17  b = 10 −2 , 0.1).Assuming that Σ is indeed dominated by a single -th harmonic, such that d ex /d ≈ d ex,m /d, we can provide a mathematical justification for several observations made in Section 3.2.
First, the oscillatory behaviour of d ex,m /d with  immediately follows from Eq. ( 16) through the factor sin  φ () .This dependence also naturally explains the quasi-periodicity of the radial oscillations of d ex /d.Indeed, it predicts the nulls of d ex,m /d to occur at   ,  = 1, .. given by the condition φ (  ) =  /.In particular, Eq. ( 10) with Ω P = Ω b gives the following WKB relation between the consecutive nodes   and  +1 of d ex /d Outside the cavity, where Ω ≪ Ω b , the distance between the consecutive nodes can be approximated in the linear WKB limit as explaining the regular nature of the d ex /d periodicity with  (this lengthscale is constant for  s = const.).The direct dependence of the radial length scale of the torque oscillations on  has been previously illustrated by a toy model in Fig. 7.
In Fig. 8 we illustrate the relation ( 17), assuming that the binary drives a two-armed ( = 2) spiral pattern, by plotting it as a continuous function of   (and approximating  +1 −   ≈ ( +1 −   )/2) in the globally isothermal case with ℎ 0 = 0.1; the horizontal dotted line shows the asymptotic limit (18).We also plot (as diamonds) the data for ( +2 −   )/2 taken from our fiducial simulations with different values of  b .One can see that the WKB prediction (17) provides a good fit to the clustering of the simulation data points, supporting our interpretation of the torque oscillations put forward in Section 3.3.In fact, given the diverse morphologies of Σ in panels (e)-(h) of Fig. 3, it may seem surprising that the agreement between the relation (17) for  = 2 and simulation data for very different  b is so good; this agreement will be explained later in Section 4. The scatter of data points around the relation ( 17) is caused by the fact that other Fourier components of Σ (i.e. ≠ 2) are also contributing to d ex /d at some level in the simulations (also, for higher  b we get only a few clearly pronounced oscillations in the inner disc).
Second, the decrease of the amplitude of d ex /d normalized by   ,b with increasing  b is due to the nonlinear effects.In the linear limit and, very importantly, in a disc with a fixed surface density, the normalization by   ,b makes the disc-perturber coupling independent of the perturber mass, see Goldreich & Tremaine (1980) & Cimerman & Rafikov (2021).This is roughly the case for the lowest4  b = 10 −3 , 10 −2 , which look strongly alike, see Fig. 5a,b.However, for higher  b the cavity size starts to change ( cav increases with  b , see Fig. 2), reducing the gas surface density near the important resonances where the density wave is launched.As a result, the relative amplitude of the torque oscillations (normalized by   ,b ) decreases as  b increases.
Third, according to Eq. ( 16) the evolution of the amplitude of d ex /d oscillations with  is determined by the behaviour of both the mode amplitude |Σ  ()| and the scaling factor   (), which characterizes the weakening of the binary potential with the distance.Figure 4 shows that |Σ  ()| initially increases with  due to the initial wave excitation but then decays as a result of non-linear dissipation, which in the absence of explicit linear dissipation (e.g.Miranda & Rafikov 2020a) is the only wave damping mechanism available (Goodman & Rafikov 2001).This figure also shows that |Σ  ()| decay is faster for higher  b as a result of stronger wave nonlinearity.For that reason the amplitude of d ex /d oscillations decreases over shorter radial interval in high  b cases.But even if |Σ  ()| were to stay unchanged, the amplitude of d ex /d oscillations would still go down as  → ∞ due to the   () behaviour alone, since  () 1/2 () ∝   as  → 0. Fourth, wave nonlinearity also explains why the radial period of d ex /d oscillations near the cavity edge tends to slightly increase with  b .Figure 8 provides some evidence for this, as the higher- b points tend to lay above the linear relation (17) for  ≲ 4 b , where the wave nonlinearity is strongest.Indeed, it is known from the studies of disc-planet interaction (Rafikov 2002a;Zhu et al. 2015;Cimerman & Rafikov 2021) that the nonlinear spiral density waves launched by massive planets have higher opening angles than the wave launched by the low-mass planets.For the same reason, in our case the more nonlinear density waves in higher  b runs travel further in  per each azimuthal wrapping than in the low- b cases, explaining the (weak) evolution of the radial period of torque oscillations with  b .Note that this effect is not captured in equations ( 17)-( 18) which are based on the linear approximation ( 9)-( 10) for the global shape of the density wave.

VORTEX-DRIVEN WAVES
In Sections 3.1 and 3.2 we pointed out several interesting anomalies in CBD behaviour as  b approaches unity: the emergence of very open spiral arms and increased irregularity of the perturbation pattern in Fig. 3h, as well as the unexpected excess of   () over  ex () in a finite radial range in Fig. 5h,j.To further investigate these for  b = 0.1, 0.3, 1, respectively.Horizontal black dotted lines mark (in both rows) the outer Lindblad resonance radii for these (, Ω P ).The dark green dashed lines mark the corotation radii of these modes (i.e.radii where ⟨Ω() ⟩  = Ω P ), and they fit well the radial locations of the cores of the vortices clearly visible in the panels (h)-(j).Going from left to right, this plot shows the transition from the binary-driven spirals dominating to the vortex-driven modes dominating.
Figure 9 illustrates our findings for different values of  b at  = 50 b .The lower panels of this figure show the non-axisymmetric component of the vortensity  − ⟨⟩  in the inner part of the CBD.The use of the non-axisymmetric component allows us to remove the dominant axisymmetric vortensity contribution and focus on any structures that may be forming near the cavity edge.The top panels of Fig. 9 show contemporaneous snapshots of the relative surface density perturbation (normalized by  b ) in a Cartesian (, ) projection (i.e. as in Fig. 6a).The use of this projection helps us visualize the pitch angle  of the spiral features (the angle between the tangent to the spiral and the local azimuthal direction) in the CBD, which is given by such that the slope of the wave crests in a Cartesian map is  tan ; in the last approximation in (19) we used the WKB expression ( 8 Vortensity maps in panels (h)-(j) reveal the reason for these changes in the appearance of Σ.For  b = 0.1 one can see three highly localized (in  and ) vortensity concentrations at  ≈ 2.2 b , at the 'wall' of the cavity, see Fig. 2.These concentrations represent fluid vortices, roughly equally spaced in  and forming  = 3 vortensity structure, explaining the three-fold symmetry of Σ perturbation inside the cavity in Fig. 3f.For  b = 0.3 one can again see three sharp vortices, although two of them are in the process of merging with each other at this time (at  −  b ≈ −0.6).They are also located at larger radii, around 3 b , near the peak of Σ in Fig. 2. Finally, for  b = 1 one can see two big vortical 'rolls' around  ≈ 3 b , which are quite extended in the azimuthal direction and form an  = 2 pattern.
These different vortical structures (we will generally refer to them as just 'vortices') are the reason behind the emergence of additional open spiral arms in panels (c)-(e).Indeed, localized concentrations of vortensity are known to launch density waves in discs (Li et al. 2001;Johnson & Gammie 2005;Mamatsashvili & Chagelishvili 2007;Heinemann & Papaloizou 2009), similar to embedded planets, and this is what happens in the present case.The vortices visible in panels (h)-(j) drive spiral density waves in the CBD -the vortexdriven modes discussed in Coleman et al. (2022a) in the context of astrophysical boundary layers.These vortices are passively advected with the disc flow and orbit at the angular frequency Ω vrt corresponding to the radii  vrt at which their peaks are located.Since  vrt >  b , one naturally has Ω vrt < Ω b .According to Eq. ( 19), this makes the vortex-driven spiral arms (with Ω P = Ω vrt ) more open, i.e. have higher , than the tightly-wrapped binary-driven spirals (for which Ω P = Ω b ).
To confirm that the open spiral arms appearing for  b ≥ 0.1 are indeed the vortex-driven density waves, in panels (c)-(e) we draw (violet dotted curves) the WKB prediction for the shape of the vortex-driven modes.We do this using Eqs.( 9)-( 10) and Ω P = Ω vrt ,  determined based on vortex characteristics from panels (h)-(j), see figure caption.We show the corotation radii of the modes, coinciding with  vrt , with dashed dark green lines and the corresponding outer Lindblad resonance locations  OLR with horizontal black dotted lines.One can see that the shape of the open vortex-driven spirals is matched by the WKB density wave predictions very well, solidifying their interpretation as vortex-driven modes.
We can now address the origin of the anomalies in the behaviour of   and  ex visible in Fig. 5j, namely the excess of   () over  ex () over some radial ranges.If the binary were the only actor responsible for the wave activity in the CBD, one would always expect   () <  ex (), which is the case in other panels in that figure.However, strong vortices present near the cavity edge change this picture since the density waves that they launch also transport angular momentum through the disc (Paardekooper et al. 2010;Coleman et al. 2022b).Since density wave launching by vortices does not involve gravity (Mamatsashvili & Chagelishvili 2007;Heinemann & Papaloizou 2009), there is no direct contribution from these modes to the gravitational d ex /d.However, the vortex-driven modes do carry positive angular momentum flux outside  vrt , providing an additional positive contribution to the   carried by the (tightlywrapped) binary-driven spirals, such that the combined   ends up exceeding  ex .Coleman et al. (2022b) also found that the vortex-driven modes excited in the boundary layers of accretion discs decay with the distance due to non-linear dissipation rather weakly, leading to a slow decay of their   (see figs. 2, 3 of that paper).We observe a similar phenomenon in panels (i),(j) of the Fig. 5, in which   decays rather slowly for  ≳ 3 b -for these two  b the vortex-driven modes dominate in the outer disc.
One may wonder why the strong, open vortex-driven spiral arms, asymmetric in the binary frame, do not provide a distinct contribution to d ex /d in Fig. 5 for  b ≥ 0.3.Indeed, Eqs. ( 15)-( 16) would seem to predict such a contribution (scaling linearly with  b , unlike the regular d ex /d) with φ () determined by the equation ( 10) with Ω P = Ω vrt .The openness of the vortex-driven spirals would then naturally translate into a longer radial wave-length of the corresponding d ex /d oscillations.
However, as mentioned earlier in Section 3.4, one has to remember that the long-term (time-averaged) behaviour of d ex /d (which is what is shown in Fig. 5) is determined only by the component of Σ that is stationary in the frame co-rotating with the binary.At the same time, since Ω vrt ≠ Ω b , the vortices, as well as the spiral pattern that they excite, are not stationary in that frame.As a result, the net contribution of the vortex-driven modes to Σ in the binary frame should average out to zero over long time intervals.We illustrate this expectation in Fig. 10, where we plot both the snapshot of Σ at  = 50 b as well as the time average (over 10  b and centred on that ) of Σ in the frame co-rotating with the binary (it is also instructive to compare with Figs.3h & 9e,j).This calculation is done for  b = 1, for which the vortices are strongest (see Fig. 9) and clearly illustrates the dramatic effect of time-averaging: it very effectively filters out the vortex-driven perturbation at radii where the local orbital period is shorter than the averaging timescale (one can still see some traces of the vortex-driven modes in panel (b) outside ≈ 5 b , where this condition is not fulfilled).As a result, the time-averaged Σ retains only the structure due to the binary-driven  = 2 modes, which are tightly wrapped.We also verified that for  b = 10 −3 time averaging does not change Σ pattern at  = 50 b as vortices play minimal role for such low  b at that time.
In the end of this exercise, when we compute d ex /d and take its mean over time longer than the local orbital period, the contribution of the vortex-driven modes to the gravitational torque averages out to zero.Thus, the left panels of the Fig. 5 show us ⟨d ex /d⟩ predominantly due to the tightly-wrapped  = 2 binary-driven spirals, with no contribution from the vortex-driven modes in the inner disc, where the averaging interval is longer than the local orbital period.This naturally explains good agreement that we observe in Figure 8 between the simulation data for very different  b and the WKB prediction (17).

Origin of the vortices
While vortices are not the main subject of this study, we still provide some discussion of their possible origin.
In our simulations CBD starts with a sharp vortensity gradient in its inner region simply because of the sharp variation of Σ(,  = 0) at the edge of the cavity, see Eq. ( 3).It is known (Ono et al. 2016(Ono et al. , 2018) ) that a radially localized and large-amplitude vortensity feature can be a natural driver of the Rossby Wave Instability (RWI, see Lovelace et al. 1999).As RWI eventually results in non-axisymmetric vorticity concentrations, it represents a natural candidate for producing nearcavity vortices in our globally isothermal setting.
However, RWI driven by the initial vortensity structure in the disc is not the whole story: if it were the case, vortices would be produced in the same way for all values of  b in Fig. 9 since the initial Σ(,  = 0) is the same.Instead, there is a clear dependence on  b , with vortices appearing earlier in higher- b systems in Fig. 9. Thus, binary mass ratio is an important parameter determining the conditions for the RWI (or some other instability) to be triggered.
The most likely reason for this  b -dependence is the nonlinear dissipation of the density waves launched by the binary in the disc: the higher  b is, the more nonlinear are the density waves (see Fig. 4), leading to faster production of vortensity at the wave shocks near the cavity and ultimately causing the RWI to be triggered earlier.This picture is similar to the excitation of vortices by the RWI at the edges of planet-induced gaps in protoplanetary discs (Koller et al. 2003;de Val-Borro et al. 2007;Lin & Papaloizou 2010), with the more massive planets driving stronger density waves, producing vortensity via shock dissipation faster and causing earlier triggering of the RWI and vortex formation (Cimerman & Rafikov 2023b).
At the same time, the initial vortensity gradient present near the cavity edge likely also plays an important role, shortening the time for vortices to appear.In that regard, the CBD situation bears similarity also with the production of vortices near the boundary layers of accretion discs as described in Coleman et al. (2022a), in which a localised vortensity structure is present from the start because of the sharp gradient of the angular velocity in the boundary layer, see Fig. 8 of that work.In the planetary gap context such a vortensity feature is initially absent and arises only due to the action of the planet (Cimerman & Rafikov 2023b;Rafikov & Cimerman 2023).
To summarize, the vortices in CBDs are likely produced by the RWI triggered by the evolving vortensity structure near the cavity, which gets established during the initial cavity clearing phase and then via a steady vortensity buildup due to the dissipation of the binary-driven density waves.Importantly, regardless of the origin of the vortices, we showed that the vortex-driven modes observed in our present work are identical to the vortex-driven modes excited in the boundary layers of accretion discs (Coleman et al. 2022a) and at the edges of the planet-induced gaps (Cimerman & Rafikov 2023b).

Vortices and lumps in CBDs
On timescales much longer than  b the CBD is often observed to develop a non-axisymmetric, eccentric cavity (MacFadyen & Milosavljević 2008).Such a cavity slowly precesses in the binary frame and typically features an overdensity -a lump -near its apocentre (Shi et al. 2012).This lump has an effect on the periodicity of accretion and is often associated with the strong  = 1 spiral in the CBD.
Recently, using 2D inviscid general-relativistic hydrodynamical simulations Mignon-Risse et al. ( 2023) have followed the emergence of a lump at the inner edge of the disc.They ascribed its origin to the RWI that naturally develops at the inner edge of the disc and produces large-scale vortices that eventually evolve into a lump.Formation of a large vortex was suggested in that work as the root cause of the perturbation symmetry transitioning from  = 2 to  = 1 in circular equal mass binaries, followed by the formation of an eccentric cavity.Rabago et al. (2023) also found vortices forming at the inner edge of a polar disc around an eccentric binary in their low-viscosity runs.
Our results on the emergence of vortices and their role in the angular momentum transport support these findings and demonstrate that neither general relativity nor high disc inclination are necessary for the vortices to appear.At the same time, our findings do not provide a clear evidence in favor or against vortices being a precursor to and the cause for the lump formation.As an alternative, Shi et al. (2012) and Noble et al. (2012) have ascribed the lump origin to the interaction of gas streams inside the cavity with the cavity edge.
One way to discriminate between the different possibilities might be to examine the sensitivity of vorticity structures in CBDs to the disc viscosity, since viscous stresses are known to suppress the RWI while lumps have been observed even in high-viscosity simulations.We provide a quick look into this issue in Fig. 11, where we show snapshots of Σ and the non-axisymmetric component of vortensity for different values of viscosity (parameterized via an effective  ∈ {10 −4 , 10 −3 , 10 −2 , 10 −1 }, see Shakura & Sunyaev 1973) in a circular equal-mass binary case.The snapshots are taken at  = 50 b , same as in Fig. 9 (right panels of Fig. 9 are identical to the left panels in Fig. 11).One can see that vortices readily form for  ≤ 10 −3 , and the disc and vortensity structure are almost the same in these runs -two strong rolls around 2.7 b launching open spiral arms into the disc, see three left panels.In higher viscosity ( ≥ 10 −2 ) runs  structure changes and while there is still a clear  = 2 periodicity in the vortensity maps, the characteristic roll structures are no longer there.The two spiral arms are still present in Σ maps, albeit less open than for lower  (most likely because high viscosity reduces the size of the cavity and moves the vortensity structures closer to the binary, increasing their Ω P and decreasing their pitch angle, see Eq. ( 19)).
In Fig. 12 we illustrate the appearance of viscous CBDs at later time,  = 300 b (we only show  ≥ 10 −4 runs).One can see that in low-viscosity discs with  ≤ 10 −3 a single dominant vortex forms by that time at  ∼ (3 − 4) b launching a very open and wide spiral arm in the disc.But in the higher  ≥ 10 −2 runs the picture hardly changes compared to  = 50 b shown in Fig. 11, with  = 2 azimuthal periodicity still clearly dominating.This difference is also clearly reflected in the shape of the central cavity as illustrated in the bottom row of Fig. 12: for low  = 10 −4 , 10 −3 the cavity is eccentric with a well-defined lump at its edge, similar to the 'asymmetric' state of D 'Orazio et al. (2013), whereas for higher  = 10 −2 , 0.1 the cavity shape conforms with the  = 2 perturbation (ellipse centred on the binary barycentre), similar to the 'point-symmetric' state of D 'Orazio et al. (2013).Thus,  = 2 symmetry breaking is clearly present in our low- discs, which exhibit vortices early on and feature a strong vortex at  = 300 b , but not in high- CBDs, at least during the run time of our viscous simulations (which is still quite short compared to many other CBD simulations evolved for thousands of orbits, e.g.Miranda et al. 2017;Siwek et al. 2023).
In conclusion of this section we also note that in all our viscous runs one can still easily see the interference of the open vortex-driven waves with the tightly-wound binary-driven waves, manifested in the form of ripples on top of the vortex-driven spiral arms.Very importantly, the oscillations of d ex /d in viscous discs remain consistent with the explanations provided in Sections 3.3, 3.4.

DISCUSSION
Our results on d ex /d oscillations and their interpretation imply that the main agent determining the radial structure and amplitude of the oscillations is the behaviour of the density waves launched by the binary into the CBD with the pattern speed Ω P = Ω b .These waves get excited close to the edge of the cavity and then freely propagate through the disc, losing their strength as they dissipate due to non-linear (or linear) processes.This implies that the torque oscillations owe their existence to coupling of the binary potential to the freely-propagating density waves in the disc.
The most obvious driver of the density waves with Ω P = Ω b in the CBD is the resonant wave excitation by the binary potential at the low order Lindblad resonances (Goldreich & Tremaine 1979, 1980;Artymowicz & Lubow 1994).But in the case of a CBD some contribution to the wave excitation may also come from the ballistic motion of the gas inside the cavity, with fluid splashback against the cavity wall exciting nontrivial response in the CBD (Farris et al. 2014).However, in our picture the actual mechanism driving the density waves in the CBD does not really matter -once the waves exist the torque density will exhibit oscillatory structure.
Note that, as mentioned in Section 4, vortices forming near the cavity edge do not contribute to the gravitational d ex /d oscillations in the time-averaged sense since their Ω P ≠ Ω b and their density waves are launched non-gravitationally.While vortices do complicate our interpretation of d ex /d behaviour, a careful analysis allows us to isolate their effect on the angular momentum transport in the CBD, see Section 4.
The primary goal of this study was to understand the spatial structure of the d ex /d oscillations, namely their radial quasiperiodicity.Equations ( 17), (18) show that this periodicity depends almost exclusively on the radial behaviour of the sound speed  s , at least in the linear regime.While all the simulations that we presented so far assume globally isothermal EoS and ℎ 0 = 0.1, this remains true for other thermodynamic assumptions as well.To show this, in Fig. 13 we show d ex /d oscillations for  b = 10 −3 but different thermodynamic setups: a lower ℎ 0 = 0.05, panel (a), and a different (locally isothermal with  = 1) EoS in panel (b).One can see that in both cases the radial periodicity of d ex /d behaves in agreement with Eqs. ( 17), (18): the (constant) radial period goes down by a factor of 2 compared to the fiducial case in panel (a), while in panel (b) the radial period decreases with .This radial period evolution is in full agreement with our theory as we additionally illustrate in Fig. 14.There we plot (similar to Fig. 8) the radial spacing between the d ex /d nodes (as a function of ) from simulations with different thermodynamic assumptions and  b = 10 −3 .One can see that in each case the simulation data points (diamonds) agree very well with the analytical WKB prediction (17), asymptotically converging to the behaviour given by Eq. ( 18).Note also that the insensitivity of our theoretical calculations in Section 3.4 to the behaviour of the background surface density ⟨Σ⟩  makes our simulations and theoretical results immune to any possible inconsistencies in setting up the initial disc profile in Section 2.2.1, e.g. for low  b .
We have also linked the evolution of the amplitude of the d ex /d oscillations to that of the individual Fourier harmonics of the CBD perturbation Σ  , see Section 3.4 and Eq. ( 16).However, understanding the radial behaviour of Σ  is beyond the scope of this study.Some information on the evolution of the wave amplitude due to the nonlinear effects can be gleaned from Fig. 4, but the full quantitative understanding of these results, including the prediction of the overall density wave amplitude as a function of  b and other parameters, will require deeper understanding of the wave excitation and ⟨Σ⟩  structure near the cavity (going beyond the simple fit (3)) in the first place.Wave amplitude evolution will also be affected by the thermodynamic assumptions about the equation of state (Miranda & Rafikov 2019) and explicit cooling mechanisms (Miranda & Rafikov 2020a,b) in the CBD.
In this work torque density oscillations have been studied in the context of a CBD with a large cavity in which the masses perturbing the disc orbit.A similar phenomenon of the torque wiggles has been explored by Cimerman & Rafikov (2023a) in smooth (i.e.not featuring gaps) protoplanetary discs in the context of disc-planet interaction, see Section 5.2.But other astrophysical discs may also feature d ex /d oscillations of similar nature.
For example, in their study of discs in cataclysmic variables (CVs) orbiting one component of the binary Ju et al. (2016) reported oscillations of d ex /d (see their Fig.10) very similar to what we and others find in CBDs; even earlier such oscillations have been reported by Savonĳe et al. (1994).Based on our current understanding, we believe that the torque density oscillations in CV discs are caused by the coupling of the gravitational potential of the secondary to the freely (inward) propagating spiral density waves launched by the sec- As  s decreases with  for the orange curve, the radial period of d ex /d oscillations also decreases, consistent with the equation (18).8) for different thermodynamic assumptions: globally isothermal EoS with ℎ 0 = 0.1 (blue) and ℎ 0 = 0.05 (orange) and locally isothermal EoS with  = 1 and ℎ 0 = 0.1 (green).Coloured points show the corresponding simulation data for the distance between d ex /d nodes, dashed and dotted curves illustrate WKB equations ( 17) and ( 18), respectively, while the dotdashed curves show the resonant torque prediction (21).
ondary at the outer edge of the disc; these waves are clearly visible in Figs. 5 and 7 of Ju et al. (2016).We leave detailed exploration of torque oscillations in CV discs to a future study.

Resonant torques and d𝑇 ex /d𝑅 oscillations
Until this work, the standard explanation for the d ex /d oscillations relied on the  = 2 eigenfunction of the resonant torque density excited at the corresponding outer Lindblad resonance where Ω 2 = (2/3)Ω b (MacFadyen & Milosavljević 2008;Farris et al. 2011;Shi et al. 2012).The general idea ascends to Meyer-Vernet & Sicardy (1987) who showed that an isolated  = 2 Lindblad resonance gives rise to a torque density in the form where  2 is the OLR radius for  = 2,  2 = 2 −2/3 (/) 2  b and (/) 2 is the disc aspect ratio evaluated in the vicinity of  2 .The dimensional, radially-independent coefficient5  2 is proportional to Σ( 2 ), which makes the determination of its exact value rather ambiguous because of the rapid variation of Σ near the resonance.For that reason, in Fig. 1 we simply choose  2 so that d 2 /d approximately matches the first peak of d ex /d.
Given the asymptotic behaviour of the Airy function Ai(−) ∝  −1/4 sin (2/3) 3/2 +  0 ,  0 = const, for  ≫ 1 (Abramowitz & Stegun 1972), Eq. ( 20) predicts the distance between the consecutive nodes of d 2 /d to obey This expression is clearly different from our prediction (18) since the dependence of  s () in the latter can be arbitrary.In particular, globally isothermal simulations presented in this work have  s () = const and the scaling of  +1 −   with   ≫  b in ( 18) is different from that in (21), as illustrated in Fig. 1.This is also very obvious in Figs. 8 and 14, where the dot-dashed curves representing Eq. ( 21) provide an inadequate fit to the simulation points for different thermodynamic assumptions, unlike our WKB predictions (17), (18).Moreover, scaling of the radial period of d ex /d with  is only a part of the story.For example, in a locally-isothermal disc with  s () ∝  −1/2 (i.e. = 1 and a radially-independent aspect ratio /) Eq. ( 17) predicts  +1 −   ∝  −1/2  , just as in Eq. ( 21).However, the two predictions are still off by a factor ≈ 2.2, see green dot-dashed and dotted curves for  = 1 in Fig. 14.Also, Eq. ( 20) implies that for  ≫  b the amplitude of d 2 /d ∝  −1/4 , whereas our result (16) for the behaviour of torque oscillations means that their amplitude should scale as  −2 |Σ 2 ()| (see Section 3.4), i.e. a very different dependence.
These observations and a poor match that d 2 /d provides for d ex /d in Figs. 1, 8 & 14 make it clear that the resonant torques cannot explain the behaviour of the torque oscillations.The reason for this naturally follows from the original work of Meyer-Vernet & Sicardy (1987), which analyzed density wave excitation only in the immediate vicinity of a Lindblad resonance.The Airy function dependence arises in the behaviour of the velocity perturbation when fluid equations are expanded near  2 and it propagates also into d 2 /d.However, by construction, Eq. ( 20) is then also valid only Both are obtained from the corresponding simulations with  b = 10 −3 .In the CBD case d ex /d has been multiplied by 30 to bring the scale of its oscillations to the level of the torque wiggles in the planetary case (see bottom panel).The top panel illustrates the tall, main peak of d ex /d near the planet, which is absent in the CBD case.Note the quasi-periodicity of the torque wiggles (orange), which has a radial length scale twice larger than that of the d ex /d oscillations in the CBD (blue).
locally, very close to  2 and is simply not applicable for  −  2 ≳  b where the conspicuous torque oscillations are found.This invalidates the idea that the mathematical structure of the resonant torque expression (20) can be the reason for the global torque oscillations.The only resonance relevant for the origin of these oscillations is the corotation of the binary and the wave pattern that it drives in the disc, since only the component of Σ stationary in the binary frame is responsible for the mean d ex /d, see Sections 3.4, 4.However, resonant torques still play a role in the genesis of d ex /d oscillations, even though indirectly: they are the primary cause of the density wave excitation in the first place.It is the subsequent coupling of the binary potential to the freely-propagating density waves that causes the oscillatory d ex /d as we described in Sections 3.3 & 3.4.The resonant torque endows density waves with the angular momentum that they carry away into the disc, and most of this angular momentum gets accumulated very close to the resonance (Meyer-Vernet & Sicardy 1987).Far from the resonance, the excitation torque density d ex /d provides an insignificant contribution to the integrated torque  ex due to its oscillatory behaviour, as can be seen in the right panels of Fig. 5.

Comparison with the planetary case
Recently Cimerman & Rafikov (2023a) described and explained the phenomenon of the so-called torque wiggles in disc-planet interaction.They showed that when a planet is embedded in a smooth disc (without a gap), the torque density profile far from the planet shows a quasi-periodic pattern of low-amplitude d ex /d oscillations which are especially coherent in the outer disc.The explanation that was provided for these torque wiggles in Cimerman & Rafikov (2023a) -coupling of the planetary potential to the freely-propagating density wave launched by the planet -is essentially the same as the explanation of d ex /d oscillations in our present work and in fact motivates it.
It is instructive to provide a side-by-side comparison of the torque oscillations and the torque wiggles, which we do in Fig. 15 by showing d ex /d profiles extracted from two simulations with  b = 10 −3 .The first one is our standard CBD simulation with a surface density profile in the form (3) featuring a central cavity (blue curves).The other one assumes a planetary setup, in which the secondary is embedded in a smooth disc, i.e. there are no radial structures around the orbit of the secondary (orange curves).Given that this  b ∼ (/) 3 is at the boundary of the linear regime, one expects gap opening to eventually occur in the latter simulation.However, we ran it for only a short period of time sufficient to measure d ex /d while the disc was not strongly perturbed.
The first obvious thing to note in Fig. 15 is the very different scale of the d ex /d in two cases -note that d ex /d in CBD has been multiplied by a factor of 30 to make it similar in scale to the torque wiggles in the outer disc.The second remarkable feature is the strong peak of d ex /d near the planetary orbit in the planetary setup, see panel (a), which is completely missing in the CBD case.Both of these differences are due to the lack of gas near the perturber in the CBD setup, which does not allow the high-order Lindblad resonances to be activated.It is well known (Goldreich & Tremaine 1980) that in the planetary setting most angular momentum is carried by the wave harmonics with azimuthal wavenumbers  ∼ /ℎ ≫ 1, which are excited close to the planet and explain the strong peak of d ex /d.Because of this peak, the relative amplitude of the torque wiggles in the planetary case is low, at the level of several per cent of the maximum |d ex /d| (Cimerman & Rafikov 2023a).On the contrary, in the CBD case, because of the central cavity, only the  = 2 outer Lindblad resonance couples to the disc sufficiently well to drive this low-order harmonic.Since the central peak is absent in this case, the torque oscillations have  (1) amplitude and are much more prominent.
Focusing next on the radial structure of d ex /d far from the planet in panel (b), one can see that outside the 'negative torque density' region (at  ∈ (1.4,1.8) b , see Dong et al. 2011;Rafikov & Petrovich 2012) planetary d ex /d shows a roughly self-similar periodic structure with two peaks and a deeper trough, with slowly decaying amplitude.At the same time, the CBD d ex /d has an almost sinusoidal shape with a modulated amplitude and the radial frequency which is twice that of the planetary d ex /d pattern.The explanation of these differences lies in the different spatial structure of the density waves in two cases: in the planetary setup there is a narrow, single-armed density wave (Rafikov 2002a;Ogilvie & Lubow 2002) which is a superposition of many azimuthal harmonics but is dominated by high  ∼ /ℎ.As shown by Cimerman & Rafikov (2023a), this results in a nontrivial shape of the torque wiggles.In the CBD case a low  = 2 harmonic dominates (see Fig. 4) since the cavity suppresses higher , ensuring a single dominant sinusoidal contribution to d ex /d as follows from Eq. ( 16).The radial frequency is twice higher in the CBD case because the density wave has two arms rather than one as in the planetary case -a similar (although not exactly the same) situation can be seen in Fig. 7 for  = 1 and  = 2. Despite these differences, the underlying cause of the quasi-periodicity of the torque wiggles and d ex /d oscillations is the same -the presence of the density wake in the disc and its coupling to the perturber's potential.Since the wake shape is set mainly by the radial profile of the sound speed  s , in both cases ⟨Σ()⟩  profile plays essentially no role in setting the radial periodicity of d ex /d.

Comparison to previous studies
As mentioned in the Introduction, a number of CBD studies have previously looked into the behaviour of d ex /d.These investigations used both 2D and 3D setups, considered standard  viscosity as well as full MHD models allowing for magnetorotational instability, and some of them included general relativistic effects.Despite these differences, all these studies confirm the large-amplitude, oscillatory nature of d ex /d in CBDs.
Unlike most other studies using explicit shear viscosity, we mainly considered the inviscid limit.This makes our work similar to Mignon-Risse et al. (2023) and Rabago et al. (2023) who also explored lowviscosity CBDs and found, similar to us, that vortices readily form at the edge of the cavity, see Section 4. Compared to many other studies (e.g.Siwek et al. 2023), our simulations are run for a rather short time.Nevertheless, our results will be helpful for understanding the outcomes of the longer-term simulations.
Working in the planetary setting, Dempsey et al. (2020) presented results on the evolution of d ex /d as the width of the gap carved by the planet increases.By examining their Fig. 5 we see a gradual transition from the radially-spaced, low-amplitude torque wiggles in the case of a narrow and shallow gap, still allowing wave excitation at high- Lindblad resonances, to the less radially-spaced, largeamplitude torque density oscillations in the case of a deep and wide gap, similar to a cavity in the CBD.This evolution can be easily understood based on our discussion in Section 5.2: as the gap depth and width increase, higher- perturbation harmonics get gradually deactivated, until  = 2 becomes the dominant one.As a result, the relative amplitude of of d ex /d oscillations increases (as the d ex /d peak near the planet is gradually eroded) and the one-armed spiral perturbation transitions to a two-armed density wave, leading to the increased radial frequency of oscillations.
In conclusion, we comment on a couple of byproducts of our work that may be useful for interpreting the results of other studies.In their investigation of CBD dynamics Mahesh et al. (2023) presented Fourier decomposition of Σ (see their Fig.8) and argued, based on the rough radial independence of the relative contribution of different harmonics, that the resonant torques are not involved in setting the cavity size in the course of disc-binary interaction.Our calculation of |Σ  | in Fig. 4 shows a very similar picture, with different harmonics maintaining roughly the same order as  varies.But our interpretation is very different: this decomposition only reflects the harmonic content of the freely-propagating density waves launched by the binary (and not the radial scaling of different multipoles of the binary potential at  ≫  b ).And resonant torques are the reason why these waves are excited in the first place, they just have little to do with the radial structure of d ex /d oscillations.Thus, we would caution against dismissing them in the CBD context.
In their study of waves in CV discs Van den Bossche et al. (2023) found the emergence of a one-armed spiral arm in the inner disc and the trend towards the dominance of  = 1 perturbation component in high Mach number (low  s ) discs.Based on our results on spiral arm merging in Figs.3e and 9b (see Sections 3.1, 3.3), we suggest that in the CV context a single arm can also arise as a result of a nonlinearly driven merger of the two spiral arms excited at the outer edge of the CV disc.This tendency would get stronger as the flow Mach number and, correspondingly, the density wave nonlinearity increase, just as observed in Van den Bossche et al. (2023).

SUMMARY
We have explored angular momentum exchange in a circumbinary disc and the properties of the excitation torque density d ex /dan important characteristic of the binary-disc gravitational coupling -focusing on its nature and radial periodicity.Using 2D, primarily inviscid, hydrodynamic simulations and analytical arguments, we found the following.
• Outside the central cavity, in the bulk of the CBD (at  ≳  b ), gravitational torque arises due to the coupling of the binary potential with the freely-propagating density waves launched near the cavity edge (similar to the origin of the torque wiggles explored by Cimerman & Rafikov (2023a) in the planetary context).
• The radial periodicity of excitation torque oscillations is set primarily by the rate at which these density waves wrap around the origin, which in turn is determined by the gas sound speed.
• The periodicity of d ex /d also depends on the symmetry properties of the waves, i.e. a number of their spiral arms.Two-armed  = 2 density waves often encountered in CBDs result in twice higher radial frequency of d ex /d oscillations than the d ex /d frequency for single-armed spirals.
• The evolution of the amplitude of d ex /d oscillations with radius is set by the density wave dissipation and the radial dependence of the binary potential.
• The local resonant Lindblad torque, often invoked to explain d ex /d oscillations, plays no direct role in their origin and periodicity.It is involved only indirectly, as the initial driver of the density waves that give rise to d ex /d oscillations.
• Angular momentum transport in the inner CBD can be affected in a non-trivial manner by the vortices that readily form (especially for high  b ) at the inner edge of the cavity in our inviscid simulations (although we observe their effect in our viscous runs as well).These vortices may be relevant for the origin of 'lumps' in CBDs.
To clarify the picture of the binary-CBD coupling, future work should aim to understand the overall amplitude and radial decay of the torque density in CBDs as a function of the binary mass ratio and disc parameters.Results of our study can also be used to understand the torque structure in other objects, e.g. in discs orbiting one of the components of the binary in cataclysmic variables, X-ray and young stellar binaries, see Section 5.
Figure1.Radial profile of the torque density d ex /d extracted from a simulation (blue) of a mass ratio  b = 10 −3 binary with a circumbinary disc, see Section 2 for the setup and details.A globally isothermal equation of state with the aspect ratio ℎ 0 = 0.1 (at  =  b , see Eq. (2)) is adopted here.Note the oscillatory, quasi-periodic behaviour of d ex /d.The salmon curve shows the resonant torque d 2 /d (depending on  via the Airy function) given by Eq. (20), with the amplitude chosen arbitrarily to match the first peak of d ex /d; see Section 3.2 for details.

Figure 2 .
Figure 2. Radial profile of the azimuthally-averaged surface density ⟨Σ⟩  for different values of  b used in this work.We show ⟨Σ⟩  time-averaged from  = 40 b to 50 b and sampled ten times per binary orbit  b , as well as the initial surface density profile (black dashed curve) given by the equation (3).

Figure 3 .
Figure 3. 2D maps at  = 50 b of the simulated (using globally isothermal EoS with ℎ 0 = 0.1) surface density (top row) and relative surface density perturbation Σ/⟨Σ⟩  − 1 (rest of figure), inversely scaled with the mass ratio  b .The central white disc marks the excised region in which the binary orbits.At this time, the binary components are aligned with the -direction.In panel (e)-(h), we do not show data for  > 8 b , in order to make it easier to follow spiral features.The black dotted curve in panel (e) shows the shape of the density wave in the WKB approximation (Eq.10) for  = 2, Ω P = Ω b ,  OLR = 1.45 b ,  2,0 = 0 and Ω() = ⟨Ω⟩  () taken from the simulation.MNRAS 000, 1-20 (2022)

Figure 4 .
Figure 4. Radial profiles of the absolute magnitude of the Fourier coefficients Σ  of the simulated density wave perturbation Σ at  = 50 b .Different panels show results for different values of  b .

Figure 5 .
Figure 5. Left: Excitation torque density d ex /d in a circumbinary disc extracted from simulations (with ℎ 0 = 0.1) for binaries with different mass ratio  b .Right: Angular momentum flux   () and cumulative (integrated) torque  ex ().All quantities are averaged over 10 orbits from  = 45 b to  = 55 b using 100 samples per orbit.Total torque and   are expressed in units of  ,b ∝  2b given by Eq. (7).See Section 3.2 for details.

Figure 6 .
Figure 6.Illustration of the origin of d ex /d oscillations using a globally isothermal simulation with ℎ 0 = 0.1 and  b = 10 −3 : (a) a map of Σ/⟨Σ⟩  − 1 illustrating the two-armed shape of the density waves in the CBD in Cartesian  −  projection, (b) d ex /d plotted on the same radial scale.Black dotted curve in panel (a) shows the linear WKB prediction (10) for the wave arm shape (with  = 2, Ω P = Ω b and some chosen  ,0), which provides an excellent fit.Vertical violet and green markers show the radii where the peaks of each spiral arm align with the binary such that  =  2 .These markers clearly show that d ex /d periodicity is directly related to the way the density waves wrap around the origin.
Figure 7 takes this point even further by showing d ex /d (rescaled by  b ) in a CBD with an artificially imposed perturbation

Figure 7 .
Figure7.Torque density (re-scaled by the binary mass ratio  b ) due to an artificially imposed surface density perturbation that assumes Σ/Σ (,  = 0) = const and has only a single azimuthal Fourier component in the form of spiral arm(s) co-rotating with the binary, with the shape given by the linear WKB relation (10): Σ ∝ Σ (,  = 0) exp{i[  −  2 ( ) − φ () ] }.The top panel shows  = 1, the bottom  = 2, for a range of  b (colours).The initial phase  ,0 in φ () is adjusted to match the location of the first peak of the d ex /d obtained from a  b = 10 −3 simulation (black dashed curve, arbitrarily rescaled for presentation purposes to match the first significant peak of  −1 b d ex /d in amplitude).The grey shaded area indicates the heavily depleted cavity region.The fact that this toy model reproduces d ex /d periodicity well for  = 2 and all  b means that it is the two-armed spiral with the pattern speed Ω P = Ω b that sets the periodicity of torque oscillations in our  b = 10 −3 simulation.

Figure 8 .
Figure 8. Radial spacing of the nodes of d ex /d -the distance between the consecutive radial locations   where d ex /d becomes zero (approximated as ( +2 −   )/2 for symmetry) -plotted as a function of   .Points of different colour correspond to our fiducial inviscid simulations (globally isothermal with ℎ 0 = 0.1) with different  b .Black dashed line shows  +1 −   computed as a continuous function of   using the analytical WKB relation (17) with  = 2, providing a good fit to the simulation data (coloured points).The dotted line reflects the asymptotic relation (18).The dot-dashed curve represents the asymptotic 'Airy' fit (21) discussed in Section 5.1.

Figure 9 .
Figure 9. Snapshots at  = 50 b of (top row) relative surface density perturbation Σ/⟨Σ⟩  , scaled by  b and (bottom row) vortensity deviation  − ⟨ ⟩  from its azimuthal average ⟨ ⟩  , plotted for different  b .Black dotted curve in panel (a) shows the linear WKB prediction (10) for the wave arm shape for  = 2 and Ω P = Ω b (with some chosen  2,0 ) providing an excellent match to the shape of the low-amplitude spiral arms driven by the  b = 10 −3 binary.In panels (c)-(e) we show WKB fits (10), (19) for vortex-driven density waves (violet dots), in which we use (, Ω P ) = (3, 0.33Ω b ), (2, 0.22Ω b ), (2, 0.23Ω b )for  b = 0.1, 0.3, 1, respectively.Horizontal black dotted lines mark (in both rows) the outer Lindblad resonance radii for these (, Ω P ).The dark green dashed lines mark the corotation radii of these modes (i.e.radii where ⟨Ω() ⟩  = Ω P ), and they fit well the radial locations of the cores of the vortices clearly visible in the panels (h)-(j).Going from left to right, this plot shows the transition from the binary-driven spirals dominating to the vortex-driven modes dominating.
) for   and assumed  ≳  OLR , Ω ≪ Ω P .The black dotted curve in panel (a) again illustrates the WKB prediction (9)-(10) for the shape of the binary-driven spiral in a  b = 10 −3 CBD assuming  = 2 and Ω P = Ω b .As was already evident in Figs.3e and 6a, it fits the shape of the spiral arms very well.The slow evolution of the two arms towards their merger into one as  increases (see Section 3.1) can also be traced in this panel.This picture hardly changes for  b = 10 −2 (panel (b)), only the merger of the two arms into one is more obvious.The vortensity maps for these values of  b (panels (f), (g)) show some non-trivial activity only inside the cavity but not in the bulk of the CBD.Things change noticeably for  b = 0.1 in panel (c) as the Σ map starts showing irregularities: while the tightly wound spiral arms are still visible in this map, their crests (still well described by the WKB prediction) are clearly perturbed by some additional quasi-regular pattern.For  b = 0.3 (panel (d)) this additional pattern becomes more pronounced in the form of a number of open spiral arms that compete in strength with the underlying original tightly wrapped spirals.Finally, in panel (e) (for  b = 1), Σ is clearly dominated by the two open spiral arms perturbed by the remnants of the original tightly-wrapped spirals, which have now become subdominant.These open spiral arms are clearly seen in Fig. 3d,h.

Figure 10 .
Figure 10.Comparison of (a) a surface density perturbation snapshot at  = 50 b and (b) time-average of the same variable in the frame co-rotating with the binary over the interval (45 − 55)  b .Data are from the globally isothermal simulation with ℎ 0 = 0.1 and  b = 1.One can see that timeaveraging in the frame co-rotating with the binary effectively eliminates the vortex-driven perturbation in the inner disc.

Figure 12 .
Figure 12.Same as Fig. 11 but at  = 300 b and only showing viscous models.We use  = 1 for the first two (a,b) and  = 2 for the second two (c,d) WKB fits (10), (19), with pattern speeds Ω P = 0.17Ω b , 0.17Ω b , 0.26Ω b , 0.34Ω b , from left to right, respectively.The bottom row of this figure illustrates the shape of the central cavity in polar coordinates in each of these runs at late times.The pink dashed circles have a radius of  c = 2 b , and provide a reference for the shape of the cavity.

Figure 13 .
Figure13.Comparison of the ex /d profiles from our fiducial globally isothermal ( = 0) runs with ℎ 0 = 0.1 (blue) and simulations with other thermodynamic assumptions (orange); all runs are for  b = 10 −3 , as in Fig.1.(a) Same EoS but ℎ 0 = 0.05, which naturally leads to doubling of the radial frequency of d ex /d oscillations, consistent with Eqs.(17), (18).(b) Same ℎ 0 = 0.1 but a different EoS -locally isothermal with  = 1, i.e.  ∝  −1 .As  s decreases with  for the orange curve, the radial period of d ex /d oscillations also decreases, consistent with the equation (18).

Figure 14 .
Figure14.Illustration of the radial periodicity of d ex /d oscillations (similar to Fig.8) for different thermodynamic assumptions: globally isothermal EoS with ℎ 0 = 0.1 (blue) and ℎ 0 = 0.05 (orange) and locally isothermal EoS with  = 1 and ℎ 0 = 0.1 (green).Coloured points show the corresponding simulation data for the distance between d ex /d nodes, dashed and dotted curves illustrate WKB equations (17) and (18), respectively, while the dotdashed curves show the resonant torque prediction (21).

Figure 15 .
Figure 15.Comparison of the excitation torque density profiles (on different scales, top vs botom) in the circumbinary setup, with central cavity (blue) and planetary setup, with no cavity or gap around the planetary orbit (orange).Both are obtained from the corresponding simulations with  b = 10 −3 .In the CBD case d ex /d has been multiplied by 30 to bring the scale of its oscillations to the level of the torque wiggles in the planetary case (see bottom panel).The top panel illustrates the tall, main peak of d ex /d near the planet, which is absent in the CBD case.Note the quasi-periodicity of the torque wiggles (orange), which has a radial length scale twice larger than that of the d ex /d oscillations in the CBD (blue).