Issue 
A&A
Volume 677, September 2023



Article Number  A10  
Number of page(s)  12  
Section  Astrophysical processes  
DOI  https://doi.org/10.1051/00046361/202346090  
Published online  24 August 2023 
Spiralwavedriven accretion in quiescent dwarf novæ
Univ. Grenoble Alpes, CNRS, IPAG, 38000 Grenoble, France
email: marc.vandenbossche@univgrenoblealpes.fr
Received:
6
February
2023
Accepted:
25
April
2023
Context. In dwarf novæ (DNe) and lowmass Xray binaries, the tidal potential excites spiral waves in the accretion disc. Spiralwavedriven accretion may be important in quiescent discs, where the angular momentum transport mechanism has yet to be identified. Previous studies were limited to unrealistically high temperatures for numerical studies or to specific regimes for analytical studies.
Aims. We perform the first numerical simulation of spiralwavedriven accretion in the cold temperature regime appropriate to quiescent discs, which have Mach numbers ≳100.
Methods. We used the new GPUaccelerated finitevolume code IDEFIX to produce global hydrodynamics 2D simulations of the accretion discs of DN systems with a sufficiently fine spatial resolution to capture the short scaleheight of cold, quiescent discs with Mach numbers ranging from 80 to 370.
Results. Running the simulations on timescales of tens of binary orbits shows transient angular momentum transport that decays as the disc relaxes from its initial conditions. We find the angular momentum parameter α drops to values of ≪10^{−2}, too weak to drive accretion in quiescence.
Key words: novae / cataclysmic variables / binaries: close / accretion / accretion disks / hydrodynamics
© The Authors 2023
Open Access article, published by EDP Sciences, under the terms of the Creative Commons Attribution License (https://creativecommons.org/licenses/by/4.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.
This article is published in open access under the Subscribe to Open model. Subscribe to A&A to support open access publication.
1. Introduction
Dwarf novæ (DNe), which are binary systems composed of a white dwarf accreting from a lowmass 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 lowmass Xray binaries composed of a black hole or neutron star accreting from a lowmass 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 gasionisation 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_{s}H, 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 & MeyerHofmeister (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 magnetorotational 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 MRIdriven 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 MRIdriven 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 MRIstable 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; RuizCarmona et al. 2020). Studies using linear (Savonije & Papaloizou 1983; Savonije et al. 1994; Xu & Goodman 2018) and nonlinear selfsimilar 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 lowtemperature 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 GPUaccelerated finitevolume 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 twodimensional 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^{1} (Lesur et al. 2023) finitevolume 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:
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:
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
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 noninertial 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 intercell 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} = 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_{1} 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_{0} is approximately 6.6 hours and a is 1.37 × 10^{11} cm (Bitner et al. 2007).
We define the Mach number at the inner edge as
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} = 0.02 a.
2.3. Algorithm
The simulations presented in this paper all use the IDEFIX code, which is a finitevolume conservative Godunov gridbased 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 manycore 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^{−6}, 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:
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 ‘onthefly’ 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 inneredge 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)
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.
Runs presented in this paper.
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
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 wavekilling zone close to the inner edge extending from r_{0} 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 onetenth of the orbital time at the inner boundary. For a wavekilling 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_{0}, 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 wavekilling zone with r_{WKZ} = 2r_{0} 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
where F_{R}(v_{ϕ}) is the intercell 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 intercell angular momentum fluxes computed by the Riemann solver in the integration loop rather than the reconstructed cellcentred values. This allows us to better capture the smallscale 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^{−2} 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 timeaveraged. 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^{−2}. 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^{−2} 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^{−3} 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 higherresolution 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 highresolution run as in the lowresolution 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^{−2} after 100 binary orbits, a value consistent with the lowerresolution runs, confirming that these low values of α are not an artefact of resolution.
4. Discussion
4.1. Onearmed 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 wavekilling zone at the inner boundary to dampen the m = 2 reflection. This suggests a less efficient, nonlinear coupling is at play in our simulations that may have been missed in lowerresolution 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 Σ_{0} can be quantified as
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_{0})≃1 at the inner edge. As σ(r)∝r^{2}, nonlinear 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 onearmed 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^{−4} and are still decaying when we end the simulation. The one exception is for q = 0.1 where transport stabilises to a larger α ≈ 10^{−3}, 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 (OtulakowskaHypka 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 (OtulakowskaHypka 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 MRIdriven transport in outburst. Our initial conditions assume an isothermal disc with a flat surface density profile. A constanttemperature 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 spiraldriven 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^{−4} 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 overstability (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 highresolution 3D hydrodynamical simulations, Manger et al. (2020) propose a scaling for the VSI α ≲ (H/R)^{2.6}, which implies α < 10^{−5} 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 spiralheated 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.

The dominant spiral mode is a onearmed 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.

Spiralwavedriven 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 Xray binaries (whose outer discs sample similar densities and temperatures).
Stable version and documentation can be found at https://github.com/idefixcode/idefix
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.univgrenoblealpes.fr), which is supported by the Grenoble research communities, as well as the JeanZay supercomputer at IDRIS http://www.idris.fr/jeanzay/. This work was granted access to the HPC resources of IDRIS under the allocation 2021A0120402231 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
 Balbus, S. A., & Hawley, J. F. 1991, ApJ, 376, 214 [Google Scholar]
 Balman, Ş., & Revnivtsev, M. 2012, A&A, 546, A112 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Bitner, M. A., Robinson, E. L., & Behr, B. B. 2007, ApJ, 662, 564 [NASA ADS] [CrossRef] [Google Scholar]
 Blondin, J. M. 2000, New Astron., 5, 53 [NASA ADS] [CrossRef] [Google Scholar]
 Chandrasekhar, S. 1961, Hydrodynamic and Hydromagnetic Stability (Oxford: Clarendon) [Google Scholar]
 Edwards, H. C., Trott, C. R., & Sunderland, D. 2014, J. Parallel Distributed Comput., 74, 3202 [CrossRef] [Google Scholar]
 Gammie, C. F., & Menou, K. 1998, ApJ, 492, L75 [NASA ADS] [CrossRef] [Google Scholar]
 Hameury, J. M. 2020, Adv. Space Res., 66, 1004 [Google Scholar]
 Harris, C. R., Millman, K. J., van der Walt, S. J., et al. 2020, Nature, 585, 357 [Google Scholar]
 Harten, A., Lax, P. D., & Van Leer, B. 1983, SIAM Rev., 25, 35 [CrossRef] [Google Scholar]
 Heemskerk, M. H. M. 1994, A&A, 288, 807 [NASA ADS] [Google Scholar]
 Hennebelle, P., Lesur, G., & Fromang, S. 2016, A&A, 590, A22 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Hirose, S., Blaes, O., Krolik, J. H., Coleman, M. S. B., & Sano, T. 2014, ApJ, 787, 1 [NASA ADS] [CrossRef] [Google Scholar]
 Hunter, J. D. 2007, Comput. Sci. Eng., 9, 90 [Google Scholar]
 Ju, W., Stone, J. M., & Zhu, Z. 2016, ApJ, 823, 81 [NASA ADS] [CrossRef] [Google Scholar]
 Kato, T. 2015, PASJ, 67, 108 [NASA ADS] [CrossRef] [Google Scholar]
 Klahr, H., & Hubbard, A. 2014, ApJ, 788, 21 [Google Scholar]
 Kley, W., Papaloizou, J. C. B., & Ogilvie, G. I. 2008, A&A, 487, 671 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Kornet, K., & Rozyczka, M. 2000, Acta Astron., 50, 163 [NASA ADS] [Google Scholar]
 Lasota, J.P. 2001, New A Rev., 45, 449 [NASA ADS] [CrossRef] [Google Scholar]
 Latter, H. N., & Papaloizou, J. C. B. 2012, MNRAS, 426, 1107 [NASA ADS] [CrossRef] [Google Scholar]
 Lesur, G., Hennebelle, P., & Fromang, S. 2015, A&A, 582, L9 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Lesur, G., Ercolano, B., Flock, M., et al. 2022, ArXiv eprints [arXiv:2203.09821] [Google Scholar]
 Lesur, G. R. J., Baghdadi, S., WafflardFernandez, G., et al. 2023, Astrophysics Source Code Library [record ascl:2306.036] [Google Scholar]
 Lovelace, R. V. E., Li, H., Colgate, S. A., & Nelson, A. F. 1999, ApJ, 513, 805 [NASA ADS] [CrossRef] [Google Scholar]
 Lubow, S. H. 1991, ApJ, 381, 259 [Google Scholar]
 Makita, M., Miyawaki, K., & Matsuda, T. 2000, MNRAS, 316, 906 [NASA ADS] [CrossRef] [Google Scholar]
 Manger, N., Klahr, H., Kley, W., & Flock, M. 2020, MNRAS, 499, 1841 [NASA ADS] [CrossRef] [Google Scholar]
 Martin, R. G., Nixon, C. J., Pringle, J. E., & Livio, M. 2019, New Astron., 70, 7 [Google Scholar]
 Masset, F. 2000, A&AS, 141, 165 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Matsuda, T., Sekino, N., Shima, E., Sawada, K., & Spruit, H. 1990, A&A, 235, 211 [NASA ADS] [Google Scholar]
 Meyer, F., & MeyerHofmeister, E. 1984, A&A, 132, 143 [NASA ADS] [Google Scholar]
 Mignone, A., Bodo, G., Massaglia, S., et al. 2007, ApJS, 170, 228 [Google Scholar]
 Mignone, A., Flock, M., Stute, M., Kolb, S. M., & Muscianisi, G. 2012, A&A, 545, A152 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Mineshige, S., & Osaki, Y. 1983, PASJ, 35, 377 [NASA ADS] [Google Scholar]
 Nelson, R. P., Gressel, O., & Umurhan, O. M. 2013, MNRAS, 435, 2610 [Google Scholar]
 Osaki, Y. 1974, PASJ, 26, 429 [NASA ADS] [Google Scholar]
 Osaki, Y. 1996, PASP, 108, 39 [NASA ADS] [CrossRef] [Google Scholar]
 OtulakowskaHypka, M., Olech, A., & Patterson, J. 2016, MNRAS, 460, 2526 [CrossRef] [Google Scholar]
 Paczynski, B. 1977, ApJ, 216, 822 [CrossRef] [Google Scholar]
 Pala, A. F., Gänsicke, B. T., Marsh, T. R., et al. 2019, MNRAS, 483, 1080 [NASA ADS] [CrossRef] [Google Scholar]
 Pjanka, P., & Stone, J. M. 2020, ApJ, 904, 90 [NASA ADS] [CrossRef] [Google Scholar]
 Reback, J., Jbrockmendel, McKinney, W., et al. 2022, https://doi.org/10.5281/zenodo.6408044 [Google Scholar]
 Rozyczka, M., & Spruit, H. C. 1993, ApJ, 417, 677 [CrossRef] [Google Scholar]
 RuizCarmona, R., Khangale, Z. N., Woudt, P. A., & Groot, P. J. 2020, MNRAS, 491, 344 [NASA ADS] [CrossRef] [Google Scholar]
 Savonije, G. J., & Papaloizou, J. C. B. 1983, MNRAS, 203, 581 [Google Scholar]
 Savonije, G. J., Papaloizou, J. C. B., & Lin, D. N. C. 1994, MNRAS, 268, 13 [NASA ADS] [CrossRef] [Google Scholar]
 Scepi, N., Lesur, G., Dubus, G., & Flock, M. 2018, A&A, 609, A77 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Shakura, N. I., & Sunyaev, R. A. 1973, A&A, 24, 337 [NASA ADS] [Google Scholar]
 Smak, J. 1971, Acta Astron., 21, 15 [NASA ADS] [Google Scholar]
 Smak, J. 1984, Acta Astron., 34, 317 [NASA ADS] [Google Scholar]
 Spruit, H. C. 1987, A&A, 184, 173 [NASA ADS] [Google Scholar]
 Stehle, R. 1999, MNRAS, 304, 687 [NASA ADS] [CrossRef] [Google Scholar]
 Toro, E. F., Spruce, M., & Speares, W. 1994, Shock Waves, 4, 25 [Google Scholar]
 Velikhov, E. P. 1959, Zhur. Eksptl’. i Teoret. Fiz., 36 [Google Scholar]
 Virtanen, P., Gommers, R., Oliphant, T. E., et al. 2020, Nat. Methods, 17, 261 [Google Scholar]
 Xu, W., & Goodman, J. 2018, MNRAS, 480, 4327 [Google Scholar]
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 ‘lowresolution’ 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 wavekilling 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 nonisothermal 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). 
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
where the accretion rate Ṁ is taken from the intercell 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^{−2} 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,
Figure D.1 shows the spiral pattern on a deprojected plot. Consistently with Fig. 2, we see that the onearmed 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: Deprojected density fluctuation , as defined in Eq. (8). 
Appendix E: RWI/KHIdriven 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: Deprojected 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 
Current usage metrics show cumulative count of Article Views (fulltext article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.
Data correspond to usage on the plateform after 2015. The current usage metrics is available 4896 hours after online publication and is updated daily on week days.
Initial download of the metrics may take a while.