Loschmidt echo for local perturbations: non-monotonic cross-over from the Fermi-golden-rule to the escape-rate regime

We address the sensitivity of quantum mechanical time evolution by considering the time decay of the Loschmidt echo (LE) (or fidelity) for local perturbations of the Hamiltonian. Within a semiclassical approach, we derive analytical expressions for the LE decay for chaotic systems for the whole range from weak to strong local perturbations and identify different decay regimes which complement those known for the case of global perturbations. For weak perturbations, a Fermi-golden-rule (FGR)-type behavior is recovered. For strong perturbations, the escape-rate regime is reached, where the LE decays exponentially with a rate independent of the perturbation strength. The transition between the FGR regime and the escape-rate regime is non-monotonic, i.e. the rate of the exponential time-decay of the LE oscillates as a function of the perturbation strength. We further perform extensive quantum mechanical calculations of the LE based on numerical wave packet evolution, which strongly support our semiclassical theory. Finally, we discuss in some detail possible experimental realizations for observing the predicted behavior of the LE.


Introduction
One of the most prominent manifestations of chaos in classical physics is the hypersensitivity of the dynamics to perturbations in the initial conditions or Hamiltonian. That is, two trajectories of a chaotic system launched from two infinitesimally close phase-space points deviate exponentially from each other; so do the trajectories starting from the same point in phase space, but evolving under slightly different Hamiltonians. In a quantum system, it is natural to consider | φ 1 |φ 2 | 2 as a measure of 'separation' of two quantum states |φ 1 and |φ 2 . The unitarity of quantum propagators renders the overlap of any two states of the same system unchanged in the course of time. Thus, quantum systems are said to be stable with respect to perturbations of the initial state. However, a perturbation of the Hamiltonian can (and usually does) result in a nontrivial time dependence of the wave function overlap, suggesting a viable approach for describing instabilities and, therefore, for quantifying chaos in quantum systems.
Peres [1] proposed to consider the overlap of the state e −iH t/h |φ 0 , resulting from an initial state |φ 0 after evolution for a time t under the Hamiltonian H , with the state e −iH t/h |φ 0 obtained from evolving the same initial state through t, but under a slightly different (perturbed) HamiltonianH . He showed that the long-time behavior of depends on whether the underlying classical dynamics is regular or chaotic. 4 de Broglie wavelength. Analytical and numerical calculations yielded a novel LE decay regime, for which M(t) ∼ exp(−2γ t) with γ being the probability (per unit time) of the corresponding classical particle to encounter the boundary deformation. γ can also be viewed as a classical escape rate from a related open billiard obtained from the original (closed) one by removing the deformation-affected boundary segment. In this work, we explore all possible strengths of a local perturbation and describe the transition from the weak to the strong perturbation regime that completes the previous picture. In particular, we show that the rate of the exponential decay of the LE oscillates as a function of the perturbation strength as the escape-rate regime is approached. The paper is organized as follows. In section 2, we develop a comprehensive semiclassical approach of the LE decay due to local Hamiltonian perturbations of increasing strength. We perform a systematic analysis of the different decay regimes and establish their relation to the previously known decay regimes in the case of a global perturbation. In section 3, we validate our semiclassical theory by numerical simulations. In section 4, we outline possible experimental realizations and focus on the possibility of introducing a local perturbation in microwave-cavities and ultra-cold atom-optics billiards. We provide concluding remarks in section 5 and point to the similarities and differences with respect to another approach to the LE for local perturbations. Technical aspects of the calculations are relegated to the appendices.

Wave-function evolution
We address the time evolution of the wave function that describes a quantum particle moving inside a classically chaotic two-dimensional billiard (corresponding to a Hamiltonian H ). We assume that initially (at time t = 0) the particle is in a coherent state Here, σ quantifies the extension of the Gaussian wave packet, whereas r 0 and p 0 are the initial mean values of the position and momentum operators, respectively. We further define the (rescaled 5 ) de Broglie wavelength of the particle as λ =h/ p 0 . In our description of the time evolution of the wave function, we rely on the semiclassical approximation (see e.g. [24]) of the wave function at a time t, φ t (r) = dr ŝ(r,r ,t) Kŝ(r, r , t)φ 0 (r ). Here is the contribution to the Van Vleck propagator associated with the classical trajectorŷ s(r, r , t) leading from point r to point r in time t. Sŝ(r, r , t) denotes the classical action integral (or the Hamilton principal function) along the pathŝ. In a hard-wall billiard  (6), together with the conditions discussed in the text, allow to represent all the trajectoriesŝ contributing to equation (4) by the single reference trajectory s.
where Lŝ(r, r ) is the length of the trajectoryŝ, and m is the mass of the particle. In equation (5), Dŝ = |det(−∂ 2 Sŝ/∂r∂r )|, and the Maslov index νŝ equals the number of caustics along the trajectoryŝ plus twice the number of particle-wall collisions (for the case of Dirichlet boundary conditions).
Since we assume that the initial wave packet is localized around r 0 within σ , only trajectories starting at points r close to r 0 are relevant for our semiclassical description. Thus we can expand the action integral Sŝ(r, r , t) in a power series in (r − r 0 ). In appendix A, we show that the power series can be terminated at the linear term, if the wave packet is narrow enough, so that the condition is satisfied, where l L = p 0 /mλ is the Lyapunov length. In equation (6), s(r, r 0 , t) is the central reference trajectory into whichŝ(r, r , t) gets uniformly deformed as r → r 0 , and p s = −∂ S s (r, r 0 , t)/∂r 0 denotes the initial momentum of the trajectory s (see figure 1). Substituting equation (6) into (5), and performing the integration in equation (4) we obtain the momentum representation of the initial wave packet. 6 We now consider a related billiard, corresponding to the perturbed HamiltonianH , that differs from the original (unperturbed) billiard by a deformation of the boundary segment B 1 of width w (see figure 2). The perturbation is thus local, and will be characterized by its extent (depending on the ratio between w and the cavity perimeter P) and its strength (that will be quantified in the sequel). In view of equation (8), the wave function describing the evolution of the particle (starting from the same initial state φ 0 ) can be written as The sum now runs over all possible trajectoriess(r, r 0 , t) of a classical particle that travels from r 0 to r in time t while bouncing off the boundary of the perturbed billiard.

Wave-function overlap for local perturbations
According to equations (8)-(10) and the definition (1) of the LE amplitude, we have where A stands for the billiard area. The shadowing theorem [21] allows us to employ the diagonal approximation (s s) in the case of a classically small perturbation 6 , thus reducing equation (11) to where is the probability distribution of the particle momentum. In billiards the action difference between the two trajectories traveling between the same initial and final points in the same time t can be written, in terms of their length difference L s , as Using the Jacobian property of the Van Vleck determinant D s , we can replace the integral over final coordinates in equation (12) by an integral over the initial momenta and obtain The dependence of L on the product pt stems from the fact that in billiards, changing the magnitude of the momentum only modifies the traveling time, but does not affect the path.  (12) identifies both of them and assigns an action difference given by equation (14). The starting momentum of the solid red trajectory belongs to the set P 1 . The third trajectory (blue solid line) hits the boundary only at B 0 and therefore is the same for both the unperturbed and perturbed systems. Hence the action difference of the corresponding trajectory pair is zero. The starting momentum of the blue trajectory belongs to the set P 0 .
We now introduce a sequence of momentum sets P n (r 0 , t), with n = 0, 1, 2, . . ., such that for any p ∈ P n the classical trajectory, starting from the phase-space point (r 0 , p), arrives at a coordinate-space point r ∈ A after time t while visiting the deformation-affected boundary segment (B 1 in figure 2) exactly n times. Thus, the trajectories with the initial momentum in P 0 (r 0 , t) only undergo collisions with the part of the boundary unaffected by the deformation (B 0 in figure 2) rendering L = 0, so that with O n (t) = P n dpW 0 (p) e i p L(r 0 ,pt)/h , for n 1.
Since we are studying billiards, the integrations over momenta can be simplified by working in polar coordinates ( p, θ ) and considering the sets n (r 0 , pt) of angles θ such that p ≡ ( p, θ) ∈ P n (r 0 , t) iff θ ∈ n . For a classically chaotic dynamics, the set 0 shrinks with increasing time t, and 0 becomes the fractal set defining the repeller of the corresponding open (scattering) problem in the limit t → ∞. Equations (17) and (18) can, respectively, be written as For long times t, where many trajectories contribute to the semiclassical expansions, the angular integrals over n can be replaced by integrals over all angles, weighted with the measures ρ n of the corresponding sets. Assuming ergodicity one readily obtains an approximation for the probability for a trajectory of length l = pt/m to visit the boundary region B 1 exactly n times: where l d is the average dwell length of paths in the related open chaotic billiard obtained from the original (closed) one by removing the boundary region B 1 . This corresponds to the classical escape rate of the open cavity that for particles with momentum p 0 is given by For a chaotic cavity with an opening w much smaller than its perimeter P, we can approximate ( [25] and references therein) l d ≈ πA/w, and therefore In our case, the escape rate γ yields a measure of the perturbation extent. The classical escape rate of an open cavity controls the fluctuations of the transmission coefficients, and therefore approximations such as (23) have been thoroughly examined in the context of quantum transport ( [26] and references therein). According to the previous discussion, we can approximate O n (t), with n = 0, 1, 2, . . ., by the averagesŌ The mean value · · · n should be taken over the set n (r 0 , pt). The chaotic nature of the dynamics will enable us to treat the averages over n in a statistical way. In view of equation (13), the θ -integral in equations (24) and (25) where I 0 denotes the modified Bessel function.
As usually assumed in the LE studies, we restrict our analysis to 'semiclassical' initial wave packets φ 0 (r) with sizes much larger than the de Broglie wavelength, λ σ.
This assumption, together with condition (7), defines the interval for the dispersion σ , where the semiclassical approach is reliable, and hence yields restrictions on the parameters of the billiard. Employing condition (27) enables us to use the asymptotic form I 0 (x) ≈ e x / √ 2π x, valid for large x, in equation (26). Thus, the probability distribution function for the magnitude of the initial momentum is given by We note that equation (28) provides a good approximation to the exact distribution function already for σ 2λ. Under this assumption the p-integrals in equations (24) and (25) are dominated by the contributions around p 0 , and we can write (in view of equation (21)) Since in equation (30) all classical quantities are evaluated for an initial momentum with magnitude p 0 , the mean values · · · n should be taken over the sets n (r 0 , p 0 t). However, for long times and a chaotic dynamics we do not expect these mean values to depend on r 0 . In the next section, we will further invoke the chaotic nature of the underlying classical dynamics in order to estimate the mean values and therefore the LE average amplitude.

Averages over trajectory distributions
For classically small perturbations the action, respectively, length difference (equation (14)) between a trajectory s (solid red segment in figure 2) and its perturbed partners (dashed segment) is given only by the contributions accumulated along the n encounters with B 1 . Differences in length arising from the free flights between collisions with the boundary (B 0 + B 1 ) are of higher order in the perturbation strength and will not be considered. We can then write where the deformation function u(ϑ j , ξ j ) is the length difference accumulated in the jth collision with B 1 , depending on the impinging angle ϑ j ∈ (−π/2, π/2) and on the coordinate ξ j ∈ (0, w) of the hitting point. The number n of collisions with B 1 is a fraction of the total number of collisions p 0 t/ml f . The mean bouncing length l f can be approximated by πA/P [27], and we suppose l d l f since P w. We note that for small perturbations L depends on p only through n.
Given the chaotic nature of the classical dynamics and the fact that the collisions with B 1 are typically separated by many collisions with B 0 , we assume {ϑ j } and {ξ j } to be random variables. Furthermore, assuming a perfect randomization of the trajectories within the billiard the probability distribution functions for {ϑ j } and {ξ j } are, respectively, Then, treating the random variables as uncorrelated we obtain where the average . . . is defined as for a function f (ϑ, ξ ). Once we specify the shape of the perturbation, the average in the right-hand side of equation (33) is readily calculated from the probability distributions of equation (32). For instance, for a piston-like deformation (see appendix B), so that the average reads where H 1 stands for the Struve H-function of the first order, and J 1 is the first-order Bessel function of the first kind. Here we note that for the case of the piston-like boundary deformation the ratio h/λ serves as a measure of the perturbation strength, whereas w (and therefore γ ) quantifies the extent of the perturbation. At this stage, however, we will keep our discussion general and do not specify the details of the local perturbation. The substitution of equation (33) into (30) yields the average LE amplitude. The latter is usually not an observable quantity, however, it will be helpful towards our semiclassical calculation of the LE. The condition of a classically small perturbation that we have adopted throughout our work, implies that u 2 w P (which for the case of the piston-like deformation is equivalent to h w P). Quantum mechanically the perturbation can be characterized by the extent w and a deformation strength For χ 1, we will be in the quantum perturbative regime [7]- [9], which will not be considered in this work. Increasing the deformation strength χ, we anticipate a richer variety of regimes than for the case of the LE under global perturbations [5,6] since the perturbation extent enters as another relevant parameter.

LE for local perturbations
According to equations (2) and (11), the semiclassical expansion for the LE contains terms involving four trajectories. The diagonal approximation, leading to equation (12) for the LE amplitude, reduces the LE to a sum over pairs of trajectories. Consequently, the semiclassical form of the LE must take into account the different possibilities for each trajectory of the pair to hit (or not) the region of the boundary where the perturbation acts. We can therefore decompose the LE as where we have introduced the non-diagonal and diagonal contributions according to and Here, the set ε p of momenta p is defined such that two trajectories starting from the phase space points (r 0 , p) and (r 0 , p ), stay 'close' to each other in phase space during time t, and thus are 'correlated' with respect to the perturbation;ε p is the momentum set complementary to ε p . We give a quantitative definition of ε p below. Following the standard notation introduced in [5], we call the diagonal term the one resulting from the identification of pairs of trajectories where the effect of the perturbation is correlated, that is, when p ∈ ε p . In the non-diagonal term, we consider the pairs of trajectories uncorrelated with respect to the perturbation, including the case where one or both orbits are unperturbed. As noted at the beginning of this section, each of the trajectories of the above pair already incorporates a diagonal approximation between a perturbed and an unperturbed trajectory with the same extreme points.

Non-diagonal contribution to the LE
Calculating the LE as an average over trajectory distributions forces us to take into account pairs of trajectories and the possible correlations among them. The correlations are particularly important for M d (t), as we show in section 2.6. On the other hand, in our discussion in section 2.4, we established that for the calculation of M nd (t) the two trajectories of the pair can be considered to be uncorrelated with respect to the perturbation, and the averages can be performed independently. Assuming, in addition, that the measure of the momentum set ε p is small compared with that of the momentum set effectively represented by the distribution W 0 (p), we write Substituting equations (29), (30) and (33) into (41) and performing straightforward algebraic operations, we arrive at our central result for the non-diagonal contribution to the LE: where Thus, for the case of the piston-like boundary deformation, the average phase factor due to a single visit of the perturbation region B 1 by the particle is given by equation (36), so that equation (43) reads In section 2.7, we study the emergence of the FGR and escape-rate regimes of the decay of M nd (as well as the transition between the two regimes) depending on the strength of the perturbation. In sections 3 and 4, we discuss their implications for numerical simulations and possible experimental observations.

Diagonal contribution to the LE
To proceed with the calculation of the diagonal contribution to the LE, equation (40), we first need to specify the set ε p such that two trajectories of time t, starting from the phase space points (r 0 , p) and (r 0 , p ∈ ε p ), stay 'correlated' during time t. As in section 2.2, it is convenient to work with polar coordinates, p = ( p, θ) and p = ( p , θ ), in which the set ε p can be defined as follows: for every p ∈ ε p one has | p − p| p and |θ − θ | θ. In turn, p and θ are subject to the requirement that the two trajectories stay 'close' to each other in phase space. Indeed, any two 'correlated' trajectories must have the same number of collisions with the billiard boundary. This condition leads to | p − p|t/m l f . Moreover, the two trajectories must also have the same number of collisions with B 1 . The spatial separation between the two trajectories at the first collision with the boundary is approximately given by |θ − θ|l f . The condition that after time t this separation is smaller that the size w of B 1 is |θ − θ|l f exp (λt) w. Here, we have used the property that for chaotic dynamics two arbitrary initially close trajectories deviate exponentially from each over with a rate given by the average Lyapunov exponent λ. Thus, we can estimate the measure of the ε p -set to be Using this quantitative description of ε p for the evaluation of equation (40), we obtain for the diagonal contribution to the LE We now argue that for boundary deformations of moderate strength the exponent in the integrand on the right-hand side of equation (46) can be neglected. The argument of the exponent is given by the total differential of the function p L(r 0 , pt, θ), and therefore its absolute value can be bounded by Here we have used that, as discussed in section 2.3, L is independent of p for a fixed number n (≈ γ t) of collisions with B 1 . Then and p θ h 13 where we have introduced the dimensionless quantity This implies that the exponent in equation (46) is smaller than unity if χ l d /l f and C l f /w. Since the ratios l d /l f and l f /w are assumed to be large, the above inequalities hold for a wide range of deformations. (Note that for the piston-like deformation, equation (35), χ = (8/3) 1/2 h/λ and C = 0, and hence the exponent is small if h λl d /l f .) Neglecting the exponent in equation (46), we obtain In the second line of this equation, we have taken into account that p p and θ 1 for times t much longer than the dwell time t d . Under the assumption (27) of a 'semiclassical' initial wave packet, the integral over p in equation (51) is dominated by the contribution around p 0 and we find We note that the exponential dependence of M d on λt does not contain the perturbation and arises from a classical probability distribution, as in the standard Lyapunov regime [5]. The dependence with respect to the perturbation extent, w, appears in the prefactor.

Decay regimes of the LE
According to equation (38) the full LE M(t) is the sum of the non-diagonal and diagonal contributions, M nd (t) and M d (t), given by equations (42) and (52), respectively. We first argue that, unlike in the case of a global Hamiltonian perturbation, the non-diagonal contribution will typically dominate over the diagonal term. The most favorable regime to observe the diagonal term would be that of perturbations satisfying λ < κγ . Therefore, the necessary requirement for resolving the diagonal term is λ < 4γ . Since λ − γ = h KS > 0 (with h KS the Kolmogorov-Sinai entropy of the chaotic repeller [28]), we see that M d (t) can possibly prevail over M nd (t) only if h KS < 3λ/4. Taking into account that h KS = λ for a closed system, we see that only relatively open cavities would allow to observe M d (t). Hence for further analysis of the LE decay we mainly focus on the non-diagonal contribution M nd (t).
According to equation (42) the LE displays different decay regimes depending on the strength of the perturbation. Thus, for weak perturbations characterized by χ 1 one can expand the phase factor e iu/λ in the Taylor series to the second order in u/λ to obtain e iu/λ ≈ 1 − χ 2 /2 and, therefore, The rate of the exponential decay given by equation (53) depends on the perturbation strength χ, in analogy to the FGR regime found for global perturbations, but is dressed with γ that provides a measure of the fraction of phase-space affected by the boundary deformation. On the other hand, in the limit of strong perturbations, χ 1, the LE decay rate is independent of the   (43)) as a function of the perturbation strength χ . We have chosen a piston-like boundary deformation (see appendix B) for which the LE decay follows equation (44) and χ = (8/3) 1/2 h/λ. The FGR decay regime (equation (53)) is recovered for χ 1 (red dashed line). As χ → ∞ the decay rate asymptotically approaches the value 2 (black dashed line) representing the escape-rate regime (equation (54)).
perturbation strength χ and is entirely determined by the extent of the deformation quantified by γ . Indeed, e iu/λ → 0 as χ → ∞ leading to M nd (t) ≈ e −2γ t for χ 1 (escape-rate).
This is the escape-rate dominated decay regime previously reported in [23]. We finally emphasize that equation (42), and therefore equations (53) and (54), hold only if the conditions (7) and (27) are satisfied, i.e.
where l L is the Lyapunov length. As a particularly interesting feature, the transition between the FGR and escape-rate regime is non-monotonic, i.e. the LE decay rate κγ , in general, oscillates as a function of the perturbation strength χ while approaching the asymptotic value 2γ . Figure 3 illustrates these distinct oscillations for the case of a chaotic billiard with the Hamiltonian perturbation generated by a piston-like boundary deformation, see appendix B. In other words, for a fixed time t and starting in a minimum, the LE can increase by orders of magnitude upon varying the perturbation strength χ . The strength of the oscillations depends on the geometry of the boundary perturbation. The oscillations are particularly pronounced for the piston-type geometry where they only very slowly merge into the escape-rate limit.

Finite-size corrections to the semiclassical limit
In deriving our approximate expression for the non-diagonal contribution to the LE (equation (42)) we made the following two important assumptions. Firstly, we restricted our discussion only to those initial wave packets whose dispersion σ is much larger than the de Broglie wavelength λ, see equation (27). In numerical simulations, however, one can only address finite σ/λ ratios, so that the theory must be improved accordingly to be capable of accounting for the results of the simulations. Secondly, in section 2, we used the simple expression, given by equation (21), for the probabilities ρ n (l) for a classical trajectory of length l to visit the deformation region n times. While this expression well approximates ρ n (l) for l l d an improved formula is needed to correctly treat the numerically assessable range of lengths, l ∼ l d . In this subsection, we address these two important issues and present an expression for M nd (t) appropriate for a quantitative comparison with the results of the numerical simulations.
The non-diagonal contribution to the LE, as given by equation (41), is given by the square of the absolute value of the sum of overlapsŌ n (t) (with n = 0, 1, 2, . . .) defined according to equations (24) and (25). Now instead of using the asymptotic form of the Bessel function I 0 (like in section 2.2), we keep the analysis general by writingŌ n (t) as an integral over the dimensionless momentum variable z = p/ p 0 : and O n (t) = 2σ 2 λ 2 ∞ 0 dz z e izu/λ n ρ n (zp 0 t/m) exp − σ 2 λ 2 (z 2 + 1) I 0 2σ 2 λ 2 z , n 1.
We now address the probability distributions ρ n (l). The central building block of our analysis here is the probability g(l) dl for a classical trajectory with the length between (L + l) and (L + l + dl) to end on the boundary deformation region B 1 subject to the condition that the previous encounter with B 1 took place at length L. Assuming ergodicity for the billiard system under consideration we approximate g by the Heaviside step function θ as where the length l 0 is the (average) minimal length that a trajectory starting on B 1 must have to return to B 1 . In view of equation (58), one readily obtains the following approximation for the survival probability: Then, the visit probabilities for n 1 are calculated as  ×g(l n − l n−1 ) ρ 0 (l n − l n−1 ) · · · g(l 2 − l 1 )ρ 0 (l 2 − l 1 ) g(l 1 ) ρ 0 (l 1 ) (60) 16 and can be shown to satisfy the following recursion relation for n 1: Note that equation (61) together with (59) simplifies to equation (21) if one puts l 0 = 0. However, for trajectories of lengths comparable to l d the minimal return length l 0 cannot be neglected and equations (59) and (61) must be used in equations (56) and (57) to yield the LE.

Numerical simulations
In order to support our semiclassical predictions, we present in this section numerical quantum mechanical calculations for a local perturbation. We use the Trotter-Suzuki algorithm [29,30] to simulate the dynamics of a Gaussian wave packet inside a desymmetrized diamond billiard (DDB). The DDB is defined as a fundamental domain of the area confined by four intersecting disks of radius R centered at the vertices of a square. We denote the length of the largest straight segment of the DDB by L (see figure 4). As proved in [31], the DDB is fully chaotic and thus has been previously considered for studying aspects of quantum chaos [32]. Our semiclassical analysis is valid for an arbitrarily shaped local perturbation acting on a region B 1 (of width w) of the boundary. A perturbation with the shape of a circular segment was used in [23]. In our present numerical simulations, we chose a piston-like perturbation ( figure 4(b)), for which analytical results can be readily obtained (see equations (35) and (44), and appendix B).
In figure 5, we present the LE decay for a DDB with L = 1000, R = 1311 (in units of the lattice spacing of the underlying tight-binding model) and a piston-like perturbation. For the present geometry, the dwell length is l d = (P/w)l f ≈ 18.7 l f , with l f being the mean free flight path. The initial wave packet has the dispersion σ = 8; its momentum direction is chosen to be parallel to the longest straight segment of the DDB (α = 0 in figure 4(b), but we have verified that the LE decay rate is independent of α. 7 The three numerical (solid-line) curves in figure 5 correspond to the following values of the deformation strength: (i) χ ≈ 1.63 (h = 3 and λ = 3), (ii) χ ≈ 1.75 (h = 3 and λ = 2.8), and (iii) χ ≈ 5.44 (h = 10 and λ = 3). These numerically obtained LE curves decay almost exponentially for times up to γ t ≈ 4 before turning over to a regime with strong irregular fluctuations around a saturation value [1]. The three curves shown illustrate the non-monotonic dependence of the decay rate with the perturbation strength.
The corresponding semiclassical LE decay curves-the three dashed-line curves in figure 5-were obtained by doing the integrals in the right-hand side of equations (56) and (57) numerically, with ρ n (l) probability distributions determined in accordance with equations (59) and (61); the minimal return length was taken to be l 0 = 0.16 l d . The infinite sum in the righthand side of equation (41) was terminated at n = 8.
The good agreement between the semiclassical predictions and the results of the full quantum-mechanical computation is evident in figure 5. The fact that the obtained LE decay deviates from the purely exponential one is entirely due to the finiteness of the σ/λ ratio: the analytical results of section 2 are recovered in the limit given by equation (55), which proves challenging in numerical simulations 8 . As we discuss in the next section, however, this limit is naturally achievable in laboratory experiments with ultra-cold atom-optics billiards, so that the latter provide a viable model system for studying the LE from local Hamiltonian perturbations.

Experimental realizations of LE with a local perturbation
Experiments on the LE are of foremost importance since they render crucial information about quantum dynamics of physical systems and their decoherence mechanisms 9 . While the examples discussed in the introduction show that the agreement between the semiclassical 8 Scaling σ → nσ also requires l L → n 2 l L to keep the ratio σ/ √ λl L unchanged. This, in turn, is equivalent to increasing the linear size of the billiard n 2 times. The increase of the number of the computational grid points is then given by the factor n 4 , and the number of time steps to propagate the wave function through the same total time multiplies by the factor n 2 . This gives the overall factor n 6 for the code running time. 9 For a recent account see section 9 in [6]. LE experiments were first performed on nuclear spins of organic molecules using NMR techniques [4,33]. The decay of the polarization was found to be quite insensitive to the coupling to external degrees of freedom or the precision of the reversal. The Gaussian decay of the experimentally measured LE is at odds with the one-body semiclassical theory, and many-body aspects of the problem have been pointed to be at the origin of such a behavior [10,34,35].
In [23], a principle experimental scheme for measuring the LE for local boundary perturbations was proposed based on a ballistic electron cavity with a small ferromagnet attached acting as the local perturbation. Such a setting provides a link between spin relaxation in a mesoscopic cavity and LE decay. Here, we discuss two further experimental settings which appear suitable for a measurement of the echo decay: microwave and cold atom cavities.
Microwave experiments allow the independent measurement of individual scattering matrix elements for the unperturbed and perturbed systems [36]. The cross-correlation of these matrix elements can then be calculated, and going into the time domain, the scattering fidelity is obtained. The latter is a good representation of the usual average fidelity amplitude when appropriate ensemble and/or energy averages are taken. Correspondingly, the LE can also be constructed from measurements of scattering matrix elements. The observation of the Lyapunov decay regime for a global perturbation in the microwave cavity can then be envisioned. However, reaching the required long time scales remains as an experimental challenge. On the other hand, the corresponding (escape-rate) decay regime for local perturbations might be easier to reach experimentally. Moreover, microwave billiards appear to be rather suited to investigate local boundary perturbations, since the piston-type deformation presented in our work can be directly realized in a microwave billiard setup. The width of the piston determines the exponent of the LE time decay. Furthermore, by moving the piston the perturbation strength χ can be directly controlled and tuned. Hence, by devising sufficiently large microwave cavities to approach the semiclassical limit it seems promising to experimentally reach both, the escape-rate regime for large χ and the non-monotonic dependence of the decay rate of the LE on the perturbation strength (see figure 3).
Studying quantum chaos in the laboratory by recreating a delta-kicked harmonic oscillator in an ion trap was proposed a decade ago by Gardiner et al [37], and a number of fruitful approaches have since then been developed and successfully realized using ultra-cold atoms confined to optical billiards [38]. For instance, the decay of quantum correlations has been measured by echo spectroscopy on ultra-cold atoms using the detuning of the trapping laser as a perturbation [39]. Below we focus on the time evolution of clouds of ultra-cold atoms in optical billiards, and show that they provide a viable system for experimental investigation of different perturbations and various regimes of the LE decay. The perturbations can be global, such as in the cases previously studied, but also local. Since the large-scale separation of system parameters, given by the condition (55), is attainable in these experimental systems we expect that a direct support for the theoretical predictions of section 2 can be obtained.
In a typical microwave echo (or Ramsey) spectroscopy experiment [39], a cloud of ultracold Rb atoms is loaded into an off-resonance optical trap. For the purpose of our study the role of the trap can be played by a hollow laser beam with the cross section corresponding to the geometry of a chaotic billiard of interest. The fabrication of such hollow laser beams, as well as the manipulation of atoms inside them, can now be performed with a high level of precision [40]. The atomic cloud, after being positioned inside the hollow beam in its focal plane and accelerated (or 'kicked') as a whole to a nonzero average momentum, is allowed to evolve freely in an effectively two-dimensional billiard, see figure 6.
The Rb atoms used in echo experiments [39] are initially prepared in a quantum state 0 equal to a direct product of an internal atomic state |s and a spatial state described by a wave function φ 0 (r), i.e. 0 = |s ⊗ φ 0 (r). The internal state evolves in a coherent superposition of the two hyperfine sub-states, denoted by |↓ and |↑ , of the ground state of rubidium. The |↓component of the total wave function of an atom experiences a laser field potential V ↓ (r) which is, in general, different from the potential V ↑ (r) exerted by the same laser on the |↑ -component. The relative difference between the two optical potentials is given by the ratio ω HF / L , wherē hω HF is the energy of the hyperfine splitting of the ground state, and L is the laser detuning from the frequency of the transition between the ground state and the first allowed excited state of the Rb atom. The application of a sequence of π/2 microwave pulses during the time evolution of the atoms, followed by a measurement of the populations of the |↓ -and |↑substates at the end of the evolution, allows one to determine the LE (corresponding to the spatial wave function φ 0 (r)) due to the difference between the potentials V ↓ (r) and V ↑ (r) as a function of the evolution time.

Perturbation laser beam
Atomic cloud

Momentum kick
Confining hollow laser beam In order to measure the LE decay due to local perturbations two different lasers have to be used. The first laser is to produce the confining hollow beam with the cross section of a desired (billiard) geometry, and has to be tuned as to exert approximately the same potential V bill (r) on the both |↓ -and |↑ -substates. The beam of the second laser plays a role of the local Hamiltonian perturbation. It has to be placed inside (and aligned with) the hollow beam of the first confining laser, and its width should be much smaller than the linear scale of the billiard (see figure 6). The frequency of this second laser (and perhaps its position inside the billiard) determines the perturbation strength χ. Altering this frequency changes the difference between the potentials δV ↓ (r) and δV ↑ (r), produced by the second laser and acting differently on the |↓ -and |↑ -substates, respectively. Thus, an echo spectroscopy experiment performed in such a system would measure the LE decay due to the difference of the atomic potentials V ↓ (r) = V bill (r) + δV ↓ (r) and V ↑ (r) = V bill (r) + δV ↑ (r); this difference is localized in an area much smaller than that of the billiard.
To date one is typically able to experimentally prepare and manipulate clouds of Rb atoms as cold as 1 µK. This temperature corresponds to the thermal speed of about 1.3 cm s −1 . At the same time, by first placing the atoms inside a far-off-resonance dipole trap then moving the trap and finally switching it off one can accelerate the atomic cloud as a whole up to 10 cm s −1 . Such a momentum kick can nowadays be easily realized in a laboratory, and does not significantly increase the temperature of the atoms. As a result one obtains a cloud of atoms moving as a whole with an average momentum that corresponds to a rescaled de Broglie wavelength λ ∼ 10 nm. The dispersion of the atomic cloud can be shrunk to σ ≈ 1 µm. Under these conditions, the number of Rb atoms composing the cloud can reach 10 5 , which is easily sufficient for Ramsey-spectroscopy-type measurements. The hollow laser beam producing the billiard confinement can reach L ≈ 1 cm in linear size. Assuming the Lyapunov scale l L of the billiard to be of the same order of magnitude as L one arrives at the following estimate for the scale separation of the system parameters: λ : σ : λl L ∼ 1 : 10 2 : 10 3 , which well satisfies the restriction given by equation (55). Under the conditions specified above, one can control the atoms for up to 5 s before the cloud gets significantly elongated in the axial direction of the hollow laser beam. Such a time allows a single atom to experience about 50 bounces with the billiard boundary, which according to our theory is sufficient for observing the LE decay regimes predicted in this work. The above considerations show that atom-optics billiards constitute promising candidates for experimentally investigating different regimes of the LE decay due to local perturbations of the Hamiltonian, and in particular put in evidence the escape-rate regime and the predicted oscillations of the LE as a function of the perturbation strength.

Conclusions
In this work, we have studied the time decay of the LE in quantum systems that are chaotic in the classical limit, due to Hamiltonian perturbations localized in coordinate space. We have provided the corresponding semiclassical theory of the LE for coherent initial states evolving in two-dimensional chaotic billiards.
In addition to the FGR decay regime, which is well known for the case of global Hamiltonian perturbations and is recovered in our theory for weak (χ 1) local perturbations, our analysis predicts a novel decay regime for strong (χ 1) perturbations that stems entirely from the local nature of the Hamiltonian perturbation, i.e. the escape-rate regime, and quantitatively describes the transition between the FGR and escape-rate regimes, as the perturbation strength is varied. In the escape-rate regime, the LE decays exponentially in time with a rate equal to twice the escape rate from an open billiard with the 'hole' at the place of the perturbation. Hence, the LE allows to mimic the decay behavior of a system without opening it. In this regime, the LE decay rate is independent of the deformation strength χ. The transition between the FGR regime and the escape-rate regime turns out to be non-monotonic: the rate of the exponential time-decay of the LE oscillates as a function of the perturbation strength.
We would like to point out that recently there has been another study [41] of the LE decay due to local perturbations. It addresses a particle moving in a two-dimensional array of point-like scatterers, and the perturbation of the Hamiltonian is achieved by slightly displacing one of the scatterers. The theoretical and experimental analysis of the system reveals a polynomial decay of the LE, namely M(t) ∼ t −2 , for long times, which cannot be compared with our findings: that work is in the perturbative regime, where the eigenstates are not significatively modified by the perturbation.
We have also performed an extensive numerical study of the LE decay to support our semiclassical theory. To this end, we have simulated the time evolution of initially coherent states in the DDB. The role of the local Hamiltonian perturbation was played by a piston-like deformation of the billiard boundary. The results of our numerical simulations exhibit strong quantitative agreement with the predictions of our theory extended to cope with initial states given by Gaussian wave packets of a dispersion comparable with the de Broglie wavelength.
While the scale separation given by equation (55) is rather challenging to satisfy numerically it can be naturally achieved in laboratory experiments with ultra-cold atoms confined to optical billiards. In this work, we have proposed a laboratory setup allowing one to investigate the LE decay from local Hamiltonian perturbations for a wide range of perturbation strengths, and to observe the predicted decay regimes. We believe that the study of the LE decay due to local perturbations provides an example of physical problems for which capabilities of laboratory experiments go beyond those of numerical simulations. Such experiments may also reveal weak-localization-type quantum corrections to the LE decay which are expected from an analysis [42] of loop contributions [43] beyond the semiclassical diagonal approximation.
For a billiard, the off-diagonal partial derivatives can be neglected compared to the diagonal ones, so that condition (A.3) can be replaced by h . (A.5) The first of the two derivatives in equation (A.5) is ∂ p 0 /∂q 0 = −m/t for a particle of mass m in a billiard. To evaluate the second derivative, we first linearize the trajectoryŝ(r, r , t) around s(r, r 0 , t), so that q ⊥ τ ≈ q ⊥ τ (q ⊥ 0 , p ⊥ 0 , τ ). Therefore, The right-hand side of equation (A.7) is given by the ratio of the two monodromy matrix elements. To facilitate our analytical presentation, we use here the monodromy matrix of the dynamics on Riemann surfaces of constant negative curvature 10 The action integral expansion (6) requires condition (A.11) to be satisfied. We finally note that in the limit λ σ , which we utilize in section 2, equation (A.11) simplifies to σ √ λl L .

Appendix B. Piston-like boundary deformation
In this appendix, we explicitly compute as an example the length-difference function u(ϑ, ξ ) of equation (31) for a piston-like local boundary deformation, see figure B.1, which is also used in our numerics. We assume that the boundary of the unperturbed billiard possesses a straight segment of length w that gets 'lifted' by the perturbation as if an imaginary 'piston' was pulled out. We denote the piston displacement by h. Assuming h much smaller than the free flight path l f of the trajectory hitting the deformation, we treat the unperturbed and perturbed trajectories as being parallel. Then the length difference u(ϑ, ξ ) accumulated due to a single collision with the deformation-affected segment of the boundary is given by an expression analogous to Bragg's diffraction formula, u(ϑ, ξ ) ≈ 2h cos ϑ. (B.1) Here, ϑ represents the collision angle as shown in figure B.1. The length difference u can be considered independent of the collision coordinate ξ for deformations such that h w. Taking into account the probability distribution function of collision angles, equation (32), we have cos ϑ = π/4 and cos 2 ϑ = 2/3. Consequently, the first two moments of the deformation function read u = π 2 h and u 2 = Similarly, the average phase factor due to a single collision of the particle with the piston is given by equation (36).