# Non-Axisymmetric Line Driven Disc Winds I - Disc Perturbations

###### Abstract

We study mass outflows driven from accretion discs by radiation pressure due to spectral lines. To investigate non-axisymmetric effects, we use the Athena++ code and develop a new module to account for radiation pressure driving. In 2D, our new simulations are consistent with previous 2D axisymmetric solutions by Proga et al. who used the Zeus 2D code. Specifically, we find that the disc winds are time dependent, characterized by a dense stream confined to relative to the disc midplane and bounded on the polar side by a less dense, fast stream. In 3D we introduce a vertical, -dependent, subsonic velocity perturbation in the disc midplane. The perturbation does not change the overall character of the solution but global outflow properties such as the mass, momentum and kinetic energy fluxes are altered by up to 100%. Non-axisymmetric density structures develop and persist mainly at the base of the wind. They are relatively small, and their densities can be a few times higher that the azimuthal average. The structure of the non-axisymmetric and axisymmetric solutions differ also in other ways. Perhaps most importantly from the observational point of view are the differences in the so called clumping factors, that serve as a proxy for emissivity due to two body processes. In particular, the spatially averaged clumping factor over the entire fast stream, while it is of a comparable value in both solutions, it varies about 10 times faster in the non-axisymmetric case.

###### keywords:

hydrodynamics - radiation: dynamics - methods: numerical - stars: winds, outflows - galaxies: active - X-rays: binaries^{†}

^{†}pagerange: Non-Axisymmetric Line Driven Disc Winds I - Disc Perturbations–LABEL:lastpage

^{†}

^{†}pubyear: 2017

## 1 Introduction

Many accretion disc systems exhibit blue shifted absorption lines, which are interpreted as evidence of outflowing gas. These outflows are likely driven from the disc, raising the question of what physical mechanism is responsible for the driving. One possibility is radiation pressure from spectral lines, due to continuum radiation from the disc scattering off resonance lines. By including line opacities, the amount of momentum transfer is enhanced by several orders of magnitude above pure electron scattering, allowing for outflows to be launched in systems with sub-Eddington luminosity (Lucy & Solomon 1970, hereafter LS70 and Castor, Abbott & Klein 1975, hereafter CAK). Originally, line driving, operating in the Sobolev approximation, was proposed as a driving mechanism for stellar winds in OB stars (LS70) and was successful in predicting their mass loss rates and outflow velocities (Pauldrach, Puls & Kudritzki 1986, Friend & Abbott 1986). Line driving has since been applied to a variety of systems, in particular those with accretion discs. Examples of such systems include high-state cataclysmic variables (CVs), young stellar objects (YSOs) and broad absorption line quasi-stellar objects (BAL QSOs).

Initially theoretical studies of line driven disc winds were limited in comparison to spherically symmetric stellar winds. Some 2D theoretical models introduced simplifying assumptions to render the axisymmetric disc wind problem analytically tracktable, such as simple scaling for the strength of the radiation field (Vitello & Shlosman 1988) or a decoupling of the radial and polar angle equations of motion (Murray et. al 1995).

These early attempts were followed by 2D simulations carried out by Pereyra, Kallman & Blondin (1997) and Pereyra (1997) in the context of CVs. These simulations found steady state solutions but were too coarse to resolve the subsonic part of the flow near the accretion disc. Proga, Stone and Drew (1998, hereafter PSD98) used non-uniform meshes to resolve the subsonic parts of the flow and found that line driven disc winds in CV and YSO systems launch complex, time dependent, outflows when disc radiation is comparable to or stronger than the stellar component. Proga, Stone and Drew (1999, hereafter PSD99) used a more efficient treatment of the radiation force which allowed for the optical depth to be calculated using the Sobolev prescription without approximating the velocity gradient tensor.

An important question is what structure can affect line driven winds? Substructure on small scales or “clumpiness” can affect the flow by altering the ionization structure of the wind, via shielding outer parts of the flow (see for example Fullerton 2011). It can further provide new time dependent, observable features, since diagnositcs such as the IR continuum and subordinate lines like H are sensitive to (Puls et al. 1993). Clumping in disc winds has been used, albeit on different length scales, to explain in part line width ratios seen in BAL QSOs (Matthews et al. 2016)

One possibility for generating clumpiness is instabilities in the flow, such as the line deshadowing instability (LDI) inherent in line-driven winds (Lucy & Solomon 1970; MacGregor, Hartmann & Raymond 1979; Carlberg 1980; Owocki & Rybicki 1984). Owocki, Castor & Rybicki (1988) found in 1D spherically symmetric simulations that small amplitude waves in the subsonic part of the flow can be amplified by the line driving into non-linear waves in the supersonic part of the flow. The LDI does not operate however in the Sobolev approximation that is used in CAK-like treatment of line driving, leading us to explore alternative mechanisms for generating clumps.

In PSD98 the base of the wind is time dependent and clumpy, and matter sometimes falls back to the disc rather than being expelled. We propose that these failed winds provide a mechanism whereby an initial non-axisymmetry in the disc induces a non-axisymmetry in the wind and some matter that fails to escape the gravitational potential falls back to the disc generating new non-axisymmetries and perpetuating the cycle. Here we study this effect by first simulating an axisymmetric, line driven disc wind to test our Athena++ radiation pressure module. We compare this to a simulation where the initial condition is perturbed and study the resulting wind outflow properties, in particlar the degree of non-axisymmetry of the solution, as well as its stability.

The structure of this paper is as follows. In Section 2, we describe our numerical methods, including the implementation of the radiation force due to electron scattering and line driving in the hydrodynamics code Athena++. In Section 3.1, we describe the results of a 3D axisymmetric disc wind, which we use as a benchmark case. We compare this run in Section 3.2 to a 3D simulation that used the same parameters, but had an initial dependent, vertical, subsonic velocity perturbation in the disc midplane at . In Section 4, we discuss possible observational implications for these results. We conclude in Section 5, where we describe possible future directions of inquiry.

## 2 Numerical Methods

We performed all numerical simulations with the publicly available MHD code Athena++ (Gardiner & Stone 2005, 2008). The basic physical setup is a gravitating central object surrounded by a thin, luminous accretion disc that provides the radiation field to drive the gas that is optically thin to the continuum. We describe our models basic equations in Section 2.1. Our simulations required developing a module for calculating the radiation force due to electron scattering and lines, described in Section 2.2. For ease in comparison with previously published results we used the parameters for PSD98 “Run 2”, which was their fiducial run for line driven disc winds in CVs. These parameters, along with the boundary conditions are described in Section 2.3.

### 2.1 Basic Equations

The basic equations for single fluid hydrodynamics driven by a radiation field are

(1a) | |||

(1b) | |||

(1c) |

where , are the fluid density and velocity respectively and is a diagonal tensor with components P the gas pressure. For the gravitational potential, we use and is the total energy where is the internal energy. The total radiation force per unit mass is . The isothermal sound speed is and the adiabatic sound speed . We take a nearly isothermal equation of state where . We can compute the temperature from the internal energy via where is the mean molecular weight and other symbols have their usual meaning.

### 2.2 Radiation Force

We assume an axisymmetric, time independent radiation field is provided by a geometrically thin accretion disc along the midplane. The frequency integrated intensity is

(2) |

where is the radial position on the disc, is the inner radius of the disc, is the Thompson cross section per unit mass, the speed of light and

(3) |

is the disc Eddington number, where is the accretion rate in the disc (e.g., Pringle 1981). In addition, we assume the central object of radius is optically thick and is not a source of driving radiation. This affects the radiation field at small radii, as the central object effectively shadows the wind from the radiation from the inner disc.

We assume all gas is optically thin to continuum radiation and every point in the wind experiences a radiation force

(4) |

which is a sum of the contributions due to electron scattering and line driving . In this continuum optically thin approximation, the radiation force due to electron scattering is

(5) |

where is the normal vector from the disc to the point in the wind, is the solid angle spanned by the disc and the integration is carried out over the entire disc. We note that the above integrand is time independent and for each gridpoint in the wind the integration can therefore be carried out before the start of the hydrodynamics simulation.

We treat the radiation due to lines using a modification of the CAK formulation where the radiation force due to lines is

(6) |

and is the so-called force multiplier. The force multiplier parametrizes how many lines are effectively available to increase the scattering coefficient. We use the Owocki, Castor & Rybicki (1988, hereafter OCR) parametrization of the line strength, where working in the Sobolev approximation, it is a function of the optical depth parameter

(7) |

where is the thermal velocity of the gas and the velocity gradient along the line of sight

(8) |

can be expressed as a sum over elements of the shear tensor where are the components of the normal vector. The OCR formulation for the force multiplier is

(9) |

where and are constants, , and is related to the maximum force multiplier via . This force multiplier parametrization follows CAK for large optical depths but saturates to as optical depth becomes very small.

Stevens and Kallamn (1990) used detailed photoionization calculations to determine the radiation force driving the stellar wind in a massive X-ray binary (MXRB). They found the ionization state depends on the SED, in particular the hardness of the X-rays. Line driven wind simulations of QSOs (Proga, Stone & Kallman 2000; Nomura et al 2016) have accounted for this in a self-consistent way by including the heating/cooling primarily due to X-rays. Our simulations apply to CVs, where the X-ray flux is softer than in QSOs. We therefore approximate the flow as isothermal and leave careful modeling of the wind ionization structure for future work.

The radiation force due to lines (see eq. 6) has an integrand that is time dependent, because it depends on the velocity gradient and the gas density through the optical depth parameter (see eq. 7). This is computationally expensive because the integration over the disc must be carried out at every time step. PSD98 assumed that the dominant contribution to the velocity gradient is along the vertical direction,

(10) |

The radiation force due to lines is then

(11) |

where we have assumed that (see appendix C in PSD98 for a derivation). Since all time dependent factors are in , the integration over the disc must only be performed at the start of the simulation. PSD99 found that this approximation gave good agreement with simulations using the full expression for Q (eq. 8) so we make use of it here. We describe our numerical integration scheme, based on PSD99, in more detail in Appendix A.

### 2.3 Simulation Parameters

We chose parameters that correspond to PSD98 Run2. The central object has mass and radius and , respectively. The sound speed , which corresponds to a hydrodynamic escape parameter at the base of the wind. For this high HEP, thermal driving, which requires no more than (Stone & Proga 2009; Dyda et al. 2017), is negligible throughout the domain, .

We impose outflow boundary conditions at the inner and outer radial boundaries and axis boundary conditions along the axis. We assume a reflection symmetry about the midplane. In the direction, we impose periodic boundary conditions, where our 3D simulation covers a range . After every full time step we reset to , to and to . We also impose that the velocity is unchanged due to resetting density by also resetting the momentum i.e. where the superscript refers to the and timestep respectively. Generally we followed the setup of PSD98 as closely as possible, given that they used the Zeus 2D code (Stone & Norman 1992).

We choose a domain size and MPI meshblocks of size . Meshblocks of size 32 is the recommended smallest size for Athena++ to achieve good scaling with number of processors. The radial domain extends over the range with logarithmic spacing . The polar angle range is and has logarithmic spacing which ensures that we have sufficient resolution near the disc midplane to resolve the density length scale of the hydrostatic equilibrium disc that forms i.e to obtain . We use linear spacing in the direction and with this resolution expect to be able to resolve the highest frequency axial modes that can be excited in this wedge.

Initially, the cells along the disc are set to have , , and the pressure . In the rest of the domain , and .

We take the disc Eddington number to be . The thermal velocity of the gas is . The line driving parameters are chosen to be , and . The disc as a source of radiation is taken to extend outside the computational domain to .

We impose a density floor , which adds matter to stay above this floor, while conserving momentum, if the density ever drops below it. In addition, we have a pressure floor , so that our nearly isothermal condition is preserved.

## 3 Results

We describe in detail two 3D simulations of line driven disc winds. The first, a 3D analogue of the “Run 2” simulation from PSD98 is used as a validation of our implementation of line driving in Athena++ and to ensure that axisymmetry is explicitly preserved in the approximation. These results are shown in Section 3.1. We then perturb the intial conditions along the midplane with a vertical, subsonic, sinusoidal, velocity perturbation. The results of this perturbed run are shown in Section 3.2 where we compare properties to the unperturbed run and also quantify its non-axisymmetric features.

### 3.1 Axisymmetric Reference Run

We performed a 3D simulation of a line driven disc wind to test that Athena++ preserves the axisymmetry of the initial condition and of the line driven wind. As a consistency check, we also compared our results to the 2D simulations of PSD98.

At the start of the simulation gas enters the simulation region from the midplane, accelerated by the pressure gradient and a radiation force that progressively ramps up over the first 100 s, corresponding to 10% of the simulation runtime. After roughly 200 s a line driven wind forms, with a wind opening angle of . Most of the mass flux is carried in a fast stream, within a roughly wide wedge, below the upper wind envelope. Above this wedge, there is a polar funnel region with density so low that it is negligible. Below this wedge, the velocity is low, and the flow turbulent enough that the mass flux is insignificant. The wind is time dependent, though the streamlines do not change significantly in time.

The structure of the wind can be further characterized by considering its properties at the outer radial boundary. We define the integrated mass flux

(12) |

In Fig. 1 (top panel), we plot the momentum density (black, solid line) and integrated mass flux (red, dashed line) as a function of at the outer boundary at . We also indicate the total mass flux from the full simulation region (see the top left corner of the top panel), assuming top/bottom symmetry. The momentum density is concentrated around , where the stream reaches the outer radius. This can also be seen from the integrated mass flux that levels off for larger values of . A spike above occurs near the midplane because of the much higher density disc and small subsonic that is positive in this snapshot. Therefore we impose a cutoff of when measuring fluxes of both simulations, because we are interested in wind properties.

In the bottom panel of Fig. 1, we plot the density and velocity structure of the wind as a function of at the outer radial boundary. The polar region is nearly void of matter with for . Then there is a fast-stream region where the velocity peaks to its maximum value and the density of the wind increases to . The density remains nearly constant at this value for between and . We define the characteristic velocity of the outflow as the velocity of the wind when this reference density is reached, roughly . We measure a wind mass loss rate to be . We find our 3D solutions both confirm the results of PSD98 and maintain axisymmetry to machine precision.

### 3.2 Non-Axisymmetric Perturbation

The main subject of this investigation is to explore non-axisymmetric effects in line driven disc wind solutions. Non-axisymmetry can be introduced by having forces that explicitly break the axisymmetry or by perturbing the axisymmetric initial conditions. Here we want to establish the level of non-axisymmetry generated in the wind by perturbing the initial conditions. We introduce an initial vertical, subsonic, sinusoidal, velocity perturbation

(13) |

in the disc midplane. We note this perturbation is the lowest mode consistent with our periodic boundary conditons in the axial direction i.e., and . We compare both the global and local properties of this perturbed solution with the axisymmetric wind from the previous section.

#### 3.2.1 Global Properties

First, we consider the global properties of the flow. In particular, we consider the mass flux

(14) |

momentum flux

(15) |

and kinetic energy flux

(16) |

across the outer radial boundary. We plot these quantities in Fig. 2, showing each of them for the unperturbed (red, solid line) and perturbed (black, solid line) solutions as well as the relative difference in fluxes between solutions (black, dashed line). At early times, the wind has not yet reached the outer boundary and fluxes are nearly zero. After an initial transient outburst at time , the solutions settle down and exhibit a relative difference of up to 100%. We note that both solutions are highly time dependent, and there is no clear trend about which solution has a more powerful wind. However, we can conclude that the perturbation has a non-negligible effect on the gross outflow properties.

#### 3.2.2 Local Deviations

To quantify departures from axisymmetry we define the relative standard deviation

(17) |

where is the averaged density distribution. This metric quantifies how density at fixed deviates, on average, in the direction. We also define the relative maximum deviation

(18) |

This provides a measure of the largest deviation from axisymmetry at a fixed . Because emission scales like , we consider over(under)-dense regions observationally interesting. To potentially emit times more (less) than the azimuthal average, the density must differ by a factor of about three from the azimuthal average. We therefore refer to over-dense regions where and under dense regions with .

We identify three qualitatively different regions in our simulations - a nearly hydrostatic disc, a wind base and the main part of the wind. To sample these regions we choose points at , and inside of each of these respective regions. These points lie on a single streamline, and are used in later analyses.

In Fig 3, we plot the averaged density and velocity at . We mark the and contours (yellow, dashed and solid contours, respectively) as well as over-dense (green contours) and under-dense (red contours) i.e. features which differ from the azimuthal average by a factor of 3. We indicate the location of our three representative points in the “disc” (“D”, orange), “base” (“B”, blue) of the wind and the main “wind” (“W”, purple), later used in more detailed analysis. For clarity we have not plotted vectors and contours where .

Qualitatively, this solution resembles the unperturbed solution. The standard deviation countours indicate that non-axisymmetries are small in the disc , but are roughly (dashed yellow contours) in most of the wind. There are some small regions at the base of the wind where (solid yellow contours).

In Fig. 4, we plot (solid) and (dashed) as a function of time for the “wind”, “base” and “disc” points in their respective color. Within the disc the effects of the perturbation are nearly constant at late time, with . This is also evident in Fig 3 from the location of the contour that follows the contours of the disc. As the wind is turned on, the deviation grows to at the base and remains at this level. In the fastest part of the wind, the deviation is , smaller than at the base by an order of magnitude.

These results are consistent with two other numerical tests we performed with line driving turned off. We perturbed a disc that is too cold to thermally launch outflows (), and found non-axisymmetries in the disc to be of the same order of magnitude as with line driving turned on. We also considered a disc that is hot enough to thermally launch outflows (), and found that non-axisymmetries were restricted to the streams at the polar/wind boundary as non-axisymmetric features in the main wind were advected from the simulation region.

We note that , indicating that the dominant contribution to the relative standard deviation is primarily from the point where the deviation is maximal. In particular, we note how both distributions exhibit the same time behaviour. We can see this as well from the over-dense and under-dense contours, which have considerable overlap with the contours. It is also important to note that there are spatially more under-dense regions than over-dense regions. This is expected since if mass in over and under-dense regions is identical then we would expect the volume of over-dense regions to be 9 times smaller than the under-dense regions.

#### 3.2.3 Clumpiness

To quantify and identify the origin of the non-axisymmetry we measure the clumping factor defined as

(19) |

where the angular brackets represent a spatial average. The average may be taken over different domains, and we outline our choices for the domain below.

We first consider azimuthal clumping,

(20) |

where the spatial average is taken over but at fixed e.g., . Azimuthal clumping is a useful metric because for axisymmetric simulations and is therefore a local measurement of deviations from axisymmetry. In Fig. 5 we plot the azimuthal clumping for the disc, base and wind point. In the disc and the wind, we find , indicating little deviation from axisymmetry. However, in the base consistent with our finding that over/under dense features are primarily at the base of the wind.

Alternatively, we may consider clumping over a surface

(21) |

where the spatial average is over and but at fixed i.e. . We take the surface clumping over the region of the fast stream . By this metric, even 2D axisymmetric winds are clumpy, see for example Fig. 2 of PSD98.

We may also consider clumping over a volume

(22) |

where the spatial average is over , and i.e. . We take the volume clumping over the region of the fast stream . For axisymmetric simulations .

In Fig. 6, we plot volume averaged clumping (top panel) for the representative time of the unperturbed (red, solid line) and perturbed (black, solid line) simulations. In the panels below, we plot the surface averaged clumping at fixed and .

We find that the unperturbed wind has the highest volume averaged clumpiness over the duration of the outflow , with , and it varies on timescales . This same quantity for the perturbed run has a smaller range . The surface averaged clumpiness exhibits the same range as the volume avergaed clumpiness, but varies on much shorter time scales .

The features in the distribution vary with as a function of time. We plot three dashed diagonal lines, indicating the motion of three troughs in the clumpiness distribution. The diagonal lines indicate how a uniform sized feature moving at constant angular velocity propagates in time. In this case, we see a dip in propagating for three orbits before dissipating. This corresponds to an approximate angular velocity Hz. This is what one might expect for features of size propagating at . This is not to say that our result is resolution dependent, rather that the timescale of variations in is set by the orbital velocity of clumps and not their outflow velocity. In the wind, the region where is maximum does in fact correspond to the part of the flow with this angular velocity (). The clumpiness of the wind is therefore dominated by material near the base of the wind,where the poloidal velocity is low. By contrast, varies on timescales for material to be ejected from the wind where . Disc perturbations can therefore cause variations of roughly a factor of 2 in , but on timescales about 10 times shorter than similar sized variations in which are due to radial variations in the flow.

#### 3.2.4 Fourier Transform

To further investigate the time variability of the non-axisymmetry, we consider the Fourier transform,

(23) |

where we apply a standard fast fourier transform (FFT) algorithm on our data. To eliminate any effects due to transient modes associated with the initial wind launching, we consider times s, as the time during which an outflow is established. In Figure 7, we plot the normalized power spectrum distribution (PSD) for density and each of the velocity components , and as a function of frequency for the disc, base and wind points. The normalization is chosen so that integrating the PSD over all frequencies equals unity.

The PSD of the disc have the most well defined features, with the PSD of each quantity having a strong peak at fundamental frequency , where is the local Keplerian frequency. We label higher multiples of this fundamental frequency via . The PSD is dominated by a zero mode. The unperturbed simulation has prominent peaks at and , while the perturbed solution has ,,, and . The PSD has no zero mode and the unperturbed solution has peaks at and the perturbed solution at , . The unperturbed PSD has peaks at , , and whereas the perturbed run at ,,,, and . The has a dominant zero mode with the perturbed run having peaks at and . The PSD being less noisy in the disc is consistent with the disc being in near hydrostatic equilibrium and maintaining the highest degree of axisymmetry.

The PSD for the base has the greatest number of modes with fundamental frequency .The unperturbed and perturbed runs have PSD peaks at all modes . The unperturbed PSD has a mode and the perturbed run has , , , and . The unperturbed PSD has a mode and the perturebd run , and . The unperturbed PSD has a , and mode and the perturbed run has ,, , and modes. The presence of very many modes is consistent with our finding that this is the most disordered part of the flow.

The wind PSD, all feature a large peak at low frequency, but unlike in the base and the disc it is hard to identify particular modes. The wind does not have any dominant modes because they are advected out of the domain.

Certain modes may be present in our results because we have only simulated a wedge with . Therefore we ran a test simulation with the full disc and concluded that this was not the case. A general feature is the perturbed wind has more modes, and these modes are at higher frequencies, than the unperturbed wind. A complete understanding of the modes present, even in the comparitively simple nearly hydrostatic disc is challenging and is left to future work.

## 4 Discussion

### 4.1 Form of Perturbation

Here we present the first study of non-axisymmetric features of line driven disc winds. We restricted ourselves to vertical, subsonic, sinusoidal velocity perturbations of the form (13), though other forms of perturbation are possible. We varied a few of the perturbation parameters, to check that our results are reasonably robust. In particular, we increased the magnitude of the perturbation from 0.1 to 0.3 . We also doubled the wavenumber of the perturbation i.e. from 8 to 16.

The overall conclusion from these experiments is that varying these parameters does not have a strong effect on the properties of the solution. Increasing the amplitude of the perturbation did not introduce any new modes. Doubling the wavenumber excited a mode in in the base and but kept other modes unchanged. This suggests the excited modes are the natural frequencies of the system and not due to the initial perturbation. After the inital mass outflow at s, increasing the magnitude or changing the wavenumber led to a change in the wind fluxes relative to our fiducial perturbation. This was a smaller difference than between our perturbed and unperturbed runs shown in Fig. 2.

Removing the linear ramp up of the line driving led to differences at early times, as the initial outburst occured at roughly s. However, after this initial discrepency the flow settles down to the same difference as seen in our other experiments.

This suggests that for this form of the perturbation, the amount by which the gross outflow properties (mass, momentum and energy fluxes) can be varied is a factor of . This is comparable to the variation we see in these quantities because of the time dependent nature of line driven disc winds.

We also performed a simulation over the entire toroidal range to ensure that simulating only a wedge would not introduce spurious modes into the solution. We did not find evidence that the wedge run had additional modes. In the case this test was possible since the runtime was under two weeks running on 32 processors. When we relax this approximation this will be computationally expensive, so this null result is encouraging in this context.

The precise form of the perturbation does not seem to have an effect on the Fourrier modes, suggesting these are the natural modes of the system. Non-axisymmetric density features grow at the base of the wind, unlike in a persistent thermal wind where they are advected away. We speculate this growth is due to gas falling back on the disc after failing to launch. This scenario will be investigated in future work.

### 4.2 Observations

In the previous section, we have primarily discussed how a disc perturbation breaks the axisymmetry of the wind solution and argued that from the point of view of simulations these are non-trivial i.e effects. However, these axisymmetries are also observationally relevant.

First we point out that axisymmetric winds can result in non-axisymmetric observables. Velocity dependent observables (Doppler shifted lines) depend on the velocity component projected onto the line of sight. Because the line of sight typically breaks the natural symmetry of the system (axial symmetry in the case of disc winds), this leads to non-axisymmetric observables. This was done for example in the context of BAL QSOs where axisymmetric disc wind solutions generated non-axisymmetric, time dependent line profiles (Proga, Rodriguez-Hidalgo & Hamann 2012). Therefore, we might expect the non-axisymmetric wind solutions we have found to further increase the time variability of these profiles.

Consider the emission measure of a particular line, that scales like . The maximum relative dispersion in , Fig. 4, shows that the density varies by . This would therefore lead to an increase in a factor of in the emisson measure at these points. These are of course over-estimates, because these are the points at which the non-axisymmetric features are strongest. However, it does suggest that if we are to calculate accurate synthetic line profiles we should use 3D wind solutions.

The effect on absorption measures of lines will be less pronounced since these vary . Therefore, we would expect these non-axisymmetric features to cause factors of at most . Observationally such variations are challenging to measure because variations of this magnitude are seen in axisymmetric simulations, as a result of the winds inherent time dependence.

A key to disentangling these effects may be the timescales over which they operate. As an example, we found the clumpiness of the wind to vary on a time scale of s, whereas the clumpiness along specific sightlines varied on s due to orbiting clumps at small radii. These timescales are too short to resolve for CV systems. However, for QSOs where the Keplerian timescale is about times longer variations due to clumpiness may be observable. Given that line emission depends sensitively on the locations of resonance points, it may be possible to see variations in the emission measure on the time scales of these clumps moving at Keplerian velocity. A detailed calculation of line profiles should be able to determine this.

## 5 Conclusion

We performed 3D simulations of line driven disc winds using the Athena++ code where the velocity gradient was calculated using the approximation. In the first 3D simulation we used axisymmetric initial conditions and found driven winds using maintain axisymmetry and produce outflows consistent with results from 2D axisymmetric simulations.

We then introduced a vertical, subsonic, non-axisymmetric velocity perturbation in the disc midplane. We compared the resulting outflow to the unperturbed simulation and found, the global properties of the outflow are largely unchanged, though fluxes may vary by up to 100 % at any given time. The disc perturbation produces azimuthal clumping factors up to a factor of in the base of the wind. Clumpiness in the overall flow does not increase but varies on timescales much shorter than for axisymmetric winds.

The disc is least affected by the perturbation where deviations from axisymmetry . In the wind base the deviation is largest with but decreases to in the fastest moving parts of the wind because any perturbations in the wind are carried out of the simulation domain. Non-axisymmetric features are, volume wise, primarily under-densities, up to a factor of 3 lower than the toroidal average. There exist also over-dense regions where densities are a factor of 3 higher than the toroidal average. These over/under dense regions exist primarily at the base of the wind, above the hydrostatic disc but upstream from the fastest moving parts of the wind.

The non-axisymmetries have potential observational effects. We find regions with difference in densities, which leads to a factor of difference in the emission from these regions. The clumpiness of the flow is also found to vary on time scales of the orbital motion, rather than the outflow velocity for axisymmetric flows. This too may be potentially observable with sufficient time resolution.

We have found that a small initial perturbation leads to non-axisymmetries in the flow. However, these asymmetries tend to saturate and not grow without bounds. This suggests that in order for non-axisymmetries to grow beyond as we have found, requires additional physics beyond what is included in this model.

We calculated the velocitiy shear in the approximation . In a later paper we will explore the case where this approximation is relaxed. This will lead to an azimuthal component of the radiation force, potentially changing the growth of non-axisymmetries beyond what we presented here.

Another approach is where we account for the shielding of the wind by matter at the base of the wind. We found that most of the non-axisymmetric clumpiness is due to matter at small radii which is launched above the disc but fails to be expelled. This matter can shield the wind from ionizing radiation. This provides a mechanism to couple the base of the wind to the broader outflow, which our optically thin treatment does not capture. This requires a more proper treatment of the radiation transfer through the wind. We also plan on including radiation from the central object.

Finally we have assumed that the radiation field is axisymmetric. This may not be the case, as spiral structures can form due to accretion from a donor star in a CV system. Alternatively, perturbations driven by the magnetorotational instabiity, can generate non-axisymmetries in the radiation field.

We conclude it is worthwhile to devote future efforts to simulating line driven disc winds in full 3D, as non-axisymmetric effects may affect the flow and its observational features.

## Acknowledgements

This work was supported by NASA under ATP grant NNX14AK44G. The authors also acknowledge useful discussions with Zhaohuan Zhu and Tim Waters.

## References

- [1]
- [2] Carlberg R. G., 1980, ApJ, 241, 1131
- [3]
- [4] Castor, J.I., Abbott, D.C., Klein, R.I., 1975, ApJ, 195, 157
- [5]
- [6] Dyda, S., Dannen, R., Waters, T., Proga, D., 2017, MNRAS, 467, 4161â4173
- [7]
- [8] Feldmeier, A., Shlosman, I., ApJ, 526, 344F
- [9]
- [10] Friend, D.B., Abbott, D.C., 1986, ApJ, 311, 701
- [11]
- [12] Fullerton A. W., 2011, in Neiner C., Wade G., Meynet G., Peters G., eds, Proc. IAU Symp. 272, Active OB Stars: Structure, Evolution, Mass Loss, and Critical Limits. Cambridge Univ. Press, Cambridge, p. 136
- [13]
- [14] Gardiner,T.A., Stone. J.M., 2005, J. Comput. Phys., 205, 509
- [15]
- [16] Gardiner,T.A., Stone. J.M., 2008, J. Comput. Phys., 227, 4123
- [17]
- [18] Gayley K. G., 1995, ApJ, 454, 410
- [19]
- [20] Higginbottom, N., Proga, D., Knigge, C., Long, K.S., Matthews, J.H., Sim, S.A., 2014, ApJ, 789:19H.
- [21]
- [22] Higginbottom, N., Proga, D., 2015, ApJ, 807:107.
- [23]
- [24] Icke, V., 1980, AJ, 85, 329
- [25]
- [26] Icke, V., 1981, ApJ, 247, 152
- [27]
- [28] Lucy L. B., Solomon P. M., 1970, ApJ, 159, 879
- [29]
- [30] MacGregor K. B., Hartmann L., Raymond J. C., 1979, ApJ, 231, 514
- [31]
- [32] Matthews,J. H., Knigge, C., Long, K.S., Sim, S.A., Higginbottom, N., Mangham, S.W., 2016, MNRAS, 458..293M
- [33]
- [34] Murray, N., Chiang, J., Grossman, S.A., Voit, G.M., 1995, ApJ, 451, 498
- [35]
- [36] Nomura, M., Ohsuga, K., Takahashi, H.R., Wada, K., Yoshida, T., 2016, ASJ, 68(1), 16
- [37]
- [38] Owocki S. P., Castor J. I., Rybicki G. B., 1988, ApJ, 335, 914
- [39]
- [40] Owocki S. P., Rybicki G. B., 1984, ApJ, 284, 337
- [41]
- [42] Pauldrach A., Puls J., Kudritzki R. P., 1986, A&A, 164, 86
- [43]
- [44] Pereyra, N.A., 1997, PhD thesis, Univ. Maryland
- [45]
- [46] Pereyra, N.A., Kallman, T.R., Blondin, J.M., 1997, ApJ, 477, 368
- [47]
- [48] Press, W.H., Teukolsy, S.A., Vetterling, W.T., Flanney, B.P., Numerical Recipes in C++, 1988, Cambridge UP.
- [49]
- [50] Pringle, J.E., 1981, ARA&A, 19, 137
- [51]
- [52] Proga, D., Stone, J., Drew, J.E., 1998, MNRAS, 295, 595
- [53]
- [54] Proga, D., Stone, J., Drew, J.E., 1999, MNRAS, 310, 476
- [55]
- [56] Proga, D., Stone, J., Kallaman, T.R., 2000, ApJ, 543, 686P
- [57]
- [58] Proga, D., Rodriguez-Hidalgo, P., Hamann, F., 2012, ASPC, 460, 171P
- [59]
- [60] Puls, J., Pauldrach, A.W.A., Kudritzki, R.P. Owocki, S.P., Najarro, F., 1993, in:Reviews in Modern Astronomy 6, Springer, Berlin, p 271
- [61]
- [62] Owocki, S.P., Castor, J.I., Rybicki, G.B., 1998, ApJ, 335, 914
- [63]
- [64] Stevens, I.R., Kallman, T.R., 1990, ApJ, 265, 321S
- [65]
- [66] Stone, J. M., Norman, M. L. 1992, ApJS, 80, 753
- [67]
- [68] Stone, J., Proga, D., ApJ, 2009, 694:205â213
- [69]
- [70] Vitello P.A.J., Shlosman, I., 1988, ApJ, 327, 680
- [71]

## Appendix A Radiation Force Module

We describe our module for computing the line driving force for a disc and stellar radiation field.

### a.1 Coordinate System

We use spherical coordinates to denote coordinates in the wind. These are the natural coordinates to use since these are the coordinates in the Athena++ simulation. The disc is centered on the origin and lies in the XY plane with coordinates and inner and outer radii and respectively. The star is also centered on the origin with radius .

Computing the radiation field requires integrating over the disc and star. PSD98 have shown that using cylindrical coordinates causes a loss of accuracy near the disc surface because of geometric foreshortening effects. PSD99 subsequently found that a spherical coordinate system centered on the wind point does not suffer from this issue so we use these coordinates as well. We choose the angular variable such that corresponds to points behind the star. We show the different coodinates in Fig. 8.

To compute the radiation force due to electron scattering (eq. 5) and line driving (eq. 6) we require the normal vector and the solid angle . The intensity profile (eq. 2) requires the disc radius

(24) |

where

(25) |

is the distance from the wind point to the disc point. It is useful to further perform the change of variables .

### a.2 Numerical Integrator

Our numerical scheme is inspired by the 3D integrator described in Press et al. (1988) but modified to integrate 2D functions. Suppose we are performing the integral

(26) |

over a region , . We approximate this integral as

(27) |

where are weights, the appropriate abscissas and we are using integration points in the x and y directions respectively. The problem consists of choosing the weights and abscissas appropiately. In general we have found Gaussian quadrature performs best in terms of rapid convergence.

We specify the function , the y boundaries and and the x boundaries and . In the first step, the algorithm chooses abscissas and weights in the interval . Then for each the abscissas and weights are chosen in the interval . The integral is then computed using (27) where we integrate over the y variable first followed by the integral over x. For both the stellar and disc integrals we take and as our integration variables. We use integration points. We artificillay break up the integral along so as to avoid numerical asymetries in the direction.

For electron scattering (eq. 5) and line driving using the approximation (eq. 11), we store the respective integration sum during the initialization step for each point in the wind. When computing the line driving force without approximation (eq. 6), we store the components of for each integration point and their corresponding weight. We need the directional vector to compute the Sobolev parameter and storing these allows us to not require calculating weights every timestep.

### a.3 Testing

We performed various tests to ensure our implementation of the radiation force was accurate to 1%. We compared our numerical integrator to the analytic expressions for electron scattering for a uniform intensity and “Newtonian” disc from Feldmeier & Shlosman (1999). We checked our implementation in the hydro by comparing to the 1D spherical CAK solution, as well as the CAK solution with finite disc correction (Pauldrach, Puls & Kudritzki 1986; Gayley 1995). Finally we compared with results from PSD98, which explored cases with both stellar and disc radiation fields with Eddington parameter covering the range .

Our tests led to us implementing two numerical features to ensure that the line driving is well behaved in Athena++, which were not necessary in the work of PSD98 and PSD99 who used the hydrodynamics code Zeus 2D (Stone & Norman 1992). We imposed a ceiling on the strength of the line driving force. Whenever the momentum is to be updated due to the contribution from line driving, we ensure that the new velocity with the Keplerian velocity at radius . We found that this prevents the formation of spuriously large velocities along the axis and is only triggered for the intial transient. However, we find late time solutions that have velocities .

To further reduce transients, we turn on the radiation field using a linear ramp over s, which corresponds to roughly 5 inner disc orbits. These numerical features that have been implemented to ensure the line driving is well behaved but they are only necessary for high disc and stellar luminosities. i.e and approximately , which we do not consider in this work. However, we describe them here for completeness as we will consider higher luminosity cases in future studies.