format: text

text:
# Accretion disc–stellar magnetosphere interaction: field line inflation and the effect on the spin-down torque.

## 1 Introduction

Accreting magnetised stars exist in the form of neutron stars in close binary systems and T Tauri stars with surrounding protostellar discs. In both cases an accretion disc is disrupted through interaction with the stellar magnetosphere and the accretion flow becomes channelled along field lines. Torques act between the star and the accreting material and they determine the spin history of the star and its equilibrium rotation period. A spin–up torque is produced by rapidly rotating material in the inner regions of the disc while a spin–down torque is produced by the more slowly rotating material in the outer parts. Calculation of the torques requires knowledge of the structure of the disc and the magnetosphere in which it is embedded.

Early work (\[Ghosh & Lamb (1979a)\]; \[Ghosh & Lamb (1979b)\]) suggested a direct dependence of the total torque on the mass accretion rate. However, recent observations of accreting neutron stars (\[Nelson et al. (1997)\]; \[Chakrabatry et al. (1997)\]) suggest that the torques can oscillate producing alternating spin–up and spin–down phases with no evidence of a consistent correlation between the torque and the mass accretion rate. ?) suggest that this may be due to a variable outer magnetosphere which causes a corresponding variation in the spin–down torque acting between the star and the disc. In this paper we investigate the structure of a magnetosphere corotating with the central star. We calculate the spin–down torque due to the star–disc interaction and we find that it is very sensitive to the amount of field line twisting that occurs through the coupling of the stellar field lines to the disc. The magnetospheric structure and the resulting torque are also sensitive to the applied boundary conditions. In the present work we neglect any torque that may arise from a stellar wind.

According to the Ghosh & Lamb model the magnetic field of the central star penetrates the accretion disc due to diffusion arising from the growth of various instabilities. In the outer parts of the disc, viscous stresses are more important than magnetic stresses so that the effects of the magnetic field on the disc structure can be ignored. In the inner parts, the flow becomes dominated by magnetic stresses which eventually lead to the truncation of the disc at an inner radius, $`R_{\mathrm{d},}`$ and the channeling of the flow to the star along stellar magnetic field lines. The twisting of the poloidal field by the vertical shear is assumed to be counteracted by the reconnection of the oppositely directed toroidal fields through the disc vertical extent. This occurs on a rapid timescale resulting to an equilibrium ratio for toroidal, $`B_\phi `$, over vertical, $`B_z`$ field components at the disc surface. Alternatively the growth of $`B_\phi `$ can be controlled by turbulent diffusion through the disc (\[Campbell (1987)\]; \[Yi (1994)\]) or reconnection of field lines in the magnetosphere (i.e. \[Livio & Pringle (1992)\]). In all these cases magnetic torques of a similar form are produced \[Wang (1995)\].

In all calculations of the magnetic torques resulting from disc-star interactions $`B_z`$ is often assumed to be approximately equal to the stellar dipole field. However, twisting magnetic field lines imply the existence of currents which will be distributed in the magnetically dominated magnetosphere which is assumed to corotate with the central star. The star-disc magnetosphere is expected to attain a force–free equilibrium in an Alfvén wave crossing time, as is the case in the solar corona. In the force–free equilibrium, the poloidal magnetic field will deviate from that of a dipole, the effect being especially significant in the case of large twisting (or $`|B_\phi /B_z|`$). We find that the magnetic torque acting on the disc may also be significantly reduced.

The related problem of the evolution of an initially potential field in response to an increasing footpoint shear through a sequence of quasi–static force–free equilibria has been addressed in the solar context. But note that, unlike in the work presented here, ideal MHD is assumed. Numerical calculations supported theoretical expectations (\[Aly (1984)\]; \[Aly (1985)\]; \[Aly (1991)\]; \[Sturrock (1991)\]) that the energy of the force–free field increases monotonically with increasing shear, approaching the energy of an open field (\[Klimchuck (1991)\]; \[Roumeliotis et al. (1994)\]; \[Mikic & Linker (1994)\]). Analytical sequences of quasi–static equilibria were also constructed using a self–similar approach (\[Newman et al. (1992)\]; \[Lynden-Bell & Boily (1994)\]; \[Wolfson (1995)\]) and it was shown that the complete opening of the initially closed field lines is possible at a finite shear. This is associated with the development of current sheets as was conjectured by ?).

?) have calculated magnetospheric force–free equilibria resulting from the twisting of an initially dipolar field anchored on a star which threads an embedded differentially rotating disc. Disc resistivity was taken to be of turbulent origin. However, because of problems associated with their numerical method these authors were not able to calculate equilibrium configurations for $`|B_\phi /B_z|`$ larger than $`2.55`$ at the surface of the disc. In most cases a maximum $`|B_\phi /B_z|1`$ was adopted. Nonetheless significant field line inflation was observed.

In the work presented here, we similarly find that moderate to large shearing of the stellar dipole field through interaction with the disc can lead to a partially open configuration, with expulsion of toroidal field, for favourable outer boundary conditions. Although the process can be inhibited by the presence of an outer boundary that confines the field, significant field line inflation may occur with a large reduction of the magnitude of the vertical field in comparison to the initial dipole value. This could have important consequences for studies of stellar spin evolution which have assumed a dipolar poloidal field at equilibrium (\[Ghosh & Lamb (1979b)\]; \[Cameron & Campbell (1993)\]; \[Yi (1995)\]; \[Ghosh (1995)\]; \[Armitage & Clarke (1996)\]).

After formulating the physical model and basic equations in Section 2, we describe our numerical method in Section 3. Using this we were able to perform calculations with a maximum value of $`|B_\phi /B_z|`$ at the disc surface of $`120`$ and with a variety of external disc radii. Results for Keplerian discs with different amounts of field twisting are given in sections 4. These calculations were performed with a field–confining, perfectly conducting outer boundary.

We also performed calculations with a partly sub–Keplerian accretion disc using the form of angular velocity, $`\mathrm{\Omega }`$ given by ?). This emulates the existence of an inner boundary layer as proposed in the Ghosh & Lamb model. Results are given in Section 4.3 for the conducting outer boundary condition. In Section 4.4 we investigate the effect of the outer boundary condition, showing that a boundary condition that allows field lines to penetrate it, leads to more open configurations.

In Section 4.5 we go on to evaluate the spin down torque acting between star and disc for the force free configurations. We show that this torque can be reduced by up to a factor of $`100,`$ in comparison to what would be obtained assuming an unmodified dipole field, in a configuration that has undergone large twisting and toroidal field expulsion.

Finally in Section 5 we summarize and discuss our results and their application.

## 2 Basic Equations

### 2.1 Force–free equilibria

We consider a magnetic field $`𝑩_{\mathrm{dp}}`$ anchored on a star that rotates with a rate $`\mathrm{\Omega }_{}.`$ We suppose that the magnetic field is axisymmetric with symmetry axis aligned with the rotation axis. A thin accretion disc orbits in the equatorial plane such that the system remains axisymmetric. We ignore the internal dynamics of the disc and assume that the material rotates with a prescribed angular velocity $`\mathrm{\Omega }.`$ For the most part this will equal or be close to the local Keplerian angular velocity $`\mathrm{\Omega }_\mathrm{K}\sqrt{GM/R^3}`$ where $`R`$ is the radial distance to the star measured in the midplane of the disc. We suppose that there is a highly conducting low density corona that corotates with the star and in which the disc is immersed.

However, the disc has non zero resistivity such that the stellar field is enabled to diffuse into the disc. Because the star/corona and disc rotate at different rates, field lines that permeate the disc are twisted in the azimuthal direction such that a toroidal component of the magnetic field is generated.

The disc is considered to be turbulent and manifesting an effective magnetic diffusivity, $`\eta .`$ This enables a steady state to be achieved in which the generation of toroidal magnetic field through differential rotation between the disc and star/corona is balanced by the effects of turbulent diffusion.

The corona responds to the twisting of the initially dipolar field by producing toroidal magnetic field and current density components. These adjust such that a force–free equilibrium is eventually attained. We expect this situation to occur when the coronal Alfvén speed sufficiently exceeds the characteristic rotational speeds. This requires a low density corona which is magnetically dominated. For the purposes of the work presented here, we assume that such a corona exists in which a force–free equilibrium is established in a time short compared to the stellar rotation period. Thus we look for steady state axisymmetric configurations involving the star, the disc and the corona.

An axisymmetric magnetic field $`𝑩`$ can be split into poloidal $`𝑩_\mathrm{p}`$ and toroidal $`B_\phi \widehat{\phi }`$ components such that:

$$𝑩_\mathrm{p}=\mathbf{}\mathbf{\times }(\mathrm{\Psi }\widehat{\phi }/(r\mathrm{sin}\theta )),$$
(1)
where $`\mathrm{\Psi }`$ is the magnetic stream function which is constant on field lines. Here we use spherical polar coordinates $`(r,\theta ,\phi )`$ based on the central star. Then

$$𝑩_\mathrm{p}=(\frac{1}{r^2\mathrm{sin}\theta }\frac{\mathrm{\Psi }}{\theta },\frac{1}{r\mathrm{sin}\theta }\frac{\mathrm{\Psi }}{r})$$
(2)
The force–free condition implies that the Lorentz force is zero everywhere, thus:

$$𝑱\mathbf{\times }𝑩=0$$
(3)
where $`𝑱=(c/4)\times 𝑩`$ is the current density (c.g.s units are used throughout). The toroidal component of equation (3) gives:

$$𝑩_\mathrm{p}\mathbf{}(r\mathrm{sin}\theta B_\phi )=0.$$
(4)
The definition of $`\mathrm{\Psi }`$ implies constancy on field lines or equivalently $`𝑩_\mathrm{p}\mathbf{}\mathrm{\Psi }=0`$. Thus (4) implies that:

$$r\mathrm{sin}\theta B_\phi =f(\mathrm{\Psi })f,$$
(5)
where $`f`$ is an arbitrary function of $`\mathrm{\Psi }.`$

Using (5) and the poloidal component of (3) we finally obtain the governing equation for $`\mathrm{\Psi }`$ in the form:

$$\mathrm{\Delta }\mathrm{\Psi }=ff^{}$$
(6)
where $`f^{}\mathrm{d}f/\mathrm{d}\mathrm{\Psi }`$ and $`\mathrm{\Delta }`$ is the differential operator defined through:

$$\mathrm{\Delta }=\frac{^2}{r^2}+\frac{1\mu ^2}{r^2}\frac{^2}{\mu ^2},$$
(7)
where $`\mu \mathrm{cos}\theta `$.

Equation (6) can be solved for $`\mathrm{\Psi }`$ if $`f`$ is known. Since $`f(\mathrm{\Psi })`$ is constant on poloidal field lines we need to calculate its value at only one point along each field line. This we do by consideration of the interaction of the field with the resistive disc.

### 2.2 Disc–stellar field interaction

The corona above the disc is assumed to corotate with the star and so a toroidal field will be generated due to the vertical shear. Neglecting the radial component of the field interior to the disc, we calculate the internal toroidal field using the toroidal component of the induction equation which in cylindrical coordinates $`(R,\phi ,z)`$ takes the form:

$$\frac{B_\phi }{t}=R𝑩_\mathrm{p}\mathrm{\Omega }+\frac{\eta }{R}\mathrm{\Delta }_\mathrm{c}(RB_\phi )+\frac{1}{R}\eta (RB_\phi )$$
(8)
where $`\mathrm{\Delta }_\mathrm{c}=^2(2/R)(/R)`$ and $`\eta =𝒟\nu `$ is the turbulent magnetic diffusivity which is assumed to be proportional to the turbulent viscosity $`\nu ,`$ the constant of proportionality being $`𝒟.`$ We use the Shakura & Sunyaev (1973) parameterization for $`\nu ,`$ $`\nu =\alpha _{\mathrm{SS}}\mathrm{\Omega }_\mathrm{K}H^2.`$ Here the viscosity parameter is $`\alpha _{\mathrm{SS}}`$ and $`H`$ is the disc semithickness. The gas velocity in the disc is taken to be $`𝒖=(0,r\mathrm{\Omega },0).`$

We consider that $`B_\phi `$ is generated by the shearing of the vertical magnetic field component arising from the $`z`$ dependence of $`\mathrm{\Omega }.`$ For simplicity we take $`\eta `$ to be a function of radial distance only. For a more general treatment see ?). Noting that the disc is thin we expect $`/z(R/H)/R`$ and so we neglect radial derivatives in equation (8). Then in a steady state we have:

$$\eta \frac{{}_{}{}^{2}B_{\phi }^{}}{z^2}=RB_z\frac{\mathrm{\Omega }}{z}$$
(9)
To find $`B_\phi `$ we need to specify $`\mathrm{\Omega }/z`$. We follow previous authors (\[Wang (1987)\]; \[Campbell (1987)\]; \[Yi (1994)\]) and take the vertical shear to be concentrated in a thin layer close to the disc surface. There the angular velocity changes between the disc midplane value and $`\mathrm{\Omega }_{}.`$ We take the vertical average of equation (9) the right hand of which gives:

$$\frac{RB_z}{H}_0^H\frac{\mathrm{\Omega }}{z}𝑑z=\frac{RB_z(\mathrm{\Omega }_{}\mathrm{\Omega })}{H},$$
(10)
where $`B_z`$ is assumed not to vary with $`z`$ thus being the disc midplane value. In the following $`\mathrm{\Omega }`$ will refer to the disc midplane value.

We assume a simple approximation for the vertical average of the dissipative term on the left-hand-side of equation (9) namely:

$$\frac{1}{H}_0^H\eta \frac{{}_{}{}^{2}B_{\phi }^{}}{z^2}𝑑z\eta \gamma \frac{B_\phi }{H^2},$$
(11)
where $`B_\phi `$ now applies to the value at the upper disc surface and $`\gamma `$ is a dimensionless parameter expected to be of order unity.

Combining equation (10) with equation (11) we obtain the following estimate for the equilibrium value of the ratio $`B_\phi /B_z`$ at the upper disc surface:

$$\frac{B_\phi }{B_z}=\pm \frac{RH}{\gamma \eta }(\mathrm{\Omega }_{}\mathrm{\Omega })$$
(12)
We have inserted the $`\pm `$ alternative to allow for different senses of rotation of the star and disc with respect to the right handed coordinate system with the negative sign corresponding to clockwise rotation. The angular velocities are then always positive. Using the turbulent diffusivity specified above and taking $`\gamma =1,`$ equation (12) gives for $`\mathrm{\Omega }=\mathrm{\Omega }_K`$:

$$\frac{B_\phi }{B_z}=\pm \frac{R}{𝒟\alpha _{\mathrm{SS}}H}\left[\left(\frac{R}{R_\mathrm{c}}\right)^{3/2}1\right]$$
(13)
From equation (13) we see that $`B_\phi `$ changes sign at the corotation radius, $`R_\mathrm{c}(GM/\mathrm{\Omega }_{}^2)^{1/3}`$. Similarly the vertically integrated $`z\phi `$ component of the magnetic torque per unit area, $`B_zB_\phi R^2/(2)`$, also changes sign such that the star transfers angular momentum to the disc outside the corotation radius while it gains angular momentum from disc material inside that radius.

In our calculations we shall assume that $`RH\mathrm{\Omega }_{}/(\gamma \eta )C=\text{constant}`$. Equation (13) then becomes (adopting the negative sign corresponding to clockwise rotation):

$$\frac{B_\phi }{B_z}=C\left[\left(\frac{R_\mathrm{c}}{R}\right)^{3/2}1\right].$$
(14)
Our assumption of $`C=\text{constant}`$ is such that for large values of $`R/R_\mathrm{c},`$ $`B_\phi /B_zC`$ which is constant. However, an important result of the calculations, namely the inflation of poloidal field lines such that $`B_z0`$ at the disc surface for large $`C`$, should not depend on having a dependence of $`C`$ on radius provided it remains large. For $`H/R=0.1`$, $`\alpha _{\mathrm{SS}}=0.01`$ and $`𝒟=\gamma =1`$ at the corotation radius, these being reasonable values for accretion discs, equation (14) gives $`|B_\phi /B_z|=C=10^3`$ for $`RR_\mathrm{c}.`$ We comment that equation (13) indicates that in a disc with constant $`H/R`$ and constant $`\alpha _{SS},`$

$$C\frac{1}{𝒟}\left(\frac{R}{R_\mathrm{c}}\right)^{3/2}.$$
(15)
Thus in this case if $`𝒟=1`$, $`C`$ would indeed increase with radius.

Our analysis is distinct from that of e.g. ?) who use a similar expression for $`|B_\phi /B_z|`$ to equation (14) (for $`RR_\mathrm{c}`$) but who assume that fast reconnection in the star–disc corona with the coronal Alfvén speed restricts $`|B_\phi /B_z|`$ to a maximum value of unity. The above discussion suggests that $`C`$ and hence $`|B_\phi /B_z|`$ at the disc surface is large. The actual magnitude of $`|B_\phi |`$ is then controlled by the requirement of a force free equilibrium in the corona. This results in a significant departure of the poloidal field from the original stellar dipole through a significant inflation and opening out of the field lines.

Equation (6) is solved for $`R_\mathrm{d}<r<R_{\mathrm{ext}}`$ where $`R_\mathrm{d}`$ and $`R_{\mathrm{ext}}`$ correspond to the inner and outer disc radii. $`R_\mathrm{d}`$ is usually taken to be small enough that the magnetic torque that is applied there overwelhms the viscous torque so that angular momentum transport is magnetically regulated.

### 2.3 The location of the inner disc radius

An accurate calculation of the inner disc radius would require the solution of the accretion disc structure problem with the inclusion of all magnetic stresses and allowing for sub–Keplerian rotation in the innermost disc region. In most studies an approximate value for $`R_\mathrm{d}`$ has been calculated based on the criterion of disc disruption at a radius where magnetic stresses start to dominate viscous stresses. In order for stellar accretion to proceed the magnetic field must not disrupt the disc exterior to the corotation radius since the magnetic torque for $`R>R_\mathrm{c}`$ imparts angular momentum to the disc gas which cannot connect to the stellar field. Assuming accretion takes place so that $`R_\mathrm{d}<R_\mathrm{c},`$ the inner disc radius can be estimated by requiring that the magnetic torque balances the rate of angular momentum advection by the flow or that:

$$\dot{M}\left[\frac{d\left(\mathrm{\Omega }R^2\right)}{dR}\right]_{R=R_\mathrm{d}}=[B_\phi B_z]_{r=R_\mathrm{d}}R_\mathrm{d}^2$$
(16)
where $`\dot{M}`$ is the mass accretion rate and the positive sign corresponds to clockwise rotation. Because of the rapid increase of the magnetic field strength inwards estimates of $`R_\mathrm{d}`$ using equation (16) (e.g. \[Wang (1995)\]; \[Yi (1995)\]) give $`R_\mathrm{d}`$ close to $`R_c`$ when the star is near its equilibrium spin rate.

However, we point out that this analysis is based on an assumed stellar dipole structure for $`B_z`$. When the effect of an increasing $`B_\phi `$ at the disc surface is taken into account we expect that both $`B_z`$ and the magnetic torques will be modified. We will return to this point below.

In our numerical calculations we have mostly set $`R_\mathrm{d}=R_\mathrm{c}`$ implicitly assuming that the star is close to its equilibrium spin rate. For completeness we also performed some calculations with a sub–Keplerian inner disc region (see Section 4.3).

## 3 Numerical solution

### 3.1 Previous work

Equation (6) has been solved by ?) using a relaxation method with a coordinate transformation $`r1/r`$ so that the solution can extend to arbitrary large radii in principle. ?) adopted equation (13) for $`B_\phi /B_z`$ with a maximum expected value of $`C3.2510^3(R_\mathrm{d}/R)^{1/8}.`$ However, in practice the numerical method would not converge unless $`C`$ was multiplied by a small parameter $`\lambda `$ whose value depends on $`R_\mathrm{c}/R_\mathrm{d}.`$ As a result the maximum value of $`|B_\phi /B_z|`$ at the disc surface that could be obtained in a converged calculation was $`2.55.`$ For this particular case $`\lambda =1.610^5`$ (A.Bardou, private communication). ?) found that the force–free equilibria consist of closed field lines whose footpoints are shifted radially outwards with respect to their position in the dipolar state and with field lines that are flattened at small $`\theta `$ as compared to dipolar ones. ?) found the same indication from a similarity solution. Numerical convergence problems were encountered also by ?) who used a relaxation method to solve equation (6) with a prescribed $`f`$ in the context of the solar corona. They found that no equilibrium solution could be calculated for $`|B_\phi /B_z|\begin{array}{c}>\\ \end{array}1`$. We comment that methods for solving equation (6) based on a simple iteration method for producing a contraction mapping (see also \[Wolfson & Verma (1991)\]) using successive previous approximations for the right hand side tend to fail to converge if $`ff^{}`$ is a sufficiently rapidly varying function of $`\mathrm{\Psi }.`$ By comparison with solutions of related algebraic problems, the lack of convergence in solving equation (6) is not necessarily related to the emergence of additional solutions through a bifurcation. On the other hand studies that employ different methods of solution do not seem to encounter convergence problems (i.e. \[Klimchuck (1991)\]; \[Roumeliotis et al. (1994)\]).

We conclude that previous studies of force–free equilibria related to star–disc interactions were restricted to modest values of $`|B_\phi /B_z|`$ at the disc surface due to restrictions enforced by the adopted numerical method.

### 3.2 Numerical method

In this paper we adopt a different method of solution of equation (6) to those indicated above. In our case equation (6) can be solved without the use of a controlling parameter $`\lambda `$ and in principle for any disc surface value of $`|B_\phi /B_z|.`$ In practice some restrictions also apply in this case but we were able to perform calculations for discs with a ratio of inner to outer radius of several decades and for a maximum value of $`|B_\phi /B_z|10^2`$.

In our method equation (6) is transformed into a parabolic PDE by moving the source term to the left of equation (6) and by equating the new expression to a multiple of the time derivative of $`\mathrm{\Psi }`$ which is expected to tend to zero close to an equilibrium configuration. We use an explicit finite–difference scheme to integrate the parabolic PDE forward in time as an initial value problem until a steady state is achieved. Following this approach rather than solving the elliptic equation (6) directly, we were able to obtain equilibria corresponding to a wide range of values of the source term.

The total magnetic stream function, hereafter denoted by $`\mathrm{\Psi }_\mathrm{t}`$, can be written as a sum of contributions arising from disc currents, $`\mathrm{\Psi }_\mathrm{d},`$ and the central dipole, $`\mathrm{\Psi }_{\mathrm{dp}}.`$ Thus $`\mathrm{\Psi }_\mathrm{t}=\mathrm{\Psi }_{\mathrm{dp}}+\mathrm{\Psi }_\mathrm{d}`$, where $`\mathrm{\Psi }_{\mathrm{dp}}`$ is given by:

$$\mathrm{\Psi }_{\mathrm{dp}}=K\frac{1\mu ^2}{r},$$
(17)
and $`K`$ is a normalizing constant. Note that $`\mathrm{\Delta }\mathrm{\Psi }_{\mathrm{dp}}=0`$ since $`f=0`$ for a dipole field. Equation (6) with $`f0`$ is then solved for $`\mathrm{\Psi }_\mathrm{d}`$ only from:

$$\mathrm{\Delta }\mathrm{\Psi }_\mathrm{d}=f(\mathrm{\Psi }_\mathrm{t})f^{}(\mathrm{\Psi }_\mathrm{t}).$$
(18)
We find it convenient to drop the subscript “$`\mathrm{d}`$” from the disc magnetic stream function and use dimensionless quantities, $`\overline{r}=r/r_{}`$, $`\overline{\mathrm{\Psi }}=\mathrm{\Psi }/\mathrm{\Psi }_0`$ and $`\overline{B}=B/B_0`$ where $`B`$ represents the magnitude of either the toroidal or the poloidal magnetic field components and we take $`\mathrm{\Psi }_0=B_0r_{}^2`$ such that $`B_0=B_{\mathrm{dp}}(r_{},0)`$, the magnitude of the dipole field on the equatorial plane. Equations (6) and (18) can then be read in dimensionless form if all quantities are replaced by their barred equivalents. For simplicity we will drop the bars from the dimensionless quantities from now on.

The computational domain is bounded by an inner radius $`r_{\mathrm{in}}`$, an outer radius $`r_{\mathrm{out}}=R_{\mathrm{ext}}`$, the disc midplane at $`\mu =0`$ and the symmetry axis, $`\mu =1.`$ The grid spacing is taken to be uniform in both $`r`$ and $`\mu .`$

At $`r=r_{\mathrm{in}}`$ and $`\mu =1`$ respectively we specify the boundary condition $`\mathrm{\Psi }=0`$. The first condition specifies only the dipole flux exists at the inner boundary while the second is a requirement of the coordinate system. On the disc midplane, $`\mu =0`$, we impose the symmetry condition, $`\mathrm{\Psi }/\mu =0`$.

We have considered two different outer boundary conditions. The first one is the Dirichlet condition: $`\mathrm{\Psi }=0`$ which forces the field due to the disc to be tangential there. Physically this corresponds to a conducting boundary which excludes the field arising from currents in the disc but not the original dipole field which may be assumed to have had infinite time to diffuse. The second outer boundary condition that we use is the Neumman condition: $`\frac{\mathrm{\Psi }}{r}=0.`$ This makes the disc field radial at the outer boundary as might be the case if there were a coronal wind there. The second boundary condition leads to much more open configurations than the first and hence to larger departures from the stellar dipole field.

We add a term $`F(\mathrm{\Psi }/\tau )`$ to the right-hand side of the dimensionless form of equation (18) where we are allowed to choose $`F`$ as an arbitrary function of position and solve the equation with a forward in time, $`\tau `$, and centered in space (FTCS) finite–difference scheme on a $`(Nr\times Nm)`$, $`(r,\mu )`$ computational grid. The finite–difference equation that we use is given at a grid point denoted with subscript $`(i,j)`$ by:

$`F_{i,j}{\displaystyle \frac{\mathrm{\Psi }_{i,j}^{n+1}\mathrm{\Psi }_{i,j}}{\mathrm{d}\tau }}`$ $`=`$ $`{\displaystyle \frac{\mathrm{\Psi }_{i+1,j}2\mathrm{\Psi }_{i,j}+\mathrm{\Psi }_{i1,j}}{(\mathrm{d}x)^2}}`$ (19)
$`+`$ $`{\displaystyle \frac{1\mu _j^2}{r_i^2}}{\displaystyle \frac{\mathrm{\Psi }_{i,j+1}2\mathrm{\Psi }_{i,j}+\mathrm{\Psi }_{i,j1}}{(\mathrm{d}\mu )^2}}+f\left(\mathrm{\Psi }_\mathrm{t}|_{i,j}\right)f^{}\left(\mathrm{\Psi }_\mathrm{t}|_{i,j}\right).`$
where $`\mathrm{d}r,`$ and $`\mathrm{d}\mu ,`$ are the uniform grid spacings in $`r,`$ and $`\mu ,`$ respectively. Variables without superscripts are at $`\tau `$ level $`n.`$ The time step is $`\mathrm{d}\tau .`$ We assume that the disc is truncated at $`R_\mathrm{d}`$ by an infinite conductor so we have $`B_\phi =0`$ for $`1<RR_\mathrm{d}.`$ For all other values of $`R`$ on the equator ($`\mu =0`$ denoted by subscript $`1`$) we use equation (12) approximated as follows:

$$B_\phi |_{i+1/2,1}=CB_z|_{i+1/2,1}\frac{\mathrm{\Omega }_{i+1/2}\mathrm{\Omega }_{}}{\mathrm{\Omega }_{}}.$$
(20)
In equation (20) $`\mathrm{\Omega }_{i+1/2}`$ is either given by the Keplerian rate or by the form proposed by ?) (see Section 4.3). $`\mathrm{\Omega }_{}`$ is calculated on the $`(i+1/2,1)`$ grid point that is closest to the prescribed value of $`R_\mathrm{c}`$. The vertical magnetic field at the disc surface is calculated at the $`(i+1/2,1)`$ points by numerical differentiation of equation (2).

Using equation (20) we calculate $`f_{i+1/2,1}=R_{i+1/2}B_\phi |_{i+1/2,1}`$ while $`(ff^{})_{i,1}`$ is calculated from:

$$(ff^{})_{i,1}=\frac{1}{2}\frac{f_{i+1/2,1}^2f_{i1/2,1}^2}{\mathrm{\Psi }_\mathrm{t}|_{i+1/2,1}\mathrm{\Psi }_\mathrm{t}|_{i1/2,1}}.$$
(21)
Since $`f`$ is a function of $`\mathrm{\Psi }_\mathrm{t}`$ only we construct a table of values of $`f`$ as function of $`\mathrm{\Psi }_\mathrm{t}|_{i,1}`$ on the equator. We then use this table to calculate $`f`$ for any required value of $`\mathrm{\Psi }_\mathrm{t}`$ in equation (19) using linear interpolation.

The main disadvantage of the method we use is the constraint on the time step coming from the requirements of numerical stability. We find that the maximum allowed timestep decreases rather steeply with $`C.`$

In order to achieve a more rapid evolution for the same computational time and always keeping in mind that we are only interested in final steady state equilibria, we have introduced a spatially varying “diffusion coefficient”, $`(1/F)`$, that multiplies the right-hand side of equation (19). This coefficient is taken to be equal to $`\mathrm{d}\tau _{i,j}/\mathrm{min}(\mathrm{d}\tau _{i,j})`$ and therefore allows for a more advanced evolution of the values at the grid points where $`\mathrm{d}\tau _{i,j},`$ the maximum allowed timestep based on local stability considerations, can be larger. Since the source term peaks near the corotation radius but it is significantly smaller elsewhere it is possible to reach the equilibrium in significantly smaller computational time in this way.

Time integrations have been performed until the right hand side of equation (19) becomes smaller than $`10^6`$ in dimensionless units although the magnetic field configuration is usually found to have converged to its final form well before this criterion is satisfied.

## 4 Numerical Results

The parameters used in the calculations presented here are summarised in Table 2. All the calculations were performed with $`r_{\mathrm{in}}=1`$. Different values were used for $`R_{\mathrm{ext}}`$ and $`C`$ and calculations were done using both outer boundary conditions on $`\mathrm{\Psi }.`$ All the calculations with $`\mathrm{\Omega }=\mathrm{\Omega }_\mathrm{K}`$ were performed with $`R_\mathrm{d}=R_\mathrm{c}`$. Our choice of $`R_\mathrm{c}=3`$ implies $`\mathrm{\Omega }_{}=0.19\mathrm{\Omega }_\mathrm{K}(R=R_{})`$ which corresponds to a rotation period of $`5\mathrm{days}`$ for $`R_{}=3R_{}`$ and $`M_{}=0.5M_{}`$. The observed rotation periods of CTTS are $`3`$$`12`$ days (\[Bouvier et al. (1997)\] and references therein).

The choice of $`R_{\mathrm{ext}}`$ was dictated mainly by the limitations of the computational method since a much larger value of $`R_{\mathrm{ext}}`$ would require increased values of $`Nr`$ and $`Nm`$ for the same resolution to be achieved as in the cases with smaller outer disc radii.

We chose $`C=3`$$`120`$ in order to examine the effect of an extended range of values of $`|B_\phi /B_z|`$ at the disc surface on the magnetospheric equilibrium. Various authors have assumed $`|B_\phi /B_z|1`$ at the disc surface in their analysis. ?) argue that anomalous resistivity in the disc would limit $`B_\phi /B_z.`$ ?) and ?) assert that fast reconnection outside the disc would lead to a similar result. These physical arguments have not so far been supplemented by detailed calculations in the context of star–disc interactions. Our maximum value of $`C`$ was chosen so that the computational times required are reasonable. Nevertheless much larger values of $`C`$ and therefore of $`B_\phi `$ on the disc surface could result in very large toroidal magnetic pressure on the disc surface and the subsequent destabilisation of the disc.

In Table 2 we have also included the parameters of calculations performed with a sub–Keplerian form of $`\mathrm{\Omega }`$. This is characterised by the “fastness” parameter $`\omega \mathrm{\Omega }_{}/\mathrm{\Omega }_\mathrm{K}(R_\mathrm{d})`$ (for more details see Section 4.3).

### 4.1 Qualitative considerations of star–disc interactions

The main result of the production of toroidal field in the star–disc corona is the presence of a feedback mechanism that reduces $`|B_z|`$ at the disc surface. In this way the magnitude of $`B_\phi `$ does not become too large to maintain a force free equilibrium when $`C`$ is very large. As $`C`$ is increased, the field becomes increasingly twisted. The result is the production of azimuthal currents as well as a toroidal field component in the corona. The azimuthal currents tend to reduce $`|B_z|`$ near the disc and hence the magnitude of $`B_\phi .`$ To see this in a qualitative way: the disc flux obeys equation (18) which, for the boundary conditions used, can be converted to an integral form:

$$\mathrm{\Psi }=_VG(r,\mu ,r^{},\mu ^{})f(\mathrm{\Psi }_\mathrm{t})f^{}(\mathrm{\Psi }_\mathrm{t})r^2𝑑\mu ^{}𝑑r^{}.$$
(22)
Here $`G`$ is an appropriate Green’s function and the stream function inside the integral, taken over the computational domain, is evaluated using the primed coordinates. Thus:

$$\mathrm{\Psi }_\mathrm{t}=_VG(r,\mu ,r^{},\mu ^{})f(\mathrm{\Psi }_\mathrm{t})f^{}(\mathrm{\Psi }_\mathrm{t})r^2𝑑\mu ^{}𝑑r^{}+\mathrm{\Psi }_{\mathrm{dp}}.$$
(23)
For small $`C,`$ $`f`$ is small and $`\mathrm{\Psi }_\mathrm{t}`$ is close to the dipole value $`\mathrm{\Psi }_{\mathrm{dp}}.`$ However, for large $`C`$ and hence $`f,`$ a solution can only be found (or the force free equilibrium can only be maintained) by otherwise reducing the magnitude of the integral. Note that if $`\mathrm{\Psi }_{\mathrm{dp}}`$ is scaled by multiplying by a constant value, both $`\mathrm{\Psi }_\mathrm{t}`$ and the integral above scale similarly. The magnitude of the integral can be reduced by having fewer field lines intersecting the disc and therefore more field lines for which $`f=0`$. If the outer boundary condition does not allow that, the magnitude of the integral can also be reduced by the field lines intersecting the disc at progressively larger radii (and smaller field magnitudes), as $`C`$ increases. In either case, the magnitude of $`B_\phi `$ is controlled and the poloidal field takes on either a more inflated or a partially open structure.

We now discuss our results for a disc with Keplerian rotation everywhere and $`R_\mathrm{d}=R_\mathrm{c}=3`$ (those cases are annotated with $`\omega =1`$ in Table 2).

### 4.2 Keplerian disc rotation profile

The calculations presented here were performed using five values of $`C`$ in the range $`3`$$`120`$ and using the Dirichlet outer boundary condition. The maximum value of $`R_{\mathrm{ext}}`$ used in these calculations decreases with $`C`$ (see Table 2). As $`C`$ increases the resolution required to resolve fully the magnetic field structure around the corotation radius increases. We have therefore restricted the external radius for $`C=30`$ or larger so that a reasonable resolution of the area around corotation can be achieved.

The maximum value of $`|B_\phi /B_z|`$ at the disc surface found in the calculations with a Keplerian disc rotation profile is $`100`$ being much larger than values considered by previous authors. Consequently we study the effect of increasing $`|B_\phi /B_z|`$ on the force–free magnetospheric field for a large region of parameter space.

Below we present results obtained with $`C=3`$ corresponding to maximum value of $`|B_\phi /B_z|=2.56`$$`2.97`$ for $`R_{\mathrm{ext}}=11`$$`80`$. When $`C=3`$ we were able to achieve the highest resolution and the largest range of $`R_{\mathrm{ext}}.`$

#### 4.2.1 Results with $`C=3`$.

All the calculations with $`C=3`$ were performed with $`Nr=200`$ and $`Nm=80`$ therefore cases with smaller $`R_{\mathrm{ext}}`$ have better resolution.

In the left panel of Fig. 2 we plot the equilibrium radial profiles of $`|B_\phi |`$ at the disc surface. Since these profiles are well resolved we can conclude that increasing the computational domain leads to smaller values of $`|B_\phi |`$. This follows from the dependence of $`B_z`$ on $`R_{\mathrm{ext}}.`$ In the right hand panel of Fig. 2 we plot the radial profile of $`B_z`$ at the disc surface for the different values of $`R_{\mathrm{ext}}`$.

The poloidal field is similar to the stellar dipole field inside corotation where $`B_\phi =0`$ but it deviates from that significantly at all other radii. Immediately outside corotation the field becomes smaller than the stellar dipole field as a result of field line inflation. At the largest radii the poloidal field eventually becomes larger than the dipole field locally. The radius at which this transition occurs moves to larger radii as $`R_{\mathrm{ext}}`$ increases indicating the influence of the outer boundary condition. The structure of the magnetic field is not radially self–similar. The departure from self–similarity in the inner part of the disc is due to the transition from $`B_\phi =0`$ inside the corotation radius to $`B_\phi `$ given by equation (14) for $`R>R_\mathrm{c}`$ while the outer part of the disc is influenced by the outer boundary condition.

Field line expansion as seen in Fig. (4) for all $`R_{\mathrm{ext}}`$ used is a consequence of field line twisting together with the assumption of a force–free equilibrium (see our discussion in Section 4.1). The inflation of field lines increases with increasing $`C`$ or field line footpoint shear until the field lines become partially or fully open (\[Aly (1985)\]; \[Wolfson & Low (1992)\]; \[Newman et al. (1992)\]; \[Lynden-Bell & Boily (1994)\]; \[Aly (1995)\]; \[Wolfson (1995)\]). In the fully open state line twisting disappears and $`B_\phi =0.`$

In Fig. 6 we plot the ratio of the final $`R_\mathrm{f}`$ to the initial or dipole $`R_\mathrm{i}`$ footpoint radius as a function of $`R_\mathrm{i}`$ for a group of field lines. The values of $`R_\mathrm{i}`$ chosen are the same as in Fig. 4 with $`R_\mathrm{i}4`$. Different curves correspond to different values of $`R_{\mathrm{ext}}`$. The ratio $`R_\mathrm{f}/R_\mathrm{i}`$ increases linearly with $`R_\mathrm{i}`$ until a turnover radius is reached where the outer boundary starts affecting the rate of inflation. For $`R_\mathrm{i}`$ larger than the turnover radius $`R_\mathrm{f}/R_\mathrm{i}`$ decreases with $`R_\mathrm{i}`$. For the cases with $`R_{\mathrm{ext}}=60`$ and $`80`$ the values of $`R_\mathrm{f}/R_\mathrm{i}`$ are similar for $`R_\mathrm{i}10`$. We therefore conclude that the maximum inflation has been reached for the innermost disc region in these cases. The maximum overall value of $`R_\mathrm{f}/R_\mathrm{i}`$ is $`3.2`$ and it corresponds to $`|B_\phi /B_z|`$ of $`2.73`$. The maximum value of $`R_\mathrm{f}/R_\mathrm{i}`$ quoted by ?) is 1.9 but this is an underestimate because of resolution problems (A.Bardou, private communication). Due to the effect of the outer boundary in our calculations it is difficult to predict how inflated the configuration would become for $`R_{ext}\mathrm{}.`$ However, we expect the field lines will remain closed within a radius of approximately $`3R_\mathrm{c}`$ for $`C=3`$.

#### 4.2.2 Results for larger values of $`C`$

We performed calculations for $`C=10`$ with $`R_{\mathrm{ext}}=60`$ and $`80.`$ Both these have $`Nr=100`$ and $`Nm=50`$ so the resolution is reduced with respect to the $`C=3`$ case. In Fig. 8 we plot the radial profiles of $`\mathrm{\Psi }_\mathrm{t}(r,0)`$ and $`B_z(r,0)`$ for $`C=10`$ and also for $`C=3`$ and $`R_{\mathrm{ext}}=80`$. For $`R_{\mathrm{ext}}=80`$ we find that $`\mathrm{\Psi }_\mathrm{t}(r,0)`$ is smaller for $`C=10`$ than for $`C=3`$ everywhere in the disc as expected. For the smaller $`R_{\mathrm{ext}}`$ value the effect of the outer boundary is stronger as it was found also for $`C=3`$. We find that $`B_z`$ decreases with $`C`$ for $`R<15`$. For larger radii the trend is reversed because of the effect of the outer boundary.

In the left panel of Fig. 10 we plot dipole and equilibrium poloidal field lines for $`C=10`$ with $`R_{\mathrm{ext}}=80`$. Field lines interior to $`2R_{ext}/3`$ are characterised by the largest flattening at small $`\theta `$ while the outer field lines are more spherical due to the effects of the outer boundary. A value of $`R_{\mathrm{ext}}`$ exceeding the outer disc radius by a factor $`10`$$`100`$ might be necessary in order to remove the effects of the outer boundary on the calculation of field line inflation (see \[Roumeliotis et al. (1994)\]).

In the right panel of Fig. 10 we plot $`R_\mathrm{f}/R_\mathrm{i}`$ as function of $`R_\mathrm{i}`$ for $`C=3`$ and $`C=10`$, each with $`R_{\mathrm{ext}}=60`$ and $`80`$. For $`R_\mathrm{i}8`$ the ratio $`R_\mathrm{f}/R_\mathrm{i}`$ is $`1.75`$ times larger for $`C=10`$ than for $`C=3`$ showing greater inflation. But for larger values of $`R_\mathrm{i}`$ this factor drops to $`1.25`$. The larger influence of the outer boundary on the equilibrium for $`C=10`$ is manifested in the smaller turnover radius observed in this case. It is not clear if the ratio of final to initial footpoint radius found for $`C=10`$ at $`R_\mathrm{i}=6`$ would continue to increase linearly with $`R_\mathrm{i}`$ if the outer boundary effect was not present. In the ?) results $`R_\mathrm{f}/R_\mathrm{i}`$ increases linearly with $`R_\mathrm{i}`$ for $`RR_\mathrm{c}.`$ In our case $`|B_\phi /B_z|`$ tends to a constant for $`RR_\mathrm{c}`$ and $`R_\mathrm{f}/R_\mathrm{i}`$ could behave similarly.

We now discuss results obtained for $`C=30`$$`120`$. These calculations were performed with $`Nr=50`$ and $`Nr=30`$ and $`R_{\mathrm{ext}}=11`$ except for one case with $`C=30`$ and $`R_{\mathrm{ext}}=21`$ which had $`Nr=120`$ and $`Nm=50`$. The disc surface equilibrium profiles of $`B_\phi `$ for $`C=3`$$`120`$ and $`R_{\mathrm{ext}}=11`$ are compared in Fig. 12. Note that the magnitude of $`B_\phi `$ varies only by a factor of $`3`$ as $`C`$ varies between $`3`$ and $`120.`$

The effect of the growing $`|B_\phi |`$ on the stream function for $`C30`$ and for $`R_{\mathrm{ext}}=11`$ is to make it larger than the dipole value, more so around $`R=R_\mathrm{c},`$ leading to an enhancement of $`B_z(r,0)`$ there. The disc surface radial profiles of $`B_z`$ for $`C=3`$$`120`$ are compared in the left panel of Fig. 14. Because the toroidal magnetic pressure increases sharply just outside $`R=R_\mathrm{c}`$ it is expected that a larger poloidal field just interior to $`R_c`$ is required for equilibrium. Thus field line deflation results. This occurs in this model because the field lines are not sheared inside corotation. For $`R>R_\mathrm{c}`$ $`B_z(r,0)`$ becomes smaller than the dipole value before increasing close to the outer boundary. Due to the increased $`|B_\phi /B_z|`$ at the disc surface, $`B_z`$ decreases with $`C`$ outside corotation as expected.

To study the dependence on the outer boundary radius for large $`C`$, we performed a calculation with $`C=30`$ and $`R_{\mathrm{ext}}=21`$ and compared the result to the case where $`R_{\mathrm{ext}}=11`$. The radial resolution is approximately the same for these cases (see Table 2). The corresponding radial profiles of $`B_z(r,0)`$ are plotted in the right panel of figure 14.

When $`R_{\mathrm{ext}}`$ is increased $`B_z(r,0)`$ values are unmodified inside $`R=10`$ so the deflation of the field lines there is independent of the outer disc radius. The exterior point at which $`B_z(r,0)`$ exceeds the dipole value moves to a larger radius for $`C=30`$ as compared to that when $`C=3`$. There is also a change to the radial dependence of $`B_z(r,0)`$ which becomes flatter at large $`R`$ for $`C30`$.

In Fig. 16 we plot the magnetic field lines for $`C=30`$ with $`R_{\mathrm{ext}}=11`$ (left panel) and $`R_{\mathrm{ext}}=21`$ (right panel). We see that the deflation of field lines for $`r<10`$ is the same in these cases. At larger radii inflation rather than deflation is observed for $`R_{\mathrm{ext}}=21`$. The field line inflation observed at $`r>10`$ for $`C=30`$ with $`R_{\mathrm{ext}}=21`$ is smaller than that found for $`C=3`$ and the same $`R_{\mathrm{ext}}`$ (compare with Fig. 4). This result is unexpected given the larger $`|B_\phi /B_z|`$ at the disc surface in the $`C=30`$ case. This is probably the result of the increased influence of the outer boundary condition on the equilibrium configuration for increasing $`C`$ values.

In the next section we will discuss results obtained using a sub–Keplerian form of $`\mathrm{\Omega }`$ and which results to a smoother radial variation of $`B_\phi /B_z`$ on the disc surface.

### 4.3 Sub–Keplerian disc rotation profile

As the magnetic field strength increases with decreasing disc radius and the magnetic stresses start to affect the radial and vertical disc equilibrium the disc is expected to ultimately corotate with the star. Interior to the corotation radius the disc will then be sub–Keplerian. In the present study we investigate the effect of sub–Keplerian rotation on $`B_\phi /B_z`$ and the equilibrium poloidal magnetic field. We modify $`\mathrm{\Omega }`$ so that it becomes sub–Keplerian for $`RR_\mathrm{c}.`$ The form of $`\mathrm{\Omega }`$ that we will use is that of ?) which is given, for $`R>R_\mathrm{d}`$, by:

$$\mathrm{\Omega }=\mathrm{\Omega }_\mathrm{K}(R_\mathrm{d})\left\{\left(\frac{R_\mathrm{d}}{R}\right)^{3/2}(1\omega )\mathrm{exp}\left[\frac{3}{2}(1\omega )^1\left(\frac{R}{R_\mathrm{d}}1\right)\right]\right\}.$$
(24)
Here $`\omega =\mathrm{\Omega }_{}/\mathrm{\Omega }_\mathrm{K}(R_\mathrm{d})`$. At $`R=R_\mathrm{d}`$ we have $`\mathrm{\Omega }=\mathrm{\Omega }_{}`$ and $`\mathrm{d}\mathrm{\Omega }/\mathrm{d}R=0`$. Interior to $`R=R_\mathrm{d}`$, the material is taken to rotate with $`\mathrm{\Omega }=\mathrm{\Omega }_{}.`$ We plot $`\mathrm{\Omega }`$ as a function of $`R`$ in the left panel of Fig. 18 for $`R_\mathrm{c}=3`$ and $`\omega =0.875`$ and $`0.2`$. The condition: $`\mathrm{\Omega }=\mathrm{\Omega }_{}`$ at $`R=R_\mathrm{c}`$ determines $`R_\mathrm{d}`$ through the chosen value of $`\omega `$.

For the largest value of $`\omega `$, $`R_\mathrm{d}`$ almost coincides with $`R_\mathrm{c}`$ and $`\mathrm{\Omega }`$ remains Keplerian until very close to $`R_\mathrm{c}`$. For smaller values of $`\omega `$ we have $`R_\mathrm{d}<R_\mathrm{c}`$ and in the case where $`\omega =0.2,`$ $`R_\mathrm{d}`$ almost reaches the stellar surface ($`R=1`$). Here as in the previous sections we have $`B_\phi =0`$ for $`R<R_\mathrm{d}`$ but a part of the differentially rotating disc may exist inside corotation and that can affect the equilibrium magnetic field. The variation of $`B_\phi /B_z`$ with $`r`$ at the disc surface can be seen in the right panel of Fig. 18 where we took $`C=30.`$

For $`\omega =0.875`$ the profile of $`B_\phi /B_z`$ is virtually unchanged with respect to the Keplerian case with $`R_\mathrm{d}=R_\mathrm{c}.`$ A larger effect on $`B_\phi /B_z`$ occurs when $`\omega =0.2.`$ In that case there is a region inside corotation where $`B_\phi /B_z`$ reverses sign.

#### 4.3.1 Results with $`\omega =0.875`$ and $`C=30`$$`120`$

These calculations have $`R_{\mathrm{ext}}=60.`$

In the left panel of Fig. 20 we plot $`B_z(r,0)`$ for $`C=30`$$`120`$ with $`\omega =0.875`$ and $`R_{\mathrm{ext}}=60`$ with the corresponding $`\mathrm{\Psi }_\mathrm{t}(r,0)`$ profiles shown in the right panel. We also plot the case $`C=3`$ with $`\mathrm{\Omega }=\mathrm{\Omega }_\mathrm{K}`$. Apart from boundary effects, $`|B_z|`$ decreases with $`C`$, the decrease being faster between $`C=60`$ and $`120`$ than between $`C=30`$ and $`60`$. These results reinforce our expectation: $`B_z0`$ for $`C\mathrm{}`$.

As seen in the right panel of Fig. 20 the disc values of $`\mathrm{\Psi }_\mathrm{t}`$ increase with $`C`$ between $`C=30`$ and $`120`$. This is contrary to theoretical prediction and to the trend observed between $`C=3`$ and $`30`$. As we will see in Section 4.4 this result depends strongly on the outer boundary condition.

The profiles obtained in the non–Keplerian case with $`\omega =0.875`$ are smoother around corotation than in the Keplerian case. This is due to the smoother radial variation of $`B_\phi /B_z`$. However, $`B_z(r,0)`$ is still larger than the dipole value there and some deflation is observed interior to $`r=R_\mathrm{c}`$. At larger radii field lines are inflated as we would expect given the increase in $`R_{\mathrm{ext}}`$.

#### 4.3.2 $`\omega =0.2`$

In this case $`R_\mathrm{d}/R_\mathrm{c}=0.34`$ and an inversion of the sign of $`B_\phi `$ may significantly affect the variation of $`B_z`$ with radius. We adopted $`C=3`$ with $`R_{\mathrm{ext}}=21`$ and $`C=30`$ with $`R_{\mathrm{ext}}=11`$ with two different resolutions (see Table 2 for details). For $`C=3`$ the difference between Keplerian and non–Keplerian cases is not significant.

In the left panel of Fig. 22 we plot the disc surface profiles of $`B_\phi `$ for $`C=30`$ and $`\omega =0.2`$ obtained with two different resolutions and for $`C=30`$ with $`\mathrm{\Omega }=\mathrm{\Omega }_\mathrm{K}`$. For $`C=30`$ the larger radial variation of $`|B_\phi |`$ in the non–Keplerian case has a more significant effect on the equilibrium. As the radial gradient of $`B_\phi ^2`$ increases for small radii a small deflation of the field lines is observed there. Eventually $`B_\phi ^2`$ decreases with radius and inflation is observed for $`R>2`$. The reverse was observed in the Keplerian case where field lines were found to be deflated everywhere in the disc (always for $`R_{\mathrm{ext}}=11`$).

In the right panel of Fig. 22 we plot the radial profile of $`B_z(r,0)`$ for the same parameters as in the left panel. We note that when $`\omega =0.2`$ the field deviates more from the dipole form around corotation, becoming larger than the dipole value locally, than in the case with $`\mathrm{\Omega }=\mathrm{\Omega }_\mathrm{K}.`$ This difference is mainly due to the larger radial variation of $`|B_\phi |`$ with radius in the non–Keplerian case. The difference observed between the profiles obtained in the high and low resolution cases is mostly due to the radial shift of the gridpoint corresponding to $`R=R_\mathrm{c}`$.

We conclude that the deflation observed in previous calculations with a Keplerian rotation profile and with $`\omega =0.875`$ is related to the restricted twisting of field lines inside the corotation radius. In the case with $`\omega =0.2`$ we have $`B_\phi 0`$ almost everywhere in the disc. Consequently field lines are much more inflated in this case.

### 4.4 A different outer boundary condition

To study the effect of the outer boundary condition on our results we have also performed calculations with the Neumann condition: $`\frac{\mathrm{\Psi }}{r}=0`$ and for $`C=3,10,60`$ and $`120`$. The parameters used in these calculations are given in Table 2. All cases have $`\omega =0.875`$ except when $`C=3`$ where a Keplerian rotation profile was used. The new boundary condition makes the disc field normal to the outer boundary. Thus we expect some of the field lines to be open at the equilibrium, at least at large disc radii far from the disc midplane.

Indeed the results show much larger inflation and opening of the field lines in this case than in the cases presented in previous sections where the Dirichlet condition was used. Also almost no deflation is observed inside corotation showing the global effect of the outer boundary condition. As seen in Fig. 24 both $`\mathrm{\Psi }_\mathrm{t}(r,0)`$ and $`B_z(r,0)`$ decrease with $`C`$ monotonically. For $`C=120`$ the equilibrium seems to have reached a limiting configuration with almost all field lines open.

Comparing the disc radial profiles of $`\mathrm{\Psi }_\mathrm{t}(r,0)`$ depicted here with those in Fig. 20 we note that $`\mathrm{\Psi }_\mathrm{t}(r,0)`$ decreases considerably when the second boundary condition is used. Also $`\mathrm{\Psi }_\mathrm{t}(r,0)`$ becomes approximately independent of $`r`$ everywhere outside the corotation area which is the limiting behaviour that we expect for large values of $`C`$. Comparing the disc radial profiles of $`B_z(r,0)`$ between figures 20 and 24 we note that for $`C=3`$ the magnetic field is similar for the two boundary conditions for $`r10`$. Further out $`B_z(r,0)`$ is smaller for the second boundary condition. This difference increases with $`C`$. We also note that $`B_z(r,0)`$ has a stronger dependence on $`r`$ when the second boundary condition is used.

In the first three panels of Fig. 26 we plot the magnetic field lines for $`C=3,10`$ and $`60`$ respectively, using the same initial (dipole) field lines as in Fig. 4 (with $`R_{\mathrm{ext}}=60`$). For $`C=3`$ the level of inflation is comparable for the two boundary conditions for $`r\begin{array}{c}<\\ \end{array}10`$. For larger initial footpoint radii field lines become much more inflated for the second outer boundary condition and eventually they open up for $`r15`$. As $`C`$ increases more field lines become open as seen by comparing the cases with $`C=10`$ and $`60`$. The case with $`C=120`$ is virtually identical to the $`C=60`$ case. Some field lines remain closed even when $`C=60`$ as can be seen by plotting more contour levels (see lower right panel of Fig. 26). The closed field lines correspond to initial radii which are very close to the inner disc radius. Consequently $`ff^{}`$ is non-zero for an increasingly smaller range of $`\mathrm{\Psi }_\mathrm{t}`$ values as $`C`$ increases. Calculations for larger values of $`C`$ could only be performed accurately using considerably larger numerical resolution. However, our present calculations have already shown that field lines at equilibrium will be approximately open for $`|B_\phi /B_z|10^2`$ if a Neumman outer boundary condition is used.

The ratio $`|B_\phi |/B_\mathrm{p}`$ is found to be less than unity for $`\mu <0.1`$. Therefore $`B_\phi 0`$ in most of the corona as expected on open field lines. In the limit of a large twist the disc field acts as to expel the poloidal magnetic field from the disc and consequently to annihilate the toroidal magnetic field, as long as the disc field is permitted to penetrate the outer boundary. In the limit of large $`C`$ the poloidal field resembles the solution $`Z_1`$ given by ?). This corresponds to a perfectly conducting disc with central hole immersed in a central dipole field where some flux escapes to infinity such that the field is non singular at the disc inner edge.

### 4.5 Modification of the spin–down and spin–up torques

For the force–free equilibria calculated above, the poloidal magnetic field differs significantly from the original stellar dipole form. As a consequence, the local magnetic torque acting on the disc may also differ significantly from what one obtains assuming $`B_zB_{\mathrm{dp}}`$. Since this torque regulates the spin evolution of the central star it is important to examine the possible effect of our results on the calculation of the total torque experienced by the star.

The total magnetic torque $`N`$ can be split into two components of opposite sign. The first arises inside corotation where angular momentum is transfered from the disc to spin up the star. The second arises outside corotation where angular momentum is transferred to the disc. In order for stellar accretion to occur we require $`R_\mathrm{d}<R_\mathrm{c}`$. In one set of calculations we assumed $`R_\mathrm{d}=R_\mathrm{c}`$ and adopted Keplerian rotation. In that case the spin–up torque would be zero. In reality we expect $`R_\mathrm{d}/R_\mathrm{c}0.91`$$`0.97`$ (i.e. \[Wang (1995)\]) giving rise to a relatively small spin–up torque. This torque is probably small compared to that arising from direct accretion onto the star.

We now calculate the total torques for the case where the disc is assumed to be Keplerian and truncated at the corotation radius.

#### 4.5.1 The Keplerian case

The total magnetic torque on the star is in general given by (e.g. \[Ghosh & Lamb (1979b)\]):

$$N=_{R_\mathrm{d}}^{R_{\mathrm{ext}}}B_\phi B_zR^2dR.$$
(25)
We adopt dimensionless parameters such that the torque is given in units of $`N_{}B_{}^2r_{}^3`$. Of course all field values are calculated on the disc surface ($`\mu =0`$). For the cases presented here the first boundary condition has been used. When $`R_\mathrm{d}=R_\mathrm{c},`$ as is considered here, $`N`$ is positive definite and only acts to spin down the star.

The total torque $`N`$ can be calculated analytically for a dipole field with $`B_\phi `$ given by equation (14) and is given by:

$$N^{\mathrm{dp}}=\frac{C}{3}\left[\frac{1}{3}+\frac{2}{3}\left(\frac{R_{\mathrm{ext}}}{R_\mathrm{c}}\right)^{9/2}\left(\frac{R_{\mathrm{ext}}}{R_\mathrm{c}}\right)^3\right]R_\mathrm{c}^3.$$
(26)
In the left hand panel of Fig. 28 we plot $`N/C`$ and $`N^{\mathrm{dp}}/C`$ as a function of $`C`$ for $`R_{\mathrm{ext}}=11`$. Although $`N^{\mathrm{dp}}/C`$ is constant with $`C`$, $`N/C`$ decreases with increasing $`C`$ as the poloidal magnetic field in the force–free equilibria decreases with increasing $`C`$. The ratio of force–free to dipole torque $`N/N^{\mathrm{dp}}`$ for $`C30`$ varies between $`0.14`$$`0.013`$ which is significantly smaller than the same ratio in the case with $`C=3`$. Although the effect of the outer boundary condition is to suppress inflation in all cases with large $`C`$ values and $`R_{\mathrm{ext}}=11`$, the effect of the twisting of the field lines on the final equilibrium is an increased reduction of the total torque with respect to its corresponding dipolar value.

In the right panel of Fig. 28 we plot $`N`$ and $`N^{\mathrm{dp}}`$ where we see that the behaviour of $`N`$ as function of $`C`$ is very different (in fact almost reversed) to that expected for a dipolar poloidal field. From this figure we see that $`N`$ increases slowly for $`C15`$ and it decreases for larger values of $`C`$. As the values of $`N`$ for $`C30`$ are between $`10`$ and $`100`$ times smaller than the dipolar values the increasing field line twisting has an important effect on the spin down of the star. The magnetic breaking timescales calculated are longer and therefore magnetic breaking models might have to be reevaluated. However, we have not taken into account here the disc region inside corotation. Also the results described here for $`C30`$ have $`R_{\mathrm{ext}}=11`$ which is relatively small and therefore the effect of the outer boundary condition is important.

#### 4.5.2 The sub–Keplerian case

In the case where $`\mathrm{\Omega }\mathrm{\Omega }_\mathrm{K}`$ the total magnetic torque includes both spin–up and spin–down components. For the case with $`\omega =0.875`$ the spin–up component is of no significance. In the following we will compare the total magnetic torques calculated using a sub–Keplerian rotation profile with $`\omega =0.875`$ to the corresponding dipolar values. Calculations with similar parameters were performed with both outer boundary conditions.

In Fig. 30 we plot $`N/C`$ and $`N^{\mathrm{dp}}/C`$ as a function of $`C`$. The values in the left panel of this figure are from calculations performed with the first outer boundary condition. Here, as for the Keplerian case, we find that $`N/C`$ decreases with $`C`$ and that the ratio $`N/N^{\mathrm{dp}}`$ ranges from $`0.56`$ for $`C=3`$ to $`0.05`$ for $`C=120`$. The larger ratios here are mainly a consequence of the use of a sub–Keplerian rotation profile. However, For $`C30`$ we have again $`N/C1/C`$ approximately. Therefore, the trend for a smaller total torque is clear and the difference between force–free and dipole cases is again significant for $`C1`$. For $`B_\phi B_z`$ the difference is not as large although we have to reserve final judgement until calculations with much larger values of $`R_{\mathrm{ext}}`$ are performed.

For calculations performed with the second outer boundary condition we found significantly smaller values of $`B_z(r,0)`$ which decreases faster with $`r`$ than when the first boundary condition is used. The steeper decline of $`B_z(r,0)`$ with $`r`$ is more noticeable in the cases with large $`C`$ values. For large $`C`$ values we therefore expect much smaller values of $`N/C`$ to arise when the second boundary condition is used. In the right panel of Fig. 30 we plot $`N/C`$ and $`N^{\mathrm{dp}}/C`$ as a function of $`C`$ for the calculations presented in Section 4.4. The ratio $`N/N^{\mathrm{dp}}`$ varies in the range $`0.45`$$`0.003`$. Although this ratio is almost unity for $`C=3,`$ as was the case for the first boundary condition, it becomes progressively smaller for larger $`C`$ values. For large values of $`C`$ we have $`N/C1/C^{3/2}`$ when the second boundary condition is used.

In summary, when $`|B_\phi /B_z|`$ on the disc surface is large the magnetic torque acting on it will be much smaller than that estimated assuming $`B_zB_{\mathrm{dp}},`$ especially when external conditions enable the opening of field lines.

## 5 Discussion

We first summarise our results. For a Keplerian rotation profile with the first boundary condition the maximum inflation observed in our calculations with $`C=3`$ is such that $`R_\mathrm{f}/R_\mathrm{i}=3.2.`$ This is a confirmation of the results of ?). The expansion of the field lines is restricted by the outer boundary in this case. However, for $`|B_\phi /B_z|3`$ at the disc surface the field lines with $`r3R_\mathrm{c}`$ do not vary much once $`R_{\mathrm{ext}}>10.`$ Thus although significantly inflated they will remain closed in the force–free equilibrium. Inflation increases with $`C`$ so that when $`|B_\phi /B_z|=10,`$ a maximum $`R_\mathrm{f}/R_\mathrm{i}5`$ is found for $`R_\mathrm{i}=8.`$

In calculations with $`C=30`$$`120`$ and $`R_{\mathrm{ext}}=11`$ a new effect was found. Inner field lines become deflated rather than inflated with respect to the initial dipolar configuration. We expect the level of deflation, and consequent increase in poloidal magnetic pressure, to depend on the rate of increase of toroidal magnetic pressure, which it has to balance around $`R=R_\mathrm{c}`$. We demonstrated, by varying $`R_{\mathrm{ext}}`$, that the deflation is not related to the outer boundary. We found that the level of field line deflation for $`r<10`$ was restricted when a smoother non Keplerian rotation profile was used. However, some level of deflation could be typical of the field structure in the region near the corotation radius. That applies especially to discs with $`\eta /\nu 1`$ and with $`R_\mathrm{d}R_\mathrm{c}`$ when a Dirichlet outer boundary condition is used.

In all the calculations presented here the equilibrium poloidal field is more nearly radial above the disc surface than the stellar dipole field. The lengthening of the field lines is accompanied by a reduction of $`|B_\phi |/B_\mathrm{p}`$ for increasing $`\mu .`$ We expect that the force–free condition constrains $`|B_\phi |/B_\mathrm{p}`$ to remain of order unity even for $`|B_\phi /B_z|1`$ at the disc surface.

We also performed calculations with a Neumman outer boundary condition that required the disc field to be radial at the outer boundary as it might occur if a wind existed. In these cases, other parameters being equal, the equilibrium field structure was more open. The ratio $`|B_\phi |/B_\mathrm{p}`$ was found to be less than unity for $`\mu <0.1`$ and $`0`$ in most of the corona as expected on open field lines. In the limit of a large twist the poloidal magnetic field tends to be expelled from the disc and consequently the toroidal magnetic field is annihilated, as long as the disc field is permitted to penetrate the outer boundary.

We also found that for large values of $`|B_\phi /B_z|`$ on the disc surface the total magnetic torque is much smaller than that estimated assuming $`B_zB_{\mathrm{dp}}`$. For $`C=120`$ the torque is only $`0.003`$ of the dipole value when the second boundary condition is used and $`0.0130.05`$ times the dipole value when the first boundary condition is used (depending on the form of $`\mathrm{\Omega }`$). Since reasonable values of $`\eta `$ for circumstellar discs tend to produce very large values of $`|B_\phi /B_z|`$ on the disc surface, we believe that the modification of the total torque discussed here will be significant if a force free equilibrium can be attained. We therefore expect significant modification of the magnetic breaking times of T Tauri stars to result from the large twisting of the magnetospheric field lines.

### 5.1 The spin–down of neutron stars

Recent observations of accreting neutron stars (\[Nelson et al. (1997)\]; \[Chakrabatry et al. (1997)\]) suggest that the total torque between star and disc can oscillate producing alternating spin–up and spin–down phases with little evidence of any correlation with the mass accretion rate. ?) have suggested that this may be due to variations in the structure of the magnetosphere.

The work presented here indicates that a large range in the value of the spin–down torque may be obtained depending on the outer boundary condition and effective disc resistivity. In cases favouring large field line inflation and open field lines the spin–down torque resulting from the disc-magnetosphere interaction becomes very small. Thus if the magnetosphere oscillates between such a state and one with smaller inflation, oscillations in the spin–down torque and hence the total torque may be produced. The change from open to closed field lines may also affect the magnetic torque that arises from a wind, if such a wind is produced. Changes in the magnetospheric configuration might occur through reconnection through the disc midplane, variations in outer conditions being more or less favourable for open field configurations, or variations in the internal disc resistivity.

## Acknowledgments

This work was supported by the European Union grant ERBFMRX-CT98-0195. V.A. acknowledges support by the State Scholarships Foundation (IKY) of the Hellenic Republic through a postgraduate studentship. The authors are grateful to Anne Bardou for useful discussions.


text:
# COSMIC VELOCITIES 2000: A REVIEW To Appear in Proceedings of the XXXVth Rencontres de Moriond: Energy Densities in the Universe

## 1 Background

The study of cosmic flows emerged as a distinct subfield of cosmology in the late 1970s, spurred by the discovery of the Tully-Fisher (TF) and Faber-Jackson relations for spiral and elliptical galaxies, respectively. With these empirical correlations between galaxy luminosity and internal velocity as distance indicators, one could measure redshift-independent distances for galaxies out to tens or even hundreds of megaparsecs. Both the TF relation and the successors of Faber-Jackson—the $`D_n`$-$`\sigma `$ and Fundamental Plane (FP) relations—yield distances with $`15`$$`20\%`$ accuracy. Although hopes for further improving the accuracy of TF and FP have proved unfounded, these distance indicators have remained the workhorses of cosmic velocity analysis for two decades. Only in recent years have newer and more accurate distance indicator methods—in particular Type Ia Supernovae and Surface Brightness Fluctuations—begun to complement (not supplant) TF and FP, as discussed further below.

The early scientific emphasis in flow studies was on determining the amplitude of Virgo infall (e.g., Tonry & Davis 1981; Aaronson et al. 1982). The Virgo cluster and its environs were then thought to dominate the local flow field. It has since been recognized that the local velocity field is more complex. There are a number of nearby attractors and voids, as well as the tidal effect of distant mass concentrations. The most sophisticated models of the local velocity field now use gravitational instability theory to predict peculiar velocities on the basis of the galaxy density field observed in redshift surveys. The redshift survey most widely used for this purpose is the IRAS redshift survey, both in its older, 1.2 Jy (Fisher et al. 1995) and newer PSCz (Saunders et al. 2000) incarnations. If one assumes IRAS galaxies trace mass and adopts the approximations of linear theory, comparison of predicted velocities with the observed ones constrains the parameter $`\beta _I=\mathrm{\Omega }_m^{0.6}/b_I.`$ Here, $`b_I`$ is the biasing parameter for IRAS galaxies, a measure of whether IRAS galaxies are more ($`b_I>1`$) or less ($`b_I<1`$) clustered than mass. (The subscript is needed because different redshift samples have different clustering amplitudes, and thus different biasing parameters.)

Measurement of $`\beta _I`$ has thus become one of the main thrusts of cosmic flow analysis in the 1990s. Early in the decade, when there was a widespread theoretical prejudice in favor of an $`\mathrm{\Omega }_m=1`$ universe, it was thought peculiar velocities might be the key to proving it. On large scales, the argument went, galaxies should trace the dominant dark matter, and the large-scale peculiar field should therefore reflect the underlying Einstein-de Sitter nature of the universe. The earliest efforts in this direction indeed seemed to bear out this suspicion, finding $`\beta _I1.3,`$ suggestive of $`\mathrm{\Omega }_m=1`$ (Dekel et al. 1993). More recently, however, cosmic flow-based estimates of $`\beta _I`$ have often, though not invariably, produced values in the $`0.4`$$`0.6`$ range, consistent with a low-density cosmology. In §3 I will summarize recent work done on this problem, and explain why I believe the low $`\beta _I`$ values are more likely to be correct.

Another aim of cosmic flow studies first arose serendipitously: efforts to detect bulk flow on very large ($`\stackrel{>}{}100h^1\mathrm{Mpc}`$) scales. This work was propelled by the discovery, in 1987, of a bulk flow stretching across the sky by the “7-Samurai” (7S) group (Dressler et al. 1987). Working with the newly discovered $`D_n`$-$`\sigma `$ relation, the 7S found that elliptical galaxies out to $`3000\mathrm{km}\mathrm{s}^1`$ redshift exhibited fairly uniform Hubble expansion from the vantage point of the Local Group (LG) barycenter. But the LG is known to move at $`630\mathrm{km}\mathrm{s}^1`$ with respect to the Cosmic Microwave Background radiation (CMB), as indicated by the CMB dipole anisotropy. Thus, the implication of the 7S findings was that the ellipticals were streaming at $`600\mathrm{km}\mathrm{s}^1`$ relative to the CMB. If there is any validity to Big Bang cosmology, though, the CMB defines an absolute standard of rest—the “cosmic rest frame,” as it were. The 7S finding thus inaugurated a long-standing puzzle in cosmology: how can large-scale, coherent bulk flows exist in a universe that seems to be so uniform on large scales?

The puzzle deepened in the following decade, with a number of groups confirming a 7S-like bulk flow, and, in several cases, finding that it continued to scales three or four times the 7S volume. The current controversy may be framed as “what is the scale of the largest bulk flow,” or, equivalently, as one of convergence scale: at what distance are the galaxies within a spherical shell finally at rest in the CMB frame? Some astronomers (this author included) have been led to wonder whether such a convergence scale existed; perhaps the CMB-defined “cosmic rest frame” was offset by $`600\mathrm{km}\mathrm{s}^1`$ from the frame in which uniform Hubble expansion is observed—a conjecture which, if correct, would call into question some of the fundamental tenets of cosmology. Fortunately—at least if you like agreement between theory and observation—it now appears that the convergence to the CMB has been detected at a distance of $`50`$$`60h^1\mathrm{Mpc}.`$ Or at least, that is what I will argue in §2, when I summarize results from newly completed surveys.

## 2 The Scale of the Largest Bulk Flows

First we should ask, from the perspective of cosmology, Why are bulk flows interesting? In particular, what does their convergence scale tell us?

The answer lies in the sensitivity of bulk flows to long-wavelength modes of the mass fluctuation power spectrum. Mass conservation in the linear regime of gravitational instability tells us that

$$𝐯=\mathrm{\Omega }_m^{0.6}\delta ,$$
(1)
where $`𝐯`$ is the peculiar velocity vector and $`\delta `$ is the mass density contrast. The corresponding equation in Fourier space is $`𝐤𝐯_k=\mathrm{\Omega }_m^{0.6}\delta _k,`$ where the subscript denotes Fourier transform. Thus, $`|𝐯_k|\delta _k/k,`$ which is to say, long-wavelength perturbations (small $`k`$) have a larger impact on large-scale peculiar velocities than they do on mass fluctuations.

One can flesh out these ideas by calculating the mean square bulk velocity on a scale $`R:`$

$$v^2(R)=\frac{\mathrm{\Omega }_m^{1.2}}{2\pi ^2}_0^{\mathrm{}}P(k)\stackrel{~}{W}^2(kR)𝑑k,$$
(2)
where $`P(k)`$ is the mass fluctuation power spectrum and $`\stackrel{~}{W}^2(kR)`$ is the Fourier transform of a top-hat window of radius $`R.`$ In Figure 1 the rms expected bulk velocity, $`V_{rms}(R)=\sqrt{v^2(R)},`$ is plotted against $`R`$ (left panel) and $`\mathrm{\Omega }_m`$ (right panel) for COBE-normalized CDM power spectra. The left panel assumes a canonical $`\mathrm{\Omega }_m=0.3,`$ $`\mathrm{\Omega }_\mathrm{\Lambda }=0.7`$ universe. For comparison, the rms density fluctuation $`\sigma __M(R)=(\delta M)/M)^2(R)`$ is also plotted. Note that $`\sigma __M`$ drops much more rapidly with scale than does $`V_{rms}.`$ This is a result of the sensitivity of large-scale bulk flow to long-wavelength modes of the power spectrum. An observational consequence is that while redshift surveys have difficulty probing the power spectrum on scales $`\stackrel{>}{}100h^1\mathrm{Mpc},`$ bulk flow studies can in principle do so.

Closer inspection of the left panel also shows why large-amplitude ($`V\stackrel{>}{}500\mathrm{km}\mathrm{s}^1`$), large-scale $`R\stackrel{>}{}60h^1\mathrm{Mpc}`$ are potentially problematic for standard structure-formation scenarios. In the CDM-type models shown here, there simply isn’t enough large-scale power to drive such flows. Stated another way, the CDM universe approaches homogeneity sufficiently rapidly with increasing scale that the coherence scale of bulk flows should be a few tens of megaparsecs at most. Moreover, the right panel shows, perhaps counterintuitively, that boosting the matter density doesn’t help. (This is a consequence of imposing COBE-normalization; normalizing the power spectrum to the cluster abundance does not substantially change our conclusions.)

However, high-amplitude bulk motions on scales $`\stackrel{>}{}100h^1\mathrm{Mpc}`$ are just what were found by three surveys conducted in the early and mid 1990s. The first of these, and the one which has achieved the greatest notoriety, was the Brightest Cluster Galaxy (BCG) survey by Lauer & Postman (1994; LP). LP used the BCGs as standard candles, thus obtaining distances to over 100 Abell clusters out to $`15,000\mathrm{km}\mathrm{s}^1`$ redshift, and found them to be coherently moving relative to the CMB at $`700\mathrm{km}\mathrm{s}^1.`$ The combination of large scale and high amplitude places the LP result far above the expected $`V_{rms}`$ values in Figure 1.

More recently, two large surveys of cluster galaxies appeared to confirm the scale and amplitude (but not the direction) of the LP result. The Streaming Motions of Abell Clusters (SMAC) survey of Hudson et al. (1999) used the FP relation to measure cluster ellipticals to about the same depth as the LP survey, and found a $`600\mathrm{km}\mathrm{s}^1`$ flow, in roughly the same direction as the that found by 7S a decade earlier. And in a Tully-Fisher survey of cluster galaxies in a shell between $`9000`$ and $`13,000\mathrm{km}\mathrm{s}^1,`$ I (Willick 1999, LP10K) found a streaming motion of about $`700\mathrm{km}\mathrm{s}^1`$ in roughly the same direction as SMAC. The above results are summarized in Table I, along with another bulk flow measurement based on more nearby galaxies from the Mark III Catalog of TF and $`D_n`$-$`\sigma `$ data (Willick et al. 1997), as measured by the POTENT algorithm (on which more in § 3).

Table I. Recent Bulk Flow Measurements

| Survey | $`R`$ ($`\mathrm{km}\mathrm{s}^1`$) | $`V_B`$ ($`\mathrm{km}\mathrm{s}^1`$) | Comments |
| --- | --- | --- | --- |
| Lauer-Postman (LP) | 15000 | 700 | BCG |
| Willick (LP10K) | 12000 | 700 | TF |
| Hudson et al. (SMAC) | 14000 | 600 | FP |
| Dekel et al. (POTENT) | 6000 | 350 | TF+$`D_n`$-$`\sigma `$ |

The four results above argue for large bulk flows, but are not fully consistent. The SMAC, LP10K, and POTENT/Mark III flows agree in direction, but the LP flow is nearly orthogonal. Also, the smaller amplitude of the smaller scale POTENT/Mark III measurement is puzzling; one would expect bulk flow amplitude to diminish monotonicaly with scale (Figure 1). Even without evidence to the contrary, then, the above results are less than convincing.

More importantly, however, 1999 saw the announcements of new survey results that contradict the findings of Table 1. The nearby Type Ia Supernova (SN Ia) data have accumulated, and as reported by Riess (2000), the sample of SN Ia distances within $`10,000\mathrm{km}\mathrm{s}^1`$ redshift show no evidence for bulk flow. This is quite important, because the SN Ia data are of a fundamentally different nature than the other distance indicators employed in cosmic flow studies. Also, as shown in the cosmological context, SN Ia have small scatter, $`0.15`$ mag. The EFAR FP survey of Colless et al. (2000) similarly finds no bulk motion on a large scale, and in particular is inconsistent with LP at better than 99% confidence. Similar results have been obtained by Dale, Giovanelli, and coworkers from their extensive TF surveys (as summarized by Dale & Giovanelli 2000). The key findings of these surveys, plus that of the Shellflow survey, to which I turn next, are summarized in Table 2.

Table II. More Recent Bulk Flow Measurements

| Survey | $`R`$ ($`\mathrm{km}\mathrm{s}^1`$) | $`V_B`$ ($`\mathrm{km}\mathrm{s}^1`$) | Comments |
| --- | --- | --- | --- |
| Riess et al. | $`10000`$ | $`0`$ | SN Ia |
| Courteau et al. (SHELLFLOW) | $`6000`$ | $`0`$ | TF |
| Colless et al. (EFAR) | $`10000`$ | $`0`$ | FP |
| Dale & Giovanelli (SFI) | $`6000`$ | $`0`$ | TF |
| Dale & Giovanelli (SCI/SCII) | $`15000`$ | $`0`$ | TF |

A group of us who were active in cosmic velocity measurement and analysis (M. Strauss, S. Courteau, M. Postman, D. Schlegel, and myself) realized in 1995 that a critical issue had not been properly addressed in the then-extant Tully-Fisher surveys: the need for extremely uniform data across the sky. We proposed for and were granted extensive NOAO time for the Shellflow project, a TF survey of a shell of 300 galaxies between $`5000`$ and $`7000\mathrm{km}\mathrm{s}^1`$ redshift. By using identical observational setups from northern and southern hemisphere NOAO telescopes we ensured that data nonuniformity could not produce spurious peculiar velocities.

In a recent paper (Courteau et al. 2000) we reported our main result: the bulk flow of the shell centered at $`6000\mathrm{km}\mathrm{s}^1`$ is $`70_{70}^{+100}\mathrm{km}\mathrm{s}^1,`$ i.e., is consistent with being at rest in the CMB frame. The residuals with respect to a fit that assumes pure Hubble flow in the CMB frame are shown in Figure 2. The points are everywhere consistent with being due to TF scatter, not coherent peculiar velocities. Inflowing and outflowing points are well mixed at all positions, indicating the absence of coherent motions.

Having worked with the Shellflow data myself, I am confident that it is of very good quality, and am convinced that systematic errors have a very small effect, if any on our results. I therefore find the conclusion of insignificant bulk flow highly persuasive. The Shellflow results refer to a specific scale, namely, a $`60h^1\mathrm{Mpc}`$ sphere. They do not directly test the LP, SMAC, and LP10K results of Table 1. However, because abundant evidence indicates that the universe approaches homogeneity monotonically with increasing scale size, I believe that the Shellflow result, if correct, is physically inconsistent with the LP, SMAC, and LP10K findings, and that the latter are therefore not to be taken at face value.

One should not, however, judge the LP, SMAC, and LP10K authors harshly (and I have an obvious reason for hoping you don’t!). Measuring distortions of the Hubble expansion that are at the level of a few percent is a challenging task, given the limitations of the distance indicator techniques we work with. False detections, if that is what they are, will most likely be seen in hindsight as inevitable products of initial efforts at a very difficult measurement.

## 3 The Value of $`\beta _I`$

In considering the scientific implications of large-scale bulk flow surveys, I made no mention of comparison with redshift surveys. That is because the full-sky redshift surveys we have are not reliable enough, at distances $`\stackrel{>}{}150h^1\mathrm{Mpc},`$ to predict what the bulk flows should be on such scales. The situation changes when we talk about the velocity field within $`50h^1\mathrm{Mpc}.`$ At these distances both the peculiar velocity and galaxy density fields are mapped with enough accuracy to do a detailed comparsion. The goals of this comparison are (1) to verify that the two maps are compatible with the gravitational instability paradigm, and (2) assuming they are, to measure $`\beta _I=\mathrm{\Omega }_m^{0.6}/b_I,`$ and if possible to go further and measure $`\mathrm{\Omega }_m`$ itself.aaaIt will be assumed in what follows that, unless otherwise specified, the redshift survey in question is one of IRAS galaxies.

The theoretical bases of the comparison are either of two forms of the linear velocity-density relation. One is the differential form, Eq. (1), which for comparison of observables is written as

$$𝐯(𝐫)=\beta _I\delta _{I,g}(𝐫),$$
(3)
where $`\delta _{I,g}(𝐫)`$ is the galaxy overdensity as determined by an IRAS redshift survey. The other is the corresponding integral form,

$$𝐯(𝐫)=\frac{\beta _I}{4\pi }d^3𝐫^{}\frac{\delta _{I,g}(𝐫^{})(𝐫^{}𝐫)}{\left|𝐫^{}𝐫\right|^3}.$$
(4)
In both Eqs. (3) and (4), the spatial positions $`𝐫`$ are assumed to be measured in velocity units—i.e., they are distances in units of the Hubble velocity. In this way, the Hubble constant itself is removed from the analysis, which is (needless to say) a useful simplification.

Although the two forms of the velocity-density relation are equivalent, they lead to two rather different analytical approaches to measuring $`\beta _I.`$ Because the controversy between the “high” and the “low” (see §1) values of $`\beta _I`$ appears to revolve around this distinction, it is worth taking a moment to understand it. To apply Eq. (3), one needs to map out the three-dimensional velocity field $`𝐯(𝐫),`$ differentiate it, and finally compare it to the galaxy density field to determine $`\beta _I.`$ Because only the radial component of $`𝐯`$ is observable, one first needs an ansatz for “three-dimensionalizing” the inherently one-dimensional velocity data. An elegant approach to this problem was developed by Bertschinger & Dekel (1989; see Dekel 1994 for a review), who argued that the large-scale velocity field should be irrotational and thus expressible as the gradient of a potential function, which could itself be computed by integrating the observed, radial velocities along rays. This algorithm, known as POTENT, thus produces a 3D velocity field, smoothed on a rather large (typically $`12h^1\mathrm{Mpc}`$) scale, which can be differentiated and then used in Eq. (3).

Eq. (4) suggests a different approach. Rather than heavily processing the velocity data, one carries out the indicated integration using the redshift survey data. One thus obtains a predicted peculiar velocity field as a function of $`\beta _I,`$ of which only the radial component, $`u(𝐫)=𝐯(𝐫)𝐫/r,`$ is used in the subsequent analysis. The predicted $`u(𝐫)`$ is compared with the observed radial peculiar velocities from the TF (FP, etc.) data sets; the final estimate of $`\beta _I`$ is that which yields the closest match predictions and observations.

The two approaches to measuring $`\beta _I`$ are known as the density-density (d-d) and velocity-velocity (v-v) comparisons. It is notable that d-d comparisons, all done using POTENT to reconstruct the 3D velocity field, have consistently produced values of $`\beta _I`$ consistent with unity—and thus, to the extent $`b_I`$ is itself not so different from one, implicit estimates of $`\mathrm{\Omega }_m`$ near unity as well. The 1993 paper by Dekel et al. was already mentioned (see §1). A more recent application of POTENT, using improved peculiar velocity data, was that of Sigad et al. (1998), who found $`\beta _I=0.89\pm 0.1.`$

Since about 1995, several v-v alternatives to the POTENT approach have been developed for measuring $`\beta _I.`$ They have differed in the way in which the IRAS (or other) redshift data are used to predict peculiar velocities, and in the way the predicted and observed peculiar velocities are compared. But as v-v methods they have more in common with one another than with POTENT; in particular, the heavy computational work is done with the redshift data, while the TF (FP, etc.) data are used only in a limited, statistical sense. These methods include the Least Action Principle approach of Shaya, Tully, & Peebles (1995), who obtain $`\beta _I0.35\pm 0.1`$), the VELMOD method of Willick et al. (1997) and Willick & Strauss (1998), who obtain $`\beta _I=0.5\pm 0.05,`$ and the ITF method (Davis, Nusser, & Willick 1996) which, applied to the Type Ia Supernova data produced a value of $`\beta _I=0.4\pm 0.1.`$

Figure 3 shows representative results from two v-v analyses. The left panel shows the VELMOD results from Willick & Strauss (1998). The statistic plotted is $`=2\mathrm{ln}P,`$ where $`P`$ is the probability of observing the TF data given the IRAS velocity model. (The calculation is done for the “inverse” TF relation, which is immune to selection biases.) The minimum of the curve occurs at the maximum likelihood value of $`\beta _I,`$ and its curvature yields the $`1\sigma `$ error, as indicated on the plot.

The right panel shows the first use of the new Surface Brightness Fluctuation (SBF) data set in $`\beta `$-measurement, as reported by Blakeslee et al. (2000). As was stated above in connection with SN Ia data, it is essential to carry out these experiments with independent distance indicators, and the SBF method is quite distinct from Tully-Fisher. Moreover, like SN Ia, SBF distances are on average about twice as accurate as TF distances, and have the additional advantage of having well-determined errors. In Figure 3, the SBF data are compared with velocities predicted by IRAS and the Optical Redshift Survey (ORS). The minimum $`\chi ^2`$ is achieved for $`\beta _I=0.44\pm 0.08,`$ consistent with the VELMOD TF results. (The ORS galaxies are more clustered than IRAS galaxies, and therefore yield a lower $`\beta `$ value.)

A truly convincing explanation for this discrepancy between the v-v and d-d $`\beta `$ values has not yet been found. It seems to me, however, that the v-v comparison is a more robust procedure. In it, the intensive data manipulation is done on the redshift survey data, which is by its nature more reliable than the distance indicator data (redshifts errors are fractionally very small; redshift-independent distance errors are always $`20\%`$). In the d-d comparison, it is the data set with much larger scatter and non-Gaussian errors that is subjected to intensive manipulation. In particular, the d-d comparison requires that noisy data first be smoothed, then integrated, and then differentiated in three dimensions. There is ample opportunity in this procedure for errors to propagate. The same noisy data in the v-v studies are left virtually untouched, save for the statistical comparison with their predicted values. This procedure is far more stable. For this reason I consider the low $`\beta `$ values derived from the v-v analyses to be more reliable.

## 4 Summary

I have argued that the two major controversies in cosmic flow analysis have been largely resolved in the last few years. With regard to bulk flows, most recent surveys show convergence to the CMB frame by a distance of $`60h^1\mathrm{Mpc}.`$ A corollary is that the observed bulk motions do not require more large-scale power than is provided by COBE-normalized CDM density fluctuation spectra. With regard to the value of $`\beta _I=\mathrm{\Omega }_m^{0.6}/b_I,`$ a number of independent analyses now suggest a low value, $`\beta _I0.4`$$`0.5.`$ If the IRAS galaxies are nearly unbiased with respect to mass, a reasonable if not airtight hypothesis, these $`\beta `$-values imply a density parameter $`\mathrm{\Omega }_m0.2`$$`0.3.`$

A word of caution is in order, however. The above conclusions represent a consensus view, not a unanimous one. The bulk flow detections listed in the first three rows of Table 1 have not been in any sense “refuted,” which is to say as far as we know there is nothing wrong with the data. New surveys, such as FP200 (see http://astro.uwaterloo.ca/$``$mjhudson/fp200 for details) will, it is hoped, settle the issue definitively. Similarly, while a majority of recent velocity-density comparisons favor low $`\beta _I,`$ the reason for the discrepant POTENT result, $`\beta _I0.9,`$ is not well understood. Tests of methods such as POTENT and VELMOD using N-body simulations are under way, and may clarify things. Moreover, the coming decade will bring much larger TF data sets from the DENIS and 2MASS infrared surveys. As always, these new data sets, if they live up to their promise, are our best hope for putting any remaining contoversy to rest.

## Acknowledgments

I am grateful to my collaborators on the Shellflow project, Stéphane Courteau, Michael Strauss, Marc Postman, and David Schlegel, as well as to Michael Hudson, Marc Davis, Avishai Dekel, John Tonry and Alan Dressler for enlightening discussions over the last several years. Special thanks go to Stéphane Courteau for organizing the highly successful Cosmic Flows 99 conference in Victoria, B.C. last summer. My research is supported by a Cottrell Scholarship of Research Corporation, NSF grant AST96-17188, and a Terman Fellowship from Stanford University.