text:
# The accretion flow in the discless intermediate polar V2400 Ophiuchi

## 1 Introduction

The magnetic cataclysmic variables can be divided into two classes. In polars (or AM Her stars) the magnetic field of the white dwarf is strong enough to lock its rotation to the binary orbit, and accretion proceeds via a stream which is deflected by the field onto a magnetic pole. In contrast, in intermediate polars (IPs or DQ Her stars) a lower-field white dwarf spins more rapidly than the orbit, and the stream feeds into an accretion disc, which then feeds field lines from its inner edge (see Warner 1995 for a comprehensive review).

The distinction is blurred, however, by the existence of several asynchronous polars, in which the spin periods differ from the orbit by $``$1 per cent (e.g. Schwope et al. 1997). There has also been a long debate on the existence of discless IPs (e.g. Hameury, King & Lasota 1986; King & Lasota 1991; Hellier 1991; Wynn & King 1992). In such systems, the stream would be expected to flip first to the upper pole, then to the lower pole, as the magnetic dipole rotates (e.g. Hellier 1991; Wynn & King 1992). Such pole flipping would occur at the frequency with which the relative geometry changes, namely the beat frequency $`\omega `$ – $`\mathrm{\Omega }`$, where $`\mathrm{\Omega }`$ and $`\omega `$ are the orbital and spin frequencies respectively.

Buckley et al. (1995; 1997) reported the first secure evidence for a discless IP, with the discovery of the Rosat source V2400 Oph (RX J1712.6–2414). They found that polarised light from the system varies at 927 s, which is interpreted as the spin period of the magnetic dipole and thus of the white dwarf. The X-rays, though, are pulsed at 1003 s, which is the beat period between the 3.42-h orbital and 927-s spin cycles; thus V2400 Oph shows the signature of pole-flipping accretion, as expected in a discless IP.

V2400 Oph is a prime opportunity to study the interaction of an accretion stream with a magnetic field in a situation where the flow is continually changing as the dipole rotates. This contrasts with most magnetic cataclysmic variables where the flow is expected to settle (at least temporarily) into a quasi-equilibrium. (In principle the asynchronous polars offer the same opportunity, but it is much harder to obtain good coverage of their $``$50-d beat cycles.) In this paper we present new X-ray observations of V2400 Oph, along with optical spectroscopy aimed at tracing the accretion flow as it connects to the field lines.

V2400 Oph is the most polarised IP, leading to a magnetic field estimate of 9–27 MG (Buckley et al. 1995; Väth 1997), the highest for any IP and overlapping with the range for low-field polars. The high field may explain why a disc has not formed in this system, whereas they do in IPs at $``$1 MG. Buckley et al. (1995) also found that the circular polarisation is always of the same sign, and so concluded that we only ever see one magnetic pole. Given, also, a lack of radial velocity motion at the orbital period, they proposed that V2400 Oph is at a low inclination, and that we only see the pole nearest us (the ‘upper’ pole). From the fact that there were no detectable variations in linear polarisation, Buckley et al. (1995) suggested that this pole was never viewed side on, and so proposed that both the inclination of the binary, and the angle between the spin and magnetic axes of the white dwarf, $`\delta `$, are low.

## 2 An RXTE observation

We report, first, on an RXTE X-ray observation obtained over the interval 2000 June 24–25. This resulted in 39 ks of data from the PCA instrument, with either 3 or 4 of the 5 PCU modules in operation at any time. A section of the resulting lightcurve is shown in Fig. 1, while the Fourier transform of the entire observation is shown in Fig. 2. Both lightcurve and transform are dominated by a pulsation at the 1003-s beat period; this has a maximum amplitude of $``$50 per cent \[(peak–trough)/peak\] and a mean amplitude of 25 per cent, as shown by the folded pulse profile in Fig. 3. There is no sign of variations at the orbital period ($`<`$ 7 per cent) or at the spin period ($`<`$ 2 per cent). A preliminary look at five other RXTE observations, of similar length and spaced at intervals of a few months, finds the same behaviour: a strong beat-cycle pulse but no modulation at the orbital or spin cycles.

The absence of an orbital modulation in the X-ray lightcurve is unexpected in a discless IP. If the stream flows to the magnetosphere, latches onto the nearest field line, and follows it onto the white dwarf (see Fig. 4), then the accretion sites will always be on the hemisphere facing the secondary, and thus their visibility should vary over the orbit, as in an AM Her star. The lack of an orbital modulation can, though, be explained if we concur with the suggestion that V2400 Oph is at so low an inclination that the upper pole is always in view and the lower pole never in view.

The next puzzle arises from the fact that the beat-cycle modulation is not total. If all the accreting material participates in the pole flipping, then, given the low inclination, the flow would flip to the far, hidden hemisphere of the white dwarf for roughly half the beat cycle, leading to 100 per cent modulation. If some of the flow does not flip, but manages to feed both poles continually, then this could explain the non-total modulation. But why, then, does this not result in a spin-cycle pulsation, as seen in all the IPs where a disc ensures continual feeding of both poles? One answer could be to invoke a field aligned exactly with the spin axis, $`\delta =0`$, but this would also remove the beat pulse, since nothing would change with beat phase. It also conflicts with the fact that we see spin-cycle variations of the emission lines (see Section 7).

A second answer might be to invoke large, extended accretion footprints, coupled with a large $`\delta `$ so that these sites would be near the white-dwarf equator. In a face-on system, any accretion sites extending across the equator would always be partly in view, and, with both poles partially visible, any pole-flipping modulation would not be total. Arguing against this, however, is the fact that the visibility of the poles would then be exquisitely sensitive to orbital phase, and we do not see any orbital modulation in the X-ray flux. Also, the fact that we see variable circular polarisation but not variable linear polarisation argues that we do not see the poles sideways on (Buckley et al. 1995).

A final possibility could invoke a field geometry that is more complex than a dipole, in which both poles were always on the near, visible hemisphere. However, this is contradicted by the fact that the circular polarisation always has the same sign, arguing that we only see one pole (Buckley et al. 1995). Thus, the absence of an X-ray spin pulse remains a puzzle, which we defer to Section 8, where we suggest a possible explanation.

## 3 Spectroscopic observations

We observed V2400 Oph with the 3.9-m AAT and the RGO spectrograph plus a TEK CCD. A 1200 lines mm-1 grating gave a resolution of 1.4Å, covering the range H$`\gamma `$ to H$`\beta `$. Observing for 5.0 h, 4.5 h & 4.5 hr on the three consecutive nights 1996 May 10–12 we obtained 900 integrations of 50 s each, thus covering $``$4 orbital cycles and $``$50 spin cycle of the star. The summed spectrum, containing He i and He ii lines in addition to the Balmer lines, is shown in Fig. 5.

As a first look at the data we computed the equivalent widths and the V/R ratios for the lines in each spectrum (V/R being the ratio of the equivalent widths on either side of the rest wavelength). The Fourier transforms of these quantities for H$`\beta `$ and He ii$`\lambda `$4686 are shown in Fig. 6.

All four Fourier transforms show significant power at the beat ($`\omega `$ – $`\mathrm{\Omega }`$) frequency, which is always stronger than the power at the spin frequency ($`\omega `$). A harmonic of the beat frequency, 2($`\omega `$ – $`\mathrm{\Omega }`$), is also prominent in the equivalent widths. Other sidebands detected include $`\omega `$ – 2$`\mathrm{\Omega }`$ in the H$`\beta `$ V/R ratios and 2$`\omega `$ – $`\mathrm{\Omega }`$ in the He ii$`\lambda `$4686 equivalent widths. There is significant power at the orbital frequency in the V/R ratios (a one-day alias being the highest peak) but not in the equivalent widths.

To help us interpret the line profiles of V2400 Oph we have computed simulations of the profiles expected in a discless IP. We first describe this model, and then present a comparison of simulated and observed profiles.

## 4 Model line profiles

The model we use to simulate the line profiles of V2400 Oph is illustrated in Fig. 4. We assume a dipole field centered at the white-dwarf centre — observations of AM Her stars often indicate the need for a more complex field (e.g. Wickramasinghe & Ferrario 2000) but we know less about the fields in IPs, so adopt the simplest approach. We assume that at some radius, $`R_{\mathrm{mag}}`$, the stream attaches to the nearest field line and follows this line, taking the gravitationally downhill direction, onto the white dwarf. (Strictly, the field strength, and thus $`R_{\mathrm{mag}}`$, would be a function of the dipole orientation, but how the material feeds onto the field lines, whether it passes through a shock, and how much the rate of feeding depends on dipole orientation are all uncertain, so again we adopt the simplest approach.)

The radius $`R_{\mathrm{mag}}`$, the angle $`\delta `$, and the inclination of the binary are all free parameters. We assume, arbitrarily, that the stream feeds a 20 swathe of azimuth. We then calculate the free-fall velocities along the field lines and project them onto the line of sight to produce line profiles. We assume that the emitting regions are optically thick, and so scale the line intensity by projected area. Also, for maximum simplicity, we assume that the field lines are not distorted by their interaction with the flow, and that the material starts infalling with near-zero velocity at $`R_{\mathrm{mag}}`$, as though it passes through a shock. We also assume that the infall time is short compared to the spin cycle, and thus do not allow for the fact that the stream will have moved on by the time the material accretes. Lastly, we rotate the dipole and stream location to mimic the spin and orbital cycles. We do not include any emission from the stream itself, further out than $`R_{\mathrm{mag}}`$, which would produce a lower-velocity S-wave on the orbital cycle. Note that the above model is very similar to that developed by Ferrario & Wickramasinghe (1999), the main difference being that they aimed for model Fourier transforms whereas we aim to compare line profiles.

Fig. 7 shows a sample set of simulated line profiles. We adopted spin and orbital cycles of 927 s and 3.42 h, a white dwarf mass of 0.80 M, and a secondary mass of 0.37 M. The inclination was 10, the dipole-offset, $`\delta `$, was 60 and the radius $`R_{\mathrm{mag}}`$ was 20 $`R_{\mathrm{wd}}`$. The four panels are for emission at different distances from the white dwarf.

A velocity variation at the orbital cycle is obvious in all panels of Fig. 7. Note that this is primarily the infall velocity of the stream, as in an AM Her star, and is not the orbital velocity of the white dwarf, which is much smaller. The spin and beat cycles are seen as the much faster variation in the simulated profiles.

## 5 Line profiles over the orbital cycle

The Fourier transforms (Fig. 6) revealed a variation of the line profiles over the orbital cycle, but when the observed profiles are folded on the orbital cycle it is hardly visible (Fig. 8). To enhance the variation we have subtracted the phase-invariant profile (for each velocity we found the phase bin with the lowest value and subtracted that from the data). In the portion of the line remaining, only 5 per cent of the original, the slight variation is discernable. Since the orbital modulation is so small, it is clear that V2400 Oph must be at a very low inclination (as already suggested by the X-ray data). But even so, the lack of variation is clearly discrepant with the model, which is calculated for $`i=10`$.

By fitting Gaussians to the phase-resolved H$`\beta `$ profiles we measure a radial-velocity variation of only 5.8 km s-1, which agrees with Buckley et al. (1995) who report a limit of $`<`$ 10 km s-1. If this is the motion of the white dwarf it implies an inclination of only 3 (assuming stellar masses of 0.7 and 0.3 M). Even this could be too high, since, as shown by the simulations, we expect the orbital motion to be dwarfed by an infall velocity of anything up to the white-dwarf escape velocity. The observations and simulations can be reconciled only if 1) the system is at an inclination of $`<`$ 1 (unlikely, since the a priori probability is only 1.5$`\times 10^4`$), or 2) much of the emission comes from material that has circularised about the white dwarf and dilutes the emission from the infalling stream. Thus the variations in the emission lines might be tracing only a part (possibly a small one) of the flow.

## 6 Line profiles over the beat cycle

Fig. 9 shows the simulated profiles folded on the beat period, again for four distances from the white dwarf. The changes in velocity as the stream changes direction to flow to a different pole are obvious. Note that, at such a low inclination, the dominant component of the velocity is the motion out of the plane following the arch of the field lines. At phase zero the stream is flowing to the upper pole, and the emission 10–20 $`R_{\mathrm{wd}}`$ out is blueshifted, since it is rising out of the plane towards us; the emission 1–10 $`R_{\mathrm{wd}}`$ out is falling back to the plane and so is redshifted.

In Fig. 10 we show the observed profiles of H$`\beta `$ and He ii$`\lambda `$4686 folded on the beat cycle. Again, we have first subtracted the phase-invariant profile (leaving only 5 and 10 per cent of the two lines, respectively), since the variation is barely discernable in the raw fold. Accompanying the data are the simulated profiles from 5–10 $`R_{\mathrm{wd}}`$ out, which give a rough match. The similarity of the H$`\beta `$ line to the simulation, with both showing changes from red shift to blue shift, suggests that we are seeing the change in velocity as the stream flips from pole to pole. In the observations, however, the flip is less abrupt than in the model, perhaps indicating that the stream feeds field lines over a larger range of azimuth than the 20 adopted in the model. The He ii$`\lambda `$4686 line does not match so well, and shows little variation in velocity. Possibly this line comes mostly from near the apex of the arch, where motion out of the plane is minimal.

Our model predicts that the line intensity is modulated at twice the beat frequency (as did the model of Ferrario & Wickramasinghe 1999), and this agrees with the presence of a 2($`\omega `$ – $`\mathrm{\Omega }`$) modulation in the equivalent widths (Fig. 6). The weaker modulation seen in the equivalent widths at $`\omega `$ – $`\mathrm{\Omega }`$ is not predicted by our model, which has symmetric upper and lower poles, but would result from any asymmetry between the poles. Another factor is that X-ray illumination of the stream will vary with beat phase, and this would also modulate the equivalent widths at the beat period.

## 7 Line profiles over the spin cycle

Fig. 11 shows the simulated profiles folded on the spin cycle. The symmetric shape is the result of the two opposite poles, but note that at any one time only one is accreting. Thus a ‘snaphot’ of the simulation would show an S-wave, snaking from red to blue over the spin cycle; after a pole-flip the inverse S-wave would be seen, and such S-waves accumulate to produce the symmetric pattern when a long sequence is folded.

Again, the dominant velocity is the motion out of the plane, so that the profiles 10–20 $`R_{\mathrm{wd}}`$ out (where material rises out of the plane) tend to show the reverse of the profiles from 1–10 $`R_{\mathrm{wd}}`$ (where material falls back). If there were no effect of inclination the profiles would repeat exactly twice per spin cycle, but the tilt to an inclination of 10 is enough to break the symmetry.

Fig. 12 shows the spin-folded H$`\beta `$ and He ii$`\lambda `$4686 profiles, again accompanied with the simulated profiles from 5–10 $`R_{\mathrm{wd}}`$, which give a good match to the He ii$`\lambda `$4686 line. In the simulation the upper pole is at the “3 o’clock” position at phase zero and rotates counterclockwise. At phases 0.2–0.3 we are seeing the upper pole furthest from us and the lower pole towards us. The material at 5–10 $`R_{\mathrm{wd}}`$ is beginning to fall back to the plane. The 10 inclination tilts this motion towards the plane of the sky, reducing the projected velocities, and the profile is thus narrowest. Half a cycle later (upper pole towards us), the inclination tilts the infall motion towards the line-of-sight, and the observed velocities are higher.

The profiles from H$`\beta `$ are noisier than those from He ii$`\lambda `$4686: after subtraction of the phase-invariant profile only 3 per cent of the original H$`\beta `$ line remained. However, the profile can be interpreted as being similar to that of He ii$`\lambda `$4686, with the exception that the blue wing seen at phases 0.4–0.8 is missing.

### 7.1 Varying the model parameters

In the previous sections we have compared observed line profiles with models calculated for emission 5–10 $`R_{\mathrm{wd}}`$ from the white dwarf. We should emphasise that we do not regard the models as fits, since the model contains too many assumptions and adjustable parameters to lead to a unique solution. Instead, we describe here the result of varying the parameters. Since the simulation calculates infall velocities, changing $`R_{\mathrm{mag}}`$ or the white-dwarf mass simply scales the velocities without altering the shape of the model profiles. Similarly, adding an injection velocity at $`R_{\mathrm{mag}}`$ would increase the velocities at all points. Increasing the 20 range of azimuth over which accretion occurs would smear out the profiles, but retain the overall shape.

However, the shape of the spin-resolved profiles is affected by the dipole offset, combined with the inclination. At zero inclination, of course, our view would not change with spin phase and there would be no spin-cycle modulation. This rules out the extreme inclinations of $`<`$ 3 discussed in Section 5 — only by adopting an inclination of at least 10 do we reproduce the observed spin-resolved profiles, in which the narrow profile at phase $``$0.2 changes to a broader, split profile at phase $``$0.7. Even with an inclination of 10, the profile is only reproduced with $`\delta `$$`>`$ 30: any lower and the change in viewing angle over the spin cycle is insufficient to reproduce the observed profiles. Of course, we can also reproduce the spin-resolved profiles using a higher inclination, but that only increases the problem of the lack of orbital modulation.

## 8 Discussion

We have shown that the emission lines in V2400 Oph vary on both the spin and the beat cycles, and that the variations are compatible with a model in which the accretion flow flips from one magnetic pole to the other over the beat cycle — the hallmark of discless accretion in an asynchronous magnetic cataclysmic variable. We have found that a simple model of this behaviour, in which a stream falls to the magnetosphere and then feeds the nearest field line, reproduces the observed line-profile changes over the spin and beat cycles. However, it should be noted that the line-profile variations involve only $``$5 per cent of the line emission, and that the bulk of the emission shows no detectable variation.

As a result of this, the variations over the orbital cycle are harder to interpret. We adopt an inclination of 10, since this is the minimum necessary to reproduce the line-profile variations over the spin cycle (Section 7). However, if all the emission came from a collimated, infalling stream, this inclination would lead to an orbital modulation of hundreds of km s-1, whereas the observed value is only 5.8 km s-1 (Section 5). This implies that the bulk of the emitting material is not in the infalling stream, but is instead symmetric about the white dwarf.

This interpretation is supported by X-ray data, in which the beat-cycle modulation is only 25 per cent deep. Given that the system is at so low an inclination that we only ever see the upper pole, this suggests that only a minority of the accreting material flips from pole to pole, and that there is a continual flow to both poles irrespective of beat-cycle phase.

What is the nature of this flow? One possibility might be a conventional accretion disc. This could feed field lines from its inner edge, forming ‘accretion curtains’ of material in the usual manner for an IP. However, why, then, do we see no X-ray spin pulse? Given that the inclination and dipole offset are sufficient to produce obvious spin-cycle variations in the emission lines, we would expect the aspect of the accretion sites to vary with spin phase, and thus absorption and scattering of X-rays in the accretion curtain would produce a spin pulse (this is best documented in AO Psc; Hellier, Cropper & Mason 1991; Hellier et al. 1996; see also the model of Kim & Beuermann 1995). Another issue is how the stream coexists with the disc, since, having different velocities, they couldn’t both occupy the orbital plane. One plausible answer is that the stream could overflow the disc (see Armitage & Livio 1998 for theoretical simulations of this). This idea has been invoked to explain the IPs showing both spin and beat pulses in the X-ray lightcurves (e.g. Hellier 1991). However, again, a spin-cycle pulsation is nearly always seen, and, taking the example of FO Aqr, the beat-cycle modulation is always weaker or even absent (Norton et al. 1992; Hellier 1993b; Beardmore et al. 1998).

In view of the difficulties in invoking a disc, we turn instead to the idea of magnetically threaded accretion flows, developed by King (1993) and Wynn & King (1995). This model treats the flow as diamagnetic blobs, which are diverted from a ballistic trajectory by the magnetic drag produced by their crossing field lines. The resulting trajectories depend on blob size and density. Simulations by King & Wynn (1999) suggest that it is possible for less-dense matter to be easily threaded and controlled by field lines (which would then act as in our simulations), while denser blobs cross field lines and circulate around the white dwarf. These blobs would lose knowledge of orbital phase, and eventually be destroyed by magnetohydrodynamic instabilities, feeding their material onto the white dwarf’s poles.

Back to the same question, why doesn’t this result in an X-ray spin pulse? We speculate that this process is much less orderly than in a disc-fed accretor. A disc might supply a steady flow over a relatively small disc–magnetosphere transition region, and so feed accretion sites that are small and localised. For example they cover only $`<`$ 0.002 of the white-dwarf surface in XY Ari (Hellier 1997b). This would lead to steady, well defined accretion curtains, able to absorb or scatter the X-ray flux and so modulate it over the spin cycle. In the diamagnetic-blob model, the accretion sites could be widely strewn and transient; they may not form coherent accretion curtains and so the spin-cycle modulation could be reduced to below the detection threshold.

### 8.1 Relation to other systems

We have detected variations of the line profiles over the spin frequency in V2400 Oph, whereas hitherto the spin frequency in this star has been seen only in circular polarimetry. This is reassuring, for the following reason. Most IPs show little or no polarised light, and in such systems an X-ray pulse is typically interpreted as the spin period. But if it was actually the beat period, how could we know? The fact that both spin and beat frequencies show up in the line profiles of a discless accretor like V2400 Oph means that we can distinguish between the spin and beat frequencies even without polarimetry. This means that we can unambiguously determine spin frequencies in the class as a whole, since nearly all well-studied IPs show both spin and beat cycles in the emission lines \[e.g. Hellier (1997a) and references therein, an exception being TV Col, where neither cycle is convincingly detected in the optical (Hellier 1993a; Augusteijn et al. 1994)\].

Secure identification of the spin and beat frequencies strengthens the classification of IPs using the idea that disc-fed accretors show X-ray pulses (predominantly) at the spin cycle while discless accretors show X-ray pulses (predominantly) at the beat cycle. Based on this, we would assign V2400 Oph as the sole secure member of the discless category, and a dozen others to the disc-fed category (e.g. Hellier 1991; 1998). An unclear case is TX Col, in which the spin pulse dominates at some epochs and the beat pulse at others (Norton et al. 1997; Wheatley 1999).

V2400 Oph has a 3.42-h orbit and is the highest-field IP at 9–27 MG. It is likely that this system is on the edge of disc formation since we can account for our observations only if there is a substantial ring-like structure of blobs circling the white dwarf. If the field were a little lower, or the orbit slightly wider (so that the angular momentum of the blobs caused them to orbit further out), the ring of blobs might turn into an accretion disc. Note that PQ Gem, the only IP with a comparable field estimate (9–21 MG; Väth et al. 1996; Potter et al. 1997) has a wider orbit of 5.19 h (Hellier 1997a). It appears to accrete though a disc, judging by the dominance of the X-ray spin pulse (Mason et al. 1992). It is thus likely that, with these two systems, we are delineating the conditions for disc formation in an IP.

The theoretical requirements for disc formation have long been debated (e.g. Hameury et al. 1986; Lamb & Melia 1988). The issue centers on the three distances $`R_{\mathrm{min}}`$ (the radius of closest approach to the white dwarf of a ballistic stream from the Lagrangian point), the larger circularisation radius, $`R_{\mathrm{circ}}`$ (where orbiting material has the same angular momentum that it had at the Lagrangian point), and $`R_{\mathrm{inner}}`$ (the radius at which the disc is magnetically disrupted). \[We distinguish $`R_{\mathrm{inner}}`$ from the $`R_{\mathrm{mag}}`$ used above since the radius is different for spherical, stream, and disc flows; see Warner (1995) for a fuller discussion and formulae for these radii.\]

It is easy to form a disc with $`R_{\mathrm{inner}}`$$`<`$$`R_{\mathrm{min}}`$, since the stream could freely orbit outside the magnetosphere, circularise at $`R_{\mathrm{circ}}`$, and form the disc. It is not possible to form a disc with $`R_{\mathrm{inner}}`$$`>`$$`R_{\mathrm{circ}}`$ since the stream has insufficient angular momentum to orbit at these radii. The intermediate case, $`R_{\mathrm{circ}}`$$`>`$$`R_{\mathrm{inner}}`$$`>`$$`R_{\mathrm{min}}`$, is less clear, but is the case applying to most IPs. For example, if we assume that PQ Gem has a red-dwarf mass of 0.5 M (appropriate for its 5.19-h orbital period), a white-dwarf mass of 0.7 M, and that its magnetosphere corotates with the Keplerian orbit at $`R_{\mathrm{inner}}`$, then, given the 833-s spin period, we have $`R_{\mathrm{inner}}`$$``$1.1$`\times 10^{10}`$ cm, which is larger than $`R_{\mathrm{min}}`$ ($``$0.6$`\times 10^{10}`$ cm) and close to $`R_{\mathrm{circ}}`$ ($``$1.1$`\times 10^{10}`$ cm).

V2400 Oph could be showing us the mechanism for disc formation under these circumstances. In this discless system a ring of diamagnetic blobs orbits inside the magnetosphere. In doing so, the blobs cross field lines and lose angular momentum to the field, and so eventually accrete. However, if the field were weaker, as in most IPs, the drag on the diamagnetic blobs would be lower, and their inward spiral would be slower. The blobs might accumulate, screening the field from each other, and so accumulate further in a runaway process leading to the formation of a disc.

## 9 Conclusions

$``$ V2400 Oph is at a low inclination, estimated as 10. We only ever see the upper magnetic pole.

$``$ There is no accretion disc in this system. The stream-fed accretion flips from pole to pole on the beat cycle, and this can be seen directly in the line profiles. However, only $``$25 per cent of the accreting material participates in this motion.

$``$ The rest of the material appears to be circling the white dwarf, perhaps in the form of diamagnetic blobs. Such blobs may have been a precursor to disc formation in other intermediate polars.

$``$ The spin-cycle variations of the emission lines are accounted for by a simple model of stream-fed accretion. The best match between model and data is found for ditances 5–10 $`R_{\mathrm{wd}}`$ from the white dwarf. The offset between the magnetic and spin axes of the white dwarf must be at least 30.


format: text

text:
# Enhancement of Ambipolar Diffusion Rates through Field Fluctuations

## 1 Introduction

In the usual paradigm of low mass star formation, molecular cloud cores are supported by magnetic fields. In order for star formation to take place, the cores must lose magnetic support, and this loss of support is generally thought to take place through the action of ambipolar diffusion (Mouschovias 1976; Shu 1983; Nakano 1984; Shu, Adams, & Lizano 1987; Lizano & Shu 1989; Ciolek & Basu 2000, 2001). This general picture has support from observations, which suggest that ion-neutral drift does indeed occur in magnetized star-forming cores (e.g., Greaves & Holland 1999).

An important issue facing this standard scenario is the time scale required for magnetic support to be removed from the cloud cores. As the observational picture comes into sharper focus, the number of observed cores without stars (e.g., Jijina, Myers, & Adams 1999) seems to be smaller than that predicted by most previous estimates from ambipolar diffusion (e.g., Ciolek & Mouschovias 1994, 1995; Lizano & Shu 1989) by a factor of 3 – 10. In other words, loss of magnetic support by diffusion appears to be too slow, with a time scale a factor of 3 – 10 times longer than suggested by the observed statistics of cloud cores. However, these previous calculations neglected a dimensionless factor that depends on the mass to flux ratio of the cores (Ciolek & Basu 2001). If the cloud cores have mass to flux ratios that approach the critical value, then the ambipolar diffusion time scale is significantly shorter than previous estimates. In particular, if the mass to flux ratio becomes supercritical, then the ambipolar diffusion time scale approaches zero. The need for this correction is bolstered by a recent compilation of Zeeman measurements of magnetic field strengths (Crutcher 1999), which suggests that many cores may have mass to flux ratios near the supercritical value. This observed sample includes only 27 cores with relatively large masses; additional measurements are necessary to clarify the observational picture.

In this study, we consider the effects of fluctuations on the mechanism of ambipolar diffusion described above. The time scale issue remains important and this work shows that ambipolar diffusion can operate more quickly in the presence of such fluctuations. In addition, because of the chaotic nature of the fluctuations, the ambipolar diffusion time scale will take on a full distribution of values for effectively “the same” initial states.

Fluctuations are expected to be present in essentially all star forming regions. Molecular clouds are observed to have substantial non-thermal contributions to the observed molecular line-widths (e.g., Larson 1981; Myers, Ladd, & Fuller 1991; Myers & Gammie 1999). These non-thermal motions are generally interpreted as arising from MHD turbulence (e.g., Arons & Max 1975; Gammie & Ostriker 1996; for further evidence that the observed linewidths are magnetic in origin, see Mouschovias & Psaltis 1995). Indeed, the size of these non-thermal motions, as indicated by the observed line-widths, are consistent with the magnitude of the Alfvén speed (e.g., Myers & Goodman 1988; Crutcher 1998, 1999; McKee & Zweibel 1995; Fatuzzo & Adams 1993). As a result, the fluctuations are often comparable in magnitude to the mean values of the fields (T. Troland, private communication).

Background fluctuations can lead to a net change in the diffusion rate because magnetic diffusion is a nonlinear process. As many authors have derived previously (e.g., see the textbook treatment of Shu 1992), and as we present below, the (dimensionless) diffusion equation takes the schematic form

$$\frac{b}{\tau }=\frac{}{\mu }\left(b^2\frac{b}{\mu }\right),$$
(1)
where $`b`$ is the magnetic field strength, $`\mu `$ is the Lagrangian mass coordinate, and we have ignored density variations. Now suppose that the magnetic field fluctuates about its mean value on a time scale that is short compared to the diffusion time (the time required for the mean value to change). We thus let $`bb(1+\xi )`$, where $`\xi `$ is the relative fluctuation amplitude. In the simplest case in which the fluctuations are spatially independent, the right hand side of equation  is thus multiplied by a cubic factor $`(1+\xi )^3`$. Although a linear correction would average out over time, this nonlinear term must always average to a value greater than unity and the corresponding diffusion time scale grows shorter by the same factor. As an example, suppose the field spends half of its time at a value of twice its mean strength and the other half of its time near zero strength. Half of the time, the effective diffusion constant is thus larger by a factor 8, whereas the other half of the time, the diffusion constant is effectively zero. In this naive example, the mean diffusion constant is thus 4 times larger due to the fluctuations. The goal of this paper is to derive a more rigorous argument for this time scale enhancement.

This effect – changing diffusion time scales because of fluctuations – is well known in mathematical subfields. Simpler problems in which random noise fields drive physical systems at different rates appear in a host of textbooks (e.g., Srinivasan & Vasudevan 1971; Soong 1973). More recently, in the context of “stochastic ratchets”, it has been shown that random fluctuations can drive a physical system to propagate “uphill”, i.e., the opposite direction of its natural propagation in the absence of fluctuations (Doering, Horsthemke, & Riordan 1994). In astrophysics, stochastic aspects of magnetic field fluctuations have been considered in the context of cosmic ray propagation (e.g., Jokipii & Parker 1969) and also in stellar atmospheres (e.g., Shore & Adelman 1976). Several previous papers have studied turbulent fluctuations in magnetically supported clouds, often by considering how the turbulence itself leads to field evolution in the absence of ambipolar diffusion (e.g., Kim 1997). The formation of cores through the dissipation of turbulence has been suggested (Myers & Lazarian 1998). Turbulence can also enhance the rate of ambipolar drift and may help explain the observed relationship between density and magnetic field strength, $`B\rho ^\kappa `$ (see Zweibel 2001). In this current work, we consider ambipolar diffusion to be the main process that forms molecular cloud cores and study how fluctuations alter its rate.

This paper is organized as follows. In §2, we reformulate the ambipolar diffusion calculation in a plane geometry, where we explicitly include fluctuations in both the magnetic field and the density field. We perform an analysis of the resulting set of equations in §3. We specialize to the limit of long wavelength fluctuations and apply the resulting formalism to astrophysical systems. We conclude, in §4, with a summary and discussion of our results. The case of short wavelength fluctuations is presented briefly in an Appendix. The most important outcome of this study is to demonstrate that fluctuations can lead to more rapid diffusion of magnetic fields in star forming regions and that the diffusion time scale takes on a distribution of values (rather than a single time scale).

## 2 Ambipolar Diffusion in the Presence of Fluctuations

In this section, we modify the standard ambipolar diffusion derivation to include the effects of fluctuations. We consider the simplest case of a planar layer of molecular cloud material. The magnetic field lines are parallel to the plane and all quantities depend only on the height $`z`$ above the midplane,

$$𝐁=B(z)\widehat{x}.$$
(2)
We note that a host of two-dimensional axisymmetric calculations have already been done (Lizano & Shu 1989; Mouschovias & Morton 1991; Ciolek & Mouschovias 1994; Basu & Mouschovias 1994; Basu 1997, 1998; Ciolek & Basu 2000). Although more idealized, this one-dimensional calculation is useful because it allows for analytical results and isolates the fluctuation effects that we are trying to elucidate (see also the discussion at the end of this section).

As derived previously (e.g., Shu 1992), with this basic configuration the equations of motion take the form

$$\frac{\rho }{t}+\frac{}{z}(\rho u)=0,$$
(3)
$$\frac{u}{t}+u\frac{u}{z}=g\frac{1}{\rho }\frac{}{z}\left(P+\frac{B^2}{8\pi }\right),$$
(4)
$$\frac{g}{z}=4\pi G\rho ,$$
(5)
and finally

$$\frac{B}{t}+\frac{}{z}(Bu)=\frac{}{z}\left(\frac{B^2}{4\pi \gamma \rho \rho _i}\frac{B}{z}\right).$$
(6)
To close the system of equations, the ion mass density $`\rho _i`$ must be specified. The ion population depends on the complex balance between the ionization rate of the neutrals (primarily through cosmic rays), the subsequent production of molecular and metal ions, and the ion-electron recombination rate in the presence of grains. A reasonable (standard) approximation for the ionic mass density in molecular cloud environments is given by

$$\rho _i=C\rho ^{1/2},$$
(7)
where $`C=3\times 10^{16}`$ cm-3/2 g1/2. However, a new observational study (Caselli et al. 2001) indicates that the ionization fraction can be significantly smaller than suggested by equation  with this value of the constant $`C`$ (see also Ciolek & Mouschovias 1998). This reduction in ionization also acts to speed up the ambipolar diffusion rate and helps alleviate the time scale problem of interest. Notice, however, that this effect is independent of our present work and can be easily incorporated by using a different value of the constant $`C`$.

Next, we introduce the fluctuations through the ansatz

$$BB(1+\xi )\mathrm{and}\rho \rho (1+\eta ).$$
(8)
The relative fluctuations $`\xi `$ and $`\eta `$ obey distributions, which ultimately determine the effects of these variations. Although we will specialize to particular distributions later on, we now keep the analysis as general as possible. In any case, the fluctuations are assumed to have zero mean so that the quantities $`B`$ and $`\rho `$ appearing in the equations of motion can be thought of as time-averaged quantities and play the same role in this generalized calculation as they do in previous treatments.

Following many previous treatments, we rewrite the problem in terms of Lagrangian coordinates so that the basic variable is the surface density of neutrals between the midplane ($`z=0`$) and a height $`z`$. This change of variables takes the form

$$\sigma _0^z\rho (z^{},t)𝑑z^{}.$$
(9)
Notice that this Lagrangian coordinate measures the distance in terms of the mean density $`\rho `$ rather than the full fluctuating density field $`\rho (1+\eta )`$. In terms of this new coordinate, the equation of continuity becomes

$$\frac{z}{\sigma }=\frac{1}{\rho }.$$
(10)
For an isothermal equation of state ($`P=a^2\rho `$), the force equation can be written in the form

$$\frac{^2z}{t^2}=4\pi G_0^z\rho (1+\eta )𝑑z^{}+\frac{a^2}{1+\eta }\frac{}{\sigma }\left[\rho (1+\eta )\right]+\frac{1}{1+\eta }\frac{}{\sigma }\left[\frac{B^2(1+\xi )^2}{8\pi }\right].$$
(11)
Since ambipolar diffusion is generally much slower than gravitational collapse, the cloud is expected to evolve near a magnetohydrodynamic equilibrium state. The inertial term in the force equation can therefore be justifiably ignored, leaving an expression which can be integrated to yield the quasi-magnetohydrodynamic equilibrium condition. In this case, however, the integration can only be carried out after a suitable averaging of the fluctuations has been done (see below). Finally, the magnetic field evolves according to a nonlinear diffusion equation of the form

$$\frac{}{t}\left[\frac{B(1+\xi )}{\rho (1+\eta )}\right]=\frac{1}{1+\eta }\frac{}{\sigma }\left\{\frac{B^2(1+\xi )^2}{4\pi \gamma C\rho ^{1/2}(1+\eta )^{3/2}}\frac{}{\sigma }\left[B(1+\xi )\right]\right\},$$
(12)
where the partial derivatives are taken with respect to constant $`t`$ or constant $`\sigma `$.

Following the notation of Shu (1983, 1992), we introduce a dimensionless surface density $`\mu `$, volume density $`p`$, magnetic field $`b`$, vertical coordinate $`y`$, and time $`\tau `$ as follows:

$$\sigma \sigma _{\mathrm{}}\mu ,$$
(13)
$$\rho \frac{2\pi G\sigma _{\mathrm{}}^2}{a^2}p,$$
(14)
$$B4\pi G^{1/2}\sigma _{\mathrm{}}b,$$
(15)
$$z\frac{a^2}{2\pi G\sigma _{\mathrm{}}}y,$$
(16)
and

$$t\left[\frac{\gamma C}{2(2\pi G)^{1/2}}\right]\left[\frac{a}{2\pi G\sigma _{\mathrm{}}}\right]\tau .$$
(17)
The closed set of equations which describe ambipolar diffusion can now be written as a continuity equation

$$\frac{y}{\mu }=\frac{1}{p},$$
(18)
a force balance equation

$$_0^yp(1+\eta )𝑑y^{}+\frac{1}{1+\eta }\frac{}{\mu }\left[p(1+\eta )\right]+\frac{1}{1+\eta }\frac{}{\mu }\left[b^2(1+\xi )^2\right]=0,$$
(19)
and a magnetic diffusion equation

$$\frac{}{\tau }\left[\frac{b(1+\xi )}{p(1+\eta )}\right]=\frac{1}{1+\eta }\frac{}{\mu }\left\{\frac{b^2(1+\xi )^2}{p^{1/2}(1+\eta )^{3/2}}\frac{}{\mu }[b(1+\xi )]\right\}.$$
(20)
These three equations, in conjuntion with the specification of the fluctuations, constitute the problem to be solved.

As is well known, this one-dimensional slab model is unrealistic in the sense that gravity “saturates”. In other words, in Lagrangian coordinates, our model has a constant gravitational field strength for a given column density. The gravitational field depends on the mass shell but not on the distance from the midplane. Hence, a Lagrangian observer would find no increase in the gravitational field strength as the slab is compressed, unlike the more realistic case of higher dimensions. This point is important because ambipolar diffusion is driven by self-gravity and the slab model lies in the regime of “weak gravity” (Mouschovias 1982). However, the enhancement of the time scale that we study here arises from the nonlinearity of the equations and will be present in any geometry. To illustrate this claim, we derive an analogous formulation for a cylindrical geometry in Appendix A. In particular, we find that the enhancement factor (due to fluctuations) takes the same form as that derived above. Although the absolute value of the ambipolar diffusion time scale depends on the geometry, the enhancement factor calculated here does not.

## 3 Analysis

In this section, we find solutions to a limited version of the problem derived in the previous section. Unfortunately, a full solution to equations  is beyond the scope of this paper and several important simplifications are necessary to make further progress.

In a complete theory, one would also derive equations of motion for the fluctuations, represented here by $`\xi `$ and $`\eta `$, and one would solve for their temporal and spatial dependences in a self-consistent manner. Such a calculation would require a theory of how MHD waves are produced, and how they cascade into turbulence. Because we have no working a priori theory of turbulence, however, our first simplification is to assume viable fluctuation distributions rather than calculate them. In other words, we adopt a semi-empirical approach. These fluctuations are observed, and we can use observations to constrain their form, but we do not explicitly calculate their behavior.

Regardless of what distributions are chosen for the fluctuations, another complication arises. The equations of motion derived above describe the evolution of one particular realization of the problem. The usual way to solve such a problem is to set up the initial state and step forward in time (either analytically or numerically) by sampling the values of the fluctuations $`\xi `$ and $`\eta `$ from their (presumably known) distributions. As time progresses, the field evolves and one obtains a solution. However, this approach – which we will adopt – does not provide a full solution to the problem. Because the fluctuations obey a distribution, a probabilistic description of the overall evolution problem is necessary. In particular, the time evolution of the magnetic field will have a distribution of possible solutions. A full solution to the entire problem (see Doering 1990) would calculate the full distribution of possible time behaviors for the magnetic field (and other fluid variables). Because this issue is compounded for the case of short wavelength fluctuations, we focus our analysis on long wavelength fluctuations. We briefly discuss the case of short wavelength fluctuations in Appendix B, which shows that additional simplifying assumptions are necessary to make progress.

### 3.1 Time Averaging

In order to describe the average behavior of the fluid fields, we introduce a method of intermediate time averaging. Previous calculations for slab models suggest that ambipolar diffusion occurs on a time scale $`\tau _{AD}1015`$ (Nakano 1984, Shu 1983, Mouschovias 1983), although this time scale is significantly shorter when the mass to flux ratio of the cores approach supercritical values (Ciolek & Basu 2001). For comparison, observations suggest that the real value is closer to $`\tau _{AD}13`$. The fluctuations themselves presumably occur within the MHD regime and thus occur on time scales $`\tau _{\mathrm{mhd}}\lambda /v_A`$, where $`\lambda `$ is the length scale of the fluctuation and $`v_A`$ is the Alfvén speed. We can write the length scale in dimensionless form through the relation

$$\lambda =\frac{a^2}{2\pi G\sigma _{\mathrm{}}}\chi ,$$
(21)
where we expect $`\chi 110`$ for the long wavelength limit considered below (see §3.3). The dimensionless time scale for which the fluctuations vary is thus given by

$$\tau _{\mathrm{mhd}}\left[\frac{\gamma C}{2(2\pi G)^{1/2}}\right]^1\frac{a}{v_A}\chi ,$$
(22)
where the dimensionless parameter $`[\gamma C/2(2\pi G)^{1/2}]`$ $``$ 8 for our adopted values of $`C`$ and $`\gamma `$. For this approach to be consistent, the fluctuation time scale must be much shorter than the ambipolar diffusion time. This requirement implies a constraint of the form $`\chi 8(v_A/a)\tau _{AD}`$, where we have used standard values for $`\gamma `$ and $`C`$ (see also Zweibel 1988). Notice that the ratio $`v_A/a`$ is typically 5 – 10, so this constraint is usually easy to satisfy (see, however, the discussion below).

To find time-averaged quantities, we must average over an intermediate time scale $`\tau _0`$ that obeys the ordering

$$\tau _{\mathrm{mhd}}\tau _0\tau _{AD}.$$
(23)
Again, we require $`\tau _{\mathrm{mhd}}\tau _{AD}`$ so that an intermediate range of time scales exists. Further, we let brackets $`\mathrm{}`$ denote time averaged quantities so that

$$p(1+\eta )\frac{1}{\tau _0}_0^{\tau _0}p(1+\eta )𝑑\tau ^{}p,$$
(24)
and

$$b(1+\xi )\frac{1}{\tau _0}_0^{\tau _0}b(1+\xi )𝑑\tau ^{}b.$$
(25)
In equating the time averaged quantities with $`p`$ and $`b`$, we are ignoring errors of $`𝒪`$ $`[\tau _0/\tau _{AD}]`$. For the reasonable choice of taking the averaging time scale to be the geometric mean $`\tau _0=\sqrt{\tau _{\mathrm{mhd}}\tau _{AD}}`$, the relative error becomes $`𝒪`$ \[$`(\tau _{\mathrm{mhd}}/\tau _{AD})^{1/2}`$\].

### 3.2 The Quasi-Equilibrium State

Next, we need to consider the equilibrium states. The time scale over which fluctuations change (roughly given by the MHD crossing time) is much shorter than the collapse time scale (roughly given by the sound crossing time) for cloud cores undergoing ambipolar diffusion. Because the fluctuation time scale is also much shorter than the ambipolar diffusion time scale, the slab is expected to evolve in a quasi-static equilibrium state supported primarily by magnetic pressure. Since the force resulting from magnetic pressure is nonlinear, however, the fluctuations can play an important role in supporting this state. The quasi-equilibrium condition is found by integrating the time-averaged force equation  to obtain

$$Kb^2+p=1\mu ^2\mathrm{where}K\frac{(1+\xi )^2}{(1+\eta )}.$$
(26)
In this geometry (used here to make the problem more tractable), an equilibrium state exists even in the absence of magnetic fields and full gravitational collapse is not possible. While this geometry does not reflect realistic conditions in molecular cloud cores, it nicely illustrates the effects of fluctuations on the ambipolar diffusion problem.

In obtaining the above form, we have assumed that the fluctuations are relatively well-behaved so that they obey constraints of the form

$$\frac{1}{1+\eta }\frac{\eta }{\mu }0\mathrm{and}(1+\xi )\frac{\xi }{\mu }0.$$
(27)
In order for these constraints to hold, the derivatives of the fluctuations must average to zero (as expected) and the derivatives must not be correlated with the fluctuations themselves. In other words, the fluctuations must be both spatially and temporally symmetric. Ambipolar diffusion, as expressed by equation , then evolves the quasi-equilibrium state from a magnetically supported configuration to a thermally supported one in the slab model.

### 3.3 Long Wavelength Fluctuations

We now consider the effects of fluctuations on the ambipolar diffusion process. As shown below, these effects vary considerably with the spectrum of the fluctuations. To start, we assume that the fluctuations in density and the magnetic field are independent of each other. We also consider the simplest case in which the fluctuations have long wavelengths. In the limiting case where the fluctuation scale $`\lambda `$ is much larger than the typical length scale $`R`$ in the cloud, the fluctuations can be considered to be spatially independent in the diffusion equation . This regime of parameter space is defined by the constraint

$$\lambda R\frac{a^2}{2\pi G\sigma _{\mathrm{}}}.$$
(28)
For the rest of this paper, we will specialize our analysis to long wavelength fluctuations that obey this constraint.

Although nature can also support short wavelength fluctuations, a rigorous treatment of their effects is made difficult by several issues: (1) In stochastic differential equations, the solutions depend on the manner in which various limits are taken (see Doering 1990). For short wavelength fluctuations, many different equivalent ways of taking the appropriate limits are possible and no unique solution exists without further specification. In this context, we would need to understand the origin of turbulent fluctuations to provide further specification. (2) In numerical treatments of the diffusion problem, the diffusion equation can be unstable on short length scales. Many numerical treatments allow for (small) errors on short size scales, but provide the correct global behavior (see Press et al. 1986). To adequately follow short wavelength fluctuations, one needs a numerical method that adequately resolves all spatial scales. In light of these difficulties, we focus on long wavelength perturbations in this paper. For purposes of illustration, however, we consider a representative approach for short wavelength fluctuations in Appendix B.

In order for the fluctuations to affect the diffusion process, the effective diffusion coefficient must sample a range of possible fluctuations during the ambipolar diffusion time $`\tau _{AD}`$ (see equation ). If the fluctuations are MHD in origin, their time scale is given roughly by the Alfvén crossing time. Sufficiently rapid fluctuations thus impose the constraint

$$N_F\frac{\tau _{AD}}{\tau _{\mathrm{mhd}}}=\left[\frac{\gamma C}{2(2\pi G)^{1/2}}\right]\frac{v_A}{a}\tau _{AD}\chi ^11,$$
(29)
where $`\tau _{AD}`$ is the (non-dimensional) time interval required for diffusion to take place and where $`\tau _{\mathrm{mhd}}`$ is given by equation . The quantity $`N_F`$ thus represents the number of times that the distributions are independently sampled during the course of the diffusion process. The time required for the field strength to decrease by a factor of e1 $`2.7`$ is usually $`\tau _{AD}=510`$. The ratio $`v_A/a`$ is typically 3 – 10, the length scale $`\chi =110`$, and the dimensionless ratio in square brackets is about 8. For reasonable values of the parameters, the left hand side of the above inequality lies in the range $`N_F`$ = 20 – 800, a wide range but always comfortably greater than unity. Furthermore, we expect the variations from case to case to differ from the expectation values by relative sizes of order $`N_F^{1/2}`$ due to incomplete sampling of the distribution functions during the diffusion process. This relative variation is thus expected to be $`0.040.2`$.

When constraints  and  are satisfied, the fluctuations are independent of the spatial derivatives appearing in the diffusion equation . Before time averaging the diffusion equation , we rewrite the equation in a specific form to simplify the averaging procedure. In particular, we write out the time derivatives, collect the terms on the right hand side, and multiply by one factor of $`(1+\eta )`$ to obtain the form

$$(1+\xi )\frac{}{\tau }\left(\frac{b}{p}\right)+\frac{b}{p}\left[\dot{\xi }\dot{\eta }\frac{(1+\xi )}{(1+\eta )}\right]=\frac{(1+\xi )^3}{(1+\eta )^{3/2}}\frac{}{\mu }\left\{\frac{b^2}{p^{1/2}}\frac{b}{\mu }\right\},$$
(30)
where the dots represent time derivatives. With the diffusion equation written in this form, the time averaging procedure removes all fluctuations except for an effective diffusion constant on the right hand side, and the diffusion equation becomes

$$\frac{}{\tau }\left(\frac{b}{p}\right)=D\frac{}{\mu }\left\{\frac{b^2}{p^{1/2}}\frac{b}{\mu }\right\}\mathrm{where}D\frac{(1+\xi )^3}{(1+\eta )^{3/2}}.$$
(31)
Notice that this effective diffusion constant $`D`$ would be unity in the absence of fluctuations.

### 3.4 Results for particular distributions

After time averaging, the diffusion problem is a rescaled version of the one solved previously. The quasi-static equilibrium state is given by equation , which contains the factor $`K`$ due to the pressure provided by fluctuations. In the diffusion equation, the effective diffusion constant is larger by the factor $`D`$. If we rescale the magnetic field strength according to $`\stackrel{~}{b}=\sqrt{K}b`$, then the quasi-static equilibrium equation  takes its standard form. The diffusion equation also takes its standard form if we rescale the time coordinate according to $`\stackrel{~}{\tau }=(D/K)\tau `$. In other words, the ambipolar diffusion process speeds up by a factor $`\mathrm{\Lambda }=D/K`$ which is determined by the distribution of the fluctuations. To make further progress, we thus have to specify the distributions.

The magnetic field fluctuations $`\xi `$ and the density fluctuations $`\eta `$ follow normalized distributions, $`f(\xi )`$ and $`g(\eta )`$, with zero mean, i.e.,

$$f(\xi )𝑑\xi =1,f(\xi )\xi 𝑑\xi =0,g(\eta )𝑑\eta =1,\mathrm{and}g(\eta )\eta 𝑑\eta =0,$$
(32)
where $`f(\xi )d\xi `$ is the probability that a fluctuation in the magnetic field has an amplitude between $`\xi `$ and $`\xi +d\xi `$ (the function $`g(\eta )`$ is defined similarly for density fluctuations). For a given choice of distributions $`f(\xi )`$ and $`g(\eta )`$, the expectation value of the diffusion constant $`D`$, the equilibrium factor $`K`$, and the time scale correction factor $`\mathrm{\Lambda }`$ = $`D/K`$ can be calculated directly.

For purposes of illustration, we first consider the fluctuations to have uniform (flat) distributions with amplitude $`A<1`$. In other words,

$$f(\xi )=\frac{1}{2A}\mathrm{for}A<\xi <A,$$
(33)
and

$$g(\eta )=\frac{1}{2A}\mathrm{for}A<\eta <A.$$
(34)
We are implicitly assuming that the density fluctuations are independent of the magnetic field fluctuations. With this particular choice for the distribution functions, the enhancement factor can be written in the form

$$\mathrm{\Lambda }=\frac{6(1+A^2)}{3+A^2}\frac{(1A)^{1/2}(1+A)^{1/2}}{\mathrm{ln}(1+A)\mathrm{ln}(1A)}.$$
(35)
The result is plotted in Figure 1, which shows that the enhancement factor becomes substantial, $`\mathrm{\Lambda }5`$, as the fluctuation amplitude $`A`$ approaches unity. In this particular treatment, the enhancement factor has a logarithmic divergence as $`A1`$.

For comparison, if we ignore density fluctuations and consider only fluctuations in the magnetic field strength, the enhancement factor simplifies to the form

$$\mathrm{\Lambda }=\frac{1+3\xi ^2}{1+\xi ^2}=\frac{1+A^2}{1+A^2/3},$$
(36)
where the first equality holds for $`\xi ^3`$ = 0 and where the second equality holds for a uniform distribution of fluctuations. For the benchmark case $`A=1`$, the magnetic field fluctuations are comparable to the total field strength and we find a modest increase in the ambipolar diffusion rate, i.e., $`\mathrm{\Lambda }`$ = 1.5. In the extreme limit $`\xi ^2\mathrm{}`$, the enhancement factor $`\mathrm{\Lambda }3`$ (for $`\xi ^3`$ = 0). Thus, as expected, fluctuations in the magnetic field strength alone lead to smaller changes than when the density field has independent variations.

We can also consider the case of density fluctuations that depend directly on the magnetic field fluctuations. For example, flux freezing arguments (Shu 1992) imply that $`B\rho ^\kappa `$, where $`\kappa `$ = 2/3 for the usual case of spherical geometry and $`\kappa `$ = 1 for a one-dimensional cloud layer. If the fields, including fluctuations about their mean values, obey such a flux freezing relation, then we have a correlation of the form $`(1+\eta )(1+\xi )^{1/\kappa }`$ and the enhancement factor becomes

$$\mathrm{\Lambda }=\frac{(1+\xi )^{33/2\kappa }}{(1+\xi )^{21/\kappa }}.$$
(37)
For our reference case of uniform distributions with amplitude $`A=1`$, we can evaluate this expression to obtain $`\mathrm{\Lambda }=(6\kappa 2)2/[(8\kappa 3)2^{1/2\kappa }]`$. For the flux freezing exponent $`\kappa =1`$ appropriate for a one-dimensional layer, for example, we find $`\mathrm{\Lambda }1.13`$.

When the fluctuations in magnetic field strength and density are correlated, the enhancement factor is generally smaller than in the absence of correlations. The special case of $`\kappa `$ = 1/2 is particularly interesting: As shown by equation , the value $`\kappa `$ = 1/2 leads to $`\mathrm{\Lambda }`$ = 1, i.e., no net enhancement of the ambipolar diffusion rate. In other words, for this particular correlation between the density and magnetic field fluctuations, the problem reduces to its old (non-fluctuating) form. This value $`\kappa =1/2`$ is found in three-dimensional calculations of magnetic cloud models (Mouschovias 1976) and is consistent with Zeeman measurements of the magnetic field strength in several molecular clouds (Crutcher 1999; Basu 2000). However, both the theoretical calculation and observational relation correspond to the mean field and not the fluctuations; one interpretation is that the relation $`B\rho ^{1/2}`$ is more of an upper envelope than a scaling law (Zweibel 2001). In any case, the correlation of the magnetic field fluctuations with the density fluctuations thus needs to be further specified.

In general, MHD disturbances exhibit a rich variety of different possible behaviors, including the various kinds of correlations between density and magnetic field fluctuations. In the linear regime, for example, pure Alfvén waves exhibit magnetic field fluctuations but have no density variations. Purely acoustic waves have density variations but no magnetic field fluctuations. When the magnetic field is perpendicular to the direction of wave propagation, then fluids can develop magnetosonic waves with correlated magnetic and density fluctuations; this configuration also supports slow (nonpropagating) modes that compress the matter (increase density) and displace the magnetic field. In real molecular clouds, MHD disturbances are a complicated and non-linear superposition of many different types of waves and other motions (one should keep in mind that these different types of waves propagate at different speeds). As a result, many different possible distributions of fluctuations are allowed by the observations. Without further constraints, the parameter space of possible fluctuations allows arbitrarily large enhancements in the ambipolar diffusion rate. This claim is substantiated by equation , which shows that $`\mathrm{\Lambda }`$ increases without bound in the limit $`A1`$ for one particular choice of distribution. To obtain bounds on $`\mathrm{\Lambda }`$, we thus need to impose additional constraints, as discussed below.

### 3.5 Application to Molecular Cloud Cores

We can scale the parameters of this theory to observed molecular clouds in order to specify the expected regime of our derived results. As mentioned in the Introduction, molecular clouds have significant non-thermal contributions to their line-widths (Larson 1981; Myers, Ladd, & Fuller 1991) and the implied non-thermal motions are often interpreted as MHD turbulence. Indeed, the amplitude of these motions is consistent with the expected Alfvén speed in these regions.

To constrain the allowed range of our fluctuations, we consider the observed line-widths to arise from the fluctuating part of the Alfvén speed in a time-averaged sense. In other words, we assume that

$$(\mathrm{\Delta }v)_{NT}^2=v_A^2_{fluc}=\frac{B^2}{4\pi \rho }\frac{(1+\xi )^2}{1+\eta }\frac{B^2}{4\pi \rho }$$
(38)
where we have subtracted off the non-fluctuating part so that the non-thermal contribution to the line-width vanishes in the absence of fluctuations. Evaluating this expression in terms of our formalism for the fluctuations, we find

$$(\mathrm{\Delta }v)_{NT}^2=2a^2\frac{b^2}{p}\left[K1\right]=2\alpha _0a^2\left[K1\right],$$
(39)
where $`a`$ is the sound speed, $`K`$ measures the amplitude of the fluctuations as defined by equation , and $`\alpha _0`$ is the initial ratio of magnetic to thermal pressure. To evaluate the ratio $`b^2/p`$, we have used our equilibrium state with no fluctuations as a reference state. Since the observed line-widths are comparable to the expected Alfvén speeds, we expect $`(\mathrm{\Delta }v)_{NT}^22\alpha _0a^2`$ and thus we expect the quantity $`(K1)`$ to be close to unity. In other words, scaling our formulation to observed molecular clouds implies that $`K2`$.

We need to apply this result to constrain the possible values of $`\mathrm{\Lambda }`$. For fluctuations obeying uniform distributions, the parameter $`K`$ = 2 for an amplitude $`A0.886`$; the corresponding value of the enhancement factor is $`\mathrm{\Lambda }2.26`$. For the special case of no density fluctuations, the enhancement factor is given by equation ; the constraint $`K=2`$ implies that $`\xi ^2`$ = 1 and hence $`\mathrm{\Lambda }`$ = 2 (in the symmetric limit where $`\xi ^3`$ = 0). In general, the inclusion of density fluctuations makes the enhancement factor larger. For typical star forming regions, our net result is that fluctuations increase the ambipolar diffusion rate by a factor of 2.

### 3.6 Comparison to Numerical Results

We now compare our analytic results with numerical simulations. Instead of working exclusively with the expectation values ($`D`$ and $`K`$), we can use the same distributions for the fluctuations $`\xi `$ and $`\eta `$ and then sample the distributions as we numerically integrate the equations of motion. Specifically, we use a time-averaged quasi-static equilibrium state (equation ) and numerically integrate the diffusion equation  by sampling the distribution functions for the fluctuations to specify the effective diffusion constant. Our basic numerical scheme is similar to that described in Shu (1983), although the time resolution must be significantly higher (and highly variable) to properly include the wide range of values for the effective diffusion constant. Each numerical simulation results in one particular realization of the time evolution. If we perform enough realizations of the problem, the average time scale should be the same as that calculated from the expectation values above.

To numerically follow the evolution of the cloud, we need to specify the starting condition. Following previous authors, we adopt the following modified standard family of initial states

$$p=\frac{1}{1+\alpha _0}\mathrm{sech}^2\left(\frac{y}{1+\alpha _0}\right)=\frac{1}{1+\alpha _0}\left(1\mu ^2\right),$$
(40)
$$b=\left[\frac{\alpha _0}{(1+\alpha _0)K}\right]^{1/2}\mathrm{sech}\left(\frac{y}{1+\alpha _0}\right)=\left[\frac{\alpha _0}{(1+\alpha _0)K}\right]^{1/2}\left(1\mu ^2\right)^{1/2},$$
(41)
$$\mu =\mathrm{tanh}\left(\frac{y}{1+\alpha _0}\right),$$
(42)
where the ratio $`\alpha _0/K`$ represents the initial magnetic to thermal pressure ratio, i.e.,

$$\frac{\alpha _0}{K}=\frac{b^2}{p}=\frac{B^2/8\pi }{a^2\rho }.$$
(43)
Through the action of ambipolar diffusion, the fluid evolves toward a final state described by the above equations in the limit $`\alpha _00`$.

In addition to recovering the expectation values for the diffusion time, the numerical treatment also finds the deviations about the mean value (the expectation value). These deviations depend on how often we sample the distributions to evaluate the effective diffusion constant. If we sample the distributions often enough, say at every numerical time step, then the numerically calculated time scale converges to the expectation values found earlier. As outlined in §3.3, however, we expect the fluctuations to vary on the MHD crossing time scale (equation ). In dimensionless time units, we thus expect the diffusion constant to take on independent new values for time intervals longer than $`\tau _X0.1a/v_A0.02`$. For the sake of definiteness, we take $`\tau _X`$ = 0.02 for all of the simulations presented here. As we show below, however, the value of $`\tau _X`$ determines the width of the distribution of possible time scales; we can estimate this width analytically and then scale our results to other choices of $`\tau _X`$.

We compare our numerical results with the analytical formulation in Figure 2. We have used uniform distributions of the fluctuations with amplitude $`A`$ (see equations ) and set $`\alpha _0=10`$ so that the magnetic pressure is ten times larger than the thermal pressure at the start of the calculation. In these simulations, the diffusion constant $`D_j`$ = $`(1+\xi )^3/(1+\eta )^{3/2}`$ is changed every time interval $`\tau _X=0.02`$ by sampling the distributions of $`\xi `$ and $`\eta `$. The resulting time scale $`\tau _\mathrm{e}`$ for the magnetic field strength to decrease by one e-folding is plotted versus the fluctuation amplitude in Figure 2. The solid curve shows the expectation value for the time scale $`\tau _\mathrm{e}=\tau _{\mathrm{e0}}/\mathrm{\Lambda }`$, where $`\tau _{\mathrm{e0}}`$ is the time scale in the absence of fluctuations. The symbols show different realizations found by numerically integrating the problem. For this set of simulations, we have selected the amplitude $`A`$ randomly as well. As expected, the numerical integrations agree with the analytic predictions of §3.3. In particular, the mean value of the distribution of time scales from the numerical experiments closely follows the expectation value.

The width of the distribution (for a given amplitude $`A`$) increases with the amplitude of the fluctuations. As as result, fluctuations not only force ambipolar diffusion to take place more rapidly, on average, but they also allow the process to sample a range of time scales. This range of possible time behaviors, as illustrated in Figure 2, is determined (in part) by the time scale $`\tau _X`$ that specifies how often the effective diffusion constant takes on different values. In the limit $`\tau _X0`$, the ambipolar diffusion time scale approaches the expectation value and the width of the time scale distribution shrinks to zero. In the opposite limit, the effective diffusion constant takes on a single (but randomly selected) value for the whole evolutionary time. In dimensionless units, the e-folding time for ambipolar diffusion is $`\tau _\mathrm{e}`$ = 5 – 10 in this regime of parameter space. With $`\tau _X`$ = 0.02, near the low end of the range estimated from MHD considerations, the effective diffusion constant takes on $`N`$ = 250 – 500 different values during an evolutionary run.

We can quantify the width of the distribution of time scales for a given distribution of fluctuations and a given amplitude. For purposes of illustration, we use uniform distributions of fluctuations with amplitude $`A`$ = 0.886, which corresponds to the case in which the fluctuations are large enough to account for the observed non-thermal line widths in molecular clouds (see §3.5). Figure 3 shows the distribution of ambipolar diffusion time scales, calculated both numerically and analytically. The numerical results were obtained by running the aforementioned code 10,000 times to build up the distribution shown by the solid curve in Figure 3. For comparison, the dashed curve shows a gaussian distribution whose width is calculated as described below. The solid spike at $`\tau 12.5`$ represents the delta-function distribution that applies in the absence of fluctuations.

The distribution of time scales $`\tau `$ can be estimated analytically. The peak of the distribution (denoted here as $`\tau `$) is determined by the enhancement factor $`\mathrm{\Lambda }`$ calculated in §3.4. For the case shown in Figure 3, $`\mathrm{\Lambda }=2.25`$ and hence $`\tau `$ = 5.54. We can also find the width of the distribution. For a given realization of the problem, the time scale $`\tau `$ is determined by the effective value $`\stackrel{~}{D}`$ of the diffusion constant, i.e.,

$$\stackrel{~}{D}=D+\frac{D}{N_F}\underset{j=1}{\overset{N_F}{}}\left[\frac{D_j}{D}1\right],$$
(44)
where $`N_F`$ is the number of independent samplings of the diffusion constant (equation ), $`D`$ is the expectation value, and the $`D_j`$ are the particular choices taken during a given run. (Notice that in the limit $`N_F\mathrm{}`$, we must have $`\stackrel{~}{D}D`$.) For simplicity, we consider $`N_F=\tau /\tau _X`$ to have a fixed value, although the exact value also varies from case to case (this complication is a higher order effect). The case-to-case variation $`\mathrm{\Delta }D`$ = $`\stackrel{~}{D}D`$ in the diffusion constant is thus proportional to the sum of a large number (here, $`N_F300`$) of random variables $`\zeta _j`$ = $`D_j/D1`$. These variables $`\zeta _j`$ are constructed to have zero mean and follow a distribution given by the convolution of $`f(\xi )`$ and $`g(\eta )`$; we denote the variance of this $`\zeta _j`$ distribution as $`\stackrel{~}{\sigma }`$. Because of the central limit theorem, the distribution of the composite variable $`\mathrm{\Delta }D`$ takes a gaussian form. Similarly, distribution of time scales takes a gaussian form and its width is given by

$$\sigma =\stackrel{~}{\sigma }\sqrt{\tau _X\tau }.$$
(45)
We can compare this formula with our numerical results for the particular case shown in Figure 3 (uniform fluctuations with $`A`$ = 0.886). The dashed curve shows a gaussian distribution with this predicted width; the distribution is in good agreement with the histogram of results (the solid curve) from the numerical simulations. More precisely, we find that the expectation value for the ambipolar diffusion time scale $`\tau _\mathrm{e}=\tau _{\mathrm{e0}}/\mathrm{\Lambda }`$ $`5.54`$ is in good agreement with that found numerically, i.e., $`\tau `$ $``$ 5.62. The chosen distributions of fluctuations show that $`\stackrel{~}{\sigma }2.51`$ and hence equation  implies that the distribution of time scales should have width $`\sigma `$ 0.835. The numerical simulations imply almost the same value ($`\sigma `$ 0.834). For all quantities that can be compared, the numerical results and the analytic predictions agree to within about one percent.

As shown by equation , the width of the time scale distribution depends on the time scale $`\tau _X`$ over which the fluctuations occur, or, equivalently, the number $`N_F`$ of independent fluctuation samples. The value chosen for our numerical simulations, $`\tau _X=0.02`$, is near the low end of the expected range. For comparison, the dotted curve in Figure 3 also shows the distribution for a $`\tau _X=0.055`$, which corresponds to $`N_F=100`$. Notice that as the fluctuation time scale $`\tau _X`$ grows even larger (and $`N_F`$ decreases), the distribution of time scales grows wider and eventually departs from a gaussian form.

## 4 Discussion

In this paper, we have explored how fluctuations in the background fields affect the rate of ambipolar diffusion. These fluctuations force the magnetic field strength and the density to sample a distribution of values, rather than take on a single value at a given point in space and time. We have used a one-dimensional molecular cloud layer as a test problem to study the effects of such fluctuations.

The first principal result of this paper is that the time scale for ambipolar diffusion is altered by these fluctuations. In particular, fluctuations drive ambipolar diffusion to take place more rapidly, with the time scale shorter by a factor $`\mathrm{\Lambda }110`$. For typical conditions in molecular clouds cores, the enhancement factor is near the lower end of this range, $`\mathrm{\Lambda }23`$, but much larger enhancements remain possible. The ambipolar diffusion time scale depends on the distribution of fluctuations. For the case of uniform distributions, for example, the ambipolar diffusion time scale varies with the amplitude $`A`$ as shown in Figure 1. In general, the time scale also depends on the shape of the distributions.

For a given distribution of fluctuations, the ambipolar diffusion time scale also varies from realization to realization (see Figure 2). Descriptions of the ambipolar diffusion process thus face an interesting complication, which is our second principal result: The time scale for loss of magnetic support takes on a distribution of values instead of a single value. Consider two identical molecular cloud regions and suppose they are laced with fluctuations following a given distribution. Because the two regions experience incomplete (and different) samplings of the fluctuation distributions as their magnetic fields diffuse outwards, they will not exhibit the same diffusion time. This feature is more general than its manifestation in this particular test problem. When a physical system contains an effectively random element – in this context through chaos and turbulence – the outcomes must be described in terms of a probability distribution. For our test problem (ambipolar diffusion in a cloud layer), the single value of the e-folding time $`\tau _\mathrm{e}`$ is replaced by a distribution of values (see Figure 3). Furthermore, the distribution of possible time scales approaches a gaussian form; the most likely value for the time scale is shorter than the case without fluctuations by the enhancement factor $`\mathrm{\Lambda }=D/K`$ (see equations ) and the width of the distribution is given by equation .

This effect on the ambipolar diffusion time scale has important implications for star formation in molecular clouds. These clouds appear to be supported by magnetic fields and the observed magnetic field strengths are commensurate with this view. However, statistics of molecular cloud cores (with and without young stellar objects) argues that the (uncorrected) time scale for ambipolar diffusion may be too long to account for the observations. This work shows that magnetic fields can diffuse more rapidly than previous estimates suggest. This speed-up, along with any other enhancements (e.g., Ciolek & Basu 2001), can help account for the observed statistics of molecular cloud cores. Another complicating issue arises: Because the ambipolar diffusion time scale takes on a distribution of values, and this distribution can be rather wide if the fluctuations change on long time scales $`\tau _X`$, some core regions will experience much faster diffusion rates than others even if they have “the same” starting conditions. In this regime of diffusion activity, the cores that actually form stars are those which evolve on the “fast” side of the distribution, whereas the cores that happen to live on the “slow” side of the distribution will fail to form new stars.

This preliminary treatment of fluctuations, including their effects on ambipolar diffusion and star formation, remains incomplete in several respects. In this paper, we have separated the calculation of the diffusion process from the determination of the fluctuations. In particular, we have assumed a priori forms for the fluctuations to study their implications. In a complete treatment, one should calculate the fluctuations and their effects in a self-consistent manner. In addition, we have focused on long wavelength fluctuations and have not considered spatial gradients in the fluctuating part of the fields. Magnetic turbulence cascades down to small scales, however, so it is possible that fields fluctuate at length scales smaller than our MHD condition. This complication should also be considered in future work. Our present treatment is limited to one-dimensional slab models so that magnetic tension is not included; two-dimensional simulations should be done in the future. Another classical problem is the heating of molecular cloud regions by ambipolar diffusion; the effects of fluctuations on this mechanism should be considered. Finally, the act of star formation provides a source of new turbulence, which drives new fluctuations and can affect the ambipolar diffusion rates of neighboring cores; this feedback effect should also be studied. In any case, fluctuations in both the magnetic and density fields introduce an effectively random element into the ambipolar diffusion process, and thereby provide a rich class of new behavior for further study.

Acknowledgements

We would like to thank Charlie Doering, Phil Myers, Steve Shore, and Frank Shu for useful discussions. We also thank the referee – Glenn Ciolek – for comments that improved the paper. MF is supported by the Hauck Foundation through Xavier University. FCA is supported by NASA through a grant from the Origins of Solar Systems Program and by Univ. Michigan through the Michigan Center for Theoretical Physics.

APPENDIX A: Cylindrical Geometry

In this Appendix, we consider the effects of fluctuations on ambipolar diffusion in a molecular cloud filament. We consider only the simplest case of magnetic field lines that are aligned with the axis of the filament and depend only on the radial coordinate $`r`$, i.e., we have

$$𝐁=B(r)\widehat{z}.$$
$`(\mathrm{A1})`$
With this basic configuration, the equations of motion take the form

$$\frac{\rho }{t}+\frac{1}{r}\frac{}{r}(r\rho u)=0,$$
$`(\mathrm{A2})`$
$$\frac{u}{t}+u\frac{u}{r}=g\frac{1}{\rho }\frac{}{r}\left(P+\frac{B^2}{8\pi }\right),$$
$`(\mathrm{A3})`$
$$\frac{1}{r}\frac{}{r}(rg)=4\pi G\rho ,$$
$`(\mathrm{A4})`$
and finally

$$\frac{B}{t}+\frac{1}{r}\frac{}{r}(rBu)=\frac{1}{4\pi \gamma C}\frac{1}{r}\frac{}{r}\left\{\frac{r}{\rho ^{3/2}}B^2\frac{B}{r}\right\},$$
$`(\mathrm{A5})`$
where we have defined $`u`$ to be the radial (and only nonvanishing) component of the velocity, and we have made use of the relationship defined by equation . In practice, these filaments will be subject to clumping instabilities in both the linear (Gehman, Adams, & Watkins 1996) and nonlinear regimes (Adams, Fatuzzo, & Watkins 1994); in this derivation, however, we neglect this issue and focus on the effects of the cylindrical geometry.

Next, we introduce the fluctuations through the ansatz given by equation  and rewrite the problem in terms of a Lagrangian description of the dynamics (e.g., see §2). For the cylindrical geometry considered here, the relevant Lagrangian coordinate is the mass per unit length $`\sigma `$ along the filament within a radius $`r`$ of the central axis, i.e.,

$$\sigma _0^r\rho (r^{},t)r^{}𝑑r^{}.$$
$`(\mathrm{A6})`$
Notice that the variable $`\sigma `$ differs from the true mass per unit length by a factor of $`2\pi `$ which has been omitted for simplicity. The original problem in the variables $`(r,t)`$ is now transformed to one in new variables $`(\sigma ,t)`$ and the derivatives transform according to

$$\frac{}{t}+u\frac{}{r}\frac{}{t}|_\sigma ,$$
$`(\mathrm{A7})`$
$$\frac{}{r}r\rho \frac{}{\sigma }|_t,$$
$`(\mathrm{A8})`$
$$u\frac{r}{t}|_\sigma ,$$
$`(\mathrm{A9})`$
With this transformation, the equation of continuity becomes

$$\frac{r}{\sigma }=\frac{1}{r\rho }.$$
$`(\mathrm{A10})`$
For an isothermal equation of state, the force equation can be written in the form

$$\frac{1}{r}\frac{^2r}{t^2}=\frac{4\pi G}{r^2}_0^r\rho (1+\eta )r^{}𝑑r^{}+\frac{a^2}{1+\eta }\frac{}{\sigma }[\rho (1+\eta )]+\frac{1}{1+\eta }\frac{}{\sigma }\left[\frac{B^2(1+\xi )^2}{8\pi }\right],$$
$`(\mathrm{A11})`$
and the nonlinear diffusion equation for the magnetic field becomes

$$\frac{}{t}\left[\frac{B(1+\xi )}{\rho (1+\eta )}\right]=\frac{1}{1+\eta }\frac{}{\sigma }\left\{\frac{r^2B^2(1+\xi )^2}{4\pi \gamma C\rho ^{1/2}(1+\eta )^{3/2}}\frac{}{\sigma }[B(1+\xi )]\right\}.$$
$`(\mathrm{A12})`$
Equations \[A10 – A12\] exhibit exactly the same form as their slab counterparts (equations ). As such, the effects of fluctuations on ambipolar diffusion in a cylindrical filament will scale exactly as for the slab geometry considered in the main text.

APPENDIX B: Short Wavelength Fluctuations

We consider the effects of short wavelength fluctuations in this Appendix. As noted in the main body of the text, this analysis is complicated by the fact that the solutions to stochastic differential equations depend on the manner in which various limits are taken (Doering 1990). The formulation presented below thus represents one possible approach to the general problem.

The most likely source of short wavelength fluctuations is MHD turbulence, which is present in most regions of molecular clouds. The MHD condition, already built into the ambipolar diffusion equations, requires that the neutral fluid remain coupled to the ions and to the magnetic field. Physically, this condition is met if the ion-neutral collision frequency $`f_{in}=\gamma \rho _i`$ exceeds the frequency associated with the MHD turbulence. The latter frequency can be approximated by $`f_{mhd}v_A/\lambda `$, where $`v_A=B/(4\pi \rho )^{1/2}`$ is the Alfvén wave speed and $`\lambda `$ is the length scale of the fluctuations. As a result, the coupling condition requires that

$$\chi >0.09\frac{b}{p},$$
$`(\mathrm{B1})`$
where $`\chi `$ is the dimensionless turbulence length scale as defined by equation .

If we assume that the fluctuations are both spatially and temporally symmetric (which means that $`\xi `$ and $`\eta `$ are not correlated with their first order derivatives), then the following relations hold:

$$F(\xi ,\eta )\frac{\eta }{\mu }0\mathrm{and}F(\xi ,\eta )\frac{\xi }{\mu }0,$$
$`(\mathrm{B2})`$
and

$$F(\xi ,\eta )\frac{\eta }{\tau }0\mathrm{and}F(\xi ,\eta )\frac{\xi }{\tau }0,$$
$`(\mathrm{B3})`$
for all well-behaved functions $`F(\xi ,\eta )`$. Under these conditions, the quasi-equilibrium state described in §3.2 remains valid.

For short wavelength fluctuations, expanding the diffusion equation using the same approach as presented in §3.3 yields the form

$$(1+\xi )\frac{}{\tau }\left(\frac{b}{p}\right)+\frac{b}{p}\left[\dot{\xi }\dot{\eta }\frac{(1+\xi )}{(1+\eta )}\right]=\frac{(1+\xi )^3}{(1+\eta )^{3/2}}\frac{}{\mu }\left\{\frac{b^2}{p^{1/2}}\frac{b}{\mu }\right\}+\left\{\frac{}{\mu }\left[\frac{(1+\xi )^3}{(1+\eta )^{3/2}}\right]\right\}\left(\frac{b^2}{p^{1/2}}\frac{b}{\mu }\right)$$

$$+\frac{(1+\xi )^2}{(1+\eta )^{3/2}}\left\{\frac{}{\mu }(1+\xi )\right\}\left\{\frac{}{\mu }\left[\frac{b^3}{p^{1/2}}\right]\right\}+\frac{b^3}{p^{1/2}}\frac{}{\mu }\left\{\frac{(1+\xi )^2}{(1+\eta )^{3/2}}\frac{}{\mu }(1+\xi )\right\}.$$
$`(\mathrm{B4})`$
We note that the second and third terms on the right hand side vanish when they are time averaged because they take the form given by equation \[$`(\mathrm{B2})`$\]. With this simplification, the time-averaged diffusion equation reduces to the form

$$\frac{}{\tau }\left(\frac{b}{p}\right)=D\frac{}{\mu }\left(\frac{b^2}{p^{1/2}}\frac{b}{\mu }\right)+G\frac{b^3}{p^{1/2}},$$
$`(\mathrm{B5})`$
where $`D`$ is defined in equation  and

$$G=\frac{}{\mu }\left[\frac{(1+\xi )^2}{(1+\eta )^{3/2}}\frac{}{\mu }(1+\xi )\right].$$
$`(\mathrm{B6})`$
For the case in which the fluctuations $`\xi `$ and $`\eta `$ are uncorrelated, the parameter $`G`$ simplifies to the form

$$G=\frac{1}{3}\frac{1}{(1+\eta )^{3/2}}\frac{^2}{\mu ^2}(1+\xi )^3.$$
$`(\mathrm{B7})`$
Similarly, for the case of perfectly correlated fluctuations, $`G`$ simplifies to the form

$$G=\frac{2}{3}\frac{^2}{\mu ^2}(1+\xi )^{3/2}.$$
$`(\mathrm{B8})`$
The relative size of the two terms on the right hand side of equation \[B5\] ultimately determines the behavior of the magnetic field diffusion. Since the fluctuations $`\xi `$ and $`\eta `$ vary on a length scale $`\chi 1`$, whereas $`b`$ and $`p`$ vary on a lengthscale $`1+\alpha _0Kv_A^2/a^21`$, a simple scaling analysis naively suggests that the second term (with coefficient $`G`$) would dominate over the first (with coefficient $`D`$). Upon closer inspection, however, we see that the derivatives of the fluctuations tend to cancel out, so that the relative sizes of $`D`$ and $`G`$ depend on the form of the fluctuations. In any case, however, this treatment does not yield an expression that can be scaled to the previous solutions with no fluctuations.

FIGURE CAPTIONS

Figure 1. Enhancement factor as a function of fluctuation amplitude. The solid curve shows the factor $`\mathrm{\Lambda }`$ by which the ambipolar diffusion process is sped up by fluctuations, as a function of their amplitude $`A`$ for uniform distributions. The dashed curves shows the enhancement factor for fluctuations in the magnetic field only, whereas the dotted curve shows the result for density fluctuations. Notice that fluctuations in the magnetic field produce a greater effect for relatively small amplitudes, but the density fluctuations are more important for larger amplitudes.

Figure 2. Comparison of numerical and analytic results. The time scale $`\tau _\mathrm{e}`$ is the e-folding time for the magnetic field strength to decrease. The time $`\tau _\mathrm{e}`$ is plotted here as a function of fluctuation amplitude $`A`$ for uniform distributions of magnetic and density fluctuations. The solid curve shows the expectation value of the e-folding time scale as calculated analytically in §3.3. The symbols show results from numerical integrations of the same problem for different samplings of the distributions (see §3.5), where we have used $`\tau _X`$ = 0.02 as the time scale over which the diffusion constant changes (see text).

Figure 3. Distribution of ambipolar diffusion times for a cloud layer with uniform fluctuations of amplitude $`A`$ = 0.886; this level of fluctuations is consistent with the non-thermal line widths observed in star forming regions. The solid curve (histogram) shows the result of 10,000 numerical simulations with different realizations of the fluctuations. The dashed curve shows the analytic prediction for the time scale distribution – a gaussian with a peak value given by the expectation value and with a width predicted by application of the central limit theorem. The dotted curve depicts a wider gaussian distribution that applies for longer fluctuation time scales (here, $`N_F`$ = 100 or $`\tau _X0.055`$). In the absence of fluctuations, the cloud maintains a single value for its ambipolar diffusion time, as shown by the delta-function spike at $`\tau _\mathrm{e}12.5`$.