text:
# The Effect of Magnetic Field Tilt and Divergence on the Mass Flux and Flow Speed in a Line-Driven Stellar Wind

## 1 Introduction

In a recent paper (ud-Doula & Owocki 2002; hereafter UO-02), we presented numerical magnetohydrodynamic (MHD) simulations of the effect of a stellar dipole magnetic field on the line-driven stellar wind from a non-rotating, hot, luminous star. We showed that the overall influence of the field on the wind depends largely on a single, dimensionless, ‘wind magnetic confinement parameter’, $`\eta _{}`$ ($`=B_{eq}^2R_{}^2/\dot{M}v_{\mathrm{}}`$), which characterizes the ratio between magnetic field energy density and wind kinetic energy density. Because the field energy declines faster than the wind energy, in the outer regions the magnetic field is always dominated by the radial wind outflow, which thus asymptotically stretches the field into a radial configuration, regardless of the strength of $`\eta _{}`$. For weak confinement $`\eta _{}<1`$, this radial opening of the field applies throughout the whole computational domain. But for stronger confinement $`\eta _{}>1`$, the magnetic field remains closed near the surface and over a limited range of latitude around the magnetic equator.

In this paper we provide further analysis and interpretation of how key wind properties, namely the surface mass flux and asymptotic flow speed (see figs. 1a,b in §2 below), are modified by the presence of a magnetic field. One central motivation is to reconcile results of our numerical MHD simulations with the previous scaling analysis done by MacGregor (1988; hereafter M-88), which had predicted that the faster-than-radial divergence of magnetic flux tubes would lead to substantial (factor 3-4) increase in the terminal velocity of a line-driven wind, compared to that of the non-magnetic, spherical-expansion case. Since the implied flow speeds ranging up to 5000-6000 km/s are never observed through blue edges of UV P Cygni lines from hot stars (Prinja & Howarth 1986; Howarth and Prinja 1987), that prediction could be a basis for inferring that hot stars lack magnetic fields with sufficient strength to substantially influence their stellar winds.

Here we develop (§3) a simplified formulation of the one-dimensional (1D) equations for a steady line-driven wind. Applying this to cases with faster-than-radial area expansion (§4), we show that the prediction of a 3-fold increase in flow speed is the consequence of certain assumptions and idealizations (point-star radiation, phenomenological flow expansion, and large line power-index $`\alpha `$) in the previous M-88 analysis; when these assumptions are relaxed (to include finite-disk correction, a dynamical flow divergence, and more realistic $`\alpha `$) the expected effect on flow speeds is consistent with the more modest $`\mathrm{}<50\%`$ increase typically found in our MHD simulations. We further show that the inferred simulation scaling of the base mass flux with the square of the cosine of the surface field tilt angle can be understood from a simple modified 1D tilted flow analysis (§5). After a brief discussion of the implications for stellar wind structure (§6), we conclude with a summary (§7) of our main results.

## 2 Summary of Key Results of 2D MHD Simulations

Figures 1a,b summarize two key results of our 2D MHD simulations, wherein we study the dynamical competition between magnetic field and wind by inserting a dipole magnetic field of varied strength (parameterized by $`\eta _{}`$) into a previously relaxed spherically symmetric CAK wind.111 The simulation results given here are analogous to those presented by UO-02 and ud-Doula (2002), but correct a minor coding error that lead the polar flow to have an artifically enhanced mass flux and an artificially reduced flow speed. As in UO-02, the model here is assumed isothermal, with a constant sound speed $`a=26.3`$ km/s, implying a ratio $`sa^2/v_{esc}^2=0.0014`$ between the gas internal energy and surface escape energy. (See §3 and the Appendices.)  For the asymptotic states when the regions of magnetic field opened by the wind have a nearly steady-state flow, both figures show results for three cases, with the degree of magnetic confinement ranging from weak ($`\eta _{}=0.1`$), through moderate ($`\eta _{}=1`$), to strong ($`\eta _{}=10`$).

In fig. 1a, the dotted curves show the latitudinal variation of the base mass flux, while the solid curves show this same flux scaled by $`\mu _B^2`$, where $`\mu _B`$ is the radial projection of a unit vector along the base magnetic field. The near constancy of this scaled flux is explained in §5, which uses a simple tilted-flow analysis of the base line-driving to show that the radial mass flux does indeed vary with the square of the field’s radial projection cosine, $`\dot{M}\mu _B^2`$.

At the magnetic poles, the radial orientation of the field means the flow is also radial, with essentially the same base mass flux as in the spherical model.222 The few percent increase is a second-order effect, associated with the finite sound speed and the resulting detachment of the flow critical point from the surface, where faster areal expansion reduces the density, and so allows driving of a somewhat greater mass flux. See Appendix B.

Fig. 1b plots the latitudinal variation of the flow speed at the maximum model radius $`R_{max}=6R_{}`$, scaled by the corresponding speed in the spherical, non-magnetic ($`B=0`$) model. The lower speeds toward the magnetic equator ($`\mu 0`$) reflect the various relative degrees of equatorial flow compression, which increases the density, and so reduces the radiative acceleration. In the highest magnetic confinement case $`\eta _{}=10`$, the increased speeds at mid-latitudes ($`0.1<|\mu |<0.8`$) result from a combination of the reduced mass flux and greater flow expansion, both of which reduce the density, and so increase the radiative acceleration.

While such latitudinal trends thus make sense qualitatively, our focus below will be on developing a semi-quantitative understanding of the wind properties at the magnetic poles, where the radial orientation of the flow makes possible a 1D analysis, if one takes proper account of the non-spherical expansion of the flow area. This was essentially the approach taken by M-88, but the factor 3 or more speed enhancements found there are substantially greater than those in fig. 1b, which even for the strongest magnetic confinement $`\eta _{}=10`$ shows a polar speed enhancement of only $`1.54`$. The next two sections develop and apply a formalism for understanding this difference.

## 3 Flow Analyses for Spherical Expansion

### 3.1 CAK Model of a Wind Driven by a Point-Source of Radiation

To provide a basis for our analysis of how the tilt and areal divergence of a magnetic field modify a wind outflow, let us first recast the basic Castor, Abbott and Klein (1975; hereafter CAK) formalism in a somewhat more tractable, physically transparent form. First, for a steady, !D stellar wind driven radially by line-scattering of radiation from a central point-source, the associated equation of motion for the radial velocity $`v`$ as a function of radius $`r`$ takes the form

$$\left[1\frac{a^2}{v^2}\right]v\frac{dv}{dr}=\frac{GM(1\mathrm{\Gamma })}{r^2}+g_{CAK}+\frac{2a^2}{r},$$
(1)
where $`GM(1\mathrm{\Gamma })/r^2`$ is the effective gravitational acceleration, and the Eddington parameter $`\mathrm{\Gamma }\kappa _eL/4\pi GMc`$ accounts for the effect of electron scattering opacity $`\kappa _e`$ interacting with the stellar luminosity $`L`$, with $`c`$ the speed oflight. The square-bracket factor with the sound speed $`a`$ on the left-hand-side accounts for gas pressure effects that allow for smooth mapping of the wind outflow onto a hydrostatic atmosphere below the sonic point, where $`v=a`$. On the right-hand-side the flow’s areal expansion also gives rise to a pressure-form term, $`2a^2/r`$, which compared to the competing gravitational term is of order $`s(a/v_{esc})^20.001`$, where $`v_{esc}\sqrt{2GM(1\mathrm{\Gamma })/R_{}}`$ is the effective escape speed from the stellar surface radius $`R_{}`$. The outward driving against gravity through the sonic point and beyond must thus be principally via the line acceleration, which under CAK’s original assumption of radially streaming radiation from a point source can be written as

$$g_{CAK}=\frac{1}{(1\alpha )}\frac{\kappa _eL\overline{Q}}{4\pi r^2c}\left(\frac{dv/dr}{\rho c\overline{Q}\kappa _e}\right)^\alpha $$
(2)
where $`\rho `$ is the density, and $`\alpha `$ is the CAK exponent. Here the line normalization $`\overline{Q}`$ is related to the usual CAK “k” parameter through $`k=\overline{Q}^{1\alpha }\left(v_{th}/c\right)^\alpha /(1\alpha )`$, but offers the advantages of being a dimensionless measure of line opacity that is independent of the assumed ion thermal speed $`v_{th}`$, with a nearly constant characteristic value of order $`\overline{Q}10^3`$ for a wide range of ionization conditions (Gayley 1995).

To proceed let us define the gravitationally scaled inertial acceleration

$$w^{}\frac{r^2vdv/dr}{GM(1\mathrm{\Gamma })}.$$
(3)
In terms of an inverse radius coordinate $`x1R_{}/r`$, note that $`w^{}=dw/dx`$, where $`wv^2/v_{esc}^2`$ represents the ratio of wind kinetic energy to the effective gravitational binding $`v_{esc}^2/2GM(1\mathrm{\Gamma })/R`$ from the stellar surface radius $`R_{}`$. We can then rewrite the equation of motion (1) in the dimensionless form

$$(1s/w)w^{}=1+C(w^{})^\alpha +2s/(1x)$$
(4)
where $`sa^2/v_{esc}^2`$, and we have eliminated the density $`\rho `$ in favor of the mass loss rate $`\dot{M}=4\pi r^2\rho v`$ for the assumed steady, spherical expansion, with the line-force constant defined by

$$C\frac{1}{1\alpha }\left[\frac{L}{\dot{M}c^2}\right]^\alpha \left[\frac{\overline{Q}\mathrm{\Gamma }}{1\mathrm{\Gamma }}\right]^{1\alpha }.$$
(5)
Note that, for fixed sets of parameters for the star ($`L`$, $`M`$, $`\mathrm{\Gamma }`$) and line-opacity ($`\alpha `$, $`\overline{Q}`$), this constant scales with the mass loss rate as $`C1/\dot{M}^\alpha `$.

The essential elements of interest in this paper can be analyzed without including the small gas pressure terms, and so we shall henceforth consider only the limit of vanishing sound speed $`a\sqrt{s}0`$. (The corrections from a small but finite $`s`$ are examined in the Appendices.) For this limit, fig. 2 illustrates a graphical solution to eqn. (4) for various values of $`C`$. For high $`\dot{M}`$ or small $`C`$ there are no solutions, while for small $`\dot{M}`$ or high $`C`$ there are two solutions. The CAK critical solution corresponds to a maximal mass loss rate, and requires a tangential intersection between line-force and combined inertia plus gravity, for which

$$\alpha C_cw_{c}^{}{}_{}{}^{\alpha 1}=1,$$
(6)
and thus

$$w_c^{}=\frac{\alpha }{1\alpha },$$
(7)
with

$$C_c=\frac{1}{\alpha ^\alpha (1\alpha )^{1\alpha }}.$$
(8)
Applying this in eqn. (5), we can then solve for the mass loss rate, yielding the standard CAK scaling

$$\dot{M}_{CAK}=\frac{L}{c^2}\frac{\alpha }{1\alpha }\left[\frac{\overline{Q}\mathrm{\Gamma }}{1\mathrm{\Gamma }}\right]^{(1\alpha )/\alpha }.$$
(9)
Moreover, since the scaled equation of motion (4) has no explicit spatial dependence, the scaled critical acceleration $`w_c^{}`$ applies throughout the wind. This can thus be trivially integrated to yield $`w=\alpha x/(1\alpha )`$, which is equivalent to the CAK velocity law

$$v(r)=v_{\mathrm{}}\left(1\frac{R_{}}{r}\right)^{1/2},$$
(10)
where the wind terminal speed $`v_{\mathrm{}}=v_{esc}\sqrt{\alpha /(1\alpha )}`$ scales with $`v_{esc}\sqrt{2GM(1\mathrm{\Gamma })/R_{}}`$, the effective escape speed from the stellar surface radius $`R_{}`$.

### 3.2 The Finite-Disk Correction Factor

The above analysis has so far been based on the idealization of radially streaming radiation, as if the star were a point-source at the origin. This was the basis of the original CAK wind solutions, although they did already identify (but did not implement) the appropriate “finite-disk correction factor” to account for the full angular extent of the star \[See CAK eqn. (50).\],

$$f_d=\frac{(1+\sigma )^{1+\alpha }(1+\sigma \mu _{}^2)^{1+\alpha }}{(1+\alpha )\sigma (1+\sigma )^\alpha (1\mu _{}^2)},$$
(11)
where $`\sigma (r/v)dv/dr1`$ and $`\mu _{}^21R_{}^2/r^2`$. When this factor is included to modify the point-star CAK line-acceleration (2), its complex dependence on radius, velocity, and velocity gradient complicates somewhat the solution of the full equation of motion (1). But full solutions have been derived independently by Friend and Abbott (1986) and Pauldrach, Puls, and Kudritzki (1986). They yield a somewhat reduced mass loss rate $`\dot{M}_{fd}\dot{M}_{CAK}/(1+\alpha )^{1/\alpha }`$, higher terminal speed $`v_{\mathrm{}}3v_{esc}`$, and flatter velocity law, roughly characterized by replacing the exponent $`1/2`$ in eqn. (10) with a somewhat higher value $`\beta 0.8`$.

If one simply assumes such a canonical “beta” velocity law, fig. (3) illustrates the specific spatial variations of $`f_d(x)`$ for various CAK exponents $`\alpha `$ and velocity exponents $`\beta `$. All these varied cases show a quite common general behavior. For the initial acceleration from a static stellar surface ($`rR_{}`$, or $`x0`$; where $`v0`$ with finite $`dv/dr`$, implying $`\sigma \mathrm{}`$), $`f_d`$ has a starting value $`f_d1/(1+\alpha )`$. It then increases outward, past unity at the point where the wind has a locally isotropic expansion (with $`dv/dr=v/r`$, so that $`\sigma =0`$), and then finally falls back asymptotically toward unity at large radii ($`x1`$), where the star is indeed well approximated by a point source.

To understand the dynamical effect of this finite-disk factor, and indeed to see how it interacts with effects of a magnetic field, it is helpful to examine how such a typical radial variation influences the scaled acceleration $`w^{}`$ through solutions to the appropriately modified equation of motion \[cf. eqn (4)\],

$$w^{}=1+fC_{}(w^{})^\alpha ,$$
(12)
where for convenience we have defined scaled quantities $`ff_d/f_d`$ and $`C_{}f_dC`$, and, for simplicity, again taken the zero-sound-speed limit $`s=0`$. Since $`f`$ increases outward from unity at the surface, the surface radius $`r=R_{}`$ ($`x=0`$) now represents a critical point (also known as a choke or throat) that fixes the required minimal value of $`C_{}=C_c=1/\alpha ^\alpha (1\alpha )^{1\alpha }`$ to allow an accelerating solution from this wind base. Accounting then for the mass loss rate scaling of $`C_{}1/\dot{M}^\alpha `$, we readily see that the net result is to reduce the maximal allowed mass loss to

$$\dot{M}_{fd}=f_d^{1/\alpha }\dot{M}_{CAK}=\dot{M}_{CAK}/(1+\alpha )^{1/\alpha }.$$
(13)
This thus provides a simple rationale for this aspect of the detailed numerical solutions of Friend and Abbott (1986) and Pauldrach et al. (1986).

As the finite-disk factor increases from its reduced base value, the lower mass loading allows for stronger CAK line-force, with now the resulting greater acceleration being “leveraged” by its effect in further increasing the line-force. From the graphical solution in fig. 2, one sees that a modest increase in the line-acceleration term \[$`C(w^{})^\alpha `$\] above the critical solution (with $`C=C_c`$) leads to a solution with a much larger scaled acceleration $`w^{}`$. \[There also arises a much shallower acceleration, but as discussed by Feldmeier et al. (2002), maintaining this beyond the critical point requires a special kind of outer boundary condition, whereas the more rapid acceleration is compatible with the standard assumption of simple expansion into a vacuum.\]

Fig. 4 illustrates this strong sensitivity to the finite-disk correction in terms of $`\sqrt{w^{}}`$, which sets the scaling of terminal speed $`v_{\mathrm{}}v_{esc}\sqrt{w^{}}`$. Note how a modest, less-than-factor-two change in $`f`$ leads to large increases in this speed scale $`\sqrt{w^{}}`$, with the trends being most dramatic for the largest $`\alpha `$. From eqn. (12), we see that in the limit of large acceleration $`w^{}1`$, the speed scaling approaches the form

$$\sqrt{w^{}}\left(C_cf\right)^{1/2(1\alpha )};w^{}1,$$
(14)
which is illustrated through the dashed curves in fig. 4. For the specific case $`\alpha =1/2`$, we can use the quadratic formula to obtain a simple expression that applies to the full range,

$$\sqrt{w^{}}=f\left[1+\sqrt{11/f^2}\right],$$
(15)
where the choice of the “plus” root gives the steeper of the two acceleration solutions.

To obtain specific values for the flow speed, we can numerically integrate the scaled equation of motion (12) accounting for the complex functional dependence of the finite disk factor (11) on radius, velocity, and velocity gradient. The left panel in fig. 5 shows how the ratio of terminal to escape speed depends on the CAK power index $`\alpha `$; the right panel plots the full radial variation of velocity for the $`\alpha =0.6`$ case. The dashed curves compare results for a similar integration when a canonical beta velocity law is used to evaluate the finite-disk factor. For $`\beta =0.75`$ the curves are indistinguishable from the full calculation, and other values of $`\beta `$ also give very good approximations.

A basic conclusion is thus that one may derive accurate finite-disk-corrected wind solutions using a beta velocity law to obtain spatially explicit approximations of the finite-disk factor, with especially good agreement obtained for $`\beta 0.75`$.

## 4 The Effect of Non-Spherical Area Expansion

### 4.1  Kopp-Holzer Expansion in a Point-Star CAK Wind

The above formulation of standard line-driven wind theory provides a convenient basis for understanding the results of models extended to include magnetic fields. As first analyzed by MacGregor (M-88), one potentially key effect regards the faster-than-radial expansion of the area in a flow tied to a rapidly diverging magnetic field. Following Kopp & Holzer (1976), it is convenient to consider an infinitesimal flow tube with a non-radial expansion, but still with a radial orientation, such as would apply, for example, to flow emanating from the poles of a surface magnetic dipole. The non-radial expansion factor can then be written as

$$h(r)\frac{A(r)}{A_{}}\frac{R_{}^2}{r^2}=\frac{B(R_{})}{B(r)}\frac{R_{}^2}{r^2},$$
(16)
where the latter equality follows from the $`𝐁=0`$ property of any magnetic field, since that then implies $`d(BA)/dr=0`$ and so $`A1/B`$ along such a radial flow tube.

In the absence of a detailed dynamical solution for the magnetic field divergence, M-88 adopted the Kopp & Holzer (1976) phenomenological non-spherical expansion factor333 We use the notation $`h`$ instead of Kopp and Holzer’s $`f`$ to distinguish from the finite-disk correction defined above.,

$$h_{KH}=\frac{h_{max}\mathrm{exp}[(rR_1)/\sigma ]+1(h_{max}1)\mathrm{exp}[(R_{}R_1)/\sigma ]}{\mathrm{exp}[(rR_1)/\sigma ]+1},$$
(17)
where $`R_1`$, $`\sigma `$, and $`h_{max}`$ are free parameters, the first two specifying the location and extent of the non-spherical expansion region, and the last fixing the total overall level of non-spherical expansion. Following M-88, we fix sample values of the two Kopp-Holzer parameters, $`R_1=1.25R_{}`$ and $`\sigma =0.1`$, and compare models with various $`h_{max}=`$ 1.5 and 2.

Applying this Kopp & Holzer expansion factor to CAK line-driven wind models computed in the idealized point-star approximation, M-88 found substantial (more than factor three) increases in wind terminal speed for even modest values of the total non-spherical expansion, $`h_{max}2`$. The above formalism provides a convenient way to understand those results, and indeed to extend them to account for the combined effects of a non-spherical expansion and a finite-disk correction. (See §4.2.) Noting from eqn. (2) that the line acceleration scales with the density $`\rho `$ as $`g_{CAK}1/\rho ^\alpha `$, we see from the steady-state mass conservation relation $`\rho 1/vAh(r)/vr^2`$ that along a radial flow with non-spherical expansion (“nse”) the line force has the form

$$g_{nse}=h(r)^\alpha g_{CAK}.$$
(18)
The relevant equation of motion thus takes the same general form as given in eqn. (12) if we simply replace the finite-disk correction with the non-spherical expansion, $`fh^\alpha `$. With this redefinition, the dependence of the scaled acceleration $`w^{}`$ on $`h^\alpha `$ is thus again given by fig. 4.

Fig. 6a illustrates results for this case, showing the spatial variation (in terms of the inverse radius coordinate $`x`$) of $`\mathrm{log}w^{}`$ (left) and wind speed ratio $`\sqrt{w}=v/v_{esc}`$ (right). (The latter is obtained by straightforward numerical integration of the acceleration $`w^{}`$ over $`x`$.) In both panels, the dotted and dashed curves represent respectively the point-star and finite-disk solutions for spherical expansion, while the dot-dash and solid curves respresent the corresponding non-spherical solutions for this standard M-88 case with $`\alpha =0.7`$ and $`h_{max}=2`$. Note that both the acceleration and velocity laws in this point-star, non-spherical expansion case (dot-dash curve) agree roughly to what is found in the corresponding spherical, finite-disk model (dashed curve). In this sense, the M-88 results for this $`h_{max}`$ can be roughly thought of as simply substituting the non-spherical expansion for the finite-disk correction.

For this case with $`h_{max}=2`$ we find in particular a terminal speed ratio $`v_{\mathrm{}}/v_{esc}=5.2`$, which, for the M-88 assumed stellar parameters that set the effective escape speed $`v_{esc}=974`$ km/s, yields a terminal speed $`v_{\mathrm{}}5145`$ km/s, or slightly above the value $`v_{\mathrm{}}=4799`$ km/s quoted in table 1 of M-88. (The modest difference is due to the neglect here gas pressure terms that were included in M-88’s more elaborate analysis; see Appendix A and B.) When compared to the terminal speed expected for point-star spherical expansion $`v_{\mathrm{}}=\sqrt{\alpha /(1\alpha )}v_{esc}=1.53\times 974=1488`$ (which is slightly below the value $`1512`$ km/s found in M-88’s solution with gas pressure included), we see that non-spherical expansion with $`h_{max}=2`$ leads to a speed increase of 3.46.

Table 1 summarizes the speed enhancement factors for this case, as well as for a series of other models to be described below.

A key point here is that the strong flow speed increases M-88 obtained from including non-spherical expansion in the point-star CAK model are quite analogous to the speed increases obtained by including a finite-disk correction. Let us next extend our analysis to include both effects.

### 4.2 Combined Effects of Non-Spherical Expansion and a Finite Disk

For a radial flow the combined modification to the line force due to finite disk and non-spherical expansion is given simply by the product of these separate444 Actually, a faster-than-radial expansion generally also implies a larger lateral velocity gradient, and thus alters the finite-disk correction factor from its usual spherical form (11). In practice, we find this to be a minor correction, except very near the stellar surface, where it leads to a slightly increased mass loss rate; see Appendix B.  correction factors,

$$g_{fd,nse}=f_dh^\alpha g_{CAK}.$$
(19)
We can thus analyze the combined dynamical effect by solving the equation of motion given again by eqn.(12), applying now the simple substitution $`f(f_d/f_{})h^\alpha `$. With this combined rescaling, fig. 4 again gives the dependence of $`\sqrt{w^{}}`$ on the combined correction factor $`f`$.

The solid curves in fig. 6a give the acceleration and velocity for the full solution of combined non-spherical expansion plus finite-disk. Compared to the above point-star case, note the acceleration is now very strong, with indeed the velocity in this case going off the scale, approaching an extaordinarily large terminal speed $`v_{\mathrm{}}=13.6v_{esc}=13,255`$ km/s! However, note that even the spherical, finite-disk model in this case also has a very large terminal speed of $`v_{\mathrm{}}=6.13v_{esc}=5971`$ km/s. As such, in relative terms of non-spherical to spherical enhancement, the finite-disk case implies a terminal speed factor of “only”2.22, which is actually significantly less than the factor 3.46 of the M-88 point-star assumption for this case.

Nonetheless, in absolute terms the implied terminal speeds well in excess of 10,000 km/s in such a model are several times higher than fastest speeds (ca. 5000 km/s) inferred from blue-edges of the P-Cygni wind lines from “normal” OB stars (Prinja & Howarth 1986; Howarth and Prinja 1987). In this sense, if finite-disk effects had been included in M-88’s original analysis, the implications of assuming hot-star magnetic fields would have seemed even more dramatic; indeed, the lack of observational evidence for such high speeds may well have been interpreted as strong evidence against existence of a dynamically significant magnetic field (i.e. with confinement parameter $`\eta _{}>0.1`$).

But now, with the benefit of the full MHD simulations of UO-02, we have a strong incentive to investigate the situation further, and indeed the simple formalism here makes it possible to readily examine the effect of varying parameters, as summarized in table 1.

#### 4.2.1 Effect of More Modest CAK Exponent $`\alpha `$

The analysis in §3, and particularly the results in figs. 4 and 5, suggest that a key parameter in influencing the wind flow speed is the CAK power-index $`\alpha `$. Thus, in comparing the M-88 results with those of the UO-02 MHD simultions, it is important to note that the former assumed $`\alpha =0.7`$ and the latter $`\alpha =0.6`$.

Fig. 6b thus provides results for the same Kopp-Holzer parameters as used in M-88 (and in fig. 6a), but now with $`\alpha =0.6`$. Note that the accelerations and speeds are all reduced. In particular, the case with non-spherical expansion combined with finite-disk now has a terminal speed of $`v_{\mathrm{}}=6.28v_{esc}=6113`$ km/s, representing a more modest factor 1.71 enhancement over the terminal speed $`v_{\mathrm{}}=3.67v_{esc}=3570`$ km/s of the corresponding spherical, finite-disk model.

#### 4.2.2 Using Dynamically Motivated Values for the Kopp-Holzer Parameters

To further facilitate comparisons with our MHD simulations results, let us next consider Kopp-Holzer parameter values that are adjusted to fit the non-spherical expansion found in these simulations. Using eqn. (16) to convert magnetic field variations into flow area, fig. 7 plots the radial variation of the polar expansion factor found in the MHD simulations assuming various magnetic confinement parameters $`\eta _{}`$. The dashed curve shows the expected result for pure dipole in limit of large confinement parameter, $`\eta _{}\mathrm{}`$. The lower, middle and upper solid curves show the MHD results for for $`\eta _{}=`$ 0.1, 1, and 10. The dotted curves plot the best-fit Kopp-Holzer parameterizations, quoted in fig. 7 of UO-02.

For the Kopp-Holzer fit to the strongest confinement case $`\eta _{}=10`$, fig. 6c shows the enhancements in acceleration and speed are again quite modest. Indeed, at the radius $`R_{max}=6R_{}`$ of the outer boundary of the UO-02 MHD simulation, the speed of this 1D model is $`v_6=5.72v_{esc}=5576`$ km/s. Relative to the speed of the spherical case at this radius, $`v_6=3.29v_{esc}`$ =3200 km/s, this respresents an enhancement factor of about $`1.56`$, which is now remarkably close to the factor $`1.54`$ found for the polar flow of the 2D MHD model (cf. fig. 1b).

Table 1 shows similarly good agreement between this modified Kopp-Holzer parameterization and the MHD simulations (fig. 1b) that show polar speed enhancements of ca. 1.12 and 1.32 for confinement parameters $`\eta _{}=`$ 0.1 and 1.

### 4.3 A Dynamically Motivated Non-Spherical Expansion Function

The overall variation of the MHD non-spherical expansion factors plotted in fig. 7 (and in fig. 7 of UO-02) suggests a simpler, more physically motivated fitting-function,

$$h_{UO}(r)=\frac{r}{R_{}}\frac{R_{}+R_c}{r+R_c},$$
(20)
which has just a single parameter, $`R_c`$, representing a “magnetic confinement radius”. This characterizes the transition between the dipole scaling $`hr`$ that applies for a strong field near the surface, and the asymptotic monopole scaling $`h`$ constant that applies at large radii, where the wind has stretched the field into a nearly radial configuration. Note that in the opposite limits $`R_c0`$ or $`R_c\mathrm{}`$ the monopole or dipole scalings apply throughout the flow, as appropriate to the limits of weak ($`\eta _{}0`$) or strong ($`\eta _{}\mathrm{}`$) magnetic confinement.

This suggests a direct scaling of $`R_c`$ with the confinement parameter $`\eta _{}`$. In general, the ratio of magnetic to wind-kinetic energy should scale as $`\eta (r)=`$ $`\eta _{}(r/R_{})^{22p}`$, where $`B(r)1/r^p`$, with $`p=3`$ for a dipole, and $`p=2`$ for a monopole. If we then assume the confinement radius occurs where this ratio has some order-unity value $`\eta _c`$, we obtain a general scaling of the form

$$\frac{R_c}{R_{}}\left[\frac{\eta _{}}{\eta _c}\right]^q,$$
(21)
where $`q=1/(2p2)`$ should be intermediate between the dipole ($`p=3`$) value $`q_d=1/4`$ and the monopole ($`p=2`$) value $`q_m=1/2`$. Taking the simple case that the order-unity constant $`\eta _c`$ is in fact unity, we find quite good fits for an exponent $`q=3/8`$, representing the “geometric mean” between the dipole and monopole scalings for $`R_c`$,

$$R_c\eta _{}^{3/8}R_{}.$$
(22)
The dot-dash curves in fig. 7 plot the radial variations of this simple, physically motivated form of the single-parameter non-spherical-expansion function (20). Note that the agreement with the corresponding MHD simulation results (solid) is almost as good as the best-fit, three-parameter, Kopp-Holzer form \[dotted; eqn. (17)\] that has less physical justification.

Fig. 6c shows the acceleration and speed when this simple scaling is applied to the strong confinement case with $`R_c=10^{3/8}R_{}=2.37R_{}`$. At the simulation outer radius $`R_{max}=6R_{}`$, the speed ehancement of this non-spherical, finite-disk model is now a factor of about $`1.51`$, which is again remarkably close to the factor $`1.54`$ found for the polar flow of the 2D MHD model (cf. fig. 1b). Table 1 shows similarly good agreement between this simpler non-spherical enhancement factor and the MHD polar speed enhancements for the other confinement parameter cases.

## 5 Mass Loss for Tilted Flow

In addition to modifying the flow divergence, a magnetic field can also alter a flow’s direction. In particular, for a strong field the initial flow off the stellar surface occurs at an angle tilted to the local radial direction. In this section we examine how such a base flow tilt should influence the local mass flux.

Since the mass loss is set very near the stellar surface, where combined factor from non-spherical divergence and finite-disk correction has its minimum, let us ignore both sphericity and magnetic divergence effects and consider the field and flow to be along a fixed direction $`\widehat{𝐬}`$, assumed to be tilted by an angle $`\theta _B`$ with respect to the local vertical direction, defined by unit vector $`\widehat{𝐳}`$. We then have $`\mu _B\mathrm{cos}\theta _B=\widehat{𝐬}\widehat{𝐳}`$, and in terms of the local horizontal direction $`\widehat{𝐱}`$, we can also write $`\widehat{𝐬}=\mu _B\widehat{𝐳}+\sqrt{1\mu _B^2}\widehat{𝐱}`$.

For our usual simplification of negligible gas pressure, the steady-state equation of motion along the fixed field direction $`\widehat{𝐬}`$ becomes \[cf. eqn. (1)\]

$$(𝐯)(\widehat{𝐬}𝐯)=g_{}(\widehat{𝐬}\widehat{𝐳})+\widehat{𝐬}𝐠_{CAK},$$
(23)
where on the right-hand-side $`𝐠_{CAK}`$ represents the vector radiative force, and we have used the fact that surface gravity $`g_{}(1\mathrm{\Gamma })GM/R_{}^2`$ is along the local vertical $`\widehat{𝐳}`$. On the left-hand-side the advective derivative of velocity along $`\widehat{𝐬}`$ can be written as

$`(𝐯)(\widehat{𝐬}𝐯)`$ $`=`$ $`v_s(\widehat{𝐬})(v_s)`$ (24)
$`=`$ $`\mu _Bv_s{\displaystyle \frac{v_s}{z}}+\sqrt{1\mu _B^2}v_s{\displaystyle \frac{v_s}{x}}`$
$`=`$ $`\mu _Bv_s{\displaystyle \frac{v_s}{z}},`$
where the last expression assumes that there are no horizontal variations in velocity.

For the line-driving, let us again assume a point-source approximation with the radiation confined to a pencil beam along the local vertical $`\widehat{𝐳}`$. In terms of an associated surface flux $`F_{}L/4\pi R_{}^2`$, the projected CAK line-force term then becomes \[cf. eqn.(2)\]

$`\widehat{𝐬}𝐠_{CAK}`$ $`=`$ $`{\displaystyle \frac{1}{1\alpha }}{\displaystyle \frac{\overline{Q}\kappa _eF_{}}{c}}(\widehat{𝐬}\widehat{𝐳})\left[{\displaystyle \frac{\widehat{𝐳}[\widehat{𝐳}𝐯]}{\overline{Q}\kappa _e\rho c}}\right]^\alpha `$ (25)
$`=`$ $`{\displaystyle \frac{\mu _B}{1\alpha }}{\displaystyle \frac{\overline{Q}\kappa _eF_{}}{c}}\left[{\displaystyle \frac{\mu _Bv_s/z}{\overline{Q}\kappa _e\rho c}}\right]^\alpha `$
$`=`$ $`{\displaystyle \frac{\mu _B}{1\alpha }}{\displaystyle \frac{\overline{Q}\kappa _eF_{}}{c}}\left[{\displaystyle \frac{\mu _Bv_sv_s/z}{\overline{Q}\kappa _e\dot{m}_sc}}\right]^\alpha `$
$`=`$ $`{\displaystyle \frac{\mu _B}{1\alpha }}{\displaystyle \frac{\overline{Q}\kappa _eF_{}}{c}}\left[{\displaystyle \frac{\mu _B^2v_sv_s/z}{\overline{Q}\kappa _e\dot{m}_rc}}\right]^\alpha ,`$
where the third equality uses mass conservation along the tilted flow to eliminate the density $`\rho `$ in favor of the mass flux in the tilted direction, $`\dot{m}_s\rho v_s`$. The final equality then uses the fact that the associated radial mass flux component $`\dot{m}_r`$ (which sets the associated mass loss through the spherical surface) is reduced by another factor $`\mu _B`$, i.e. $`\dot{m}_r=\mu _B\dot{m}_s`$. After dividing through by $`\mu _B`$ and making the modified definition \[cf. eqn.(3)\],

$$w^{}\frac{v_sdv_s/dz}{g_{}},$$
(26)
the momentum equation (23) can again be recast in the familiar form \[cf. eqn. (4)\],

$$w^{}=1+C\mu _B^{2\alpha }(w^{})^\alpha ,$$
(27)
where $`C`$ is defined in eqn. (5), and is related to the maximal value of the radial mass flux. Equation (27) is mathematically identical to eqn. (4) under the substitution $`CC\mu _B^{2\alpha }`$. As such, following the same procedure as in §2.1, we find the surface mass flux scales as

$$\dot{m}_r=\mu _B^2\dot{m}_{CAK}$$
(28)
where $`\dot{m}_{CAK}\dot{M}_{CAK}/4\pi R_{}^2`$ is the analogous surface mass flux in the CAK model for spherical mass loss.

The $`\mu _B^2`$ reduction of the radial surface mass flux predicted by this simplified planar, 1D flow analysis is actually in quite good agreement with results of our full 2D MHD simulations (Ud-Doula 2002; UO-02). Indeed, the solid curves of fig. 1a above show that the mass flux scaled by $`\mu _B^2`$ is roughly constant in latitude, and equal to the mass flux in the spherical case. Physically, one factor of $`\mu _B`$ can be attributed simply to geometric projection of the tilted flow onto the vertical normal to the star’s spherical surface, while the other can be seen to result from the tilted flow acceleration having less effect in desaturating the optically thick absorption of the radially streaming radiation.

This lower initial mass flux combined with the faster-than-radial area expansion leads to a lower density, which through the CAK line-force leads to a greater acceleration and thus a faster terminal speed. In cases with a large magnetic confinement parameter (i.e. $`\eta _{}=10`$), the simulations thus tend to show the greatest terminal speeds at latitudes just away from the magnetic equator, which originate from the open field lines that are furthest from the magnetic pole, and for which the initial flow tilt at the surface wind base is greatest. (See fig. 1b here and right panel of fig. 8 in UO-02).

## 6 Discussion

The mass flux and speed variations arising in such magnetic models of line-driven winds provide a potentially attractive mechanism for explaining the extensive wind structure and variability commonly inferred from observations of UV P-Cygni lines from OB stars (e.g., Kaper et al. 1996; Prinja et al. 1998). If the models here were extended to include even a modest stellar rotation along an axis that is tilted relative to the magnetic dipole axis, then the speed variations with magnetic latitude would naturally lead to Co-rotating Interaction Regions (CIRs), much as occur between faster and slower streams in the solar wind. Following initial suggestions by Mullan (1984), such CIRs now represent a popular general paradigm for explaining the commonly observed “Discrete Absorption Components” (DACs) in P-Cygni absorption troughs. But in previous hydrodynamical simulations by Cranmer and Owocki (1996), the speed variations leading to CIR formation had to be induced somewhat artificially, by assuming a modulation of the base wind driving, as might result from, e.g., bright or dark spots on the stellar surface. Since the analyses here (and the UO-02 MHD simulations) show that even modest magnetic confinement parameters can lead to quite substantial latitudinal variations in flow speed, a tilted dipole field seems a much more likely mechanism for inducing such wind CIRs.

Note that such CIRs would also tend to limit the highest speeds occuring from the most rapidly diverging flow, simply from the truncation by the interaction with the slower speed flow. As such the cases with the highest overall flow speeds would be those with a strong, rotation-aligned dipole, for which the high speed polar wind wouldn’t be brought into interaction with the slower, equatorial wind. Such cases may be relatively rare, but are possible. The flow speeds would be quite high, perhaps up to twice that of a spherical wind. But the analysis here shows that they needn’t be absurdly high, at least for moderate values of the CAK power index, i.e. $`\alpha 0.6`$.

## 7 Concluding Summary

This paper analyzes the effects of the flow tilt and non-spherical expansion associated with magnetic channeling on the mass flux and flow speeds of a line-driven stellar wind. Our main results are summarized as follows:

1. A faster-than-radial expansion leads to an enhancement in wind speed, but the relative corrections are typically about 50% when compared to a finite-disk-corrected model with a moderate CAK power-index $`\alpha =0.6`$, much less than the factor three or more inferred by M-88 in their analysis of a point-star model with $`\alpha =0.7`$. The analysis here is in good quantitative agreement with results from numerical MHD simulations.
2. The non-radial expansion obtained from our numerical MHD models differs from the heuristic form used by Kopp and Holzer (1976), with the rapid expansion beginning right at the wind base, and extending over a quite large radial distance. We propose a simpler, physically motivated fit function \[eqn. (20)\] controlled by a single parameter, $`R_c`$, representing a characteristic radius for wind magnetic confinement, and scaling as $`R_c\eta _{}^{3/8}`$ with the magnetic confinement parameter $`\eta _{}`$ ($`B_{eq}^2R_{}^2/\dot{M}v_{\mathrm{}}`$).
3. The radial mass flux at the stellar surface with a tilted magnetic field is reduced by the square of the cosine of the tilt-angle ($`\mu _B^2`$) compared to a non-magnetic, spherical wind. Physically, one factor of $`\mu _B`$ stems from the geometric projection of the tilted flow onto a surface normal, while the other arises from the reduced desaturation of absorption of radially streaming radiation by acceleration along this tilted flow.
4. A perturbation analysis (in Appendix A and B) shows that the corrections from a small, but non-zero gas pressure scale with the ratio of sound-speed to escape-speed, $`a/v_{esc}`$. Relative to a zero-sound-speed, finite-disk-corrected spherical wind, typical increases in the mass loss rate are 10-20%, with comparable relative decreases in the wind terminal speed. In non-spherical models, the additional gas pressure correction is typically just 1-2%.

Overall, the results confirm that even modest magnetic fields can have a substantial influence in line-driven winds, with faster areal divergence enhancing the wind acceleration and flow speed. But even the largest inferred enhancements still generally allow for speeds that are within the range of observational constraints. As such, strong magnetic fields in hot-star winds should not be precluded, as they might have been if previous analyses implying large speed enhancements had not been reduced by the more complete study here.

Acknowledgements. This research was supported in part by NSF grant AST-0097983, NASA grants NAG5-11886 and NAG5-11095, and the NASA Space Grant College program at the University of Delaware. SPO acknowledges support of a PPARC visiting fellowship and thanks J. Brown of the University of Glasgow and A. Willis of University College London for their hospitality during his sabbatical-year visits. We thank R. Townsend and D. Cohen for helpful discussions and comments.

## Appendix A  Non-Zero-Sound-Speed Correction for a Spherical, Finite-Disk Model

Let us examine here the effect of a finite gas pressure associated with a non-zero sound speed. In the zero-sound-speed limit considered in §§3-5, the minimal value of the finite-disk-correction factor, $`f_d=1/(1+\alpha )`$, occuring near the stellar surface, $`rR_{}`$, sets the critical conditions that determine the maximal allowed steady mass loss rate, $`\dot{M}_{fd}`$. But with a finite sound speed, the reduced effective inertia $`(1s/w)w^{}`$ in the base region near the sonic point shifts this critical location for most difficult driving slightly away from the surface wind base. But then, because of the somewhat larger finite-disk factor, this in turn allows a somewhat larger maximal mass loss rate.

To analyze these effects, let us recast the equation of motion in the form,

$$F(x,w,w^{})(1s/w)w^{}+1\frac{f}{(1+\delta m_s)^\alpha }C_c(w^{})^\alpha 2s/(1x)=0,$$
(A1)
where $`s(a/v_{esc})^2`$ is typically of order $`10^3`$, and

$$\delta m_s\frac{\dot{M}_{fd,fs}}{\dot{M}_{fd}}1$$
(A2)
is the fractional correction of the mass loss rate due to this small, but non-zero $`s`$. The new critical point where this mass loss is fixed now occurs somewhat away from the surface, $`x_c>0`$, but in the usual cases that $`s1`$, we expect both $`\delta m_s1`$ and $`x_c1`$.

In the narrow region from the sonic point to this critical point, the driving is dominated by the need to overcome gravity. This implies an outward acceleration that scales as the inward gravity, namely $`w^{}(x)w_c^{}=\alpha /(1\alpha )`$, which in turn implies that the velocity in this initial region nearly follows the standard CAK (“beta=1/2”) solution, $`w(x)w_c^{}x`$. Using this, we obtain a spatially explicit form for the finite-disk factor, which upon expansion to first order in $`xx_c1`$ becomes

$$f(x)1+4x.$$
(A3)
We then find the equation of motion (A1) takes the spatially explicit form

$$F(x)=\frac{(1\alpha )s}{x}+4x\alpha \delta m_s+2s=0,$$
(A4)
where we have also expanded $`1/(1+\delta m_s)^\alpha `$ and kept terms to first order in the smallness parameter $`s`$. The critical point is located where this function has a spatial minimum,

$$0=\left[\frac{dF}{dx}\right]_c=4\frac{(1\alpha )s}{x_c^2},$$
(A5)
which solves to

$$x_c\sqrt{(1\alpha )s}/2.$$
(A6)
Applying this in eqn. (A4) and keeping only the leading term of order $`\sqrt{s}`$, we find the mass loss correction is given by

$$\delta m_s\frac{4\sqrt{1\alpha }}{\alpha }\frac{a}{v_{esc}}.$$
(A7)
Beyond this critical point, the increase in the finite-disk factor means that solutions to the equation of motion quickly approach their large gradient limit, $`w^{1\alpha }fC^c/(1+\delta m_s)^\alpha `$. (See dashed curves in fig. 4.) Since the wind terminal speed scales as $`v_{\mathrm{}}\sqrt{w^{}}`$, this implies that the fractional speed corrections from a finite sound speed should scale as

$$\delta v_{\mathrm{},s}\frac{v_{\mathrm{},s}}{v_{\mathrm{},s=0}}1(1+\delta m_s)^{\alpha /2(1\alpha )}1\frac{\alpha \delta m_s}{2(1\alpha )}\frac{2}{\sqrt{1\alpha }}\frac{a}{v_{esc}}.$$
(A8)
For a typical case with $`\alpha =1/2`$ and $`s(a/v_{esc})^2=0.001`$, we find $`x_c0.012`$, $`\delta m_s0.18`$, and $`\delta v_{\mathrm{},s}0.09`$. For the same $`s`$ but $`\alpha =0.6`$, we get $`\delta m_s0.13`$ and $`\delta v_{\mathrm{},s}0.10`$, while for $`\alpha =0.7`$, we find $`\delta m_s0.10`$ and $`\delta v_{\mathrm{},s}0.12`$. Applying these corrections to the scalings from the zero-sound-speed calculations in §3 generally gives very good agreement (to about a percent) with the full solutions.

It is worth noting here that, even though the corrections for finite gas pressure are modest, they are not, perhaps, as small as one might have expected. That is, even though the ratio of internal energy to escape energy is of order $`s=(a/v_{esc})^2`$, the corrections scale with $`\sqrt{s}=a/v_{esc}`$. Moreover, this scaling is multiplied by a large coefficient, which stems from the relatively steep increase in the leveraged finite disk factor – e.g. scaling as $`f^{1/\alpha }(x)1+8x`$ for $`\alpha =1/2`$, equivalent to an $`r^8`$ dependence on the radius $`r`$. As a result, even though the reduced inertia of a mildly supersonic flow only leads to a relatively modest outward shift of the critical point from the stellar surface, the stronger finite-disk factor leads to a stronger radiative driving, and thus allows a larger mass loss. Finally, that larger mass loss also implies a comparable relative reduction in the wind terminal speed.

## Appendix B  Non-Zero-Sound-Speed Correction for a Non-Spherical, Finite-Disk Model

In the case of non-spherical expansion, the analysis follows much as above, except that now we must include the non-spherical expansion factor $`h(x)`$. In the region near and below the critical point $`xx_c1`$ this expands to

$$h(x)1+h_o^{}x,$$
(B1)
where, for the form defined in eqn. (20),

$$h_o^{}=\frac{R_c}{R_{}+R_c}1/(1+\eta _{}^{3/8}).$$
(B2)
This non-spherical expansion implies a stronger lateral velocity gradient, and for our assumed velocity law scaling ($`wx`$) this can be shown to modify the finite-disk correction expansion (A3) near the surface to

$$f_h(x)1+(4+h_o^{})x.$$
(B3)
Including then also the usual effect of this non-spherical expansion on the density, we make the replacement $`ff_hh^\alpha 1+[4+(1+\alpha )h_o^{}]x`$ and thereby find that equations (A4) and (A5) still apply if we simply make the substitution $`44+(1+\alpha )h_o^{}`$. The critical point thus now occurs at

$$x_c\sqrt{\frac{(1\alpha )s}{4+(1+\alpha )h_o^{}}},$$
(B4)
while the mass loss correction is now given by

$$\delta m_s\frac{2}{\alpha }\sqrt{(1\alpha )[4+(1+\alpha )h_o^{}]}\frac{a}{v_{esc}},$$
(B5)
and the speed correction by

$$\delta v_{\mathrm{},s}\frac{4}{\sqrt{(1\alpha )[4+(1+\alpha )h_o^{}]}}\frac{a}{v_{esc}}.$$
(B6)
For the weak, moderate, and strong magnetic confinement cases ($`\eta _{}=`$ 0.1, 1, and 10) we find $`h_o^{}`$ 0.3, 0.5, 0.7. For the typical $`s`$ of order $`10^3`$, we find the extra changes from non-spherical divergence range from one to two percent in both the mass loss rate and flow speed. In particular, applying the parameters $`s=0.0014`$ and $`\alpha =0.6`$ used in the UO-02 MHD simulations, we estimate extra mass loss enhancements for the weak, moderate, and strong-field cases to be respectively 0.9%, 1.5%, and 2.1%. These values are roughly consistent with the polar mass fluxes plotted in fig. 1a.


format: text

text: # 1 X-ray Properties of Nearby Abell Cluster Sample


text:
# Distribution of epicenters in the Olami-Feder-Christensen model

## Abstract

We found that the well established Olami-Feder-Christensen (OFC) model for the dynamics of earthquakes is able to reproduce a new striking property of real earthquake data. Recently, Abe and Suzuki found that the epicenters of earthquakes could be connected to generate a graph, with properties of a scale-free network of the Barabási-Albert type. However, only the non conservative version of the Olami-Feder-Christensen model is able to reproduce this behavior. The conservative version, instead, behaves like a random graph. Those findings, besides indicating the robustness of the model to describe earthquake dynamics, reinforce that conservative and nonconservative versions of the OFC model are qualitatively different, and propose a completely new dynamical mechanism that, without an explicit rule of preferential attachment, is able to generate a free scale network. The preferential attachment is in this case a “by-product” of the long term correlations built by the self-organized critical state. We believe that the detailed study of the properties of this network can reveal new aspects of the dynamics of the OFC model, contributing to the understanding of self-organized criticality in non conserving models.

self-organized criticality; earthquake; complex networks, complex systems

The concept of self-organized criticality was introduced by Bak, Tang and Wiesenfeld Bak  in 1987, as a possible explanation of scale invariance in nature. To illustrate their basic ideas, they presented a cellular automaton model, the sandpile model, so called because of a possible analogy between its dynamical rules and the movement of sand or snow in avalanches. Since this seminal work, a great number of cellular automata and coupled map models have been investigated, in an attempt to elucidate the essential mechanisms hidden in such a wide class of different non-linear phenomena whose statistics of events (or avalanches) are governed by power-laws. However, up to now, one still lacks from a general theoretical framework for self-organized criticality. Success in analytical investigations have been achieved in many models. For a revision see, for instance, book ; turcotte .

In this context, a model that has been widely studied in the literature is the Olami-Feder-Christensen (OFC) model for the dynamics of earthquakes. The original OFC model, introduced in 1992 OFC , is a two-dimensional coupled map model defined on a square lattice, whose dynamical rules were inspired in a spring-block model proposed to describe the dynamics of earthquakes. Earthquakes, in the real world, are associated with many power-laws, the most known of them being the Gutenberg-Ritcher law for the distribution of avalanche energies. The OFC model assigns - to each site of a square lattice - a real variable $`z_{i,j}`$ (energy or tension), initially chosen at random in the interval $`[0,z_c)`$, where $`z_c`$ is a threshold value. $`z_{i,j}`$ increases slowly throughout the lattice and each time that, for a given site, $`z_{i,j}`$ exceeds $`z_c`$, the system relaxes. A fraction $`\alpha z_{i,j}`$ of the tension of site $`(i,j)`$ is then distributed to each of its nearest neighbors. As a consequence, the tension of some of its neighbors may also exceed $`z_c`$, generating an ‘avalanche’ that will only stop when $`z_{i,j}<z_c`$ again for all sites of the lattice. We have assumed open boundaries in our simulations.

Within the OFC model there is a dissipation parameter $`\alpha `$. If $`\alpha =0.25`$ the total tension in the lattice, $`z_{i,j}`$, is conserved during the avalanching process, in the bulk of the lattice (there is always dissipation in the boundaries). But if $`\alpha <0.25`$ there is some dissipation also in the bulk of the system. Because of those facts, this model has been widely studied in literature: it is, at the same time, a prototype of self-organization in systems with non-conservative relaxation rules (the existence of SOC in the non-conservative models is, up to now, not well understood Prado ; Chris ; Lise ; Miller ) and also a paradigm of the success of SOC ideas, since it is able to reproduce important aspects of the dynamics of earthquakes.

Recently, Abe and Suzuki Abe1  observed a new power-law in the statistics of earthquakes. They analyzed earthquake data from both the district of southern California and Japan, connecting their epicenters in order to generate a graph. Each area analyzed was divided into small cubic cells; they associated to each of these cells a node every time an earthquake started inside it. The epicenters of two successive earthquakes were linked, defining an edge. In this way the data has been mapped into a complex growing graph that behaves like a scale-free network of the Barabási-Albert type Barabasi . The degree distribution of the graph decays as a power-law. The clustering coefficient and the diameter of a cluster were also calculated, showing small-world network properties Watts . These features have revealed a novel aspect of earthquakes as a complex critical phenomenon.

In this paper we studied the Olami-Feder-Christensen model to see if it could also predict this new striking behavior. We found that the non conserving version of the model reproduces the behavior of experimental data, even for a very small degree of non conservation. The degree distribution of the evolving network formed by its epicenters is scale free. However, the conservative version of the model has a qualitatively different behavior, more similar to a random graph, whose degree distribution is Poisson, indicating that most of the nodes have the same degree and - although random - the corresponding network is much more homogeneous. These results are in agreement with some recent observations, reinforcing that conservative and non conservative versions of the OFC model are quite different. Hergarten and Neugebauer (2002) Hergarten  studying the efficiency of the OFC model to predict foreshocks and aftershocks, de Carvalho and Prado (2003)Josue , studying the transient behavior of the OFC model and Miller and Boulter (2003) Miller03  studying the distribution of values at which supercritical sites topple have also reported qualitatively different behaviors between the conservative and non conservative OFC model.

In a complex graph, the edges are not distributed in an regular way and not all nodes have the same number of edges. One possible way to characterize complex networks is through its distribution function $`P(k)`$, which gives the probability that a random selected node has exactly $`k`$ edges. $`k`$ is called the degree of the node. In a random graph, since the edges are placed randomly among the nodes, the majority of nodes have approximately the same degree, close to the average connectivity $`k`$, and the distribution $`P(k)`$ is a Poisson distribution with a peak at $`P(k)`$. Most complex networks, however, have a distribution function $`P(k)`$ that deviates significantly from a Poisson distribution. In particular, for a large number of networks, associated with a wide class of systems, ranging from the world wide web to metabolic networks, $`P(k)`$ has a power-law tail, $`P(k)k^\gamma `$. Such networks are called scale free Barabasi , and have called the attention of many researchers in the last years.

We simulated the OFC model in a square lattice, building graphs with a procedure very similar to what has been employed by Abe and Suzuki. Each site that gives birth to a new avalanche is an epicenter; each epicenter defines a node, and every node is then connected to the node where the next epicenter occurs, establishing a link or edge between them. After many avalanches this procedure generates a complex network (or graph), and we have studied some of its statistical properties.

After eliminating a transient of at least $`10^6`$ events, we calculate numerically the distribution function $`P(k)`$ for the graph constructed from the time sequence of epicenters in the OFC model, for different values of $`\alpha `$ and different lattice sizes. As the first and last sites are the only ones with an odd number of edges, they were eliminated. Our results for the distribution $`P(k)`$ can be seen in figure 1. It is clear that, if $`\alpha <0.25`$ (figure 1a), the distribution is scale-free for some decades, with an exponent $`\gamma `$ that varies linearly with $`\alpha `$ (see figure 2), at least for values of $`\alpha `$ not too far from the conservative regime. The network grows toward the inside of the lattice, with the most connected sites in the borders and the most inner sites being the last ones being added to it (see figure 3a). The complex structure, however, is not a boundary effect. If we take out the border sites and adjusts the scale, we see that the same spatial structure is reproduced (figure 3b). Because one needs a growing network to observe the scale free-behavior Barabasi , after a certain number of events, as a consequence of the finite size of the lattice, most of the sites of the lattice have already become part of the network. At this point the scale free behavior starts to break.

If the system is conservative, however, the distribution function $`P(k)`$ has a well defined peak, indicating a higher degree of homogeneity among the nodes (figure 1b). Figures 4a and 4b, that shall be compared with figure 3, shows the spatial distribution of connectivities in the lattice. As expected, it is much more homogeneous. This homogeneous behavior is not destroyed if we vary the statistic of events.

Finally, our findings seems also to be robust with respect to the cell size. If we increase the size of the cell, defining, for instance, four adjacent sites of the lattice as a unique cell, there is no change in the results, not even in the exponent $`\gamma `$ that characterizes the degree distribution $`P(k)`$, as shown in figure 5.

In conclusion, we have shown that the non conservative version of the Olami-Feder-Christensen model is able to reproduce the scale free network associated to the dynamics of the epicenters observed on real earthquake data. The conservative version of the model displays a qualitatively different behavior, being more close to a random graph. The smallest degree of non conservation seems to be enough to change the behavior of the model, since for $`\alpha =0.249`$ we see that $`P(k)`$ has already a well defined power law behavior for some decades. Those findings, besides giving an indication of the robustness of this model to reproduce the dynamics of earthquakes, reproducing the experimental findings of Abe and Suzuki, present a completely new dynamical mechanism to generate a free scale network. There is no explicit rule of preferential attachment, and the preferential attachment observed in the network is a signature of the model dynamics. We hope that a complete study of the properties of the network can help to solve some still controversial aspects of the Olami-Feder-Christensen model and of self-organized critical behavior, and can be interesting and useful even if a more detailed study of earthquake data comes to show in the future that the the results reported by Abe and Suzuki are not universal.

Figure Captions

Figure 1: Normalized degree distribution $`P(k/I)`$ for different values of $`\alpha `$. $`I`$ is the total number of epecenters. (a) non conservative regime: the results show a free scale network behavior in all cases. The curves for $`\alpha <0.249`$ have been shifted upwards along the $`y`$ axis for clarity, otherwise they would all coincide. In all cases $`L=200`$ and the number of epicenters registered is $`10^5`$. (b) Conservative regime: the degree distribution is similar to a random graph; in this case we have $`L=200`$ and an statistics of $`10^6`$ events. Lowering the statistics does not change this behavior.

Figure 2: Exponent $`\gamma `$, that characterizes the power law behavior of $`P(k)`$, for different values of $`\alpha `$. $`\gamma `$ seems to increase linearly with $`\alpha `$. In all cases $`L=200`$ the number of epicenters is $`10^5`$.

Figure 3: Spatial distribution of node degrees in the non conservative case, for $`\alpha =0.249`$, $`L=200`$ and $`10^5`$ events. Sites associated with nodes of higher degree are darker and, as one can see, are closer to the boundaries. Figure (b) is a blow up of (a). The $`20`$ sites closer to the boundaries have not been not plotted and the scale has been changed in order to show the details. We can see that the structure of the network is reproduced and is not a boundary effect.

Figure 4: Spatial distribution of node degrees for the conservative case. Sites associated with nodes of higher degree are darker, $`L=200`$ and the number of epicenters is $`10^6`$. (a) The same scale of figure 3a has been used. (b) The scale has been changed to reveal details of the structure of the network that, in this case, is much more homogeneous and quite different than the one observed in the non conservative regime.

Figure 5: The normalized degree distribution $`P(k/I)`$ for $`\alpha =0.249`$, $`L=200`$ and $`10^5`$ events, for different cell sizes. $`I`$ is the total number of epecenters. (a) $`L=200`$ and each site of the lattice defines a cell. (b) $`L=400`$ and each four adjacent sites are in the same cell. The curve has been shifted upwards in the axis $`y`$ for clarity.


text:
# Polar Ring Galaxies and the Tully-Fisher relation: implications for the dark halo shape

## 1. Introduction

Polar Ring Galaxies (PRGs) are peculiar objects composed of a central spheroidal component, the host galaxy, surrounded by an outer ring, made up of gas, stars and dust, which orbits nearly perpendicular to the plane of the gas-poor central galaxy (PRGC; Whitmore et al. 1990). Previous papers (Arnaboldi et al. 1995, 1997; Iodice et al. 2002a, 2002b, 2002c) found that even where the morphology of the host galaxy resembles that of an early-type system, PRGs show many similarities with late-type galaxies. The PRGs are characterized by a large amount of neutral hydrogen (H i), always associated with the polar structure (Schechter et al. 1984; van Gorkom et al. 1987; Arnaboldi et al. 1997), and by a gas-to-total luminosity ratio in the B-band typical of late-type galaxies.

By exploring the properties of the host galaxy and ring in the optical and NIR, for a sample of PRGs, Iodice et al. (2002a, 2002b, 2002c) found that the connection with spirals is tighter. The Tully-Fisher relation (TF) is the most important scaling relation for disks (Tully & Fisher, 1977): this is an empirical relationship between the disk rotational velocity ($`V_{rot}`$) and its absolute luminosity ($`L`$), where $`LV_{rot}^4`$, approximately. In the past few years, several studies have asserted the validity of the TF relation for some classes of disk galaxies which show different photometric and kinematical properties with respect to ‘classical’, high-surface-brightness spiral galaxies (Matthews, van Driel & Gallagher 1998a, 1998b; McGaugh et al. 2000; Chung et al. 2002). These latest developments indicate that the TF relation is probing a very close liaison between the dark halo parameters and the total quantity of baryons in galaxies: the dark halo, which is responsible for the H i linewidth and the flat rotation curve in the outer regions of a disk, is tuned to the total amount of baryons in the luminous component.

In PRGs, the H i linewidth ($`\mathrm{\Delta }V`$) measures the dynamics along the meridian plane, which is dominated by the dark matter, while the baryons are nearly equally distributed between the host galaxy and the polar ring. We wish to investigate the position of the PRGs in the $`\mathrm{log}(\mathrm{\Delta }V)L`$ plane, and study via N-body simulations of 3-D systems whether the dark halo shape may influence their position in the $`\mathrm{log}(\mathrm{\Delta }V)L`$ plane, with respect to the TF relation of bright disks. The question of the dark halo shape is important i) to constrain dark matter models, through cosmological simulations (Navarro, Frenk & White 1996, 1997; Bullock et al. 2001) which predict the distribution of the halo shapes and the universal radial dependence of the dark matter distribution; ii) to give hints on the nature of dark matter (see Combes 2003 as a review); and furthermore iii) the dark halo properties in PRGs can give important constraints on the formation scenarios for these peculiar objects, which is still an open issue (see Iodice et al. 2002a, 2002c, and references therein).

## 2. Observations

New near-infrared J, H and Kn images are available for a sample of PRGs (Iodice et al. 2002a, 2002b, 2002c), which are selected from the PRGC, and for ESO 235-G58. The B band magnitude is known for many PRGs (Van Driel et al. 2002, 2003; Gallagher et al. 2002), and the H i integrated line profile data were obtained from several published hydrogen observations of PRGs, carried out by e.g. Richter et al. (1994), van Gorkom et al. (1987), van Driel and collaborators (2000, 2002, 2003), with several radio telescopes.

To derive the TF relation for normal disk galaxies, we will use the very large and detailed dataset available in the I-band from Giovanelli et al. (1997). We estimate the B band magnitude for spiral disks in the sample from Giovanelli et al. using the observational relation between morphological type index and the $`BI`$ colors (de Jong 1996).

### 2.1. PRGs and the TF Relation for Spiral Galaxies

In Fig. 1 we show the K-band TF relation for a sample of spiral galaxies studied by Verheijen (2001). The values of $`M_K`$ and $`\mathrm{log}(\mathrm{\Delta }V_{20})`$ for PRGs are also shown on this plot. Our and Verheijen’s data sets have very similar photometric properties and limiting magnitudes. We see that five PRGs lie near to the high-velocity boundary of the TF relation, or show larger velocities (for a given luminosity) than disk galaxies. Only the PRG AM2020-504 shows a lower rotation velocity for its K-band absolute magnitude.

In the B-band, the tendency of PRGs to have larger velocities with respect to the TF is confirmed for a larger sample of PRGs compared with data for 787 disk galaxies (Giovanelli et al. 1997). In Fig. 1, we see that two PRGs lie at lower velocities than those predicted by the TF for spiral galaxies, two objects lie on the TF relation, and twelve PRGs either lie on the high velocity boundary of the TF relation or show much larger velocities. Van Driel et al. (2003) also find that most kinematically confirmed PRGs show larger H i profile widths than bright spiral galaxies, at a given luminosity.

## 3. The PRG positions in the $`\mathrm{log}(\mathrm{\Delta }V)L`$ plane and their implications for the dark halo

Observing higher or lower velocities with respect to the linear TF of disk galaxies is relevant for the discussion on the dark halo shape. Via analytical models and simple assumptions about the mass distribution, either luminous or dark, we can estimate where the PRGs ought to lie in the $`\mathrm{log}(\mathrm{\Delta }V)M_{K/B}`$ plane (Fig. 1), with respect to the TF relation for disk galaxies. If there were no dark matter, and the gravitational potential were oblate in the same sense as the flattened host galaxy, the polar ring would acquire an eccentric shape. When the polar ring and the host galaxy are both seen edge-on, which is close to being the case for most of our PRGs, the net effect will be that the line-of-sight (LOS) polar ring velocities are reduced. In the logarithmic, scale-free, potential case, a simple formula gives the expected velocity ratio between the major and the minor axis components as:

$$ϵ_v=1\frac{v_{minor}}{v_{major}}=ϵ_\rho 2ϵ_\mathrm{\Phi }$$
(1)
from Gerhard & Vietri (1986), where $`ϵ_\rho `$ is the flattening (1-axis ratio) of the density distribution and $`ϵ_\varphi `$ is the potential flattening. Therefore we would expect PRGs to have on average lower velocities with respect to what would be measured in the equatorial plane.

This implies that when the polar structure is eccentric, the observed LOS velocities in $`\mathrm{\Delta }V_{20}`$ are the smallest, i.e. those from the particles in the polar regions, see Fig. 2. Thus the observed $`\mathrm{\Delta }V_{20}`$ depends on both the mean velocity along the ring, and the ring eccentricity. On the contrary, Fig. 1 shows that the majority of PRGs have larger velocities than expected in the $`\mathrm{log}(\mathrm{\Delta }V)M_{K/B}`$ plane. Therefore we need to investigate how these velocities can be produced, and how they may depend on the intrinsic properties of the dark galaxy halo. We computed a series of N-body models of the formation of PRGs (Bournaud & Combes 2003). Our 3-D and 2-D simulations of PRGs show different positions in the TF plane depending on the shape of their dark halos. When the halo is oblate and flattened towards the host galaxy, the observed velocity in the polar ring are then smaller, and PRGs lie on the left-hand of the TF relation for bright disks. When the dark halo is flattened towards the polar ring plane, the observed velocities as larger, shifting the PRGs to the right-hand side of the diagram as in Fig. 1 (see also Iodice et al. 2003).

## 4. Discussion

In the $`\mathrm{log}(\mathrm{\Delta }V)M_{K/B}`$ planes shown in Fig. 1, most PRGs have larger H i rotation velocities than standard spiral galaxies, at a given K or B-band luminosity of the stellar component. Our N-body simulations have suggested that a likely explanation for this effect is a flat dark halo, whose main plane is aligned with the host galaxy meridian plane, and prevents the polar ring to become eccentric. The question arises if other effects, i.e. non-homogeneities in the TF relations for spirals, caused by bar and/or non edge on disks, or larger M/L ratios, can produce similar results and therefore be alternative explanations for the high velocities observed in PRGs.

Can the offset between the TF relation for bright spirals and PRGs be caused by the PRGs being less luminous than spirals at a given velocity? No, it cannot: as shown in Fig. 2 the host galaxies in 5 PRGs fall on the TF relation for bright spirals, which indicates 1) that the M/L ratio of this component is different from those of standard S0s, and 2) a luminous-to-dark matter content ratio similar to those of standard bright disks. Gerhard et al. (2001) showed that elliptical galaxies follow a TF relation in the $`\mathrm{log}(\mathrm{\Delta }V)M_{}`$, where $`M_{}`$ is the total mass in the luminous component, which is shallower than the relation for spiral galaxies, even when the maximal $`M/L_B`$ is adopted to compute the total stellar masses. This led Gerhard and collaborators to infer that elliptical galaxies have slightly lower baryonic mass than spiral galaxies of the same circular velocities, and that their dark halos are denser than halos of spiral galaxy with the same $`L_B`$. How much more massive must the dark halo be, to account for the velocities observed in polar rings? The observed value of $`\mathrm{\Delta }V_{20}`$ depends largely on the dark halo shape and the ring eccentricity, while it varies only as the square root of the total mass (and depends even less on the dark mass). Thus, a large amount of dark matter is needed when the halo is not polar: it ranges from factor 2 for a spherical halo, to a factor 3.5/4 for an oblate/prolate halo. Such a massive halo would cause a large offset of the host galaxy from the TF relation of bright disks, which is not observed.

#### Acknowledgments.

M.A. wishes to thank the SOC for such a stimulating Symposium, and gratefully acknowledge the IAU for financial support.

## References

Arnaboldi, M., et al. 1995, Planet. Space Sci., 43, 1377

Arnaboldi, M., et al. 1997, AJ, 113, 585

Bekki, K. 1998, ApJ, 499, 635

Bournaud, F. & Combes, F. 2003, A&A, 401, 817

Bullock, J. S., et al. 2001, MNRAS, 321, 559

Chung, A., et al. 2002, AJ, 123, 2387

Combes, F. 2003, NewA Reviews, 46, 755

de Jong, R. S. 1996, A&A, 313, 377

Gallagher, J. S., et al. 2002, ApJ, 568, 199

Giovanelli, R., et al. 1997, AJ, 113, 22

Gerhard, O. & Vietri, M. 1986, MNRAS, 223, 377

Gerhard, O., et al. 2001, AJ, 121, 1936

Iodice, E., et al. 2002a, AJ, 123, 195

Iodice, E., et al. 2002b, A&A, 391, 103

Iodice, E., et al. 2002c, A&A, 391, 117

Iodice, E., et al. 2003, ApJ, 585, 730

Matthews, L.D., van Driel, W. & Gallagher, J.S. 1998a, AJ, 116, 1169

Matthews, L.D., van Driel, W. & Gallagher, J.S. 1998b, AJ, 116, 2196

Mathieu, A., Merrifield, M.R. & Kuijken, K. 2002, MNRAS, 330, 251

McGaugh, S. S., et al. 2000, ApJ, 533, 99

Navarro, J.F., Frenk, C.S. & White, S.D. 1996, ApJ, 462, 563

Navarro, J.F., Frenk, C.S. & White, S.D. 1997, ApJ, 490, 493

Reshetnikov, V. & Sotnikova, N. 1997, A&A, 325, 933

Richter, O.-G., Sackett, P. D. & Sparke, L. S. 1994, AJ, 107, 99

Schechter, P.L., et al. 1984, MNRAS, 210, 497

Tully, R. B. & Fisher, J. R. 1977, A&A, 54, 661

van Driel, W., et al. 2000, A&A, 141, 385

van Driel, W., et al. 2002, A&A, 386, 140

van Driel, W., et al. 2003, A&A, in prep

van Gorkom, J. H., Schechter, P. L. & Kristian, J. 1987, ApJ, 314, 457

Verheijen, M. A. 2001, ApJ, 563, 694

Whitmore, B. C., Lucas, R. A. & McElroy, D. B. 1990, AJ, 100, 1489 (PRGC)