Recovery of the fidelity amplitude for the Gaussian ensembles

Using supersymmetry techniques analytical expressions for the average of the fidelity amplitude f_epsilon(tau)=are obtained, where H_epsilon=H_0+(sqrt{epsilon}/(2 pi) )*V, and H_0 and H_epsilon are taken from the Gaussian unitary ensemble (GUE) or the Gaussian orthogonal ensemble (GOE), respectively. As long as the perturbation strength is small compared to the mean level spacing, a Gaussian decay of the fidelity amplitude is observed, whereas for stronger perturbations a change to a single-exponential decay takes place, in accordance with results from literature. Close to the Heisenberg time tau=1, however, a partial revival of the fidelity is found, which hitherto remained unnoticed. Random matrix simulations have been performed for the three Gaussian ensembles. For the case of the GOE and the GUE they are in perfect agreement with the analytical results.


Introduction
The concept of fidelity has been developed as a tool to characterize the stability of a quantum-mechanical system against perturbations [1]. Originally fidelity was introduced as the squared modulus of the overlap integral of a wave packet with itself after the development forth and back under the influence of two slightly different Hamiltonians. Let H 0 be the unperturbed Hamiltonian and the perturbed one. This somewhat unusual definition of the perturbation strength ǫ has been applied for later convenience. Then the fidelity is given by where ψ(0) is the wave function at the beginning, often chosen as a Gaussian wave packet with minimum uncertainty. It is assumed that H 0 has mean level spacing of one, and thus τ is given in units of the Heisenberg time. The variance of the off-diagonal elements of V is chosen to be one. In all what follows it is assumed that ǫ is of the order of one thus guaranteeing that the shift of the levels due to the parameter variation is of the order of the mean level spacing.
Depending on the strength of the perturbation one can discriminate roughly three regimes. In the perturbative regime, where the strength of the perturbation is small compared to the mean level spacing, the decay of the fidelity is Gaussian. As soon as the strength of the perturbation becomes of the order of the mean level spacing, a cross-over to exponential decay is observed, with a decay constant obtained from Fermi's golden rule [2,3]. For very strong perturbations the decay becomes independent of the strength of the perturbation. Here the decay is still exponential, but now the decay constant is given by the classical Lyapunov exponent [4]. It has been proposed by Pastawski et al [5] to look for such a behaviour in a spin-echo experiment on isolated spins coupled weakly to a bath of surrounding spins [6].
A paper of Gorin et al [7] is of particular relevance for the present work. The authors calculated the Gaussian average of the fidelity amplitude in the regime of small perturbations using the linear-response approximation, where C(τ ) is given by and 1 − b 2,β (τ ) is the spectral form factor. β is the universality index, i. e. β = 1 for the Gaussian orthogonal ensemble (GOE), β = 2 for the Gaussian unitary ensemble (GUE), and β = 4 for the Gaussian symplectic ensemble (GSE). For an explicit calculation knowledge of the spectral form factor is thus sufficient. By an exponentiation of the above formula, the authors were able to describe the cross-over from Gaussian to exponential decay with increasing perturbation strength quantitatively. The Lyapunov regime is non-universal and thus not accessible in a random matrix model. It is obvious that the linear-response approximation must break down for large perturbations. In the present work supersymmetry techniques are applied to calculate the ensemble average of the fidelity decay. Since this calculation is non-perturbative, the results hold for arbitrary values of the perturbation strength. We shall see that the calculation reveals an important generic feature, which is unaccessible by any perturbative approach.

The model
In the present paper we shall discuss the ensemble average of the fidelity amplitude, since for this quantity the calculation is much easier than for the originally introduced quantity (2). Since both the unperturbed Hamiltonian and the perturbation are taken from the Gaussian ensembles, the choice of the initial wave packet ψ(0) is irrelevant. The ensemble average may thus be written as where it is assumed that the Hamiltonian has been truncated to a finite rank N. The Hamiltonian introduced in equation (1) has the disadvantage that the mean density of states changes with ǫ. The more it is somewhat inconvenient for the present calculation that the variances of the matrix elements of H 0 and V differ. We therefore adopt a slightly different parameter variation, It is assumed that the matrix elements of H 0 and H 1 have the same variance, have zero average, (H 0 ) ij = (H 1 ) ij = 0, are uncorrelated, (H 0 ) ij (H 1 ) kl = 0, and are Gaussian distributed, where From random matrix theory it is known that the corresponding ensemble averaged density of states is given by Wigner's semi-circle law with a value of one in the centre of the circle, (see e. g. reference [8]). The resulting variances of the matrix elements of H 0 and H 1 are given by It follows for the GUE A similar expression is obtained for the GOE. Since only the difference of φ 2 and φ 1 enters expression (12), we may assume without loss of generality φ 1 = −φ 2 = φ/2. Ansatz (7) for the parameter variation is obtained from equation (1) by means of the substitutions where φ is thus of O( 1 √ N ), and ǫ is given in the limit of large N by ǫ = 4Nφ 2 .
Only terms up to O(φ 2 ) will survive the limit N → ∞ as we shall see later. The details of the parameter dependence are irrelevant. With all these substitutions equation (6) may be transformed into where with E ± = E ± ıη, and the abbreviations c = cos(φ/2), s = sin(φ/2). Using standard supersymmetry techniques [9], this can be written as where x = (x 1 , ξ 1 , . . . , x N , ξ N ) T , y = (y 1 , η 1 , . . . , y N , η N ) T , and We adopt the usual convention and use latin letters for commuting, and greek ones for anticommuting variables, respectively. Equation (18) is still true for all Gaussian ensembles, but now we have to discriminate between the GUE and the GOE.

The GUE case
Using equation (8), the calculation of the Gaussian average over H 0 and H 1 is elementary.
The result for H 0 may be expressed as where Whenever supermatrices are involved, traces and determinants are to be interpreted as super traces and determinants, respectively, in the definition of reference [9]. In short hand notation equation (21) may be written as where and 1 2 is the two-dimensional unit matrix. Introducing the notation where each S ij is a 2 × 2 matrix, and the indices 'A', 'R' refer to the 'advanced' and 'retarded' components, respectively, the sum entering equation (18) may concisely be written as where Next, a Hubbard-Stratonovich transformation is applied to equation (20), where U is the supermatrix with the 2 × 2 components For the integrals in equation (27) to be well-defined, the u ij integrations have to be performed from −∞ to ∞, whereas theū ij integrations are from −ı∞ to ı∞. (In literature usually an additional factor of ı is introduced in the lower right corner of the matrix (29) to avoid integrations along the imaginary axis.) In the same way we obtain where and with Collecting the results we obtain from equation (18) Now the x, y integrations can be performed resulting in Introducing the notation E 1/2 =Ē ± E/2, this can be written as where 1 4 is the 4 × 4 unit matrix, and L has been given in equation (23). Substituting theV integrations can be performed yielding where The second equation (39) has been obtained by an integration by parts. Equation (39) is still exact, but now the limit N → ∞ is performed. Since λ β is of O(N) (see equation (9) This suggests to diagonalize U, and perform the integrations over the elements of the diagonal matrix U D by means of the saddle point technique. The saddle points are obtained from the zeros of g ′ (u), whence follows The plus and the minus sign belong to the advanced saddle point u A , and the retarded one u R , respectively. The matrix U D at the saddle point is thus given by The matrix T diagonalizing U may be parameterized as where B and C are 2 × 2 supermatrices. Inserting equations (43) and (44) into equation [41), we obtain for the matrix U at the saddle point where In the last equation we used expression (10) for the mean density of states ρ.
We are now left with where only terms in φ surviving the N → ∞ limit have been taken. The brackets denote the average over the angular variables entering the matrix T , taken at the saddle point. Using equation (45) we obtain for the quantities entering on the right hand side of equation (47) Tr (U AA − U RR ) = 2ı∆ Tr BC , Tr The matrices B and C are diagonalized by means of the transformation where and (see e. g. Chapter 10 of reference [10]). Inserting these expressions into equations (48) to (50), we obtain, where expression (9) for λ β and expression (15) for ǫ were used, and For the calculation of the average (47) over the angular variables the 'surface volume' element is needed, (see again reference [10]). The integral over the anticommuting variables is easily performed. Only σ P and σ Q depend on the variables α, α * , and β, β * , respectively, and the corresponding integrals reduce to whence follows Collecting the results, we obtain from equation (47) The t, t * integration is over the whole plane, whereas thet,t * integration is restricted to the unit circlett * ≤ 1. Introducing polar variables, we obtain Inserting this result into equation (16), and introducingĒ = (E 1 + E 2 )/2 and E = E 1 − E 2 as new integration variables, we get, fixing the constant of proportionality by the condition f ǫ (0) = 1, TheĒ integration is nothing but an energy average. Restricting the discussion to the band centre, we may discard this average and obtain The integral is easily performed with the result where and s ′ (x) denotes its derivative. Equation (66) is the central result of this section. It gives an analytic expression for the GUE average of the fidelity amplitude for arbitrary perturbation strengths. f ǫ (τ ) and its first derivative are continuous, but the second derivative shows a discontinuity at τ = 1. A similar situation is known for the spectral form factor, where, however, for the GUE already the first derivative is discontinuous.
The solid lines in Figure 1 show the GUE fidelity amplitude for different values of the perturbation strength ǫ. For ǫ ≪ 1 the fidelity decay is predominantly Gaussian.
For ǫ = 1 and small times τ an exponential decay is found, with a cross-over to Gaussian behaviour at τ = 1, both observations in accordance with results known from literature. For ǫ ≫ 1 the fidelity decay is exponential for short times. Close to τ = 1, however, there is a conspicuous partial revival of the fidelity which had not been reported before, as it seems. For still longer times the decay becomes Gaussian again.
In the limit of small perturbations equation (66) reduces to This is in complete accordance with the results obtained by Gorin et al [7].

The GOE case
The first steps in the calculation of the ensemble average of the fidelity amplitude for the GOE are the same as for the GUE. Equation (22) for S remains correct, but now S is an 8 × 8 matrix with z n given by In taking the adjoint of z n one has to consider that the complex conjugate of the complex conjugate of an antisymmetric variable is defined as (α * ) * = −α, see the appendix of reference [9].
Up to equation (40) the further procedure is nothing but a step-by-step repetition of the calculation for the GUE case. The main problem for the GOE case arises from the diagonalization of matrix U, see equation (41). Not all of the matrix elements of S are different, as is evident from its definition, with the consequence that S obeys a number of symmetries which are inherited by the matrix U. The matrices T have to be chosen such that all symmetries are conserved. It is a highly non-trivial task to find the best parameterization for the matrix elements of T obeying these constraints. Fortunately this problem has already been solved by Verbaarschot, Weidenmüller, and Zirnbauer in their disseminating work [9]. We just cite their essential results: Equations (51) still hold, but now B D and C D are equal and given by The parameterization of the matrices P and Q is complicated, and is given in the appendices of reference [9]. For the present purpose it is sufficient to note that the angular averages over the matrices σ P = P σP −1 and σ Q = QσQ −1 (see equations (57) and (58)) up to a constant factor again yield the unit matrix, as can be shown by explicit calculation. All formulas of section 3 can thus be applied directly to the GOE situation. The surface volume element for the only remaining variables t 1 , t 2 ,t is given by and the integrations are from 0 to ∞ for t 1 and t 2 , and from 0 to 1 fort. Collecting the results, and proceeding in exactly the same way as for the GOE case, we finally end up with Substituting u = (x + y)/2 and v = (x − y)/2, we obtain where the constant of proportionality again was fixed by the condition f ǫ (0) = 1. For ǫ = 0 the right hand side of equation (74) must be one by construction, but it is not straightforward to show this explicitly. Since the corresponding calculation may be of some interest, it is reproduced in Appendix A. Equation (74) gives an explicit expression for the fidelity amplitude for the GOE case. It is not yet suited directly for a numerical integration, since the integrand contains a number of singularities. But it is not difficult to remove them by suitable substitutions of integration variables. This is done in Appendix B.
The dashed lines in Figure 1 show the results of the calculation for the same ǫ parameters as before. We notice that the partial recovery of the fidelity close to τ = 1 is still present for large ǫ values, but is considerably less pronounced than for the GUE case.

Numerical simulations
In this section we present random matrix simulations to affirm the analytical findings for the Gaussian orthogonal and unitary ensembles. Further we show numerical results for the Gaussian symplectic ensemble which has not been treated analytically.
In our simulations the Hamiltonians H 0 and H 1 are random matrices of dimension N × N with variances of the diagonal and off-diagonal elements given by equation (11). To calculate the fidelity amplitude, we write expression (6) as where In the numerical simulations the trace in equation (75) was restricted to 20 percent of the eigenvalues in the centre of the spectrum where the mean level density is still about constant. The average was taken over up to 8000 random matrices for H 0 , and for each of them over 50 random matrices for H 1 . For larger values of the perturbation strength ǫ it became more and more important to choose the dimension N of the matrices large enough to avoid finite-size effects. N = 500 proved to be sufficient for ǫ ≤ 10.
The results for the three Gaussian ensembles are shown in Figure 2. For the GOE and the GUE the numerical simulations are in perfect agreement with the analytical result for all ǫ values shown. For comparison, the fidelity amplitudes in the linear response approximation [7] (see equations (3) and (4)) are shown as well. For small perturbation strengths and small values of τ , the linear response result is a good approximation, but the limits of its validity are also clearly illustrated. In particular, it does not show any indication of the recovery near τ = 1.

Discussion
This work extends the results by Gorin et al [7] to the regime of strong perturbations using supersymmetry techniques. An intuitive explanation for the surprising recovery of the fidelity amplitude at the Heisenberg time can be given in terms of the Brownianmotion model for the eigenvalues of random matrices. The behaviour of the fidelity amplitude has its direct analogue in the Debye-Waller factor of solid state physics (see reference [11]). It is stressed that our result is generic and not restricted to random matrix systems. For instance, the fidelity recovery has recently been observed in a spin-chain model by Pineda et al [12].
The results of the present work may be easily extended to all situations, where the Gaussian averages (see equation (18)) lead to expressions allowing a subsequent Hubbard-Stratonovich transformation. This is, e. g., the case, if H 0 is taken from the GOE, and H 1 is from the GUE, or is purely imaginary antisymmetric. Perturbations, where the diagonal is zero, are of particular interest, since in such a situation the decay of fidelity freezes [13,14].
There might be still another application of the formulas derived in this paper. If f ǫ (τ ) is expanded into a power series of ǫ, the linear term can be expressed in terms of the spectral form factor, i. e. the Fourier transform of the two-point correlation function [7]. In a similar way the coefficients of the nth power of ǫ depend on all k-point correlation functions up to k = n. f ǫ (τ ) may thus be used as a generating function to obtain these terms in a simple way.