text:
# Equivalence of two approaches for the inhomogeneous density in the canonical ensemble

## Abstract

In this article we show that the inhomogeneous density obtained from a density-functional theory of classical fluids in the canonical ensemble (CE), recently presented by White et al \[Phys. Rev. Lett. 84, 1220 (2000)\], is equivalent to first order to the result of the series expansion of the CE inhomogeneous density introduced by González et al \[Phys. Rev. Lett. 79, 2466 (1997)\].

A statistical mechanics ensemble is a collection of identical systems under the same external conditions. Although the choice of a particular ensemble for studying a concrete system should be guided by the conditions in which the system is found, one can choose —due to mathematical or computational convenience— any ensemble for analyzing the equilibrium properties of the system. This way of proceeding, based on the equivalence of the ensembles in the thermodynamic limit, is only justified for systems with a very large number of particles. For small systems, however, the ensembles are no longer equivalent and the external conditions must determine the choice of ensemble.

In this context, the use of density-functional theory (DFT) for the study of classical inhomogeneous fluids has been usually limited to the grand canonical ensemble (GCE), where the temperature $`T`$ and the chemical potential $`\mu `$ are fixed by an external reservoir. A large variety of inhomogeneous situations has been successfully studied by means of DFT in the GCE . These situations include fluids confined in narrow pores or capillaries, or even spherical cavities , which are implicitly assumed to be open, i.e., allowing exchange of particles with a reservoir. This assumption is crucial for situations with a small number of particles where, depending on the choice of ensemble, important differences may arise in the equilibrium microscopic structure of the system . If one wishes to investigate the properties of a small closed system at temperature $`T`$, the study must be performed in such a way that one obtains results in the canonical ensemble (CE) because the number of particles $`N`$ is fixed. In DFT this goal can be achieved by means of two different approaches. On one hand the DFT could be formulated in the canonical ensemble , with a minimum free-energy principle with fixed $`T`$ and $`N`$, and an appropriate CE functional. Very recently, this approach has been explicitly realized  by considering an approximate expression for the CE functional. On the other hand, one can perform a conventional DFT study in the GCE and then relate the obtained properties to those of the CE. This approach was followed in Refs.  where the CE density profile of a hard-sphere fluid in a small spherical cavity was calculated by means of a series expansion in terms of the corresponding GCE profile. The aim of the present paper is to show that these two approaches yield equivalent results to order $`1/N`$. For clarity we start with a brief summary of the main results of the two approaches.

The first approach is based on the following series expansion of the CE density profile $`\rho _\text{c}(𝐫)`$ in terms of its corresponding GCE density profile $`\rho _{\text{gc}}(𝐫)`$:

$$\rho _\text{c}(𝐫)=\rho _{\text{gc}}(𝐫)\frac{1}{2}\mathrm{\Delta }^2(N)\frac{^2}{N^2}\rho _{\text{gc}}(𝐫)+O(\frac{1}{N^2})$$
(1)
where the grand canonical profile is obtained for a chemical potential $`\mu `$ that leads to an average number of particles $`N`$ equal to $`N`$ —the fixed integer number of particles in the CE, and $`\mathrm{\Delta }^2(N)N^2N^2`$ is the mean square fluctuation of the number of particles in the GCE. Higher-order terms in the above expansion also depend on fluctuations and on variations of the GCE profile w.r.t. $`N`$ . In DFT, the grand canonical profile $`\rho _{\text{gc}}`$ is the solution of the usual GCE Euler-Lagrange equation with chemical potential $`\mu `$ and external potential $`V_{\text{ext}}(𝐫)`$,

$$\frac{\delta _{\text{gc}}[\rho ]}{\delta \rho (𝐫)}|_{\rho =\rho _{\text{gc}}}+V_{\text{ext}}(𝐫)=\mu ,$$
(2)
where $`_{\text{gc}}[\rho ]`$ is the GCE free-energy functional. For a chemical potential leading to a given $`N`$, Eq. (2) can be rewritten as :

$$\rho _{\text{gc}}(𝐫)=N\mathrm{exp}\left[\beta V_{\text{ext}}(𝐫)+c^{(1)}(𝐫;[\rho _{\text{gc}}])\right]/𝑑𝐫\mathrm{exp}\left[\beta V_{\text{ext}}(𝐫)+c^{(1)}(𝐫;[\rho _{\text{gc}}])\right],$$
(3)
where $`\beta =1/k_\text{B}T`$ is the inverse temperature and $`c^{(1)}`$ is the one-body direct correlation function

$$c^{(1)}(𝐫;[\rho ])=\beta \frac{\delta (_{\text{gc}}[\rho ]_{\text{gc-id}}[\rho ])}{\delta \rho (𝐫)},$$
(4)
being $`_{\text{gc-id}}`$ the usual ideal-gas free-energy. This correlation function is the first member of the direct correlation hierarchy

$$c^{(n)}(𝐫_1,\mathrm{},𝐫_n;[\rho ])=\beta \frac{\delta ^n(_{\text{gc}}[\rho ]_{\text{gc-id}}[\rho ])}{\delta \rho (𝐫_1)\mathrm{}\delta \rho (𝐫_n)}.$$
(5)
From Eq. (3), in DFT it is possible to obtain density profiles normalized for a given $`N`$. This allows for obtaining approximate CE profiles using (1) where the derivatives w.r.t. $`N`$ are calculated numerically. (We note that, using the thermodynamic identity $`\mathrm{\Delta }^2(N)=N/(\beta \mu )`$, the mean square fluctuation can also be expressed as a derivative w.r.t. $`N`$.) This procedure was used in  to obtain the CE density profile of a hard-sphere fluid confined in a hard spherical cavity.

The second approach consists of an approximate expression for the free energy functional in the CE. On the basis of the standard saddle point relation between the CE Helmholtz free energy and the GCE grand potential  the following approximation for the CE free-energy functional $`_\text{c}`$ was proposed in Ref. :

$$\beta _\text{c}[\rho ]\beta _{\text{gc}}[\rho ]+\frac{1}{2}\mathrm{log}2\pi \mathrm{\Delta }^2(N;[\rho ])$$
(6)
where the functional dependence of the GCE mean square fluctuation $`\mathrm{\Delta }^2(N)`$ is made explicit. Since we are now working in the canonical ensemble, the equilibrium density profile $`\rho _\text{c}(𝐫)`$ is obtained by minimizing the functional $`_\text{c}[\rho ]+𝑑𝐫\rho (𝐫)V_{\text{ext}}(𝐫)`$ subject to the constraint

$$𝑑𝐫\rho _\text{c}(𝐫)=N.$$
(7)
Using the Lagrange multiplier technique one obtains

$$\frac{\delta _\text{c}[\rho ]}{\delta \rho (𝐫)}|_{\rho =\rho _\text{c}}+V_{\text{ext}}(𝐫)=\lambda $$
(8)
where the Lagrange multiplier $`\lambda `$ must be calculated from the constraint (7). This equation can be re-expressed as

$`\rho _\text{c}(𝐫)`$ $`=`$ $`N\mathrm{exp}\left[\beta V_{\text{ext}}(𝐫)+c^{(1)}(𝐫;[\rho _\text{c}])+\xi (𝐫;[\rho _\text{c}])\right]/{\displaystyle 𝑑𝐫}`$ (10)
$`\times \mathrm{exp}\left[\beta V_{\text{ext}}(𝐫)+c^{(1)}(𝐫;[\rho _\text{c}])+\xi (𝐫;[\rho _\text{c}])\right]`$
where

$$\xi (𝐫;[\rho ])\beta \frac{\delta (_\text{c}[\rho ]_{\text{gc}}[\rho ])}{\delta \rho (𝐫)},$$
(11)
which, for the saddle-point (SP) approximation (6), becomes

$$\xi (𝐫;[\rho ])\xi _{\text{sp}}(𝐫;[\rho ])=\frac{1}{2}\frac{1}{\mathrm{\Delta }^2(N;[\rho ])}\frac{\delta \mathrm{\Delta }^2(N;[\rho ])}{\delta \rho (𝐫)}.$$
(12)
In order to calculate $`\xi _{\text{sp}}(𝐫;[\rho ])`$ it is important to express the mean square fluctuation $`\mathrm{\Delta }^2(N;[\rho ])`$ as a functional of the density. This can be done conveniently in the GCE by means of the density-density correlation function

$$𝒢(𝐫_1,𝐫_2)=\beta ^1\frac{\delta \rho (𝐫_1)}{\delta (\mu V_{\text{ext}}(𝐫_2))},$$
(13)
since this function normalizes to the mean square fluctuation, i.e.,

$$\mathrm{\Delta }^2(N)=𝑑𝐫_1𝑑𝐫_2𝒢(𝐫_1,𝐫_2).$$
(14)
In addition, taking into account that $`𝒢`$ is the functional inverse of the second derivative of the GCE free-energy

$`𝒢^1(𝐫_1,𝐫_2)`$ $`=`$ $`\beta {\displaystyle \frac{\delta (\mu V_{\text{ext}}(𝐫_1))}{\delta \rho (𝐫_2)}}=\beta {\displaystyle \frac{\delta ^2_{\text{gc}}[\rho ]}{\delta \rho (𝐫_1)\rho (𝐫_2)}}`$ (15)
$`=`$ $`{\displaystyle \frac{1}{\rho (𝐫_1)}}\delta (𝐫_1𝐫_2)c^{(2)}(𝐫_1,𝐫_2),`$ (16)
and satisfies the Ornstein-Zernike relation

$$𝑑𝐫_2𝒢^1(𝐫_1,𝐫_2)𝒢(𝐫_2,𝐫_3)=\delta (𝐫_1𝐫_3),$$
(17)
one obtains

$`\mathrm{\Delta }^2(N;[\rho ])`$ $`=`$ $`{\displaystyle 𝑑𝐫\mathrm{\Gamma }(𝐫)}`$ (18)
$`{\displaystyle \frac{\delta }{\delta \rho (𝐫)}}\mathrm{\Delta }^2(N;[\rho ])`$ $`=`$ $`{\displaystyle 𝑑𝐫_1𝑑𝐫_2\frac{\delta 𝒢^1(𝐫,𝐫_1)}{\delta \rho (𝐫_2)}\mathrm{\Gamma }(𝐫_1)\mathrm{\Gamma }(𝐫_2)}`$ (19)
$`=`$ $`\left({\displaystyle \frac{\mathrm{\Gamma }(𝐫)}{\rho (𝐫)}}\right)^2+{\displaystyle 𝑑𝐫_1𝑑𝐫_2c^{(3)}(𝐫,𝐫_1,𝐫_2)\mathrm{\Gamma }(𝐫_1)\mathrm{\Gamma }(𝐫_2)}`$ (20)
where

$$\mathrm{\Gamma }(𝐫)𝑑𝐫_1𝒢(𝐫,𝐫_1),$$
(21)
is obtained from the following averaged Ornstein-Zernike relation

$$\mathrm{\Gamma }(𝐫)=\rho (𝐫)+\rho (𝐫)𝑑𝐫_1\mathrm{\Gamma }(𝐫_1)c^{(2)}(𝐫,𝐫_1).$$
(22)
In deriving (19) we have considered the functional derivative w.r.t. density of the Ornstein-Zernike relation (17) and exploited the fact that $`𝒢`$ and $`𝒢^1`$ are functional inverses. We note that the key difference between the GCE result (3) and the CE density (10) is the term $`\xi (𝐫;[\rho _\text{c}])`$. We also note that, in this approach, using (10) one directly obtains the CE profile while, in the previous approach, the result (3) of GCE-DFT had to be inserted into (1) in order to obtain an approximation for the CE equilibrium density. In what follows we shall show that both approaches, agree to first order, though they yield different results due to higher-order contributions in the saddle point approach. We first derive some useful relations and then we show the equivalence to first order of the approaches.

Our starting point is the well-known result of GCE density functional theory that, for given intermolecular potential, chemical potential $`\mu `$, and temperature $`T`$, only one external potential can determine a specified equilibrium density profile . Thus there must exist an external potential $`\overline{V}_{\text{ext}}(𝐫)`$ so that its corresponding GCE equilibrium density $`\rho _{\text{gc}}(𝐫;[\overline{V}_{\text{ext}}])`$ (here the functional dependence of $`\rho _{\text{gc}}`$ is made explicit) is equal to the CE result $`\rho _\text{c}`$. Performing a functional expansion of $`\rho _{\text{gc}}(𝐫;[\overline{V}_{\text{ext}}])`$ about $`V_{\text{ext}}`$ and using definition (13) we obtain:

$`\beta ^1\mathrm{\Delta }\rho (𝐫)`$ $`=`$ $`{\displaystyle 𝑑𝐫_1𝒢(𝐫,𝐫_1)\mathrm{\Delta }V_{\text{ext}}(𝐫_1)}`$ (24)
$`+{\displaystyle \frac{1}{2!}}{\displaystyle 𝑑𝐫_1𝑑𝐫_2\frac{\delta 𝒢(𝐫,𝐫_1)}{\delta (\mu V_{\text{ext}}(𝐫_2))}\mathrm{\Delta }V_{\text{ext}}(𝐫_1)\mathrm{\Delta }V_{\text{ext}}(𝐫_2)}+\mathrm{},`$
where

$$\mathrm{\Delta }\rho (𝐫)=\rho _{\text{gc}}(𝐫;[\overline{V}_{\text{ext}}])\rho _{\text{gc}}(𝐫;[V_{\text{ext}}])=\rho _\text{c}(𝐫)\rho _{\text{gc}}(𝐫)$$
(25)
and

$$\mathrm{\Delta }V_{\text{ext}}(𝐫)=V_{\text{ext}}(𝐫)\overline{V}_{\text{ext}}(𝐫).$$
(26)
Therefore, Eq. (24) provides a link between $`\mathrm{\Delta }\rho `$ and $`\mathrm{\Delta }V_{\text{ext}}`$ via a functional expansion where the coefficients belong to the standard distribution function hierarchy.

At this point we would like to emphasize the role played by $`\overline{V}_{\text{ext}}`$ in the present work. $`\overline{V}_{\text{ext}}`$ is the external potential that, at chemical potential $`\mu `$, yields the canonical profile $`\rho _\text{c}`$ in a GCE approach and this implies that $`\rho _\text{c}`$ is the solution of the following GCE Euler-Lagrange equation

$$\frac{\delta _{\text{gc}}[\rho ]}{\delta \rho (𝐫)}|_{\rho =\rho _\text{c}}+\overline{V}_{\text{ext}}(𝐫)=\mu .$$
(27)
This fact makes meaningful the use of functionals of $`\rho _\text{c}`$ like $`𝒢^1(𝐫_1,𝐫_2;[\rho _\text{c}])=`$ $`\beta \delta (\mu \overline{V}_{\text{ext}}(𝐫_1))/\delta \rho _\text{c}(𝐫_2)`$, its inverse, $`𝒢(𝐫_1,𝐫_2;[\rho _\text{c}])=`$ $`𝒢(𝐫_1,𝐫_2;[\overline{V}_{\text{ext}}[\rho _\text{c}]])=`$ $`\beta ^1\delta \rho _\text{c}(𝐫_1)/\delta (\mu \overline{V}_{\text{ext}}(𝐫_2))`$, or the mean square fluctuation $`\mathrm{\Delta }^2(N;[\rho _\text{c}])`$. On the other hand, since $`\rho _\text{c}(𝐫)`$ is the equilibrium density profile in the canonical ensemble, it is the solution of the CE Euler-Lagrange equation (8) where, in comparison with this GCE equation, the free energy is $`_\text{c}`$ —the CE functional, the external potential is $`V_{\text{ext}}`$, and the Lagrange multiplier $`\lambda `$ is used in the place of the chemical potential $`\mu `$. From Eqs. (8) and (27) and definition (11) we obtain

$$\mathrm{\Delta }V_{\text{ext}}(𝐫)=\lambda \mu +\beta ^1\xi (𝐫;[\rho _\text{c}]).$$
(28)
In the uniform limit, where $`\rho _{\text{gc}}(𝐫)\rho _0N/V=N/V`$ and also $`\rho _\text{c}(𝐫)\rho _0`$, from (2) and (27) one has $`\mathrm{\Delta }V_{\text{ext}}(𝐫)=0`$ and thus

$$\mu \lambda =\beta ^1\xi (\rho _0),$$
(29)
and (28) can be rewritten as

$$\beta \mathrm{\Delta }V_{\text{ext}}(𝐫)=\xi (𝐫;[\rho _\text{c}])\xi (\rho _0).$$
(30)
Using the saddle-point approximation $`\xi _{\text{sp}}`$ \[Eq. (12)\], this expression could be employed in (24) to obtain an approximation for the difference $`\mathrm{\Delta }\rho `$. Conversely, since Eq. (1) gives an approximation for $`\mathrm{\Delta }\rho `$, an expansion inverse to (24) would provide a way to obtain $`\xi `$. This inverse expression can be easily derived by substituting the functional expansion of $`\delta _{\text{gc}}[\rho ]/\delta \rho (𝐫)`$ about $`\rho _{\text{gc}}`$ in Eq. (27). Using definition (15), we obtain

$`\beta \mathrm{\Delta }V_{\text{ext}}(𝐫)`$ $`=`$ $`{\displaystyle 𝑑𝐫_1𝒢^1(𝐫,𝐫_1)\mathrm{\Delta }\rho (𝐫_1)}`$ (32)
$`+{\displaystyle \frac{1}{2!}}{\displaystyle 𝑑𝐫_1𝑑𝐫_2\frac{\delta 𝒢^1(𝐫,𝐫_1)}{\delta \rho _{\text{gc}}(𝐫_2)}\mathrm{\Delta }\rho (𝐫_1)\mathrm{\Delta }\rho (𝐫_2)}+\mathrm{}`$
where we have exploited the fact that $`\rho _{\text{gc}}(𝐫)`$ is the solution of the usual GCE Euler-Lagrange equation (2).

Expansions (24) and (32) are asymptotically exact relations linking $`\mathrm{\Delta }\rho `$ and $`\mathrm{\Delta }V_{\text{ext}}`$. However, these expansions need to be truncated in order to become suitable for practical applications. In particular, to first order, Eq. (24) becomes

$$\beta ^1\mathrm{\Delta }\rho (𝐫)𝑑𝐫_1𝒢(𝐫,𝐫_1)\mathrm{\Delta }V_{\text{ext}}(𝐫_1),$$
(33)
and Eq. (32) reduces to

$$\beta \mathrm{\Delta }V_{\text{ext}}(𝐫)𝑑𝐫_1𝒢^1(𝐫,𝐫_1)\mathrm{\Delta }\rho (𝐫_1).$$
(34)
Approximations (33) and (34) are, by virtue of Eq. (17), equivalent equations; this fact shows the consistency of the truncation of the expansions. Either (33) or (34) provide a simple (first order) relation between the differences $`\mathrm{\Delta }\rho `$ and $`\mathrm{\Delta }V_{\text{ext}}`$. Using these equations, we shall show that the approximation (1) for $`\mathrm{\Delta }\rho `$ is equivalent to first order to the saddle point approximation $`\xi _{\text{sp}}`$ \[Eq. (12)\]. By considering the derivative of the GCE Euler-Lagrange equation (2) w.r.t. $`N`$ at fixed $`V_{\text{ext}}`$, one obtains the exact relation

$$𝑑𝐫_1𝒢^1(𝐫,𝐫_1)\frac{\rho _{\text{gc}}(𝐫_1)}{N}=\frac{1}{\mathrm{\Delta }^2(N)}$$
(35)
where the chain rule for functional differentiation together with definition (15) and the identity $`N/(\beta \mu )=\mathrm{\Delta }^2(N)`$ have been used. This equation can be rewritten, via Eqs. (17) and (21), as

$$\frac{\rho _{\text{gc}}(𝐫)}{N}=\frac{\mathrm{\Gamma }(𝐫)}{\mathrm{\Delta }^2(N)}.$$
(36)
Differentiating this equation w.r.t. $`N`$ and using Eq. (1), we obtain

$$\mathrm{\Delta }\rho (𝐫)\frac{1}{2}\left(\frac{\mathrm{\Gamma }(𝐫)}{N}\frac{\mathrm{\Gamma }(𝐫)}{\mathrm{\Delta }^2(N)}\frac{\mathrm{\Delta }^2(N)}{N}\right),$$
(37)
which inserted into (34) yields the following approximation for $`\mathrm{\Delta }V_{\text{ext}}`$:

$$\beta \mathrm{\Delta }V_{\text{ext}}(𝐫)\frac{1}{2\mathrm{\Delta }^2(N)}\frac{\mathrm{\Delta }^2(N)}{N}\frac{1}{2}𝑑𝐫_1𝒢^1(𝐫,𝐫_1)\frac{\mathrm{\Gamma }(𝐫_1)}{N}$$
(38)
where we have used the identity

$$𝑑𝐫_1𝒢^1(𝐫,𝐫_1)\mathrm{\Gamma }(𝐫_1)=1$$
(39)
which follows from Eqs. (35) and (36) \[or, equivalently, from (17) and (21)\]. Considering the derivative of Eq. (39) w.r.t. $`N`$, and using the chain rule and (36), Eq. (38) can be re-expressed as

$`\beta \mathrm{\Delta }V_{\text{ext}}(𝐫)`$ $``$ $`{\displaystyle \frac{1}{2\mathrm{\Delta }^2(N)}}\left({\displaystyle \frac{\mathrm{\Delta }^2(N)}{N}}+{\displaystyle 𝑑𝐫_1𝑑𝐫_2\frac{\delta 𝒢^1(𝐫,𝐫_1)}{\delta \rho _{\text{gc}}(𝐫_2)}\mathrm{\Gamma }(𝐫_1)\mathrm{\Gamma }(𝐫_2)}\right)`$ (40)
$`=`$ $`{\displaystyle \frac{1}{2\mathrm{\Delta }^2(N)}}\left({\displaystyle \frac{\mathrm{\Delta }^2(N)}{N}}{\displaystyle \frac{\delta \mathrm{\Delta }^2(N)}{\delta \rho _{\text{gc}}(𝐫)}}\right),`$ (41)
where in the last equality we have used (19). Comparing (41) with (30) we obtain

$$\xi (𝐫;[\rho _\text{c}])\frac{1}{2\mathrm{\Delta }^2(N)}\frac{\delta \mathrm{\Delta }^2(N)}{\delta \rho _{\text{gc}}(𝐫)}=\xi _{\text{sp}}(𝐫;[\rho _{\text{gc}}])$$
(42)
and

$$\xi (\rho _0)\frac{1}{2\mathrm{\Delta }^2(N)}\frac{\mathrm{\Delta }^2(N)}{N}=\xi _{\text{sp}}(\rho _0).$$
(43)
where $`\xi _{\text{sp}}`$ was defined in (12). Finally, we note that in (42) it is shown that $`\xi [\rho _\text{c}]`$ is approximately equal to $`\xi _{\text{sp}}[\rho _{\text{gc}}]`$, i.e., $`\xi _{\text{sp}}`$ evaluated at $`\rho _{\text{gc}}`$ instead of $`\rho _\text{c}`$. This approximate equality also holds for $`\xi _{\text{sp}}[\rho _\text{c}]`$ as we show in what follows. Expanding $`\xi _{\text{sp}}`$ about $`\rho _{\text{gc}}`$ we obtain

$$\xi _{\text{sp}}(𝐫;[\rho _\text{c}])=\xi _{\text{sp}}(𝐫;[\rho _{\text{gc}}])+𝑑𝐫_1\frac{\delta \xi _{\text{sp}}(𝐫;[\rho _{\text{gc}}])}{\delta \rho _{\text{gc}}(𝐫_1)}\mathrm{\Delta }\rho +\mathrm{}$$
(44)
and, taking into account that $`\xi `$ is already a $`O(\mathrm{\Delta }\rho )`$ quantity,

$$\xi _{\text{sp}}(𝐫;[\rho _\text{c}])\xi _{\text{sp}}(𝐫;[\rho _{\text{gc}}])+O(\mathrm{\Delta }\rho )^2$$
(45)
This proves the equivalence to first order of the two approaches for obtaining an approximate density profile in the canonical ensemble.

In summary, we have shown that two different approaches for obtaining the density profile of a fluid in the canonical ensemble are equivalent to first order. The demonstration was based on considering an external potential $`\overline{V}_{\text{ext}}`$ for which the equilibrium density in the grand canonical ensemble is precisely the canonical ensemble result. Using this external potential we have been able to work in the framework of the grand canonical ensemble where approximations similar to those carried out in this work are commonly encountered.

The proof of the equivalence gives additional support to the saddle point approximation for the CE free-energy functional introduced in . This approximation allows for a CE-DFT treatment of fluids confined in a closed cavity with excellent agreement with simulation data. However, the SP free-energy functional was proposed on the basis of the well-known SP relation between the equilibrium CE free-energy and the grand potential of a homogeneous fluid, assuming that this relation would be also a good approximation for inhomogeneous fluids. This assumption is thus reinforced by the results of the present paper which are valid for any inhomogeneous situation. As a final remark, we would like to stress the fact that our demonstration has focused on the approximate CE density rather than on the free-energy functionals, and the equivalence between the two approaches must be understood in this sense.

We thank financial support by the Comisión Interministerial de Ciencia y Tecnología of Spain under Grant PB 98-0261.


format: text

text:
# Excitation of H2 and HD in Shocks and PDRs

## 1. Introduction

Shocks and photodissociation regions show strong H2 emission which is detectable in the visible, near- and mid-infrared. Among the important H2 lines observed are those of pure rotational transitions in the vibrational ground state, at wavelengths between about 2.4 and 28 $`\mu `$m. This spectral range has been difficult to access from the ground, and only recently has the Infrared Space Observatory (ISO) opened up this window to sensitive spectroscopy. Much attention has focused on the pure rotational emission, because it provides important clues about the physical state of the emitting gas. ISO and ground-based observations were able to trace the populations of levels at high energy, up to near the dissociation limit at 4.5 eV above the ground state. Excitation diagrams such as those shown in Figs. 1 or 4 indicate that in shocks and PDRs (see Wright, this volume), the higher energy levels typically show remarkably high excitation temperatures, of order 2000 to 3000 K. Could this emission from high levels actually be thermally excited ? Although in principle possible, the kinetic temperatures and gas densities required for their thermalization seem uncomfortably high. What then could be the cause of this high excitation ? In the following, we will briefly discuss several possible mechanisms.

## 2. H2 in Shocks: the Orion Outflow

OMC-1, the molecular cloud behind the Orion M42 Nebula, is the best-studied massive star-forming region. It embeds a young stellar object which produces a spectacular outflow that shocks the surrounding molecular gas, producing strong H2 infrared line emission fragmented in a multitude of clumps and arcs (Schultz et al. 1998). Peak 1 (Beckwith et al. 1978) is the brighter of the two H2 emission concentrations of the outflow. Toward Peak 1 we obtained 2.4–45 $`\mu `$m spectra with ISO, which show over 60 $`\mathrm{H}_2`$ ro-vibrational and pure rotational lines, the latter ranging from 0–0 S(1) to 0–0 S(25) (Rosenthal et al. 1999, 2000). The extinction-corrected level column density distribution derived from the observed line intensities is shown in Fig. 1 (left).

### 2.1. Excitation of high-energy H2 levels toward Peak 1

Only a fraction $``$ 10-3 of the total column density 1.9 $`\times 10^{21}`$ cm-2 of warm H2 toward Peak 1 resides in high energy levels, $`E(v,J)/k`$$`>`$ 10,000 K. For these levels the distribution of column densities implies an excitation temperature, $`T_{ex}`$, of about 3000 K; $`T_{ex}`$ is defined by a fit to the column densities of the form $`N(v,J)/g`$$``$$`\mathrm{exp}(E(v,J)/kT_{ex})`$. If the observed $`N(v,J)`$ is decomposed as a sum of LTE distributions with a range of excitation temperatures, then we find that about 99% of all warm H2 has $`T_{ex}`$$`<`$ 2000 K, and 90% has $`T_{ex}`$$`<`$ 1000 K (see Rosenthal et al. 2000). If the distribution of excitation temperatures were to reflect the actual distribution of kinetic temperatures, such a steep distribution of temperatures would be difficult to reconcile with the expected smooth temperature profile of planar C-type shocks. In C-type shocks the gas temperature changes smoothly, and a large fraction of the warm gas is near the maximum temperature reached in the shock (Timmermann 1996). Even with a distribution of shocks of different velocities and a corresponding range of peak temperatures, an excitation temperature distribution like that found toward Peak 1 is difficult to construct, unless one were to combine shocks with very disparate filling factors, e.g., several percent at 40 km s-1, the rest at 20 km s-1. Even in bow shocks, which are often cited to match the observed excitation (Smith 1991), the velocity changes slowly with distance from the apex, and the observed distribution of temperatures would hardly be expected.

In dissociative J-type shocks (e.g. Hollenbach & McKee 1989), the molecules are destroyed in the shock, and they reform in a postshock layer where the temperature has dropped much below 3000 K, somewhat dependent on the H2 formation rate at higher temperatures, which is very uncertain (e.g. Bertoldi 1997). Dissociative J-shocks can therefore also not account for the high excitation H2 we observe.

Even if temperatures of 3000 K or more can be reached in nondissociative shocks, the higher H2 levels would remain sub-thermally excited unless the gas density is high enough that the collisional excitation and deexcitation rates are comparable to those for radiative decay. A “critical” gas density can be estimated for a given level as that density for which the total collisional deexcitation rate of the level equals its total radiative decay rate. We find that even at kinetic temperatures of 3000 K, gas densities well above $`10^6`$ cm-3 would be needed to maintain the higher $`v`$=0 levels at populations resembling LTE. Since such high densities may not prevail in the shocked gas of the Orion outflow, we explore mechanisms other than thermal excitation that might account for high populations in the upper energy states.

### 2.2. Time-dependent C-shocks

When a high velocity outflow strikes dense molecular gas, and a C-type shock is first established, J-type shocks can temporarily form within the C-shock (Chièze et al. 1998; Flower & Pineau des Forêts 1999). In such an embedded J-shock, a small column of molecular gas is briefly heated to high temperatures without being dissociated, and if the gas density is sufficiently high, this would produce a high-excitation tail in the column density distribution. The lifetime of the embedded J-shock is small, and the high-excitation populations are transient, unless shocks were continuously reforming. Embedded J-shocks may also form when a C-shock encounters a dense clump. Since the time scale for the fading of the high-excitation tail in newly forming C-shocks is somewhat smaller than the age of the Orion outflow, we should expect spatial variations in the upper level excitation. Whether such gradients exist is unknown, but it would provide an interesting test for the transient J-shock model.

### 2.3. H2 formation pumping

H2 formation is one other mechanism to populate the high-energy H2 levels. Molecular hydrogen is believed to form on the surfaces of dust grains. Some fraction of the 4.5 eV released during the formation of an H2 molecule is required to leave the grain, and the remainder is split between translation, rotation, and vibration of the new molecule.

Toward Peak 1 we found that a fraction 6.8 $`\times 10^4`$ of the total warm H2 column is in states with energy $`E(v,J)/k`$$``$ 10,000 K. Could H2 formation in a steady state account for such a fraction of molecules in highly excited states ? We can try a very simple estimate. The pumping rate due to formation pumping is equal to the H2 formation rate, $`n(\mathrm{H})n_\mathrm{H}R_{gr}`$, where we adopt $`R_{gr}`$$``$ 5 $`\times 10^{17}`$ cm3 s-1 as the H2 formation rate coefficient per hydrogen nucleus. The radiative decay rate we estimate from a characteristic radiative lifetime of a level at $``$ 2 eV, which is $``$ 106 sec, and we prolong this by an average number of $``$ 5 transitions for the molecule to reach its ground state. Thereby defining an effective $`A`$-coefficient for the decay rate of highly excited states, $`A_x`$$``$ 2 $`\times 10^7`$ s-1, the population balance for the excited states writes

$$R_{gr}n_\mathrm{H}n(\mathrm{H})=n_x(\mathrm{H}_2)A_x,$$
(1)
which yields an excited H2 fraction

$$\frac{n_x(\mathrm{H}_2)}{n(\mathrm{H}_2)}=\frac{n(\mathrm{H})}{n(\mathrm{H}_2)}\frac{n_\mathrm{H}R_{gr}}{A_x}=5\times 10^4\left(\frac{n_\mathrm{H}}{10^6\mathrm{cm}^3}\right)\left(\frac{n(\mathrm{H})}{2n(\mathrm{H}_2)}\right),$$
(2)
that agrees well with the observed value. This simple estimate shows that H2 formation could in principle account for some of the high-excitation level populations, provided the density is high enough, the atomic fraction is not small, and the formation rate coefficient in the warm shocked gas is comparable to the value implied at $`100`$ K from Copernicus observations, $`R_{gr}`$$``$ 3 $`\times 10^{17}`$ cm3 s-1 (Jura 1975).

Can H2 formation also produce the observed shape of the high energy level populations ?

The energy required to desorb from the dust grain depends on whether the atoms are bound by van der Waals forces (physisorbed, $``$ 0.03 eV), such as on water ice, or by chemical bonds (chemisorbed, $``$ 2 eV), like on a silicate surface. The exact level distribution of newly formed H2 is unknown, but it could very well contribute to the observed excitation at intermediate energies, $`E(v,J)`$$``$ 1–3 eV. Experiments and theoretical calculations are yet inconclusive regarding the excitation of a newly formed, desorbed hydrogen molecule. Hunter & Watson (1978) consider physisorption and find a partition between rotational, vibrational, and translational energy, $`E_r`$$``$ 0.8 eV, $`E_v`$$``$ 3.3 eV, $`E_t`$$``$ 0.8 eV. Duley & Williams (1986) assume the hydrogen to be physisorbed on silicate surfaces and find $`E_r`$$`<`$ 0.01 eV, $`E_v`$$``$ 3.3 eV, $`E_t`$$``$ 0.2 eV. Duley & Williams (1993) find that on PAH and graphite surfaces, $`E_r`$$`<`$ 0.01 eV, $`E_v`$$`<`$ 10-3 eV, whereas Duley & Williams (1993) and Leonas & Pjarnpuu (1981) find that on H2O ice and polymeric hydrocarbon surfaces, $`E_v`$$``$ 4 eV. Evidently, there is no conclusive agreement yet, except for a trend to expect that newly formed molecules must be vibrationally hot and rotationally cold.

Using the PDR codes of Draine & Bertoldi (1996), we computed the non-LTE level populations of partially dissociated gas of fixed density and temperature, and assuming that the H2 formation distribution follows

$$\varphi (v,J)=\varphi (0,0)g_n(J)g_fe^{E(v,J)/kT_f},$$
(3)
where the nuclear degeneracy factor, $`g_n(J)`$=1 for even $`J`$, 3 for odd $`J`$. The formation temperature, $`T_f`$, is a free parameter, and so is the formation degeneracy, $`g_f`$. The exponential form of the level distribution (Black & van Dishoeck 1987) is somewhat arbitrary, but in lack of observational constraints, seems reasonable. Fig. 2 shows the computed level populations. The high excitation tail in the population of levels with energy $`>`$ 5000 K is due to formation pumping, whereas the lower levels are populated through collisional excitation. The difference in the level population when adopting $`g_f`$=2$`J`$+1 or $`g_f`$=1+$`v`$ is illustrated in Fig. 2 (left). When emphasizing vibrational excitation, the individual vibrational level populations show a splitting similar to that seen for fluorescent excitation. When the weight lies on rotational excitation, the level populations remain tighter, more like what we observe toward Peak 1. Fig. 2 (right) shows how the level populations depend on the value of the formation temperature, $`T_f`$, where we adopt $`g_f`$=2$`J`$+1.

No single temperature gas slab is able to reproduce the level distributions observed toward Peak 1. However, the sum of two components with hydrogen nuclei density $`n_\mathrm{H}`$ = 106 cm-3, atomic fraction $`n(\mathrm{H})/n_\mathrm{H}`$ = 0.5, and respective temperatures of 200 K and 800 K, with column densities $`N_{\mathrm{H}_2}`$ = 1.2 $`\times 10^{22}`$ cm-2 and 1.2 $`\times 10^{21}`$ cm-2, can reproduce the observed level column distribution very well (Fig. 1, right). In order to produce a narrow distribution of level columns in the Boltzmann plot, we would have to assume pronounced rotational excitation of the newly-formed H2. This is somewhat contrary to the expectation of vibrationally hot, but rotationally cold new molecules.

### 2.4. Non-thermal collisions

An even more important pumping mechanism for the high-excitation levels could be non-thermal collisions between molecules and ions in a magnetic shock. In magnetic C-type shocks, which are believed to be responsible for most of the emission in Peak 1, the gas is accelerated through fast inelastic collisions. In a magnetic precursor the ions, which are tied to the magnetic field, collide with the undisturbed preshock gas at relative velocities comparable to the shock speed. Such non-thermal ion-molecule collisions lead to the acceleration of the molecules and to their internal excitation. High-velocity molecules subsequently collide with other molecules, resulting in a cascade of collisions during which the relative kinetic energy is in part converted to internal excitation of the molecules (O’Brien & Drury 1996). In sufficiently fast C-shocks, the ion–H2 and H2–H2 collisions can lead to a significant collisional dissociation rate. The molecules dissociated in a partially dissociative, steady-state shock reform further downstream, so that across such a shock the H2 dissociation rate equals the H2 reformation rate. For every collisionally dissociated molecule there is likely to be a much larger number of inelastic collisions which did not lead to dissociation, but to the excitation of the molecule into high rovibrational states. The high-excitation H2 level column densities thereby created should then be larger than those caused by H2 formation alone.

We conclude that non-thermal collisions in partially dissociative shocks could pump the high-excitation states in the H2 electronic ground state to the levels observed. No detailed shock models are available yet which account for this process.

### 2.5. 0–0 S(25)

The $`J`$=27 level observed through the 0–0 S(25) line (Fig. 3) appears overpopulated by a factor of seven (see Fig. 1) over what would be expected from an extrapolation of the lower energy level populations. The $`J`$=27 level lies 3.6 eV above ground, only 0.9 eV from the dissociation limit. H2 that is newly formed on grains is less likely to populate rotational states so high, as some fraction of the formation energy is lost to overcome the grain surface potential, and some goes into vibrational excitation and translational kinetic energy.

Unless we misidentified the 0–0 S(25) line, it appears that a different mechanism may be populating this level, and possibly other high levels. The gas phase formation of H2 via H-, e.g., might be able to leave the new molecule in such a high rotational state.

## 3. Excitation of H2 in PDRs

A “photodissociation region” is the interface separating a region which is predominantly molecular, and a region where the ultraviolet energy density is sufficiently high that hydrogen is mostly atomic. Young O and B stars located near molecular clouds produce conspicuous PDRs. A number of general reviews of PDRs have appeared recently (Hollenbach & Tielens 1997, 1999; Sternberg, Yan, & Dalgarno 1998; Walmsley 1998), and the current status of PDR modeling with regard to ISO observations was discussed by Draine & Bertoldi (1999).

The structure of a PDR is determined primarily by the attenuation of the far-ultraviolet (6–13.6 eV) radiation field, as one moves from the ionization front into the PDR. The dominant process is the photodissociation of H2, the rate for which is determined by both H2 self-shielding (Draine & Bertoldi 1996) and attenuation by dust.

The $`(v,J)`$ ro-vibrational excited states of H2 are populated by inelastic collisions, by UV pumping, and by formation on grains. At the typical densities $`n_\mathrm{H}`$$`\stackrel{>}{}`$ 104 cm-3 of bright PDRs, collisions maintain the ($`v`$=0, $`J`$) levels of H2 in approximate thermal equilibrium for $`J`$$`\stackrel{<}{}`$ 9. Therefore measurements of the quadrupole emission line intensities with ISO provide a good indicator of the gas temperature. ISO observations toward the PDRs in, e.g., S140 (Timmermann et al. 1996), NGC 7023 (Fuente et al. 1999), or NGC 2023 (Draine & Bertoldi 2000), provided unequivocal evidence for gas temperatures in the 500–1000 K range in a portion of the photodissociation region (PDR) where the H2 fraction is appreciable. These temperatures were higher than expected from current models of the heating and cooling processes in PDRs, and therefore require reconsideration of the physics of the gas and dust in PDRs.

The higher energy H2 levels are not populated by collisions, but by the decay of electronically excited states, which were pumped through the absorption of FUV photons. Because most molecules are excited from their ro-vibrational ground states, most subsequent transitions change the rotational quantum number, $`J`$, only by one or two, and the decay back into the electronic ground state populates states within a limited range of $`J`$. The cascade down the vibrational states to the ground states can then significantly populate states up $`J`$=10 or beyond. Because $`J`$=0, 1 pose a reflective boundary for the cascade down, the lowest $`J`$ levels within each vibrational state build up enhanced populations, which results in the characteristic population distribution of fluorescently excited H2 displayed in Fig. 4.

As an example we compare the level populations observed toward a bright filament in the NGC 7023 PDR with our best fit PDR model, which adopts $`n_\mathrm{H}`$ = 105 cm-3 and an incident FUV flux $`\chi `$ = 5000 times 1.2 $`\times 10^7`$ cm-2 s-1. The observations shown include our ISO as well as ground-based spectroscopy. Our PDR model includes all basic collisional and radiative processes, including detailed thermal balance (Draine & Bertoldi 1996, 1999, 2000). The agreement between the model and the observations is remarkably good considering the simplicity of a constant density, planar slab model.

Black & van Dishoeck (1987), Le Bourlot et al. (1995), and Bertoldi (1997) questioned whether processes other than collisions and fluorescence could affect the high energy or high rotational state populations in PDRs. Burton et al. (1992) suggested that H2 formation might be responsible for an apparent excess in the $`v`$=4 levels observed in the NGC 2023 PDR.

Would we expect formation pumping to have a significant effect on the high-excitation level populations in PDRs ? Of all photo-excitations of H2 by FUV photons, about 10 to 15% lead to dissociation, the remaining to pumping of high-excitation levels in the electronic ground state. In equilibrium, dissociation is balanced by reformation of H2, so that the ratio of the rates of formation pumping and fluorescent pumping of the high-excitation levels in the electronic ground state is $``$ 15/85. Overall, fluorescent pumping should therefore dominate over formation pumping by a factor five. Unless the level distribution, $`\varphi (v,J)`$, of newly formed H2 is strongly concentrated toward a small number of states (see Black & van Dishoeck  for models with pumping of specific levels), there will be little chance to detect signs of it in the level distributions of H2 in PDRs.

However, we should emphasize that the value of the H2 formation rate coefficient, $`R_{gr}`$, does strongly affect the total emission from a PDR, since it governs its molecular abundance structure, and thereby the fraction of the incident FUV energy that is reemitted by H2, as opposed to dust.

## 4. HD Emission from the Orion Outflow

The ideal molecule for measuring the interstellar deuterium abundance is HD, whose low rotational and ro-vibrational transition lines lie in the near- and mid-infrared. We used ISO and the United Kingdom Infrared Telescope (UKIRT) to search for HD lines toward Peak 1. With ISO we there detected one faint HD line, 0–0 R(5) at 19.4 $`\mu `$m (Fig. 5, Bertoldi et al. 1999), and with UKIRT we subsequently found the 1–0 R(5) line at 2.46 $`\mu `$m, providing the first observation of HD emission with a ground-based telescope, and the first measurement of the rotation-vibrational excitation of HD in the ISM (Ramsay Howat et al., in prep.).

Surprisingly, we found that the HD excitation traced by the two observed transitions exactly mirrors the excitation of H2, despite the fact that the much faster radiative transitions of HD place the critical densities of its vibrationally excited states significantly above those of H2, and much higher than the estimated gas density in the outflow region. The cause for the similarity of the H2 and HD excitation is unresolved. It may imply a strong role of non-thermal excitation mechanisms, such as the direct coupling of the HD excitation to that of H2 through the exchange reaction $`\mathrm{H}_2`$ \+ D $``$ HD + H. Current shock models are unable to explain the observed excitation of HD.

The similarity of the HD and the well-traced H2 excitation permits the yet most accurate determination of the \[HD\]/\[H2\] abundance ratio in the warm, shocked gas, and with a 30% correction for chemical depletion of HD relative to H2 (due to an asymmetry in the above-mentioned exchange reaction – see Bertoldi et al. , Ramsay Howat et al. in prep.) we are able to estimate a deuterium abundance \[D\]/\[H\] = (5.1 $`\pm `$ 1.5) $`\times 10^6`$.

This value is low compared with previous D and H absorption measurements towards nearby stars, which found values higher by a factor of two to three. But it is still consistent with two other measurement of the deuterium abundance towards Orion (Wright et al. 1999; Jenkins et al. 1999). This indicates that in the Orion region, either stars have burnt up a large fraction of the primordial deuterium, or more likely, that there is significant trapping of deuterium on dust grains. Future observations from space and from the SOFIA airborne observatory will eventually tell.

#### Acknowledgments.

We are thankful to G. Pineau des Fôrets, J. Black, E. van Dishoeck, and C. Wright for valuable input, and to the SWS Data Center at MPE, especially to E. Wieprecht, for their support. BTD was supported in part by NSF grant AST-9619429.

## References

Beckwith, S., Persson, S., Neugebauer, G., & Becklin, E. 1978, ApJ, 223, 464

Bertoldi, F. 1997, in Proceedings of the first ISO workshop on Analytical Spectroscopy, eds. A.M. Heras et al. (Noordwijk: ESA Publications Division, ESA SP-419), 67

Bertoldi, F., Draine, B.T., Drapatz, S., Rosenthal, D., Timmermann, R., & Wright, C.M. 1999, in The Physics and Chemistry of the Interstellar Medium, eds. V. Ossenkopf et al. (GCA-Verlag Herdecke)

Bertoldi, F., Timmermann, R., Rosenthal, D., Drapatz, S., & Wright, C.M. 1999, A&A, 346, 267

Black, J. & van Dishoeck, E.F. 1987, ApJ, 322, 412

Burton, M., Bulmer, M., Moorhouse, A., Geballe, T.R., & Brand, P.W.J.L. 1992, MNRAS, 257, 1

Chièze, J.-P., Pineau des Forêts, G., & Flower, D. 1998, MNRAS, 295, 672

Draine, B.T. & Bertoldi, F. 1996, ApJ, 468, 269

———— 1999, in The Universe as seen by ISO, eds. P. Cox & M.F. Kessler (ESA SP-427)

———— 2000, in H2 in Space, eds. F. Combes & G. Pineau des Forêts, Cambridge University Press

Duley, W.W. & Williams, D.A. 1986, MNRAS, 223, 177

———— 1993, MNRAS, 260, 37

Flower, D. & Pineau des Forêts, G. 1999, MNRAS, 308, 271

Fuente, A., Martín-Pintado, J., Rodriguez-Fernandez, N.J., de Vicente, P., & Kunze, D. 1999, ApJ, 518, L45

Hollenbach, D.J. & McKee, C.F. 1989, ApJ, 342, 306

Hollenbach, D.J. & Tielens, A.G.G.M. 1997, ARA&A, 35, 179

———— 1999, Rev. Modern Phys., 71, 173

Hunter, D.A. & Watson, W.D. 1978, ApJ, 226, 477

Jenkins, E., Tripp, T., Wozniack, P., et al. 1999, ApJ, 520, 182

Jura, M. 1975, ApJ, 197, 575

Leonas, V.B. & Pjarnpuu, A.A. 1981, Sov. Ast. Let., 7, 19

Le Bourlot, J., Pineau des Forêts, G., Roueff, E., Dalgarno, A., & Gredel, R. 1995, ApJ, 449, 178

Martini, P., Sellgren, K., & Hora, J.L. 1997, ApJ, 484, 296

O’Brien, I. & Drury, L. 1996, MNRAS, 280, 550

Rosenthal, D., Bertoldi, F., Drapatz, S., & Timmermann, R. 1999, in The Universe as seen by ISO, eds. P. Cox & M.F. Kessler (ESA SP-427), 561

Rosenthal, D., Bertoldi, F., & Drapatz, S. 2000, A&A, submitted

Smith, M. 1991, MNRAS, 253, 175

Schultz, A.S.B, Colgan, S.W.J., Erickson, E.F., et al. 1999, ApJ, 511, 282

Sternberg, A., Yan, M., & Dalgarno, A. 1998, in Molecules in Astrophysics: Probes and Processes, ed. E.F. van Dishoeck (Kluwer: Dordrecht), 141

Timmermann, R., Bertoldi, F., Wright, C.M., Drapatz, S., Draine, B.T., Haser, L., & Sternberg, A. 1996, A&A, 315, L281

Timmermann, R. 1996, ApJ, 456, 631

Walmsley, C.M. 1998, Astr. Lett. & Comm., 37, 1

Wright, C.M., van Dishoeck, E.F., Cox, P., Sidher, S., & Kessler, M.F. 1999, ApJ, 515, L29


text:
# 1 Introduction

## 1 Introduction

The Ising model is one of the most studied and basic models for phase transitions. In this article, the current status of experimental studies characterizing two classic models of second-order phase transitions in short-range interaction systems in the presence of quenched disorder, the random-exchange Ising model (REIM) and the random-field Ising model (RFIM), is presented. The discussion concentrates on experiments in dilute, insulating, anisotropic antiferromagnets, the systems that have yielded the best understood data for these two models. The REIM is realized in zero magnetic field and the RFIM with a field applied along the spin-ordering direction. The REIM is rather well characterized experimentally, theoretically and through computer simulations. The $`d=2`$ RFIM is fairly well characterized, although the scaling behavior of scattering near the destroyed phase transition is still being investigated. The understanding of the RFIM for $`d=3`$ is not as complete, though significant progress has been made in the past few years, and it is this model that will be the main focus of this short review. The early history of the $`d=3`$ RFIM was fraught with controversial interpretations of the data, a result of severe nonequilibrium effects. Nevertheless, some experimental groups realized from the start that underlying the observed, complicated behavior is a new kind of phase transition. Efforts to characterize the new critical behavior were thwarted by the severe nonequilibrium effects. These nonequilibrium effects have recently been overcome by going to sufficiently high magnetic concentration and a complete characterization of the universal $`d=3`$ RFIM critical behavior is possible and underway. The most recent static critical behavior will be compared to results from computer simulations and theory. In addition to these low-field behaviors, much has been learned about the high-field limit of the RFIM. An overview will be given of the phase diagram and the different behaviors observed.

Experiments have been performed on REIM and RFIM systems for more than two decades. Since the experiments performed some time ago have been reviewed previously, they will be included here only as needed to give a perspective on the current physical understanding of the models. Theory and simulation results will be included as needed for the interpretation of the experiments. Another classic model of ordering in the presence of disorder, the spin glass, will be covered in another review in this Ising Colloquium and so will not be discussed here in detail, although some spin-glass-like behaviors do occur at low magnetic concentrations and at high magnetic fields.

Table 1 shows the most frequently measured static critical behaviors associated with a phase transition. We will make reference to the universal parameters defined in Table 1 as needed.

## 2 Experiments on Pure $`d=2`$ and $`d=3`$ Anisotropic Antiferromagnets

Observations of asymptotic static critical behavior in the pure $`d=2`$ and $`d=3`$ Ising antiferromagnets are very well documented. The magnetic specific heat ($`C_m`$) critical behavior has been characterized using the optical linear birefringence techniques on the $`d=2`$ $`\mathrm{Rb}_2\mathrm{CoF}_4`$ system and the $`d=3`$ $`\mathrm{FeF}_2`$ system. The measured pure Ising critical exponents $`\alpha `$ and the amplitude ratios $`A^+/A^{}`$ are in superb agreement with very many theoretical and simulation results. The birefringence technique is particularly useful and more accurate than pulsed specific heat techniques since it is insensitive to the large phonon contributions which are particularly difficult to handle for $`d=2`$.

The critical behavior of the staggered susceptibility and correlation length have been determined with neutron scattering for $`\mathrm{K}_2\mathrm{CoF}_4`$ for $`d=2`$ and in $`\mathrm{FeF}_2`$ for $`d=3`$. In general, the scattering line shapes for the pure and REIM systems away from the Bragg scattering point follow the scaling behavior of the spin-spin correlation function

$$\chi _s(q)=A^\pm \kappa ^{\eta 2}f(q/\kappa ),$$
(1)
where $`\gamma =\nu (2\eta )`$. For both $`d=2`$ and $`d=3`$, the scaling functions used in data analysis are approximate ones that differ significantly from the mean-field (MF) Lorentzian

$$f(q/\kappa )=\frac{1}{1+(q/\kappa )^2},$$
(2)
as can be seen in Fig. 1 where various scaling functions are compared. The deviation from the Lorentzian is more pronounced for $`d=2`$ and for $`T<T_c(H)`$ in both dimensions.

The order parameter critical behavior has been determined using neutron scattering in the $`d=2`$ compounds. The Mössbauer technique was used in the study of the $`d=3`$ system.

The results for the pure Ising model in $`d=2`$ and $`d=3`$ are summarized in Table 2. Note that the Rushbrooke scaling relation

$$2\beta +\gamma +\alpha 2$$
(3)
is satisfied as an equality for both cases. Included in Table 2 are the results from a few theoretical and simulation studies. No attempt is made to review the vast literature on the pure Ising models.

## 3 Random-Exchange Experiments in Dilute $`d=2`$ and $`d=3`$ Anisotropic Antiferromagnets

The REIM is realized in dilute, anisotropic insulating antiferromagnets when the site dilution does not result in strongly frustrated bonds (that would lead to spin-glass behavior). Random-exchange phase transitions are observed in $`d=2`$ and $`d=3`$ systems and these appear to be in good accord with theory and simulations.

The $`d=2`$ REIM $`C_m`$ critical behavior was observed using the birefringence technique on the magnetically dilute antiferromagnet $`\mathrm{Rb}_2\mathrm{Co}_{0.85}\mathrm{Mg}_{0.15}\mathrm{F}_4`$. The approximately logarithmic divergence is compatible with theoretical predictions. The scattering critical behavior of the compound $`\mathrm{Rb}_2\mathrm{Co}_\mathrm{x}\mathrm{Mg}_{1\mathrm{x}}\mathrm{F}_4`$ was analyzed using approximate scattering line shapes developed for the pure $`d=2`$ Ising model. The successful analysis using these line shapes suggests that the correct line shape is close to the pure one. The static critical behavior of the $`d=2`$ Ising model is quite well characterized by experiments and theory as shown in Table 3.

The $`d=3`$ REIM is similarly well characterized with birefringence, neutron scattering and Mössbauer experiments employing $`\mathrm{Fe}_\mathrm{x}\mathrm{Zn}_{1\mathrm{x}}\mathrm{F}_2`$ with the results shown in Table 3 along with some theoretical and simulation results with which they agree very well. The critical behavior of the specific heat of the $`d=3`$ REIM, measured with birefringence techniques, is shown in Fig. 2. Monte Carlo simulations based on the $`\mathrm{Fe}_\mathrm{x}\mathrm{Zn}_{1\mathrm{x}}\mathrm{F}_2`$ system are shown in Fig. 3. The birefringence technique yields a negative specific heat exponent $`\alpha `$ as predicted, consistent with a universality class different from the pure Ising model where $`\alpha `$ is positive. Note that, just as in the pure case, the birefringence technique is consistent with pulsed specific heat techniques, though the latter technique suffers from greater concentration gradient sensitivity and the large phonon specific heat component. The critical behaviors of the staggered susceptibility and correlation length were determined from neutron scattering experiments. The order parameter critical behavior was determined from Mössbauer studies. The scattering line shape scaling functions are not known from theory and were therefore determined directly from the scattering data in neutron scattering experiments using $`\mathrm{Fe}_{0.93}\mathrm{Zn}_{0.07}\mathrm{F}_2`$. The results shown in Fig. 1 clearly indicate that the scaling functions are fairly close to those of the pure $`d=3`$ case.

The REIM universal static critical parameters for $`d=2`$ and $`d=3`$ are shown in Table 3 along with theoretical and simulation results. In both dimensions the agreement is excellent. Note that the Rushbrooke scaling relation (Eq. 3) is satisfied as an equality for the REIM.

## 4 The $`d=3`$ Magnetic Percolation Threshold Concentration

As the magnetic percolation threshold concentration, $`x_p`$, is approached from above in zero field, the equilibrium phase transition is expected to approach zero temperature. For a particular magnetic structure, $`x_p`$ depends on what interactions exist between the different neighboring spins. For example, for $`\mathrm{Fe}_\mathrm{x}\mathrm{Zn}_{1\mathrm{x}}\mathrm{F}_2`$, $`x_p=0.245`$ provided only the dominant interaction between the body-center and corner spins is considered. Except very close to $`x_p`$, the much smaller interactions can be ignored. Close to $`x_p`$, however, the smaller interactions may drastically affect the behavior and even prevent ordering above $`x_p`$ when they frustrate the dominate interaction. Both the extremely slow dynamics near percolation and the sensitivity to tiny frustrating interactions can cause the system to exhibit spin-glass-like behavior.

A good deal of effort has focused on the properties of the system $`\mathrm{Fe}_\mathrm{x}\mathrm{Zn}_{1\mathrm{x}}\mathrm{F}_2`$ for $`x`$ near $`x_p`$. This system has a small frustrating interaction. The spin-glass-like properties were first elucidated in experiments by Montenegro, et al.. For $`H>0`$ and $`x=x_p`$, there exists a boundary that resembles a de Almeida-Thouless boundary with curvature $`TT_oH^{2/\varphi }`$ where $`\varphi =3.4`$, a typical spin-glass value. It was shown with neutron scattering that there is no antiferromagnetic long-range ordering below this boundary. Much of the behavior is very reminiscent of a canonical spin glass. The detailed behavior of this sample has been studied experimentally and extensively modeled in local mean-field and Monte Carlo simulations.

## 5 $`d=2`$ Random-Field Behavior

Scaling arguments for Zeeman and domain wall energies by Imry and Ma as well as considerations by Binder leave little doubt that the $`d=2`$ Ising transition is destroyed by the introduction of arbitrarily small random fields. Birefringence and neutron Bragg scattering experiments bear this out; no sharp phase transition is observed in equilibrium, though the rounded transition exhibits the expected scaling behavior. The equilibrium region is separated from a lower temperature region of strong hysteresis observed in the difference between data obtained upon heating after cooling in zero field to low temperatures and then applying the field (ZFC) and upon simply cooling the sample in the field (FC). The boundary separating these regions is time-scale dependent. The domain dynamics induced with the application of a magnetic field as well as those remaining after the field is removed at low temperatures have been studied experimentally and theoretically.

## 6 $`d=3`$ RFIM behavior for $`x_p<x<x_e`$ at low fields.

The behavior for concentrations between $`x_p`$ and the percolation threshold concentration for vacancies, $`x_e=1x_p`$, with a relatively small applied field occupied the bulk of early experimental efforts. Much of the controversy over interpretations of experimental data involved this region of concentration and fields. As a result of the equivalence of the dilute anisotropic antiferromagnet in small fields and the random-field ferromagnet often studied theoretically, it was believed that concentrations near $`x=0.5`$ would yield strong random-field effects in reasonably small fields and would be the best realizations of the $`d=3`$ RFIM for phase transition studies. At the time of the first experiments, it was generally believed, based on many theoretical arguments, that no phase transition would be observed. Indeed, early neutron scattering FC experiments, by Yoshizawa, et al., seemed to bear this out. In particular, a resolution-limited Gaussian Bragg peak does not occur upon FC, although subsequent experiments show that the samples retain long-range order below the phase boundary if ZFC. In contrast, the first $`C_m`$ studies yielded compelling evidence for a fundamentally new phase transition governing the behavior. The phase boundary $`TT_c(H)H^{2/\varphi }`$ behaves as predicted with $`\varphi =1.42\pm 0.03`$ for random-exchange to random-field crossover.

Despite the sharp $`C_m`$ peak, many have argued against it as evidence of the existence of a phase transition. The most recent of these is the “trompe l’oeil transition” phenomenological model. Among the assumptions of this model are that the birefringence and $`C_m`$ experiments do not yield the same behavior, the uniform magnetization is reflected by the square of the staggered magnetization, and conventional scaling is inoperative. This phenomenological model was shown to be inconsistent when all available data are considered. In contrast, more proven and conventional techniques of analyzing the experimental data, as described in this review, have been very fruitful in providing consistent results in a meaningful scaling context.

It is clear that the phase transition underlying the behavior is obscured by nonequilibrium behavior below the boundary, $`T_{eq}(H)`$, lying just above the phase transition and scaling in the same manner, albeit with a slightly larger amplitude. The equilibrium behavior above $`T_{eq}(H)`$ can be used to extrapolate the scattering data to the obscured phase transition boundary. When done carefully, the boundary determined in this way coincides with that determined via $`C_m`$ experiments, which are much less sensitive to the non-equilibrium behavior that distorts the neutron scattering data. For a long time, the nonequilibrium experiments represented the best random-field results available. One of the particularly interesting predictions of the change in the critical behavior induced by the random fields is that the order parameter critical exponent $`\beta `$ should decrease from $`0.35`$ to a value near zero. This can only be measured below $`T_c(H)`$, i.e. in the nonequilibrium region, so it was not clear what would be observed. Experiments on thin films were made for $`x=0.52`$, well below $`x_e0.76`$. Not surprisingly, the results were peculiar. The curvature of the Bragg intensity versus $`T`$ was such that it would require $`\beta >>0.5`$, which is hard to justify theoretically. Magnetic x-ray scattering data showed similar behavior near surfaces in bulk samples, though they were interpreted under the “trompe l’oeil” phenomenology and the Bragg scattering was not separated from the fluctuation scattering. For some years it appeared that the problems of metastability below $`T_{eq}(H)`$ could not be avoided, i.e. that they were intrinsic to the random-field behavior as realized in dilute antiferromagnets. However, insight into the origins of the metastable domains finally led to experiments at high magnetic concentration as a way to avoid the nonequilibrium behavior, as discussed later.

The metastable domains formed upon FC are themselves quite interesting and their dynamics were studied in some detail both experimentally and through simulations.

## 7 RFIM behavior for $`x_p<x<x_e`$ at high fields.

The general phase diagram features, shown in Fig. 4 for the RFIM at high fields were investigated in pioneering pulsed-field magnetization measurements in $`\mathrm{Fe}_\mathrm{x}\mathrm{Zn}_{1\mathrm{x}}\mathrm{F}_2`$. Low temperature single spin flips and the phase boundary are shown in Fig. 5 and, interestingly, the behavior of the upper phase boundary appears to be different for $`x<x_e`$ and $`x>x_e`$. This is consistent with the differentiation of the behaviors observed in neutron scattering experiments above and below $`x_e`$.

In recent years, it has become clear that weak RFIM (small applied field) and strong RFIM regimes exist for $`x_p<x<x_e`$. The $`\mathrm{Fe}_{0.31}\mathrm{Zn}_{0.69}\mathrm{F}_2`$ system exhibits typical low-field behavior for $`H<1.5`$ T, with $`T_NT_c(H)`$ and $`T_NT_{eq}(H)`$ scaling as $`H^{2/\varphi }`$ with $`\varphi 1.4`$. At larger fields, however, the curvature for $`T_NT_{eq}(H)`$ changes to $`\varphi 3.4`$, a value close to that observed in spin glasses. Some of these features were suggested qualitatively. While the lower region has been shown to have antiferromagnetic long-range order upon ZFC, no long-range antiferromagnetic order is observed at the higher fields below $`T_{eq}(H)`$. Instead, spin-glass-like behavior is observed. This is clearly the same type of behavior observed for all fields at the percolation threshold concentration. The same type of distinctive low and high field behaviors have been observed for $`x`$ as large as $`0.60`$. The two regions are separated by an equilibrium boundary, as observed for $`x=0.56`$ (Fig. 6) and $`0.60`$, and it appears to decrease towards the $`H=0`$ boundary at finite temperature well below $`T_c(H)`$ and approach the phase transition line at finite field, a point separating the sharp transition observed at low field and the more glassy transition at higher fields. This also is consistent with specific heat peaks that are very sharp at low fields and quite rounded at high fields.

The distinction between high- and low-field behavior behavior is observed as well in ac susceptibility experiments. At low fields, there exists a single peak which seems to be associated with extremely slow dynamics, either from activated dynamics or at least power-law behavior with a very large dynamic exponent. There is little hysteresis between the ZFC and FC procedures at low $`H`$. At larger fields the peak splits in the ZFC procedure only, with a sharp peak at slightly lower temperatures than the broader peak. This splitting appears to be associated with the upper region of the phase diagram corresponding to the spin-glass-like behavior. These effects have not been investigated for $`x>x_p`$. The high field region for these concentrations is still an open area for research.

## 8 $`d=3`$ RFIM Equilibrium Critical Behavior for $`x>x_e`$

It was, of course, realized very early that the metastable domain walls at low magnetic concentrations took advantage of vacancies. What was not fully appreciated was that the Imry-Ma domain wall energy argument is not applicable when domain walls can to a great extent pass through vacancies, avoiding the energy cost of breaking magnetic bonds. With sufficient vacancies, i.e. for $`x<x_e`$, domain walls can take advantage of vacancies to such an extent that the domain wall energy can be insignificant. Interestingly, every experiment that could have detected low temperature hysteresis, particularly neutron scattering and capacitance experiments, was done for $`x0.72`$, which is below $`x_e=0.76`$. Higher concentrations were avoided since the generated random fields are small, resulting in relatively very narrow asymptotic random-field critical regions around $`T_c(H)`$. Nevertheless, concentrations well above $`x_e`$ are necessary to study the equilibrium critical behavior and require high magnetic fields and very fine temperature resolution.

Figure 2 shows the specific heat data for $`\mathrm{Fe}_{0.93}\mathrm{Zn}_{0.07}\mathrm{F}_2`$, measured with optical linear birefringence. The specific heat was also measured to demonstrate that, in agreement with theory and contrary to the so-called ‘trompe l’oeil’ phenomenology, the data from the both techniques yield the same critical parameters. An important advantage of the birefringence technique is its relative insensitivity to the concentration gradients that tend to smear the transition. Interestingly, the critical behavior appears to be very similar to that of lower concentrations where metastable domains dominate the scattering behavior. The specific heat was studied using Monte Carlo simulations based on the $`\mathrm{Fe}_{0.93}\mathrm{Zn}_{0.07}\mathrm{F}_2`$ system and the results are shown in Fig. 3. Although the simulations are not of sufficient resolution to extract the critical exponent, the similarities in the shapes of both the REIM and RFIM indicate the same qualitative change from the asymmetric cusp at $`H=0`$ to the nearly symmetric peaks at $`H>0`$. This result, however, is not seen in all MC simulations.

The neutron scattering experiments on $`\mathrm{Fe}_{0.93}\mathrm{Zn}_{0.07}\mathrm{F}_2`$ show no hysteresis below the phase transition, in stark contrast with samples with $`x<x_e`$. There is no evidence that domains form upon ZFC or FC in this concentration range and the line shapes are independent of the thermal cycling procedure, implying equilibrium conditions.

Neutron scattering experiments on this sample were difficult to analyze since the RFIM line shape is not known from theory. In general, two different scaling functions are involved with the form

$$\chi _s(q)=A^\pm \kappa ^{\eta 2}f(q/\kappa )+B^\pm \kappa ^{\overline{\eta }4}g(q/\kappa ).$$
(4)
with two independent sets of critical behavior exponents. However, simulations and high temperature series expansions strongly suggest a simpler scenario. The predictions are that the exponents are simply related, $`\overline{\eta }`$ being twice $`\eta `$, and the new scaling function $`g(q/\kappa )`$ is, to a good approximation, the square of $`f(q/\kappa )`$.

Not only were the universal RFIM critical parameters obtained in this experimental study, but the scaling analysis yielded the spin-spin correlation scaling function $`f(q/\kappa )`$. This scaling function is compared to several other known spin-spin correlation scaling functions in Fig. 1. Note that the $`d=3`$ RFIM one seems the furthest away from the MF Lorentzian of all the examples for $`T>T_c(H)`$, whereas for $`T>T_c(H)`$ the pure $`d=2`$ case is further than the $`d=3`$ RFIM, though both are very far from MF behavior.

The critical parameters for the RFIM with $`x>x_e`$ are shown in Table 4. Certainly more effort is needed to complete the experimental entries and to find reconciliation between the simulation and experimental results. Note that some sets of exponents from the simulations violate the Rushbrooke scaling relation in Eq. 3.

## 9 The Vacancy Percolation Threshold Concentration

The regions of low temperature nonequilibrium behavior and equilibrium behavior have been shown to be separated at relatively small $`H`$ by a nearly vertical sharp boundary at $`xx_e=0.755`$ in $`\mathrm{Fe}_\mathrm{x}\mathrm{Zn}_{1\mathrm{x}}\mathrm{F}_2`$ using Monte Carlo studies. Figure 7 shows simulations on lattices of size $`2\times L^3`$ with $`L=64`$ modeled closely after $`\mathrm{Fe}_\mathrm{x}\mathrm{Zn}_{1\mathrm{x}}\mathrm{F}_2`$. Hysteresis is observed upon FC and ZFC for $`x<x_e`$ but not above. The hysteresis for $`x<x_e`$ increases for larger lattices or slower thermal cycling, showing that it is not simply an artifact of the simulations not being run long enough. The concentration dividing equilibrium and nonequilibrium behavior is very close to or equal to the vacancy percolation threshold concentration $`x_e=0.755`$. Apparently, the percolation vacancy structure facilitates the formation of domain walls.

To further investigate this boundary, recent experiments have been done on a sample with $`x=0.76`$, just above $`x_e`$. No evidence of domains has been observed for small $`H`$. Since earlier experiments for $`x=0.72`$ in $`\mathrm{Fe}_\mathrm{x}\mathrm{Zn}_{1\mathrm{x}}\mathrm{F}_2`$ gave clear evidence for domain formation, including a reversal of the Bragg intensity curvature just below $`T_c(H)`$ upon ZFC, the boundary must be $`0.72<x_e<0.76`$, in agreement with the MC simulations. There have been no theoretical studies reported explaining the existence or nature of this boundary.

## 10 The Current Situation and Outlook

The recent measurement of the equilibrium critical behavior in the random-field Ising model has side-stepped the great difficulties encountered in the interpretation of data below the transition that are obfuscated by nonequilibrium phenomena in many studies at lower concentrations. Certainly, there is wide agreement that a phase transition exists. Although experiments for $`x>x_e`$ are much more difficult since the random-field region is very narrow, they are being done. Interestingly, the experimental results are not in agreement with much of the theory and simulation results, unlike the REIM and pure Ising model. A reliable characterization of some aspects of the $`d=3`$ RFIM universality class behavior remains to be completed. The phase diagram of $`\mathrm{Fe}_\mathrm{x}\mathrm{Zn}_{1\mathrm{x}}\mathrm{F}_2`$ has proven to be quite rich in detail. An important area of the phase diagram for which a good understanding is being developed is the large random field limit for $`x_p<x<x_e`$. With the the progress being made along these two lines of inquiry, a rather complete experimental characterization of the RFIM in dilute antiferromagnets seems near at hand.

Recent work has been supported by Department of Energy Grant No. DE-FG03-87ER45324.


text:
# Disorder and Two-Particle Interaction in Low-Dimensional Quantum Systems

## 1 Introduction

Until 1994, the theoretical and experimental research on transport in disordered systems clearly supported the scaling hypothesis of localization for non-interacting electrons . The systems studied usually fell into the predicted universality classes, and, if they didn’t, then they could be shown not to be generic . However, real electrons of course interact , and their interaction is of relevance for the transport properties of disordered systems , especially in 2D and 1D where screening  is less efficient than in 3D. The influence of weak interactions has been investigated extensively using perturbation theory and the perturbative renormalization group (RG) . One of the key results is that the lower critical dimension of the MIT is $`d_c^{}=2`$ as it is for non-interacting electrons. The application of the perturbative RG in 1D  has lead to the prediction that all thermodynamic states remain localized in the presence of repulsive many-body interactions.

Due to the persistent current problem  and new experiments on 2D electron systems  which show striking signatures of a metal-insulator transition, these theoretical considerations received a lot of renewed attention. In order to theoretically study the effects of the interplay between disorder and interactions, one should in principle solve a problem with an exponentially growing number of states in the Hilbert space with increasing system size. At present, this can be achieved only for a few particles in 1D  and very few particles in 2D . However, in 1994 Shepelyansky  proposed to simply look at two interacting particles (TIP) in a random environment. In particular, he suggested that the two particles would form pairs even for repulsive interactions such that the TIP pairs would have a larger localization length than the two single particles (SP) separately. Thus the interaction would lead to an enhanced possibility of transport . The perhaps even more surprising part of the prediction is that the TIP pairs will have a localization length $`\lambda _2`$ such that at pair energy $`E=0`$

$$\lambda _2U^2\lambda _{1}^{}{}_{}{}^{2},$$
(1)
where $`U`$ represents the onsite interaction strength and $`\lambda _1`$ is the SP localization length. Since $`\lambda _1105/W^2`$  in 1D, this implies large values of $`\lambda _2`$ for small disorders $`W`$.

The first numerical studies devoted to the TIP problem used the transfer-matrix method (TMM)  to investigate the proposed enhancement of the pair localization length $`\lambda _2`$ . The TMM of  contained an additional artificial infinitely-long-ranged interaction that tends to mask the onsite interaction . The TMM of  avoids this problem, but is restricted to small system sizes and results for localization lengths of $`\lambda _2300`$ had been deduced on systems of size $`M=100`$. Therefore, two of us studied the TIP problem by a different TMM  at large system size $`M300`$ and found that (i) the enhancement $`\lambda _2/\lambda _1`$ decreases with increasing $`M`$, (ii) the behavior of $`\lambda _2`$ for $`U=0`$ is equal to $`\lambda _1`$ in the limit $`M\mathrm{}`$ only, and (iii) for $`U0`$ the enhancement $`\lambda _2/\lambda _1`$ also vanishes completely in this limit. Consequently, we concluded  that the TMM applied to the TIP problem in 1D measures an enhancement of the localisation length which is due to the finiteness of the systems considered. Although Ref.  has been criticized , we emphasize that subsequent publications have shown  that there are no variants of TMM that reproduce Eq. (1). Furthermore, in a later numerical approach , based on Green function methods, Song and v. Oppen argue that our extrapolations for $`M\mathrm{}`$ were off by $`11\%`$ only, whereas the original TMM of  deviated by about a factor of $`3`$ . Thus while our criticized TMM results are valid, various other numerical investigations by other groups  as well as ourselves  convincingly demonstrated some enhancement. The reason for the failure of the TMM approach of  has been explained by Song and v. Oppen  by arguing that the TMM measures a localization length $`\lambda _f<\lambda _2`$ due to the cigar-shape geometry  of the TIP states.

Reliable numerical approaches to the TIP problem are nowadays based on the computation of the decay of the Green function . Other direct numerical approaches to the TIP problem have been based on the time evolution of wave packets , exact diagonalization , variants of level statistics  and analysis of multifractal properties , perturbative methods  and mappings to effective models . In these investigations an enhancement of $`\lambda _2`$ compared to $`\lambda _1`$ has been found as remarked above but the quantitative results tend to differ both from the analytical prediction in Eq. (1), and, albeit less, from each other. Furthermore, a check of the functional dependence of $`\lambda _2`$ on $`\lambda _1`$ is numerically very expensive since it requires very large system sizes $`M\lambda _2\lambda _1`$. Extensions of the original arguments have been proposed for TIP in 2D  and 3D , for TIP close to a Fermi sea , and for long-range interactions in 1D .

The basic idea leading to the prediction (1) is based on looking at the interaction matrix element between two eigenstates $`\psi _{kl}=\psi _k\psi _l`$ and $`\psi _{nm}=\psi _n\psi _m`$ of the non-interacting system . Here $`\psi _k,\psi _l,\psi _n,\psi _m`$ denote SP eigenstates localized with localization length $`\lambda _1`$ around sites $`k,l,n,m`$. For an onsite interaction  $`U_{j=1}^Nn_jn_j`$ (with $`n_{j\tau }`$ denoting the number operator at site $`j`$ for spin $`\tau `$) only states with $`|kl|\lambda _1`$, $`|nm|\lambda _1`$, $`|kn|\lambda _1`$, $`|lm|\lambda _1`$ will give significant contributions to the interaction matrix element

$`u`$ $`=`$ $`\psi _{kl}|U|\psi _{nm}`$ (2)
$`=`$ $`U{\displaystyle \underset{j=1}{\overset{N}{}}}\psi _k^{}(j)\psi _l^{}(j)\psi _n(j)\psi _m(j).`$
These conditions are illustrated in Fig. 1.

If one assumes  that the SP state is given as

$$\psi _k(j)\frac{1}{\sqrt{\lambda _1}}\mathrm{exp}\left(\frac{|jk|}{\lambda _1}+i\theta (j)\right)$$
(3)
with $`\theta (j)`$ a random phase, one finds  that the typical interaction matrix element has a magnitude of

$$u\lambda _1^{3/2}$$
(4)
since it is the sum of $`\lambda _1`$ random contributions of magnitude $`\lambda _1^2`$. Shepelyansky next calculated the decay rate $`\mathrm{\Gamma }`$ of a non-interacting eigenstate by means of Fermi’s golden rule $`\mathrm{\Gamma }U^2/\lambda _1t`$ . Since the typical hopping distance is of the order of $`\lambda _1`$ the diffusion constant is $`DU^2\lambda _1/t`$. Within a time $`\tau `$ the particle pair visits $`NU\lambda _1^{3/2}t^{1/2}\tau ^{1/2}`$ states. Diffusion stops when the level spacing of the visited states is of the order of the frequency resolution $`1/\tau `$. This determines the cut-off time $`\tau ^{}`$ and the corresponding pair-localization length is obtained as $`\lambda _2\sqrt{D\tau ^{}}(U/t)^2\lambda _1^2`$ in agreement with Eq. (1). Applicability of Fermi’s golden rule requires $`\mathrm{\Gamma }t/\lambda _1^2`$ which is equivalent to $`U^2\lambda _1/t^21`$. This is exactly the condition for an enhancement of $`\lambda _2`$ compared to $`\lambda _1`$. Alternatively, the model may be mapped to a random-matrix model (RMM) with entries chosen according to Eq. (4) .

## 2 Numerical results for the random-matrix model of TIP

The arguments presented in the last section are of a qualitatively nature and Eq. (1) must be checked for quantitative accuracy. Even before testing (1), it is already worthwhile to check the validity of (4) and the subsequent arguments or the RMM approach . In Ref. , it had been shown that the assumption of a Gaussian distribution of the matrix elements $`u`$ was oversimplified. The distribution showed long tails making the arithmetic average unsuitable to characterize the typical value. In Ref.  we have paid special attention to the exact dependence of $`u`$ on $`\lambda _1`$ and system size. To this end, we diagonalized the 1D Anderson model for a given $`M`$ and $`W`$ and computed $`u`$ by averaging over all suitable states and many disorder configurations. We showed that due to the strongly non-Gaussian distribution of $`u`$, one should rather use the logarithmic average than the arithmetic average as the typical value for the computation of $`u(\lambda _1)`$. But whereas the arithmetic average  gives $`u_{\mathrm{abs}}\lambda _1^{1.5}`$, the typical value obeys $`u_{\mathrm{typ}}\lambda _1^{1.95}`$. Following the arguments above, this would imply $`\lambda _2\lambda _1^{1.1}`$, i.e., a very small enhancement. We emphasize that this result does not mean that there is no enhancement of the localization length. Rather, the results of Ref.  indicate that the arguments of Ref.  capture the physics, but only in a somewhat simplified form. One step towards a better agreement between the analytical and the numerical approaches is to take into account the energy denominators in the computation of $`u`$, e.g., to consider only interaction matrix elements for states whose energy spacings are of the order of $`U`$ or smaller . In this case we find that there is a slight decrease in the value of the typical exponent and correspondingly a slight increase in TIP delocalization yielding $`\lambda _2\lambda _1^{1.4\pm 0.2}`$. This suggests that higher orders in perturbation theory than the first order RMM approach  are important. Furthermore, the exponent $`1.4\pm 0.2`$ is in reasonable agreement with previous results in the literature .

## 3 RMM approach for toy models

In this section we show that a naive application the RMM approach may give qualitatively incorrect results even if the RMM contains the correct dependence of the matrix elements on the disorder strength. To this end we consider two toy models which seem to be closely related to the TIP problem. For these models, viz. Anderson models of localization with additional perturbing random potentials, we show that mapping onto RMMs and estimating the localization length by Fermi’s golden rule leads to an incorrect enhancement of the localization length.

### 3.1 2D Anderson model with perturbation on a line

The first example is set up to lead to the same RMM as the TIP problem. It consists of the usual 2D Anderson model perturbed by an additional weak random potential of strength $`U`$ at the diagonal $`x=y`$ in real space. Since the perturbation increases the width of the disorder distribution at the diagonal we expect it to decrease the localization length. We map the model onto an RMM following the arguments for the TIP problem sketched in Sec. 1. Again, the eigenstates of the unperturbed system are localized with a localization length $`\lambda _1`$ and approximately given by

$`\psi _n(x,y){\displaystyle \frac{1}{\lambda _1}}\mathrm{exp}\left[{\displaystyle \frac{|𝐫𝐫_n|}{\lambda _1}}+i\theta _n(𝐫)\right]`$ (5)
where $`𝐫=(x,y)^T`$ is the coordinate vector of the particle and $`\theta `$ is again assumed to be a random phase. The Hamiltonian of the 2D perturbed Anderson model differs from the TIP Hamiltonian in two points: (i) the diagonal elements (given by the random potential) are independent random numbers instead of being correlated as in the TIP problem and (ii) the perturbing potential $`U(x,x)[U,U]`$ at each diagonal site is random instead of having a definite sign and modulus $`U`$ as in the TIP problem. However, none of these points enters the mapping procedure outlined in Sec. 1. Thus, we find that the perturbation couples each state close to the diagonal ($`|x_ny_n|<\lambda _1`$) to $`𝒪(\lambda _1^2)`$ other such states. The interaction matrix element is again a sum of $`𝒪(\lambda _1)`$ terms of magnitude $`U/\lambda _1^2`$ and random phases and as before $`uU\lambda _1^{3/2}`$. Consequently, our toy model is mapped onto exactly the same RMM as TIP in a random potential. Therefore, the resulting localization length along the diagonal is also given by Eq. (1). We thus arrive at the surprising conclusion, that adding a weak random potential at the diagonal of a 2D Anderson model leads to an enormous enhancement of the localization length along this diagonal, in contradiction to the expectation expressed above, viz. that increasing disorder leads to stronger localization.

As for the TIP case  we now numerically check whether the relation $`uU\lambda _1^{3/2}`$ between the coupling matrix element $`u`$ and the localization length $`\lambda _1`$ of the unperturbed system is correctly described by the RMM. Since in 2D a simple analytic formula for the dependence of $`\lambda _1`$ on the disorder $`W`$ does not exist, we first compute estimates $`\lambda _1(M)`$ for quasi-1D strips of finite strip width $`M`$ with $`1\%`$ accuracy by TMM. In Fig. 2, we show data of $`\lambda _1(M)`$ as a function of $`W`$. In the following, we take $`\lambda _1(50)`$ to compute the coupling matrix elements.

Next, we calculate both the arithmetic average $`u_{\mathrm{abs}}=|u|`$ and the logarithmic average $`u_{\mathrm{typ}}=\mathrm{exp}[\mathrm{log}(|u|)]`$ for different values of $`W`$ and various $`M\times M`$ squares. Disorder averaging is over $`20`$ samples and we study $`u_{\mathrm{abs}}`$ and $`u_{\mathrm{typ}}`$ as functions of $`\lambda _1(M)`$. We emphasize that instead of the well-known extrapolations of $`\lambda _1(M)`$ to infinite system size by means of FSS , we take the finite-size approximants $`\lambda _1(M)`$ on purpose, since we compute $`\lambda _2`$ also for comparable finite sizes only.

As for the TIP model  the distribution $`P_\mathrm{o}(u)`$ of the (off-diagonal) coupling matrix elements is strongly non-Gaussian, suggesting that $`u_{\mathrm{typ}}`$ rather than $`u_{\mathrm{abs}}`$ is the relevant quantity. The results for $`u_{\mathrm{abs}}`$ and $`u_{\mathrm{typ}}`$ are presented in Fig. 3.

The dependence of $`u_{\mathrm{abs}}`$ on $`\lambda _1(M)`$ for $`2\lambda _1(M)12`$ follows $`u_{\mathrm{abs}}\lambda _1(M)^{1.6\pm 0.1}`$ in agreement with the RMM value of 3/2 and with Ref. . Furthermore, here we also have $`u_{\mathrm{typ}}\lambda _1(M)^{1.5\pm 0.1}`$. We note that the change of the slopes of $`u_{\mathrm{abs}}`$ and $`u_{\mathrm{typ}}`$ at $`\lambda _1(M)M/2`$ is entirely due to the finite sample sizes .

Consequently, in contrast to the TIP problem the RMM model for the 2D perturbed Anderson model of localization contains the correct dependence of the coupling matrix elements on the localization length of the unperturbed system, but still it leads to an incorrect enhancement of the localization length along the diagonal.

### 3.2 1D Anderson model with perturbation

An even more striking contradiction can be obtained for an 1D Anderson model of localization. The eigenstates are again given by Eq. (3) with $`\lambda _1`$ known from second order perturbation theory  and numerical calculations  to vary as $`\lambda _1t^2/W^2`$ for small disorder. We now add a weak random potential of strength $`U`$ at all sites. Since the result is obviously an 1D Anderson model with a slightly higher disorder strength $`W_u>W`$ the localization length will be reduced, $`\lambda (U)t^2/W_u^2`$.

The mapping onto an RMM can be performed in complete analogy to the TIP problem and the 2D Anderson model discussed above. The perturbing potential leads to transitions between the unperturbed eigenstates $`\psi _n`$. Each such state is now coupled to $`𝒪(\lambda _1)`$ other states by coupling matrix elements $`\psi _n|U|\psi _n^{}`$ with magnitude $`uU\lambda _1^{1/2}`$ since we sum over $`\lambda _1`$ contributions with magnitude $`U/\lambda _1`$ and supposedly random phases.

The application of Fermi’s golden rule in this 1D case leads to a diffusion constant $`DU^2\lambda _1^2/t`$. The number of states visited within a time $`\tau `$ is $`NU\lambda _1t^{1/2}\tau ^{1/2}`$. Again, diffusion stops at a time $`\tau ^{}`$ when the level spacing of the states visited equals the frequency resolution. This gives $`\tau ^{}U^2\lambda _1^2/t^3`$. The localization length $`\lambda `$ of the perturbed system thus reads $`\lambda \sqrt{D\tau ^{}}U^2\lambda _1^2`$ as in Eq. (1), in contradiction to the correct result.

Again we numerically check the relation between $`u_{\mathrm{abs}}`$ and $`u_{\mathrm{typ}}`$ and the unperturbed localization length $`\lambda _1`$. In Fig. 4, we show results obtained for chains with various lengths and $`50`$ disorder configurations for each $`W`$. $`\lambda _1`$ is computed by TMM.

For $`10\lambda _1250`$, $`u_{\mathrm{abs}}`$ varies as $`\lambda _1^{0.48\pm 0.10}`$ as predicted above. $`u_{\mathrm{typ}}`$ varies as $`\lambda _1^{0.59\pm 0.10}`$. Both variations are compatible with the RMM value of $`1/2`$ for the exponent. Again we need $`\lambda _1<M/2`$ in order to suppress finite size effects.

Consequently, although the RMM model for the 1D perturbed Anderson model of localization contains the correct dependence of the coupling matrix elements on the localization length of the unperturbed system, it still leads to an incorrect enhancement of the localization length.

## 4 Application of the block-scaling picture to toy models

Let us now discuss the relation of these results to Imry’s block-scaling picture (BSP)  for the TIP problem. In this approach one considers blocks of linear size $`\lambda _1`$ and calculates the dimensionless pair conductance on that scale,

$$g_2\frac{u^2}{\mathrm{\Delta }^2},$$
(6)
where $`u`$ represents the typical interaction-induced coupling matrix element between states in neighboring blocks and $`\mathrm{\Delta }t/\lambda _1^2`$ is the level spacing within the block. If the typical coupling matrix element depends on $`\lambda _1`$ as $`uU\lambda _1^\alpha `$ the pair conductance obeys

$$g_2(U/t)^2\lambda _1^{42\alpha }.$$
(7)
For the 2D Anderson model with perturbation considered above, the BSP can be applied analogously. Again, we consider blocks of linear size $`\lambda _1`$ and compute the typical perturbation-induced matrix elements between these blocks. We then find that according to the BSP the conductance of a 2D Anderson model with additional weak perturbing potential along the diagonal is given by Eq. (6). Using $`\alpha =1.5\pm 0.1`$ as obtained above from the numerical data for $`u_{\mathrm{abs}}`$ and $`u_{\mathrm{typ}}`$, we then have $`g_2(U/t)^2\lambda _1`$. Thus we conclude that the BSP does not work for our 2D toy model, because it yields the same unphysical result as the RMM approach of section 3.1.

Let us also apply the BSP to the 1D toy example. The level spacing in a 1D block of size $`\lambda _1`$ is $`\mathrm{\Delta }t/\lambda _1`$, and the coupling matrix element between states in neighboring blocks is $`uU\lambda _1^{1/2}`$. Thus, the conductance of the perturbed system on a scale $`\lambda _1`$ is obtained as $`g_2(U/t)^2\lambda _1`$. For large $`\lambda _1`$ this again contradicts the correct result, viz. a decrease of the conductance compared to the unperturbed system. Thus, the BSP applied to the two toy models introduced in the present work gives the same qualitatively incorrect results for the localization properties as the RMM. This is not surprising since the only ingredients of the BSP are the intra-block level spacing $`\mathrm{\Delta }t/\lambda _1^2`$ and the inter-block coupling matrix elements $`u`$ which also enter the RMM.

## 5 Conclusions

We have presented two toy models which seem to be closely related to the TIP problem. For these toy models the usual analytical arguments given to support the delocalization of TIP, viz. the RMM and the BSP do not work. However, the large-scale numerical simulations  have convincingly shown that an enhancement of the two-particle localization length due to the interaction exists, even though the detailed results are more complicated than the original prediction (1). This leads, of course, to the question, under which conditions the RMM mapping and the BSP give the correct result, at least qualitatively. While a general answer to this question is not known, it has been suggested  that the difference between the TIP and our toy models is an additional symmetry in the TIP problem.

In summary, the two examples suggest that additional physical insight is needed before applying the RMM. In addition, we expect that taking into account the energies of the states as in Ref.  for TIP will result in a reduced enhancement, i.e., a smaller value of $`\alpha `$, in the analytical predictions. This will in turn give a better agreement with the numerically determined dependence of the TIP localization length on $`\lambda _1`$.


text:
# Dynamical models for dusty disk galaxies

## 1. Introduction

The knowledge of the dynamics of disk galaxies is essential in order to understand their structure and history. Unfortunately, disk galaxies are difficult systems to model dynamically, for several reasons. One of them is the presence of a large amount of interstellar dust, which obscures the light along the lines-of-sight. Using extended radiative transfer models it is nowadays possible to recover quite accurately the three-dimensional light and dust distribution in disk galaxies (Kylafis & Bahcall 1987, Xilouris et al. 1999). But also the observed kinematics are affected by dust obscuration. Indeed, each element along a line-of-sight carries its own kinematic information, and the projected kinematics are a weighted mean of all these contributions. We adopt the technique outlined in Baes et al. (2000a,b) in order to investigate in detail the effects of dust extinction on the mean projected velocity $`\overline{v}_p`$ and the projected velocity dispersion $`\sigma _p`$.

## 2. The modelling and results

We adopt a galaxy model which consists of a double exponential disk and a de Vaucouleurs bulge. We construct a dynamical model (i.e. a potential and a phase-space distribution function) for this galaxy. We choose a potential that gives rise to a flat rotation curve and represents a halo-disk structure (Batsleer & Dejonghe 1994). Using the Quadratic Programming modelling procedure (Dejonghe 1989) we then construct a two-integral distribution function that is consistent with the light density. We add a double exponential dust disk to this model. Finally, the dust-affected $`\overline{v}_p`$ and $`\sigma _p`$ can be calculated for various values of the inclination and optical depth.

For galaxies which are face-on or moderately inclined, the effects of dust extinction on $`\overline{v}_p`$ and $`\sigma _p`$ are negligibly small. In the edge-on case, the dust-affected $`\overline{v}_p`$-profile tends to apparent solid body rotation, as we only see the stars moving on the outer near edge of the disk. In meanwhile, the projected dispersion decreases drastically as a function of optical depth for the inner lines-of-light, as dust obscuration strongly reduces the contribution of the high random motions of the bulge stars. Both effects are critically dependent on inclination, and they are already much weaker for galaxies which are only a few degrees from exactly edge-on (see also Bosma et al. 1992).

## 3. Conclusion: dynamical modelling of disk galaxies

From our results it is clear that the effects of dust obscuration on $`\overline{v}_p`$ and $`\sigma _p`$ are negligible for moderately inclined galaxies. Hence it is quite safe to neglect dust extinction in the interpretation of projected kinematics. This leads us to propose the following strategy to construct dynamical models for disk galaxies. Intermediately inclined disks are the best choice, as spectra at different position angles will then show different projections of the velocity ellipsoid.

First, one should determine the three-dimensional light distribution of the galaxy, using deprojection techniques which take the dust into account. The accuracy of the results can be tested by comparing models in different wavebands with the galactic extinction curve (Xilouris et al. 1999) or by comparing the derived extinction profile with FIR/submm emission (Alton et al. 2000). Then, a set of potentials which are consistent with the rotation curve and the light distribution need to be determined. For each potential a three-integral model can be constructed. Input for the fit should be the light density and the projected kinematics along (at least) both major and minor axes. The goodness of fit of the different models can then be used to constrain the set of possible potentials, which will reveal the mass distribution in the galaxy. The velocity field can then be analysed, in particular the behaviour of the velocity ellipsoid. This can shed a light on the mechanism responsible for the dynamical history of the disk (Jenkins & Binney 1990, Gerssen et al. 1997, 2000).

## References

Alton P. B., Xilouris E. M., Bianchi S., Davies J., Kylafis N., 2000, A&A, 356, 795

Baes M., Dejonghe H., 2000, MNRAS, 313, 153

Baes M., Dejonghe H., De Rijcke S., 2000, MNRAS, in press (astro-ph/0006110)

Batsleer P., Dejonghe H., 1994, A&A, 287, 43

Bosma A., Byun Y., Freeman K. C., Athanassoula E., 1992, ApJ, 400, L21

Dejonghe H., 1989, ApJ, 343, 113

Gerssen J., Kuijken K., Merrifield M. R., 1997, MNRAS, 288, 618

Gerssen J., Kuijken K., Merrifield M. R., 2000, astro-ph/0005454

Jenkins A., Binney J., 1990, MNRAS, 245, 305

Kylafis N. D., Bahcall J. N., 1987, ApJ, 317, 637

Xilouris E. M., Byun Y. I., Kylafis N. D., Paleologou E. V., Papamastorakis J., 1999, A&A, 344, 868