text:
# Stability of accretion discs threaded by a strong magnetic field

## 1 Introduction

In the formation and evolution of Young Stellar Objects magnetic fields are thought to play a key role. Magnetic fields, beside rotation and thermal pressure, stabilize molecular clouds from gravitational infall. By magnetic braking, thermal cooling or diffusive magnetic processes like ambipolar diffusion a supercritical cloud can form. The subsequent collapse, which proceeds preferentially along the rotation axis and the magnetic field lines, results in a dense central object, i.e. a protostar and an accretion disc (Mestel 1965; Shu et al. 1993; see Tomisaka 1995 for a numerical study). After the first, dynamical collapse we expect the disc around the central object to be still threaded by a fraction of the interstellar magnetic flux of the original cloud.

Though microscopic diffusion processes are important for the evolution of magnetic fields in star forming regions, effective diffusion of matter across field lines is also possible through instabilities driven by the magnetic field itself. In the molecular cloud cores where the magnetic field is thought to disengage itself from the matter, the gravitational and magnetic energy densities are believed to be of comparable magnitude (see McKee et al. 1993 for a review). This is just the condition under which magnetic instabilities can operate at rates competitive with the gravitational collapse rate. It would be somewhat of a coincidence if microscopic diffusion were always the dominant process, and effective diffusion by magnetically-driven instabilities never played a role in the entire contraction process from cloud to star. This is one important reason to study the possible effects of non–axisymmetric magnetic instabilities in the process of star formation.

Another reason is provided by the possible connection of magnetic fields in discs with outflows and jets. A remnant of the cloud core’s magnetic flux, anchored in the accretion disc, represents a poloidal magnetic field component with the right geometry to launch and accelerate disc gas along the magnetic field lines (Bisnovatyi-Kogan and Ruzmaikin 1976; Blandford & Payne 1982). P Cygni line profiles in the spectra of Young Stellar Objects might find their explanation in such disc winds (Edwards, Ray & Mundt 1993).

In addition to their role in accelerating outflows, poloidal magnetic fields anchored in the accretion disc may also play an important role in the collimation of bipolar outflows and jets (Blandford 1993; Spruit, Fogglizio & Stehle 1997). We refer to Hughes (1991) for an extensive introduction into the observations and physics of beams and jets. For a review and tutorial introduction to the magnetic acceleration model see Spruit 1996; for some issues of current interest see Ogilvie and Livio (1998), Cao and Spruit (1994, 2000).

The ability of the disc to produce a magnetically accelerated outflow depends rather critically on the strength and radial distribution of the poloidal field at its surface. This distribution, in turn, is determined by the rate at which an advected large scale magnetic field in the disc is able to diffuse outward against the inward drift velocity compressing it.

Van Ballegooijen (1989) studied the radial transport of magnetic field lines assuming an isotropic turbulent viscosity $`\nu `$, related to an equally isotropic magnetic diffusivity $`\eta `$ by a constant ratio, the magnetic Prandtl number. Along the same lines Lubow, Papaloizou & Pringle (1994) studied the magnetic field dragging by turbulent processes in an accretion disc which was initially threaded by an externally imposed magnetic field. The conclusion from these studies is that inward dragging of an external field is difficult to achieve, if the assumption of an isotropic magnetic Prandtl number of order unity holds. It follows that one should expect to see, in this case, only internally generated fields like those obtained in numerical simulations (Stone et al. 1996; Brandenburg et al. 1995).

If the initial magnetic field is sufficiently strong to contribute to support against gravity, however, as appears to be the case in cloud cores, it is likely to suppress these dynamo processes, since a weak-field instability (Balbus & Hawley 1992) is thought to be an essential ingredient in this kind of dynamo. In this case, which we shall call here the strong field case, the most plausible source of turbulence in the disc are non–axisymmetric instabilities caused by the strong magnetic field itself.

One well-known instability is interchange, a local instability driven by the magnetic field energy, and which plays an important role in controlled fusion devices and in solar magnetic fields. The example of solar prominences shows that these instabilities can be quite efficient in transporting mass across field lines (cf. Priest 1982). Linear interchange instability has been studied previously by Spruit & Taam (1990) for uniformly rotating discs, and for differentially rotating discs by Spruit, Stehle & Papaloizou (1996) and Lubow & Spruit (1996), in a shearing sheet approximation. The instability appears when the ratio $`B_z/\mathrm{\Sigma }`$ of the vertical magnetic field strength to the surface mass density decreases with distance from the centre of the disc. Its behavior is similar to convection, with the gradient of $`B_z/\mathrm{\Sigma }`$ replacing the entropy gradient.

As in the case of solar prominences, it is the magnetic tension force due to the curvature of the field lines in the $`rz`$-plane that drives the instability in geometrically thin discs (Anzer 1967; Spruit and Taam 1992; for a detailed analysis of the magnetic forces in thin discs see Ogilvie 1997). Since non–axisymmetry is also crucial, this makes the magnetic field three-dimensional. Our study thus differs conceptually from the numerical study by Kaisig Tajima and Lovelace (1992), where a 2-dimensional cylindrical configuration with vertical magnetic field lines was assumed. In that case, the only magnetic force is the magnetic pressure gradient.

A full three-dimensional numerical treatment of the problem is made difficult by the very large range of characteristic speeds expected in the problem. Inside the disc, the Alvfén speed is not larger than the sound speed, but in the low density regions outside the disc it can easily approach the speed of light. We show here how this difficulty can be turned into an advantage, such that only a part of the problem needs to be done in 3 dimensions, while the numerically most demanding parts can be done using only the two dimensions in the disc plane.

In the above, we have introduced the magnetic disc model with the example of protostellar discs. The physics addressed by our calculations, however, is equally applicable to discs in compact binaries or active galactic nuclei.

### 1.1 poloidal fields vs. dynamo-generated fields

The magnetic problem we study is in several ways complementary to that of the simulations done by Stone et al. (1996) or Brandenburg et al. (1995), and Hawley (2000). There the magnetic fields are generated locally by a dynamo process, which can already start from a very weak intitial field. If a strong field is present initially, such that the magnetic energy density is comparable to the thermal energy density, the dynamo process is suppressed. This limit can be written in the form

$$\frac{B^2}{4\pi \mathrm{\Sigma }}c_\mathrm{s}\mathrm{\Omega },$$
(1)
where $`\mathrm{\Sigma }`$ is the surface mass density and $`c_\mathrm{s}`$ the sound speed.

A small scale dynamo process like the magnetic turbulence seen in these simulations does not create a net flux of field lines through the disc (Hawley 2000), and the overall field structure is therefore at least of quadrupole order at large distance from the disc. For advected magnetic field lines, as they are advected radially and anchored in the disc, the total magnetic flux through the disc is non–zero. The global structure of the advected magnetic field, as seen from large distances, is thus close to a dipole magnetic field distribution. The radial force exerted by such a field is predominantly the tension force, whose magnitude integrated over the disc thickness is $`B_rB_z/2\pi `$. A poloidal field can in principle exist up to strengths such that this radial force starts contributing significantly to the support against gravity. This limit can be written as

$$\frac{B^2}{4\pi \mathrm{\Sigma }}\mathrm{\Omega }^2r.$$
(2)
This limit on $`B^2`$ is a factor $`\mathrm{\Omega }r/c_\mathrm{s}`$ larger than the strength (1) at which dynamo-generated fields are suppressed. Once strong poloidal fields exist in a disc, they suppress the magnetic turbulence that would make them diffuse out of the disc by van Ballegooijen’s argument. There is thus a large range in parameter space where a poloidal field in a disc would be affected only by its own internal instabilities. These are the subject of the present investigation.

## 2 Two-dimensional discs with three-dimensional magnetospheres

We neglect the self–gravity of the disc so that the central star is the only cause of gravitational force acting on the disc. Following Spruit & Taam (1990, hereafter ST, see also Tagger et al. 1990 for a similar discussion) we assume a geometrically thin disc with an internal velocity field that does not depend on the vertical coordinate (perpendicular to the disc). The equations of motion can then be integrated across the disc, resulting in a two-dimensional problem confined to the plane of the disc. This plane can in principle be of arbitrary time-dependent shape, as in ST, but we limit the calculations here to the case of a plane disc, without displacements of the disc surface in the vertical direction (i.e. without corrugations or ‘bending modes’. For results on such modes, in the same approximation for the magnetic field, see Agapitou et al. 1997).

One may wonder when the approximation of height-independent velocity fields in the disc is justified. The assumption clearly eliminates the possibility of dynamo-generation of magnetic fields, but is appropriate for strong fields in the sense discussed in the preceding section, namely fields which contribute, more than the pressure force, to support against gravity. For such strong fields, winding-up of field lines inside the disc can be ignored, because the small amount of differential rotation encountered by the field line as it crosses from one side of the disc to the other is not enough to bend field lines significantly. Instead, the differential rotation across the disc will adjust, due to the magnetic forces, such that the rotation rate is constant along a field line. In the thin-disc limit, the amount of differential rotation across the disc vanishes as $`H/r`$, justifying the assumption made (for a discussion of the thin limit of magnetic discs, from a more mathematical point of view, see Ogilvie 1997).

Because of the thin disc limit taken, the dominant magnetic force in the problem is the curvature force due to the bend of the field lines crossing the disc (see ST). Outside the disc, the magnetic forces dominate over fluid forces on account of the low gas density. We call this region of low plasma-$`\beta `$ the ‘magnetosphere’ of the disc (not to be confused with the magnetosphere of an accreting magnetic star). We simplify the physics by assuming the field in this magnetosphere to be current free so that it is derivable from a magnetic potential $`\psi `$. This is justified if the matter density in the magnetosphere is sufficiently low. The magnetospheric field in then approxmately force free. The case of the solar corona shows that such a field is in practice also close to a potential (current free) field. This is due to the fact that at low $`\beta `$, ‘forced’ processes (e.g. Parker 1979, for recent numerical simulations see Galsgaard and Nordlund 1997) are fast, and keep the degree of twisting in the field low.

In this approximation, all currents are confined to the plane of the disc, and are proportional to the jump in the tangential field components across the disc. These tangential components ($`B_\varphi ,B_r`$) are determined by the normal field component through the solution of the (3-dimensional) potential problem in the magnetosphere, for which the normal component $`B_z`$ provides the boundary condition.

The magnetic forces acting on the fluid in the disc plane are given by the difference in the magnetic stress acting on the upper and lower surfaces of the disc. They are proportional to the product of the tangential and normal components of the magnetic field at the disc surface (see ST for details).

In the computations, closed inner and outer boundaries are used for the disc ($`v_r=0`$, so that no mass or magnetic flux crosses these boundaries. Thus the total magnetic flux through the disc and the total disc mass are constant in time.

Apart from the addition of the magnetic force term, the hydrodynamical problem is the same as in ordinary two-dimensional disc hydrodynamics. We use an Eulerian grid with the van Leer (1977) scheme for upwind differencing. An outline of the model and its basic assumptions has been given previously in Stehle (1997).

## 3 Equations

To describe the non–axisymmetric disc response resulting from large scale magnetic fields, we adopt a thin–disc approximation in which the vertical velocity vanishes. We use a cylindrical coordinate system ($`r,\varphi ,z`$) defined with the origin at the central mass $`M`$. The mass surface–density is $`\mathrm{\Sigma }=_{\mathrm{}}^+\mathrm{}\rho dz`$ with $`\rho =\rho (r,\varphi ,z)`$ the gas density.

In the thin disc limit $`H/r1`$ where $`H`$ is the disc thickness the only contribution of the Lorentz force per unit surface area is due the magnetic tension of the field lines (ST). This is because for magnetic fields varying on a length scale $`LH`$, the curvature force is larger than the magnetic pressure gradient by a factor $`L/H`$. Then the radial component of the equation of motion for the disc, assumed to be inviscid, reads

$`{\displaystyle \frac{v_\mathrm{r}}{t}}`$ $`+`$ $`v_\mathrm{r}{\displaystyle \frac{v_\mathrm{r}}{r}}+{\displaystyle \frac{v_\varphi }{r}}{\displaystyle \frac{v_\mathrm{r}}{\varphi }}{\displaystyle \frac{v_\varphi ^2}{r}}`$
$`=`$ $`{\displaystyle \frac{1}{\mathrm{\Sigma }}}{\displaystyle \frac{P}{r}}+{\displaystyle \frac{B_\mathrm{z}[B_\mathrm{r}]}{4\pi \mathrm{\Sigma }}}g,`$
where $`v_\mathrm{r}`$ and $`v_\varphi =r\mathrm{\Omega }`$ are the radial and azimuthal components of the velocity vector and $`g`$is the acceleration of gravity due to the central star. $`B_\mathrm{z}`$ is the vertical magnetic field component at disc midplane and $`[B_\mathrm{r}]=B_\mathrm{r}^+B_\mathrm{r}^{}=2B_\mathrm{r}^+`$ the jump of the field vector from above ($`B_\mathrm{r}^+`$) to below ($`B_\mathrm{r}^{}`$) the disc plane. Neglecting any possible disc warps it is assumed that the field geometry is antisymmetric with respect to the disc plane (i.e. $`B_\mathrm{z}^+=B_\mathrm{z}^{}`$ and $`B_\mathrm{r}^+=B_\mathrm{r}^{}`$, see ST). For use in what follows we define the magnetic acceleration in radial direction,

$$g_\mathrm{m}=B_\mathrm{z}B_\mathrm{r}^+/2\pi \mathrm{\Sigma }.$$
(4)
The azimuthal equation of motion is

$`{\displaystyle \frac{v_\varphi }{t}}`$ $`+`$ $`v_\mathrm{r}{\displaystyle \frac{v_\varphi }{r}}+{\displaystyle \frac{v_\varphi }{r}}{\displaystyle \frac{v_\varphi }{\varphi }}+{\displaystyle \frac{v_\mathrm{r}v_\varphi }{r}}`$
$`=`$ $`{\displaystyle \frac{1}{\mathrm{\Sigma }}}{\displaystyle \frac{P}{r\varphi }}+{\displaystyle \frac{B_\mathrm{z}[B_\varphi ]}{4\pi \mathrm{\Sigma }}}.`$
The continuity equation is

$$\frac{\mathrm{\Sigma }}{t}+\frac{1}{r}\frac{}{r}(\mathrm{\Sigma }rv_\mathrm{r})+\frac{1}{r^2}\frac{}{\varphi }(\mathrm{\Sigma }rv_\varphi )=0.$$
(6)
A similar equation holds for the vertical component of the magnetic field in the limit of complete flux freezing:

$$\frac{B_\mathrm{z}}{t}+\frac{1}{r}\frac{}{r}(B_\mathrm{z}rv_\mathrm{r})+\frac{1}{r^2}\frac{}{\varphi }(B_\mathrm{z}rv_\varphi )=0,$$
(7)
which expresses that the magnetic flux density $`B_z`$ is conserved.

The gas pressure is computed from the vertically integrated internal energy $`e`$ assuming an ideal gas for the equation of state, $`P=(\gamma 1)e`$, with ratio of specific heats $`\gamma `$. For the computations reported here, the value $`\gamma =1.4`$ was used. The adiabatic evolution of the internal energy is given by

$`{\displaystyle \frac{e}{t}}`$ $`+`$ $`{\displaystyle \frac{1}{r}}{\displaystyle \frac{}{r}}(rev_\mathrm{r})+{\displaystyle \frac{1}{r^2}}{\displaystyle \frac{}{\varphi }}(rev_\varphi )`$
$`=`$ $`P\left[{\displaystyle \frac{1}{r}}{\displaystyle \frac{}{r}}(rv_\mathrm{r})+{\displaystyle \frac{}{r\varphi }}(v_\varphi )\right].`$
To close our set of equations we have to determine the inclination of the magnetic field lines to the surface of the accretion disc, i.e. we need to derive $`B_\mathrm{r}^+(r,\varphi )`$ and $`B_\varphi ^+(r,\varphi )`$. Their values determine the magnetic forces in the equations of motion. By our assumption of a potential field in the magnetosphere, electric currents only exist within the plane of the accretion disc, and the vertical component of the electric current vanishes everywhere.

The current is therefore of the form

$$𝐣_\mathrm{h}(r,\varphi ,z)=𝐣_\mathrm{h}(r,\varphi )\delta (z)\text{and}j_\mathrm{z}(r,\varphi ,z)=0$$
(9)
where $`\delta (z)`$ is the Dirac–delta function. A subscript h denotes vectors parallel to the plane of the accretion disc. Disk winds and ionized particles in the disc magnetosphere, neglected in our model, will certainly contribute to current flows in the disc magnetosphere. The calculation of the 3–dimensional disc magnetosphere in force-free approximation is, however, beyond the present computational feasibility. The assumption of a potential field for the structure of the disc magnetosphere is equivalent to assuming that the magnetosphere is sufficiently close to its minimum energy state.

In the thin disc limit, the magnetic field at the disc surface $`(B_\mathrm{r}^+,B_\varphi ^+,B_\mathrm{z})`$ is connected to the disc currents by

$$B_\mathrm{r}^+=\frac{2\pi }{c}j_\varphi ,B_\varphi ^+=\frac{2\pi }{c}j_\mathrm{r}$$
(10)
and

$`B_\mathrm{z}(r,\varphi )`$ $`=`$ $`{\displaystyle \frac{1}{c}}{\displaystyle _0^{2\pi }}{\displaystyle _{r_{\mathrm{in}}}^{r_{\mathrm{out}}}}{\displaystyle \frac{_\mathrm{r}^{}(r^{}j_\varphi )_\varphi ^{}j_\mathrm{r}}{R}}dr^{}d\varphi ^{}`$
$``$ $`{\displaystyle \frac{1}{c}}{\displaystyle _{\delta F_{disc}}}{\displaystyle \frac{1}{R}}𝒋\mathbf{(}𝒓^{\mathbf{}}\mathbf{)}d𝒍^{\mathbf{}}`$
where the second term on the right hand side (RHS) of eq. (3) is the current flow at the disc–boundaries and $`R`$ is given by

$$R^2=r^2+r^22rr^{}\mathrm{cos}(\varphi \varphi ^{}).$$
(12)
Eq. (3) is derived from a partial integration of the vector potential field $`𝑨`$ in Coulomb gauge so that $`𝑩=\mathbf{}\times 𝑨`$, and (Landau & Lifshitz 1975):

$$𝑨(r,\varphi )=\frac{1}{c}_0^{2\pi }_{r_{\mathrm{in}}}^{r_{\mathrm{out}}}\frac{𝒋(r^{},\varphi ^{})}{R}dr^{}d\varphi ^{}.$$
(13)
Given the magnetic flux distribution in the disc $`B_\mathrm{z}(r,\varphi )`$ eq. (3) has to be inverted to yield the currents $`𝒋_\mathrm{h}`$ in the disc. The inversion is unique by applying the fact that currents do not accumulate, $`\mathrm{div}𝒋=0`$:

$$\frac{1}{r}\frac{}{r}(rj_\mathrm{r})+\frac{1}{r}\frac{}{\varphi }j_\varphi =0.$$
(14)
This connects the radial $`j_\mathrm{r}`$ and azimuthal $`j_\varphi `$ component of the current, and has to be solved together with the inversion of eq. (3).

As our accretion disc has a central hole, the solution of the potential field problem with $`B_\mathrm{z}(r,\varphi )`$ given in the disc is no longer unique, since the domain space is multiply connected. To any solution of the inversion problem (3), an arbitrary multiple of the solution for $`B_z=0`$ can be added. This special solution consists of closed field lines wrapping through the hole of the disc and around its outer edge without crossing the disc itself (like the windings of a ring-core coil). We thus have additionally to specify the number of field lines which pass through the central hole, i.e. the total magnetic flux through the hole of the disc (see also Lubow, Papaloizou & Pringle 1994)

$$\mathrm{\Psi }=_0^{2\pi }_0^{\mathrm{r}_{\mathrm{in}}}rB_\mathrm{z}(r,\varphi ,z=0)drd\varphi .$$
(15)
A Fourier Transform of eq. (15) shows that only the axisymmetric component contributes to the magnetic flux through the disc hole whereas all other components cancel. The degree of freedom introduced by the presence of a hole thus enters only into the computation of the axisymmetric component of the field.

Eqs. (6), (7) and (3) specify the time evolution of $`\mathrm{\Sigma }`$, $`B_\mathrm{z}`$ and $`e`$ when the velocities $`v_\mathrm{r}`$ and $`v_\varphi `$ in the plane of the accretion disc are known. These follow from Euler’s equation (3) and (3) which includes thermal, gravitational and magnetic forces. To solve for the magnetic forces we invert equation (3) with the differential constraint (14) for the unknown currents, and an assumed value for the magnetic flux though the hole. Eq. (10) gives the relation between $`B_\mathrm{r}^+`$, $`B_\varphi ^+`$ and $`j_\varphi `$, $`j_\mathrm{r}`$ and thus the magnetic forces are determined. In the next chapter we show how we solve the hydrodynamic equations and the magnetospheric potential problem numerically.

## 4 Numerical Method

We solve the equations on an Eulerian grid with equidistant spacing in the $`r`$ and $`\varphi `$ coordinates. The inner disc rim is at $`r_{\mathrm{in}}=0.1r_{\mathrm{out}}`$ with $`r_{\mathrm{out}}`$ the radius of the outer edge of the disc. The number of grid points and grid spacing in radial direction are $`n_r`$ and $`\mathrm{\Delta }r=(r_{\mathrm{out}}r_{\mathrm{in}})/n_\mathrm{r}`$. We use a staggered grid such that the scalar quantities $`\mathrm{\Sigma },B_\mathrm{z},e`$ are defined at the cell centres and the vector quantities $`(v_\mathrm{r},v_\varphi )`$, $`(B_\mathrm{r},B_\varphi )`$ and $`(j_\varphi ,j_\mathrm{r})`$ at the cell boundaries. The equations are used in dimensionless form, as follows. As unit of length we take the outer disc radius, as the unit of time the inverse of the Keplerian rotation rate at the outer edge. Thus, the Keplerian velocity at the outer edge is unity, and the time for one Keplerian orbital period of the outer edge is $`T=2\pi `$. For all the calculations presented here we use $`n_\mathrm{r}\times n_\varphi =156\times 128`$ grid cells in $`r`$ and $`\varphi `$ direction respectively. We also performed some calculations on a smaller grid of $`64\times 64`$ grid cells. Comparing these models to calculations done with a higher resolution we find only differences on the order of the applied grid spacing. We are thus convinced that the models presented here are resolved sufficiently.

### 4.1 The hydrodynamic part

The hydrodynamic part of the calculations is done with a natural extension to the scheme described in Stehle & Spruit (1999).

The equations are written in conservative form. Terms in the equations are divided into advection and source terms. The advection from one grid cell to the other is done with the upwind–differencing scheme of van Leer (1977). Thermal pressure, gravitational forces and compressional heating are calculated following Stone & Norman (1992). Viscous processes and radiative cooling from the surface of the accretion disc are neglected.

The induction equation, i.e. eq. (7), is an additional equation compared with the hydrodynamic case. It has the same forms as the continuity equation, and is treated numerically in the same way.

The magnetic force is a new term in the equations of motion. Contrary to the pressure force, which is calculated by a local derivative, the magnetic forces are given by the bending of the magnetic fields where they pass through the accretion disc. The inclination of the magnetic field lines to the accretion disc depends, however, on the global magnetic field structure in the magnetosphere and only when this is known, the magnetic forces be calculated. We introduce the numerical method to calculate the magnetic forces in Section (4.2).

The time step is controlled by the Courant–Friedrich–Levy condition. In addition to the sound speed, there is a magnetic wave speed in the problem. The wave is compressive, but since the restoring force results from the change in the external potential upon compression of the field lines rather than of the magnetic pressure itself, it is a dispersive wave. The phase velocity is (ST, Tagger et al. 1990)

$$\frac{\omega }{k}=\frac{B_\mathrm{z}}{\sqrt{2\pi \mathrm{\Sigma }k}},$$
(16)
where $`k=2\pi /\lambda `$ is the wavenumber and $`\lambda `$ the wavelength. The group speed, which carries the wave information, is a factor 2 lower. The highest magnetic wave speed in the discretized equations is therefore obtained for the highest wavenumber that can be represented by the grid. By the Nyquist theorem, this is $`k=\pi /\mathrm{\Delta }r`$. The magnetic wave speed $`v_\mathrm{m}`$ that enters the Courant condition is thus

$$v_\mathrm{m}=\frac{B_\mathrm{z}}{2\pi \sqrt{2\mathrm{\Sigma }/\mathrm{\Delta }r}}.$$
(17)
A Courant factor of 0.75 was found to be sufficient for numerical stability. In most cases, however, we find the time step to be controlled by the azimuthal velocity $`v_\varphi (r_{\mathrm{in}})`$ at the inner accretion disc rim. The additional magnetic wave constrains the time step only for models where the magnetic acceleration $`g_\mathrm{m}`$ approaches the acceleration of gravity $`g`$.

The boundary conditions used are solid boundaries at the inner and outer edges of the disc, i.e. $`v_r=P/r=0`$. These are sufficient for the present calculations in which only the short-term evolution of the disc is followed. For the longer-term evolution, one would want to use conditions that allow accretion to take place through the boundaries. This is beyond the scope of the present study.

### 4.2 The solver for the disc magnetosphere

The magnetic forces are calculated at each time step from the magnetic flux distribution $`B_z(r,\varphi )`$ in the disc. This involves two steps. First, eq. (3) is inverted to obtain the currents $`𝐣_\mathrm{h}`$ from the flux distribution $`B_z`$. In the second step, the forces $`B_zB_r^+/2\pi `$, $`B_zB_\varphi ^+/2\pi `$ are computed using eqs. (10), and added to the hydrodynamic forces.

The inversion of eq. (3) has two complications. Because of charge conservation the current is not an arbitrary function of $`r`$ and $`\varphi `$, but must satisfy eq. (14), i.e. $`\mathrm{div}𝐣=0`$. This condition can be used to eliminate one of the current components $`j_r`$ and $`j_\varphi `$ in favor of the other. Eq. (3) can then be read as an integral equation determining one of the current components in terms of the vertical field component.

Second, the region on which the computations are done is not simply connected. The central hole in the computational domain generates an additional parameter in the potential problem, namely the net magnetic flux through the hole $`\mathrm{\Psi }`$. For the present exploratory calculations, we use fixed boundary conditions. Hence $`v_r=0`$ at the boundaries, and with the induction equation (7) no magnetic flux enters or leaves through the edges of the disc. This also implies that $`\mathrm{\Psi }`$ is constant in time. More general boundary conditions are possible that take account of the advection of magnetic flux into the hole. These will be needed in more realistic calculations of discs accreting on to protostars, for example.

#### 4.2.1 Fourier decomposition in $`\varphi `$

Since the equations for the potential problem are homogeneous in the azimuth $`\varphi `$, Fourier transforms can be used in this direction. Since they are also linear, the Fourier components do not mix, and one can solve for each of the Fourier components separately. If the number of azimuthal grid points or Fourier modes is $`n_\varphi `$, this reduces the computing effort required by a factor of the order $`n_\varphi `$ compared with straightforward discretization in $`\varphi `$. By using Fourier decomposition, the computing effort for the potential problem can be kept at a level comparable to the hydrodynamic parts of the calculation.

Thus we write the magnetic flux distribution as

$$B_\mathrm{z}(r,\varphi )=B_\mathrm{z}^0+\underset{m=1}{\overset{\mathrm{n}_\varphi /2}{}}\left(B_\mathrm{z}^{m,\mathrm{s}}\mathrm{sin}(m\varphi )+B_\mathrm{z}^{m,\mathrm{c}}\mathrm{cos}(m\varphi )\right),$$
(18)
where $`B_\mathrm{z}^0(r)`$, $`B_\mathrm{z}^{m,\mathrm{s}}(r)`$ and $`B_\mathrm{z}^{m,\mathrm{c}}(r)`$ are only functions of $`r`$. Similar equations hold for the currents $`j_r,j_\varphi `$.

The Fourier amplitudes are given by

$$B_\mathrm{z}^{m,\mathrm{c}}(r)=\frac{1}{\pi }_0^{2\pi }B_\mathrm{z}(r,\varphi )\mathrm{cos}(m\varphi )d\varphi ,$$
(19)
and similarly for $`B_z^{m,\mathrm{s}}`$. Using (3) this becomes

$`B_\mathrm{z}^{m,\mathrm{c}}(r)`$ $`=`$ $`{\displaystyle \frac{1}{\pi c}}{\displaystyle _{r_{\mathrm{i}n}}^{r_{\mathrm{o}ut}}}{\displaystyle _0^{2\pi }}{\displaystyle _0^{2\pi }}dr^{}d\varphi d\varphi ^{}`$
$`\mathrm{cos}(m\varphi ){\displaystyle \frac{_r^{}[r^{}j_\varphi (r^{},\varphi ^{})]_\varphi ^{}j_r(r^{},\varphi ^{})}{(r^2+r_{}^{}{}_{}{}^{2}2rr^{}\mathrm{cos}(\varphi \varphi ^{}))^{1/2}}}.`$
This can be written as

$`B_\mathrm{z}^{m,c}(r)`$ $`=`$ $`{\displaystyle \frac{1}{\pi cr}}{\displaystyle _{r_{\mathrm{i}n}}^{r_{\mathrm{o}ut}}}{\displaystyle _0^{2\pi }}dr^{}d\varphi ^{}`$
$`\mathrm{cos}(m\varphi ^{})[_r^{}(r^{}j_\varphi )_\varphi ^{}j_r]K_m(r^{}/r),`$
where

$$K_m(x)=_0^{2\pi }\frac{\mathrm{cos}(m\varphi )}{(1+x^22x\mathrm{cos}(\varphi ^{}))^{1/2}}d\varphi ^{}.$$
(22)
We evaluate this function numerically. For $`m=0`$ it can be expressed in terms of the complete elliptical function of the second kind $`F(x)`$ (Gradstein & Ryzhik 1981),

$$K_0(x)=4\chi F(\chi )\text{with}\chi =\mathrm{min}(x,1/x).$$
(23)
Substituting the Fourier expansions of $`j_\mathrm{r}`$ and $`j_\varphi `$ and integrating over $`\varphi ^{}`$, we get

$$B_z^{m,c}(r)=\frac{1}{cr}_{r_{\mathrm{i}n}}^{r_{\mathrm{o}ut}}[_r^{}(r^{}j_\varphi ^{m,\mathrm{c}})mj_r^{m,\mathrm{s}}]K_m(r^{}/r)dr^{}.$$
(24)
For $`m0`$ the charge conservation condition (14) can be used to eliminate $`j_\varphi `$, which yields

$`B_z^{m,c}(r)`$ $`=`$ $`^{m,c}(r,j_r)={\displaystyle \frac{1}{cr}}{\displaystyle _{r_{\mathrm{i}n}}^{r_{\mathrm{o}ut}}}dr^{}`$
$`\left[{\displaystyle \frac{1}{m}}_r^{}(r^{}_r^{}r^{}j_r^{m,s})m_\varphi ^{}j_r^{m,s}\right]K_m(r^{}/r),`$
where $`^{m,\mathrm{c}}(r,j_\mathrm{r})`$ is an abbreviation for the integral operator. A similar equation holds for the $`B_\mathrm{z}^{m,\mathrm{s}}`$ component.

The axially symmetric component reads

$$B_z^0(r)=^0(r,j_\varphi )=\frac{1}{c}_{r_{\mathrm{i}n}}^{r_{\mathrm{o}ut}}_r^{}(r^{}j_\varphi ^0)K_0(r^{}/r)\mathrm{d}r^{}.$$
(26)
The flux through the hole in the middle of the disc is

$`\mathrm{\Psi }`$ $`=`$ $`_\mathrm{\Psi }(j_\varphi )=2\pi {\displaystyle _0^{r_{\mathrm{i}n}}}rB_\mathrm{z}^0dr=`$
$`=`$ $`{\displaystyle \frac{2\pi }{c}}{\displaystyle _0^{r_{\mathrm{in}}}}{\displaystyle _{r_{\mathrm{i}n}}^{r_{\mathrm{o}ut}}}_r^{}(r^{}j_\varphi ^0)K_0(r^{}/r)\mathrm{d}r^{}\mathrm{d}r.`$
Eqs. (4.2.1) and (26) are integral equations for the $`(n_\varphi 1)`$ Fourier amplitudes $`j_r^{m,\mathrm{c}}(r)`$ and $`j_r^{m,\mathrm{s}}(r)`$ of the radial current, and the azimuthal current distribution $`j_\varphi ^0(r)`$. To solve these, we take finite differences in $`r`$, which turns each of the eqs. (4.2.1) and (26) into a set of linear algebraic equations. The matrices involved are fixed in time and need to be inverted only once. If the number of grid points in $`r`$ is $`n_\mathrm{r}`$, the computing effort for the potential problem is therefore of the order $`n_\mathrm{r}^2n_\varphi `$ per time step, or $`n_\mathrm{r}`$ per grid point and time step. Since the number of operations per grid point and time step in the hydrodynamic part of the calculation is a substantial, but fixed number independent of $`n_\mathrm{r}`$, the computing expense for the potential problem does not dominate the overall expense except for very large numbers of radial grid points. In fact, for a grid of $`n_\mathrm{r}\times n_\varphi =156\times 128`$ grid cells the magnetic part of the calculation takes about 50% of the CPU–time.

#### 4.2.2 Discretization in $`r`$

The discretization in $`r`$ of the integrals in (4.2.1) & (26) is different for the non–axisymmetric components (4.2.1) and the axisymmetric component (26). The Fourier amplitudes $`B_\mathrm{z}^m`$ of the field are naturally defined at the same radial positions as the values of the field $`B_\mathrm{z}`$ itself, i.e. at the centres $`r_i`$ of the cells defined in the hydrodynamic part of the scheme. To discretize (26), the currents $`j_\varphi ^0`$ are defined at the boundaries $`r_{i+1/2}(r_i+r_{i+1})/2`$ between these cells. Because of the topology of the domain used (a disc with a hole), there is one more of such boundaries than there are cell centres. Since the boundary conditions do not impose constraints on the azimuthal current at the boundaries, there is then also one more current $`j_\varphi ^0`$ than there are field amplitudes $`B^0`$. This additional degree of freedom is balanced by the the hole-flux condition (15), which arises from the same topological property.

#### 4.2.3 The axisymmetric component

To evaluate the integral in (26) we interpolate the current linearly between the values $`J_j^0`$ at its grid points:

$$j_{\varphi ,j}^0(r)=(x1)J_j^0+xJ_{j+1}^0,$$
(28)
where $`x=(rr_{j+1/2})/(r_{j+3/2}r_{j+1/2})`$. Inserting this into (26), $`B_z^0`$ is a linear function of the currents $`J_j^0`$, with coefficients $`B_{ij}^0`$:

$$B_z^0(r_i)=\underset{j}{}B_{ij}^0J_j^0.$$
(29)
Evaluation of these coefficients involves integrals of the type

$$I^0(r)=_{r_\mathrm{a}}^{r_\mathrm{b}}(\alpha +\beta r^{})K_0(r^{}/r_i)dr^{}.$$
(30)
The elliptic function involved in $`K_0`$ is evaluated by polynomial approximation (Abramowitz and Stegun 1984), the integral by Bulirsch’s algorithm (Press et al. 1995).

Similarly, assigning the index $`i=0`$ to the hole-flux $`\mathrm{\Psi }`$, (4.2.1) can be written as

$$\mathrm{\Psi }=\underset{j}{}B_{0j}^0J_j^0.$$
(31)
The coefficients $`B_{ij}^0`$ then form a square matrix of dimension $`n_r+1`$, relating a vector of magnetic variables, $`\mathrm{b}^0=(\mathrm{\Psi },B_\mathrm{z}^0(r_1),\mathrm{}B_\mathrm{z}^0(r_{n_r}))`$, to the currents:

$$b_i^0=\underset{j=0}{\overset{n_r}{}}B_{ij}^0J_j^0.$$
(32)
Since we use a fixed Eulerian grid, the matrix elements $`B_{ij}^0`$ are fixed and the inversion of the matrix can be done once and for all for a given computational grid. This inversion is done by the LU decomposition method (Press et al. 1995).

#### 4.2.4 The non–axisymmetric components

The procedure for the non–axisymmetric components is very similar to that for the axisymmetric component, except that $`\mathrm{\Psi }`$ does not appear because it is already determined by the axisymmetric component. For $`m0`$ the currents $`J_\mathrm{r}^m`$, like the $`B_\mathrm{z}^m`$, are defined at the cell centres, so that there is an equal number of each.

Because we have used a partial integration in deriving (4.2.1), and have used charge conservation to eliminate $`j_\varphi `$, the integrand in (4.2.1) contains a second derivative with respect to $`r`$. In order to evaluate it at the same order of accuracy as the axisymmetric coefficients, a third order interpolation is needed. We choose cubic spline interpolation between neighbouring grid–cells. The coefficients then involve expressions of the type

$$𝒦_m^l=_0^xx^lK_m(x)dx,$$
(33)
with $`l=1,2,3`$. Coefficients up to $`m=64`$ were evaluated with an accuracy of $`10^7`$.

#### 4.2.5 Test calculations

To test the accuracy of our numerical solution of the potential field problem, we apply the inversion of eq. (3) to a known field configuration. We proceed by first specifying $`𝒋_\mathrm{i}(r)`$ analytically. Then we derive $`B_\mathrm{z}(r)`$ from a numerical integration of eq. (3). The integration yields $`B_\mathrm{z}(r_\mathrm{i})`$ at the grid cell centres $`r_\mathrm{i}`$ with an accuracy of $`10^6`$. These values are used to invert $`B_\mathrm{z}`$ with our magnetosphere–solver to reconstruct the disc currents $`𝒋_\mathrm{f}(r_{\mathrm{i}+1/2})`$ and field strengths which are then compared with the original function values at these points.

An example is shown in Fig. (1). This is an $`m=1`$ mode with the current distribution

$$j_{\mathrm{r},1}^{\mathrm{m}=1}=(rr_{\mathrm{in}})^2(r_{\mathrm{out}}r)^2$$
(34)
with $`r_{\mathrm{in}}=0.01r_{\mathrm{out}}`$. The reconstructed current distribution matches the input distribution within an accuracy of 1%, for the 100-point grid used.

Additional tests of accuracy of the code as a whole were done by comparing results at different resolutions.

## 5 Uniformly rotating discs

The linear stability of uniformly rotating discs with magnetic fields of the type considered here was studied by ST. The condition for instability of the interchange type is

$$a\frac{B_\mathrm{r}^+B_\mathrm{z}}{2\pi \mathrm{\Sigma }\mathrm{\Omega }^2}\frac{\mathrm{d}}{\mathrm{d}r}\mathrm{ln}\left(\frac{B_\mathrm{z}}{\mathrm{\Sigma }}\right)<0.$$
(35)
Disks where $`a0`$ everywhere are stable to interchange instability. The derivation of this condition does not take account of possible global instabilities.

In this section we study two examples, one of a disc which is stable (model 1) and one that is unstable (model 2) according to condition (35).

Uniformly rotating discs are set up numerically by modifying the gravitational potential $`\mathrm{\Phi }(r)`$ such that magnetic, centrifugal and gravitational forces are just balanced for the case $`v_\varphi (r)=`$const. While this case is of limited astrophysical interest, it serves to test the agreement with the predictions from the linear theory, to check for possible global modes not covered by condition (35), and to get an impression of the nonlinear development of interchange instability in the present case.

We study the discs in a frame of reference which corotates with the disc. The corresponding non–inertial terms are added to the equations in section 2 (see also Stehle & Spruit 1999). Both models presented here are advanced in time with zero magnetic flux through the central hole, i.e. with $`\mathrm{\Psi }=\text{const.}=0`$.

The uniformly rotating case as defined above is a one-parameter family, governed by the ratio $`c_\mathrm{u}`$ of magnetic to the centrifugal forces:

$$c_\mathrm{u}=\left(\frac{g_\mathrm{m}}{g}\right)_{\mathrm{max}}=\left(\frac{rB_\mathrm{z}B_\mathrm{r}}{2\pi \mathrm{\Sigma }v_\varphi ^2}\right).$$
(36)
Other parameters such as the amplitudes of the central potential, the magnetic field strength and the surface density can be scaled out of the equations. In cases where $`c_\mathrm{u}1`$ magnetic forces are unimportant and the disc rotates freely. In the case $`c_\mathrm{u}1`$ magnetic forces dominate and rotation can be neglected.

In the results shown, the gas pressure included in the calculations for numerical reasons (section 3) has only a small influence.

### 5.1 Model 1: $`\mathrm{d}\left(B_\mathrm{z}/\mathrm{\Sigma }\right)/\mathrm{d}r=0`$

According to eq. (35) discs with $`\mathrm{d}\left(B_\mathrm{z}/\mathrm{\Sigma }\right)/\mathrm{d}r=0`$ are expected to be stable against the interchange instability. The initial density and magnetic field distributions are specified as Gaussian humps in $`r`$ \[$`\mathrm{\Sigma }(r)\mathrm{exp}((r0.5r_{\mathrm{out}})^2/\mathrm{\Delta }^2)`$\], with a maximum at $`r=0.5r_{\mathrm{out}}`$ and width $`\mathrm{\Delta }=0.1r_{\mathrm{out}}`$.

Models with four different field strength are followed: a weakly magnetized disc ($`c_\mathrm{u}\mathrm{9\hspace{0.17em}10}^4`$), and one with a high magnetic support (model 1a, $`c_\mathrm{u}=225`$). Intermediate cases are chosen with $`c_\mathrm{u}=0.09`$ and $`c_\mathrm{u}=9`$ . We perturb the initial stationary, axisymmetric models with a low amplitude ($`10^7v_\varphi (r_{\mathrm{out}})`$) point–to–point noise in $`v_\mathrm{r}`$.

We find the extreme cases $`c_\mathrm{u}0`$ and $`c_\mathrm{u}\mathrm{}`$ to be stable. The total kinetic energy in the radial velocity component. (i.e. $`E_{\mathrm{kin},\mathrm{r}}=\mathrm{\Sigma }v_\mathrm{r}^2dF`$) was constant for the whole calculation of $``$15 disc orbits.

In the intermediate cases $`c_\mathrm{u}1`$ a very weak form instability was observed, with characteristics different from an interchange. The energy in the radial motions $`E_{\mathrm{kin},\mathrm{r}}`$ increased during the first 15 orbits by a factor $`100`$$`1000`$ . A global disc pattern is excited, showing spiral arms with $`m20`$. The wave saturates at a strength $`(B_\varphi /B_\mathrm{z})10^6`$$`10^4`$. The waves are of low amplitude and do not measurably transport angular momentum. After a linear rise time of some orbits the disc is found to be stationary again, i.e. $`_\mathrm{t}\mathrm{\Sigma }(r,\varphi )=0`$, even though it is now slightly non–axisymmetric.

The exact nature origin of this weak instability is not entirely clear at the moment. In any case, the amplitude of the motions observed is low compared to those of the instabilities described below, and are not relevant for actual accretion discs where the magnetic support is always less than gravity, $`c_\mathrm{u}<1`$.

### 5.2 Model 2: $`\mathrm{d}\left(B_\mathrm{z}/\mathrm{\Sigma }\right)/\mathrm{d}r<0`$

Next we study accretion discs where condition (35) predicts the interchange instability to be present.

We choose $`\mathrm{\Sigma }=`$const. and $`B_\mathrm{z}(r)`$ decreasing with $`r`$ as a power law

$$B_\mathrm{z}(r)=B_{\mathrm{z},0}\left(\frac{r}{r_{\mathrm{out}}}\right)^{4/5}.$$
(37)
With this choice, the instability parameter $`a`$ has a minimum at $`r0.2r_{\mathrm{out}}`$. We expect the instability to first manifest itself near this radius.

As in model 1 we perturb the initial axisymmetric model with a low amplitude point–to–point noise in $`v_\mathrm{r}`$.

Instability is found in all models of this sequence. Fig. (2) shows the time evolution of $`E_{\mathrm{kin},\mathrm{r}}\frac{1}{2}\mathrm{\Sigma }v_r^2dF`$ for these models. The model parameters for this sequence are shown in Table 1. At first $`E_{\mathrm{kin},\mathrm{r}}`$ increases exponentially. $`E_{\mathrm{kin},\mathrm{r}}`$ saturates after several orbits and the subsequent evolution is highly non–linear.

The instability causes a significant redistribution of the disc mass. This is illustrated in Fig. (3), which shows the evolution of the surface mass density for model 2c. It is seen that the instability first operates at $`r0.2r_{\mathrm{out}}`$ as expected from the local minimum of $`a`$. A high mode number $`m15`$ dominates at first. The pattern of motions resembles that of convective cells or the plumes of a Rayleigh-Taylor instability. The influence of rotation is evident in the asymmetry of the plumes. With time the instability is seen to spread over the whole disc. The small instability cells merge and grow in size as they drift to larger radii, as is characteristic of Rayleigh-Taylor instability.

We identify this instability with the interchange instability as the disc pattern looks very similar to convective cells as predicted by Spruit et al. (1995), and starts at the point where the linear analysis predicts the disc to be most unstable.

We conclude that the instability in the uniformly rotating case sets in as predicted from linear theory and has the nonlinear development of an interchange instability.

## 6 Non–uniformly rotating discs

We now study accretion discs revolving around the gravitational field of a point mass. According to the linear analysis of Spruit et al. (1995) and Lubow & Spruit (1995), shear due to differential rotation acts as a stabilizing factor on the interchange instability. This analysis predicts that instability appears only in regions of the disc where magnetic forces contribute significantly to support against gravity. The predicted linear growth is algebraic (a power law of time) rather than exponential.

We study the evolution of two different initial setups. First we follow a model where the magnetic field decreases as a power of radius (model 3) and then a case where it decreases exponentially (model 4). The initial structure of the models is summarized in Tab. (2) and (3). The initial $`v_\varphi (r)`$ is found from the radial force balance between magnetic, gravitational and centrifugal forces. We then perturb $`v_\mathrm{r}`$ by point–to–point low amplitude noise and follow the subsequent disc evolution numerically. As before all models are calculated with zero magnetic flux through the central hole of the disc, $`\mathrm{\Psi }=0`$.

### 6.1 Model 3: $`B_\mathrm{z}r^{5/4}`$ and $`\mathrm{\Sigma }r^{3/2}`$

For the initial state in model sequence 3we choose a surface density varying as $`\mathrm{\Sigma }r^{3/2}`$. In order to contain the effects of the instability within the computational domain as much as possible, we choose a magnetic field distribution with strength vanishing towards the boundaries (Fig. 4).

The degree of support against gravity by the magnetic field, as measured by the ratio $`g_\mathrm{m}/g_\mathrm{g}`$, increases with radius and has a maximum at $`0.85r_{\mathrm{out}}`$ (Fig. 4d). The parameters of the model sequence are summarized in Tab. (2). The parameter $``$ specifies the mean magnetic flux per surface–mass:

$$=\frac{_{\mathrm{disc}}B_\mathrm{z}dF}{_{\mathrm{disc}}\mathrm{\Sigma }dF}.$$
(38)
and can be used to compare model calculations.

Instead of interchange type instability, we find, in all three cases that the initial set up is unstable to a global, non–axisymmetric instability. The wave pattern of the instability can be traced from one edge of the disc to the other (see Fig. 7). Initially the kinetic energy $`E_{\mathrm{kin},\mathrm{r}}`$ growth exponentially on a dynamical time-scale. The initial growth rate, as given by $`\gamma _{\mathrm{Ekin},\mathrm{r}}=\mathrm{d}\mathrm{ln}E_{\mathrm{kin},r}/\mathrm{d}t`$, is largest for the disc with the highest magnetic support.

Fig. (5) shows the time evolution of the power

$$P_\mathrm{m}=\frac{_{r_{\mathrm{in}}}^{r_{\mathrm{out}}}\left((B_z^{m,\mathrm{c}})^2+(B_z^{m,\mathrm{s}})^2\right)^{1/2}rdr}{_{r_{\mathrm{in}}}^{r_{\mathrm{out}}}B_z^0rdr}$$
(39)
in Fourier mode $`m`$ of the field strength $`B_\mathrm{z}`$, integrated over the whole disc, and relative to the axisymmetric component $`m=0`$. In all three calculations we find that the modes $`m=1`$–4 grow equally fast. The relative power saturates nearly at the same level, independent of the degree of magnetic support. The relative power in the first 4 modes is similar. For the calculation with the highest magnetic field strength, however, (model 3a) the $`m=1`$ component dominates for a period of about 15 orbits around $`t/T40`$ (the numerical time unit $`T=1/\mathrm{\Omega }_{\mathrm{out}}`$, where $`\mathrm{\Omega }_{\mathrm{out}}`$ is the orbital frequency at the outer edge).

Fig. (6) shows the evolution of the mass in the inner disc ($`r_{\mathrm{in}}r0.3r_{\mathrm{out}}`$) in units of the total disc mass $`M_{\mathrm{disc}}`$. Mass piles up in the inner part of the disc, the faster the higher the magnetic support. This is accompanied by an outward transport of angular momentum by the magnetic instability. For model 3a we find a roughly linear increase of the inner disc mass with time corresponding to an accretion rate $`\dot{M}_{\mathrm{disc}}(r<0.3r_{\mathrm{out}})=\mathrm{2\hspace{0.33em}10}^4M_{\mathrm{disc}}/T`$, but superposed on this trend is a much more active ‘outburst’ around $`t40`$. During this active episode the accretion rate is about 10 times higher (i.e. for $`35<t/T<40`$). Fig. (5) shows that during the outburst $`P_1`$ is larger by at least a factor 3–5 compared to the other modes, and by a factor $`10`$ larger than in the preceding phase. Figs. (7) show that a strong m=1 spiral wave, traveling outward from a crescent-shaped disturbance in the inner disc, is present during the outburst.

The time-scale for mass accretion in model 3a is longer, by a factor $`10^3`$$`10^4`$, than the dynamical time-scale. For discs with less magnetic support the mass accretion time-scales are so long that we have been unable to follow their evolution beyond the initial development of the instabilities.

The evolution of the disc pattern as seen in $`B_\mathrm{z}`$ is shown in Fig. (7). The corresponding images of the surface density $`\mathrm{\Sigma }`$ are found in Stehle (1997). It is seen that the global instability starts with rather high mode numbers, $`m5`$…8 (Fig. 7a). Subsequently the waves are wound up (Fig. 7b) and it is only later that the modes $`m=1`$…4 become dominant. After about 30 orbits of the outer disc edge, a prominent $`m=1`$ spiral arm develops. It causes mass accretion rates 10 times higher than during the preceding phase. The relative strength in the $`m=`$1 component is largest near the inner disc edge. At $`r0.2r_{\mathrm{out}}`$ a prominent crescent-shaped field strength enhancement shows up (best seen in Fig. 7c). It is accompanied by a similar enhancement in the surface density (Stehle 1997). The crescent rotates approximately with the local orbital rate. An $`m=1`$ wave travels outward from this rotating crescent. The maximum of the crescent moves in to smaller radii (compare Figs. 7c and e), where the rotation rates are higher. The outward traveling wave becomes correspondingly more tightly wound.

This behavior is very reminiscent of a purely hydrodynamic form of global instability observed in hot, partially pressure-supported discs (Blaes and Hawley 1988; Różyczka & Spruit 1993). This instability takes place only when the degree of support against gravity by pressure becomes noticeable, and it also takes the form of a crescent rotating at the local orbital rate. It generates shock waves traveling outwards and inwards. These waves cause the mass in the disc to spread, while the angular momentum lost by the crescent causes it to spiral in towards the centre, behaving much like a solid object in doing so. Though the waves in the present calculation are rather different from hydrodynamic shock waves, we suspect that the same mechanism is at work. The peculiar behavior of the crescent mode and the fact that it appears only at certain phases suggests that it is a basically nonlinear phenomenon, not related directly to the linear global modes of the system. The nature of this phenomenon warrants further study.

With time, mass piles up near to the inner and outer edge of the accretion disc. The effect of the instability is thus much like that of viscous spreading, but it must be stressed that this is due to a global transport of angular momentum by the spiral wave, which can not be reduced to the action of a local viscosity. After the approximately 10 outer disc orbits during which the $`m=1`$ component dominated, its relative strength compared to the other modes decreases again and the phase of high mass accretion rates is finished (Fig. 7f).

At the end of the outburst the disc density distribution has completely changed. It now decreases approximately exponentially with radius rather than as a power law. The same is true for the $`B_\mathrm{z}`$ distribution (Fig. 8).

The dissipation taking place during the redistribution of mass and magnetic flux causes the disc to heat up, so that the gas pressure increases over its initial low value. At the end of the calculation the Mach number of the orbital motion in the inner disc has decreased to values of 5–10.

For the calculations done with magnetic support less than in model 3a the mass accretion rates are too low for significant redistribution of mass to occur over the 50 orbits we were able to follow. It is thus unclear if discs with less magnetic support also show outbursts, or if a threshold in the magnetic field strength exists below which the outburst mechanism cannot operate.

Since the outburst was a transient but obviously very effective transporter of mass and angular momentum, one wonders what caused its decline. After the outburst the degree of support of the disc against gravity by the magnetic forces has decreased substantially. Since the amplitude of all global waves seen in our results increases sharply with the degree of magnetic support, it is possible that the outburst declined because the crescent instability operates only at sufficiently high degrees of support. This view is consistent with the properties of the purely hydrodynamic crescent instability.

Another possibility is that an exponential dependence on radius is perhaps a more stable configuration, towards which the instability tends to develop. To test this possibility, we investigate in the next model a sequence of discs where $`\mathrm{\Sigma }`$ and $`B_\mathrm{z}`$ decrease exponentially with radius.

### 6.2 Model 4: $`B_\mathrm{z}\mathrm{exp}(r)`$ and $`\mathrm{\Sigma }\mathrm{exp}(r)`$

The final disc structure of model 3a motivates us to study the stability of accretion discs where the initial magnetic field and surface mass density decrease exponentially with radius. In this sequence we study two such cases, which differ only in the initial field strength. In model 4a we choose a field with magnetic support $`(g_\mathrm{m}/g)_{\mathrm{max}}=0.36`$ at maximum and in model 4b one with $`(g_\mathrm{m}/g)_{\mathrm{max}}=0.16`$. The initial structure of model 4a is shown in Fig. (9).

We take $`B_z/\mathrm{\Sigma }=`$ const. for the whole disc. $`\mathrm{\Sigma }(r)`$ and $`B_z(r)`$ decline linearly to zero at the inner edge. The initial axisymmetric distribution is again perturbed in $`v_r`$ by a low amplitude point–to–point noise.

Values for the ratio $``$ are listed in Tab. (2) for models 3 and in Tab. (3) for models 4. A comparison of these values shows that the total magnetic flux through the disc in models 4 is higher than in models 3, for the same degree of magnetic support. This difference is partly caused by the rapid decline of the magnetic field strength towards the disc edges in models 3 and partly from the fact that $`\mathrm{d}(B_\mathrm{z}/\mathrm{\Sigma })/\mathrm{d}r>0`$ for the initial distribution in model 3.

Fig. (10) shows the time evolution of the relative power $`P_\mathrm{m}`$ in the Fourier modes $`m=`$1–4, for models 4a and 4b. In both models the relative power $`P_\mathrm{m}`$ again increases initially on a short, dynamical, time-scale and the mode strength at which they saturate is similar to what we observe in model sequence 3. As expected from model 3a at times after the outburst, $`P_{\mathrm{m}=1}`$ is comparable to $`P_{\mathrm{m}=3}`$ and $`P_{\mathrm{m}=4}`$ but significantly weaker than $`P_{\mathrm{m}=2}`$. The $`m=2`$ mode appears to be the dominant mode during most time of our integration.

Even though $``$ is significantly higher than in model 3, the mass accretion rates at the inner part of the disc are small compared to model 3. The time-scale to clear the disc mass completely at these rates is now about $`10^5`$ outer edge orbits. This can be seen in Fig. (11), where we plot the evolution of the mass in the innermost part of the disc, $`M(r<0.3r_{\mathrm{out}})`$ in units of the total disc mass $`M_{\mathrm{disc}}`$ (see also Tab. 3). It is also interesting to note that the mass accretion is now much less dependent on the strength of the magnetic field. The low values for $`\dot{M}(r<0.3r_{\mathrm{out}})`$ are accompanied by a low level of power in the $`m=1`$ mode. This agrees with model 3, where the highest mass accretion rates are also found at a time where $`P_{\mathrm{m}=1}`$ is large.

In Fig. (12) we show the radial velocities of model 4a at $`t/T=51.1`$, the time where we stopped the integration, Fig. (13) shows the azimuthally averaged surface density $`<\mathrm{\Sigma }>`$ at that time.

The pattern shows tightly wound spiral arms which can be traced from one edge of the accretion disc to the other. The radial velocities are rather small, of the order $`10^2`$ of the orbital velocity at the outer edge of the disc. The mass redistribution has been strongest for $`r0.7r_{\mathrm{out}}`$ where it has taken only a few binary orbits. This is just the region where the magnetic support was strongest. In comparison, the mass redistribution in the inner regions is small.

In both sequences 3 and 4 the redistribution of mass and magnetic flux appears to be closely related to the degree of support against gravity by the magnetic field configuration. It appears in regions where $`g_\mathrm{m}/g`$ exceeds 5–10 per cent. The comparison also shows that, unlike interchange instability, the global instability that causes this redistribution is not directly related to the flux-to-mass distribution $`B/\mathrm{\Sigma }`$.

## 7 Conclusions and discussion

We have studied, by numerical simulation, the stability of accretion discs threaded by strong large scale magnetic fields, assuming that the discs are geometrically thin. The simulations solve the MHD equations for a vertically averaged accretion disc in the $`(r,\varphi )`$–plane. The disc magnetosphere is calculated in potential field approximation, i.e. we treat the disc as a current sheet and assume the magnetic field outside the disc to be current free. This approximation allows us to follow the evolution of a 3-dimensional field configuration with the computational effort of a 2-dimensional simulation.

Previous analytical studies by Spruit & Taam (1990), predicted that uniformly rotating discs would show an interchange-type instability, a local instability that appears at any location where the magnetic field $`B_z`$ decreases with radius more rapidly than the surface density $`\mathrm{\Sigma }`$. To test this prediction, as well as the stability of the numerical method, we first computed discs with initially uniform rotation. The results agree with the analytic stability condition and growth rates. No additional, non–local, forms of instability were found in the simulations of initially uniformly rotating discs. The nonlinear development of the instability agrees with that expected of an interchange instability like convection or Rayleigh-Taylor.

Linear analysis (Spruit et al. 1995; Lubow & Spruit 1995) predicts that differentially rotating discs with approximately Keplerian rotation are much more stable than uniformly rotating discs. It predicts that instability of the interchange type occurs only when the local shear rate is less than the growth rate of the instability in a uniformly rotating but otherwise identical disc. For smooth distributions of $`B_z`$ and $`\mathrm{\Sigma }`$ with $`r`$ this is equivalent to the condition that the magnetic forces contribute significantly to the support of the disc against gravity (Spruit et al. 1995). Such discs are strongly magnetized in the sense that $`v_\mathrm{A}\mathrm{\Omega }H`$, where $`v_\mathrm{A}`$ is the Alfén speed at the midplane of the disc and $`H`$ the disc thickness. The importance of this prediction is that it suggests that even quite strong poloidal magnetic fields might still be stable in an accretion disc.

To see if this prediction holds up in a full numerical simulation, we have done a sequence of calculations for discs with approximately Keplerian rotation, in which the magnetic contribution to support against gravity ranges from a few % to 36%. At these values, the linear results predict a weak form of interchange instability, with perturbations slowly growing as a power of $`t`$. In contrast with the uniformly rotating case, however, no evidence of this form of instability was found in the simulations. Instead, a new, global, exponentially growing form of instability appears. It takes the form of tightly wound spiral arms which can be traced from one edge of the accretion disc to the other. This instability was found both in cases where the magnetic field and the surface mass decrease with radius as a power law (model sequence 3) and where they decrease exponentially (model sequence 4). Its presence seems to depend primarily on the degree of magnetic support of the disc; its amplitude is a steep function of this quantity.

The instability acts similar to viscous spreading in the sense that an outward angular momentum transport takes place which allows mass accretion from the outer to the inner regions of the disc. Mass accretion is of the order of $`10^6`$$`10^4M_{\mathrm{disc}}/T`$ where $`M_{\mathrm{disc}}`$ is the total disc mass and T the Kepler time-scale $`1/\mathrm{\Omega }`$ at the outer edge of the accretion disc.

In the case with the highest degree of magnetic support an additional, more violent form of instability is observed. It takes place during a limited period of $`10`$ disc orbits. The mass accretion rate in this episode is 10 times higher than the long term average in the simulation. During this time the disc is perturbed by a strong one-armed spiral wave excited by a density enhancement rotating at the local orbital frequency at the location of highest magnetic support. At the end of the outburst the surface density and the magnetic field distribution in the disc decline exponentially with disc radius and the subsequent mass accretion is less efficient.

The time-scale for the redistribution can be rather large compared to the dynamical time-scale. Especially for discs with low magnetic support this time-scale is so large that, extrapolating the results from the simulation, the final configuration would be obtained only after some 1000–100000 disc orbits. Only when the magnetic support of the disc is of the order of the gravitational or centrifugal forces, does the time-scale become small enough that the redistribution of the disc mass proceeds on time-scales of several 10 disc orbits.

With the present results we can not establish whether the instabilities found can also lead to instability and redistribution of mass and angular momentum at lower degrees of magnetic support, or is limited to stronger fields. The difference between these possibilities is of obvious importance for accretion on longer time-scales. If the instabilities found here generally require magnetic support exceeding a few per cent, quite strong poloidal fields might exist in the inner regions of accretion discs.

In the protostellar context, such strong fields may also be relevant for the problem of binary formation (cf. Sigalotti and Klapp 2000)

ACKNOWLEDGMENTS

The work reported here was done as part of the research network ’accretion on to compact objects and protostars’ supported by EC grant ERB–CHRX–CT93–0329. R.S. has been supported by the ’Deutscher Akademischer Austauschdienst im Rahmen des 2.Hochschulsonderprogramms’. R.S. wants to thank J.Papaloizou for many stimulating discussions.


format: text

text:
# Interpolation in non-positively curved Kähler manifolds

## 1. Introduction

Let $`(X,\omega )`$ be a Kähler manifold, $`(L,h)X`$ be a holomorphic hermitian line bundle on $`X`$, and $`\mathrm{\Lambda }`$ a discrete set of points in $`X`$. Consider

$$_{\omega ,h}:=\{fH^0(X,L)/_X||f||_h^2dV_\omega <+\mathrm{}\}$$

and

$$l_{\mathrm{\Lambda },h}^2:=\{a\mathrm{\Gamma }(\mathrm{\Lambda },L)/\underset{p\mathrm{\Lambda }}{}||a(p)||_h^2<+\mathrm{}\}$$

We say that the discrete set $`\mathrm{\Lambda }`$ is interpolating for $`_{\omega ,h}`$, if for every $`a`$ in $`l_{\mathrm{\Lambda },h}^2`$ there is a section $`f`$ in $`_{\omega ,h}`$ such that $`f_{|\mathrm{\Lambda }}=a`$.

In the complex plane with its usual flat metric, there is a complete characterisation of interpolating sets for pluri-sub-harmonic weights comparable to $`|z|^2`$.

Theorem\[Berndtsson and Ortega Cerdà , Ortega-Cerdà and Seip \]

Consider the trivial line bundle on $``$ with metric $`h_\mathrm{\Phi }:=||^2e^\mathrm{\Phi }`$ where $`\mathrm{\Phi }`$ is a pluri-sub-harmonic function on $``$ whose Laplacian is uniformly bounded. A set $`\mathrm{\Lambda }`$ in $``$ is interpolating for $`_{h_\mathrm{\Phi }}`$ if and only if it is uniformly separated and

$$\rho >0,\epsilon >0z,\mathrm{\Delta }\mathrm{\Phi }(z)\frac{cardB(z,\rho )\mathrm{\Lambda }}{\rho ^2}+\epsilon .$$

Here $`B(z,\rho )`$ is the Euclidean ball of centre $`z`$ and radius $`\rho `$.

A simply-connected complete Riemannian manifold with non-positive sectional curvature is called a Cartan-Hadamard manifold. An hermitian complex manifold is called Cartan-Hadamard if its underlying Riemannian space is Cartan-Hadamard.

We prove

###### Theorem 1.

Let $`(X,\omega )`$ be a Kähler Cartan-Hadamard manifold with sectional curvature bounded from below by $`k^2`$. Let $`\mathrm{\Lambda }`$ be a discrete set of points in $`X`$ and $`(L,h)X`$ be an holomorphic hermitian line bundle on $`X`$ almost analytic around $`\mathrm{\Lambda }`$ (see definition 4).

If $`\mathrm{\Lambda }`$ is a uniformly $`\omega `$-separated set such that $`\rho >0,\epsilon >0,zX`$ ,

$$ic_h(L)+ricci(\omega )n\frac{cardB_\omega (z,\rho )\mathrm{\Lambda }}{\rho ^2}\left(1+k\rho \mathrm{coth}k\rho \right)\omega +\epsilon \omega $$

then it is interpolating for $`_{\omega ,h}`$.

On the other hand, Lindholm  obtained on $`^n`$ necessary conditions for interpolation in terms of the Monge-Ampère measure $`(i\overline{}\mathrm{\Phi })^n`$ of the weight. For $`n>1`$, there is therefore a gap between necessary and sufficient conditions expressed in terms of eigenvalues of the curvature.

Seip also obtained a complete description of interpolating sets in the unit disk with its hyperbolic measure for special weights $`A\mathrm{log}(1|z|^2)`$ in the trivial line bundle in terms of the hyperbolic lower density for $`\mathrm{\Lambda }`$. In the spirit of , we define for $`(X,\omega )`$ a Kähler Cartan-Hadamard manifold with sectional curvature bounded from above by $`1/\kappa ^2`$ , $`\mathrm{\Lambda }`$ a discrete set of points in $`X`$, and $`xX`$

$$D_{\mathrm{\Lambda },\kappa }(x):=\underset{\genfrac{}{}{0pt}{}{p\mathrm{\Lambda }}{d_\omega (x,p)1}}{}\mathrm{log}\mathrm{tanh}^2\left(\frac{d_\omega (x,p)}{2\kappa }\right).$$

We prove

###### Theorem 2.

Let $`(X,\omega )`$ be a Kähler Cartan-Hadamard manifold with sectional curvature bounded from above by $`1/\kappa ^2`$ and from below by $`k^2`$. Assume $`\mathrm{\Lambda }`$ is uniformly $`\omega `$-separated, and $`(L,h)`$ is almost analytic around $`\mathrm{\Lambda }`$. Then if

$$\underset{X}{sup}D_{\lambda ,\kappa }<+\mathrm{}\text{ and }ic_h(L)+ricci(\omega )\epsilon \omega $$

for some positive $`\epsilon `$, the set $`\mathrm{\Lambda }`$ is interpolating for $`_{\omega ,h}`$.

In the next section, we recall some comparison theorem in Riemannian geometry from  and . In section 2, we prove that under some assumption relating the complex structure and the metric structure uniformly on $`\mathrm{\Lambda }`$, the uniform separation is necessary. The section 4 is devoted to prove theorem 1 and the section 5 to prove theorem 2.

Note that using the extension theorem of Ohsawa-Takegoshi-Manivel , one can get interpolation results on Stein manifolds, expressed in terms of an equation for the set $`\mathrm{\Lambda }`$.

In the whole work, we will denote by $`C`$ any constant which only depends on quantities fixed before.

## 2. Comparison theorems in Riemannian geometry

In a Cartan-Hadamard manifold the exponential map is at each point a diffeomorphism. There is no conjugate points. The distance between two points is achieved by a unique geodesic. Hence, the distance function from any point is smooth. Recall that by a result of Wu (see  page 100) a Kähler Cartan-Hadamard manifold is Stein. As a consequence of these facts, any holomorphic line bundle on a Kähler Cartan-Hadamard manifold is trivialisable. Nevertheless, it is more convenient to state the results, especially the admissibility condition, in terms of hermitian line bundles instead of pluri-sub-harmonic weights.

By comparison theorem of Rauch (see for example  page 76), lower bounds on the sectional curvature yields upper bounds on the Jacobi vector fields and hence on the norm of the differential of the exponential map. By other forms of the comparison theorem (see  page 118-119), upper bounds on sectional curvature yields lower bounds on volume.

We will state the comparison theorem for the Hessian (see  Theorem A page 19) in a form suitable for our purpose. The proof of it relies on a precise form of the formula for second order variation of arc length, and on a relation between real and complex Hessian on Kähler manifolds (see  page 102).

Theorem Let $`M`$ be a Kähler Cartan-Hadamard manifold of dimension $`n`$ and constant sectional curvature $`k`$. Let $`f`$ be a non-decreasing function on $`[0,+\mathrm{}[`$. Let $`0`$ be a point in $`M`$. Assume $`i\overline{}f(d_M(0,))=g(d_M(0,))\omega _M`$ for a function $`g`$ on $`[0,+\mathrm{}[`$.

Let $`X`$ a Kähler Cartan-Hadamard manifold of dimension $`n`$ and sectional curvature bounded from below (resp. from above) by $`k`$. Let $`x_0X`$. Then

$$i\overline{}f(d_X(x_0,))g(d_X(x_0,))\omega _X.$$

resp.

$$i\overline{}f(d_X(x_0,))g(d_X(x_0,))\omega _X.$$

## 3. On uniform separation

We generalize the condition on the weight posed in earlier works on interpolation.

###### Definition 1.

Fix $`M_2`$ in $`[0,+\mathrm{}[`$. An hermitian line bundle $`(L,h)U`$ on an open set $`U`$ of $`^n`$ has an $`M_2`$-regular frame $`e`$ on $`U`$ if the weight function $`\mathrm{\Phi }:=\mathrm{log}h(e)`$ can be written as the sum

$$\underset{1\alpha N}{}|\sigma _\alpha |^2+\mathrm{\Phi }_{def}$$

where the $`\sigma _\alpha `$ are holomorphic functions on $`U`$ with oscillation $`sup_{x,yU}|\sigma _\alpha (x)\sigma _\alpha (y)|^2`$ bounded by $`M_2`$ and all second order real derivatives of the deformation weight $`\mathrm{\Phi }_{def}`$ are bounded by $`M_2`$ and $`M_2`$ on $`U`$.

###### Definition 2.

An hermitian line bundle $`(L,h)X`$ is said to be regular around $`\mathrm{\Lambda }`$ if there exists positive numbers $`r,\lambda ,\mu `$ and $`M_2`$ such that for all $`p\mathrm{\Lambda }`$, there exists a coordinates chart $`\mathrm{\Psi }_p:(U,0)^n(V,p)X`$ centered at $`p`$ such that

* $`\mathrm{\Psi }_p(U)B_\omega (p,r)`$
* $`(\mathrm{\Psi }_p^{}L,h)`$ has an $`M_2`$-regular frame $`e`$ on $`\mathrm{\Psi }_p^1(B_\omega (p,r))`$.
* The differential $`d\mathrm{\Psi }_p:(TU,\mathrm{\Psi }_p^{}\omega (p))(TV,\omega )`$ satisfies $`\mu d\mathrm{\Psi }_p\lambda `$ on $`\mathrm{\Psi }_p^1(B_\omega (p,r))`$.

Remark first that (iii) implies

$$\mu d_{eucl}d_\omega \lambda d_{eucl}\text{ on }B_\omega (p,r)$$

where $`d_{eucl}`$ is computed with the constant metric $`\mathrm{\Psi }_p^{}\omega (p)`$ on $`U`$.

In this section, we will prove the following

###### Theorem 3.

Let $`X,L`$ and $`\mathrm{\Lambda }`$ as in the introduction. Assume $`X`$ is a Kähler Cartan-Hadamard manifold and that $`(L,h)X`$ is regular around $`\mathrm{\Lambda }`$ .

Then, if $`\mathrm{\Lambda }`$ is interpolating for $`_{\omega ,h}`$, it is uniformly $`\omega `$-separated.

### 3.1. On regularity

It can be proved from the fact that $``$ is interpolating in $``$ for the weight $`|z|^2`$ that the non-uniformly separated set $`\{\sqrt{n},n\}`$ is interpolating in $``$ for the weight $`|z^2|^2`$. This explains our restriction on the metric $`h`$.

### 3.2. Mean value estimates

According to  (Theorem B page 43), for all $`p`$ in a Kähler Cartan-Hadamard manifold, for all holomorphic function $`f`$ on $`B_\omega (p,r)`$ and all $`xB_\omega (p,r/2)`$,

$$|f(x)|^2\frac{1}{vol_{eucl}(B_{eucl}(r/2))}_{B_\omega (x,r/2)}|f|^2𝑑V_\omega .$$

We now assume the hypothesis of theorem 3. The following estimates relies on the simple inequalities for $`x,z`$ in $`B_\omega (p,r)`$

$`|\sigma _\alpha (x)|^2+2\mathrm{}\left(\overline{\sigma _\alpha (x)}(\sigma _\alpha (z)\sigma _\alpha (x))\right)`$ $`=`$ $`|\sigma _\alpha (z)|^2|\sigma _\alpha (z)\sigma _\alpha (x)|^2`$
$``$ $`|\sigma _\alpha (z)|^2M_2`$
and with the coordinates $`(z_k)`$ provided by $`\mathrm{\Psi }_p`$,

$`\mathrm{\Phi }(x)+2\mathrm{}\left({\displaystyle \frac{\mathrm{\Phi }}{z_k}(x)(z_kz_k(x))}\right)`$ $``$ $`\mathrm{\Phi }(z){\displaystyle \frac{1}{2}}M_2(2n)^2|zx|_{eucl}^2`$
$``$ $`\mathrm{\Phi }(z){\displaystyle \frac{1}{2}}M_2(2n)^2{\displaystyle \frac{r^2}{\mu ^2}}.`$
Then, for all holomorphic section $`f`$ of $`L`$ on $`B_\omega (p,r)`$ and all $`x`$ in $`B_\omega (p,r/2)`$,

$`f(x)_h^2=|f(x)|^2e^{\mathrm{\Phi }(x)}`$
$``$ $`C{\displaystyle _{B_\omega (x,\frac{r}{2})}}\left|fe^{\underset{1\alpha N}{}\overline{\sigma _\alpha (x)}(\sigma _\alpha (z)\sigma _\alpha (x))\underset{1kn}{}\frac{\mathrm{\Phi }_{def}}{z_k}(x)(z_kz_k(x))}\right|^2e^{\mathrm{\Phi }(x)}𝑑V_\omega `$
$``$ $`C{\displaystyle _{B_\omega (x,r/2)}}f(x)_h^2𝑑V_\omega .`$
### 3.3. Gradient estimates

Let $`f`$ be a holomorphic section of $`L`$ on $`B_\omega (p,r)`$ and $`x`$ in $`B_\omega (p,r/4)`$. Denote by $`s`$ any coordinate among the $`z_kz_k(x)`$ centered at $`x`$. Then

$`{\displaystyle \frac{}{s}}(f(x)_h^2)_{|s=0}`$
$`=`$ $`\overline{f(x)}e^{\mathrm{\Phi }(x)}\left({\displaystyle \frac{f}{s}}_{|s=0}{\displaystyle \frac{\mathrm{\Phi }}{s}}_{|s=0}\right)`$
$`=`$ $`\overline{f(x)}e^{\mathrm{\Phi }(x)}{\displaystyle \frac{}{s}}\left(f(s)e^{\underset{1\alpha N}{}\overline{\sigma _\alpha (x)}(\sigma _\alpha (s)\sigma _\alpha (x))\frac{\mathrm{\Phi }_{def}}{s}(0)s}\right)_{|s=0}`$
$`=`$ $`\overline{f(x)}e^{\frac{\mathrm{\Phi }(x)}{2}}{\displaystyle \underset{|s|=\frac{r}{4\lambda }}{}}f(s)e^{\underset{1\alpha N}{}\overline{\sigma _\alpha (x)}(\sigma _\alpha (s)\sigma _\alpha (x))\frac{\mathrm{\Phi }_{def}}{s}(0)s}e^{\frac{\mathrm{\Phi }(x)}{2}}{\displaystyle \frac{ds}{s}}.`$
Now, note that thanks to the regularity assumption

$$\left|f(s)e^{\underset{1\alpha N}{}\overline{\sigma _\alpha (x)}(\sigma _\alpha (z)\sigma _\alpha (x))\frac{\mathrm{\Phi }_{def}}{s}(0)s}\right|^2e^{\mathrm{\Phi }(x)}Cf(s)_h^2.$$

By the mean value estimate, we get

$$\left|\frac{}{s}(f(x)_h^2)_{|s=0}\right|C_{B_\omega (p,r)}f_h^2𝑑V_\omega .$$

### 3.4. Uniform separation

Choose an interpolating set $`\mathrm{\Lambda }`$ for $`_{\omega ,h}`$. Assume it is not uniformly $`\omega `$-separated. Then, there exists a subset $`\mathrm{\Lambda }^{}`$ (still interpolating) non-uniformly separated but such that any intersection of three balls among the $`B_\omega (p_i,r)`$ is empty. Hence, by the mean value estimate, the evaluation map

$$\begin{array}{cccc}ev_\mathrm{\Lambda }^{}:& _{\omega ,h}& & l_{\mathrm{\Lambda }^{},h}^2\\ & f& & f_{|\mathrm{\Lambda }^{}}\end{array}$$

is well defined and continuous. By the open mapping theorem, $`\mathrm{\Lambda }^{}`$ is stably interpolating, i.e.

$$K,al_{\mathrm{\Lambda }^{},h}^2,f_{\omega ,h}/f_{|\mathrm{\Lambda }^{}}=a\text{ and }_Xf_h^2𝑑V_\omega K\underset{p\mathrm{\Lambda }^{}}{}|a(p)|_h^2<+\mathrm{}.$$

For all $`p_0`$ in $`\mathrm{\Lambda }^{}`$ denote by $`f^{p_0}`$ any stable solution for $`a0`$ at $`p_0`$ and $`a=0`$ on $`\mathrm{\Lambda }^{}\{p_0\}`$. By the gradient estimate, we infer if $`p\mathrm{\Lambda }^{}B_\omega (p_0,r/4)`$,

$$f^{p_0}(p_0)_h^2=\left|f^{p_0}(p_0)_h^2f^{p_0}(p)_h^2\right||pp_0|_X||f||_h^2dV_\omega \frac{K}{\mu }d_\omega (p,p_0)f^{p_0}(p_0)_h^2$$

which contradicts the assumption that $`\mathrm{\Lambda }^{}`$ is not uniformly separated.

## 4. Interpolation on non-positively curved manifolds

We will prove theorem 1 in this section. We first have to describe the admissible metrics.

###### Definition 3.

An hermitian line bundle $`(L,h)U`$ on an open set $`U`$ of $`^n`$ has an $`M_2`$-almost analytic frame $`e`$ on $`U`$ if the weight function $`\mathrm{\Phi }:=\mathrm{log}h(e)`$ can be written as the sum

$$\underset{1\alpha N}{}|\sigma _\alpha |^2+\mathrm{\Phi }_{def}$$

where the $`\sigma _\alpha `$ are holomorphic functions on $`U`$ and all second order real derivatives of the deformation weight $`\mathrm{\Phi }_{def}`$ are bounded by $`M_2`$ and $`M_2`$ on $`U`$.

Note that nothing is now assumed on the oscillation of the holomorphic functions $`\sigma _\alpha `$.

###### Definition 4.

An hermitian line bundle $`(L,h)X`$ is said to be almost analytic around $`\mathrm{\Lambda }`$ if there exists positive numbers $`r_0`$ and $`M_2`$ such that for all $`p\mathrm{\Lambda }`$, there exists a coordinates chart $`\mathrm{\Psi }_p:(U,0)^n(V,p)X`$ centered at $`p`$ such that

* $`\mathrm{\Psi }_p(U)B_\omega (p,r_0)`$
* $`(L,h)`$ has a $`M_2`$-almost analytic frame $`e`$ on $`\mathrm{\Psi }_p^1(B_\omega (p,r_0))`$.
* The differential $`d\mathrm{\Psi }_p:(TU,\mathrm{\Psi }_p^{}\omega (p))(TV,\omega )`$ satisfies $`\mu d\mathrm{\Psi }_p`$
on $`\mathrm{\Psi }_p^1(B_\omega (p,r_0))`$.

### 4.1. Local choice

We now make the assumptions of theorem 1. Denote by $`2\delta _0:=inf_{\genfrac{}{}{0pt}{}{p,q\mathrm{\Lambda }}{pq}}(d_\omega (p,q),r_0)`$. Let $`al_{\mathrm{\Lambda },h}^2`$ and $`p\mathrm{\Lambda }`$. On the ball $`B_\omega (p,\delta _0)`$, consider the holomorphic function

$$f_p:=a(p)\mathrm{exp}\left(\underset{1kn}{}\frac{\mathrm{\Phi }_{def}}{z_k}(p)(z_kz_k(p))+\underset{1\alpha N}{}\overline{\sigma _\alpha (p)}(\sigma _\alpha \sigma _\alpha (p))\right).$$

It is in fact a normal frame at $`p`$ for $`(L,h)`$. The following estimates relies on the simple inequalities for $`z`$ in $`B_\omega (p,\delta _0)`$ which may be regarded as higher dimensional analogues of Riesz formula.

$`|\sigma _\alpha (z)|^2+2\mathrm{}\left(\overline{\sigma _\alpha (p)}(\sigma _\alpha (z)\sigma _\alpha (p))\right)`$ $`=`$ $`|\sigma _\alpha (p)|^2|\sigma _\alpha (z)\sigma _\alpha (p)|^2`$
$``$ $`|\sigma _\alpha (p)|^2.`$
and

$`\mathrm{\Phi }_{def}(z)+2\mathrm{}\left({\displaystyle \underset{1kn}{}}{\displaystyle \frac{\mathrm{\Phi }_{def}}{z_k}}(p)(z_kz_k(p))\right)`$ $``$ $`\mathrm{\Phi }_{def}(p)+{\displaystyle \frac{1}{2}}M_2(2n)^2|zp|_{eucl}^2`$
$``$ $`\mathrm{\Phi }_{def}(p)+{\displaystyle \frac{1}{2}}M_2(2n)^2{\displaystyle \frac{\delta _0^2}{\mu ^2}}`$
Then, for all $`zB_\omega (p,\delta _0)`$,

$`f_p(z)_h^2=|f_p(z)|^2e^{\mathrm{\Phi }(z)}`$ $``$ $`Ca(p)_h^2`$
where $`C`$ is a constant independent of $`p\mathrm{\Lambda }`$.

### 4.2. Gluing

Let $`\chi `$ be a cut-off function defined on $`^+`$ decreasing, equal to 1 on $`[0,1/4]`$, and to 0 on $`[1,+\mathrm{}[`$. Then,

$$F:=\underset{p\mathrm{\Lambda }}{}f_p\chi \left(\frac{d_\omega (p,)^2}{\delta _0^2}\right)$$

is a smooth solution of the interpolation problem, holomorphic around $`\mathrm{\Lambda }`$, with the growth condition

$`{\displaystyle _X}F_h^2𝑑V_\omega `$ $``$ $`{\displaystyle \underset{p\mathrm{\Lambda }}{}}{\displaystyle _{B(p,\delta _0)}}Ca(p)_h^2𝑑V_\omega `$
$``$ $`C{\displaystyle \underset{p\mathrm{\Lambda }}{}}a(p)_h^2<+\mathrm{}.`$
Here, we used Rauch theorem to infer that the balls $`B(p,\delta _0)`$ have volume bounded by a constant independent of $`p`$.

### 4.3. Finding holomorphic sections

We will apply the $`L^2`$-estimates existence theorem with a singular weight (see for example  theorem 5.1) in order to keep the right values on $`\mathrm{\Lambda }`$. Consider the auxiliary weight

$$v(z):=n\underset{q\mathrm{\Lambda }}{}\left(1\frac{d_\omega (q,z)^2}{\rho ^2}+\mathrm{log}\frac{d_\omega (q,z)^2}{\rho ^2}\right)1_{B_\omega (q,\rho )}(z)$$

Remark first that,

$`{\displaystyle _X}\overline{}F_h^2e^v𝑑V_\omega `$ $`=`$ $`{\displaystyle _X}{\displaystyle \underset{p\mathrm{\Lambda }}{}}\left|\overline{}\chi \left({\displaystyle \frac{d_\omega (p,z)^2}{\delta _0^2}}\right)\right|^2f_p_h^2e^vdV_\omega `$
The term $`|\overline{}\chi \left(\frac{d_\omega (p,z)^2}{\delta _0^2}\right)|^2`$ depends on the differential of the exponential map. Hence it is bounded by a constant independent of the point $`p`$ by Rauch theorem. Now,

$`{\displaystyle _X}\overline{}F_h^2e^v𝑑V_\omega `$ $``$ $`C{\displaystyle \underset{p\mathrm{\Lambda }}{}}a(p)_h^2{\displaystyle _{B_\omega (p,\delta _0)B_\omega (p,\delta _0/2)}}e^v𝑑V_\omega .`$
We have to estimate $`v`$ around $`p\mathrm{\Lambda }`$. Choose $`q`$ in $`\mathrm{\Lambda }`$. For $`z(B_\omega (p,\delta _0)B_\omega (p,\delta _0/2))B_\omega (q,\rho )`$, $`d_\omega (q,z)d_\omega (p,q)d_\omega (z,p)2\delta _0\delta _0/2`$ if $`qp`$, and $`d_\omega (q,z)\delta _0/2`$ if $`q=p`$. Hence,

$$1\frac{d_\omega (q,z)^2}{\rho ^2}+\mathrm{log}\frac{d_\omega (q,z)^2}{\rho ^2}\mathrm{log}\frac{d_\omega (q,z)^2}{\rho ^2}\mathrm{log}\frac{\delta _0^2}{4\rho ^2}$$

Recalling that the number of $`q`$ in $`\mathrm{\Lambda }`$ such that $`zB_\omega (q,\rho )`$ is uniformly bounded because $`\mathrm{\Lambda }`$ is assumed to be uniformly separated and the volume of $`B_\omega (z,\rho )`$ is bounded from above by the volume of a ball radius $`\rho `$ in the hyperbolic space of dimension $`n`$ and of sectional curvature $`k^2`$, one infers that

$$_X\overline{}F_h^2e^v𝑑x<+\mathrm{}.$$

Now, the computation of the curvature is simple because the derivatives of

$$\left(1\frac{d_\omega (p,z)^2}{\rho ^2}+\mathrm{log}\frac{d_\omega (p,z)^2}{\rho ^2}\right)$$

vanishes along $`B_\omega (p,a)`$.

$`i\overline{}v`$ $`=`$ $`n{\displaystyle \underset{p\mathrm{\Lambda }}{}}\left({\displaystyle \frac{i\overline{}d_\omega (p,z)^2}{\rho ^2}}+i\overline{}\mathrm{log}d_\omega (p,z)^2\right)1_{B(p,a)}(z).`$
By the comparison theorem with the model $`M=(,||^2)`$ and $`f=\mathrm{log}`$ we get that $`\mathrm{log}d_\omega (z,p)^2`$ is a pluri-sub-harmonic function. On the other hand, using the hyperbolic space of curvature $`k^2`$ as model, with $`f=x^2`$, we get (see  page 35)

$$i\overline{}d_\omega (p,z)^2(1+kd_\omega (p,z)\mathrm{coth}kd_\omega (p,z))\omega (1+k\rho \mathrm{coth}k\rho )\omega $$

on $`B_\omega (p,\rho )`$. Hence,

$`ic_h(L)+i\overline{}v+ricci(\omega )`$
$``$ $`ic_h(L)+ricci(\omega )n{\displaystyle \frac{cardB_\omega (z,\rho )\mathrm{\Lambda }}{\rho ^2}}(1+k\rho \mathrm{coth}k\rho )\omega `$
$``$ $`\epsilon \omega .`$
Solving the equation $`\overline{}G=\overline{}F`$ on $`X`$ endowed with the complete metric $`\omega `$ for the line bundle $`L`$ endowed with the metric $`he^v`$, we get a smooth section $`G`$ of $`L`$ on $`X`$ such that $`FG`$ is holomorphic and $`_XG_h^2e^v𝑑V_\omega `$ is finite. Now, because $`v`$ is non-positive, $`_XG_h^2𝑑V_\omega `$ is also finite. Furthermore, the estimate of $`v`$ around $`p\mathrm{\Lambda }`$ leads to

$$_{B_\omega (p,\delta _0)}\frac{G_h^2}{zp^{2n}}C_XG_h^2e^v𝑑x<+\mathrm{},$$

so that $`G`$ vanishes on $`\mathrm{\Lambda }`$.

### 4.4. Comments on the factor $`\rho ^2`$

The first remark is an homogeneity property. The quantity

$$\frac{cardB_\omega (z,\rho )\mathrm{\Lambda }}{\rho ^2}\omega $$

which appears in theorem 1 applied for flat manifolds is equal to the quantity $`{\displaystyle \frac{cardB_{t\omega }(z,\sqrt{t}\rho )\mathrm{\Lambda }}{(\sqrt{t}\rho )^2}}t\omega `$ for all positive $`t`$.

The second remark is a restriction property, which, together with the criterion for interpolation in $``$ yields a reason for the factor $`\rho ^2`$.

###### Proposition 1.

Consider $`^n`$ with its usual flat metric $`\omega `$ and $`(L,h)^n`$ an hermitian line bundle with $`M_2`$-regular frame on a neighbourhood of $`^{n1}`$ in $`^n`$.

If $`\mathrm{\Lambda }`$ is interpolating for $`_{\omega ,h}`$ then $`\mathrm{\Lambda }^{n1}`$ is interpolating for $`_{\omega _{|^{n1}},h_{|^{n1}}}`$.

Proof Let $`al_{\mathrm{\Lambda }^{n1},h}^2`$. There exists a function $`f`$ in $`_{\omega ,h}`$ achieving $`f(p)=a(p)`$ for $`p\mathrm{\Lambda }^{n1}`$ and $`f(p)=0`$ for $`p\mathrm{\Lambda }^{n1}`$. To check the growth condition,

$`|f(z^{},0)|^2e^{\mathrm{\Phi }(z^{},0)}`$
$``$ $`{\displaystyle \frac{1}{\pi \delta _0^2}}{\displaystyle \underset{z_nD(0,\delta _0)}{}}\left|f(z^{},z_n)e^{\underset{1\alpha N}{}\overline{\sigma _\alpha (z^{},0)}(\sigma _\alpha (z)\sigma _\alpha (z^{},0))\frac{\mathrm{\Phi }_{def}}{z_n}(z^{},0)z_n}\right|^2e^{\mathrm{\Phi }(z^{},0)}`$
$``$ $`Ce^{1/2M_22^2\delta _0^2+NM_2}{\displaystyle \underset{z_nD(0,\delta _0)}{}}|f(z^{},z_n)|^2e^{\mathrm{\Phi }(z^{},z_n)}𝑑x_n𝑑y_n.`$
We can now apply Fubini theorem because the measure on $`^n`$ is a product measure. We get

$`_{^{n1}}|f(z^{},0)|^2e^{\mathrm{\Phi }(z^{},0)}C_^n|f(z)|^2e^{\mathrm{\Phi }(z)}<+\mathrm{}.`$

## 5. Interpolation on negatively curved manifolds

We will now prove theorem 2.

The ball $`𝔹(\kappa )`$ of radius $`\kappa `$ in $`^{2n}`$ endowed with the metric $`{\displaystyle \frac{4|dz|^2}{(1\frac{|z|^2}{\kappa ^2})^2}}`$ has constant sectional curvature $`1/\kappa ^2`$. The distance between $`0`$ and $`z`$ is $`2\kappa \text{arctanh}\frac{|z|}{\kappa }`$. Hence,

$$i\overline{}\mathrm{log}\mathrm{tanh}^2\frac{d_{𝔹(\kappa )}(0,)}{2\kappa }=i\overline{}\mathrm{log}|z|^20.$$

Hence, by comparison theorem, for every $`pX`$,

$$i\overline{}\mathrm{log}\mathrm{tanh}^2\frac{d_\omega (p,)}{2\kappa }0.$$

Moreover, $`v_p:=n\mathrm{log}\mathrm{tanh}^2\frac{d_X(p,)}{2\kappa }`$ is negative and has a pole of order $`2n`$ at $`p`$. Consider $`v:=_{p\mathrm{\Lambda }}v_p`$. By assumption on the density,

$$v(z)\underset{X}{sup}D_{\mathrm{\Lambda },\kappa }+\underset{\genfrac{}{}{0pt}{}{p\mathrm{\Lambda }}{d_\omega (z,p)1}}{}v_p(z)$$

Just as we estimated the auxiliary weight in the previous proof, we can estimate $`v`$ on $`B_\omega (p,\delta _0)B_\omega (p,\delta _0/2)`$. The rest of the proof follows the same lines as the proof of theorem 1.

The function $`v`$ may be compared to a pluri-complex Green function having poles on $`\mathrm{\Lambda }`$.