text:
# Simulation-Based Investigation of a Model for the Interaction Between Stellar Magnetospheres and Circumstellar Accretion Disks

## 1 Introduction

The formation of outflows from young stellar objects (YSO’s) has been a subject of interest for some time (see summaries by Königl & Pudritz, 2000; Reipurth & Bally, 2001). These outflows tend to consist of a highly collimated “jet” and a poorly collimated outflow, often referred to as the “disk wind” (Kwan & Tademaru, 1988, 1995). Jets have outflow velocities of $`50`$ to 500 km s-1 and may, in some cases, exhibit spectral characteristics of a hot plasma cooling down (Bacciotti, 1997). The disk wind component of the outflow typically has directionally averaged speeds of 0 – 50 km s-1 (Hirth et al., 1997; Solf & Böhm, 1993; Solf, 1997) and has spectral characteristics of a cooler plasma. Some systems contain a zero velocity component, which may be due to an oblique shock interaction between the disk wind and the disk.

As observational techniques and capabilities have improved, new insights have been gained concerning the outflow properties very close to the star. For example, jet collimation can occur very close to the source star-disk system. Observations of DG Tauri using Solf’s method (Solf, 1989) place the collimation region within 30 – 40 AU, which is the resolution of the observations of Solf & Böhm (1993). Recent measurements by Bacciotti et al. (2000) of the spatial distribution of the line emission for the H$`\alpha `$, \[N II\] 6583, \[O I\] 6363, and \[S II\] 6731 lines, at different distances from DG Tauri, imply an even closer collimation distance for this system of 14 AU. Also, Hubble Space Telescope observations of the HH 30 system (Burrows et al., 1996) imply that the jet must be fully collimated within 30 AU. For a general discussion of the collimation of YSO jets, see Eislöffel et al. (2000).

Observations of HH 30 also imply that the jet formation mechanism is episodic, producing outbursts (or “knots” along the jet axis) every $`2.5`$ years (Burrows et al., 1996). A probably related result is the observation of DG Tau by Solf (1997), who showed that the complicated \[O I\] 6300 line profile from the immediate environment of DG Tau changes drastically within 1 year. Larger scale observations of other systems such as HH 111 imply a knot ejection timescale of $`15`$ years (Reipurth et al., 1992), but there is no consensus as to whether or not the timescales of HH 30 and HH 111 are due to the same physical principle.

One result of so many recent observational insights is that significant progress has been made in understanding the mechanism which generates YSO outflows. While hydrodynamic means have been investigated (e.g., a stellar wind which is collimated via shock focusing; see, e.g., the early paper by Königl, 1982; Frank & Mellema, 1996), the models that have met the most success have relied on magnetic fields to launch and collimate the outflow. These models fall into two categories: those that rely solely on a magnetic field which threads the disk to centrifugally fling the wind (and eventually collimate into a jet; Pudritz & Ouyed, 1997; Camenzind, 1997; Kuwabara et al., 2000, see, e.g.,), and those that depend on the interaction between the stellar magnetic field and the accretion disk (e.g., Lovelace et al., 1995, 1999; Hayashi et al., 1996; Miller & Stone, 1997; Goodson et al., 1997, 1999; Goodson & Winglee, 1999). In the model of Goodson et al. (1997, henceforth GWB), Goodson et al. (1999, henceforth GBW), and Goodson & Winglee (1999, henceforth GW), a stellar dipole field couples the star to the disk, differential rotation between the star and the disk rapidly inflates the magnetic field via helicity injection, and an outflow burst is produced with both jet and disk wind properties. Magnetic reconnection allows the process to repeat. Hayashi et al. (1996) and GWB identified this basic mechanism within a few stellar radii, and GBW illustrated the resulting outflow out to several AU with a magnetohydrodynamic (MHD) simulation in a single case study. GW laid out some of the physics driving key characteristics of the outflow and star-disk interaction, allowing some predictions to be made beyond their case study.

The work presented in GBW and GW implied that the interaction between the stellar magnetosphere and the innermost regions of the accretion disk could produce an outflow that is qualitatively similar to those observed. The outflow consists of a knotty jet (collimated within 1 AU) that is morphologically similar to the observations of HH 30 (Burrows et al., 1996). However, the size scale is an order of magnitude smaller than HH 30, and the timescale for knot formation is an order of magnitude shorter. In addition to the apparent mismatch between observation and simulation results, issues that must be addressed are:

1. The dipole field from the star initially threads the disk everywhere. Thus, it is not clear if the outflow is a result of the interaction between the star and the disk (as claimed in GW), or if it is due to the magnetic field which threads the disk everywhere at $`t`$ = 0.
2. The results of GBW represent a single case study. The sensitivity of the jet formation mechanism to key attributes of the star-disk system was not directly established.

To address these issues, we repeated the basic simulation of GBW for various disk densities, stellar magnetic field strengths, and magnetic field configurations. Section 2 contains the details of the parametric cases and of the outflows that result in each case. The parametric simulations also allow us to evaluate the disk oscillation model of GW in §3. A semi-analytical approach, including a formulation of the extraction of angular momentum from the disk, leads to a modification and enables testing of the modified model. A summary and discussion follows, in §4.

## 2 Simulations and Outflow Results

We used the 2.5-D (axisymmetric) MHD simulation code employed by GWB, GBW, and GW to examine, parametrically, the effect of varying the accretion disk density, the stellar dipole field strength, and the initial magnetic field configuration. The simulation code employs a two-step Lax-Wendroff (finite difference) scheme (Richtmyer & Morton, 1967) to simultaneously solve the following time-dependent, MHD equations:

$`{\displaystyle \frac{\rho }{t}}`$ $`=`$ $`\mathbf{}(\rho 𝒗)`$ (1)
$`\rho {\displaystyle \frac{(𝒗)}{t}}`$ $`=`$ $`\rho (𝒗\mathbf{})𝒗\mathbf{}P{\displaystyle \frac{GM_{}\rho }{(𝒓^2+𝒛^2)}}\widehat{𝑹}+{\displaystyle \frac{1}{\mathrm{c}}}(𝑱\times 𝑩)`$ (2)
$`{\displaystyle \frac{e}{t}}`$ $`=`$ $`\mathbf{}[𝒗(e+P)]\left[{\displaystyle \frac{GM_{}\rho }{(𝒓^2+𝒛^2)}}\widehat{𝑹}\right]𝒗+𝑱𝑬`$ (3)
$`{\displaystyle \frac{𝑩}{t}}`$ $`=`$ $`\mathrm{c}(\mathbf{}\times 𝑬)`$ (4)
and uses

$`𝑬`$ $`=`$ $`{\displaystyle \frac{1}{\mathrm{c}}}(𝒗\times 𝑩)+\eta 𝑱`$ (5)
$`𝑱`$ $`=`$ $`{\displaystyle \frac{\mathrm{c}}{4\pi }}(\mathbf{}\times 𝑩)`$ (6)
$`e`$ $`=`$ $`{\displaystyle \frac{1}{2}}\rho |𝒗|^2+{\displaystyle \frac{P}{\gamma 1}}`$ (7)
where $`\rho `$ is the density, $`𝒗`$ the velocity, $`P`$ the scalar gas pressure, $`e`$ the internal energy density, $`𝑩`$ the magnetic field, $`𝑱`$ the volume current, $`𝑬`$ the electric field, c the speed of light, and $`\gamma `$ the ratio of specific heats (we used $`\gamma `$ = 5/3 in our simulations). In the momentum and energy equations (2 and 3), we include a source term for gravity in which $`G`$ is the gravitational constant, $`M_{}`$ is the stellar mass, $`𝒓`$ and $`𝒛`$ are the usual cylindrical coordinates, and $`𝑹`$ is the spherical distance from the center of the star ($`|𝑹|^2=|𝒓|^2+|𝒛|^2`$). The last term in equation 5 allows us to explicitly model the Ohmic diffusion of the magnetic field by specifying the value of the resistivity, $`\eta `$. As discussed by GW, the numerical diffusivity intrinsic to the simulations is measured to be $`5\times 10^{16}`$ cm2 s-1. In order to properly capture the physics of magnetic diffusion, we use $`\eta =10^{17}`$ cm2 s-1 in all of our simulations.

We use a “nested box” scheme, and for the work presented here, all simulations used 5 nested boxes, each with half the resolution (and so twice the physical size) as the one interior to it. Each box consists of a 401 $`\times `$ 100 (in the cylindrical $`r`$ and $`z`$ direction, respectively) grid. As in GBW and GW, the grid spacing is 0.1 $`R_{\mathrm{}}`$ in the innermost box, but the outermost box dimensions are 3.2 AU by 0.8 AU. The innermost (in $`𝑹`$) boundary consists of a 15 gridpoint sphere (the “surface of the star”), where $`\rho `$, $`P`$, $`𝑩`$, and $`𝒗`$ are held fixed at the initial values. We enforce reflection symmetry across the equatorial plane. Along the rotation axis, we require that $`B_\varphi =B_r=v_\varphi =v_r=0`$, and all other quantities are equal to the value at $`r=dr`$ (where $`dr`$ is the grid spacing). We adopt outflow conditions (the spatial derivative across the boundary is zero for all quantities) at the outermost boundary of the outermost box.

The initial conditions for the baseline case of this study is parametrically identical to the case illustrated by GBW. They describe the initial conditions for this case in considerable detail (see their table 1), which includes a disk similar to a classic $`\alpha `$ model accretion disk (Shakura & Sunyaev, 1973) with $`\dot{M}_{\mathrm{acc}}`$ = 10-7 $`M_{\mathrm{}}`$ yr-1 and $`\alpha `$ = 0.2 ($`r_{\mathrm{au}}`$)-1.5 (where $`r_{\mathrm{au}}`$ is the usual cylindrical coordinate in units of AU). The stellar magnetic field is dipolar, with a field strength of 1800 G at the pole of the 1 $`M_{}`$, 1.5 $`R_{}`$ central star. Initially, the stellar corona is in hydrostatic equilibrium, and the density falls of as $`𝑹^4`$. The star rotates with a period of 1.8 days (which is fairly close to the shortest observed rotation period of CTTS; see, e.g., Bouvier et al., 1986), and the disk is initially truncated at 8.5 $`R_{}`$ and has a midplane density of $`2.2\times 10^{10}`$ g cm-3 at 10 $`R_{}`$. The dipole field threads the disk everywhere at the start of the simulation.

To save computational time, this work focuses more on the outflow generation mechanism and the integrated outflow properties instead of the larger scale outflow morphology. The simulation cases were all carried out for a period of at least 160 days of real time (requiring 30-60 days of computing time per case).

Table 1 serves as a quick reference for all simulation cases in this study. The three cases which examine variations of initial disk density are the “low density,” “baseline,” and “high density” cases. The three which examine variations of stellar field strength are the “weak field,” “baseline,” and “strong field” cases. The relative factors by which the disk density and the magnetic field strength were varied were chosen based on energy balance. This is because the importance of the magnetic field is measured by the ratio of the magnetic energy to the kinetic energy in the disk. Since the disk kinetic energy is proportional to the disk density, and the magnetic energy scales as the square of the stellar field strength, we varied the initial disk density by factors of 4 and the initial magnetic field strength by factors of 2.

We also investigated simulations with a “primordial,” vertical magnetic field in the disk. So that the stellar magnetosphere was truncated (at the region of disk and stellar field equality) near the star, the simulations use stellar dipole strength equal to the weak field case, plus a vertical field of 0.07 Gauss everywhere in the simulation region. These parameters give a stellar magnetosphere that is truncated at 29 $`R_{}`$ in the equatorial plane. The stellar magnetosphere has an open or closed topology, depending on the polarity of the primordial field. The vertical field in these simulations cases is relatively weak (in that the stellar dipole is dominant for an extended region; the ratio of magnetic to kinetic energy near the disk inner edge is $`10^5`$) but is sufficient to assess the effect of the dipole field that initially threads the disk at all radii in the other cases. In order to initially exclude the disk from the stellar magnetosphere, we truncated the disk at 29 $`R_{}`$ for these models. Therefore, our models with various magnetic field topologies comprise three cases: One with a downward-directed ($`z`$) vertical field (“closed magnetosphere”), one without a vertical field (“truncated disk”), and one with an upward-directed field (“open magnetosphere”). All three have a stellar dipole field with the same strength as the weak field case, and all have a disk truncation radius at 29 $`R_{}`$ (which is 3.4 times further out than all the other simulations).

### 2.1 Outflow Morphology

Figures 13 illustrate the general outflow morphology for all cases. In these figures, the density grayscale (logarithmic) and poloidal velocity vectors represent data in the outermost box for a single “snapshot” in time. Data from the inner 100 $`\times `$ 100 grid points of the outer box are shown, reflected about the rotation axis to better show the jet. The grayscale emphasizes material in the outflow, but does not capture the relatively high density distribution near the star and disk inner edge. For a detailed look at the complex and time-dependent density distribution near the star, see figures 3 and 4 of GBW. Note that different times are represented for each case. This is because the launching mechanism is essentially periodic, with the period being different for each case. The innermost regions of the disk oscillate radially, and each oscillation triggers a large scale reconnection event which alters the shape of the outflow in a periodic way. For the purpose of comparison, we show the various cases at a similar phase (but not at the same time) of their outflow.

As evident in figure 1, disk density variations of a factor of 16 have no major effect on the outflow morphology. The basic outflow, in all cases, consists of a wide angle component (the “disk wind”) and a jet. Also note that even at this close-in scale (on the order of an AU), the general characteristics of both the jets and disk winds are similar – knots in the jet are apparent in each case (most clearly in the high-density case), as are the “wispy” structures between the disk and jet (associated with magnetic reconnection; see GBW).

While the general morphologies are remarkably similar, some differences between the outflows in figure 1 are apparent. The low density case has a higher outflow speed, a lower outflow density, and a larger spacing between knots in the jet that the other cases. Conversely, the high density disk case has a higher density, slower outflow, and the spacing between the knots in the jet is reduced considerably. This reduction in knot formation period is related to a reduction in the period of radial disk oscillations (discussed further in §3).

Figure 2 represents the simulation cases which examined variations in stellar magnetic field strength. Again, the morphologies of all runs are similar, indicating that variations of a factor of 4 in magnetic field strength produce no significant variation in large scale morphology. Each outflow produces a knotty jet with a disk wind and a wispy structure in the region between the disk and jet. However, differences between the outflows in figure 2 are evident. The outflow is slower in the weak field case and faster in the strong field case. Also, the structure of the disk wind in the weak field case has a more complicated appearance than the other two cases.

Figure 3 illustrates the outflows for the cases with different large scale magnetic field topologies. Insets in each of the three panels show the topology in the initial conditions. Once again, note that the general morphologies of the outflow are all similar. Recalling from table 1 that these three topologies are variants of the weak field case, it is interesting to note that all of these cases exhibit the same general characteristics of the weak field case (fig. 2). All produce outflows with jets and disk winds (including wispy structures) with velocities similar to the weak field case, and all three cases in figure 3 produce outflows with a complicated structure (relative to the other cases) on the one AU scale. This figure is particularly noteworthy because, in the bottom two panels, the disk is threaded everywhere by a vertical field, so it has much more magnetic energy in the disk (at large radii) than the cases with only a stellar dipole. The fact that the outflows in these cases are not significantly different from the top panel, supports the underlying premise of the model: It is the interaction between the stellar magnetosphere and the surrounding accretion disk that produces the outflow, and the effect of the dipole field lines initially threading the disk at large radii is insignificant.

The primary conclusion to be drawn from figures 1 \- 3 is that the outflow formation mechanism is robust. Also, to the extent these parameters were varied, the jet morphology produced by this mechanism is independent of the density of the accretion disk, the strength of the stellar magnetic field, and the large scale magnetic field topology. As seen in these figures, however, the magnitude and timescales of outflow characteristics do depend on the magnetic field strength and the density of the accretion disk, and the remainder of this paper addresses these dependencies in a quantitative way.

### 2.2 Sensitivity of Outflow Properties

Figure 4 shows the variation in radial outflow velocities for all cases. For the data in this figure we first calculated the radial velocity ($`𝒗\widehat{𝑹}`$) on a sphere with a radius of 0.67 AU (90 grid points in the outer box) centered on the star. Due to the periodic nature of the system, and because we were most interested in the global outflow properties, we used a time average value. Since it takes a while for the outflow to clear out the entire simulation region, our average only includes data from the interval of 50 – 160 days. We thus obtain the average velocity as a function of angle from the pole. The data in figure 4 represent the maximum of the time-averaged velocity, which occurred at $`40`$ degrees for all cases. The leftmost panel shows the velocities as a function of the disk variation, the middle as a function of the initial magnetic field energy, and the right panel as a function of the vertical field strength. Note that the entire range of velocities (100 to 200 km s-1) is well within the range of observed outflow velocities (50 to 500 km s-1) and corresponds to the velocities in most of the observed HH objects (see Raga et al., 1996).

In panel a), the radial outflow velocity decreases with increasing disk density, while in b), the outflow velocity increases with increasing magnetic field energy. Finally, panel c) shows that the outflow velocity is roughly independent of the large-scale magnetic field topology. It is important to note that the velocities in panel c) are similar to the weak field case. Remember that the simulations represented in c) have their disk initially truncated 3.3 times further out than the weak field case, and two of them have an additional vertical “primordial” field. We see from panel b) that the outflow velocity depends on the magnetic energy (from the stellar dipole), so the fact that panel c) shows essentially no trend suggests that it is only the magnetic energy in the stellar magnetosphere that is important for accelerating the outflow.

The data in panels a) and b) can be fit with a power law ($`v_\mathrm{w}\rho _\mathrm{d}^n`$ and $`v_\mathrm{w}E_{\mathrm{mag}}^n`$, respectively, where $`\rho _\mathrm{d}`$ is the initial disk density and $`E_{\mathrm{mag}}`$ is the initial magnetic energy associated with the stellar dipole field). A least squares minimization fit gives $`n`$ = $`0.2`$ and 0.3 for panel a) and b), respectively. We note that there is not necessarily a physical justification for fitting all of the data with a power law, but we will use it to quantify the dependence of various system properties on the parameters we varied.

Figure 5 shows the variation in mass outflow rates, $`\dot{M}_\mathrm{w}`$, for the three parametric studies. To obtain these points (and the ones for $`L_\mathrm{w}`$ and $`\dot{J}_\mathrm{w}`$, discussed below), we first calculated the mass flux (also energy flux and angular momentum flux) through a sphere at 0.67 AU. The mass outflow rate (also $`L_\mathrm{w}`$ and $`\dot{J}_\mathrm{w}`$) is given by an integration over the surface of the sphere and a time average over the interval of 50 – 160 days. In an effort to exclude the disk (which can be inflowing or outflowing, but with very high mass relative to the jet and disk wind), the data shown represent an integration over the surface up to an opening angle of 75 degrees. To get an idea of the relative uncertainties (due to the fact that we are representing a time-dependent property by a time-averaged value) in the data, we calculated the cumulative average of $`\dot{M}_\mathrm{w}`$ during the 50 – 160 day interval. We take three times the standard deviation of this cumulative average over the interval of 85 – 160 days as our uncertainty. This quantity reflects the uncertainty due to the periodicity of $`\dot{M}_\mathrm{w}`$ and due to the settling down of the system to a fixed period (discussed further in §3).

All of the data in figure 5 have $`\dot{M}_\mathrm{w}10^7`$$`10^8`$ $`M_{\mathrm{}}`$ yr-1, consistent with observations of YSO outflows (Hartmann, 1997). Also, $`\dot{M}_\mathrm{w}`$ increases with both disk density and magnetic field energy. Once again, a power law nicely quantifies this dependence with $`n`$ = 0.2 for both panel a) and b) (within uncertainty). The variation of mass outflow rate with the different field topologies, panel c), shows no clear trend.

Figure 6 shows the variation in outflow luminosity, $`L_\mathrm{w}`$, for the parametric cases studied. The energy flux is the sum of the kinetic, magnetic, and the thermal energy fluxes.

$$\mathrm{\Phi }_{\mathrm{eng}}=\left(\frac{1}{2}\rho v^2+\frac{P}{\gamma 1}+\frac{B^2}{8\pi }\right)(𝒗\widehat{𝑹}),$$
(8)
and the data in figure 6 represent the time-averaged, integrated flux (as for $`\dot{M}_\mathrm{w}`$ above). The format of the figure is the same as figures 4 and 5. The disk density has no obvious influence on $`L_\mathrm{w}`$, and the best fit power law is consistent with $`n`$ = 0. The outflow luminosity does, however, depend on the initial magnetic field energy, suggesting a power law with $`n`$ = 0.4. Finally, the addition of a weak vertical field has no significant impact on outflow luminosity.

An analysis of the angular momentum carried in the outflow, $`\dot{J}_\mathrm{w}`$, is deferred to §3.1, where we use a semi-analytical formulation to understand the controlling parameters of the angular momentum extracted from the inner edge of the accretion disk.

### 2.3 Large Scale Magnetic Field Evolution

In the jet launching process described in GBW (also in Hayashi et al., 1996, and, to some degree, in Lovelace et al., 1999), the large scale magnetic field undergoes rapid inflation. This inflation results from the differential rotation between the star and the disk which twists the magnetic field, injecting helicity, and increasing the overall magnetic energy. The resulting magnetic field has characteristics that differ significantly from those of a dipole. Figure 7 illustrates the initial magnetic field expansion of the baseline case, demonstrating the rapid change from a dipole to a more complicated and open geometry.

The basic process of magnetospheric inflation can be summarized as follows (see GWB and GBW for details): The star is initially connected to the disk via the dipole magnetic field. Differential rotation between the star and disk induces an azimuthal component of the magnetic field, increasing the magnetic pressure ($`|𝑩|^2/8\pi `$) in the twisted region, which in turn inflates the poloidal component of the magnetic field.

The expanded magnetic field has many properties that are different from the simple dipole. Figure 8 illustrates the most significant differences for the discussion at hand. The figure shows the total magnetic field strength, $`|𝑩|`$, on a 45 degree line in the meridional plane for the baseline case. The lower solid line shows the original dipole field. The top solid line shows the total field, $`|𝑩|`$(tot) = $`(B_r^2+B_\varphi ^2+B_z^2)^{1/2}`$, and the middle solid line shows the poloidal field, $`|𝑩|`$(pol) = $`(B_r^2+B_z^2)^{1/2}`$, both at $`t`$ = 61.6 days.

Note that the dipole field decreases with distance as $`R^3`$, the evolved total field decreases as $`R^n`$ with $`n`$ near 1.7, and the evolved poloidal field decreases with $`n`$ roughly equal to 2.0. Also note that the evolved field has a “dipole-like” ($`n3`$) region close to the star ($`R5R_{}`$), with a clear transition to an “expansion-dominated” ($`n`$ 1.5-2) region beyond. The total and poloidal field strengths are almost identical in the dipole-like region, indicating that $`B_\varphi `$ is negligible there. A piece-wise continuous fit was made of the expanded field (along a $`45^{}`$ slice), of the form:

$$|𝑩|=\mathrm{Max}[B_{01}\left(\frac{R_{}}{R}\right)^{n_1},B_{02}\left(\frac{R_{}}{R}\right)^{n_2}]$$
(9)
In this fit, $`n_1`$ is the power law that applies in the interior, dipole-like region, and $`n_2`$ is the power law that applies in the exterior, expanded region. $`B_{01}`$ and $`B_{02}`$ are the field strengths that one would obtain at the surface of the star, if the power law could be extrapolated all the way to the stellar surface (this is roughly the true surface value for the dipole-like fit, but not for the expanded field fit). Both the total and poloidal field are fit in this way. The dashed lines in figure 8 are the fits to the total and poloidal field at $`t=61.6`$ days.

Table 2 lists the values of the magnetic field fit parameters for all of the simulations. Since the values are time-dependent (periodic), and it takes some time for the flow to be established in the entire simulation grid, table 2 data have been averaged during the interval of 50 – 160 days. To get an idea of the uncertainty of these values (due to the episodic nature of the system), we first obtain the cumulative average during this interval. The uncertainty listed in the table is the standard deviation of this cumulative average during the interval of 85 – 160 days. For each case, column 2 and 3 contain the $`B_{01}`$ (in Gauss) and $`n_1`$ values. Since the total and poloidal field is nearly identical in this interior fit, we list only the value for the poloidal field. The last 4 columns contain similar values for the exterior portion of the fit for both the total and poloidal field components.

When $`n_2`$(pol) is near 2.0, the poloidal field is essentially open ($`n=2`$ for a purely radial field). This is roughly true for all cases, though the ones with weaker dipole fields (row 4 and 6 – 8) have $`n_2`$(pol) significantly less than than 2. The primary conclusion from this data is that the total magnetic field strength will behave as $`R^n`$ where $`n`$ is between 1.5 and 2.0. The cases with additional vertical field in the disk (last two rows of table 2) have the largest magnetic to kinetic energy ratios (far from the disk inner edge) of all cases. However, since the magnetic field in the outflow originates in the stellar magnetosphere, the field in the evolved outflow of these cases is identical to the weak field and truncated disk case.

### 2.4 Summary of Parametric Results

The most striking result of this parametric study is that every case examined produces an outflow similar to that of GBW. For a variation of disk density of a factor of 16, a variation in stellar magnetic field of a factor of 4, and variations in initial disk truncation points and large scale magnetic field topologies, all of the simulations produce outflows with jets and disk winds. The only differences are in the details of the star-disk interaction and in the magnitude of the outflow properties.

Table 3 summarizes the sensitivities of the outflow properties and morphology to increasing the disk density and increasing the stellar magnetic energy. The listed numerical values are the spectral index ($`n`$) from power laws that represent various data (e.g., the dotted lines in figures 4 \- 6). The least squares power law fits were weighted by the uncertainty (as described in §2.2), and the errors in table 3 are the formal 1-sigma errors of the fit. Note that some of outflow properties listed in table 3 are discussed in the following section.

The main result of table 3 is that variation of key system parameters do indeed cause variations in outflow properties, but the relationship is weak (i.e., $`|n|`$ is small). The star-disk-field system is self-regulating, due to the interaction between the disk and the outflow. The system adjusts to mitigate the impact of any single parameter change. Also, the large scale field topology does not significantly affect any of the outflow properties. This result implies that it is indeed the interaction between the stellar magnetosphere and the accretion disk that drives the outflow, as opposed to the existence (via artificial initial conditions as in GBW) of a dipole field that threads the disk everywhere at $`t`$ = 0.

## 3 Evaluation of the Disk Oscillation Model of GW

As noted in GBW and GW, the inner edge of the disk undergoes radial oscillations. These oscillations are actually a flux exchange process – on each oscillation, the stellar magnetosphere diffuses into the innermost regions of the disk and differential rotation inflates the magnetic field to an open configuration. The inflated field topology spins the disk down via magnetocentrifugal launching of disk plasma (see, e.g., Blandford & Payne, 1982). As it spins down and moves inward toward the star, the disk inner edge is stripped away from the rest of the disk. At the point of “closest approach”, a reconnection event alters the magnetic field topology, and the spun down portion of the disk flows along the newly reconnected field lines, accreting onto the star.

Figure 9 illustrates the disk oscillations for the simulation cases where the accretion disk density is varied. The solid line is the radial location of the inner edge of the disk as a function of time. The dashed line shows the location of the Keplerian co-rotation point (where the stellar solid body rotation rate is equal to the Keplerian rotation rate). The figure clearly demonstrates the dependence of disk oscillations on the disk density. As disk density is decreased, the amplitude and period of the oscillations increase. Note that it takes time for the system to settle down to a constant period (the physical explanation for this is discussed in §3.2). In most cases, the system becomes quasi-periodic after about 50 days, but the low density case seems only to be settling down near the end of the 160 day simulation. Much of our uncertainty in all quantities is due to this effect, and we have calculated our uncertainties to reflect this.

As the average disk location approaches the Keplerian co-rotation point, the difference between the disk and stellar magnetosphere rotation rates is reduced, so the magnetosphere will expand more slowly. This is true in the high density case (bottom panel), but here the disk oscillations are actually faster than in the other cases. Also, even in the high density case, the oscillation period is much longer than the inflation timescale (which is roughly the time for the disk inner edge and star to rotate to a differential angle of $`\pi `$; see Uzdensky, 2002; GW). This supports the basic tenet of GW that the oscillation period is driven by the spin-down rate of the inner edge of the disk (via magnetic torques) rather than the inflation rate of the magnetic field.

Ideally, the parametric simulation results presented here would allow for the disk oscillation model of GW to be tested. Unfortunately, the data in §2.3 require us to relax their assumption that the expanded magnetic field has reached a power law configuration characterized by $`|𝑩|=B_{}(R_{}/R)^2`$. While this may be true for the poloidal field of many cases (see table 2), it is more generally true that the expanded field reaches a configuration with $`|𝑩|=B_0(R_{}/R)^n`$, with $`n`$ between 1.5 and 2.0 (and $`B_0`$ not necessarily equal to the stellar surface value). We will also discard the assumption of GW that the accretion disk has reached a secular steady state, where the accretion rate onto the star (regulated by the interaction between the stellar magnetosphere and the inner edge of the disk) equals the accretion rate given by the $`\alpha `$-disk description of the disk. While this secular steady state does exist near the end of simulations, the accretion rate in the outer portion of the simulated disk is unlikely to be equal to the $`\alpha `$-disk model. This is because the physics of the Shakura-Sunyaev $`\alpha `$-disk model are not fully captured in the simulations. However, we will show in the following section that the extraction of angular momentum from the inner edge of the disk can be simply understood and determined from simulation data. In light of this, we present a modified disk oscillation model (in §3.2), in which the angular momentum extraction can be a prescribed value – thus removing the assumption of secular stead state in the disk.

We should point out that it is our ultimate goal to be able to determine the properties of the outflow ($`\dot{M}_\mathrm{w}`$, $`v_\mathrm{w}`$, $`L_\mathrm{w}`$, oscillation period, etc.) generated by the star-disk interaction a priori from the physical parameters of the system ($`R_{}`$, $`B_{}`$, disk accretion rate, etc.). Unfortunately, the complexity of the model, which includes a relatively large number of parameters, makes this difficult. The following work in this section is thus meant to be an early step in that direction. As such, this semi-analytical work relies on some quantities derived from the outcome of the simulations (e.g., $`B_0`$). The work presented in §3.1 and §3.2, therefore, is useful for testing our understanding of the most relevant physics controlling the oscillations, and thus the accretion/ejection mechanism. It will also serve as a check for self-consistency in the simulations.

### 3.1 Extraction of Angular Momentum from the Disk

Consider an inner region of the disk that is threaded by a magnetic field. From the analytical work of Weber & Davis (1967), we note that the angular momentum transported in the wind of a magnetic rotator is $`\dot{J}_\mathrm{w}=\dot{M}_\mathrm{w}\mathrm{\Omega }_{\mathrm{fp}}R_\mathrm{A}^2`$ (where $`\mathrm{\Omega }_{\mathrm{fp}}`$ is the angular velocity of the footpoint of the flux tube, and $`R_\mathrm{A}`$ is the Alfvèn radius). Assuming nearly rigid rotation of plasma out to $`R_\mathrm{A}`$ (for a discussion of this assumption, see, e.g., Michel, 1969), plasma will be magnetocentrifugally launched along field lines via the Blandford-Payne mechanism (Blandford & Payne, 1982), provided the field lines make an angle of more than 30 degrees from vertical. This angle is typically 45 to 55 degrees in our simulations. For a simple approach, we’ll assume that the velocity of the wind at the Alfvèn radius is $`v_{R_\mathrm{A}}=\mathrm{\Omega }_{\mathrm{fp}}R_\mathrm{A}`$. This is strictly true for $`R_\mathrm{A}R_\mathrm{d}`$ (where $`R_\mathrm{d}`$ is the position of the field footpoint in the disk). The Alfvèn radius, then, can be estimated by examining the energy balance at that location:

$$\frac{1}{2}\rho _{R_\mathrm{A}}v_{R_\mathrm{A}}^2=\frac{B_{R_\mathrm{A}}^2}{8\pi },$$
(10)
When expressed in terms of outflow properties, and assuming the magnetic field behaves as $`|𝑩|=B_0(R_{}/R)^n`$, this becomes

$$\frac{1}{2}\frac{\dot{M}_\mathrm{w}\mathrm{\Omega }_{\mathrm{fp}}}{4\pi fR_\mathrm{A}}=\frac{B_0^2}{8\pi }\left(\frac{R_{}}{R_\mathrm{A}}\right)^{2n}$$
(11)
Here, we have assumed that $`\dot{M}_\mathrm{w}=4\pi R_\mathrm{A}^2\rho _{R_\mathrm{A}}v_{R_\mathrm{A}}`$. Note that this assumption is not valid for an anisotropic wind, so we introduced a filling factor, $`f`$, by which the mass outflow rate is reduced. Reducing equation 11 and solving for $`R_\mathrm{A}`$ gives

$$R_\mathrm{A}=\left(\frac{fB_0^2R_{}^{2n}}{\dot{M}_\mathrm{w}\mathrm{\Omega }_{\mathrm{fp}}}\right)^{\frac{1}{2n1}}$$
(12)
and so

$$\dot{J}_\mathrm{w}=(\dot{M}_\mathrm{w}\mathrm{\Omega }_{\mathrm{fp}})^{1\frac{2}{2n1}}(fB_0^2R_{}^{2n})^{\frac{2}{2n1}}$$
(13)
The dependence of $`\dot{J}_\mathrm{w}`$ on various outflow properties is interesting. The strength of the dependence relies only on $`n`$, the magnetic field power law fall off. Note that if $`n`$ = 1.5, equation 13 becomes

$$\dot{J}_\mathrm{w}=fB_0^2R_{}^3,$$
(14)
and $`\dot{J}_\mathrm{w}`$ is independent of both $`\dot{M}_\mathrm{w}`$ and $`\mathrm{\Omega }_{\mathrm{fp}}`$. In general, for any $`n`$ between 1.5 and 2.0, the exponent on the factor $`(\dot{M}_\mathrm{w}\mathrm{\Omega }_{\mathrm{fp}})`$ in equation 13 is substantially smaller than on the factor $`(fB_0^2R_{}^{2n})`$. Therefore, the extraction of angular momentum should be most sensitive to the magnetic field strength.

The simulations verify this conclusion. Figure 10 shows the angular momentum outflow rate for all cases. The data (squares) are calculated as described in §2.2, and shown in the same format as figure 4. It is clear that the angular momentum carried in the simulation outflow is more dependent on magnetic energy than on disk density, which influences $`\dot{M}_\mathrm{w}`$ and $`\mathrm{\Omega }_{\mathrm{fp}}`$ (see figures 5 and 9). In fact, panel a) of figure 10 shows no obvious dependence of $`\dot{J}_\mathrm{w}`$ on disk density at all (also see the fifth row of table 3). This also appears to be true for the cases with various field topologies shown in panel c), where there is no clear trend, especially considering the larger uncertainties.

The simulations also offer a quantitative test of equation 13. We used the time-averaged simulation values for $`\dot{M}_\mathrm{w}`$, $`\mathrm{\Omega }_{\mathrm{fp}}`$, $`B_0`$, and $`n`$ (the magnetic field power law fall off) for each case and calculated $`\dot{J}_\mathrm{w}`$ using equation 13. Since the inner edge of the disk spends most of the time beyond 5 $`R_{}`$ (see figure 9), the value of $`B_0`$ and $`n`$ to use comes from the outer portion of the piece-wise fit to the magnetic field strength in the simulations. Also, since angular momentum is added to the wind via the poloidal field anchored in the disk, the calculated values of $`\dot{J}_\mathrm{w}`$ use $`B_{02}`$(pol) and $`n_2`$(pol). We chose a single value of the filling factor ($`f`$ = 0.43), so that the calculated values match the baseline result. These calculated values (triangles in figure 10) can be directly compared to the simulated values (squares), and the agreement between them is remarkable. The largest disagreement is (for the closed magnetosphere case) at the 50% level.

In the context of equation 13 and the assumptions at the beginning of this section, it seems clear that variations in the location of the inner edge of the disk (which is the launching point for much of the outflow) and in $`\dot{M}_\mathrm{w}`$ (which is related to the conditions at the footpoint of the flux tube) have little relative impact on the angular momentum ejected from the system. The primary reason for this is the restructuring of the initial magnetic field. As evident in table 2 (also see figs 7 and 8), the magnetic field rapidly evolves into a state that follows a radial power law with $`n`$ between 1.5 and 2.0, increasing the importance of $`B_0`$ relative to $`\dot{M}_\mathrm{w}`$ and $`\mathrm{\Omega }_{\mathrm{fp}}`$.

In addition to this, the magnetic field in the outflow seems to be restructured in such a way that it “softens” the effect (on $`\dot{J}_\mathrm{w}`$) of the parameters varied. For example, the $`B_{02}`$(tot) values for the weak field, baseline, and strong field cases are 550, 750, and 1180 Gauss, respectively. The entire range in these values is only a factor of 2, even though the entire range of stellar surface dipole field strengths (which is the only simulation parameter changed among these cases) is a factor of 4. Thus, a given variation in the stellar surface magnetic field strength cannot be used in conjunction with equation 13 alone to determine the resulting variation in $`\dot{J}_\mathrm{w}`$. The self-regulating nature of the system indeed makes prediction difficult.

### 3.2 Modified Disk Oscillation Model

In their oscillation model, GW assumed that the magnetic field associated with the stellar magnetosphere would diffuse into the disk radially at a given radial velocity ($`v_{\mathrm{diff}}`$). The differing azimuthal velocities of the inner edge of the disk and the stellar magnetosphere (which rotates at the angular rotation rate of the star) leads to a rapid inflation of the magnetosphere (note that $`v_\varphi v_r`$ and $`v_{\mathrm{diff}}`$). The inflated magnetic configuration removes angular momentum from the inner edge of the disk (via magnetocentrifugal launching of disk plasma), spinning it down. The inner edge of the disk will eventually achieve the condition where its inward radial velocity exceeds the velocity at which the magnetic field is diffusing out into the disk ($`|v_r||v_{\mathrm{diff}}|`$). At this point, the inner edge of the disk separates from the rest of the disk, and spirals inward to the star, where a fraction of it is ultimately accreted by the star, and a fraction is expelled in the outflow. Thus, the diffusive velocity ($`v_{\mathrm{diff}}`$) determines how much of the disk separates for each oscillation, and so sets the mass accretion rate.

There are two main assumptions in the disk oscillation model of GW that must be modified. They calculated the radial velocity of the disk inner edge (as it spins down; their eq. 9) as a function of $`r`$, $`v_{\mathrm{diff}}`$, $`R_\mathrm{d}`$ (the radial position of the disk inner edge), $`\mathrm{\Sigma }(R_\mathrm{d})`$ (the surface density there), and $`B_{}`$. We first modify their assumption that $`|𝑩|=B_{}(R_{}/R)^2`$ to the more general case, $`|𝑩|=B_0(R_{}/R)^n`$, and get

$$v_r=B_0^2R_{}^{2n}r^{22n}\frac{\{GM_{}[2\pi R_\mathrm{d}\mathrm{\Sigma }(R_\mathrm{d})]^2rB_0^2R_{}^{2n}\pi R_\mathrm{d}\mathrm{\Sigma }(R_\mathrm{d})r^{42n}\}^{\frac{1}{2}}}{GM_{}[2\pi R_\mathrm{d}\mathrm{\Sigma }(R_\mathrm{d})]^2B_0^2R_{}^{2n}\pi R_\mathrm{d}\mathrm{\Sigma }(R_\mathrm{d})(42n)r^{32n}}$$
(15)
Note that when $`n`$ = 2, equation 15 reduces to GW’s equation 9.

The second assumption of GW to discard is that the accretion rate from the disk inner edge (regulated by the interaction of star’s magnetosphere and the disk) is the same as the accretion rate given analytically by the $`\alpha `$-disk description. Instead, we assume that the angular momentum carried in the outflow (see §3.1 and fig. 10) is equal to the amount extracted from the disk inner edge.

Figures 11 and 12 contain the simulated disk oscillation amplitudes and periods (squares and dotted lines; in the same format as fig. 4). Since it takes some time for the oscillations to settle down (see fig. 9), the simulated amplitude of the disk oscillations plotted in figure 11 are the maximum separation of the star and disk inner edge 160 days into the simulation. The simulated periods (fig. 12) are also taken at this time. As a measure of the uncertainty (represented by error bars in the figures), we use the difference between the amplitude at 160 days and the value at 85 days. For the period, we use the same method, but the error bars have been doubled to show up better on the plot.

The “settling down” of the disk oscillations is a result of the initial conditions. Equation 15 indicates that, as $`R_\mathrm{d}`$ increases, $`v_r`$ decreases (due to a decrease in the magnetic torque and an increase of the amount of angular momentum at that radius), and a ring from the disk inner edge has further to fall, so the period increases. Also, GW showed (their fig. 7) that the accretion rate decreases for increasing $`R_\mathrm{d}`$. If the initial state of the system (in the simulations) is such that the accretion rate regulated by the star-disk interaction ($`\dot{M}_{\mathrm{inner}}`$) is greater than the accretion rate further out in the disk (at $`r>R_\mathrm{d}`$; $`\dot{M}_{\mathrm{outer}}`$), the inner region of the disk will be cleared out (thus increasing $`R_\mathrm{d}`$) until the two accretion rates are equal.

Notice that the cases in panels a) of figures 11 and 12 show a relatively strong variation in amplitude and period, even though the angular momentum extraction rate is roughly equal among them (see fig. 10). This is explained by the fact that, for a denser disk, $`\dot{M}_{\mathrm{outer}}`$ is larger, so the final state of the system will have a smaller amplitude and, therefore, a shorter period.

The previous section showed that angular momentum extraction depends most strongly on the magnetic field strength. This seems to imply that, for stronger magnetic fields, the oscillation periods should be shorter (since a ring of disk material would spiral in faster). This would indeed be the case if $`R_\mathrm{d}`$ were a fixed quantity. However, when the magnetosphere-disk interaction can extract more angular momentum from a given radius, $`\dot{M}_{\mathrm{inner}}`$ increases, so the final state of the system will actually have a larger $`R_\mathrm{d}`$ and a longer period. This is evident in panels b) of figures 1012.

To calculate the period and amplitude of disk oscillations, following GW, we used an iterative process in which we first pick a radial location for the disk inner edge ($`R_\mathrm{d}`$) and integrate equation 15 to obtain a spin down time. Equation 15 requires the case-specific values of $`B_{02}`$(tot) and $`n_2`$(tot) listed in table 2. This method gives a simultaneous prediction of the oscillation amplitude and period. As opposed to GW, we adopt the $`R_\mathrm{d}`$ value that satisfies the requirement that the angular momentum extracted from the ring of material (as it falls toward the star) divided by the oscillation period (taken as twice the spin down time) equals the angular momentum extraction rate in the simulations (displayed in fig. 10). As in GW, we need an estimate of $`v_{\mathrm{diff}}`$, in order to calculate the mass of the ring of material that is stripped off of the disk inner edge. We use $`v_{\mathrm{diff}}\eta /\beta h`$ ($`h`$ is the initial disk scale height and depends on $`R_\mathrm{d}`$), which differs from GW’s approximation by the factor $`\beta `$. This factor is necessary because the average disk scale height during the simulations is larger than the initial value. This has two causes: a) the magnetic interaction increases the scale height above the initial analytical value, and b) the numerical diffusion of disk material is not negligible ($`h`$ is typically only a few grid points at $`R_\mathrm{d}`$). For each case, then, we used our explit value of $`\eta =10^{17}`$ cm2 s-1, and chose $`\beta `$ to best fit the disk oscillation period. In all cases, $`\beta `$ ranged from 3 to 9.

The predictions of the modified model presented here (triangles in fig.’s 11 and 12) match the simulated results (squares) well for the cases with various stellar magnetic field strengths and disk densities, but are consistently too small for the cases with various field topologies. Remember that the the cases in panel c) of figures 11 and 12 have a disk which is initially truncated further out than all other cases. Since their initial truncation radius is further out than the final disk oscillation amplitude of the weak field case (to which they are otherwise similar), these cases have an initial state in which $`\dot{M}_{\mathrm{inner}}<\dot{M}_{\mathrm{outer}}`$. The mismatch between the data an predictions for these models may be due to a significant alteration of the disk density profile from its initial state.

The agreement of the predictions of this modified model with the data implies that the fundamental assumptions of the model are valid (at least in the simulations). Those assumptions are: 1) The disk oscillations are associated with a flux interchange process where magnetic flux diffuses into the inner-most edge of the disk. 2) The disk oscillations are driven by spin-down torques associated with the magnetic field that threads the disk. 3) The spin-down rate is driven by the surface density of the accretion disk and the magnetic field topology. 4) What appear to be large scale disk oscillation are actually the spin-down of a ring of material from the inner most edge of the disk.

## 4 Summary and Discussion

The mechanism presented by GBW and GW for the launching of jets is promising, but it is difficult to make a priori predictions, given the large number of parameters in the star-disk system. One has to consider the effects of the mass of the star and disk, the magnetic field strength, and the stellar rotation rate, to name a few. In addition, the disk responds to the conditions in its environment, so one must also consider the self-consistent interaction between the disk, outflow, and changing magnetic field topology. The complexity of this model prompted us to carry out a parameter survey using magnetohydrodynamic simulations. Using the case presented by GBW as a baseline, we examined cases with various initial disk densities, various stellar dipole field strengths, and cases with additional “primordial” field associated with the disk. The results of this survey are:

1. The fundamental mechanism for jet formation is robust.
2. The star-disk-outflow system is self-regulating.
3. The mechanism is independent of magnetic flux that may exist in the disk.
4. We have improved the disk oscillation model of GW.

The robust nature of the mechanism is best demonstrated by examination of figures 13. For variations of initial disk density, variations in stellar magnetic field strength, and variations in initial disk truncation points and large scale magnetic field topologies, all of the simulations produced outflows with similar morphologies. All outflows consisted of a highly collimated, knotty jet and a less collimated “disk wind” (including a “wispy” structure associated with magnetic reconnection). Also, the outflows were always episodic with the period set by the oscillations of the inner edge of the disk. The magnitude and timescales of outflow characteristics, however, do depend on the magnetic field strength and the density of the accretion disk, and we were able to address these dependencies in a quantitative way.

As seen in figures 46, 1012, and in table 3, variation of key system parameters do indeed cause variations in outflow properties, but the relationship is fairly weak (the strongest power law relating a system property to the initial conditions has $`n=0.44`$). This self-regulating nature of the system suggests that the self-consistent interaction between the disk, outflow, and stellar magnetosphere “softens” the effect of any single parameter change. For example, we found in §3.1 that the large scale magnetic field structure rapidly evolves to a configuration that removes angular momentum from the disk at a rate that depends most strongly on the evolved field and much more weakly on other parameters. Also, when the stellar magnetic field is changed, the outer region of the evolved field (which is important for the extraction of angular momentum) has characteristics which only partially reflect the surface value change.

The robust and self-regulating nature of this system imply that, if such a mechanism does occur in astrophysical systems, it should indeed be ubiquitous, and may actually be the prime mechanism for the formation of non-relativistic astrophysical jets. Unfortunately the self-regulating nature also means that observable (or deduceable) properties of outflows, such as the mass outflow rate or the rate of angular momentum expelled from the system, will contain a limited amount of information about the properties of the star-disk systems.

Figures 46 and 1012 also demonstrated with some certainty that the system properties examined do not depend on the amount or direction of magnetic flux that initially threads the disk. For each system property we examined, there is at least a weak dependence on the initial stellar magnetic field energy, evident in panel b) of each of these figures. However, as seen in the c) panels, for the cases with additional vertical field (and therefore additional magnetic energy in the disk), there were no convincing trends. Implicitly, therefore, it is the interaction between the stellar magnetosphere and the inner edge of the accretion disk that drives the outflow, as opposed to the existence (via artificial initial conditions) of a dipole field that threads the disk everywhere at $`t`$ = 0. Note, however, that for stronger vertical fields, magnetocentrifugal launching of disk plasma may increase at larger radii (leading to an increased collimation of the material launched from the disk inner edge). Magnetized disk models (e.g., Pudritz & Ouyed, 1997; Camenzind, 1997; Kuwabara et al., 2000) will become important, independent from and in addition to the mechanism presented here. A system with a dynamically important stellar dipole field and disk-associated field might produce an outflow with an episodic/unsteady core (closest to the axis of rotation) and a “cocoon” of material originating from the disk.

We also evaluated and modified the disk oscillations model of GW (§3). Using a semi-analytical method, we tested and confirmed many of their basic assumptions, including the idea that the oscillation period is set by spin down time of disk inner edge, not by the inflation of magnetosphere. Therefore, we demonstrated that the basic physical mechanism for disk oscillations is understood.

There is still work to be done with the accretion/ejection mechanism presented here in order to decide whether such a process occurs in real astrophysical systems. Future work is planned to determine the importance of other properties of the system, though the qualitative mechanism is valid for a wide range of parameters. In particular, the stellar rotation rate certainly influences the details of the outflow characteristics, but the inflation of the magnetosphere will always occur, even in the extreme case of a non-rotating star, due to the differential rotation in the accretion disk. Also, a more realistic treatment of the accretion disk, possibly including the effects of the disk viscosity or the magnetorotational instability (Balbus & Hawley, 1991), would affect the star-disk interaction through the accretion rate ($`\dot{M}_{\mathrm{outer}}`$), and will thus change the period and amplitude of the disk oscillations. Future work will also examine the system under true three-dimensional conditions. The effect of radiative cooling/heating may not be negligible, especially considering that the coupling of the field to the gas is dependent on the specific ionization state and composition of the circumstellar matter.

The disk oscillation process, and therefore the jet launching mechanism, is fundamentally driven by a diffusive process which is poorly understood. Recall that the mechanism requires magnetic reconnection (a diffusive process with physics that are partially understood; see Priest & Forbes, 2000) for accretion and re-expansion. GBW point out that this model is not particular sensitive to the physics of magnetic reconnection, since regions of oppositely directed magnetic field are driven together by the disk oscillations. However, this is only a partial solution to the dilemma, since the disk oscillation mechanism itself is initiated by the diffusion of stellar magnetospheric flux into the disk inner edge.

The timescale for knot formation predicted by this mechanism is a factor of $`15`$ different than observed for HH 30, and the rudimentary scaling laws offered by this parametric study do not offer an obvious solution to the discrepancy. The timescale predicted by this work (several tens of days) cannot presently be resolved by direct imaging (there would be several tens of knots in a typical optical resolution element), so one cannot yet rule it out. However, even if there exists a short timescale present in astrophysical systems that can be explained by this mechanism, there is then no detailed theoretical prediction of the observed longer timescale. A primary focus of future work with this model will be the timescale for disk oscillations/knot formation.

Another criticism of this work (as well as many other models; see §1) addresses the assumption that the central star has a strong dipole magnetic field (dominating over other multipolar components which may exist). Detailed observations of the magnetic fields in T Tauri stars (see, e.g., Safier, 1998; Johns-Krull et al., 1999) suggest that strong, ordered fields may not dominate; rather, the stellar field topology appears to be time-variable and non-axisymmetric. Whether such a field can participate in the mechanism described here (e.g., via non-axisymmetric magnetic flux tubes connecting the star and disk) is not known.

We thank Frank Shu, Chris Johns-Krull, and Bo Reipurth for useful discussion about the model and observed outflow properties. We also thank the anonymous referee for helpful suggestions that improved the paper. This research was supported by NSF grant AST 97-29096.


format: text

text:
# Solid state effects on exciton states and optical properties of PPV

## Abstract

We perform ab initio calculations of optical properties for a typical semiconductor conjugated polymer, poly-para-phenylenevinylene, in both isolated chain and crystalline packing. In order to obtain results for excitonic energies and real-space wavefunctions we explicitly include electron-hole interaction within the density-matrix formalism. We find that the details of crystalline arrangement crucially affect the optical properties, leading to a richer exciton structure and opening non-radiative decay channels. This has implications for the optical activity and optoelectronic applications of polymer films.

Ordered films of organic conjugated polymers are of strategic relevance for novel optoelectronic devices frie+99nat . In addition, they offer an ideal scenario for the study of electronic and excitonic confinement, being composed of quasi-one-dimensional (1D) systems, arranged however in a three-dimensional (3D) crystalline environment. Long linear chains can indeed be thought of as 1D systems with the highest degree of 1D confinement for a polyatomic system ($``$ 2-5 Å). Since the main optoelectronic characteristics derive from the mobile $`\pi `$-electrons, delocalized along the chain backbone, and non bonding to neighbouring chains, the vast majority of studies adopt the single-chain model: in fact, complete ab initio theoretical studies of the optical properties have been performed for isolated chains, highlighting the strong confinement expected for such systems rohl-loui99prl ; ruin+01sm . Recently, very simplified models of polymers in a “crystalline medium” have been approached through ab initio theory hors+99prl ; however the effect of crystal structure or side chains is completely neglected. This is also usual practice in the intepretation of experimental data: e.g. data on poly-para-phenylenevinylene (PPV) and its alkylated and metoxilated derivatives (MH-PPV, MEH-PPV, etc…) are usually bundled together as representative of PPV, based on the supposition that details of the 3D structure are not relevant. The picture emerging from this analysis is far from clear, however: quoted exciton binding energies differ by an order of magnitude leng+94prl ; chan+97prb ; alva98-01all ; camp+96prl ; bart-bass97prl ; mose+01sm\_a , there is an on-going controversy about the existence of charge-transfer excitons or excimers conw96sm ; wu-conw97prbrc ; yan+95prl ; camp+96prl ; frol+00prl , and about conditions for efficient light emission. The question then arises if, on the contrary, solid state effects cannot be negleted and crystallization entails also structure-specific 3D interchain coupling effects: recent data on oligothiophene crystals seem to suggest that this might be the case mucc+00prb ; buss+02tobe .

Here we address these issues through a full ab initio calculation of optical spectra and real-space exciton wavefunctions for PPV, in both the isolated chain and the crystalline phase. In the single PPV chain, we find that the lowest singlet exciton extends over a few monomers along the chain and is optically active, with a binding energy of about 0.6$`÷`$0.7 eV. In the crystal we find evidence of significant interchain coupling, with decrease of the first optically active exciton binding energy to about 0.2$`÷`$0.3 eV alva98-01all ; camp+96prl ; bart-bass97prl ; as a consequence of both interchain interaction and detailed crystalline structure, a dark split-off singlet exciton appears at lower energy, which introduces a new path for non-radiative emission and electron-hole dissociation, and is expected to decrease drastically the photoluminescence efficiency yan+95prl ; furthermore, we find a much richer exciton structure, with charge-transfer excitons close enough to the first direct exciton camp+96prl  to account for early electron-hole dissociation and apparently-vanishing binding energies mose+01sm\_a ; fitting all the pieces together, our results can explain the controversies on the optical behavior of these crystals.

The exciton binding energy, $`E_b`$, is a central quantity in the photophysics of polymer semiconductors because it is related to the spatial extension of the excitation and to the probability of radiative emission/absorption and electric-field-induced dissociation of the exciton. It is defined as the energy required to separate a bound electron-hole pair, the exciton, into a free electron and a free hole conw96sm ; alva98-01all ; camp+96prl ; bart-bass97prl ; hill+00cpl . For PPV, the experimental values that are currently reported range from 0.2$`÷`$0.4 eV alva98-01all ; camp+96prl ; bart-bass97prl  to around 1 eV leng+94prl ; chan+97prb , but energies of less than 0.1 eV are also proposed mose+01sm\_a . These values correspond to spatial extensions varying from a single monomer (Frenkel-like excitons) up to many unit cells (Wannier-like excitons).

From the theoretical point of view, calculating the properties of excitons in extended polymeric systems is much more difficult than in conventional inorganic-semiconductor quantum wires ross-moli96all , because it requires to combine an accurate treatment of the Coulomb interaction with a realistic microscopic description of the single-particle electronic structure beyond effective mass approximations. Standard quantum-chemistry schemes belj+99jcp  for treating excitations in finite molecular systems are intrinsically not suitable for crystals. To this end, we have implemented a first-principles density-matrix approach where the correlated electron-hole spectrum is obtained directly from the interband polarizations $`d_\nu ^{}c_\mu ^{}`$, resulting from the creation of an electron in conduction state $`\mu `$, with energy $`ϵ_\mu ^e`$ and a hole in valence state $`\nu `$ with energy $`ϵ_\nu ^h`$ ($`c_\mu ^{}`$ and $`d_\nu ^{}`$ are the Fermionic creation operators; the labels $`\mu =(c,𝐤)`$ and $`\nu =(v,𝐤)`$ stand for both the band index and the wavevector). The polarization eigenstates are obtained by direct diagonalization of the two-body Schrödinger equation:

$$(\overline{ϵ}_\mu ^e\overline{ϵ}_\nu ^h)A_{\mu \nu }+e^2\underset{\mu ^{}\nu ^{}}{}(2\delta _SV_{\mu \nu ^{},\nu \mu ^{}}W_{\mu \mu ^{},\nu \nu ^{}})A_{\mu ^{}\nu ^{}}=E_xA_{\mu \nu }$$

Here the first term on the right-hand side accounts for the uncorrelated electron and hole quasi-particle energies, the second for electron-hole screened direct ($`W`$) and unscreened exchange ($`V`$) Coulomb interaction ($`\delta _S=1`$ for singlets, 0 otherwise). The interaction kernels are written in terms of the single-particle wavefunctions $`\psi `$:

$$W_{\mu \mu ^{},\nu \nu ^{}}=\psi _\mu ^e(𝐫_1)\psi _\mu ^{}^e(𝐫_1)w(𝐫_1,𝐫_2)\psi _\nu ^h(𝐫_2)\psi _\nu ^{}^h(𝐫_2)𝑑𝐫_1𝑑𝐫_2$$

$$V_{\mu \nu ^{},\nu \mu ^{}}=\psi _\mu ^e(𝐫_1)\psi _\nu ^{}^h(𝐫_1)v(𝐫_1,𝐫_2)\psi _\nu ^h(𝐫_2)\psi _\mu ^{}^e(𝐫_2)𝑑𝐫_1𝑑𝐫_2$$

where $`v(𝐫_1,𝐫_2)=|𝐫_1𝐫_2|^1`$ is the bare Coulomb potential while $`w(𝐫_1,𝐫_2)=\epsilon ^1(𝐫_1,𝐫)v(𝐫,𝐫_2)𝑑𝐫`$ is the screened one. The excitonic eigenenergies $`E_x`$ and eigenstates $`𝑨`$ are provided by the solution of the above eigenvalue problem; hence, this description hohe01prb  closely resembles the earlier one obtained within a Green’s function scheme, which comprises the solution of the Bethe-Salpeter equation sham-rice66pr ; hank-sham79prl .

The first step is therefore the calculation of non-interacting electron and hole quasiparticle energies $`\overline{ϵ}`$ and wavefunctions $`\psi `$ that we obtain through a density-functional theory hamiltonian in the local density approximation kohn-sham65pr ; pwscf  (with norm-conserving pseudopotentials and a plane-wave basis with 50 Ry energy cutoff), and by subsequently adding a rigid shift to the conduction band energies to account for the self-energy correction (taken from Ref. rohl-loui99prl, ). Moreover, in order to make the full calculation computationally feasible for complex systems as polymer crystals, we approximate the full effective screening function by a constant diagonal dielectric tensor $`\underset{¯}{\underset{¯}{\epsilon }}`$, whose components are calculated ab initio within linear response pwscf . For the case of the isolated polymer chain, the assumption of space homogeneity is indeed very drastic, so we introduce an effective screening tensor, which is defined within the volume $`\mathrm{\Omega }_{\text{eff}}`$ where the interaction actually takes place, while $`\underset{¯}{\underset{¯}{\epsilon }}=1`$ outside $`\mathrm{\Omega }_{\text{eff}}`$. The components of $`\underset{¯}{\underset{¯}{\epsilon }}`$ are again calculated ab initio within linear response and $`\mathrm{\Omega }_{\text{eff}}`$ is determined selfconsistently with $`\underset{¯}{\underset{¯}{\epsilon }}`$ to encompass the region of significant electron-hole correlation function.

In order to better understand the role of interchain coupling for the crystal case, we first consider the isolated chain as reference system. We show in Fig. 1(a) the results for the absorption spectrum of the isolated chain rohl-loui99prl ; ruin+01sm , which is highly polarized along the chain ($`z`$-direction); for $`x`$ and $`y`$ polarizations absorption close to band gap is negligible. We find an optically active $`{}_{}{}^{1}B_{u}^{}`$ singlet exciton with energy approximately 0.6$`÷`$0.7 eV below the single-particle band-gap energy. This lowest singlet state extends over a few unit cells along the chain rohl-loui99prl ; ruin+01sm , the large binding energy arising essentially from the strong lateral confinement. To gain insight in the nature of this state, in Fig. 2 we plot its electron and hole contributions: $`\rho ^h(𝐫_h)=_\mu |_\nu A_{\mu \nu }\psi _\nu ^h(𝐫_h)|^2`$ and $`\rho ^e(𝐫_e)=_\nu |_\mu A_{\mu \nu }\psi _\mu ^e(𝐫_e)|^2`$. Both have their maxima on the vinylene chains but in different bonding configurations, and thus vinylene atoms are expected to undergo the largest photoinduced lattice relaxation; this suggests the assignment of the phonon replica observed in PPV photoluminescence to the vinylene C-C stretch mode (around 0.2 eV tian+91jcp ; saka+92jpc ; mart+99prb ; capa-cald02tobe ) and explains the difference in lineshape eckh+89jcp  with poly-para-phenylene (PPP), where the vinylene chains are not present.

We point out that the exciton signature is accompanied by suppression of the typical 1D van Hove singularity ($`1/\sqrt{E}`$) ruin+01sm . This is a phenomenon resulting from dimensionality, also found in inorganic semiconductor quantum wires ross-moli96all ; ogaw-taka91prb . It is relevant for the experimental determination of the exciton binding energy: the energy separation between the exciton and the onset of the 1D absorption continuum cannot be extracted directly from a one-photon optical spectrum, since the oscillator strength associated to the latter is vanishingly small ross-moli96all ; a separate measurement is required to extract the quasiparticle gap leng+94prl ; alva98-01all ; camp+96prl ; rina+01prb .

We now turn to the PPV crystal, that we calculate in a realistic herring-bone structure chen+92pol ; gome-conw93prb , by keeping an orthorombic cell with $`a=8.07`$ Å, $`b=5.08`$ Å, $`c=6.54`$ Å and the setting angle $`\mathrm{\Phi }=52^{}`$ (see inset in Fig. 3). The influence of crystalline structure can be felt already on the single-particle bands: due to the specific crystal symmetry, with two chains per unit cell, all bands are doubled gome-conw93prb ; the splitting between doublets gives rise to the splitting between excitonic states (the so-called Davydov splitting davy62 ), and is strongly dependent on interchain interaction. This feature alone points to a clear difference between PPV and other substituted PPV derivatives, such as MEH-PPV, that cannot pack in the herring-bone structure, have only one chain per cell, and thus would not show Davydov splittings. Note that this very important feature cannot be seen either from structureless crystal models hors+99prl  or simple molecular aggregate models without translational symmetry corn+01sm . Other general features are that dispersion along the $`\mathrm{\Gamma }`$-Z direction (along the chain) is of the order of $``$ 2.2 eV, as for the isolated chain, and dispersion perpendicular to the chain is vanishingly small.

We show in Fig. 1(b) the correlated spectrum for the crystal. It is again characterized by a strong singlet exciton peak (marked DA in the figure) and quenching of the oscillator strength at the onset of the single-particle continuum. While the single-particle band gap is reduced in the crystal gome-conw93prb , the optical gap in the correlated spectrum is close to the single-chain result. Hence, the binding energy of this exciton is reduced to about 0.2 eV, consistent with the increased screening in the crystal hors+99prl . The first absorption peak is still strongly polarized along z, confirming the strong optical anisotropy of the material.

Now we turn our attention to the other remarkable features of the optical behavior: (a) the single-particle band-splitting results in an exciton Davydov splitting of $``$ 120 meV, with the corresponding introduction of an optically inactive (dark) state below the first active exciton (DD arrow in the figure); (b) other excitonic states are also introduced above the DA peak, with no oscillator strength along z, but just barely active on the perpendicular polarization (CT1 and CT2 in the figure).

The corresponding exciton wavefunctions for DA and CT1 are shown in Fig. 3, where we plot the probability density distributions $`|\mathrm{\Psi }^{exc}(𝐫_e,𝐫_h)|^2`$. Note that the exciton wavefunction $`\mathrm{\Psi }^{exc}(𝐫_e,𝐫_h)`$, expressing the probability amplitude of finding the electron in $`𝐫_e`$ and the hole in $`𝐫_h`$, is computed within our formalism directly in terms of the interband polarization $`𝑨`$ as $`\mathrm{\Psi }^{exc}(𝐫_e,𝐫_h)=_{\mu \nu }A_{\mu \nu }\psi _\mu ^e(𝐫_e)\psi _\nu ^h(𝐫_h)`$. The DA state is remarkably confined on a single chain, only slightly more extended in $`z`$-direction than in the isolated chain. (The DD state has the same distribution as the DA, by symmetry: they are both direct as electrons and holes are on the same chain.) The CT1 state, on the other hand, shows very small on-chain contribution, with high probability of finding electron and hole on different chains, and can be characterized as an excimer; same comments apply to the CT2 state.

In summary, we have found that the lowest correlated states of PPV crystal do involve more than one chain and are very dependent on crystal structure buss+02tobe ; however, they retain a quasi-one-dimensional character, and the consequent suppression of 1D singularity precludes the use of one-photon experiments to extract exciton binding energies. Furthermore, we have found that, due to interchain coupling and symmetry effects, the lowest electronic excitation is a dark exciton, that is optically inactive but provides a non-radiative decay path. This can strongly quench the photoluminescence with respect to non-interacting chains. In view of the impact of crystallization on the optical properties, we studied also the influence of the packing density: without changing the symmetry, we increased the interchain distance in the $`\pi `$-stack direction $`b`$ (see inset in Fig. 3); we find that an increase of $``$ 20 % leads to an increase in the binding energy (to $`0.4`$ eV) and a drastic decrease in the Davydov splitting by an order of magnitude, thus restoring radiative efficiency. These results suggest that the spread in photoluminescence (PL) and photoconductivity (PC) data observed by different groups may arise from differences in the aggregation state of the polymer. Moreover, since real films have regions with different chain packings, the dominant contributions to PL and PC can originate from different regions of the same sample (the least and the most closely packed, respectively). Another solid state effect is the appearance of bound charge-transfer excitations close above the first exciton peak: this not only can enhance photoconductivity, but also can be the origin of the infra-red active modes seen on two-photon experiments mose+01sm\_a ; mose+00cpl .

Our results show that solid state effects are indeed very important for the final optical structure of a polymer film. This conclusion could only be reached because of the possibility of calculating first-principles optical spectra of realistic three-dimensional polymer structures.

We are grateful to Fausto Rossi for many relevant suggestions in the early stages of this work. We also thank Ulrich Hohenester, Lucia Reining, Eric K. Chang, and Guido Goldoni for useful discussion. This work was supported in part by MURST (Cofin-99) - Italy, by FAPESP - Brazil, and by a CNR-CNPq bilateral project.


text:
# Global versus Local Ferromagnetism in a Model for Diluted Magnetic Semiconductors Studied with Monte Carlo Techniques

## Abstract

A model recently introduced for diluted magnetic semiconductors by Berciu and Bhatt (PRL 87, 107203 (2001)) is studied with a Monte Carlo technique, and the results are compared to Hartree-Fock calculations. For doping rates close to the experimentally observed metal-insulator transition, a picture dominated by ferromagnetic droplets formed below a $`T^{}`$ scale emerges. The moments of these droplets align as the temperature is lowered below a critical value $`T_C`$$`<`$$`T^{}`$. Our Monte Carlo investigations provide critical temperatures considerably smaller than Hartree-Fock predictions. Disorder does not seem to enhance ferromagnetism substantially. The inhomogeneous droplet state should be strongly susceptible to changes in doping and external fields.

PACS numbers: 75.50.Pp, 71.30.+h, 75.40.Mg

Recent advances in experimental techniques have allowed for the introduction of magnetic elements into semiconducting hosts beyond the solubility limit. This has stimulated the exciting research area of diluted magnetic semiconductors (DMS) with the prospect of manipulating the charge as well as the spin degrees of freedom, and its possible technological applications in such novel fields as “spintronics”. Prominent among these new compounds is the III-V system Ga1-xMnxAs, where ferromagnetic (FM) transition temperatures as high as 110 K were obtained. Although some of the basic ingredients of the physics of these materials such as the local antiferromagnetic (AF) exchange between the Mn spins and the charge carriers are known, the origin of such fairly high critical temperatures is still in question. The ferromagnetism itself can be understood as being carrier-induced, in a similar fashion as the FM state in the double-exchange (DE) model for manganites at intermediate doping. However, alternative descriptions based on Ruderman-Kittel-Kasuya-Yosida (RKKY) interactions between the impurity spins are claimed to lead to transition temperatures comparable to those observed. Yet it is not obvious whether RKKY can apply for systems with small Fermi energies as is the case for compensated semiconductors. Additionally, the effects of disorder that accompanies chemical doping was neglected in this approach. Recent annealing studies have shown the relevance of disorder and defects in DMS. It was proposed , based on Hartree-Fock (HF) calculations, that disorder in the dopant position and hopping amplitudes has a strong effect on the magnetic properties of these compounds, and should be incorporated in trying to understand the behavior of DMS.

The above described subtle interplay between localizing and delocalizing tendencies in these models often cannot be described correctly in a mean-field picture, but it requires nearly exact methods such as Monte Carlo (MC) techniques. A well-known example are the manganites in the low-doping regime, where mean-field approximations and exact solutions lead to different conclusions about the nature of the metal-insulator transition (MIT). The cluster formation and percolative picture that has emerged for these systems cannot be captured properly with mean-field approximations. Therefore, here the model introduced in Ref.6 for DMS in the low-doping limit describing electrons in an impurity band is investigated using MC methods, and the results will be compared to mean-field calculations.

Although the system under study is actually hole-doped, it is described in terms of an electron-doped material. Therefore, the acceptor level is treated as a donor level lying approximately 100 meV below the conduction band. Furthermore, only the impurity band itself will be considered, whereas the valence (conduction) band states are neglected. This should be a good approximation as long as the Fermi energy of the “electrons” is considerably smaller than the gap between the donor state and the lowest conduction level. The simplest such model  takes into account the hopping between random Mn sites and the AF interaction between the impurity spins and the charge carriers. The disorder potential and Coulomb interactions are neglected and the model is given by

$`H_{DMS}`$ $`=`$ $`{\displaystyle \underset{i,j,\alpha }{}}(t_{ij}\mu \delta _{ij})(c_{i\alpha }^{}c_{j\alpha }+\text{H.c.})`$ (2)
$`+{\displaystyle \underset{i,j,\alpha ,\beta }{}}J_{ij}(c_{j\alpha }^{}{\displaystyle \frac{1}{2}}𝝈_{\alpha ,\beta }c_{j\beta })𝐒_i,`$
where $`c_{i\alpha }^{}`$ is the creation operator for an electron at the impurity site $`i`$ with spin $`\alpha `$, and $`t_{ij}`$ is the hopping amplitude between sites $`i,j`$. For two sites at distance $`r`$=$`|ij|`$ an exponential form, $`t(r)`$=$`2(1+r/a_B)\mathrm{exp}(r/a_B)`$ Ry, was assumed, with $`a_B`$ the Bohr radius associated with the impurity site and Ry its corresponding binding energy $`E_b`$. $`𝐒_i`$ represents the (assumed classical) spin at the Mn ($`3d^5`$) site with $`|𝐒_i|`$=$`5/2`$, $`𝝈`$ denotes the vector of Pauli matrices, $`\mu `$ the chemical potential, and $`J_{ij}`$ describes the non-local interaction between impurity sites and mobile holes. Owing to the localized nature of the Mn spins, it can be expressed as $`J_{ij}`$=$`J\mathrm{exp}(2r/a_B)`$, with $`J`$ an AF coupling constant. $`J`$ is chosen as the energy unit and periodic boundary conditions (PBC) are imposed. The effects of an external field can easily be introduced via a Zeeman term.

Following Ref.6, we assume the hole binding energy $`E_b`$=$`112.4`$meV (1Ry) and, additionally, we have $`J`$=15meV. It was checked that our conclusions below do not depend crucially on the exact values of the exchange coupling and, therefore, the results should be typical for the systems considered. The lattice constant $`a`$ of the GaAs lattice is known to be $`5.65`$Å, and $`a_B`$=7.8Å. With these values, the various nonlocal couplings $`J_{ij}`$ and hoppings $`t_{ij}`$ can easily be calculated once a random distribution of the impurity ions within the host lattice (of size $`L^3`$) has been chosen. In addition, these systems are presumably strongly compensated through As antisite donors resulting in a relatively small effective hole concentration of 10% of the Mn ions.

The previous mean-field study of $`H_{DMS}`$ resulted in some interesting conclusions about the magnetic behavior of DMS. In particular, an increase in $`T_C`$ by a factor of $``$2 was reported when a disorder Mn distribution was studied, compared with a periodic distribution. This is a challenging counterintuitive prediction that deserves to be tested with techniques beyond the HF approach, even more so as it is known that critical temperatures typically are overestimated by mean-field methods, as it happens, e.g., in the DE model. Fortunately, Hamiltonian (2) is especially suited for MC calculations, and the procedure used here has been extensively discussed in the study of the manganites.

The MC calculations were performed on a variety of lattice sizes, but the bulk of our work was done studying systems with $`N`$=80 impurity sites and a doping level of $`x`$$``$0.02, which lies in the neighborhood of the realistic MIT. Other lattice sizes ($`N`$=50-100), with slightly different values of $`x`$, yield similar results and lead to the same conclusions. The density $`p`$ of charge carriers is fixed at $`p`$=$`0.1`$ (with respect to the number of Mn dopants), unless otherwise noted. The random placement of Mn ions in the GaAs host lattice creates regions of high and low impurity density, respectively. Some of the results presented below are for only one such configuration of Mn impurities, but it was confirmed that the results for other configurations are within 10$`\%`$ of the data shown. Furthermore, for the configuration selected the HF critical temperature $`T_C^{HF}`$ is the same as for larger systems ($`N`$=200, 300) with similar values of $`x,p`$ and $`J`$. Most interesting is the magnetization behavior, with and without applied magnetic field. It is analyzed by measuring the magnetization per site, $`m`$=$`M`$/$`N`$ ($`M`$=$`\sqrt{_{ij}𝐒_i𝐒_j}`$), of the Mn spins, neglecting the contribution of the charge carriers. The latter are of much less importance as they have a spin of $`1/2`$ only, and also a density that is only a fraction of that of the Mn spins. Fig.1 compares $`m`$ for the MC as well as for the HF solution, with both results obtained for the same configuration of Mn impurities.

In this particular case $`T_C^{HF}`$$``$0.85 (similar to earlier results on larger systems ) for the impurity configuration considered, whereas the MC solution suggests an enhancement of FM correlations at $``$0.3($`\pm `$0.1). Similar values and comparable deviations are found for other configurations and system sizes as well (Fig.2(a)), which supports the notion that mean-field treatments generally overestimate $`T_C`$’s. At this same temperature, $`M`$ (scaled by $`|𝐒|`$) starts to deviate from the totally disordered case $`M`$=1 as well and also a pronounced change in d$`M`$/d$`T`$ is observed (Fig.2)(b), signalling the development of FM correlations. Furthermore, the concept of a disorder-enhanced $`T_C`$ does not appear in agreement with MC. Fig.1 points out the influence of disorder as it also includes HF and MC results for the periodic case. From the magnetization data as well as the magnetic fluctuation, $`\chi `$$``$$`M^2M^2`$, a critical behavior around $`T_C^{per}`$$``$0.25$`\pm `$0.05 can be deduced for the periodic system, which is slightly below the critical region for the disordered case. This appears in agreement with previous HF results , but the reduction in $`T_C`$ is much less pronounced. It will be argued below, however, that the actual transition temperature is lower than $`T_C^{per}`$. A thorough analysis reveals that the disordered system breaks apart into occupied and almost empty domains (in both approaches) and the main contribution to $`m`$ is entirely due to a single, large cluster, with average particle density $`n^{cl}`$$`>`$$`p`$. This leads to an increased temperature where FM correlations start building up compared to the periodic case since the effective local magnetic field is enhanced in this region. Single sites or very small clusters are virtually unoccupied and do not order at any finite temperature. For this reason the magnetization is never fully saturated in the MC approach as in the HF description, but reaches a plateau at about 15% of the possible value for a sample with the given Mn doping rate. If $`x`$ is selected considerably lower, it would be very unlikely for such a cluster to exist and there should be no finite saturation magnetization, as it is observed for samples with $`x`$$``$0.005. The conclusions above are also supported by analyzing the spin correlations between different clusters A, B (Fig.3).

This reveals that each cluster turns FM as the temperature is reduced, but the overall magnetization will remain low as the magnetization vectors for different domains point in different directions, whereas in HF they are required to be aligned along a chosen $`S_z`$-axis. In other words, the idea that delocalization of charge carriers leads to a spin alignment of the Mn ions across the whole sample, appears to be an artifact of the approximation employed. It also leads to the unphysical observation of totally aligned impurity spins in the extremely diluted limit $`x`$$`<`$0.001 at very low temperatures. But it is interesting to note from Fig.3 that the spin correlations between sites of a given cluster and all other sites assume a finite value at temperatures below $``$0.08 (10K). In addition, finite correlations develop between A, B at T$``$0.2. The exact value of this temperature likely depends on the relative position between two specific domains, but is in good agreement with the value $`T_C`$=0.08 quoted above, where all sites were taken into account. Additionally, $`M`$/$`|𝐒|`$ is proportional to the number of spins below this temperature as it should be for a FM ordered state (Fig.2(b)). In any case, these values are lower than $`T_C^{per}`$ and they suggest two different temperature scales in the problem: $`T^{}`$ ($``$ 0.3), where isolated magnetic domains are formed, and a smaller true critical temperature $`T_C`$ ($``$ 0.1) where long-range order is established. Below $`T^{}`$, spins in the vicinity of a given cluster start to align, i.e. the magnetic domain grows in size as the gain in kinetic energy outweighs the loss in entropy. Further delocalization of the charge carriers then causes the alignment of these domains resulting in finite bulk magnetization. The associated $`T_C`$ is considerably lower than the one predicted by mean-field theory. This idea is also qualitatively depicted in Fig.2(c). In this scenario disorder enhances $`T^{}`$, but reduces $`T_C`$.

These percolation-like effects are also visible once a magnetic field is applied. As shown in the inset of Fig.1, small fields of $`h`$$``$0.04 (5-10T) are sufficient to introduce a robust change in the total magnetization of the sample at low temperatures. From Fig.3, it can be concluded that this happens mostly through the mutual alignment of already preformed FM areas, whereas the magnetization of a given domain remains comparatively unaffected by those small fields. The susceptibilities should be even stronger for smaller values of $`x`$, especially in the insulating regime. A pronounced magnetoresistance effect was indeed observed in lightly doped, insulating InAs and GaAs systems in a similar temperature regime. The situation described here is reminiscent of the colossal effects that can be observed in transition metal oxides. In that context it was conjectured that disorder by chemical doping leads to cluster formation and concomitant large responses to small changes in various parameters.

The possible percolative nature of the magnetic transition in DMS has been discussed recently, and it was proposed that the relevant analog is the problem of randomly placed spheres which overlap and, thus, form an infinite connected area once the radius of the spheres is larger than the percolation radius $`r_{per}`$. In $`d`$=3, $`r_{per}`$=1.33$`r_o`$, with $`r_o`$ ($``$$`x^{1/3}`$) the interspin distance (and percolation radius) in the periodic case. The MC data ($`r_{per}`$$``$2-3 lattice units) agree well with this prediction. In this scenario an increase in $`x`$ will lead to a smaller $`r_{per}`$ and, therefore, to an increase in $`T_C`$. At some point, however, $`r_{per}`$ will be smaller than the typical cluster size at $`T^{}`$, which makes $`T^{}`$ the maximum critical temperature possible. This temperature, although slightly lower, is in the range of the observed $`T_C^{max}`$ and small changes in the parameters of $`H_{DMS}`$ might be sufficient to result in higher transition temperatures.

This can be demonstrated by altering the compensation rate. Fig.4 shows $`m`$ for the same impurity configuration as above, but with $`p`$=$`0.20`$ and $`p`$=$`0.05`$, respectively. As expected, a lowered compensation rate leads to enhanced magnetism. The analysis of the spin-correlations (inset Fig.4) shows that this is mostly due to an increased intercluster spin correlation, whereas the $`m`$ vs. $`T`$ curve for a given domain (not shown) is largely unchanged from the one at $`p`$=$`0.1`$. Thus, the increased number of charge carriers mainly leads to larger droplets, favoring the creation of a coherent FM state. Additionally, the $`m`$ vs. $`T`$ curve acquires a more Brillouin-like look. The opposite is true for the case of even stronger compensation, $`p`$=$`0.05`$, where a small reduction in $`T_C`$ is observed.

In summary, a model for DMS has been studied with MC methods at low temperatures. For the doping rates considered, the system can be described in terms of ferromagnetic droplets. A change in parameters such as temperature or external fields can lead to the mutual alignment of these domains, while $`T_C`$ is considerably reduced compared to a mean-field analysis. The HF predictions of a disorder-enhanced transition point are not in agreement with nearly exact MC methods. The predicted $`T_C`$ is somewhat lower than measured experimentally, but this value depends considerably on parameters that are not very well known. The inhomogeneous state is expected to show large susceptibilities, e.g., a pronounced magnetoresistance effect.

Discussions with M. Berciu, especially about the implementation of the Hartree-Fock method for (2), are greatfully acknowledged. Part of this work (E.D.) was supported by NSF grant DMR-0122523.