text:
# Numerical methods for non-LTE line radiative transfer: Performance and convergence characteristics

## 1 Introduction

Molecular lines are excellent probes of the physical and chemical conditions in interstellar clouds, protostellar envelopes, circumstellar shells around late-type stars, photon-dominated regions etc. Furthermore, molecular line transitions play a key role in probing the properties of galaxies and their evolution. The interpretation of such lines requires the use of line radiative transfer programs which can calculate accurately the non-LTE (local thermodynamic equilibrium) level populations and the resulting output spectra. See Black (2000) for a recent review.

It is known from stellar atmosphere research that subtle errors in radiative transfer algorithms can lead to significantly incorrect results (Mihalas 1978). A particularly well-known problem of this kind is an insufficiently stringent convergence criterion at high optical depths. In the absence of a-posteriori checks on numerical results, the best way to validate the methods is by the use of various techniques, if possible with many independent codes of different types. Once the reliability and the behavior of the codes has been established and understood, they can be safely used within the limits of parameter space within which tests have been carried out.

In this paper, the setup and results of a large-scale campaign to compare line radiative transfer programs are described. The test problems are spherically symmetric, and can be modeled by a 1-dimensional radiative transfer code. Codes capable of handling more than one dimension could be compared in a similar way in the future.

The main aim of this paper is to examine and compare the various radiative transfer methods that are currently used in the astrophysical community for modeling the shape and strength of molecular lines emerging from interstellar and circumstellar matter. A series of problems is chosen that represent the difficulties that may be encountered in comparison with data from present-day and future ground-based telescopes, such as the Sub-Millimeter Array (SMA), the Atacama Large Millimeter Array (ALMA) and future airborne and space-borne missions such as SOFIA and the Herschel Space Observatory. The test problems and their solutions presented here are available to the community via a web-site, allowing users to run the same problems, check the accuracy of their codes and thus speed up the further development of their own radiative transfer codes. It is hoped that this will stimulate a more widespread use of these codes for the interpretation of molecular line observations.

We present two test problems, each at two different optical depths (i.e. in total four problems). The first problem is a fictive 2-level molecule in a spherical envelope with a powerlaw density with no systematic velocity and a constant temperature. This problem is meant to tune all codes before examining the more complex, realistic problem. The main test problem is based on the 1-dimensional (1-D) inside-out collapse model, where the level populations are computed for HCO+ at various abundances. The HCO+ ion is chosen as a representative example of the molecule which samples gas with a large range in densities and is readily observable in a variety of astrophysical regions. Both a high and a low optical depth model, representative of the main isotope HCO+ and the less abundant isotope H13CO+, are used to check the convergence properties of each code. All test problems and their results can be found at the webpage http://www.strw.leidenuniv.nl/$``$radtrans.

## 2 Molecular line radiative transfer

The radiative transfer problem is represented by an equation describing the emission, absorption and movement of photons along a straight line in a medium:

$$\frac{dI_\nu }{ds}=j_\nu \alpha _\nu I_\nu .$$
(1)
with the notation adopted from Rybicki & Lightman (1979). Equation (1) is the differential description of the intensity $`I_\nu `$ along a photon path ($`ds`$) at frequency $`\nu `$, where $`j_\nu `$ \[erg s-1 cm-3\] and $`\alpha _\nu `$ \[cm-1\] denote the emission and absorption coefficients. Another way of writing Eq. (1) is:

$$\frac{dI_\nu }{d\tau _\nu }=S_\nu I_\nu ,$$
(2)
where $`S_\nu =j_\nu /\alpha _\nu `$ is referred to as the source function (the emissivity of the medium per unit optical depth), and the optical depth $`\tau _\nu `$ is defined in differential form as $`d\tau _\nu =\alpha _\nu ds`$. Equation (2) can be written in integral form, which is the form that is most often used in radiative transfer codes:

$$I_\nu (\tau )=_0^\tau S_\nu (\tau ^{})e^{\tau ^{}\tau }𝑑\tau ^{},$$
(3)
where $`\tau `$ is the optical depth between the point where $`I_\nu `$ is evaluated and spatial infinity along the line (i.e., $`s=\mathrm{}`$). This integral is evaluated along all possible straight lines through the medium. In practice, these will be a discrete sample of lines covering space and direction as well as possible.

For the problem of molecular line transfer the emission and absorption coefficients are determined by the transition rates between the various rotational and/or vibrational levels of the molecule, and the population of these levels. For a transition from level $`i`$ to level $`j`$ (where the energy of level $`i`$ is greater than that of level $`j`$) the emission and absorption coefficients are given by:

$`j_{ij}(\nu )`$ $`=`$ $`n_iA_{ij}\varphi _{ij}(\nu ),`$ (4)
$`\alpha _{ij}(\nu )`$ $`=`$ $`(n_jB_{ji}n_iB_{ij})\varphi _{ij}(\nu ),`$ (5)
where $`n_i`$ and $`n_j`$ \[cm-3\] are the population densities of the upper ($`i`$) and lower ($`j`$) level, and $`A_{ij}`$, $`B_{ij}`$ and $`B_{ji}`$ are the Einstein coefficients. The function $`\varphi _{ij}(\nu )`$ represents the line profile for this transition, which is properly described as a Voigt profile, a combination of a (micro-turbulent) Gaussian and intrinsic Lorentzian line broadening.

The source function for a particular transition $`S_{ij}`$ is independent of velocity if one assumes complete frequency redistribution, i.e., the frequency deviation from line center of absorbed and emitted photons are uncorrelated. The source function then becomes:

$$S_{ij}=\frac{j_{ij}(\nu )}{\alpha _{ij}(\nu )}=\frac{n_iA_{ij}}{n_jB_{ji}n_iB_{ij}}.$$
(6)
The relative level populations $`n_i`$ are determined from the statistical equilibrium equation:

$$\begin{array}{c}\hfill \underset{j>l}{}[n_jA_{jl}+(n_jB_{jl}n_lB_{lj})\overline{J}_{jl}]\\ \hfill \underset{j<l}{}[n_lA_{lj}+(n_lB_{lj}n_jB_{jl})\overline{J}_{lj}]\\ \hfill +\underset{j}{}[n_jC_{jl}n_lC_{lj}]=0,\end{array}$$
(7)
where $`C_{ij}=n_{\mathrm{col}}K_{ij}`$ with $`K_{ij}`$ the collisional rate coefficients in cm3 s-1 and $`n_{\mathrm{col}}`$ the density of collision partners, taken here to be H2 in $`J`$=0. $`\overline{J}_{jl}`$ is the integrated mean intensity over the line profile:

$$\overline{J}_{ij}=\frac{1}{4\pi }I_\nu (\mathrm{\Omega })\varphi (\nu )𝑑\mathrm{\Omega }𝑑\nu .$$
(8)
The symbol $`\mathrm{\Omega }`$ represents the spatial direction in which the intensity $`I_\nu (\mathrm{\Omega })`$ is measured.

A useful concept for the analysis of radiative transfer calculations is the excitation temperature $`T_{\mathrm{ex}}`$ of the transition between level $`i`$ and level $`j`$, given by

$$T_{\mathrm{ex}}=\frac{h\nu _{ij}}{k}\left[\mathrm{ln}\left(\frac{g_j}{g_i}\frac{n_i}{n_j}\right)\right]^1,$$
(9)
where $`k`$ is Boltzmann’s constant and $`g_i`$ the statistical weight of level $`i`$. The energy difference between the two levels is given by $`h\nu _{ij}`$. In local thermodynamic equilibrium (LTE), $`T_{\mathrm{ex}}`$ equals the local gas temperature, while if $`T_{\mathrm{ex}}`$ higher or lower, the excitation is super- or sub-thermal. In addition, the intensity is proportional to $`T_{\mathrm{ex}}`$ in the optically thin limit.

Equations (3)-(8) form a coupled set of equations. The way in which the quantities depend on each other is depicted in Fig. 1. The intensities are found by integrating the source function (Eq. 3). The source function and extinction coefficients depend on the level populations (Eqs. 4–5). These in turn depend on the intensities (Eqs. 7–8). To solve this set of equations one must determine the radiation field and the level populations simultaneously. Since the radiation field couples the level populations at different spatial positions to each other through the transfer integral (Eq. 3), the only way of solving this system directly seems to be through complete linearization and solving a huge matrix equation involving all level populations at all spatial positions. In practice, the evaluation of the matrix elements and the inversion of the problem consumes far too much CPU-time as well as memory and is therefore beyond current computing capabilities. An alternative and much simpler way is to iteratively evaluate all equations, following the arrows, as illustrated in Fig. 1. This method is called “Lambda Iteration” and includes formal integral codes as well as Monte Carlo methods. Though simple, its principal disadvantage is that it converges slowly at high optical depths. Many radiative transfer codes therefore use a hybrid scheme: the direct inversion of a simplified subset of the equations, and iteration to solve the remaining problem. In this paper we have used codes of the Lambda Iteration type (including Monte Carlo methods), of the “Complete Linearization” type, and hybrid schemes such as “Approximate/Accelerated Lambda Iteration” and “Accelerated Monte Carlo”.

## 3 Methods and codes

In this paper four methods and eight different codes are compared. Any practical method discretizes the problem to determine the level populations at each position. A short description is given in this section of each of these codes. The first method (LI) is not used, but given as an introduction for the other methods. In Table LABEL:tab:\_codes, the different codes and principal authors are indicated. Detailed descriptions of each of the codes can be found in the references given.

### 3.1 Lambda Iteration (LI)

The Lambda Iteration method is the most basic method of all. It involves the iterative evaluation of level populations and intensities until the system has converged. The name “Lambda Iteration” originates from the fact that the process of iteration can be mathematically written into a formalism involving a “Lambda Operator”. This Lambda Operator represents the entire procedure of computing the line-weighted mean intensities $`\overline{J}_{ij}`$ from the source function. It involves the formal integral Eq. (3) along all possible lines through the medium, and includes the angle-frequency integrations of Eq. (8) to obtain the $`\overline{J}`$. The Lambda Operator is defined as:

$$\overline{J}_{ij}=\mathrm{\Lambda }_{ij}[S_{ij}],$$
(10)
and is therefore a global operator. The $`\mathrm{\Lambda }_\nu `$ operator is a matrix operator connecting all points and all levels to each other. To solve the radiative transfer numerically, the problem has to be discretized both in space and frequency and initial level populations have to be assumed. With these level populations, the Lambda Operator is constructed by solving the radiative transfer for $`I_\nu `$ and $`\overline{J}`$, which is inserted in Eq. (7) to calculate the new level populations. The procedure is repeated until the relative change in the level populations or mean intensities between two successive steps falls below a desired convergence criterium. The time to reach convergence of the solution is proportional to $`\tau ^2`$ for $`\tau 1`$ and is therefore extremely slow for highly optically thick regions. The very small local changes in the level populations in successive iterations could then be easily mistaken for convergence.

### 3.2 Monte Carlo (MC)

The Monte Carlo (MC) method is based on the simulation of basic physical processes with the aid of random numbers. This makes the formalisms of these codes relatively simple and intuitive as one has to deal only with the basic formulae. But one must take proper care of the statistics and be sure that all regions in space and frequency are well sampled by the randomly distributed photon packages. The Monte Carlo method for molecular line transfer has been described by Bernes (1979) for a spherically symmetric cloud with a uniform density. A major advantage of MC codes is the possibility of non-regular grids, in particular for multi-dimensional problems. In fact, MC methods are straightforward to extend from 1-D to 2-D (Hogerheijde & van der Tak 2000) and even 3-D (Park & Hong 1998, Juvela 1997). One of the major disadvantages of the method is the long CPU time which is needed to lower the random error intrinsic in the method. This error decreases inversely proportional to the square root of the number of simulated photons. In addition, the method suffers from similar convergence problems as Lambda Iteration, as it is formally a variant of the Lambda Iteration method.

Thus far, the Monte Carlo method has been implemented in two ways. In one set of codes, the radiation field is randomly sampled by discrete photons that are emitted and absorbed in the material. In the other set, the radiation field is sampled by random directions along which the radiative transfer is evaluated.

### 3.3 Approximated/Accelerated Lambda Iteration (ALI)

This scheme is similar to the Lambda iteration scheme, with the difference that the equations are pre-conditioned to speed up the convergence. The simplest (but still highly efficient) way to implement the ALI scheme is by splitting the Lambda Operator into the local self-coupling contribution and the remainder:

$$\mathrm{\Lambda }_{ij}=(\mathrm{\Lambda }_{ij}\mathrm{\Lambda }_{ij}^{})+\mathrm{\Lambda }_{ij}^{}$$
(11)
with $`\mathrm{\Lambda }_{ij}`$ the lambda operator between the $`i`$ and $`j`$ levels and $`\mathrm{\Lambda }_{ij}^{}`$ the diagonal (local), or tridiagonal (local+nearest-neighbor), part of the full Lambda operator. In general more reliable convergence properties are obtained with tridiagonal or higher bandwidths (Hauschildt et al. 1994). Since an approximate operator $`\mathrm{\Lambda }_{ij}^{}`$ of this kind is easily invertible, and since its matrix elements are relatively easy to calculate, one can use this operator for the matrix inversion, and use the iteration for solving the remainder (non-local part) of the transfer problem. Since the self-coupling of a cell at high optical depth is the bottle-neck that slows down convergence in Lambda Iteration codes, this removal of the local self-coupling by direct inversion of $`\mathrm{\Lambda }_{ij}^{}`$ can speed up the convergence drastically. A full description of the method was given by, e.g., Rybicki & Hummer (1991, 1992).

In addition to this operator splitting technique, one can apply certain iteration-improvement schemes such as the Ng-acceleration method (Ng (1974)). These are methods that can improve the convergence of any linearly converging iteration scheme. In the Ng-method one uses the previous results to estimate the convergence behavior of the problem. After every four iteration steps an extrapolation can then be performed towards the expected convergence. The number of iterations is not a strict requirement but four is typically found to give a reliable and significant acceleration. This scheme is very effective as can be seen in Figs. 1 to 3 in Rybicki & Hummer (1991), where the number of iterations is plotted versus the convergence.

### 3.4 Accelerated Monte Carlo (AMC)

The difference in Accelerated Monte Carlo and Monte Carlo methods can be similarly described as the ALI compared to LI method. Most important is the splitting of the mean intensity in an external field and a local contribution. Formally,

$`\overline{J}`$ $`=`$ $`(\mathrm{\Lambda }\mathrm{\Lambda }^{})[S_{ul}^{}(\overline{J})]+\mathrm{\Lambda }^{}[S_{ul}(\overline{J})]`$ (12)
$`=`$ $`\overline{J}^{external}+\overline{J}^{local}`$ (13)
$`=`$ $`{\displaystyle \frac{1}{N}}{\displaystyle \underset{i}{}}I_{0,i}e^{\tau _i}+{\displaystyle \frac{1}{N}}{\displaystyle \underset{i}{}}S_{ul}[1e^\tau ]`$ (14)
with S$`{}_{ul}{}^{}(\overline{J})`$ the results from the previous iteration and $`N`$ the number of photons. This is an important issue as in the standard Monte Carlo approach, most time is taken by photons trapped in an optically thick cell, where due to absorption and randomly directed emission the photons follow a random walk path through the cell, with one step for each iteration.

The splitting can be done in different ways. Juvela (1997) implements a diagonal lambda operator in Eqs. (12)-(14) by counting explicitly those photons which were emitted and absorbed within the same cell. A reference field described by Bernes et al. (1979) is used to decrease the random fluctuations caused by the Monte Carlo sampling. Hogerheijde & van der Tak (2000) perform sub-iterations to calculate the local contribution in each cell, thus ensuring that in each cell the local field and populations are always consistent with the external field due to all other cells. A third possibility of acceleration, adopted by Schöier (2000), is the use of core-saturation (Rybicki 1972, Hartstein & Liseau 1998), where optically thick line center photons are replaced by the local source function.

### 3.5 Local Linearization (MULTI type codes)

This method was developed by Scharmer & Carlsson (1985) and Harper (1994) to produce the MULTI and SMULTI line radiation transport codes respectively. This approach perturbs Eqs. (2) to (8) linearly, neglecting second order or higher terms. The linearization greatly reduces the effect of high optical depth terms and allows rapid convergence. The MULTI and SMULTI codes use the Olson diagonal approximate operator (Olson, Auer, & Buchler (1986)) scheme to save on storage. This operator uses the diagonal of the $`\mathrm{\Lambda }`$ matrix when computing the solution to the linearized equations. Although this adds in approximations and delays convergence, it greatly reduces the required storage capacity and so for many problems it allows the problem to be tackled in the first place. Olson’s scheme simply assumes that the changes in the intensity are related to the changes in the source and opacity terms (the level populations). The main difference to the other schemes is that the solution to the linearized equations returns changes, $`\delta n`$ and $`\delta J`$ to the level populations and the mean intensities respectively. There is always the possibility that a converged solution contains unphysical negative populations, and this provides a useful indicator for poor sampling and errors. In the version used here a ALI-type scheme to complement the MULTI is used, providing numerical stability but also slower convergence in optically thick media.

### 3.6 Radiative transfer codes

Each code described in this paper uses one or more of the above techniques to calculate a converged set of level populations. In this section we list the codes which use each of the convergence methods listed above. Then for each code we describe how convergence is accelerated, how each code samples the volume under consideration, and state the convergence criteria used by each code.

1. Monte Carlo (MC):

* F. Schöier (Schöier (2000)); The rate of convergence depends on the number of model photons and the iterative procedure. In the problems presented here, the counters of stimulated emission are reset after each set of 5 iterations. The number of iterations needed for convergence in the “classical” Monte Carlo scheme is of the same order as the maximum optical depth in the model. The core saturation method is included to speed up the convergence in the high optical depth case. To reach an accuracy of $``$10-2 in the derived level populations, $`10^510^6`$ model photons per iteration are generally needed in the optically thick case.
2. Accelerated Lambda Iteration (ALI):

* V. Ossenkopf (Ossenkopf, Trojan, & Stutzki (2001)); All discretizations are done by the code such that sufficiently small differences between the grid points are guaranteed. This concerns the radial shells, the grid of rays, the frequency grid, and the number of levels used. The convergence is measured in terms of radiative energy densities, not level populations, which are extrapolated from the Auer acceleration scheme. Locally, the code may take into account turbulent sub-clumping using the statistical description of Martin, Sanders & Hills (1984), but this is not used here.
    * S. Doty; the input data are interpolated onto the spatial grid of the problem with a variable method (log/linear/combined). This, coupled with extreme care in the ray-tracing integration, helps to ensure high accuracy. A local approximate lambda operator and Ng acceleration are used to solve both the line and continuum transfer. Convergence is determined by changes in both level populations and the local net heating rate, with small ($`10^4`$) changes required for the entire spatial grid for both accelerated and non-accelerated iterations.
    * C.P. Dullemond (Dullemond & Turolla (2000)); The formal integration proceeds via a short-characteristics method using the three-point quadrature formula of Olson & Kunasz (1987). The discrete local angles $`\mu _i`$ are chosen using the tangent-ray prescription. This means that successive short characteristics match up and no interpolations are required. At every grid point there are 41 angular points in $`\mu `$ with three extra around $`\mu 0`$ to prevent undersampling. Using the ALI method with a local operator and Ng acceleration, the system is iterated towards a solution with a convergence criterion $`10^6`$ in level population.
    * H. Wiesemeyer (Wiesemeyer 1997) ; The rate of convergence is improved by applying the method of minimization of residuals (Auer 1987). Other acceleration methods, such as vector and matrix extrapolation, are currently being implemented. The solution to the equation of radiative transfer is evaluated by various numerical methods (either quadrature rules, or multi-value methods), according to the accuracy required by the problem to be solved. The code uses long characteristics, to preserve an isotropic distribution of rays at any point of the model volume, following the multivariate quadrature rule of Steinacker et al. (1996). Its performance is thus rather optimized to solve smooth problems in more than one dimension. All discretizations are performed by the code.
3. Accelerated Monte Carlo (AMC):

* M. Juvela (Juvela 1997); This code uses the reference field method to reduce random fluctuations. The velocity is discretized into 50 channels and the number of rays per iteration was taken to be 350. The ray generation is weighted such that the same number of rays was shot through each annulus. The random number generators are reset after each iteration, making it possible to use Ng acceleration on every third iteration starting with the fifth iteration. Convergence was tested only for the eight lowest levels down to 1$`\times 10^3`$ in level population.
    * M. Hogerheijde & F. van der Tak (Hogerheijde & van der Tak (2000)); In this code, one cell at a time is considered, with $`N`$ rays entering the cell from random directions. The radiative transfer is followed along each of these rays, starting with the CMB field at the edge of the cloud. In this way it is possible to calculate separately the local contribution to the radiation field in the target cell, allowing significant reduction of the computing time for optically thick cells. For a first order estimate of the radiation field, the same set of random numbers is used to describe the ray directions for each iteration, thus resembling a fixed-ray (ALI) code. This process yields a solution free of random noise but with possible inadequate sampling of directions and velocities. In the second stage a different set of random numbers is used for each iteration, providing true random sampling of the radiation field. The number of rays is increased for each cell until convergence is reached, ensuring proper angular sampling of the radiation field everywhere. The equations of statistical equilibrium are solved in each iteration to a fractional error of $`10^6`$ in all levels except the highest. The user specifies a signal-to-noise ratio $`S`$; the MC noise is reduced by increasing the number of photons such that the fractional error of the levels between iterations is smaller than 1/$`S`$. For the test problem, $`S`$=100 was used.
4. Complete Linearization (MULTI-type):

* J. Yates (Rawlings & Yates (2001)); In the problems presented, the SMMOL code used 400 rays through the cloud to compute the intensities at each radial grid-point in the spherical cloud using a finite difference method to compute the intensities. The grid sampling along each ray can adapt simply to take account of large changes that affect the optical depth, typically caused by velocity changes along the ray. In this calculation, typically 100 grid-points were used to sample the full absorption profile. The SMMOL codes uses a convergence criterion of $`1\times 10^4`$ for both the level populations and the mean intensities.

## 4 Description of the test models

### 4.1 Models 1a/1b: a simple 2-level molecule

The first two test problems (problems 1a/1b) have a simple and cleanly defined setup without velocity gradients and using a constant temperature and line width. Complicating factors such as a different treatment of the spatial grid and frequency sampling are thus kept to a minimum, so that every code should in principle be able to solve these problems down to their specified accuracy.

The simple 2-level molecule test setup consists of a spherically symmetric cloud with a power law hydrogen number density specified by

$$n_{\text{H}_2}(r)=n_{\text{H}_2}(r_0)\left(\frac{r}{r_0}\right)^\alpha \text{cm}^3,$$
(15)
where $`r`$ is the radius in cm. We take $`n_{\text{H}_2}(r_0)=2.0\times 10^7\text{cm}^3`$, $`r_0=1.0\times 10^{15}\text{cm}`$ and $`\alpha =2`$. The kinetic temperature of this problem is taken to be constant: $`T_{\text{kin}}(r)=20\text{K}`$. The abundance of the molecule (which we will specify below) is constant as well: $`X_{\text{mol}}(r)n_{\text{mol}}(r)/n_{\text{H}_2}(r)=X_{\text{mol}}`$, and the systematic velocity is zero everywhere. The spherically symmetric cloud ranges from $`r_{\text{in}}=1.0\times 10^{15}\text{cm}`$ to $`r_{\text{out}}=7.8\times 10^{18}\text{cm}`$. For $`r<r_{\text{in}}`$ and for $`r>r_{\text{out}}`$ the density is assumed to be zero. The incoming radiation at the outer boundary is the $`T=2.728\text{K}`$ microwave background radiation.

We choose a fictive 2-level molecule which is specified by

$`E_2E_1`$ $`=6.0\text{cm}^1=8.63244\text{K}`$ (16)
$`g_2/g_1`$ $`=3.0`$ (17)
$`A_{21}`$ $`=1.0\times 10^4\text{s}^1`$ (18)
$`K_{21}`$ $`=2.0\times 10^{10}\text{cm}^3\text{s}^1`$ (19)
in which the downward collision rate is $`C_{21}=n_{\text{H}_2}K_{21}\text{s}^1`$. The total (thermal+turbulent) line width $`a`$ is given by $`a=0.150\text{km}\text{s}^1`$. The line profile is assumed to be a Gaussian:

$$\varphi (\nu )=\frac{c}{a\nu _0\sqrt{\pi }}\mathrm{exp}\left(\frac{c^2(\nu \nu _0)^2}{a^2\nu _0^2}\right)$$
(20)
where $`c`$ is the speed of light in units of $`\text{km}\text{s}^1`$, and $`\nu _0`$ is the frequency at line center.

We solve the problem for two different abundances:

Problem 1a: $`X_{\text{mol}}`$ $`=1.0\times 10^8`$
Problem 1b: $`X_{\text{mol}}`$ $`=1.0\times 10^6`$
Problem 1a is a simple case with a moderate optical depth ($`\tau 60`$ from star to infinity at line center). Problem 1b has much higher optical depth ($`\tau 4800`$) and is therefore numerically stiffer.

This problem was originally described and worked out by Dullemond & Ossenkopf (see Dullemond 1999).

### 4.2 Models 2a/2b: A collapsing cloud in HCO+

The second two problems (problems 2a/2b) are examples of problems typically encountered in the field of sub-millimeter molecular line modeling. These problems have much of the complicated physics included, such as velocity- and temperature gradients, non-constant line widths, multiple levels etc. Moreover, the problem is based on a numerical model provided externally, and therefore introduces the complication of a given intermediate-resolution spatial grid. The resulting spread of the results therefore gives a good indication of how accurate the predicted line strengths and shapes for a typical every-day-life model calculation are.

The model on which these “realistic” test problems are based was constructed by Rawlings et al. (1992, 1999) to analyze HCO+ data for an infalling envelope around a protostar. The prototypical example of such a ‘class 0’ young stellar object is B335. The model describes how the cloud core collapses from the inside-out. Starting from a nearly isothermal sphere in pressure balance, a perturbation triggers the center of the cloud to collapse. This sends out a rarefaction wave at the local speed of sound, which causes the outer parts of the cloud to collapse as well. The model is similar to the analytical inside-out collapse model by Shu (1977), but includes more realistic physics. No rotation is assumed, which makes it ideal for a 1-D spherical comparison.

Figure 2 shows the structure of the cloud at a particular time during the collapse phase. This is the input model for our test cases. The collapse can clearly be seen in the radial velocity, which is 0 for radii greater than 1017 cm, and directed towards the source for smaller radii. The density profile is given by a power-law of the form $`n(r)=n_0(r/r_0)^m`$, where $`m=1.5`$ inside the collapsing sphere ($`r<`$ 1017 cm) and $`m=2.0`$ outside. The model parameters are specified at 50 radii logarithmically spaced between $`10^{16}`$ and $`4.6\times 10^{17}`$ cm.

The infall of a cloud has a significant influence on the molecular line emission. The lines will be skewed to the blue, or double peaked with a stronger blue-shifted component, depending on the velocity field and optical depth of the line (e.g. Myers et al. 2000). This feature can be intuitively understood because the foreground gas has a lower excitation temperature than the background gas. This only holds when the foreground gas is optically thick enough in the line. Otherwise a single-peaked lineshape, centered on the source velocity, results.

For the comparison presented here, the ion HCO+ is used, where an abundance compared to H2 is adopted of 1$`\times `$10-9 for problem 2a, and 1$`\times `$10-8 for problem 2b. This causes the optical depths $`\tau `$ for the lowest 4 transitions to range from 0.1 to 10 for problem 2a, and from 1 to 100 for problem 2b. A common file of Einstein-$`A`$ coefficients and collisional rate coefficients of HCO+ with H2 in $`J`$=0 (Monteiro (1985), Green (1975)) has been used in all models. These can be downloaded from http://www.strw.leidenuniv.nl/$``$radtrans. 21 levels of the HCO+ ion are included in the model. Since the rate coefficients are temperature dependent, all downward (upper$``$lower) rate coefficients are given for a number of temperatures (10, 20, 30, 40 K). At each temperature, a linear interpolation is performed to calculate the local downward collisional rate coefficients. Performing this interpolation accurately is extremely important because the exponential relation between the excitation and de-excitation coefficients can introduce large deviations. The upward (lower$``$upper) collisional rate coefficient is subsequently calculated using the detailed balance relation,

$$\frac{C_{lu}}{C_{ul}}=\frac{g_l}{g_u}e^{E_{ul}/kT}.$$
(21)
The fact that the molecular data, in particular the collisional rate coefficients, are not known with infinite precision introduces an additional uncertainty in the solution of the radiative transfer equations. The use of different collisional rate coefficients has an impact on the level populations and thereby on the strength of the predicted emission lines. However, the spread in available rate coefficients in the literature has no effect on the comparison presented here, as all the codes use the same molecular data as input. On a webpage http://www.strw.leidenuniv.nl/$``$moldata molecular files for a large number of molecular species will be made available in the future (Schöier et al., in preparation), including a comparison of different rates from the literature and their effect on the model results.

With these sets of problems the different codes are tested on the following features:

* Convergence of the code
* Radiative transfer in optically thin lines
* Radiative transfer in optically thick lines
* Sampling of the radiation transfer in the presence of a velocity field
* Incorporation of the cosmic microwave background radiation ($`T_{\mathrm{CMB}}`$=2.728 K), which influences the lower levels near the outside of the cloud.

## 5 Results

### 5.1 Results of problems 1a/1b

Test problems 1a and 1b are relatively simple in the sense that no velocity gradients and temperature variations are present. Also, since the problem setup is specified analytically, the grid resolution can be chosen to be as high as required for the particular code the participant is using. By virtue of the zero systematic velocity, the treatment of the inner boundary is not much of an issue, because the line is perfectly thermalized at that radius. These test problems therefore purely test the radiative transfer as such.

The results for the problems 1a and 1b are shown in Fig. 3.

The interior of the cloud is completely thermalized ($`T_{\mathrm{ex}}=T_{\mathrm{kin}}`$). The decline of $`n_2`$ at larger radii is due to the departure from LTE. At the outer radii the population again saturates to a constant value, which is due to the microwave background radiation ($`T_{\mathrm{ex}}=T_{\mathrm{CMB}}`$). All codes agree reasonably well, although small differences remain distinguishable. The relative differences between the codes are shown in Fig.4 in comparison to the mean.

If we consider the ’mean’ solution to be the correct one, then the typical errors are of the order of a few percent or less. The errors are clearly greater for problem 1b than for problem 1a. This is to be expected, since problem 1b has higher optical depth and is therefore numerically more challenging. For problem 1a, the codes show a small spread with a standard deviation of 0.2 %. i.e. they agree with each other to within $`210^3`$ relative difference, comparable to the accuracy criteria of the codes, clearly visible in the noise for the Monte Carlo results. For problem 1b the best agreement reached between the codes shows a standard deviation with a local maximum of 2% $`(\mathrm{\Delta }n_2/n_2210^2)`$. To achieve maximum accuracy, the codes used a spatial resolution of $`(R_{i+1}R_i)/R_i0.05`$ for this problem, but this high spatial resolution did not reduce the spread between the codes to the $`10^3`$ level.

For problem 1b (the most difficult of the two) an ALI code with Ng acceleration with spatial resolution $`(R_{i+1}R_i)/R_i0.04`$ typically takes about 160 iterations to meet the converge criterion of $`\mathrm{\Delta }n_2/n_2<10^7`$. The same problem, but for a less stringent convergence criterion ($`\mathrm{\Delta }n_2/n_2<10^4`$) and for a lower spatial resolution ($`(R_{i+1}R_i)/R_i0.2`$) typically takes about 40 iterations. The relative error introduced by the lower spatial resolution and lower convergence criterion for this problem is at most 1%, provided that second-order (or higher) integration of the transfer equation is used. The use of first-order integration introduces errors of the order of $`12\%`$, even for the high-resolution runs. For test problem 1a, at low resolution and with a weak convergence criterion ($`\mathrm{\Delta }n_2/n_2<10^4`$), an ALI code typically uses 18 iterations.

For the AMC-type codes with Ng acceleration (Juvela 1997), problem 1b typically requires about 75 iterations, with about 120 photon packages (rays) per iteration. It was found that the application of Ng acceleration at the very start of the iteration had a negative influence on the convergence rate. Good convergence was achieved by applying Ng only after about 25 iterations. As in the case of the ALI type codes, problem 1a was much more benign: 15 iterations with about 120 photon packages (rays) per iteration.

It should be said that the AMC-type codes allow considerable freedom in how to disentangle the number of iterations and the number of photon packages per cell per iteration. Therefore the numbers given here may vary from code to code, i.e., in most codes the number of photon packages per iteration per cell is used. For the accelerated Monte Carlo codes, acceleration of the convergence only operates when some of the cells are optically thick. If all cells are optically thin but the entire source optically thick, the local radiation field becomes negligible compared to the overall radiation field and separating local and global contributions no longer effectively speeds up convergences.

The SMULTI code converged to below 10-4 in 30 iterations and to below 10-5 in 50 iterations for the second model. For the first model the code converged in 10 iterations.

### 5.2 Results of problems 2a/2b

#### 5.2.1 Level populations

The results for the test problems 2a and 2b are shown in Fig. 5. The direct comparison is shown for two different levels of HCO+, $`J`$=1 and $`J`$=4. These levels were chosen because they represent the emitting levels of two easily and often observed lines of the molecule probing different density regimes. The $`J`$=1–0 line lies at millimeter wavelengths at 89 GHz and the critical density of the $`J`$=1 level is $`3.4\times 10^4`$ cm-3. The $`J`$=4–3 line occurs in the sub-millimeter at 356 GHz, with a critical density for exciting the $`J`$=4 level of $`1.8\times 10^6`$ cm-3. The different model results are plotted on top of each other, together with one additional calculation where the CMB radiation field was deliberately ignored. In the level populations, the effect of the CMB is visible as a lack of excitation of the $`J`$=1 level in the outer regions. The figure immediately shows that in the outer region, the excitation at line center of the $`J`$=1–0 line is controlled by the CMB field and not by the local temperature and density. The $`J`$=1 level lies only 4.3 K above ground, close to the $`T_{\mathrm{CMB}}`$ (2.728 K) temperature, so that the 1–0 transition can be effectively excited by the peak of the CMB radiation. The density in the outer regions of the collapsing cloud is too low for the molecule to be excited to LTE. The $`J`$=4 level population is less affected by the CMB radiation field.

To show a more quantitative measure of the accuracy of the results, the level populations are plotted versus their difference to the mean of all the results (Fig. 6). Only results including the CMB radiation field are used to calculate the mean. At the inner boundary, the solutions diverge into two main groups. This is a result from the inner boundary condition adopted by the different authors. The density and temperature of the test problem were not specified within the inner radius and were either chosen as an empty sphere or a non-rotating sphere of constant density. The different solutions near the inner boundary have a relatively small error (Fig. 6) showing that in this case the inner boundary has little influence on the overall solution and can be ignored in the specification of the problem.

The calculated populations of the $`J`$=1 level show a standard deviation ranging from 2% close to the star down to 1.5% for problem 2b and below 1% in problem 2a. The $`J`$=4 level has a larger spread with a standard deviation of 8% close to the inner boundary up to 12% at $`r`$=2$`\times 10^{17}`$ cm and down to 2% at the outer boundary for the high $`\tau `$ (problem 2b) case. The HCO+ lines have high critical densities and the $`J`$=4 level population is far from LTE. In addition, the optical depth is high,making the calculation for the 4–3 line particularly difficult. In problem 2a the $`J=4`$ level rises from 3% close to the inner boundary up to 4% in the transition region and down to 1% close to the outer boundary. Lacking an analytical solution, we take the average of the numerical results as the best estimate of the exact solution. The deviations from the average therefore give an estimate of the implicit and explicit approximations of the codes. Fitting models to observational data with criteria better than these deviations will not lead to more accurate estimates of the physical parameters.

For all the problems, the level populations were also calculated using LI by Dullemond and Wiesemeyer. The populations agreed within the errors as shown in Figs. 4 & 6 for the low $`\tau `$ cases. However in both the 1b and 2b problems the solution was not reached using similar convergence criteria as the ALI calculations. Even though the LI method itself should in principle find the same solution, the number of iterations and more strict convergence criteria make it a numerically too costly task to solve.

The problem is particularly sensitive to the gridding of the physical parameters. A comparison of two runs, where one had twice as many gridcells, showed a significant decrease of the error compared to the mean value for one of the Monte Carlo code (C). This shows the importance of a correct gridding of the problem.

#### 5.2.2 Excitation temperature

As a second comparison, the excitation temperature (Fig. 7) is plotted for all solutions, which may be a more intuitive manner to show the results. When level populations are in LTE, the excitation temperature equals the local kinetic temperature. Comparing the results with Fig. 2b, the excitation temperature of the 1–0 line is close to LTE near the inner boundary but becomes subthermal at larger radii. The $`J`$=4–3 excitation temperature is always well below its LTE value, which is due to its much higher critical density. As Fig. 7 shows, the excitation temperature never drops below 2.728 K, except in the case where the CMB radiation was ignored. In the inner region ($`r<5\times 10^{16}`$ cm), the excitation of the $`J`$=1 level is dominated by collisions, while in the outer parts, it is dominated by radiation. The transition region between these two extremes shows the largest differences in the results (Fig. 6a).

#### 5.2.3 Line profiles

Figure 8 compares the results in terms of the line profiles. The calculated level populations at each position in the cloud are used to compute the velocity profiles of the selected lines using a program which calculates the sky brightness distribution. To ensure an equivalent emitting mass, the inner sphere is assumed to be empty for all models. The profiles are calculated by constructing a plane through the origin of the cloud perpendicular to the line of sight, with a spatial resolution small enough to sample the physical distributions. A ray-tracing calculation is performed through this plane from $`\mathrm{}`$ to $`+\mathrm{}`$, keeping track of the intensity in the different velocity bins. For this calculation the program SKY was used, part of the RATRAN code (http://talisker.as.arizona.edu/$``$michiel/) By applying a single common code for these calculations, no secondary uncertainties are introduced. The resulting sky brightness distribution is convolved with a beam of 14$`\mathrm{}`$ for the $`J`$=4–3 transition, appropriate for the James Clerk Maxwell Telescope (JCMT) at this frequency. For the $`J`$=1–0 transition, the beam size is taken to be 29$`\mathrm{}`$ representative of the IRAM 30m telescope. The distance is assumed to be 250 pc.

Comparison of Fig. 8a and c shows that increasing the abundance of HCO+ by a factor of 10 changes the $`J`$=1–0 peak intensity by only a factor of 2.5. The $`J`$=1–0 line becomes more self-absorbed, indicating that the line has indeed become more optically thick. The $`J`$=4–3 emission line is optically thin for the low abundance, as seen from its single-peaked profile in Fig. 8b. The more optically thick model shows the characteristic asymmetric line structure, skewed to the blue. The differences in the level populations are visible in the line profiles. In the $`J`$=1–0 line profiles (Fig. 8c) all solutions lie on top of each other in the line wings and only a slight difference can be seen at line center. The error is largest in the region with low velocities represented by the larger errors in the center of the line profile. The integrated intensity profiles differ by a few % for low $`\tau `$ and 6-7% for the high $`\tau `$ case. The error is in general larger for the $`J`$=4 level, which shows up as a larger error in the integrated lines for the $`J`$=4–3 transition. For low $`\tau `$ the error is $``$ 2%, but for high optical depths deviations up to 12% are found for the most outlying case. Although substantial, these differences are generally less than the current calibration uncertainties of the observational data, which are typically 20% to 30%.

The differences in the profiles shown here give a rough indication of the quality of fitting necessary to interpret observational data. It should be noted that some differences are likely due to different implementations of the cloud model. One can fit a model precisely to an observational data set using a particular radiative transfer code; however the results will not necessarily be the same if a different code is used. The derivation of physical parameters is therefore always limited in accuracy to the error bars given by the code itself in comparison to other codes.

## 6 Discussion and conclusion

In spite of the fact that the spread in our solutions lies within the accuracy limits of current-day (sub-)millimeter instruments, it is important to understand how this spread comes about. The first two test cases were defined analytically, and avoided difficulties of gridding as much as possible. The small spread of these results, even at high optical depths, seems to show that in principle the radiative transfer is done correctly by all codes. The larger spread in the HCO+ collapsing cloud problem can therefore be attributed to the problem of gridding, both in space and in frequency. The coarse spatial gridding in the setup of the problem, and the presence of strong velocity gradients makes it very likely that different interpretations of the sub-grid behavior of temperature, density and velocity lead to slightly different results.

The main results of the comparison of radiative transfer codes for molecular lines can be summarized as follows.

* All relevant methods currently used in molecular astrophysics agree to a few % for optical depths up to $`\tau `$4800. At high optical depth and in transition layers between collision-dominated and radiation-dominated excitation, relative differences up to 2 % can arise for well defined “simple” problems. In practice, for more complex models, the relative differences are higher due to gridding and geometry problems, reaching relative differences up to 12 % locally for the high optical depth model described in this paper (problem 2b) . However, in most practical applications in molecular astrophysics the optical depths are lower than in our test problems, and therefore this 12% error can be regarded as an upper limit.
* The error bars on the line profiles are generally much less than 10%, and therefore lie within the calibration errors of typical (sub-)mm observations. Only in the worst case these errors may be comparable to the observational uncertainty.
* The choice of gridding is of extreme importance, and is one of the major causes of the deviations in the problems presented. Special care should be taken when constructing models for specific problems.
* Abundances and excitation temperatures derived from lines which are formed in the transition layer should be interpreted with caution.
* The direct comparison of results of different programs speeds up significantly the debugging process of new programs.

The work presented here is one step in the direction of standardization of radiative transfer computations in molecular rotational lines. Suggestions for future work in this area include: (1) establishment of a database for collisional rate coefficients, accessible to the community through the WWW; (2) a campaign to have inelastic rate coefficients measured or calculated for those molecules for which these data are unavailable, and (3) comparison of 2D radiative transfer codes, with a test problem based on, e.g., the rotating flattened collapse described by Terebey et al. (1984) or the “sheet” models of protostellar collapse by Hartmann et al. (1996). These developments will be vital to interpreting the high-quality data which e.g. HIFI (the heterodyne instrument onboard the Herschel Space Observatory) and ALMA will provide.

###### Acknowledgements.

We are grateful to the anonymous referee for suggestions and comments that improved the paper considerably. We thank Ewine van Dishoeck for useful advice and discussions during this work. The comparison of line radiative transfer codes was started during a workshop in Leiden in May 1999 in the Lorentz Center. Financial support for the workshop from the Lorentz Center, the Netherlands Research School for Astronomy (NOVA), the Leidsch Kerkhoven Bosscha Fund (LKBF), the Netherlands Organization for Scientific Research (NWO) and the UK Particle Physics and Astronomy Research Council (PPARC) is gratefully acknowledged. SDD acknowledges support from the Research Corporation.


format: text

text:
# Exact integral equation for the renormalized Fermi surface

(August 27, 2002)

## Abstract

The true Fermi surface of a fermionic many-body system can be viewed as a fixed point manifold of the renormalization group (RG). Within the framework of the exact functional RG we show that the fixed point condition implies an exact integral equation for the counterterm which is needed for a self-consistent calculation of the Fermi surface. In the simplest approximation, our integral equation reduces to the self-consistent Hartree-Fock equation for the counterterm.

In his authoritative book on interacting Fermi systems Nozières wrote fourty years ago Nozieres64 : “In practice, we shall never try to calculate the Fermi surface, which is much too difficult.” What is the reason for this difficulty? Formally, the Fermi surface of an interacting Fermi system is defined as the set of all wavevectors $`𝐤_F`$ satisfying Luttinger60

$$ϵ_{𝐤_F}\mu +\mathrm{\Sigma }(𝐤_F,i0)=0,$$
(1)
where $`ϵ_𝐤`$ is the energy dispersion in the absence of interactions, $`\mu `$ is the chemical potential, and $`\mathrm{\Sigma }(𝐤,\omega )`$ is the exact self-energy of the interacting system footnotereal . For simplicity we assume an infinite and spin-rotationally invariant system at zero temperature, so that $`\mathrm{\Sigma }(𝐤,\omega )`$ is independent of the spin. Unfortunately, the function $`\mathrm{\Sigma }(𝐤_F,i0)`$ in Eq. (1) is not known a priori, so that the calculation of the true Fermi surface requires the solution of the many-body problem.

For weak interactions, one might try to determine the Fermi surface perturbatively by simply calculating $`\mathrm{\Sigma }(𝐤,i0)`$ in powers of the interaction and substituting the result into Eq. (1). However, in general the perturbation series contains anomalous terms Kohn60  with unphysical singularities, which are generated because the ground state of the non-interacting system evolves into an excited state of the interacting system when the interaction is adiabatically switched on. As discussed by Nozières Nozieres64 , this artificial level crossing can be avoided by introducing counterterms which are determined by the requirement that the Fermi surface remains fixed as the interaction is adiabatically switched on. This intuitive idea can be implemented perturbatively as follows Nozieres64 ; Feldman96 : Suppose we would like to know the true Fermi surface of a system with Hamiltonian $`H=H_0+H_1`$, where $`H_1`$ describes some general two-body interaction and the non-interacting part is given by

$$H_0=\underset{𝐤,\sigma }{}ϵ_𝐤c_{𝐤\sigma }^{}c_{𝐤\sigma }.$$
(2)
Here $`c_{𝐤\sigma }`$ are the usual annihilation operators of fermions with momentum $`𝐤`$ and spin $`\sigma `$. An expansion in powers of $`H_1`$ leads to Feynman diagrams where vertices corresponding to $`H_1`$ are connected by propagators $`G_0(𝐤,\omega )=[\omega ϵ_𝐤+\mu ]^1`$. These are singular for $`\omega =0`$ and $`ϵ_𝐤=\mu `$, which is not the true Fermi surface defined in Eq. (1). If the perturbative expansion is truncated at a finite order, this leads to the unphysical divergencies mentioned above Kohn60 . To avoid these, we add the counterterm $`_{𝐤\sigma }\mathrm{\Sigma }(𝐤_F,i0)c_{𝐤\sigma }^{}c_{𝐤\sigma }`$ to $`H_0`$ and subtract it again from $`H_1`$, writing $`H=H_0^{}+H_1^{}`$, with

$`H_0^{}={\displaystyle \underset{𝐤\sigma }{}}[ϵ_𝐤+\mathrm{\Sigma }(𝐤_F,i0)]c_{𝐤\sigma }^{}c_{𝐤\sigma },`$ (3)
and $`H_1^{}=H_1_{𝐤\sigma }\mathrm{\Sigma }(𝐤_F,i0)c_{𝐤\sigma }^{}c_{𝐤\sigma }`$. Here $`𝐤_F`$ is the wavevector closest to $`𝐤`$ lying on the Fermi surface, see Fig. 1.

Using Eq. (1), the corresponding free propagator can then be written as $`G_0^{}(𝐤,\omega )=[\omega ϵ_𝐤+ϵ_{𝐤_F}]^1`$, which by construction is singular on the true Fermi surface. But how do we find the counterterm necessary to calculate $`ϵ_{𝐤_F}`$? Following the usual strategy adopted in field theory ZinnJustin89 , we may expand the irreducible self-energy associated with the modified interaction $`H_1^{}`$ perturbatively in powers of $`H_1^{}`$ and require that, order by order in perturbation theory, the corrections vanish when we set $`\omega =0`$ and $`𝐤=𝐤_F`$. This renormalized perturbation theory leads to complicated integral equations for $`\mathrm{\Sigma }(𝐤_F,i0)`$, which must be solved numerically Neumayr02 .

A few years ago Anderson Anderson93  critically discussed a simplified version of the renormalized perturbation theory outlined above, where $`\mathrm{\Sigma }(𝐤_F,i0)`$ is replaced by a constant $`\delta \mu `$. He correctly pointed out that in general it is not allowed to ignore the momentum-dependence of the counterterms, and argued that in two dimensions the effective interaction between electrons with opposite spin and wavevectors on the true Fermi surface depends in a subtle way on the boundary conditions which cannot be adequately taken into account perturbatively. To gain a better understanding of this problem, it should be useful to have an algorithm for calculating the Fermi surface which does not rely on the perturbative expansion of counterterms in powers of the interaction. In this work we show that such an algorithm follows in a straightforward way from the exact functional renormalization group (RG) approach described in Ref.Kopietz01 .

In the past decade several authors have used Wilsonian RG methods to study interacting Fermi systems Kopietz01 ; Benfatto90 ; Polchinski92 ; Shankar94 ; Zanchi96 ; Salmhofer98 ; Halboth00 ; Binz02 , using different versions of the RG. In particular, the RG transformations used in Refs. Zanchi96 ; Salmhofer98 ; Halboth00 ; Binz02  exclusively focus on the mode elimination step. Although such a procedure is sufficient if one considers only the one-loop flow of marginal couplings, in general a complete Wilsonian RG transformation consists not only of the mode elimination, but includes also the rescaling of momenta, frequencies and fields Ma76 . While the field rescaling (i.e. the wavefunction renormalization) becomes only important beyond the one-loop approximation, the RG flow of all relevant and irrelevant couplings is determined in an essential way by the rescaling step even at the one-loop level. We emphasize that the rescaling of momenta and frequencies is more than a trivial mathematical change of variables – it is crucial to detect possible fixed points of the RG and to calculate critical exponents Ma76 . Moreover, as will be shown below, the length and time rescaling is very helpful for a derivation of the self-consistent Hartree-Fock approximation within the framework of the exact RG. The importance of rescaling in the presence of a Fermi surface has been emphasized by Polchinski Polchinski92  and by Shankar Shankar94 . In Ref.Kopietz01  we have shown how the exact functional RG approach developed previously Zanchi96 ; Salmhofer98 ; Halboth00  can be modified to include the rescaling step. Here, we shall show that the rescaling is crucial to obtain a non-perturbative definition of the Fermi surface as a fixed point manifold of the RG. We obtain an explicit algorithm for calculating the Fermi surface which does not rely on the iterative procedure of fixing the counterterms order by order in perturbation theory. The fixed point property of the Fermi surface has also been emphasized by Ferraz Ferraz02 , who recently discussed the Fermi surface renormalization in a special two-dimensional system using the field theoretical RG.

The exact functional RG can be formulated in terms of an infinite hierarchy of coupled differential equations for the irreducible $`2n`$-point vertices $`\mathrm{\Gamma }_\xi ^{(2n)}(K_1^{},\mathrm{},K_n^{};K_n,\mathrm{},K_1)`$, where $`K=(\sigma ,𝐤,i\omega )`$ is a collective label for spin projection $`\sigma `$, momentum $`𝐤`$, Matsubara frequency $`i\omega `$. Here $`\xi `$ is an infrared cutoff with units of energy which regularizes the singularity of the free propagator. To scale wavevectors toward the Fermi surface, it is convenient to perform a non-linear coordinate transformation in momentum space, $`𝐤=\widehat{𝐧}k_F(\widehat{𝐧})+\widehat{𝐯}_F\xi q/|𝐯_F|`$, and eliminate $`𝐤`$ in favor of the dimensionless variable $`q`$ and the unit vector $`\widehat{𝐧}`$. Here $`k_F(\widehat{𝐧})`$ is the length of $`𝐤_F`$ parameterized by $`\widehat{𝐧}`$, and $`𝐯_F=_𝐤ϵ_𝐤|_{𝐤_F}`$ is the Fermi velocity of the non-interacting system at the true $`𝐤_F`$, see Fig. 1. The corresponding unit vector is denoted by $`\widehat{𝐯}_F=𝐯_F/|𝐯_F|`$. Geometrically, $`q=𝐯_F𝐩/\xi =𝐯_F(𝐤\widehat{𝐧}k_F(\widehat{𝐧}))/\xi `$ measures the distance of a given $`𝐤`$ from the Fermi surface in units of the infrared cutoff. We also define rescaled frequencies $`ϵ=\omega /\xi `$, and label the degrees of freedom by $`Q=(\sigma ,\widehat{𝐧},q,iϵ)`$ instead of $`K`$. To implement the scaling toward the Fermi surface Polchinski92  we consider the RG flow of the following rescaled vertices Kopietz01

$`\stackrel{~}{\mathrm{\Gamma }}_t^{(2n)}(Q_1^{},\mathrm{},Q_n^{};Q_n,\mathrm{},Q_1)=`$ (4)
$`\nu _0^{n1}\xi ^{n2}\left[Z_t^{\widehat{𝐧}_1^{}}\mathrm{}Z_t^{\widehat{𝐧}_n^{}}Z_t^{\widehat{𝐧}_n}\mathrm{}Z_t^{\widehat{𝐧}_1}\right]^{1/2}`$
$`\times \mathrm{\Gamma }_\xi ^{(2n)}(K_1^{},\mathrm{},K_n^{};K_n,\mathrm{},K_1).`$
Here $`t=\mathrm{ln}(\xi /\xi _0)`$ is a logarithmic flow parameter, $`\nu _0=\frac{d𝐤}{(2\pi )^D}\delta (ϵ_𝐤ϵ_{𝐤_F})`$ is the density of states (per spin projection) of the non-interacting system at the true Fermi surface, and

$$Z_t^{\widehat{𝐧}}=\left[1\frac{\mathrm{\Sigma }_\xi (𝐤_F,\omega +i0)}{\omega }|_{\omega =0}\right]^1$$
(5)
is the wavefunction renormalization factor, where $`\mathrm{\Sigma }_\xi (𝐤,i\omega )`$ is the irreducible self-energy of the system with infrared cutoff $`\xi =\xi _0e^t`$. To determine the true Fermi surface self-consistently it turns out to be useful to subtract the counterterm $`\mathrm{\Sigma }(𝐤_F,i0)`$ from the irreducible two-point vertexKopietz01 , defining

$$\stackrel{~}{\mathrm{\Gamma }}_t^{(2)}(Q)=\frac{Z_t^{\widehat{𝐧}}}{\xi }\mathrm{\Gamma }_\xi ^{(2)}(K)=\frac{Z_t^{\widehat{𝐧}}}{\xi }\left[\mathrm{\Sigma }_\xi (K)\mathrm{\Sigma }(𝐤_F,i0)\right].$$
(6)
The subtracted two-point vertex satisfies the exact flow equation Kopietz01

$`_t\stackrel{~}{\mathrm{\Gamma }}_t^{(2)}(Q)`$ $`=`$ $`(1\eta _t^{\widehat{𝐧}}q_qϵ_ϵ)\stackrel{~}{\mathrm{\Gamma }}_t^{(2)}(Q)`$ (7)
$``$ $`{\displaystyle _Q^{}}\dot{G}_t(Q^{})\stackrel{~}{\mathrm{\Gamma }}_t^{(4)}(Q,Q^{};Q^{},Q),`$
where $`\eta _t^{\widehat{𝐧}}=_t\mathrm{ln}Z_t^{\widehat{𝐧}}`$ is the flowing anomalous dimension (which vanishes for large $`t`$ if the system is a Fermi liquid). The integration measure in Eq. (7) is defined by

$$_Q=\underset{\sigma }{}\frac{dS_{\widehat{𝐧}}}{S_D}𝑑qJ(\widehat{𝐧},q)\frac{dϵ}{2\pi },$$
(8)
where $`dS_{\widehat{𝐧}}`$ is a surface element and $`S_D`$ is the surface area of the unit sphere in $`D`$ dimensions, and $`J(\widehat{𝐧},q)`$ is the Jacobian associated with the transformation $`𝐤(\widehat{𝐧},q)`$ divided by $`\nu _0\xi ^2`$. The function $`\dot{G}_t(Q)`$ in Eq. (7) is defined by

$$\dot{G}_t(Q)=\frac{\delta (\stackrel{~}{\mathrm{\Omega }}_Q1)}{Z_t^{\widehat{𝐧}}\left[iϵ\xi _t^{\widehat{𝐧}}(q)\right]+\stackrel{~}{\mathrm{\Gamma }}_t^{(2)}(Q)},$$
(9)
where $`\xi _t^{\widehat{𝐧}}(q)=(ϵ_{𝐤_t}ϵ_{𝐤_F})/\xi `$, with $`𝐤_t=\widehat{𝐧}k_F(\widehat{𝐧})+\widehat{𝐯}_F\xi _0e^tq/|𝐯_F|`$. Here $`\stackrel{~}{\mathrm{\Omega }}_Q`$ is some function that measures the distance from the Fermi surface, for example $`\stackrel{~}{\mathrm{\Omega }}_Q=|\xi _t^{\widehat{𝐧}}(q)||q|`$. The flow equation of the rescaled irreducible four-point vertex $`\stackrel{~}{\mathrm{\Gamma }}_t^{(4)}(Q_1^{},Q_2^{};Q_2,Q_1)`$ is not explicitly needed in this work; it can be found in Ref.Kopietz01 .

The shape of the Fermi surface of the interacting system is determined by the RG flow of the couplings

$$\stackrel{~}{\mu }_t^{\widehat{𝐧}}\stackrel{~}{\mathrm{\Gamma }}_t^{(2)}(\sigma ,\widehat{𝐧},q=0,iϵ=i0).$$
(10)
From Eq. (7) we see that these couplings satisfy the exact flow equation

$$_t\stackrel{~}{\mu }_t^{\widehat{𝐧}}=(1\eta _t^{\widehat{𝐧}})\stackrel{~}{\mu }_t^{\widehat{𝐧}}+\dot{\mathrm{\Gamma }}_t^{(2)}(\widehat{𝐧}),$$
(11)
with

$$\dot{\mathrm{\Gamma }}_t^{(2)}(\widehat{𝐧})=_Q^{}\dot{G}_t(Q^{})\stackrel{~}{\mathrm{\Gamma }}_t^{(4)}(Q_0,Q^{};Q^{},Q_0),$$
(12)
where $`Q_0=(\sigma ,\widehat{𝐧},q=0,iϵ=i0)`$. Note that in dimensions $`D>1`$ there are infinitely many couplings $`\stackrel{~}{\mu }_t^{\widehat{𝐧}}`$, labelled by the unit vector $`\widehat{𝐧}`$. Obviously, each $`\stackrel{~}{\mu }_t^{\widehat{𝐧}}`$ with $`\eta _{\mathrm{}}^{\widehat{𝐧}}=lim_t\mathrm{}\eta _t^{\widehat{𝐧}}<1`$ is relevant, so that some fine tuning of the bare couplings $`\stackrel{~}{\mu }_0^{\widehat{𝐧}}`$ is necessary to force $`\stackrel{~}{\mu }_t^{\widehat{𝐧}}`$ to flow into a fixed point of the RG. Because a finite limit $`\stackrel{~}{\mu }_{\mathrm{}}^{\widehat{𝐧}}=lim_t\mathrm{}\stackrel{~}{\mu }_t^{\widehat{𝐧}}`$ means that we have found the true Fermi surface of the interacting system Kopietz01 , we conclude that the detailed shape of the Fermi surface is sensitive to the numerical values of the bare couplings of the theory. Polchinski Polchinski92  pointed out that such a fine tuning of the bare couplings is unnatural, because physical effects which depend on the precise shape of the Fermi surface are not protected against small perturbations.

To further elucidate the relation between the relevant couplings $`\stackrel{~}{\mu }_t^{\widehat{𝐧}}`$ and the shape of the Fermi surface, let us explicitly derive from Eq. (11) a self-consistency condition for the Fermi surface. Therefore it is useful to transform Eq. (11) into an integral equation,

$$\stackrel{~}{\mu }_t^{\widehat{𝐧}}=e^{t_0^t𝑑\tau \eta _\tau ^{\widehat{𝐧}}}\left[\stackrel{~}{\mu }_0^{\widehat{𝐧}}+_0^t𝑑t^{}e^{t^{}+_0^t^{}𝑑\tau \eta _\tau ^{\widehat{𝐧}}}\dot{\mathrm{\Gamma }}_t^{}^{(2)}(\widehat{𝐧})\right].$$
(13)
Suppose now that we have adjusted the bare couplings such that for $`t\mathrm{}`$ the flowing couplings $`\stackrel{~}{\mu }_t^{\widehat{𝐧}}`$ indeed approach finite fixed point values. Assuming that the associated anomalous dimensions $`\eta _{\mathrm{}}^{\widehat{𝐧}}`$ are smaller than unityfootnoteeta , we conclude from Eq. (13) that the limit $`\stackrel{~}{\mu }_{\mathrm{}}^{\widehat{𝐧}}=lim_t\mathrm{}\stackrel{~}{\mu }_t^{\widehat{𝐧}}`$ can only be finite if the initial values $`\stackrel{~}{\mu }_0^{\widehat{𝐧}}`$ are chosen such that

$`\stackrel{~}{\mu }_0^{\widehat{𝐧}}`$ $`=`$ $`{\displaystyle _0^{\mathrm{}}}𝑑te^{t+_0^t𝑑\tau \eta _\tau ^{\widehat{𝐧}}}\dot{\mathrm{\Gamma }}_t^{(2)}(\widehat{𝐧})`$
$`=`$ $`{\displaystyle _0^{\mathrm{}}}𝑑te^{t+_0^t𝑑\tau \eta _\tau ^{\widehat{𝐧}}}{\displaystyle _Q^{}}\dot{G}_t(Q^{})\stackrel{~}{\mathrm{\Gamma }}_t^{(4)}(Q_0,Q^{};Q^{},Q_0).`$
This is an implicit equation for $`\stackrel{~}{\mu }_0^{\widehat{𝐧}}`$, relating it to the values of the two-point vertex and the four-point vertex on the entire RG trajectory. Keeping in mind that the right-hand side of Eq. (LABEL:eq:integral2) implicitly depends on $`\stackrel{~}{\mu }_t^{\widehat{𝐧}}`$ and that according to Eq. (6)

$$\mathrm{\Sigma }(𝐤_F,i0)\mathrm{\Sigma }_{\xi _0}(𝐤_F,i0)=\xi _0\stackrel{~}{\mu }_0^{\widehat{𝐧}}/Z_0^{\widehat{𝐧}},$$
(15)
it is obvious that Eq. (LABEL:eq:integral2) can be regarded as an integral equation for the counterterm $`\mathrm{\Sigma }(𝐤_F,i0)`$, the solution of which yields the true shape of the Fermi surface. At this point it is instructive to transform Eq. (LABEL:eq:integral2) back to unrescaled variables, choosing for simplicity the initial conditions $`\mathrm{\Sigma }_{\xi _0}(𝐤_F,i0)=0`$ and $`Z_0^{\widehat{𝐧}}=1`$. Using the above definitions we find that Eq. (LABEL:eq:integral2) is equivalent with

$`\mathrm{\Sigma }(𝐤_F,i0)`$ $`=`$ $`{\displaystyle \underset{\sigma ^{}}{}}{\displaystyle \frac{d𝐤^{}}{(2\pi )^D}\frac{d\omega ^{}}{2\pi }\frac{\mathrm{\Theta }(\xi _0\xi _𝐤^{})}{i\omega ^{}ϵ_𝐤^{}+\mu \mathrm{\Sigma }_{\xi _𝐤^{}}(K^{})}}`$ (16)
$`\times \mathrm{\Gamma }_{\xi _𝐤^{}}^{(4)}(𝐤_F,i0,\sigma ,𝐤^{},i\omega ^{},\sigma ^{};𝐤^{},i\omega ^{},\sigma ^{},𝐤_F,i0,\sigma ),`$
where $`\xi _𝐤^{}=|ϵ_𝐤^{}ϵ_{𝐤_F^{}}|`$. Note that the right-hand side of Eq. (16) involves the flowing self-energy and four-point vertex at the scales $`\xi =\xi _𝐤^{}`$ which depend on the distance from the true Fermi surface. We emphasize that the exact integral equation (16) and the equivalent rescaled equation (LABEL:eq:integral2) fix the counterterm $`\mathrm{\Sigma }(𝐤_F,i0)`$ from the requirement that for $`t\mathrm{}`$ all couplings approach finite fixed point values.

If the system is a Luttinger liquid, then the rescaled version (LABEL:eq:integral2) of the self-consistency equation is more convenient than the unrescaled version (16), because for a Luttinger liquid the marginal part of the four-point vertex vertex $`\mathrm{\Gamma }_\xi ^{(4)}`$ without wavefunction renormalization factors diverges for $`\xi 0`$, while the rescaled four-point vertex $`\stackrel{~}{\mathrm{\Gamma }}_t^{(4)}`$ defined in Eq.(4) approaches for $`t\mathrm{}`$ a finite limit. In this case the divergence of the unrescaled $`\mathrm{\Gamma }_\xi ^{(4)}`$ is canceled by the vanishing wavefunction renormalization factors at the Luttinger liquid fixed point Busche02 .

Given a solution $`\mathrm{\Sigma }(𝐤_F,i0)`$ of Eq. (16), we may calculate the compressibility of the system by substituting the result for $`\mathrm{\Sigma }(𝐤_F,i0)`$ into Eq. (1) and solving for $`k_F(\widehat{𝐧})`$, which implicitly depends on the chemical potential $`\mu `$. According to the Luttinger theorem Luttinger60  the density $`n(\mu )`$ of the system is determined by the volume enclosed by the Fermi surface,

$$n(\mu )=\underset{\sigma }{}\frac{d𝐤}{(2\pi )^D}\mathrm{\Theta }(k_F(\widehat{𝐧})|𝐤|),$$
(17)
so that we obtain for the compressibility $`\chi _n`$

$$n^2\chi _n=\frac{n}{\mu }=\underset{\sigma }{}\frac{dS_{\widehat{𝐧}}}{(2\pi )^D}k_F^{D1}(\widehat{𝐧})\frac{k_F(\widehat{𝐧})}{\mu }.$$
(18)
Note that the compressibility is a functional of the shape of the true Fermi surface of the many-body system, i.e. the compressibility characterizes the RG fixed point. Because the rescaling of momenta and frequencies is crucial to obtain fixed points of the RG Ma76 , the compressibility and other uniform susceptibilities are not so easy to calculate within alternative versions of the RG which omit the rescaling Salmhofer98 ; Halboth00 ; Binz02 .

From our point of view, the Fermi surface is a fixed point manifold of the RG, so that it is meaningless to talk about the RG flow of the Fermi surface. However, it is possible to define a “flowing Fermi surface” $`𝐤_{F,t}`$ via Kopietz01

$$ϵ_{𝐤_{F,t}}\mu +\mathrm{\Sigma }_{\xi _0e^t}(𝐤_{F,t},i0)=0,$$
(19)
which by construction approaches the true Fermi surface for $`t\mathrm{}`$. If we choose the initial conditions at $`t=0`$ (corresponding to $`\xi =\xi _0`$) such that $`\mathrm{\Sigma }_{\xi _0}(𝐤_{F,0},i0)=0`$, then $`𝐤_{F,0}`$ is the Fermi surface of the non-interacting system at the same chemical potential as the interacting system. This corresponds to a non-interacting system at the density $`n_0(\mu )=_\sigma \frac{d𝐤}{(2\pi )^D}\mathrm{\Theta }(k_{F,0}(\widehat{𝐧})|𝐤|)`$, which is in general different from the density $`n(\mu )`$ of the interacting system given in Eq. (17). Note that in Fermi liquid theory one usually works at constant density. However, as discussed by Nozières Nozieres237 , for a self-consistent calculation of the Fermi surface it is more convenient to determine the counterterm at constant chemical potential $`\mu `$, and calculate the corresponding density afterwards footnotedens . The practical advantages of such a procedure have been recognized previously in Ref. Gonzalez00 . Moreover, in the field-theoretical approach advanced by Ferraz Ferraz02  the chemical potential is also a RG invariant.

In the simplest approximation Eqs. (LABEL:eq:integral2) and (16) reduce to the Hartree-Fock self-consistency equation for the counterterm $`\mathrm{\Sigma }(𝐤_F,i0)`$. To see this, we approximate the flowing four-point vertex in Eq. (16) as follows

$`\mathrm{\Gamma }_{\xi _𝐤^{}}^{(4)}(𝐤_F,i0,\sigma ,𝐤^{},i\omega ^{},\sigma ^{};𝐤^{},i\omega ^{},\sigma ^{},𝐤_F,i0,\sigma )`$ (20)
$`\mathrm{\Gamma }_{\xi =0}^{(4)}(𝐤_F,i0,\sigma ,𝐤_F^{},i0,\sigma ^{};𝐤_F^{},i0,\sigma ^{},𝐤_F,i0,\sigma )`$
$`\mathrm{\Gamma }_0^{(4)}(𝐤_F,\sigma ;𝐤_F^{},\sigma ^{}),`$
i.e. we project all momenta onto the Fermi surface, ignore the frequency dependence, and replace the flowing vertex by its fixed point value. Moreover, at this level of approximation we may also replace $`\mathrm{\Sigma }_{\xi _𝐤^{}}(K^{})\mathrm{\Sigma }(𝐤_F^{},i0)`$ on the right-hand side of Eq. (16). For $`\xi _0\mathrm{}`$ we then obtain the Hartree-Fock self-consistency equation

$`\mathrm{\Sigma }(𝐤_F,i0)`$ $`=`$ $`{\displaystyle \underset{\sigma ^{}}{}}{\displaystyle \frac{d𝐤^{}}{(2\pi )^D}\mathrm{\Gamma }_0^{(4)}(𝐤_F,\sigma ;𝐤_F^{},\sigma ^{})}`$ (21)
$`\times `$ $`\mathrm{\Theta }\left(\mu ϵ_𝐤^{}\mathrm{\Sigma }(𝐤_F^{},i0)\right).`$
For a given interaction vertex $`\mathrm{\Gamma }_0^{(4)}(𝐤_F,\sigma ;𝐤_F^{},\sigma ^{})`$, Eq. (21) can be used to study possible Fermi surface instabilities such as the Pomeranchuk instability Pomenanchuk58  at the Hartree-Fock level, which should always be the first step before more elaborate methods are used. A trivial generalization of Eq. (21) with a spin-dependent counterterm $`\mathrm{\Sigma }(𝐤_F,i0,\sigma )`$ leads to the self-consistent Hartree-Fock equation for spontaneous ferromagnetism, implying the usual Stoner instability for strong enough interactions in the spin-channel. More generally, if we allow for other types of symmetry breaking, within the same approximations as above the exact RG fixed point equation for the two-point vertex can be reduced to the Hartree-Fock self-consistency equation for the corresponding order parameter. Note that no approximation has been made in deriving Eqs. (LABEL:eq:integral2) and (16), so that these exact fixed point equations (or generalizations thereof for other types of symmetry breaking) can serve as a starting point for a systematic calculation of corrections to the Hartree-Fock approximation.

In summary, in this work we have shown how the true Fermi surface of an interacting Fermi system can be defined self-consistently as a fixed point property of the RG. Our main result are the two equivalent integral equations (LABEL:eq:integral2) and (16), which determine the counterterm $`\mathrm{\Sigma }(𝐤_F,i0)`$ necessary to calculate the Fermi surface from Eq. (1). In the simplest approximation, Eq. (16) reduces to the Hartree-Fock self-consistency condition for the counterterm. However, systematic improvements are possible.

Very recently Dusuel and Douçot Dusuel02  presented a detailed analysis of the Fermi surface deformations in quasi one-dimensional electronic systems, using perturbation theory and the RG method. They realized that some “slight modification” of the Wilson-Polchinski RG approach is necessary in order to use this approach for a self-consistent calculation of the Fermi surface, but admitted that a practical implementation of such a modification remains to be attempted. We have shown here that the necessary modification of the functional RG used in Refs. Zanchi96 ; Salmhofer98 ; Halboth00 ; Binz02  is simply the usual rescaling step Kopietz01 , which is an essential part of the “orthodox” Wilsonian RG Ma76 . We conclude that the functional RG approach for interacting Fermi systems in the form advanced in Ref. Kopietz01  provides an elegant solution to the difficult Nozieres64  problem of self-consistently constructing the true Fermi surface.

We thank Lorenz Bartosch, Tom Busche, Alvaro Ferraz, Volker Meden, Walter Metzner for discussions. This work was partially supported by the DFG via Forschergruppe FOR 412, Project No. KO 1442/5-1.


text:
# Another simple proof of a theorem of Chandler Davis

## Abstract.

In 1957, Chandler Davis proved that unitarily invariant convex functions on the space of hermitian matrices are precisely those which are convex and symmetrically invariant on the set of diagonal matrices. We give a simple perturbation theoretic proof of this result. (Davis’ argument was also very short, though based on completely different ideas).

###### Key words and phrases:

convexity, symmetric, hermitian, orthogonally invariant  The author was partially supported by the Department of Mathematical Sciences of the National Science Foundation. He would like to thank Omar Hijab for helpful conversations.

Consider an orthogonally invariant function $`f`$ defined on the set of $`n\times n`$ symmetric matrices. Such a function has to factor through the spectrum:

(1)
$$f(M)=g\lambda (M),$$

where $`g`$ is a symmetric function:

(2)
$$g(\lambda _1,\mathrm{},\lambda _n)=g(\lambda _{\sigma (1)},\mathrm{},\lambda _{\sigma (n)}),$$

for any permutation $`\sigma .`$

In the sequel we shall further assume that $`f`$ is a $`C^2`$ convex function, and under this assumption we shall show that such functions are precisely those decomposing as per Eq. (1), with convex $`g.`$ The argument and the statement are identical for unitarily invariant functions of Hermitian matrices; in that setting the theorem was proved in \[Davis57\], by a completely different argument (Davis made no regularity assumption, but this is easily dispensed with (see Remark 0.2)).

To show this, let $`M=P+tQ,`$ and let $`\stackrel{~}{f}_{P,Q}(t)=f(M).`$ It is enough to show that for any symmetric $`P,Q,`$

$$\frac{d^2\stackrel{~}{f}_{P,Q}}{dt^2}(0)>0.$$

We compute (dropping the subscript, since from now on $`P,Q`$ do not vary):

(3)
$$\frac{d\stackrel{~}{f}}{dt}=\underset{i=1}{\overset{n}{}}\frac{g}{\lambda _i}\dot{\lambda _i}.$$

(4)
$$\frac{d^2\stackrel{~}{f}}{dt^2}=\underset{1i,jn}{}\frac{^2g}{\lambda _i\lambda _j}\dot{\lambda }_i\dot{\lambda _j}+\underset{i=1}{\overset{n}{}}\frac{g}{\lambda _i}\ddot{\lambda }_i.$$

The first sum is positive, since it equals

$$\dot{\lambda }^tH(g)\dot{\lambda },$$

and the Hessian $`H(g)`$ is positive definite by assumption. It now suffices to show that the second sum (which can be written as $`g\ddot{\lambda }`$) is non-negative.

By continuity, it is sufficient to prove this for $`P`$ whose spectrum is simple. By orthogonal invariance, we can compute in a basis where $`P`$ is diagonal. According to \[Kato80, page 81\], in that case

$$\ddot{\lambda }_i=\underset{ji}{}(\lambda _j\lambda _i)^1Q_{ij}^2,$$

and so

$$f\ddot{\lambda }=\underset{j>i}{}Q_{ij}^2\frac{\frac{g}{\lambda _j}\frac{g}{\lambda _i}}{\lambda _j\lambda _i}.$$

The result then follows from the Lemma below.

###### Lemma 0.1.

Let $`f`$ be a convex function such that $`f(x,y)=f(y,x),`$ where $`yx.`$ Then

$$\frac{\frac{f}{x}\frac{f}{y}}{xy}>0.$$

###### Proof.

The conclusion of the lemma is obviously equivalent to:

$$\left(\frac{f}{x}\frac{f}{y}\right)\left(xy\right)>0.$$

To show this, consider

$$h(t)=f((1t)x+ty,tx+(1t)y).$$

The function $`h(t)`$ is convex, and $`h(0)=h(1).`$ This obviously implies that $`\dot{h}(0)0.`$ Since

$$\dot{h}(0)=\left(\frac{f}{x}\frac{f}{y}\right)(yx),$$

the conclusion follows. ∎

###### Remark 0.2.

Since smooth symmetric convex functions are dense in the set of all symmetric functions, we have actually shown the result for all convex symmetric functions, as follows:

* We start with a symmetric locally integrable function $`g`$ on the set of diagonal matrices.
* The function $`g`$ can be extended (by orthogonal invariance) to the set of all symmetric matrices, to obtain a function $`f.`$
* Convolving $`f`$ with a $`C^{\mathrm{}}`$ Gaussian approximation $`D_n`$ to the delta function, then averaging over the orthogonal orbits we obtain a convex orthogonally invariant function $`f_n,`$ such that the $`ff_n1/n,`$ where $``$ denotes the $`sup`$-norm. The restriction of $`f_n`$ to the diagonal matrices is smooth and convex (since convolving with a positive kernel and averaging both preserve convexity) function $`g_n.`$
* We apply the smooth argument above to $`g_n`$ and $`f_n,`$ to show that $`f_n`$ is a convex function.
* Since $`f`$ is a sup-norm limit of convex functions $`f_n,`$ the result follows.

###### Remark 0.3.

The computation above can also be extended to certain infinite dimensional situations. In particular, it is not hard to see that if $`P`$ has a compact resolvent and $`Q`$ is *bounded* then if the function $`g`$ is the the spectral zeta function or the logarithm of the determinant of the operator. In particular, the regularized $`\mathrm{log}det`$ of the Laplacian on a Riemann surface is convex with respect to bounded perturbations, as are special values of the spectral zeta function. This is very much *not* true for arbitrary perturbations: in particular, when $`P=Q=\mathrm{\Delta },`$ the convexity h$`\mathrm{log}det`$ depends on the sign of the spectral zeta function at $`0,`$ which is the same as the sign of the Euler characteristic, so that $`\mathrm{log}det`$ is *not* convex with respect to the “diagonal” perturbation of the Laplacian whenever the Riemann surfvace has negative Euler characteristic.


text:
# Size and Location Control of Si Nanocrystals at Ion Beam Synthesis in Thin SiO2 Films

thanks: Tel. +49 351 260 3148; Fax. +49 351 260 3285

\]http://www.fz-rossendorf.de

(March 0 d , 2024)

## Abstract

Binary collision simulations of high-fluence $`1\text{ }\mathrm{keV}`$ Si+ ion implantation into $`8\text{ }\mathrm{nm}`$ thick SiO2 films on (001)Si were combined with kinetic Monte Carlo simulations of Si nanocrystal (NC) formation by phase separation during annealing. For nonvolatile memory applications, these simulations help to control size and location of NCs. For low concentrations of implanted Si, NCs form via nucleation, growth and Ostwald ripening, whereas for high concentrations Si separates by spinodal decomposition. In both regimes, NCs form above a thin NC free oxide layer at the SiO2/Si interface. This, self-adjusted layer has just a thickness appropriate for NC charging by direct electron tunneling. Only in the nucleation regime the width of the tunneling oxide and the mean NC diameter remain constant during a long annealing period. This behavior originates from the competition of Ostwald ripening and Si loss to the Si/SiO2 interface. The process simulations predict that, for nonvolatile memories, the technological demands on NC synthesis are fulfilled best in the nucleation regime.

nanocrystals, phase separation, nonvolatile memory, nucleation, spinodal decomposition

Recently, nonvolatile memory concepts based on nanocrystals (NCs) embedded in the gate oxide of MOS transistors have attracted much interest Tiwari et al. (1996). For that aim, NCs have been synthesized by a variety of techniques like chemical vapor deposition Takahashi et al. (2000), ion implantation Normand et al. (2001); von Borany et al. (1999), and Si aerosol deposition Ostraat et al. (2001). Ion implantation followed by thermally activated precipitation of the implanted impurity atoms is most compatible with current silicon technology. By low-energy Si+ ion implantation into thin SiO2 layers on (001)Si, NCs of Si were formed a few nanometers above the Si/SiO2 interface Normand et al. (2001). This allows charging of the NCs by direct electron tunneling, which is a prerequisite for high endurance and low operation voltages De Salvo et al. (2001). Further optimization of location and size of ion beam synthesized NCs for memory application requires a deeper understanding of the mechanisms involved, which determine (i) the built-up of Si supersaturation by high-fluence ion implantation and (ii) NC formation by phase separation. For that aim, process simulations were divided into two steps. The Si implantation was studied using the binary collision code TRIDYN Möller and Eckstein (1984), which includes dynamic target changes. The phase separation of Si from SiO2 during subsequent annealing has been simulated with a kinetic lattice Monte-Carlo code, which describes the thermally activated processes.

The TRIDYN depth profiles are shown in Fig. 1 for $`1\text{ }\mathrm{keV}`$ Si+ ion implantation into SiO2. TRIDYN takes into account dynamic target changes due to ion deposition, ion erosion and ion beam mixing. The input parameters required by the simulation include the displacement and surface binding energies of target atoms. The displacement energies $`E_d`$ for both, Si and O, were assumed to be $`8\text{ }\mathrm{eV}`$. This value proved to yield satisfactory agreement between earlier TRIDYN simulations of ion mixing and experiments Möller (1986) and is also consistent with the choice of Sigmund and Gras-Marti Sigmund and Gras-Marti (1981) in their theoretical treatment of ion mixing. For the present problem, simulations with $`E_d`$ varied by factors of 0.5 and 2 did not show any significant differences of the Si deposition profiles. The surface binding energies of Si and O are assumed to vary linearly with the surface composition in a way that they balance the enthalpies of sublimation and decomposition of Si and SiO2, respectively Möller and Posselt (2001); Kelly (1980).

The Si profile broadens by ion beam mixing, sputtering and swelling. For comparison, a Si implantation profile calculated by TRIM Ziegler et al. (1985) has been added to Fig. 1, which is much sharper than the corresponding TRIDYN profile. Accordingly, the Si peak concentration is significantly overestimated by TRIM. Additionally, sputtering has led to a Si enrichment at the target surface.

During subsequent annealing, Si implanted into the thin SiO2 separates from the oxide phase. In general, the process of phase separation is expected to be a sequence of physical mechanisms like nucleation, growth, and Ostwald ripening of Si precipitates or, at higher Si concentrations, spinodal decomposition and interfacial energy minimization of the Si/SiO2 mixture, respectively. However, these mechanisms are the result of a variety of elementary events (like bond breaking, diffusional jumps of atoms, chemical reactions etc.) that occur in random sequence. Here, they are studied by a kinetic 3D lattice Monte Carlo (KMC) method, which is discussed elsewhere in detail Strobel et al. (2001).

The kinetics of Si atoms is described in a solid host matrix (SiO2), which is the background or ”system’s vacuum”. Thereby, an underlying fcc lattice has been assumed, which is the most isotropic lattice. Within this host dissolved Si diffuses and can form precipitates. (The lattice spacing was chosen such that the correct atomic Si density is obtained.) Applying the classical lattice gas model with attractive Si-Si interaction, the energetics is determined by the nearest-neighbor Ising model. The Metropolis algorithm Metropolis et al. (1953) is used to describe the kinetics of the system.

Si dissolved in the matrix performs diffusional jumps from one lattice site to another with the probability $`1/\tau \mathrm{exp}\left\{E_A/k_BT\right\}`$, where $`E_A`$ is the activation energy of diffusion, $`1/\tau `$ is the attempt frequency and $`k_BT`$ has its usual meaning. Due to a time scale normalization, the internal time unit of the KMC simulation is a Monte Carlo step (MCS) given by $`\tau \mathrm{exp}\{E_A/k_BT\}`$, thus depending on the temperature. Si dimers, trimers and larger agglomerates can form and exhibit a binding energy, which is the product of the Si-Si bond strength $`E_B`$ by the number of Si-Si bonds. A Si atom having $`n_i`$ Si neighbors jumps to an empty neighboring site having $`n_f`$ Si neighbors with a reduced probability $`1/\tau \mathrm{exp}\left\{\left[E_A+(n_in_f)E_B\right]/k_BT\right\}`$ if $`n_f>n_i`$. Otherwise, the diffusional jump probability remains valid.

In principle, the bond strength $`E_B`$ for our KMC simulation can be determined from the solubility of Si in SiO2 via the detailed balance of Si attachment/detachment at the Si/SiO2 interface. (It has to be taken into account that the coordination number in the fcc lattice is 12 instead of 4 in the case of Si lattice.) However, the diffusivity and the solubility of Si in SiO2 are largely unknown. Thus, a direct relation of simulation and experimental temperatures is difficult, which holds for the time scale too. Nevertheless, the path of systems’s evolution towards equilibrium and the regimes of phase separation predicted by KMC simulations may improve the process understanding substantially.

In Fig. 2, snapshots of KMC simulations of phase separation are shown for $`E_B/k_BT=2.0`$. Thereby, the TRIDYN profiles of Fig. 1 were used as initial Si distributions. The size of the simulation cell is $`56\times 56\times \mathrm{\Delta }z\text{ }\mathrm{nm}^3`$, where $`\mathrm{\Delta }z`$ is $`8\text{ }\mathrm{nm}`$ plus swelling due to implantation. Besides the image at the lower right corner, a quarter of the simulation cell is shown only. The simulation cell borders on a fixed (001) layer of the Si substrate to account for the Si/SiO2 interface. Periodic boundary conditions were applied in the oxide plane, whereas reflecting boundary conditions were assumed at the surface.

Two regimes of phase separation are predicted by the KMC simulations. A ”nucleation and growth” regime is observed for $`2\times 10^{15}\text{ }\mathrm{Si}^+\mathrm{cm}^2`$. Si NCs form by nucleation and grow further at the expense of Si supersaturation (Fig. 2 (a), 100 kMCS). Later on, NCs grow by Ostwald ripening and finally dissolve by Si loss to the SiO2/Si interface. Above $`2\times 10^{15}\text{ }\mathrm{Si}^+\mathrm{cm}^2`$ a ”spinodal decomposition” regime is identified, where due to the vanishing nucleation barrier LeGoues et al. (1984) non-spherical, elongated Si structures are formed (Fig. 2 (b), 100 kMCS). At even higher Si concentrations ($`1\times 10^{16}\text{ }\mathrm{cm}^2`$), above the percolation threshold, the phase separated Si becomes laterally connected (Fig. 2 (c), 100 kMCS). This network of Si does not decay into droplets during longer annealing (Fig. 2 (c), 4000 kMCS). An electrical charge brought to this network can spread over several tens of nanometer, i.e. the phase separated Si behaves like a floating gate in a conventional metal-oxide-silicon transistor.

Below the percolation threshold at $`5\times 10^{15}\text{ }\mathrm{cm}^2`$, the initially non-spherical Si structures evolve into spherical NCs, which can hardly be distinguished at this late stage from that formed by nucleation and growth (see Fig. 2 (a),(b), 4000 kMCS). The SiO2/Si interface acts as an effective sink for Si in both regimes, which results in a zone denuded of NCs. However, a more detailed consideration reveals some differences between the nucleation regime and the spinodal decomposition regime. As shown in Fig. 3 a), for the nucleation regime ($`2\times 10^{15}\text{ }\mathrm{cm}^2`$) the width of the denuded zone is constant over long annealing times. A similar dependence is observed for the NC size (Fig. 3 (b)). The competition between Ostwald ripening and Si loss to the Si/SiO2 interface keeps the mean NC diameter constant over a long period of annealing. The Si loss manifests itself in a rapid decrease in NC density (Fig. 3 (b)). On the other hand, for the spinodal decomposition regime($`5\times 10^{15}\text{ }\mathrm{cm}^2`$), the interface minimization of the non-spherical Si structures leads to a narrowing of the denuded zone (Fig. 3 (a)). Moreover, Ostwald ripening is more effective than Si loss to the interface and, hence, up to $`4000\text{ }\mathrm{kMCS}`$ the mean NC size is increasing as can be seen in Fig. 3 (b). In contrast to the nucleation regime, the NCs become larger and dissolve slower ( Fig. 3 (b)). Of course, a long lasting annealing results also in this regime in a complete dissolution of NCs.

It should be emphasized that NCs form behind a zone completely denuded of NCs (see cross section views of Fig, 2). Additionally, the distance of the NCs from the interface is small enough to allow their charging by direct electron tunneling. This self-alignment of NCs as well as their degradation free charging/decharging is crucial for application in nonvolatile memories De Salvo et al. (2001). Varying ion implantation energy and annealing temperature gives additional control over the width of the depleted zone and the NC size, which will be described elsewhere Müller and Heinig . Using Si implantation profiles predicted by TRIM instead of TRIDYN lead to a strong overestimation of the percolation behavior during phase separation, which will be shown there too.

In conclusion, two regimes of Si NC formation by 1 keV Si+ ion implantation into thin SiO2 and subsequent annealing have been found by process simulation. Below a Si fluence of $`2\times 10^{15}\text{ }\mathrm{cm}^2`$, NCs form by nucleation and growth, while at higher fluences spinodal decomposition occurs. At $`1\times 10^{16}\text{ }\mathrm{cm}^2`$, percolation leads to a spatially connected 2D pattern of Si in SiO2. The present KMC simulations predict that NC fabrication for nonvolatile memory applications should be performed preferably in the nucleation regime. In this regime, the width of the denuded zone does not depend critically on the annealing time and/or temperature. The synthesized NCs are small ($`\mathrm{2..3}\text{ }\mathrm{nm}`$ diameter) and of high area density ($`>10^{12}`$ cm-2).

This work was sponsored by the European Community under the auspices of the GROWTH project GRD1-2000-25619.