text:
# BLACK HOLE MAGNETOSPHERES AROUND THIN DISKS DRIVING INWARD AND OUTWARD WINDS

## 1 INTRODUCTION

It is widely believed that a supermassive black hole surrounded by an accretion disk can work as a central engine of active galactic nuclei. If magnetic field exists around the black hole and the disk, it is important to understand the magnetospheric structure for explaining the highly energetic phenomena. The magnetohydrodynamical (MHD) approach in general relativity may be necessary to study various features of the black hole magnetospheres. Unfortunately, even if the magnetosphere is assumed to be stationary and axisymmetric, any global models describing both inward winds in the strong-gravity region near black hole and outward winds in the far distant region have not been constructed, in spite of a large number of works devoted to the MHD problems. This is mainly because it is very difficult to analyze the highly nonlinear Grad-Shafranov equation for the magnetic stream function which defines the poloidal field lines (see e.g., Beskin (1997) for review).

To avoid the mathematical difficulty in the full MHD system, one may consider vacuum solutions of the magnetic stream function as a practical approximation applicable to the magnetically dominated models. In particular, the well-known examples such as split-monopole, uniform and paraboloidal fields in Schwarzschild geometry (Michel, 1973; Wald, 1974; Blandford & Znajek, 1977) have been used as the zeroth-order approximation for the perturbation method in slow-rotation and force-free fields, from which the first-order perturbation has been derived for giving a weak toroidal magnetic field and discussing the Blandford-Znajek process (Blandford & Znajek, 1977). Another important MHD problem to be studied is fluid motion along poloidal field lines. For examples, the existence of accretion flows onto the black hole passing through some MHD critical points has been shown under a fixed shape of field lines (Takahashi et al., 1990). The vacuum solutions allowing a field configuration of astrophysical interest will become a useful tool in an attempt to discuss the trans-critical motion of outgoing and ingoing magnetized winds, except in the region where the effect of fluid inertia should dominate to change crucially the field structure.

Of course, even for the vacuum fields, some current distribution should be assumed to exist in a restricted region of the magnetosphere. In this paper we focus on the magnetic field generated by a toroidal current distribution in a thin disk around the Schwarzschild black hole. The source of the above-mentioned fields with split-monopole and paraboloidal structures can be also disk currents. However, such models contain only open field lines and fail to describe a magnetic connection of the disk with the black hole. This will be partially due to the absence of the disk inner edge separated from the event horizon. The magnetic connection by closed field lines is expected to exist and have astrophysically important effects (Punsly, 1991a; Gruzinov, 1999; Blandford, 2000; Li, 2000c; Krolik, 2000). Hence, our purpose here is to present explicitly a global magnetospheric model with both closed and open field lines threading a thin disk and connecting magnetically the inner and outer parts, respectively, with the black hole and the far distant region, which is relevant to the problems of inward and outward disk-driven winds. Though the vacuum model would miss some important features of more realistic magnetospheres, our modeling could be a preliminary step for understanding the various MHD processes in further detail.

To obtain general models based on a disk current distribution, it may be convenient to consider linear superposition of the magnetic field generated by a single ring current located at a fixed radius on the equatorial plane which has been found in Kerr geometry (Petterson, 1975). Assuming the current distribution from an inner edge to an outer one, Li (2000a) has derived a formal expression of the magnetic stream function in terms of the infinite sum of the multipole fields, without revealing the global shape of the field lines. Because the each multipole field of which the structure has been discussed by Ghosh (2000) still fails to give closed field lines regular at the event horizon, the infinite sum is essential to understanding the global shape.

In §2 we would like to develop a different approach to arrive at the infinite sum formula for the magnetic stream function. In §3 we succeed in reducing it to an analytical simple form by locating the disk inner edge at a radius much larger than the horizon radius and the outer edge at infinite radius. The global structure of closed and open field lines is shown, and it is found that only the very narrow region of the disk separated from the inner edge can be magnetically connected with the black hole, together with the result that the field lines become nearly cylindrical at the event horizon and conical at infinity. In §4 we also consider a modification of the global shape due to an addition of a uniform external field. Though collimated field lines appear at the far distant region, it is found that for a strong uniform field the magnetic connection between the black hole and the disk breaks down. Finally we discuss the MHD effect due to slow rotation of the magnetosphere (namely, angular momentum transfer by magnetically dominated inward winds from the disk to the black hole), which could lead to radial accretion of disk plasma and gradual disruption of the closed field lines (Punsly & Coroniti, 1989). In the following we use the geometrical units with $`G=c=1`$.

## 2 THE MAGNETIC STREAM FUNCTION

Let us consider a stationary, axisymmetric magnetic field in Schwarzschild background of the metric

$$ds^2=\left(1\frac{2M}{r}\right)dt^2+\left(1\frac{2M}{r}\right)^1dr^2+r^2\left(d\theta ^2+\mathrm{sin}^2\theta d\varphi ^2\right),$$
(1)
where $`M`$ is the mass of the black hole. The poloidal components of the magnetic field may be given by

$$\stackrel{}{B}_p=\frac{\stackrel{}{}\mathrm{\Psi }\times \stackrel{}{\varphi }}{2\pi r\mathrm{sin}\theta },$$
(2)
where $`\stackrel{}{\varphi }`$ is a unit vector in the azimuthal direction, and $`\mathrm{\Psi }(r,\theta )`$ is the so-called magnetic stream function defining the poloidal field lines by $`\mathrm{\Psi }=`$ constant. If the vacuum field is assumed, the Maxwell equations reduce to the simple form

$$r^2\frac{}{r}\left[\left(1\frac{2M}{r}\right)\frac{\mathrm{\Psi }}{r}\right]+\mathrm{sin}\theta \frac{}{\theta }\left(\frac{1}{\mathrm{sin}\theta }\frac{\mathrm{\Psi }}{\theta }\right)=0,$$
(3)
which allows the separable solutions

$$\mathrm{\Psi }=R_\nu (r)\mathrm{\Theta }_\nu (\theta ).$$
(4)
with the separation constant $`\nu `$

For the analysis of these functions it is convenient to introduce the variables

$$w=\frac{r}{M}1,x=\mathrm{cos}\theta .$$
(5)
Then, from the regularity at the polar axis $`\theta =0,\pi `$ we obtain the angular part written by

$$\mathrm{\Theta }_\nu =_{|x|}^1P_\nu (x^{})𝑑x^{},$$
(6)
where $`P_\nu `$ is the Legendre function of the first kind. Though in this paper we focus on the field with equatorial symmetry, the parameter $`\nu `$ is not limited to a positive odd integer for allowing a discontinuous change of the radial component of $`\stackrel{}{B}_p`$ on the upper and lower surfaces of the equatorial disk.

To determine the radial part $`R_\nu `$, the condition for regularity at the event horizon $`r=2M`$ and at infinity $`r\mathrm{}`$ becomes important. We use the Legendre function $`P_\nu (w)`$ to give the solution regular at $`w=1`$ of the form

$$R_\nu =(w+1)P_\nu (w)_1^wP_\nu (w^{})𝑑w^{},$$
(7)
where we assume to be $`\nu =(1/2)+ik`$ with arbitrary positive $`k`$, because the oscillatory real function $`P_\nu (w)`$ has the amplitude decreasing in proportion to $`1/\sqrt{w}`$ at large $`w`$, and the linear superposition

$$\mathrm{\Psi }(r,\theta )=_0^{\mathrm{}}f(k)R_\nu (r)\mathrm{\Theta }_\nu (\theta )𝑑k$$
(8)
will be able to describe a global magnetic field satisfying the boundary conditions.

Now we calculate the coefficient $`f(k)`$ in the linear superposition corresponding to the case of a disk current with the inner edge at $`w=w_0>1`$. Recall that in previous works the solution was written by an infinite sum of the multipole fields with $`\nu `$ of natural numbers. Then, in the distant region $`w>w_0`$ the Legendre function $`P_\nu (w)`$ should be replaced by $`Q_\nu (w)`$ with a change of the coefficient due to the condition for continuity of $`\mathrm{\Psi }(r,\theta )`$ at $`w=w_0`$. We can avoid these slightly complicated steps by using the integral form (8). Though the toroidal surface current in the thin disk is proportional to $`r^2(\mathrm{\Psi }/\theta )`$ at $`\theta =\pi /2`$, here we rather focus on the equality

$$\frac{}{r}\left(\frac{\mathrm{\Psi }}{\theta }\right)=\frac{w+1}{M}\frac{d}{dw}h(w),$$
(9)
which is evaluated at the upper surface of the equatorial plane. The remarkable point is that the function $`h(w)`$ is given by

$$h(w)=_0^{\mathrm{}}g(k)P_\nu (w)𝑑k,$$
(10)
where we have

$$g(k)=P_\nu (0)f(k),$$
(11)
and the inversion formula (Erdélyi et al., 1953) holds as follows,

$$g(k)=k\mathrm{tanh}(\pi k)_1^{\mathrm{}}P_\nu (w)h(w)𝑑w.$$
(12)
Then, if the function $`h(w)`$ is specified according to the current distribution, we can derive the coefficient $`f(k)`$ giving the magnetic stream function $`\mathrm{\Psi }(r,\theta )`$ valid in the whole region.

Further, to obtain the explicit form of $`f(k)`$ through the well-known integral formula involving two Legendre functions, let us assume in the disk region $`w>w_0`$ to be

$$h(w)=A_0Q_s(w).$$
(13)
Though $`A_0`$ may be an arbitrary real parameter, we limit $`s`$ to be positive for assuring that $`\mathrm{\Psi }/\theta `$ does not diverge at large $`w`$, where the Legendre function $`Q_s(w)`$ of the second kind falls off as $`w^{s1}`$. (If any different model of the current distribution becomes necessary, the linear superposition of $`Q_s`$ with various $`s`$ in the range $`0<s<\mathrm{}`$ is a possible procedure to construct the solution $`\mathrm{\Psi }`$. However, we would like to emphasize that this simple example for the current distribution can be useful for revealing interesting features of the magnetospheric structure, as will be shown in the following.) Because in the inner region $`1w<w_0`$, where the field lines should vertically thread the equatorial plane, we have

$$h(w)=0,$$
(14)
it is easy to arrive at the result

$$g(k)=\frac{A_0k\mathrm{tanh}(\pi k)}{k^2+[s+(1/2)]^2}F_\nu ,$$
(15)
where

$$F_\nu =(w_0^21)\left[P_\nu (w_0)\frac{dQ_s(w_0)}{dw_0}Q_s(w_0)\frac{dP_\nu (w_0)}{dw_0}\right].$$
(16)
Strictly speaking, even if the condition (14) is assumed, equation (9) only implies that $`\mathrm{\Psi }/\theta `$ is constant at the upper or lower surface of the equatorial plane in the inner range $`1w<w_0`$. Hence, we must confirm the validity of equation (15), by calculating the value on the event horizon $`w=1`$, which is given by

$$\frac{\mathrm{\Psi }}{\theta }=_0^{\mathrm{}}2g(k)𝑑k,$$
(17)
at the upper surface of the equatorial plane. In fact, the right-hand side of equation (17) is shown to vanish by virtue of the equality

$$Q_s(w_0)=_0^{\mathrm{}}\frac{k\mathrm{tanh}(\pi k)}{k^2+[s+(1/2)]^2}P_\nu (w_0)𝑑k$$
(18)
which is a result of the inversion formula (12).

We have obtained the magnetic stream function of the integral form (8) in which the coefficient $`f(k)=g(k)/P_\nu (0)`$ is given by equation (15). To calculate the integral with respect to $`k`$ with the application of the residue theorem, it will be useful to consider the integration contour along the real axis ($`\mathrm{}<k<\mathrm{}`$) and the large semi-circle in the complex $`k`$-plane, based on the relation between two Legendre functions written by

$$P_\nu (z)=\frac{i}{\pi }\mathrm{coth}(\pi k)\left[Q_\nu (z)Q_{\nu 1}(z)\right].$$
(19)
If $`w<w_0`$ we replace $`P_\nu (w_0)`$ in the integrand according to equation (19), and the large semi-circle is chosen to be in the lower (or upper) half of the complex k-plane for the integral of the term involving $`Q_\nu `$ (or $`Q_{\nu 1}`$). Then, by taking account of the contribution only from the poles present at $`k=\pm i[s+(1/2)]`$ and at $`k=\pm i[2n+(3/2)]`$ where $`P_\nu (0)`$ vanishes for $`n=0,1,2,\mathrm{}`$, we obtain the result valid in the inner region $`1w<w_0`$ as follows,

$$\mathrm{\Psi }=\underset{n=0}{\overset{\mathrm{}}{}}\alpha _nG_{2n+1}R_{2n+1}(r)\mathrm{\Theta }_{2n+1}(\theta ),$$
(20)
where

$$\alpha _n=\frac{A_0(1)^n(4n+3)}{(2n+s+2)(2n+1s)}\frac{(2n+1)!!}{(2n)!!}$$
(21)
and

$$G_{2n+1}=(w_0^21)\left[Q_{2n+1}(w_0)\frac{dQ_s(w_0)}{dw_0}Q_s(w_0)\frac{dQ_{2n+1}(w_0)}{dw_0}\right].$$
(22)
In the outer region $`w>w_0`$, however, we must use equation (19) for $`P_\nu (w)`$ \[instead of $`P_\nu (w_0)`$\] with the same choice of the large semi-circle, and the contribution from the poles is divided into the two parts

$$\mathrm{\Psi }=\mathrm{\Psi }_1+\mathrm{\Psi }_2$$
(23)
with the forms written by

$$\mathrm{\Psi }_1=\underset{n=0}{\overset{\mathrm{}}{}}\alpha _n\mathrm{\Theta }_{2n+1}(\theta )\left[F_{2n+1}\stackrel{~}{R}_{2n+1}(r)G_{2n+1}_1^{w_0}P_{2n+1}(w^{})𝑑w^{}\right],$$
(24)
and

$$\mathrm{\Psi }_2=\frac{A_0}{\sqrt{\pi }}\mathrm{\Gamma }\left(\frac{1s}{2}\right)\mathrm{\Gamma }\left(1+\frac{s}{2}\right)\mathrm{\Theta }_s(\theta )\stackrel{~}{R}_s(r)$$
(25)
where the function $`\stackrel{~}{R}_\nu (r)`$ for any $`\nu `$ is defined by

$$\stackrel{~}{R}_\nu (r)=(w+1)Q_\nu (w)_{w_0}^wQ_\nu (w^{})𝑑w^{},$$
(26)
and $`\mathrm{\Gamma }(z)`$ is the gamma function.

It is easy to check that the continuity of these expressions (20) and (23) for $`\mathrm{\Psi }`$ at the boundary $`w=w_0`$ holds, because through the application of the expansion in terms of the Legendre polynomials to the function $`u(x)`$ defined by $`u=P_s(x)`$ for $`x>0`$ and $`u=P_s(|x|)`$ for $`x<0`$ we obtain

$$\underset{n=0}{\overset{\mathrm{}}{}}\frac{(1)^n(4n+3)}{(2n+s+2)(2n+1s)}\frac{(2n+1)!!}{(2n)!!}P_{2n+1}(x)=\frac{1}{\sqrt{\pi }}\mathrm{\Gamma }\left(\frac{1s}{2}\right)\mathrm{\Gamma }\left(1+\frac{s}{2}\right)u(x).$$
(27)
Further, equation (27) is useful to show that for a positive $`s`$ the asymptotic behavior of $`\mathrm{\Psi }`$ at $`w\mathrm{}`$ becomes

$$\mathrm{\Psi }A_0_{w_0}^{\mathrm{}}Q_s(w^{})𝑑w^{}\times v(x),$$
(28)
where the dependence on $`x`$ is given by

$$v(x)=\underset{n=0}{\overset{\mathrm{}}{}}\frac{(1)^n(4n+3)}{(2n+1)(2n+2)}\frac{(2n+1)!!}{(2n)!!}_{|x|}^1P_{2n+1}(x^{})𝑑x^{}.$$
(29)
By substituting the value of $`s=0`$ into the equality (27), we can calculate the summation to be $`v=1|x|`$. Hence, we conclude that this black hole magnetosphere based on a disk current model has the asymptotically conical shape of magnetic field lines. This asymptotic structure is a consequence of the disk currents extending to infinite distance. Though such a current distribution is unphysical, the existence of open field lines will be useful to discuss the problem of outgoing winds from the disk to infinity in our modeling.

## 3 THE GLOBAL STRUCTURE OF FIELD LINES

Though the asymptotic behavior of $`\mathrm{\Psi }`$ at infinity is clear, the expression written by the infinite sum of the multipole fields is still complicated for understanding the features of the field line structure near the inner edge and the event horizon. Hence, let us assume that the inner edge is far distant from the event horizon. For example, the inner edge of the disk around the Schwarzschild black hole may be located at the innermost stable circular orbit $`r=6M`$ on the equatorial plane. Then, we have $`w_0=5`$, for which the approximation $`w_01`$ will be roughly allowed.

In the following we focus on equation (20), because the multipole sum should arrive at the same result even if equation (23) is used. By virtue of the approximation $`w_01`$, we obtain the dominant behavior of the Legendre function such that $`Q_s(w_0)w_0^{s1}`$. Except in the region near the event horizon where $`w`$ is of order unity, we can calculate equation (20) under the assumption $`1ww_0`$, which leads to the approximate relation

$$P_{2n+1}(w)Q_{2n+1}(w_0)\frac{2n+1}{(2n+2)(4n+3)}\left(\frac{w}{w_0}\right)^{2n+2}.$$
(30)
Of course, for the angular part involving the integral of $`P_{2n+1}(x)`$ we cannot use such an approximation. Hence, our key step is to rewrite the Legendre polynomial through the integral formula

$$P_{2n+1}(\mathrm{cos}\theta )=\frac{1}{\pi }_0^\pi (\mathrm{cos}\theta +i\mathrm{sin}\theta \mathrm{cos}\phi )^{2n+1}𝑑\phi ,$$
(31)
which allows us to obtain

$$\mathrm{\Psi }A_0Q_s(w_0)w_{|x|}^1𝑑x^{}\left[\frac{1}{\pi }_0^\pi 𝑑\phi K(y)\right].$$
(32)
The variable $`y`$ and the function $`K`$ are defined by

$$y=\frac{w}{w_0}\left[x^{}+i\mathrm{cos}\phi \sqrt{1(x^{})^2}\right],$$
(33)
and

$$K=\underset{n=0}{\overset{\mathrm{}}{}}\frac{(1)^n(2n+1)}{2n+2+s}\frac{(2n+1)!!}{(2n+2)!!}y^{2n+1},$$
(34)
for which the infinite sum becomes possible to give

$$K=\frac{1}{y}\left(1\frac{1}{\sqrt{1+y^2}}\right)(s+1)y^{s1}_0^y(y^{})^{s1}\left[1\frac{1}{\sqrt{1+(y^{})^2}}\right]𝑑y^{}.$$
(35)
This expression of $`K(y)`$ is also valid in the outer region $`w>w_0`$, and it is easy to verify that at infinity where $`y1`$ the magnetic stream function given by equation (32) can keep the same asymptotic form as equation (28).

Next, let us give the magnetic stream function near the event horizon. Because $`w`$ is of order unity, the term with $`n=0`$ dominates in the multipole sum (20) under the approximation $`w_01`$, and we have

$$\mathrm{\Psi }\frac{A_0Q_s(w_0)}{4(s+2)w_0}(1x^2)(w+1)^2,$$
(36)
which shows the cylindrical magnetic field written by $`\mathrm{\Psi }=\pi B_0r^2\mathrm{sin}^2\theta `$ with the strength

$$B_0=\frac{A_0Q_s(w_0)}{4\pi M^2(s+2)w_0}.$$
(37)
This is in accordance with the result claimed for general multipole expansions by King, Lasota, & Kundt (1975). However, if $`w_0`$ is of order unity, the dominance of the cylindrical field in the sum (20) apparently breaks down even near the event horizon, and a more complicated structure of field lines will become possible there.

We hereafter fix $`B_0`$ to be positive in equation (37). Then, for $`s>0`$, the magnetic stream function becomes negative at infinity, where we have

$$\mathrm{\Psi }\left(1+\frac{2}{s}\right)4\pi B_0(Mw_0)^2\left(1|x|\right).$$
(38)
Because this change of $`\mathrm{\Psi }`$ occurs for all $`\theta `$, a closed field line of $`\mathrm{\Psi }=0`$ connecting the equatorial plane with the polar axis should exist in the magnetosphere. This critical field line can thread the black hole along the polar axis $`|x|=1`$ and divide the magnetosphere into the inner region with closed field lines and the outer region with open field lines. Note that if equation (32) is written by the sum of the multipole fields we obtain the $`r`$-dependence at $`x=0`$ as follows,

$$\mathrm{\Psi }B_0(Mw)^2\underset{n=0}{\overset{\mathrm{}}{}}\frac{2(s+2)}{2n+2+s}\left[\frac{(2n+1)!!}{(2n+2)!!}\right]^2\left(\frac{w}{w_0}\right)^{2n},$$
(39)
which claims $`\mathrm{\Psi }`$ to be positive everywhere on the inner equatorial plane $`w<w_0`$ between the disk and the black hole. Hence, the critical field line of $`\mathrm{\Psi }=0`$ should thread the disk at some point which we denote by $`w=w_c`$ (or $`r=r_c`$). For example, if we consider the model of $`s=2`$ with $`K`$ of the form

$$K=\frac{3}{y^3}\left(\sqrt{1+y^2}1\right)\frac{1}{y}\left(\frac{1}{2}+\frac{1}{\sqrt{1+y^2}}\right),$$
(40)
the numerical estimation of equation (32) gives $`w_c1.9w_0`$. (The critical field line also reaches to the polar axis at $`w1.8w_0`$, along which $`w`$ is approximately constant.)

To give the magnetic stream function approximately valid in the whole range from $`w1`$ to $`ww_0`$, we modify equation (32) which reduces to $`\mathrm{\Psi }\pi B_0(Mw)^2(1x^2)`$ in the limit $`ww_0`$. We note that $`\mathrm{\Psi }`$ given by equation (32) is smoothly matched to the cylindrical field (36) which is valid near the event horizon, only by multiplying the factor $`(w+1)^2/w^2`$ which becomes equal to unity in the region $`w1`$. Hence, the modified form of the magnetic stream function should be

$$\mathrm{\Psi }=4(s+2)B_0(w+1)^2(w_0/w)M^2_{|x|}^1𝑑x^{}\left[_0^\pi 𝑑\phi K(y)\right]$$
(41)
with $`K`$ given by equation (35).

We have arrived at the main result (41) available as a model for describing the global structure of the black hole magnetosphere. The numerical example of magnetic field lines for $`s=2`$ and $`w_0=10`$ is drawn in Figures 1 and 2, which are useful to see the magnetic connection between the disk and the black hole in the inner region and the conical shape at large $`r`$, respectively. The magnetic spot in the disk connected with the black hole is found to be a very narrow region. This is because $`\mathrm{\Psi }B_0M^2`$ near the event horizon, while $`\mathrm{\Psi }B_0(Mw_0)^2`$ on the disk unless $`ww_c`$. If the width of the magnetic spot is denoted by $`\mathrm{\Delta }w_c`$ (or $`\mathrm{\Delta }r_c`$), at $`w=w_c`$ we obtain

$$B_0M^2\left(\frac{\mathrm{\Psi }}{w}\right)\mathrm{\Delta }w_c,$$
(42)
which leads to $`\mathrm{\Delta }w_c1/w_01`$ (i.e., $`\mathrm{\Delta }r_cM`$). Because the area $`2\pi r_c\mathrm{\Delta }r_c`$ of the magnetic spot is just of order of the horizon area $`16\pi M^2`$, no significant amplification of the magnetic field occurs near the event horizon. Steady inflows of the disk plasma to the black hole may occur along the magnetic field lines, starting from the magnetic spot of the disk which is apart from the inner edge. Then, the inflows reach to the polar region at a distance roughly equal to $`Mw_0`$ and finally fall to the black hole along the cylindrical field lines.

Plasma outflows from the outer part of the disk may also emanate along the open field lines of $`\mathrm{\Psi }<0`$ extending to infinity. It should be remarked, however, that the co-existence of two topologically distinct zones in the magnetosphere is not due to a general relativistic effect of the black hole, because such a global structure is also possible for an appropriate superposition of a dipole and a uniform field, giving closed and open field lines, respectively. The role of the black hole is rather to prohibit the presence of a dipole field at the event horizon. Hence, the closed field lines should be sustained by disk currents only, which is the point essential to our modeling.

In astrophysical magnetospheres various MHD processes are expected to become important for producing highly energetic phenomena. For example, in the region very close to the event horizon, strong gravity forces the plasma to accrete in radial direction and to bend the cylindrical field lines (Punsly, 1991b; Hirotani et al., 1992). The magnetic field (41) generated by disk currents provides a starting point to analyze such MHD effects of ambient and disk plasma on the magnetospheric structure, which will be discussed in the next section.

## 4 DISCUSSION

If the plasma inertia is taken into account even under the condition of magnetic domination, the asymptotic collimation from the conical shape of field lines should occur in logarithmic scales of the cylindrical radius (Chiueh, Li, & Begelman, 1991; Tomimatsu, 1994). Apparently such an inertial back-reaction of the plasma on the poloidal field does not work in the vacuum fields. Then, we present a model with collimated structure, by adding the uniform field (which is a solution of the vacuum Maxwell equations)

$$\mathrm{\Psi }_u=\pi B_ur^2\mathrm{sin}^2\theta $$
(43)
to $`\mathrm{\Psi }`$ given by equation (41). Let us denote the new magnetic stream function by $`\mathrm{\Psi }^{}(=\mathrm{\Psi }+\mathrm{\Psi }_u)`$. Because $`\mathrm{\Psi }`$ is negative in the far distant region, $`B_u`$ should be chosen to be negative for allowing a smooth collimation of the conical field lines. If the uniform field is weak, namely $`B_0+B_u>0`$, the value of $`\mathrm{\Psi }^{}`$ remains positive at the event horizon to keep the magnetic connection between the disk and the black hole (Fig. 3). However, if $`B_0+B_u`$ becomes negative, the field lines threading the black hole can reach to infinity only (Fig. 4).

Of course, this addition of a uniform field is just a mathematical procedure, and the strong collimation shown in Figures 3 and 4 is not a consequence of the physical confinement mechanisms by self-generated toroidal fields and/or by external pressure with a boundary at finite cylindrical radius. It remains quite uncertain whether the details of the field-line shape shown in these figures become approximately valid in magnetically dominated wind regions or not. Nevertheless, it is interesting to note that the vacuum models can successfully describe the transition of the structure of field lines threading the black hole (i.e., from magnetic connection with the disk to that with the remote load), which may occur in more realistic magnetospheres as an observable change of the astrophysical activity.

To discuss the magnetic activity by using the vacuum field $`\mathrm{\Psi }^{}`$ (or $`\mathrm{\Psi }`$) as a background field perturbed by magnetically dominated winds from the disk, let us consider the magnetosphere slowly rotating with the angular velocity $`\mathrm{\Omega }_F`$ of a magnetic field line and the hole’s one $`\mathrm{\Omega }_H`$ (i.e., $`r\mathrm{sin}\theta \mathrm{\Omega }_F1`$ and $`M\mathrm{\Omega }_Ha/4M1`$). If a magnetic field line threading the black hole is closed in the disk at $`r=r_d`$ (see Figs. 1 and 3), the angular velocity $`\mathrm{\Omega }_F`$ is nearly equal to the local Keplerian one, i.e.,

$$\mathrm{\Omega }_F\mathrm{\Omega }_K\left(M/r_d^3\right)^{1/2},$$
(44)
and we have $`r_dM\chi w_0`$, where $`\chi `$ is a numerical factor in the range $`1<\chi 1.9`$ for the $`s=2`$ model. Now the approximation of slow rotation is assured by virtue of the estimation such that $`r_d\mathrm{\Omega }_F1/\sqrt{\beta w_0}1`$. However, the inner edge in realistic accretion disks cannot be too distant from the event horizon (namely, $`w_0`$ cannot be too large), and from the relation

$$\mathrm{\Omega }_H/\mathrm{\Omega }_FM\mathrm{\Omega }_H(\beta w_0)^{3/2},$$
(45)
we expect to be $`\mathrm{\Omega }_H<\mathrm{\Omega }_F`$ as a typical case for slowly rotating black holes. Hence, the inward winds from the magnetic spot carry angular momentum from the disk to the black hole, and the disk plasma which remains in the magnetic spot without flowing into the magnetic tube will be torqued and fall inward. This means that the disk currents located in the inner region and sustaining the loop field lines become unstable for radial accretion, and the closed loop field disconnected from the disk must disappear owing to the infinite redshift effect at the event horizon (Punsly & Coroniti, 1989).

In the magnetically-dominated wind region the toroidal magnetic field $`B_T`$ is approximately given by

$$B_T4\pi \eta L/\alpha \varpi ,$$
(46)
where $`\alpha `$ and $`\varpi `$ are the lapse function and the cylindrical radius, respectively, in Kerr geometry. By virtue of the approximation of slow rotation we obtain $`\alpha =\sqrt{1(2M/r)}`$ and $`\varpi =r\mathrm{sin}\theta `$. Along a poloidal field line the particle flux $`\eta `$ per unit flux tube and the total angular momentum $`L`$ of magnetized fluid are conserved, and the angular momentum flux per unit area carried by the winds towards the black hole is given by $`\eta LB_p`$, which is a consequence of exerting the magnetic torque on the disk plasma. (Now $`L`$ is positive, while $`\eta `$ is negative for ingoing winds.) Then, the characteristic timescale $`t_1`$ to extract the plasma angular momentum in the disk would be

$$t_1h\rho _dr_d^2\mathrm{\Omega }_K/(\eta LB_p)(r_d/h)/\beta _e\mathrm{\Omega }_K,$$
(47)
where $`\rho _d`$ is the mass density in the disk with the thickness $`h`$, and $`\beta _e`$ is the ratio of the external magnetic stress $`B_TB_p`$ at the disk surface to the local pressure in the disk. One may also consider the internal torque transporting angular momentum radially outward through the disk, for which the accretion timescale $`t_2`$ is estimated to be

$$t_2(r_d/h)^2/\beta _i\mathrm{\Omega }_K,$$
(48)
and the ratio $`\beta _i`$ of the internal stress to the local pressure is usually supposed to be less than unity (e.g., $`\beta _i0.1`$). In our modeling the external torque responsible for carrying angular momentum into the winds should be dominant if compared with the internal torque, and we require the condition $`\beta _e\beta _ih/r_d`$, which becomes valid for $`hr_d`$, unless the turbulent internal disk field is supposed to be much stronger than any ordered vertical field at the disk surface. Hence, for $`\beta _e\beta _i`$, the duration of the closed loop structure could be quite longer than the dynamical timescale $`\mathrm{\Omega }_K^1`$ particularly in a thin disk.

In significant long term accretion of plasma beyond the timescale $`t_1`$ one must consider the change of the closed field-line topology, unless dynamo action in the inner disk region can efficiently work to recover the closed loops. The process of disruption of the closed loops would be quasi-stationary, and our models could give a qualitative picture of the structure change as follows: Initially strong disk currents can sustain nearly circular loops of field lines shown in Figure 1 (i.e., $`B_0|B_u|`$). As the disk currents become weaker (i.e., as the ratio $`B_0/|B_u|`$ decreases) by virtue of the radial accretion, the loops begin to collapse towards the black hole (especially in vertical direction), as was shown in Figure 3. Finally, only the loops disconnected from the horizon remains near the inner edge (see Fig. 4 for $`B_0|B_u|`$). The loop field completely annihilates as a set of O-points, if the radial accretion from the inner edge continues still more (Punsly & Coroniti, 1990b). During this process, a part of electromagnetic energy stored in loop field lines could be outward transferred as a Poynting flux to infinity. We expect that the magnetic field (41) to become a convenient initial state for analytical and numerical MHD calculations to study the problem of non-steady energy release in plasma accretion toward the black hole (see, e.g., Koide et al. (2000) for numerical investigations, in which the initial magnetic field is chosen to be the Wald solution in Kerr background).

In the slowly rotating magnetospheres the closed field line plays the role of a path for the energy and angular momentum transfer from the disk to the black hole. \[ Li (2000c) claimed that the transfer from the black hole to the disk occurs only for $`a/M>0.36`$.\] However, the spin-up of the black hole should stop when the structure shown in Figure 4 has been established. Then, the energy extraction from the rotating black hole becomes possible as a mechanism of energy release. The angular velocity $`\mathrm{\Omega }_F`$ of the open field lines connecting the black hole with the remote load may satisfy the condition $`\mathrm{\Omega }_F\mathrm{\Omega }_H/2`$, which allows the twisted field with the negative toroidal component $`B_T`$ (i.e., $`\eta <0`$ and $`L<0`$) to bring negative angular momentum and energy to the black hole. The efficiency of this Blandford-Znajek mechanism is still controversial, because the validity of the condition $`\mathrm{\Omega }_F\mathrm{\Omega }_H/2`$ is based on the assumption that a causal contact between the black hole and the remote plasma holds. In plasma accretion onto the black hole the poloidal inward velocity exceeds the fast magnetosonic speed near the horizon, and no information could propagate outward across some critical surface (Punsly & Coroniti, 1990a). In fact, the black hole-driven winds satisfying $`\mathrm{\Omega }_F\mathrm{\Omega }_H`$ have been found to exist (Punsly, 1998). For such winds the main mechanism of energy extraction from the black hole could be the screw instability, which makes the twisted magnetic tubes unwind and causes a sudden transport of the magnetic free-energy into the plasma’s kinetic energy. For the cylindrical magnetic field shown in Figure 4 we can use the instability condition known as the Kruskal-Shafranov criterion (Kadomtsev, 1966) of the form

$$|B_T|/2\pi \delta >|B_p|/S,$$
(49)
where $`\delta `$ and $`S`$ are the width and the length of the cylindrical magnetic tube. This has been proved by Gruzinov (1999) in the force-free magnetosphere, but it should be noted that the analysis has been limited to the case such that $`\mathrm{\Omega }_F=0`$ and $`SM`$ (i.e., the approximation of no rotation and no gravity). Then, a naive application (Li, 2000b) of the Kruskal-Shafranov criterion to the magnetosphere where $`B_T`$ is roughly estimated to be $`r\mathrm{sin}\theta \mathrm{\Omega }_FB_p`$ seems to be very dubious. Here, the toroidal field is estimated by equation (46), which remains valid even for $`\mathrm{\Omega }_F=0`$, and considering the region well away from the event horizon (i.e., $`\alpha 1`$), we rewrite the instability condition into the form

$$S>|B_p|\delta ^2/(2|\eta L|),$$
(50)
of which the right-hand side is constant if the magnetic flux $`|B_p|\delta ^2`$ in the cylinder is conserved. This means that if a larger angular momentum flux $`|\eta LB_p|\delta ^2`$ is carried away from the black hole by magnetically dominated winds in the cylinder with the width $`\delta M`$, the twisted field lines become screw-unstable in a shorter scale $`S`$. After the disruption of the screw-unstable structure of magnetic field, the rotating black hole would begin to twist again the magnetic field lines, and the energy release from the black hole could quasi-periodically occur.

As was first proposed by Gruzinov (1999), the closed field lines connecting the black hole with the disk may be also screw-unstable. However, for the field lines closed in the disk the toroidal component will be mainly generated by disk rotation with the angular velocity $`\mathrm{\Omega }_K\mathrm{\Omega }_F`$. Then, we must discuss the instability condition by including the effect of $`\mathrm{\Omega }_F`$, which is a problem studied in future works (Matsuoka, Tomimatsu, & Takahashi, 2001).

Finally it should be emphasized that the closed field-line topology presented in this paper is a very transient structure in rapidly rotating plasma-filled magnetospheres. The angular velocity $`\mathrm{\Omega }_F\mathrm{\Omega }_K`$ of closed loops connecting the black hole with the disk will be smaller than the hole’s angular velocity $`\mathrm{\Omega }_H`$, and the energy and angular momentum transfer from the black hole to the disk could spin up the disk plasma (Li, 2000c). Then, as a consequence of this energy extraction from the black hole, outward flows of plasma from the disk would be generated, and owing to the frozen-in condition the field lines threading the black hole are eventually disconnected from the disk. The closed field lines anchored in plasma orbiting in the ergosphere would also open up to infinity by the strong dynamo action, as was conjectured to occur in field-aligned pulsars (Punsly & Coroniti, 1990b). Further, the timescale $`t_2`$ for radial accretion due to internal torque may become shorter than the timescales to open up the loops closed in the disk. Then, the closed loops accrete toward the black hole and annihilate owing to the infinite redshift effect (Punsly & Coroniti, 1990b). These are astrophysically interesting processes of rapid energy release, and we could observe transient flaring states of black hole magnetospheres, if the repeated formation of closed magnetic loops due to strong disk currents in the inner region occurs.

We are grateful to T. Matsuoka for valuable discussion and to an anonymous referee for strong criticism which helped us to improve the paper .


format: text

text:
# Absence of a Slater Transition in The Two-Dimensional Hubbard Model

## Introduction

In this letter we report the study of the metal-insulator transition (MIT) and its relation to antiferromagnetism (AFM) in the two-dimensional (2D) Hubbard model at half-filling. Results from numerical simulations  have convincingly shown that the ground state of the model is an AF insulator, with the Néel temperature constrained to be $`T_N=0`$ by the Mermin-Wagner theorem. However, the nature of the MIT is less clear, as there are two conflicting opinions concerning the MIT and its relation to AFM. The first opinion  is that the strong and weak coupling MI transitions are very different. When the local Coulomb repulsion parameter $`U`$ is larger than the non-interacting band width $`W`$, the ground state is an insulator with a large charge gap $`\mathrm{\Delta }U`$. The MIT occurs well before the onset of magnetism at the temperature $`T_gU`$, and the spin and charge degrees of freedom are decoupled. The superexchange interaction couples the spins with the exchange constant $`\left|J\right|4t^2/U`$, and the spins govern the low-energy physics. This type of MIT, which is purely due to local correlations, is called a Mott transition. In weak coupling, a spin density wave (SDW) instability develops at $`T=0`$ because of the nesting of the Fermi surface. The MIT is the direct consequence of the Brillouin zone folding generated by magnetic ordering. This type of MIT will be referred to as a Slater transition. It is believed that this regime can be well described by the usual many-body weak-coupling approaches. The second opinion, is due to Anderson who has argued that the strong-coupling behavior presented above occurs for both strong and weak coupling so that a Mott gap is present for all $`U>0`$ as in one dimension. As the temperature falls, local moments develop first because of the MIT and then they order so that AFM is the consequence of the MIT, not the converse. Because of this strong interaction, there is no adiabatic continuity between the non-interacting and interacting eigenstates. The conventional renormalizable many-body perturbation theory cannot describe this physics.

The nature of the MIT has been addressed in the Hubbard model in infinite dimensions by applying the dynamical mean-field approximation (DMFA). The unfrustrated model is an example of the first scenario presented above with a Slater AF insulator in weak coupling, and a Mott MIT at strong coupling. The one-dimensional model is an example of the second scenario, since there is a gap in the charge excitation spectrum for all non-zero $`U`$, but no long range magnetic order.

In 2D, the Mermin-Wagner theorem precludes any AFM transition at any finite temperature. Hence, if there is a MIT at finite $`T`$, it cannot be attributed to the Slater mechanism. Finite size simulations (FSS) employing the quantum Monte Carlo (QMC) technique have been applied to study this model quite extensively during the last decade. But the difficulty with FSS is that at low temperatures, the correlation length is greater than the lattice size. Thus the effects of correlations are overestimated for smaller clusters because they are artificially closer to criticality than a system in the thermodynamic limit. This tendency is only reduced by increasing the cluster size, which moves the system in the direction of the thermodynamic limit. The spurious AFM gap that opens at finite $`T`$ in the FSS renders the disentanglement of the MI and the AF transitions practically impossible. The situation is radically different in the dynamical cluster approximation (DCA), which is an extension of the DMFA. In the DCA, the system is already in the thermodynamic limit; however, the DCA restricts correlations to within the cluster length. As the cluster size increases, longer range correlations are progressively included. The FSS and DCA results are complimentary in that the DCA (FSS) always underestimates (overestimates) the effects of non-local correlations. Thus, the DCA will always systematically underestimate the gap formation inherent in the Anderson picture of the MIT. Therefore, if a gap is found in the DCA for a finite cluster, it will persist in the thermodynamic limit.

Using the DCA, we show that a MIT occurs in the 2D Hubbard model at half-filling and finite temperature at both weak and strong coupling regimes. This result is to be contrasted with the weak-coupling approaches which predict that a gap will only exist at $`T=0`$ as a consequence of AFM. We argue that at weak couplings the symbiosis of short-range AFM correlations and local moment formation is the key mechanism for the opening of a gap. The local moments forms at a relatively high $`T`$ leading to a short-range AFM order, that in turn enhances the moments which enhance the AF order, etc. This leads to the destruction of the Fermi liquid state found in the DMFA studies . We discuss the validity of our conclusions drawn from finite clusters at the thermodynamic limit.

## Method

A detailed discussion of the DCA formalism was given in previous publications . The compact part of the free energy is coarse-grained in reciprocal space, projecting the problem onto a finite-sized cluster of $`N_c`$ points embedded in a self-consistently determined host. When $`N_c=1`$, the DCA recovers the DMFA, and when $`N_c\mathrm{}`$, it becomes exact. The cluster problem is solved with a generalized version of the Hirsch-Fye QMC algorithm, and the spectra are analytically continued with the maximum entropy method. We use the DCA to study the single band Hubbard model in a two-dimensional square lattice at half-filling. The model is characterized by $`W`$ and $`U`$. We set $`W=2`$ and we vary $`U`$ from $`U=W/4`$ to $`U=2W`$.

## Results

Jarrell , Georges and Krauth  and Rozenberg et al.  originally obtained a numerically exact solution from the DMFA in the limit of infinite dimensions using QMC. The analysis of these mean-field equations revealed the existence of a transition between a paramagnetic state and an AFM state. The paramagnetic solution was shown to be a Fermi liquid for small values of the interaction $`U`$ and an insulator for large $`U`$. Georges and Krauth identified this transition by observing the behavior of the double occupancy which is given by $`D=n_{}n_{}`$, where $`n_\sigma `$ is the density of spin $`\sigma `$, and of the single-particle Green’s function $`G(\tau )`$. $`D(T)`$ displays a minimum in the Fermi liquid regime. But when the system is an insulator for $`U>U_c`$, $`D`$ always increases with T. We observe the same qualitative behavior as in when $`N_c=1`$.

We will perform the same kind of analysis to study the effect of short-range AFM correlations that are absent in the DMFA. The behavior of $`D`$ as function of $`N_c`$ depends strongly upon $`U`$. In the strong coupling regime, $`D(T)`$ is essentially independent of $`N_c`$ (c.f inset to Fig. 1). This was expected, since for large $`U`$, local fluctuations dominate; they are already captured in $`N_c=1`$. In the weak coupling regime, the minimum found for $`N_c=1`$ flattens progressively as $`N_c`$ increases from $`8`$ to $`16`$. At $`N_c=36`$, there is a downturn in $`D`$ at low $`T`$. By localizing the moments at low temperatures, the system gains free energy by taking advantage of the short range magnetic order. For the value $`U=1`$, we are already deep in the weak-coupling regime of the model. Hence, we expect this behavior to be generic in this region. This finding is similar to a recent finite system QMC simulation .

The density of states (DOS) $`\rho (\omega )`$ (shown in Fig. 2) confirms the destruction of the Fermi liquid quasi-particle peak by short-range AFM correlations. With increasing $`N_c`$, the gap opens fully, and the Hubbard side bands become more pronounced, consistent with the suppression of $`D`$ shown in Fig. 1. Since we work within the paramagnetic equations, the opening of the gap in the weak-coupling regime cannot be attributed to the Slater mechanism.

In Fig. 3 we show a tentative phase diagram obtained with $`N_c=1,8,16,64`$. The boundary line between the metallic and the insulating phases is defined as the point where a true gap opens to the system within our numerical accuracy; the criterion for the gap opening is roughly $`\rho (\omega =0)<1\times 10^2`$. Hence for temperatures just above this line, $`\rho (\omega )`$ displays a pseudogap. As expected, for $`U>W`$, the boundary line is almost independent of the cluster size. For $`U<W`$ however, the $`N_c=1`$ line ends at $`U=U_cW`$. For $`U<U_c`$, a Fermi liquid is stable at low $`T`$. For larger clusters, we don’t find a quasiparticle peak, except when $`UU_c`$, and $`U_c`$ seems to approach zero as $`N_c`$ increases. This seems to indicate that there is no weak-coupling fixed point in the two dimensional Hubbard model at half-filling. However, since an exponential behavior of the gap is expected in this limit, the numerical cost becomes too high at very small couplings; the smallest value of $`U`$ we were able to observe a pseudogap is $`U=0.5`$. For instance, it was necessary that $`N_c>36`$ in order to see a gap at $`U=1`$ (as shown in Fig. 2). We believe that this gap is due to a synergism between short ranged order and moment formation as evidenced by the suppression of $`D`$ (Fig. 1) and the enhancement of the Hubbard side bands (Fig. 2) seen when the gap opens.

It has been argued that the SDW picture might capture the description of the insulating state of the cuprates in the weak to intermediate coupling regime . We show in the inset to Fig. 3 the gap $`\mathrm{\Delta }`$ and the mean-field transition temperature, which to a very good approximation is also the pseudogap temperature $`T^{}`$ characterized by a peak in the bulk susceptibility and the onset of the pseudogap in the DOS. $`T^{}`$ defines the temperature where short-ranged AF correlations emerge. If the gap were due only to these correlations, then we expect $`\mathrm{\Delta }T^{}`$. For example, in the SDW approximation $`\mathrm{\Delta }/T^{}=1.76`$. We find that $`\mathrm{\Delta }/T^{}`$ strongly deviates from this prediction. In fact, we find that $`\mathrm{\Delta }T^{}>T_g`$ hence, SDW ordering is not the mechanism that leads to the gap formation.

We now show that the conclusions drawn from finite clusters are valid in the thermodynamic limit. For this, we exploit the complementarity between the DCA and FSS . In Fig. 4, we show the imaginary-time Green’s function $`G(\tau )`$ at the Fermi point $`X=(\pi ,0)`$. This quantity has a more rapid decay from it maximum at $`G(\beta /2)`$ when the effects of the correlations are stronger so that the gap is more pronounced. In finite systems, the decay is sharper for smaller lattices while in the DCA the opposite occurs. This behavior marks the fundamental difference between the FSS and the DCA. At low temperatures, in FSS, the correlation length is greater than the lattice size. Thus the effect of correlations are overestimated for smaller clusters because these systems are artificially closer to criticality than a system in the thermodynamic limit. This tendency is reduced by increasing the cluster size. The situation is radically different in the DCA where the system is already in the thermodynamic limit. The DCA restricts correlations to within the cluster length. As the cluster size increases, long range correlations are progressively included. Thus, the effects of the correlations increase with cluster size. Since the two techniques must become identical in the thermodynamic limit, the $`G(\tau )`$ curve in the this limit is bracketed by the FSS and DCA curves.

This complementarity is also seen in the spectra. In the weak-coupling regime, the FSS and the DCA start from two different physical limits. A finite system is always an insulator as soon as $`T`$ is less than the energy separation between the ground state and the first excited state. Hence finite systems always show gaps or pseudogaps if $`T`$ is sufficiently low. On the contrary, for $`N_c=1`$, the DCA is identical to the DMFA which at weak couplings yields a local Fermi metal. As the cluster size increases, one expects the two techniques, that become identical in the thermodynamic limit, to converge towards this limit from complementary directions. In the strong coupling regime, the DMFA has a gap too. But this it is always smaller than the one found in the FSS for the same set of parameters. Hence the convergence towards the thermodynamic limit in the DCA will be from smaller gaps to larger, while the converse occurs in FSS. The DOS shown in Fig. 5 supports these conclusions. The finite size gap in the FSS decreases when the cluster size goes from $`N_c=16`$ to $`64`$. In the DCA, there is a pseudogap for $`N_c=16`$ that turns into a true gap when the cluster size is increased to $`64`$. Since by construction the DCA underestimates the gap, we can affirm that at this temperature, the gap exists in the thermodynamic limit. Its actual value is bracketed by the FSS and the DCA.

## Conclusion

We have presented well-controlled results on the MIT in the 2D Hubbard model obtained by the DCA. We find that the double occupancy and the single-particle DOS show evidence of a gap for both small and large values of $`U`$ at finite $`T`$. Since, a finite-temperature gap persists well into the weak coupling regime, and there is no long range order until $`T=0`$, the Slater mechanism is likely not responsible for the metal-insulator transition in the two dimensional Hubbard model. Instead, the gap is due to a synergism between short ranged order and moment formation, as evidenced by the concomitant opening of the gap and the suppression of the double occupancy. Since the charge fluctuations are suppressed at finite $`T`$ by the MIT and the spins do not order until T=0, spin-charge separation also persists into the weak coupling regime. As the DCA systematically underestimates the gap formation, these conclusions are valid in the thermodynamic limit. The resulting phase diagram is consistent with Anderson’s view that the effective Hamiltonian for the 2D Hubbard model at half-filling for all $`U>0`$ and $`\mathrm{\Delta }T`$ is the 2D Heisenberg Hamiltonian. Finally, our results put some constraints on theories of high temperature superconductivity. Since theories based on weak-coupling expansions are unable to describe these results, their prediction in the doped case are questionable.

## Acknowledgments

Part of this work was performed during the visit of S.M. to the University of Southern California, we would like to thank G. Bickers for the warm hospitality and for useful conversations. We gratefully acknowledge discussions with G. Baskaran which motivated this work, we thank C. Huscroft for his help on parallel programming, and we thank M. Ma for helpful discussions. This work was supported by the National Science Foundation grants DMR-0073308 and PHY94-07194 and the Ohio Supercomputing Center.


text:
# 1 Introduction

## 1 Introduction

Since the original formulation of the Sen’s conjecture  there has been significant progress in understanding the nonperturbative aspects of the string field theory. Initially the existence of translationally invariant vacuum with conjectured energy density was established numerically  to a rather high accuracy by the level expansion method  in the Witten’s cubic string field theory . More recently the Sen’s conjecture has been proved rigorously in the framework of background independent string field theory . Nevertheless it seems worth continuing to look for the exact tachyon condensate in the original cubic string field theory since it can teach us many things .

Various insights into the nature of the tachyon condensate has already been obtained in . One particular suggestion for the exact form of the condensate based on the noncommutative geometric nature of the string field action was made in . Independently of this recent suggestion we will show in this letter, that we can obtain some new exact information about the tachyon condensate. We point out that some anomalous symmetries  of the string field theory vertex can be used to derive an infinite set of identities for any string field which solves the equations of motion. We will see that at level $`n`$ we get $`n`$ additional constraints on the coefficients of the string field which are in reasonable agreement with the explicit results from level truncation scheme . It would be very interesting if one could find even further symmetries which would then fix all the coefficients completely.

## 2 Anomalous symmetries

The string field theory action as given by  takes the form of noncommutative Chern-Simons action

$$S[\mathrm{\Psi }]=\frac{1}{\alpha ^{}g_o^2}\left(\frac{1}{2}\mathrm{\Psi },Q\mathrm{\Psi }+\frac{1}{3}\mathrm{\Psi },\mathrm{\Psi }\mathrm{\Psi }\right)$$
(1)
with the noncommutative multiplication defined by

$$\mathrm{\Psi }_1\mathrm{\Psi }_2=bpz\left(V|\mathrm{\Psi }_1\mathrm{\Psi }_2\right),$$
(2)
where $`bpz`$ denotes the $`bpz`$ conjugation in conformal field theory and the vertex $`V|`$ was reviewed in the oscillator formulation  and studied in a background independent manner in .

From  we know that the vertex $`V|`$ satisfies certain identities. For us will be important in particular the following ones for $`n`$ even

$`V|{\displaystyle \underset{i=1}{\overset{3}{}}}(L_n^{(i)}L_n^{(i)})`$ $`=`$ $`3k_n^xV|,`$
$`V|{\displaystyle \underset{i=1}{\overset{3}{}}}(J_n^{(i)}+J_n^{(i)})`$ $`=`$ $`3(h_n^{gh}+3\delta _{n,0})V|,`$ (3)
where $`L_n`$ and $`J_n`$ denote matter Virasoro and ghost current generators respectively. The constants $`k_n^x`$ and $`h_n^{gh}`$ take for $`n`$ even the following values

$`k_n^x`$ $`=`$ $`{\displaystyle \frac{135}{27}}{\displaystyle \frac{n}{2}}(1)^{\frac{n}{2}},`$
$`h_n^{gh}`$ $`=`$ $`(1)^{\frac{n}{2}}.`$ (4)
For $`n`$ odd there would be extra signs between the generators in (2) and the right hand side would vanish. We are not interested in this case since it will not lead to any information about the tachyon condensate. Note that the additional term on the right hand side of the second equation in (2) accounts for the nontensor character of the ghost number current.

Let us study now the variation of the action (1) under the infinitesimal variations of the string field

$`\delta \mathrm{\Psi }`$ $`=`$ $`(L_nL_nk_n^x)\mathrm{\Psi },`$
$`\delta \mathrm{\Psi }`$ $`=`$ $`(J_n+J_nh_n^{gh}3\delta _{n,0})\mathrm{\Psi }`$ (5)
respectively. Under these variations the cubic term in the action is obviously invariant due to the invariance of the vertex (2). On the other hand we know that the total action should also be invariant as long as $`\mathrm{\Psi }`$ satisfies equations of motion. Combining these two facts we get from the kinetic term

$`\mathrm{\Psi }|[Q,L_n]|\mathrm{\Psi }`$ $`=`$ $`k_n^x\mathrm{\Psi }|Q|\mathrm{\Psi },`$
$`\mathrm{\Psi }|[Q,J_n]|\mathrm{\Psi }`$ $`=`$ $`h_n^{gh}\mathrm{\Psi }|Q|\mathrm{\Psi }.`$ (6)
Let us note that both commutators on the left hand side are modes of conformal primary fields, the latter being minus the BRST current $`J^B`$.

## 3 Explicit checks

To compare the above formulas with the results obtained in level expansion scheme in  one should first of all impose the Siegel gauge condition $`b_0|\mathrm{\Psi }=0`$ on the string field and simplify the commutators. For the first equation of (2) one has simply

$$[Q,L_n]=nc_0L_n+\mathrm{}$$
(7)
where the dots stand for terms which do not contribute. For the second equation one can use a little trick. Write the left hand side as

$$\mathrm{\Psi }|[Q,J_n]|\mathrm{\Psi }=\mathrm{\Psi }|\{J_n^B,b_0\}c_0|\mathrm{\Psi }$$
(8)
where we used the facts that $`[Q,J_n]=J_n^B`$ and $`b_0|\mathrm{\Psi }=0`$. The anticommutator can be easily evaluated using the operator product expansion (see e.g. ). Both formulas (2) thus simplify in the Siegel gauge to

$`\mathrm{\Psi }|c_0L_n|\mathrm{\Psi }`$ $`=`$ $`{\displaystyle \frac{1}{n}}k_n^x\mathrm{\Psi }|c_0L_0^{tot}|\mathrm{\Psi },`$
$`\mathrm{\Psi }|c_0(nJ_n+L_n^{tot})|\mathrm{\Psi }`$ $`=`$ $`h_n^{gh}\mathrm{\Psi }|c_0L_0^{tot}|\mathrm{\Psi },`$ (9)
where $`L_n^{tot}`$ denotes the total Virasoro generator. These identities can be easily checked for the numerical values obtained in . Let us define $`r_n^{L,J}`$ to be the ratio of the left and right hand sides of the first or second equation of (3) respectively. Then inserting for simplicity the values for the string field coefficients from  obtained at the level (4,8) we get the following results

$`r_2^L`$ $`=`$ $`1.069,r_4^L=1.044,`$
$`r_2^J`$ $`=`$ $`1.004,r_4^J=0.939.`$
We see that the above identities are preserved within $`7\%`$. This can be compared with the value of the potential which is for the same values about $`1.4\%`$ away from the expected value. This discrepancy in the errors by a factor of five does not necessarily mean that there are mistakes neither in the derivation nor in the numerical evaluation. In fact we know that the convergence properties of the level truncation approximation depends rather strongly on what kind of calculation we are doing. In an unpublished work we have studied the properties of the string field algebra unity $`|I`$ in the level truncation using the universal recursive methods of . Keeping only terms up to level 8 in the unity $`|I`$ and during the whole calculation we got for example

$`L_2|0|I`$ $`=`$ $`0.990L_2|0+0.108L_2^{tot}|00.196L_2^{tot}L_2^{tot}|0+\mathrm{},`$
$`L_2^{tot}|0|I`$ $`=`$ $`0.990L_2^{tot}|0+0.009L_2^{tot}L_2^{tot}|0+\mathrm{},`$ (10)
where the dots stand for terms which are relatively smaller or of higher levels where one can understand bigger errors. Looking at these values one might wonder whether after all the string algebra unity is unity also for the state $`L_2|0`$. The experience from calculations at lower levels where the errors are much bigger suggests that it really converges, hopefully to the correct state. The fact that calculations involving matter Virasoro generators converge much more slowly can be easily traced back to the presence of the Virasoro anomaly.


text:
# Extended Abstract – Model-Based Debugging of Java ProgramsThis work was partially supported by the Austrian Science Fund project P12344-INF and project N Z29-INF. 11footnote 1In M. Ducassé (ed), proceedings of the Fourth International Workshop on Automated Debugging (AADEBUG 2000), August 2000, Munich. COmputer Research Repository (http://www.acm.org/corr/), cs.SE/0011027; whole proceedings: cs.SE/0010035.

## 1 Introduction

Model-based reasoning is a central concept in current research into intelligent diagnostic systems. It is based on the assumption that sources of incorrect behavior in technical devices can be located and identified via the existence of a model describing the basic properties of components of a certain application domain. When actual data concerning the misbehavior of a system composed from such components is available, a domain-independent diagnosis engine can be used to infer which parts of the system contribute to the observed behavior. Model-based Diagnosis provides a set of proven algorithms and methods for searching faults  and identifying points of measurement .

This paper describes the application of the model-based approach to the debugging of Java programs written in a subset of Java. We show how a simple dependency model can be derived from a program, demonstrate the use of the model for debugging and reducing the required user interactions, give a comparison of the functional dependency model with program slicing , and finally discuss some current research issues.

## 2 Model-Based Diagnosis

The model-based approach is based on the notion of providing a representation of the correct behavior of a technical system. By describing the structure of a system and the function of its components, it is possible to ask for the reasons why the desired behavior was not achieved. In the diagnosis community, the model-based approach has achieved wide recognition due to its advantages:

* once an adequate model has been developed for a particular domain, it can be used to diagnose different actual systems of that domain
* the model can be used to search for single or multiple faults in the system without alteration
* different diagnosis algorithms can be used for a given model
* the existence of a clear formal basis for judging and computing diagnoses

Using the standard consistency-based view as defined by Reiter , a diagnosis system can be seen formally as a tuple $`(SD,COMP)`$ where $`SD`$ is a logical theory sentence modeling the behavior of the given system (in our case the program to be debugged), and $`COMP`$ a set of components, i.e., statements. A diagnosis system together with a set of observations $`OBS`$, i.e., a test-case, forms a diagnosis problem. A diagnosis $`\mathrm{\Delta }`$, i.e., a bug candidate, is a subset of $`COMP`$, with the property that the assumption that all statements in $`\mathrm{\Delta }`$ are incorrect, and the rest of the statements is correct, should be consistent with $`SD`$ and $`OBS`$. Formally, $`\mathrm{\Delta }`$ is a diagnosis iff $`SDOBS\{\neg AB(C)|CCOMP\mathrm{\Delta }\}\{AB(C)|C\mathrm{\Delta }\}`$ is consistent. A component not working as expected, i.e., a statement containing a bug, is represented by the predicate $`AB(C)`$.

The basis for this is that an incorrect output value (where the incorrectness can be observed directly or derived from observations of other signals) cannot be produced by a correctly functioning component with correct inputs. Therefore, to make a system with observed incorrect behavior consistent with the description and avoid a contradiction, some subset of its components must be assumed to work incorrectly. In practical terms, one is interested in finding minimal diagnoses, i.e., a minimal set of components whose malfunction explains the misbehavior of the system (otherwise, one could explain every error by simply assuming every component to be malfunctioning). Basic properties of the approach as well as algorithms for efficient computation of diagnoses are described in .

Starting from straightforward work that used a logic program directly as system descriptions , in the last years the use of model-based reasoning (MBR) for debugging of software has been examined in a wider context . All of the approaches have in common that they use a model derived from a program for locating (or, rarely, correcting) a bug. They differ in the considered programming language (ranging from purely logical languages to hardware description languages – in particular VHDL , functional, and finally imperative languages), and the type of model (qualitative , dependency-, or value-based models). The purpose of previous research was to show the applicability of MBR in the software domain by introducing models, often for special purpose languages. Our current work deals with the extension and application of these principles to a mainstream language (Java). This paper presents first results of an implemented debugger prototype using different example programs. The JADE debugger currently implements a functional dependency model that extends our earlier work . The granularity of the debugger, i.e., the elements of a program that are considered to be faulty or not, is currently set to the statement level (instead of individual expression level) for efficiency reasons.

The JADE debugger combines the standard operation modes of diagnosis systems and standard debuggers. First, the program is converted into the dependency representation which is compiled into a logical model. Once the program has been executed actual observations of its behavior can be provided. This behavior together with the logical model is used by the diagnosis engine to compute bug candidates that map back to positions (and statements) within the program to be debugged (see Figure 1).

## 3 Modeling for Debugging

We first give an overview of computing the dependencies for Java programs and then describe the derivation of the system description which is used for diagnosis.

### Dependencies

For the sake of brevity we omit the discussion of object-oriented features such as dynamic binding in this paper and concentrate on the ”basic” imperative features of the language. A more detailed technical description of the basic idea behind the conversion algorithm (excluding external side effects and method calls) can be found in , and a discussion about references and side effects in . In order to compute dependencies, we must consider that the variables occurring in the methods change their values during program execution. This is handled by assigning a unique index to all locations where a variable occurs as target of an assignment (we refer to each such location as an occurrence of the variable. It is the various variable occurrences that dependencies are computed for.

Let $`x`$ and $`y`$ be indexed variable occurrences of a given method $`m`$. We say that $`x`$ depends on $`y`$ iff the value of $`y`$ determines the value of $`x`$ for at least one input vector. This definition is based on earlier work on dependencies in software debugging, e.g., . Beside debugging dependencies are used for verification (see ). Formally, we define a functional dependency as a pair $`(x,M_x)`$, where $`x`$ is a variable occurrence and $`M_x`$ is a set of variable occurrences such that $`x`$ depends on every $`yM_x`$. We can now compute all functional dependencies of a particular statement by determining the functional dependencies of all variables used within the statement.

Whereas the variables correspond to the ports of traditional diagnosis components, the natural choice for components given the abstract nature of the dependency-based representation are statements. The set of diagnosis components can be viewed as a diagnosis system where the connections are formed by the variable occurrences inside the components, i.e., components $`c_i`$ and $`c_j`$ are connected iff one component establishes the functional dependency $`(v_i,M)`$ and the other, $`(w_j,\{\mathrm{},v_i,\mathrm{}\})`$. Variable occurrences $`v_0`$ are inputs of the whole diagnosis system. A variable occurrence $`v_i`$ is an output iff there is no other occurrence $`v_j`$ such that $`j>i`$. Since during conversion indices are always increased, the resulting diagnosis system, i.e., the graph representation, is acyclic.

Computing functional dependencies for Java programs requires compiling each method declared for a class. A method $`m`$ of a Java program is converted by sequentially converting its statements into diagnosis components.

We illustrate the computation of functional dependencies using the example of Figure 2. The functional dependencies for s1=a\*c; are $`\{(s1_2,\{a_1,c_1\})\}`$ because the value of s1 is given by the product of the values of a and c. In summary we obtain the following dependency sets (indices are ignored) for the 5 statements: $`fd(C_4)=\{(s1,\{a,c\})\}`$, $`fd(C_5)=\{(s2,\{b,d\})\}`$, $`fd(C_6)=\{(s3,\{c,e\})\}`$, $`fd(C_7)=\{(f,\{s1,s2\})\}`$, $`fd(C_8)=\{(g,\{s2,s3\})\}`$ where $`C_i`$ denotes the statement in line $`i`$.

### The System Description

After computing all dependencies, we map them to a logical representation, which can be directly used for model-based debugging. For this purpose we assume that the statements are given as a set $`COMP`$, and that for all statements, the functional dependencies are defined. The set of functional dependencies for a statement $`st`$ is written as $`fd(st)`$.

Functional dependencies describe behavior implicitly by describing influences between variables. Instead of speaking about real values, we can only speak about whether a value $`v`$ is correct (written as $`ok(v)`$) or not (written $`nok(v)`$). We can further write that if a statement $`s`$ is assumed to be correct (i.e., $`\neg AB(s)`$ holds) and all input variables have a correct value then the value of variables used as target in an assignment statement must be correct. Formally, the system description is given by:

$$_{(o,M)fd(C)}\left[\neg AB(C)\underset{xM}{}ok(x)ok(o)\right]SD$$

where $`CCOMP`$ is a statement. In addition, we know that it is impossible that a variable value is known to be correct and incorrect at the same time. Therefore, we have to add the rule $`ok(v)nok(v)`$ to the model $`SD`$, for each variable occurrence $`v`$ in the program. The described model can be used together with a standard MBD algorithm for computing bug locations.

For software debugging, the observations required for diagnosis are given by the specified behavior, in our case the expected input/output vectors. By comparing the specified output with the computed output, we can classify the correctness of variables. Variables $`v`$ that are assumed to have the correct value lead to the observation $`ok(v)`$. Variables with an incorrect value are represented by $`nok(v)`$.

In Figure 2(b) the evaluation trace for the call test(3,2,2,3,3) for the Java program given in Figure 2(a) is given. The trace only presents the lines of code which are involved in the current evaluation, and the new environments created. To distinguish different local variables they are indexed with the name of the method where they are declared. In this case there is no return value. From the dependencies computed above, we get the logical model $`SD`$:

| $`\neg AB\left(C_4\right)ok\left(a\right)ok\left(c\right)ok\left(s1\right)`$ | $`\neg AB\left(C_5\right)ok\left(b\right)ok\left(d\right)ok\left(s2\right)`$ |
| --- | --- |
| $`\neg AB\left(C_6\right)ok\left(c\right)ok\left(e\right)ok\left(s3\right)`$ | $`\neg AB\left(C_7\right)ok\left(s1\right)ok\left(s2\right)ok\left(f\right)`$ |
| $`\neg AB\left(C_8\right)ok\left(s2\right)ok\left(s3\right)ok\left(g\right)`$ |  |

| $`ok\left(a\right)nok\left(a\right)`$ | $`ok\left(b\right)nok\left(b\right)`$ |  | $`ok\left(s3\right)nok\left(s3\right)`$ |
| --- | --- | --- | --- |

In this example we assume that the method call test(3,2,2,3,3) should lead to values f=12 and g=0, i.e., that line 8 should be g=s2-s3 instead of g=s2+s3. For this case, we get observations $`OBS`$:

$$ok(a)ok(b)ok(c)ok(d)ok(e)ok(f)nok(g)$$

Using $`SDOBS`$ we get 3 diagnoses, each pinpointing a single possible bug location: $`\{C_5\},\{C_6\},\{C_8\}`$. The other statements can be ignored in this case. Using the measurement selection algorithm from  we can compute the optimal next question to be presented to the user in order to distinguish between the 3 candidates.

## 4 Diagnosing with the Dependency Model

The JADE debugger with dependency model can be proven to be complete with regard to bugs that do not alter the dependency structure of the program, since all statements that may cause a wrong value are considered and therefore are diagnosis candidates. However, discrimination capability can be low. Consider the example program from Figure 2, together with the specified values $`f=12`$ and $`g=0`$. In this case, the debugger returns the candidate $`\{C5\}`$ which could be eliminated when using a value-based model. Now assume $`C5`$ is incorrect, all other statements are correct, and apply the test case from Figure 2. From $`C4`$ and $`C6`$ and the input values we derive $`s1_{test}=s3_{test}=6`$. Using these values together with $`C7`$ and $`f_{test}=12`$, we get $`s2_{test}=6`$. Now using this value together with the assumption $`C8`$ is correct leads to $`g_{test}=12`$, contradicting our specified value $`g_{test}=0`$. Hence, $`\{C5\}`$ is not a diagnosis w.r.t the value-based model, illustrating the (unsurprising) fact that a model based purely on dependencies is too weak to discriminate between all possible program errors.

Concerning performance, for smaller programs and interactive debugger use, diagnosis times should be in the single second range although longer times are acceptable for very large programs. It is obvious that searching for all single bugs using our model is restricted to $`O(n^2)`$, where $`n`$ denotes the number of diagnosis components, i.e., in our case, statements. Using empirical results from  we can expect that computing all single bugs for Java methods with several hundred statements should be done in less than 1 second, a result that is consistent with the experience from the VHDL domain .

Like program slicing , our dependency model is based on static analysis of the code, i.e., it is computed using the program structure and does not use the runtime program behavior for fault localization. A program slice is defined as the part of a program possibly influencing the value of given variables and not occurring within the program after a given position. The slice for our running example for variables $`\{g\}`$ and position 8 comprises the lines 8, 6, and 5. This result is equal to the one obtained by our dependency model and the question arises about the differences between both approaches. Our dependency model is more hierarchically organized, e.g., formally a conditional statement is viewed as a single diagnosis component and sub-divided only after being identified as faulty. To allow a comparison of slicing with MBD using the functional dependency model, we assume an appropriate mapping. In this respect we obtain similar results from both techniques except in the cases where several variables have a faulty value after program execution. In this situation the model-based approach tries to minimize the source of the misbehavior leading to fewer solutions, while slicing does not. Given the example, assume that line 5 is faulty, leading to wrong values for variables $`f`$ and $`g`$. The slice for $`f`$ and $`g`$ is the whole program, while our dependency model would deliver only line 5 as single fault candidate.

## 5 Empirical results

The following experiments show the results of the JADE debugger using various Java methods from our example library, modified at randomly selected statements:
Example 1: The $`adder`$ method implements a binary full adder mapping three inputs to two outputs.
Example 2: The $`library`$ method is part of a small application that creates a sample library and then computes the author who has published the most books of all authors whose books can be found in the library (methods involving object-oriented language structures, such as multiple objects, instance method calls, class & instance variables, etc…)
Example 3: Sorting methods: various sorting methods, providing the full complement of control statements: loops, selection statements, and method calls.

Table 1 shows (from left to right) the tested method (in which a single error has been installed), the total number of statements in each method, the index of the buggy statement within the method (which can directly be used to compare the outcomes of the Jade  debugger tests with ”manual” use of a debugger where the user steps through the code sequentially until the erroneous line is found), and finally the number of user interactions which are needed to exactly locate the bug, classified by type:
Setup: verify the system output connections, i.e. all output variables of the method.
Query: verify the value of a variable at a particular statement
Loop: debug the condition of a loop or selection statement
Exprs: find the smallest sub-expression of a particular statement. This allows further debugging of method and constructor calls.
Iter.: determine the first loop iteration in which the error occurs

The ”Total” column from the right shows the total number of user interactions, and the ”Total2” column the sum of all setup and query interactions. The latter figure determines the debugger’s performance at statement level and can therefore be compared with Column 4. These two columns give a simple comparison of the user interactions needed to find the exact location of a faulty statement with the two different debugging strategies, i.e. model-based vs. unsupported debugging.

Locating a bug in the $`adder`$ method with a traditional debugging tool requires 10 user interactions versus 4.43 with the JADE debugger. This fairly drastic difference in favor of the model-based technique is due to the simple block structure of the method with the model-based debugger exploiting its knowledge about the underlying functional dependency structure of the block. The longer the block, the greater the advantage. For the same reason, the faulty version of the $`library`$ method requires 28 user interactions in the traditional approach compared to 6 using the JADE debugger. To locate a bug in a sorting algorithm a traditional debugger on average needs 5.45 steps versus 4.12 for the JADE debugger, the advantage being less clear due to the complex control structures.

Overall, 24 tests are compiled in the table with 9.04 user interactions on average required in the traditional case and 4.37 for the JADE debugger, clearly highlighting the potential abilities of model-based debugging tools.

Note that as with techniques like program slicing, the bug location efficiency of the model-based approach does not depend on the actual location of the error in the code. The number of user interactions it takes to find the bug is therefore much better assessable.

## 6 Discussion and Conclusion

One of the main advantages of the model-based approach is the ability to incorporate multiple models in the same formal and computational framework. The dependency-based representation described in this paper has the advantage of being simple and computationally efficient to solve, thus providing an answer to the important question of scalability of the approach (our use of a dependency based mechanism in the VHDL domain was usable for very large programs ). Its disadvantage of low discrimination can be counteracted by combining it with a more detailed (i.e., value-based) model, e.g., by focusing in on smaller program parts isolated by the dependency model, and diagnosing them with the more detailed model, which is currently undergoing implementation . The value-based model, unlike the dependency model, is dynamic and based on the evaluation trace, similar to the manner in which Dynamic Program Slicing  extends classical Program Slicing. The value-based model records, in a detailed manner, which program parts actually contribute to the faulty behavior. By propagating value assumptions forward and backward via the standard diagnosis algorithms it provides better discrimination. In  an approach for combining program slicing and algorithmic debugging for debugging of procedural programs was introduced. The ideas of using test specifications and test results can be easily incorporated into our approach. In contrast to  we can change the debugging performance by changing the underlying model without changing the underlying algorithms.

The JADE debugger is a prototype system for research purposes and for demonstrating the underlying model-based techniques. Its main application interface can be seen in Figure 3. The system is not yet applicable in a real production environment. The incorporated Java evaluator that is used for computing variable values to be presented to the user during debugging implements only the basic Java functionality and some important classes and is far away from being JDK compliant. In addition, JADE assumes that the Java source code for all involved classes is available. Both problems have to be tackled for a production version of JADE.

We have described the application of model-based technology in the building of a prototype intelligent debugger for Java programs. The approach is based on automatically building a formal internal model of the executed program, which a generic diagnosis engine then uses, together with observations of incorrect program output, for identifying possible sources of the error in the program. This combines modeling flexibility with the ability to reuse standard algorithms like measurement selection. Following an approach that we have followed in the domain of hardware design languages, the currently used model is purely dependency based for simplicity and quick diagnosis, and will be combined with a more detailed value-based model for improved discrimination capabilities.


text:
# First results from the BOOMERanG experiment

## Introduction

Acoustic oscillations in the primeval plasma before recombination ($`z1100`$) produce a delicate pattern of horizon and sub-horizon size structures in the Cosmic Microwave Background. In the framework of inflationary adiabatic perturbations the angular power spectrum of these structures features a harmonic series of peaks, with the first peak at multipole $`\mathrm{}200`$.

Here we present a high confidence measurement of the first peak resulting from the analysis of the multiband, high quality image of about 3$`\%`$ of the sky obtained in the long duration flight of BOOMERanG deB00 ; Lan00 . We show that the BOOMERanG image is a faithful representation of the CMB at angular scales smaller than 10o, and we shortly discuss the methods we used to make sure that instrumental and astrophysical contaminations are negligible in our results. We present the measured power spectrum, and we discuss the cosmological implications of this measurement.

## The instrument and the data set

The BOOMERanG payload and the important subsystems have been described in Mau97 ; Mas98 ; Mas99 ; Pia00 ; Cri00 . It is a scanning telescope measuring simultaneously eight pixels in the sky. Four pixels feature multiband photometers (150, 240 and 410 GHz), two pixels have single-mode, diffraction limited detectors at 150 GHz and two pixels have single-mode, diffraction limited detectors at 90 GHz. The NEP of these detectors is below 200 $`\mu K_{CMB}\sqrt{s}`$ at 90, 150, and 240 GHz, and the angular resolution ranges from 10 to 18 arcmin FWHM.

The continuous scan of the sky ($`\pm 30^o`$ in azimuth at $`1^o/s`$ or $`2^o/s`$ ) allows a wide sky coverage. Moreover, it encodes different spherical harmonic components of the sky brightness into different frequencies within a detection bandwidth ($``$ 0.2-20 Hz) free from 1/f noise and drifts. The instrument is flown aboard a stratospheric balloon at 38 km of altitude, to avoid the bulk of atmospheric emission and noise. Using a long duration balloon flight carried out by NASA-NSBF around Antarctica we mapped $`1800`$ square degrees in the region of the sky where Galactic contamination is minimal Sch98 . The instrument collected 57 million samples for each detector during the scans. A pointing solution matrix (accurate to 3’ rms) has been constructed using the data from the attitude sensors (sun sensors, laser gyroscopes, differential GPS). The beam response has been calibrated at ground, using a tethered thermal source in the far field of the telescope. This calibration has been checked in flight observing compact HII regions. The responsivity has been calibrated against the CMB dipole, which is visible along the scans as a $``$ 3 mK linear drift. Systematic effects limit the precision of this calibration to 10$`\%`$. For further analysis the data have been high-passed in the time domain, so that 1/f noise and drifts are reduced. In the same process, all the structures in the sky larger than 10o are effectively removed. A time-time correlation matrix has been constructed using an iterative method which effectively separates signal and noise in the data Pru00 . The maximum likelihood map and its covariance have been obtained from the data vector, the time-time correlation matrix and the pointing matrix as described in Bor99 . We use the HEALPIX Gor98  pixelization with 14’ pixel side as a compromise between coverage of high multipoles and computation speed.

## The maps and the power spectra

The maps in the four frequency channels have been published in deB00 . In figure 1 we increase the contrast and show the central part of a map obtained using a naive coadding in pixels of data from three of the 150 GHz channels. The 90, 150, 240 GHz maps all present the same degree-size structures, with amplitude ratios characteristic of CMB anisotropies (see figure 2), thus demonstrating that the bulk of the anisotropy we measure is cosmological in origin.

Dust contamination is evident only at lower Galactic latitudes, and mainly in the 240 GHz channel. At high Galactic latitudes the 410 GHz channel is a good monitor of low level emission from interstellar dust. We find a good correlation between this map and the IRAS/DIRBE map extrapolated at 410 GHz (model 8 in ref. Sch98 ) and filtered in the time domain in the same way as our data. This correlation allows us to estimate low level contributions from IRAS-correlated dust to the CMB map shown in fig.1. We find (see deB00 ; Mas00  for details) that this dust can produce at most 2$`\%`$ of the mean square anisotropy plotted in the 150 GHz map. We also studied the effect of point sources in the observed region. We used the PMN survey and the WOMBAT Wom00  software to extrapolate the radio sources flux to 150 GHz. The $`rms`$ brightness fluctuation due to all the sources in our field corresponds to $`7\mu K_{CMB}`$, and has to be compared to the total measured sky fluctuation $`80\mu K_{CMB}`$. Its contribution to the angular power spectrum has been computed and is reported in fig.3. The maximum likelihood power spectrum of the central part of the map (1$`\%`$ of the sky) is computed as described in Bor99 ; deB00  and is plotted in fig.3. A distinct peak at multipole $`\mathrm{}200`$ is evident, and significant power is detected at all the multipoles up to the highest observed ($`\mathrm{}=625`$).

In figure 4 we plot the rms signal in the center part of the maps at 90, 150, 240 and 410 GHz, computed as $`\mathrm{\Delta }T_{rms}=\sqrt{_{\mathrm{}}(2\mathrm{}+1)c_{\mathrm{}}/4\pi }`$. The $`\mathrm{\Delta }T_{rms}`$ have been converted into brightness fluctuations using the measured spectral efficiency of the different freqeuency channels. The $`rms`$ of the 410 GHz brightness has been estimated from the component correlated with the IRAS/DIRBE maps; the amplitude of the thermal dust spectrum has been adjusted to fit such $`rms`$. The amplitudes of the other theoretical spectra for Galactic emission have been scaled from DeO00  taking into account the different window functions of BOOMERanG and COBE and assuming an $`\mathrm{}^3`$ power spectrum. Once again it is evident that the data from the three channels more sensitive to the CMB are fit very well by the derivative of a 2.726K blackbody, while all the other reasonable spectra for millimeter wave emission originated in our galaxy cannot fit our measurements.

## Cosmological interpretation of the data

Having been able to exclude all the suspect sources of local contamination in the data of the power spectrum shown in figure 3, we can reasonably assume in the following that the detected power is in the CMB and comes from structures at the recombination epoch ($`z1100`$). In the framework of adiabatic inflationary models, acoustic oscillations in the primeval plasma at horizon and sub-horizon scales produce a distinctive series of peaks in the power spectrum of the CMB. The BOOMERanG data support this picture. The position, amplitude and width of the peak evident in fig.3 are consistent with the general adiabatic inflationary scenario Hu97 , while the simplest models based on topological defects do not fit the data as well xxx . Using a quadratic fit we find that the peak is located at multipole $`\mathrm{}=(197\pm 6)`$ (1$`\sigma `$). This result can shift a little bit if one assumes skewed power spectra similar to the adiabatic scenario ones Kno00 . This position is consistent with a flat geometry of space, but it is not univocally related to the density parameter $`\mathrm{\Omega }`$ if we allow for a non vanishing cosmological constant Wei00 . The continuous line in fig.3, a very good fit to the measured data, has parameters ($`\mathrm{\Omega }`$, $`\mathrm{\Omega }_b`$, $`\mathrm{\Omega }_\mathrm{\Lambda }`$, $`h`$, $`n_s`$) = (1.035, 0.06, 0.425, 0.75, 1.0) where $`\mathrm{\Omega }_b`$, $`\mathrm{\Omega }_\mathrm{\Lambda }`$, $`\mathrm{\Omega }`$ are the density parameters for baryons, cosmological constant, and for the total mass-energy density respectively; $`h`$ is the dimensionless Hubble constant, $`n_s`$ is the spectral index of the spectrum of primordial density perturbations. A full Bayesian analysis has been carried out Lan00  in order to measure the cosmological parameters and compute significant confidence intervals. In order to do this, we need to specify the prior distribution assumed for each of the parameters. This is especially important in the case of CMB power spectrum measurements, since an important geometrical degeneracy is present Bon99 , so different combinations of the parameters produce very similar power spectra. Given a flat model with $`\mathrm{\Omega }_\mathrm{\Lambda }0.7`$, a closed model producing the same power spectrum can be found decreasing $`\mathrm{\Omega }_\mathrm{\Lambda }`$ and $`h`$ and increasing $`\mathrm{\Omega }_b`$, and an open model can be found doing the reverse. The marginalized likelihood curve for $`\mathrm{\Omega }`$ will then depend on the density of models found in each direction along that particular path. Our 95$`\%`$ confidence intervals for $`\mathrm{\Omega }`$ range from $`(0.881.12)`$ to $`(0.971.35)`$ depending on the assumed priors and parametrizations deB00 ; Lan00 . This strongly suggests a flat geometry of the Universe, and at least implies, with $`95\%`$ confidence, a curvature length larger than 2 to 3 times the Hubble length today. This measurement has very important cosmological consequences. Several independent measurements point to a matter density parameter significantly below unity Don99 ; Bah00 ; Bla99 ; Jus00 ; Wit00 . As a consequence a different form of energy is needed to fill the gap between the total energy density of the Universe and the energy density in matter. A cosmological constant is a good candidate, and has the appealing feature of explaining the recent evidence coming from the observation of distant supernovae Rie98 ; Per99  for an acceleration in the expansion rate of the Universe . The fine tuning problem characteristic of the cosmological constant Wei89  can be overcome invoking quintessence Ost95 ; Cal98 ; Arm00 ; Ame00 ; Bal200 . In the same adiabatic perturbations framework, we constrain $`n_s`$ in the range $`(0.711.01)`$ to $`(0.861.24)`$ (95$`\%`$ confidence) again depending on the priors Lan00 . Our constraints on $`\mathrm{\Omega }`$ and $`n_s`$ are consistent and support the simplest inflationary scenarios Lin81 ; Alb82 ; Wat00 . The third parameter constrained by BOOMERanG is the density of baryons. In the power spectrum this parameter controls the relative amplitude of the first peak to the second one (see e.g. Hu97 ). From the power spectrum of fig.4, this ratio is larger than expected for a standard primordial nucleosynthesis ($`\mathrm{\Omega }_bh^2=(0.019\pm 0.002)`$, 1$`\sigma `$ Tyt00 ), thus suggesting a physical density of baryons somewhat higher. Depending on the assumed priors, we obtain 95$`\%`$ intervals for $`\mathrm{\Omega }_bh^2`$ ranging from $`(0.0190.045)`$ to $`(0.0260.048)`$ Lan00 . This should be considered as a spectacular agreement, since we are comparing the density of baryons inferred from the physics of acoustic oscillations $``$ 300000 years after the big bang to the density of baryons inferred from the physics of nuclear reactions a few minutes after the big bang. Moreover, the measurement techniques are completely orthogonal and subject to completely different systematic effects. Should more precise measurements imply a more significant disagreement, new physics will be needed, and several hypothesis have been proposed already aa ; bb ; cc ; dd ; ee ; ff ; gg . One way to improve the $`\mathrm{\Omega }_b`$ consistency is to force the $`n_s`$ value to be smaller. $`n_s0.9`$ is not the best marginalized value for $`n_s`$, but the fit to the data is very good and still consistent with inflation (see e.g. Kin00 ). In fig.5 we plot the BOOMERanG constraints in the ($`\mathrm{\Omega }_m`$,$`\mathrm{\Omega }_\mathrm{\Lambda }`$) plane. Here $`\mathrm{\Omega }_m=\mathrm{\Omega }_b+\mathrm{\Omega }_c`$, with $`\mathrm{\Omega }_c`$ density parameter for non relativistic dark matter. It is evident that the data are fully consistent with a flat geometry of the Universe, represented by the $`\mathrm{\Omega }_m+\mathrm{\Omega }_\mathrm{\Lambda }=1`$ line; however, the CMB anisotropy data alone cannot distinguish between the different contributions to the total mass-energy density. This degeneracy is efficiently removed once additional data are considered, like the data from distant supernovae or the data from the large scale distribution of galaxies. The combination of our data with any of the two is quite powerful, and suggests a ’concordance’ model with $`\mathrm{\Omega }_\mathrm{\Lambda }0.7`$ and $`\mathrm{\Omega }_m0.3`$ (see Lan00  for details). In ref. Teg00  all the available CMB power spectrum data and the power spectrum of the distribution of IRAS galaxies from the redshift survey data are fit simultaneously, confirming the results above.

## Conclusions

BOOMERanG has produced a faithful, wide map of CMB anisotropy at angular scales smaller than 10o. The power spectrum derived from the central part of the map ($`1\%`$ of the sky) is consistent with the inflationary cosmological model with gaussian adiabatic density fluctuations. A careful analysis of the power spectrum data, rigorously taking into account the priors and the degeneracies related to these parameters, has been carried out. We conclude that the geometry of our universe is close to flat, and that the primordial density flucutations have a nearly Harrison-Zeldovich power spectrum. The density of baryons inferred from these measurements is slightly higher than the standard nucleosynthesis value. Combining the BOOMERanG data with data from the observations of distant supernovae or from large scale structure studies we obtain significant determinations of both $`\mathrm{\Omega }_m=(0.48\pm 0.13)`$ and $`\mathrm{\Omega }_\mathrm{\Lambda }=(0.66\pm 0.07)`$ (1$`\sigma `$). All these results come from a single detector, and use the central part of the observed region. We are currently working on the combined analysis of all the 12 detectors sensitive to CMB.

## Acknowledgments

The BOOMERanG project has been supported by PNRA, Universitá “La Sapienza”, and ASI in Italy, by NSF and NASA in the USA, and by PPARC in the UK. We would like to thank the entire staff of the NSBF, and the US Antarctic Program personnel in McMurdo for their excellent preflight support and a marvelous LDB flight. DoE/NERSC provided the supercomputing facilities. Web sites: (http://oberon.roma1.infn.it/boomerang) and (http://www.physics.ucsb.edu/$``$boomerang).