State-of-the-art simulations of line-driven accretion disc winds: realistic radiation-hydrodynamics leads to weaker outflows

Disc winds are a common feature in accreting astrophysical systems on all scales. In active galactic nuclei (AGN) and accreting white dwarfs (AWDs), specifically, radiation pressure mediated by spectral lines is a promising mechanism for driving these outflows. Previous hydrodynamical simulations have largely supported this idea, but relied on highly approximate treatments of ionization and radiative transfer. Given the sensitivity of line driving to the ionization state and radiation field in the outflow, here we present a new method for carrying out 2.5D radiation-hydrodynamic simulations that takes full account of the frequency-dependent radiative transfer through the wind, the corresponding ionization state and the resulting radiative accelerations. Applying our method to AWDs, we find that it is much harder to drive a powerful line-driven outflow when the interaction between matter and radiation is treated self-consistently. This conclusion is robust to changes in the adopted system parameters. The fundamental difficulty is that discs luminous enough to drive such a wind are also hot enough to over-ionize it. As a result, the mass-loss rates in our simulations are much lower than those found in earlier, more approximate calculations. We also show that the ultraviolet spectra produced by our simulations do not match those observed in AWDs. We conclude that, unless the over-ionization problem can be mitigated (e.g. by sub-grid clumping or a softer-than-expected radiation field), line driving may not be a promising mechanism for powering the outflows from AWDs. These conclusions are likely to have significant implications for disc winds in AGN also.


INTRODUCTION
Signatures of outflowing gas are present in a wide range of astrophysical settings including young stellar objects (YSOs) (e.g.Güdel et al. 2018;Moscadelli et al. 2019), hot stars (e.g.Crowther et al. 2006;Haucke et al. 2018), accreting white dwarfs (AWDs) (e.g.Froning et al. 2012), X-ray binaries (XRBs) (e.g.Tetarenko et al. 2018;Díaz Trigo & Boirin 2016) and quasars/active galactic nuclei (AGN) (e.g.Gibson et al. 2009).In some cases, the driving mechanism is clear.For example, the stellar winds in hot OB-stars (10, 000 K ≲  eff,OB ≲ 50, 000 K) are almost certainly driven by the radiation pressure that is produced when UV photons are scattered by resonance lines in the partially ionized stellar atmosphere.The theory of these 'line-driven' winds has been developed over the course of several decades, beginning with the seminal paper by Castor et al. (1975, hereafter referred to as CAK).
Line-driving has also been suggested as a possible mechanism ★ E-mail: C.Knigge@soton.ac.uk for powering the disc winds in AWDs and AGN, since the disc temperatures in these systems are comparable to  eff,OB .In XRBs, line-driving is unlikely to be dominant, because the ionization state of the outflows seen in these system is too high. 1 Here, either thermal (where gas is heated to the point where its thermal velocity exceeds the local escape velocity) or magnetic driving are more promising mechanisms (Miller et al. 2016;Higginbottom et al. 2017;Tomaru et al. 2020).
Radiation-driven winds arise in the environments of massive luminous objects whenever the outward force associated with the scattering and/or absorption of photons exceeds the inward pull of gravity (Eddington 1916).The famous Eddington limit is associated with the simplest case of a fully ionized, hydrogen-only plasma, in which case the radiative force is solely due to electron scattering.Line driving is an extension of this mechanism to partially ionized media, where interactions between photons and bound electrons can dramatically increase the net radiative force.The ratio of this radiative force relative to the electron scattering case is often referred to as the force multiplier, M, a nomenclature first introduced by CAK.Force multipliers can reach values as high as M ≃ 2000 if the ionization state of the gas is optimal for line driving (Gayley 1995).
Hydrodynamic simulations of line-driven disc winds have been carried out both for AGN-scale systems (Proga et al. 2000;Proga & Kallman 2004, PK;Nomura et al. 2016;Nomura & Ohsuga 2017;Nomura et al. 2020) and for AWDs (Proga et al. 1998, hereafter PSD;Pereyra et al. 1997Pereyra et al. , 2000;;Dyda & Proga 2018a;Dyda & Proga 2018b, hereafter DP).A crucial aspect of all such simulations is the spatial distribution of the force multiplier throughout the outflow.The force multiplier depends directly on the ionization state of the wind, which is expected to be set by radiative processes.Of course, the radiation field itself is also modified as a result of its interactions with the outflowing material, so matter and radiation are intricately and non-linearly coupled in line-driven flows.
In previous hydrodynamic simulations of disc winds, these complex interactions between matter and radiation have been simplified significantly.Usually, this has been achieved by considering an isothermal outflow whose ionization state is either completely fixed (e.g.PSD; Pereyra et al. 1997Pereyra et al. , 2000;;Dyda & Proga 2018a; DP; all in the AWD setting) or estimated via an ionization parameter, , assuming a radiation field modified only by grey attenuation (e.g.Proga et al. 2000;Proga & Kallman 2004, PK;Nomura et al. 2016;Nomura & Ohsuga 2017;Nomura et al. 2020; all in the AGN setting).With these assumptions, the force multiplier is then estimated via the standard CAK - parameterization, modulo a slight modification first introduced by (Owocki et al. 1988) designed to limit the maximum force muliplier to M max (c.f.Equation 4 in PSD).Here,  is the so-called "optical depth parameter", which depends on the local velocity gradient of the flow.The powerlaw index  measures the relative importance of optically thin and optically thick lines.Finally, the normalization constant  encodes the overall efficiency of line-driving.For their AWD model, PSD adopted  = 0.2 and  = 0.6, based upon typical values for winds in OB stars (Gayley 1995).For their AGN simulations, Proga et al. (2000) retain  = 0.6, but use an analytical approximation for  (), based on the calculations of Stevens & Kallman (1990).The maximum force multiplier is typically set to M max ≃ 4400.
All of these simulations generated line-driven winds that were broadly capable of accounting for the observed outflow signatures in these systems, with typical mass-loss rates on the order of  wind ∼ 10 −3  acc ).That said, Drew & Proga (2000) and Proga et al. (2002) already noted some tension between simulations and observations of disc winds in AWDs, with simulations requiring higher-than-expected accretion rates (by a factor of 2-3) to generate enough wind mass-loss and account for the observed strengths of the wind-formed lines.In any case, it is not certain that the assumptions and approximations that had to be made at the time are entirely benign.For example, Sim et al. (2010) and Higginbottom et al. (2014) carried out detailed radiative transfer calculations for the AGN simulation presented by PK and found that line-driving failed in the simulated outflow, because multi-dimensional scatter-ing effects caused the wind to become over-ionized2 .Similarly, the prescription for the force multiplier adopted in these simulations implicitly assumed that -just as in OB-stars -the ionization state in AWDs is near optimal for line driving.Moreover, in generating synthetic spectral lines for the simulated outflows, the ionic species associated with these relevant transitions were assumed to be dominant.Given all this, it has been a long-standing ambition to carry out new radiation-hydrodynamics (RHD) simulations of line-driven disc winds in which multi-frequency radiative transfer and ionization are treated in detail.However, the computational challenges posed by this task are significant.
Here we describe our first attempt to overcome the above challenges and conduct Monte Carlo radiation-hydrodynamics simulations of a disc wind system in which the radiative transfer and force multiplier are calculated realistically and self-consistently.We focus on AWDs as our initial application because, relative to AGN, they can be modelled with a smaller spatial dynamic range and a simpler spectral energy distribution (SED).The remainder of this paper is organized as follows: In Section 2, we present the method we have developed to carry out 2.5D, multi-frequency RHD simulations with a full treatment of ionization and line-driving forces.In Section 3, we present the results of several simulations, including some that can be directly compared to those presented by PSD.Finally, in Section 4, we discuss our results, explore the limitations of our calculations and summarize our conclusions.

METHOD
The fundamentals of our coupled Monte Carlo radiative transfer and hydrodynamics method are as outlined by Higginbottom et al. (2020).Briefly, we use an operator-splitting approach in which the hydrodynamic equations, are solved using the publicly available Godunov-style code pluto 4.4 (Mignone et al. 2007).Here , v,  and Φ are respectively the density, the velocity vector, the thermal pressure and the gravitational potential.The radiative force, g rad , is calculated by using our Monte Carlo ionization and radiative transfer code python 3 (Long & Knigge 2002;Higginbottom et al. 2013;Matthews et al. 2015).
The two aspects of our computational approach -hydrodynamics and ionization/radiative transfer -are described in more detail in §2.2 and §2.3, respectively.Note that we use an isothermal equation of state in pluto, so that we do not have to solve an energy equation.
Our overall iterative numerical procedure can be summarized as follows.After an initialization stage, which we describe below in §2.1, we start by running the hydrodynamic calculation for a short amount of simulation time, Δ RAD (see §2.3), to evolve the density and velocity in the wind.Once this is done, we pass the new velocity, density and temperature fields to python to run a Monte Carlo radiation transport simulation and calculate the ionization state of the gas and the direction-dependent flux in each cell.We then use a stand-alone code (see §2.4) to calculate force multiplier tables using the new ionization state and direction-dependent fluxes coming from python.These force multiplier tables are passed to pluto along with the direction-dependent flux in order to compute g rad .To close the loop, we run the hydrodynamics simulation for another Δ RAD .
We repeat this entire process many times, until the wind reaches a quasi-steady state.All of our simulations employ pluto's standard hydrodynamics module, with linear reconstruction, and use the Harten, Lax, Van Leer (HLL) approximate Riemann Solver.We use 2 nd order Runge Kutta time-stepping and adopt a Courant-Friedrichs-Levy number of 0.4 in our calculations.All of our python calculations utilize the socalled "simple atom" approach described by Long & Knigge (2002), and we adopt solar abundances throughout.In the rest of this section, we provide more details on each step of these RHD simulations.

System Parameters and Initial Conditions
In all of our models, the accretor is taken to be a WD with mass  WD = 0.6  ⊙ and radius  WD = 8.7 × 10 8 cm.We employ a spherical polar coordinate systems and cover one quadrant of the environment close to the WD for our hydrodynamical calculations.Thus our simulation domain extends from  min = 0 to  max = /2 in polar angle and from  min =  WD to  max ≃ 10  WD radially.We use a stretched grid with fixed zone size ratios (as in PSD), with  +1 /  = 1.07 and  +1 /  = 0.95, to ensure that we achieve the finest resolution in the wind launching region, i.e. close to the accretor and close to the disc plane.The number of grid points in the  and  directions are   = 128 and   = 96, respectively.
For our fiducial model, we focus on the simplest astrophysically interesting case: a system containing a luminous accretion disc in which there is no significant radiation associated with the central object.We adopt an accretion rate of  acc =  × 10 −8 M ⊙ yr −1 for this model, which is very much at the high end of the plausible range for "normal" wind-driving AWDs, such as nova-like variables (Howell & Mason 2018).This choice is intended to maximize the chance of generating a robust outflow; it also allows a direct comparison to identical models in PSD and DP.
We also present results for several additional models, characterized by (i) the omission of radiative transfer and ionization effects; (ii) the inclusion of a luminous boundary layer (BL); (iii) a 10× higher accretion rate; (iv) a 3× and 10× lower accretion rate; (v) a flat radial disc temperature profile.Information regarding all of these models is provided in Table 1.
We assume a geometrically thin disc, so all of the disc radiation emanates from the system mid-plane at  = /2.Since the disc represents the mass reservoir for the outflow, we fix the density at  = /2 and hold that value constant throughout the simulation.This density,  0 , should ideally be set to a value that corresponds to the upper parts of the disc atmosphere, just below the wind launching region.For our simulations, it needs to be set high enough to ensure that the critical and sonic points of the outflow are inside the simulation domain, but low enough to prevent the hydrostatic regions within the domain from becoming completely optically thick.We typically set  0 = 10 −9 g cm −3 to achieve this goal.We checked that multiplying and dividing  0 by a factor of 2 did not change the results of our simulations significantly.Hence, it seems that, as long as  0 is set to properly capture the sonic point while avoiding a highly optically thick base of the wind, our results are not very sensitive to the choice of  0 .
The initial configuration of the simulation in pluto corresponds to a hydrostatic disc atmosphere.Thus the velocities in each grid cell are set to the Keplerian velocity at a given distance from the central source, while the densities are initialized to Here,  iso = √︁ /  is the isothermal sound speed and  = 0.6 is the mean molecular weight.We currently keep the temperature fixed at  = 40, 000 K throughout the simulation (so that  iso = 24 km s −1 ).We have verified that this choice of wind temperature does not have a significant effect on our results (see §2.3 and §4).
To initialize the radiative force, we set up a hydrostatic disc atmosphere with a very low density of 10 −20 g cm −3 at its base.This ensures that the entire domain is optically thin.We then use python to calculate the corresponding ionization state and direction-dependent fluxes and follow the procedure described at the beginning of §2 to initialize the radiation force in pluto for the first time step.

Hydrodynamics
We have implemented radiative driving in the hydrodynamics by making use of pluto's built-in support for "body forces".These are user-defined external forces that act on the gas in each grid cell, through a g rad term on the right-hand side of the momentum conservation equation.More specifically, during each hydrodynamic timestep, Δ HD , the net radiative acceleration in each cell is estimated as a the vector sum, over  n = 36 directions that uniformly sample all possible vectors in the  = 0 plane.Here,   and ì  UV, are, respectively, the radiative acceleration and UV flux in direction .Note that we performed a convergence test by doubling the directional resolution to  n = 72 and find no significant difference in our final results.Indeed, the character of the wind is nearly identical with a change in mass-loss rates within 1% between the two simulations.
The quantity M (  ) in Equation 5 is the force multiplier, which depends on the ionization state of the material in the cell (see Section 2.4).In the present context, the key point is that the force multiplier is a function of the dimensionless optical depth parameter,   , which is given by In this expression, d  represents an infinitesimal step along the direction n , so that |d(ì  • n )/d  | is the gradient of the projected velocity component along this direction.Thus M (  ) is both direction dependent and sensitive to the instantaneous local velocity field.These superficially simple expressions capture the complex nonlinear coupling between matter and radiation in line-driven flows: the radiation field and ionization state depend on the outflow dynamics (which control the density and velocity fields), but the dynamics also depend on the radiation field and ionization state (which control  UV, and M).In our approach, pluto is responsible for updating  and ì , while python is responsible for updating (in each cell)  UV, and generating a look-up table for M (  ) across a wide range of   .
In practice, we do not run a python simulation after each hydrodynamic time-step, but only after a fixed interval, Δ RAD ≫ Δ HD , of simulation time.This amounts to the assumption that the ionization state of, and radiation field in, each cell are roughly constant over timescales of Δ RAD .Typically, The characteristic velocity and mass-loss rate may be underestimated for this model, because much of the mass loss takes place in the outer disc, near the edge of the computational domain.
Table 1.Parameters adopted in the simulations and some derived quantities.For each simulation, we provide the model designation (column "Model"), a brief description of its key features ("Comments"), the accretion rate ("  acc "), a brief description of the radial disc temperature distribution ("T d,visc (R)"), a brief description of the prescription used for M () ("Force Multiplier"), the ratio of boundary layer to disc luminosity ("L BL /L disc "), the boundary layer temperature ("T BL "), the Eddington fraction Γ ("L tot /L Edd "), the equivalent models in PSD and DP ("PSD" and "DP"), the wind mass-loss rate ("  wind ") and the velocity velocity of the fast parts of the wind ("  ").
our simulations; we have checked that our results are not sensitive to Δ RAD in this regime.
After each call to python, we update ì  UV, in each cell (see Section 2.3) and also construct a new look-up table for M (  ) for the cell, based on its current ionization state (see Section 2.4).However,  and ì  in Equation 6 are updated during each hydrodynamic time-step, so the line optical depth parameters, force multipliers and radiative accelerations do evolve on this time-scale as well.
As noted above, the density at the mid-plane ( = /2) is kept fixed and acts as a mass reservoir.We also impose a a density floor of 10 −20 g cm −3 throughout the simulation domain.The mid-plane and inner radial boundary conditions (BCs) are both reflective, whilst the  = 0 BC is axisymmetric.An outflow BC is used for the outer radial edge of the grid.

Ionization and Radiative Transfer
Each time python is called after Δ RAD of simulation time, it reads in the latest snapshot of the density and velocity fields calculated by pluto, as well as the wind ionization state produced by the previous call to python.Photon packets are then injected into the outflow by the accretion disc and the wind itself.These photon packets are tracked as they make their way through the wind, where they can undergo attenuation due to bound-free and free-free opacity, as well as electron scattering and bound-bound interactions.
The photon packets passing through a cell are used to construct estimators that describe the radiation field inside the cell.These are used to update the ionization state of the wind and to calculate the radiative accelerations (see below).
Unless otherwise indicated, we assume the disc initially has an effective temperature distribution that follows the standard Shakura-Sunyaev (1973) form, where  represents the cylindrical radius and Here,  acc is the accretion rate through the disc.For the purpose of generating photon packets, the disc is split into concentric annuli, each of which is initially taken to radiate as a blackbody with effective temperature  ,eff =  ,visc ().This then fixes both the number and the frequency distribution of the photon packets generated by a given annulus.
The plasma in the simulation domain itself also produces photons (via free-free, free-bound and bound-bound processes).Since the temperature of this material is kept constant in space and time, the simulation does not enforce strict energy conservation: the gas in each cell can, in principle, emit more or less energy than it absorbs.This is not a serious problem in practice, since typically the plasma is not a significant energy source in our simulations.We have nevertheless carried out two tests to check the validity of this approximation.First, as discussed in more detail in Section 4, we have run a much more detailed python simulation for the final state of our fiducial wind model, in which the thermal state of the wind is iterated to convergence along with the ionization state.The resulting temperatures are comparable to the  = 40, 000 K we assume in our RHD simulations.Second, we have repeated our fiducial RHD simulation with different gas temperatures ( = 30, 000 K and  = 50, 000 K).The differences between these runs are relatively minor and do not affect our main results and conclusions.
As photons propagate through the wind, they are used to update two sets of estimators related to the radiation field: (i) the angleaveraged mean intensity,   , and (ii) the direction-dependent UV flux, ì  UV, .As before, the index  here denotes one of  n = 36 directions chosen to fairly sample all possible vectors in the  = 0 plane.The UV band is defined as running from 7.4 × 10 14 Hz to 3 × 10 16 Hz.The flux estimator, ì  UV, , is used by pluto to calculate the radiative accelerations (c.f.Equation 5), while the mean intensity estimator,   , is used to calculate the ionization state of the gas and the corresponding relationship between the force multiplier and the optical depth parameter (see Section 2.4 below).
Once the photon transport phase is concluded, python calculates the ionization state of the gas based on the local density, temperature and frequency-dependent intensity,   .The construction of the   estimator implicitly takes account of any attenuation of photon packets in the wind, as well as any re-emission from the wind itself.Crucially, it also takes account of photons arriving in a cell from other parts of the wind, including via scattering events.This scattered and/or reprocessed component of the radiation field makes it much harder for the inner flow to successfully shield material further out.In the context of line-driven disc winds in AGN, it has already been demonstrated that including this component is critically important (Sim et al. 2010;Higginbottom et al. 2014). 4n each call to python, we carry out at least two iterations of the entire photon generation, radiation transport and ionization equilibrium calculation.This allows us to deal with the effect of irradiation on the disc temperature distribution.Some photon packets generated by the central source / boundary layer -as well as some photon packets that have undergone reprocessing in the outflow -will inevitably hit the disc surface.In python, we can treat such photons in one of three ways: (i) we can discard them; (ii) we can assume that they are absorbed and reprocessed; (iii) we can assume that they are specularly reflected.In the simulations described here, we adopt option (ii).Thus, after the first iteration, we update the effective temperature distribution across the disc such that  4 ,eff () =  4 ,visc () +  irr (), where  irr () is the irradiating flux incident on the disc at radius .

Creating Look-Up Tables for the Force Multiplier
After each call to python, we follow the method outlined by Parkin & Sim (2013) to create updated look-up tables for the local force multiplier, based on the SED and ionization parameter in each cell.We do not use python directly for this purpose because the calculation of accurate force multipliers requires a larger line list than the one used in python (which is sufficient for obtaining the accurate estimate of   we need to calculate the ionization state).Instead, we pass the ionization state and the local estimate of   to a stand-alone code that has access to over 450,000 lines (see Parkin & Sim 2013, for details).For each transition, the code first calculates the quantity where  and  refer to the lower and upper level of the relevant transition, respectively.Similarly,  , and  , are the usual Einstein coefficients for absorption and stimulated emission, respectively.The number densities of ions in these levels,   and   , are calculated using the ion densities supplied by python and assuming that the excited state populations follow a Boltzmann distribution (which is equivalent to assuming that the excited levels are in local thermodynamic equilibrium).The force multiplier is then calculated from a weighted sum of contributions from all the lines.In principle, the weighting factor of each transition depends on the specific flux  , at the wavelength of that transition (see e.g.Parkin & Sim 2013).However, to avoid the computational expense of finding  , for each direction bin in every grid cell, we make the assumption that this weighting factor can be replaced with its angle-averaged value.
That is, we approximate where Δ  is the Doppler width of the line, given by where  0 is the central frequency of the line.Figure 1 shows the form of M () for three representative cells in different regions of the simulation, with the form of M () adopted by PSD and DP shown for comparison.While the shape of the curves is rather similar, the figure shows that the self-consistent calculation of M () leads to variation of the form of M () and in general a lower normalisation compared to the idealised CAK approximation.As we shall see in the next section, these fairly modest differences in M () lead to real and substantive changes to the character and strength of the resulting outflow in a self-consistent simulation.

Fiducial Model
Figure 2 shows a snapshot of the density distribution for our fiducial model (HK22D), while Figure 3 shows the mass-loss rate through the outer boundary of this model as a function of simulation time.
The density distribution in Figure 2 corresponds to  = 821 s of simulation time.This is ≃1/5 of the sound crossing time-scale, but Figure 3 shows that the flow has converged to a quasi-steady state at this point.Figure 4 provides details on the angular distribution of outflow properties in the fiducial model, taken from the same snapshot as Figure 2. Roughly speaking, the main parts of the outflow are concentrated within a 10°wedge at polar angles between 50°and 60°.
The most important feature of the simulation is the total massloss rate,  wind ≃ 5 × 10 −14 M ⊙ yr −1 , i.e.  wind < 10 −6  acc .The mass loss rate is roughly 100× lower than the value seen in the corresponding simulations reported by PSD and DP, which did not include any ionization and radiative transfer calculations.The lower mass-loss rate we find is associated with both a lower density in the outflowing wind and a reduced velocity.As discussed in more detail below, this is not just a theoretical issue.As a direct consequence of the reduced mass-loss rate in our simulation, it also fails to reproduce the wind signatures that are commonly observed in AWDs.
Why is the outflow in our simulation so much weaker and slower than in earlier calculations?There are two major contributing factors: First, even relatively luminous AWDs -such as the system described by our fiducial model -are at best marginally capable of powering a line-driven wind (Proga 2005).For the physical conditions typically found in such outflows, the maximum force multiplier that can be achieved is M max ≃ 2000 (Gayley 1995).In order to drive strong line-driven disc winds, accreting systems thus require Eddington ratios of at least Γ = / Edd ≃ 5 × 10 −4 .As shown in Table 1, our fiducial model is characterised by Γ = 10 −3 , only just above the threshold for efficient line-driving, even under optimal ionization conditions.As a result, successful line-driving in AWDs requires near-optimal conditions.Against this background, the second key factor is that the ionization conditions are not near-optimal, contrary to what had previously been assumed.PSD and DP adopted the standard CAK approximation for the force multiplier, M () =  −  , with parameters based on models of line-driven winds from hot stars, including a hard upper limit of M max = 4400.With this approximation, M () ≃ M max throughout much of the wind volume, and M ≳ 1000 in the acceleration zone (see below and Figure 5).This allows a strong wind to be driven from the disc.By contrast, our simulations show that the physical conditions in the outflow are actually not all that conducive to line driving, so the system struggles to launch a fast and powerful wind.
The underlying problem is demonstrated quantitatively in Figure 6.In hot stars, the single element that contributes most to the linedriving force is Iron (e.g.Vink et al. 2001;Noebauer & Sim 2015).There, Iron is mainly found in Fe iv and Fe v, and these ionization stages also dominate the bound-bound opacities.By contrast, the upper left panel of Figure 6 shows the Iron ionization state throughout our fiducial model.The characteristic ionization level is much higher than in hot stars, and the ionization stages required for efficient line driving are largely absent.
But why is the ionization state so high?In hindsight, this should not come as a surprise, because the maximum temperature in the disc is  ,max ≃ 74, 000 K. Figure 7 shows the mean intensity in three cells that lie in different regimes of the outflow (the location of these cells is highlighted in Figure 6).There is significant energy in the SED above the Fe v ionization threshold, and so the gas is efficiently ionized beyond this stage.
In order to verify this interpretation, we have tried to replicate the results of PSD and DP by adopting their CAK approximation for the force multiplier, with  = 0.59 and  = 0.6. 5This simulationmodel HK22Df in Table 1 -is much faster computationally, because there is no ionization step and radiative transfer. 6Figure 5 allows a comparison of the resulting distribution of the force multiplier to that found in our fiducial model.As expected, M is much higher when the CAK approximation is adopted.This over-estimate of the driving force artificially increases the mass-loss rate and wind velocity, bringing them in line with the higher values previously found by PSD and DP.All of this demonstrates how important it is to take full account of ionization and radiation transfer in dynamical simulations of line-driven winds.
A major benefit of our RHD approach is that python is also capable of producing high quality synthetic spectra for any given simulation.Figure 8 shows a set of UV spectra generated from the snapshot of our fiducial model.There are wind-formed spectral features, but they are all associated with high-ionization transitions.Observationally, the strongest UV lines in wind-driving AWDs tend to be N v 1240 Å, Si iv 1400 Å and C iv 1550 Å.Of these, only N v 1240Å is associated with a reasonably strong feature in the synthetic spectra.Thus, as expected based on the very low mass-loss rate, this simulation is not a good match to observations.

Modified Models
Our fiducial model fails to launch a powerful wind because (a) it is barely luminous enough even under optimal conditions, and (b) it is too highly ionized to generate such conditions.Could reasonable changes in the adopted system parameters overcome these challenges?There are two obvious levers.First, we can try to make the system more luminous.The natural ways to do this are to account for the possible contribution of a boundary layer between the WD and the accretion disc or to simply increase the accretion rate.Both of these will, however, produce even harder SEDs, potentially exacerbating the ionization problem.Second, we can try to soften the SED.The only natural way to accomplish this is to reduce the accretion rate, although this will clearly make the system even less luminous.That said, it is well known that disc SEDs are fairly poorly described by sums of blackbodies, and energy dissipation/reprocessing in the upper disc atmosphere and base of the wind might conceivably lead to a softer SED than expected in the Shakura-Sunyaev picture (e.g.Knigge et al. 1998;Shaviv & Wehrse 1991;Hubeny & Long 2021).This idea can be crudely approximated by adopting a flatter-thanstandard effective temperature distribution across the disc.In the rest of this section, we will explore each of these modifications.

Including a Luminous Boundary Layer
In the standard Shakura-Sunyaev picture, the accretion disc releases 50% of the total available accretion luminosity.If the accretor is a non-rotating object with a physical surface (i.e.not a black hole), the rest is radiated away by a narrow boundary layer that connects the 5 These parameters are designed to produce equivalent force multipliers to those in PSD and DP.The numerical values for  are different, however, because we use a different value for the thermal velocity. 6This means that, as in PSD and DP, the UV flux is not attenuated in the wind, in contrast to our fiducial model with radiative transfer.surface to the inner disc.At the accretion rates of interest here, the boundary layer in AWDs is expected to be optically thick and radiate approximately as a blackbody.In order to obtain the largest plausible luminosity boost with the smallest plausible effect on ionization, we optimistically assume that the BL covers the entire surface of the WD.This yields the lowest reasonable effective temperature for a luminous BL, set by  BL =  disc = 4 2 WD  4 BL .For our model parameters, this leads to  BL ≃ 1.1 × 10 5 K.
When we add this boundary layer SED contribution to the simulation (HK22C in Table 1), its negative effect on the ionization state of the outflow completely dominates its positive effect on the luminosity of the system.In fact, there is no organized fast outflow at all in the final state of this simulation.

Increasing the Accretion Rate
Another way to increase the luminosity of the system is to increase the assumed accretion rate.As noted in Section 2.1, the value  acc =  × 10 −8 M ⊙ yr −1 adopted in our fiducial model is actually already a high value for AWDs.We have nevertheless run a simulation -HK22F in Table 1 -in which we increased the accretion rate to  acc =  × 10 −7 M ⊙ yr −1 .We stress that this is not a realistic accretion rate for "normal" AWDs.In fact, it is near or above the threshold where steady nuclear burning would be induced on the WD surface, a situation thought to be found in supersoft X-ray binaries (van den Heuvel et al. 1992).
Ignoring these concerns for the moment, we find that this high- acc model does produce a notably stronger line-driven wind, with  wind ≃ 10 −12 M ⊙ yr −1 .However, this mass-loss rate is still a factor of ≃ 30 less than obtained by PSD.The characteristic velocity in our simulation,  ≃ 1000 km s −1 , is also much lower than that found by PSD,  ≃ 7000 km s −1 .Finally, the outflow remains very highly ionized and produces no strong UV signatures except O vi 1032 Å.

Decreasing the Accretion Rate
As noted above, we can try to prevent the outflow from becoming over-ionized by lowering the accretion rate, since this will soften the SED.For example, in our fiducial model, the maximum effective temperature in the accretion disc is  ,max ≃ 74, 000 K. Reducing the accretion rate by an order of magnitude, to a value of  acc =  × 10 −9 M ⊙ yr −1 , lowers this to  ,max ≃ 41, 000 K. However, this change also reduces the overall luminosity by a factor of 10, bringing the Eddington ratio down to Γ ≃ 10 −4 .
Our RHD simulation for this case is labelled HK22B1 in Table 1.It does not produce an outflow.This should not come as a surprise, since Γ < M −1 max for this model.In fact, the matching outflow in PSD (their Model 1) also failed to produce a supersonic outflow, even though their use of the CAK approximation assumes near-optimal ionization conditions for line-driving.
In order to test if there is a "Goldilocks zone" for line-driving around  acc ≃ 10 −8 M ⊙ yr −1 , we have also run a model for this accretion rate.This is listed as HK22B2 in Table 1 and is characterized by Γ ≃ M −1 max .It, too, does not produce an outflow.

Flattening the disc Temperature Distribution
As a final test, we consider a disc that generates the same luminosity as that in our fiducial model, but with a flat effective temperature distribution.For our system parameters, and accounting for the finite radial extent of our simulation domain ( max = 10  WD ), the implied constant disc temperature is then  ,visc ≃ 40, 000 K. The corresponding simulation is listed as HK22Ds in Table 1.Since the maximum disc temperature is greatly reduced, this model produces a markedly softer SED than our fiducial model, at the same luminosity.Since the ionization state of the wind is now more conducive to line-driving, the mass-loss rate increases by about an order of magnitude to  wind ≃ 4 × 10 −13 M ⊙ yr −1 .However, the characteristic wind velocity in this case is also much lower,  wind ≃ 500 km s −1 .This happens because, with fractionally more luminosity emerging at large disc radii, the wind is now driven preferentially from the outer disc.The lower wind speed therefore primarily reflects the reduced Keplerian and dynamical velocities in this region. 7igure 9 shows the synthetic UV spectra for this model for a range of inclination angles.These spectra do show some of the features seen in observations, notably the N v and C iv resonance features at 1240 Å and 1550 Å, respectively.However, as a result of the low wind velocity in this model, the line profiles are quite narrow and only slightly blue-shifted with respect to their rest wavelengths.There is also no Si iv 1400 Å absorption feature, because the ionization state of the wind is still too high to produce this line.This contrasts with observations of AWDs, in which Si iv 1400 Å is commonly observed (e.g.Hartley et al. 2002).

DISCUSSION AND CONCLUSIONS
We have presented new RHD simulations of line-driven accretion disc winds that include a detailed treatment of ionization and radiative transfer.Focusing on the outflows from AWDs, our main finding is that the physical conditions in these systems are far less conducive to efficient line-driving than implicitly assumed in previous calculations.This is entirely because these outflows turn out to be more highly ionized than those of OB stars.Earlier simulations had adopted approximate force multipliers based on these hot single stars, causing them to systematically overestimate the driving forces.
The overall outflow characteristics are highly sensitive to these issues, because AWDs are only marginally luminous enough to drive a powerful line-driven wind, even under near-optimal ionization conditions.Perhaps most importantly, the mass-loss rates in our simulations are significantly lower than previous estimates, by two orders of magnitude for our fiducial model.In line with this, the corresponding synthetic spectra do not resemble those observed.All of these results raise serious questions about the viability of line-driving as the primary mechanism for generating the winds observed in AWDs.
Fundamentally, our fiducial AWD model struggles to launch a powerful line-driven disc wind for two reasons: its relatively low Eddington ratio (Γ = / Edd ≃ 10 −3 ) and the strongly ionizing radiation field produced by its hot disc ( ,max ≃ 74, 000 K).
We have therefore also explored whether simple modifications to the model might generate a stronger outflow that is more in line with observations.Including a luminous boundary layer in order to provide more driving photons fails as a strategy: it over-ionizes the outflow even more and causes the outflow to stall completely.Increasing the accretion rate does generate a higher mass-loss rate, but the resulting outflow remains highly ionized and fails to produce the observed wind signatures.This type of model is also physically unrealistic, because accretion at such high rates would likely trigger  steady nuclear burning on the WD, making the system a supersoft source, rather than a "normal" wind-driving AWD (i.e. a nova-like cataclysmic variable or a dwarf nova in outburst).Softening the SED by reducing the accretion rate reduces the Eddington ratio even more and again chokes the wind completely.).This may be compared to the same plot for our fiducial model in Figure 8. Spectra are shown for a range of inclinations and normalized such that the flux levels correspond to a system observed at 100 pc.The positions of the Lyman limit and several key UV resonance lines are marked with vertical grey lines.
The only modification that is marginally successful in both increasing the mass-loss rate and producing more of the observed spectral signatures is to adopt a flat effective temperature profile for the accretion disc.This may be viewed as an ad-hoc attempt to describe the putative reprocessing effect of the disc-wind transition region on the SED.This change leads to an order-of-magnitude increase in  wind and the appearance of C iv 1550 Å in the synthetic spectra.However, since this outflow is driven primarily from larger radii, its characteristic speed is significantly lower than implied by observations.The synthetic spectra also presents hardly any wind-formed Si iv 1400 Å, a relatively low-ionization wind signature commonly seen in AWDs.
Overall, our RHD simulations suggest a rather bleak outlook for line-driving in AWDs: systems luminous enough to drive a wind are also hot enough to over-ionize it.Drew & Proga (2000) already highlighted the low Eddington ratios found in AWDs as a serious problem for models of line-driven winds in these systems.As also noted by Drew & Proga (2000), the only obvious loophole to their conclusion was that the empirically inferred accretion rates (and hence luminosities) of AWDs are systematically underestimated.Our results now close this loophole, since even if AWD were accreting at higher rates than thought they would over-ionize the wind, hence reducing the efficiency of the line-driving mechanism, so as to prevent the wind from becoming massive enough to explain observations.Might these conclusions change if we relax one or more of the assumptions and approximations made in our calculations?One such assumption is that the accreting material is characterized by solar abundances.The efficiency of line-driving does depend on composition: increasing the abundance(s) of atomic species that contribute significantly to the net line force -e.g. by increasing the overall metallicity -will tend to increase the strength of the resulting outflow (CAK).However, most "normal" AWDs appear to exhibit roughly solar or slightly sub-solar abundances (e.g.Harrison 2016;Pala et al. 2017).The main exception to this is a sub-set of AWDs in which the accreting material appears to have undergone some degree of CNO processing (Gänsicke et al. 2003).Even in these systems, the line-driving efficiency is unlikely to be significantly different from the solar-abundance case (c.f.Abbott 1982). 8 Perhaps the most serious remaining approximation is that we treat the wind as isothermal.We are not yet able to include heating and cooling in our models self-consistently, but we expect to be able to in the not too distant future.This will be particularly important for modelling line-driven disc winds in AGN.The SED and geometry of the radiation field are more complex in these systems and may give rise to significant temperature gradients in the outflow.
In the meantime, we can take a snapshot of our fiducial model and run python on it until thermal equilibrium is achieved.Figure 10 shows the resulting temperature structure.As expected, the base of the wind heats up to roughly the same temperature as the underlying thin disc, and the gas in the acceleration zone settles at temperatures not too far from the 40, 000 K assumed in the RHD simulations.Since the effect of the wind temperature on the ionization state is relatively modest, we believe that relaxing the isothermal approximation in our RHD simulations of AWDs is unlikely to fundamentally change our current conclusions.
A more promising way to make line-driving more efficient in AWDs is to relax the assumption that the disc radiates as an ensemble of blackbodies that follow the Shakura-Sunyaev effective temperature profile.Among all the modifications we have tried so far, our (extremely crude) first attempt to implement this idea -by imposing a constant effective temperature across the disc -is the only one that was marginally successful.We have also carried out some experiments with models in which the disc is represented by an ensemble of stellar atmospheres, rather than an ensemble of blackbodies.Unfortunately, first indications are that the overall disc SED in this case actually has a more pronounced high energy tail, leading to an even higher ionization state -and thus a weaker wind -than obtained with the blackbody approximation.Nevertheless, if the outflows from AWDs are due to line-driving, a good description of the disc SED is likely to be the most critical ingredient for future RHD simulations.
Another potentially impactful assumption of our work is that the flow remains smooth on all scales smaller than the grid resolution of our simulations.This may be important because if the flow were clumped on small scales this would effectively increase the density of the material which enhances recombination.Thus small-scale inhomogeneity could help to alleviate the overionization problem.This is a promising avenue for further study since clumping is widely thought to operate in line-driven stellar winds (see e.g.Fullerton 2011, and references thererin) and has also been invoked in disc wind models where it can lead to better agreement with observations (see e.g.Matthews et al. 2016).
There is, of course, no guarantee that line-driving is responsible for generating the outflows from AWDs.Magnetic fields, in particular, can also drive accretion disc winds, via centrifugal forces (e.g.Blandford & Payne 1982;Pelletier & Pudritz 1992) and/or magnetic pressure gradients (e.g Uchida & Shibata 1985;Stone & Norman 1994).In fact, there is at least one significant empirical strike against 8 We note that another distinct class of AWDs may also provide an important test for line-driven disc wind models: AM CVn stars.In these ultra-short period systems ( orb ≲ 60 min), the donor is itself a fully or partially degenerate dwarf.As a result, the accreting material is severely Hydrogen deficient and may show other abundance anomalies also.Interestingly, the prototype of this class displays clear blue-shifted absorption and P-Cygni profiles in its ultraviolet spectrum (Wade et al. 2007).line-driving in AWDs: multi-epoch UV observations of two winddriving systems do not exhibit the expected strong correlation between wind activity and continuum level (Hartley et al. 2002).In any case, our simulations clearly show that the nature of the wind-driving mechanism in AWDs is far from settled.
The challenge posed by our simulations for line driving in AWDs is likely to extend to AGN as well.Sim et al. (2010) and Higginbottom et al. (2014) already showed explicitly that the line-driven AGN winds in the simulations carried out by PK would, in reality, be overionized.In these simulations, the ionization state of the wind was assumed to be dominated by a compact luminous X-ray source that is coincident with an accreting 10 8 M ⊙ black hole.Our work here suggests that the same fundamental issue may affect even AGN where the SED is dominated by a relatively cool accretion disc (e.g.X-ray weak quasars with higher black hole masses).The characteristic disc temperatures in such systems are quite similar to those found in luminous AWDs (e.g.Matthews 2016), so achieving high force multipliers is likely to be difficult.At sufficiently high Eddington ratios, modest force multipliers may be enough to power a radiativelydriven disc wind regardless.However, the viability of line driving across the full spectrum of AGN needs to be urgently reassessed with RHD simulations of the kind we have presented here.In order to test and validate the method presented in Section 2, we perform a 2D simulation of a OB, line-driven stellar wind and compare our results with the 1D calculations of (Noebauer & Sim 2015, hereafter NS15).We adopt the same physical parameters as in NS15, i.e. the mass of the central star is  ★ = 52.5⊙ , the luminosity of the star is  ★ = 10 6  ⊙ and its surface effective temperature (which is also set to be the isothermal temperature of the wind) is  eff = 4.2 × 10 4 .The physical domain is extending from 1 to 10  ★ in radius, with  ★ = 1.317 × 10 12 cm, and from 0 to /2 in latitude.We use the same logarithmically stretched grid as in our disc simulations but with a higher resolution in radius of   = 512 points so as to resolve the sonic point and a slightly lower grid resolution of   = 32 to save computational time.We find that the use of a stretched grid in the  direction is actually quite important for preserving spherical symmetry throughout the simulation.The ideal strategy is to use a cell ratio in the  direction that preserves a constant solid angle across the grid.Otherwise, different level of Monte-Carlo noise between cells can introduce strong asymmetries in the simulation.The inner radial boundary condition is set to a hydrostatic density profile for the density and so as to conserve the radial flux of mass for the radial velocity when the flow is outgoing (the radial velocity is set to zero otherwise).The outer radial boundary condition is an outflow boundary condition that prevents radial inflow from the outer edge of the simulation.The latitudinal and azimuthal velocities use outflow conditions at both radial boundaries.The latitudinal boundary condition is outflowing for all quantities at the pole and the midplane.
This stellar wind simulation enables us to test our discretization over direction of the flux (see Section 2.2) in a different geometry.It allows us, in particular, to test if our implementation does not introduce any undesired asymmetry in the flow.In this spirit, we use the same number of directions  n = 36 as in our disc case.Finally, we use Δ RAD = 200Δ HD for the coupling between PYTHON and PLUTO.
Figure A1 shows a snapshot 2D density map from our stellar wind simulation once the mass loss rate has converged to a constant value.The stellar wind is clumpy with density contrasts that can be as large as one order of magnitude.This clumpiness cannot be due to the line-driving instability (Owocki & Rybicki 1984) since we use the Sobolev approximation in python.Determining whether the clumps are due to a hydrodynamic instability or a radiative/ionization instability is out of the scope of this paper.Despite the clumpiness of the wind, the density structure remains spherically symmetric, with the clumps being distributed isotropically, as expected for a stellar wind.Note that the structure looks slightly more homogeneous near the −pole.This is due to the larger cells in the -direction near the pole.
We also plot on Figure A2 the latitudinally averaged profiles of the density and the radial velocity as a function of radius.The solid lines show the results of our 2D simulations while the dashed lines show the results from the 1D simulation of NS15, which despite being 1D include a correction for the finite cone effect and the inclusion of scattering in the wind.It is remarkable that despite the presence of large clumps in our simulation, our density and velocity profiles agree quite well (within a factor of 2 for the density and 1.5 for the velocity profiles) with the 1D simulation of NS15.We believe that the quantitative differences between the 1D and 2D simulations could come from the approximated way the finite cone effect is accounted for in 1D and/or the impact of the 2D clumpy structure of the wind.
To conclude, the conservation of the spherical symmetry of the stellar wind (despite its clumpiness) and the agreement between the density and velocity profiles between our 2D simulations and the 1D simulation of NS15 validates the method we presented in Section 2. The fact that we can reproduce the results of NS15 in hot stellar winds also demonstrates that the main result of our paper, which is that AWD winds are overionized and so cannot reproduce the observed mass-loss rate, is not due to our method but due to improved selfconsistent treatment of the interaction between matter and radiation compared to previous studies.

APPENDIX B: LOGARITHMIC DENSITY MAP
In order to facilitate comparisons between Figure 5

Figure 1 .
Figure 1.The force multiplier M as a function of the dimensionless optical depth parameter  for three representative cells in the simulation domain (the locations of the cells are shown in Figure6by colour coded symbols).For comparison, the dot-dashed green line shows the form of M () for the CAK  −  approximation (equation 1) adopted by PSD and DP (see section 3.1 and footnote 4 for an explanation of the value of  = 0.59).

Figure 2 .
Figure 2. The density and velocity structure of the fiducial model.The density is plotted as a colour-map, with velocity vectors overlaid.The grey lines represent streamlines launched from the  = 86.6°surface.The solid black line shows the location of the Mach 1 surface.The black, red and blue dots show the position of three representative cells in the simulation used for analysis in Figure 6 and Figure 7. See also Figure B1 in the §B for a version of that figure with logarithmic scales for the axis.A line-driven disc wind is produced, but the flow is significantly weaker and slower than in earlier calculations.

Figure 3 .
Figure 3.The variation of the mass-loss rate through the outer boundary (  wind ) as a function of simulation time in the fiducial model.The vertical dashed line shows the time-stamp used to generate the figures illustrating the model.The high density disc region is excluded in calculating  wind .

Figure 4 .
Figure 4. Characteristic physical properties of the fiducial model as a function of polar angle, .In the upper panel, we show the product of radial velocity and density,   , at the outer radial boundary (solid line, left scale) and the cumulative mass-loss rate,  wind (  ′ <  ) = ∫  0

Figure 5 .Figure 6 .
Figure 5.The distribution of the force multiplier, M, throughout the fiducial model (left panel, model HK22D) and through an equivalent model with no explicit treatment of ionization and radiative transfer (right panel, model HK22Df; a fixed k −  formulation is adopted in this model [see text for details]).Since M is direction-dependent, the quantity shown here is the force multiplier in the direction of maximum acceleration.The black lines indicate the Mach = 1 surfaces.

Figure 7 .
Figure 7.The estimator of the mean intensity,   , in three representative cells of the snapshot from our fiducial model.The black, blue and red lines correspond to the corresponding cells highlighted in Figure The thick grey line is a 70,000 K blackbody that is shown for reference.The vertical dashed lines mark the important C iv and Fe v ionization edges.

Figure 8 .Figure 9 .
Figure8.Synthetic UV spectra generated from a snapshot of the fiducial model for a range of inclinations.The spectra are normalized such that the flux levels correspond to a system observed at 100 pc.The positions of the Lyman limit and several key UV resonance lines are marked with vertical grey lines.

Figure 10 .
Figure10.The final wind temperatures obtained for a snapshot of the fiducial model after running the stand-alone version of python to convergence on it (such that both ionization/recombination and heating/cooling are in balance).As usual, the black line marks the Mach = 1 surface.The red contour marks  = 40, 000 K, the constant value adopted in our RHD simulations.

Figure A1 .
Figure A1.2D density map for our stellar wind model.The wind is quite clumpy but we recover a spherically, symmetric wind, validating our 2D implementation of the radiative fluxes and forces presented in §2.
, Figure 6, Figure 10 and the density structure shown in Figure 2, we provide in Figure B1 a map of the density with the axes on a logarithmic scale.This paper has been typeset from a T E X/L A T E X file prepared by the author.

Figure A2 .
Figure A2.Density (top figure) and radial velocity (bottom figure) profiles as a function of radius for our stellar wind model.The solid black line show the (latitudinally-averaged) results of our 2D simulation and the green dashed line show the results from the 1D simulation of NS15 including a finite cone effect and full scattering in the wind.As expected, the results of our 2D simulation are in agreement with the results of NS15.This validates our 2D implementation of the radiative fluxes and forces presented in §2.

Figure B1 .
Figure B1.Density map for our fiducial model with logarithmic axes.The solid black line shows the location of the Mach 1 surface.The black, red and blue dots show the position of three representative cells used for analysis in Figure 6 and Figure 7.