text:
# Microscopic versus mesoscopic local density of states in one-dimensional localization

\[

## Abstract

We calculate the probability distribution of the local density of states $`\nu `$ in a disordered one-dimensional conductor or single-mode waveguide, attached at one end to an electron or photon reservoir. We show that this distribution does not display a log-normal tail for small $`\nu `$, but diverges instead $`\nu ^{1/2}`$. The log-normal tail appears if $`\nu `$ is averaged over rapid oscillations on the scale of the wavelength. There is no such qualitative distinction between microscopic and mesoscopic densities of states if the levels are broadened by inelastic scattering or absorption, rather than by coupling to a reservoir.

\]

Localization of wave functions by disorder can be seen in the fluctuations of the density of states, provided the system is probed on a sufficiently short length scale . The local density of states (LDOS) of electrons can be probed using the tunnel resistance of a point contact  or the Knight shift in nuclear magnetic resonance , while the LDOS of photons determines the rate of spontaneous emission from an atomic transition . (In the photonic case one can study the effects of localization independently from those of interactions.)

For each length scale $`\delta `$ characteristic for the resolution of the probe, one can introduce a corresponding LDOS $`\nu _\delta `$. It is necessary that $`\delta `$ is less than the localization length, in order to be able to see the effects of localization — the hallmark  being the appearance of logarithmically normal tails $`\mathrm{exp}(\mathrm{constant}\times \mathrm{ln}^2\nu _\delta )`$ in the probability distribution $`P(\nu _\delta )`$.

Much of our present understanding  of this problem in a wire geometry builds on the one-dimensional (1D) solution of Altshuler and Prigodin . In the simplest case one has a single-mode wire which is closed at one end and attached at the other end to an electron reservoir. (The optical analogue would be a single-mode waveguide that can radiate into free space from one end.) In 1D the localization length equals twice the mean free path $`l`$, which is assumed to be large compared to the wavelength $`\lambda `$. One can then distinguish the microscopic LDOS $`\nu =\nu _\delta `$ for $`\delta \lambda `$, and the mesoscopic LDOS $`\stackrel{~}{\nu }=\nu _\delta `$ for $`\lambda \delta l`$. While $`\nu `$ oscillates rapidly on the scale of the wavelength, $`\stackrel{~}{\nu }`$ only contains the slowly varying envelope of these oscillations. Altshuler and Prigodin calculated the distribution $`P(\stackrel{~}{\nu })`$ and surmised that $`P(\nu )`$ would have the same log-normal tails. We will demonstrate that this is not the case for the small-$`\nu `$ asymptotics.

The calculation of Ref.  was based on the Berezinskii diagram technique , which reconstructs the probability distribution from its moments. (An alternative approach , using the method of supersymmetry, also proceeds via the moments.) An altogether different scattering approach has been proposed by Gasparian, Christen, and Büttiker , and more recently by Pustilnik . We have pursued this approach and arrive at a relation between $`\nu `$, $`\stackrel{~}{\nu }`$ and reflection coefficients. This allows a direct calculation of the distributions. We find that $`P(\nu )`$ and $`P(\stackrel{~}{\nu })`$ have the same log-normal tail for large densities, but the asymptotics for small $`\nu `$ and $`\stackrel{~}{\nu }`$ is completely different. The strong fluctuations of $`\nu `$ on the scale of the wavelength lead to a divergence $`P(\nu )\nu ^{1/2}`$ for $`\nu 0`$, while the distribution of the envelope vanishes, $`P(\stackrel{~}{\nu })0`$ for $`\stackrel{~}{\nu }0`$. This qualitative difference between microscopic and mesoscopic LDOS is a feature of an open system. Both $`P(\nu )`$ and $`P(\stackrel{~}{\nu })`$ vanish for small densities if the wire is closed at both ends and the levels are broadened by inelastic scatterers (for electrons) or absorption (for photons).

We consider a 1D wire and relate the microscopic and mesoscopic LDOS at energy $`E`$ and at a point $`x=0`$ to the reflection amplitudes $`r_\mathrm{R}`$, $`r_\mathrm{L}`$ from parts of the wire to the right and to the left of this point. The Hamiltonian is $`H=(\mathrm{}^2/2m)^2/x^2+V(x)`$ for non-interacting electrons. (For photons of a single polarization we would consider the differential operator of the scalar wave equation.) We will put $`\mathrm{}=1`$ for convenience of notation. We start from the relation between the LDOS and the retarded Green function,

$`\nu =\pi ^1\mathrm{Im}G(0),`$ (1)
$`(E+i\eta H)G(x)=\delta (x),`$ (2)
with $`\eta `$ a positive infinitesimal. We assume weak disorder ($`kl1`$, with $`k=2\pi /\lambda `$ the wavenumber), so that we can expand the Green function in scattering states in a small interval around $`x=0`$,

$$G(x)=c_\mathrm{L}(e^{ikx}+r_\mathrm{L}e^{ikx})\theta (x)+c_\mathrm{R}(e^{ikx}+r_\mathrm{R}e^{ikx})\theta (x).$$
(3)
(The function $`\theta (x)=1`$ for $`x>0`$ and $`0`$ for $`x<0`$.) The coefficients $`c_\mathrm{L}`$ and $`c_\mathrm{R}`$ are related by the requirement that the Green function be continuous at $`x=0`$,

$$c_\mathrm{L}(1+r_\mathrm{L})=c_\mathrm{R}(1+r_\mathrm{R}).$$
(4)
Substitution of Eq. (3) into Eq. (2) gives a second relation between $`c_\mathrm{L}`$ and $`c_\mathrm{R}`$, from which we deduce that

$$G(0)=\frac{(1+r_\mathrm{L})(1+r_\mathrm{R})}{iv(1r_\mathrm{R}r_\mathrm{L})},$$
(5)
with $`v`$ the velocity. Using Eq. (1) we arrive at the key relation between the microscopic LDOS and the reflection coefficients,

$$\nu =(\pi v)^1\mathrm{Re}(1+r_\mathrm{L})(1r_\mathrm{R}r_\mathrm{L})^1(1+r_\mathrm{R}).$$
(6)
In order to perform the local spatial average that gives the mesoscopic LDOS $`\stackrel{~}{\nu }`$, we use that the reflection coefficients oscillate on the scale of the wavelength. If we shift $`x_0`$ slightly away from the origin to a point $`x^{}`$, one has $`r_\mathrm{L}e^{2ikx^{}}r_\mathrm{L}`$ and $`r_\mathrm{R}e^{2ikx^{}}r_\mathrm{R}`$. The product $`r_\mathrm{R}r_\mathrm{L}`$, however, does not display these oscillations—only this combination should be retained. Hence

$$\stackrel{~}{\nu }=(\pi v)^1\mathrm{Re}(1+r_\mathrm{R}r_\mathrm{L})(1r_\mathrm{R}r_\mathrm{L})^1.$$
(7)
In what follows we will measure $`\nu `$ and $`\stackrel{~}{\nu }`$ in units of $`\nu _0=(\pi v)^1`$, which is the macroscopic density of states and the ensemble average of $`\nu `$, $`\stackrel{~}{\nu }`$.

Let us now demonstrate the power of the two simple relations (6) and (7). We take the wire open at the left end and study the density at a distance $`L`$ from this opening. At the right end the wire is assumed to be closed, giving rise to a reflection coefficient $`r_\mathrm{R}=\mathrm{exp}(i\varphi _\mathrm{R})`$ with uniformly distributed phase $`\varphi _\mathrm{R}`$ in the interval $`(0,2\pi )`$. The reflection coefficient $`r_\mathrm{L}=\sqrt{R}\mathrm{exp}(i\varphi _\mathrm{L})`$ is parameterized through the uniformly distributed phase $`\varphi _\mathrm{L}`$ and the reflection probability $`R`$ in the interval $`(0,1)`$. The ratio

$$u=(1+R)(1R)^1$$
(8)
has the probability distribution

$$\rho (u)=\frac{e^{s/4}}{\sqrt{\pi }(2s)^{3/2}}\underset{\mathrm{arcosh}u}{\overset{\mathrm{}}{}}𝑑z\frac{ze^{z^2/4s}}{(\mathrm{cosh}zu)^{1/2}},$$
(9)
with $`s=L/l`$ and $`l`$ the mean free path for backscattering. The mesoscopic LDOS (7) can be written in terms of the variables $`u`$ and $`\varphi =\varphi _\mathrm{L}+\varphi _\mathrm{R}`$,

$$\stackrel{~}{\nu }=\left(u\sqrt{u^21}\mathrm{cos}\varphi \right)^1.$$
(10)
Averaging first over $`\varphi `$ we find

$$P_{\mathrm{open}}(\stackrel{~}{\nu })=\frac{\stackrel{~}{\nu }^{3/2}}{\pi \sqrt{2}}\underset{a}{\overset{\mathrm{}}{}}𝑑u\frac{\rho (u)}{\sqrt{ua}},a=\frac{1}{2}(\stackrel{~}{\nu }+\stackrel{~}{\nu }^1).$$
(11)
The subsequent integration with Eq. (9) yields

$$P_{\mathrm{open}}(\stackrel{~}{\nu })=\frac{\stackrel{~}{\nu }^{3/2}e^{s/4}}{2\sqrt{\pi s}}\mathrm{exp}\left(\frac{1}{4s}\mathrm{ln}^2\stackrel{~}{\nu }\right).$$
(12)
The distribution function (12) is the celebrated result of Altshuler and Prigodin . It displays log-normal tails for both large and small values of $`\stackrel{~}{\nu }`$. Indeed, the two tails are linked by the functional relation

$$P(1/\stackrel{~}{\nu })=\stackrel{~}{\nu }^3P(\stackrel{~}{\nu }).$$
(13)
This relation follows directly from Eq. (11) and hence requires only a uniformly distributed phase $`\varphi `$, regardless of the distribution function $`\rho (u)`$ of the reflection probability. As we will now show, such a relation does not hold in general for the microscopic LDOS $`\nu `$, and the asymptotics of its distribution function for small and large values of $`\nu `$ can be entirely different.

The calculation is facilitated by the fact that $`\nu `$ is related to $`\stackrel{~}{\nu }`$ by

$$\nu =2\stackrel{~}{\nu }\mathrm{cos}^2(\varphi _\mathrm{R}/2)\text{if}|r_\mathrm{R}|=1.$$
(14)
Moreover, $`\stackrel{~}{\nu }`$ is statistically independent of $`\varphi _\mathrm{R}`$ because the latter enters $`\stackrel{~}{\nu }`$ only in combination with $`\varphi _\mathrm{L}`$, which itself is uniformly distributed. The distribution of the microscopic LDOS hence follows directly from Eq. (12),

$$P_{\mathrm{open}}(\nu )=\frac{\nu ^{3/2}e^{s/4}}{\pi \sqrt{2\pi s}}_0^1\frac{dt}{\sqrt{1t}}\mathrm{exp}\left(\frac{1}{4s}\mathrm{ln}^2\nu \right),$$
(15)
where we substituted $`t=\mathrm{cos}^2(\varphi _\mathrm{R}/2)`$. The asymptotic behavior is

$`P_{\mathrm{open}}(\nu )={\displaystyle \frac{\mathrm{exp}(3s/4)}{2^{1/2}\pi }}\nu ^{1/2},\nu e^s,`$ (17)
$`P_{\mathrm{open}}(\nu )={\displaystyle \frac{2^{1/2}\mathrm{exp}(s/4)}{s^{1/2}\pi ^{3/2}}}\nu ^{3/2},e^s\nu e^s,`$ (18)
$`P_{\mathrm{open}}(\nu )={\displaystyle \frac{\mathrm{exp}[s/4\mathrm{ln}^2(\nu /2)/4s]}{\pi \nu ^{3/2}\mathrm{ln}^{1/2}(\nu /2)}},\nu e^s.`$ (19)
In the second and third region this is similar to the behavior of $`P_{\mathrm{open}}(\stackrel{~}{\nu })`$ in Eq. (12). In the region of the smallest densities, however, $`P_{\mathrm{open}}(\nu )`$ is not log-normal like $`P_{\mathrm{open}}(\stackrel{~}{\nu })`$ but diverges $`\nu ^{1/2}`$.

The different tails arise from two qualitatively different mechanisms that produce small values of $`\nu `$ and $`\stackrel{~}{\nu }`$. For the mesoscopic LDOS this requires remoteness of $`E`$ from the eigenvalues of wave functions localized within a localization length around $`x_0`$. As a consequence, $`P(\stackrel{~}{\nu })`$ is intimately linked to the distribution function of resonance widths . Small values of the microscopic LDOS $`\nu `$ are attained at nodes of the wave function which solves the wave equation with open boundary conditions, independent of the energy. The nodes are completely determined by the small-scale structure of the wave function, which is a real standing wave $`\mathrm{cos}(kx+\alpha )`$ with random phase $`\alpha `$ . \[We recognize the square of this wave amplitude in Eq. (14).\] The resulting $`\nu ^{1/2}`$ divergence of the probability distribution has the same origin as in the Porter-Thomas distribution for chaotic wave functions .

The two distributions for the open wire are plotted in Fig. 1, together with the result of a numerical simulation in which the Green function inside the wire is calculated recursively . The comparison of theory and numerics is free of any adjustable parameter—the velocity was taken from the dispersion relation, and the mean free path was obtained from the disorder strength within the Born approximation.

We now show that this qualitative difference between the microscopic and mesoscopic LDOS is absent in a closed wire. If the wire is decoupled from the reservoir we need another source of level broadening to regularize the delta functions in the LDOS. Following Ref. , we will retain a finite imaginary part $`\eta `$ of the energy, corresponding to spatially uniform absorption (for photons) or inelastic scattering (for electrons), with rate $`2\eta `$. Eqs. (6) and (7) still hold provided $`\eta E`$. The reflection coefficients can be written as

$$r_{\mathrm{R},\mathrm{L}}=\sqrt{R_{\mathrm{R},\mathrm{L}}}e^{i\varphi _{\mathrm{R},\mathrm{L}}},$$
(20)
where $`\varphi _\mathrm{R}`$ and $`\varphi _\mathrm{L}`$ are uniformly distributed phases and $`R_\mathrm{R}`$, $`R_\mathrm{L}`$ are independent reflection probabilities. In an infinitely long wire they have the same distribution

$$\rho (R)=\frac{\omega e^\omega }{(1R)^2}\mathrm{exp}\left[\omega (1R)^1\right],\omega =4\eta l/v.$$
(21)
After elimination of the phases the distribution of the mesoscopic LDOS takes again the form (11), where $`u`$ now stands for the combination

$$u=(1+R_\mathrm{R}R_\mathrm{L})(1R_\mathrm{R}R_\mathrm{L})^1.$$
(22)
Eq. (21) implies for $`u`$ the distribution

$$\rho (u)=\omega ^2\left(1\frac{}{\omega }\right)e^{\omega (u1)}K_0\left(\omega \sqrt{u^21}\right).$$
(23)
The resulting distribution function of the mesoscopic LDOS is

$`P_{\mathrm{closed}}(\stackrel{~}{\nu })={\displaystyle \frac{\omega ^2\stackrel{~}{\nu }^{3/2}}{\pi \sqrt{2}}}{\displaystyle _a^{\mathrm{}}}{\displaystyle \frac{du}{\sqrt{ua}}}e^{\omega (u1)}`$ (24)
$`\times \left[uK_0\left(\omega \sqrt{u^21}\right)+\sqrt{u^21}K_1\left(\omega \sqrt{u^21}\right)\right],`$ (25)
with $`a`$ defined in Eq. (11). It vanishes for small densities as

$$P_{\mathrm{closed}}(\stackrel{~}{\nu })=2^{1/2}\omega \stackrel{~}{\nu }^2\mathrm{exp}(\omega /\stackrel{~}{\nu }),\stackrel{~}{\nu }\omega .$$
(26)
This should be compared with the known distribution

$$P_{\mathrm{closed}}(\nu )=\sqrt{\frac{2\omega }{\pi }}\nu ^{3/2}\mathrm{exp}\left[\omega \frac{1}{2}\omega (\nu +\nu ^1)\right]$$
(27)
of the microscopic LDOS. In contrast to the open wire, both distributions vanish for $`\nu `$, $`\stackrel{~}{\nu }0`$. This is illustrated in Fig. 2, which compares the analytical predictions to numerical data obtained by diagonalization of a Hamiltonian. The comparison is again free of any adjustable parameter.

We note in passing that the asymptotic behavior (26) differs from the asymptotic behavior

$$P_{\mathrm{closed}}(\stackrel{~}{\nu })\frac{1}{4}(\pi \omega )^{1/2}\stackrel{~}{\nu }^{3/2}\mathrm{exp}(\pi ^2\omega /16\stackrel{~}{\nu }),$$
(28)
given in Ref.  for $`\omega 1`$. There the distribution function was reconstructed from the leading asymptotics of the moments $`lim_{\omega 0}\stackrel{~}{\nu }^n=\omega ^{1n}n!/(2n1)`$. This would be a valid procedure if the distribution would depend only on the product $`\omega \stackrel{~}{\nu }`$ in the limit $`\omega 0`$, which it does not. The subleading terms of the moments have to be included for $`\stackrel{~}{\nu }\omega `$. Indeed, our distribution function has the same leading asymptotics of the moments, but has a different functional form. This illustrates the potential pitfalls of the restoration procedure which are circumvented by our direct method.

In conclusion, we have given exact results for the distributions of the local densities of states in one-dimensional localization, contrasting the microscopic length scale (below the wavelength) and mesoscopic length scale (between the wavelength and the mean free path). Contrary to expectations in the literature, the log-normal asymptotics at small densities applies only to the mesoscopic LDOS $`\stackrel{~}{\nu }`$, while the distribution of the microscopic LDOS $`\nu `$ diverges $`\nu ^{1/2}`$ for $`\nu 0`$. This is of physical significance because many of the local probes act on atomic degrees of freedom and hence measure $`\nu `$ rather than $`\stackrel{~}{\nu }`$. The strong length scale dependence of the LDOS disappears if the electrons (or photons) are scattered inelastically (or absorbed) before reaching the reservoir. Both $`P(\nu )`$ and $`P(\stackrel{~}{\nu })`$ then have an exponential cutoff at small densities.

It is an interesting open problem whether the qualitative distinction between $`\nu `$ and $`\stackrel{~}{\nu }`$ in an open wire carries over to the quasi-one dimensional geometry with $`N>1`$ modes. An analytic theory could build on the multi-channel generalization of Eq. (6),

$$\nu =\mathrm{Re}\mathrm{Tr}\widehat{M}(1+\widehat{r}_\mathrm{L})(1\widehat{r}_\mathrm{R}\widehat{r}_\mathrm{L})^1(1+\widehat{r}_\mathrm{R}).$$
(29)
Now $`\widehat{r}_\mathrm{L}`$ and $`\widehat{r}_\mathrm{R}`$ are $`N\times N`$ reflection matrices and the matrix $`\widehat{M}_{nm}=2(\pi A)^1(v_nv_m)^{1/2}\mathrm{sin}(\stackrel{}{q}_n\stackrel{}{r}_0)\mathrm{sin}(\stackrel{}{q}_m\stackrel{}{r}_0)`$ contains the weights of the $`N`$ scattering states with transversal momentum $`\stackrel{}{q}_n`$ and longitudinal velocity $`v_n`$ at the transversal position $`\stackrel{}{r}_0`$ on the cross section of the wire (area $`A`$).

Another promising direction for future research is to study what happens to the LDOS if the wire is coupled to a superconductor rather than to a normal electron reservoir . The convenient expressions for $`\nu `$ and $`\stackrel{~}{\nu }`$ in terms of the reflection matrices from two independent parts of the wire, derived in this paper, can be directly generalized to include Andreev reflection at the interface.

We acknowledge discussions with A. D. Mirlin and M. Pustilnik. This work was supported by the Dutch Science Foundation NWO/FOM, by the NSF under grant no. DMR 0086509, and by the Sloan foundation.


format: text

text:
# Mass Determinations of Star Clusters

## 1. Introduction

Open and globular clusters are respectively located, because of their local definition, in the plane and in the halo of our Galaxy. This review will not make any formal distinction between these two kinds of star clusters since (i) the above definitions apply only to our Galaxy and (ii) we do not know if there are any genuine differences in their formation mechanisms. Consequently, we may use the words “open” and “globular”, but we shall essentially mean “star clusters”, whether they are light or massive, or young or old. This approach is also justified by the fact that, in simulations of the Galactic globular cluster system, the dynamical evolution of an initial mass distribution, of either Gaussian or power-law type, leads always to a predicted distribution consistent with observations: light star clusters dissolve rather quickly, while heavy ones survive longer (Baumgardt 1998, 2001). See also Zhang & Fall (1999) in the case of the populations of star clusters of the Antennae Galaxies.

Individual masses of star clusters are not easy to measure. As an example, we shall quote the four different mass estimates given for the giant Galactic globular cluster $`\omega `$ Centauri by Ogorodnikov et al. (1976). From a few stellar radial velocities $``$tot $``$ 7 $`\times `$ $`10^5_{}`$, from the same radial velocities corrected for the effect of the internal rotation of the cluster $``$tot $``$ 3 $`\times `$ $`10^6_{}`$, from low-quality proper motions measured on photographic plates $``$tot $``$ 3 $`\times `$ $`10^7_{}`$, and from gravitational focusing from the same proper motions $``$tot $``$ 2.8 $`\times `$ $`10^8_{}`$. These four different approaches provide results differing by 3 orders of magnitude! Fortunately, things have improved since then, but even the best mass determinations of star clusters remain rather uncertain, typically by a factor of 2. It is interesting to realize that the age of the Universe, intuitively more difficult to determine than the mass of a nearby star cluster, is more accurately known since the Hubble constant $`H_{}`$ is now estimated at a better than 20 % level.

## 2. The Giant Galactic Clusters $`\omega `$ Centauri and 47 Tucanae

Significant improvements in the quality and numbers of radial velocities have allowed the determinations of reliable mass estimates through the use of various dynamical models. On the one hand, there is the parametric approach. For example, Gunn & Griffin (1979) developed multi-mass models whose distribution functions $`f`$ depend on the stellar energy per unit mass $`\epsilon `$ and the specific angular momentum $`l`$. Such models are spherical and have a radial anisotropic velocity dispersion ($`\overline{v_r^2}`$ $``$ $`\overline{v_\theta ^2}`$ = $`\overline{v_\varphi ^2}`$). They are called King-Michie models and associate the lowered Maxwellian of the King model with the anisotropy factor of the Eddington models:

$$f(\epsilon ,l)[\mathrm{exp}(2j^2\epsilon )\mathrm{exp}(2j^2\epsilon _t)]\mathrm{exp}(j^2l^2/r_a^2)$$
$`(1)`$
On the other hand, one never knows which of the assumptions underlying such a model are adhered to by the real system and which are not (Dejonghe & Merritt 1992). These arguments suggest that it might be profitable to interpret kinematical data from globular clusters in an entirely different manner, placing much stronger demands on the data and making fewer ad hoc assumptions about the distribution function $`f`$ as well as the gravitational potential $`\mathrm{\Phi }`$. Ideally, the unknown functions should be generated non-parametrically from the data, in an approach pioneered by Merritt (1993, 1996). We provide hereafter the results of two studies, parametric and non-parametric, respectively, of the globular cluster $`\omega `$ Centauri, both studies using exactly the same observational data, viz., the surface brightness profile and 469 stellar radial velocities.

$``$ Parametric: a simultaneous fit of these radial velocities and of the surface brightness profile to a multi-mass King-Michie dynamical model provides mean estimates of the total mass for $`\omega `$ Cen equal to $``$tot = 5.1 $`\pm `$ 0.6 $`\times `$ $`10^6`$$``$, with a corresponding mean mass-to-light ratio $``$$`/L_V`$ = 4.1 (Meylan et al. 1995).

$``$ Non-parametric: the potential and mass distribution infered in this method provide a total mass for $`\omega `$ Cen equal to $``$tot = 2.9 $`\pm `$ 0.4 $`\times `$ $`10^6`$$``$, with a corresponding mean mass-to-light ratio $``$$`/L_V`$ = 2.3 (Merritt et al. 1997).

## 3. Universal Mass-to-Light Ratio?

With a similar parametric approach applied to the observational constraints obtained for NGC 1835, an old Large Magellanic Cloud globular cluster, King-Michie models give $``$tot = 1.0 $`\pm `$ 0.3 $`\times `$ $`10^6`$$``$, corresponding to a mean mass-to-light ratio $``$$`/L_V`$ = 3.4 $`\pm `$ 1.0 (Meylan 1988, Dubath & Meylan 1990).

These studies show that when the same kind of dynamical models (King-Michie) constrained by the same kind of observations (surface brightness profile and central value of the projected velocity dispersion) are applied to an old and bright Magellanic globular cluster, viz., NGC 1835, the results seem similar to those obtained in the case of Galactic globular clusters. Consequently, the rich old globular clusters in the Magellanic clouds could be quite similar in mass and $``$$`/L_V`$ to the rich globular clusters in the Galaxy.

## 4. The Final Word from Thousands of Stellar Space Motions

Recently, the amount of data related to stellar motions, viz., radial velocities and proper motions, has increased significantly. In a pioneering ground-based study of $`\omega `$ Centauri, van Leeuwen et al. (2000) measured the individual proper motions of 7853 probable member stars, from photographic plates with epochs ranging from 1931 through 1935 and 1978 through 1983. An internal proper motion dispersion of 1.0 to 1.2 mas yr-1, equivalent to 25 to 29 km s-1 for a distance of 5.1 kpc, is found for members near the cluster center. This dispersion decreases to 0.3 mas yr-1, equivalent to 7.5 km s-1 in the outer regions. The full dynamical interpretation of these proper motions, combined with a few thousand radial velocities, is in preparation by the same group.

Another group, using HST/WFPC2 images, has obtained slightly better measurements of proper motions for about 15,000 stars in the core of 47 Tucanae, this within a time baseline of only 4 years, with 3 epochs separated by 2 years. See Anderson, King, & Meylan 1998 for a progress report after the second epoch. These data, combined with the radial velocities of about 5,000 stars, will provide, as in the case of $`\omega `$ Centauri, an insight into the dynamics of the core of 47 Tucanae, with fundamental by-products such as cluster distance and photometry of variable stars and binaries.

With their new proper-motion techniques and software (Anderson & King 2000) applied to the HST/WFPC2 archive data for the first epoch and to their own observations for the second epoch, the members of the same group will soon have available similar sets of proper motions for about ten among of the richest and most nearby Galactic globular clusters.

The above studies, with thousands of proper motions and radial velocities constraining dynamical models with three integrals of the motion as well as non-parametric ones, will allow a significant step forward in our understanding of the internal dynamics of massive star clusters.

## 5. Mayall II $``$ G1, a Giant Globular Cluster in M31

Mayall II $``$ G1 is one of the brightest globular clusters belonging to M31, the Andromeda galaxy. Observations with HST/WFPC2 provide photometric data for the $`I`$ vs. $`VI`$ and $`V`$ vs. $`VI`$ color-magnitude diagrams. They reach stars with magnitudes fainter than $`V`$ = 27 mag, with a well populated red horizontal branch at about $`V`$ = 25.3 mag (Meylan et al. 2001). From model fitting, that study determines a rather high mean metallicity of \[Fe/H\] = – 0.95 $`\pm `$ 0.09, somewhat similar to 47 Tucanae. In order to determine the true measurement errors, Meylan et al. (2001) have carried out artificial star experiments. They find a larger spread in $`VI`$ than can be explained by the measurement errors. They attribute this to an intrinsic metallicity dispersion among the stars of G1, which may be the consequence of self-enrichment during the early stellar/dynamical evolutionary phases of this cluster. So far, only $`\omega `$ Centauri, the giant Galactic globular cluster, has been known to exhibit such an intrinsic metallicity dispersion. This is, a phenomenon certainly related to the deep potential well of each of these two star clusters, which are massive enough to retain the gas expelled by the first generations of very massive stars.

The structural parameters of G1 are deduced from the same HST/WFPC2 data. Its surface brightness profile provides its core radius $`r_c`$ = 0.14′′ = 0.52 pc, its tidal radius $`r_t`$ $``$ 54′′ = 200 pc, and its concentration $`c`$ = log ($`r_t/r_c`$$``$ 2.5. Such a high concentration indicates the probable collapse of the core of G1. KECK/HIRES observations provide the central velocity dispersion $`\sigma _{obs}`$ = 25.1 km s-1, with $`\sigma _p`$(0) = 27.8 km s-1 once aperture corrected.

Three estimates of the total mass of this globular cluster can be obtained. The King-model mass is $``$K = 15 $`\times `$ $`10^6`$$``$ with $``$$`/L_V`$ $``$ 7.5, and the Virial mass is $``$Vir = 7.3 $`\times `$ $`10^6`$$``$ with $``$$`/L_V`$ $``$ 3.6. The King-Michie model fitted simultaneously to the surface brightness profile and the central velocity dispersion value provides mass estimates ranging from $``$KM = 14 $`\times `$ $`10^6`$$``$ to 17 $`\times `$ $`10^6`$$``$ (Meylan et al. 2001).

The spread between the three mass determination values listed in Table 1 give a better idea than their individual formal (smaller) errors about their true uncertainties. The masses of both clusters are known to about a factor of two. Although not very precise, all of these mass estimates make G1 more than twice as massive as $`\omega `$ Centauri, the most massive Galactic globular cluster. G1 is unique in M31 because of its projected location 40 kpc away from the center of the galaxy, but there are at least three other bright globular clusters in this galaxy which have velocity dispersions $`\sigma _{obs}`$ larger than 20 km s-1, implying rather large masses.

## 6. On the Origin of the Most Massive Globular Clusters

Such large masses are related to the metallicity spread the origin of which is still unknown. It may come either (i) from metallicity self-enrichment in a massive globular cluster, (ii) from primordial metallicity inhomogeneity in a binary proto-cluster cloud, followed by early merger, or (iii) from the fact that the present globular cluster is merely the remaining core of a previously larger entity, e.g., originally a dwarf galaxy subsequently pruned by dynamical evolution.

Although $`\omega `$ Centauri is the best studied globular cluster, because of its size and relative proximity ($``$ 5.1 kpc), many conundrums remain: viz., (i) the metallicity spread among stars (Freeman & Norris 1981), (ii) its double Main Sequence (Anderson 1997), (iii) the different kinematics between metal-rich and -poor stars (Norris et al. 1997), and (iv) a correlation between metallicity and age, implying that this cluster enriched itself over a timescale $``$ 3 Gyr (Hughes & Wallerstein 2000 and Hilker & Richtler 2000).

Let us consider (see Table 2 and Meylan et al. 2001) the following four parameters relative to G1: the central velocity dispersion $`\sigma _p`$(0) = 28 km s-1, the integrated absolute visual magnitude $`M_V`$ = – 10.94 mag, the core radius $`r_c`$ = 0.52 pc, and the central surface brightness $`\mu `$(0,$`V`$) = 13.47 mag arcsec-2. The positions of G1 in the different diagrams defined by Kormendy (1985), using the above four parameters, always put it on the sequence defined by globular clusters, and definitely away from the other sequences defined by elliptical galaxies, bulges, and dwarf spheroidal galaxies (Fig. 1). The same is true for $`\omega `$ Centauri (Meylan et al. 2001).

Little is known about the positions, in these diagrams, of the nuclei of nucleated dwarf elliptical galaxies, which could be the progenitors of the most massive, if not all, globular clusters (Zinnecker 1988, Freeman 1993). The above four parameters are known only for the nucleus of one dwarf elliptical, viz., NGC 205, and their values put this object, in Kormendy’s diagram, close to G1, right on the sequence of globular clusters (see Table 2 and Fig. 1). This result does not prove by itself that all massive globular clusters are the remnant cores of nucleated dwarf galaxies.

At the moment, only the anti-correlation of metallicity with age recently observed in $`\omega `$ Centauri suggests that this cluster enriched itself over a time scale of about 3 Gyr (Hughes & Wallerstein 2000 and Hilker & Richtler 2000). This contradicts the general idea that all the stars in a globular cluster are coeval, and may favor the origin of $`\omega `$ Centauri as being the remaining core of a larger entity, e.g., of a former nucleated dwarf elliptical galaxy. In any case, by the mere fact that their large masses imply complicated stellar and dynamical evolutions, the very massive globular clusters may blur the former clear (or simplistic) difference between globular clusters and dwarf galaxies.

## 7. Mass Estimates of Young Star Clusters

So far, we have discussed only mass estimates related to old, rich star clusters. HST has triggered numerous studies of young, bright starburst clusters which may be quite massive. See, e.g., Holtzmann et al. (1992) in the case of NGC 1275 clusters, and Schweizer & Seitzer (1993) and Whitmore et al. (1993) in the case NGC 7252. Conti & Vacca (1994) use HST/FOC UV imaging of various bright knots in the nuclear starburst region of the Wolf-Rayet galaxy H2 2-10. The luminosities of the knots are compared to those predicted from stellar population synthesis models for ages between 1 and 10 Myr, from which they estimate the masses of the knots to be between 105 and 106 $``$. In a more direct way, Böker et al. (1999) measure a velocity dispersion $`\sigma _{obs}`$ = 33 $`\pm `$ 3 km s-1 in the nuclear star cluster of the face-on giant Scd spiral galaxy IC 342 and deduce, by fitting the central surface brightness profile, a cluster mass $``$tot = 6 $`\times `$ $`10^6`$$``$. They infer a best-fitting cluster age in the range 10 - 60 $`\times `$ $`10^6`$yr.

It is worth emphasizing that mass estimates of young star clusters are more difficult to obtain and more uncertain than those of older clusters. This is due to the frequent lack of velocity dispersion measurement and because of the large uncertainties inherent in the use of stellar population synthesis models for young stellar systems whose luminosities evolve quickly (see, e.g., Bruzual 2001).

## 8. The Slow Destruction of Star Clusters

Numerical works have demonstrated that continual two-body relaxation within globular clusters combined with weak tidal encounters between globular clusters and the Galactic disk and/or bulge will lead, around each globular cluster, to the development of both a halo of unbound stars and tidal tails. All clusters observed, which do not suffer from strong observational biases, present tidal tails (Grillmair et al. 1995, Leon et al. 2000). These tidal tails exhibit projected directions preferentially aligned with the cluster orbit or towards the Galactic center, betraying their recent dynamical evolution through disk and/or bulge shocking. See also Odenkirchen et al. (2001).

Recent theoretical work corroborate these observations. In N-body simulations of globular clusters moving in the Galactic potential well, Combes et al. (1999) observe that once the particles (stars) are unbound, they slowly drift along the globular cluster path and form two huge tidal tails. A cluster is always surrounded by two giant tidal tails and debris, in permanence along its orbit. The length of these tidal tails is of the order 5 tidal radii or more. The orientation of these tidal tails is the signature of the last disk crossing and can constrain strongly the cluster orbit and the Galactic model. Each disk/bulge crossing may extract up to about 1% of the total mass of the star cluster, leading slowly but surely to its complete evaporation. The lighter the cluster and the stronger the tidal shocks, the faster the destruction process.

## 9. Do Globular Clusters Have Dark Halo?

Are globular clusters the most massive stellar systems without non-baryonic dark matter? Since both the velocity dispersion profile and the rotation curve of $`\omega `$ Centauri decrease towards the outer parts of this cluster, we may conclude that there is no dynamical evidence of any massive halo made of non-baryonic dark matter.

In order to investigate further this point among the most poorly studied Galactic clusters, Côté et al. (2001) have obtained KECK/HIRES high-accuracy radial velocities of about 20 stars in six outer-halo globular clusters, located between 20 and 100 kpc from the center of our Galaxy. The velocity dispersions range between $`\sigma _{obs}`$ = 1 and 5 km s-1, with corresponding $``$$`/L_V`$ values between 1 and 4. This is exactly what is expected from the nearby well studied globular clusters such as $`\omega `$ Centauri and 47 Tucanae.

However there is one single, although conspicuous, exception in the sample of Côté et al. (2001): Pal 13 exhibits a velocity dispersion $`\sigma _{obs}`$ = 2.6 $`\pm `$ 0.3 km s-1. With its low and uncertain total luminosity, Pal 13 has a corresponding $``$$`/L_V`$ in the range 10 $`<`$ $``$$`/L_V`$ $`<`$ 40. This is quite unique for a globular cluster. Simulations show that such a result would be mimicked, with difficulties, only by an uncomfortably large fraction of binary stars. Binaries exist in globular clusters, but not in large quantities. We should apply Ockham’s razor before invoking the presence of a massive dark halo. Such a high $``$$`/L_V`$ value could be an indication of a cluster in the late phases of its dynamical evolution, when a large fraction of its total mass is made of white dwarfs, as predicted by Vesperini & Heggie (1997). It could also be explained by a velocity dispersion inflated because of lack of Virial equilibrium: Pal 13 seems to be in the advanced stages of tidal disruption, a status supported by the recent observations of extra-tidal extension around this cluster (Siegel et al. 2001 and Côté et al. 2001). These two explanations do not exclude each other.

## 10. Is Tidal Disruption an Ubiquitous Phenomenon?

Leon et al. (2000) show that any Galactic globular cluster observed which does not suffer from strong observational biases (e.g., intervening absorption by Galactic cirrus along the line-of-sight) displays a pair of tidal tails. The tidal disturbance may also be unveiled through its effect on the surface brightness profile. An isolated cluster suffering no tidal shock will have a normal King-model profile characterised by a core surrounded by an envelope with a steep profile. Any tidally-perturbed cluster will have an envelope profile departing from the King-model profile as if there were a very strong background of stars: the stronger the tidal effect, the higher the level of the background (i.e., the higher the departure from King-model profile), which is made of stars escaping the cluster gravitational attraction.

When observing a tidally perturbed globular cluster, there is a high probability of measuring radial velocities of escaping stars, which are no more in Virial equilibrium with the cluster. Such stars would inflate the measured velocity dispersion. The stronger the tidal disturbance, the higher the inflated velocity dispersion, and the larger the mass-to-light ratio. This is exactly what is observed in Fig. 2, from Majewksi et al. (2001), which displays the surface brightness profiles of two globular clusters (NGC 288 with $``$$`/L_V`$ = 2.1 and Pal 13 $``$$`/L_V`$ $``$ 20) and five dSph galaxies (Sculptor with $``$$`/L_V`$ = 3.0, Leo I with $``$$`/L_V`$ = 5.6, Leo II with $``$$`/L_V`$ = 17, Carina with $``$$`/L_V`$ = 31, and Ursa Major with $``$$`/L_V`$ = 79). From the bottom panel (NGC 288) to the top one (Ursa Major) we observe a sequence of increasingly disturbed profiles, with higher and higher departures from the King models, with the most disturbed profile being for Ursa Major. Interestingly, this sequence of increasingly disturbed profiles correspond to a sequence of increasing mass-to-light ratios, as if the $``$$`/L_V`$ values were a direct measure of the intensity of the tidal shocks. Is the evidence for dark matter in dSph galaxies evaporating with their stars?

### Acknowledgments.

It is a pleasure to thank Steve Majewksi (Virginia) and my collaborators Patrick Côté (Rutgers) and George Djorgovski (Caltech) for allowing me to present results in advance of publication. I also thanks Jennifer Lotz (JHU), Roeland van der Marel (STScI), and Brad Whitmore (STScI) for interesting discussions and information about NGC 205.

## References

Anderson J., 1997, Ph.D. thesis, University of California, Berkeley

Anderson J., King I.R., Meylan G., 1998, BAAS, 30, 1347

Anderson J., King I.R., 2000, PASP, 112, 1360

Baumgardt H., 1998, A&A, 330, 480

Baumgardt H., 2001, in ASP Conf. Ser. Vol. ???, Modes of Star Formation and the Origin of Field Populations, ed. E.K. Grebel and W. Brandner (San Francisco: ASP), in press

Böker T., van der Marel R.P., Vacca W.D., 1999, AJ, 118, 831

Bruzual G.A., 2001, in XI Canary Islands Winter School of Astrophysics, Galaxies at High Redshift, ed. I. Pérez-Fournon and F. Sánchez, (Cambridge: Cambridge Contemporary Astrophysics), in press

Combes F., Leon S., Meylan G., 1999, A&A, 352, 149

Conti P.S., Vacca W.D., 1994, ApJ, 423, L97

Côté P., Djorgovski S.G., Meylan G., Castro S., McCarthy J.K., 2001, AJ, submitted

Dejonghe H., Merritt D., 1992, ApJ, 391, 531

Dubath P., Meylan G., Mayor M., Magain P., 1990, A&A, 239, 142

Freeman K.C., 1993, in IAU Symp. 153, Galactic Bulges, ed. H. Dejonghe & H.J. Habing (Dordrecht: Kluwer), 263

Freeman K.C., Norris J.E., 1981, ARA&A, 19, 319

Grillmair C. J., Freeman K. C., Irwin M., 1995, AJ, 109, 2553

Gunn J.E., Griffin R.F., 1979, AJ, 84, 752

Hilker M., Richtler T., 2000, A&A, 362, 895

Holtzmann J.A., et al. 1992, AJ, 103, 691

Hughes J., Wallerstein G., 2000, AJ, 119, 1225

Kormendy J., 1985, ApJ, 295, 73

Leon S., Meylan G., Combes F., 2000, A&A, 359, 907

Majewksi S.R., et al., 2001, in preparation

Meylan G., 1988, ApJ, 331, 718

Meylan G., Heggie D.C., 1997, A&AR, 8, 1-143

Meylan G., Mayor M., Duquennoy A., Dubath P., 1995, A&A, 303, 761

Meylan G., Sarajedini A., Jablonka P., et al., 2001, AJ, in press

Merritt D., 1993, ApJ, 413, 79

Merritt D., 1996, AJ, 112, 1085

Merritt D., Meylan G., Mayor M., 1997, AJ, 114, 1074

Norris J.E., Freeman K.C., Mayor M., Seitzer P., 1997, ApJ, 487, L187

Odenkirchen M., Grebel E.K., et al., 2001, ApJ, 548, L165

Ogorodnikov K.F., Nezhinskii E.M., Osipkov L.P., 1976, Sov. Astr. Lett. 2, 57

Schweizer F., Seitzer P., 1993, ApJ, 417, L29

Siegel M.H., Majewski S.R., Cudworth K.M., Takamiya M., 2001, AJ, 121, 935

van Leeuwen F., Le Poole R.S., Reijns R.A., Freeman K.C., de Zeeuw P.T., 2000, A&A, 360, 472

Vesperini E., Heggie D.C., 1997, MNRAS, 289, 898

Whitmore B.C., Schweizer F., Leitherer C., Borne K., Robert C., 1993, AJ, 106, 1354

Zhang Q., Fall S.M., 1999, ApJ, 527, L81

Zinnecker H., Keable C.J., Dunlop J.S., Cannon R.D., Griffiths W.K., 1988, in IAU Symp. 126, Globular Cluster Systems in Galaxies, ed. J.E. Grindlay, A.G.D. Philip, (Dordrecht: Kluwer), 603

## Discussion

Armandroff:  You raised the possibility of tidal disruption leading to the large velocity dispersion in Ursa Minor. Models for the disruption of dSphs (Piatek & Pryor, 1995, AJ, 109, 1071; Oh et al., 1995, ApJ, 442, 142) predict a kinematic signature that resembles a rotation curve, as opposed to simply an increased dispersion. Large samples of radial velocities in Draco and Ursa Minor have not shown evidence for the kinematic signature of disruption. Do you have evidence for Ursa Minor or other dSphs matching the modeling?

Meylan:  The two papers you mention present results of N-body simulations which reach conclusions different from those obtained from other, more recent N-body simulations (Kroupa, 1997, New Astron, 2, 139; Klessen & Kroupa, 1998, ApJ, 498, 143). Kroupa’s work shows that dSphs are obtained under the extreme hypothesis that dSph progenitors are not dominated by dark matter but are significantly shaped by tides through many periastron passages; Piatek & Pryor simulated only one passage, and Kroupa’s results are completely consistent with their findings for one passage only. There are now in Draco photometric observations of stars beyond its measured tidal boundary (Piatek et al., 2001 AJ, 121, 841). It will be essential to obtain, for a few of these dSphs, samples of thousands of stellar radial velocities. At the moment, I simply find the correlation between tidal disturbance intensity and $``$$`/L`$ values rather intriguing. This definitely calls for more studies.

Da Costa:  Comment: Integrated spectrum of nucleus of NGC 205 is A-type, so its stellar population is very different from globular clusters like G1. Question: What is the core radius of the NGC 205 nucleus and does if fit a King model? (different from say M32 with its power-law cusp).

Meylan:  About the comment: Yes, the stellar populations in the core of NGC 205 are younger than those in G1. This means that the position, in the panels defined by Kormendy (see Fig. 1 above), of a dynamical system older than about 1 Gyr does not depend strongly on its age. Answer to the question: The core radius of the nucleus NGC 205 is 100 mas = 0.35 pc (Jones et al., 1996, ApJ, 466, 742), similar to a very dense globular cluster easily fitted by a King model. The star cluster in the center of NGC 205 is very different from the nucleus of M32, which has, e.g., a much larger velocity dispersion $`\sigma _{obs}`$ = 150 km s-1(Joseph et al., 2001, ApJ, 550, 668). See also Meylan et al. (2001).

Grillmair:  Two comments: 1) The onset of tidal effects (e.g., break from normal profile to power-law tidal tail, velocity dispersion increase) can depend significantly on the orbital phase of the cluster (strongest effect at apogalacticon). 2) Ben Moore (1996, ApJ, 461, L13) put an upper limit of $``$$`/L`$ of NGC 7089 $``$ M2, based on existence of tidal tails.

Wenderoth:  Around 15 years ago, the idea that $`\omega `$ Centauri is a merger was proposed. Since then, is there more evidence to support this idea or to reject it?

Meylan:  There was a paper by Icke & Alcaino, 1988, A&A, 204, 115, suggesting that $`\omega `$ Centauri could be the result of a merger, this in order to explain (i) its spread in metallicity, (ii) its strong flattening and (iii) its large mass. There is now more evidence about the intrinsic complexity of this globular cluster, as mentioned above in Section 6, although none of these observational facts supports the merger scenario over the other two alternatives. If merger there was, it must have been very early in the life of the two merging proto-cluster clouds, since the mere encounter of two current globular clusters would not induce a merger, except in the case of orbital-parameter adjustment with vanishingly small probabilities.


text:
# Time Reversal for Classical Waves in Random Media

## Version française abrégée

Dans cette note, nous nous intéressons à la propagation d’ondes classiques en milieu aléatoire et à leurs propriétés de refocalisation après renversement temporel. Dans les expériences conduites par M. Fink  et ses collaborateurs, une source localisée émet des ondes acoustiques à travers un milieu hétérogène. Le signal est enregistré en temps par des capteurs localisés en espace et positionnés loin de la source, puis renversé en temps et enfin retransmis dans le milieu (ce qui est arrivé en dernier repart en premier).

Dans ces expériences, il a été observé le phénomène suivant, qui à première vue pourrait paraître contre-intuitif : le signal réémis refocalise à l’emplacement de la source initiale d’autant mieux que le milieu est hétérogène, alors que la qualité de la refocalisation est très mauvaise en milieu homogène. L’explication réside dans la multiplicité des “chemins” menant de la source aux capteurs. Plus le milieu est aléatoire, plus il donne d’informations sur la source aux capteurs, pourvu que les enregistrements soient suffisamment longs.

Plusieurs travaux sur ce phénomène existent dans la littérature physique . Les résultats mathématiques existants portent sur le cas mono-dimensionnel  et le régime de faisceaux étroits . Dans cette note, nous proposons une théorie qui explique le phénomène de refocalisation pour des ondes classiques à hautes fréquences en milieu multi-dimensionnel. Notre approche, s’inspirant de , se fonde sur la transformée de Wigner, définie en (4), qui permet l’analyse des corrélations de noyaux de Green en des points voisins. Nous pouvons écrire le signal retransmis après renversement du temps en fonction de cette transformée de Wigner (5).

Selon les caractéristiques du milieu aléatoire, la transformée de Wigner satisfait différentes équations limites quand la fréquence caractéristique des ondes tend vers l’infini. Nous analysons les deux cas où nous obtenons à la limite soit une équation de Liouville (limite semiclassique), soit une équation de transport de type Boltzmann (limite de couplage faible). Dans le premier cas, nous sommes capables d’obtenir une dérivation rigoureuse mathématiquement. La formule de refocalisation est donnée en (8)-(12). Dans le second cas, la dérivation de la formule de reconstruction (8)-(14) est formelle en certains endroits. La théorie permettant d’obtenir l’équation de transport limite n’est complète que pour l’équation (scalaire) de Schrödinger . De plus elle nécessite de moyenner sur les réalisations du milieu aléatoire, bien que les expériences physiques montrent que ce soit superflu.

## 1. Refocusing in time-reversal

Time-reversal experiments  can be briefly described as follows. A point source signal is sent through an inhomogeneous medium and is recorded by a spatially localized array of receivers-transducers. The signal is subsequently reversed in time and reemitted back into the medium; that is, the part of the signal that was recorded last is sent back first. The refocusing of the reemitted signal is then observed on the spot of the original source.

Refocusing in a homogeneous medium is poor when only a few receivers are used to record the signal. The most striking and somewhat counterintuitive observation is that inhomogeneities in the medium enhance the refocusing effects. This is because waves recorded at the time-reversal array have sampled a larger part of the medium than in the homogeneous case and carry more information about the source location. Furthermore, the refocused signal in a random medium is self-averaging, that is, independent of the realization of the random medium.

Many time-reversal experiments have been conducted in the recent past and studied theoretically . Mathematical studies have concentrated so far on the one space dimension case  and on the narrow beam approximation . In this note we present an analysis of the enhanced refocusing property in genuinely three dimensional heterogeneous media in the regimes of random geometrical acoustics and radiative transport. Our main tool is the asymptotic analysis of the equation for the Wigner transform . It turns out that in the high frequency limit the refocused signal may be written as a convolution with a kernel $`F`$ expressible in terms of the limit Wigner distribution. More precisely, this kernel is expressed (9) in terms of the Fourier transform in $`𝐤`$ of the energy densities $`a(t,𝐱,𝐤)`$ in the phase space (6). The latter are smooth in $`𝐤`$ under evolution in random media but not in a homogeneous domain. Therefore the spatial localization of the kernel $`F`$ is improved in a random medium compared to the homogeneous case. In the geometric acoustics limit, the Wigner transform asymptotically satisfies a Liouville or a Fokker-Plank equation similar to the one used to describe the refocusing of narrow beam signals . In the radiative transfer regime, the Wigner distribution satisfies a linear Boltzmann equation. The rigorous passage to this limit is not complete for classical waves  and our results in this part remain formal.

## 2. Reemitted signal and the Wigner transform

The propagation of many classical waves is described by systems of $`m`$ first-order differential equations of the form (repeated indices are summed over throughout the paper)

$$A(𝐱)\frac{𝐮(t,𝐱)}{t}+D^j\frac{𝐮(t,𝐱)}{x_j}=0,𝐱^d,d=2\text{ or }3$$
(1)
with initial conditions $`𝐮(0,𝐱)=𝐮_0(𝐱)`$. In this paper, we address time-reversal for acoustic waves with $`d=3`$ to fix notation. The corresponding theory for other types of waves will be presented elsewhere . In this context, $`𝐮`$ is the 4-vector composed of the three components of the velocity field $`𝐯`$ and the scalar pressure field $`p`$. The matrix $`A(𝐱)=\mathrm{Diag}(\rho (𝐱),\rho (𝐱),\rho (𝐱),\kappa (𝐱))`$, where $`\rho (𝐱)`$ is the density of the underlying medium and $`\kappa (𝐱)`$ its compressibility. The $`4\times 4`$ matrices $`D^j`$ have entries $`D_{ml}^j=\delta _{m4}\delta _{lj}+\delta _{l4}\delta _{mj}`$. The solution of (1) is $`𝐮(t,𝐱)=_^dG(t,𝐱;𝐲)𝐮_0(𝐲)𝑑𝐲.`$ Here $`G(t,𝐱;𝐲)`$ is the matrix-valued Green’s propagator from $`𝐲`$ to $`𝐱`$. It solves (1) with the initial condition $`G(0,𝐱;𝐲)=I\delta (𝐱𝐲)`$ where $`I`$ is the $`4\times 4`$ identity matrix. Our time reversal setting is as follows. At time $`t`$ the wave field is truncated in a region $`\mathrm{\Omega }`$, where the transducers-receivers are located. Then the direction of the acoustic velocity field is reversed and the signal is re-emitted into the medium. The back-propagated signal emanating from the source term $`𝐮_0`$ at time $`0`$ can be written using the Green’s propagator as

$$𝐮^B(𝐱)=_{^{3d}}\mathrm{\Gamma }G(t,𝐱;𝐲)\mathrm{\Gamma }G(t,𝐲^{};𝐳)\chi _\mathrm{\Omega }(𝐲)\chi _\mathrm{\Omega }(𝐲^{})f(𝐲𝐲^{})𝐮_0(𝐳)𝑑𝐲𝑑𝐲^{}𝑑𝐳.$$
(2)
Here $`\mathrm{\Gamma }`$ is a matrix that models the time reversal process. It is given by $`\mathrm{\Gamma }=\mathrm{Diag}(1,1,1,1)`$, so that the velocity field $`𝐯`$ is replaced by $`𝐯`$ and the pressure field $`p`$ remains unchanged. The function $`\chi _\mathrm{\Omega }(𝐲)`$ may equal $`1`$ on $`\mathrm{\Omega }`$ and $`0`$ outside, or could be a more general function modeling amplification of the time-reversed signal by various receivers. We also allow for some blurring of the recorded signal before it is transmitted back. The convolution with a filter $`f`$ describes this blurring. Introducing the adjoint Green’s matrix $`G_{}`$, solution of

$$\begin{array}{c}\frac{G_{}(t,𝐱;𝐲)}{t}+\frac{}{x_j}(G_{}(t,𝐱;𝐲))D^jA^1(𝐱)=0,\hfill \end{array}$$
(3)
with initial condition $`G_{}(0,𝐱;𝐲)=\delta (𝐱𝐲)A^1(𝐱)`$, we observe that $`\mathrm{\Gamma }G_{}(t,𝐱;𝐲)A(𝐱)\mathrm{\Gamma }=G(t,𝐲;𝐱)`$.

We now rescale our problem. Refocusing is expected to be important in the close vicinity of the source location. We take the initial source as a localized function $`𝐮_0(𝐱)=𝐒\left(\frac{𝐱𝐱_0}{\epsilon }\right)`$ of finite amplitude and rescale the filter accordingly: $`\frac{1}{\epsilon ^d}f(\frac{𝐲𝐲^{}}{\epsilon })`$. Here $`\epsilon 1`$ is a small parameter that measures the ratio of the width of the initial pulse to the propagation distance $`L`$ between the source location and the array of receivers-transducers. An observation point $`𝐱`$ is close to $`𝐱_0`$ and we write it as $`𝐱=𝐱_0+\epsilon 𝝃`$. After a change of variables, equation (2) becomes

$$𝐮_\epsilon ^B(𝝃;𝐱_0)=\underset{^{3d}}{}\mathrm{\Gamma }G(t,𝐱_0+\epsilon 𝝃;𝐲)G_{}(t,𝐱_0+\epsilon 𝐳;𝐲^{})A(𝐱_0+\epsilon 𝐳)\mathrm{\Gamma }𝐒(𝐳)\chi _\mathrm{\Omega }(𝐲)\chi _\mathrm{\Omega }(𝐲^{})f(\frac{𝐲𝐲^{}}{\epsilon })𝑑𝐲𝑑𝐲^{}𝑑𝐳.$$

The Wigner transform is a natural tool in the analysis of the correlation function of wave fields at neighboring points . In our context, we define it as

$$W_\epsilon (t,𝐱,𝐤)=_{^{2d}}[_^de^{i𝐤𝐳}G(t,𝐱\frac{\epsilon 𝐳}{2};𝐲)G_{}(t,𝐱+\frac{\epsilon 𝐳}{2};𝐲^{})\frac{d𝐳}{(2\pi )^d}]\chi _\mathrm{\Omega }(𝐲)\chi _\mathrm{\Omega }(𝐲^{})f(\frac{𝐲𝐲^{}}{\epsilon })𝑑𝐲𝑑𝐲^{}.$$
(4)
This allows us to recast the expression for the retransmitted signal as follows

$$𝐮_\epsilon ^B(𝝃;𝐱_0)=_{^{2d}}\mathrm{\Gamma }W_\epsilon (t,𝐱_0+\epsilon \frac{𝝃+𝐳}{2},𝐤)e^{i𝐤(𝐳𝝃)}A(𝐱_0+\epsilon 𝐳)\mathrm{\Gamma }𝐒(𝐳)𝑑𝐳𝑑𝐤.$$
(5)
Let us define the space $`𝒜(^{2d})`$ as the subset of $`𝒮^{}(^{2d})`$ of matrix-valued distributions $`\eta (𝐱,𝐤)`$ such that $`_^dsup_𝐱\widehat{\eta }(𝐱,𝐲)d𝐲`$ is bounded. Here $`\widehat{\eta }(𝐱,𝐲)`$ is the inverse Fourier transform of $`\eta (𝐱,𝐤)`$ in the second variable only. We denote by $`𝒜^{}`$ the dual space to $`𝒜`$. Then we have
Lemma 1. The Wigner transform $`W_\epsilon (t,𝐱,𝐤)`$ is bounded in $`𝒞^0(0,T;𝒜^{}(^{2d}))`$ independent of $`\epsilon `$ provided that $`\widehat{f}(𝐤)L^1(^d)`$. As a consequence, it converges weakly along a subsequence $`\epsilon _k0`$ to a distribution $`W(t,𝐱,𝐤)𝒞^0(0,T;𝒜^{}(^{2d}))`$.

The proof of this lemma is obtained by rewriting $`W_\epsilon `$ as

$$W_\epsilon (t,𝐱,𝐤)=_^d\widehat{f}(𝐪)_^de^{i𝐤𝐳}\stackrel{~}{G}_\epsilon (t,𝐱\frac{\epsilon 𝐳}{2};𝐪)\stackrel{~}{G}_\epsilon (t,𝐱+\frac{\epsilon 𝐳}{2};𝐪)\frac{d𝐳}{(2\pi )^d}𝑑𝐪,$$

where $`\stackrel{~}{G}_\epsilon (t,𝐱;𝐪)=_^dG(t,𝐱;𝐲)\chi _\mathrm{\Omega }(𝐲)e^{i𝐪𝐱/\epsilon }𝑑𝐲`$ solves (1) with initial data $`\chi _\mathrm{\Omega }(𝐱)e^{i𝐪𝐱/\epsilon }I`$, and $`\stackrel{~}{G}_\epsilon (t,𝐱;𝐪)=_^dG_{}(t,𝐱;𝐲)\chi _\mathrm{\Omega }(𝐲)e^{i𝐪𝐱/\epsilon }𝑑𝐲`$ solves (3) with initial data $`\chi _\mathrm{\Omega }(𝐱)e^{i𝐪𝐱/\epsilon }A^1(𝐱)`$. The functions $`\stackrel{~}{G}_\epsilon `$ and $`\stackrel{~}{G}_\epsilon `$ are uniformly (in $`\epsilon `$ and $`𝐪`$) bounded in $`L^2`$, and hence  $`_^de^{i𝐤𝐳}\stackrel{~}{G}(t,𝐱\frac{\epsilon 𝐳}{2};𝐪)\stackrel{~}{G}_{}(t,𝐱+\frac{\epsilon 𝐳}{2};𝐪)\frac{d𝐳}{(2\pi )^d}`$ is uniformly bounded in $`𝒞^0(0,T;𝒜^{}(^{2d}))`$. This implies the result of the Lemma.
The Wigner distribution at time $`t=0`$ is given by $`W(0,𝐱,𝐤)=|\chi _\mathrm{\Omega }(𝐱)|^2\widehat{f}(𝐤)A^1(𝐱)`$. The dispersion matrix $`L(𝐱,𝐤)=A^1(𝐱)k_jD^j`$ has a double eigenvalue $`\omega _0=0`$ that corresponds to vortical modes, and simple eigenvalues $`\omega _{1,2}=\pm c(𝐱)|𝐤|`$, $`c(𝐱)=1/\sqrt{\rho (𝐱)\kappa (𝐱)}`$. As was shown in  the limit Wigner distribution may be decomposed as

$$W(t,𝐱,𝐤)=\underset{j=1}{\overset{2}{}}a_{ij}^0(t,𝐱,𝐤)𝐛_i^0𝐛_j^0+a_1(t,𝐱,𝐤)𝐛^1𝐛^1+a_2(t,𝐱,𝐤)𝐛^2𝐛^2.$$
(6)
Here $`𝐛`$ are eigenvectors of the matrix $`L`$. If the source $`𝐒(𝐳)`$ has no vortical waves, that is, if $`\widehat{𝐒}(𝐤)`$ is parallel to $`𝐤`$, then the limit of the back-propagated signal is

$`𝐮^B(𝝃;𝐱_0)`$ $`=`$ $`{\displaystyle \underset{m=1}{\overset{2}{}}}{\displaystyle _{^{2d}}}a_m(t,𝐱_0,𝐤)\mathrm{\Gamma }𝐛^m(𝐱_0,𝐤)𝐛^m(𝐱_0,𝐤)e^{i𝐤(𝐳𝝃)}A(𝐱_0)\mathrm{\Gamma }𝐒(𝐳)𝑑𝐳𝑑𝐤`$ (7)
$`=`$ $`{\displaystyle _^d}F(t,𝝃𝐳;𝐱_0)𝐒(𝐳)𝑑𝐳=(F(t,;𝐱_0)𝐒)(𝝃),`$ (8)
with $`a_1(0,𝐱,𝐤)=a_2(0,𝐱,𝐤)=|\chi _\mathrm{\Omega }(𝐱)|^2\widehat{f}(𝐤)`$. The quality of the refocusing of the back-propagated signal is determined by the decay properties in $`𝝃`$ of the kernel

$$F(t,𝝃;𝐱_0)=\underset{m=1}{\overset{2}{}}_^da_m(t,𝐱_0,𝐤)\mathrm{\Gamma }𝐛^m(𝐱_0,𝐤)𝐛^m(𝐱_0,𝐤)e^{i𝐤𝝃}A(𝐱_0)\mathrm{\Gamma }𝑑𝐤.$$
(9)
In the limit where $`f=\delta `$ and $`\mathrm{\Omega }=^d`$, we verify that $`𝐮^B(𝝃;𝐱_0)=𝐒(𝝃)`$. When all the information is propagated back, the refocusing is perfect, as is physically expected. However, when the medium is homogeneous and $`c(𝐱)=c_0`$, the amplitudes $`a_{1,2}(t,𝐱_0,𝐤)=|\chi _\mathrm{\Omega }(𝐱_0c_0\widehat{𝐤}t)|^2\widehat{f}(𝐤)`$ become more and more singular in $`𝐤`$ as time grows (their gradient in $`𝐤`$ grows linearly in time). The corresponding kernel $`F=F_H`$ given by (9) being a Fourier transform in $`𝐤`$ of a function with growing gradients decays more and more slowly (in $`𝝃`$) as time grows. This leads to poorer and poorer refocusing.

## 3. Refocusing via the Liouville and Fokker-Plank equations

In this section, we analyze the random geometric acoustics limit, where the correlation length $`\delta `$ of the random medium is small compared to the propagation distance $`L`$ but large relative to the wavelength $`\epsilon `$. The mode $`a=a^1`$ satisfies the Liouville equation

$$_ta^\delta +_𝐤\omega _\delta _𝐱a^\delta _𝐱\omega _\delta _𝐤a^\delta =0,$$
(10)
where $`\omega (𝐱,𝐤)=(c_0+\sqrt{\delta }c_1(𝐱/\delta ))|𝐤|`$. The Liouville equation is obtained thanks to the a priori estimate of Lemma 1. We assume that $`c_1(𝐱)`$ is a mean-zero stationary random process with a smooth rapidly decaying correlation function $`R(𝐱)`$. In the limit $`\delta 0`$ we have $`𝔼\left\{a^\delta \right\}a(t,𝐱,𝐤)`$, where $`a`$ solves the Fokker-Plank equation

$$\frac{a}{t}=\frac{}{k_m}\left(k^2D_{jm}(\widehat{𝐤})\frac{a}{k_j}\right)c_0\widehat{𝐤}_𝐱a$$
(11)
with the initial data $`a(0,𝐱,𝐤)=|\chi _\mathrm{\Omega }(𝐱)|^2\widehat{f}(𝐤)`$. Here, $`k=|𝐤|`$, and the diffusion matrix is

$$D_{jm}(\widehat{𝐤})=\frac{1}{2}_{\mathrm{}}^{\mathrm{}}\frac{^2R}{z_jz_m}(c_0\widehat{𝐤}s)𝑑s,$$

with $`\widehat{𝐤}=𝐤/k`$. The operator on the right side of (11) is degenerate because $`D_{ij}(𝐤)k_j=0`$. However, it is hypoelliptic on $`X=_𝐱^3\times 𝕊_𝐤^2`$ and thus possesses a smooth Green’s function $`G_L(t,𝐱,𝐤;𝐱^{},𝐤^{})`$ on $`X`$ so that

$$a(t,𝐱,𝐤)=_{^{2d}}G_L(t,𝐱,𝐤;𝐲,\widehat{𝐤}^{})|\chi _\mathrm{\Omega }(𝐲)|^2\widehat{f}(\widehat{𝐤}^{},k^{})𝑑𝐲𝑑𝐤^{}.$$

Hence we can introduce the new filter $`F_L=F_L^1+F_L^2`$ with

$$F_L^1(t,𝝃;𝐱_0)=_^de^{i𝐤𝝃}a(t,𝐱_0,𝐤)\mathrm{\Gamma }𝐛^1(𝐤)𝐛^1(𝐤)A_0\mathrm{\Gamma }𝑑𝐤,$$
(12)
and $`F_L^2`$ defined similarly. In the limit $`\epsilon 0`$ and $`\delta 0`$ we obtain the reconstruction (8) for $`𝔼\left\{u^B\right\}`$ with $`F=F_L`$. Now the function $`a(t,𝐱,𝐤)`$ is smoother in $`𝐤`$ than in the homogeneous case and thus the filter $`F_L`$ decays faster than $`F_H`$. This produces sharper refocusing and eliminates spurious Fresnel zones. The process $`a^\delta (t,𝐱,𝐤)`$ does not converge to the deterministic function $`a(t,𝐱,𝐤)`$ in probability pointwise in $`𝐱`$ and $`𝐤`$. However, for any test function $`\varphi (𝐱,𝐤)𝒮(^{2d})`$ the process $`\varphi ^\delta (t)=(a^\delta ,\varphi )`$ converges to $`(a,\varphi )`$ in probability. That means that the refocused signal $`u^B`$ is deterministic. This self-averaging of the refocused pulse follows from slight modifications of the analysis in .

## 4. Refocusing via the transport and diffusion equations

The radiative transport regime arises when the random medium has weak fluctuations at the wavelength scale. This case is more complicated than the one treated in Section 3 because the waves now fully interact with the medium. Propagation of wave correlations is described in terms of a Boltzmann-type radiative transport equation . The random fluctuations of the density and compressibility are assumed to be of the form $`\rho _\epsilon (𝐱)=\rho _0+\sqrt{\epsilon }\rho _1(\frac{𝐱}{\epsilon }),`$ $`\kappa _\epsilon (𝐱)=\kappa _0+\sqrt{\epsilon }\kappa _1(\frac{𝐱}{\epsilon }),`$ where $`\rho _0`$, $`\kappa _0`$ are constants, and $`\rho _1(𝐱)`$, $`\kappa _1(𝐱)`$ are mean-zero stationary random processes. A rigorous derivation of the transport equation for the Wigner distribution is only available for Schrödinger equations , Based on the formal analysis in  for hyperbolic systems, we can still conjecture that $`a_{1,2}(t,𝐱,𝐤)`$ solve a linear transport equation, at least after ensemble averaging over realizations of the random fluctuations. It is moreover known that for long distances of propagation and large times, the transport solution becomes independent of the direction $`\widehat{𝐤}`$ and is approximated by that of a diffusion equation:

$$_ta(t,𝐱,k)D(k)\mathrm{\Delta }_𝐱a(t,𝐱,k)=0,$$
(13)
where $`D(k)`$ is a diffusion coefficient that depends on the power spectrum of the random fluctuations (see  for instance). When the filter $`f`$ is radially symmetrical, the solution of (13) may be written as $`a_{1,2}(t,𝐱,k)=\psi (t,𝐱,k)\widehat{f}(k)`$. Here $`\psi (t,𝐱,k)`$ is the solution of (13) with initial data $`\psi (0,𝐱,k)=|\chi _\mathrm{\Omega }(𝐱)|^2`$. We then obtain that the new filter $`F_B`$ is scalar and is given by

$$F_B(t,𝝃;𝐱_0)=_^de^{i𝐤𝝃}\psi (t,𝐱_0,k)\widehat{f}(k)𝑑𝐤.$$
(14)
The retransmitted field (in the weak sense and after ensemble averaging) is given now by (8) with $`F=F_BI`$. Note that similarly to the random geometrical optics case treated in Section 3 the amplitude $`a(t,𝐱,𝐤)`$ is smoother in $`𝐤`$ than in the homogeneous case and the filter $`F_B`$ decays faster than $`F_H`$, thus produces a sharper refocusing than in a homogeneous medium.

Acknowledgement. We thank G. Papanicolaou and K. Solna for fruitful discussions. G. Bal was supported by NSF grant DMS-0072008 and L. Ryzhik by NSF grant DMS-9971742.


text:
# Effects of Transport Memory and Nonlinear Damping in a Generalized Fisher’s Equation

## I Introduction

Fisher’s equation murray  describes the dynamics of a field $`u(x,t)`$ subject to diffusive transport and logistic growth:

$$\frac{u}{t}=D\frac{^2u}{x^2}+ku(1u/K).$$
(1)
This kind of reaction-diffusion equation is relevant in chemical kinetics as well as in ecological contexts where $`u`$ diffuses via a diffusion constant $`D`$ and grows with a linear growth rate $`k`$ while the environment imposes a carrying capacity $`K`$. Equation (1) belongs to a family of single component models of broad applicability. Fisher proposed it as a deterministic model of the spread of a favored gene in a population fisher36 . Other systems in which Eq. (1) plays a significant role include flame propagation, neutron flux in a nuclear reactor and the dynamics of defects in nematic liquid crystals murray ; mikhailov ; canosa ; saarloos ; borzi . Work has been done on the propagation of wave-like structures in unbounded systems borzi ; hassan , as well as the formation of spatial structures in finite geometries izus .

Diffusion is typically a limit of a more coherent motion interrupted by scattering events which is valid when the scattering events are extremely frequent. A useful manner of describing situations in which the frequency of scattering events is not infinitely large is through the incorporation of memories with finite decay times kenkrereineker . An analytical investigation of reaction-diffusion phenomena based on the introduction of exponential memory in the context of the logistic nonlinearity of the Fisher’s equation, and a piecewise linearity approximation, was recently given by Manne, Hurd and Kenkre (MHK) MHK . A number of interesting results appeared in that analysis, including a generalization of the relation between the minimum speed of shock fronts and the system parameters, and the appearance of a speed limit beyond which the fronts would exhibit spatial oscillations.

As stated in Ref. MHK , the generalization of Fisher’s equation is not unique, the one employed in Ref. MHK  having been chosen for analytic tractability. A more natural generalization results in a more complex structure of the equation, involves a nonlinearity in the damping, and is the subject of the present paper. Our investigations show that, in the light of the present generalization, some of the interesting predictions of Ref. MHK  are modified drastically, while others are only slightly changed.

Reference MHK  first generalized the linear part of Eq. (1) by incorporating an exponential memory, converted the resulting integro-differential equation via differentiation into a differential (specifically the telegrapher’s) equation, and finally added the nonlinear logistic term to obtain the starting point of the analysis. Here, however, we take as our point of departure,

$$\frac{u}{t}=D_0^t\varphi (t\tau )\frac{^2u}{x^2}𝑑\tau +kf(u),$$
(2)
which reduces to Fisher’s equation if the nonlinearity is logistic, and if the memory function $`\varphi (t)`$ is a $`\delta `$-function in time. Generally, the decay time of $`\varphi (t)`$ is a measure of the time between scattering events: we will take $`\varphi (t)`$ to have a simple exponential form $`\varphi (t)=\alpha e^{\alpha t}`$, where $`1/\alpha `$ represents the scattering time.

Transformation of Eq. (2) into a differential equation is trivially done by differentiating once with respect to time:

$$\frac{^2u}{t^2}+[\alpha kf^{}(u)]\frac{u}{t}=v^2\frac{^2u}{x^2}+\alpha kf(u),$$
(3)
where, as in MHK , we put $`D\alpha =v^2`$, the physical meaning of $`v`$ being the speed dictated by the medium in the absence of scattering. This is the speed at which the underlying quasiparticle, whose number density or probability density is described by $`u(x,t)`$, moves ballistically (coherently) in between scattering events. If the coherent motion is interrupted too often by scattering, the motion looks diffusive and one returns to the Fisher limit.

Equation (3) becomes, in the absence of the nonlinearity, $`k=0`$, the well-known telegrapher’s equation morse  suggested by Lord Kelvin for the description of propagation of transatlantic telegraph signals. Equation (3) differs from the starting point of the MHK analysis, viz.,

$$\frac{^2u}{t^2}+\alpha \frac{u}{t}=v^2\frac{^2u}{x^2}+\alpha kf(u),$$
(4)
in a significant way. The former shows that the damping coefficient multiplying $`u/t`$ is not constant, and can even be negative.

Whereas Eq. (4) is obviously simpler in form than Eq. (3), the latter could be argued to stem from a more natural generalization of the Fisher equation, involving a simple replacement of a $`\delta `$-function memory by a finite decay time memory in the transport term $`D^2u/x^2`$ in (1). We are interested in comparing the physical consequences of the two generalizations in the context of traveling wave solutions. We do this in two steps. In Section II, we give a general analysis of traveling wave solutions without approximating the equation in any way. We investigate the minimum speed of wavefronts, the existence of spatially oscillatory fronts given in MHK, and give a qualitative description of all the kinds of front shapes that the system can support as a function of the parameters. In a certain regime in parameter space, a dynamical stabilization of the unstable state of the logistic equation is observed, with fronts of state $`u=0`$ invading the state $`u=K`$. In Section III, we replace the nonlinearity by a piecewise linear approximate form and obtain analytic solutions following the method of MHK. We end with concluding remarks in Section IV.

## II Nonlinear analysis

We look for waves moving in the direction of increasing $`x`$: $`u(x,t)=KU(xct)=KU(z)`$, where $`c`$ is the speed of the nonlinear wave, generally different from the speed of linear waves $`v`$, dictated by the medium. We are also renormalizing with respect to the carrying capacity $`K`$ to simplify the reaction term in the logistic case. With this ansatz, we obtain from (3) the following ordinary differential equation:

$$(v^2c^2)U^{\prime \prime }+c[\alpha kf^{}(U)]U^{}+\alpha kf(U)=0.$$
(5)
For a logistic reaction term: $`f(U)=U(1U)`$ and $`f^{}(U)=12U`$, so that Eq. (5) becomes:

$$mU^{\prime \prime }+c(\alpha k+2kU)U^{}+\alpha kU(1U)=0,$$
(6)
where we have followed the notation of MHK, $`m=v^2c^2`$, to emphasize the formal similarity between (6) and the equation of motion of a damped oscillator of mass $`m`$, subject to the nonlinear force $`\alpha kU(1U)`$.

A convenient way to analyze the solutions of (6) is to write it as a first order system:

$$\begin{array}{ccc}\hfill U^{}& =V\hfill & f(U,V)\hfill \\ \hfill V^{}& =\frac{1}{m}[(k\alpha 2kU)cV\alpha kU(1U)]\hfill & g(U,V).\hfill \end{array}$$
(7)
The system (7) has two equilibria: $`(U^{},V^{})=(0,0)`$ and $`(U^{},V^{})=(1,0)`$. We can analyze the character of these by looking at the linear behavior in a neighborhood of the equilibria.

At the equilibrium $`(0,0)`$ we have the eigenvalues:

$$\lambda _{}=\frac{c(k\alpha )\sqrt{(k\alpha )^2c^24k\alpha m}}{2m}.$$
(8)
If we assume that $`U`$ is a density or a concentration, solutions where $`U`$ oscillates below $`U=0`$ are not allowed. This imposes a condition on the eigenvalues $`\lambda _{}`$ to be real, from which we obtain:

$$cc_{min}=v\frac{1}{\sqrt{1+\frac{1}{4}(y1/y)^2}},$$
(9)
where $`y=\sqrt{\alpha /k}`$. This relation states that there is a minimum value of the speed of the nonlinear waves. It is to be compared to the known result in the context of Fisher’s equation: $`c_{min}=2\sqrt{kD}`$ (see for example murray ). If we make the formal identification $`D\alpha =v^2`$, as in (3), we have:

$$cc_{min}=2v\sqrt{1/y}.$$
(10)
Equation (9) is also to be compared to the previously found result of Ref. MHK :

$$cc_{min}^{MHK}=v\frac{1}{\sqrt{1+\frac{1}{4}y^2}},$$
(11)
Figure 1 displays a comparison of the minimum speeds $`c_{min}`$ as a function of the system parameter $`\alpha /k`$ in the three cases: purely diffusive (Fisher limit), the MHK generalization, and the present generalization. It can be seen that the two generalizations approach asymptotically the behavior of the purely diffusive situation for large values of $`\alpha /k`$. This is to be expected since $`\alpha \mathrm{}`$, $`v\mathrm{}`$, $`v^2/\alpha =D`$ is the diffusive limit of Eq. (2)\]. At low values of the ratio $`\alpha /k`$, however, there is a notable difference. The finite correlation time generalizations allow waves with lower speed than those present in Fisher’s equation. Moreover, the present generalization predicts a sharply different behavior in the region $`\alpha <k`$, due to the anti-damping present in Eq. (6). This growing branch of the $`c_{min}`$ function turns out to have drastic consequences in the nature of the traveling fronts when the system is far from the diffusive regime. On the one hand, it can be seen that waves with speed $`c=v`$ are obtained at a finite value of the memory constant $`\alpha `$, namely $`\alpha =k`$. When approaching the ballistic limit beyond this point, traveling waves are allowed with speeds lower than $`v`$. The nature of these will be analyzed in the following discussion.

Around the equilibrium $`(1,0)`$ the eigenvalues are:

$$\mu _{}=\frac{c(k+\alpha )\sqrt{(k+\alpha )^2c^2+4k\alpha m}}{2m}.$$
(12)
Oscillations around this equilibrium are in principle possible, as observed in MHK , since there is no impediment to the solution growing beyond $`U=1`$. The condition for oscillations is found from the radicand of the eigenvalues $`\mu _{}`$:

$$(k+\alpha )^2c^2+4k\alpha (v^2c^2)<0,$$
(13)
which, after some manipulation, becomes

$$(\alpha k)^2c^2+4k\alpha v^2<0.$$
(14)
Both terms in this relation are positive, implying that the relation cannot be satisfied by the wave speed for any set of parameters.

A similar analysis can be carried out with the equation (4), corresponding to MHK model. The character of the equilibria for both generalizations is resumed in Tables 1 and 2. Note that there is a $`c_{osc}`$ in Table 2 Case (7), not present in the results resumed in Table 1 because Eq. (14) cannot be satisfied. Different combinations of parameters allow a variety of traveling fronts connecting one equilibrium to the other. Cases (1) to (7) in Table 2 are the fully nonlinear equivalents of the piecewise linear situations analyzed in MHK . The situations resumed in Table 1 are discussed below.

Equation. (6) can be exploited directly to understand the nature of the solutions of the traveling wave problem, without obtaining analytical solutions via simplifications. We will base the following discussion on the mechanical interpretation of Eq. (6), which can be reinterpreted, following MHK, as describing the motion of a particle of mass $`m`$ in a nonlinear potential, subject to a state-dependent damping that can be positive or negative. Since the mass of the particle can also be positive or negative, depending on whether the velocity $`c`$ is lower than $`v`$ or not, we analyze the two cases separately.

### II.1 $`c<v`$

In this case a particle of mass $`m>0`$ is moving in a potential $`\varphi (U)=\alpha k(U^2/2U^3/3)`$ \[see Fig. 2 (a)\]. This potential has a minimum at $`U=0`$, but we have ruled out solutions that oscillate around it based on the positivity of $`U`$ \[this are cases (1) and (4) in Table 1\]. Another possible solution is an overdamped trajectory that connects the equilibrium at $`U=1`$ with that at $`U=0`$. This corresponds to a traveling front of the state 1 invading the system at state 0. The damping coefficient is $`\gamma =c(\alpha k+2kU)`$. We can see that if $`\alpha >k`$, then $`\gamma >0`$ for all values of $`U`$ along this trajectory. This is a solution that connects the saddle with the stable node shown in the case (2) of Table 1.

Now, if $`\alpha <k`$, the motion is damped on part of the trajectory \[when $`U>(1\alpha /k)/2`$, see Fig. 2 (c)\] and anti-damped on the rest, specifically near $`U=0`$. This means that the equilibrium at $`U=0`$ is unstable, and there is a solution connecting it to the equilibrium at $`U=1`$. It is a trajectory connecting the unstable node at $`(0,0)`$ to the stable manifold of the saddle at $`(1,0)`$ \[see Table 1, case (5)\]. This trajectory corresponds to a front of state 0 invading the state 1.

### II.2 $`v<c`$

In this situation $`m`$ is negative, so we multiply Eq. (6) by $`1`$, to obtain:

$$|m|U^{\prime \prime }c(\alpha k+2kU)U^{}=\alpha kU(1U).$$
(15)
Now the particle of mass $`|m|`$ moves in the potential $`\varphi (U)=\alpha k(U^2/2U^3/3)`$ \[refer to Fig. 2 (b)\], where there is a stable equilibrium at $`U=1`$. The sign of the damping coefficient is also reversed, and the motion is always anti-damped when $`\alpha >k`$. When $`\alpha <k`$, it is always anti-damped in the vicinity of $`U=1`$. It can be seen in Table 1 \[cases (3) and (6)\] that the equilibrium at $`(0,0)`$ is always a saddle, and the equilibrium at $`(1,0)`$ is always an unstable node. There is a trajectory that connects the unstable node to the stable manifold of the saddle, corresponding to a front of state 1 invading the state 0.

We note that, contrary to the generalized telegrapher’s equation studied in MHK , we do not have here oscillating wave shapes (corresponding to the unstable spiral at $`(1,0)`$ shown in Table 2). Also to be remarked are the solutions where the state $`U=0`$—which is an unstable state of the logistic equation— invades the state $`U=1`$. This stabilization of the null state has been made possible by interplay of memory and reaction, through the state-dependent damping coefficient of Eq. (3).

A phase diagram summarizing these results is shown in Fig. 3. A qualitative indication of the nature of the front wave is given, together with a reference to the cases enumerated in Table 1.

## III Piecewise linearization

Explicit solutions corresponding to the front waves described above cannot be found in the fully nonlinear situation. However, a piecewise linearization of Eq. (3) can be made, just as in MHK , taking the following as a reaction function:

$$f(U)=\{\begin{array}{cc}U/a,\hfill & Ua,\hfill \\ (bU)/(ba),\hfill & Ua,\hfill \end{array}$$
(16)
where $`a<b`$. With this reaction term, which incidentally, generalizes somewhat the logistic form (which corresponds to $`a=b`$), the oscillator equations in the traveling waves ansatz (5) become:

$`mU^{\prime \prime }+2\gamma _1U^{}+k_1^2U`$ $`=`$ $`0,Ua,`$ (17)
$`mU^{\prime \prime }+2\gamma _2U^{}+k_2^2(bU)`$ $`=`$ $`0,Ua.`$ (18)
where

$$\gamma _1=\frac{c}{2}\left(\alpha \frac{k}{a}\right),\gamma _2=\frac{c}{2}\left(\alpha \frac{k}{ba}\right),$$
(19)
$$k_1^2=\frac{\alpha k}{a},k_2^2=\frac{\alpha k}{ba}.$$
(20)
Suppose that we are looking for a front solution that interpolates from $`U=b`$ as $`z\mathrm{}`$ to $`U=0`$ as $`z\mathrm{}`$. Without loss of generality, we take $`U=a`$ at $`z=0`$. We look for solutions of two oscillators: one to the right of $`z=0`$, and one to the left of $`z=0`$. We define the following variables:

$$\begin{array}{cccc}U_R(z)\hfill & =& U(z),\hfill & z0,\hfill \\ U_L(z)\hfill & =& bU(z),\hfill & z0.\hfill \end{array}$$
(21)
Equations (17) and (18) become:

$`mU_R^{\prime \prime }+2\gamma _1U_R^{}+k_1^2U_R`$ $`=`$ $`0,`$ (22)
$`mU_L^{\prime \prime }+2\gamma _2U_L^{}+k_1^2U_L`$ $`=`$ $`0.`$ (23)
The wave front shape, which is a solution to Eqs. (17) and (18), can be easily found in terms of exponentials, by matching solutions of Eqs. (22) and (23) and their derivatives on either side of $`z=0`$. In some cases, one of the solutions of Eqs. (22) or (23) contains a growing exponential, and the corresponding factor has to be set to zero to avoid unbounded growth. Examples of all the cases described in the previous section are shown in Fig. 4, within this piecewise linear scheme. The solutions have been arbitrarily shifted in the variable $`z=xct`$ to avoid superposition of the curves. In the model linearized according to (16), the regimes of damping and anti-damping of Eq. (5) are separated by the condition $`\alpha =2k`$, which is used in the figure caption to classify the curves. (The condition separating the damping from the anti-damping regimes is $`\alpha =k/a`$ in the piecewise linear model, instead of the $`\alpha =k`$ of the fully nonlinear model.) The case numbers corresponding to Table 1 are also shown in the caption.

## IV Conclusions

Reaction-diffusion systems in which the transport process is wave-like at short times and diffusive at long times are the focus of the present investigation. This passage of the character of the motion from coherent to incoherent is a general feature of all physical systems and may be represented by a memory function whose decay time (or correlation time) represents the demarcation. We have analyzed here a generalization of Fisher’s reaction-diffusion equation which has a logistic nonlinearity describing the reaction process. Our study has centered on traveling wave solutions. From arguments without approximation we have found a generalization of the known Fisher’s equation result regarding the minimum speed of the traveling waves. While the generalization converges to the Fisher’s equation result in the limit of diffusive transport (memory that decays infinitely fast), sharp differences occur in the wave-limit (see Fig. 1): new kinds of solutions are found to be possible involving “inverse” fronts, in which the state $`U=0`$ invades the state $`U=1`$ (see Fig. 3). We have also found explicit analytic solutions (typical cases plotted in Fig. 4) via the piecewise linearization introduced for this problem elsewhere MHK  and pointed out a number of differences in the predictions arising in our present analysis. Some of these differences are expected but at least one of them is surprising: solutions which oscillate spatially when they have a speed above a certain value (called $`c_{osc}`$ in Ref. MHK ; see also Table 2) predicted in Ref. MHK  are found to disappear in the present analysis which can be argued to be a more natural generalization of Fisher’s equation to include wave-like transport. There are a number of issues, such as the stability of the solutions against perturbations introduced into the system, which will be studied in future investigations.

We mention in passing the problem of what kind of initial conditions eventually evolve into a traveling front. The general problem can be very difficult, but a simplified treatment is possible as follows. Consider that the initial state has the behavior $`u(x,0)Ae^{ax}`$ for $`x\mathrm{}`$. Correspondingly, we suppose that the leading edge of the traveling wave has the form:

$$u(x,t)=Ae^{a(xct)},$$
(24)
with $`a`$ and $`A`$ arbitrary.

We substitute this into the differential equation, and suppose that, at the leading edge, $`u0`$:

$$(ac)^2u+\alpha acukacu+o(u^2)=v^2a^2u+\alpha ku+o(u^2).$$
(25)
Disregarding the terms in $`u^2`$ we arrive at a dispersion relation between $`c`$ and $`a`$:

$$c=\frac{k\alpha }{2a}\pm \sqrt{\frac{(\alpha +k)^2}{4a^2}+v^2}.$$
(26)
The negative sign before the square root in this expression leads to negative $`c`$, and can be disregarded. The expression with the positive sign gives the speed of the front wave that will eventually develop from an initial condition that has an exponential decay of coefficient $`a`$.

###### Acknowledgements.

This work was supported in part by the Los Alamos National Laboratory via a grant made to the University of New Mexico (Consortium of the Americas for Interdisciplinary Science) and by the National Science Foundation’s Division of Materials Research via grant No. DMR0097204. G. A. thanks the support of the Consortium of the Americas for Interdisciplinary Science and the hospitality of the University of New Mexico.