Decoherence factor as a convolution: an interplay between a Gaussian and an exponential coherence loss

This paper identifies and investigates nature of the transition between Gaussian and exponential forms of decoherence. We show that the decoherence factor (that controls the time dependence of the suppression of the off-diagonal terms when the density matrix is expressed in the pointer basis representation) can be described by the convolution of Gaussian and exponential functions, their contributions modulated by the strength of the system-environment interaction. In the strong and weak coupling limits, decoherence reduces to the familiar Gaussian and exponential forms, respectively. The mechanism is demonstrated with two paradigmatic examples of decoherence -- a spin-bath model and the quantum Brownian motion.


I. INTRODUCTION
Quantum systems lose quantum coherence when coupled to their environments. Decoherence results in decay of the off-diagonal terms between system's pointer states-states that are immune to decoherence-in its density matrix [1][2][3]. It is an essential ingredient responsible for the quantum-to-classical transition [4,5] and the main obstacle in the development of quantum information technologies [6][7][8].
Decoherence is caused by the flow of quantum information from the system S into its environment E. Decoherence turns the initial superpositions of pointer states into their mixture. Decoherence factor is often presumed to decay exponentially, but this is not always the case. For instance, power law decay may appear [9][10][11][12][13][14]. Nonetheless, exponential decay is widely regarded as typical in various settings [15][16][17][18][19][20].
The interplay between Gaussian and exponential decoherence has not been-as yet-understood or even investigated. We show that decoherence factor can be described by a convolution f * g of an exponential f (t) = e −Γ|t|/2 and a Gaussian g(t) = e −σ 2 t 2 /2 : The Gaussian contribution originates from the spectrum: It is the Fourier transform of the spectral density (density of states) of the interaction Hamiltonian. Its origin is easiest to understand for pure decoherence, when the self-Hamiltonians of e.g. the spin-1 2 S and the N spin-1 2 E can be neglected [25,26]. The 2 N eigenvalues of the interaction Hamiltonian are then obtained by summing the couplings of S to the environment spins with the sign (+ or -) given by their relative orientations. Even for modest N the distribution of the eigenvalues approaches quickly a Gaussian. Gaussian spectrum arises in more general circumstances for physical systems with local interactions [27].
The exponential factor accounts for the response of the environment to the system-environment interactions: According the Fermi's golden rule, the environment decays from its initial state due to the coupling with the system with the rate Γ determined by the coupling strength.
The Fourier transform of the decoherence factor is then simply a product of two Fourier images. One is the energy spectrum; the other 'Breit-Wigner -like' represents the decay of the overlap between the eigenstates of the environment with and without coupling to the system. The overlap function will be defined rigorously in the next section.
These two contributions are responsible for the interplay between exponential and Gaussian time-dependence of the decoherence factor. In the limit of weak (Γ σ) and strong (Γ ∼ σ) coupling, the convolution-hence, the decoherence factor-exhibits predominantly exponential or Gaussian decay, respectively. Below we derive the convolution form of the decoherence factor in the spin-bath model, and then apply it to another paradigmatic model of decoherence, the quantum Brownian motion.

II. SPIN-BATH MODEL
The model consists of a single spin linearly interacting with a generic environment. The composite SE can be described by the Hamiltonian: where σ z is the Pauli operator, H I,E act on the environment, and λ is a dimensionless coupling strength. This Hamiltonian is relevant for a large class of pure decoherence processes of single-spin systems, and has been investigated in various contexts [25,26,28,29] including quantum phase transitions [30,31]. Consider the composite SE that starts as a product: and the central spin is prepared in a pure state |ψ S (0) = a|0 + b|1 . For the sake of simplicity, in the following analysis, the environment initial state ρ E is assumed to be an eigenstate |n of H E with energy E n . The results can be generalized (e.g., to a thermal state) by averaging with respect to the specific ensemble. Decoherence of the central spin is reflected in the decoherence factor r(t) that suppresses off-diagonal element of its reduced density matrix ρ(t): Upper and lower panels correspond to different initial energy eigenstates and coupling strengths, i.e., λ = 1 and 2, respectively. The first column is the plots of the full overlap functions, with a zoom-in shown on the second column and the third column (semi-log scale). Red solid (dashed) curves are the best fit to the Lorentzian (exponential) function. Inset is the histogram for the (Gaussian) spectral density.
Simple algebra shows that it is given by: where we have introduced an effective "perturbation" H P using the Baker-Campbell-Hausdorff formula, which appears in the same order as the coupling H I (see Appendix). The second line of the above equation can be interpreted as the survival probability of the initial state after an imperfect echo loop, i.e., a forward evolution followed by a perturbed backward evolution. This quantity is also known as the Loschmidt echo [32][33][34]. Denote the eigenenerges of the "perturbed" Hamiltonian, H E + λH P , as E k , and the corresponding eigenstate as |k = n C k n |n . The decoherence factor (5) further reduces to a Fourier transformation, Here, the discrete sum is replaced by a continuous integral with a "weighting function"-the global spectral density η(E). F (E k , E n ) ≡ |C k n | 2 is an averaged smooth function in continuous energy E k of the squared overlap |C k n | 2 between the perturbed and unperturbed eigenfunctions, which we call the overlap. The local density of states (LDOS), a.k.a the strength function, is defined as Both the LDOS and the overlap have peaks around the resonance energy E n . The LDOS was introduced and intensively studied in the field of nuclear physics [35] and encountered elsewhere (see e.g. [36]). In Ref. [35], using random matrix theory, it was derived to generally take the form where Γ(E) = 2πV 2 η(E) is given by Fermi's golden rule; V 2 is the average of the square of the off-diagonal elements of the perturbation matrix λH P , when written in the basis of the unperturbed Hamiltonian H E . It is known [37][38][39] that in the weak perturbation regime, where the width of the peak of the LDOS is much smaller than the width of the spectral density η(E), Γ(E) ≈ 2πV 2 η(E n ) can be treated as a constant. In this case, the LDOS reduces to a Lorentzian (Breit-Winger) distribution, whose Fourier transformation yields an exponential decay in time. On the other hand, for large perturbations, the LDOS can be approximated by a Gaussian [38,39], with a width bounded by the bandwidth of the spectral density.
To provide a unified treatment of the decoherence factor, including the cross-over regime, we first note that the LDOS (8) is valid for an arbitrary perturbation strength [35]. Comparing with definition (7), we identify the overlap function that takes the form This is not a precise Lorentzian distribution since the 'width' term is energy-dependent. Nevertheless, we will show that the overlap is insensitive to this energy dependence, and has an approximate (unnormalized) Lorentzian form, even for large perturbations that make the LDOS (7) far from a Lorentzian shape.
To establish this we perform a Taylor expansion of the energy-dependent width Γ(E) around the peak position E n . The LDOS (8) is of a Lorentzian form only at zeroth order of the expansion. However, the overlap F (E, E n ) has an exact (unnormalized) Lorentzian form up to the second order. Both first and second-order energy-dependent terms can be absorbed into the quadratic term (E − E n ) 2 . In this case, the effective width Γ eff of the overlap is where η(E) and its derivatives are all evaluated at E n . This expansion also shifts the effective peak position of the overlap by a small energy residue, The corresponding Fourier transform of the overlap (9) is then proportional to an exponential decay e −Γ eff |t|/2 multiplied by a phase, which, after absorbing the phase factor e iEnt in the decoherence factor (6), eventually contributes to a net oscillation e iErt .
In physical systems, the Hamiltonian is in general given by few-body interactions or local terms, whose spectral densities satisfy universal Gaussian distributions [27]. This, upon Fourier transform, gives rise to a Gaussian decay in the time domain, e −σ 2 t 2 /2 , where σ is the bandwidth (standard deviation) of the spectral density. Here, without the loss of generality, we assume the spectral density has a maximum at E 0 = 0. A non-zero E 0 only contributes a phase factor e iE0t to the decoherence factor.
The decoherence factor (6), up to an overall slow oscillation e iErt , is the Fourier transform of the product of two spectral images in the energy domain-the Lorentzian overlap and the Gaussian spectral density. The Fourier transform of such a product is a convolution of the two individual Fourier transforms of these two spectral images. They are exponential and Gaussian functions in the time domain, as derived above. Consequently, the decoherence factor is where Erfc is the complementary error function. The full derivation is presented in the appendix. For weak (Γ σ) and strong (Γ ∼ σ) couplings, this expression reduces to exponential and Gaussian decays, respectively. When we define the decoherence time as the time scale for the O(1) decay of the decoherence factor, in the exponential and Gaussian regime, the decoherence time scales as ∼ 1/Γ and ∼ 1/σ, respectively.
Recall also that λ quantifies the interaction strength between the system and the environment. According to Fermi's golden rule, for λ 1, Γ ∝ λ 2 . For the other extreme case, when the interaction eventually dominates the environment internal dynamics as studied in Refs. [25,26], the bandwidth of the Hamiltonian, σ ∝ λ. We conclude that the decoherence time scales as ∼ λ −2 and ∼ λ −1 in the weak and strong coupling limit, respectively.
It is worth stressing that, the above derivation of the convolution description of the decoherence factor is based on two explicit assumptions, that is, i) The energy spectrum of a typical physical Hamiltonian (with local interactions) takes a universal Gaussian form, and ii) the LDOS take a Lorentzian-like form (8). These two conditions can be derived using random matrix theory [27,35] under mild assumptions, and has been tested in numerous contexts. In this sense, we claim the convolution description of the decoherence factor of the spin-bath model is typical. Hence, for a physical system with local interactions, that can be well described by the spin-bath model, the decoherence factor is expected to take the form as a Gaussian-exponential convolution. On the other hand, there are exceptions for these conditions as well. For instance, one can get a strict exponential decay of the decoherence factor by manually tuning the energy spectrum of the bath in to an exponential form [25]. In the following discussion of the quantum Brownian motion, we will see another example, where at low temperatures, the non-Gaussian nature of the edge of the environment spectrum modifies the general behavior of the decoherence factor.
To illustrate the above results, we numerically simulate a specific bath term consisting many interacting spin-1/2 particles. More precisely, the linear coupling term and the bath internal Hamiltonian in Eq. (2) are; Decay of the decoherence factor at various coupling strengths λ in semi-log scale, with the corresponding overlap function (second column). The initial state of the environment is prepared as a pure energy eigenstate in the middle of the spectrum. Triangles are the numerical data, with the corresponding convolution functions (11) plotted in red curves. To guide the eyes, we also plotted the best fit to Gaussian (dotted) and exponential (dashed) functions.
where α = x, y, z, and σ α E,i are Pauli operators vectors for the environment spins while a i and b α ij are random numbers drawn from standard normal distribution.
In Fig. 1 (inset), the spectral density is shown to be well described by a Gaussian distribution. The claim for the Lorentzian distribution of the overlap is also justified. Note that the asymptotic tail of the overlap becomes exponential. This agrees with Wigner's derivation [40,41] for the tail of the strength function of banded matrices. The matrix of the physical perturbation represented in the basis of the unperturbed Hamiltonian is a banded matrix, i.e., its off-diagonal elements typically decay exponentially fast with the distance from the diagonal [42].
To visualize the change of the time dependence of the decoherence factor with the growth of the coupling strength (and hence, with the width of the overlap), here we focus on the case where the environment is initially prepared in an energy eigenstate. The results can be generalized to a thermal environment by averaging over a thermal distribution. This will be addressed in the study of quantum Brownian motion. Figure 2 depicts the decoherence factor at various coupling strengths, demonstrating the transition between exponential and Gaussian decays, with a crossover well-described by their convolution.

III. QUANTUM BROWNIAN MOTION
Encouraged by the above study of the spin-bath model, we expect that the convolution paradigm is robust and applies to other decoherence problems. We now test this conjecture with the familiar and extensively studied quantum Brownian motion (QBM) model [43][44][45][46]. Surprisingly, the possibility of the exponential to Gaussian transition has been so far overlooked.
QBM describes a single harmonic oscillator with a unit mass and a unit bare frequency, linearly interacting with a thermal bath of non-interacting harmonic oscillators. The interaction Hamiltonian between the system oscillator and

the bath oscillators is
where c k are the coupling strength. This model can be solved exactly [46,47]. It is well-known that, to study the reduced dynamics of the system, one can cast the effect of the bath oscillators into a single function, the spectral density function where m k and ω k are the masses and frequencies of the bath oscillators. Here, we consider a Ohmic environment, as it is widely studied in the literature. Its spectral density function is given by Above, Λ is a high-energy cutoff, and γ 0 quantifies the interaction strength. (Compared to Eq. (14), it can be seen that γ 0 is quadratic in the linear coupling strength between the system and environment oscillators). We provide detailed solution of the QBM in the appendix, where it is clear how the spectral density function enters the dynamics of the system oscillator. We assume the system and the bath are initially uncorrelated. The bath is in a thermal state with temperature T , and the system oscillator is prepared in a "Schrödinger cat" superposition of two Gaussian wavepackets located at different positions, i.e., ψ(t = 0) = ψ 1 + ψ 2 , with where N is a normalization factor. Due to couplings to the bath, the system oscillator evolves into a mixed state, which can be written as where ρ 1,2 are the reduced states resulting from the evolution of ψ 1,2 alone, and ρ int is the interference term between them [2,48]. Decoherence of the system is reflected in the decay of ρ int . To quantify this, we transform each term of the reduced state (17) into a corresponding Wigner representation in the phase space, W 1,2 and W int , and study the decay of the peak-to-peak ratio between the Wigner functions (see Appendix C.1): This quantity is known to be a good measure of decoherence for QBM [47]. Figure 3 depicts the decay of r W at various values of γ 0 . Both the exponential and Gaussian regime are accurately described by the convolution paradigm. For QBM, we define the decoherence time τ D through r W (τ D ) = 1/e. Since γ 0 in the spectral density is quadratic in the coupling strength between the system and environment oscillators [46], the decoherence time τ D is thus expected to scale as ∼ 1/γ 0 and ∼ 1/γ 1/2 0 , in the exponential and Gaussian regime, respectively. Both of these scaling behaviors are verified in the inset of Fig. 3. We present in the appendix the full solution of the model and the rationale for the convolution description.
As discussed before, the convolution behavior of the decoherence factor is accurate when the spectrum of the environment is well described by a Gaussian distribution, and the overlap is given by the Lorentzian-like distribution (9). These conditions may not be fulfilled when the state of the environment is close to the edge of the spectrum, e.g., the temperature is extremely low. In this case, the overlap cannot extend below the ground state energy. This low-energy cutoff may modify the shape of the overlap, and hence reduce the accuracy of our analysis in terms of convolution functions. To test this situation, we simulated r W (t) at zero temperature, shown in Fig. 4. It can be seen that at early times, r W (t) can still be accurately described by the convolution function (Fig. 4 left). This further demonstrated the robustness of the convolution paradigm. The low energy cutoff of the overlap mainly affects the long time behavior of its Fourier transform. Indeed, at long times, as shown in Fig. 4 (right), the decay of r W (t) is suppressed by a power-law.
It is worth stressing that, the range of γ 0 shown in Fig. 3 (main plot) is small, so that the time dependence appears to be Markovian [47]. We have developed an exact numerical treatment, aimed at exploring of the exponential-to-Gaussian transition through a wide range of coupling strength. We have focused on the Markovian regime in this paper because the exponential-to-Gaussian crossover turned out to occur already in this regime. The question is nevertheless whether the results presented above still hold when the coupling is strong. Our theory predicts that, for coupling strengths stronger than that in the crossover regime, we should observe Gaussian decays of the decoherence factor. The exact numerical treatment allows us to confirm that. As expected, we observed nearly perfect Gaussian decay for larger couplings (blue line in the inset of Fig. 3). Now the puzzling fact is that our main result predicts exponential-to-Gaussian transition from weak to strong couplings, respectively, with a crossover in between. However, for the QBM the crossovers have all been observed already in the weak coupling (Markovian) regime. This looks like a disagreement. The answer to this apparent problem is that the QBM dynamics do not quite fit into the spin-bath model. The appearance of the convolution and of the exponential-to-Gaussian transition is due to the presence of an effective "echo" operator (5) (detailed in Appendix C.3). As a result, this transition is not completely determined by the coupling strength γ 0 . Rather, the effective coupling strength is set by γ 0 together with the relevant energy scale of the system. Let us briefly discuss the mechanism here. Suppose the coupling term between the system S and the bath B is H I = H S I ⊗ H B I . In the truncated Hilbert space below the high energy cut-off, we can write the eigenvalues of H S I as {λ k } k=1,2... . As shown in the appendix, λ k serves as the effective "coupling strength" that controls the exponential-to-Gaussian transition. Now, λ k becomes proportional to the true coupling strength of the QBM, γ 2 0 . It also depends on the energy scale of the system, so that at high temperatures, λ k can extend to large values. Therefore, even though γ 0 is in the Markovian regime, the effective strength λ k can be in the strong regime.
The above argument also suggests that, at lower temperatures, we should expect that the crossover regime happens at larger values of γ 0 . This is indeed confirmed by comparing Figs. 3 and 4. Quantitative analysis of this interplay deserves further investigation.

IV. DISCUSSION
We have provided a unified and robust description of the interplay between the exponential to Gaussian decoherence in terms of the convolution functions. The relative import of these two types of decoherence is controlled by the interaction strength between the system and the environment as was seen in the two paradigmatic decoherence models, the spin-bath and the quantum Brownian motion.
The convolution description of the decoherence factor is explicitly derived using the standard spin-bath model, and therefore is not expected to work in any circumstance of decoherence. However, because of its simplicity, it is robust and can be used to model a broad class of systems. Our result straightforwardly applies to such scenarios. The quantum Brownian motion (QBM) does not fit into the paradigm of the spin-bath model. However, as elaborated in the previous section, the QBM exhibits similar structure that is responsible for the convolution description of the decoherence factor, namely, the appearance of the echo operator in (5). This is responsible for the observed accurate description of decoherence using convolution functions in the QBM model. Moreover, we expect that the convolution paradigm applies to a broader context than decoherence, whenever the same imperfect echo process is involved. This may explain the Gaussian to exponential transition in information scrambling [49], decay processes [50,51], as well as recent experiments [52] showing similar crossover behaviors (although study of the convolution paradigm in fields such as nuclear physics where LDOS has very different origins than in decoherence is beyond the scope of our discussion).
Decoherence is the main obstacle in implementing quantum computing. The time dependence of the decoherence factor will affect the error correction strategy (e.g., selection of the time intervals between the error correcting operations). We hope our result will provide valuable input into the hardware design in overcoming decoherence.
In practice, the interaction H I , as well as the effective perturbation H P , are much smaller compared to the environment Hamiltonian H E . This allows to expand echo operator as a power series. Intuitively, the effective perturbation H P should be twice as large as the interaction H I , as it only appears in the forward evolution loop. This can be verified by expanding the logarithm of the echo operator. Using the Baker-Campbell-Hausdorff formula, up to first order, we have the convolution theorem assertsF where the convolution * is defined as For convenience, we adapt the Fourier transform for the angular frequency 2πx, namely, with the corresponding inverse Fourier transform normalized by a pre-factor 1/2π. Under this convention, the convolution theorem readsF This allows to evaluate the decoherence factor in the main text as Assuming that the global spectrum density has a Gaussian shape with a band-width (standard deviation) σ, its Fourier transform is a Gaussian decayF [η(E)] = e −σ 2 t 2 /2 . (B10) As derived in the main text, the overlap takes the form with an effective width Γ eff and a small shift E r of the peak position E n . Its Fourier transform is given bŷ Here we omitted the normalization factor since the decoherence factor always start from unity. Evaluating the convolution of the above two Fourier transforms of the overlap and the spectral density, we arrive at the final expression of the decoherence factor Note that the above solution gives the unnormalized magnitude of the decoherence factor, which starts from unity at the initial time. The decoherence factor also has a residue oscillation e −iErt induced by the small shift E r of the resonance peak of the overlap function.

The solution
In this section, we collect the essential results of the solution that allow us to perform numerical simulations. More details of the derivations can be found in Refs. [46,47].
The total reduced density matrix (C6) can be solved using path integral [46] ρ(x, y, t) = dx dy J(x, y, t; x , y , t )ρ(x , y , t ), where the propagator, written in terms of the variables X = x + y and Y = x − y, is given by Here, the functions a ij (t) and b i (t) can be determined by the equations with boundary conditions u 1 (0) = u 2 (t) = 1 and u 1 (t) = u 2 (0) = 0, and µ(s), η(s) are the noise and dissipation kernels determined by the spectral density: We solve the above integrodifferential equation for u(t) with a shooting method incorporating the imposed boundary conditions. With these solutions, we can then compute the peak-to-peak ratio (C9) between the Wigner functions: where Here, the peak values for W 1,2 are the same, so we do not distinguish between them.

Convolution in QBM
As shown in the main text, the convolution function fits the decoherence factor of QBM very well. In this section, we explore the mechanism behind this. Instead of solving the reduced dynamics of the system harmonic oscillator alone, here we need to include the bath degrees of freedom into consideration.
The bath is prepared in a thermal state. However, for simplicity, in the following discussion we assume the bath is initially in a pure energy eigenstate. The result can be generalized to thermal states by doing a simple thermal average. Under this condition, the initial total wavefunction of the system oscillator and the bath is which results in two branches of the total wavefunction at any time instant t, i.e., The reduced states of the system oscillator from each branch, e.g., ρ 1,2 = Tr|Ψ 1,2,B Ψ 1(2),B |, are mixed states, whose Wigner representations (the first two terms in Eq. (C8)) in the phase space take Gaussian forms. As shown in Ref. [47], during time evolution, these two Gaussians are always separated and are approximately orthogonal to each other, as are the system's reduced density matrices ρ 1,2 (t) -they are (approximately) supported on orthogonal subspaces of the total Hilbert space. Therefore, we can decompose the total wavefunction at any time instant into: In this expansion, the basis states for the first branch, {ψ k 1 (t)} k , are orthogonal to those of the second branch, {ψ k 2 (t)} k . Here, we truncate the infinite dimensional Hilbert space to a finite (but large) N -dimensional Hilbert space within the relevant energy scale. By choosing a sufficiently large N , this can approximate the exact wave function to any degree of accuracy.
The above assumptions are specifically made for the QBM model and are justified by the known solutions. From now on, we will keep the discussion in a more formal and abstract manner beyond the QBM model. We write the linear interaction in a generic form As a matter of convention, we also choose the basis {ψ k 1,2 (t)} k as an eigenbasis of the interaction Hamiltonian H S Ithis is always possible since for any basis of the truncated Hilbert space we can diagonalize H S I and use the resulting eigenbasis to expand the wavefunction into the form (C27). The corresponding eigenvalues of the basis are denoted as {λ k 1(2) (t)} k . As introduced in the previous sections, the reduced state of the system oscillator can be factorized into where ρ 1 (ρ 2 ) is reduced from the first (second) branch of the total wavefunction, and ρ int comes from the interference part. In our 2N -dimensional truncated total Hilbert space, ρ int has matrix element With this, we can build up the evolution trajectories for the bath wavefunctions {ψ k B1 (t)} k and {ψ k B2 (t)} k . Consider the evolution of the total wavefunction from the given state Ψ(t) at time t to time t + ∆t with an infinitesimal ∆t.
where H 0 is the internal Hamiltonian (non-interacting part) of the total system. Hence, we conclude that the trajectory ψ k B1 (t) (and similarly for the second branch ψ k B2 (t)) can be interpreted as an evolution with a time-dependent Hamiltonian, i.e., where the time ordering ← − P is imposed. Plugging the above equation into Eq. (C30), we get the matrix element of the interference part of the reduced density matrix, Note that in the above equation we used a boundary condition imposed by the initial state of the bath, i.e., ψ k B1,2 (0) = ψ B . We also omitted the time ordering operator, but the integral is time ordered and is interpreted as an evolution generated by the Hamiltonian H 0 with small time-dependent perturbations λ 1,2 (t)H B I . Hence, the interference part of the reduced density matrix is related to the echo operators, which then contribute to the convolution function in the same manner as in the qubit model discussed in the main text. It is worth emphasizing that the above derivation is by no means a method for solving the reduced density matrix from the original model