format: text

text:
# Persistence of Manifolds in Nonequilibrium Critical Dynamics

## Abstract

We study the persistence $`P(t)`$ of the magnetization of a $`d^{}`$ dimensional manifold (i.e. the probability that the manifold magnetization does not flip up to time $`t`$, starting from a random initial condition) in a $`d`$-dimensional spin system at its critical point. We show analytically that there are three distinct late time decay forms for $`P(t)`$ : exponential, stretched exponential and power law, depending on a single parameter $`\zeta =(D2+\eta )/z`$ where $`D=dd^{}`$ and $`\eta `$, $`z`$ are standard critical exponents. In particular, our theory predicts that the persistence of a line magnetization decays as a power law in the $`d=2`$ Ising model at its critical point. For the $`d=3`$ critical Ising model, the persistence of the plane magnetization decays as a power law, while that of a line magnetization decays as a stretched exponential. Numerical results are consistent with these analytical predictions.

PACS numbers: 05.70.Ln, 05.50.+q, 05.70.Jk

Following a rapid quench of a spin system from infinite temperature to zero temperature, domains of competing ground states form and grow with time. In the standard nonconserved dynamics a given spin inside the sample flips only when a domain wall passes through it, which happens rather rarely. As a result, persistence, i.e. the probability that the spin remains unflipped up to time $`t`$ decays slowly as a power law, $`P(t)t^\theta `$ at late times. In contrast, if the spin system is quenched to its critical temperature $`T_c`$, any given spin fluctuates rapidly due to the finite temperature and the persistence of a single spin has an exponential tail at late times. However, for a quench to $`T_c`$, the persistence of the total magnetization of the sample (as opposed to that of a single spin) again decays as a power law, $`P(t)t^{\theta _g}`$, where the exponent $`\theta _g`$ has been argued to be a new nonequilibrium critical exponent. The global persistence has since been studied in a wide variety of systems.

In a $`d`$-dimensional sample, a single spin is a $`0`$-dimensional manifold. On the other hand, the global magnetization corresponds to summing over all the spins of the full sample which is a $`d`$-dimensional manifold. To interpolate between these two limits, it is natural to study the persistence of the magnetization of a $`d^{}`$-dimensional manifold with $`0d^{}d`$. In the two limits $`d^{}=0`$ and $`d^{}=d`$, the persistence has very different asymptotic decay, respectively exponential and power law. The following question then naturally arises: if one tunes the manifold dimension $`d^{}`$ from $`d^{}=0`$ to $`d^{}=d`$, how does the asymptotic behavior of persistence changes from an exponential ($`d^{}=0`$) to a power law decay ($`d^{}=d`$)? Does this change occur abruptly at some intermediate value of $`d^{}`$, or is there an intermediate regime of $`d^{}`$ where the behavior is neither exponential nor power law, but something in between? In this Letter we address this interesting issue and we show analytically that indeed there is a novel intermediate regime of $`d^{}`$ where the persistence has a stretched exponential tail. Our main results can be summarized in terms of a single number $`\zeta =(D2+\eta )/z`$, where $`D=dd^{}`$ is the co-dimension of the manifold and $`\eta `$ and $`z`$ are the standard critical exponents: $`z`$ is the dynamical exponent that describes the temporal growth of the correlation length $`\xi (t)t^{1/z}`$ and $`\eta `$ describes the power law decay of the equilibrium spin-spin correlation function at $`T_c`$, $`\varphi (0)\varphi (r)r^{(d2+\eta )}`$ for large separation $`r`$. Depending on the value of $`\zeta `$, the persistence of the magnetization of a $`d^{}`$-dimensional manifold at $`T_c`$ has the asymptotic behavior,

$`P(t)`$ $``$ $`t^{\theta (d^{},d)},\zeta <0`$ (1)
$``$ $`\mathrm{exp}\left(a_1t^\zeta \right),\mathrm{\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}0}\zeta 1`$ (2)
$``$ $`\mathrm{exp}\left(b_1t\right),\zeta >1,`$ (3)
where the exponent $`\theta (d^{},d)`$ depends on $`d^{}`$ and $`d`$ and $`a_1`$, $`b_1`$ are constants. Strictly speaking, in the intermediate regime $`0\zeta 1`$ we will show that, $`\mathrm{exp}\left(a_2t^\zeta \mathrm{ln}t\right)P(t)\mathrm{exp}\left(a_1t^\zeta \right)`$ for large $`t`$. We will derive the results in Eq. (3) within the mean-field theory (valid for $`d>4`$), in the $`n\mathrm{}`$ limit of the $`O(n)`$ model, followed by a general scaling theory.

Two specific applications of our general results are as follows: (i) Consider the persistence of the line magnetization ($`d^{}=1`$) in the $`d=2`$ Ising model at $`T_c`$. Using $`d^{}=1`$, $`d=2`$, $`\eta =1/4`$ and $`z2.172`$, one gets $`\zeta =(dd^{}2+\eta )/z0.3453<0`$. Hence Eq. (3) predicts a power law decay for the persistence of the line magnetization; (ii) For the $`d=3`$ Ising model, using $`\eta 0.032`$ and $`z2`$, one finds that while for the plane magnetization ($`d^{}=2`$), $`\zeta 0.484<0`$, for the line magnetization ($`d^{}=1`$), $`0<\zeta 0.016<1`$. Our Eq. (3) then predicts a power law decay for the persistence of the plane magnetization, but a stretched exponential decay for the line magnetization. Numerical simulations for small samples in $`2`$ and $`3`$ dimensions are consistent with these analytical predictions.

Our starting point is the standard Langevin equation for the vector order parameter $`\stackrel{}{\varphi }=(\varphi _1,\mathrm{},\varphi _n)`$,

$$_t\varphi _i=^2\varphi _ir\varphi _i(u/n)\varphi _i\underset{j=1}{\overset{n}{}}\varphi _j^2+\eta _i,$$
(4)
where $`\stackrel{}{\eta }(𝐱,t)`$ is a Gaussian white noise with mean zero and correlator $`\eta _i(𝐱,t)\eta _j(𝐱^{},t^{})=2\delta _{i,j}\delta ^d(𝐱𝐱^{})\delta (tt^{})`$. The magnetization of a $`d^{}`$-dimensional manifold will be defined by the vector field, $`\psi _i(x_{d^{}+1},\mathrm{},x_d,t)=\varphi _i(\stackrel{}{x},t)_{i=1}^d^{}\frac{dx_i}{\sqrt{L}}`$, obtained by integrating the order parameter over the $`d^{}`$ directions. Here $`L`$ denotes the length of the sample along any direction. For a vector order parameter, we define the persistence $`P(t)`$ of the manifold to be the probability that any given component of the manifold magnetization, say $`\psi _1`$, does not change sign up to time $`t`$. Since all the components of the spin are equivalent, henceforth we will drop the subscript $`i`$ of $`\psi _i`$ for convenience. Note that the magnetization $`\psi `$ is a field over the remaining $`D=dd^{}`$ dimensional space whose coordinates $`(x_{d^{}+1},\mathrm{},x_d)`$ will be relabelled by the vector $`\stackrel{}{r}=(r_1,r_2,\mathrm{}r_D)`$ for convenience. The observation that allows us to make analytical predictions for the persistence of the magnetization $`\psi (\stackrel{}{r},t)`$ is that it is a Gaussian variable at all finite times. This follows simply from the fact that $`\psi (\stackrel{}{r},t)`$ is a sum of $`L^d^{}`$ random variables which are correlated but only over a finite correlation length $`\xi (t)t^{1/z}`$. Thus in the thermodynamic limit when $`t^{1/z}L`$, the central limit theorem asserts that $`\psi (\stackrel{}{r},t)`$ is a Gaussian field. Hence its persistence is determined by the autocorrelation function $`C(t_1,t_2)=\psi (\stackrel{}{r},t_1)\psi (\stackrel{}{r},t_2)`$.

We start with the mean-field theory, valid for $`d4`$, where we set $`u=0`$, and also $`r=0`$ (at the critical point) in Eq. (4). Next we integrate the Langevin equation over the $`d^{}`$ space directions and then solve the resulting linear equation in the Fourier space. Defining $`\stackrel{~}{\psi }(\stackrel{}{k},t)=\psi (\stackrel{}{r},t)e^{i\stackrel{}{k}\stackrel{}{r}}d^Dr`$, it is easy to compute the two-time correlation function,

$`\stackrel{~}{\psi }(\stackrel{}{k},t_1)\stackrel{~}{\psi }(\stackrel{}{k},t_2)`$ $`=`$ $`\mathrm{\Delta }(\stackrel{}{k})e^{k^2(t_1+t_2)}`$ (5)
$`+`$ $`{\displaystyle \frac{1}{k^2}}\left(e^{k^2|t_1t_2|}e^{k^2(t_1+t_2)}\right),`$ (6)
where $`\mathrm{\Delta }(\stackrel{}{k})=\stackrel{~}{\psi }(\stackrel{}{k},0)\stackrel{~}{\psi }(\stackrel{}{k},0)`$ is taken to be a constant (for uncorrelated initial condition). At late times, the initial condition dependent term becomes negligible and it is sufficient to retain only the second term on the right-hand side of Eq. (6). The autocorrelation function is then obtained by integrating over the $`\stackrel{}{k}`$ space, $`C(t_1,t_2)=d^Dke^{k^2a^2}\stackrel{~}{\psi }(\stackrel{}{k},t_1)\stackrel{~}{\psi }(\stackrel{}{k},t_2)`$ where we have introduced a soft ultraviolet cut-off $`a`$. One finds

$$C(t_1,t_2)=\beta A\left[(t_1+t_2+a^2)^{2\beta }(|t_1t_2|+a^2)^{2\beta }\right],$$
(7)
where $`2\beta =(2D)/2`$ and $`A`$ is an unimportant constant.

Consider first the case $`D<2`$, i.e., $`\beta >0`$. In this case there is no need for the ultraviolet cut-off since $`\psi ^2(\stackrel{}{r},t)=C(t,t)`$ does not diverge at any finite time even if one puts $`a=0`$. Setting $`a=0`$ in Eq. (7) we find that the correlation function $`C(t_1,t_2)`$ is still non-stationary. However, one can render it stationary by employing a well known trick, where one introduces a normalized Gaussian process, $`X=\psi /\sqrt{\psi ^2}`$ which, when observed in the logarithmic time $`T=\mathrm{ln}t`$ becomes a stationary Gaussian variable with the correlator, $`C(T)=X(0)X(T)=\mathrm{cosh}(T/2)^{2\beta }|\mathrm{sinh}(T/2)|^{2\beta }`$. Interestingly, the same Gaussian correlator also appears in the context of the persistence of the Edwards-Wilkinson type rough interfaces. It is well known that for such a stationary Gaussian correlator (decaying exponentially for large $`T`$), the persistence of the process also decays exponentially for large $`T`$, $`P(T)\mathrm{exp}(\theta T)`$. In terms of the real time $`t=e^T`$, this indicates a power law decay, $`P(t)t^\theta `$, with persistence exponent $`\theta `$. For this particular correlator, the exponent $`\theta `$ has been studied in great detail in the context of the interface problem and the exponent is known to depend continuously on the roughness parameter $`\beta =(2D)/4>0`$.

For $`D>2`$, on the other hand, one needs the ultraviolet cut-off $`a`$ explicitly in order to keep $`\psi ^2(\stackrel{}{r},t)=C(t,t)`$ finite. In that case, the appropriate scaling limit is $`t_1,t_2\mathrm{}`$, but keeping their difference $`|t_1t_2|`$ fixed. In this limit, Eq. (7) reduces to a stationary correlator in the original time variable, $`C(t_1,t_2)(|t_1t_2|+a^2)^{(D2)/2}`$ that decays as a power law for large $`|t_1t_2|`$. To calculate the persistence of Gaussian stationary processes with an algebraically decaying correlator is nontrivial. However, there exists a powerful theorem due to Newell and Rosenblatt, which states that if the stationary correlator decays as $`C(t)t^\alpha `$ with $`\alpha >0`$ for large time difference $`t=|t_1t_2|`$, then the persistence $`P(t)`$ (probability of no zero crossing between $`t_1`$ and $`t_2`$) of such a process has the following asymptotic behaviors: (i) $`P(t)\mathrm{exp}(K_1t)`$ if $`\alpha >1`$ and (ii) $`\mathrm{exp}(K_2t^\alpha \mathrm{ln}t)P(t)\mathrm{exp}(K_3t^\alpha )`$ if $`0<\alpha <1`$, where the $`K_i`$’s are constants. Applying this theorem to our problem, we find that $`P(t)\mathrm{exp}(K_1t)`$ for $`D>4`$ and $`\mathrm{exp}(K_2t^{(D2)/2}\mathrm{ln}t)P(t)\mathrm{exp}(K_3t^{(D2)/2})`$ for $`2<D<4`$. In the borderline case $`D=4`$, there will be an additional logarithmic correction.

Combining these results for $`D<2`$ and $`D>2`$ and noting that $`z=2`$ and $`\eta =0`$ within the mean-field theory, we find that the explicit exact results for the mean-field theory derived above are just the special cases of the general result in Eq. (3) provided one uses the mean-field value $`\zeta =(D2)/2`$ in Eq. (3).

The mean-field theory is valid for $`d4`$. In order to access the physically relevant dimensions $`d4`$, we now consider another solvable limit, namely the $`n\mathrm{}`$ limit where Eq. (4) becomes

$$_t\varphi _i=^2\varphi _i[r+S(t)]\varphi _i+\eta _i,$$
(8)
where $`S(t)=u\varphi _i^2`$ has to be determined self-consistently. The critical point corresponds to $`r+S(\mathrm{})=0`$. This self-consistent determination of $`S(t)`$ can be done using standard techniques and one finds that at late times, $`S(t)S(\mathrm{})(4d)/4t`$ for $`2<d4`$. Substituting this result into Eq. (8), summing over the $`d^{}`$ directions, solving the resulting equation in the Fourier space and finally integrating over the $`\stackrel{}{k}`$ space (as in the mean-field theory), we finally arrive at the following autocorrelation function for the manifold magnetization $`\psi (\stackrel{}{r},t)`$ in $`2<d4`$,

$$C(t_1,t_2)=A_1(t_1t_2)^{(4d)/4}_0^{t_1}\frac{dt^{}t_{}^{}{}_{}{}^{(d4)/2}}{(t_1+t_22t^{}+a^2)^{D/2}},$$
(9)
where $`t_1t_2`$, $`A_1`$ is a constant and $`a`$ represents the soft ultraviolet cut-off as before.

For $`D<2`$, as in the mean-field theory, one can set the cut-off $`a=0`$ and the resulting nonstationary correlator in Eq. (9) can be reduced to a stationary correlator for the normalized process $`X=\psi /\sqrt{\psi ^2}`$ in the logarithmic time $`T=\mathrm{ln}t`$,

$$A(T)=\left[\mathrm{cosh}(T/2)\right]^{\mu D/2}\frac{B[\mu ,2\beta ,2/(1+e^T)]}{B[\mu ,2\beta ]},$$
(10)
where $`\mu =(d2)/2`$, $`2\beta =(2D)/2`$, $`B[m,n]`$ is the standard Beta function and $`B[m,n,x]=_0^x𝑑yy^{m1}(1y)^{n1}`$. Since the stationary correlator in Eq. (10) decays exponentially for large $`T`$, $`A(T)\mathrm{exp}[(d+D2)T/4]`$, one concludes  that the corresponding persistence also decays exponentially for large $`T`$, $`P(T)\mathrm{exp}(\theta T)`$ and hence as a power law in the original $`t=e^T`$ variable, $`P(t)t^\theta `$. Determining the exponent $`\theta `$ analytically is still a challenging task. However one can make progress in the limit of small co-dimension $`D0`$. Note that for $`D=0`$, i.e. for the global persistence, Eq. (10) becomes a pure exponential $`A(T)=\mathrm{exp}[(d2)T/4]`$ for all $`T`$, indicating that the process is Markovian. One then finds $`P(T)\mathrm{exp}(\theta _0T)`$ with $`\theta _0=(d2)/4`$. For $`D`$ nonzero but small, one can expand the correlator in Eq. (10) around the Markov process ($`D=0`$) and then use a perturbation theory result to calculate $`\theta `$ to first order in $`D`$. We get $`\theta =\theta _0+D\theta _0^2I_d/\pi +O(D^2)`$ for all $`2<d4`$, where $`I_d`$ is a complicated $`d`$-dependent integral. For special values of $`d`$, this integral simplifies. For example, for $`d=4`$, we get, $`\theta =1/2+(2\sqrt{2}1)D/4+O(D^2)`$ and for $`d=3`$, $`\theta =1/4+0.183615\mathrm{}D+O(D^2)`$.

For $`D>2`$, it is evident from Eq. (9) that one needs to keep a nonzero cut-off $`a`$ in order that the integral does not diverge at $`t_1=t_2`$. In this case, for large $`t_1`$, the dominant contribution to the integral comes from the regime $`t^{}t_1`$. It is easy to see that in the limit when $`t_1,t_2`$ are both large with their difference $`|t_1t_2|`$ fixed, the autocorrelator in Eq. (9) reduces to a stationary one, $`C(t_1,t_2)B_1(|t_1t_2|+a^2)^{(D2)/2}`$, where $`B_1`$ is an unimportant constant. One can then invoke the Newell-Rosenblatt theorem once again to conclude that for large $`t`$, the persistence $`P(t)\mathrm{exp}(\kappa _1t)`$ if $`D>4`$ and $`\mathrm{exp}(\kappa _2t^{(D2)/2}\mathrm{ln}t)P(t)\mathrm{exp}(\kappa _3t^{(D2)/2})`$ for $`2<D<4`$, where the $`\kappa _i`$’s are constants. Combining these results we thus find that the $`n\mathrm{}`$ results for the persistence of manifold magnetization in $`2<d4`$ are also compatible with our general results in Eq. (3) on noting that $`\zeta =(D2)/2`$ since $`\eta =0`$ and $`z=2`$ within the large $`n`$ limit.

Taking hints from the two solvable cases above we now construct a general scaling theory valid for all $`d2`$. The two-point correlation function of the order parameter, at the critical point, has the generic scaling form, $`\varphi (\mathrm{𝟎},t_1)\varphi (𝐱,t_2)x^{(d2+\eta )}F(xt_1^{1/z},t_2/t_1)`$ for large distance $`x`$ and large times $`t_1`$, $`t_2`$, where $`\eta `$ and $`z`$ are the standard critical exponents defined in the introduction. This indicates that in the Fourier space, $`\stackrel{~}{\varphi }(𝐊,t_1)\stackrel{~}{\varphi }(𝐊,t_2)K^{(2\eta )}G(Kt_1^{1/z},t_2/t_1)`$, where $`𝐊`$ is a $`d`$-dimensional vector conjugate to $`𝐱`$. The manifold magnetization $`\psi `$ is obtained by summing the order parameter $`\varphi `$ over $`d^{}`$ directions. This is equivalent to putting $`K_i=0`$ along the $`i=1,\mathrm{},d^{}`$ directions. One then obtains the scaling behavior of the two point correlator of the manifold magnetization, $`\stackrel{~}{\psi }(\stackrel{}{k},t_1)\stackrel{~}{\psi }(\stackrel{}{k},t_2)k^{(2\eta )}g(kt_1^{1/z},t_2/t_1)`$, where $`\stackrel{}{k}`$ is now a $`D=dd^{}`$ dimensional vector. For example, within the mean-field theory, the scaling function $`g(x,y)=\mathrm{exp}[x^2|1y|]\mathrm{exp}[x^2(1+y)]`$, as evident from Eq. (6) after dropping the $`\mathrm{\Delta }`$ dependent term at late times. The autocorrelation function $`C(t_1,t_2)=\psi (\stackrel{}{r},t_1)\psi (\stackrel{}{r},t_2)`$ is then obtained by integrating over $`\stackrel{}{k}`$,

$$C(t_1,t_2)=\frac{d^Dk}{k^{2\eta }}g(kt_1^{1/z},t_2/t_1)e^{k^2a^2},$$
(11)
where $`a`$ is the soft ultraviolet cut-off as before.

Consider first the case when $`D2+\eta <0`$. One can then set the cut-off $`a=0`$ (since the integral in Eq. (11) is convergent at the upper limit), and one obtains $`C(t_1,t_2)t_1^{(D2+\eta )/z}f(t_2/t_1)`$ in the limit $`t_1,t_2\mathrm{}`$ with $`t_2/t_1`$ arbitrary. Note that the function $`f(x)x^{\lambda _c/z}`$ for large $`x`$ such that $`C(t_1,t_2)t_2^{\lambda _c/z}`$ for $`t_2t_1`$ where $`\lambda _c`$ is the standard autocorrelation exponent. This nonstationary Gaussian correlator can then be reduced, as before, to a stationary one for the normalized variable $`X=\psi /\sqrt{\psi ^2}`$ in the logarithmic time, $`T=\mathrm{ln}t`$ and one gets, $`A(T)=X(0)X(T)=\mathrm{exp}[(D2+\eta )T/2z]f(e^T)/f(1)`$. Since $`A(T)\mathrm{exp}\left[(\lambda _c(D2+\eta )/2)T/z\right]`$ for large $`T`$, it follows, as before, that the persistence $`P(T)\mathrm{exp}(\theta T)`$ for large $`T`$. This means that the persistence decays as a power law in the original time variable $`t=e^T`$, $`P(t)t^\theta `$ for large $`t`$.

In the complementary case $`D2+\eta >0`$, the integral in Eq. (11) is, for $`t_1=t_2`$, divergent near the upper limit without the cut-off. Hence one needs to keep $`a`$ nonzero and then the appropriate scaling limit is obtained by taking $`t_1,t_2`$ both large keeping their difference $`|t_1t_2|`$ fixed but arbitrary. Then one can replace the scaling function $`g(kt_1^{1/z},t_2/t_1)`$ in Eq. (11) by another function $`g_1(k|t_1t_2|^{1/z})`$ of a single scaling variable, as in the two previous solvable cases. Performing the integral, one then finds $`C(t_1,t_2)|t_1t_2|^{(D2+\eta )/z}`$ for $`|t_1t_2|a^2`$. This correlator is stationary and decays as a power law. Invoking the Newell-Rosenblatt theorem once more, we find that $`P(t)`$ decays exponentially for $`(D2+\eta )/z>1`$ and as a stretched exponential for $`0<(D2+\eta )/z<1`$. Combining this with the result for $`D2+\eta <0`$ outlined in the previous paragraph, gives our general result in Eq. (3) on defining $`\zeta =(D2+\eta )/z`$.

To test the analytical predictions (i) and (ii), we have done preliminary Monte Carlo simulations for $`d=2`$ and $`d=3`$ Ising model for small size lattices with periodic boundary conditions. The results are summarized in Fig. 1. The systems were evolved using heat-bath Monte Carlo dynamics at their bulk critical couplings: $`K_c=[\mathrm{ln}(1+\sqrt{2})]/2`$ in $`d=2`$ and $`K_c0.221656`$ in $`d=3`$. For $`d=2`$, the data for the persistence of the line magnetization for system sizes $`L=63`$, $`95`$ and $`127`$ (in each case the data was averaged over $`1000`$ samples) shows a power law decay $`P(t)t^\theta `$ for $`tL^z`$ and crosses over to a faster decay for $`tL^z`$. A fit to the linear part of the data on the log-log plot gives an estimate $`\theta 0.72`$. Similarly the persistence for the plane magnetization in $`d=3`$ shows a power law decay, $`P(t)t^\theta `$ with $`\theta 0.88`$, estimated from small lattice sizes $`L=15`$ and $`31`$. The estimates of these exponents are only rough and may shift a little with larger lattices. The persistence of the line magnetization in $`d=3`$, in contrast, has a much faster decay. Note that our theory predicts a decay, $`P(t)\mathrm{exp}(a_1t^\zeta )`$ where $`\zeta 0.016`$. Such a small stretching exponent is difficult to determine from the small size lattices due to the strong finite-size effects at late times. All we can say is that the present data for the line magnetization in $`d=3`$ is consistent with a decay of persistence that is faster than a power law but slower than an exponential. More extensive simulations and a somewhat sophisticated finite size scaling analysis are required to pin down the precise value of the stretching exponent.

In summary, we have shown that the persistence of submanifolds in a critical system decays with time in a manner that depends on the dimensions $`d`$ and $`d^{}`$ of the system and manifold respectively. The crossover between the power-law decay observed for the ‘global’ persistence ($`d^{}=d`$) and the exponential decay of the ‘local’ persistence ($`d^{}=0`$), occurs via an intermediate regime of stretched exponential decay with stretched exponent $`\zeta =(dd^{}2+\eta )/z`$, for $`0\zeta 1`$. This latter behavior is predicted for the line magnetization of the $`d=3`$ critical Ising model. Numerical simulations are consistent with the analytical predictions.


text:
# Chapter 0 Dynamical Modeling of SAURON Galaxies

## Chapter 0 Dynamical Modeling of SAURON Galaxies

### Abstract

We describe our program for the dynamical modeling of early-type galaxies observed with the panoramic integral-field spectrograph SAURON. We are using Schwarzschild’s numerical orbit superposition method to reproduce in detail all kinematical and photometric observables, and recover the intrinsic orbital structure of the galaxies. Since catastrophes are the most prominent features in the orbital observables, two-dimensional kinematical coverage is essential to constrain the dynamical models.

### 1 Introduction

We have observed a representative sample of 72 E, S0 and Sa galaxies (de Zeeuw et al. 2002) with the panoramic integral-field spectrograph (IFS) SAURON (Bacon et al. 2001) mounted on the 4.2-m William Herschel telescope. The galaxies were observed to one effective radius, with an effective spatial sampling of 0.8<sup>′′</sup>. For these objects we extracted the stellar kinematics including the $`h_3`$ and $`h_4`$ Gauss–Hermite moments. Gas kinematics and line-strength distribution were measured as well.

We are in the process of constructing dynamical models for the E, S0 galaxies of the sample for which the kinematics and photometry are consistent with axisymmetry. Not to bias our conclusions the models need to be able to reproduce any general orbital distribution. For this we adopted Schwarzschild’s (1979) numerical orbit-superposition method, which is able to fit all kinematical and photometric observations (Rix et al. 1997; van der Marel et al. 1998; Cretton et al. 1999). A similar approach was adopted by other groups to measure the black hole (BH) masses in galaxy nuclei (e.g. Gebhardt et al. 2003).

The implementation of Schwarzschild’s method by van der Marel was adapted (Cappellari et al. 2002) for use with the Multi-Gaussian Expansion parametrization of the galaxy surface brightness (Emsellem et al. 1994; Cappellari 2002) and was applied to IFS data (Verolme et al. 2002). There we showed that IFS observations are important to constrain model parameters such as the galaxy inclination, the stellar mass-to-light ratio and the BH mass. Triaxial Schwarzschild models have also been developed (see Verolme et al. 2003), and will be used to study the SAURON kinematics. The goal is to derive intrinsic shapes, nuclear BH masses and orbital structure, to understand which models are preferred by the galaxies, and to set constraints on galaxy formation scenarios.

### 2 Catastrophes in Orbital Observables

Schwarzschild’s dynamical modeling method consists of finding a positive linear combination of a representative set of orbital building-blocks so as to best fit the kinematical and photometric observables. Here we discuss the nature of orbital observables to better understand the results from the models.

In a three-dimensional time-independent potential, orbits which conserve three integrals of motion are called regular and can be reduced to translations on three-tori in a six-dimensional phase space (e.g. Binney & Tremaine 1987). If there are no resonance conditions between the three frequencies of the motion on the torus, the trajectory of the system covers the torus surface with uniform density after an infinite amount of time. This property can be used to derive accurate time-averaged projected quantities for an orbit, without actually integrating the trajectory for an infinite time. It can be done from the knowledge of the transformation from the torus coordinates to the configuration space coordinates. Although the transformation is generally not known analytically, approximate methods exist for finding it from a short numerical integration of the orbit (e.g. Copin, Zhao, & de Zeeuw 2000).

At an even higher level, one can study the qualitative nature of orbital observables without even performing a single numerical integration of the orbit, using general results from catastrophe theory (e.g. Arnold 1992). When a manifold (e.g. a torus; Fig. 1) is projected onto a two-dimensional (2D) space, its projection is characterized by catastrophes. In 2D space only two stable catastrophes can appear, the fold and the cusp. The fold catastrophe is a curve on the plane, defined as the location where two of the inverse images in the projection merge and disappear. The cusp is a single point, lying at the intersection of folds, where three inverse images coalesce into a single one. If the density is constant on the manifold, the limiting surface density on the projection is $`\mathrm{\Sigma }(𝐱)|𝐱|^{1/2}`$ near the fold catastrophe and is $`\mathrm{\Sigma }(𝐱)|𝐱|^{2/3}`$ near the cusp. These singularities are integrable, as the mass is conserved by projection, and reduce to finite values in realistic situations (e.g. integrated over a pixel or convolved with a PSF), the density being higher at the cusps.

All the general features of the projection of a smooth surface on the plane can be applied to the projection of the uniform-density orbital torus. Fold and in particular cusp catastrophes will be the most prominent features in the projection of any regular orbit on the sky plane. This is precisely what one observes by numerically integrating a regular orbit in an axisymmetric potential, and computing orbital observables in a Monte Carlo fashion. Fig. 2 shows the projected surface brightness and velocity field of a characteristic regular orbit. Its luminosity-weighted projected velocity is significantly different from zero only at the cusp positions, and similar behavior is observed for higher order moments. In an axisymmetric potential regular orbits are the dominant component, which means that catastrophes are important. One can draw the following conclusions, which we will discuss in more detail elsewhere:

* A model composed mainly of regular orbits has a large flexibility in fitting general surface brightness distributions and complex kinematics. The model observables at a given position depend strongly only on the weights assigned to the orbits having a cusp at that position. Each orbit can be used to optimize the fit at the position of its cusps, and complex features (including noise in the data) can be reproduced.
* 2D kinematical coverage is essential to constrain the orbital structure in a galaxy from the observables. Orbital cusps can appear anywhere on the projected image of the galaxy, and the weight of an orbit can be tightly constrained only when its cusp falls within an observed kinematical aperture;

In the case of triaxial galaxy potentials, regular non-resonant triaxial orbits will show the same general features as the axisymmetric case (Fig. 3). Also here 2D kinematical coverage is crucial. However chaos can be significant in triaxial potentials with central singularities (Gerhard & Binney 1985). Chaotic orbits do not lie on manifolds in phase space and their projections are not dominated by catastrophes.

### 3 Dynamical Modeling of NGC 4473

As an example of the models we are constructing, we show in Fig. 4 a data-model comparison for the best-fitting axisymmetric Schwarzschild dynamical model to the SAURON observations of the E5 galaxy NGC 4473. This galaxy was chosen because of its peculiar velocity dispersion field, remaining essentially constant along the galaxy major axis. For an optimal extraction of the kinematics the SAURON datacube was first binned spatially to a constant S/N using the Voronoi 2D-binning method by Cappellari & Copin (2003). Our Schwarzschild modeling software is able to take the irregular shape of the bins precisely and efficiently into account.

As expected from the above considerations, the model is able to reproduce the mean velocity $`V`$, the velocity dispersion $`\sigma `$, and the Gauss–Hermite moments $`h_3`$ and $`h_4`$, well within the errors, over the whole observed field. From the analysis of the model orbital structure we find that two components, of different angular momentum, are needed to fit the observations and explain the observed peculiar kinematics.

The analysis of the orbital structure of the galaxies in the sample we are modeling, will provide, as a function of morphological type and luminosity, information on the intrinsic shapes, the mass distribution and orbital structure, including decoupled components. The knowledge of the DF, in connection with the stellar population analysis we are performing, provides clues to the galaxy formation processes and the link between the formation of the galaxy and the formation of its central BH.

### References


text:
# Semiclassical evaluation of quantum fidelity

## Abstract

We present a numerically feasible semiclassical (SC) method to evaluate quantum fidelity decay (Loschmidt echo, FD) in a classically chaotic system. It was thought that such evaluation would be intractable, but instead we show that a uniform SC expression not only is tractable but it gives remarkably accurate numerical results for the standard map in both the Fermi-golden-rule and Lyapunov regimes. Because it allows Monte Carlo evaluation, the uniform expression is accurate at times when there are $`10^{70}`$ semiclassical contributions. Remarkably, it also explicitly contains the “building blocks” of analytical theories of recent literature, and thus permits a direct test of the approximations made by other authors in these regimes, rather than an a posteriori comparison with numerical results. We explain in more detail the extended validity of the classical perturbation approximation (CPA) and show that within this approximation, the so-called “diagonal approximation” is automatic and does not require ensemble averaging.

uniform semiclassical approximation, fidelity decay, Loschmidt echo

The question of stability of quantum motion, originally formulated by Peres PERES , has recently attracted much interest, due to its relevance to quantum computation and decoherence in complex systems. Peres defined stability in terms of quantum fidelity $`M(t)`$, the overlap at time $`t`$ of two states, which were identical at time $`t=0`$, but afterwards propagated in slightly different dynamical systems, described by Hamiltonians $`H^0`$ and $`H^V=H^0+V`$,

$$M(t)=\left|\psi \left|\mathrm{exp}\left(iH^Vt\right)\mathrm{exp}\left(iH^0t\right)\right|\psi \right|^2.$$
(1)
This quantity is also called Loschmidt echo, because it can be interpreted as an overlap of a state propagated forward for time $`t`$ with $`H^0`$ and then backward for time $`t`$ with $`H^V`$, with the initial state. We consider $`H^0`$ to be strongly chaotic, although our method is not limited to this case. Even with this restriction, the decay of fidelity has a surprisingly rich behavior: Most surprising recently was the derivation in Ref. JALABERT  that for certain range of perturbations the decay rate is independent of the perturbation strength.

The Loschmidt echo is physically realizable, for example in NMR spin echo experiments, where back-propagation under a slightly different Hamiltonian is feasible NMR\_ECHO1 ; NMR\_ECHO2 ; NMR\_ECHO3 . There are other examples, which often go unnoticed. An example is neutron scattering, where the scattering kernel can be written as in Eq. (1), with $`H^V`$ a momentum boosted version of $`H^0`$. Many numerical investigations of FD have been undertaken in various systems JACQUOD ; CERRUTI ; CUCCHIETTI ; CUCCHIETTI2 ; CUCCHIETTI3 ; CUCCHIETTI4 ; uniformCERRUTI ; PROSEN1 ; PROSEN2 ; PROSEN3 ; ZDINARIC ; JACQUOD2 ; JACQUOD3 ; SILVESTROV ; BENENTI1 ; BENENTI2 ; ECKHARDT ; WISNIACKI1 ; WISNIACKI2 ; WISNIACKI3 ; KOTTOS ; EMERSON ; WEINSTEIN ; ADAMOV ; WANG ; GEORGEOT ; VANICEK . Depending on the strength of perturbation, there exist at least four qualitatively different regimes of the decay in chaotic systems CUCCHIETTI : As the perturbation increases, these regimes are perturbative (PT), Fermi-golden-rule (FGR), Lyapunov (L), and the strong SC regime.

In the PT regime, in which the characteristic matrix element of the perturbation is smaller than the mean level spacing $`\mathrm{\Delta }`$, the decay can be described by a combination of perturbation theory and random-matrix theory (RMT), and is Gaussian CERRUTI ; CUCCHIETTI ,

$$M_{PT}\left(t\right)\mathrm{exp}\left(\overline{V^2}t^2/\mathrm{}^2\right).$$
(2)
For intermediate perturbation strengths, the decay follows the Fermi golden rule JACQUOD  and is exponential,

$$M_{FGR}\left(t\right)\mathrm{exp}\left(\mathrm{\Gamma }t/\mathrm{}\right)$$
(3)
where $`\mathrm{\Gamma }=2\pi \overline{V^2}/\mathrm{\Delta }`$. In Ref. CUCCHIETTI  it is shown that this FGR decay is equivalent to the exponential decay derived semiclassically in Refs. JALABERT ; CERRUTI . In other words, $`\mathrm{\Gamma }=2K/\mathrm{}`$ where $`K`$ is the classical action diffusion constant,

$$K=_0^{\mathrm{}}𝑑tV\left[𝐫\left(t\right)\right]V\left[𝐫\left(0\right)\right].$$

In the Lyapunov regime, derived in Ref. JALABERT , FD actually does not depend on the strength of perturbation, but only on the Lyapunov exponent $`\lambda `$ of the chaotic system,

$$M_L\left(t\right)\mathrm{exp}\left(\lambda t\right).$$
(4)
We are able to find a numerically feasible uniform VANICEK ; VANICEK1 ; VANICEK2 ; VANICEK3  SC method to evaluate FD in the FGR and Lyapunov regimes. As a result, we can directly test all approximations made in the derivation of results (3) and (4) from Refs. JALABERT ; CERRUTI . The method starts with a SC approach based on the CPA JALABERT ; CERRUTI , and ends with a form of initial value representation (IVR) ivr ; IVR2  which makes the numerical calculation manageable and the SC approximation itself more accurate.

Following notation of Ref. JALABERT , we want to find FD for an initial Gaussian wave packet

$$\psi (𝐫;0)=\left(\pi \sigma ^2\right)^{d/4}\mathrm{exp}\left[\frac{i}{\mathrm{}}𝐩_0\left(𝐫𝐫_0\right)\frac{\left(𝐫𝐫_0\right)^2}{2\sigma ^2}\right].$$

It is centered at $`𝐫_0`$ with dispersion $`\sigma `$ and has an average momentum $`𝐩_0.`$ We propagate this state with a SC Van Vleck-Gutzwiller propagator vanVleck

$$K^{sc}(𝐫^{\prime \prime },𝐫^{};t)=\underset{j}{}\left(2\pi i\mathrm{}\right)^{d/2}C_j^{1/2}\mathrm{exp}\left(\frac{i}{\mathrm{}}S_ji\frac{\pi }{2}\nu _j\right).$$

Here $`C_j=\left|det\left(^2S_j/𝐫^{\prime \prime }𝐫^{}\right)\right|`$ is the absolute value of the Van Vleck determinant, $`S_j(𝐫^{\prime \prime },𝐫^{};t)`$ is the action along the $`j`$th trajectory connecting $`𝐫^{}`$ with $`𝐫^{\prime \prime }`$,

$$S_j(𝐫^{\prime \prime },𝐫^{};t)=_0^t𝑑t^{}L[𝐫\left(t^{}\right),\dot{𝐫}\left(t^{}\right),t^{}]$$

and $`\nu _j`$ is the Maslov index.

Expanding each contribution about a central trajectory TOMSOVIC , the overlap amplitude of the semiclassically propagated states becomes JALABERT

$`O(t)=\psi ^{sc,V}\left(t\right)|\psi ^{sc}\left(t\right)=\left(\sigma ^2/\pi \mathrm{}^2\right)^{d/2}{\displaystyle d^dr}`$ (5)
$`\times {\displaystyle \underset{j,j^{}}{}}\left(C_j^VC_j^{}\right)^{1/2}\mathrm{exp}[{\displaystyle \frac{i}{\mathrm{}}}(S_j^VS_j^{}){\displaystyle \frac{i\pi }{2}}(\nu _j^V\nu _j^{})]`$
$`\times \mathrm{exp}\left\{\left[\left(𝐩_j^{}𝐩_0\right)^2+\left(𝐩_j^{}^{}𝐩_0\right)^2\right]\sigma ^2/2\mathrm{}^2\right\}`$
where $`S_j=S_j(𝐫,𝐫_0;t)`$ and superscript $`V`$ denotes quantities in the perturbed system. At this point, two crucial approximations are made in Refs. JALABERT ; CERRUTI : First, only the diagonal terms $`j=j^{}`$ are considered. Ref. JALABERT  claims that these are the only terms surviving the average over impurities in disordered systems. Below we show that this is not a separate approximation, but that it follows from the CPA and does not require any ensemble averaging. CPA, the second approximation used in Refs. JALABERT ; CERRUTI  is based on an apparently hopeless assumption that the perturbation does not affect trajectories (i.e., $`C_j^VC_j`$ and $`\nu _j^V\nu _j`$) but only affects the actions, through

$$\mathrm{\Delta }S_j=S_j^VS_j=_0^t𝑑t^{}V\left[𝐫_j\left(t^{}\right)\right].$$
(6)
Of course this assumption is wrong for individual trajectories which deviate exponentially with time. The reason why the approximation works in quantum mechanics is subtle: The first step to understanding why it yields accurate wave functions lies in the structural stability of the manifolds, as pointed out in Ref. CERRUTI . Assuming that perturbation does not cause a bifurcation and does not significantly change the stable manifold, the evolved manifolds almost exactly overlap whereas the same initial points deviate exponentially by sliding along the manifold CERRUTI .

The second step goes as follows: consider trajectories $`A(t)`$, $`A^V(t)`$ under the flow $`H^0`$, $`H^V`$, respectively. Let $`A(0)=A^V(0)`$ be a point on the Lagrangian manifold supporting the wave function at $`t=0`$. While $`A^V(t)`$ exponentially diverges from $`A(t)`$, if the evolved manifolds (almost) exactly overlap, we can find a point $`B(0)`$ on the manifold at $`t=0`$ such that $`B^V(t)`$ (almost) coincides with $`A(t)`$. Because of the exponential sensitivity to the initial conditions, point $`B(0)`$ will be exponentially close to $`A(0)`$. Trajectories $`A(t)`$ and $`B(t)`$ remain exponentially close for all times, so if we use these particular trajectories to find $`\psi (t)`$ and $`\psi ^V(t)`$, respectively, the CPA will be justified.

The diagonal approximation and CPA enormously simplify expression (5) for the overlap amplitude:

$`O(t)=\left(\sigma ^2/\pi \mathrm{}^2\right)^{d/2}`$ (7)
$`\times {\displaystyle }d^dr{\displaystyle \underset{j}{}}C_j\mathrm{exp}[i\mathrm{\Delta }S_j/\mathrm{}(𝐩_j^{}𝐩_0)^2\sigma ^2/\mathrm{}^2].`$
At this point, both Refs. JALABERT ; CERRUTI  resort to statistical arguments to obtain an analytical result. Expression (7) for the overlap would be very difficult to implement numerically for three reasons: First, in chaotic systems there is an exponentially growing number of contributing trajectories. Second, the accuracy would be compromised by proliferating caustic singularities in the Van Vleck determinant $`C_j`$ whenever $`𝐫/𝐩_j^{}=0.`$ Finally, for each trajectory we would have to perform a computationally expensive root-search to find initial $`𝐩_j^{}`$ that satisfies $`𝐫(𝐫_0,𝐩_j^{},t)=𝐫`$. However, there exists a beautiful and simple way to eliminate the exponential number of contributions, caustic singularities, and the root search, all at the same time. All three problems can be solved if we evaluate the overlap (5) in the initial momentum instead of final position representation. Exactly one point on the evolved manifold corresponds to each initial momentum, so no summation is necessary. The new “Van Vleck determinant” is exactly 1, so there will be no Maslov indices either. With all these simplifications, the SC evaluation becomes tractable; in principle, it yields the same result that an arduous evaluation of (7) would:

$`O(t)=\left(\sigma ^2/\pi \mathrm{}^2\right)^{d/2}{\displaystyle d^dp^{}}`$ (8)
$`\times \mathrm{exp}\left[i\mathrm{\Delta }S[𝐫(𝐫_0,𝐩^{},t),𝐫_0,t]/\mathrm{}\left(𝐩^{}𝐩_0\right)^2\sigma ^2/\mathrm{}^2\right].`$
The only assumption required to derive (8) is the validity of CPA, in the extended sense described above. Ensemble averaging used in Ref. JALABERT  is unnecessary: result (8) works for pure states. Expression (8) is a special form of IVR ivr ; IVR2 . In general, IVR avoids the singularities and the root search, but at a cost of replacing a sum over classically allowed paths by an integral over all initial momenta. In our case, it is even better, since we also eliminated the integral over final position $`𝐫`$. We remark that (8) can also be obtained by changing the integration variable in (7) from final $`𝐫`$ to initial $`𝐩^{}`$, but our derivation avoids the intermediate step (7) that requires making diagonal approximation in (5). We note the unique property of IVR: in this representation, FD is only due to dephasing. In other representations, the decay can also have a component due to the decay of classical overlaps.

We chose to test our method on the standard map used in Ref. CERRUTI ,

$`q_{j+1}`$ $`=q_j+p_j(\mathrm{mod}1),`$
$`p_{j+1}`$ $`=p_j{\displaystyle \frac{k}{2\pi }}\mathrm{sin}\left(2\pi q_{j+1}\right)(\mathrm{mod}1).`$
Perturbation is effected by replacing the parameter $`k`$ by $`k+ϵ`$. Choice of an $`n`$-dimensional Hilbert space for the quantized map fixes the effective Planck constant to be $`\mathrm{}=\left(2\pi n\right)^1`$. We note that results of exact quantum and SC computations which we present below are for initial position eigenstate with $`q_0=0.5`$ rather than a wavepacket.

In previous numerical experiments analytical predictions of Gaussian or exponential decay have been compared to an exact quantum calculation: see, e.g., Refs. JACQUOD ; CERRUTI ; CUCCHIETTI ; BENENTI1 ; SILVESTROV . While we also have an exact quantum benchmark (FFT) with which to compare the expressions for various regimes, we reiterate that it would be hard from a mere comparison of final results for $`M(t)`$ to determine the source of errors. We proceed by discussing how the uniform method helps to analyze various regimes of decay. In the PT regime (see Fig. 1), we do not expect any SC approach to work very well except for short times (much shorter than the Heisenberg time $`t_H=h/\mathrm{\Delta }`$). The RMT analytical result $`M_{PT}`$ from Ref. CERRUTI  gives an excellent agreement in this case. The inset shows, however, that before the Gaussian decay $`M_{PT}`$ sets in at the Heisenberg time, the uniform expression follows $`M_{exact}`$ much better.

As the perturbation $`ϵ`$ increases, we enter the regimes with exponential decay of fidelity. If the perturbation is strong quantum mechanically, but does not significantly change the stable manifold, CPA may be used. Even within the CPA, there are two types of decay, discussed already in Ref. JALABERT . First, there is decay related to dephasing of trajectories with uncorrelated actions. Second, there is decay related to dephasing of very near trajectories with correlated actions. For smaller perturbations, the first type of decay is slower and dominates the behavior of fidelity: this happens in the FGR regime. For larger perturbations, dephasing of uncorrelated trajectories is so fast that the quantum overlap is determined by the fraction of near trajectories that have remained in phase. This is the case in the Lyapunov regime. Transition from the PT to FGR regime occurs for $`ϵ^232\pi ^2n^3[1+2J_2(k)]^1`$ CERRUTI  when most of the overlap has decayed before Heisenberg time. Transition from the FGR to Lyapunov regime occurs for $`ϵ^28\pi ^2\lambda n^2[1+2J_2(k)]^1`$ when the FGR decay rate is larger than $`\lambda `$.

Using Eq. (8), fidelity can be written as a weighted average of terms $`\mathrm{exp}\left[i(\mathrm{\Delta }S^{}\mathrm{\Delta }S^{\prime \prime })/\mathrm{}\right]`$,

$`M_{unif}\left(t\right)=\left({\displaystyle \frac{\sigma ^2}{\pi \mathrm{}^2}}\right)^d{\displaystyle d^dp^{}d^dp^{\prime \prime }\mathrm{exp}\left[\frac{i}{\mathrm{}}(\mathrm{\Delta }S^{}\mathrm{\Delta }S^{\prime \prime })\right]}`$
$`\times \mathrm{exp}\left\{\left[\left(𝐩^{}𝐩_0\right)^2+\left(𝐩^{\prime \prime }𝐩_0\right)^2\right]\sigma ^2/\mathrm{}^2\right\}`$ (9)
where $`\mathrm{\Delta }S^{\prime \prime }`$ corresponds to a trajectory with initial momentum $`𝐩^{\prime \prime }`$. Assuming the averaging window (i.e. the momentum width of the wavepacket) is large enough, we can make the replacement

$$\mathrm{exp}\left[i(\mathrm{\Delta }S^{}\mathrm{\Delta }S^{\prime \prime })/\mathrm{}\right]\mathrm{exp}\left[i(\mathrm{\Delta }S^{}\mathrm{\Delta }S^{\prime \prime })/\mathrm{}\right]$$
(10)
in Eq. (9) where averaging is over all initial momenta $`p^{}`$, $`p^{\prime \prime }`$. In the FGR regime where dephasing is determined by uncorrelated trajectories, a further simplification

$$\mathrm{exp}\left[i(\mathrm{\Delta }S^{}\mathrm{\Delta }S^{\prime \prime })/\mathrm{}\right]e^{i\mathrm{\Delta }S^{}/\mathrm{}}e^{i\mathrm{\Delta }S^{\prime \prime }/\mathrm{}}$$
(11)
is possible. Due to the central limit theorem, in chaotic systems distribution of $`\mathrm{\Delta }S`$ approaches a Gaussian and

$$\mathrm{exp}\left(i\mathrm{\Delta }S/\mathrm{}\right)=\mathrm{exp}\left[i\mathrm{\Delta }S/\mathrm{}\sigma _{\mathrm{\Delta }S}^2/2\mathrm{}^2\right]$$
(12)
where $`\sigma _{\mathrm{\Delta }S}^2=2Kt`$ is the action variance at time $`t`$. Applying approximations (10), (11), (12) in Eq. (9), confirms Eq. (3) for the FGR decay JALABERT ; CERRUTI . Fig. 2 shows FD in the FGR regime. In the inset, histogram of action differences is compared with a Gaussian fit, confirming assumption (12). It is apparent that $`M_{unif}`$ matches $`M_{exact}`$ better than the $`M_{FGR}`$ since $`M_{unif}`$ takes into account the precise initial conditions without the averaging assumption (10) and since $`M_{FGR}`$ uses an analytic result for $`K`$, which is only approximate CERRUTI . Careful inspection of the short time regime (not shown) reveals that $`M_{unif}`$ agrees with $`M_{exact}`$, since unlike $`M_{FGR}`$, $`M_{unif}`$ does not depend on the central limit theorem which guarantees the Gaussian assumption (12) at later times. Finally, we would like to point out that the uniform expression is very accurate at time $`t120`$ when there are approximately $`10^{70}`$ semiclassical contributions in the sum (7)!

In the Lyapunov regime, FD is determined by dephasing of near trajectories with correlated actions JALABERT , invalidating simplification (11). Now the action difference $`\mathrm{\Delta }S^{}\mathrm{\Delta }S^{\prime \prime }`$ depends on the initial momenta $`p^{}`$, $`p^{\prime \prime }`$. Using reasoning similar to Ref. JALABERT  or statistical arguments for a random walk with an exponentially increasing time step RANDOM\_WALK , it can be shown that the action difference is also Gaussian distributed, with zero average and variance

$`\left[\mathrm{\Delta }S\left(p^{}\right)\mathrm{\Delta }S\left(p^{\prime \prime }\right)\right]^2`$ $`(D/2\lambda )e^{2\lambda t}\left(p^{}p^{\prime \prime }\right)^2,`$ (13)
$`D`$ $`=2{\displaystyle _0^{\mathrm{}}}𝑑tV^{}\left(0\right)V^{}\left(t\right).`$
We can therefore make the replacement

$$\mathrm{exp}\left[\frac{i}{\mathrm{}}(\mathrm{\Delta }S^{}\mathrm{\Delta }S^{\prime \prime })\right]\mathrm{exp}\left[\frac{D}{4\lambda \mathrm{}^2}e^{2\lambda t}\left(p^{}p^{\prime \prime }\right)^2\right]$$

in Eq. (9) and (10) to find

$$M_L\left(t\right)\left(1+e^{2\lambda t}D/2\lambda \sigma ^2\right)^{1/2}\left(2\lambda \sigma ^2/D\right)^{1/2}e^{\lambda t},$$

confirming Eq. (4). For the precise definition of $`\lambda `$, see Ref. SILVESTROV (one has to be careful about the averaging process). Fig. 3 displays $`M\left(t\right)`$ in the Lyapunov regime. It shows that while $`M_L`$ gives an accurate average decay only for $`\lambda t1`$, $`M_{unif}`$ correctly follows the behavior of $`M_{exact}`$ even for short times $`t\lambda ^1`$. The inset shows the variance of $`\mathrm{\Delta }S\left(p^{}\right)\mathrm{\Delta }S\left(p^{\prime \prime }\right)`$ as a function of $`p^{}p^{\prime \prime }`$ at a fixed time and justifies the assumption made in Ref. JALABERT  in derivation of perturbation independent decay: for near trajectories, the variance grows quadratically with $`p^{}p^{\prime \prime }`$ (fitted line gives an exponent 2.003), in accordance with Eq. (13), while for distant trajectories, in accordance with the derivation of the FGR regime, the variance is independent of $`p^{}p^{\prime \prime }`$,

$$\left[\mathrm{\Delta }S\left(p^{}\right)\mathrm{\Delta }S\left(p^{\prime \prime }\right)\right]^2=2\sigma _{\mathrm{\Delta }S}^2=4Kt.$$
(14)
The time dependence of $`\left[\mathrm{\Delta }S\left(p^{}\right)\mathrm{\Delta }S\left(p^{\prime \prime }\right)\right]^2`$ for fixed $`p^{}p^{\prime \prime }`$ is shown in Fig. 4. Part a) shows that for short times when trajectories are still correlated, this dependence is exponential, in agreement with Eq. (13). Part b) shows that for longer times, when correlation is lost, the dependence is linear, as expected from Eq. (14).

To conclude, we have explicitly evaluated SC expressions which were thought to be intractable numerically, yielding remarkably accurate results for FD in the FGR and Lyapunov regimes. We provided a more detailed explanation why CPA works and employed our method to test other approximations used in Refs. JALABERT ; CERRUTI .

This research was supported by the National Science Foundation under Grant No. NSF-CHE-0073544, by ITAMP at the Harvard-Smithsonian Center for Astrophysics, Harvard University, and by the Mathematical Sciences Research Institute at Berkeley. We would like to acknowledge helpful discussions with S. Tomsovic. E.J.H. acknowledges the hospitality of the Max Planck Institute for Complex Systems, Dresden and the Humboldt Foundation for hospitality and support.