text:
# An accelerated Monte Carlo method to solve two-dimensional radiative transfer and molecular excitation

## 1 Introduction

The dense and cool material in the interstellar medium of galaxies plays an important role in the life cycle of stars, from the earliest phases of star formation to the shells around evolved stars and the gas and dust tori around active galactic nuclei. Line emission from atoms and molecules, and continuum emission from dust particles, at radio, (sub) millimetre and infrared wavelengths are indispensable tools in the study of a wide variety of astrophysical problems. This is illustrated by the large number of infrared and submillimetre observatories planned for the near future, such as the Smithsonian Millimeter Array (SMA), the Atacama Large Millimeter Array (ALMA), the Far-Infrared and Submillimetre Space Telescope (FIRST) and the Stratospheric Observatory for Infrared Astronomy (SOFIA).

An essential step in the interpretation of the data from these instruments is the comparison with predicted emission from models. This paper presents a numerical method to solve the radiative transfer and molecular excitation in spherically symmetric and cylindrically symmetric source models. At the comparatively low densities of interstellar gas, the excitation of many molecules is out of local thermodynamic equilibrium (LTE), and the transfer of line (and continuum) radiation plays a significant role in determining the molecular excitation Leung and Liszt (1976); Black (2000). Geometry thus becomes an important element, and the high spatial resolution of current and future instruments often demands that at least two-dimensional (axisymmetric) source structures are considered. In the implementation of our method discussed in this paper, we have limited the source structure to spherical and cylindrical symmetries. The large and often multidimensional parameter space further requires a fast and reliable method, which needs to be easily applicable to many different astrophysical problems.

This need for reliable and flexible tools calls for the use of Monte Carlo techniques, where the directions of integration are chosen randomly. This approach was first explored by Bernes Bernes (1979) for non-LTE molecular excitation; later, Choi et al. Choi et al. (1995), Juvela Juvela (1997) and Park & Hong Park and Hong (1998) augmented it and expanded it to multiple dimensions. However, Monte Carlo methods can be quite slow, especially at large optical depths ($`\tau 100`$), which has limited their use so far. We will show that this problem can be overcome by using a technique inspired on Accelerated Lambda Iteration: the local radiation field and excitation are solved self-consistently and separated from the overall radiative transfer problem (see §3.4). The greatest virtue of our code is its ability to deal with a wide variety of source models for many atomic and molecular species, with or without a dust continuum. Although for any individual problem a somewhat more efficient method could be constructed (§ 3.5), the Monte Carlo approach frees the user from having to fine-tune such a method and allows the user to focus on the astrophysics of the problem at hand.

This paper does not discuss the influence of radiative transfer on the source structure, through the thermal balance, ionization and chemistry (Takahashi et al., 1983; Ceccarelli et al., 1996; Doty and Neufeld, 1997, for example). However, our code is suited to become part of an iterative scheme to obtain self-consistent solutions for the source structure including radiative transfer and molecular excitation.

Throughout this paper, examples from studies of star formation will serve to illustrate the various topics – and to stress the link with the analysis of observations. Sect. 2 introduces a simple, spherically symmetric model of the core of an interstellar cloud, collapsing to form a star. Sect. 3 then discusses the coupled problem of radiative transfer and molecular excitation. It introduces the two most commonly used solution methods, and shows that these are closely related. This opens the possibility of constructing a hybrid method which combines the benefits of both; the implementation of this combined approach in our code is deferred to the Appendix. The paper continues by exploring the capabilities of our code through a number of astrophysically relevant examples, based on extensions of the simple one-dimensional model of Sect. 2. A brief summary concludes the paper in Sect. 5.

## 2 An illustrative, one-dimensional model

The formation of stars occurs in dense condensations within interstellar molecular clouds, which collapse under the influence of their own gravity. A widely used theoretical description of this process, constructed by Shu Shu (1977), starts with the singular isothermal sphere,

$$\rho (r,t=0)=\frac{a^2}{2\pi G}r^2.$$
(1)
Here, $`\rho `$ is the density as function of radius $`r`$ and time $`t`$, $`a`$ is the isothermal sound speed, and $`G`$ is the gravitational constant.

At $`t=0`$ collapse starts at the center ($`r=0`$). After a time $`t`$, all regions $`r<at`$ are collapsing, with speed $`v(r,t)`$ increasing from 0 at $`r=at`$ to free-fall, $`vr^{0.5}`$, well within this ‘collapse expansion wave’ ($`rat`$). Shu Shu (1977) constructed a solution for the density and velocity field of the collapsing core which is self-similar in the spatial coordinate $`xat`$. The density follows a power-law behaviour as function of radius, with $`\rho r^{1.5}`$ for $`rat`$, $`\rho r^1`$ just inside $`r=at`$, and the undisturbed $`\rho r^2`$ outside $`r=at`$ (Fig. 1).

Many authors have tested this model against observations of cloud cores and envelopes around young stellar objects (YSOs), e.g. Zhou et al. Zhou et al. (1993), Choi et al. Choi et al. (1995), Ward-Thompson et al. Ward-Thompson et al. (1996) and Hogerheijde & Sandell Hogerheijde and Sandell (2000). Especially the spectral-line signature of collapse (Fig. 1d) has received much attention as a probe of ongoing collapse, although this signature is shared by all collapse models and is not unique to the particular model described here. The exact line shape, however, depends quantitatively on the adopted model.

The interpretation of this signature needs non-LTE radiative transfer. Both collisional and radiative processes can excite molecules, and for each transition a *critical density* can be defined where the two are of equal importance. At lower densities radiation dominates, while at higher densities collisions drive the level populations to thermodynamic equilibrium. The large range of densities of star forming cores ensures that many molecules and transitions will go through the entire range of excitation conditions, while line emission will have a significant impact on the excitation at the intensities and opacities expected for typical abundances of many species, not only locally but throughout the envelope (Fig. 1c).

In the following we will use this model to illustrate our method of solving the coupled problem of radiative transfer and excitation. In particular, we will consider emission lines of HCO+ and H13CO+, which are readily observed and often used as tracers of dense gas. The strong $`J=10`$, $`32`$ and $`43`$ lines at millimetre wavelengths have critical densities of $`2\times 10^5`$, $`4\times 10^6`$, and $`1\times 10^7`$ cm-3, using the molecular data in Table 1. We assume an abundance of HCO+/H$`{}_{2}{}^{}=5\times 10^9`$ and an isotopic ratio of 1:65 for H13CO+: HCO+. The sound speed of the adopted model is $`a=0.24`$ km s-1, its age $`t=1\times 10^5`$ yr, and its outer radius 8000 AU. The total mass of the model is 0.73 M. The kinetic temperature follows $`T_{\mathrm{kin}}=30\mathrm{K}(r/1000\mathrm{AU})^{0.4}`$, appropriate for a centrally heated envelope at a luminosity of $`2`$ L (Adams et al., 1987, e.g.). The turbulent line width of 0.2 km s-1 is smaller than the systematic velocities except in the outermost part (Fig. 1b).

## 3 Solving radiative transfer and molecular excitation

### 3.1 The coupled problem

The equation of radiative transport reads, in the notation of Rybicki & Lightman Rybicki and Lightman (1979),

$$\frac{dI_\nu }{ds}=\alpha _\nu I_\nu +j_\nu ,$$
(2)
or, equivalently,

$$\frac{dI_\nu }{d\tau _\nu }=I_\nu +S_\nu .$$
(3)
Here, $`I_\nu `$ is the intensity at frequency $`\nu `$ along a particular line of sight parameterized by $`ds`$, $`\alpha _\nu `$ is the absorption coefficient in units cm-1, and $`j_\nu `$ the emission coefficient with units erg s-1 cm-3 Hz-1 sr-1. The second form of the equation is a useful change of variables, with the source function $`S_\nu j_\nu /\alpha _\nu `$ and the optical depth $`d\tau _\nu \alpha _\nu ds`$. We consider both molecules and dust particles as sources of emission and absorption ($`j_\nu =j_\nu (\mathrm{dust})+j_\nu (\mathrm{gas})`$; $`\alpha _\nu =\alpha _\nu (\mathrm{dust})+\alpha _\nu (\mathrm{gas})`$), but ignore scattering. Although not impossible to include in our code, scattering effects are usually negligible at wavelengths longer than mid-infrared.

When $`\alpha _\nu `$ and $`j_\nu `$ are known at each position in the source, the distribution of the emission on the sky simply follows from ray tracing. However, in many cases, $`\alpha _\nu `$ and $`j_\nu `$ will depend on the local mean intensity of the radiation field

$$J_\nu \frac{1}{4\pi }I_\nu 𝑑\mathrm{\Omega }.$$
(4)
Here, $`J_\nu `$ is the average intensity received from all solid angles $`d\mathrm{\Omega }`$, and $`I_\nu `$ is the solution of Eq. (2) along each direction under consideration. The latter integration extends formally to infinity, but in practice only to the edge of the source with any incident isotropic radiation field like the cosmic microwave background (CMB) as boundary condition.

For thermal continuum emission from dust, $`j_\nu (\mathrm{dust})`$ and $`\alpha _\nu (\mathrm{dust})`$ are simply given by

$$j_\nu (\mathrm{dust})=\alpha _\nu (\mathrm{dust})B_\nu (T_{\mathrm{dust}}),$$
(5)
where $`B_\nu `$ is the Planck function at the dust temperature $`T_{\mathrm{dust}}`$, and

$$\alpha _\nu (\mathrm{dust})=\kappa _\nu \rho _{\mathrm{dust}},$$
(6)
where $`\kappa _\nu `$ is the dust opacity in cm-2 per unit (dust) mass and $`\rho _{\mathrm{dust}}`$ is the mass density of dust. Our code can use any description of $`\kappa _\nu `$ (Ossenkopf and Henning, 1994; Pollack et al., 1994; Draine and Lee, 1984; Mathis et al., 1977, e.g.).

In the case of emission and absorption in a spectral line, $`\alpha _\nu ^{ul}(\mathrm{gas})`$ and $`j_\nu ^{ul}(\mathrm{gas})`$ are determined by absorption and emission between radiatively coupled levels $`u`$ and $`l`$ with populations (in cm-3) $`n_u`$ and $`n_l`$. The energy difference between levels $`\mathrm{\Delta }E=E_uE_l`$ corresponds to the rest frequency of the transition, $`\nu _0=\mathrm{\Delta }E/h`$, where $`h`$ is Planck’s constant. The emission and absorption coefficients between levels $`u`$ and $`l`$ are strongly peaked around $`\nu _0`$ with a frequency dependence described by a line-profile function $`\varphi (\nu )`$,

$`j_\nu ^{ul}(\mathrm{gas})`$ $`=`$ $`{\displaystyle \frac{h\nu _0}{4\pi }}n_uA_{ul}\varphi (\nu ),`$ (7)
$`\alpha _\nu ^{ul}(\mathrm{gas})`$ $`=`$ $`{\displaystyle \frac{h\nu _0}{4\pi }}(n_lB_{lu}n_uB_{ul})\varphi (\nu ).`$ (8)
The Einstein $`A_{ul}`$, $`B_{lu}`$, and $`B_{ul}`$ coefficients determine the transition probabilities for spontaneous emission, absorption, and stimulated emission, respectively, and depend on molecule. In most interstellar clouds the line profile is dominated by Doppler broadening due to the turbulent velocity field

$$\varphi (\nu )=\frac{c}{b\nu _0\sqrt{\pi }}\mathrm{exp}\left(\frac{c^2(\nu \nu _0)^2}{\nu _0^2b^2}\right),$$
(9)
where the turbulence is assumed to be Gaussian with a full width $`b`$. In the presence of a systematic velocity field, the line profile is angle-dependent and the projection of the local velocity vector onto the photon propagation direction enters $`(\nu \nu _0)`$.

Together, collisions and radiation determine the level populations through the equation of statistical equilibrium,

$$\begin{array}{c}n_l\left[_{k<l}A_{lk}+_{kl}(B_{lk}J_\nu +C_{lk})\right]=\hfill \\ _{k>l}n_kA_{kl}+_{kl}n_k(B_{kl}J_\nu +C_{kl}).\hfill \end{array}$$
(10)
The collision rates $`C_{kl}`$ depend on the density and the collisional rate coefficients of molecular hydrogen and other collision partners, and on temperature through the detailed balance of the up- and downward coefficients. Eq. (10) can be easily solved through matrix inversion for each position in the source provided the radiation field $`J_\nu `$ is known. However, $`J_\nu `$ contains contributions by the CMB, dust and spectral lines, and since the spectral line term depends on the level populations through Eqs. (2), (7) and (8), the problem must be solved iteratively. Starting with an initial guess for the level populations, $`J_\nu `$ is calculated, statistical equilibrium is solved and new level populations are obtained; through the Monte Carlo integration, the new populations yield a new value for $`J_\nu `$, after which the populations are updated; etc., until the radiation field and the populations have converged on a consistent solution.

When the physical conditions do not vary much over the model, an approximate value of $`J_\nu `$ can be found from the local conditions alone. This idea is the basis of the Large Velocity Gradient method, the Sobolev method, microturbulence, or the escape probability formalism (Sobolev, 1960; Goldreich and Kwan, 1974; Leung and Liszt, 1976; de Jong et al., 1980, e.g.). Also, in specific geometries, the integration over all solid angles and along the full length of the line of sight of Eqs. (3) and (4) can be greatly reduced, making the problem tractable. This sort of technique has most application in stellar and planetary atmospheres; the Eddington approximation is an example.

However, in many astrophysical situations including the example of § 2, such simplifications cannot be made, and Eqs. (3) and (4) need to be fully solved to get $`J_\nu `$. Compared to the relative ease with which statistical equilibrium can be solved (Eq. 10), obtaining $`J_\nu `$ becomes the central issue. Direct integration of Eqs. (3) and (4) with, e.g., Romberg’s method, is infeasible for realistic models, but based on a finite set of directions a good approximation of $`J_\nu `$ can be obtained. The next two sections describe two different methods to choose this set and construct $`J_\nu `$ in this way.

### 3.2 Constructing $`J_\nu `$ and the $`\mathrm{\Lambda }`$-operator

For computational purposes, source models are divided into discrete grid cells, each with constant properties (density, temperature, molecular abundance, turbulent line width, etc.). It is also assumed that the molecular excitation can be represented by a single value in each cell, which requires instantaneous spatial and velocity mixing of the gas. Appropriate source models have small enough cells that the assumption of constant excitation is valid. The systematic velocity field is the only quantity that is a vector field rather than a scalar field, and in our code it is allowed to vary in a continuous way within each cell. We divide the integration along a ray into subunits within a cell to track the variation of the velocity projected on the ray.

Such a gridded source model lends itself easily to the construction of a finite set of integration paths to build up $`J_\nu `$. The average radiation field can be thought of as the sum of the emission received in cell $`i`$ from each of the other cells $`j`$ after propagation through the intervening cells and weighted with the solid angle subtended by each of these cells $`j`$ as seen from cell $`i`$. The combination of radiative transfer and statistical equilibrium can be written as

$$J_\nu =\mathrm{\Lambda }[S_{ul}(J_\nu )].$$
(11)
This equation states that the radiation field is given by an operator $`\mathrm{\Lambda }`$ acting on the source function $`S_{ul}`$, which depends on the level populations and hence $`J_\nu `$ (Eqs. 7, 8, 10). Considering the narrow frequency interval around the transition $`u`$$`l`$, we have replaced $`S_\nu `$ by $`S_{ul}[j_{\nu _0}(\mathrm{dust})+j_\nu ^{ul}(\mathrm{gas})𝑑\nu ]/[\alpha _{\nu _0}(\mathrm{dust})+\alpha _\nu ^{ul}(\mathrm{gas})𝑑\nu ]`$. This corresponds to the assumption of instanteous redistribution of excitation mentioned above. In a gridded source model, one can think of $`\mathrm{\Lambda }`$ as a matrix describing how the radiation field in cell $`i`$ depends on the excitation in all other cells. The elements in the matrix then represent the radiative coupling between cell pairs.

Eq. (11) can be solved iteratively, where an updated value of $`J_\nu `$ is obtained by having $`\mathrm{\Lambda }`$ operate on the previous populations, $`S_{ul}^{}`$,

$$J_\nu =\mathrm{\Lambda }[S_{ul}^{}(J_\nu )].$$
(12)
Since $`S_{ul}^{}`$ is already known, this only involves matrix multiplication, compared to the much more expensive matrix inversion required to solve Eq. (11). Because of this elegant notation, iterative schemes for non-LTE excitation and radiative transfer are commonly referred to as $`\mathrm{\Lambda }`$-iteration, even if no $`\mathrm{\Lambda }`$-operator is ever actually constructed. These methods share the use of the same set of rays throughout the calculation, in contrast to Monte Carlo methods, which use random rays as discussed in §§3.3 and 3.5.

Constructing the $`\mathrm{\Lambda }`$-operator in multidimensional source models is taxing on computer memory because of all the possible connections to keep track of simultaneously. Techniques exist to reduce the number of elements Dullemond and Turolla (2000), but these are complex and may require some fine-tuning for individual source geometries. Alternatively, computer *memory* can be exchanged for computing *time* by solving the problem serially, calculating the radiation field in each of the cells due to the other cells one at a time.

### 3.3 The Monte Carlo method

One way of solving Eq. (12) is to directly sum the contribution from all other cells to the radiation field in each of the individual cells. This corresponds to replacing the integral in Eq. (4) by a summation. With a judicially chosen fixed set of directions or rays, as most $`\mathrm{\Lambda }`$-iteration codes do, a good approximation of $`J_\nu `$ can be found in this way (Phillips, 1999, e.g.). However, this procedure requires care, since the necessary angular sampling depends, in principle, on characteristics of the excitation solution of the problem at hand.

Since our aim is to construct a method that can be applied to many different source models without too much fine-tuning, we adopt the Monte Carlo approach to obtain $`J_\nu `$. Analogous to the Monte Carlo method to solve the definite integral of a function \[see chapter 7 of Press et al. Press et al. (1992) for a discussion of Monte Carlo integration, and further references\], Eq. (4) can be approximated by the summation over a *random set* of directions. This has the advantage that all directions are sampled to sufficient detail: if too few directions are included, subsequent realizations will give different estimates of $`J_\nu `$ (see § 3.5 for further discussion of this issue).

Originally Bernes (1979), the Monte Carlo approach was phrased in terms of randomly generated ‘photon packages’, which are followed as they travel through the source and which together approximate the radiation field. Fig. 2 illustrates that a formulation in terms of randomly chosen directions from each cell yields an equivalent estimate of $`J_\nu `$. The only difference is the direction of integration in Eq. (2). Where the former approach follows the photons as they propagate through the cells, the latter backtracks the paths of the photons incident on each cell. As the next section will discuss, this latter approach lends itself better to improvements in its convergence characteristics. Treatment of non-isotropic scattering is more complicated in this approach, and since scattering is not important at the wavelengths of interest here, $`10\mu `$m, scattering is not included in the code. Implementations of the Monte Carlo method more appropriate for scattering are available in the literature Wood et al. (1996a, b); Wolf et al. (1999).

### 3.4 Convergence and acceleration

Besides estimating $`J_\nu `$, an important aspect of non-LTE radiative transfer is convergence towards the correct solution in a reasonable amount of time. Since the solution is not a priori known, convergence is often gauged from the difference between subsequent iterative solutions. This relies on the assumption that when $`J_\nu `$ and the populations are far from their equilibrium solution, corrections in each iteration are large. Large optical depth can be a major obstacle to this behaviour: emission passing through an opaque cell will rapidly lose all memory of its incident intensity and quickly tend toward the local source function. The distance over which the information about changes in excitation can travel is one mean free path per iteration, so that the required number of iterations grows $`\tau ^2`$ characteristic of random walk. This effect makes it hard to determine if the process has converged.

Accelerated or Approximated Lambda Iteration (Rybicki and Hummer, 1991, ALI), circumvents this problem by defining an approximate operator $`\mathrm{\Lambda }^{}`$ such that

$$J_\nu =\left(\mathrm{\Lambda }\mathrm{\Lambda }^{}\right)[S_{ul}^{}(J_\nu )]+\mathrm{\Lambda }^{}[S_{ul}(J_\nu )].$$
(13)
An appropriate choice for $`\mathrm{\Lambda }^{}`$ is one which is easily invertible and which steers rapidly toward convergence. This occurs if $`J_\nu `$ is dominated by the second term on the right hand side of the equation, where $`\mathrm{\Lambda }^{}`$ works on the *current* source function as opposed to the solution from the previous iteration.

After several attempts Scharmer (1981), Olson, Auer, & Buchler Olson et al. (1986) found that a good choice for $`\mathrm{\Lambda }^{}`$ is the diagonal, or sometimes tri-diagonal, part of the full operator $`\mathrm{\Lambda }`$. This choice for $`\mathrm{\Lambda }^{}`$ describes the radiation field generated locally by the material in each cell, and its direct neighbours in the case of the tri-diagonal matrix. Eq. (13) then gives $`J_\nu `$ as the sum of the field incident on each cell due to the previous solution for the excitation {$`(\mathrm{\Lambda }\mathrm{\Lambda }^{})[S_{ul}^{}]`$}, and a *self-consistent* solution of the *local* excitation and radiation field {$`\mathrm{\Lambda }^{}[S_{ul}]`$}. In opaque cells, the radiation field is close to the local source function, and Eq. (13) converges significantly faster than Eq. (12); for optically thin cells, both formalisms converge equally fast.

Formulating the Monte Carlo method in terms of randomly generated photon packages traveling through the source does not easily permit separation of the locally generated field and the incident field for each cell. However, such a separation is possible when $`J_\nu `$ is constructed by summation over a set of rays, which each start at a random position within the cell and point in a random direction. For ray $`i`$, call the incident radiation on the cell $`I_{0,i}`$ and the distance to the boundary of the cell $`ds_i`$. The current level populations translate this distance into an opacity $`d\tau _i`$, and give the source function $`S_{ul}`$. The average radiation field from $`N`$ rays then follows from Eqs. (3) and (7),

$`J_\nu `$ $`=`$ $`{\displaystyle \frac{1}{N}}{\displaystyle \underset{i}{}}I_{0,i}\mathrm{e}^{\tau _i}+{\displaystyle \frac{1}{N}}{\displaystyle \underset{i}{}}S_{ul}[1\mathrm{e}^{\tau _i}]`$ (14)
$`=`$ $`J_\nu ^{\mathrm{external}}+J_\nu ^{\mathrm{local}}.`$ (15)
Here, $`S_\nu `$ and $`d\tau _i`$ contain both line and continuum terms, and $`I_{0,i}`$ includes the CMB. The radiation field is now the sum of the external ($`J_\nu ^{\mathrm{external}}`$) and internal ($`J_\nu ^{\mathrm{local}}`$) terms. Since the external term is evaluated using populations from the previous Monte Carlo iteration (through $`\tau _i`$ and $`S_{ul}`$), this scheme is akin to accelerated $`\mathrm{\Lambda }`$-iteration. Within Monte Carlo iterations, sub-iterations are used to find the locally self-consistent solution of $`S_{ul}`$ and $`\tau _i`$ for given $`J_\nu ^{\mathrm{external}}`$.

The main computational cost of this strategy lies in following a sufficient number of rays out of each cell through the source. Iteration on Eq. (14) is relatively cheap and speeds up convergence considerably in the presence of opaque cells. Fig. 3 illustrates this, by showing the evolution of the fractional error of the solution of the simple problem posed in §2 for optically thick HCO+ and thin H13CO+ excitation (for a fixed set of directions – see below).

Population inversions require careful treatment in radiative transfer codes, since the associated opacity is negative and the intensity grows exponentially. In general, an equilibrium will be reached where the increased radiation field quenches the maser. Since iterative schemes solve the radiative transfer before deriving a new solution for the excitation, the radiation field can grow too fast if population inversions are present. Our code handles negative opacities by limiting the intensity to a fixed maximum which is much larger than any realistic field. Proper treatment requires that the grid is well chosen, so that masing regions are subdivided into small cells where the radiation field remains finite. Our code can deal with the small population inversions that occur in many problems including the model presented in §2. However, to model astrophysical masers, specialized codes are required (Spaans and van Langevelde, 1992, e.g.).

### 3.5 The role of variance in Monte Carlo calculations

Because the Monte Carlo method estimates $`J_\nu `$ from a randomly chosen set of directions, the result has a variance, $`\sigma `$, which depends on the number $`N`$ of included directions as $`\sigma 1/\sqrt{N}`$. As explained above (§ 3.3), this variance is a strength rather than a weakness of the Monte Carlo method. Since it is not a priori known how many directions are required for a fiducial estimate of $`J_\nu `$, this method automatically arrives at an appropriate sampling by increasing $`N`$ until the variance drops below a predefined value.

The variance of a solution is usually estimated from the largest relative difference between subsequent iterations. In our implementation (see appendix), the number $`N`$ of rays making up $`J_\nu `$ in a particular cell is doubled each time the variance in that cell exceeds a certain value; the variance is evaluated using the largest relative difference between three subsequent solutions with the same $`N`$. This cell-specific scheme ensures that the radiation field is sufficiently sampled everywhere, and at the same time prevents oversampling of cells which are close to LTE and/or weakly coupled to other regions.

The variance as estimated from the difference between subsequent solutions only reflects the noise if random fluctuations dominate the difference. There will be systematic differences between subsequent solutions if these are still far from convergence. Therefore, many Monte Carlo methods consist of two stages. In the first stage, a fixed number of photons will yield a roughly converged solution; in the second stage, the number of photons is increased until the noise becomes sufficiently small.

In our implementation, this first stage consists of iterations with a fixed number of directions making up $`J_\nu `$ in each cell, $`N_0`$, which depends on the model. The directions are randomly distributed, but in each iteration, the *same* set of random directions is used by resetting the random number generator each iteration. Without random fluctuations in $`J_\nu `$, the difference between subsequent solutions purely reflects the evolution toward convergence. The first stage is considered converged when this ‘noise’ is a factor of ten smaller than the user-specified level.

For a sufficiently large $`N_0`$ (typically a few hundred), the excitation in each cell now is close to the final solution, except for imperfections in the sampling of $`J_\nu `$. In the second stage, each iteration uses a *different* set of random directions to estimate $`J_\nu `$: the random number generator is no longer reset. Based on the resulting variance, the number of rays in each cell is doubled each iteration, until the noise on the level populations in each cell is below a given value. If $`N_0`$ was initially insufficient, the variance will contain a significant contribution from systematic differences between iterations. Even though this will slow down the code by artificially increasing the number of rays in these cells as the code over-compensates the variance, ultimately the code will still converge to the correct solution.

The separation between local and incident radiation fields in our method (§3.4) keeps the system responsive to changes even in the presence of significant optical depth. This accelerates the convergence, but also increases the noise level. The literature mentions several methods to reduce the noise of Monte Carlo methods, e.g., with a reference field Bernes (1979); Choi et al. (1995) or quasi-random instead of pseudo-random numbers Juvela (1997). These schemes are useful when assumptions about the solution are available, but may slow down convergence if the initial guess is far off. Since the ‘first stage’ described above and the presence of noise prevents Monte Carlo methods from ‘false convergence’, we have not included any noise reduction techniques in our code, to keep it as widely applicable as possible.

### 3.6 Implementation and performance characteristics

Appendix A describes the structure of the program in detail, and provides a reference to a web site where we have made the source code of its spherically symmetric version publicly available. To test its performance, Appendix B compares results obtained with our code to those of other codes.

The main part of the program deals with calculating the excitation through the source model. In a separate part, comparison to observations proceeds by integrating Eq. (2) on a grid of lines of sight for a source at a specified inclination angle and distance. The resulting maps of the sky brightness may be convolved with a beam pattern for comparison with single-dish data, or Fourier transformed for comparison with interferometer data.

Appendix B describes tests of the validity of the results of the program. We have also tested up to what optical depth the program can be used, and found that this depends on source model. These tests were done on a Sun Ultrasparc 2 computer with 256 Mb internal memory and a 296 MHz processor. For a simple, homogeneous HCO+ model with $`n=10^4`$ cm-3 and $`T=30`$ K, the code produces accurate results within an hour for values of $`N`$(HCO+) up to $`10^{17}`$ cm-2, corresponding to $`\tau 20,000`$ in the lowest four rotational transitions. Higher$`J`$ lines are less optically thick under these physical conditions. For such opaque models, ‘local’ approximations fail badly, because the excitation drops sharply at the edge of the model (Bernes 1979; Appendix A).

For a power-law, Shu-type model, performance is somewhat slower. The dense and warm region fills only a small volume, while its radiation has a significant influence on the excitation further out, and modeling this effect requires a large number of rays. We have used the specific model from the Leiden workshop on radiative transfer (Appendix B) for various values of the HCO+ abundance. Within a few hours, accurate results are obtained for values of HCO+/H2 up to $`10^6`$, corresponding to $`\tau =6002000`$ in the lowest four lines.

These test cases should bracket the range of one-dimensional models of interest. For two-dimensional models, the limitations of present-day hardware are much more prohibitive. As a test case, we have used the flattened infalling envelope model from § 4.1 for various HCO+ abundances. Within 24 hours, the above machine provides a converged solution for HCO+/H2 up to $`10^8`$, corresponding to a maximum optical depth of $`100`$. Realistic models often have higher opacities, and call for the use of parallel computers. However, as faster computers are rapidly becoming available, we expect that these limitations will become less relevant in the near future. For both one- and two-dimensional models, the second, ray-tracing part of the code to create maps of the sky brightness takes only a fraction of the computer resources of the first part.

### 3.7 Alternative accelerators

Another method to tackle slow convergence in the presence of large opacities is *core saturation* Rybicki (1972); Doty and Neufeld (1997), where photons in the optically thick line center are replaced by the source function and no longer followed, while photons in the still optically thin line wings which carry most information are more accurately followed. This scheme has been implemented in a Monte Carlo program by Hartstein & Liseau Hartstein and Liseau (1998), but involves a choice where to separate the line core from the line wings. Since the effectiveness of the method depend on this choice, we have not included core saturation in our program.

A completely different approach to accelerating convergence is to extrapolate the behaviour of the populations from the last few iterations. This so-called Ng acceleration Ng (1974) is not implemented in our code, because extrapolating from an inherently noisy estimate may be dangerous.

## 4 Astrophysically relevant examples

A first example of the application of our code has already been given in §§2 and 3.4. This model is a spherically symmetric (one-dimensional) model; many authors have already illustrated the capability of Monte Carlo and other methods in solving one-dimensional problems (Bernes, 1979; Zhou, 1995; Choi et al., 1995, e.g.). This section presents a number of astrophysically relevant examples, again drawn from star formation studies, to illustrate the distinguishing properties of our code – in addition to accelerated convergence: the capability to calculate axisymmetric source models with a wide range of spatial scales, and the effects of dust continuum on radiative transfer and excitation.

The models presented in the following sections all include continuum radiation fields from dust. For these star-forming regions, we have chosen the model of Ossenkopf & Henning Ossenkopf and Henning (1994) for the dust emissivity, which includes grain growth for a period of $`10^5`$ yr at an ambient density of $`10^6`$ cm-3 and thin ice mantles. Except for the calculations in §4.3 where we specifically examine the effect of dust on the excitation including infrared transitions, only (sub) millimetre transitions were included in the excitation calculations which are not significantly influenced by the relatively weak continuum field.

### 4.1 A young stellar object with rotation

Observations of nearby young stellar objects often show flattened structures rather than the spherical symmetry of models like that of Shu Shu (1977), presumably caused by ordered magnetic fields and/or rotation (Hogerheijde et al., 1998, e.g.). These mechanisms probably influence the accretion behaviour, and may give rise to a rotationally supported circumstellar disk. To test these ideas against observations, cylindrically symmetric source models need to be considered. This section examines a model of protostellar collapse that includes flattening due to rotation following the description of Terebey, Shu, & Cassen Terebey et al. (1984), and its appearance in aperture synthesis maps.

The model of Terebey et al. Terebey et al. (1984) treats rotation as a small perturbation to the solution of Shu Shu (1977) for a collapsing envelope, joined smoothly to the description of a circumstellar disk by Cassen & Moosman Cassen and Moosman (1981). In addition to the sound speed and age, which we take identical to the values of §2 of $`a=0.24`$ km s-1 and $`t=1\times 10^5`$ yr, this model is parameterized by a rotation rate $`\mathrm{\Omega }`$. This gives rise to a centrifugal radius $`R_cam_0^3t^3\mathrm{\Omega }^2/16`$, within which the infalling material accretes onto a thin disk. Here, $`m_0=0.975`$ is a numerical constant. We choose $`\mathrm{\Omega }=5.9\times 10^{13}`$ s-1, so that $`R_c=800`$ AU. We assume that inside $`R_c`$ the material accretes onto a thin disk, and that most molecules rapidly freeze out onto dust grains (cf. §4.2). Effectively, we assume the region within $`R_c`$ to be empty for this calculation. Fig. 4 (top left) shows the adopted density structure. All other parameters are similar to the model of §2.

To follow the power-law behaviour of the density in the model, a total of $`15\times 15`$ cells are spaced exponentially in the $`R`$ and $`z`$ directions. To reach a final signal-to-noise ratio of 10, with 300 rays initially making up the radiation field in each of the cells, the Monte Carlo code requires 5 hours CPU time on a UltraSparc 10 to converge on the HCO+ solution. For comparison, the optically thin and more readily thermalized 13CO problem takes only 10 minutes. The resulting excitation temperature distribution (Fig. 4; middle panels) does not deviate much from that of the spherically symmetric model of §2, apart from the flattened distribution of the material at the center: rotation is only a small perturbation on the overall source structure. As a result, the appearance is mostly unaffected in single-dish observations which do not resolve scales comparable to $`R_c`$ where flattening is important.

Millimetre interferometers, on the other hand, can resolve these scales at the distances of the nearest star-forming regions, $`140`$ pc. Fig. 4 (right panels) shows the integrated emission in HCO+ $`J`$=1–0 and 3–2 after sampling at the same spatial scales as real interferometer observations and subsequent image reconstruction. Delays between 100–1000 ns were used, corresponding to angular scales of $`2^{\prime \prime }`$$`40^{\prime \prime }`$ for the 1–0 line and $`0\stackrel{}{.}6`$$`13^{\prime \prime }`$ for 3–2. Hence, emission on scales $`6000`$ AU ($`2000`$ AU at 3–2) is filtered out. This results in a reconstructed (‘cleaned’) image that is dominated by the central, flattened regions of the envelope when the object is seen edge-on. Aperture synthesis observations of embedded protostars in Taurus show similar structures (Ohashi et al., 1997; Hogerheijde et al., 1998, e.g.).

### 4.2 A circumstellar disk

Planetary systems form from the disks that surround many young stars Beckwith (1999); Mannings et al. (2000). Observational characterization of these disks is of prime importance to increase our understanding of the processes that shape planetary systems. Here, we present simulations of molecular line observations of a circumstellar disk around a T Tauri star as obtained with current and planned millimetre-interferometric facilities.

The physical structure of the model disk is that of a passive accretion disk in vertical hydrostatic equilibrium as described by Chiang & Goldreich Chiang and Goldreich (1997). This description includes the effect of ‘backwarming’ of the disk by thermal radiation of a thin, flared surface layer that intercepts the stellar light. The total amount of material in the superheated surface layer is too small to be detectable, but the overall effect of increased mid-plane temperature is significant. The effective temperature of the central star is 4000 K and its luminosity is 1.5 L. The outer radius of the disk is 700 AU, with a total mass of 0.04 M.

The chemical structure of the disk follows Aikawa & Herbst Aikawa and Herbst (1999), who describe the radial and vertical composition of a flared disk. Freezing out of molecules onto dust grains is one of the dominant processes influencing the gas-phase composition in disks, and strongly depends on temperature and density. In the dense and cold midplane, many molecules will be depleted onto grains. However, close to the star where temperatures are raised, and away from the midplane where densities are lower and depletion time scales longer, gas-phase abundances will be significant. In addition, ultraviolet radiation and X-rays may penetrate the upper layers of the disk, photodissociating molecules and increasing the abundance of dissociation products like CN. Fig. 5 shows the distribution of the density, temperature, and abundances of CO, HCO+, HCN and CN in the adopted model. We have used the results presented in Aikawa & Herbst (Aikawa and Herbst, 1999, their Figs. 6 and 7; high ionization case), and parameterized the abundances as function of scale height.

For the Monte Carlo calculations, a gridding is adopted that follows the radial power-law density profile in 14 exponentially distributed cells and the vertical Gaussian profile in 13 cells linearly distributed over 3 scale heights. Convergence to a signal-to-noise ratio of 10 requires approximately 6 hours CPU time per model on an UltraSparc 10 workstation, starting with 100 rays per cell and limiting the spatial dynamic range to 36 (i.e., the smallest cell measures 10 AU on the side). The resulting excitation and emission depends on the competing effects of increased abundance and decreased density with distance from the midplane. Fig. 6 shows the excitation temperature of selected transitions and molecules. Fig. 7 shows a number of representative simulated observations, at resolutions of $`2^{\prime \prime }`$, $`0\stackrel{}{.}5`$, and $`0\stackrel{}{.}2`$ attainable with current and planned (sub) millimetre interferometers. Van Zadelhoff et al. (in prep.) present a simpler analysis of similar models.

### 4.3 A high mass young stellar object

The above examples referred to the formation of stars with masses of up to few M and luminosities $`100`$ L. Stars of higher mass spend their first $`10^4`$ yr in envelopes of $`100`$ M. With their luminosities of $`10^4`$$`10^5`$ L, these stars heat significant parts of their envelopes to several hundred K, shifting the peak of the Planck function to the wavelengths of the vibrational transitions of many molecules. Stars of lower mass and luminosity only heat small regions to a few hundred K, and the impact on the excitation is correspondingly smaller. As an example, Figure 8 shows two models of the HCN submillimetre line emission, with and without pumping through the bending (“$`\nu _2`$”) mode at $`14.02`$ $`\mu `$m. For computational convenience, only energy levels up to $`J`$=10 within the first vibrationally excited and ground states have been included, including l-type doubling in the excited state. This doubling occurs due to the two possible orientations of the rotational and vibrational motions with respect to each other. Collisional rate coefficients between rotational levels are from Green (1994, priv. comm., see http://www.giss.nasa.gov/data/mcrates); within vibrational levels, they were set to $`10^{12}`$ cm3s-1. The source model is that of the young high-mass star GL 2136 by van der Tak et al. van der Tak et al. (2000). Based on its luminosity of $`7\times 10^4`$ L and dust mass of $`100`$ M, the star has heated a region of radius $`3000`$ AU to $`100`$ K, making pumping through the $`14`$ $`\mu `$m HCN bending mode important. We have assumed that $`T_{\mathrm{gas}}=T_{\mathrm{dust}}`$, as is true for high density regions. The dust emissivity, from Ossenkopf & Henning Ossenkopf and Henning (1994), is the same as in the previous sections. As seen in Fig. 8, the effect of dust is especially strong for the rotational lines within the $`\nu _2`$=1 band, which occur at frequencies slightly offset from the ground state transitions. These lines have indeed been detected towards similar objects (Ziurys and Turner1986, , e.g.).

The shells around evolved stars is another area where inclusion of infrared pumping by dust is essential to understand the rotational line emission (Ryde et al., 1999, e.g.). Many molecules that are commonly observed through rotational lines at millimetre wavelengths have ro-vibrational bands in the mid-infrared, and can be pumped by warm dust. In a few cases, pumping through *rotational* lines at far-infrared wavelengths is important as well, for example CS Hauschildt et al. (1993) and all hydrides, most notably water Hartstein and Liseau (1998).

## 5 Conclusion

This paper describes a computer code to solve the coupled problem of radiative transfer and molecular excitation for spherically and cylindrically symmetric source geometries. It is based on the Monte Carlo method, but incorporates elements from accelerated lambda iteration which greatly improve convergence in the presence of significant optical depth. In particular, the code separates excitation due to the local environment from the response to the radiation field after propagation through the source. This approach combines the flexibility of a Monte Carlo method with the reliability of accelerated lambda iteration. We expect our code to be a valuable tool in the interpretation of the wealth of data that current and future instruments operating from the millimetre to the infrared will yield. A number of examples (§4) already illustrates the applications to problems in star formation studies.

###### Acknowledgements.

The authors wish to thank Ewine van Dishoeck for many stimulating discussion on the topic of this paper, her careful reading of the manuscript, and making available a data base with molecular parameters for this code; David Jansen for maintaining this data base; Lee Mundy for assistance with the construction of an earlier version of the part of the code calculating the sky brightness distributions; Marco Spaans for discussions and hospitality at the Astronomy Department of the Johns Hopkins University where an initial version of the code was written; Minho Choi and Neal Evans for the use of their Monte Carlo program for testing purposes; the organizers (Gerd-Jan van Zadelhoff, Kees Dullemond and Jeremy Yates) and participants of the May 1999 workshop on Radiative Transfer in Molecular Lines at the Lorentz Center of Leiden University; and George Rybicki for suggesting that iterating on the level populations may be of benefit in a Monte Carlo code. The research of M. R. H. is supported by the Miller Institute for Basic Research in Science; that of FvdT by NWO grant 614-41-003.

## Appendix A The code

Our code is applicable to a wide range of astrophysical problems involving (sub) millimetre and infrared ($`10`$ $`\mu `$m) spectral-line and continuum observations. The one-dimensional (spherically symmetrical) version of our code is publicly available for all interested researchers via http://astro.berkeley.edu/$``$michiel. The two-dimensional (cylindrically symmetric) code is available on a collaborative basis through the authors (see the same website for contact information). This appendix gives a concise description of the implementation of the accelerated Monte Carlo method.

Fig. 9 gives an overview of the structure of our code, which consists of two parts. The first runs the Monte Carlo simulation solving the radiative transfer and molecular excitation. The second part uses this solution to calculate the emission that would be observed from this source above the atmosphere and with perfect spatial and velocity resolution, given a source distance and, for cylindrically symmetric models, inclination. This latter part can also be used to calculate the continuum radiation emitted by the source. Its output format is that used by the MIRIAD package \[Multichannel Image Reconstruction, Image Analysis, and Display; Sault et al. Sault et al. (1995)\]. This package, designed to analyse interferometer interferometric spectral line data, includes many processing options such as convolution with a single-dish beam and modeling of aperture-synthesis visibilities, as well as a wide variety of imaging capabilities. MIRIAD also allows easy conversion to the ubiquitous FITS format and portability to other software packages.

Both parts of the code are controlled by UNIX C-shell scripts that extract information from the provided input and compile an executable code. In this way, the size of several arrays containing the source model, the collisional rate coefficients, etc., can be adjusted to the required size, minimizing memory requirements. The source code is written in FORTRAN-77.

Following the flow chart of Fig. 9, the following steps describe the Monte Carlo part of the code in more detail.

1. The code starts by reading a list of keywords, detailing the required signal-to-noise ratio on the level populations, the initial number of photons in each cell ($`N_0`$), and pointers to the source model, the systematic velocity field (if any), the description of the dust emissivity, and the molecular energy levels and collisional rate coefficients. The velocity field can be defined simply through the source model with each grid cell moving at a constant speed, or it can be a constantly varying function over each cell. The source model can be a series of concentric shells covering a region from the origin to a maximum radius, or a series of stacked cylinders fully covering a region out to a maximum radius and height.

Collisional rate coefficients are available for many astrophysically interesting species and common collision partners such as H2 in the $`J=0`$ and in the $`J=1`$ levels, e-, and He. Our code currently allows for two simultaneous collision partners, e.g., H2 and e-, each with its own density and temperature. For molecular ions such as HCO+, excitation due to collisions with electrons can be significant compared to collisions with H2 at fractional ionization ($`10^5`$). Often, listed rate coefficients are equivalent rates per H2 molecule including contributions from He at cosmic abundance. The results of our code, and any non-LTE calculation, sensitively depend on the quality of the rate coefficients. Recently, Black Black (2000) discussed the need for good rate coefficients and the effects of other implicit assumptions of radiative transfer codes.
2. In the first stage of the calculation, the radiation field is based on $`N_0`$ rays per cell, each starting at a random position equally distributed over the cell volume, pointing in a random direction, and at a random frequency within 4.3 times the local line width around the local systematic velocity vector. The value of 4.3 corresponds to the width where the line profile has dropped to less than 1% of its peak. In this stage, in each iteration the *same* series of random numbers is used, so that there are no random fluctuations in the coverage of the radiation field.
3. For each ray, the distance $`ds`$ from the ray’s origin to the nearest boundary of the cell along its randomly chosen direction is calculated. The incident radiation $`I_0`$ along the ray then follows from integrating Eq. (2) in a stepwise manner from cell edge to cell edge, attenuating the contribution from each cell by all intervening cells, with the cosmic microwave background as a boundary condition. The only quantity that changes when stepping through a cell is the direction, and possibly the magnitude, of the systematic velocity vector, which enters Eq. 2 through the line profile function $`\varphi (\nu )`$. Changes of $`\varphi (\nu )`$ within cells are tracked by subdividing the integration into small steps as needed.
4. Armed with the set of $`(ds,\varphi (\nu ),I_0)`$ for each ray, the radiation field $`J_\nu `$ in the cell follows from Eq. (14). A consistent solution of this equation and the level populations (Eq. 10) quickly follows from iteration to a relative accuracy of $`10^6`$ in the populations. Limitations on masering are discussed in Sect. 3.4.
5. The first stage of the code repeats items 2–4 until the largest relative fractional difference between the populations in all cells of three subsequent solutions is ten times better than ultimately required. Since the angular sampling is the same in each iteration, these differences are free of random noise but might not adequately sample all directions and frequencies.
6. The second stage of the code proceeds along similar lines as the first stage, but with a *different* set of random numbers in each iteration. The only other difference is, that each time the maximum fractional error in the populations in a cell exceeds the requested accuracy, the number of rays in that cell is doubled. This stage lasts until all cells comply with the required accuracy, after which the solutions are written out to a file.

The second part of the program calculates the emission distribution on the sky for a given source distance and inclination by simple ray tracing. The output from the Monte Carlo code forms the input for this ray-tracing code. Since it uses much of the same code as the Monte Carlo part, geometry and radiative transfer being the same, it is not further discussed here.

## Appendix B Comparison with other codes

This section describes two cases to test our code against well-documented calculations with Monte Carlo codes from the literature. For further tests, we refer the reader to the web-page collecting a number of standard test cases, which has resulted from the 1999 workshop on Radiative Transfer in Molecular Lines at the Lorentz Center of Leiden University (http://www.strw.leidenuniv.nl/$``$radtrans).

### B.1 Bernes’ CO cloud

In his seminal paper on Monte Carlo methods for radiative transfer and molecular excitation, Bernes Bernes (1979) presents a constant-density, constant-temperature, optically thick cloud model. The density of the cloud, $`n_{\mathrm{H}_2}=2000`$ cm-3, is below the critical density of the CO transitions, and the excitation is dominated by radiative trapping. The excitation temperatures of the CO transitions drop off rapidly in the outer regions of the cloud. This necessitates fine sampling of these regions. Figure 10 shows that our code reproduces the original results within the accuracy of our and Bernes’ calculations. This simple model forms a critical test for the code’s ability to correctly handle excitation by radiative trapping. The total run time for the model was approximately 5 minutes on a UltraSparc 10 workstation, using the same collisional rate coefficients as Bernes (Table 1).

### B.2 Model for B 335 by Choi et al.

Another critical element of any radiative transfer code is its ability to correctly deal with systematic velocity fields. The inside-out collapse model as outlined in §2 is well suited for such a test, because of its wide range in velocities from zero to many times the turbulent line width combined with significant optical depth. As a test case, we calculate the populations and the emergent spectrum of several CS and H2CO lines, following the model for the young stellar object B 335 of Choi et al. Choi et al. (1995). This model is similar to that of § 2, with $`a=0.23`$ km s-1 and $`t=1.3\times 10^5`$ yr. The turbulent line width is 0.12 km s-1. Only the temperature structure is different from § 2: Choi et al. use continuum observations to constrain the temperature distribution, which we follow as closely as possible from their Fig. 3.

Fig. 11 compares the resulting excitation temperature distribution with the results of Choi et al. The agreement is very good for CS, where we used the same collisional rate coefficients (Table 1). For H2CO the agreement is less favourable, but we were unable to use the exact same rate coefficients. Simple tests show that the deviation is comparable to what can be expected from the difference in the molecular data. This variation corresponds to a 10% difference in the emergent line profiles.

Fig. 11 also plots the excitation temperatures obtained for the same model but using the cylindrically symmetric code. Both codes clearly give consistent answers; the small ‘wiggles’ in the excitation temperatures as function of radius in the output of the cylindrically symmetric calculation can be attributed to geometrical defects when trying to fit a sphere in a series of stacked cylinders.


format: text

text:
# Casimir energy for a scalar field with a frequency dependent boundary condition

## I Introduction

The Casimir effect  arises as a distortion of the vacuum energy of quantized fields due to the presence of boundaries (or nontrivial topologies) in the quantization domain. This effect, which has a quantum nature associated with the zero-point oscillations in the vacuum state, is significant in diverse areas of physics, from statistical physics to elementary particle physics and cosmology.

In particular, in the last years there has been a great interest in the Casimir energy of electromagnetic fields in the presence of dielectric media, due to Schwinger’s suggestion  that it could play a role in the explanation of the phenomenon of sonoluminescence .

The results obtained on this subject by different groups through several calculation techniques (as Green’s functions methods, van der Waals forces, $`\zeta `$-function methods and asymptotic developments for the density of states - see references  among others) are rather controversial, and some basic issues remain to be clarified. In this respect, it is our aim to contribute to the understanding of the problem by studying a simplified model, which incorporates a frequency cut-off in the boundary conditions at the separation between media, to emulate the behavior of real dielectrics.

It should be mentioned that authors introducing a cut-of in the wave number to describe the behavior of real dielectrics agree with Schwinger’s explanation of the phenomenon (See, for example, ). On the other hand in Ref. , where the permittivity of a non-magnetic medium is modeled by means of a (angular momentum dependent) one-absorption-frequency Sellmeir relation, and a $`\zeta `$-function technique is employed to sum-up the contributions of the proper oscillation modes of the electromagnetic field, it is concluded that the vacuum energy of a bubble embedded in this material has the wrong behavior with the radius (besides its absolute value being far too small) to be relevant to sonoluminescence.

Candelas  has been the first to remark the importance of dispersion in connection with the Casimir effect. Later, Brevik and collaborators  have considered several models with spherical and cylindrical geometries, showing that dispersive effects may reverse the direction of the Casimir force acting on the boundary between media. In Ref.  it is studied the case of spherically symmetric media satisfying the relation $`ϵ(\omega )\mu (\omega )=1`$ (of special interest in QCD), both with a sharp frequency cut-off on $`\mu (\omega )`$ and with a behavior analogous to the one-absorption-frequency Sellmeir’s formula for dielectrics. In Ref. , for non-magnetic dispersive media with spherical geometry, and based upon the Minkowski energy-momentum tensor, the Casimir surface force is worked out when the permittivity $`ϵ(\omega )`$ presents a step along the imaginary frequency axis. In the limit of perfect conductivity, a non-dispersive term is recovered which is in agreement with Boyer result , while an attractive dispersive contribution is also found.

In what follows we consider a simple model of a scalar field subject to frequency dependent local boundary conditions on the surface of a sphere. Thus our main goal is to establish a method for calculating the change of the Casimir energy of the field when the radius of the sphere is varied, in a situation where the boundary conditions impose a physical frequency cut-off $`\mathrm{\Omega }`$.

To this end we consider the very simple case of a scalar field whose modes corresponding to eigenfrequencies $`\omega \mathrm{\Omega }`$ are confined to the interior of a sphere of radius $`R`$, satisfying local homogeneous boundary conditions.

On the other hand, we will assume that the boundary is completely transparent for those modes with $`\omega >\mathrm{\Omega }`$. Therefore, their contribution to the difference of Casimir energies for two different values of $`R`$ will cancel out, no matter the regularization employed for its definition. Consequently, we will subtract these contributions, which amounts to a redefinition of the reference energy level in an $`R`$-independent way.

For the evaluation of the vacuum energy of the low frequency modes we will employ asymptotic expansions in an incomplete-$`\zeta `$ summation technique, to be discussed in the following. This approach will allow for the identification of the volume, surface, …, dominant terms in the Casimir energy.

As a final exercise we will show that, with reasonable values for the frequency cut-off and the bubble radius, the amount of Casimir energy available in this model is comparable with the energy emitted during each cycle of a typical sonoluminesce experiment . In view of the controversy existing on the subject, it would be of great interest to apply this method to a similar model for the case of the electromagnetic field in the presence of dielectric media, calculation which will be presented in a forthcoming paper .

## II The model and its incomplete $`\zeta `$-function

Let us consider a free scalar field in $`𝐑^3`$ satisfying at the surface of a sphere of radius $`R`$, local boundary conditions which depend on the frequency $`\omega `$ of the field modes.

We will make the assumption that the boundary is completely transparent for the modes of frequencies grater than a cut-off $`\mathrm{\Omega }`$, while for $`\omega \mathrm{\Omega }`$ the modes satisfy Dirichlet boundary conditions,

$$\begin{array}{c}\left(\mathrm{}+\frac{\omega ^2}{c^2}\right)\psi _\omega (\stackrel{}{r})=0,\mathrm{for}r<R,\\ \\ \psi _\omega (\stackrel{}{r})|_{r=R}=0,\end{array}$$
(1)
being confined to the interior of the sphere.

Writing $`\psi _\omega (\stackrel{}{r})=f_l(r)Y_l^m(\theta ,\phi )`$, we get for the radial function

$$\left(\frac{d^2}{dr^2}+\frac{2}{r}\frac{d}{dr}\frac{l(l+1)}{r^2}+\frac{\omega ^2}{c^2}\right)f_l(r)=0,\mathrm{for}r<R,$$
(2)
where the eigenfrequencies are determined by imposing the condition

$$f_l(r)|_{r=R}=0.$$
(3)
The solutions regular at the origin are given by $`f_l(r)=\sqrt{{\displaystyle \frac{\pi }{2z}}}J_\nu (z)`$, with $`\nu =l+1/2`$ and $`z=\omega _{\nu ,n}r/c`$, where the eigenfrequencies are

$$\omega _{\nu ,n}=\frac{c}{R}j_{\nu ,n},$$
(4)
being $`j_{\nu ,n}`$ the $`n`$-th zero of the Bessel function $`J_\nu (z)`$.

We will be interested in the difference between the vacuum energies of two situations differing in the value of $`R`$. Then, we can disregard the contributions of those modes with frequencies $`\omega >\mathrm{\Omega }`$ because, being independent of the position of the boundary, their contributions cancel out (whatever the regularization employed in defining the vacuum energy would be). This simply amounts to an $`R`$-independent subtraction, which is nothing but a redefinition of the zero energy level.

Therefore, we should evaluate the (finite) sum

$$E(R)=\underset{\nu =1/2}{\overset{\nu _0}{}}2\nu \underset{n=1}{\overset{N_\nu }{}}\frac{1}{2}\mathrm{}\omega _{\nu ,n},$$
(5)
where $`N_\nu `$ is the number of positive zeroes of $`J_\nu (z)`$ less than or equal to $`x=\mathrm{\Omega }R/c`$, the factor $`2\nu =2l+1`$ is the eigenvalue degeneracy, and $`\nu _0`$ is the maximum value of $`\nu `$ for which $`N_\nu 1`$.

We are interested in an analytic, rather than numeric, evaluation of eq. (5). Although this is a finite sum, we will employ a summation method based on the evaluation of an incomplete $`\zeta `$-function, an approach which could be applied in more complex situations. We can employ the following representation:

$$\underset{n=1}{\overset{N_\nu }{}}j_{\nu ,n}=\underset{n=1}{\overset{N_\nu }{}}j_{\nu ,n}^s|_{s=1},$$
(6)
where the sum in the right hand side obviously exists for any $`s𝐂`$ ***Notice that the sum in the r.h.s. of eq. (6) evaluated at $`s=0`$ gives $`N_\nu `$, the number of eigenfrequencies contributing to the Casimir energy of the field (after the subtraction made to define it) for a given value of the angular momentum $`l=\nu 1/2`$..

Since $`J_\nu (z)`$, for $`\nu >1`$, has only real zeros, and its non-vanishing zeros are all simple , we can employ the Cauchy theorem to represent the sum in the r.h.s. of (6) as an integral on the complex plane,

$$\underset{n=1}{\overset{N_\nu }{}}j_{\nu ,n}^s=\frac{1}{2\pi ı}_Cz^s\frac{J_\nu ^{}(z)}{J_\nu (z)}𝑑z,$$
(7)
where the curve $`C`$ encircles counterclockwise the $`N_\nu `$ first positive zeros of $`J_\nu (z)`$.

For $`\mathrm{}(s)`$ large enough, the contour $`C`$ can be deformed into two straight vertical lines, one crossing the horizontal axis at $`\mathrm{}(z)=x`$ and the other at $`\mathrm{}(z)=0^+`$. Indeed, expressing the integrand in terms of the modified Bessel function

$$I_\nu \left(w\right)=e^{ı\frac{\pi }{2}\nu }J_\nu \left(e^{i\frac{\pi }{2}}w\right)$$
(8)
(valid for $`\pi <\mathrm{arg}(w)\pi /2`$), and taking into account its asymptotic behavior for large arguments , it is easily seen that, for $`0<xj_{\nu ,n},n`$, the integral

$$\zeta _\nu (s,x)\frac{1}{2\pi ı}_{xı\mathrm{}}^{x+ı\mathrm{}}z^s\frac{J_\nu ^{}(z)}{J_\nu (z)}𝑑z,$$
(9)
converges absolutely and uniformly in the open half-line $`s>1`$, from which it can be meromorphically extended to the whole complex $`s`$-plane.

Therefore, for $`s>1`$,

$$\underset{n=1}{\overset{N_\nu }{}}j_{\nu ,n}^s=\zeta _\nu (s,0^+)\zeta _\nu (s,x).$$
(10)
And, since the left hand side of (10) is holomorphic in $`s`$, the singularities of $`\zeta _\nu (s,x)`$ must be independent of $`x`$.

On the other hand, for $`y>0`$ ,

$$\begin{array}{c}I_\nu \left(yix\right)=e^{ı\pi \nu }I_\nu \left(y+ix\right)\\ \\ I_\nu \left(y+ix\right)=\left(I_\nu \left(yix\right)\right)^{}.\end{array}$$
(11)
So, for real $`s>1`$ we can write

$`\zeta _\nu (s,x)=`$

$$\mathrm{}\left\{\frac{x^{1s}}{\pi }e^{i{\displaystyle \frac{\pi }{2}}\left(s+1\right)}_0^{\mathrm{}}(yı)^s\frac{I_\nu ^{}\left(x(yı)\right)}{I_\nu \left(x(yı)\right)}𝑑y\right\}.$$
(12)
In order to construct the analytic extension of $`\zeta _\nu (s,x)`$ to $`s1`$, we subtract and add to the integrand in (12) the first few terms obtained from the uniform asymptotic (Debye) expansion  of the Bessel functions,

$$\frac{I_\nu ^{}\left(\nu t\right)}{I_\nu \left(\nu t\right)}=\frac{1}{\nu }D_\nu (t)+𝒪(\nu ^3),$$
(13)
where

$$\begin{array}{c}D_\nu (t)=\nu D^{(1)}(t)+D^{(0)}(t)+\nu ^1D^{(1)}(t)=\\ \\ \frac{\nu \sqrt{1+t^2}}{t}\frac{t}{2\left(1+t^2\right)}+\frac{4tt^3}{8\nu \left(1+t^2\right)^{\frac{5}{2}}}.\end{array}$$
(14)
valid for large $`\nu `$ with fixed $`t`$. We will see that this approximation allows for the identification of the volume, surface, … contributions to the vacuum energy.

Changing the integration variable in eq. (12) to $`tz(yı)`$, with $`z=x/\nu >0`$, we get

$`\zeta _\nu (s,x)=`$

$$\mathrm{}\left\{\frac{\nu ^s}{\pi }e^{i{\displaystyle \frac{\pi }{2}}\left(s+1\right)}_{iz}^{\mathrm{}iz}t^s\frac{d\left(\mathrm{ln}I_\nu \left(\nu t\right)\right)}{dt}𝑑t\right\}.$$
(15)
So, we must consider the integral

$$\begin{array}{c}_{iz}^{\mathrm{}iz}t^s\frac{d\left(\mathrm{ln}I_\nu \left(\nu t\right)\right)}{dt}𝑑t=_{iz}^{\mathrm{}iz}t^sD_\nu (t)𝑑t+\\ \\ +_{iz}^{\mathrm{}iz}t^s\left\{\frac{d\left(\mathrm{ln}I_\nu \left(\nu t\right)\right)}{dt}D_\nu (t)\right\}𝑑t.\end{array}$$
(16)
The second integral in the right hand side of this equation converges for $`s>2`$, since the integrand can be estimated by means of the next ($`𝒪(\nu ^3)`$) term in the Debye expansion (eq. (13)), which behaves as $`𝒪(t^3)`$ for large $`|t|`$. It can be numerically evaluated at $`s=1`$. This will not be done in this paper.

In the following, we will consider only the first integral in the right hand side of (16), retaining the first few terms of its expansion in powers of $`\nu ^1`$ consistent with the approximation made in eq. (13).

Notice that the integrand is an algebraic function, having singularities at $`t=0,\pm ı`$, and behaving as $`𝒪(t^0)`$ for large $`|t|`$. This integral converges absolutely and uniformly for $`s>1`$, where it defines an analytic function to be meromophically extended to the region of interest of the parameter $`s`$. As we will see, this extension reveals the singularities of $`\zeta (s,x)`$ as simple poles, whose residues are independent of $`x`$ (a necessary condition to give a finite result in (10) for any $`s`$).

In fact, by virtue of the analyticity of the integrand, we can change the path of integration to write

$`{\displaystyle _{iz}^{\mathrm{}iz}}t^sD_\nu (t)𝑑t=`$

$$\begin{array}{c}_{iz}^1t^s\left(\frac{\nu \sqrt{1+t^2}}{t}\frac{t}{2\left(1+t^2\right)}+\frac{4tt^3}{8\nu \left(1+t^2\right)^{\frac{5}{2}}}\right)𝑑t\\ \\ +_1^{\mathrm{}}t^s\{\frac{\nu \sqrt{1+t^2}}{t}\frac{t}{2\left(1+t^2\right)}+\frac{4tt^3}{8\nu \left(1+t^2\right)^{\frac{5}{2}}}\\ \\ [\nu (1+\frac{1}{2t^2})\frac{1}{2t}\frac{1}{8\nu t^2}]\}dt+\\ \\ +_1^{\mathrm{}}t^s\left[\nu \left(1+\frac{1}{2t^2}\right)\frac{1}{2t}\frac{1}{8\nu t^2}\right]𝑑t.\end{array}$$
(17)
The first integral in the r.h.s. of eq. (17), containing the whole dependence on $`x=\nu z`$, is holomorphic in $`s`$ and can be directly evaluated at the required value of this parameter. On the half-line $`(1,\mathrm{})`$ we have subtracted and added the first terms in the series expansion of $`D_\nu (t)`$ for large $`t`$, thus making the second integral to converge for $`s>2`$. The third one must be evaluated for $`s>1`$ and then analytically continued to the relevant values of $`s`$. This can be exactly done, its contribution to $`\zeta _\nu (s,x)`$ in eq. (15) being the real part of

$$\begin{array}{c}\frac{e^{{\displaystyle \frac{i}{2}}\pi \left(1+s\right)}}{\nu ^{1+s}\mathrm{\hspace{0.17em}8}\pi }\left(\frac{8\nu ^2}{1s}+\frac{4\nu }{s}+\frac{14\nu ^2}{1+s}\right)\end{array}$$
(18)
This expression has simple poles at $`s=0,\pm 1`$, which are the only singularities of $`\zeta _\nu (s,x)`$ for $`\mathrm{}(s)>2`$. Notice that the residues of $`\zeta _\nu (s,x)`$ are independent of $`x`$,

$$\begin{array}{c}Res\zeta _\nu (s,x)|_{s=1}=\frac{1}{\pi },\\ \\ Res\zeta _\nu (s,x)|_{s=0}=0,\\ \\ Res\zeta _\nu (s,x)|_{s=1}=\frac{14\nu ^2}{8\pi },\end{array}$$
(19)
and in agreement with the results in  (where $`\zeta _\nu (s,0^+)`$ is studied).

For example, for $`\zeta _\nu (s,x)`$ around $`s=1`$ and for $`\nu <x`$ (which will be needed in Section IV to evaluate the vacuum energy), one straightforwardly obtains the Laurent expansion

$$\begin{array}{c}\zeta _\nu (s,x>\nu )=\frac{14\nu ^2}{8\pi \left(1+s\right)}+\\ \\ +[\frac{\left(4\nu ^21\right)}{8\pi }[\mathrm{log}\left(\frac{\nu }{2}\right)+\mathrm{log}(z+\sqrt{z^21})]\frac{\nu ^2}{4\pi }+\\ \\ \frac{\nu ^2}{2\pi }z\sqrt{z^21}\frac{3z8z^3}{24\pi \left(z^21\right)^{3/2}}\frac{1}{3\pi }+𝒪(\nu ^1)]+\\ \\ +𝒪(s+1),\end{array}$$
(20)
with fixed $`z=x/\nu 1`$.

On the other hand, for $`x0^+`$ a similar calculation leads to

$$\begin{array}{c}\zeta _\nu (s,x=0^+)=\frac{14\nu ^2}{8\pi \left(1+s\right)}+[\frac{\left(4\nu ^21\right)}{8\pi }\mathrm{log}\left(\frac{\nu }{2}\right)\\ \\ \frac{\nu ^2}{4\pi }\frac{\nu }{4}\frac{1}{3\pi }+𝒪(\nu ^1)]+𝒪(s+1).\end{array}$$
(21)
In the following Section we will evaluate, as a function of $`\nu `$, the number of modes contributing in eq. (10), and in Section IV, their contributions to the vacuum energy.

## III The number of contributing modes

In this Section we address ourselves to the determination of $`\nu _0`$ in (5), the maximum value of $`\nu `$ for which $`N_\nu 1`$. Although in the simple case under study the zeros of $`J_\nu (w)`$ are well known , we prefer to establish a criterium which can be applied in more general situations.

First, notice that

$$N_\nu (x)\underset{n=1}{\overset{N_\nu }{}}j_{\nu ,n}^s|_{s=0}=\left[\zeta _\nu (s,0^+)\zeta _\nu (s,x)\right]|_{s=0}$$
(22)
is a discontinuous function of $`x`$, having a step of height $`1`$ at each positive zero $`j_{\nu ,n}`$ of the Bessel function $`J_\nu (w)`$.

Then, $`\nu _0`$ can be determined from the condition

$$N_{\nu _0}(x)=N_{\nu _0}(j_{\nu _0,1}+0)=1,$$
(23)
with $`N_{\nu _0}(j_{\nu _0,1}0)=0`$.

Taking into account eq. (19) and the fact that the second and third integrals in the r.h.s. of eq. (17) are real, it is straightforward to obtain from eqs. (15-17) that

$$\begin{array}{c}\zeta _\nu (s=0,x)=\frac{\nu }{2}\frac{1}{4}+\\ \\ +\mathrm{}\{\frac{ı\nu }{\pi }(\sqrt{1+e^{ı\pi }z^2}\mathrm{log}(1+\sqrt{1+e^{ı\pi }z^2}))\\ \\ +\frac{ı}{4\pi }\mathrm{log}(1+e^{ı\pi }z^2)+\\ \\ +\frac{ı(2+3z^2)}{24\nu \pi \left(1+e^{ı\pi }z^2\right)^{3/2}}\}+𝒪(\nu ^2),\end{array}$$
(24)
where we have taken $`z=x/\nu 1`$. In particular, for $`x0^+`$,

$$\zeta _\nu (s=0,x=0^+)=\frac{\nu }{2}\frac{1}{4}+𝒪(\nu ^2),$$
(25)
in coincidence with .

Now, taking the difference between (25) and (24) we get a smooth approximation, $`\stackrel{~}{N}_\nu (x)+𝒪(\nu ^1)`$, to the step function $`N_\nu (x)`$ in (22). It is easily seen that, for $`\nu >x`$, $`\stackrel{~}{N}_\nu (x)=0`$ while, for $`\nu <x`$, we have

$$\begin{array}{c}\stackrel{~}{N}_\nu (x)=\frac{\nu }{\pi }\left(\sqrt{z^21}\mathrm{arctan}(\sqrt{z^21})\right)\\ \\ \frac{1}{4}\frac{2+3z^2}{24\nu \pi \left(z^21\right)^{3/2}},\end{array}$$
(26)
with $`z=x/\nu `$.

Let us now determine the value $`\nu _0`$ for which

$$\stackrel{~}{N}_{\nu _0}(x)=1/2.$$
(27)
To this end, we propose an expansion of the form

$$\sqrt{z_0^21}=\epsilon _1\nu _0^{1/3}+\epsilon _3\nu _0^{3/3}+𝒪(\nu _0^{5/3}),$$
(28)
which makes sense for $`\nu _01`$ and $`z_0=x/\nu _01`$. Replacing in (26) and imposing (27), the coefficients $`\epsilon _k`$ can be determined order by order in $`\nu _0^{1/3}`$, to get

$$x=\nu _0+1.857\nu _0^{1/3}+1.034\nu _0^{1/3}+𝒪(\nu _0^1)$$
(29)
or, inverting this development,

$$\nu _0=x1.857x^{1/3}+0.1155x^{1/3}+𝒪(x^1).$$
(30)
(Notice that $`\nu _0<x`$.)

Equation (29) is in excellent agreement with the expression of the first non-vanishing zero of $`J_{\nu _0}(w)`$, for a large order $`\nu _0`$ : $`j_{\nu _0,1}=\nu _0+1.8557\nu _0^{1/3}+1.03315\nu _0^{1/3}+𝒪(\nu _0^1)`$.

## IV The dominant contributions to the vacuum energy

In this Section we evaluate the first contributions to the vacuum energy obtained from the Debye expansion employed in Section II. As we will see, this allows to determine the volume, surface, …, terms in the Casimir energy of the scalar field.

According to the results in the previous Section, we are interested in the Laurent expansion of $`\zeta _\nu (s,x>\nu )`$ and $`\zeta _\nu (s,x=0^+)`$ around $`s=1`$, given in eqs. (20) and (21). As already remarked, the singular parts cancel out in the difference on the r.h.s. of eq. (10) (see eq. (19)). For the difference of the finite parts we get

$$\begin{array}{c}zF(x,\nu )\left[\zeta _\nu (s,0^+)\zeta _\nu (s,x)\right]_{s=1}=\\ \\ =\frac{\nu ^2}{2\pi }\left(z\sqrt{1+z^2}\mathrm{log}(z+\sqrt{1+z^2})\right)\\ \\ \frac{\nu }{4}+\frac{3z8z^3}{24\pi \left(1+z^2\right)^{\frac{3}{2}}}+\frac{1}{8\pi }\mathrm{log}(z+\sqrt{1+z^2})\\ \\ +𝒪\left(\nu ^1\right).\end{array}$$
(31)
This is a good approximation as long as $`\nu 1`$ and $`z=x/\nu 1`$.

Our aim is now to evaluate the sum in eq. (5),

$$\begin{array}{c}E(R=xc/\mathrm{\Omega })=\\ \\ =\frac{\mathrm{}\mathrm{\Omega }}{x}\underset{\nu =1/2}{\overset{\nu _0}{}}\nu \left[\zeta _\nu (1,x)\zeta _\nu (1,0^+)\right]=\\ \\ =\mathrm{}\mathrm{\Omega }\underset{\nu =1/2}{\overset{\nu _0}{}}F(x,\nu ),\end{array}$$
(32)
with $`\nu _0`$ given in (30).

The function $`F(x,\nu )`$ is non-negative and has a pronounced maximum at $`\nu x/2`$ (i.e. $`z2`$). Thus, the use of the approximation in eq. (31) is consistent if $`x1`$.

From eq. (31), it is not difficult to see that the successive terms in the Euler - Maclaurin summation formula ,

$$\begin{array}{c}\underset{\nu =1/2}{\overset{\nu _k}{}}F(x,\nu )=_{1/2}^{\nu _0}F(x,\nu )𝑑\nu +\\ \\ +\frac{1}{2}\left\{F(x,\nu _0)+F(x,1/2)\right\}+\\ \\ +\frac{1}{12}_\nu \left[F(x,\nu )\right]|_{\nu =1/2}^{\nu =\nu _0}+\mathrm{},\end{array}$$
(33)
are of increasing order in $`x^1`$. So, retaining the first few terms consistent with the approximation made, we get for the vacuum energy

$$\begin{array}{c}\frac{E(R)}{\mathrm{}\mathrm{\Omega }}=\\ \\ x^3\left(\frac{\left(\sqrt{1+z_0^2}\left(5+2z_0^2\right)\right)}{24\pi z_0^3}+\frac{2z_0^43\mathrm{log}(z_0+\sqrt{1+z_0^2})}{24\pi z_0^4}\right)\\ \\ x^2\left(\frac{\sqrt{1+z_0^2}}{4\pi z_0^2}+\frac{\pi +3\mathrm{log}(z_0+\sqrt{1+z_0^2})}{12\pi z_0^3}\right)\\ \\ x\left(\frac{\left(27+17z_0^2\right)}{48\pi z_0\sqrt{1+z_0^2}}+\frac{6\pi 10z_0^2+69\mathrm{log}(z_0+\sqrt{1+z_0^2})}{48\pi z_0^2}\right)\\ \\ +𝒪(x^0),\end{array}$$
(34)
where $`z_0=x/\nu _0`$, $`x=R\mathrm{\Omega }/c`$ and $`\nu _0`$ is given in eq. (30). In this equation one can recognize volume, surface and curvature contributions to $`E(R)`$. Notice that we could have retained any number of terms of the asymptotic expansion in eq. (13), and performed the same steps as in the present calculation, to get the Casimir energy to any order in $`x^1`$.

Finally replacing in eq. (34) $`z_0x/\nu _0`$ and $`\nu _0`$ by its expression we get

$$E(R)=\mathrm{}\mathrm{\Omega }\left[\frac{x^3}{12\pi }\frac{x^2}{12}0.1343x^{\frac{4}{3}}+𝒪(x)\right],$$
(35)
where one can see that the volume and surface terms are dominant for $`x>>1`$. Notice that also non-integer powers of the radius $`R`$ appear as a consequence of the relation between $`\nu _0`$ and $`x`$, eq. (30).

## V Conclusions

In equations (34-35) we have derived the dominant contributions to the vacuum energy of a scalar field in a model with a frequency dependent boundary condition, consisting in the confinement of the modes with low frequency (up to a physical cut-off $`\mathrm{\Omega }`$) to the interior of a sphere of radius $`R`$.

These modes are subject to Dirichlet boundary conditions at the surface of the sphere, while those with frequency higher than $`\mathrm{\Omega }`$ are free, being the boundary completely transparent to them. This characteristic of the model allows for the subtraction of the contribution of the high frequency modes to the vacuum energy, which amounts (independently of the regularization employed to define it) to an $`R`$-independent redefinition of the zero energy level, having no consequences on the evaluation of energy differences.

In so doing, we have represented the sum over the eigenfrequencies up to the cut-off $`\mathrm{\Omega }`$ in terms of an incomplete $`\zeta `$-function associated with the Laplacian operator in the sphere with Dirichlet boundary conditions (see eq. (10)). The function $`\zeta (s,x)`$ as in eq. (12) is well defined only for $`\mathrm{}(s)>1`$. So, it was analytically continued from $`s>1`$ to the relevant values of this parameter ($`s=0`$, needed to evaluate the maximum angular momentum $`l_0=\nu _01/2`$ giving rise to eigenfrequencies less than or equal to $`\mathrm{\Omega }`$, and $`s=1`$, necessary to evaluate the contribution to the Casimir energy of the modes with angular momentum $`l=\nu 1/2`$) by approximating the behavior of the integrand in eq. (12) employing the Debye asymptotic expansion of the modified Bessel functions appearing in its expression.

This procedure has lead to a meromorphic function having simple poles with (exactly evaluated) cut-off independent residues (see eq. (19)), a necessary condition to have a finite result for the sum in eq. (10).

The finite part of $`\zeta (s,x)`$ has been evaluated up to terms of a given order in $`\nu ^1`$. Although we have retained only the first terms in this asymptotic expansion, one can follow the same steps to determine the finite part of $`\zeta (s,x)`$ at $`s=1`$ or $`0`$ up to any required precision in $`\nu ^1`$.

Finally, the application of the Euler - Maclaurin summation formula has lead to an expression for the Casimir energy of the model in which one can recognize volume, surface and curvature contributions (see eq. (34)).

For a cut-off corresponding to $`x=R\mathrm{\Omega }/c>>1`$, the dominant terms in the vacuum energy, eq. (35), are proportional to the volume ($`V=4\pi R^3/3`$) and area ($`S=4\pi R^2`$) of the sphere,

$$\begin{array}{c}\frac{E(R)}{\mathrm{}\mathrm{\Omega }}=V\frac{\mathrm{\Omega }^3}{16\pi ^2c^3}+\xi S\frac{\mathrm{\Omega }^2}{12\pi ^2c^2}+\mathrm{}\end{array},$$
(36)
with $`\xi =\pi /4`$.

It is worthwile to remark that, for a similar model where the low frequency modes of the scalar field are subject to Neumann (rather than Dirichlet) boundary condition, we get the same expression for the dominant terms with $`\xi =+\pi /4`$.

These two dominant terms are in complete agreement with those obtained from the expansion of the density of states in powers of the inverse wavelengthIndeed, for scalar fields subject to local homogeneous boundary conditions the density of states is modified by finite volume effects . The first correction in the asymptotic expansion for large wavelength $`k`$ is given by
$$\begin{array}{c}\underset{n}{}V\frac{d^3k}{(2\pi )^3}+S\xi \frac{d^3k}{(2\pi )^3k}+\mathrm{}\end{array},$$
(37)  where the coefficient $`\xi `$ takes the value $`\xi =\pi /4`$, $`\xi =+\pi /4`$ and $`\xi =+\pi /4`$ for Dirichlet, Neumann and Robin boundary conditions respectively. Then, introducing the dispersion relation $`\omega (k)=ck`$ and a cut-off in the wavelength given by $`K=\mathrm{\Omega }/c`$, it is easy to get eq. (36) for the vacuum energy. For junction boundary conditions, as those appearing in problems with dielectrics, the coefficient $`\xi `$ has more involved expressions, which can be difficult to evaluate . . The relation between $`\nu _0`$ and $`x`$, eq. (30) (or, equivalently, the expression of the first zero of the Bessel function $`J_{\nu _0}(w)`$ in terms of the order $`\nu _0`$ ) introduces also non-integer powers of the radius $`R`$ (see eq. (35)).

As a final exercise, we can use eq. (35) in a very schematic model pretending to mimic the phenomenon of sonoluminescence. We will adopt the values of the radius and emitted energy corresponding to a typical sonoluminescent bubble , and estimate the cut-off $`\mathrm{\Omega }`$ needed to produce this amount of energy. To this end, we will simply take the difference of the low frequency contribution to the Casimir energy of the scalar field for two different values of the bubble radius.

If the bubble collapses from an initial radius $`R=4\times 10^5m`$ to a final radius of one tenth this value, and the emitted energy is $`=1.2\times 10^{12}Joule`$, by imposing the equality

$$\frac{R}{\mathrm{}c}\left(E(R)E(R/10)\right)=1.516\times 10^9,$$
(38)
and taking into account eq. (35), it follows that $`x=490`$, justifying the use of the approximation obtained.

This implies that $`\mathrm{\Omega }=3.675\times 10^{15}1/sec`$, which corresponds to a cut-off in wavelengths in the ultraviolet of $`\mathrm{\Lambda }=5.129\times 10^7m=5129\AA `$, not far from the region where the refraction index of water becomes essentially 1 . This strongly suggests to consider a similar model for the case of the electromagnetic field in the presence of dielectric media, calculation which will be presented elsewhere .

Acknowledgements: The authors thank E.M. Santangelo for useful comments. This work was partially supported by Fundaciones Andes - Antorchas under Contract Nr. C-13398/6, ANPCyT (Argentina) - PICT’97 Nr. 00039, CONICET (Argentina) - PIP Nr. 0459, and FONDECyT (Chile) - grants Nr. 1980577 and 7980011.