format: text

text:
# Modeling the Spectral Signatures of Accretion Disk Winds: A New Monte Carlo Approach

## 1 Introduction

Non-magnetic cataclysmic variables stars (CVs) are semi-detached binary systems in which a late-type companion star loses mass to a white dwarf (WD) primary via Roche-lobe overflow onto an accretion disk. The discovery of winds in high-state CVs – from the presence of blue-shifted absorption troughs and P-Cygni profiles in the ultraviolet (UV) resonance lines – dates back about two decades (Heap et al., 1978; Córdova & Mason, 1982). However, most of the fundamental theoretical questions concerning the nature and origin of this mass loss are only now beginning to be addressed (Drew & Proga, 2000).

Recent work suggests that the impact of the outflowing material on observations of CVs may have been underappreciated, in that its presence may leave more numerous and subtle observational signatures than just the shapes of the “classical” UV wind lines. For example, it has been suggested that a powerful disk wind may be responsible for the peculiar behavior of the SW Sex stars (Honeycutt, Schlegel, & Kaitchuck, 1986; Dhillon & Rutten, 1995), and that the often single-peaked optical emission line profiles seen in these systems (and also in other nova-like variables \[NLs\]) are the direct consequence of a wind-induced velocity gradient in the line-emitting material (Murray & Chiang, 1996). Similarly, Knigge et al. (1997) suggested that virtually all of the lines in the UV spectrum of a typical high-state CV, Z Cam, are formed in the outflow, either in the supersonic portion of the wind or in a lower velocity portion of the wind near the interface with the disk photosphere. Based on an analysis of low resolution (1140-8950 Å) spectra of the eclipsing novalike variable UX UMa, Knigge et al. (1998) then suggested that continuum emission from an optically thin disk wind and/or atmosphere provided one possible explanation for the absence of a clear Balmer jump and the flatter than expected UV continua of high-state CVs (Wade, 1984, 1988). Optical spectra of at least one CV, BZ Cam, shows P-Cygni features in the hydrogen and helium lines that are most easily interpreted in terms of a highly variable wind (Ringwald & Naylor, 1998). In addition, the EUVE spectrum of the dwarf nova (DN) U Gem in outburst appears to be dominated by strong wind features to such an extent that any attempt at a detailed spectral analysis must account for the outflowing material (Long et al., 1996). Finally, there is also direct evidence that CV winds act as partially ionized absorbers of (soft) X-rays in high-state, non-magnetic systems (Baskill, Wheatley, & Osborne, 2001).

The first attempts to interpret UV line profiles in high-state disk-dominated CVs were quite naturally based on wind modeling theory that was being developed to interpret the UV spectra of early-type stars (Mauche & Raymond, 1987; Drew, 1987; Vitello & Shlosman, 1988). These studies, along with a burgeoning amount of observational analysis, indicated the mass loss rates in CV winds are a substantial fraction ($``$10%) of the accretion rate, that the acceleration length scale is substantial (20-100 $`R_{wd}`$), and that the winds are strongly affected by rotation, presumably reflecting an origin in the rapidly rotating inner disk.

The evident importance of rotation led to a second generation of radiative transfer codes that provided for azimuthally symmetric flows. In the CV wind code described in Shlosman & Vitello (1993), hereafter SV93, the wind ionization structure is computed first and the radiative transfer then calculated using the Sobolev (SA) and disconnected approximations. This method is comparatively fast, but does not treat multiple scattering correctly in detail. It was applied by Vitello & Shlosman (1993) to show that, depending on the wind geometry, wind mass loss rates of order 1-10% of the accretion rate could account for the C IV line profiles of RW Sex, RW Tri and V Sge . By contrast, the Monte Carlo (MC) code described by Knigge, Woods, & Drew (1995), hereafter KWD95, simulates the transfer of resonance line photons through an outflow exactly, but does not perform an ionization calculation. It was used by Knigge & Drew (1997) to model high resolution time-resolved spectra of the C IV region in UX UMa to show that a relatively dense, slowly outflowing transition region between the disk photosphere and the fast wind was necessary to understand the profiles in that system.

Both of these radiative transfer codes have limitations. Neither method solves self-consistently for the thermal balance in the wind (constant wind temperature is assumed); neither is designed to allow for the thermal creation or destruction of line photons in the wind (but see Knigge et al., 1997, for an approximate treatment in the MC context). Finally, both of these codes are intended for modeling of individual spectral lines; as a result it is usually unclear whether a wind model intended to reproduce the C IV profile has any predictive power for other wind signatures.

Here, we present a new, hybrid SA/MC method designed to provide improved spectral modeling of azimuthally symmetric winds. Various aspects of our method build on the work of Lucy and collaborators (Mazzali & Lucy, 1993; Lucy, 1999a, b) who have to-date considered only spherically symmetric situations. In addition to combining the respective strengths of the SA and MC methods in the (conservative) transfer of resonance lines, our method solves self-consistently for the thermal and ionization balance in the outflow and accounts for the destruction and creation of both line and continuum photons in the wind. As a result, our method can be used to synthesize any part of the emergent, wind-affected spectrum, allowing detailed analyses of previously neglected wind signatures.

## 2 Method

### 2.1 Basics

The idea behind our Monte Carlo algorithm is simple: a large number of suitably randomized and distributed energy packets (“photons”) are generated to represent the radiation field of the accretion disk, the WD, the boundary layer (BL) and the wind itself. These photons are followed as they make their way through the wind, where they can be scattered by resonance lines, scattered by electrons, absorbed in photoionization events or destroyed via collisional de-excitations. After accounting for these interactions, all photons escaping towards the observer are binned in wavelength to yield the wind-attenuated spectrum of the system as viewed from the observer’s viewing angle.

Not surprisingly, the implementation of this simple idea is more complex. Each complete run of our code actually consists of two separate sets of Monte Carlo simulations. In the first, the ionization and temperature structure of the outflow is calculated by iteration. To this end, an initial wind structure is assumed first, and a family of photons representing the entire system luminosity is generated and followed through the outflow. Photons scatter when the scattering optical depth

$$\tau _{scat}ln(1p)$$
(1)
where p is a random number between 0 and 1 generated with a uniform distribution. The characteristics of all photons that pass through each grid cell are stored and used to obtain better estimates for the ionization state and temperature of that cell. A new iteration is then started by generating a new family of photons. After a number of iterations, the temperature and ionization structure of the outflow converges. This converged structure is used as a fixed input to the second set of MC simulations. In this set, the detailed emergent spectrum is calculated for a user-defined wavelength range. This can be as narrow as a single line profile or as wide as the entire spectral energy distribution. The wind structure is not (and, generally speaking, cannot be) refined any further at this stage, since only photons in a restricted frequency range are generated.

In the following sections, we will provide a more detailed look at the inner workings of our code. We will begin with a description of the two kinematic wind models we have implemented to date.

### 2.2 Wind Models

The underlying geometry of our Monte Carlo code is cylindrical. The center of the WD defines the origin of the grid. The code is designed to handle winds that are symmetric with respect to reflection through the disk plane, as well as azimuthally symmetric with respect to the disk axis. As part of the initialization process, a grid of points with logarithmic (or linear) spacing is established, onto which the model of the wind is impressed. Linear interpolation is used to assure that variables, e.g. velocity and velocity gradients, are continuous throughout the wind. Provision is made for defining regions in which the wind exists and regions in which it does not. These features make it straightforward to support various wind model geometries. Despite the fact that the grid is cylindrical, one can accommodate models in which the wind occupies a biconical region emerging at arbitrary angles from the disk, or alternatively one that is approximately spherical. In the discussion that follows, we will use terminology specific to CVs, but the code is not limited to that astrophysical situation.

We currently implement the two kinematic descriptions of CV winds used by SV93 and by KWD95 in our code. We selected these because they are the two prescriptions that have been used most extensively to model C IV line profiles. Both are biconical wind flows, and both assume that material flows out from the disk along streamlines. Both separate the velocity into poloidal ($`v_l=\sqrt{v_\rho ^2+v_z^2}`$) and azimuthal components $`v_\varphi `$. Both assume that the outflowing material shares the Keplerian speed of the accretion flow at the footpoint of a streamline and that the material in the wind conserves specific angular momentum on its way out of the system. As a result the azimuthal velocity is large close to the disk, and is largest at the inner edge of the disk where the Keplerian velocity of the disk is greatest.

In the SV93 prescription (see Fig. Modeling the Spectral Signatures of Accretion Disk Winds: A New Monte Carlo Approach), the wind emerges from the disk between $`r_{min}`$ and $`r_{max}`$ along streamlines whose orientation with respect to the system are described by an angle111We note in passing that in a few cases we have consolidated for clarity within this article the notation of the two models. We trust that this will not be confusing for those who refer back to SV93 and KWD95

$$\theta =\theta _{min}+(\theta _{max}\theta _{min})x^\gamma $$
(2)
where

$$x=\frac{r_or_{min}}{r_{max}r_{min}}$$
(3)
and $`r_o`$ refers to the position of the footpoint of a streamline. A highly collimated flow will result if the difference between $`\theta _{min}`$ and $`\theta _{max}`$ is small while $`\gamma `$ controls the concentration of the stream lines toward either of the two wind boundaries.

The poloidal velocity along the streamlines is defined to be

$$v_l=v_o+(v_{\mathrm{}}(r_o)v_o)\left(\frac{(l/R_v)^\alpha }{(l/R_v)^\alpha +1}\right)$$
(4)
The scale length $`R_v`$ and the exponent $`\alpha `$ control the acceleration of the wind between a fixed velocity $`v_o`$, normally 6 $`\mathrm{km}\mathrm{s}^1`$, at the base of the wind and the terminal velocity $`v_{\mathrm{}}(r_o)`$. The terminal velocity of each streamline varies depending on the location of the streamline in the inner and outer disk, being characterized as a fixed multiple of the escape velocity at the footpoint of the streamline. Thus the poloidal velocity is greatest for stream lines that originate from the inner regions of the disk, since the gravitational potential that must be overcome is greatest there.

The mass loss per unit surface area $`\delta \dot{m}/\delta A`$ of the disk is controlled by a parameter $`\lambda `$ such that

$$\frac{\delta \dot{m}}{\delta A}\dot{m}_{wind}r_o^\lambda cos(\theta (r_o))$$
(5)
With this prescription, the overall mass loss rate declines with radius if $`\lambda `$ is somewhat less than -2.

To use the SV93 prescription, therefore, one must provide the basic parameters of the system, the mass of the WD, the accretion rate, the inner and outer radius of the disk, and in addition, for the wind $`\dot{m}_{wind}`$, $`r_{min}`$, $`r_{max}`$, $`\theta _{min}`$, $`\theta _{max}`$, $`\gamma `$ ,$`R_\nu `$, $`\alpha `$, $`\lambda `$, and the multiple of the escape velocity to be used for $`v_{\mathrm{}}`$.

The KWD95 characterization of the wind has somewhat fewer parameters. In it, as indicated in Fig. Modeling the Spectral Signatures of Accretion Disk Winds: A New Monte Carlo Approach, the outflow in each hemisphere of the disk system is described as if it originated from a point a distance D on the other side of the disk plane. The inner and outer boundary of the wind is determined by the WD radius and the outer radius of the disk. For large D the wind will be highly collimated along the axis of the disk, and for very small D it will be highly collimated in a thin region with large opening angle. The least collimation occurs when

$$D=\sqrt{R_{disk}R_{wd}}.$$
(6)
For example, if $`R_{disk}=50R_{wd}`$, then the wind with the least collimation angle has a value of D of $`7.07R_{wd}`$ and extends from $`\theta _{min}`$ of 8° to $`\theta _{max}`$ of 82°. On the other hand, if D were 10$`\times `$ larger, then $`\theta _{min}`$ would be 0.8° and $`\theta _{max}`$ would be 35°.222Rather than use D, the distance along the z axis from which streamlines originate, KWD95 characterized their models in terms of the ratio $`D/D_{min}`$ where $`D_{min}`$ corresponds to the point at which the solid angle of the disk is maximized from an observer located at $`D_{min}`$. This is approximately $`0.5\sqrt{R_{disk}R_{wd}}`$, with the exact value of the coefficient depending on the ratio of $`R_{disk}`$:$`R_{wd}`$. We have chosen to deprecate the use of $`D/D_{min}`$, except when it is necessary to refer back to models calculated by KWD95.

In their description, the mass loss rate per unit surface area

$$\frac{\delta \dot{m}}{\delta A}T(R)^{4\alpha }$$
(7)
For $`\alpha =0`$, mass loss rate per unit area is uniform, while for $`\alpha =1`$ mass loss is proportional to the luminous flux of the disk.

Like SV93, KWD95 adopt a velocity law reminiscent of those used to describe stellar winds, but modified slightly to reflect the fact the wind is arising from the disk. More specifically, they define the poloidal velocity

$$v_l=c_s+(v_{\mathrm{}}c_s)\left(1\frac{R_v}{l+R_v}\right)^\beta $$
(8)
where $`c_s`$ is the speed at the base of the flow, $`R_v`$ is the scale length for velocity law, and $`\beta `$ is the exponent that determines the number of scale lengths over which the wind accelerates. Like SV93, $`v_{\mathrm{}}`$ is defined as a multiple of the escape velocity at the footpoint of each stream line.

KWD95 assume the poloidal velocity at the base of the wind is just the sound speed of the disk surface, defined as

$$c_s(R)=10\sqrt{\frac{T_{eff}(R)}{10^4K}}\mathrm{km}\mathrm{s}^1.$$
(9)
where $`T_{eff}`$ is given by equation 10 below (Frank, King & Raine, 1985).

Thus for the KWD95 prescription, one must define the basic parameters of the system as was the case for SV93, and then $`\dot{m}_{wind}`$, D, $`\alpha `$, $`R_s`$, $`\beta `$, and the multiple of the escape velocity to be used for $`v_{\mathrm{}}`$. The KWD95 prescription has fewer parameters, mostly due to a less flexible prescription for controlling the direction of streamlines and for limiting the region on the disk from which the wind emerges. Depending on one’s perspective, the differences in flexibility between the two prescriptions can be taken as an advantage or a disadvantage. In practice, the differences between them are relatively modest, and one can generally recast a wind described in terms of one description into the other.

### 2.3 Radiation Sources

The accretion disk, the WD, a boundary layer and the wind are all sources of radiation in the code. The disk and the WD are treated as opaque. Irradiation of the disk and WD by each other and/or by the wind (“backwarming”) is not included in the current version of the code. The wind itself emits both radiation due resonantly scattered photons and radiation due to the fact that the wind has been heated by other sources of emission.

The companion star is not included as a radiation source. Including the companion star as a radiation source would require one to adopt a full 3-d grid for the wind, since in that case the radiation field and consequently the physical conditions in the wind would vary with azimuth. The companion is included as a dark absorber (located outside the wind) so that the code can be used to model line profiles through an eclipse.

The disk, with an inner radius of $`R_{wd}`$ and outer radius $`R_D`$, is assumed to be geometrically thin and optically thick and to follow the standard radial effective temperature distribution (see, e.g. Wade, 1984).

$$T_{eff}(R)=T_{}\left(\frac{R_{wd}}{R}\right)^{3/4}\left(1\sqrt{\frac{R_{wd}}{R}}\right)^{1/4},$$
(10)
where $`R_{wd}`$ is the WD radius. $`T_{}`$ is a function of the mass accretion rate, $`\dot{M}_{acc}`$, and the WD mass and radius, or more specifically,

$$T_{}=\left(\frac{3GM_{wd}\dot{M}_{acc}}{8\sigma \pi R_{wd}^3}\right)^{1/4},$$
(11)
with $`\sigma `$ being the Stefan-Boltzmann constant. The disk emits as either an ensemble of blackbodies (BBs) or an ensemble of Kurucz (1993) and/or Hubeny (1988) model stellar spectra (SAs; see, e.g. Wade, 1984; Long et al., 1994). When assumed to emit as a stellar atmosphere, the prescription for the effective gravity suggested by Herter et al. (1979) is utilized:

$$g(r)=g_o\left(\frac{r}{R_{wd}}\right)^{15/8}\left[1\left(\frac{r}{R_{wd}}\right)^{1/2}\right]^{1/8}$$
(12)
where

$$g_o=5.96\times 10^5\left(\frac{M_{wd}}{M_{\mathrm{}}}\right)^{5/8}\left(\frac{R_{wd}}{10^9cm}\right)^{15/8}\left(\frac{\dot{M}}{10^{16}gs^1}\right)^{1/8}cms^1.$$
(13)
The broadening effect of the disk’s Keplerian rotation on the emitted spectrum is properly included.

The WD is taken to be spherical and to have a temperature of $`T_{wd}`$. Its total luminosity is $`L_{wd}=4\pi \sigma R_{wd}^2T_{wd}^4`$. Its spectrum can be generated as a blackbody, or, more commonly, from a stellar atmosphere.

Radiation due to the boundary layer (BL) between the disk and the WD is included in a simple fashion, and is specified in terms of a BL temperature, $`T_{BL}`$, and luminosity, $`L_{BL}`$. The size and extent of the boundary layer is not well known on observational or theoretical grounds. As a first approximation, we have assumed the BL to radiate from the entire surface of the WD. Usually, the BL is modeled as a BB, although it, too, can be described by a stellar model atmosphere, provided its temperature is low enough.

The wind both absorbs and radiates in our code, as will be discussed more fully in Section 2.5. Each wind cell emits free-free, free-bound and bound-bound radiation, in accordance with its current ionization state and temperature.

All photons (energy packets) in the code are assigned a direction ($`\theta ,\varphi `$), weight ($`w_p`$) and frequency ($`\nu `$) at the time of emission. Since the disk, WD, and BL all are assumed to be optically thick, they emit photons into the “outward” hemisphere only. Correspondingly, the probability, $`p(\theta ,\varphi )`$, that one of these photons is emitted in a direction ($`\theta ,\varphi )`$) relative to the local outward normal is proportional to

$$p(\theta ,\varphi )d\theta d\varphi =\eta (\theta )\mathrm{sin}\theta \mathrm{cos}\theta d\theta \times d\varphi $$
(14)
where $`\eta (\theta )`$ is the limb darkening. For “solid” surfaces, we assume linear limb darkening, viz.

$$\eta (\theta )=a(1+b\mathrm{cos}\theta ).$$
(15)
Here, the parameter $`b`$ defines the strength of the limb-darkening, which we take in the usual Eddington approximation to be 3/2. The parameter $`a`$ is a normalization constant that can be included either explicitly or implicitly. (The only important point in the MC approach is that the total probability must sum to unity.) The Eddington approximation is known to provide a fairly good description of limb darkening on the sun (see, e.g. Mihalas, 1978). The same approximation was used by SV93 and KWD95 in their CV wind codes, and by us (see, e.g. Knigge et al., 1997) in fitting the disk continua of high state CVs.

Each photon bundle is initially given a weight ($`w_p`$) that corresponds to a fraction of the physical luminosity, or, more usually, the band-limited luminosity (within wavelength or frequency limits) of the system. As a result, both the normalization and spectral shape of models generated with the code have meaning. In the simplest case when all photon bundles have the same weight, $`w_p=L_{tot}^{}/N_{tot}`$, where $`L_{tot}^{}`$ and $`N_{tot}`$ are the sum of luminosities and number of photons emitted by the disk, WD, boundary layer and wind in one complete iteration. The number of disk, WD, BL, and wind photons emitted in one complete iteration is given by $`N_{disk}=N_{tot}L_{disk}/L_{tot}^{}`$, $`N_{wd}=N_{tot}L_{wd}/L_{tot}^{}`$, $`N_{BL}=N_{tot}L_{BL}/L_{tot}^{}`$ and $`N_{wind}=N_{tot}L_{wind}/L_{tot}^{}`$, respectively. Note that $`L_{tot}^{}`$ as defined is not the luminosity that would be measured by external observers, but rather the sum of the luminosities of the individual radiation sources. Some of the the energy contained in these photons will be absorbed by the wind. Indeed, the wind must (and does ignoring a small amount of adiabatic cooling) absorb as much energy as it emits. As one would hope, the actual luminosity of the system as calculated by summing the weights of all photons that emerge from the system is quite close to the sum of the $`L_{disk}`$, $`L_{wd}`$, and $`L_{BL}`$. The only energy that is not properly accounted for currently in the code is due to the small fraction of photons that are “lost” due to impact on one of the “solid” surfaces before emerging from the system.

In order to ensure that the temperature and ionization structure of the wind converges, a statistically significant number of ionizing photons must pass through each grid cell during the ionization/thermal balance stage of the MC calculation (see section 2.4 below). In the absence of a luminous boundary layer, the fraction of the luminosity produced in photons with energies in energies above the He II edge at 54.4 eV is small, since the effective temperature of a typical CV disk is of order 50,000 K. For this reason, we split frequency space into multiple bands, and require that the numbers of photon bundles generated in each band exceed certain minimum fractions. For the calculations discussed here for example we have four bands, with breaks at 13.6, 24.6 and 54.4 eV corresponding to the ionization potentials of H I, He I and He II, and we require that at least 10% of the photon bundles are in the three higher energy bands. This requires a straightforward adjustment of the photon weights in each band, with the largest change in the high energy band where the number of photons was increased over that which would be produced via uniform sampling.

Cumulative distribution functions (CDFs), defined by

$$P(x)=_{x_{min}}^xp(x^{})𝑑x^{},$$
(16)
where $`p(x^{})`$ is the probability density for a process, are used extensively in generating and following photons within the code. CDFs allow one to decide the source of the photon, the location of the photon on that source, the direction of the photon, and the frequency of the photon, where the photon is scattered, etc. For example, for generating the frequency of a photon packet from a specific SED, the CDF is simply

$$P(\nu )=\frac{_{\nu _{min}}^\nu L_\nu (\nu ^{})𝑑\nu ^{}}{_{\nu _{min}}^{\nu _{max}}L_\nu (\nu ^{})𝑑\nu ^{}}.$$
(17)
By definition, all CDFs increase monotonically from zero to one and can thus be easily inverted to yield appropriately distributed random frequencies. In our code, we create an array of dimension $`n`$ that contains the frequencies at which $`P(\nu )`$ reaches values of $`1/n`$. We then generate a random number $`p`$ from a uniform distribution from 0 to 1, locate the array elements that bracket $`p`$, and then linearly interpolate between the associated frequencies to get the “exact” frequency of the photon to be generated. It is also straightforward to extend this approach to assure that certain specific frequencies, e.g. photoionization edges, are included in such arrays. The exact way in which CDFs are implemented depends on the process in question. In some cases, e.g. for determining the optical depth of scattering, the integral can be calculated in closed form. In others, such as in generating photons from a BB, the CDF is constructed numerically (in units of $`h\nu /kT`$) and used throughout the subsequent calculation.

Wind photons are the most computationally intensive, because a small number of photons is generated in each cell and because the emission processes, especially line emission, are particularly complex. Our approach is to store the bound-free, free-free, and bound-bound luminosities within the appropriate frequency limits of each cell each time we calculate a new ionization structure for the wind. In addition, we store a very coarse CDF for the bound-bound emission, namely the frequencies at which the CDF reaches 0.1, 0.2, etc. We also store the band-limited luminosities for each separate recombination process. When a wind photon must be generated, based on the fraction of the $`L_{tot}`$ comprised by the wind, we first select the cell and type – bound-bound, bound-free, free-free – of photon based on the corresponding band-limited luminosities of each process in the individual cell. If a free-free or free-bound photon must be generated, we use the process described by equation 17. For line photons, where linear interpolation is clearly inappropriate, the analogous sum over a frequency ordered list of lines is

$$P(\nu _i)=\frac{_{i^{}=o,i}L_{\nu _i^{}}}{L_{lines}},$$
(18)
where here $`L_{\nu _i^{}}`$ refers to the luminosity of a specific line within a specific cell in the grid, and $`L_{lines}`$ is the band-limited luminosity due to lines in that same cell.

### 2.4 Ionization and Thermal Balance

Given the appropriate atomic data, an arbitrary set of elements can be included in our MC calculations (see Section 2.7). Ionization balance in the code is carried out approximately, usually using the formalism developed by Mazzali & Lucy (1993). This method is basically an improvement on the standard nebular approximation (e.g. Mihalas, 1978, Eqn 5-46).

Specifically, Mazzali & Lucy (1993) suggest the ionization formula

$$\frac{n_{j+1}n_e}{n_j}=W[\xi +W(1\xi )]\left(\frac{T_e}{T_R}\right)^{1/2}\left(\frac{n_{j+1}n_e}{n_j}\right)_{T_R}^{},$$
(19)
which, in principle, accounts for ionizations from and recombinations to all levels of an ion. In this equation, the last (“starred”) term on the right hand side refers to the Saha abundances, $`W`$ is an effective dilution factor, $`\xi `$ is the fraction of recombinations going directly to the ground state, and $`T_R`$ and $`T_e`$ are the radiation and electron temperature respectively. 333For diagnostic purposes, the code also implements several other approaches to ionization structure of the wind, including an option to fix the abundances of individual ions and an option to force LTE ionization populations.

Following Mazzali & Lucy (1993), we define the radiation temperature $`T_R`$ by computing the specific intensity-weighted ($`J_\nu `$-weighted) mean photon frequency, $`\overline{\nu }`$, for each cell in each iteration. For a black-body radiation field described by temperature $`T_R`$, this mean frequency satisfies the relation

$$h\overline{\nu }=3.832kT_R,$$
(20)
and so

$$T_R=\frac{h\overline{\nu }}{3.832k}.$$
(21)
Both geometric dilution and absorption will cause the local intensity, $`J`$, to have a value smaller than the black-body intensity at the local radiation temperature, $`B(T_R)`$. Here, we have dropped the usual subscript $`\nu `$ on $`J`$ and $`B`$ to denote that we are concerned with frequency integrated quantities here. We therefore define the effective dilution factor $`W`$ by demanding that

$$J=WB(T_R)=W(\sigma /\pi )T_R^4.$$
(22)
Note that in our MC approach, the frequency integrated intensity contributed by a single photon (energy packet) to the radiation field in a given cell, $`\mathrm{\Delta }J`$, is given by

$$4\pi \mathrm{\Delta }J=w_pV^1l.$$
(23)
where $`w_p`$ is the average weight of the photon during its passage through the cell, $`l`$ is the distance the photon travels through the cell, and $`V`$ is the cell volume.

The quantity $`\xi `$ in equation 19 is the fraction of recombinations going directly to the ground state. In principle, all recombination processes should be accounted for in calculating $`\xi `$. In practice, we estimate $`\xi `$ by using only radiative rates calculated in the hydrogenic approximation. More specifically, we use the analytic fits of Pequignot, Petitjean, & Boisson (1991), which in turn reflect the results of Martin (1988). Physically, the first appearance of $`\xi `$ in the square brackets in equation 19 accounts for recombinations to levels other than the ground state, whereas the second term involving $`\xi `$ accounts for photoionizations from excited levels.

There is one important difference between the calculation of the ionization balance in our code and the prescription of Mazzali & Lucy (1993). They use the ad hoc approximation $`T_e=0.9T_R`$ to estimate the local electron temperature. We calculate $`T_e`$ self-consistently for each cell by considering the relevant heating and cooling processes explicitly, as discussed below. It is worth emphasizing in this context that the thermal and ionization equilibria are inextricably interlinked: the strength of the heating and cooling processes (and hence $`T_e`$) depend on the ionization state, but the ionization state itself also depends on the local electron temperature. Indeed, it is this non-linear inter-dependence that requires the use of iterative techniques in radiative transfer codes of all kinds, including our own MC code. The Mazzali & Lucy (1993) approximation for $`T_e`$ reduces the complexity of this inter-dependence, but (quite intentionally) does not remove it completely: in their formalism, $`T_e`$ depends on $`T_R`$, which depends on the ionization state, which in turn depends on $`T_e`$.

The ionization and temperature structure of the wind is calculated iteratively in an initial series of Monte Carlo calculations. A guess at the radiation temperature and electron temperature is made, and this is used to calculate a trial ionization structure. A family of typically 100,000 photons is generated from all emission sources (disk, star, boundary layer, and wind) and propagated through the wind. Next, the radiation temperature in each cell is calculated, and the luminosity emitted by each cell is compared to the luminosity absorbed by that cell. The current estimate of the electron temperature is then adjusted, so as to reduce the difference between emitted and absorbed luminosities. Finally the ionization structure is updated, using this new estimate of $`T_e`$. At this point, a new iteration can begin. The whole process is repeated until both thermal and ionization structures have converged (in the sense described below).

In the process of calculating the ionization structure, we limit the maximum change in the electron temperature of each cell in one iteration. This is necessary not only because the interaction of the heating of a cell is strongly dependent on the ionization concentrations, which in turn reflects the electron temperature, but also because of statistical fluctuations in the number and frequency of photons passing through a given cell. The number of photons that pass through a cell depends on the effective solid angle of that cell, which varies, and can be relatively small for cells lie along the cones that define the edges of the wind. Therefore, we adaptively limit the fractional electron temperature change in a given cycle based on the history of the cell. The maximal temperature change allowed increases if the desired temperature adjustment is the same direction as a previous iteration or if the amplitude of the change is decreasing. Cells of this type are said to be “converging”. Conversely, if the desired temperature adjustment has changed sign from that of the previous iteration and is larger in amplitude, the maximum allowed adjustment for the cell is decreased. Limiting the rate of change of $`T_e`$ help to stabilize the progression of a cell towards its equilibrium structure. It smooths out statistical fluctuations and prevents excessive changes in the ionization structure of major heating or cooling constituents.

The code does not contain a formal convergence criterion; instead we fix the number of ionization cycles at the outset and monitor the results. In most cases, of order 20 iterations of 100,000 photons each are required to reach the point at which essentially all of the cells are “converging”, and in which $`>90\%`$ of the cells show temperature fluctuations of less than 5% and have heating and cooling balanced within 10%. At this point, the average temperature variation of the entire wind will no longer show a definite trend in direction in subsequent iterations, nor will the magnitude of the adjustment exceed 100 K. If this is not the case, we have the ability to continue the ionization calculation from its current point (or, more commonly, to investigate why it is not converging).

This simplified calculation of the ionization equilibrium is clearly one of the limitations of our MC approach. We have therefore compared the ionization equilibrium calculated by our code to that predicted by cloudy (Ferland, 2000) for the simple case of an optically thin, solar abundance gas cloud with a H density of $`10^{12}\mathrm{cm}^3`$ irradiated by a black body radiation field with dilution factors W ranging from 1 to 10-6 and radiation temperatures $`T_r`$ ranging from 2,000 to 1000,000 K. Fig. Modeling the Spectral Signatures of Accretion Disk Winds: A New Monte Carlo Approach shows one typical comparison for W = 0.01. The overall shapes of the ionization curves, here illustrated by He, N, C, and O are similar. The temperatures at which an ionization state predominates are not precisely the same, but typically the discrepancy in temperature is less that 15%. In this particular comparison, C IV, N V, and O VI all peak at somewhat lower temperatures in our modified nebular approximation than in the cloudy calculation. On the whole, however, given the relative crudeness of our approximation, the degree of agreement is quite reassuring.

### 2.5 Radiative transfer in the wind

The fundamental equation of radiative transfer is:

$$\frac{dI_\nu }{ds}=\kappa _\nu I_\nu +ϵ_\nu .$$
(24)
Here $`dI_\nu /ds`$ is the rate of change of the specific intensity $`I_\nu `$ in the direction s, $`\kappa _\nu `$ is the opacity and $`ϵ_\nu `$ is the emissivity of the wind. In our MC code, the second term on the right hand side of equation 24 (the wind emissivity) is accounted for by allowing each wind cell to generate photons in accordance with the cell’s electron temperature and ionization state. The first term on the right hand side (absorption) is implemented by distinguishing between “pure absorption” and “scattering” processes. Free-free absorption and photoionization are treated as “pure absorption” processes. These processes continuously remove energy from a photon bundle as it makes its way through the wind, but they do not change the direction of the bundle. By contrast, electron and resonance line scattering are implemented as scattering processes. Photon bundles scatter once their opacity due to these processes exceeds the threshold defined in equation 1. When the bundle scatters, the direction, frequency and energy of the bundle all change discontinuously.

#### 2.5.1 Line scattering and emission

Line transfer in our MC code is calculated in the Sobolev approximation. In this approximation a photon travels until it encounters a resonant surface, where the photon frequency $`\nu =\nu _{lu}(1v/c)`$ where $`\nu _{lu}`$ is the frequency associated with a particular transition and $`v`$ is the velocity of the wind along the line of sight of the photon.

In the Sobolev limit, the optical depth through that surface, allowing for stimulated emission, is given by

$$\delta \tau =(n_l\frac{g_l}{g_u}n_u)\left(\frac{h\nu }{4\pi }B_{lu}\right)\frac{\lambda }{|dv/ds|}=(n_l\frac{g_l}{g_u}n_u)\left(\frac{\pi e^2}{mc}f_{lu}\right)\frac{\lambda }{|dv/ds|}$$
(25)
where $`n_l`$ and $`n_u`$ refer to the density of ions in the lower and upper state of a particular transition, and $`dv/ds`$ is the velocity gradient in the direction of photon travel (see, e.g. Rybicki and Hummer, 1978).

All levels are calculated in a two-level approximation. The code does not incorporate a complete model of the levels of each ion. Specifically, for the purposes of calculating the level populations, we use an on-the-spot approximation for each level, i.e.

$$\frac{n_u}{n_l}=\frac{c_{lu}+B_{lu}WB_\nu (T_r)}{c_{ul}+A_{ul}+B_{ul}WB_\nu (T_r)}$$
(26)
where $`A_{ul}`$, $`B_{ul}`$, and $`B_{lu}`$ are the Einstein coefficients associated with a particular transition, $`c_{lu}`$ and $`c_{ul}`$ are collisional excitation and de-excitation rates. In the large J limit, the usual situation near the disk, this becomes

$$\frac{n_u}{n_l}=\frac{g_u}{g_l}\frac{W}{e^{h\nu /kT_r}+W1}$$
(27)
The lower level population $`n_1`$ is calculated from the diluted black-body approximation described by equation 37.

When a photon bundle is scattered, it loses energy to the wind, and is re-emitted at the frequency of the resonance transition in the co-moving local frame. The optical depth through the resonant surface can be very large, and therefore in a real CV wind a photon will scatter multiple times before exiting the region where scattering is occurring. In addition, although individual atoms scatter isotropically (at least under the assumption of complete redistribution), the fact that the optical depth is not the same in all directions will cause photons to emerge from a resonant surface anisotropically. Here, we use an elementary escape probability formalism to estimate both how much energy is lost to electrons at each scattering site, and the anisotropy. Specifically, we assume that each time a photon scatters, collisional de-excitation reduces the weight of a photon bundle by a fraction

$$q=\frac{c_{ul}(1e^{h\nu /kT_e})}{A_{ul}+c_{ul}(1e^{h\nu /kT_e})}$$
(28)
where $`c_{ul}`$ is the collisional de-excitation rate for an ion in the upper to lower state. Following Rybicki and Hummer (1983), we assume that the local escape probability, the probability that an emitted photon will escape the region in a single scatter in a given direction $`\stackrel{}{n}`$, is given by

$$P_{esc}=\frac{1e^{\tau (\stackrel{}{n})}}{\tau (\stackrel{}{n})}$$
(29)
where $`\tau (\stackrel{}{n})`$ is the optical depth through the entire resonant surface in that direction.

With this definition, it is straightforward to show that the fraction of the energy of the photon bundle that will be absorbed at a single scattering site is

$$Q=\frac{q}{q+<P_{esc}>(1q)}$$
(30)
where $`<P_{esc}>`$ is the angle-averaged local escape probability. In terms of photons weights, this means simply that following a scatter

$$w_{p,out}=(1Q)w_{p,in}.$$
(31)
In addition to scattering line photons, the wind is allowed to emit wind photons via thermal processes. The effective emissivity of the plasma is, however, reduced from that expected for an optically thin plasma. For complete redistribution, the two-level source function is given by

$$S=(1q)<J>+qB(T_e)$$
(32)
where q is given by equation 28. By analogy, in extending this to multiple scattering within a single resonant region, it is clear that the two-level source function can be written

$$S=\frac{ϵ}{\kappa }=(1Q)J+QB(T_e)$$
(33)
where Q is defined by equation 30, and that the effective thermal emissivity is $`Q\kappa B(T_e)`$. Thus the effective thermal contribution due to lines is given by

$$4\pi ϵ_{bb,eff}=Qn_uA_{ul}h\nu $$
(34)
where the sum is over all line transitions. As is the case for scattered photons, thermally emitted line photons are emitted anisotropically in accordance with equation 29.

#### 2.5.2 Photoionization and Recombination

As a photon bundle travels through a cell it loses energy as a result of continuum processes. The final weight of the photon after accounting for continuum absorption is $`w_{p,f}=w_{p,i}e^{\tau _{cont}}`$, where $`\tau _{cont}`$ is the optical depth due to all continuum processes.

The free-bound optical depth experienced by an energy packet with frequency $`\nu >\nu _{thresh}`$ as it travels a distance $`l`$ through a grid cell is

$$\tau _{phot}=\kappa _{phot}l=n_{jk}\sigma (\nu )_{jk}l.$$
(35)
where $`\sigma (\nu )_{jk}`$ is the photoionization cross section and $`n_{jk}`$ is the number density of ions of a particular ion j in a particular level k. Of the energy absorbed by each bound-free transition, a fraction $`\nu _t/\nu `$ is spent overcoming the binding energy ($`h\nu _t`$) of the level of the ion, but the remaining part, $`(1\nu _t/\nu )`$, is taken up as kinetic energy by the freed electron and thus goes into heating the electron pool. Consequently, the thermal energy gained by a grid cell (per unit time and per transition) due to one photoionizing energy packet passing through it is

$$G_{phot}=w_{p,i}(1e^{\tau _{cont}})\frac{\tau _{phot}}{\tau _{cont}}(1\frac{\nu _t}{\nu }).$$
(36)
Photoionization can occur from ground or excited states. For photoionization calculations, $`n_{jk}`$ is calculated assuming

$$n_{jk}=n_jWg_ke^{E_k/kT_r}/z_j(T_r)$$
(37)
where W is the usual dilution factor for that cell, $`g_k`$ is the statistical weight of level k and $`z_j`$ is the partition function for the ion calculated similarly. This is of course not a proper non-LTE calculation of the level populations, but is intended as a first-order correction to LTE.

The free-bound emissivity to the jth level of the ith ionization state of an ion is calculated from detailed balance and given by:

$$4\pi ϵ(\nu )_{jk,fb}=h\nu \left(\frac{2\pi mk}{h^2}\right)^{3/2}\left(\frac{8\pi }{c^2}\right)\frac{g_{jk}}{g_eg_{j+1}}T_e^{3/2}\nu ^2\sigma (\nu )_{jk}e^{h(\nu \nu _t)/kT_e}$$
(38)
where $`g_{j+1}`$ is the ground state multiplicity of the recombining ion, and $`g_{jk}`$ is the multiplicity of the excited jth level of the ion after recombination. The total rate at which energy is lost from a grid cell due to radiative recombinations is therefore given by

$$L_{fb}=V\underset{ions}{}\underset{levels}{}_{\nu _{thresh}}^{\mathrm{}}ϵ_{fb}(\nu )𝑑\nu $$
(39)
where V is the volume associated with the cell.

#### 2.5.3 Free-free heating and cooling

The angle-integrated, free-free volume emissivity due to an ion j of charge $`Z`$ is

$$4\pi ϵ_{ff}(\nu )=n_en_j\frac{8}{3}\left(\frac{2\pi }{3}\right)^{1/2}\frac{Z^2e^6}{m^2c^3}\left(\frac{m}{kT_e}\right)^{1/2}g_{ff}e^{h\nu /kT_e}.$$
(40)
Consequently, the rate at which energy is lost via free-free radiation from a grid cell in the MC simulation is just

$$L_{ff}=V\underset{ions}{}_0^{\mathrm{}}ϵ_{ff}(\nu )𝑑\nu .$$
(41)
By Kirchhoff’s law, the free-free absorption coefficient is

$$\kappa _{ff}=\frac{ϵ_{ff}}{B_\nu (T_e)}.$$
(42)
Given this, the heating rate due to free-free absorption may be calculated analogously to that due to photoionization. Unlike the case of photoionization heating, all of this lost weight represents energy that is (temporarily) being added to the free electron pool. Thus the thermal energy (per unit time) gained by a grid cell due to free-free absorption as an energy packet travels through the cell is

$$G_{ff}=w_{p,i}(1e^{\tau _{cont}})\frac{\tau _{ff}}{\tau _{cont}}.$$
(43)
### 2.6 Extraction of photons at specific inclination angles

In determining the ionization structure of the wind we use a traditional “live-or-die” method. Photon bundles move through the wind until they escape the system (“live”) or until they encounter an opaque surface (“die”). In spherically symmetric situations, the live-or-die method is perfectly adequate for constructing detailed spectra as well since all escaping photons can be used. If the live-or-die method is applied to our situation however, one can use only those photon bundles that escape in a certain narrow inclination range to construct the spectrum for that inclination. One can use different photons to construct spectra at a variety of inclinations in a single calculation, but since most photons escape at intermediate inclination angles, the photon noise in the simulated spectra will vary significantly.

To address this problem and to allow the calculation of spectra in eclipsing systems, we use a variation of the “viewpoint” technique developed by Woods (1995) and described by KWD95. Specifically, we do create and follow photons throughout the grid in a standard live-or-die fashion. But we use these photons to identify where in the grid photons scatter (or in some sense to create a source function in the wind). At each place where a photon originates or scatters, we calculate the relative probability for a photon to be emitted or to scatter in the specific direction of the observer. This probability is used to increase or decrease the weight of a photon to be extracted in that direction. For isotropic scattering in the wind, the weight is unchanged. However, for resonance scattering, the weight is modified by $`P_{esc}/<P_{esc}>`$, the probability that the photon will scatter toward a specific observer divided by the angle-averaged escape probability. A similar correction is applied for photons from the disk or WD photosphere, which, as described earlier, emit anisotropically as a result of limb darkening.

This modification of the weight accounts for local effects in photon generation or scattering. To account for non-local effects, the optical depth depth $`\tau _{tot}`$ is then calculated by extracting the photon along the line of sight to the observer. Assuming it does not encounter a hard surface – the disk, the WD, or the secondary – the photon weight is then reduced by $`e^{\tau _{tot}}`$ and the photon is allowed to contribute to the spectrum. In this manner, we are able to create simulated spectra with similar amounts of photon noise regardless of the inclination angle. (It is similar and not identical because there are differences in the numbers of photons that hit a solid surface.) To verify our implementation of the viewpoint technique, we have conducted various tests with the code that confirm that spectra calculated using the live-or-die method and the viewpoint technique produce the same spectral shapes.

### 2.7 Atomic data

Our MC code incorporates a fairly flexible scheme of reading external data files to acquire atomic data on elements, ions, level configurations, photoionization cross-sections and oscillator strengths. For the sample models discussed in the following sections, we have incorporated atomic data for H, He, C, N, O, Ne, Na, Mg, Al, Si, S, Ar, Ca, Fe, and Ni, using the abundances specified by Verner, Barthel, & Tytler (1994). The basic data on ions for these elements, principally ionization potential and ground state multiplicities were obtained from Verner, Verner, & Ferland (1996). For H, He, C, N, and O, topbase (Cunto, Mendoza, Ochsenbein, & Zeippen, 1993) is the source of both ground and excited state photoionization cross-sections; for higher Z elements we have used analytic approximations for ground-state populations (only) due to Verner, Ferland, Korista, & Yakovlev (1996). For lines, we have used either the list of Verner, Verner, & Ferland (1996), which is a list of $``$5000 transitions that are ground-state connected, or a subset of the Kurucz & Bell (1995) line list. In the latter case, we abstract $``$55,000 lines from the entire line list based on the oscillator strength and approximate temperature where the associated levels are expected to be populated. We construct the associated levels for each ion from information contained in the line lists, using a method similar to that described by Lucy (1999b).

## 3 Comparison to KWD95 and SV93

As indicated earlier, the two codes most widely applied to CV spectra are those of KWD95 and SV93. Therefore it is important to make a connection to these two codes in similar situations.

To illustrate the line shapes produced by their code, KWD95 calculated spectra for two rotating bi-conical flows with the velocity law described by equation 8. The two models differed only in the degree of collimation. Basically, they assumed a 0.8$`M_{\mathrm{}}`$ WD, accreting at a rate of $`10^9M_{\mathrm{}}\mathrm{yr}^1`$ from a disk with a radius of 50 $`R_{wd}`$. They assumed that the wind had a mass loss rate of 1% of the accretion rate, a velocity scale length of 10 $`R_{wd}`$, and a terminal velocity 1.5$`\times `$ that of the local escape velocity. For the weakly collimated model, the took D to be 3.57 $`R_{wd}`$, which produced a wind with opening angles ranging from 16° to 86°; for the tightly collimated model they took D to be 57.12 $`R_{wd}`$, implying opening angles of 1° to 41°.444These correspond to models with $`d/d_{min}`$ of 1 and 16 in the terminology employed by KWD95. The complete set of parameters for this model are summarized in Table 1. As noted earlier, the code described by KWD95 was designed to model a few lines in detail. It is also a Monte Carlo code, but it does not use the Sobolev approximation and instead attempts to follow the radiative transfer exactly. It does not include an ionization calculation of the wind, and so the ionization state must be specified. For the fiducial models described by KWD95, C IV was assumed to be the dominant ionization state of C everywhere. In addition, thermal emission from the wind was not included.

In an attempt to make a direct comparison between our new code and the code developed by KWD95, we have calculated model spectra for the parameters used by KWD95 to generate their Fig. 9. For this calculation, we fixed the ionization structure and temperature of the wind, and suppressed wind emission. We used a 100 x 100 grid555It is perhaps obvious, but each grid cell describes cylindrical region above and below the disk, since we assume the wind model is symmetric through the disk plane. and a flight of 106 photons to generate CIV line profiles for the same inclination angles shown by KWD95, i.e. 30, 60, 75 and 85°.

The results of this comparison are shown in Fig. Modeling the Spectral Signatures of Accretion Disk Winds: A New Monte Carlo Approach. The profiles calculated with our new code are rendered in black, while the profiles calculated with the current version of KWD95’s code are overplotted in grey. At low inclination angles, when the observer is looking through the oncoming wind at the disk, the profiles for the weak and strongly collimated winds exhibit strong absorption at wavelengths shorter than the rest wavelengths of 1548.2 and 1550.8 Å of the C IV doublet. There are also emission wings near the rest wavelengths arising from scattering in the wind. At higher inclinations the emission portion of the profile increases in importance in both wind models. This is due to the decreasing projected area of the disk and exacerbated by limb-darkening. The fact that the weakly collimated wind has broader, but less deep, absorption profiles at low inclination and a higher peak line-to-continuum ratio at large inclinations is due to the fact that it is a more efficient scatterer. It is more efficient, because velocity gradients in the dense portion of the wind are more varied, and the “size” of the wind is larger than in the denser, but tightly collimated model.

The overall structure of the profiles is very similar, with all of the same trends appearing in both sets of calculations. With the exception of the highest inclination profiles (i=85°), both codes produce similar line to continuum ratios, similar overall widths, and the same narrow absorption features. At 85°, the main differences are in the central portion of the profile where KWD95 produces more flux relative to the continuum than does our new code. We have carried out extensive tests to verify that there are no discrepancies in the velocity and density structures in the two codes. We can isolate unscattered and scattered emission in the two codes and we know that all of the profile differences arise in the scattered light. We suspect, but cannot prove, that the remaining differences are due to the simplifications inherent in the Sobolev approximation. We know, for example, that in our Sobolev code the scattered light profiles are sensitive to the anisotropy of the scattering. In principle at least, KWD95 should be expected to produce a more exact version of the profile in multiple scattering situations.

We have also carried out a comparison to the fiducial model used by SV93, using the parameters listed in Table 2. Briefly, this model is for a system with a 40,000 K WD with a mass of 0.8 $`M_{\mathrm{}}`$ and radius of $`6\times 10^8\mathrm{cm}`$. The disk extends from the WD surface to 34 $`R_{wd}`$ and is presumed to be in a steady-state with an $`\dot{m}_{disk}`$ of $`10^8M_{\mathrm{}}\mathrm{yr}^1`$. The wind is constrained to emanate from the disk in a broad annulus between 4 and 12 $`R_{wd}`$. At the inner wind edge, streamlines are inclined at 12° with respect to the disk normal, whereas at the outer wind edge, the streamline inclination is 65 °. The total mass loss rate in the wind is 10% of the accretion rate. The scale length for defining the velocity of the wind is 100 Rwd, so the wind does not reach terminal velocity until the wind is well outside the system. The terminal velocity along each stream line is 3 times the escape velocity at the footpoint of the wind, so the most rapidly moving portions of the wind are located at its inner edge.

For this particular calculation, we employed the 100 x 100 grid, and used 20 cycles of $`10^6\mathrm{p}`$hotons to calculate the ionization structure; at this point, 96% of the cells were converged by our criteria. Consistent with SV93, all of the sources of continuum radiation were simulated in terms of appropriately-weighted black bodies. We then generated $`2\times 10^6`$ photons to create simulated CIV line profiles at 10 degree intervals; these are plotted as black lines in Fig. Modeling the Spectral Signatures of Accretion Disk Winds: A New Monte Carlo Approach. For this calculation, C IV was represented as a singlet and emission from the wind was suppressed in an attempt to mimic SV93’s calculation as closely as possible. In this case, the line profiles at very low inclination show emission profiles. Unlike the KWD95 prescription, the SV93 fiducial model does not have a wind covering the innermost portion of the disk (where much of the UV emission arises) and as a result there is little or no absorption at very low inclination. When the inclination reaches 30°, and the observer is viewing the system from within the wind cone, absorption profiles develop. The blue edge of the absorption is sharp, but moves redward as inclination increases. This is due to the fact that the fastest stream lines, which are at the inner edge of the wind, intersect the line of sight with increasing obliquity. The emission wing grows increasingly strong with inclination for the same reasons as in the KWD95 model above. When the inclination reaches 70°, which is beyond $`\theta _{max}`$, the profile is completely in emission.

Also shown in grey in Fig. Modeling the Spectral Signatures of Accretion Disk Winds: A New Monte Carlo Approach are the profiles calculated by SV93 for this model, as digitized from their Fig. 5. As was the case for comparisons with KWD95, the same general features are seen in both calculations, and there is quite good agreement in the emission line portion of the profiles; this may reflect the fact that the SV93 escape probability assumptions are rather similar to our own. There is, however, a difference in the shapes of the absorption profiles at intermediate inclination angles (40-60°). The blue edges of the lines lie at about the same wavelengths, and the edges are sharp in both sets of models. However, the SV93 profiles exhibit maximum absorption at about 1540 Å, with the absorption depths declining towards the blue edge. These differences are almost surely due to the details of the ionization structure calculated in the two cases. The ionization fraction of C IV that we calculate in the inner 50 $`R_{wd}`$ portion of the wind is shown in Fig. Modeling the Spectral Signatures of Accretion Disk Winds: A New Monte Carlo Approach, and should be compared with Fig. 4 of SV93. However, in our calculation, the fractional abundance of CIV exceeds 0.5 over a significant portion of the wind at all cylindrical radii from the WD; by contrast, in the SV93 calculation, the CIV density declines significantly beyond 25 $`R_{wd}`$ in the radial direction and also beyond the inner 5 $`R_{wd}`$ above the disk. As a result, the optical depth at the blue edge of the C IV profile is much higher in our calculation. It is not clear at this stage which ionization structure is more correct. The SV93 ionization calculation is not described in detail and we have adopted a rather simplified approach to ionization ourselves. If our ionization calculation is more accurate, it suggests that somewhat lower mass loss rates will be required to model the line spectra of CVs.

Overall, we conclude that the agreement between our calculations and those of SV93 is reassuring, certainly for us, and possibly for SV93 as well. The differences in the profile shapes we calculate are small compared to those that will result from current uncertainties in the basic parameters of CV winds, such as $`\dot{m}_{wind}`$, the degree of collimation, and the velocity law.

## 4 A typical calculation

Our primary motivation in developing a new radiative transfer code was to enable us to simulate spectra over a wide wavelength range. Fig. Modeling the Spectral Signatures of Accretion Disk Winds: A New Monte Carlo Approach illustrates the results a simulation covering 1000 Å in the FUV waveband for the SV93 fiducial model view at an inclination of 62.5°. For this and other examples described below, $`10^6`$ photon bundles were created and followed through the wind to create the final spectrum. Since both disk and WD spectral distributions have been formed from appropriately weighted black bodies, all of the features in the spectrum are due to the effects of the wind, including the sharp cutoff at 912 Å. In this simulation, C IV and N V show well-developed P Cygni profiles. O VI is not evident. A number of lower ionization potential lines of ions, such as C III, N III, and Si IV, appear only in absorption. These ions are present only near the disk; as a result, we see them only in absorption along the line of sight to the disk. Ly$`\alpha `$ and Ly$`\beta `$ are also seen in absorption, even though the ionization fraction of neutral hydrogen is quite small.

In Fig. Modeling the Spectral Signatures of Accretion Disk Winds: A New Monte Carlo Approach, the same model is shown for inclinations ranging from 10° to 80°, but in this case the disk and WD are simulated in terms of appropriately weighted stellar atmospheres. At 10°, when the observer views the system from “inside” the inner wind cone, the effects of the wind on the spectrum are quite subtle, and most of the features are photospheric with the exception of the narrow emission profiles of C IV. At angles of 27.5, 45, and 62.5°, the observer is viewing the disk through larger and larger columns of wind material, so more and more ions become apparent. Often these ions deepen features in the photospheric spectrum, since, not surprisingly, the ionization state of the wind immediately above the disk is not that different from the photosphere. At 62.5°, the line of sight skims over the outer lower temperature portion of the disk, so the lower ionization lines appear in the spectrum. The profiles of some lines (e.g. C IV and N V) already exhibit strong emission at this inclination. This reflects that fact that the ionization stages corresponding to these lines are highly extended in the wind. As a result, relatively large portions of the line forming regions are not viewed projected against the bright disk, and these regions contribute only scattered emission. A comparison of the spectrum in Fig. Modeling the Spectral Signatures of Accretion Disk Winds: A New Monte Carlo Approach to the 62.5° panel in Fig. Modeling the Spectral Signatures of Accretion Disk Winds: A New Monte Carlo Approach shows how important an accurate “photospheric” model of the disk and WD remains for modeling the FUV spectra of CVs. Finally, at 80°, when the observer views the system from “outside” the wind cone, the spectrum is dominated by wind emission lines.

## 5 Modeling the Spectrum of a Real CV

In the absence of an agreed physical model of the wind of a DN in outburst, any description of a bi-conical flow necessarily has a large number of parameters. It is far from clear how to explore the parameter space, especially since the luminosity, geometry and SED of the boundary layer is unknown in almost all CVs.

However, an example of the type of “fit” that can be obtained is shown in Fig. Modeling the Spectral Signatures of Accretion Disk Winds: A New Monte Carlo Approach, which shows data obtained with HUT during a normal outburst of Z Cam compared to a model spectrum generated using a wind geometry nearly identical to that described by Knigge et al. (1997). Qualitatively, the model spectrum and the data resemble one another quite well. In overlaying the two spectra, we have assumed the value of E(B-V) and NH used by Knigge et al. (1997), and simply renormalized. The strong resonance lines of C IV, N V, and O VI have approximately the right profiles, and many of the features of the lower ionization lines can be seen in both the model and in the data. The specific model shown here has the same amount of collimation, the same acceleration length and terminal velocity, the same $`\dot{m}_{wind}`$ of $`1\times 10^{10}M_{\mathrm{}}\mathrm{yr}^1`$, and a similar value of $`6\times 10^9M_{\mathrm{}}`$ for $`\dot{m}_{disk}`$ (well within the range derived from fitting the continuum spectrum). The disk continuum was generated from an ensemble of stellar atmospheres, and the Kurucz & Bell (1995) line list was utilized in the wind. The model shown is for an inclination of 69°, which is near the maximum of range allowed for Z Cam. The main effect of reducing the inclination angle slightly would be to decrease the emission wing of C IV. The only real change from the model described by Knigge et al. (1997) is that the model shown here is has the parameter $`\alpha `$ controlling the mass loss rate as a function of radius is 1 instead of 0.4. In our model the mass loss rate per unit area is proportional to the surface flux, which concentrates the wind closer to the WD. This helped to reduce the strength of several low ionization lines, which otherwise produced narrow deep absorption lines not seen in the data.

As noted above, we did not systematically search for a best fit to the Z Cam spectrum. Therefore, given our limited experience at this point, we cannot rule out that different wind structures might produce similar emergent spectra. We did, however, generate other models around the one described here. Small changes, in for example the collimation of the wind or the acceleration length, do improve individual lines, but do not qualitatively change one’s impression of the fit overall. One interesting fact is that it is very difficult to fit both the region below 1000 Å and the resonance lines when model atmospheres are used to generate the disk spectrum. This is because both the wind and the the turnover in the underlying disk continuum depress the overall flux in this region. This is clearly an area we will have to explore in our attempts to fit FUSE spectra in the future.

## 6 Summary

The inherent advantage of a Monte Carlo approach in simulating the spectra of cataclysmic variables is the ease with which one can model complex geometries. In the Monte Carlo approach, axially symmetric winds are not significantly more difficult to model than spherical winds. In the case of our code, we have now progressed to the point of spectral verisimilitude, where it is difficult to distinguish a simulation from data. The wind prescriptions required to produce verisimilitude are not that different from those which were developed to model C IV line profiles.

On the other hand, preliminary comparisons of the models and data suggest that considerable work will be required to reproduce the spectra of some dwarf novae. Nevertheless, we are already having some successes, such as the model shown here for Z Cam and our model of the FUSE spectrum of SS Cygni (Froning et al. 2001). We are now embarking upon an effort to fully explore the parameter space inherent in the existing models and to compare these models to a variety of HUT, FUSE and HST spectra of high state CVs. This should allow us to determine whether there is a single class of models, broadly or narrowly collimated for example, that can be used to approximate the wind geometry in the majority of CVs.

This work was supported by NASA through grants associated with our analysis of HST and FUSE spectra of dwarf novae. Specifically we acknowledge support from grants G0-7362 and GO-8279 from the Space Telescope Science Institute, which is operated by AURA, Inc., under NASA contract NAS5-26555, as well from grants NAG5-9283 and NAG5-10338, which are associated with our FUSE projects. We have also benefitted from many useful conversations from our colleagues, most notably Janet Drew, Ivan Hubeny, and Tom Marsh. This work was completed during a delightful visit by one of us (ksl) to the Department of Physics & Astronomy of the University of Southampton.


text:
# A semi-quantum version of the game of Life

## I Introduction

John Conway’s game of Life gar71  is a well known two dimensional cellular automaton where cells are arranged in a square grid and have binary values generally known as dead or alive. The status of the cells change in a discrete fashion each “generation” depending upon the number of neighboring cells that are alive, the general idea being that a cell dies if there is either overcrowding or isolation. There are many different rules that can be applied for birth or survival of a cell and a number of these give rise to interesting properties such as still lives (stable patterns), oscillators (patterns that periodically repeat), spaceships or gliders (fixed shapes that move across the Life universe), glider guns, and so on gar71 ; ber82 ; gar83 . Conway’s original rules are one of the few that are balanced between survival and extinction of the life “organisms.” In this version a dead (or empty) cell becomes alive if it has exactly three living neighbors, while an alive cell survives if and only if it has two or three living neighbors. Much literature on the game of Life and its implications exists. For recent discussion on the possibilities of this and other cellular automata the interested reader is referred to wol02 .

The recent interest in quantum games mey99 ; eis99 ; mar00 ; ben00a ; ben00b ; john01 ; du01 ; iqbal01 ; iqbal02 ; flitney02a ; flitney02b  suggests the possibility of applying the idea of superposition of states in quantum mechanics to the game of Life. Unfortunately Conway’s Life is irreversible while, in the absence of a measurement, quantum mechanics is reversible. In particular, operators that represent measurable quantities must be unitary. A full quantum Life would be problematic given the known difficulties of quantum cellular automata mey96 . Interesting behavior can still be obtained in a semi-quantum mechanical Life by representing the cells by classical sine-wave oscillators with a period equal to one generation, an amplitude between zero and one, and a variable phase. The amplitude of the oscillation represents the coefficient of the alive state so that the square of the amplitude gives the probability of finding the cell in the alive state when a measurement of the “health” of the cell is taken. If the initial state of the system contains at least one cell that is in a superposition of eigenstates the neighboring cells will be influenced according to the coefficients of the respective eigenstates, propagating the superposition to the surrounding region.

If the coefficients of the superpositions are restricted to positive real numbers we simply obtain a probabilistic Life universe (i.e. a classical mixture of states) where we do not expect to see qualitatively new phenomena. By allowing the coefficients to be complex, that is, by allowing phase differences between the oscillators, qualitatively new phenomena, for example interference effects, may arise.

## II Classical Mixture of States

To represent the state of a cell we introduce the following notation:

$$|\psi =a|\text{alive}+b|\text{dead},$$
(1)
subject to the normalization condition

$$|a|^2+|b|^2=1.$$
(2)
$`|a|^2`$ and $`|b|^2`$ represent the probabilities of measuring the cell as alive or dead, respectively. If the values of $`a`$ and $`b`$ are restricted to non-negative real numbers the cell is a classical mixture of the alive and dead states. In our model $`|a|`$ is the amplitude of the oscillator. Restricting $`a`$ to non-negative real numbers corresponds to the oscillators all being in phase.

The birth, death and survival operators have the following effects

$`B|\psi `$ $`=`$ $`|\text{alive}`$ (3)
$`D|\psi `$ $`=`$ $`|\text{dead}`$
$`S|\psi `$ $`=`$ $`|\psi .`$
A cell can be represented by the vector

$$\left(\begin{array}{c}a\\ b\end{array}\right).$$

The $`B`$ and $`D`$ operators are not unitary. Indeed they can be represented in matrix form as

$`B`$ $`=`$ $`\left(\begin{array}{cc}1& 1\\ 0& 0\end{array}\right),`$ (6)
$`D`$ $`=`$ $`\left(\begin{array}{cc}0& 0\\ 1& 1\end{array}\right).`$ (9)
A new generation is obtained by determining the number of living neighbors each cell has and then applying the appropriate operator to that cell. The number of living neighbors in our model is the amplitude of the superposition of the oscillators representing the surrounding eight cells. This process is carried out on all cells effectively simultaneously. When the cells are permitted to take a superposition of states, the number of living neighbors need not be an integer. Thus a mixture of the $`B`$, $`D`$ and $`S`$ operators may need to be applied. For consistency with standard Life the following conditions will be imposed upon the operators that produce the next generation:

* If there are an integer number of living neighbors the operator applied must be the same as that in standard Life.
* The operator that is applied to a cell must continuously change from one of the basic forms to another as the sum of the $`a`$ coefficients from the neighboring cells changes from one integer to another.
* The operators can only depend upon this sum and not on the individual coefficients.

If the sum of the $`a`$ coefficients of the surrounding eight cells is

$$A=\underset{i=1}{\overset{8}{}}a_i$$
(10)
then the following set of operators, depending upon the value of $`A`$, is the simplest that has the required properties

$`0A1;`$ $`G_0=`$ $`D,`$ (11)
$`1<A2;`$ $`G_1=`$ $`(\sqrt{2}+1)(2A)D+(A1)S,`$
$`2<A3;`$ $`G_2=`$ $`(\sqrt{2}+1)(3A)S+(A2)B,`$
$`3<A<4;`$ $`G_3=`$ $`(\sqrt{2}+1)(4A)B+(A3)D,`$
$`A4;`$ $`G_4=`$ $`D.`$
For integer values of $`A`$, the $`G`$ operators are the same as the basic operators of standard Life, as required. For non-integer values in the range $`(1,4)`$, the operators are a linear combination of the standard operators. The factors of $`\sqrt{2}+1`$ have been inserted to give more appropriate behavior in the middle of each range. For example, consider the case where $`A=3+1/\sqrt{2}`$, a value that may represent three neighboring cells that are alive and one the has a probability of one half of being alive. The operator in this case is

$$G=\frac{1}{\sqrt{2}}B+\frac{1}{\sqrt{2}}D.$$
(12)
Applying this to either a cell in the alive, $`\left(\begin{array}{cc}1& \\ 0& \end{array}\right)`$ or dead, $`\left(\begin{array}{cc}0& \\ 1& \end{array}\right)`$ states will produce the state

$$|\psi =\frac{1}{\sqrt{2}}|\text{alive}+\frac{1}{\sqrt{2}}|\text{dead}$$
(13)
which represents a cell with a 50% probability of being alive. That is, $`G`$ is an equal combination of the birth and death operators, as we might have expected given the possibility that $`A`$ represents an equal probability of three or four living neighbors. Of course the same value of $`A`$ may have been obtained by other combinations of neighbors that do not lie half way between three and four living neighbors, but one of our requirements is that the operators can only depend on the sum of the $`a`$ coefficients of the neighboring cells and not on how the sum was obtained.

The new state of a cell is obtained by calculating $`A`$, applying the matrix $`G`$ corresponding to the appropriate operator:

$$\left(\begin{array}{c}a^{}\\ b^{}\end{array}\right)=G\left(\begin{array}{c}a\\ b\end{array}\right),$$
(14)
and then normalizing the resulting state so that $`|a^{}|^2+|b^{}|^2=1`$. It is this process of normalization that means that multiplying the matrix by a constant has no effect. Hence, for example, $`G_2`$ for $`A=3`$ has the same effect as $`G_3`$ in the limit as $`A3`$, despite differing by the constant factor $`(\sqrt{2}+1)`$.

## III Semi-quantum Life

To get qualitatively different behavior from classical Life we need to introduce a phase associated with the coefficients, that is, a phase difference between the oscillators. We require the following features from this version of Life:

* It must smoothly approach the classical mixture of states as all the phases are taken to zero.
* Interference, that is the partial or complete cancellation between cells of different phases, must be possible.
* The overall phase of the Universe must not be measurable. That is, multiplying all cells by $`e^{i\varphi }`$ for some real $`\varphi `$ will have no measurable consequences.
* The symmetry between $`(B,|\text{alive})`$ and $`(D,|\text{dead})`$ that is a feature of the original game of Life should be retained. That is, if the state of all cells is reversed ($`|\text{alive}|\text{dead}`$) and the operation of the $`B`$ and $`D`$ operators is reversed the system will behave in the same manner.

In order to incorporate complex coefficients while keeping the above properties, the basic operators are modified in the following way:

$`B|\text{dead}`$ $`=`$ $`e^{i\varphi }|\text{alive},`$ (15)
$`B|\text{alive}`$ $`=`$ $`|\text{alive},`$
$`D|\text{alive}`$ $`=`$ $`e^{i\varphi }|\text{dead},`$
$`D|\text{dead}`$ $`=`$ $`|\text{dead},`$
$`S|\psi `$ $`=`$ $`|\psi ,`$
where the superposition of the surrounding oscillators is

$$\alpha =\underset{i=1}{\overset{8}{}}a_i=Ae^{i\varphi },$$
(16)
$`A`$ and $`\varphi `$ being real positive numbers. That is, the birth and death operators are modified so that the new alive or dead state has the phase of the sum of the surrounding cells. The operation of the $`B`$ and $`D`$ operators on the state $`\left(\begin{array}{cc}a& \\ b& \end{array}\right)`$ can be written as

$`B\left(\begin{array}{cc}a& \\ b& \end{array}\right)`$ $`=`$ $`\left(\begin{array}{cc}a+|b|e^{i\varphi }& \\ 0& \end{array}\right),`$ (21)
$`D\left(\begin{array}{cc}a& \\ b& \end{array}\right)`$ $`=`$ $`\left(\begin{array}{cc}0& \\ |a|e^{i\varphi }+b& \end{array}\right),`$ (26)
with $`S`$ leaving the cell unchanged. The presence of the modulus signs means that the $`B`$ and $`D`$ operators are not linear. The modulus of the sum of the neighboring cells, $`A`$, determines which operators apply, in the same way as before (see Eqn. (11)). The addition of the phase factors for the cells allows for interference effects since the coefficients of alive cells may not always reinforce in taking the sum, $`\alpha =a_i`$. A cell with $`a=1`$ still has a unit probability of being measured in the alive state but its effect on the sum will cancel that of a cell with $`a=1`$. It will be noted that the phase of the dead state has no effects. It is retained in order to maintain the alive $``$ dead symmetry.

A useful notation to represent semi-quantum Life is to use an arrow whose length represents the amplitude of the $`a`$ coefficient and whose angle with the horizontal is a measure of the phase of $`a`$. That is, the arrow represents the phasor of the oscillator at the beginning of that generation. For example

$``$ $`=`$ $`\left(\begin{array}{c}1\\ 0\end{array}\right),`$ (29)
$``$ $`=`$ $`e^{i\pi /4}\left(\begin{array}{c}1/\sqrt{2}\\ 1/\sqrt{2}\end{array}\right)=\left(\begin{array}{c}(1+i)/\sqrt{2}\\ (1+i)/\sqrt{2}\end{array}\right),`$ (34)
etc. Then $`\alpha `$ is the vector sum of the arrows. This notation includes no information about the $`b`$ coefficient. The magnitude of this coefficient can be determined from $`a`$ and the normalization condition. As noted previously, the phase of the $`b`$ coefficient has no effect on the future progression of the game so it is not necessary to represent this.

## IV Results and Conclusions

The above rules have been implemented in Mathematica math . Some simple results are indicated in the attached figures. Figure 1 shows an elementary example of interference in a block still life with one cell out of phase with its neighbors. Figures 2 and 3 show the sort of structures, a wick and a loop that can be obtained from a line of cells of alternating phase, that is of units of $``$’s. The stability of the line of $``$’s results from the fact that while each cell in the line has exactly two living neighbors, the cells above or below this line have a net of zero (or one at a corner) living neighbors, due to the canceling effect of the opposite phases. No new births around the line will occur unlike the case where all the cells are in phase. The equivalent structure in Conway’s Life is a diagonal line of live cells, but in our case we are able to make the loop any shape. Figure 4 shows a stable boundary that results from the appropriate adjustment of the phase differences between the cells. The angles have been chosen so that each cell in the line has between two and three living neighbors, while the empty cells above and below the line have either two or four living neighbors and so remain life-less. Wicks and boundaries are known in standard Life but generally require more complex structures.

To summarize, John Conway’s game of Life is a two dimensional cellular automaton where the new state of a cell is determined by the sum of the neighboring states that are in one particular state generally referred to as “alive”. In semi-quantum Life cells may be in a superposition of the alive and dead states with the coefficient of the alive state being represented by an oscillator. The equivalent of evaluating the number of living neighbors of a cell is to take the superposition of the oscillators of the surrounding states. The amplitude of this superposition will determine which operator(s) to apply to the central cell to determine its new state, while the phase gives the phase of any new state produced. Such a system is able to reproduce some of the aspects of quantum mechanics such as interference.

Obviously this paper just touches on some of the results that can be obtained with this new scheme. More detailed possibilities will not emerge until a systematic search of the semi-quantum Life possibilities is undertaken but it can be seen that some new effects and structures occur and that some of the known effects in Conway’s Life can occur in a simpler manner.