text:
# 1 Introduction

## 1
Introduction

In 1948 Casimir made the remarkable prediction that there is an attractive force between a pair uncharged parallel plane perfect conductors . Furthermore, he argued that this force arises solely from a shift in the energy of the vacuum state of the quantized electromagnetic field. An early attempt by Sparnaay  to observe this force was inconclusive, but in recent years several new experiments  have been performed which seem to give good agreement with Casimir’s prediction. (To be more precise, most of these experiments actually measure the force between a plate and a sphere and incorporate a theoretical correction to compare to Casimir’s result. Of the recent experiments, only that of Bressi et al  uses two parallel plates.)

If the energy of the vacuum state is zero in the limit of infinite plate separation, then the attractive force found by Casimir would seem to imply a negative vacuum energy at finite separation. In fact, Brown and Maclay  showed that for perfectly conducting plates, one has a constant negative vacuum energy density. This conclusion is of great theoretical interest, because negative energy density has the potential to cause some rather bizarre effects in gravity theory. (See, for example, Ref.  and references therein.) However, questions have been raised as to whether the negative energy density will still arise in a more realistic treatment in which the plates are not perfect conductors . In particular, Helfer and Lang  calculated the energy density outside of a single half-space filled with a nondispersive dielectric material and obtained a positive result. They interpreted this as a positive self-energy density associated with a single plate which would add to the negative interaction energy density between a pair of plates. Helfer and Lang conjecture that the net Casimir energy density might be positive when the self energy is accounted for. If this conjecture is correct, then the situation would be analogous to that of the energy density in classical electrostatics. A pair of oppositely charged particles have a negative interaction energy, but the net energy density, which is proportional to the square of the electric field, is always positive.

However, the Helfer and Lang calculation does not include dispersion, which is essential in a realistic treatment. Numerous authors, beginning with Lifshitz , have studied the effects of dispersion upon Casimir forces. However, these authors have been concerned with the force or the total energy, and not the local energy density. The purpose of this paper is to present a calculation of the Casimir energy density in a model in which dispersion is included. For this purpose, we will use the methods of source theory developed by Schwinger and coworkers . This is a method based upon the calculation of Green’s functions which is especially well suited to dissipative materials, and was used by Schwinger et al  to rederive the results of Lifshitz. Milonni and Shih  have used conventional quantum electrodynamics to reproduce some of the results of source theory. There has also been considerable interest in recent years in quantization of the electromagnetic field inside dissipative materials using operator methods . The relation between the results of the latter set of authors and those of Schwinger et al has not yet been clarified.

The outline of this paper is as follows: In Sect. 2 we review the source theory approach as applied to parallel interfaces of dielectric media. In Sect. 3 we compute the expectation values of the squares of the electric and magnetic fields in the vacuum region outside a half-space filled with a uniform dispersive dielectric. We extend this calculation to the case of two parallel dielectric half-spaces and also discuss the energy density in Sect. 4. Conclusions are given in Sect. 5.

## 2 Green-Function Approach for Multilayer Dielectrics

This section is a review of the formalism of Schwinger et al . One begins by writing the Maxwell equations for the macroscopic electromagnetic fields produced by an external polarization source $`𝐏`$, which formally describes the zero point fluctuations of the fields333 Heaviside-Lorentz units with $`c=\mathrm{}=1`$ will be used in this paper. Also, it is assumed that the magnetic permeability is unity.

$$𝐁=\times 𝐀,$$

$$𝐄=\dot{𝐀}\varphi ,$$

$$\times 𝐁=ϵ\dot{𝐄}+\dot{𝐏},$$
(1)
$$(ϵ𝐄+𝐏)=0,$$

where $`ϵ`$ is the dielectric constant of the medium. The wave equation for the electric field resulting from the Maxwell equations is

$$\times (\times 𝐄)ϵ\ddot{𝐄}=\ddot{𝐏},$$
(2)
By assuming a linear relation between sources and fields, the electric field can be written as a spacetime integral

$$𝐄(x)=d^4x^{}\stackrel{}{𝚪}(x,x^{})𝐏(x^{}),$$
(3)
where $`x=(t,𝐫)`$, $`x^{}=(t^{},𝐫^{})`$ and $`\stackrel{}{𝚪}`$ is a Green’s dyadic, which satisfies (2) with a $`\delta `$-function source. Let

$$\stackrel{}{𝚪}(𝐫,𝐫^{},\omega )=_{\mathrm{}}^{\mathrm{}}𝑑\tau e^{i\omega \tau }\stackrel{}{𝚪}(x,x^{}),$$
(4)
where $`\tau =tt^{}`$. From (2) and (3), it follows that $`\stackrel{}{𝚪}(𝐫,𝐫^{},\omega )`$ satisfies the following equation:

$$\times (\times \stackrel{}{𝚪})+\omega ^2ϵ\stackrel{}{𝚪}=\omega ^2\stackrel{}{\mathrm{𝟏}}\delta (𝐫𝐫^{}).$$
(5)
So far, the discussion has been purely classical. At this point, Schwinger et al  use source theory to identify the Green’s dyadic $`\stackrel{}{𝚪}`$ with an “effective product of electric fields”

$$\frac{i}{\mathrm{}}E_j(𝐫)E_k(𝐫^{})=\mathrm{\Gamma }_{jk}(𝐫,𝐫^{},\omega ).$$
(6)
We can interpret this as the Fourier transform of the electric field correlation function. From the Maxwell equation $`\times 𝐄=\dot{𝐁}`$, one finds the corresponding expression for the magnetic field:

$$\frac{i}{\mathrm{}}B_j(𝐫)B_k(𝐫^{})=ϵ_{jlm}ϵ_{knp}(_l_n^{}/\omega ^2)\mathrm{\Gamma }_{mp}(𝐫,𝐫^{},\omega ).$$
(7)
Note that $`\mathrm{}`$ makes its first appearance in these expressions. These expressions can be identified with the vacuum expectation values of products of field operators, which appear in the more conventional field theory approach to quantization of the electromagnetic field. From now onward, we revert to units in which $`\mathrm{}=1`$. In order to calculate the field correlation functions, one needs to find the Green’s function $`𝚪`$ occurring in (5). This amounts to solving a classical boundary value problem.

The interfaces between the media are chosen to be perpendicular to the $`z`$ direction, so for now it will only matter that the dielectric constant changes in the $`z`$ direction only. Therefore, it is convenient to introduce a transverse spatial Fourier transform

$$\stackrel{}{𝚪}(𝐫,𝐫^{},\omega )=𝑑𝐤_{}\frac{1}{(2\pi )^2}e^{i𝐤_{}(𝐫𝐫^{})_{}}\stackrel{}{𝚪}(z,z^{},𝐤_{},\omega ),$$
(8)
where the vector $`𝐤_{}`$ can be chosen to point along the $`+x`$ axis ($`k=\left|𝐤_{}\right|`$).

Some components of $`\stackrel{}{𝚪}`$are found to be

$`\mathrm{\Gamma }_{xx}`$ $`=`$ $`{\displaystyle \frac{1}{ϵ}}\delta (zz^{})+{\displaystyle \frac{1}{ϵ}}{\displaystyle \frac{}{z}}{\displaystyle \frac{1}{ϵ^{}}}{\displaystyle \frac{}{z^{}}}g^B,`$ (9a)
$`\mathrm{\Gamma }_{yy}`$ $`=`$ $`\omega ^2g^E,`$ (9b)
$`\mathrm{\Gamma }_{zz}`$ $`=`$ $`{\displaystyle \frac{1}{ϵ}}\delta (zz^{})+{\displaystyle \frac{k^2}{ϵϵ^{}}}g^B,`$ (9c)
$`\mathrm{\Gamma }_{xz}`$ $`=`$ $`i{\displaystyle \frac{k}{ϵϵ^{}}}{\displaystyle \frac{}{z}}g^B,`$ (9d)
$`\mathrm{\Gamma }_{zx}`$ $`=`$ $`i{\displaystyle \frac{k}{ϵϵ^{}}}{\displaystyle \frac{}{z^{}}}g^B,`$ (9e)
where $`ϵ^{}=ϵ(z^{})`$, and $`g^E`$, the “transverse electric”, and $`g^B`$, the “transverse magnetic” Green’s functions satisfy

$`\left[{\displaystyle \frac{^2}{z^2}}+k^2\omega ^2ϵ\right]g^E(z,z^{})`$ $`=`$ $`\delta (zz^{}),`$ (10a)
$`\left[{\displaystyle \frac{}{z}}{\displaystyle \frac{1}{ϵ}}{\displaystyle \frac{}{z}}+{\displaystyle \frac{k^2}{ϵ}}\omega ^2\right]g^B(z,z^{})`$ $`=`$ $`\delta (zz^{}).`$ (10b)
By introducing the quantity

$$\kappa ^2=k^2\omega ^2ϵ,$$
(11)
(2) can be written as:

$`\left[{\displaystyle \frac{^2}{z^2}}+\kappa ^2\right]g^E(z,z^{^{}})`$ $`=`$ $`\delta (zz^{}),`$ (12a)
$`\left[{\displaystyle \frac{}{z}}{\displaystyle \frac{1}{ϵ}}{\displaystyle \frac{}{z}}+{\displaystyle \frac{\kappa ^2}{ϵ}}\right]g^B(z,z^{})`$ $`=`$ $`\delta (zz^{}).`$ (12b)
So, in order to find the field correlation functions as defined in (6) and (7) in a given situation, one needs to solve these equations with the appropriate boundary conditions. We consider here two cases.

## 3 One Interface Case

We now specialize the above discussion to a situation in which the inhomogeneity of the dielectric constant is due to a plane interface separating a dielectric substance from a vacuum:

$`z>0:ϵ(z)`$ $`=`$ $`1,`$
$`z<0:ϵ(z)`$ $``$ $`ϵ_d.`$ (13)
Here $`ϵ_d`$ is a function of frequency, but not of position.

### 3.1 Boundary Conditions

In solving (12a) and (12b), we use the following boundary conditions. At $`z=z^{}`$, $`g`$ is continuous but the derivative is discontinuous at this point :

$$\frac{g}{z}|_{zz_{}^{}}^{zz_+^{}}=1.$$
(14)
At the boundary ($`z=0`$) we use the conditions for continuity of $`E_x`$, $`E_y`$, $`ϵE_z`$, and $`B_i`$. The first three, as seen from (6), imply the continuity of $`\mathrm{\Gamma }_{xx}`$, $`\mathrm{\Gamma }_{yy}`$, and $`ϵ\mathrm{\Gamma }_{zz}`$ and subsequently, from (2), the continuity of $`g^E`$, $`g^B`$, and

$$\frac{1}{ϵ}\frac{}{z}\frac{1}{ϵ^{}}\frac{}{z^{}}g^B.$$

The continuity of $`B_x`$implies that of $`_z_z^{}\mathrm{\Gamma }_{yy}`$, as seen from Eq. (24), which is given below. From this, using (6) and (9b), we deduce the continuity of $`g^E/z`$.

The solutions $`g^E`$ and $`g^B`$ in the vacuum region have the form

$`g^E`$ $`=`$ $`{\displaystyle \frac{e^{\kappa _0zz^{}}+re^{\kappa _0(z+z^{})}}{2\kappa _0}},`$ (15a)
$`g^B`$ $`=`$ $`{\displaystyle \frac{e^{\kappa _0zz^{}}+r^{}e^{\kappa _0(z+z^{})}}{2\kappa _0}},`$ (15b)
where

$`r`$ $``$ $`{\displaystyle \frac{\kappa _0\kappa _1}{\kappa _0+\kappa _1}}`$ (16a)
$`r^{}`$ $``$ $`{\displaystyle \frac{\kappa _0ϵ_d\kappa _1}{\kappa _0ϵ_d+\kappa _1}}.`$ (16b)
Here $`\kappa _0`$ and $`\kappa _1`$ represent the quantity $`\kappa `$ as defined in (11) for the vacuum region $`(ϵ=1)`$ , and for the dielectric half-space region $`(ϵ=ϵ_d)`$ , respectively, and $`r`$ and $`r^{}`$ can be identified as reflection coefficients for two polarization states, $``$ and $``$ respectively, corresponding to electric field vector being perpendicular or parallel to the plane of incidence of an linearly polarized electromagnetic wave .

### 3.2 The Electric Field

Using (6), we write the formal expectation value of the square of the electric field at coincident points as

$$E^2_f=i_{\mathrm{}}^{\mathrm{}}𝑑\omega \frac{1}{2\pi }_0^{\mathrm{}}𝑑kk\frac{1}{2\pi }\mathrm{\Gamma }_{kk}=\frac{i}{2\pi ^2}_0^{\mathrm{}}𝑑\omega _0^{\mathrm{}}𝑑kk\mathrm{\Gamma }_{kk}.$$
(17)
In the second step, we assumed that the integrand is an even function of $`\omega `$. By complex rotation $`(\omega i\zeta )`$, this becomes:

$$E^2_f=\frac{1}{2\pi ^2}_0^{\mathrm{}}𝑑\zeta _0^{\mathrm{}}𝑑kk\mathrm{\Gamma }_{kk}.$$
(18)
Note from (11) that $`\kappa ^2>0`$ when $`\omega `$ is imaginary. By means of (2), all of the components of $`\stackrel{}{𝚪}`$ in a given region can be written in terms of $`\mathrm{\Gamma }_{xx}`$ and $`\mathrm{\Gamma }_{yy}`$:

$`\mathrm{\Gamma }_{xz}(z,z^{})`$ $`=`$ $`{\displaystyle \frac{ik}{\kappa ^2}}{\displaystyle \frac{}{z^{}}}\mathrm{\Gamma }_{xx}(z,z^{}),`$
$`\mathrm{\Gamma }_{zx}(z,z^{})`$ $`=`$ $`{\displaystyle \frac{ik}{\kappa ^2}}{\displaystyle \frac{}{z}}\mathrm{\Gamma }_{xx}(z,z^{}),`$ (19)
$`\mathrm{\Gamma }_{zz}(z,z^{})`$ $`=`$ $`{\displaystyle \frac{k^2}{(\kappa ^2)}}{\displaystyle \frac{}{z}}{\displaystyle \frac{}{z^{}}}\mathrm{\Gamma }_{xx}(z,z^{})+{\displaystyle \frac{\omega ^2}{\kappa ^2}}\delta (zz^{}).`$
By taking the limit $`zz^{}`$, and thus omitting the delta function, $`\mathrm{\Gamma }_{kk}`$ becomes

$$\mathrm{\Gamma }_{kk}=\mathrm{\Gamma }_{xx}+\mathrm{\Gamma }_{yy}+\frac{k^2}{(\kappa ^2)^2}_z_z^{}\mathrm{\Gamma }_{xx},$$
(20)
or by (2), using $`ϵ=1`$,

$`\mathrm{\Gamma }_{kk}`$ $`=`$ $`\omega ^2g^E+_z_z^{}g^B+{\displaystyle \frac{k^2}{(\kappa ^2)^2}}_z_z^{}(_z_z^{}g^B)`$ (21)
$`=`$ $`\omega ^2g^E+(k^2+_z_z^{})g^B.`$
Using (11) and (3.1), this becomes

$$\mathrm{\Gamma }_{kk}=\frac{\omega ^2}{\kappa }+\frac{1}{2\kappa }\left[\omega ^2r+\left(2k^2\omega ^2\right)r^{}\right]e^{2\kappa z}.$$
(22)
Equation (18) gives a formal expectation value only, because the integral is divergent. However, the divergence comes only from the $`\omega ^2/\kappa `$ term in $`\mathrm{\Gamma }_{kk}`$ and is independent of $`z`$. It is the usual empty space vacuum divergence. We will henceforth drop this term and denote the resulting finite expectation value by $`E^2`$. The renormalization results in a quantity which vanishes at large distances from the interface: $`E^20`$ as $`z\mathrm{}`$, which amounts to finding the difference in $`E^2`$ with the boundary and without it. Thus we find

$$E^2=\frac{1}{4\pi ^2}_0^{\mathrm{}}𝑑\zeta _0^{\mathrm{}}𝑑k\frac{k}{\kappa }\left[\zeta ^2r+\left(2k^2+\zeta ^2\right)r^{}\right]e^{2\kappa z}.$$
(23)
### 3.3 The Magnetic Field

Now we compute the expectation value of the magnetic field. Using Eq. (7), we find

$`iB_x(𝐫)B_x(𝐫^{})`$ $`=`$ $`{\displaystyle \frac{1}{\omega ^2}}\left(_z_z^{}\mathrm{\Gamma }_{yy}_y_z^{}\mathrm{\Gamma }_{zy}_z_y^{}\mathrm{\Gamma }_{yz}+_y_y^{}\mathrm{\Gamma }_{zz}\right),`$
$`iB_y(𝐫)B_y(𝐫^{})`$ $`=`$ $`{\displaystyle \frac{1}{\omega ^2}}\left(_z_z^{}\mathrm{\Gamma }_{xx}_x_z^{}\mathrm{\Gamma }_{zx}_z_x^{}\mathrm{\Gamma }_{xz}+_x_x^{}\mathrm{\Gamma }_{zz}\right),`$ (24)
$`iB_z(𝐫)B_z(𝐫^{})`$ $`=`$ $`{\displaystyle \frac{1}{\omega ^2}}\left(_y_y^{}\mathrm{\Gamma }_{xx}_y_x^{}\mathrm{\Gamma }_{xy}_x_y^{}\mathrm{\Gamma }_{yx}+_x_x^{}\mathrm{\Gamma }_{yy}\right).`$
From the definition of $`𝐤_{}`$, it follows that all derivatives in $`y`$ vanish, so we can write the sum of the above terms as

$`iB_i(𝐫)B_i(𝐫^{})`$ $`={\displaystyle \frac{1}{\omega ^2}}(_z_z^{}\mathrm{\Gamma }_{yy}+_z_z^{}\mathrm{\Gamma }_{xx}`$
$`_x_z^{}\mathrm{\Gamma }_{zx}_z_x^{}\mathrm{\Gamma }_{xz}+_x_x^{}\mathrm{\Gamma }_{zz}+_x_x^{}\mathrm{\Gamma }_{yy}).`$ (25)
Using (19), we have

$`_x_x^{}\mathrm{\Gamma }_{zz}(𝐫,𝐫^{},\omega )`$ $`=`$ $`{\displaystyle \frac{d𝐤_{}}{(2\pi )^2}\frac{(k^2)^2}{(\kappa ^2)^2}_z_z^{}\mathrm{\Gamma }_{xx}},`$
$`_z_x^{}\mathrm{\Gamma }_{xz}(𝐫,𝐫^{},\omega )`$ $`=`$ $`{\displaystyle \frac{d𝐤_{}}{(2\pi )^2}\frac{k^2}{\kappa ^2}_z_z^{}\mathrm{\Gamma }_{xx,}}`$ (26)
$`_x_z^{}\mathrm{\Gamma }_{zx}(𝐫,𝐫^{},\omega )`$ $`=`$ $`{\displaystyle \frac{d𝐤_{}}{(2\pi )^2}\frac{k^2}{\kappa ^2}_z_z^{}\mathrm{\Gamma }_{xx}}.`$
This leads to

$$iB(x)B(x^{})=\frac{d\omega }{2\pi }\frac{d𝐤_{}}{(2\pi )^2}\left[\frac{1}{\omega ^2}\left(k^2+_z_z^{}\right)\mathrm{\Gamma }_{yy}(z,z^{})+\frac{\omega ^2}{\kappa ^4}_z_z^{}\mathrm{\Gamma }_{xx}(z,z^{})\right].$$
(27)
Using (2), this becomes

$$iB(x)B(x^{})=\frac{d\omega }{2\pi }\frac{d𝐤_{}}{(2\pi )^2}\left[\left(k^2+_z_z^{}\right)g^E(z,z^{})+\omega ^2g^B(z,z^{})\right].$$
(28)
Following the same procedure as used above in calculating $`E^2`$, we find the finite mean squared magnetic field to be

$$B^2=\frac{1}{4\pi ^2}_0^{\mathrm{}}𝑑\zeta _0^{\mathrm{}}𝑑k\frac{k}{\kappa }\left[\left(2k^2+\zeta ^2\right)r\zeta ^2r^{}\right]e^{2\kappa z}.$$
(29)
Note from (23) and (29) that $`E^2B^2`$ under interchange of $`r`$ and $`r^{}`$. Now, the mean energy density can be calculated as

$$U=\frac{1}{2}\left(E^2+B^2\right).$$
(30)
Using (23) and (29), this becomes

$$U=\frac{1}{4\pi ^2}_0^{\mathrm{}}𝑑\zeta _0^{\mathrm{}}𝑑k\frac{k^3}{\kappa }\left(r+r^{}\right)e^{2\kappa z}.$$
(31)
We can write $`U`$ in a form more convenient for numerical calculation by introducing polar coordinates $`u`$ and $`\theta `$ ($`\zeta =u\mathrm{cos}\theta ,`$ $`k=u\mathrm{sin}\theta `$):

$$U=\frac{1}{4\pi ^2}_0^{\mathrm{}}𝑑uu^3_0^{\frac{\pi }{2}}𝑑\theta (\mathrm{sin}\theta )^3\left(r+r^{}\right)e^{2uz}.$$
(32)
We use the Drude model for the dielectric function

$$ϵ_d(\omega )=1\frac{\omega _p^2}{\omega ^2},$$
(33)
where $`\omega _p`$ is the plasma frequency. From (16a), (16b), and (33), we find

$`r`$ $`=`$ $`{\displaystyle \frac{u\sqrt{u^2+\omega _p^2}}{u+\sqrt{u^2+\omega _p^2}}},`$ (34a)
$`r^{}`$ $`=`$ $`{\displaystyle \frac{u^2(\mathrm{cos}\theta )^2+\omega _p^2u(\mathrm{cos}\theta )^2\sqrt{u^2+\omega _p^2}}{u^2(\mathrm{cos}\theta )^2+\omega _p^2+u(\mathrm{cos}\theta )^2\sqrt{u^2+\omega _p^2}}}.`$ (34b)
By a substitution ($`\mathrm{cos}(\theta )t`$), U becomes

$$U=\frac{1}{4\pi ^2}_0^{\mathrm{}}𝑑uu^3\left[_0^1𝑑t(1t^2)\left(r+r^{}\right)\right]e^{2uz}.$$
(35)
By the same coordinate transform, (23) and (29) become

$`E^2`$ $`=`$ $`{\displaystyle \frac{1}{4\pi ^2}}{\displaystyle _0^{\mathrm{}}}𝑑uu^3\left\{{\displaystyle _0^1}𝑑t\left[t^2r+(2t^2)r^{}\right]\right\}e^{2uz},`$ (36a)
$`B^2`$ $`=`$ $`{\displaystyle \frac{1}{4\pi ^2}}{\displaystyle _0^{\mathrm{}}}𝑑uu^3\left\{{\displaystyle _0^1}𝑑t\left[(2t^2)rt^2r^{}\right]\right\}e^{2uz}.`$ (36b)
The plot for $`E^2`$ and $`B^2,`$ as well as $`U`$ is shown in Figure 1. As we can see from the figure, the energy density is positive. Now we consider some limiting cases.

### 3.4 The Fields near the Interface

To see how $`U`$ behaves for small $`z`$ (large $`u`$), we first perform the $`t`$ integration in (35), which can be done analytically, and then Taylor expand the resulting expression in the brackets in powers of $`u^1`$. That is, we are expanding all of the integrand except for the exponential factor. To the leading order we find:

$$U\frac{\sqrt{2}\omega _p}{64\pi }\frac{1}{z^3}.$$
(37)
The asymptotic behavior of the mean squared fields (36a), (36b) in this limit is

$`E^2`$ $``$ $`{\displaystyle \frac{\sqrt{2}\omega _p}{32\pi }}{\displaystyle \frac{1}{z^3}},`$ (38a)
$`B^2`$ $``$ $`{\displaystyle \frac{5\omega _p^2}{96\pi }}{\displaystyle \frac{1}{z^2}}.`$ (38b)
We see that $`E^2`$ dominates over $`B^2`$, so that $`U1/2E^2;`$ this is due fact that the leading order in the expression in braces in (36a) is $`u^1`$ as compared to $`u^2`$ in (36b). If we compare these expressions to ones that would result if dispersion were not included in the calculation, it can be seen from (3.3) that in this case $`E^2z^4`$, and same for $`B^2`$ (see also ). As seen from (3.4), the inclusion of dispersion in the calculation reduces the power in $`z`$ up to two orders, but it does not remove the singularity of the results at $`z=0`$, as might be naively expected.

After more careful consideration, it is not surprising that dispersion alone is insufficient to render the results finite at the interface. The integrals for $`E^2`$ and $`B^2`$ at $`z=0`$ diverge quartically in a frequency cutoff. However, in general dielectric functions go as

$$ϵ(\omega )1+O(\omega ^2)$$
(38c)
as $`\omega \mathrm{}`$. Thus the reflection coefficients will go to zero no faster than $`\omega ^2`$, leaving the integrals quadratically divergent. This argument explains why $`B^2z^2`$ for small $`z`$, but understanding the behavior of $`E^2`$ requires examining the dependence of the reflection coefficients $`r`$ and $`r^{}`$ upon the transverse momentum $`k`$. In fact the contribution of the coefficient $`r`$, which describes modes with the polarization vector perpendicular to the plane of incidence, does go as $`z^2`$. This coefficient depends only upon frequency, and falls as $`\omega ^2`$ for large $`\omega `$, as can be seen from (34a). The coefficient $`r^{}`$ describes modes with the polarization vector parallel to the plane of incidence, and goes to one as $`\theta \pi /2`$ (corresponding to grazing incidence) for all frequencies. It is this behavior which leads to the $`z^3`$ singularity in $`E^2`$ and hence in $`U`$. (The role of cutoffs for the quantized electromagnetic field in dielectrics has been discussed in more detail by Candelas . Barton  has recently emphasized the fact that dispersion alone will not remove all divergences.)

The divergence of $`U`$ is not considered to be physical, but as resulting from the idealization of the wall as a perfectly smooth surface. One way of removing this singularity is to allow the position of the boundary to fluctuate . It seems plausible that such effects as surface roughness, or the atomic nature of matter on small scales can also introduce a physical cutoff that makes the mean squared fields and the energy density finite everywhere.

### 3.5 Case of a Perfect Conductor

Now we consider the limit $`ϵ\mathrm{}`$. In this limit, as seen from (11), (16a) and (16b), $`r1`$, and $`r^{}1`$ . Equation (31) implies that $`U`$ becomes zero, as expected, and (36a) and (36b) give:

$`E^2`$ $``$ $`{\displaystyle \frac{3}{16\pi ^2}}{\displaystyle \frac{1}{z^4}},`$ (39)
$`B^2`$ $``$ $`{\displaystyle \frac{3}{16\pi ^2}}{\displaystyle \frac{1}{z^4}}.`$ (40)
These well-known results are consistent with the asymptotic Casimir-Polder potential :

$$V_{CP}\frac{3}{32\pi ^2}\frac{\alpha _0}{z^4}=\frac{1}{2}\alpha _0E^2,$$
(41)
where $`\alpha _0`$ is the static polarizability of an atom near the interface.

## 4 The energy density between dielectrics

In this section we calculate the energy density in a vacuum region of width $`a`$ between two dielectric half-spaces. We define the dielectric constant as:

$`0<z<a:ϵ`$ $`=`$ $`1,`$
$`z<0\mathrm{and}z>a:ϵ`$ $`=`$ $`ϵ_d.`$ (42)
In the vacuum region, $`g^E`$ occurring in (12a) has the form

$`g^E`$ $`={\displaystyle \frac{1}{2\kappa _0}}\{{\displaystyle \frac{\mathrm{exp}(\kappa _0zz^{})+r^2\mathrm{exp}(2\kappa _0a)\mathrm{exp}(\kappa _0zz^{})}{r^2\mathrm{exp}(2\kappa _0a)1}}+`$
$`+{\displaystyle \frac{r\left[\mathrm{exp}(\kappa _0(z+z^{}))+\mathrm{exp}(2\kappa _0a)\mathrm{exp}(\kappa (z+z^{}))\right]}{r^2\mathrm{exp}(2\kappa _0a)1}}\},`$ (43)
and $`g^B`$ has the same form as $`g^{E,}`$ with $`r`$ and $`r^{}`$ interchanged, where $`r`$ and $`r^{}`$ are defined in (16a) and (16b). We can now calculate the electric and magnetic fields in this region.

### 4.1 The Electric Field

The first and the second term on the right hand side of (21) for the present case, in the limit $`zz^{}`$, are

$`\omega ^2g^E`$ $`=`$ $`{\displaystyle \frac{\omega ^2}{2\kappa }}{\displaystyle \frac{1+r^2e^{2\kappa a}}{1r^2e^{2\kappa a}}}+{\displaystyle \frac{\omega ^2}{2\kappa }}{\displaystyle \frac{r(e^{2\kappa z}+e^{2\kappa a}e^{2\kappa z})}{1r^2e^{2\kappa a}}}`$ (44)
$`(k^2+_z_z^{})g^B`$ $`=`$ $`{\displaystyle \frac{\omega ^2}{2\kappa }}{\displaystyle \frac{1+r^{}{}_{}{}^{2}e_{}^{2\kappa a}}{1r^{}{}_{}{}^{2}e_{}^{2\kappa a}}}+{\displaystyle \frac{k^2+\kappa ^2}{2\kappa }}{\displaystyle \frac{r^{}(e^{2\kappa z}+e^{2\kappa a}e^{2\kappa z})}{1r^{}{}_{}{}^{2}e_{}^{2\kappa a}}},`$ (45)
so that $`\mathrm{\Gamma }_{kk}`$ becomes

$`\mathrm{\Gamma }_{kk}`$ $`={\displaystyle \frac{\omega ^2}{\kappa }}+{\displaystyle \frac{\omega ^2}{\kappa }}\left({\displaystyle \frac{r^2}{e^{2\kappa a}r^2}}+{\displaystyle \frac{r^{}^2}{e^{2\kappa a}r^{}^2}}\right)+`$
$`+{\displaystyle \frac{1}{\kappa }}\left[\omega ^2{\displaystyle \frac{r}{1r^2e^{2\kappa a}}}+(2k^2\omega ^2){\displaystyle \frac{r^{}}{1r^{}{}_{}{}^{2}e_{}^{2\kappa a}}}\right]e^{\kappa a}\mathrm{cosh}\left[\kappa (2za)\right].`$ (46)
We again drop the first term on the right hand side of the above expression. After introducing polar coordinates and using (11) and (46), (18) we find

$`E^2`$ $`={\displaystyle \frac{1}{2\pi ^2}}{\displaystyle _0^{\mathrm{}}}duu^3{\displaystyle _0^{\frac{\pi }{2}}}d\theta \mathrm{sin}(\theta )\{\mathrm{cos}^2(\theta )({\displaystyle \frac{r^2}{r^2e^{2ua}}}+{\displaystyle \frac{r^{}^2}{r^{}{}_{}{}^{2}e^{2ua}}})+`$
$`+[\mathrm{cos}^2(\theta ){\displaystyle \frac{r}{1r^2e^{2ua}}}+(1+\mathrm{sin}^2(\theta )){\displaystyle \frac{r^{}}{1r^{}{}_{}{}^{2}e_{}^{2ua}}}]e^{ua}\mathrm{cosh}\left[u(2za)\right]\},`$ (47)
or with ($`\mathrm{cos}(\theta )t`$)

$`E^2`$ $`={\displaystyle \frac{1}{2\pi ^2}}{\displaystyle _0^{\mathrm{}}}duu^3{\displaystyle _0^1}dt\{t^2({\displaystyle \frac{r^2}{r^2e^{2ua}}}+{\displaystyle \frac{r^{}^2}{r^{}{}_{}{}^{2}e^{2ua}}})+`$
$`+[t^2{\displaystyle \frac{r}{1r^2e^{2ua}}}+(2t^2){\displaystyle \frac{r^{}}{1r^{}{}_{}{}^{2}e_{}^{2ua}}}]e^{ua}\mathrm{cosh}\left[u(2za)\right]\}.`$ (48)
### 4.2 The Magnetic Field

In the same way, (28) leads to

$`B^2`$ $`={\displaystyle \frac{1}{2\pi ^2}}{\displaystyle _0^{\mathrm{}}}duu^3{\displaystyle _0^1}dt\{t^2({\displaystyle \frac{r^2}{r^2e^{2ua}}}+{\displaystyle \frac{r^{}^2}{r^{}{}_{}{}^{2}e^{2ua}}})+`$
$`+[(2t^2){\displaystyle \frac{r}{1r^2e^{2ua}}}t^2{\displaystyle \frac{r^{}}{1r^{}{}_{}{}^{2}e_{}^{2ua}}}]e^{ua}\mathrm{cosh}\left[u(2za)\right]\}.`$ (49)
This expression differs from the one for $`E^2`$ only in the $`z`$\- dependent term with $`rr^{}`$. A plot of $`E^2`$ , $`B^2`$ and $`U`$ is shown in Figure 2 for $`\omega _pa=200.`$ It can be seen in this figure that significant, but not complete cancellation occurs between $`E^2`$ and $`B^2.`$

### 4.3 The Energy Density

Now using (30), (48), and (49), the energy density in the vacuum region can be calculated as

$`U`$ $`={\displaystyle \frac{1}{2\pi ^2}}{\displaystyle _0^{\mathrm{}}}duu^3{\displaystyle _0^1}dt\{t^2({\displaystyle \frac{r^2}{r^2e^{2ua}}}+{\displaystyle \frac{r^{}^2}{r^{}{}_{}{}^{2}e^{2ua}}})+`$
$`+(1t^2)[{\displaystyle \frac{r}{1r^2e^{2ua}}}+{\displaystyle \frac{r^{}}{1r^{}{}_{}{}^{2}e_{}^{2ua}}}]e^{ua}\mathrm{cosh}\left[u(2za)\right]\}`$ (50)
By analyzing this expression we can make some conclusions about the sign of $`U`$. First we note that it is position dependent, and we also note that the first term is always negative and the second term is always positive. The overall sign of $`U`$ depends on the choice of $`a`$ and $`\omega _p`$. As $`\omega _pa`$ grows, $`U`$ at the midpoint decreases, becoming negative for $`\omega _pa100`$, as seen in Figure 3 and Figure 4.

In Figure 4 we see how the energy density at the center of the vacuum region changes as the product $`\omega _pa`$ increases. It can be seen both in Figure 3 and Figure 4 that U approaches the value given in (52) as $`\omega _pa`$ becomes large. The separation at which $`U`$ becomes negative at the center is

$$a>a_c=\frac{99}{\omega _p}=1.3\mu \mathrm{m}\left(\frac{14.8e\mathrm{V}}{\omega _p}\right),$$
(51)
where $`14.8e\mathrm{V}`$ is the plasma frequency of aluminum.

### 4.4 A Perfect Conductor Case

In the limit $`\omega _p\mathrm{}`$, $`r1`$ and $`r^{}1.`$ Then only the first ($`z`$-independent) term survives in (50), and we get the familiar result :

$$U=\frac{\pi ^2}{720a^4}.$$
(52)
It is also interesting to examine the expressions for $`E^2`$ and $`B^2`$ in this limit. After performing the $`t`$ integration, (48) can be written as

$$E^2=\frac{1}{2\pi ^2}_0^{\mathrm{}}𝑑uu^3\left[\frac{2}{3}\frac{1}{1e^{2ua}}+\left(\frac{e^{2u(az)}}{1e^{2ua}}+\frac{e^{2uz}}{1e^{2ua}}\right)\right],$$
(53)
and after performing the $`u`$ integration as

$$E^2=\frac{\pi ^2}{720a^4}+\frac{1}{32\pi ^2a^4}\left[\psi ^{(3)}\left(\frac{z}{a}\right)+\psi ^{(3)}\left(1\frac{z}{a}\right)\right].$$
(54)
Here $`\psi ^{(3)}=\frac{d^4}{dz^4}\mathrm{ln}\mathrm{\Gamma }(z)`$ is the polygamma function of order three. It satisfies the reflection formula

$$\psi ^{(3)}\left(\frac{z}{a}\right)+\psi ^{(3)}\left(1\frac{z}{a}\right)=\pi \frac{d^3}{d(\frac{z}{a})^3}\mathrm{cot}\left(\pi \frac{z}{a}\right).$$
(55)
This yields

$$E^2=\frac{\pi ^2}{720a^4}+\frac{\pi ^2}{16a^4}\frac{1+2\mathrm{cos}^2(\pi \frac{z}{a})}{\mathrm{sin}^4(\pi \frac{z}{a})}.$$
(56)
In the same way, one finds

$$B^2=\frac{\pi ^2}{720a^4}\frac{\pi ^2}{16a^4}\frac{1+2\mathrm{cos}^2(\pi \frac{z}{a})}{\mathrm{sin}^4(\pi \frac{z}{a})}.$$
(57)
These results are in agreement with those given by previous authors .

### 4.5 Energy density near the Boundary

To see how $`U`$ grows near the interface we note that for small $`z`$ (large $`u`$), (50) becomes

$$U\frac{1}{4\pi ^2}_0^{\mathrm{}}𝑑uu^3_0^1𝑑t\left(1t^2\right)\left(r+r^{}\right)e^{2uz},$$
(58)
so it reduces to (35), the solution for one interface case. In the $`z0`$ limit, this expression reduces to (37). By the same reasoning, the expressions for $`E^2`$ and $`B^2`$, (48) and (49) reduce to (36a) and (36b), respectively, in the small $`z`$ limit.

## 5 Conclusion

In this paper we have examined the mean squared fields and the energy density in the region between a pair of half-spaces filled with dispersive media. We found that these quantities diverge at the boundaries of the media, despite the inclusion of dispersion in the calculation. This divergence indicates a breakdown of a continuum description in which the dielectric function changes discontinuously at the boundary. It also shows that there is a positive self-energy density in the region outside of a single plate. Nonetheless, we have also found that it is possible for the net energy density in the region between the plates to be negative, depending upon the plate separation and the plasma frequency of the material involved. The existence of an attractive Casimir force is not an indicator of whether the energy density at the center of the plates is actually negative or not.

We have found that the energy density at the center becomes negative when $`\omega _pa>100`$. Thus for fixed plasma frequency $`\omega _p`$, the energy density always becomes negative for sufficiently large separation $`a`$. Of course, in this limit the magnitude of the energy density is also becoming small. Similarly, for fixed $`a`$, the energy density becomes negative for sufficiently large $`\omega _p`$. In the limit that $`\omega _p\mathrm{}`$ our results approach the constant negative energy density of the perfectly conducting plates. It should not come as a surprise that there is a regime of negative energy density. The calculation assuming perfect conductivity does have a region of validity so long as $`\omega _p`$ is large and one is not too close to one plate. Qualitatively similar behavior has recently been found for the vacuum energy density near a domain wall .

In this paper, we assumed a particular form for the dielectric function, (33), given by the collisionless Drude model. This is a good model for many metals especially alkali metals, and is the generic form for the dielectric function of all materials at high frequencies. Thus taking a different form for $`ϵ_d`$ would change the details of our results, especially far from an interface, but should lead to the same limiting forms near an interface. We have also assumed zero temperature throughout this paper. For systems at room temperature, this should be a good approximation when the separations are of the order of a few micrometers. More generally, one can ignore thermal effects at distances small compared to $`1/(kT)`$.

Contrary to the view expressed by Lamoreaux , the appearance of negative energy density in a quantum field theory is very natural. One can easily find quantum states of the free quantized electromagnetic field in empty space which have local negative energy densities. A squeezed vacuum state is an example . The energy density of a quantized field has to be defined as a difference between that in empty Minkowski spacetime, and that in a given state and is no longer positive definite, as it was for a classical field. Apart from coupling to gravity, which produces extremely small effects, no clear way has been found to directly observe the local energy density. In certain limits, the negative energy density in a squeezed vacuum state has been shown theoretically  to produce an effect on the magnetic moment of a spin system. Whether this effect could ever be observed, and whether negative Casimir energy density can produce similar effects is unknown.

Acknowledgement: We would like to thank J.-T. Hsiang for valuable discussions. This work was supported in part by the National Science Foundation under Grant PHY-9800965.


format: text

text:
# Responses of quark condensates to the chemical potential

## Abstract

The responses of quark condensates to the chemical potential, as a function of temperature $`T`$ and chemical potential $`\mu `$, are calculated within the Nambu–Jona-Lasinio (NJL) model. We compare our results with those from the recent lattice QCD simulations \[QCD-TARO Collaboration, Nucl. Phys. B (Proc. Suppl.) 106, 462 (2002)\]. The NJL model and lattice calculations show qualitatively similar behavior, and they will be complimentary ways to study hadrons at finite density. The behavior above $`T_c`$ requires more elaborated analyses.

The variation of quark condensates in medium plays a key role to understand the behavior of hadron masses and chiral symmetry restoration br91 ; hl92hkl93 . Recently, we calculated for the first time the second order response of the quark condensate to the chemical potential $`\frac{d^2\overline{q}q}{d\mu ^2}`$ at $`\mu `$ = 0 using lattice QCD taro-lattice  following the method in taro-prd . It was found that the response is negative both below and above $`T_c`$. One of interesting results is that the response to the isoscalar chemical potential ($`\mu _S`$ = $`\mu _u`$ = $`\mu _d`$) is almost the same as that to the isovector chemical potential ($`\mu _V`$ = $`\mu _u`$ = – $`\mu _d`$), where $`\mu _u`$ ($`\mu _d`$) is the $`u`$ ($`d`$) quark chemical potential. It would be interesting if we can check this result within effective models of QCD. In this work we present an Nambu–Jona-Lasinio (NJL) model njl ; vw91k92  calculation of $`\frac{d\overline{q}q}{d\mu }`$ and $`\frac{d^2\overline{q}q}{d\mu ^2}`$ at $`\mu `$ = 0 and compare our results with those from the lattice QCD simulations taro-lattice .

An SU(2) NJL model Lagrangian will be enough for that purpose. However, we are also interested in the responses of the s-quark condensate and a comparison with those of the u,d-quark condensates will be useful for future studies in the lattice calculations. Thus we used an SU(3) NJL model in our calculations and found that the effects of the flavor mixing are negligible. We present only the results for the u,d-quark condensates in this paper, and we will make the comparison in a forthcoming paper.

First, let us consider an SU(3) NJL model Lagrangian hk94 :

$`L`$ $`=`$ $`\overline{q}(i\gamma dm)q+{\displaystyle \frac{1}{2}}g_S{\displaystyle \underset{a=0}{\overset{8}{}}}\left[(\overline{q}\lambda _aq)^2+(\overline{q}i\lambda _a\gamma _5q)^2\right]`$ (1)
$`+`$ $`g_D[\mathrm{det}\overline{q}_i(1\gamma _5)q_j+h.c.],`$
where $`\lambda _a`$ are the Gell-Mann matrices and $`m`$ is a mass matrix for current quarks, $`m`$=diag($`m_u`$, $`m_d`$, $`m_s`$). We take the following parameters in hk94 :

$`\mathrm{\Lambda }=631.4\mathrm{MeV},g_S\mathrm{\Lambda }^2=3.67,g_D\mathrm{\Lambda }^5=9.29`$
$`m_u=m_d=5.5\mathrm{MeV},m_s=135.7\mathrm{MeV},`$ (2)
where $`\mathrm{\Lambda }`$ is the momentum cut-off. The third term in Eq.(1) is a reflection of the axial anomaly, and causes a mixing in flavors. For example, the constituent quark masses are given as follows.

$`M_u`$ $`=`$ $`m_u2g_S\overline{u}u2g_D\overline{d}d\overline{s}s,`$
$`M_d`$ $`=`$ $`m_d2g_S\overline{d}d2g_D\overline{u}u\overline{s}s,`$
$`M_s`$ $`=`$ $`m_s2g_S\overline{s}s2g_D\overline{u}u\overline{d}d,`$ (3)
where $``$ means the statistical average. In this work we concentrate mostly on the Case II in hk94 , where only $`g_D`$ has a temperature dependence

$`g_D(T)=g_D(T=0)\mathrm{exp}[(T/T_0)^2],`$ (4)
while other coupling constants and the cut-off are independent of $`T`$ and chemical potential (or density). Here, we set $`T_0`$ = 0.1 GeV taking into account the restoration of $`U_A(1)`$ symmetry as in hk94 .

In the mean-field approximation the above Lagrangian leads to the following gap equation hk94 :

$`Q_i2N_c{\displaystyle \underset{p}{}}\left({\displaystyle \frac{M_i}{E_{ip}}}f(E_{ip})\right)=\overline{q}_iq_i,`$ (5)
where the index $`i`$ denotes the $`u`$, $`d`$, and $`s`$ quarks. $`N_c`$ is the number of colors and $`M_i`$ is the constituent quark mass, and $`E_{ip}=\sqrt{M_i^2+p^2}`$. $`f(E_{ip})=1n_{ip}\overline{n}_{ip}`$, where $`n_{ip}`$ and $`\overline{n}_{ip}`$ are the distribution functions of the $`i`$th quark and antiquark, respectively, i.e.,

$`n_{ip}={\displaystyle \frac{1}{1+\mathrm{exp}((E_{ip}\mu _i)/T)}},`$
$`\overline{n}_{ip}={\displaystyle \frac{1}{1+\mathrm{exp}((E_{ip}+\mu _i)/T)}}.`$ (6)
The l.h.s. of Eq.(5) is a function of $`\overline{u}u`$, $`\overline{d}d`$, $`\overline{s}s`$, $`\mu _i`$, and $`T`$. Then, we obtain the first order response of quark condensates $`\frac{d\overline{u}u}{d\mu }`$, $`\frac{d\overline{d}d}{d\mu }`$, and $`\frac{d\overline{s}s}{d\mu }`$ by differentiating both sides with respect to $`\mu `$ at a fixed $`T`$, i.e., we solve the following equations

$`{\displaystyle \frac{dQ_u(\overline{u}u,\overline{d}d,\overline{s}s,\mu _u)}{d\mu }}`$ $`=`$ $`{\displaystyle \frac{d\overline{u}u}{d\mu }},`$
$`{\displaystyle \frac{dQ_d(\overline{u}u,\overline{d}d,\overline{s}s,\mu _d)}{d\mu }}`$ $`=`$ $`{\displaystyle \frac{d\overline{d}d}{d\mu }},`$
$`{\displaystyle \frac{dQ_s(\overline{u}u,\overline{d}d,\overline{s}s,\mu _s)}{d\mu }}`$ $`=`$ $`{\displaystyle \frac{d\overline{s}s}{d\mu }}.`$ (7)
Here, we consider two types of $`\mu `$ following the notation in the lattice QCD simulations taro-prd ; taro-lattice . One is the isoscalar chemical potential $`\mu _S`$ = $`\mu _u`$ \+ $`\mu _d`$, and the other is the isovector one $`\mu _V`$ = $`\mu _u`$$`\mu _d`$.

Let us comment on the definition of chemical potential used in the lattice calculations and ours. There will be no difference between the two definitions if we choose the same u- and d-quark mass. In the NJL model calculations we chose rather general definitions for the isoscalar and isovector chemical potentials, i.e., $`\mu _S`$ = $`\mu _u+\mu _d`$ and $`\mu _V`$ = $`\mu _u\mu _d`$ to show the quark mass dependence by taking different u- and d-quark masses.

Using Eq. (7) one can show that both $`\frac{d\overline{q}q}{d\mu _S}`$ and $`\frac{d\overline{q}q}{d\mu _V}`$ are zero at $`\mu `$=0. This is the same as the lattice QCD calculation taro-lattice . At finite chemical potential, the absolute value of the quark condensate decreases with increasing chemical potential and the variation is proportional to the chemical potential in the present NJL model mc .

Next, consider the second order response $`\frac{d^2\overline{q}q}{d\mu ^2}`$. This is obtained by differentiating each equation in Eq. (7) with respect to $`\mu _S`$ (or $`\mu _V`$) again. In Fig 1 we show $`\frac{d^2|\overline{q}q|}{d\mu ^2}`$ for the $`u`$ and $`d`$ quark condensates at $`\mu `$ = 0. Within the present NJL model the second order response to the isoscalar chemical potential $`\mu _S`$ is the same as that to the isovector chemical potential $`\mu _V`$. This is consistent with the result in taro-lattice . In fact, there are a few different terms between the isoscalar and the isovector cases in the lattice QCD simulations. However, it is found that the contribution of those terms is negligible as shown in the figure 3 of Ref.taro-lattice .

In the lattice calculations the same u- and d-quark mass was taken taro-lattice . In this paper, however, we would like to show the quark mass dependence of the response in the figure. Our main interest is to compare the responses below and above the Mott temperature depending on the quark mass. The present results may be useful for future lattice QCD simulations which use different u- and d-quark masses.

We take $`m_u`$ = 4 MeV and $`m_d`$ = 7 MeV as well as $`m_u`$ = $`m_d`$ = 5.5 MeV to show the quark mass dependence of $`\frac{d^2|\overline{q}q|}{d\mu ^2}`$. Although the cut-off and the coupling constants should be modified, we use the same coupling constants and the cut-off for both cases and study the behavior of $`\frac{d^2\overline{q}q}{d\mu ^2}`$. Below the pion Mott temperature $`T_{m_\pi }`$ ($``$ 148 MeV in this work), the response of the $`u`$ ($`m_u`$ = 4 MeV) quark condensate is larger than that of the $`d`$ ($`m_d`$ = 7 MeV) quark condensate, while above the Mott temperature the behavior of the response is opposite (Here, $`T_{m_\pi }`$ is determined as a temperature at which the sum of the $`u`$ and $`d`$ constituent quark masses equals to the pion mass, i.e., $`M_u+M_d=m_\pi `$). This behavior results from the quark mass dependence of the first order response at finite chemical potential shown in Fig. 2, where we take $`\mu _u`$ = $`\mu _d`$ = 0.02 and 0.04 GeV, respectively. That behavior in Fig. 2 is not altered even in the Case I hk94 , where all the couplings and the cut-off are independent of temperature and chemical potential. For example, in the Case I and at $`\mu `$ = 0.02 GeV, the pion Mott temperature is about 203 MeV and the extremum is located at around 186 MeV.

In Fig. 1, we also present the relative variation $`\frac{d^2\overline{q}q}{d\mu ^2}`$/$`\overline{q}q`$ at $`\mu `$=0 for comparison. Now, $`\frac{d^2\overline{q}q}{d\mu ^2}`$/$`\overline{q}q`$ has the extremum at the Mott temperature and the (absolute value of) variation decreases with increasing the quark mass. We find similar behavior in the lattice data, although $`\overline{q}q`$ is very small over $`T_c`$ in this case.

As we know, we can not detect a quark condensate itself in experiments. However, variations of quark condensates, such as shown in Fig. 2, are necessary to predict responses of hadron masses to the chemical potential, e.g., $`\frac{dm}{d\mu }`$ in mc . A more detailed analysis on the extremum point and its effects on hadron masses is in progress.

In summary, we have calculated the first and second order responses of quark condensates to the chemical potential within the NJL model and found that they are consistent with those from the recent lattice QCD simulations.

###### Acknowledgements.

We would like to thank Prof. Su H. Lee for valuable comments. Two of the authors (A.N. and T.T.) were supported by the Grant-in-Aide for Scientific Research by the Ministry of Education and Culture, Japan (11440080, 12554008, 13135216, and 13740164), and S.C. was supported by the Korean Ministry of Education under the BK21 program and in part by the Japan Society for the Promotion of Science (JSPS). This paper is based on a contribution to the Prof. Osamu Miyamura memorial symposium held at Institute for Nonlinear Sciences and Applied Mathematics (INSAM), Hiroshima University, Nov. 16–17, 2001.


text:
# Vertical structure models of T Tauri and Herbig Ae/Be disks

## 1 Introduction

The dust continuum emission often observed from T Tauri stars and Herbig Ae/Be stars is widely believed to originate from a dust/gas disk surrounding these stars (see e.g. Beckwith & Sargent,1996). These disks presumably have mass between $`10^4M_{}`$ and few$`\times 10^1M_{}`$, and are thought to be the disks of dust and gas from which the stars were formed. At early stages of the disk’s evolution, the disks may still accrete, and lose mass to the central star (Calvet, Hartmann, & Strom Calvet et al. (2000)). The energy released during this process is then the main source of power for the dust continuum emission from the disk. At later stages of the disk’s evolution the mass accretion rate drops, and the outer regions of the disk start to become dominated by stellar irradiation instead of viscous energy release. As the accretion rate drops even further, more and more of the disk becomes “passive”, and irradiation eventually will be the only source of heating for the disk. As was discussed by Kenyon & Hartmann (1987), such a passive disk can adopt a flaring shape, i.e. a shape in which the ratio of the vertical surface height $`H_s(R)`$ over $`R`$ increases with $`R`$, the distance from the central star. In this way the disk’s surface always “sees” the stellar surface, and therefore intercepts some fraction of the stellar radiation at all radii.

The detailed vertical structure of such a disk can be computed by solving the equations of vertical pressure balance coupled to the equations of radiative transfer. Such computations have been done by e.g. D’Alessio et al. (1998,1999), and Bell et al. (1997,1999). Their models include detailed physics, including viscous dissipation, heating by cosmic rays, and a study of the effects of self-gravity. A basic weakness of these models is their simplified treatment of radiative transfer, which uses the frequency-integrated moment equations in the Eddington approximation with Planck- and Rosseland-mean opacities. The reason for this simplification is that solving the full angle- and frequency-dependent radiative transfer equations is a challenging technical problem. It is known that straightforward methods for radiative transfer (e.g. Lambda Iteration and Monte Carlo methods) converge extremely slowly at high optical depth, and one is never completely sure that a true solution has been reached.

In the field of stellar atmospheres these problems are well known, and many different sophisticated radiative transfer algorithms have been developed over the years. Such algorithms include Accelerated Lambda Iteration, Complete Linearization methods, and Variable Eddington Factor methods (for a review of radiative transfer techniques, see e.g. Kudritzki & Hummer 1990; Hubeny 1992; Mihalas 1978).

In particular the latter method has already been used successfully in application to accretion disk models. Hubeny (1990) showed that when the Eddington factors and mean opacities are known, the disk temperature at every optical depth can be expressed by a simple analytic formula. He presented an analytical solution for the case in which the Eddington approximation is adopted and the Planck- and Rosseland mean opacities are used. Malbet & Bertout (1991) and Malbet, Lachaume & Monin (2001) also presented disk equations based on the variable Eddington factor method, and present solutions to the full angle-dependent transfer problem. Their solutions, however, are based on the assumption of a grey opacity, and therefore fall short of our aims. A full frequency-angle dependent vertical structure model has been presented by Sincell & Krolik (1997) for X-ray irradiated accretion disks around active galactic nuclei. Technically, that work comes close to the kind of calculations we wish to make in this paper, but in addition to the different physics involved, their calculations only solve the upper irradiated skin of the disk.

It is the goal of this paper to present complete vertical structure models for irradiated passive flaring disks based on full angle-frequency-dependent radiative transfer and vertical hydrostatic equilibrium. We use the method of Variable Eddington Factors as our main radiative transfer algorithm, and we present a variant of this algorithm that works fast and is stable under all circumstances. An Accelerated Lambda Iteration (ALI) algorithm is used to check the results and make sure that a true solution to the transfer equation is obtained.

## 2 The disk model

The structure equations for a passive irradiated disk can be grouped in three sets. First we have the equations describing the transfer of primary (stellar) photons as they move from the star radially outwards, and eventually get absorbed by the dust grains in the upper layers of the disk. The energy absorbed in these surface layers will be re-emitted at infrared wavelengths. Half of this radiation will travel upwards and escape from the disk. The other half will move downwards into the disk, where it will once again be absorbed and re-emitted. Since the disk’s optical depth can be rather high, in particular at short wavelengths, this process can repeat itself many times. In this way the energy diffuses all the way down to the equatorial plane of the disk. This process of radiative diffusion can be well described by the equations for plane-parallel 1-D vertical radiative transfer. It is here that most treatments in the literature use the moment equations with the Eddington approximation and mean opacities (henceforth MEMO method). In this paper we treat this problem in a fully angle-frequency dependent way. Once this problem has been solved and the temperature stratification has been found, we can proceed to solve the third and final equation: the equation of vertical pressure balance.

These three coupled sets of equations are solved iteratively: we move from stage 1 (primary stellar radiation), to stage 2 (diffuse radiation field) to stage 3 (vertical pressure balance), and back to stage 1. This is repeated until the relative difference in the density between successive iterations drops below the convergence criterion, generally taken to be $`10^2`$.

### 2.1 Stage 1: Extinction of primary photons

The impinging of primary (stellar) radiation onto the surface of the disk is modeled using a “grazing angle” recipe, similar to that used by D’Alessio et al. (1998) and Chiang & Goldreich (1997, henceforth CG97). Deliberately we do not use full 2-D/3-D ray-tracing, since this procedure can lead to numerical instabilities (Dullemond in prep.). The “grazing angle” recipe models the irradiation using vertical plane-parallel radiative transfer, with the radiation entering the disk under an angle $`\beta (R)`$ with respect to the disk’s surface. At each radius $`R`$ and each vertical height $`z`$ we evaluate the local primary stellar flux:

$$F_\nu (R,z)=\frac{L_\nu }{4\pi R^2}\mathrm{exp}\left(\tau _\nu (R,z)/\beta (R)\right),$$
(1)
where $`L_\nu `$ is the stellar luminosity and the vertical optical depth $`\tau _\nu (R,z)`$ is defined as

$$\tau _\nu (R,z)=_z^{\mathrm{}}\rho (R,z)\kappa _\nu 𝑑z,$$
(2)
with $`\rho (R,z)`$ the dust density (g cm-3) and $`\kappa _\nu `$ the dust absorption opacity (cm2 g-1). We ignore scattering, and we ignore the higher order geometrical effects which start to play a role when $`z/R1`$. The grazing angle $`\beta (R)`$ is defined as $`\beta (R)=0.4R_{}/R+Rd(H_s/R)/dR`$, where $`H_s`$ is the surface height of the disk and $`R_{}`$ the stellar radius (see CG97). We define $`H_s`$ as the height above the midplane where 63% (i.e., $`1\mathrm{exp}(1)`$) of the integrated stellar radiation has been absorbed. Although this definition is somewhat arbitrary, it is not critical for the results presented in this paper.

It is convenient to express $`\beta (R)`$ as

$$\beta (R)=0.4\frac{R_{}}{R}+\xi (R)\frac{H_s}{R},$$
(3)
where $`\xi `$ is the flaring index defined as:

$$\xi \frac{d\mathrm{lg}(H_s/R)}{\mathrm{lg}(R)}.$$
(4)
This flaring index is a number of order $`\xi 2/7`$ (CG97). Its value is computed self-consistently during the iteration procedure. Usually the iteration is started with a guess for $`\xi `$ (for instance $`\xi =2/7`$), and updated after each few iteration steps. In order to avoid numerical instabilities, $`\xi `$ is always evaluated two radial gridpoints away from the point where it is used (see appendix of Chiang et al. 2001).

It is important to note that a proper self-consistent computation of the flaring index $`\xi `$ is crucial if one wants to achieve energy conservation. Assuming it to be some fixed value is guaranteed to result in disks that emit more (or less) radiation than they receive.

### 2.2 Stage 2: Vertical radiative transfer

Once the function $`F_\nu (R,z)`$ is known, one can compute the amount of absorbed primary radiation per volume element:

$$q(R,z)=_0^{\mathrm{}}\rho (R,z)\kappa _\nu F_\nu (R,z)𝑑\nu $$
(5)
This energy will then be re-emitted as infrared radiation, half of which will diffuse towards lower $`z`$ into the disk interior. The transfer of this re-emitted IR radiation through the disk can be approximated as a 1-D slab geometry transfer problem of high optical depth.

The intensity of this diffuse infrared radiation field is denoted by $`I_{\mu ,\nu }`$ and obeys the following radiative transfer equation:

$$\mu \frac{dI_{\mu ,\nu }}{dz}=\rho \kappa _\nu (B_\nu (T)I_{\mu ,\nu }),$$
(6)
where $`\kappa _\nu `$ is the opacity per unit of matter, $`\rho `$ is the material density and $`B_\nu (T)`$ is the Planck function. This differential equation has to be integrated along $`z`$ for all $`\nu `$ and $`\mu `$, in order to obtain the intensity function $`I_{\mu ,\nu }(z)`$. The dust temperature $`T`$ is determined by assuming thermodynamic equilibrium with the radiation field:

$$_0^{\mathrm{}}\rho \kappa _\nu B_\nu (T)𝑑\nu =_0^{\mathrm{}}\rho \kappa _\nu J_\nu 𝑑\nu +\frac{q}{4\pi },$$
(7)
where the mean intensity $`J_\nu `$ is defined as

$$J_\nu (z)=\frac{1}{2}_1^1I_{\mu ,\nu }(z)𝑑\mu .$$
(8)
Equations (6,7,8) form a set of coupled integro-differential equations in the $`\mu `$, $`\nu `$, $`z`$ space. A straightforward way to solve them would be to use the method of Lambda Iteration (LI): first evaluate $`I_{\mu ,\nu }(z)`$ (Eq. 6), then $`J_\nu (z)`$ (Eq. 8) and finally $`T(z)`$ (Eq. 7), and iterate this procedure until convergence is reached. However, it is well known that this method converges very slowly at high optical depths (see e.g. Rybicki & Hummer 1991).

Instead, we use the method of variable Eddington factors (VEF, see appendix A for a description of the method). This method estimates the shape of the radiation field by using a moment equation, but eventually reaches a solution that solves the full frequency- and angle-dependent transfer problem. It converges generally within $``$10 iterations, and is therefore much faster and more reliable than the ALI method.

### 2.3 Stage 3: vertical hydrostatic equilibrium

Once the temperature at all locations in the disk is known, one can integrate the pressure balance equation to find the density structure. The vertical pressure balance equation is:

$$\frac{dP(R,z)}{dz}=\rho (R,z)\frac{GM}{R^3}z.$$
(9)
Since $`P=k\rho T/\mu m_u`$ (with $`\mu =2.3`$ for a H2, He mixture), and the derivatives of $`T`$ are known, this equation can be cast into the form of an integral in $`\rho `$. This integration is started at $`z=0`$, taking $`\rho (R,0)`$ at first as an arbitrary value. One can then compute the resulting vertical surface density:

$$\mathrm{\Sigma }=2_0^{\mathrm{}}\rho (R,z)𝑑z,$$
(10)
and then renormalize the density profile to the actual surface density. This new density structure is then compared to the density structure of the previous iteration. If the convergence criterion is not met, the iteration proceeds by going back to Stage 1. This iterative procedure has proven to be stable and to converge quickly, typically after about five to eight iterations.

## 3 Resulting vertical structure

As our example case we take a star with $`T_{\mathrm{eff}}=3000\text{K}`$, $`R_{}=2.0R_{}`$ and $`M_{}=0.5M_{}`$. Our disk has an inner radius of $`R_{\mathrm{in}}=3R_{}`$, an outer radius of $`R_{\mathrm{out}}=300\text{AU}`$ and a gas+dust surface density $`\mathrm{\Sigma }=\mathrm{\Sigma }_0(R/\text{AU})^1`$ with $`\mathrm{\Sigma }_0`$ the surface density at $`1\text{AU}`$. We assume that the gas and the dust are always mixed in the mass ratio $`100:1`$, and we take for the dust opacity the opacity of astronomical silicate of Draine & Lee (1984) for grains of 0.1 $`\mu `$m size.

In Fig. 1 the vertical structure of the disk at a distance of $`1\text{AU}`$ from the star is shown. For this calculation we took $`\mathrm{\Sigma }_0=10^3\text{g}/\text{cm}^2`$, which corresponds to a visual optical depth through the disk $`\tau _V=2.3\times 10^4`$. To better illustrate the results, we compare then to those obtained for the same parameters using the more approximated MEMO method of solution. The left panel shows the vertical temperature profile. One sees that the full transfer model has a much smoother structure than the MEMO solution, and that the midplane temperature is significantly lower, at least for the high optical depth case shown here. The main effect of the lower midplane temperature is a slightly smaller pressure scale height at the equator ($`H_p/R=0.022`$ for the VEF models versus $`H_p/R=0.028`$ for the MEMO one, as expected from the factor 1.57 difference in the midplane temperature between the two models), so that the disk is a bit more compressed toward the equatorial plane.

The right panel of Fig. 1 shows the run of the density as function of $`z/R`$. One can see how the two solutions differ only by a small factor near the midplane, since the pressure scale heights are not very different. At higher $`z/R`$ the differences become much larger (an order of magnitude at $`z/R0.0750.1`$). This can be understood as a slight vertical shift in $`z`$ of the steep density profile.

Fig. 2 illustrates some technical aspects of the radiation transfer solutions, which clarify the results just shown. Firstly, one can see (right panel) how the Eddington factor $`f(z)`$ drops below $`1/3`$ in the upper layers of the disk. This is because in these optically thin layers, most of the intensity is more or less parallel to the disk instead of pointing out of the disk. This is an effect that is opposite to what happens for instance in stellar winds, where most of the radiation eventually gets beamed outwards, and the value of $`f`$ becomes close to $`1`$.

It is also clear from Fig. 2 (left panel) that as a result of the use of the Rosseland mean opacity in the MEMO equations, the mean intensity in the MEMO method remains virtually constant, whereas in reality the mean intensity drops strongly towards the equatorial plane. This results in a non-isothermal disk structure, as opposed to what would be expected on the basis of the usual diffusion approximation considerations. According to these considerations one would expect to have a perfectly isothermal disk interior since there are no sources of heat in the disk. Indeed, if a grey opacity is used in the full transfer model as a test, the disk interior becomes isothermal as expected. But as soon as non-grey dust opacities are included, the isothermality is broken. The reason is that some radiation is able to leak out of the disk at long wavelengths, where the opacity is much smaller than at the wavelength of the bulk of the radiation. The slight positive temperature gradient causes a downward diffusive flux at shorter wavelengths that exactly cancels this energy leakage, so that the total flux remains zero deep within the disk.

In Fig.3 the vertical temperature structure is shown for different values of the radius. Near the star, where the disk optical depth is very large, the equatorial plane temperature is usually below that of the MEMO model, whereas at large radii, where the disk becomes optically thin to its own diffuse radiation, the MEMO models predict lower temperature values than what we find. The latter effect is related to the fact that the MEMO method treats the diffuse radiation field in a frequency-integrated way, and has therefore no information on how the energy is distributed over frequency. In other words: the MEMO method uses mean opacities based on the local dust temperature, whereas the local radiation temperature may be different in the optically thin case.

In Fig.3 it can also be seen that at small radii the model-predicted temperature for large $`z/R`$ is exactly constant. In reality we expect a small decrease of $`T`$ for increasing $`z/R`$ since the distance to the star increases as $`\sqrt{R^2+z^2}`$. But since higher order geometrical effects are ignored in our model this is not taken into account. Fortunately the gas+dust density at large $`z/R`$ is very small, and this effect has no practical influence on the resulting spectra.

Finally, we show in Fig.4 the radial dependence of the equatorial plane temperature (left panel) and of the pressure scale and surface height of the disk for the VEF and MEMO solutions (right panel). As was already known from Fig.3, the midplane temperature is lower than in the MEMO approximation at most radial intervals, but becomes higher at the outer edge. However, at the very inner edge it is again the same in both approaches. This is because at these radii the vertical optical depth of the disk becomes so large that even at millimeter wavelengths there is no possibility for the gas to cool, and the disk becomes isothermal. Fig.4, right panel, shows the radial profile of the pressure scale height, which follows approximately that of the midplane temperature, as discussed, and of the surface height. It is interesting to note how the MEMO approximation and the full transfer model give almost equal values of $`H_s`$. This is a result of the higher temperatures at higher elevations above the midplane. This can be seen from the fact that the upper part of the temperature rise, which is where the direct stellar radiation is absorbed, is located at about the same $`z`$ in both models.

With these differences between the simple (MEMO) and the complex (VEF) method, one may wonder whether a compromise between the two methods may give results that are close enough to the real one. Even with the efficiency of the VEF method, solving the full angle-dependent radiative transfer problem is a time-consuming process. It may therefore not be very practical for e.g. the computation of the evolution and/or the dynamics of the disk, where the radiative transfer problem needs to be solved thousands of times. In Fig.5 we show the same slice as in Fig.1. Here we added the computed temperature profiles for two alternative methods: one in which the full $`\mu `$-dependence is solved, but with Planck- and Rosseland mean opacity replacing the $`\kappa _J`$ and $`\kappa _H`$, and one in which the Eddington approximation is used, but with the full frequency-dependent opacities.

One sees that the former does not improve much on the MEMO method. But the latter in fact gets very close to the ’right’ answer. This shows that most of the problems with the MEMO method originate from the wrong treatment of mean opacities, while the use of the Eddington approximation iself ($`f_\nu =1/3`$ and $`H(z=z_{\mathrm{max}})=J(z=z_{\mathrm{max}})/\sqrt{3}`$) does not cause major problems. We conclude that an approximate method, in which frequency-dependent radiative transfer is done in the Eddington approximation, is in fact a relatively accurate method. This method avoids the CPU-time consuming full frequency-angle dependent radiative transfer, and is therefore faster than the full VEF method. Nevertheless, in this paper we adopt the complete VEF method, since it is still fast enough for our purpose, and it is exact, at least in the 1-D approximation made here.

### 3.1 Effect on the spectral energy distribution

The new vertical structure models have a spectral energy distribution (SED) that is not very different from the ones predicted by the vertical structures from the MEMO approximation. In fact, they do not even differ much from the even simpler semi-analytic treatment of the Chiang & Goldreich model, if certain modifications are made to the latter (see Chiang et al. 2001 and Dullemond, Dominik & Natta 2001, henceforth DDN01).

This may be rather surprising, since the disk structure is so drastically different, with equatorial temperatures differing up to 70%. The explanation of this lies in the fact that the total emitted flux of the disk is prescribed by energy balance. At each radius it is expected that the disk emits both an optically thin component and an optically thick one. Each component carries half of the emitted flux, which equals the total absorbed flux. This means for instance that the optically thick component (at far-IR wavelengths) must have an effective temperature that is the same in all models. In the full transfer models the effective temperature is therefore the same as the effective temperature of the Chiang & Goldreich model, even though the spectrum may not be a perfect blackbody.

There are however some differences that are worth noting. The full transfer models predict higher fluxes in the mid and far-infrared than MEMO models. The difference is of the order of 30% in the wavelength range between 50 and 500 $`\mu `$m for a T Tauri star, and between 50 and 100 $`\mu `$m for a hotter Herbig Ae star. However, the shape of the SED is not very different, even if the full transfer models fall slightly more rapidly towards the far infrared. The SEDs of the CG97/DDN01 models are somewhat flatter in the same range of wavelengths, since they are composed of two equally strong components (see the case of the HAe star, where the two “bumps” are clearly visible). This is a direct consequence of the discrete two-layered structure of the CG97/DDN01 model. In the full transfer model the surface layer has a continuous temperature gradient, which smoothly matches to the photosphere. This produces the smoother SEDs at far-IR wavelengths. The bumps in the SEDs of the CG97/DDN01 model are therefore an artifact of the adopted two-layers simplification.

The region aroung the 10 micron feature are virtually not affected by the full treatment of radiative transfer. The emission in this region is dominated in all models by emission from the optically thin surface layers, whose properties are quite independent of what happens deep below in the disk, and is therefore much less affected by the different methods of solving the radiative transfer.

### 3.2 Effects on molecular line intensities

Knowledge of the correct vertical structure of disks (temperature and density) is of great importance for predicting the abundances of different molecular species and the intensity of their lines. In a separate paper we discuss the differences between various models in detail, in the context of molecular line observations from T Tauri and Herbig Ae stars, for which the SED is sufficiently well known (Van Zadelhoff & Dullemond in prep). Here we confine ourselves to a demonstration of the effects on the intensity of the lines of CO and its two main isotopomers 13CO and C18O, assuming that gas and dust temperatures are the same. The CO molecule can be used as a temperature tracer, and is therefore suited to probe the differences between the two models.

A first insight of how line intensities depend on the disk model can be obtained from Fig. 7, which shows how the temperature varies as a function of the gas column density (for our template T Tauri star model and $`R`$=220 AU). If we assume for simplicity that the line optical depth is proportional to the gas column density, and that the brightness temperature in the line is roughly equal to the gas temperature where the optical depth in the line is $`\tau 2/3`$, one can immediately see that different models will predict the same brightness temperature for lines that are very optically thick, but that MEMO models will underestimate significantly the brightness temperature of more optically thin lines, which form more deeply in the disk where the predicted temperatures differ.

For a true comparison of models to the observations a number of effects have to be taken into account, including NLTE effects, size of the telescope beam, and chemical abundances of the molecular species. The NLTE level populations have been calculate using a 2D Monte Carlo code (Hogerheijde & van der Tak 2000). The reason for calculating the full NLTE level populations instead of LTE is due to the low densities in the upper layer of the disks (see van Zadelhoff et al. 2001). From these populations, the emitted line profiles were constructed, convolving the disk with a beam of 4.0″, the apparent size of the source on the sky at a distance of 150 pc. In all cases the disk is assumed to be seen under an inclination of 60°. The abundance ratios with respect to H2 are assumed to be constant (\[CO\]/\[H2\]=8.$`\times 10^5`$) in the entire disk. The abundances of the two isotopomers follow the isotopic ratios of \[12C\]/\[13C\]=60 and \[12C\]/\[18C\]=500, similar to the solar neighborhood.

The molecular line emission results are shown in Fig. 8 and summarized in Table 1, in which the line emission integrated over velocity is given. The three isotopomers, due to their differences in abundance become optically thick at different heights in the disk. CO reaches the $`\tau `$=1 plane close to the disk surface, where temperature and density are very similar in the two solutions. The line profiles shows only small changes between models with differences up to 5% in the integrated line intensities. The 13CO becomes optically thick in the intermediate layer showing slightly larger differences between the two models (up to 15%). C18O has a optical depth of a few, and its lines form in the region of the disk where the VEF and MEMO solutions differ more. The differences between the two models are $`30`$%. The CO isotopomers can be used to measure the vertical temperature profile in disks (e.g. Guilloteau & Dutrey 1994), and in this context the differences between different radiation transfer methods become significant.

For each isotope, the difference between the more optically thick 3-2 and the less optically thick 6-5 follows the same pattern, as seen in Fig. 8. However, the difference in the ratio of integrated intensities between models is small when compared with the typical observational errors (10-20 %).

These calculations were performed assuming constant abundances throughout the disk. In reality this will not be the case due to dissociation of molecules in the upper layers by stellar and interstellar radiation and freeze-out of molecules on the grain surface at $`T<20`$K (Aikawa et al. 2002). This latter process would enhance the difference between the two models as more freeze-out is expected from the colder MEMO model.

## 4 Discussion

The model described in this paper is relatively elementary in the sense that all kinds of detailed microphysics are ignored. But within the limited physics with which it is defined, and the geometrical simplifications made, it is a reasonably accurate solution. The 1-D vertical radiative transfer and the vertical pressure balance are solved without any approximations other than the discretization itself. In this respect it is a step forwards compared to the models of passive disks published in the literature so far. The models worked out in this paper can be downloaded from a website111http://www.mpa-garching.mpg.de/PUBLICATIONS/DATA/ radtrans/diskmodel/, which also features a couple of one-annulus test setups for the radiative transfer problem in these disks.

The advantage of ignoring all the complex and detailed microphysics is that we have a well-defined disk model with very few parameters. And all the features of this model can be explained in terms of two pieces of physics only: radiative transfer and hydrostatics. This simple model can then serve as a basic model to which more physics can be added later, such as the inclusion of dust settling and coagulation, the inclusion of dust scattering, a treatment of internal heat processes such as viscous dissipation, etc. A study of the effect of dust scattering on the disk’s structure and the outcoming spectrum is under way (Dullemond & Natta in prep.).

There are however a few uncertainties concerning the global structure of the disk, even within the simple problem definition used in this paper. These have to do with the possibility of self-shadowing. The present calculations start from the Ansatz that the disk has a flaring shape, so that the star’s radiation can illuminate the disk at every radius. It then finds a solution that is consistent with this Ansatz. But in reality perhaps part of the disk might be non-flaring and reside in the shadow of the inner regions of the disk. This could be the result of the history of formation of the disk, or it may have to do with the outer parts of the disk becoming too optically thin to sustain a positive derivative of $`H_s/R`$. But a self-shadowed region in the disk could also be the end-product of an instability (Dullemond 2000; Chiang 2000). In these cases the Ansatz of positive flaring index is not confirmed, and the solutions might be different from the smooth flaring disk solutions presented here. In order to investigate the possibility of such alternative disk solutions, one must include full 2D radiative transfer and full 2D hydrostatics into the model calculations. Due to the enormous complexity of this problem, in particular at high optical depths, we do not venture into this direction at present.

The models presented in this paper describe the middle and outer parts of a passive flaring circumstellar disk. The very inner parts have a different structure than the one described here. Very close to the star the dust has evaporated and one is left with a pure gas disk. If the gas is optically thin, the inner rim of the dusty part of the disk is directly exposed to the radiation of the star. It will therefore be much hotter than the rest of the disk, that is only irradiated on the surface. This inner rim will therefore puff up and consistute a new component in the spectral energy distribution of the disk (Natta et al. 2001). A physical description of this inner rim, and the effects it has on the complete disk structure is described by Dullemond, Dominik & Natta (2001). But a full 2D model for this inner rim again requires 2D radiative transfer, and is beyond the scope of this paper.

## 5 Conclusion

We have presented in this paper the results of calculations of the vertical structure of irradiated circumstellar disks in hydrostatic equilibrium. The main difference with previous models (D’Alessio et al., Malbet et al. etc.) is the full angle- and frequency-dependent treatment of radiative transfer. It turns out that the frequency-dependent treatment of the problem is of crucial importance for obtaining the correct vertical structure (temperature and density) of the disk. Over most of the disk the midplane temperature is significantly lower than the value predicted by a simpler radiative transfer method that uses Planck and Rosseland mean opacities. The correct treatment of the angular dependence is of lesser importance. In applications where computation time is critical, the use of the Eddington approximation is reasonable, as long as the full frequency dependence of the radiation field is taken into account.

If the disk model is used as input for further calculations, the need to implement a reliable radiation transfer method is obvious. This is, for example, the case when computing the expected intensity and profile of molecular lines, as shown in §4 for the lines of the CO isotopomers, which form at different depth in the disk. As more realistic chemistry and line transfer models are used, the differences between different disk models is likely to become even more relevant (see van Zadelhoff and Dullemond, in prep.), and the VEF models should be adopted to constrain physical and chemical parameters in circumstellar disk environments.

The SED is also affected by the change in vertical structure, but the differences are not very large. In general the SED compares reasonably well with the results from both the MEMO approach and the CG97/DDN01 model. This means that the model of Chiang & Goldreich (1997), with improvements described by Chiang et al. (2001) and Dullemond, Dominik & Natta (2001), is a fairly robust model to compute the SED of T Tauri stars and Herbig Ae/Be stars.

###### Acknowledgements.

We are grateful to E. Krügel for assisting us with the testing of our radiative transfer program and to M. Hogerheijde and F. van der Tak for use of their 2D Monte Carlo code for line transfer. CPD acknowledges support from the European Commission under TMR grant ERBFMRX-CT98-0195 (‘Accretion onto black holes, compact objects and prototars’). Astrochemistry in Leiden is supported through a Spinoza grant from the Netherlands Organization for Scientific Research (NWO).

## Appendix A Radiative transfer with Variable Eddington Factors

The idea at the basis of the variable Eddington factors method can be traced back to the sixties (see e.g. Mihalas & Mihalas 1984 and references therein). For accretion disk theory the equations for this method were presented by Hubeny (1990) and Malbet & Bertout (1991). Here we present a slight variation of this approach. The advantage of our method is that it can be used up to any optical depth even in cases where the opacity is non-grey and the internal energy dissipation (and therefore the net flux) within the disk is zero.

Define the mean intensity $`J_\nu `$, the Eddington flux $`H_\nu `$ and the second moment of radiation $`K_\nu `$ as follows:

$`J_\nu `$ $`=`$ $`{\displaystyle \frac{1}{2}}{\displaystyle _1^1}I_{\mu ,\nu }𝑑\mu `$ (11)
$`H_\nu `$ $`=`$ $`{\displaystyle \frac{1}{2}}{\displaystyle _1^1}I_{\mu ,\nu }\mu 𝑑\mu `$ (12)
$`K_\nu `$ $`=`$ $`{\displaystyle \frac{1}{2}}{\displaystyle _1^1}I_{\mu ,\nu }\mu ^2𝑑\mu .`$ (13)
By multiplying the transfer equation (6) by resp. $`\mu ^k`$, and integrating over $`d\mu `$ one obtains the $`k`$-th moment equation of radiative transfer. The first two moment equations are:

$`{\displaystyle \frac{dH_\nu }{dz}}`$ $`=`$ $`\rho \kappa _\nu (B_\nu (T)J_\nu )`$ (14)
$`{\displaystyle \frac{dK_\nu }{dz}}`$ $`=`$ $`\rho \kappa _\nu H_\nu .`$ (15)
One can write the second moment $`K_\nu `$ as a dimensionless factor times the mean intensity $`J_\nu `$:

$$K_\nu =f_\nu J_\nu ,$$
(16)
where $`f`$ is the Eddington factor. For isotropic radiation the Eddington factor is equal to $`f=1/3`$. For purely beamed radiation ($`I_{\mu ,\nu }=0`$ for $`\mu 1`$) this factor is $`f=1`$. If one makes an assumption for this value, then Eqs.(14,15) form a closed set of equations, which can be solved subject to boundary conditions. Often the assumption is made that the radiation field will be more or less isotropic, and therefore that $`f_\nu =1/3`$. This is the Eddington approximation, and stands at the basis of many approximate transfer codes.

But if one had some way of computing (rather than guessing) the Eddington factor, then the moment equations are not an approximation anymore, and in fact will be exactly equivalent to the full transfer equation. One way of doing so is to iteratively switch between the moment equations and the real transfer equation. Start with a guess for $`f_\nu (z)`$ (take it for instance $`1/3`$), and solve the moment equations (14,15) to find the temperature. Then integrate the formal transfer equation (6) using the current temperature. Using the resulting radiation field $`I_{\mu ,\nu }(z)`$ one can then compute the Eddington factor $`f_\nu (z)`$. Once this is done, one starts all over again, using the newly computed Eddington factor. This process is repeated until a converged solution for $`T(z)`$ is reached. Since in this procedure the $`f_\nu (z)`$ is computed on the fly, one can be sure that this solution is a real solution, and not an approximate one. It is important to note that if this solution is inserted into the full set of transfer equations (6,7,8), then one will see that these are solved as well. In other words: the moment equations were used to find a solution to the real transfer equations.

The advantage of the moment equations (14,15) is that they can be solved directly, using a two-point boundary value approach. In fact, the frequency-integrated moment equations are even simpler to solve directly. The frequency-integrated version of Eq.(14) is:

$$\frac{dH}{dz}=\rho \left(\kappa _P(T)\frac{\sigma }{\pi }T^4\kappa _JJ\right),$$
(17)
where $`\kappa _P(T)`$ is the Planck mean opacity and

$$H=_0^{\mathrm{}}H_\nu 𝑑\nu .$$
(18)
Using the concept of energy conservation, one can replace the right-hand-side of Eq.(17) with the source term:

$$\frac{dH}{dz}=\frac{q}{4\pi }$$
(19)
This equation can be trivially integrated from the equator (starting with $`H=0`$) up to $`z=z_{\mathrm{max}}`$. At that point we can compute the value of the frequency integrated mean intensity $`J`$:

$$J(z=z_{\mathrm{max}})=H(z=z_{\mathrm{max}})/\psi ,$$
(20)
where

$$J=_0^{\mathrm{}}J_\nu 𝑑\nu ,$$
(21)
and $`\psi `$ is the ratio of $`H/J`$ as computed from the full radiative transfer. Using the frequency-integrated version of Eq.(15),

$$\frac{d(fJ)}{dz}=\rho _0^{\mathrm{}}\kappa _\nu H_\nu 𝑑\nu ,$$
(22)
(where $`f`$ is the $`J`$-mean of $`f_\nu `$), and starting from the value of $`J(z=z_{\mathrm{max}})`$ computed above, one can integrate back towards the equator to find $`J(z)`$. The temperature then follows from

$$\kappa _P(T)\frac{\sigma }{\pi }T^4=\kappa _JJ+\frac{q}{4\pi },$$
(23)
where $`\kappa _J`$ is the $`J_\nu `$-mean of the opacity computed from the full transfer.

By iterating on the above procedure, one can quickly find solutions to the transfer equation. An even quicker convergence can be reached by applying a linear convergence amplifier like Ng’s algorithm (Ng 1974).

Typically one needs about 70 gridpoints in $`z`$, 40 points in $`\mu `$ and, dependent on the kind of opacity table and the wavelength range, about 60 points in $`\nu `$. It is important to make sure that there are enough $`\mu `$ gridpoints near $`\mu =0`$, because in plan-parallel geometry most of the radiation is at small values of $`\mu `$. We use a logaritmic grid in $`\mu `$ spanned between $`\mu _00.01`$ and $`1`$.

It is important to note why, in Eq.(22), we did not use a flux-weighted mean opacity $`\kappa _H`$ (or Rosseland opacity) multiplied by the frequency-integrated Eddington flux $`H`$, as was suggested for instance by Malbet & Bertout (1991). The reason is that in a disk without internal heat production, the total flux $`H`$ is zero. The frequency-dependent flux $`H_\nu `$, on the other hand, is non-zero, consisting of downward diffusive flux at short wavelengths and an equal amount of upward flux at long wavelengths. For non-grey opacities one therefore has a virtually zero $`H`$ close to the equator, but a very non-zero $`_0^{\mathrm{}}\kappa _\nu H_\nu 𝑑\nu `$ at the same location. An evaluation of $`\kappa _H=_0^{\mathrm{}}\kappa _\nu H_\nu 𝑑\nu /_0^{\mathrm{}}H_\nu 𝑑\nu `$ would therefore yield a diverging number, and the algorithm becomes unstable, whereas the equation in the form of Eq.(22) yields a stable algorithm.

The method described in this section has proven to work well, and reaches a solution to an accuracy of $`10^4`$ in temperature within 11 iterations and to $`10^8`$ within 22 iterations, independent of the optical depth. On a Pentium III with 850 MHz clock frequency a typical radiative transfer problem is solved in about 5 seconds.

### A.1 Testing

In order to verify the reliability of the code, a few tests are performed. First of all we do internal consistency checks. Since the VEF method computes the $`J`$ and $`H`$ from both the full transfer as well as the moment equations, one can a-posteriori check whether the two are the same. For all our models this turns out to be indeed the case. Then we compare the solution to the results from an independent ALI code. For those cases in which the ALI in fact converges to a satisfiable degree, we find that the solution are indeed very close to each other.

To further convince ourselves of the validity of the solutions provided by the VEF method, we applied the VEF method to a problem with grey opacities. For such a problem it is known that the MEMO method should work reasonably well. And as expected, the VEF method is consistent with the MEMO results.

### A.2 Convergence

In order to demonstrate the advantage of the VEF method over methods like LI and ALI, we show here the convergence history for a test problem. Consider a circumstellar disk which has $`\mathrm{\Sigma }=10^3\text{g}/\text{cm}^2`$ in dust mass at 1 AU from a central star of $`T_{\mathrm{eff}}=3000\text{K}`$ and $`R_{}=2.0R_{}`$. The vertical pressure scale height of the disk is assumed to be $`H_p=0.028\text{AU}`$, and the density profile is assumed to be Gaussian. The flaring index equals $`0.2`$, and the opacity table used is for astronomical silicate (Draine & Lee 1984). In Fig. 9 the convergence history for the temperature is plotted using three different methods: LI with Ng, ALI with Ng and the VEF method with Ng. The convergence history is only shown for the first 100 iterations, within which only the VEF method converges. After about 400 iterations, the ALI+Ng method in fact also converges to within $`1\%`$ of the solution, but the LI+Ng method does not converge even after 1000 iterations.

It should be noted that the method of Accelerated Lambda Iteration has not been widely used for dust continuum transfer. It has been developed mainly for line transfer, for which it has proven to be reasonably effective. We adapted the algorithm for use with dust continuum transfer, and chose as our approximate operator the diagonal of the full lambda operator.