© The Authors 2023

1. Introduction

Dwarf novæ (DNe), which are binary systems composed of a white dwarf accreting from a low-mass companion that fills its Roche lobe, display outbursts due to changes in the flow of matter through the accretion disc surrounding the white dwarf (Smak 1971; Osaki 1974). Like low-mass X-ray binaries composed of a black hole or neutron star accreting from a low-mass companion, DNe cycle between an outburst state and a quiescent state, albeit on a much shorter recurrence timescale of the order of a month.

The disc instability model (DIM; e.g., Lasota 2001; Hameury 2020 for reviews) proposes a mechanism to explain the luminosity variation of these systems, which is related to a gas-ionisation and opacity hysteresis cycle. During quiescence, the disc gas is mostly neutral and marginally optically thick (τ ≈ 1). During this phase, matter from the companion accumulates until a surface density threshold is reached in the disc. At this point, the strong dependence of the opacity on temperature as a result of hydrogen ionisation leads to a thermal runaway that heats the disc. Eventually, the disc ends up hot, fully ionised, and optically thick: the outburst state. In this state, the disc has a higher accretion rate onto the white dwarf than the mass fed at the outer disc from the secondary. The disc density therefore decreases until its density reaches a lower threshold where hydrogen starts to recombine. The typical recurrence time for these systems is around 40 days in total, with about a month for quiescence and a week for the outburst phase (Osaki 1996).

The transport of matter is driven by the redistribution of angular momentum in the disc, historically described as a turbulent viscosity ν = αc_sH, where c_s is the local sound speed and H the local scale height Shakura & Sunyaev (1973). The parameter α depends on the physics that governs the transport of angular momentum. Assuming a radiatively efficient thin disc, it is possible to derive the value of α required to reproduce the outburst cycles. Mineshige & Osaki (1983), Meyer & Meyer-Hofmeister (1984), Smak (1984), and Martin et al. (2019), amongst others, derived an α parameter value of 0.1–0.3 during outburst and 0.01 during quiescence, based on observations.

The magneto-rotational instability (MRI, Velikhov 1959, Chandrasekhar 1961, Balbus & Hawley 1991) is now widely accepted as the source of turbulent angular momentum transport during the hot, outburst phase (Latter & Papaloizou 2012, Hirose et al. 2014). However, the MRI fails to describe the angular momentum transport process during quiescence. During this colder phase, the plasma is too weakly ionised, and so MRI-driven turbulence is expected to be largely reduced (Gammie & Menou 1998). This conjecture has been confirmed through local shearing box models (Scepi et al. 2018), which show that the molecular resistivity of such cold plasma suppresses MRI-driven angular momentum transport. Consequently, we do not yet understand the angular momentum transport mechanism that acts during the quiescent phase.

An alternative route to angular momentum transport in MRI-stable discs is spiral shock waves. Spiral waves are known to be excited by the tidal potential of the binary, and patterns observed in the orbital light curves or Doppler tomograms of cataclysmic variable (CV) accretion discs are consistent with their presence (e.g., Pala et al. 2019; Ruiz-Carmona et al. 2020). Studies using linear (Savonije & Papaloizou 1983; Savonije et al. 1994; Xu & Goodman 2018) and non-linear self-similar analysis (Spruit 1987; Hennebelle et al. 2016) find that the angular momentum transport driven by spiral shocks is strongly correlated with temperature (typically with α ∝ T^3/2). This arises because the temperature controls the disc Mach number (see Sect. 2), which determines the opening angle of spiral shocks and, in turn, the angular momentum flux. Therefore, in low-temperature quiescent discs with Mach numbers of the order of several hundred we expect spiral shock to drive very weak angular momentum transport.

However, the spiral scenario was recently revived by global simulations of close binary systems (Ju et al. 2016; Pjanka & Stone 2020) that emphasise the important role of spiral density waves in the dynamics of these objects. Due to the high computational cost of global models, these simulations focused on hot discs, with Mach numbers of the order of a few tens. Simulating high Mach numbers is challenging, as this requires high resolution to accurately capture the shallow angle between the tightly wound spiral and the Keplerian flow (Matsuda et al. 1990; Rozyczka & Spruit 1993; Blondin 2000; Makita et al. 2000; Kley et al. 2008). Simulations are also typically limited in the radial range and time span over which they follow the evolution of the spiral waves. It has therefore not yet been feasible to determine from direct global simulations whether or not spiral shock waves are a viable angular momentum transport mechanism in quiescence, when the Mach number can reach 900 at the inner edge of the disc (from cold branch temperature estimates of the DIM).

In this paper, we propose to tackle this question using a new GPU-accelerated finite-volume code, IDEFIX, which allows us to explore these quiescent regimes. We present the first 2D global hydrodynamics model with a sufficiently fine resolution to resolve the spiral waves in a cold quiescent disc. The code is also fast enough for us to explore timescales well beyond transients due to relaxation from the initial state. In Sect. 2, we review the analytical formulation and present the code used to solve the equations. In Sect. 3, we show two-dimensional simulations of isothermal discs with realistically high Mach numbers. We discuss the consequences of our findings in Sect. 4 and present conclusions in Sect. 5.

2. Methods

In this paper, we use the IDEFIX¹ (Lesur et al. 2023) finite-volume astrophysics code to solve the hydrodynamic Euler equations in polar coordinates (R, ϕ). More details on the code are given in Sect. 2.3.

2.1. Analytical formulation

We solve the Euler equations for hydrodynamics:

(1)

(2)

for a fluid of density ρ, velocity v, with pressure p and in a gravitational potential Ψ. The above system is closed with the ideal gas equation of state:

(3)

The runs presented here are isothermal: the temperature and sound speed c_s are constant over time and space. Here, since we only present 2D simulations, the volume density ρ and the pressure p are vertically integrated. We write these vertically integrated quantities Σ and P, respectively.

We work in a rotating frame centred on the white dwarf in which the secondary star is fixed. For our binary systems, the tidal potential Ψ in the rotating frame is given by

(4)

where M_WD and M_s are the masses of the white dwarf and secondary star, respectively. G is the gravitational constant, and a is the binary separation. The third term of this potential accounts for the non-inertial reference frame. The inertial forces are added by the solver. In practice, the centrifugal force is added as a radial source term with v_ϕ being the total azimuthal velocity (i.e. taking both the rotation of the reference frame and a possible advection velocity). The Coriolis force is included as a modified inter-cell momentum flux to guarantee angular momentum conservation at machine precision (Mignone et al. 2012).

2.2. Grid and units

For all simulations in this paper, we use a logarithmic grid in the radial direction, such that ΔR ∝ R and a uniform grid in the azimuthal direction. The former allows a finer resolution of the inner edge of the disc, where the spiral shock structures are smaller (Spruit 1987). We set the inner boundary of the disc at r₀ = 0.01 a, which corresponds to the expected radius of the white dwarf. The outer radius of the integration domain is chosen to be the radial position of the Lagrange point L₁ of the binary system. With these definitions, the outer radius changes with the mass ratio q = M_s/M_WD. However, we choose the resolution such that going to a further radius only appends points to the grid, with no modification of the inner region (this concerns the runs Ma = 250 from q = 0.3 to q = 0.1 and q = 0.7).

In this paper, we choose the time unit as the binary period . The length unit is the binary separation a. Following Ju et al. (2016), these units are scaled on the dwarf nova system SS Cygni, where T₀ is approximately 6.6 hours and a is 1.37 × 10¹¹ cm (Bitner et al. 2007).

We define the Mach number at the inner edge as

(5)

where is the Keplerian angular frequency around the white dwarf. We note that this definition depends on the chosen inner boundary of the simulation; for example, in an isothermal setup . This must be taken into account when comparing with previous work. In particular, Ju et al. (2016) set their inner boundary further out at r₀ = 0.02 a.

2.3. Algorithm

The simulations presented in this paper all use the IDEFIX code, which is a finite-volume conservative Godunov grid-based code with a structure similar to Pluto (Mignone et al. 2007). The code is written in C++17 and uses the Kokkos portability library (Edwards et al. 2014) for many-core shared memory parallelisation and MPI for distributed memory parallelisation, allowing for high performance on most available architectures.

We use the HLLC Riemann problem solver (Harten et al. 1983; Toro et al. 1994). We also use the Fargo algorithm (Masset 2000) with a Keplerian advection velocity to speed up the integration, as implemented in Mignone et al. (2012). We verified that we are able to reproduce the work of Ju et al. (2016) with this code (see details in Appendix B).

2.4. Boundary conditions and density floor

In the azimuthal direction, we use periodic boundary conditions while we use ‘Keplerian outflow’ boundary conditions in the radial direction. These last boundary conditions are identical to standard ‘outflow’ boundaries, with the exception of the azimuthal velocity, which is set to its Keplerian value in the ghost zones. We additionally use a density floor to limit the density contrast in the simulation. If the density of one cell drops below the threshold Σ_floor = 10⁻⁶, then this value is used as density instead. The density floor is implemented in a total momentum conserving fashion, so as to prevent the injection of angular momentum by the procedure.

2.5. Initial conditions

We take the initial profile of the disc to have a uniform surface density Σ = 1. The initial velocity profile is Keplerian with and no initial radial velocity. This state is not an equilibrium state, because the initial velocity only takes into account the gravitational influence of the central white dwarf, not the companion.

We pretruncate the disc at the outer edge to shorten the time needed for the outer disc to settle down. The initial truncation radius is chosen to be slightly larger than the expected r_max of Paczynski (1977). The truncation is performed using the following mask applied to the initial density profile:

(6)

where r_T is the chosen truncation radius and δ_T = 0.01 is the chosen truncation width. We note that this truncation is not perfect as it is axisymmetric, unlike the predicted last stable disc orbit (Paczynski 1977). At larger radii, we enforce the previously described density floor. We typically wait for the system to relax over a few binary orbital periods before measuring any quantity.

2.6. Averaging methods

In this work, the dynamical timescale is ≃ 1/10 000th of the binary period, while we are following the evolution of the system for dozens of binary orbits. To capture the evolution in a computationally efficient way, we used an ‘on-the-fly’ time and azimuthal averaging scheme: we sample the simulation every 1/10 000th of a binary orbit and average over 100 of these samples, producing an output each 1/100th of a binary orbit. This is sufficient to capture the rapid inner-edge dynamics.

3. Results

In this section, we present globally isothermal models with realistically low temperatures (high Mach numbers).

3.1. Resolution

We chose the resolution such that the expected spiral wave is well resolved everywhere in the disc. Linear theory predicts that the wavelength at the inner edge of the disc is (Savonije et al. 1994)

(7)

where m is the azimuthal wave number. For the m = 2 mode that is expected to be dominant (Savonije et al. 1994), a resolution of N_R = 1000 gives a minimum of 2.5 radial points per wavelength for our maximum Mach number Ma = 370 (Fig. A.1). We consider this sufficient to resolve the entire spiral in the disc. We use about as many points in the azimuthal direction as in the radial direction to keep a cell aspect ratio of , that is, of order unity. Therefore, the resolution of our isothermal runs is 1081 radial points and 1024 azimuthal points (Table 1). The run with Ma = 250 and q = 0.1 has a larger radial extension and consequently larger N_R to keep the grid spacing unchanged in the inner disc. The opposite is true for the run with q = 0.7.

3.2. Spiral waves

As expected from previous works (Sect. 1), we see that spiral waves develop in the accretion disc due to the tidal potential (Fig. 1). To get a better grasp of the excited modes, we show the Fourier decomposition of the density in Fig. 2. The plotted values correspond to the amplitude of the Fourier expansion of the relative density fluctuation smoothed over 50 radial points (moving average), defined from

(8)

Fig. 1. Density maps after ten binary orbits for isothermal runs with different Mach numbers. (a) Ma = 80. (b) Ma = 140. (c) Ma = 250. (d) Ma = 370.

where is the azimuthal average of the surface density. Figure 2 shows strong m = 1 and m = 2 spiral modes are present in all runs. The m = 2 mode predicted by linear theory dominates only at low Mach numbers and in the outer disc. The simulations show that the m = 2 mode is initially excited at the outer boundary, propagates through the disc, and is reflected when it reaches the inner edge. The strong m = 1 mode becomes visible in density maps after this reflection together with the formation of an eccentric inner cavity.

Fig. 2. Azimuthal spectrum of the relative density fluctuation distribution for modes m = 1 (red solid line) and m = 2 (green dashed line) after ten binary orbits. (a) Ma = 80. (b) Ma = 140. (c) Ma = 250. (d) Ma = 370.

To investigate the impact of the chosen inner boundary conditions on the appearance of the m = 1 mode, we performed an additional run at Ma = 250 and q = 0.3 implementing a wave-killing zone close to the inner edge extending from r₀ to r_WKZ. In this zone, we relax all physical quantities to the initial Keplerian state of the disc, with constant surface density and no radial velocity. The associated relaxation timescale is set to one-tenth of the orbital time at the inner boundary. For a wave-killing zone to be effective, one needs to choose its spatial extension to be a few wavelengths of the mode one wishes to neutralise (see Fig. A.1). In our case, we first take r_WKZ = 1.2r₀, so that the m = 2 spiral mode is cancelled at the inner edge. However, this had no significant impact either on the dynamic or on the eccentric cavity at the inner edge. Using a more extended wave-killing zone with r_WKZ = 2r₀ yields the same results. We conclude that there is more at play than a mere reflection exciting the m = 1 mode, most likely because we do not fall within the linear theory regime (see Sect. 4).

3.3. Measured α

We measure the effective angular momentum transport parameter α (Shakura & Sunyaev 1973) in our simulations as

(9)

where F_R(v_ϕ) is the inter-cell angular momentum flux in the radial direction R, which can be written , where is the departure from the Keplerian velocity. The averages are computed over the whole azimuthal domain ϕ ∈ [0,2π], and over a few binary orbits. The exact time span over which the average is carried out will always be specified in the following. We use the inter-cell angular momentum fluxes computed by the Riemann solver in the integration loop rather than the reconstructed cell-centred values. This allows us to better capture the small-scale dynamics. We compare our method to measure α with that of Ju et al. (2016) in Appendix C.

Figure 3 shows the time evolution of α in the disc. In all simulations, we first see a transient regime at early times (t ≲ 20) with a high α that slowly decays. At later times (t ≳ 40), the transport parameter in the disc drops to α ≪ 10⁻² in all simulations, regardless of the mass ratio q or the Mach number. The relaxation timescale appears to be dependent on the Mach number. For example, the simulation with Ma = 370 shows a slower relaxation than the Ma = 250 run (bottom panel of Fig. 3). The state reached after 80 orbits may still not be a steady state and α is still decreasing. We did not go beyond this stage, as α is already much lower than required by the DIM.

Fig. 3. Measured angular momentum transport parameter. Top panel: measured angular momentum transport parameter α azimuthally averaged and smoothed over three binary orbits for the run Ma = 250 and q = 0.3. Here we plot the absolute value of the transport parameter, and note that its sign changes in the spiral waves (oscillations at R ≳ 0.1) and is mostly positive everywhere else. Bottom panel: azimuthally and radially (over 0.02 ≤ R ≤ 0.1) averaged α, smoothed over three binary orbits for three runs.

We caution that regions of very high α at the inner and outer edges are unphysical regions where the density floor is triggered because of the growth of an eccentric inner cavity (Sect. 3.2) or the tidal truncation of the outer disc. We also note that matter accumulates at the outer radius of the disc at lower temperatures. This is a result of relaxation of the disc truncation from its initial state, which assumes a Keplerian velocity profile and an analytical guess for the outer radius. This residual matter is not efficiently redistributed by the spiral waves given the weaker angular momentum transport at lower temperatures.

Figure 4 shows the radial distribution of α at late times for simulations with different Mach numbers, keeping q = 0.3. The plotted values are both azimuthally- and time-averaged. The strong oscillations around R = 0.2 − 0.3 reflect the effect of the tidal torque on the spiral waves. Here, unlike in Ju et al. (2016), the torque term of the angular momentum conservation equation does not exactly balance the Reynolds stress term. As a consequence, the measured accretion rate and the measured α are not exactly zero, but rather oscillate, which creates regions of slightly negative accretion rate. The efficiency of transport decreases with disc temperature, and in all cases we get α ≪ 10⁻². We do not find evidence for a radial steepening of α with higher Mach numbers (e.g., Blondin 2000).

Fig. 4. Measured angular momentum transport parameter α at late times (between t = 80 and t = 89.9) for different Mach numbers with q = 0.3, averaged over time and azimuth.

Figure 5 shows how α changes with varying mass ratio q, keeping Ma = 250. Again, we obtain α ≪ 10⁻² at late times. However, unlike the runs at q = 0.3 and q = 0.7, the transport appears to stabilise at a higher value of α ≈ 2 × 10⁻³ for q = 0.1. This is consistent with a possible coupling with the 3:1 resonance, which is inside the disc for q < 0.3 (Lubow 1991).

Fig. 5. Measured angular momentum transport parameter α at late times (between t = 80 and t = 89.9) for different mass ratios with Ma = 250, averaged over time and azimuth.

3.4. Convergence study

In order to make sure that the resolution of our runs is sufficiently fine, we performed a higher-resolution run of our coolest setup. For this simulation, we increase the resolution by a factor of four, with a grid of size 4096 × 4096. This resolution is much finer than what is required to resolve the m = 2 spiral wave from linear theory. However, the theoretical estimate may not hold given that we are far from the linear regime (see Sect. 4). We observe the same spiral wave structure in the high-resolution run as in the low-resolution run, albeit with a slower decrease in α. This is likely a consequence of the decreased numerical diffusivity at higher resolution. The relaxation timescales that we observe in our runs (Fig. 3) are therefore partly controlled by the numerical diffusivity. We reach α ≪ 10⁻² after 100 binary orbits, a value consistent with the lower-resolution runs, confirming that these low values of α are not an artefact of resolution.

4. Discussion

4.1. One-armed spirals

Our simulations show that the m = 1 modes are excited concomitantly with m = 2 spiral waves. Previous works indicated that m = 1 modes could be excited in a close binary system provided (1) that the mass ratio is low enough (q < 0.3) for the 3:1 resonance radius to lie inside the disc and (2) that the disc is viscous, so as to couple the m = 2 and m = 3 modes to excite the m = 1 modes (Heemskerk 1994; Stehle 1999; Kornet & Rozyczka 2000; Kley et al. 2008). Here, we satisfy neither of these criteria: we do not include viscosity and our results hold for both q = 0.3 and q = 0.7. We do not observe a strong global eccentricity growth of the disc as Kley et al. (2008) do: the eccentricity mode is limited to the inner disc region, keeping the overall shape of the outer disc circular (see Appendix D). The emergence of this inner eccentric cavity may be of relevance to observations that suggest the presence of a truncated disc (e.g., Balman & Revnivtsev 2012).

Our results also hold when including a wave-killing zone at the inner boundary to dampen the m = 2 reflection. This suggests a less efficient, non-linear coupling is at play in our simulations that may have been missed in lower-resolution simulations with lower Mach numbers. Following the analytical development from Savonije et al. (1994), we can show that linear theory is expected to fail for the realistic temperature regimes that we explore here. From the general solution derived in this latter paper, the amplitude of the density perturbation of the solution Σ′ compared to the unperturbed density profile Σ₀ can be quantified as

(10)

where ω is the binary angular frequency. For a quiescent disc temperature of 3000 K and with typical masses for the white dwarf and its companion star, we obtain σ(r₀)≃1 at the inner edge. As σ(r)∝r², non-linear effects increase at larger radii. The solution derived by Savonije et al. (1994) is based on a perturbative WKB solution where σ ≪ 1 and with the assumption that r ≪ a. The emergence of a one-armed spiral wave could not have been predicted by linear theory, as the regime probed here strongly deviates from this approximation.

4.2. Angular momentum transport

The initial high value of the α parameter mentioned in Sect. 3.3 is likely a consequence of our initial condition, which assumes circular Keplerian orbits for the gas. This is clearly not an equilibrium when a secondary star is present. Relaxation from these conditions takes more than 20 orbits. This is longer than for example the time over which Ju et al. (2016) ran their simulations, meaning that their values of α may still be impacted by their initial conditions (which are similar to ours). However, the relaxation time is only a few binary orbits at the lower Mach numbers of their simulations.

On long timescales, our simulations show that α is always much smaller than the typical α ≈ 0.01 in quiescence required by the comparison of DIM models to observed light curves. The values we reach are of the order of α ≈ 10⁻⁴ and are still decaying when we end the simulation. The one exception is for q = 0.1 where transport stabilises to a larger α ≈ 10⁻³, probably because the 3:1 resonance allows waves inside the disc to be excited. This can only apply to systems with q < 0.3, such as the WZ Sge subclass (Kato 2015), while DNe show mass ratios in the range of 0.1 ≲ q ≲ 0.8 (Otulakowska-Hypka et al. 2016).

In our simulations with the highest Mach number, the relaxation time is comparable to the total integration time of 100 binary orbits, which is on the order of the recurrence time between two normal outbursts in DNe systems (Otulakowska-Hypka et al. 2016). Spiral waves therefore might still play a role in transporting angular momentum in quiescence if the relaxation to low α is slow enough. There are two caveats that make this scenario unlikely. First, our convergence run indicates that the relaxation time is likely to depend on the numerical resolution, and so deducing a reliable decay timescale from the simulation is not straightforward. Second, this would depend critically on how the disc relaxes from MRI-driven transport in outburst. Our initial conditions assume an isothermal disc with a flat surface density profile. A constant-temperature disc is a reasonable approximation for a DIM quiescent disc. However, the surface density profile is expected to increase with radius (Σ ∝ R) when the disc enters quiescence due to the dependence on radius of the critical densities associated with the temperature hysteresis (Lasota 2001). This may introduce differences with our flat Σ profile and, more generally, induce a complex interaction in the time evolution of the α and Σ density profiles.

4.3. RWI/KHI spiral excitation

We also verified that the spiral-driven accretion that we observe is indeed a consequence of the presence of a secondary star. Lesur et al. (2015) showed that an infalling envelope around the disc could trigger a mixture of a Rossby Wave instability (Lovelace et al. 1999) and centrifugal instability at the disc outer edge, creating spiral perturbations that propagate inwards. The outer edge of our disc fulfils the instability criteria of these two instabilities, in part because of the density floor that induced an inflow from the ‘void’ surrounding the disc, but also because of the sharp truncation caused by the secondary.

To check the impact of these edge instabilities, we simulated the same discs without the companion, keeping the disc truncated density profile. We find that in this case, these instabilities only drive transport α ≪ 10⁻⁴ in the disc bulk and are therefore of negligible importance in the present work (see Fig. E.1).

However, it is possible that the tidal potential of the secondary star selects and excites a specific mode of the edge instabilities. This could be a possible mechanism to explain the unexpected apparition of m = 1 spiral modes.

4.4. Expected instabilities in 3D

Several instabilities expected in 3D are absent from our 2D approach, among which the vertical shear instability (VSI, Nelson et al. 2013) and the convective over-stability (COV, Klahr & Hubbard 2014). The existence of these two instabilities is mostly controlled by the dimensionless cooling timescale Ω_Kτ_cool where τ_cool measures the typical thermal relaxation timescale resulting from the heating/cooling equilibrium. In the case of CVs, we expect Ω_Kτ_cool ≪ 1, a regime in which we expect the VSI to be present but not baroclinic instabilities (Lesur et al. 2022).

The question is then how much angular momentum transport is to be expected from the VSI. Because the VSI is driven by the vertical shear of the flow, it depends strongly on the disc thickness. Using high-resolution 3D hydrodynamical simulations, Manger et al. (2020) propose a scaling for the VSI α ≲ (H/R)^2.6, which implies α < 10⁻⁵ in a quiescent DNe disc subject to VSI turbulence. This means that VSI turbulence, if present, is probably negligible as an angular momentum transport process in these systems.

4.5. Accretion stream impact and thermodynamics

The contribution of the matter flux coming from the companion star (called accretion stream) is expected to be the strong heating at the ‘hot spot’ where it meets the disc. As a consequence, we did not include it in this isothermal study. Modelling the thermodynamics of the disc would allow us to know the heating efficiency of the spiral shocks presented in this work. Previous works such as that of Ju et al. (2016) showed that without any cooling, and using γ = 5/3, the disc heats up quickly. This suggests that some parts of the spiral-heated disc might still be unstable to MRI even during quiescence.

5. Conclusion

We present the first 2D hydrodynamics simulations of a DN disc in realistic quiescence temperature regimes. Our main results can be summarised as follows.

1. The dominant spiral mode is a one-armed m = 1 mode once the disc has settled. Linear theory, which predicts an m = 2 mode, fails at low temperatures (high Mach number), as expected.

2. Spiral-wave-driven accretion in realistic temperature regimes is much lower than the values of α required by the DIM. Running the simulations clearly shows that, regardless of the initial transport, there is decays as the disc relaxes from its initial conditions.

We conclude that hydrodynamic spiral waves appear unable to explain accretion in quiescent dwarf novæ and low mass X-ray binaries (whose outer discs sample similar densities and temperatures).

Acknowledgments

The authors would like to thank the anonymous referee for their comments that helped them improve this work. GL acknowledges support from the European Research Council (ERC) under the European Union Horizon 2020 research and innovation program (Grant agreement No. 815559 (MHDiscs)). The simulations and data reduction used in this article were performed using the GRICAD infrastructure (https://gricad.univ-grenoble-alpes.fr), which is supported by the Grenoble research communities, as well as the Jean-Zay supercomputer at IDRIS http://www.idris.fr/jean-zay/. This work was granted access to the HPC resources of IDRIS under the allocation 2021-A0120402231 made by GENCI. Data reduction for this paper was performed by making extensive use of the SciPy (Virtanen et al. 2020), NumPy (Harris et al. 2020), pandas (Reback et al. 2022) and Matplotlib (Hunter 2007) Python libraries. IDEFIX uses the KOKKOS portability tool (Edwards et al. 2014).

References

Appendix A: Linear theory wavelength

The wavelength from linear theory is shortest at the inner edge of the disc (Eq. 7). Figure A.1 shows how this wavelength compares with the grid spacing for our ‘low-resolution’ runs. Our resolution (with 1081 radial points) is sufficiently fine up to much higher Mach than explored in this work. In practice, the m = 1 mode that develops has a much wider wavelength. This plot also shows that the width of the wave-killing zones implemented at the inner boundary is much wider than the linear theory m = 2 wavelength by several orders of magnitude (see Sect. 3.2).

Fig. A.1. Grid spacing versus analytical wavelength for the spiral mode m = 2 as a function of the radius.

Appendix B: Reproducibility with IDEFIX

We verified that we could reproduce the results of Ju et al. (2016) with the code IDEFIX. Figures B.1, B.2, and B.3 can be compared with their Figure 12. Although we only present isothermal simulations in the main text, we also compared non-isothermal runs with an adiabatic index of γ = 1.1. The general agreement is very good, with small differences likely due to variations in implementation between IDEFIX and Athena++.

Fig. B.1. Density profile after four binary orbits, averaged over four binary orbits (t = 25 to t = 50 using the units of Ju et al. 2016).

Fig. B.2. Same as figure B.1 but for the accretion rate.

Fig. B.3. Same as figure B.1 but for the effective α.

As the simulations presented in this paper do not use the same time units (by a factor of 2π), the units in this section are converted to the same units as in Ju et al. (2016). We also note that: (1) the α parameter plotted on Fig. B.3 is computed in the same way as in Ju et al. (2016), i.e. from the radial mass flux; (2) our simulation box extends down to R = 0.02 but the plot only extends down to R = 0.04 to ease comparison with Ju et al. (2016); and (3) the Mach numbers cannot be directly compared because of differences in the location of the inner boundary (where it is defined).

Appendix C: Definition(s) of α

The angular momentum transport parameter can be defined in several ways, and not all definitions are equivalent, depending on the exact system. Ju et al. (2016) argue that for binary systems such as the one studied here, this parameter should be defined from the accretion rate rather than from the Reynolds (and, in MHD, Maxwell) stress. It has been claimed that this effective transport parameter α_eff better reflects the influence of the tidal torque of the companion star. Ju et al. (2016) take

(C.1)

where the accretion rate Ṁ is taken from the inter-cell fluxes computed by the Riemann solver at each time step. We chose to measure α from the angular momentum flux (Eq. 9). Both definitions can be obtained from the angular momentum balance equation, using several approximations. In the case of the definition from the radial mass flux, there is a clear contribution from the tidal torque. In our simulations, the disc is clearly not in steady state, and so a definition of the transport parameter based on the accretion rate Ṁ cannot easily be decomposed into a Reynolds stress term (e.g. our α) plus a tidal torque term, as this would require Ṁ to be constant with respect to the radius. In practice, there is indeed a slight difference in the values computed using these two methods; however, these differences do not change the conclusions of this paper. As shown in Figure C.1, regardless of the chosen definition, we obtain α ≪ 10⁻² at late times.

Fig. C.1. Angular momentum transport parameter α averaged over azimuth and ten binary orbits for the run Ma = 250 and q = 0.3. The definition from fluxes (this paper) is plotted with a solid red line and the definition from radial mass flux (Ju et al. 2016) is plotted as a green dashed line.

Appendix D: Pitch angle measurements

The pitch angle θ of a spiral R(ϕ) is defined as . From the Fourier decomposition of Eq. (8), we can define the spiral associated to each mode m as (up to a constant phase). Therefore,

(D.1)

Figure D.1 shows the spiral pattern on a deprojected plot. Consistently with Fig. 2, we see that the one-armed spiral wave dominates the dynamic. The left panel shows that pitch angle of the spiral wave is consistently < 5° in the disc. The spiral waves are tightly wound. This explains how, even though the m = 1 mode dominates, the disc remains circular in shape overall, rather than very eccentric.

Fig. D.1. Density profile properties of the run Ma = 250 after 38.8 binary orbits. Left panel: Pitch angle of the m = 1 and m = 2 spiral modes in the disc in the density profile. Right panel: De-projected density fluctuation , as defined in Eq. (8).

Appendix E: RWI/KHI-driven accretion

As discussed above, we verified that the mixture of RWI/KHI and centrifugal instability drives negligible accretion compared to the spiral shocks triggered by the tidal potential of the secondary star. The setup used to produce Fig. E.1 is the same as for the run with Ma = 250 and q = 0.1. The mass of the white dwarf is the same but the secondary mass M_s is set to zero in the potential (Eq. 4). This setup is very similar to the setup of Lesur et al. (2015).

Fig. E.1. Measured angular momentum transport parameter. Top panel: Measured angular momentum transport parameter α azimuthally averaged and smoothed over three binary orbits for a run with Ma = 250 and no companion star. Here we plot the absolute value of the transport parameter. Bottom panel: Azimuthally and radially (over 0.02 ≤ R ≤ 0.1) averaged α, smoothed over three binary orbits for three runs. Compared to Fig. 3, the α scale of the bottom panel has been extended to show lower values of the transport parameter.

All Tables

All Figures

	Fig. 1. Density maps after ten binary orbits for isothermal runs with different Mach numbers. (a) Ma = 80. (b) Ma = 140. (c) Ma = 250. (d) Ma = 370.
In the text

	Fig. 2. Azimuthal spectrum of the relative density fluctuation distribution for modes m = 1 (red solid line) and m = 2 (green dashed line) after ten binary orbits. (a) Ma = 80. (b) Ma = 140. (c) Ma = 250. (d) Ma = 370.
In the text

Fig. 3. Measured angular momentum transport parameter. Top panel: measured angular momentum transport parameter α azimuthally averaged and smoothed over three binary orbits for the run Ma = 250 and q = 0.3. Here we plot the absolute value of the transport parameter, and note that its sign changes in the spiral waves (oscillations at R ≳ 0.1) and is mostly positive everywhere else. Bottom panel: azimuthally and radially (over 0.02 ≤ R ≤ 0.1) averaged α, smoothed over three binary orbits for three runs.

In the text

	Fig. 4. Measured angular momentum transport parameter α at late times (between t = 80 and t = 89.9) for different Mach numbers with q = 0.3, averaged over time and azimuth.
In the text

	Fig. 5. Measured angular momentum transport parameter α at late times (between t = 80 and t = 89.9) for different mass ratios with Ma = 250, averaged over time and azimuth.
In the text

	Fig. A.1. Grid spacing versus analytical wavelength for the spiral mode m = 2 as a function of the radius.
In the text

	Fig. B.1. Density profile after four binary orbits, averaged over four binary orbits (t = 25 to t = 50 using the units of Ju et al. 2016).
In the text

	Fig. B.2. Same as figure B.1 but for the accretion rate.
In the text

	Fig. B.3. Same as figure B.1 but for the effective α.
In the text

	Fig. C.1. Angular momentum transport parameter α averaged over azimuth and ten binary orbits for the run Ma = 250 and q = 0.3. The definition from fluxes (this paper) is plotted with a solid red line and the definition from radial mass flux (Ju et al. 2016) is plotted as a green dashed line.
In the text

	Fig. D.1. Density profile properties of the run Ma = 250 after 38.8 binary orbits. Left panel: Pitch angle of the m = 1 and m = 2 spiral modes in the disc in the density profile. Right panel: De-projected density fluctuation , as defined in Eq. (8).
In the text

Fig. E.1. Measured angular momentum transport parameter. Top panel: Measured angular momentum transport parameter α azimuthally averaged and smoothed over three binary orbits for a run with Ma = 250 and no companion star. Here we plot the absolute value of the transport parameter. Bottom panel: Azimuthally and radially (over 0.02 ≤ R ≤ 0.1) averaged α, smoothed over three binary orbits for three runs. Compared to Fig. 3, the α scale of the bottom panel has been extended to show lower values of the transport parameter.

In the text