ar
X
iv
:q
ua
nt
-p
h/
03
02
19
2v
2
1
7
Se
p
20
03
Semiclassical evaluation of quantum fidelity
Jiˇr´ı Van´ıcˇek1, 2 and Eric J. Heller1, 3
1Department of Physics, Harvard University, Cambridge, Massachusetts 02138
2Mathematical Sciences Research Institute, Berkeley, California 94720
3Department of Chemistry, Harvard University, Cambridge, Massachusetts 02138
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 Lya-
punov regimes. Because it allows Monte Carlo evaluation, the uniform expression is accurate at
times when there are 1070 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 approxima-
tion” is automatic and does not require ensemble averaging.
PACS numbers: 05.45.Mt, 03.65.Sq, 03.65.Yz
The question of stability of quantum motion, originally
formulated by Peres [1], has recently attracted much in-
terest, 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 sys-
tems, described by Hamiltonians H0 and HV = H0 +V ,
M(t) =
∣
∣
〈
ψ
∣
∣exp
(
iHV t
)
exp
(
−iH0t
)∣
∣ψ
〉∣
∣
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 H0 and then backward for time
t with HV , with the initial state. We consider H0 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. [2] that for certain
range of perturbations the decay rate is independent of
the perturbation strength.
The Loschmidt echo is physically realizable, for ex-
ample in NMR spin echo experiments, where back-
propagation under a slightly different Hamiltonian is fea-
sible [3, 4, 5]. 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 HV a
momentum boosted version of H0. Many numerical in-
vestigations of FD have been undertaken in various sys-
tems [6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20,
21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32]. Depending
on the strength of perturbation, there exist at least four
qualitatively different regimes of the decay in chaotic sys-
tems [6]: 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 ∆, the decay can be described by a combi-
nation of perturbation theory and random-matrix theory
(RMT), and is Gaussian [6, 7],
MPT (t) ≈ exp
(
−V 2 t2/~2
)
. (2)
For intermediate perturbation strengths, the decay fol-
lows the Fermi golden rule [8] and is exponential,
MFGR (t) ≈ exp (−Γt/~) (3)
where Γ = 2πV 2/∆. In Ref. [6] it is shown that this
FGR decay is equivalent to the exponential decay derived
semiclassically in Refs. [2, 7]. In other words, Γ = 2K/~
where K is the classical action diffusion constant,
K =
∫ ∞
0
dt 〈V [r (t)]V [r (0)]〉 .
In the Lyapunov regime, derived in Ref. [2], FD actu-
ally does not depend on the strength of perturbation, but
only on the Lyapunov exponent λ of the chaotic system,
ML (t) ∼ exp (−λt) . (4)
We are able to find a numerically feasible uniform
[9, 33, 34, 35] 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. [2, 7]. The method starts with a
SC approach based on the CPA [2, 7], and ends with a
form of initial value representation (IVR) [36, 37] which
makes the numerical calculation manageable and the SC
approximation itself more accurate.
Following notation of Ref. [2], we want to find FD for
an initial Gaussian wave packet
ψ (r; 0) =
(
πσ2
)−d/4
exp
[
i
~p0 · (r − r0)−
(r− r0)2
2σ2
]
.

2It is centered at r0 with dispersion σ and has an average
momentum p0. We propagate this state with a SC Van
Vleck-Gutzwiller propagator [38]
Ksc(r′′, r′; t) =
∑
j
(2πi~)−d/2C1/2j exp
( i
~Sj − i
π
2
νj
)
.
Here Cj =
∣
∣det
(
∂2Sj/∂r′′∂r′
)∣
∣ is the absolute value of
the Van Vleck determinant, Sj (r′′, r′; t) is the action
along the jth trajectory connecting r′ with r′′,
Sj (r′′, r′; t) =
∫ t
0
dt′L [r (t′) , r˙ (t′) ,t′]
and νj is the Maslov index.
Expanding each contribution about a central trajec-
tory [39], the overlap amplitude of the semiclassically
propagated states becomes [2]
O(t) =
〈
ψsc,V (t) |ψsc (t)
〉
=
(
σ2/π~2
)d/2
∫
ddr (5)
×
∑
j,j′
(
CVj Cj′
)1/2
exp
[ i
~
(
SVj − Sj′
)
− iπ
2
(
νVj − νj′
)
]
× exp
{
−
[
(
p′j − p0
)2
+
(
p′j′ − p0
)2
]
σ2/2~2
}
where Sj = Sj (r, r0; t) and superscript V denotes quan-
tities in the perturbed system. At this point, two crucial
approximations are made in Refs. [2, 7]: First, only the
diagonal terms j = j′ are considered. Ref. [2] 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 averag-
ing. CPA, the second approximation used in Refs. [2, 7]
is based on an apparently hopeless assumption that the
perturbation does not affect trajectories (i.e., CVj ≈ Cj
and νVj ≈ νj) but only affects the actions, through
∆Sj = SVj − Sj = −
∫ t
0
dt′ V [rj (t′)] . (6)
Of course this assumption is wrong for individual tra-
jectories which deviate exponentially with time. The rea-
son 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. [7]. 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 mani-
fold [7].
The second step goes as follows: consider trajectories
A(t), AV (t) under the flow H0, HV , respectively. Let
A(0) = AV (0) be a point on the Lagrangian manifold
supporting the wave function at t = 0. While AV (t)
exponentially diverges fromA(t), if the evolved manifolds
(almost) exactly overlap, we can find a point B(0) on the
manifold at t = 0 such thatBV (t) (almost) coincides with
A(t). Because of the exponential sensitivity to the initial
conditions, pointB(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
ψ(t) and ψV (t), respectively, the CPA will be justified.
The diagonal approximation and CPA enormously sim-
plify expression (5) for the overlap amplitude:
O(t) =
(
σ2/π~2
)d/2
(7)
×
∫
ddr
∑
j
Cj exp
[
i∆Sj/~ −
(
p′j − p0
)2 σ2/~2
]
.
At this point, both Refs. [2, 7] resort to statistical ar-
guments to obtain an analytical result. Expression (7)
for the overlap would be very difficult to implement nu-
merically for three reasons: First, in chaotic systems
there is an exponentially growing number of contribut-
ing trajectories. Second, the accuracy would be com-
promised by proliferating caustic singularities in the Van
Vleck determinant Cj whenever ∂r/∂p′j = 0. Finally, for
each trajectory we would have to perform a computation-
ally expensive root-search to find initial p′j that satisfies
r
(
r0,p′j, t
)
= r. However, there exists a beautiful and
simple way to eliminate the exponential number of con-
tributions, 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 de-
terminant” 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) =
(
σ2/π~2
)d/2
∫
ddp′ (8)
× exp
[
i∆S [r (r0,p′, t) , r0, t] /~ − (p′ − p0)2 σ2/~2
]
.
The only assumption required to derive (8) is the validity
of CPA, in the extended sense described above. Ensem-
ble averaging used in Ref. [2] is unnecessary: result (8)
works for pure states. Expression (8) is a special form
of IVR [36, 37]. 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 ini-
tial momenta. In our case, it is even better, since we also
eliminated the integral over final position r. We remark
that (8) can also be obtained by changing the integration
variable in (7) from final r to initial p′, but our deriva-
tion avoids the intermediate step (7) that requires mak-
ing 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.

3We chose to test our method on the standard map used
in Ref. [7],
qj+1 = qj + pj (mod 1),
pj+1 = pj −
k
2π sin (2πqj+1) (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
~ = (2πn)−1. We note that results of exact quantum
and SC computations which we present below are for
initial position eigenstate with q0 = 0.5 rather than a
wavepacket.
In previous numerical experiments analytical predic-
tions of Gaussian or exponential decay have been com-
pared to an exact quantum calculation: see, e.g., Refs.
[6, 7, 8, 10, 11]. While we also have an exact quantum
benchmark (FFT) with which to compare the expres-
sions 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 tH = h/∆). The
RMT analytical result MPT from Ref. [7] gives an ex-
cellent agreement in this case. The inset shows, how-
ever, that before the Gaussian decay MPT sets in at the
Heisenberg time, the uniform expression follows Mexact
much better.
0 1000 2000 3000 4000
1
0.05
0.2
0.1
0.5
FGR
exact
PT
uniform
0 25 50
0.999
1
M
t
M
t
FIG. 1: Fidelity in the perturbative regime (k = 18, λ ≈ 2.21,
ǫ = 10−4, tH ≈ n = 350). Inset: detail for short times.
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. [2]. First, there is decay related to de-
phasing of trajectories with uncorrelated actions. Sec-
ond, there is decay related to dephasing of very near
trajectories with correlated actions. For smaller per-
turbations, the first type of decay is slower and domi-
nates the behavior of fidelity: this happens in the FGR
regime. For larger perturbations, dephasing of uncor-
related 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 Lya-
punov regime. Transition from the PT to FGR regime
occurs for ǫ2 ≈ 32π2n−3[1 + 2J2(k)]−1 [7] when most of
the overlap has decayed before Heisenberg time. Tran-
sition from the FGR to Lyapunov regime occurs for
ǫ2 ≈ 8π2λn−2[1+ 2J2(k)]−1 when the FGR decay rate is
larger than λ.
Using Eq. (8), fidelity can be written as a weighted
average of terms exp [i(∆S′−∆S′′)/~],
Munif (t) =
( σ2
π~2
)d ∫
ddp′
∫
ddp′′ exp
[ i
~ (∆S
′−∆S′′)
]
× exp
{
−
[
(p′ − p0)2 + (p′′ − p0)2
]
σ2/~2
}
(9)
where ∆S′′ corresponds to a trajectory with initial mo-
mentum p′′. Assuming the averaging window (i.e. the
momentum width of the wavepacket) is large enough, we
can make the replacement
exp [i(∆S′−∆S′′)/~] ≈ 〈exp [i(∆S′−∆S′′)/~]〉 (10)
in Eq. (9) where averaging is over all initial momenta p′,
p′′. In the FGR regime where dephasing is determined
by uncorrelated trajectories, a further simplification
〈exp [i(∆S′−∆S′′)/~]〉 ≈
〈
ei∆S′/~
〉〈
e−i∆S′′/~
〉
(11)
is possible. Due to the central limit theorem, in chaotic
systems distribution of ∆S approaches a Gaussian and
〈exp (i∆S/~)〉 = exp
[
i 〈∆S〉 /~ − σ2∆S/2~2
]
(12)
where σ2∆S = 2Kt is the action variance at time t. Apply-
ing approximations (10), (11), (12) in Eq. (9), confirms
Eq. (3) for the FGR decay [2, 7]. 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 Munif matches Mexact better
than the MFGR since Munif takes into account the pre-
cise initial conditions without the averaging assumption
(10) and sinceMFGR uses an analytic result forK, which
is only approximate [7]. Careful inspection of the short
time regime (not shown) reveals that Munif agrees with
Mexact, since unlike MFGR, Munif 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 accu-
rate at time t ≈ 120 when there are approximately 1070
semiclassical contributions in the sum (7)!
In the Lyapunov regime, FD is determined by dephas-
ing of near trajectories with correlated actions [2], in-
validating simplification (11). Now the action difference
∆S′−∆S′′ depends on the initial momenta p′, p′′. Using
reasoning similar to Ref. [2] or statistical arguments for a
random walk with an exponentially increasing time step

4P(

)S
50 100 150 200
10-4
10-3
10-2
10-1
exact
FGR
uniform
PT
Lyap.
ergodic
M
t
∆S
∆
-5 0 5
0
0.1
FIG. 2: Fidelity in the FGR regime (k = 18, λ ≈ 2.21,
ǫ = 5× 10−4, n = 3500). Horizontal dashed line (“ergodic”)
is the limit of FD due to the finite size of Hilbert space. Inset:
Histogram of action differences compared to a Gaussian fit.
[40], it can be shown that the action difference is also
Gaussian distributed, with zero average and variance
〈
[∆S (p′)−∆S (p′′)]2
〉
≈ (D/2λ)e2λt (p′ − p′′)2 , (13)
D = 2
∫ ∞
0
dt 〈V ′ (0)V ′ (t)〉 .
We can therefore make the replacement
〈
exp
[ i
~(∆S
′−∆S′′)
]〉
≈ exp
[
− D
4λ~2 e
2λt (p′ − p′′)2
]
in Eq. (9) and (10) to find
ML (t) ≈
(
1 + e2λtD/2λσ2
)−1/2 ≈
(
2λσ2/D
)1/2 e−λt,
confirming Eq. (4). For the precise definition of λ, see
Ref. [10](one has to be careful about the averaging pro-
cess). Fig. 3 displays M (t) in the Lyapunov regime. It
shows that while ML gives an accurate average decay
only for λt ≫ 1, Munif correctly follows the behavior of
Mexact even for short times t ∼ λ−1. The inset shows
the variance of ∆S (p′)−∆S (p′′) as a function of p′− p′′
at a fixed time and justifies the assumption made in Ref.
[2] in derivation of perturbation independent decay: for
near trajectories, the variance grows quadratically with
p′ − p′′ (fitted line gives an exponent 2.003), in accor-
dance with Eq. (13), while for distant trajectories, in
accordance with the derivation of the FGR regime, the
variance is independent of p′ − p′′,
〈
[∆S (p′)−∆S (p′′)]2
〉
= 2σ2∆S = 4Kt. (14)
The time dependence of
〈
[∆S (p′)−∆S (p′′)]2
〉
for fixed
p′ − p′′ is shown in Fig. 4. Part a) shows that for short
times when trajectories are still correlated, this depen-
dence 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).
0 5 10 15 20
10-6
10-4
10-2
100
10-15 10-10 10-5 100
10-20
10-10
100
t
M
V
ar
ia
nc
e
p’’−p’
FIG. 3: Fidelity in the Lyapunov regime (k = 7, λ ≈ 1.28,
ǫ = 5 × 10−4, n = 105). Meaning of lines same as in Fig. 2.
Inset: Variance of ∆S (p′′) −∆S (p′) as a function of p′′ − p′
at time t = 7. Dots are numerically calculated, dashed line
is the horizontal asymptote 2σ2∆S, solid line is a linear fit for
small p′′ − p′′, in agreement with Eq. (13).
0 25 50 75 100
0
1000
2000
3000
0 10 20 30
10-20
10-10
1
(
)
S’
’−
∆
S’
2
∆
t t
b.a.
FIG. 4: Variance of ∆S (p′′)−∆S (p′) as a function of t for
p′′ − p′ = 10−11: a) exponential dependence for short times,
b) linear dependence for long times.
To conclude, we have explicitly evaluated SC expres-
sions which were thought to be intractable numerically,
yielding remarkably accurate results for FD in the FGR
and Lyapunov regimes. We provided a more detailed ex-
planation why CPA works and employed our method to
test other approximations used in Refs. [2, 7].
This research was supported by the National Science
Foundation under Grant No. NSF-CHE-0073544, by
ITAMP at the Harvard-Smithsonian Center for Astro-
physics, 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 Foun-
dation for hospitality and support.

5[1] A. Peres, Phys. Rev. A 30, 1610 (1984).
[2] R.A. Jalabert and H.M. Pastawski, Phys. Rev. Lett. 86,
2490 (2001).
[3] W.K. Rhim et al., Phys. Rev. Lett. 25, 218 (1970).
[4] H.M. Pastawski et al., Phys. Rev. Lett. 75, 4310 (1995).
[5] H.M. Pastawski et al., Physica 283A, 166 (2000).
[6] F.M. Cucchietti et al., Phys. Rev. E 65, 046209 (2002).
[7] N.R. Cerruti and S. Tomsovic, Phys. Rev. Lett. 88,
054103 (2002).
[8] P. Jacquod, P.G. Silvestrov, and C.W.J. Beenakker,
Phys. Rev. E 64, 055203 (2001).
[9] J. Van´ıcˇek, Uniform semiclassical approximations and
their applications, doctoral thesis, Harvard University,
Cambridge, MA (2003).
[10] P.G. Silvestrov, J. Tworzydlo, and C.W.J. Beenakker,
Phys. Rev. E 67, 025204(R) (2003).
[11] G. Benenti and G. Casati, Phys. Rev. E 65, 066205
(2002).
[12] F.M. Cucchietti, H.M. Pastawski, and D.A. Wisniacki,
Phys. Rev. E 65, 045206(R)(2002).
[13] F.M. Cucchietti, D.A.R. Dalvit, J.P. Paz, and W.H.
Zurek, quant-ph/0306142.
[14] F.M. Cucchietti, H.M. Pastawski, and R.A. Jalabert,
cond-mat/0307752.
[15] T. Prosen, T.H. Seligman, and M. Zˇdinaricˇ,
quant-ph/0304104.
[16] T. Prosen and M. Zˇdinaricˇ, quant-ph/0306097.
[17] M. Zˇdinaricˇ and T. Prosen, J. Phys. A 36, 2463 (2003).
[18] T. Prosen, Phys. Rev. E 65, 036208 (2002).
[19] P. Jacquod, I. Adagideli, and C.W.J. Beenakker, Phys.
Rev. Lett. 89, 154103 (2002).
[20] P. Jacquod, I. Adagideli, and C.W.J. Beenakker, Euro-
phys. Lett. 61, 729 (2003).
[21] G. Benenti, G. Casati, and G. Veble, Phys. Rev. E 67,
055202(R) (2003).
[22] B. Eckhardt, J. Phys. A 36, 371 (2003).
[23] D.A. Wisniacki, E.G. Vergini, H.M. Pastawski, and F.M.
Cucchietti, Phys. Rev. E 65, 055206(R) (2002).
[24] D.A. Wisniacki and D. Cohen, Phys. Rev. E 66, 046209
(2002).
[25] D.A. Wisniacki, Phys. Rev. E 67, 016205 (2003).
[26] T. Kottos and D. Cohen, Europhys. Lett. 61, 431 (2003).
[27] J. Emerson, Y.S. Weinstein, S.Lloyd, and D.G. Cory,
Phys. Rev. Lett. 89, 284102 (2002).
[28] Y.S. Weinstein, J. Emerson, S.Lloyd, and D.G. Cory,
quant-ph/0201064 (2002).
[29] Y. Adamov, I. V. Gornyi, and A.D. Mirlin, Phys. Rev. E
67, 056217 (2003).
[30] W. Wang and B. Li, nlin.CD/0208013.
[31] B. Georgeot and D.L. Shepelyansky, Eur. Phys. J. D 19,
263 (2002).
[32] N.R. Cerruti and S. Tomsovic, J. Phys. A 36, 3451
(2003).
[33] J. Van´ıcˇek and E.J. Heller, Phys. Rev. E 64, 026215
(2001).
[34] J. Van´ıcˇek and E.J. Heller, Phys. Rev. E 67, 016211
(2003).
[35] J. Van´ıcˇek and D. Cohen, J. Phys. A: Math. Gen. 36,
9591 (2003).
[36] W.H. Miller, J. Chem. Phys., 53, 3578 (1970).
[37] W.H. Miller, J. Phys. Chem. A 105, 2942 (2001).
[38] J.H. Van Vleck, Proc. Natl. Acad. Sci. 14, 178 (1928).
[39] S. Tomsovic and E.J. Heller, Phys. Rev. E 47, 282 (1993).
[40] J. Van´ıcˇek and E.J. Heller (not published).