text:
# Status of MiniBooNE

## 1 LSND EVIDENCE

The Liquid Scintillator Neutrino Detector ran at the Los Alamos Neutron Science Center over the years 1993 to 1998. Using the accelerator’s 800 MeV proton beam, the experiment studied neutrinos coming from the decay $`\mu ^+e^+\nu _e\overline{\nu }_\mu `$ of muons at rest and from the decay $`\pi ^+\mu ^+\nu _\mu `$ of pions in flight. The muon decay at rest (DAR) process allowed a search for $`\overline{\nu }_\mu \overline{\nu }_e`$ appearance via the reaction $`\overline{\nu }_epe^+n`$; neutrinos from pion decay in flight (DIF) provided a search for $`\nu _\mu \nu _e`$ appearance via the reaction $`\nu _eCe^{}N`$.

The LSND detector was a cylindrical tank of 167 tons of mineral oil doped with a small amount of scintillator viewed by 1220 photomultiplier tubes. It was positioned approximately 30 m from the beam target area. For the 1993-95 running period the primary proton target was a water target, and in 1996-98 the target was made mostly of tungsten. The neutrino yield from the water target was higher, so that about 59% of the total DAR flux and 62% of the total DIF flux came from the 1993-95 period.

The DAR process yields neutrinos with energies up to the 52.8 MeV decay endpoint, and the DIF process produced $`\nu _\mu `$ up to 300 MeV. The detector uses Cherenkov light for electron particle identification. The reaction $`\overline{\nu }_epe^+n`$ is observed by correlating the electron track with the 2.2 MeV gamma ray that follows from neutron capture on free protons in the liquid, $`npd\gamma `$.

LSND has presented results periodically. Analysis of 1993-95 DAR data produced an excess above backgrounds of $`51.0_{19.5}^{+20.2}\pm 8.0`$ events, corresponding to an oscillation probability of $`(0.31\pm 0.13)`$% (combined statistical and systematic errors) . The DIF data from the same period yielded an excess of $`18.1\pm 6.6\pm 4.0`$ events, and an oscillation probability of $`(0.26\pm 0.11)`$% .

In 2001, LSND presented final results using all of its data, combining the $`\overline{\nu }_\mu \overline{\nu }_e`$ and $`\nu _\mu \nu _e`$ searches into a single analysis, and employing new event reconstruction with better spatial resolution . A total excess of $`87.9\pm 22.4\pm 6.0`$ events was observed in the $`\overline{\nu }_\mu \overline{\nu }_e`$ search, corresponding to an oscillation probability of $`(0.264\pm 0.081)`$%. The common event selection was optimized for the DAR region below 60 MeV, and was therefore less effective than the previous DIF analysis in removing background events above 60 MeV. This analysis found no significant signal in the DIF energy region, where the observed excess was $`8.1\pm 12.2\pm 1.7`$ events.

The Karlsruhe Rutherford Medium Energy Neutrino experiment at the ISIS facility of the Rutherford Lab also searched for $`\overline{\nu }_\mu \overline{\nu }_e`$ oscillations using $`\overline{\nu }_\mu `$ from $`\mu ^+`$ decay at rest . The experiment employed a 56 ton segmented liquid scintillator detector located about 18 m from the beam target. The 1997-2001 run, known at KARMEN 2, did not observe an excess of events, finding 15 events while expecting a background of $`15.8\pm 0.5`$ events. However, KARMEN 2 only had the sensitivity to rule out a portion of the LSND allowed parameter space. A comparison of the favored LSND oscillation parameters with the KARMEN 2 limit is shown in Fig. 1. A combined statistical analysis of LSND and KARMEN 2 finds regions of oscillation parameters compatible with both experiments .

The LSND neutrino mass difference, $`\mathrm{\Delta }m_{LSND}^2>0.1`$ eV2, is distinct from those indicated by oscillations of solar and atmospheric neutrinos . This poses the problem that three neutrino masses cannot explain the three $`\mathrm{\Delta }m^2`$’s. Proposed solutions include consideration of a sterile neutrino, supersymmetry, $`CPT`$ violation, or lepton number violating muon decay . In addition, the relatively high neutrino mass would contribute to the question of dark matter. It remains important to confirm or refute the LSND evidence, which MiniBooNE sets out to do.

## 2 MINIBOONE OVERVIEW

At MiniBooNE we have begun collecting data to search for $`\nu _\mu \nu _e`$ oscillations. The $`\nu _\mu `$ beam energy is in the range 0.5 – 1.5 GeV with a small intrinsic $`\nu _e`$ component and we will search for an excess of electron neutrino events in a detector located approximately 500 m from the neutrino source. The baseline to neutrino energy ratio is thereby similar to that of LSND, $`L/E1`$ m/MeV, while neutrino energies are more than an order of magnitude higher.

The MiniBooNE neutrino beam is initiated by a primary beam of 8 GeV protons from the Fermilab Booster. The Booster is a reliable, high intensity machine, expected to run at least $`2\times 10^7`$ s per year, while delivering $`5\times 10^{12}`$ protons per 1.6 $`\mu `$s pulse at a rate of 5 Hz to MiniBooNE. By this proton accounting, one nominal year for the experiment is expected to supply $`5\times 10^{20}`$ protons on target. The Booster has the capacity to provide protons for several Fermilab efforts, and the Booster continues to deliver protons to the TeVatron collider program.

A secondary beam is produced when the 8 GeV protons strike a beryllium target positioned inside a magnetic horn. At present, positively charged particles (mostly $`\pi ^+`$’s) from the target are focused forward by the single horn into a 50 m decay tunnel. The $`\nu _\mu `$ beam is produced from the decay of these secondary particles. At a later time, the polarity of the horn can be changed in order to focus $`\pi ^{}`$’s to generate a $`\overline{\nu }_\mu `$ beam. Decay lengths of 50 and 25 m are possible through the use of two steel and concrete beam absorbers. One absorber is permanently positioned at the end of the decay tunnel. The intermediate absorber can either be lowered into the decay tunnel or raised out of the way. This ability to vary the decay length will provide a check of experimental systematics associated with $`\nu _e`$ contamination.

The MiniBooNE neutrino detector consists of 800 tons of pure mineral oil contained in a 40-foot (12.2 m) diameter spherical tank. A structure in the tank supports phototubes and optically isolates the most outer 35 cm of oil from the rest, turning the outer oil into a veto region that should stay quiet while a neutrino produces light only in the inner, main region. The main region is viewed by 1280 20-cm phototubes, providing 10% photocathode coverage of the 445 ton fiducial volume. The veto region contains 240 20-cm phototubes mounted in pairs on the tank wall. A schematic of the detector is shown in Fig. 2. The center of the detector is positioned about 490 m from the end of the decay tunnel, and about 6 m below ground, corresponding to the level of the neutrino beam. An overburden of about 3 m of soil shields the detector enclosure.

The detector records the arrival time and total charge for each hit phototube. From this information, the track position and direction are determined. Electrons from $`\nu _e`$ interactions are identified via their characteristic Cherenkov and scintillation light signatures (the undoped oil scintillates modestly). Backgrounds to the oscillation search will be due to $`\nu _e`$ contamination in the beam and to the misidentification as electrons of muons and $`\pi ^0`$’s produced in the detector. A display of one of MiniBooNE’s first candidate neutrino events is shown in Fig. 3.

## 3 BEAMLINE STATUS

Intense activity over summer 2002 focused on preparing the new 8 GeV proton beamline coming off of the Fermilab Booster to the MiniBooNE production target and secondary beam focusing horn. MiniBooNE started receiving proton beam on target in late August.

The focusing horn was designed to run at 5 Hz and for 200 Mpulses — a higher pulse rate and more pulses than any previous horn. In fall 2001 the horn was successfully bench-tested with more than 10 Mpulses making it already the world’s most pulsed horn.

In summer 2002 the proton beam was commissioned in a few stages. In June protons were sent through the new beamline to a stop ahead of the MiniBooNE target area. Next the target and horn were installed and then removed in a dry run of procedures to be used for removing a damaged, radioactive horn. In July, without the target and horn in place, protons were sent through the target location, where a temporary multiwire beam monitor recorded the beam profile. The multiwire measurements were used to commission the final focus and to begin calibration of permanent beam position monitors (BPM’s). In August the target and horn were reinstalled. The final configuration was modified slightly to include a multiwire just upstream of the target, which was used to complete the calibration of the target BPM’s.

High intensity Booster operation is restricted by radiation limits, and beam losses will have to be reduced before MiniBooNE can reach its full beam intensity. In September the Booster was typically supplying $`4\times 10^{12}`$ protons per pulse. (Losses were observed to grow disproportionately when trying to run with $`5\times 10^{12}`$ protons per pulse.) By mid-September the pulses were arriving on target at 1 Hz. The experimental priority is to increase the repetition rate to 5 Hz (whereas $`4\times 10^{12}`$ protons per pulse is adequate). Operation above 1 Hz will become possible with the planned installation of a new extraction magnet in the winter. However, under current run conditions 5 Hz operation would exceed the radiation limits for beamline elements.

The Booster now supplies protons to both the TeVatron collider and to MiniBooNE. The start of MiniBooNE running has had negligible effect on collider operation.

## 4 DETECTOR STATUS

Installation of phototubes and of all related in-tank calibration and monitoring hardware was completed in September 2001. Filling the tank with oil started in December and completed in May 2002. The data acquisition electronics were commissioned over the same period, and phototube signals produced by cosmic rays and by a laser calibration system indicated the rising oil level.

The laser calibration system provides pulses to four light-diffusing flasks at various locations in the detector. The system is being used to determine phototube time offsets and time slewing corrections. Energy calibrations are being performed using electrons from the decay of stopped cosmic ray muons. In addition, a cosmic ray calibration system is being commissioned, which consists of a hodoscope muon tracker above the detector with seven scintillator cubes (7.6 and 10 cm on a side) located under the tracker inside the detector. The system will provide the entering position and direction of cosmic ray muons. Cosmic rays that stop in a cube are a sample of muons with well-known path length in the detector, and will be used to calibrate the position, energy, and direction determination of the reconstruction algorithms.

## 5 PROSPECTS

As the first neutrinos were delivered to MiniBooNE the detector was ready to record them. It now remains to be seen how reliable and intense the beam will be. The beam intensity came up quickly, and within about three weeks of the start of the run, $`4\times 10^{12}`$ protons per pulse at a rate of 1 Hz were being delivered. The 1 Hz operation is currently limited by beamline hardware, which will be replaced this winter to allow 5 Hz operation. Work is ongoing to control beam losses so that 5 Hz beam operation is possible within radiation limits.

The experiment is currently running with a $`\nu _\mu `$ beam. If oscillations occur as indicated by LSND, with $`10^{21}`$ protons on target MiniBooNE will observe an excess of several hundred electron events that will be significantly above background. The option will exist during the run to reverse the horn polarity to provide a $`\overline{\nu }_\mu `$ beam. The $`\overline{\nu }_\mu `$ beam will have a reduced flux, and the antineutrino cross section is smaller, but the backgrounds are also expected to be smaller, so that the experiment would retain significant sensitivity. MiniBooNE’s expected 90% C.L. exclusion contours using a $`\nu _\mu `$ and a $`\overline{\nu }_\mu `$ beam are shown in Fig. 1.


format: text

text:
# Numerical simulations of mass outflows driven from accretion disks by radiation and magnetic forces.

## 1 Introduction

Powerful mass outflows from accretion disks are observed in many astrophysical environments such as active galactic nuclei (AGN); many types of interacting binary stars, e.g., non-magnetic cataclysmic variables (nMCVs); and young stellar objects (YSOs). Magnetic fields, the radiation force and thermal expansion have been suggested as mechanisms that can drive disk winds. These three mechanisms have been studied extensively using analytic as well as numerical methods. As a result of these studies, theoretical models have been developed that allow us to estimate under what physical conditions each of these mechanisms is efficient in launching, accelerating and collimating disk outflows. The notion that one universal physical mechanism can explain all outflows from accretion disks is very appealing. However, theoretical studies and observational results indicate that none of the three mechanisms alone is sufficient and probably such a single mechanism does not exist. Therefore it make sense to consider a hybrid model in which more than one mechanism is involved.

Magnetically driven winds from disks are the favored explanation for the outflows in many astrophysical environments. Blandford & Payne (1982) (see also Pelletier & Pudritz 1992) showed that the centrifugal force can drive a wind from the disk if the poloidal component of the magnetic field, $`𝐁_𝐩`$ makes an angle of $`>30^o`$ with respect to the normal to the disk surface. Generally, centrifugally-driven MHD disk winds (magnetocentrifugal winds for short) require the presence of a sufficiently strong, large-scale, ordered magnetic field threading the disk with a poloidal component at least comparable to the toroidal magnetic field, $`|B_\varphi /B_p|<1`$ (e.g., Cannizzo & Pudritz 1988, Pelletier & Pudritz 1992). Several groups have studied numerically axisymmetric outflows using the Blandford & Payne mechanism (e.g., Ustyugova et al. 1995, 1999; Romanova et al. 1997; Ouyed & Pudritz 1997a, 1997b, 1999; Krasnopolsky, Li & Blandford 1999; Kato, Kudoh & Shibata 2002). An important feature of magnetocentrifugal winds is that they require some assistance to flow freely and steadily from the surface of the disk, to pass through a slow magnetosonic surface (e.g., Blandford & Payne 1982). The numerical studies mentioned above do not resolve the vertical structure of the disk but treat it as a boundary surface through which mass is loaded on to the magnetic field lines at a specified rate.

Winds from disks can driven by the magnetic pressure. In particular, the toroidal magnetic field can quickly builds up due to the differential rotation of the disk so that $`|B_\varphi /B_p|1`$. In such a case, the magnetic pressure of the toroidal field can give rise to a self-starting wind (e.g., Uchida & Shibata 1985; Pudritz & Norman 1986; Shibata & Uchida 1986; Stone & Norman 1994; Contopoulos 1995; Kudoh & Shibata 1997; Ouyed & Pudritz 1997b). To produce a steady outflow driven by the magnetic pressure a steady supply of advected toroidal magnetic flux at the wind base is needed, otherwise the outflow is likely to be transient (e.g., K$`\ddot{\mathrm{o}}`$nigl 1993, Contopoulos 1995, Ouyed & Pudritz 1997b). It is still not clear whether the differential rotation of the disk can produce such a supply of the toroidal magnetic flux to match the escape of magnetic flux in the wind and even if it does whether such a system will be stable (e.g., Contopoulos 1995, Ouyed & Pudritz 1997b and references therein).

One of the reasons for favoring magnetic fields as an explanation for mass outflows from accretion disks is the fact that magnetic fields are very likely crucial for the existence of all accretion disks. The magnetorotational instability (MRI) (Balbus & Hawley 1991; and earlier by Velikov 1959 and Chandrasekhar 1960) has been shown to be a very robust and universal mechanism to produce turbulence and the transport of angular momentum in disks at all radii (Balbus & Hawley 1998). It is therefore likely that magnetic fields control mass accretion inside the disk and play a key role in producing a mass outflow from the disk. However, it has been demonstrated observationally and theoretically that accretion disks are capable of losing mass also via a radiation-driven wind, provided the disk luminousity in ultraviolet (UV) is high enough.

Radiation-driven disk winds have been extensively modeled recently (e.g., Pereyra, Kallman & Blondin 1997; Proga, Stone & Drew 1998, hereafter PSD 98; Proga 1999; Proga, Stone & Drew 1999, hereafter PSD 99; Feldmeier & Shlosman 1999; Feldmeier, Shlosman & Vitello 1999; Proga, Stone & Kallman 2000; Pereyra, Kallman & Blondin 2000; Proga & Kallman 2002). These studies showed that radiation pressure due to spectral lines can drive winds from luminous disks. This result has been expected (e.g., Vitello & Shlosman 1988). These studies, in particular those by PSD 98, also showed some unexpected results. For example, the flow is unsteady in cases where the disk luminosity dominates the driving radiation field. Despite the complex structure of the unsteady disk wind, the time-averaged mass loss rate and terminal velocity scale with luminosity, as do steady flows obtained where the radiation is dominated by the central object. In the most favorable conditions (i.e., high UV flux and low X-ray flux) the radiation force due to spectral lines (the line force) can exceed the radiation force due to electron scattering by a factor as high as $`2000`$ (e.g., Castor Abbott, & Klein 1995, hereafter CAK; Abbott 1982; Gayley 1995). Thus systems with UV luminosity, $`L_{UV}`$ as low as a few $`10^4`$ of their Eddington limit, $`L_{Edd}`$, can produce a powerful high velocity wind.

Generally, one can argue that in all accretion disks, with $`L_{UV}>`$ a few $`10^4L_{Edd}`$ mass outflows have been observed (Proga 2002). For example, accretion disks around: massive black holes, white dwarfs (as in AGN and nMCVs with $`\mathrm{\Gamma }_{UV}0.001`$) and low mass young stellar objects (as in FU Ori stars with $`\mathrm{\Gamma }_{UV}`$ a few $`\times `$ 0.01) show powerful fast winds. Systems that have too low UV luminosities to drive a wind include accretion disks around neutron stars and low mass black holes as in low mass X-ray binaries (LMXBs) and galactic black holes. These systems indeed do not show outflows similar to those observed in nMCVs, AGN and FU Ori. However, outflows, even in systems which appear to be luminous enough to produces radiation-driven disk winds, cannot be fully explained by just line driving.

For example, Drew & Proga (1999) applied results from pure line-driven (LD for short) disk wind models to nMCVs. In particular, they compared mass loss rates predicted by the models with observational constraints. Drew & Proga (1999) concluded that either mass accretion rates in high-state nMCVs are higher than presently thought by a factor of 2-3 or that radiation pressure alone is not quite sufficient to drive the observed hypersonic flows. The difficulty in accounting for the mass loss rate in a pure LD disk wind model for nMCVs is simply a reflection of the fact that the nMCV luminosities just barely satisfy the basic requirement, i.e., $`L_{UV}<7\times 10^4L_{Edd}`$. Synthetic line profiles computed based on pure LD wind models confirmed Drew & Proga’s conclusion (Proga et al. 2002). Our study has been partially motivated by this conclusion because if indeed radiation pressure alone does not suffice to drive the observed hypersonic flow then an obvious candidate to assist radiation pressure in these cases is MHD (e.g., Drew & Proga 1999).

In this paper, we study how magnetic fields can change disk winds driven by the line force for a given disk luminosity. We assume in our models that the transport of angular momentum in the disk is dominated by local disk viscosity, for instance due to the MRI in weakly magnetized disks (Balbus & Hawley 1998). Here we add magnetic fields to the PSD 99 model and solve self-consistently the full set of ideal MHD equations.

The outline of this paper is as follows. We describe our calculations in Section 2. We present our results and discuss their perceived limitations in Section 3. The paper ends in Section 4, with our conclusions.

## 2 Method

### 2.1 Equations and Numerical Techniques

To calculate the structure and evolution of a wind from a disk, we solve the equations of ideal MHD

$$\frac{D\rho }{Dt}+\rho 𝐯=0,$$
(1)
$$\rho \frac{D𝐯}{Dt}=𝐏\rho \mathrm{\Phi }+\frac{1}{4\pi }(\times 𝐁)\times 𝐁+\rho 𝐅^{rad}$$
(2)
$$\frac{𝐁}{t}=\times (𝐯\times 𝐁).$$
(3)
Here the convective derivative $`D/Dt`$ is equivalent to $`/t+𝐯`$. The dependent quantities $`\rho `$, $`𝐯`$, and $`P`$ are gas mass density, velocity, and scalar isotropic gas pressure, respectively, and $`𝐁`$ is the magnetic field. The gas in the wind is isothermal with a sound speed $`c_s`$. We calculate gravitational acceleration using the Newtonian potential $`\mathrm{\Phi }`$. The term $`𝐅^{rad}`$ in the equation of motion (2) is the total radiation force per unit mass. We solve these equations in spherical polar coordinates $`(r,\theta ,\varphi )`$, assuming axial symmetry about the rotational axis of the accretion disk ($`\theta =0^{}`$)

The geometry and assumptions needed to compute the radiation field from the disk and central object are as in PSD 99 (see also PSD 98). The disk is flat, Keplerian, geometrically-thin and optically-thick. We specify the radiation field of the disk by assuming that the local disk intensity follows the radial profile of the so-called $`\alpha `$-disk (Shakura & Sunyaev 1973), and therefore depends only on the mass accretion rate in the disk, $`\dot{M}_a`$, and the mass and radius of the central object, $`M_{}`$ and $`r_{}`$. In particular, the disk luminosity, $`L_DGM_{}\dot{M}_a/2r_{}`$. In models where the central object radiates, we take into account the irradiation of the disk, assuming that the disk re-emits all absorbed energy locally and isotropically. We express the central object luminosity $`L_{}`$ in units of the disk luminosity $`L_{}=xL_D`$.

Our numerical algorithm for evaluating the line force is described in PSD 99. Here we briefly describe the key elements of our calculations of the line force. We use the CAK force multiplier to calculate the line-driving force. In this approximation, a general form for this force at a point defined by the position vector $`𝐫`$ is

$$𝐅^{rad,l}(𝐫)=_\mathrm{\Omega }M(t)\left(\widehat{n}\frac{\sigma _eI(𝐫,\widehat{n})d\mathrm{\Omega }}{c}\right),$$
(4)
where $`I`$ is the frequency-integrated continuum intensity in the direction defined by the unit vector $`\widehat{n}`$, and $`\mathrm{\Omega }`$ is the solid angle subtended by the disk and central object at the point W. The term in brackets is the electron-scattering radiation force, where $`\sigma _e`$ is the mass-scattering coefficient for free electrons, and $`M(t)`$ is the force multiplier – the numerical factor which parameterizes by how much spectral lines increase the scattering coefficient. In the Sobolev approximation, $`M(t)`$ is a function of the optical depth parameter

$$t=\frac{\sigma _e\rho v_{th}}{\left|dv_l/dl\right|},$$
(5)
where $`v_{th}`$ is the thermal velocity and $`dv_l/dl`$ is the velocity gradient along $`\widehat{n}`$. The velocity gradient can be written as

$$\frac{dv_l}{dl}=Q\underset{i,j}{}\frac{1}{2}\left(\frac{\delta v_i}{\delta r_j}+\frac{\delta v_j}{\delta r_i}\right)n_in_j=\underset{i,j}{}e_{ij}n_in_j,$$
(6)
where $`e_{ij}`$ is the symmetric rate-of-strain tensor. Expressions for the components of $`e_{ij}`$ in spherical polar coordinates are given in Batchelor (1967).

We adopt the CAK analytical expression for the force multiplier as modified by Owocki, Castor & Rybicki (1988, see also PSD98)

$$M(t)=kt^\alpha \left[\frac{(1+\tau _{max})^{(1\alpha )}1}{\tau _{max}^{(1\alpha )}}\right]$$
(7)
where $`k`$ is proportional to the total number of lines, $`\alpha `$ is the ratio of optically thick to optically-thin lines, $`\tau _{max}=t\eta _{max}`$, and $`\eta _{max}`$ is a parameter related to the opacity of the most optically thick lines. The term in square brackets is the Owocki, Castor & Rybicki correction for the saturation of $`M(t)`$ as the wind becomes optically thin even in the strongest lines, i.e.,

$$\underset{\tau _{max}0}{lim}M(t)=M_{max}=k(1\alpha )\eta _{max}^\alpha .$$

We discretize the $`r\theta `$ domain into zones in our calculation of the wind structure. Our resolution in the $`r`$ and $`\theta `$ directions is sufficiently high to ensure that the subsonic portion of the model outflow is sampled by at least a few grid points in both $`r`$ and $`\theta `$. This requirement and the nature of the problem combine to demand an increasingly fine mesh toward the disk plane: here the density declines dramatically with height, and, moreover, the velocity in the wind increases rapidly. Our numerical resolution consists of 100 zones in each of the $`r`$ and $`\theta `$ directions, with fixed zone size ratios, $`dr_{k+1}/dr_k=d\theta _l/d\theta _{l+1}=1.05`$.

The boundary conditions for the hydrodynamic dependent quantities are specified as follows. At $`\theta =0`$, we apply an axis-of-symmetry boundary condition. For the outer radial boundary, we apply an outflow boundary condition. For the inner radial boundary $`r=r_{}`$ and for $`\theta =90^o`$, we apply reflecting boundary conditions.

Our simulations begin with a vertical, independent of $`r`$ magnetic field configuration and a Keplerian disk embedded in the rotating ambient medium of low density. The ambient medium is initially outflowing with the escape velocity at $`r=r_{}`$ in the direction perpendicular to the disk midplane. Below we give the details of our initial conditions and conditions in the first grid zone above the equatorial plane follows.

We proceed with setting the initial conditions in the following way. We start with adopting the rotational velocity $`v_\varphi =\sqrt{\frac{GM_{}}{r\mathrm{sin}\theta }}`$ for $`r\mathrm{sin}\theta >r_{}`$ and $`v_\varphi =0`$ elsewhere. Thus the gas above the disk rotates on cylinders with the disk Keplerian velocity whereas the gas above the non-rotating central object has zero rotational velocity.

Our initial density profile is given by the condition of hydrostatic equilibrium in the latitudinal direction for a gas with a given initial rotational velocity. To ensure an exact numerical equilibrium initially, we first initialize the density in the first grid zone above the equatorial plane, $`\rho _0`$. We assume $`\rho _0`$ is radius independent. Then we integrate the latitudinal equation of motion from $`\theta =90^{}`$ to $`\theta =0^{}`$ to compute the pressure using the numerical difference formula in our code. Finally, we compute the density from the isothermal equation of state. To reduce the problems caused by very high Alfv$`\stackrel{´}{\mathrm{e}}`$nic velocities in regions of very low density (i.e., to prevent the time step from being prohibitively small), we set a lower limit to the density on the grid as $`\rho _{min}(r)=10^{15}(r_o/r)^2\mathrm{g}\mathrm{cm}^3`$ and enforce it at all times in all models.

For the initial poloidal velocity, we adopt $`v_r=\sqrt{GM/r}\mathrm{sin}\theta `$ and $`v_\theta =\sqrt{GM/r}\mathrm{cos}\theta `$ for the region where the lower limit to the density is used (see above). This choice for the poloidal velocity is motivated by practical concerns (i.e., again to prevent the time step from being prohibitively small) and is particularly useful for models with no radiation force due to the central object. For $`x=0`$, the gas above the disk collapses onto the disk and the time step becomes so small that the simulations practically stop.

To represent steady conditions in the photosphere at the base of the wind, during the evolution of each model we continued to apply the constraint that in the first zone above the equatorial plane the density is fixed at $`\rho =\rho _0`$ at all times. During the evolution of our standard models, $`\rho _0`$ was fixed at $`10^4\mathrm{g}\mathrm{cm}^3`$.

There are several differences between our initial and boundary conditions for the hydrodynamic dependent quantities and those we used in previous work (PSD 98 and PSD 99). The most important difference is in our treatment of the wind base. In PSD 99, to represent steady conditions in the photosphere at the base of the wind, during the evolution of each model we continued to apply the constraints that in the first zone above the equatorial plane the radial velocity $`v_r=0`$, the rotational velocity $`v_\varphi `$ remains Keplerian, and the density is fixed at $`\rho =\rho _0`$ at all times. Physically, $`\rho _0`$ is analogous to the density in the photosphere of the disk at the base of the wind provided $`\rho _0`$ is relatively low. In PSD 99 and PSD 98, the interior of the disk itself was treated as negligibly thin and was excluded from the models (for a disk temperature of $`10^4`$ K at $`r=2r_{}`$, the disk scale height $`H`$ is $`H/r_{}10^3`$). In PSD 98 and PSD 99, the arbitrary value for $`\rho _0`$ was fixed typically at $`10^9\mathrm{g}\mathrm{cm}^3`$. As discussed in PSD 98, the gross properties of LD winds are unaffected by the value of $`\rho _0`$ provided it is large enough that the acceleration of the wind up to the sonic point is resolved with at least a few grid points. This technique, when applied to calculations of spherically symmetric LD winds from stars, produces a solution that relaxes to the appropriate CAK solution within a few dynamical crossing times. However, applying these constraints in MHD simulations can cause very strong evolution of the gas close to the equator and consequently the evolution of the mass outflow. We have performed many tests and found that, when the above constraints are used, the disk gets disrupted very quickly (within a couple orbits at $`r=r_{}`$) by MRI, making it impossible to study disk winds. Inclusion of magnetic fields in the model puts new constraints on $`\rho _0`$. For strong magnetic fields and low $`\rho _0`$ (i.e., a small value of the plasma parameter $`\beta 8\pi P/B^2`$ on the equatorial plane), the MRI or magnetic braking can very quickly reduce the density near the disk midplane. In particular, we observe a dramatic drop in the density, of several orders of magnitude, between the first and second zone above the equatorial plane. This change in the density profile makes impossible to study outflows from a ‘steady state’ disk. By choosing a large $`\rho _0`$ for a given magnetic field, we can reduce the dynamical importance of the magnetic field in the region close to the midplane so that the base of the wind can remain in a steady state for a long time. For these reasons, our value of $`\rho _0`$ is $`\rho _0=10^4\mathrm{g}\mathrm{cm}^3`$ and during the evolution of our standard models we allow all the dependent quantities (including $`\rho `$, $`v_r`$, and $`v_\varphi `$) to float everywhere on the grid. For the adopted value of $`\rho _0`$, the plasma parameter on the equator, $`\beta _08\pi P/B^2=8\pi c_s^2\rho _0/B^2`$ is very high for all standard models (i.e., $`2\times 10^5\beta _02\times 10^{11}`$). In Section 3.3, we discuss how our results depend on the conditions along the equatorial plane, in particular what difference it makes if $`\rho _010^9\mathrm{g}\mathrm{cm}^3`$, and $`\rho `$, $`v_r`$, and $`v_\varphi `$ are set as in PSD 99.

We are left with describing our initial magnetic configuration and the boundary conditions for the magnetic field. We assume a force-free configuration (Lorentz force $`(\times 𝐁)\times 𝐁=\mathrm{𝟎}`$) by simply setting $`(\times 𝐁)=\mathrm{𝟎}`$. We consider one straightforward initial magnetic configuration satisfying these constraints: a uniform vertical field configuration defined by the magnetic potential $`𝐀=(A_r=0,A_\theta =0,A_\varphi =Ar\mathrm{sin}\theta )`$. We scale the magnitude of the magnetic field using a parameter, $`\beta _w^{}32\pi ^2c_s^2\rho _w/B^2`$ defined for a fiducial wind density of $`\rho _w=10^{15}\mathrm{g}\mathrm{cm}^3`$:

$$A=2\pi \sqrt{(2c_s^2\rho _w/\beta _w^{})}.$$
(8)
The boundary conditions for the magnetic field are: at $`\theta =0`$, we apply an axis-of-symmetry boundary condition; for the outer radial boundary, we apply an outflow boundary condition. For the inner radial boundary $`r=r_{}`$ we apply reflecting boundary conditions while for $`\theta =90^o`$ we apply an equatorial-symmetry boundary condition (Stone & Norman 1992b). In our standard models we allow all three components of the magnetic field to float.

In reality, gas near the disk midplane (inside the disk) is turbulent because of the MRI but in near hydrostatic equilibrium. We would like to stress that, although we include the region very close to the midplane in our standard models, we do not claim that we model the disk interior. To do the latter we would need to add physical processes such as magnetic field dissipation and radiative transfer appropriate to optically thick disk gas, and to solve the equation of energy. We would also need to resolve better the disk so we could capture, for example, the fastest growing modes of the MRI. The most unstable wavelength of the MRI, $`\lambda 2\pi v_A/\mathrm{\Omega }`$, increases with height in the disk for a given angular velocity, $`\mathrm{\Omega }`$, because the Alfv$`\stackrel{´}{\mathrm{e}}`$n speed increases with height. In general, near the equator, we deal with a magnetized stratified disk where magnetic field is generated either by the MRI or magnetic braking. It is computationally prohibitive to resolve adequately a thin disk in global calculations even in two dimensions. Nevertheless, as we will discuss in section 3.3, our simulations are consistent with high resolution local simulations of magnetized stratified disks (e.g., Stone et al. 1996; Miller and Stone 2000).

To solve eqs. (1)-(3), we use the ZEUS-2D code described by Stone & Norman (1992a, 1992b).

### 2.2 Model Parameters

As in PSD 99 and Proga (2000), we calculate disk winds with model parameters suitable for a typical nMCV (see Table 1 in PSD 99, Table 1 in Proga 2000 and our Table 1). We vary the disk and central object luminosity and the strength of magnetic field (i.e., $`\beta _w^{}`$). We hold all other parameters fixed: $`M_{}=0.6\mathrm{M}_{}`$, $`r_{}=8.7\times 10^8\mathrm{cm}`$, $`c_s=14\mathrm{km}\mathrm{s}^1`$, $`k=0.2`$, $`\alpha =0.6`$ and $`M_{max}=4400`$ (see Table 1 in PSD 98). Nevertheless we can use our results to predict the wind properties for other parameters and systems – such as AGN and YSOs – by applying the dimensionless form of the hydrodynamic equations and the scaling formulae as discussed in PSD 98 and Proga (1999). We note that by adopting the uniform vertical field configuration, our model has two free parameters more than the pure LD disk wind model of PSD 99: the plasma parameter $`\beta _w^{}`$ and the initial toroidal magnetic field, $`B_\varphi `$. In this paper, we hold $`B_\varphi =0`$ but vary $`\beta _w^{}`$.

Summarizing, our model has three crucial independent parameters: the Eddington factor corresponding to the disk luminosity $`\mathrm{\Gamma }_DL_D/L_{Edd}=\frac{\sigma _e\dot{M}_a}{8\pi cr_{}}`$; the central object luminosity expressed in terms of the disk luminosity, $`x=L_D/L_{}`$, and the strength of the magnetic field expressed in terms of the plasma parameter for a fiducial density of $`\rho _w=10^{15}\mathrm{g}\mathrm{cm}^3`$, $`\beta _w^{}`$. As in PSD 99, we assume that all the radiation is emitted in the ultraviolet and does not evolve. The model also has other parameters that can be calculated self-consistently, in principle, for a given set of $`\mathrm{\Gamma }_D`$, $`x`$, $`\beta _w^{}`$ and the spectral distribution of the radiation field. These parameters are: the sound speed $`c_s`$, and the parameters of the CAK force multiplier, $`k`$, $`\alpha `$ and $`M_{max}`$. For the pure LD case, the wind solution does not depend on the value of the sound speed (see CAK for LD stellar winds and Proga 1999 for LD disk winds). The basic requirement for the strength of the disk radiation to drive a wind is:

$$\mathrm{\Gamma }_D[1+M_{max}]>1$$
(9)
(e.g., PSD 98; Proga 1999; see also Proga & Kallman 2002 and Proga 2002 for a discussion of the generalized version of this requirement). Therefore, $`k`$ and $`\alpha `$ do not matter as much as the value of $`M_{max}`$.

## 3 Results

Pure LD winds from a disk fall into two categories: 1) intrinsically unsteady with large fluctuations in density and velocity, and 2) steady with smooth density and velocity (PSD 98 and PSD 99). The type of flow is set by the geometry of the radiation field, parametrized by $`x`$: if the radiation field is dominated by the disk ($`x<1`$) then the flow is unsteady, and if the radiation is dominated by the central object ($`x>1`$) then the flow is steady. The geometry of the radiation field also determines the geometry of the flow: the wind becomes more polar as $`x`$ decreases. However, the mass-loss rate and terminal velocity are insensitive to geometry and depend more on the total system luminosity, $`L_D+L_{}`$. Regardless of the type of flow, pure LD winds consist of a dense, slow outflow that is bounded on the polar side by a high-velocity stream. The mass-loss rate is mostly due to the fast stream.

In Proga (2000), we recalculated some of the PSD 99 models to check how inclusion of the magnetocentrifugal force, corresponding to purely poloidal $`𝐁`$, will change LD disk winds. We found that flows which conserve specific angular velocity have a larger mass loss rate than their counterparts with purely LD flow, which conserve specific angular momentum. The difference in the mass loss rate between winds conserving specific angular momentum and those conserving angular velocity can be several orders of magnitude for low disk luminosities but vanishes for high disk luminosities. Winds which conserve angular velocity have much higher velocities than angular momentum conserving winds. In Proga (2000), we also found that fixing the wind geometry stabilizes winds which are unsteady when the geometry is derived self-consistently. Additionally, as expected, the inclination angle $`i`$ between the poloidal velocity and the normal to the disk midplane is important. Non-zero inclination angles allow the magnetocentrifugal force to increase the mass loss rate for low luminosities, and increase the wind velocity for all luminosities.

In this paper, we also recalculate some of the PSD 99 flows, but this time we check, by solving the full set of ideal MHD equations and allowing $`B_\varphi 0`$, how self-consistent inclusion of magnetic fields to the PSD 99 model will change LD disk winds. We summarize the properties of PSD 99 models and our new simulations in Table 1. Columns (2) to (4) give the input parameters that we varied: the mass accretion rate $`\dot{M}_a`$, the relative luminosity of the central object, $`x`$, and the plasma parameter $`\beta _w^{}`$, respectively. Column (5) lists the final time at which we stopped each simulation (all times here are in units of $`\tau =\sqrt{r_{}^3/GM_{}}=2.88\mathrm{sec}`$). Columns (6) to (8) give some the gross properties of the disk wind: the wind mass loss rate, $`\dot{M}_w`$, the wind velocity at the outer radial boundary, $`v_r(10r_{})`$, and the wind half-opening angle, $`\omega `$, respectively. We measure $`\omega `$ from the equator to the upper envelope of the wind. Table 1 also contains comments regarding some runs \[column (9)\] and explains our convention of labeling our runs.

### 3.1 Outflow from a luminous magnetized accretion disk

In this section we describe the properties and behavior of our model MHD-LD C0D in which $`\dot{M}_a=10^8\mathrm{M}_{}\mathrm{yr}^1`$, $`x=0`$, $`\beta _w^{}=2\times 10^2`$. This model is a rerun of the fiducial ‘$`x=0`$’ unsteady model discussed in detail in PSD 98 and PSD 99.

Figure 1 presents a sequence of maps showing density, velocity field and toroidal magnetic field (left, middle and right panels) from model C0D, plotted in the $`r,z`$ plane. The length of the arrows in the middle panels is proportional to $`(v_r^2+v_\theta ^2)^{1/2}`$. To show better the evolution of the wind with lower velocities we use the maximum lenght of the arrows in regions of high velocity, i.e., $`(v_r^2+v_\theta ^2)^{1/2}200\mathrm{km}\mathrm{s}^1`$. We also suppress velocity vectors in regions of low density (i.e., $`\rho `$ less than $`10^{15}\mathrm{g}\mathrm{cm}^3`$). The pattern of the direction of the arrows is an indication of the shape of the instantaneous streamlines. The solid lines in the middle panels mark the location where the poloidal Alfv$`\stackrel{´}{\mathrm{e}}`$n speed, $`v_{Ap}B_p/\sqrt{4\pi \rho }`$, equals the poloidal fluid speed, $`v_p`$ (i.e. the Alfv$`\stackrel{´}{\mathrm{e}}`$n surface). As in the pure LD case, after $`10`$ time units disk material fills the grid for $`\theta >30^o`$ and remains in that region for the rest of the run. In the early phases of the evolution the MHD-LD disk wind resembles its pure LD counterpart except that it is steadier. However, in the late phases the flow undergoes a dramatic change not seen in the LD case: slow and dense material rises from the disk at large radii. In the slow wind, the density increases whereas the velocity decreases with time until the wind settles to a time-averaged steady state. These changes in the slow wind occur on relatively a long time scale and are caused by the evolution of the magnetic field. In particular, $`B_\varphi `$ is generated near the base of the wind as indicated by the rise of the $`B_\varphi `$ contours in Figure 1 (the strength of $`B_\varphi `$ decreases with $`z`$). The gradient of $`B_\varphi ^2`$ drives the slow wind from the disk at large radii. At the end of the simulation, the LD wind is replaced by a wind driven primarily by the magnetic pressure.

Figure 2 presents the density and poloidal velocity field in the wind at the end of our simulations at 680 $`\tau `$ (top left and middle panels). The dashed and solid line in the top right panel of Figure 2 shows the contours of the angular velocity $`\mathrm{\Omega }`$ and specific angular momentum of fluid, $`Lr\mathrm{sin}\theta v_\varphi `$, respectively. The bottom left, middle and right panels of Figure 2 show the contours of the $`\beta `$ plasma parameter, the poloidal magnetic field and the contours of the toroidal magnetic field, $`B_\varphi `$, respectively. The solid lines in the bottom middle panel show the location where the strength of the poloidal magnetic field equals the toroidal magnetic field, $`B_p=B_\varphi `$.

Comparing our model C0D with its pure LD counterpart, we find that model C0D has a relatively smooth density distribution and is relatively steady. For example, the pure LD outflow is intrinsically unsteady and characterized by large amplitude velocity and density fluctuations. Infall as well as outflow from a disk can occur in different regions of the wind at the same time. However, in model C0D the poloidal velocity is mostly organized and there is no inflow onto the disk in the wind domain. Model C0D has not reached a steady state even after $`400\tau `$. However, we observe that this model is close to reaching such a state; for example, the mass flux density settles to some time-averaged maximum.

As in the pure LD case, we find that the wind consists of a dense, slow outflow that is bounded on the polar side by a high-velocity stream. However, there is an important difference, namely that the dense, slow outflow is significantly denser and somewhat faster in our model C0D than the pure LD case. Finally, as one could have expected for an MHD wind, the gas pressure dominates over the magnetic pressure in the region very close to the mid-plane ($`\beta 1`$) whereas the opposite is true in the wind domain ($`\beta 1`$, see bottom left panel).

The above differences in the two wind solutions appear to result from the Lorentz force due to the gradient of the toroidal magnetic field. Our simulations start with zero $`B_\varphi `$ but this situation changes very quickly as toroidal field is generated by rotation. We note that $`B_\varphi ^2`$ is higher than the line force in the wind domain (expect for the fast stream) by a couple of orders of magnitude. Comparing the contours of $`B_\varphi `$ (bottom right panel of Fig. 2) and the poloidal velocity field (top middle panel), we find that the velocity field is normal to the $`B_\varphi `$ contours (i.e., $`v_p`$ is parallel to the $`B_\varphi ^2`$ gradient) in the slow wind (lower right-hand corner of each panel). In contrast to the slow wind, in the fast stream where the gas is driven by the line force and the magnetic pressure is unimportant, $`v_p`$ is tangent to the $`B_\varphi `$ contours.

The quantities presented in Figures 1 and 2 show that the outflow in model C0D is not a magnetocentrifugal wind. There are several diagnostics of magnetocentrifugal driving. For example, in the acceleration zone of a steady state magnetocentrifugal wind $`B_p>B_\varphi `$. However in model C0D, the location where $`B_p=B_\varphi `$ is near the disk and $`B_p<B_\varphi `$ in most of the acceleration zone (see lower middle panel in Figure 2). Additionally, in a magnetocentrifugal wind the total angular momentum per unit mass, $`lv_\varphi r\mathrm{sin}\theta r\mathrm{sin}\theta B_\varphi B_p/(4\pi \rho v_p)`$, the second term due to the twisted magnetic fields is comparable to the first term. Subsequently, the wind is corotating with the underlying disk up to approximately the Alfv$`\stackrel{´}{\mathrm{e}}`$n point. In model C0D, the Alfv$`\stackrel{´}{\mathrm{e}}`$nic surface is very close to the disk (see top middle panel) and the total specific angular momentum is mostly due to the fluid even near the wind base. Thus the wind is corotating with the disk only over a relatively small length scale. The top right panel of Figure 2 is helpful to distinguish between a magnetocentrifugal wind and a fluid angular momentum conserving wind. In the former, the angular velocity is conserved along the streamlines below the Alfv$`\stackrel{´}{\mathrm{e}}`$n point while in the latter, $`L`$ is conserved along the streamlines. Comparing the top right panel and the top middle panel of Figure 2, we clearly see that the contours of $`L`$, not of $`\mathrm{\Omega }`$, are aligned with the streamlines represented by the arrows of the poloidal velocity field.

Next we consider the angular dependence of the flow at large radii. Figure 3 shows the angular dependence of density, radial velocity, mass flux density, and accumulated mass loss rate at $`r=r_o=10r_{}`$ at the end of the simulation of model C0D. The accumulated mass loss rate is given by:

$$d\dot{m}(\theta )=4\pi r_o^2_{0^o}^\theta \rho v_r\mathrm{sin}\theta d\theta .$$
(10)
The gas density is a very strong function of angle for $`\theta `$ between $`90^o`$ and $`25^o`$. Between the disk mid-plane at $`\theta =90^o`$ and $`\theta 85^o`$, $`\rho `$ drops by $`7`$ orders of magnitude, as expected for a density profile determined by hydrostatic equilibrium. For $`25^o<\theta <85^o`$, the wind domain, $`\rho `$ varies between $`10^{16}`$ and $`10^{11}\mathrm{g}\mathrm{cm}^3`$. For $`\theta <25^o`$, density again decreases rapidly, but this time to so low a value that it becomes necessary to replace it by the numerical lower limit $`\rho _{min}`$. The region with $`\rho \rho _{min}`$ is not relevant to our analysis as it has no effect on the disk flow. The radial velocity at $`10r_{}`$ increases gradually from zero at the equator to $`100\mathrm{km}\mathrm{s}^1`$ at $`\theta 60^o`$, then it drops to nearly zero at $`\theta 60^o`$. Over the angular range $`65^o>\theta >25^o`$, $`v_r`$ increases from $`0`$ up to $`1200\mathrm{km}\mathrm{s}^1`$.

The cumulative mass loss rate is negligible for $`\theta <35^o`$ because of the very low prevailing gas density. Beginning at $`\theta >35^o`$, $`d\dot{m}`$ increases to $`3\times 10^{13}\mathrm{g}\mathrm{s}^1`$ at $`\theta 55^o`$. This increase of $`d\dot{m}`$ is due to the fast stream. Then, in the slow dense outflow, the cumulative mass loss rate increases to $`7\times 10^{14}\mathrm{g}\mathrm{s}^1`$ at $`\theta 86^{}`$. For even higher $`\theta `$, in the region close to the disk plane, where the gas density starts to rise very sharply and where the motion is subsonic and typically more complex, the cumulative mass loss rate is subject to enormous fluctuations (some of which may even be negative). In the example shown in Figure 3, the total mass loss rate through the outer boundary, $`\dot{M}_{tot}=`$ $`d\dot{m}(90^o)`$ reaches $`10^{19}\mathrm{g}\mathrm{s}^1`$! This value of $`\dot{M}_{tot}=`$ is most certainly dominated by the contribution from the slow-moving region very close to the disk mid-plane – a contribution that is very markedly time-dependent. We ignore the value of $`\dot{M}_{tot}`$ as it is related more to subsonic oscillations of the ’disk’ rather than to the supersonic wind we model. In the remaining part of the paper, for the wind mass loss rate we use the value of the cumulative mass loss rate at $`\theta =82^{}`$.

We note that the increase of the mass loss rate in the slow wind is totally due to magnetic fields and there is no enhancement of the line force in this region. In fact, the opposite is true – the line force in the dense slow region in model C0D is significantly reduced by the action of the magnetic fields. The latter increases the density of the outflow which in turn reduces the line force. The slight increase in the velocity and its gradient in the slow wind is far from compensating the reduction of the line force due to the increase of the density.

We conclude that magnetic fields can change qualitatively and quantitatively a radiation-driven disk wind. In particular, the magnetic pressure can dominate the driving of the wind and reduce the role of the line force. In model C0D, we find that a disk loses mass via a LD wind in the inner part of the wind (the fast stream) and via a MHD-driven wind in the outer part of the wind (the dense slow wind). We expect that increasing the relative strength of the magnetic pressure to the radiation pressure (e.g., by reducing $`\beta _w^{}`$ and therefore $`\beta _0`$) should lead to the entire wind domain being driven by the magnetic force. Conversely, decreasing the magnetic pressure should lead to the wind being driven by the line force. Therefore we consider next a limited parameter survey to check whether our expectation is correct and to see how the wind solution changes quantitatively with the strength of the radiation and magnetic fields.

### 3.2 Parameter survey

We consider only the parameter space of our models that will define the major trends in disk wind behavior. We focus on a survey of how the mass loss rate, outflow velocity and geometry change with disk luminosity, relative central object luminosity and strength of magnetic field. In Section 3.3 we will consider how our results depend on our treatment of the conditions in the first grid zone above the equatorial plane.

In Figure 4 we show (a) the wind mass loss rate, $`\dot{M}_w`$, as a function of the total Eddington factor, $`(1+x)\mathrm{\Gamma }_D`$ (note that $`(1+x)\mathrm{\Gamma }_D`$ is proportional $`\dot{M}_a`$ for a given $`x`$) for various $`\beta _w^{}`$ and (b) $`\dot{M}_w`$, as a function of $`\beta _w^{}`$, for various $`(1+x)\mathrm{\Gamma }_D`$. The top panel of Figure 4 also shows $`\dot{M}_w`$ as a function of $`(1+x)\mathrm{\Gamma }_D`$ as predicted by the pure LD wind model with and without taking into account the fact that the force multiplier has a maximum value. The thick solid line corresponds to the prediction where $`M(t)`$ can be arbitrarily high whereas the thin solid line corresponds to the prediction where $`M(t)`$ reaches maximum at $`M_{max}=4400`$ (see Proga 1999). In the top panel of Figure 4, it can be seen that $`\dot{M}_w`$ is a very strong function of $`(1+x)\mathrm{\Gamma }_D`$ for $`\beta _w^{}=\mathrm{}`$. As shown in PSD 98 and Proga (1999), two dimensional models of LD disk winds predict mass loss rates (as well as velocities) very similar to those predicted by the original CAK formulae when the stellar Eddington factor is replaced with the the total Eddington factor and the formulae are corrected for the fact that for small Eddington factors the force multiplier reaches its maximum and consequently the radiation force cannot exceed gravity. Note that the drop in $`\dot{M}_w`$ occurs for $`(1+x)\mathrm{\Gamma }_D0.0002=1/M_{max}`$.

Motivated by the conclusion from section 3.1, we have performed a few simulations using different values of $`\beta _w^{}`$. We would like to estimate a $`\beta _w^{}`$ range for which the line force dominates and magnetic fields play a small role, and a range for which the line force is unimportant and the wind is totally controlled by magnetic fields.

We here focus on simulations for a fiducial ‘$`x=0`$’ model with $`\dot{M}_a=10^8\mathrm{M}_{}\mathrm{yr}^1`$ and a fiducial ‘$`x=1`$’ steady state model with $`\dot{M}_a=\pi \times 10^8\mathrm{M}_{}\mathrm{yr}^1`$ and various $`\beta _w^{}`$. For both series of simulations (i.e., series ‘C0’ corresponding to runs labeled C0A, C0B, C0D, C0E, and C0F and series ‘D1’ corresponding to runs labeled D1A, D1B, D1D, D1E, D1F, and D1G, see Table 1), we find that indeed for sufficiently high $`\beta _w^{}`$ the wind mass loss rate and characteristic velocity are very similar for those in the pure LD counterparts. For example, we note that for models of series D1, the symbols corresponding to $`\dot{M}_w`$ with $`\beta _w^{}=2,2\times 10^2`$, and $`2\times 10^3`$ overlap with the the symbol corresponding to the pure LD model D1A (see also the bottom panel of Figure 4). However, there is a qualitative difference between the MHD-LD winds and the LD wind even for high $`\beta _w^{}`$ cases: the MHD-LD winds are less unsteady than the LD winds. For high $`\beta _w^{}`$, the Lorentz force can be too weak to drive a wind but if $`\beta 1`$ in the wind (as in the cases we explored), the gradient of the magnetic pressure can reduce the density fluctuation near the wind base and in the wind.

The bottom panel of Figure 4 shows more clearly the point we made above that for high $`\beta _w^{}`$, $`\dot{M}_w`$ is as for pure LD cases and does not change with $`\beta _w^{}`$. However, starting from a certain value of $`\beta _w^{}`$, $`\dot{M}_w`$ becomes a strong function of $`\beta _w^{}`$. The value of $`\beta _w^{}`$ at which $`\dot{M}_w`$ starts to increase with decreasing $`\beta _w^{}`$ is higher for the lower luminosity system (the solid line) than for the higher luminosity system (the dashed line). Generally, the strong $`\dot{M}_w`$ dependence on $`(1+x)\mathrm{\Gamma }_D`$ almost disappears when strong magnetic fields are added to the model and this is clear in the top panel, where $`\dot{M}_w`$ becomes insensitive to $`\mathrm{\Gamma }_D`$ as $`\beta _w^{}`$ decreases.

As Table 1 shows, the wind geometry is also sensitive to $`\beta _w^{}`$: the half-opening angle of the wind $`\omega `$, increases with decreasing $`\beta _w^{}`$ so the wind becomes more polar as $`\beta _w^{}`$ decreases. We note that even the ‘$`x=1`$’ models become polar for $`\beta _w^{},2\times 10^4`$ (see models D1F and D1G).

Simulations for fixed $`\beta _w^{}`$ and $`x`$, but with varying the system luminosity, help us to check whether the magnetic fields can drive a wind for disk luminosities too low to the drive a LD wind. These simulations also help us to check whether the line force can ‘regain’ control over the wind if the system luminosity increases.

For the system luminosity of model B1A, the line force is too weak to produce a supersonic outflow (e.g., PSD 99). However, in model B1F with the same radiation field as in model B1A, there is a strong robust disk outflow. In the models B1F and D1F with the same low $`\beta _w^{}`$ but different system luminosities, the wind is mostly driven by MHD. However in model E1F, with the system luminosity higher than in model D1F by a factor of 10 the wind is LD and very similar to that with zero magnetic field (compare model E1F with model E1A).

We conclude that MHD driving is robust and does not require the line force, e.g., the MHD driving does not require the line force to launch a wind from the disk. We find a negative feedback between the magnetic field and the line force, i.e., the higher the magnetic field the lower the line force and vice versa. We also conclude that although MHD driving can produce a strong wind it does so by driving a relatively dense wind with relatively low velocities. The latter do not seem to depend on the escape velocity from the radius from which the wind is launched. We note that $`v_r`$ of the slow wind depends very weakly on $`\theta `$ at $`r_o`$. On the other hand, the line force drives a relatively fast wind with a velocity sensitive to the launching radius and consequently to $`\theta `$ (PSD 98).

### 3.3 Dependence of wind evolution and properties on the treatment of the wind base

Our model for mass outflows from accretion disks is a hybrid of an LD model and an MHD-driven flow. We have extensively tested our LD disk wind model and applied it to several systems (see references in Section 1). In this section we present a brief review of our test runs designed to check the MHD part of our model. In particular, we have performed a few tests aimed at reproducing qualitatively results already published on MHD disk winds. Additionally, we have explored a parameter region of our model ($`0.1<\beta _0<10`$) for which we expect to resolve the fastest growing modes of the MRI.

There have been many numerical studies of MHD disk winds (our list of references in the introduction is far from complete). For simplicity of our presentation we will reference in more detail to the work of only two groups: Stone & Norman (1994) and Ouyed & Pudritz (1997a; 1997b; 1999). The former included the disk as well as the wind in their simulations whereas the latter included only the wind and treated the disk as the lower boundary of the wind. Both approaches have been commonly adopted in the literature.

Clearly, inclusion of the disk structure in calculations of disk outflow is highly desirable. However, the physics of the disk is very complex and its proper modeling is very demanding. On the other hand, one would hope that it is possible to capture the key elements of a disk wind without modeling the disk interior. Work by Ouyed & Pudritz (1997a; 1997b; 1999) is an example of studies of magnetocentrifugal disk winds while PSD 98, PSD 99 is an example of studies of LD disk winds where the disk interior was not included. However there is an important difference between modeling a magnetocentrifugal disk wind and an LD disk wind, namely the treatment of the lower boundary for the disk. In numerical simulations of magnetocentrifugal disk winds the mass flux density from the disk is given whereas in numerical simulations of LD disk winds the mass flux density is a result. The key reason for this difference is the location of the critical surface at which the mass flux density is determined. In MHD simulations, this surface corresponds to the slow magneto-sonic surface which is located inside the disk, below the photosphere. Thus unless the disk is included in simulations, one must assume the mass flux density from the lower boundary. On the other hand, the critical surface in LD disk winds is in the supersonic part of the wind above the photosphere. Therefore, including the photosphere of the disk but not the whole disk suffices to capture to the transition between the sub-critical and super-critical parts of the LD wind and subsequently the mass flux density is determined by the physics of the flow (i.e., Feldmeier & Shlosman 1999).

We have performed test simulations of pure MHD and MHD-LD disk winds following the approach used in studies of magnetocentrifugal winds. Our test runs included those where we set $`v_r=0`$, $`v_\varphi =v_{Keplerian}`$, $`\rho =\rho _0`$ and $`B_\varphi =0`$ in the first grid zone above the equatorial plane at all times. This setting is very similar to the one used in numerical simulations of magnetocentrifugal disk winds (Ouyed & Pudritz 1997a; 1997b; 1999) with the exception that they also set the velocity in the direction perpendicular to the midplane, $`v_z`$, to some small subsonic velocity (i.e., the mass flux density $`\rho v_z`$ is fixed). In the test runs, we allow $`v_\theta `$ to float and let the mass flux density be fixed by the solution to the problem.

In short, we found such an approach to modeling MHD and MHD-LD disk winds unsatisfactory. For example, in pure MHD cases insisting on a Keplerian flow in the first grid zone above the equatorial plane at all times does not allow proper modeling of magnetic braking of the disk. By setting $`v_r=0`$, we prevent or at least significantly reduce the collapse of the disk and dragging of magnetic field line (e.g., Ouyed & Pudritz 1997a). However, when we allowed all variables to float we could successfully reproduce the evolution of a collapsing disk and an MHD outflow. In particular, in strong field cases, the disk undergoes dramatic evolution on a short time scale (i.e., of order of $`40\tau `$; in other words, a couple of orbital periods at $`r_{}`$), consistent with the time scale for magnetic braking of an aligned rotator (Mouschovias & Paleologou 1980).

Similarly, we found that the evolution of the disk in a pure MHD case depends on the conditions in the first grid zone above the equatorial plane for weak magnetic fields. We performed several test simulations with weak magnetic field, so that the disk was unstable to MRI. In particular, we adopted $`\rho _0=10^9\mathrm{g}\mathrm{cm}^3`$ and $`\beta _w^{}=8\times 10^7`$ yielding $`\beta _0=6.4`$, the parameters for which the fastest growing MRI mode is resolved for our stratified disk. For quantities in the first grid zone above the equatorial plane set as in magnetocentrifugal models, we found that the disk evolved very quickly. Initially, we observed a characteristic exponential growth of the magnetic field, which causes the disk to separate vertically into horizontal planes. This so-called channel solution is typical of the development of the axisymmetric phase of the MRI (Hawley & Balbus 1991; Goodman & Xu 1994). As expected, the gas streaming is directed outward in the channel closest to the equator and inwards in the channel higher up from the equator. However, after about three orbital periods the disk is destroyed. We find many similarities in the behavior of the disk in our simulations and in previous simulations. For example, Stone et al. (1996) and Miller & Stone (2000) reported that in their three-dimensional simulations of a stratified disk (a local shearing box approximation), an initially uniform vertical magnetic field leads to a strong radial streaming in the early phases of the evolution. However, in the late phases, high magnetic pressures disrupt the disk and the entire computational domain is magnetically dominated, i.e., $`\beta <1`$ even near the equatorial plane. The high magnetic pressures are created by strong large-scale radial streaming that occurs at large heights because in a stratified disk, most unstable wavelength increases with height. As in those local three dimensional simulations, we observed that strong radial streaming occurs in the horizontal planes into which our disk separates. Miller & Stone (2000) concluded that enforcing Keplerian rotation in the boundary conditions is appropriate for following the long-term evolution. We have arrived at the same conclusion, that it is simply impossible to model on a sufficiently long-time scale an outflow from a disk that is quickly disrupted. We found that when we do not enforce Keplerian rotation, more precisely when we allow all the fluid and magnetic field quantities to float, then the disk evolves more slowly with radial streaming. However when we add the line force and assume too low a density the disk can be totally disrupted because of the LD wind.

Motivated by the results from the above test simulations, for our standard models, we decided to allow all quantities to float in the first grid zone above the equatorial plane so that they can evolve self-consistently in all models in Tabel 1. To ensure we retain a disk for long enough to acquire a reasonably settled outflow, we increased $`\rho _0`$ from $`10^9\mathrm{g}\mathrm{cm}^3`$ as in PSD 99 to $`10^4\mathrm{g}\mathrm{cm}^3`$. The latter change means that we do not resolve the fastest growing MRI mode in the disk. To do so we would have to increase the strength of the magnetic field for a given resolution or increase the resolution. Unfortunately, we cannot afford either of these modifications because they would make our already computationally demanding simulation even more demanding. Therefore, in this first attempt to model MHD-LD disk winds we resolve the wind from the disk but do not resolve the disk itself. Thus we cannot resolve the fastest growing MRI mode within the disk scale height, but only most unstable modes at larger heights.

We finish with a remark that in the literature most of the numerical simulations of MHD disk winds were stopped after one or two or at most a few orbital periods. The reason for this is obvious but it is disappointing that one cannot perform yet global simulations of magnetized disks, treating properly accretion as well as outflow and choosing physical parameters comparable to those in real systems.

### 3.4 Limitations of models and future work

The most important limitation of our model is an inadequate spatial resolution for modeling the MRI inside the disk (see discussion above). As is fitting for a first exploration of MHD-LD wind models from disks, we aim to examine only the parameter space of our models that will define the major trends in disk wind behavior. Therefore our priority is to set up the simulation in such a way that the base of the wind is relatively stable and corresponds to a steady state accretion disk. Obviously the problem with modeling MRI disappears when the magnetic field is strong. In the context of MHD winds, a strong magnetic field case for which the disk is MRI-stable is also physically interesting.

We have explored MHD-LD models where initially $`\beta <1`$ everywhere on the grid. However, we found then Alfv$`\stackrel{´}{\mathrm{e}}`$n speeds so high that the resulting time step was extremely small. We emphasize that our choice of $`\beta _w^{}`$ is limited by the constraints on the density in the computational domain, i.e., the need to have a relatively large contrast between the density near the equatorial plane and that high above the plane. Increasing arbitrarily the lower limit, $`\rho _{min}`$, to reduce the density contrast and subsequently reduce the Alfv$`\stackrel{´}{\mathrm{e}}`$n speed is not suitable to modeling LD winds. The line force is very sensitive to the density and our wind solution can then include gas with a spurious density set by the lower limit. Thus we were very cautious in choosing $`\rho _{min}`$, so that the region with $`\rho \rho _{min}`$ has no effect on the disk flow.

We note that other numerical simulations explored a relatively low density contrast between the disk and the ambient gas (e.g., in Stone & Norman 1994; Ouyed & Pudritz 1997a, 1997b, 1999; and Krasnopolsky, Li & Blandford 1999 the density contrast is $`<10^3`$). Here we deal with density contrasts several orders of magnitude higher.

The fact that our results strongly depend on the magnetic field points to a need to explore different configurations for the initial magnetic field and to move from two-dimensional axisymmetric simulations to fully three-dimensional simulations. We are interested in the long-time evolution of the flow. Therefore three-dimensional simulations are required as there exist no self-sustained axisymmetric dynamos. Thus, contrary to the outflows from stars, simulations of outflows from magnetized disks – with or without radiation pressure – should include the disks themselves, not just the disk photosphere, and should be performed in three dimensions.

## 4 Conclusions

We have studied winds from accretion disks with magnetic fields and the radiation force due to lines. We use numerical methods to solve the two-dimensional, time-dependent equations of ideal MHD. We have accounted for the radiation force using a generalized multidimensional formulation of the Sobolev approximation. For the initial conditions, we have considered uniform vertical magnetic fields and geometrically thin, optically thick Keplerian disks. We allow the magnetic field and the disk to evolve. In particular, we do not enforce Keplerian rotation in the first grid zone above the equatorial plane. Although the gas near the equatorial plane departs only slightly from Keplerian rotation in our self-consistent calculations, its evolution is notably different when we enforce Keplerian rotation. We find that the magnetic fields very quickly start deviating from purely vertical due to the MRI. This leads to fast growth of the toroidal magnetic field as field lines wind up due to the disk rotation. As a result, the toroidal field dominates over the poloidal field above the disk and the gradient of the former drives a slow and dense disk outflow, which conserves specific angular momentum. Depending on the strength of the magnetic field relative to the system luminosity, the disk wind can be radiation- or MHD driven. For example, for our model parameters $`\dot{M}_a=10^8\mathrm{M}_{}\mathrm{yr}^1`$ and $`x=0`$, the wind is radiation-driven for $`\beta _w^{}>1`$ and, as its pure LD counterpart, consists of a dense, slow outflow that is bounded on the polar side by a high-velocity stream. The mass-loss rate is mostly due to the fast stream. As the magnetic field strength increases (i.e., $`\beta _w^{}<10^1`$), first the slow part of the flow is affected. In particular, the slow wind becomes denser and faster than its pure LD counterpart by a factor of $`>100`$ and $`>3`$, respectively. Consequently, the dense wind begins to dominate the mass-loss rate. In very strong magnetic field (i.e., $`\beta _w^{}<10^3`$) or pure MHD cases, the wind consists of only a dense, slow outflow without the presence of the distinctive fast stream so typical of pure LD winds. Our simulations indicate that winds launched by magnetic fields are likely to remain dominated by the fields downstream because of their relatively high densities. The radiation force due to lines may not be able to change a dense MHD wind because the line force strongly decreases with the density.

Our results show that, as expected, a hybrid model predicts mass-loss rates higher than those predicted by a pure LD model. We plan to compute synthetic line profile based on our simulations and check whether our MHD-LD models can resolve the problem of nMCV winds. However even now we can say that our MHD-LD models may not solve the problem because to explain nMCV winds we need a model that predicts not only a higher $`\dot{M}_w`$ than a LD wind model but also $`\dot{M}_w`$ must be mostly due to a fast wind not a dense slow wind as we find our models (Proga et al. 2002).

ACKNOWLEDGMENTS: We thank J.M. Stone, M.C. Begelman, J.E. Drew, P.J. Armitage, A. Feldmeier, and M. Ruszkowski for useful discussions. In particular, we thank A. Feldmeier for sharing with us his results and conclusions on MHD-LD disk winds. We also thank an anonymous referee for comments that helped us clarify our presentation. We acknowledge support from NASA under LTSA grant NAG5-11736. We also acknowledge support provided by NASA through grant AR-09532 from the Space Telescope Science Institute, which is operated by the Association of Universities for Research in Astronomy, Inc., under NASA contract NAS5-26555. Computations were partially supported by NSF grant AST-9876887.

## REFERENCES

Abbott, D.C. 1982, ApJ, 259, 282

Balbus, S.A., & Hawley, J.F. 1991, ApJ, 376, 214

Balbus, S.A., & Hawley, J.F. 1998, Rev. Mod. Phys., 70, 1

Batchelor G.K. 1967, An Introduction to Fluid Mechanics (Cambridge: Cambridge University Press)

Blandford R.D., Payne D.G. 1982, MNRAS, 199, 883

Cannizzo J. K., Pudritz R.E. 1988, ApJ, 327, 840

Castor J.I., Abbott D.C., Klein R.I. 1975, ApJ, 195, 157 (CAK)

Contopoulos J. 1995, ApJ, 450, 616

Chandrasekhar, S. 1960, Proc. Nat. Acad. Sci., 46, 253

Drew J.E., Proga D., in “Cataclysmic Variables”, Symposium in Honour of Brian Warner, Oxford 1999, ed. by P. Charles, A. King, D. O’Donoghue, in press

Feldmeier, A., & Shlosman, I. 1999, ApJ, 526, 344

Feldmeier, A., Shlosman, I., &, Vitello, P. 1999, ApJ, 526, 357

Gayley, K.G. 1995, ApJ, 454, 410

Goodman J., Xu G. 1994 ApJ, 432, 213

Hawley J.F, Balbus S.A. 1991, ApJ, 376, 214

Kato, S.X., Kudoh T., & Shibata K. 2002, ApJ, 565, 1035

K$`\ddot{\mathrm{o}}`$nigl A. 1993, in “Astrophysical Jets”, ed. by D.P. O’Dea (Cambridge: Cambridge Univ. Press), 239

Krasnopolsky R., Li Z.-Y., Blandford R., 1999, ApJ, 526, 631

Kudoh, T., Matsumoto, R., & Shibata, K. 1998, ApJ, 508, 186

Kudoh T., Shibata K. 1997, ApJ, 474, 362

Li, Z.-Y. 1995, ApJ, 444, 848

Li, Z.-Y. 1996, ApJ, 465, 855

Miller, K.A., Stone, J.M. 2000, ApJ, 534, 398

Mouschovias, T. Ch., Paleologou, E.V. 1980, ApJ, 237, 877

Ogilvie G.I., Livio M. 1998, ApJ, 499, 329

Ouyed R., Pudritz R.E. 1997a, ApJ, 482, 717

Ouyed R., Pudritz R.E. 1997b, ApJ, 484, 794

Ouyed R., Pudritz R.E. 1999, MNRAS, 309, 233

Owocki S.P., Cranmer S.R., Gayley K.G. 1996, ApJ, 472, L115

Pelletier G., Pudritz R.E. 1992, ApJ, 394, 117

Pereyra N.A., Kallman T.R. Blondin J.M. 1997, ApJ, 477, 368

Pereyra N.A., Kallman T.R. Blondin J.M. 2000, ApJ, 532, 563

Proga D. 1999, MNRAS, 304, 938

Proga D. 2000, ApJ, 538, 684

Proga D. 2002, Mass Outflow in Active Galactic Nuclei: New Perspectives, ASP Conference Proceedings, Vol. 255. Edited by D. M. Crenshaw, S. B. Kraemer, and I. M. George. ISBN: 1-58381-095-1. San Francisco: Astronomical Society of the Pacific, 2002., p.309

Proga, D., & Kallman, T.R. 2002, ApJ, 565, 455

Proga, D., Kallman, T.R., Drew, J.E., & Hartley, L.E., 2002, ApJ, 572, 382

Proga D., Stone J.M., Drew J.E. 1998, MNRAS, 295, 595 (PSD 98)

Proga D., Stone J.M., Drew J.E. 1999, MNRAS, 310, 476 (PSD 99)

Proga D., Stone J.M., Kallman T.R. 2000, ApJ, 543, 686

Pudritz R.E., Norman C.A. 1986, ApJ, 301, 571

Romanova, M.M., Ustyugova, G.V., Koldoba, A.V., Chechetkin, V.M., & Lovelace, R.V.E. 1997, ApJ, 482, 708

Shakura N.I., Sunyaev R.A. 1973 A&A, 24, 337

Shibata K., Uchida Y. 1986, PASJ, 38, 631

Stone, J.M., Hawley, J.F., Gammie, C.F., & Balbus, Steven A. 1996 ApJ, 463, 656

Stone J.M., Norman M.L. 1992a, ApJS, 80, 753

Stone J.M., Norman M.L. 1992b, ApJS, 80, 791

Stone J.M., Norman M.L. 1994, ApJ, 433, 746

Uchida Y., Shibata K. 1985, PASJ, 37, 515

Ustyugova G.V., Koldoba A.V., Romanova M.M. Chechetkin V.M. Lovelace R.V.E. 1995, ApJ, 439, 39L

Ustyugova G.V., Koldoba A.V., Romanova M.M. Chechetkin V.M. Lovelace R.V.E. 1999, ApJ, 516, 221

Velikov, E.P. 1959, JETP, 36, 1398

Vitello P.A.J., Shlosman I., 1988 ApJ, 327, 680

Fig. 1. A sequence of density maps (left), velocity fields (middle) and contours of the toroidal magnetic field (right) from run C0D after 37.8, 166.2, 294.7, 423.1, and 551.6 $`\tau `$ (time increases from top to bottom). Run C0D is the example of a MHD-LD flow discussed in detail in section 3.1. For clarity, we suppress velocity vectors for regions with very low density (i.e., $`\rho <10^{15}\mathrm{g}\mathrm{cm}^3`$). The solid line overplotted on the velocity maps marks the the Alfv$`\stackrel{´}{\mathrm{e}}`$nic surface (i.e., location where $`v_{Ap}=v_p`$, middle panels). The $`B_\varphi `$ contours are for -30,-20,-10,-5.,-1, 0, and 5. Dotted lines denote negative values of $`B_\varphi `$.

Fig. 2. Two-dimensional structure of several quantities from model C0D after 680 $`\tau `$. The contours of $`L`$ are for 1.0, 1.5, 2 ,2.5 , and 3.0 while for $`\mathrm{\Omega }`$ are for 0.025, 0.05, 0.1, 0.2, 0.4, and 0.8. Both $`L`$ and $`\mathrm{\Omega }`$ are in units of the specific angular momentum and angular velocity on the Keplerian disk at $`r_{}`$. The contours of $`log\beta `$ are for -3, -2.5 , -2.0, and -1.0. The solid lines overplotted on the poloidal magnetic field (bottom middle panel) show the location where the strength of the poloidal magnetic field equals the toroidal magnetic field, $`B_p=B_\varphi `$. The contours for $`B_\varphi `$ are as in Figure 1.

Fig. 3. Quantities at the outer boundary in model C0D after 680 $`\tau `$. The ordinate on the left hand side of each panel refers to the solid line, while the ordinate on the right hand side refers to the dotted line. Of particular note is the continuous strong increase of the accumulated mass loss rate with increasing polar angle. This is associated with the fact that the mass loss rate is dominated by a slow dense wind originating at large radii.

Fig. 4. Model mass loss rates as functions of the total Eddington factor (top panel) and the initial plasma parameter, $`\beta _w^{}`$ (bottom panel). In the top panel, open circles are for pure LD models, and the other shapes correspond to different values of $`\beta _w^{}`$ ($`\beta _w^{}=2\times 10^0`$ crosses, $`\beta _w^{}=2\times 10^2`$ asterisks, $`\beta _w^{}=2\times 10^3`$ diamonds, $`\beta _w^{}=2\times 10^4`$ triangles, and $`\beta _w^{}=2\times 10^6`$ squares). Table 1 specifies other model parameters. The thick solid line represents results in the stellar CAK case while the thin solid line represents the stellar CAK case corrected for a finite value of the maximum force multiplier $`M_{max}`$ (see the text for more detail). The alternative ordinate on the right hand side of the lower panel is the dimensionless wind mass loss rate parameter $`\dot{M}_w^{}`$ defined in PSD 98 (see equation 22 in PSD 98). The solid line in bottom panel represents the mass loss rate for models of series C0 while the dashed line represents the mass loss rate for models of series D1 (see section 3.2).

Table 1. Summary of results for disc winds

| Run | $`\dot{M}_a`$ | $`x`$ | $`\beta _w^{}`$ | $`t_f`$ | $`\dot{M}_w`$ | $`v_r(10r_{})`$ | $`\omega `$ | comments |
| --- | --- | --- | --- | --- | --- | --- | --- | --- |
|  | (M yr<sup>-1</sup> ) |  |  | $`\tau `$ | (M yr<sup>-1</sup>) | $`(\mathrm{km}\mathrm{s}^1)`$ | degrees |  |
| pure LD |  |  |  |  |  |  |  |  |
| C0A | $`10^8`$ | 0 | $`\mathrm{}`$ | 800 | $`5.5\times 10^{14}`$ | 900 | 50 | run A in PSD 99 |
| B1A | $`\pi \times 10^9`$ | 1 | $`\mathrm{}`$ | 1300 |  |  |  | no supersonic outflow |
| D1A | $`\pi \times 10^8`$ | 1 | $`\mathrm{}`$ | 1000 | $`2.1\times 10^{11}`$ | 3500 | 32 | run C in PSD 99 |
| E1A | $`\pi \times 10^7`$ | 1 | $`\mathrm{}`$ | 1000 | $`1.6\times 10^9`$ | 7000 | 32 |  |
| pure MHD |  |  |  |  |  |  |  |  |
| A0C | $`0`$ | 0 | $`2\times 10^1`$ | 1500 | $`6.4\times 10^{13}`$ | 100 | 65 |  |
| LD-MHD |  |  |  |  |  |  |  |  |
| C0B | $`10^8`$ | 0 | $`2\times 10^0`$ | 400 | $`6.4\times 10^{14}`$ | 900 | 55 |  |
| C0D | $`10^8`$ | 0 | $`2\times 10^2`$ | 680 | $`>7.1\times 10^{12}`$ | 1000 | 60 |  |
| C0E | $`10^8`$ | 0 | $`2\times 10^3`$ | 450 | $`>3.2\times 10^{11}`$ | 1000 | 65 |  |
| C0F | $`10^8`$ | 0 | $`2\times 10^4`$ | 170 | $`>1.0\times 10^{10}`$ | 1300 | 65 |  |
| B1F | $`\pi \times 10^9`$ | 1 | $`2\times 10^4`$ | 300 | $`1.3\times 10^{10}`$ | 600 | 70 |  |
| D1B | $`\pi \times 10^8`$ | 1 | $`2\times 10^0`$ | 250 | $`2.4\times 10^{11}`$ | 3500 | 32 |  |
| D1D | $`\pi \times 10^8`$ | 1 | $`2\times 10^2`$ | 460 | $`2.7\times 10^{11}`$ | 3500 | 35 |  |
| D1E | $`\pi \times 10^8`$ | 1 | $`2\times 10^3`$ | 670 | $`>1.0\times 10^{10}`$ | 3500 | 40 |  |
| D1F | $`\pi \times 10^8`$ | 1 | $`2\times 10^4`$ | 420 | $`>4.8\times 10^{10}`$ | 3500 | 60 |  |
| D1G | $`\pi \times 10^8`$ | 1 | $`2\times 10^6`$ | 120 | $`>1.9\times 10^{09}`$ | 3500 | 75 |  |
| E1F | $`\pi \times 10^7`$ | 1 | $`2\times 10^4`$ | 260 | $`1.6\times 10^9`$ | 7000 | 32 |  |

$``$ We use the following convention to label our runs: the first character in the name refers to $`\dot{M}_a`$, (i.e., A, B, C, D, and E are for 0, $`\pi \times 10^9`$, $`10^8`$, $`\pi \times 10^8`$, $`\pi \times 10^7`$, respectively). The second character refers to $`x`$ (i.e., 0 and 1, are for 0 and 1, respectively) and finally the third character refers to $`\beta _w^{}`$ (i.e., A, B, C, D, E, F, and G are for $`\mathrm{}`$, $`2\times 10^0`$, $`2\times 10^1`$, $`2\times 10^2`$, $`2\times 10^3`$, $`2\times 10^4`$, and $`2\times 10^6`$, respectively).


text:
# A reinterpretation of Volcano Ranch lateral distribution measurements to infer the mass composition of cosmic rays

## 1 Introduction

Over its lifetime, the Volcano Ranch array  collected data on the lateral distribution of air showers produced by cosmic rays of energies above $`10^{17}`$ eV. This lateral distribution is parameterized by a so-called lateral distribution function. The steepness of this function, given by $`\eta `$ , is sensitive to the depth of maximum ($`X_{\mathrm{max}}`$) of the shower, and therefore to the primary composition and to the character of the initial hadronic interactions.

Our analysis of precise Volcano Ranch measurements of $`\eta `$ is the first analysis of Volcano Ranch data using modern Monte Carlo tools. To simulate the development of the air shower we use the aires  code (version 2.4.0), which contains the hadronic interaction generators qgsjet98  and sibyll2.1 . The results of these simulations are convolved with a simulation of the detector response carried out using geant4 .

## 2 Volcano Ranch

The pioneering Volcano Ranch array was an array of scintillation counters in operation from 1959-1974 at the MIT Volcano Ranch station located near Albuquerque, New Mexico. One of its many distinctions was the detection of the first cosmic ray event with an energy estimated at $`10^{20}`$ eV . Over its lifetime, the array existed in three distinct configurations. The first configuration consisted of twenty scintillator detectors of surface area 3.26 $`\mathrm{m}^2`$. Nineteen detectors were spaced 442 m apart. A twentieth detector was placed at various locations and used for the measurement of the density of particles. A second configuration had a larger spacing of 884 $`\mathrm{m}`$.

It is the third configuration which is most significant for our analysis. In this configuration, the number of detectors was quadrupled by splitting up each of the twenty detectors into eighty detectors of surface area 0.815 $`\mathrm{m}^2`$ spaced 147 m apart. This configuration, with many more detectors spaced closer together, allowed for very precise measurements of the lateral distribution of signals in the detectors. The steepness of the lateral distribution, and fluctuations in it, may lead to an estimate of the primary composition.

## 3 Simulation of detector response

To simulate the detector response of the array to the ground particles, we utilized the general-purpose simulation toolkit geant4. Our procedure follows the prescription found in , where the detector response to electrons, gamma, and muons is simulated in the energy range 0.1 to 1.0e5 MeV and for five bins per decade of energy. This detector response is convolved with the results of aires air shower simulations to obtain scintillator yield in minimum ionizing particles per square meter $`(\mathrm{mips}/\mathrm{m}^2)`$.

## 4 Lateral distribution function

For Volcano Ranch data, a generalized version of the Nishimura-Kamata-Greisen (NKG) formula was used by Linsley to describe the lateral distribution of particles at ground in $`\mathrm{mips}/\mathrm{m}^2`$ . This lateral distribution function is given as

$$S_{vr}(r)=\frac{N}{r_m^2}C(\alpha ,\eta )\left(\frac{r}{r_m}\right)^\alpha \left(1+\frac{r}{r_m}\right)^{(\eta \alpha )}$$
(1)
normalized to $`N`$ with

$$C=\frac{\mathrm{\Gamma }(\eta \alpha )}{2\pi \mathrm{\Gamma }(2\alpha )\mathrm{\Gamma }(\eta 2)}.$$
(2)
Here $`r_m`$ is the Moliere radius, which is $``$ 100 m at the Volcano Ranch elevation.

From a subset of showers detected by the array in the third configuration, the following form of $`\eta `$ as a function of zenith angle and shower size was found to fit the data :

$$\eta (\theta ,N)=b_0+b_1(\mathrm{sec}\theta 1)+b_2log_{10}(\frac{N}{10^8})$$
(3)
with $`b_0=3.88\pm 0.054`$, $`b_1=0.64\pm 0.07`$, and $`b_2=0.07\pm 0.03`$.

## 5 Preliminary comparison of lateral distribution

Here we compare the Volcano Ranch lateral distribution with modern Monte Carlo simulations. In Figure 1 we show how well the simulations (aires/qgsjet98 with the detector response included) reproduce the average lateral distribution of 707 showers of estimated size $`N=10^8`$, measured by the array in its first configuration .

As a further comparison, for fixed bins of zenith angle, we determine the number of particles at ground level from a fit to the lateral distribution function (with $`\alpha =1`$) and compare to the average functional form of $`\eta `$ given by Equation 3. The results of this comparison for $`\mathrm{sec}\theta =1.11.2`$ can be seen for proton and iron showers using aires/qgsjet98 and aires/sibyll2.1 in Figures 2 and 3. One can see that the average form of $`\eta `$ over a realistic range of mass composition, from proton to iron, is reproduced by the simulations. At this stage, any attempt to draw a firm conclusion about composition or the suitability of interaction models is premature; we have limited Monte Carlo statistics and are only comparing to an average functional form of $`\eta `$. What we are content to take from this stage of the analysis is encouragement from the fact that the simulation can reproduce the data. A more detailed analysis will follow in a later work.

## 6 Towards an estimate of primary composition

A more robust estimator of primary composition comes from an analysis of the measured fluctuations of $`\eta `$, which were found to be significantly greater than that expected from measurement error alone . From our knowledge of the characteristics of showers produced by nucleons and nuclei, one would expect that measured fluctuations in the distribution of $`\eta `$ would be smaller for a heavier composition. This can be understood if we consider an air shower produced by a nucleus with $`A`$ nucleons as a superposition of $`A`$ showers each with energy $`E/A`$. Thus in iron-initiated showers the average fluctuations in shower development are reduced.

At Volcano Ranch, measurements of $`\eta `$ were made on a shower-to-shower basis for fixed bins of zenith angle . Our procedure is to simulate showers in these bins of zenith angle, find $`\eta `$ from a fit to the lateral distribution function, and compare to data. The results of this preliminary comparison can be seen in Figures 4 and 5.

## 7 Conclusions and future plans

Using modern Monte Carlo tools previously unavailable during its operation we have been able to reproduce lateral distribution measurements taken at the Volcano Ranch array. With an eye towards estimating the primary composition from these measurements, our ability to reproduce them gives us some confidence that our analysis procedure is correct and indicates that current hadronic interaction models can describe the data fairly well. A preliminary analysis of the distribution of the steepness of the lateral distribution function, $`\eta `$, indicates that indeed something can be said about primary composition from an analysis of the data, but that further analysis is needed to make any stronger statement. However, it is encouraging that the data from Volcano Ranch may be used to say something useful about mass composition thirty years after its closure.