Benchmark problems for water radiative transfer

Introduction

On 22-24 March 2004, a workshop was held at the Lorentz Center of Leiden University (The Netherlands), to verify existing computer programs for radiative transfer in interstellar water. This web page describes the test problems and shows the results.

There are three test problems: a static two-level problem, an expanding two-level problem, and an AGB star with 45-level ortho- and para-water. The first two problems have analytic solutions and were contributed by David Neufeld (JHU). The third problem, by Jeremy Yates (UCL), should approach real astrophysics. The test campaign was coordinated by Floris van der Tak (MPIfR).

The programs participating in the benchmark fall into two categories: grid-based (ALI) codes and Monte Carlo codes. The codes labeled MALI, MULTI and Doty are ALI codes, while FLS, AMC and Bernes are Monte Carlo codes. Not all codes have produced results for each test case (yet).

Problem set 1: Two-level molecule

The first two problems use a fictive molecule that only retains the lowest two energy levels of ortho-water. The cloud has a uniform temperature (40 K) and density (10⁴ cm^-3) within radii of 0.001 pc and 0.1 pc. Here is a full description of the static two-level problem and its analytic solution. Note that for the benchmark campaign, a different value for the line width was adopted. To see how this change affects the analytic solution, click here.

The participants sent numerical results in the form of the population of the upper level as a function of radius. Here are those results for a water abundance of 10^-10, and here are those for an abundance of 10^-8.

The level populations can be converted into excitation temperatures, defined through the Boltzmann equation. Here are the results of the static problem for low abundance in excitation temperature form, and here are those for the high water abundance. These plots show that radiative deexcitation dominates over collisional excitation, because the excitation temperature is much smaller than the kinetic temperature.

The dash-dotted line is the excitation temperature calculated using the `local' approximation (with an escape probability code, comparable to the `large velocity gradient' method). Clearly, under these conditions (high optical depth and fast radiative decay), the correct solution is only found by programs that resolve the cloud spatially and kinematically.

From these results, the expected sky brightness was calculated by a ray tracer, and the result convolved with a beam much larger than the cloud. Here are the line profiles for low abundance and for high abundance.

To compare with the analytic solution, the line profiles were integrated over velocity. As seen from this plot, all codes recover about 80% of the theoretically expected emission (the dashed line) at low abundances. At higher abundances, the MULTI and MALI codes are off by an additional factor of about 2.

Problem set 2: An expanding sphere

The second set of problems also considers the fictive two-level molecule, under the same physical conditions. However, the gas cloud is now isotropically expanding at a rate of 100 km/s/pc. This document describes the problem in detail, and gives its analytic solution. (Note a typo in Equation 1: the Sobolev optical depth is proportional to the cube, not the square of the wavelength. Also, the numerical coefficient in the expression for t just below Equation 1 should be 1.51066e-08, a factor 1.8 higher than reported. Thanks to Eric Chung for pointing this out; see this Web page for details).

The numerical results for the expanding sphere are here as populations for low abundance, here as populations for high abundance, here as excitation temperatures for low abundance, and here as excitation temperatures for high abundance.

The excitation temperature should monotonically increase towards the center of the cloud. Some programs indicate a decrease, which is an artifact of gridding the velocity field along with the temperature and the density. This gridding removes the isotropic character of the expansion and gives the innermost cells too little optical depth.

At small radii, the analytic solution indicates an upper level population of 5.613 10^-4 for the high abundance (corresponding to an excitation temperature of 3.57 K). For the low abundance, a population of 3.250 10^-4 (or 3.33 K excitation temperature) is expected. The results computed with AMC approach these values closely, although the Monte Carlo noise is not negligible. Other Monte Carlo codes show less noise, but seem to systematically underestimate the excitation.

Problem set 3: An AGB star

The last set of problems considers a model with temperature and density gradients as they may occur in the envelopes of AGB stars. The radiative transfer now needs to be solved for ortho- and para-water (at a 3:1 ratio) with 45 energy levels each. In addition, the envelope contains warm dust, whose continuum radiation must be taken into account.

In going from the inner radius (2 AU) to the outer radius (25 AU), the density drops from 5 10¹⁰ to 9 10⁷ cm^-3. The kinetic temperature drops from 2000 to 67 K and the dust temperature from 2000 to 240 K. The dust opacity is parametrized as a broken power law. This document gives a detailed description of the AGB test problems. The benchmark campaign used a water abundance of 10^-6.

There is no analytical solution for this test problem. Here are numerical solutions of the population of the ground and first excited states of para-water, and of the ground and first excited states of ortho-water, and the excitation temperatures of the ortho ground state line (557 GHz) and the para ground state line (1113 GHz).

Conclusions

Some conclusions of this benchmark campaign are:

1. For the two-level problem, agreement between codes is within 5% for high abundances, and much better for low abundances. For the AGB star problem, agreement is within 20%.

2. Velocity fields need to be programmed in functional form, not just gridded along with density and temperature.

3. To ensure proper convergence of stiff problems, results must be calculated from `cold' (optically thin) and `warm' (LTE) starting conditions, and convergence declared when the two agree. Merely declaring convergence when differences between successive iterations are small leads to poorly defined results.