Orbital evolution of close binary systems: comparing viscous and wind-driven circumbinary disc models

Previous work has shown that interactions between a central binary system and a circumbinary disc (CBD) can lead to the binary orbit either shrinking or expanding, depending on the properties of the disc. In this work, we perform two-dimensional hydrodynamical simulations of CBDs surrounding equal mass binary systems that are on fixed circular orbits, using the Athena++ code in Cartesian coordinates. Previous studies have focused on discs where viscosity drives angular momentum transport. The aim of this work is to examine how the evolution of a binary system changes when angular momentum is extracted from the disc by a magnetised wind. In this proof-of-concept study, we mimic the effects of a magnetic field by applying an external torque that results in a prescribed radial mass flux through the disc. For three different values of the radial mass flux, we compare how the binary system evolves when the disc is either viscous or wind-driven. In all cases considered, our simulations predict that the binary orbit should shrink faster by a factor of a few when surrounded by a wind-driven circumbinary disc compared to a corresponding viscous circumbinary disc. In-spiral timescales of $\sim 10^6$-$10^7$yr are obtained for circular binaries surrounded by CBDs with masses typical of protoplanetary discs, indicating that significant orbital shrinkage can occur through binary-disc interactions during Class I/II pre-main sequence phases.


INTRODUCTION
About half of main sequence Sun-like stars in the Solar neighbourhood are members of multiple systems (Abt 1983;Duquennoy & Mayor 1991), and their properties provide important hints about their formation and evolution.Approximately 34% are in binaries, 8% are in triples and 3% are in higher order systems, and binary orbital periods can be approximated by a log-Normal distribution with a peak at  b ∼ 10 5 days, corresponding to a separation of  b ∼ 50 au for an assumed total binary mass of 1.5 M ⊙ (Raghavan et al. 2010).The distributions of binary eccentricities,  b , and mass ratios,  b , are essentially flat except for an excess of 'stellar twins' for  b > 0.95.Moe et al. (2019) showed the binary fraction of Solar-type stars is strongly anticorrelated with stellar metallicity for periods < 10 5 days, but for longer periods the anticorrelation disappears.From this they suggest that close systems form via disc fragmentation, where opacity regulates gravitational instability, and wider systems originate from the fragmentation of optically thin turbulent molecular cloud cores.
Pre-main sequence binary systems are also common (Ghez et al. 1993;Leinert et al. 1993;Mathieu 1994).Detections of several young embedded systems on close orbits (< 30 au) have been made among Class 0/I sources (Tobin et al. 2016), with some showing evidence for circumbinary material (Tobin et al. 2022).Among Class II sources, several binary systems surrounded by Keplerian circumbinary discs (CBDs) have been detected (e.g.Mathieu et al. 1997;Dutrey et al. ★ E-mail: g.a.turpin@qmul.ac.uk 1994; Czekala et al. 2019).In the case of the short-period spectroscopic binaries AK Sco, V4046 Sgr, DQ Tau and UZ Tau, there is strong evidence that the binary orbits and disc midplanes on large scales are closely aligned with each other (Czekala et al. 2019), such that they might be viewed as precursors to transiting circumbinary planet systems similar to Kepler-16b, -34b and -35b (Doyle et al. 2011;Welsh et al. 2012).Recent analyses of close orbiting systems ( b < 10 au) using spectroscopy from the APOGEE survey indicates that the pre-main sequence binary fraction is consistent with that observed in the field (Kounkel et al. 2019), in agreement with earlier studies (Mathieu 1994).Furthermore, the period distribution of Class II/III binaries in the range 2 -10,000 days ( b = 0.05-10 au) is almost identical to the field star population (Kounkel et al. 2019), showing that many of the properties of close binary systems are established in the pre-main sequence phase.
Understanding the formation and early evolution of close binary systems remains an active area of research (see the recent review by Offner et al. 2023).Fragmentation of turbulent molecular cloud cores is expected to produce binaries with separations of ∼ few 100 au to ∼ few 1000 au (Offner et al. 2010), and fragmentation of massive protostellar discs is expected to occur on length scales ≳ 50 au (Kratter et al. 2010).Hence, the formation of close binary systems with  b < 10 au requires substantial migration to occur.Moe & Kratter (2018) considered the formation of very close binaries ( b ≤ 10 days) via Kozai-Lidov oscillations and dynamical instabilities in triple systems, combined with dissipation through stellar tides, and showed that significant additional dissipation from circumstellar/circumbinary material would be required to produce sufficient numbers of close systems to be in agreement with the observed frequency in main sequence and pre-main sequence systems.Simulations that produce clusters of stars from the fragmentation of turbulent molecular clouds successfully produce close binary systems with separations < 10 au through a combination of dynamical interactions in multiple systems, accretion, and dissipation of energy through interactions with circumbinary discs and nearby cloud material (Bate et al. 2002;Bate 2012).These simulations evolve for time scales on the order of the free fall time (∼ 10 5 yr) and are limited to resolving interactions on scales ≳ 1 au, so the ability of this scenario to produce the shortest period binaries is yet to be demonstrated.
Observations and simulations show that close pre-main sequence binary systems are expected to be surrounded by circumbinary discs, and energy and angular momentum exchange with these discs, leading to orbital shrinkage, is one means by which the closest binary systems could arise from initially more widely spaced systems (Lin & Papaloizou 1979;Artymowicz et al. 1991).However, a number of outstanding questions remain about the outcome of interactions between pre-main sequence binaries and circumbinary discs: during which stage (Class 0/I/II) of the pre-main sequence lifetime does most of the orbital evolution occur?; what role do disc-binary interactions have in shaping the distribution of binary eccentricities?; how does the evolution depend on the initial properties of the central binary and the circumbinary disc?In this paper, we start to address this latter question and examine how disc-binary interactions change when angular momentum transport in the circumbinary disc arises through the launching of a magnetised wind instead of via turbulent viscosity.In this initial proof-of-concept study, we use 2D hydrodynamical simulations combined with a prescribed external torque to mimic the effect of a magnetic torque, similar to that adopted in recent studies of disc-planet interactions (Lega et al. 2022;Nelson et al. 2023).
Circumbinary discs are thought to play an important role in the evolution of binary black holes as well as in pre-main sequence stellar systems, and may provide a solution to the 'final parsec problem' (Begelman et al. 1980;Milosavljević & Merritt 2003).Hence they have been the subject of numerous recent studies.Previous hydrodynamical simulations of CBDs have adopted a range of numerical schemes and codes: SPH (Bonnell & Bate 1994a,b;Ragusa et al. 2021); polar meshes using finite-difference ZEUS-type codes (Pierens & Nelson 2007;Mutter et al. 2017;Pierens et al. 2020;Coleman et al. 2022) and energy-conserving finite-volume codes such as PLUTO (Thun et al. 2017a;Miranda et al. 2016;Sudarshan et al. 2022;Penzlin et al. 2022); unstructured meshes using AREPO (Muñoz et al. 2019(Muñoz et al. , 2020)); Lagrangian moving meshes using the Disco code (Duffell et al. 2020;D'Orazio & Duffell 2021); Cartesian meshes using Athena++ and MARA (Shi et al. 2012;Moody et al. 2019;Tiede et al. 2020;Zrake et al. 2021).In this work we use the Athena++ code with refined Cartesian meshes.
Hydrodynamical simulations of CBDs show that a common outcome is the formation of a tidally truncated eccentric inner cavity surrounded by an eccentric, precessing circumbinary disc.Gas accretes through the CBD towards the cavity, and then forms circumstellar discs (CSDs, sometime also referred to as mini-discs) around the central stars after streamers of gas are pulled off the inner edge of the CBD by the individual stars.In general, accretion onto the stars is found to be pulsed due to the orbital phasing associated with the generation of the streamers.A number of studies have adopted relatively thick and viscous discs models (ℎ = 0.1, where ℎ ≡ / is the disc aspect ratio, and  = 0.1) due to their rapid equilibration properties (e.g.Miranda et al. 2016;Muñoz et al. 2019;Moody et al. 2019).For equal mass binaries on fixed circular orbits, it is found that interaction between the stars and the surrounding gas causes the binary orbit to expand rather than shrink, primarily because the CSDs exert a positive torque that is larger in magnitude than the negative torque exerted by the CBD.As the disc thickness is lowered to a value of ℎ ≲ 0.05, however, the binaries are found to in-spiral (Heath & Nixon 2020;Tiede et al. 2020;Dittmann & Ryan 2022;Penzlin et al. 2022).Decreasing the viscosity can also lead to in-spiralling binaries (Penzlin et al. 2022).Muñoz et al. (2020) studied circular binaries for a range of mass ratios (0.1 ≤  b ≤ 1.0) and found binaries with a mass ratio of  b ≥ 0.3 lead to expanding binaries, even when reducing the viscosity.Studies have shown the eccentricity of binaries can also play a role in the expansion or contraction of the binary orbit (Muñoz & Lai 2016;Muñoz et al. 2019Muñoz et al. , 2020;;D'Orazio & Duffell 2021;Zrake et al. 2021).
The previously described numerical work has largely been carried out in a regime where angular momentum transport in the CBD and CSDs is assumed to arise because of turbulent viscosity, parameterised through the alpha prescription (Shakura & Sunyaev 1973), and presumably originating from the magnetorotational instability (MRI) (Balbus & Hawley 1991).So far, no work has been carried out to investigate how the interaction changes if angular momentum transport and accretion are instead driven by a magnetised wind (e.g.Bai & Stone 2013), as is believed to occur in protoplanetary discs where the ionisation fraction is low.In this paper, we study and compare the two leading theories of angular momentum transport in accretion discs, viscous and wind-driven.A key point is that viscous discs evolve diffusively, whereas mass transport in a wind-driven disc is purely advective, such that we expect there to be significant differences in disc structure between the two cases even when the mass flux being transported through the disc is the same (Lega et al. 2022;Nelson et al. 2023).We measure the gravitational and accretional torques to determine whether the binary separation increases or decreases, and at what speed.
The paper is structured as follows.In Section 2 we outline the numerical methods, basic equations, and physical models.In Section 3 we discuss the basic theory required to understand how torques are calculated and how they relate to the migration rate of a binary.In Section 4 we calibrate the numerical viscosity that arises from the use of a Cartesian grid, and in Section 5 we perform a comparative study with previous work carried out by Tiede et al. (2020).We present our viscous and wind-driven disc results in Sections 6 and 7, respectively.In Section 8 we make a comparison between the results obtained from the two different types of disc models, and we discuss how the morphology of the disc effects the torque in Section 9. Finally, we summarise the work and draw conclusions in Section 10.

NUMERICAL METHOD AND PROBLEM SETUPS
In this work we present three sets of simulations.The first uses viscous ring spreading experiments to calibrate the level of numerical viscosity that arises because we simulate rotating flows on a Cartesian mesh.The second involves a suite of simulations of CBDs with different Mach numbers (or equivalently disc aspect ratios) to compare the results of our setup with those of Tiede et al. (2020) and Penzlin et al. (2022).The third involves a comparison between CBD simulations that adopt different prescriptions for the angular momentum transport mechanism acting in the disc: viscous versus wind-driven.In the following sections we describe the numerical method used to perform these different simulation suites and the initial and boundary conditions that were implemented.

Basic equations and numerical setup
All simulations presented here were performed using Athena++ (Stone et al. 2020) in 2D Cartesian coordinates.Athena++ uses a finite volume Godunov scheme to solve the following vertically integrated equations of hydrodynamics where Σ is the surface density of the disc, v is the velocity of the gas,  is the gas pressure, T vis is the viscous stress tensor, and F  is the gravitational force per unit area.The terms Σ sink and Σ sink v represent the accretion of mass and momentum by sink particles that are implemented in our CBD runs, as described below.The isothermal sound speed is defined according to where ℎ is the local aspect ratio in the disc and  is the gravitational potential.The aspect ratio is related to the Mach number, M, through the expression M = 1/ℎ.By default, Athena++ treats the gas as being adiabatic, and we implement an effective locally isothermal equation of state by imposing instantaneous cooling at each timestep, following Moody et al. (2019).
The viscous ring spreading simulations described below are performed with a single central star located at the origin of the coordinate system.When binary stars are implemented they are of equal mass and remain on fixed, circular orbits in all simulations, with a semimajor axis of  b = 1.The binary centre of mass is located at the origin of the coordinate system.The gravitational potential of the binary is given by where  is the gravitational constant, the total binary mass  b =  1 +  2 = 1,  1 and  2 are the distances to each star, and  s is the gravitational softening length, chosen to be  s = 0.05 b .Accretion of mass and momentum onto sink particles surrounding the stellar components is implemented through the Σ sink term in eqns (1) and (2), which is given by where  sink =   = 0.05 b and  −1 sink = 8Ω b .This accretion prescription is equivalent to that used by Tiede et al. (2020).
For all runs involving CBDs, the extent of the domain runs between ,  ∈ [−32, 32].We use a base resolution of   ×   = 2400 × 2400 cells, but implement a static mesh refinement region,   ,   ∈ [−6.4, 6.4] where the resolution is increased by a factor of two and which contains the binary and the cavity region of the circumbinary disc.This equates to a cell size of 0.01 3 b within the refinement region, which is similar to the resolution of 0.0156 b within the cavity used in Tiede et al. (2020).
Our boundary conditions reinforce the initial conditions at each timestep unless stated otherwise.

Setup for numerical viscosity calibration
Before carrying out inviscid wind-driven CBD simulations, it is first important to calculate the numerical viscosity of the setup.We must ensure that we are not being dominated by effects from numerical viscosity, and that the dominant driver of mass accretion through the disc is due to the wind-driven torque.A pressureless, viscously spreading ring follows the analytical solution (Lynden-Bell & Pringle 1974) where  ring is the mass of the disc,  = / 0 and  = 12  −2 0 are the normalised distance and time quantities respectively,  is the kinematic viscosity coefficient and  1 4 is the modified Bessel function of order 1/4.The surface density of the ring is initialised with this analytical solution, at time  = 0.018, and the initial velocity is Keplerian.The aspect ratio is set to ℎ = 0.005, which allows the ring to act as a pressureless fluid while maintaining numerical stability.We adopt outflow conditions at the boundaries.
We first evolve a viscous ring with a kinematic viscosity of  = 10 −5 , to ensure the spreading of the ring correctly follows the analytical solution, and this is shown in Figure 1.To calculate the numerical viscosity, we run simulations with  = 0 for different numbers of grid cells to see how the numerical viscosity varies with numerical resolution.Any spreading of the ring will now be due to the numerical effects, allowing us to calculate the numerical viscosity of the simulations by fitting the results to the analytical solution eqn.(6).Further details of the setup, and and a comparison between the simulation codes PLUTO, FARGO and Athena++, may be found in Joseph et al. (2023).The results of this viscous ring problem are presented in Section 4.

Setup for Tiede et al. (2020) comparison
In order to calibrate the code against previous simulations of CBDs, we perform simulations of varying Mach number with equivalent initial conditions to Tiede et al. (2020), allowing us to make a direct comparison to their work.
The grid setup is as described in Section 2.1, where the initial condition is a finite disc in quasi-steady state with an inner cavity that has been mildly depleted: The disc has a peak density at   = 4 b and a surface density floor of Σ floor = 10 −6 .The calculations performed by Tiede et al. (2020) used a constant kinematic viscosity coefficient of  = where Ω b = 1 is the orbital angular frequency of the binary.The Mach number is varied from 10 to 40 with an interval of 5 (equivalent to varying ℎ between ℎ = 0.1 and ℎ = 0.025), resulting in a suite of 7 simulations.We present the results from this calibration and comparison in Section 5.

Setup for wind-driven versus viscous CBDs
The grid setup for these runs has been described in Section 2.1.To mimic an inviscid disc where angular momentum transport and accretion is due to a magnetically driven wind, we apply an azimuthal acceleration on the gas, which for a disc around a single star would take the form (Lega et al. 2022;Nelson et al. 2023) Figure 1.1D surface density evolution of a viscously spreading ring with a kinematic viscosity of  = 1 × 10 −5 .The crosses are the analytical solution at viscous time , using Equation 6.
Here,  is the radial flux of gas through the disc, whose value we choose at the start of a simulation and keep constant throughout.This breaking force replicates in a simple manner the torque the disc would feel due to a centrifugally-driven magnetised wind.
To adapt this equation for binaries we define three regions of the domain, moving from the outer regions to the inner regions of the CBD-plus-binary system: an outer CBD that orbits around the centre of mass of the binary and in which the radial mass flux is directed towards the centre of mass; an intermediate region that falls inside the cavity and where the gas is treated as being ballistic (we refer to this as the ballistic zone); CSDs/mini-discs that surround each star and within which the radial accretion flow needs to be directed towards the individual stars.We implement the prescribed acceleration as follows: In the above,  1,2 and  1,2 are the mass of and the radius from the individual binary components and  COM is the radius from the centre of mass of the binary.Figure 2 shows how this binary wind torque prescription varies with position in the computational domain.
Given that our choice about the size of the ballistic zone is somewhat arbitrary and might affect the results, we have undertaken a study of what happens if different values are adopted.This is presented in appendix A, and significant changes in the results are only found to arise when  COM ≥ 3 b .
We perform three simulations defined by three values of  = [2.36× 10 −3 , 6.28 × 10 −4 , 1.23 × 10 −4 ], where these values assume Σ 0 = 1 (which normalises the total disc mass, as defined below).The largest  is equivalent to the mass flux rate of a viscous disc with ℎ = 0.05 and  = 0.1 at  = 1.0.We then choose the intermediate  to be 2 × 10 −4 and the lowest value to be 5 times smaller than this.
In order to study how a wind-driven CBD influences the evolution of the binary orbit compared to a viscously accreting disc, we must also produce equivalent viscous disc models.As  is independent of radius in the wind-driven case, we must ensure we have the same in the viscous case.To have a mass flux that is independent of radius for a Keplerian, viscous disc with a power-law surface density profile, we require Σ() = Σ 0   and () =  0   , where  = −.For the viscous models we implement an alpha viscosity prescription, where  =    (Shakura & Sunyaev 1973).Since the alpha viscosity profile has () ∝  1/2 , our resulting surface density profile must be Σ() = Σ 0  −1/2 .The alpha values that give rise to the chosen  mass flux rates are  = [0.1, 2.6 × 10 −2 , 5.2 × 10 −3 ].The remaining setup description in this section applies to both viscous and winddriven disc models.
We adopt a typical aspect ratio for a protoplanetary disc of ℎ = 0.05.The initial velocity profile takes into consideration the required mass flux through the disc due to the wind or equivalent viscosity, This ensures the disc does not have to relax under the action of the disc torque and allows a quasi-steady state to be reached more quickly.As the boundaries reinforce the initial conditions, the gas flow into the domain will also have the required radial velocity.
We implement an initial relaxation phase consisting of two stages.Over the first 1000 orbits of each simulation we allow the disc to relax with a single star of mass  * = 1.0 at the centre of the domain.This allows the initial disc profile to relax to the mass flux due to either the wind-driven torque or alpha viscosity.We then introduce a low mass secondary, where mass is linearly redistributed from the initial star to the secondary over the next 500 orbits, to reach the required mass ratio of  2 / 1 = 1.0.During this mass redistribution phase the total mass of the binary is kept constant, with  1 +  2 = 1.0, the centre of the domain remains the centre of mass of the two stars and the binary semi-major axis remains constant with  b = 1.0.We will refer to these 500 orbits as the mass ramping phase, and it is introduced to ensure numerical stability of the low viscosity and inviscid simulations, preventing strong shocks which could lead to negative densities and pressures.

TORQUES AND ORBITAL EVOLUTION
Before presenting our results, we first outline how the torques on the binary components are calculated and how they relate to the migration of the binary.The analysis below is very similar to that which has been presented previously, for example see Miranda et al. (2016), Muñoz et al. (2019), Dittmann & Ryan (2021) and Penzlin et al. (2022).We start with the orbital angular momentum of the binary, which can be expressed as where  b is the orbital eccentricity and  b is the reduced mass If we consider circular binaries ( b = 0), differentiate  b with respect to time and then divide the result by  b , we obtain showing that a large negative torque,  b , and/or a large positive accretion rate,  b , lead to rapid shrinkage of the orbit, as expected.
We now want to obtain an expression that allows us to determine when the evolution of a binary transitions from out-spiralling to in-spiralling.Rewriting eqn.( 16) as we obtain Defining the accretion eigenvalue 19), and rearranging terms, gives Hence, a critical value for the accretion eigenvalue is given by ℓ 0 = 3/8 √︁   b  b , below which  b / b will be negative and the semimajor axis of the binary will shrink.
To calculate the accretion eigenvalue, ℓ 0 , we only require the mass accretion rate and the torque exerted by the surrounding gas on the binary.The total torque on the binary due to the gas has two main contributions, gravitational and accretional torques,   =  grav +  acc .These two components are defined as where, The superscripts (1), ( 2) and () in the above refer to the binary components, where  ∈ {1, 2}.The gravitational and accretional torques as well as the mass accretion rate are calculated and saved to file over 600 times per binary orbit, with the accretion eigenvalue being calculated in post production analysis.

NUMERICAL VISCOSITY STUDY
The results of the ring spreading test for an input value of the kinematic viscosity  = 10 −5 is shown in Fig. 1, and it is clear that the simulation follows the evolution predicted by the analytical solution.
The results from our suite of inviscid runs performed at different resolutions are shown in Fig. 3, where the values of the effective viscosity  generated by the numerical diffusion was obtained by fitting the numerical results to eqn. 6 (see also Joseph et al. 2023).We report a value for the numerical viscosity of  = 2.9 × 10 −8 for a resolution 2400 × 2400.For a physical disc with realistic parameters  (ℎ ∼ 0.05), this translates to an  value at  = 1 of  = 1.2 × 10 −5 .The numerical viscosity obtained by Athena++ is found to be very similar to that obtained by the PLUTO code using a similar numerical setup (Joseph et al. 2023), resulting in a numerical viscosity that scales as  num ≃ Δ 2.07 , We note that a resolution of 2400 × 2400 is equivalent to the base grid in our domain containing the inviscid, wind-driven CBD, and that the static mesh refinement region has a grid spacing that is smaller by a factor of two compared to the base grid.Using the relationship shown in Fig. 3, we estimate the refinement region should have a numerical viscosity of  = 7.5×10 −9 .The numerical viscosity in both the base and refinement regions of the domain is therefore small enough to ensure that it is never the dominant driver of gas flow in any of our inviscid, wind-driven disc models.
The adoption of a cartesian grid means that the global angular momentum in the disc should not be conserved to machine precision.
We have examined the level to which the angular momentum is conserved during the highest resolution ring spreading simulation,    and we find that the percentage change is just ∼ 4 × 10 −5 %.Further details are provided in appendix B.

COMPARISON WITH TIEDE ET AL (2020)
In this section we present a comparison to the suite of runs performed by Tiede et al. (2020), which varies Mach number from 10 to 40 for a disc with a constant kinematic viscosity  = √ 2 × 10 −3  2  Ω  .The setup for these runs is described in Section 2.3.
Figure 4 shows images of the surface density of the inner 5 b of the domain after 1000 binary orbits for four chosen Mach numbers.As Mach number increases, the cavity becomes deeper and we see significantly more nonlinear substructure in the form of spiral density waves at the cavity edge and accretion streamers within the cavity.
In Fig. 5 we plot the total mass accretion rate onto both binary components,  b =  1 +  2 , for the four Mach numbers in which surface density snapshots have been presented.As reported in Tiede et al. (2020), we also find modest suppression to the mass accretion onto the binary with increasing Mach number.Note these simulations are for finite disc models, and so it is expected that the mass accretion rates will decline with time, as observed.
The time evolution of the accretion eigenvalue, ℓ 0 , is presented in Fig. 6 as 30-orbit averaged values, ⟨ℓ 0 ⟩ 30 = ⟨  b ⟩ 30 /⟨  b ⟩ 30 .Short term variability in ℓ 0 increases for higher Mach number, however there is no long term evolution.This relationship is also noticed in Tiede et al. (2020).The dashed line is the value of the accretion eigenvalue averaged between orbits 500 and 1000.We compile these averaged values for all 7 runs in Fig. 7, comparing to the results in Tiede et al. (2020) as well as the results from Penzlin et al. (2022), who also carried out a similar study using a polar grid instead of a Cartesian one.We also note that Muñoz et al. (2019) performed a simulation equivalent to the Mach 10 run using an irregular adaptive mesh with AREPO.We plot this single result as a star.
Figure 7 shows that our results are in good agreement with Muñoz et al. (2019) for the Mach 10 case, and as we move to higher Mach numbers our results align with those of Penzlin et al. (2022) while showing the same overall trend with Mach number observed by Tiede et al. (2020).Due to the highly nonlinear nature of the problem, and how changes in resolution and grid set up can alter the value for ℓ 0 , it is not surprising that different studies using different codes and numerical setups produce somewhat varying results.The dotted line indicates the critical value, ℓ 0 = 3/8 √︁      for which circular, equal mass binaries will have  b / b = 0. Comparing all three studies, we see that equal mass binaries are predicted to transition from outspiralling to in-spiralling within the range of Mach 17-25, for discs with a kinematic viscosity of  =

VISCOUS DISC RESULTS
If we are to answer the question is an inviscid wind-driven CBD more effective at driving orbital shrinkage of a binary system than a viscous CBD?, we must first perform simulations in the viscous regime with equivalent mass flux rates through the disc.As described in Section 2.4, these runs are performed with a Mach number M = 20 (disc aspect ratio ℎ = 0.05), with three different mass flux rates  = [2.36× 10 −3 , 6.28 × 10 −4 , 1.23 × 10 −4 ], which correspond to alpha values of  = [0.1,2.6 × 10 −2 , 5.2 × 10 −3 ] evaluated at  =  b .We present the results from the three alpha viscosity simulations in this section.
Snapshots of the surface density distributions at different times near the ends of the simulations for each value of the mass flux are shown in Fig. 8.As expected, the binary carves out a low-density central cavity region where the semi-major axis of the edge of the cavity is ∼ 3 − 5 b (Artymowicz & Lubow 1994).This central cavity region is asymmetric and contains accretion streamers that transport mass to each binary component.The cavity wall is eccentric in all three simulations, which is expected for circumbinary discs where the central binaries are on circular orbits (Papaloizou et al. 2001;Pierens & Nelson 2013;Thun et al. 2017b).In the outer disc, the binary excites spiral density waves which propagate outwards.Figure 9 shows the evolution of the eccentricities and semi-major axes of the cavity edges over the last 1000 binary orbits of each simulation.These two quantities are calculated by drawing a line from the cell with the highest density to the centre of mass.We define the apocentre of the ellipse to be on this line and to coincide with where the density is 10% of the maximum.The line is then extended through the binary centre of mass to the opposite side of the cavity up to the cell with a value of 10% of the maximum density, which we define as the pericentre of the ellipse.The resulting ellipses are plotted in Fig. 8.We see that the run with the largest mass flux has the smallest eccentricity and semi-major axis, the run with the intermediate value of the mass flux has the largest values of eccentricity and semi-major axis, and the run with the smallest value of the mass flux has values of the eccentricity and semi-major axis that are slightly smaller than the intermediate mass flux case.Hence, the behaviour of the cavity is a non-monotonic function of mass flux/viscosity in our runs, although it should be noted that there are large fluctuations in the quantities of interest, and the values for the lowest and intermediate mass flux cases are similar to each other.This relationship contrasts with Dittmann & Ryan (2022), who found, for a Mach 20 disc with circular binaries, a monotonic relationship with cavity edge eccentricity and semimajor axis against viscosity.Although we note their simulations used constant kinematic viscosities rather than the alpha prescription.
Figure 10      gas flows through the disc and onto the stars, unimpeded by the tidal torques due to the binaries.In the case with the lowest mass flux, the accretion rate onto the sinks appears to be larger than the mass flow through the disc, suggesting that the disc may be undergoing long term secular evolution because of the smaller viscosity and hence longer viscous time scale in this run.The opening of a gap or cavity can induce a local change in the gas flow rate in a viscous disc, as highlighted recently by Nelson et al. (2023) in the context of a planet opening a gap in a disc, and this possibly explains the behaviour observed in the bottom panel of Fig. 10.Further discussion of the temporal evolution of the accretion rate onto the binaries is given in Section 8 below, where we compare the viscous and the wind-driven runs.
We plot 30-orbit averaged accretion eigenvalues, ⟨ℓ 0 ⟩ 30 = ⟨  b /  b ⟩ 30 , for all three viscous models in Fig. 11.This averaging is different to the one used in Section 5 in order to be consistent with the work carried out by Penzlin et al. (2022), with whom we compare results to in Section 8.1.The horizontal dashed line is the averaged accretion eigenvalue over the last 1000 binary orbits of the simulation in each case.For the mass flux rates of  = [2.36×10−3 , 6.28×10 −4 , 1.23×10 −4 ], we calculate an average accretion eigenvalue of ℓ 0 = [0.25,0.35, -0.01] √   respectively.All of these values are below the critical value,  0 crit = 3/8, for circular equal mass binaries and therefore correspond to in-spiralling binaries.Initially, there does not appear to be a clear relationship between mass flux and ℓ 0 , however we note that Penzlin et al. (2022) conducted a study comparing different viscosities and their respective ℓ 0 values across multiple disc aspect ratios and binary mass ratios.Their work found that ℓ 0 does not vary monotonically with viscosity, and there is a maximum in ℓ 0 near  ≈ 0.03.Considering the intermediate mass flux simulation in this work corresponds to an alpha viscosity of  = 0.026, we report a similar result.

WIND-DRIVEN DISC RESULTS
In this section we present the results from the wind-driven disc models described in Section 2.4, with mass flux rates of  = [2.36× 10 −3 , 6.28 × 10 −4 , 1.23 × 10 −4 ].In Fig. 12 we plot surface density snapshots taken near the end of the simulation for each value of .
Again the binaries open a low-density eccentric cavity region with semi-major axis ∼ 3 − 5 b into which accretion streamers penetrate after being pulled from the cavity wall by the binary components.
We show the evolution of the eccentricity and semi-major axis of the cavities in Fig. 13.These values are calculated using the same method described in Section 6.We see that the model with the intermediate value of the mass flux has the largest and most eccentric cavity, whereas the model with the smallest mass flux has the smallest and least eccentric cavity, which is not what is seen with the viscous models.Furthermore, comparing each of the corresponding viscous and wind-driven models, we do not see a clear pattern of behaviour in the structures of the cavities.For the highest mass flux, the winddriven model produces a wider and more eccentric cavity, whereas the converse is true for the other two mass fluxes we consider.We caution that the lowest right-hand panel in Fig. 9 and the middle and lowest right-hand panels in Fig. 13 show evidence of residual secular evolution, so it is possible the cavity structures have not fully converged by the end of the runs.
We plot the total binary mass accretion rates,  b =  1 +  2 , for the three wind-driven simulations in Fig. 14.As in Fig. 10, the panels descend with mass flux rates from top to bottom where the horizontal dashed line shows the supplied  through the CBD and the dotted line indicates the total mass accretion rate onto the sink cells.Unlike the viscous runs, a quasi-steady state is not reached immediately in any of the wind driven simulations, instead taking approximately 3000, 5000 and 5000 binary orbits from top to bottom respectively.On longer time scales, however, the mass fluxes through the disc and the mass accretion rates onto the sinks come into decent agreement for each of the runs, indicating that systems are close to being in a steady state where mass flows through the discs and onto the stars at similar rates.This behaviour is quite different to that described by Nelson et al. (2023) for the case of an accreting giant planet in a wind-driven disc.There, the tidal torque of the planet was able to act against the wind-driven torque and effectively halt gas accretion into the gap and onto the planet.In the case of a binary, it appears that the development of a highly dynamic and eccentric disc cavity prevents this mode of evolution from occurring, and instead gas can accrete freely through the cavity and onto the stars.Further discussion of the behaviour of the accretion rates is provided below in Section 8.
Studying Fig. 15 further, we observe the same pattern in all three panels.Each simulations show an initial transient phase where the value of ℓ 0 is steadily decreasing over the period immediately after the binary mass ramping phase.This is followed by a very sharp drop off and recovery.After the recovery period, the accretion eigenvalue eventually levels out to maintain a steady average.The time taken to reach this final steady state increases for decreasing supplied .The level of variability in ℓ 0 over the course of the simulation also increases for decreasing .Further analysis linking the shape of this plot to the morphology of the disc is carried out in Section 9.
We leave the key comparisons and analysis between the viscous and wind-driven models to the following section.

COMPARISON BETWEEN RUNS
Here we compile, compare and analyse the accretion rates and the accretion eigenvalue results from both the viscous and wind-driven disc models.We also investigate the contributions to the accretion eigenvalue and the binary orbital migration rate from four regions in the disc.We label these as follows: inner region, ballistic region, cavity region and outer region.The inner region represents the area within 0.5 b from each individual binary component, therefore containing the mini-discs/CSDs around each star.The ballistic region covers the gas within 1.0 b from the binary centre of mass that is not part of the inner region, corresponding to the area where the wind-driven acceleration is equal to zero, as described in Section 2.4 and indicated in Fig. 2. The cavity region is defined between 1.0 b to 2.5 b in radius from the binary centre of mass, approximately representing the gas up to the outer edge of the cavity (although note that in the case of an eccentric cavity this may not be strictly true).Finally, the outer region corresponds to gas outside the radius of 2.5 b from the binary centre of mass, and therefore contains the main bulk of the outer circumbinary disc.

Viscous and wind-driven disc comparison
Figure 16 compares the accretion eigenvalue averaged over the last 1000 orbits of each simulation for viscous (blue) and wind-driven (orange) simulations.It is clear that a wind-driven accretion disc model produces a significantly lower accretion eigenvalue.We also plot results from Penzlin et al. (2022) in Fig. 16, who performed similar calculations to our viscous simulations.Although the values for ℓ 0 in this work are lower than those in Penzlin et al. (2022), the shape of the plot when moving from higher to lower mass fluxes (or viscosities) remains the same.Both sets of results show a maximum in ℓ 0 near  = 2 × 10 −4 .The slight discrepancy in values between the two results may originate from the different grid setups.Penzlin et al. (2022) carry out simulations on a polar grid with an inner boundary at radius  = 1.0 b and with the binary components orbiting interior to the boundary.The solid circle plotted in Fig. 16 corresponds to the Mach 20 simulation (ℎ = 0.05) we presented in Section 5 when comparing with Tiede et al. (2020).Although the initial conditions between the discs differ (finite disc versus power-law surface density profile), the extracted value for the average accretion eigenvalue is almost identical.

Binary in-spiral timescales
We use the inverse of eqn.(20) to estimate the time for the binary components to merge in each simulation, and we plot the results in Fig. 17.For these calculations we assume an initial separation of 1 AU, a total binary mass of  b =  ⊙ and a total disc mass of  disc = 0.05  b .
The merger times are all in the range ∼ 10 6 -10 7 yr, indicating that binaries that remain on near circular orbits will undergo significant orbital shrinkage within the typical lifetimes of protoplanetary discs, which also fall in the range ∼ 10 6 -10 7 yr (Haisch et al. 2001).Hence, disc-induced shrinkage of binary orbits from an initially more separated state appears to be a viable mechanism for forming close binary systems during the Class I/II phases of pre-main sequence binaries.In future work, we will consider the evolution of live binaries whose orbits evolve under the direct influence of CBDs, in order to see if the above result continues to hold when the eccentricities and semi-major axes of binaries evolve in tandem.
We note that a larger-in-magnitude negative value of ℓ 0 does not always result in faster migration.Equation ( 16) shows that the net torque on the binary and the mass accretion rate both provide important contributions to the rate of migration.So, although ℓ 0 can be used to determine whether a binary orbit should expand or contract, its magnitude should not be automatically equated with the migration speed without also considering the mass accretion rate.Larger-in-magnitude and negative ℓ 0 values can be offset by smaller mass accretion rates (see eqn. 20).An example of this can be observed when comparing ℓ 0 (Fig. 16) and the migration time (Fig. 17) for the smallest  wind-driven model.The mass accretion rate in this simulation is smaller than in the runs with larger mass flux values, so even though ℓ 0 has the most negative value the migration timescale is estimated to be longer than for the other wind-driven disc models.

Contributions to ℓ 0 and 𝑎/𝑎
Figure 18 plots the contributions to ℓ 0 from the different regions of the disc described at the beginning of Section 8.2 for both the viscous and wind-driven disc runs.The effect of accretion onto the stars is also included, and here we note that the relevant points plotted in Fig. 18 combine the effects of both mass accretion and the accretion of momentum.All contributions to ℓ 0 from accretion are within the the tight range of 0.21 and 0.23 across both sets of simulations, hence providing a contribution below ℓ 0 crit = 3/8 that would contribute to the in-spiral of a binary, and we discuss the relative contributions from mass and momentum accretion later in this section.There does not appear to be a strong relation between the accretion contribution and mass flux or disc type (viscous or wind-driven).
The inner region contains the mini-discs/CSDs that surround each individual star, defined on the domain by  1,2 < 0.5 b where  1,2 is the distance from each star, respectively.This region has the largest positive contribution to ℓ 0 and is therefore acting to drive the stars to wider separations.Once again, there is not a strong relationship between mass flux rate and the contribution from the inner regions to ℓ 0 .We find that each wind-driven value is just 5-10% lower than the viscous counterpart with the same mass flux through the disc although, as we discuss below, this small difference does not necessarily reflect the magnitudes of the torques exerted by the mini-discs because ℓ 0 also depends on the mass accretion rate onto the binary, which also varies between the runs.
The ballistic zone contains very little gas across all simulations, and therefore has minimal impact on the orbital evolution of the binaries.The contribution to ℓ 0 due to this small amount of gas across all six runs ranges between values of ℓ 0 ballistic ≈ 0.04 − 0.07.
The region of the disc that produces the strongest negative torque on the binaries, and which therefore is most responsible for causing binary orbits to shrink, is known to be the cavity region (Tiede et al. 2020).Circular binaries are expected to carve out an eccentric cavity with a semi-major axis of  ≈ 2.5 b (Artymowicz & Lubow 1994).We therefore calculate the ℓ 0 contribution from the gas outside the ballistic zone up to a radius of  = 2.5 b .One must be careful not to think of this region as containing the entire cavity, as the eccentricity of the disc causes the cavity wall to extend interior and exterior to a circle of radius  = 2.5 b .It is clear from Fig. 18 that this part of the disc is indeed responsible for the strongest negative torque on the binaries for all simulations.
Figure 18 shows the value of ℓ 0 in the outer regions of the disc,  > 2.5 b .However, as previously mentioned, the outer edge of the cavity will extend into this circular region.Gravitational torque maps of the disc show that the torque from the gas that orbits significantly outside of the vicinity of the cavity wall (i.e. ≳ 7.0 b ) cancels out and so makes very little contribution to ℓ 0 .Thus, it is appropriate to combine the cavity and outer disc regions and analyse their collective effect on the orbital evolution, instead of focusing separately on what we call the cavity and outer regions.It is immediately clear that this combination contributes most heavily to changes in ℓ 0 across the whole suite of simulations and dictates the shape of the relationship between ℓ 0 and  in Fig. 16.
Since the combination of the cavity and outer regions corresponds to all of the gas outside of radius  b , centred on the centre of mass, and that this combination is the most dominant factor responsible for changes to ℓ 0 , understanding the morphology of the disc is key to understanding the influence of this region, and this is considered below in Section 9.
To clarify which effects are responsible for causing the binary orbits to shrink, and to understand why the wind-driven discs are predicted to drive more rapid orbital shrinkage than the equivalent viscous discs, we consider eqn 16 and plot the different contributions to this expression in Fig. 19.The top panel shows the torque contributions to / from each region of the disc described above for both the viscous and wind-driven runs, and here the accretion contribution includes only the accretion of momentum.The bottom panel shows the total torque contribution, the contribution from the mass accretion rate onto the binaries and the sum of these.Note the units associated with the quantities plotted in Fig. 19 are dimensionless code units, as here we are only interested in the relative contributions to /.
Looking at the top panel of Fig. 19, we see the inner regions containing the mini-discs provide a positive torque contribution to /, and that the torque in the viscous case is always significantly larger than in the case of the wind-driven discs.This could be attributed to the fact that the mini-discs in the wind-driven simulations are slightly smaller and less dense than in the viscous case, in part because these models are less effective at driving gas into the cavity region compared to their viscous counterparts (see the discussion below about the cavity structures and Fig. 26), and in part because the wind-driven torque simply removes angular momentum from the mini-discs, and hence causes them to shrink in radius.Viscosity, however, redistributes the angular momentum and causes the minidiscs to expand to their tidal truncation radius.We see also that the accretion of momentum provides a positive torque contribution, and again this torque is larger in the viscous models because of the somewhat higher accretion rates.The outer disc and cavity regions always provide a negative torque contribution, and we see that the viscous Turning now to the lower panel of Fig. 19, we see that in the viscous cases the total torque contribution to / is actually positive in all runs, and it is only the contribution from the mass accretion onto the binaries that ultimately causes the orbits to shrink.In the wind-driven cases, however, both the total torque contribution and the mass accretion act to shrink the binary orbit, explaining why the simulations presented here predict that the in-spiral of binaries should occur more rapidly in wind-driven discs when compared to equivalent viscous discs.

DISC MORPHOLOGY EVOLUTION
We now examine the morphology of the disc in more detail and study how it influences the mass accretion rate, accretion eigenvalue and migration behaviour.
Since the surface density is initialised with a power law profile without an already depleted cavity region, the gravitational torque due to the binary is solely responsible for clearing the cavity.This clearing occurs during the 500 binary orbit mass ramping phase, where the binary mass ratio increases from 0 to 1.When the stars have reached their final masses, in contrast to Figs. 8 and 12, the cavity edge is circular and not eccentric.A transition to an eccentric cavity edge must therefore occur sometime later in the simulation.For the highest  wind-driven case, this transition is seen to occur after between 2000 and 3000 binary orbits, as shown in Fig. 20.The transition from a circular to eccentric cavity edge is accompanied by significant changes in the mass accretion rate and accretion eigenvalue, as seen in Figs. 14 and 15, respectively.This behaviour is observed across all simulations except for the largest  viscous run.In this model the cavity edge becomes eccentric during the initial 500 orbit mass ramping period when the binary is first introduced.The change in disc structure therefore marks the end of the initial transient phase and a quasi-steady state in the measured parameters is reached shortly after.
During the transition from a circular to an eccentric cavity edge, a density lump is also formed during all of the viscous and winddriven simulations.To highlight this feature, Fig. 21 shows snapshots of the surface density for the largest  viscous case covering a time interval of 8 orbits.This lump is observed to orbit around the edge of the cavity with the same period as the local gas, and for our Mach 20 discs the orbital period of the lump is approximately 7 binary orbit periods.This feature has also been observed and commented upon extensively in previous works (e.g.Miranda et al. 2016;Dittmann & Ryan 2022;Wang et al. 2022).At a certain orbital phase, the spiral density waves excited by the binary appear to congregate and impact the same region of the cavity wall.This creates an accumulation of gas at this position and leads to formation of a density lump.The lump then orbits around the cavity edge until it is broken up as it approaches the pericentre of the eccentric cavity.Some of the gas associated with the dispersing lump is directed down the accretion streamers towards the binary components, whilst what is left rejoins the over-dense region at the apocentre as the next density lump is being formed.In all runs this cycle starts as soon as the disc becomes eccentric and continues for the rest of the simulation.The lump becomes less prominent for lower mass flux rates with there also being a reduction in size and density when moving from viscous to inviscid wind-driven discs.
The orbiting density lump has an effect on the short term variability of the mass accretion rate and accretion eigenvalue.As the lump is most prominent in the largest  viscous case, we choose this simulation for further analysis.Figure 22 shows the mass accretion rate and the accretion eigenvalue over a period of 50 binary orbits.The red line in the mass accretion panel is an orbit averaged value, shown in order to reduce the noise and highlight meaningful changes in the mass accretion rate, and the times associated with the vertical dashed lines correspond to the surface density snapshots in Fig. 21.As mentioned above, when the density lump is broken up some gas is directed down the accretion streamers towards the binary components.This can be seen in the mass accretion rate (especially in the averaged value indicated by the red line), where a sharp increase in the amount of mass being accreted occurs just after the density lump has directed material down the accretion streamers.The mass accretion rate then steadily decreases over a period of ∼7 binary orbits until the density lump directs mass down the streamers again.We note this is has been observed in many previous works (e.g.Miranda et al. 2016;Dittmann & Ryan 2022;Wang et al. 2022), and is commonly referred to as the lump-modulated accretion variability.This modulation has typically been shown to have a period of 5 binary orbits, however these works have used larger values for both the viscosity ( = 0.1) and disc aspect ratio (ℎ = 0.1).Fig. 23 plots the mass accretion rate for the Mach 10 run with  = 0.1 from Section 5, showing a modulation period of 5 binary orbits and therefore agreeing with previous works.This is to be contrasted with the 7 binary orbit periods associated with the density lump in the Mach 20 disc models.Clearly the orbital period of the lump depends on the structure of the CBD and hence on the disc parameters.Looking to the bottom panel of Fig. 22, we notice an inverted behaviour in the accretion eigenvalue when we look at the upper envelope of the (not averaged) value of ℓ 0 .Instead of sharply increasing and steadily declining, the envelope increases and suddenly drops (the higher frequency oscillations in this panel are due to the binary's orientation with respect to the eccentric disc, resulting in two oscillations in ℓ 0 per orbit).As the density lump approaches the pericentre, by definition the lump is closer to the stellar components and therefore produces a larger torque on them.The lump is then broken up, dispersing its gas content over a larger area and reducing the torque on the binary.Since the accretion eigenvalue is defined as ℓ 0 =  b /  b , material being directed down the accretion streams onto the binaries during the breakup results in a further decrease in the value of ℓ 0 , and explains why the magnitude of the accretion eigenvalue is the inverse of the mass accretion rate, at least when we look at the envelope associated with the time evolution of ℓ 0 .
The orbit averaged mass accretion rate increases by around 30-40% due to the break up of the orbiting density lump, however the instantaneously measured mass accretion rate can vary by as much as two orders of magnitude.This raises the question of whether this dramatic variation in accretion rate could be observed in the light curves of young binary systems that have circular orbits as considered here.Furthermore, from simple geometric considerations, we would expect an equal mass binary on a circular orbit within the eccentric cavity of a CBD to show variation in the accretion rate with a period equal to half the binary orbit, since this is the time between successive conjunctions between one of the binary components and the pericentre of the disc, which is approximately when accretion streamers are pulled from the cavity wall.In the discussion above, we noted that the time variation of ℓ 0 showed a distinctive signal at half the binary orbit period, primarily because the torque due to the cavity-plus-outer regions of the CBD varies on this timescale.However, closer inspection of the accretion rate on short time scales also shows the expected signal at half the binary orbit period, as illustrated by Fig. 24.To investigate further, we Fourier analysed the accretion rate time series for both viscous and wind-driven disc models, and the resulting power spectra are shown in Fig. 25.As expected, we see a strong signal at low frequency corresponding to ∼ 7 binary orbit periods because of the effect of the density lump.In addition, however, we also see a prominent signal popping up at a higher frequency equal to 2, corresponding to half the binary orbit period.We note that Muñoz & Lai (2016) report accretion rates onto the individual components from their simulation of a binary on a circular orbit (see their Fig.3), and they note there is a half period lag between the accretion rate onto each star reaching repeated maxima, which corresponds to the half orbital period modulation of the total accretion rate we have discussed above.
This modulation of the accretion rate is one phenomenon that could allow comparison between observations of pre-main sequence binaries on circular orbits and the simulations presented here.The well-known systems AK Sco, UZ Tau and DQ Tau are on quite eccentric orbits (see table 3 in Czekala et al. 2019), and so cannot be compared with our simulations.V4046 Sgr, however, consists of nearly equal mass components on circular orbits, and evidence of variation in H emission occurring with a period equal to half of the binary orbital period has been reported by Quast et al. (2000),  consistent with the time variation of the accretion rate obtained in our simulations.So far, however, no variation on the longer time scale of ∼ 7 binary orbits associated with growth and break-up of the density lump have been reported for observations of V4046 Sgr, and identification of such a signal in observational data would provide confirmation that such a feature does indeed arise in close binary systems surrounded by CBDs.Identification of the precise period associated with this feature would also provide strong constraints on the physical properties of the disc in the cavity region, since we have shown that the orbital period of the lump depends on the disc thickness/Mach number.
Finally, in Fig. 26 we show azimuthally averaged surface density profiles which are also averaged in time from surface density snapshots with intervals of 100 binary orbits over the final 1000 binary orbits of each simulation.These demonstrate that the cavity has a steeper outer edge and is deeper for the wind-driven models, as expected because the diffusive behaviour of the viscous models is absent in the wind-driven models.This figure also shows the expected behaviour as a function of the mass flux through the disc, with the lowest  runs showing the deepest cavities.The shapes of cavity outer edges and the depths of cavities near the binaries are consistent with expectations about the torque contributions to / shown in Fig. 19, where it was noted that the torques in the viscous discs were larger in magnitude than in the corresponding wind-driven discs.

SUMMARY AND CONCLUSIONS
We have presented a suite of simulations that examine the interaction between a central binary system and a circumbinary disc, with the primary aim of determining how the interaction changes when angular momentum is transported in the disc by viscosity compared to when it is transported by an external torque arising from a magnetised wind.In this initial proof-of-concept study, we consider equal-mass binary components on fixed circular orbits, and 2D disc models defined on a Cartesian refined-mesh where the magnetised wind torque is mimicked by an external torque prescription.Our main focus is on pre-main sequence binaries, although our results may also have applicability to the evolution of binary black holes.
We begin by calibrating the numerical viscosity, and find that at a fiducial radius of  = 1 the numerical viscosity in our production runs corresponds to  ∼ 10 −5 .Hence it is too small to have a significant impact on the results of our simulations.
We next calibrate the code by running circumbinary disc simulations as a function of the Mach number (or equivalently the aspect ratio) in the disc.We obtain results that are in good agreement with previous studies (Muñoz et al. 2019;Tiede et al. 2020;Penzlin et al. 2022), with the disc developing an inner cavity that is surrounded by a precessing, eccentric circumbinary disc.Furthermore, we find that there is a transition in behaviour when increasing the Mach number above ∼ 20, such that binary orbits are predicted to shrink rather than expand due to interaction with the surrounding gas (Tiede et al. 2020;Penzlin et al. 2022).
The main result of this paper is that for disc models with Mach number M = 20 (ℎ = 0.05), simulations that mimic wind-driven accretion through the disc predict that binary orbits should shrink faster compared to equivalent viscous disc models (i.e.those with the same global mass flux through the disc).This behaviour oc-curs primarily because the positive torque experienced by the binary due to mini-discs surrounding each component is relatively large in the viscous models, slowing the in-spiral rate.Using reasonable estimates for the masses of the binary systems and CBDs that we simulate, the estimated in-spiral times for pre-main sequence binaries on circular orbits with initial separations of ∼ 1 AU are in the range ∼ 10 6 -10 7 yr, indicating that orbital shrinkage through interaction with CBDs of initially well separated binary systems is a viable mechanism for forming short period binaries.
Significant variations in the mass accretion rates and the values of ℓ 0 are observed on short and intermediate timescales for all runs.These arise because of the orbital phasing of the binary components relative to the eccentric precessing disc, and because all runs show the development of a density lump that corotates with disc material at the edge of the cavity once the disc becomes noticeably eccentric, in agreement with previous work (Miranda et al. 2016;Dittmann & Ryan 2022;Wang et al. 2022).The break up of this lump, as it approaches the pericentre of the disc cavity edge, leads to additional gas being pulled into accretion streamers that feed the mini-discs around the individual stars, causing modulation of the accretion rate.The close pre-main sequence binary system V4046 Sgr consists of nearly equal mass components on circular orbits (e.g.Czekala et al. 2019), and hence is a good system to compare with our simulations.Although no variation in accretion rate has been reported that corresponds to the evolution of a density lump, there is evidence of variation in H emission occurring with a period equal to half of the binary orbital period (Quast et al. 2000), consistent with the time variation of the accretion rate obtained in our simulations due to the binary orbiting within a slowly precessing eccentric cavity.
In spite of the simplified treatment of magnetic effects in this proof-of-concept study, it nonetheless lays the groundwork for future simulations that examine the orbital evolution of pre-main sequence binaries due to interactions with their circumbinary discs.In particular, it will be important to examine the evolution of binary orbits using live binaries, since the interaction between disc and binary is bound to change as the orbit itself evolves, especially if the orbital eccentricity grows.In addition, it will also be important to move towards much more sophisticated models of wind-driven discs that actually simulate explicitly the evolution of the magnetic field and its effect on the disc.Recent work that examines the evolution of gas giant planets embedded in magnetised discs shows that the magnetic field and its torque behaves differently near the edge of the tidallyinduced gap compared to torque prescriptions similar to that used in this paper (Lega et al. 2022;Nelson et al. 2023;Aoyama & Bai 2023;Wafflard-Fernandez & Lesur 2023), and an open question is how the magnetic field and its associated torque will evolve in a circumbinary disc?These and other issues will be the subject of future work.

APPENDIX A: BALLISTIC ZONE TESTS
Since the ballistic zone we define for the wind driven disc cases is not a particularly physical representation of the torque experienced by a wind-driven circumbinary disc, we varied the radius where the wind torque is set to 0 in order to test if it has an effect on the evolution of the system.We ran this test on the largest  = 2.35 × 10 −3 case for 4000 binary orbits and choose  bal = [1.5, 2.0, 3.0] b .We plot the 30-orbit averaged mass accretion rate and accretion eigenvalue in Figs.A1 and  A2.These figures also include the first 4000 orbits of the result from Section 7. The averaged accretion eigenvalue between orbits 3750 and 4000 have been calculated to be ℓ 0 = [−0.13,−0.11, −0.14] for ballistic radii of  bal = [1.0,1.5, 2.0] b respectively.So although increasing the ballistic zone radius also increases the time to reach a quasi-steady state, both mass accretion and ℓ 0 show little change in the converged values between  bal = 1.0, 1.5 and 2.0  b .When the ballistic zone has a radius of  bal = 3.0, an eccentric cavity does not form within the simulation time, and the mass accretion rate and ℓ 0 do not reach the same final values.We do not plot ℓ 0 for this case as the fluctuations are so great (between -20 and 10) that it obscures the rest of the plot.From this test we conclude that as long as the size of the ballistic zone is not too large, i.e.. > 3.0 b , the results of the simulations are reasonably well converged.

APPENDIX B: ANGULAR MOMENTUM CONSERVATION FOR VISCOUS RING
The use of a cartesian grid in the simulations means that the angular momentum of an orbiting disc or ring is not expected to be conserved to machine precision.Here we examine the global angular momentum conservation during a simulation of the viscous ring problem described in Section 4. We use a resolution of 2400 × 2400 cells.The angular momentum, J = r × p, where p is the linear momentum, is calculated for each cell between the radii 0.2 and 1.8 and is then summed over all cells.Figure B1 plots the percentage change in the total angular momentum relative to the initial conditions, and we note that a very small change of 4 × 10 −5 % occurs over the duration simulation, which ran for ∼ 450 Ω −1 , where Ω is the angular velocity measured at a radius of unity.The observed increase in angular momentum up until  ≈ 100 is probably due to residual pressure effects  causing the ring to expand slightly, since the time over which it occurs is approximately equivalent to the ring width divided by the sound speed.However, once this phase is over the angular momentum is very well conserved.

Figure 2 .
Figure 2. Magnitude of the wind torque on the inner 5  of the domain, set by the regions defined in Equation (10).

Figure 3 .
Figure 3. Numerical viscosities for different grid resolutions calculated using Athena++.Values taken from our recent comparison between different simulation codes, presented in Joseph et al. (2023), are shown as blue crosses.The red cross is the calculated numerical viscosity for the base grid resolution in this work.The fit shows that the numerical viscosity scales as   = Δ 2.07 .

Figure 4 .
Figure 4. Surface density snapshots after 1000 binary orbits for four values of Mach number.The dashed circle at the centre of the domain represents the binaries orbit.The cavity becomes deeper and the number of spiral density waves increases at higher Mach numbers.

Figure 6 .
Figure 6.30-orbit averaged accretion eigenvalue, ℓ 0 , over the full 1000 orbits.The dashed horizontal line indicates the value of ℓ 0 averaged between orbit 500 and 1000.Short term variability increases with Mach number, however there is no secular evolution.

Figure 7 .
Figure7.Averaged accretion eigenvalue, ℓ 0 , against Mach number.The horizontal dashed line is the critical value, ℓ 0  , below which the binaries semi-major axis contract.ℓ 0 decreases with higher Mach number, crossing the critical value between Mach 15-25.Interactions between the disc and the binary therefore drive the binaries to closer separations for discs with higher Mach number.
plots the mass accretion rate onto the binaries,  b =  1 +  2 , for the full duration of the simulation in each viscous disc case.The panels are defined by their equivalent mass flux rate through the disc, with the highest  in the top panel and the lowest in the bottom.The horizontal dashed lines indicate the values of the supplied  for each simulation, as described in Section 2.4.The

Figure 9 .
Figure 9. Cavity eccentricity (left) and semi-major axis (right) evolution for the set of viscous disc simulations.From top to bottom the panels show results from the highest, middle, and lowest mass flux rate simulations.

Figure 10 .
Figure 10.Binary mass accretion rate,  b =  1 +  2 , for each viscous disc simulation.The horizontal dashed line indicates the supplied  and the dotted line indicates the  into the sinks.The two highest  models reach a quasi-steady state almost immediately after the initial 1500 orbit setup period, whereas the lowest  model experiences an initial transient phase of approximately 3000 binary orbits.

Figure 11 .
Figure 11.Accretion eigenvalue for the suite of viscous disc simulations.The horizontal dashed line is averaged value over the last 1000 orbits of the simulation.An initial transient phase is seen again in the lowest  model.

Figure 13 .
Figure 13.Cavity eccentricity (left) and semi-major axis (right) evolution for the wind-driven disc simulations.

Figure 14 .
Figure 14.Same as Figure 10 but for the inviscid wind-driven simulations.

Figure 15 .
Figure 15.Same as Figure 11 but for the inviscid wind-driven disc simulation.

Figure 16 .
Figure 16.Averaged accretion eigenvalue, ℓ 0 , for all simulations in both the viscous (orange) and inviscid wind-driven (blue) disc regimes.For comparison, results from Penzlin et al. (2022) are plotted in yellow.The solid circle corresponds to the Mach 20 case presented in Section 5, and the dotted line corresponds to ℓ 0 = 3/8.For all mass flux values the inviscid wind model produces an accretion eigenvalue that is significantly less than the viscous counterpart.

Figure 17 .
Figure 17.Migration time calculated using the averaged accretion eigenvalues and mass accretion rates in Equation (20).To calculate these values we have used  b =  ⊙ and  b = 1AU.The disc mass has been scaled to be 5% of the total binary mass.

Figure 18 .
Figure 18.Contributions to the accretion eigenvalue, ℓ 0 , from the regions defined in Section 8 for the viscous (top) and the wind-driven (bottom) disc models.

Figure 19 .
Figure 19.The top panel shows the torque contributions to / in eqn.16 from each region of the disc.Dotted lines: viscous disc runs.Solid lines: wind-driven disc runs.The bottom panel shows the total torque contribution to /, the contribution from mass accretion onto the binary, and the sum of these.Dotted lines: viscous disc runs.Solid lines: wind-driven disc runs

Figure 20 .
Figure 20.Surface density snapshots for the largest  inviscid wind simulation.The circumbinary discs cavity edge can be seen to transition from circular to eccentric.This transition is accompanied by variability in the mass accretion rate and the accretion eigenvalue.

Figure 21 .Figure 22 .
Figure 21.Surface density snapshots of the largest  viscous disc case.A density lump is created during the transition from a circular to eccentric cavity edge and orbits the cavity.It undergoes a cycle of constant creation (top left panel) and destruction (bottom left panel), with some material from the breakup being directed down the accretion streamers (bottom right panel)

Figure 24 .
Figure24.Binary mass accretion rate shown over a shorter time interval for the largest  wind-driven run, in order to illustrate the fact that the accretion rate varies with a period equal to half the binary orbital period.

Figure 25 .
Figure 25.Power Spectrum of the mass accretion.There is a clear peak at 2 indicating accretion variation due to the motion of the binary and a large peak at 1/7 from the modulation due to the lump.

Figure 26 .
Figure 26.1D surface density profiles.The vertical dashed lines at a radius of  = 1 and 2.5 b approximately represents the outer edge of the circumstellar discs and the inner edge of the circumbinary disc respectively.

Figure A2 .
Figure A2.30-orbit accretion eigenvalue for different ballistic zone radii.The ballistic zone with radius 3 b is not plotted as the ℓ 0 variability is very large (between -20 and 10) throughout the full 4000 orbits.

Figure B1 .
Figure B1.Percentage change in the total angular momentum measured between radii  = 0.2-1.8during a simulation of the viscous ring problem.Time is shown in units of Ω −1 at  = 1.