UvA-DARE

Decoherence is the main obstacle to quantum computation. The decoherence rate per qubit is typically assumed to be constant. It is known, however, that quantum registers coupling to a single reservoir can show a decoher-ence rate per qubit that increases linearly with the number of qubits. This eﬀect has been referred to as superdecoherence, and has been suggested to pose a threat to the scalability of quantum computation. Here, we show that su-perdecoherence is absent when the spectrum of the single reservoir is continuous, rather than discrete. The reason of this absence, is that, as the number of qubits is increased, a quantum register inevitably becomes susceptible to an ever narrower bandwidth of frequencies in the reservoir. Furthermore, we show that for superdecoherence to occur in a reservoir with a discrete spectrum, one of the frequencies in the reservoir has to coincide exactly with the frequency the quantum register is most susceptible to. We thus fully resolve the conditions that determine the presence or absence of su-perdecoherence. We conclude that superdeco-herence is easily avoidable in practical realizations of quantum computers.


Introduction
In principle, quantum computers can solve problems that are intractable on any classical computer.The largest obstacle to bringing this in practice is decoherence [1], and it is essential to understand the sources and effects of decoherence under practical circumstances encountered in actual quantum computers.As we inch towards full-scale quantum computing, where we are already facing systems with on the order of a hundred qubits [2][3][4], the system size dependence of decoherence becomes of increasing importance.
Decoherence is commonly studied in a simplified spin-boson model, where only the dephasing effects of the bosonic bath are taken into account [5][6][7][8][9][10][11][12][13][14].Henceforth we will refer to this model as simply 'the dephasing model'.This model is exactly solvable, and at the same time broadly relevant because dephasing times are typically much shorter than relaxation times [7,15,12].It should be noted, how-ever, that there are situations where it does not accurately describe the decoherence process because of non-perturbative effects [15,12].If, in the dephasing model, each qubit is assumed to couple to its own, independent reservoir, the decoherence rate per qubit is constant.If, on the other hand, the qubits couple to single reservoir, the decoherence rate per qubit scales linearly with the number of qubits for certain states [6,16,17,7,8,18,19].This effect has been referred to as superdecoherence, in analogy with superradiance.
Superdecoherence has been predicted in Refs.[7,15], and has been observed experimentally in an ion-trap quantum computer [20].Although some states suffer superdecoherence, the probability of running into such a state during the course of an actual algorithm may be extremely small [17].Additionally, if the decoherence is dominated by relaxation, rather than dephasing, is has been shown that superdecoherence does not occur for the Greenberger-Horne-Zeilinger (GHZ) and the Hadamard state [21].Also the particular model of solid-state qubits coupling to a single phonon reservoir has been shown not to give rise to superdecoherence [22].The latter approach focuses on a specific setting of the dephasing model: the geometry of the quantum register is assumed to be a linear array, and the phonon reservoir is assumed to be three-dimensional and thermal, with a continuous spectrum and a linear dispersion relation.Therefore, it is unable to reveal the general underlying physical reasons for the absence of superdecoherence.The reason why superdecoherence emerges in other settings of the single-reservoir dephasing model remained unknown.
Here, we fully resolve the physical conditions that determine the presence or absence of superdecoherence in the dephasing model, with all qubits coupling to a single bath.We do not make any assumptions about the geometry of the quantum register, the dimension d, the reservoir dispersion relation, or the directional dependence of the spin-boson interaction.For the reservoir state, we assume a very general initial condition that applies to practically relevant situations.In this general setting, we find that the (im)possibility of superdecoherence due to a single reservoir is completely determined by the boundedness of the spectral density and the occupation density of the reservoir.
The spectral density is the density of modes at a given frequency.If the reservoir admits only a discrete set of frequencies, such as the electromagnetic field in an ideal cavity, the spectral density is given by a sum of delta functions, and is hence unbounded.If, on the other hand, the reservoir admits a continuum of frequencies, such the electromagnetic field in an imperfect cavity or free space, the reservoir spectral density is a bounded function of frequency.
The occupation density, on the other hand, tells us to what extent a given mode in the reservoir is exited.It is typically a bounded function of the mode frequency.However, if only a single frequency is excited, the occupation density is described by a delta function which is centered at that frequency.This is the case when the bosonic field is the electromagnetic field, and a mode is excited by a laser with vanishing spectral bandwidth.In contrast, if this laser has a nonzero spectral bandwidth, also the occupation density remains bounded.
Here, we prove that superdecoherence is absent in single-reservoir dephasing if both the reservoir spectral density and the reservoir occupation density are bounded.We call these reservoirs continuous because, in this case, both spectra have continuous support.An important physical quantity in the proof is the dephasing susceptibility, which we define as the only part of the decoherence rate that depends on the system.It is closely related to, but different from, the so-called array factor, which arises in classical antenna arrays [23], quantum antenna arrays [24], and interdigital transducers that couple to surface acoustic waves [25].The dephasing susceptibility captures the extent to which a reservoir frequency contributes to the dephasing process if this frequency is present in the reservoir.
The reason for the presence of superdecoherence in some situations, is that there may be frequencies for which the dephasing susceptibility scales quadratically with the number of qubits.If one of these frequencies coincides with a frequency for which either the reservoir spectral density or the reservoir occupation density diverges, superdecoherence is exhibited.This is because, in this specific case, the decoherence rate scales with the system size in the same way as the peak of the dephasing susceptibility.This explains why superdecoherence is exhibited when either the spectral density or the occupation density is unbounded.
The reason for the absence of superdecoherence in continuous reservoirs (i.e.bounded spectral density and occupation density) is that peaks in the dephasing susceptibility inevitably become narrower as the system size is increased.Specifically, we show that, if the dephasing susceptibility has a peak whose height scales as the square of the number of qubits, the width of this peak must scale inversely with the number of qubits.That is to say, the quantum register may be increasingly susceptible to a given reservoir frequency a analogue, where a linear array of L = 10 classical dipoles, with lattice spacing a, is placed in the electromagnetic field.As a whole, the array couples strongly to the mode with wave number k = 0 (not shown).The array does not couple at all to modes with wave number k = ±2π/(aL) (shown in blue).This is because, for these modes, all potential energies arising form the dipole-field interaction cancel exactly.(Bottom) The modulus squared of the coupling strength in the classical analogue, as a function of the wave number k, depicting two modes that do not couple to the array in blue.The array mainly couples to modes in a bandwidth less than ∆k = 4π/(aL), which goes as ∼ 1/(aL).
as the system size grows, but the bandwidth of this susceptibility must at the same time decrease.This effect mitigates the total decoherence rate, and the net effect is that superdecoherence is suppressed.

A classical analogue
The cause of the inverse scaling of the bandwidth of the susceptibility, which is responsible for the absence of superdecoherence in continuous reservoirs, can be sketched with a classical analogue.We leave the treatment of the quantum dephasing susceptibility for Sec. 3. Consider L identical, classical, noninteracting electric dipoles in a linear array with spacing a, as depicted in Fig. 1 (top).(This geometry is chosen for explanatory reasons.Our results concerning the quantum dephasing susceptibility hold for general register geometries.)In the initial state of the array, all dipoles point upwards.For simplicity, consider only the electromagnetic modes whose momentum is colinear with the array and are polarized in the direction of the dipole moments.The dipoles couple to the electromagnetic field, giving an initial potential energy V = C 1 L =1 E , where C 1 is some constant, E is the electric field at the th dipole, and k the wave number.In terms of the Fourier transform E(k) := L =1 e ikr E , where r = a( − 1) is the position of the th qubit, the initial poten- the coupling strength between the array and the mode with wave number k. See Fig. 1 (bottom) for a plot of |f (k)|2 .From the previous expression for f (k), and from the plot, we can see the array couples most strongly to the electromagnetic field mode with wave number k = 0. We can also see that the array does not couple at all to modes with wave number ±2π/(aL).In real space, this is because, for this wave number, all potential energies cancel exactly [also see Fig. 1 (top)].Thus, the bandwidth of modes to which the array couples strongly is at most ∆k = 4π/(aL), which scales inversely with the length of the array.

Spin-boson dephasing
In this section, we introduce the model of spin-boson dephasing, following references [6][7][8].First, we consider the case of a single qubit coupling to a bosonic reservoir, and extend this to multiple qubits, each of which couples to its own, independent, bosonic reservoir.In both of these cases, superdecoherence cannot occur under any circumstance.Subsequently, this situation is contrasted with the scenario where all qubits couple to a single bosonic reservoir, in which case superdecoherence may in fact occur.We make some generalizations concerning the initial reservoir state, the details of which can be found in Appendix A. We use units where c = = k B = 1.

Single qubit
Consider a single qubit ('the system'), with an internal Hamiltonian H S = ∆J z , that is placed in a bosonic reservoir.Here ∆ is the level spacing and J z the spin-z operator.We work in the computational basis, where this operator is diagonal, and has eigenstates |1/2 and |−1/2 .The internal Hamiltonian of the reservoir is given by k a k the number operator of a bosonic mode with wave vector k.Here a k (a † k ) is the bosonic annihilation (creation) operator of the mode with wave vector k.The sum is over all k that are admitted by the reservoir.The set of ks that are admitted by the reservoir depends on the physical details of the reservoir.The reservoir couples to the qubit via the interaction term , with g k the coupling strength between the qubit and the mode with wave vector k.There are many explicit physical settings that may lead to this interaction term [26], but here, we do not assume such a specific setting.Since the only system operator in the interaction term is J z , H SB causes dephasing only.Putting all terms together, the dephasing model of a single qubit reads In this and the following sections, we assume that the overall system-reservoir state is a product state, ρ(0) ⊗ ρ B (0).Here ρ(0) (no subscript) is a general initial system state, and ρ B (0) is the initial reservoir state.The latter is assumed to be a product state of single-mode states, ρ B (0) = k ρ B,k (0), with ρ B,k (0) the initial state of the mode with wave vector k.The state ρ B,k (0) is assumed to be a displaced thermal state, that is, where α k is the displacement (which can be any complex number), N k the number operator, T k the (kdependent) temperature, Z the normalization, and D the displacement operator.(In App.A we show displacement is irrelevant in the dephasing process, so we do not give an expression for D here.)Possible ρ B,k (0) admitted by this parameterization include the regular single-mode thermal states (T k ≥ 0 and α k = 0), the coherent states (T k = 0, |α k | ≥ 0), and the vacuum state (T k = 0, α k = 0).We call a reservoir completely thermal if the overall initial reservoir state ρ B (0) equals the regular thermal density matrix with temperature T , that is, if ρ B (0) = e −ω k N k /T /Z .In our parameterization of initial reservoir states, this is the specific case where α k = 0 and T k = T for all k.Our form of the initial reservoir state is a generalization of that used in references [5-8, 11, 14, 17], where the assumption is that the initial reservoir state is completely thermal.
It can be shown that the absolute value of the i, jth entry (with i, j ∈ {−1/2, 1/2}) of the system density matrix, after time t, is given by where Γ i−j (t) is the decoherence function (refs. [6-8], 1 App.A).In the current model, dephasing is the only decoherence mechanism.Therefore the decoherence rate can be defined as 1/T 2 , were T 2 is the dephasing time, here defined as the smallest time t for which Γ i−j (t) = 1.
In general, the decoherence function only depends on the difference d = i − j. 2 It is given by Here γ d = |d| and Under the current assumptions on the initial reservoir state, the occupation number Nk of the mode  k is given by the Bose-Einstein distribution with (k-dependent) temperature T k . 3That is, This need not be an isotropic function on k-space.For the specific case of the completely thermal reservoir (i.e.T k = T and α k = 0 for all k), the occupation number is in fact isotropic, and depends on the mode energy only, We do not assume any particular dispersion relation, nor the reservoir to be completely thermal, unless stated otherwise.

Independent reservoirs
Now consider L copies of the system-reservoir combination described in the previous subsection.This setting is known as independent dephasing.The overall Hamiltonian reads H ind L = (H 1 ) ⊗L .This is depicted schematically in Fig. 2

(left).
We denote states in the computational basis of the L-qubit quantum register by |i ≡ |i 1 , . . ., i L .It can be shown that, under the evolution by H ind L , the absolute value of the (i, j)th entry of the system density matrix equals |ρ ij (t)| = e −Γ d (t) |ρ ij (0)|, with d the difference vector d = i − j and the decoherence function.Here, we have singled out the factor γ d = L =1 |d | for later reference.This factor is the only part of the decoherence function that depends on L, and it is at most proportional to L. Thus, for independent dephasing, the decoherence function scales at worst linearly with the system size, That is, the decoherence rate per qubit is at most constant in the system size. 3The initial state of the reservoir may still be a general thermal displaced state.Displacement of a mode does affect the expectation value of its number operator, but only the thermal part contributes to Γ d (t).See Appendix A for details.

A single reservoir
Now consider the situation where all qubits couple to a single reservoir, as is depicted schematically in Fig. (2).Again, k runs over all wave vectors that are supported by the reservoir.Now, the coupling constant g k depends on both the wave vector and the qubit location.If the reservoir consist of plane-wave modes, g k = g k e ik•r .For single-reservoir dephasing, it can be shown that (Refs.[6][7][8], App.A) the density matrix equals as before, but now where r m := r − r m is the vector pointing from the location of qubit to that of qubit m.In contrast to the situation of independent dephasing, γ d (k) now depends on k and contains a double sum over the qubit indices.The summand of γ d (k) can at most equal unity, which is attained, for example, if d = 1 for all , and k = 0. Thus, if indeed k = 0 is admitted by the reservoir, at worst.The possibility of quadratic, rather than linear scaling of the decoherence function with L is called superdecoherence.The decoherence rate per qubit can thus scale with the system size, which is problematic for error correction [27,28].

The continuum limit
We may write Eq. ( 7) in a more meaningful form, starting by introducing D = k δ(k − k ), so that we may replace the sum by an integral, where d is the dimension of the reservoir.Here D(k) is a density of states on k-space, currently describing a discrete set of modes.Note that D is unbounded at those modes, and vanishes elsewhere.In the continuum limit, the peaks merge into a bounded and continuous density of states on k-space.We then have where D is unbounded for discrete reservoirs, and bounded in the continuum limit.In the continuum limit, Nk becomes an occupation density rather than an occupation number.Note D(k) is different from the usual density of states, because the latter is a function of frequency only.For the electromagnetic field in free space, without boundary conditions, D is proportional to a constant with length dimension d.Equation ( 9) is the most general form of the decoherence function in the dephasing model because it can describe both discrete and continuous reservoirs.We will work with this form from now on.
One feature of Eq. ( 9) (and the preceding, less general forms) is that we can easily separate the vacuum contributions form those that are due to reservoir excitations.That is, we may write with where For the dephasing susceptibility to be well-defined, the integral in Eq. ( 9) has to converge.This is guaranteed by a high frequency cutoff.Physically, this arises because, as a function of ω k , either D goes to zero, or the coupling strength g k goes to zero, or a combination of both.Here, we assume that after some cutoff frequency ω c , the product At this point, this cutoff does not impose any restriction on the physical systems described because ω c can be arbitrarily large.Even in continuous reservoirs, it is possible in theory that a single mode k is excited, but no modes in its neighborhood (in k-space).Then, the occupation density is unbounded at that mode, Nk ∝ δ(k − k ).We call a reservoir continuous if, in contrast, both D(k)|g k | 2 and Nk are bounded functions of k.
A common assumption [6][7][8], that we will only make occasionally, is that D(k)|g k | 2 and Nk are isotropic, and that the reservoir dispersion relation is linear.For a linear dispersion relation, ω k = v|k| for some constant v. Working in units where v = 1 for notational convenience, we may then transform to spherical coordinates and write Here Ω is the d − 1 dimensional solid angle, and called the spectral density of the reservoir.A common form is [5-9, 29-33, 14] with α d a constant with length dimension d − 1, and ω c the cutoff frequency.This expression is often extended to include even non-integer d, which may be encountered in reservoirs with fractal properties [9].Depending on the dimension, these reservoirs are called subohmic (d < 1), Ohmic (d = 1), or superohmic (d > 1).In this manuscript, we do not assume isotropy, unless stated otherwise, and we will manly work with the general form of the decoherence function [Eq.( 9)].
In the following sections, we study the qualitative system-size scaling of the decoherence function.For completeness, however, in Appendix C we show explicit solution for Γ (vac) L , and derive simplified approximate solutions in the regimes t 1 and t → ∞.

Dephasing susceptibility
In this section, we identify γ d (k) as an important physical quantity and derive some of its properties, especially regarding its system size dependence.Namely, γ d (k) is determined solely by the system, and it is the only part of the decoherence function that depends on the system.So it fully captures the influence of the system on the decoherence function.The function γ d (k) weighs the severity of the influence of the mode k if this mode was to be 'offered' by the reservoir, and depends on the system geometry and the index (i, j).We call it the dephasing susceptibility of the reservoir.In Appendix D, we show how this susceptibility relates to the dynamical fidelity susceptibility of decoherence-free subspaces, which we introduced in previous work [34].
To illustrate the qualitative behavior of the dephasing susceptibility, we first consider the array model.It consists of a linear array of L noninteracting qubits with spacing a that couple to a single reservoir with dimension d = 1.Two system states we consider are Dephasing susceptibility Both states are of the form (|i + |j )/ √ 2, and thus have only a single nonzero matrix element in the upper right triangle of their density matrix.That is, in the computational basis, and similarly for the density matrix associated with GHZ .The difference vectors d = i − j belonging to these off-diagonal matrix elements are Thus, with Eq. ( 8), we find for L even.Here, we write γ GHZ instead of γ d GHZ for conciseness, and similarly for γ GHZ .Plots of γ GHZ (k) for various L can be found in Fig. 3. 4 These closed form formulas are ill-defined when the denominator vanishes.The original form [Eq. ( 8)] does not have this anomaly.It is to be understood that at these points, the closed form formulas are determined by their limit values.Then the resulting functions are smooth.
The dephasing susceptibility γ GHZ is dominated by the peak at ak = π, whose height is L 2 .Depending on the reservoir, this may result in superdecoherence.From Eq. ( 7), we see that if the reservoir is discrete and supports the mode ak = π, the decoherence function scales as L 2 , even in the vacuum.We stress that, as shown by this simple example, superdecoherence is possible even when the coupling constants g k depend on the the qubit location.Hence permutation symmetry of the Hamiltonian is not a prerequisite for superdecoherence.
Consider the two points around the peak where γ GHZ = 0.The previous equations about the minima show that the distance between these points equals ∆k = 4π/(aL).Thus, the bandwidth of modes the off-diagonal matrix element of GHZ is most susceptible to scales inversely with the system size.This decreasing bandwidth is shown by the dephasing susceptibility in general.This is because, mathematically, γ d (k) is the spectral density of the difference vector d.That is, Here the sine vanishes because it is antisymmetric under exchange of and m.If the qubits are placed on a lattice, the dephasing susceptibility is periodic in k.
The summand in the original definition of γ d (k) [Eq.( 8)] is at most unity.This is achieved, for example, when d n = 1 for all n ∈ {1, . . ., L} and k = 0. Thus, Nevertheless, the integral of the dephasing susceptibility over one reciprocal unit cell C is bounded by 2πL/V , where V is the volume one real-space unit cell.This follows directly from the fact that the dephasing susceptibility is the spectral density function of d, and Parseval's theorem, This shows that if the dephasing susceptibility has a peak of height L 2 , the width of that peak must scale as 1/L.As we show in the following section, this relation causes a mitigation of the dephasing process, causing the absence of superdecoherence in continuous reservoirs.
A related question about the dephasing susceptibility is how large γ i−j (k) is typically if we fix k and L and vary (i, j).In Appendix B, we show that the distribution of γ i−j (k) over (i, j) is approximated by a Gaussian, with a standard deviation that is at most L/(2π).This means γ i−j (k) is typically on the order of L and that there are few i − j such that γ i−j (k) ≈ L 2 .

Asymptotic system size scaling
In this section, we derive our main results, which are upper bounds on the system size scaling of the decoherence function.An important quantity herein is the decoherence function, because this is the only factor in the integrand of the decoherence function that depends on the system size.In turn, the dephasing susceptibility depends on L because L is the length of the vector d.The exact scaling of Γ d with L depends on how entries are added to d as the L is increased.In principle, this can be done according to any prescription.
For example, we could consider the dephasing associated with d, where d increases in length by adding a random number for every qubit we add.A more physically relevant situation, is for example to consider the coherence of the state |GHZ or GHZ , as a function of the system size.The results in this section hold for any description, unless stated otherwise, but some descriptions may arise more naturally than others.
To tidy up notation, and to highlight L dependence, we will now write γ L (k) instead of γ d (k) and likewise Γ L (t) instead of Γ d (t).At the same time, we use γ GHZ (ω) and γ GHZ (ω) for the dephasing susceptibilities of the off-diagonal matrix elements of |GHZ and GHZ , respectively.Likewise, we write Γ GHZ and Γ GHZ .
The starting point of our derivation is the most general form of the decoherence function, in which the vacuum contributions are separated from the excitation contributions [see Eq. ( 10)].Both contributions are of the form of Eq. (11).Assume t = t 0 is fixed.The mathematical property that the integral of γ L over one reciprocal unit cell scales linearly with the number of qubits [Eq.(25)] ensures that also Γ (vac/ex) L scales linearly with the number of qubits, provided that ξ (vac/ex) is bounded.This is shown as follows.Assume ξ (vac/ex) (t 0 , k) is bounded.The integral Γ (vac/ex) L equals a sum of integrals, where each domain of integration is one reciprocal unit cell C, Each term is upper bounded by the integral of γ L (k) over a single reciprocal unit cell after the integral is rescaled by the maximum of ξ (vac/ex) (t 0 , k) on that unit cell, By Eq. ( 25), The high-frequency cutoff [Eq.( 14)] ensures the sum converges, no matter the value of the cutoff ω c .Thus, we obtain the main mathematical result of this paper: if ξ (vac/ex) (t 0 , k) is bounded, then The relevant physical question then, is when (12)].The temporal factor τ (t 0 , ω k ) is a bounded function of ω k for every t 0 .The remaining factor D(k)|g k | 2 is bounded for continuous reservoirs (see Sec. 2).Therefore, in continuous reservoirs, This says that in continuous reservoirs, vacuum fluctuations cannot cause superdecoherence.Now consider the excitation contribution and Nk are bounded.By Eq. ( 26) we have, in that case, Together with Eq. ( 27), this says there is no superdecoherence in continuous reservoirs.
Conversely, we can consider the situations in which ξ (vac/ex) is unbounded.First, consider ξ (vac) .It is unbounded if the reservoir is discrete, that is, if Even though, in this case, the conditions of Eq. ( 26) are not satisfied, this does not lead to superdecoherence per se.It is clear that Γ (vac) L scales superlinearly with L only when one of the modes in D coincides exactly with a mode to which the matrix element is superlinearly susceptible.This is also illustrated by Fig. 3 and Eq. ( 7): there is superdecoherence in the array model when the state is GHZ , and π/a ∈ D. If, in the array model, π/a / ∈ D, but instead π/a + δ ∈ D, with 0 < |δ| 1, there is no superdecoherence.Note this in an asymptotic statement, and that, in the latter situation (π/a / ∈ D, π/a + δ ∈ D), and for finite L, linear scaling of Γ L with L only occurs after 1/L is approximately smaller than |δ|.Further discussion on finite-size effects can be found in Sec. 5.
Secondly, consider ξ (ex) .It is unbounded if the reservoir is discrete, like in the previous paragraph.It may additionally be unbounded if Nk is unbounded.This happens when a mode k is excited but no modes in its neighborhood are excited.Again, this does not need to lead to superdecoherence per se.It is only when k coincides exactly with a mode the matrix element is highly susceptible to that superlinear scaling of Γ (ex) i−j is possible.

Completely thermal reservoirs
If the reservoir has a continuous spectrum, and the initial reservoir state is completely thermal, ξ (ex) (t 0 , k) [Eq.(13)] is possibly unbounded because N th ω [Eq.( 5)] has an algebraic divergence at the origin.In this subsection, we show this nevertheless does not lead to superdecoherence (i.e. it does not lead to quadratic scaling of the decoherence function with the system size).However, superlinear scaling may be obtained, but only in subohmic reservoirs at nonzero temperature.This is shown as follows.Consider ξ (ex) with Nk = N (th) ω .Note that τ [Eq.( 3)] is constant to first order at the origin, so that it cannot contribute to the divergence.Thus, ξ (ex) is bounded near the origin if D(k)|g k | 2 goes to zero fast enough near the origin.In the remainder of this subsection, we will assume the isotropic setting of Eq. ( 15), with J(ω) as in Eq. ( 16).Then the condition for bounded ξ (ex) becomes d ≥ 1.This means there is no superlinear scaling of the decoherence function for Ohmic (d = 1) and superohmic (d > 1) continuous thermal reservoirs.

Subohmic thermal reservoirs
For subohmic reservoirs (d < 1), ξ (ex) in fact diverges at the origin.Here, we show how this can only lead to superlinear scaling of Γ (ex) L with L when γ L scales superlinearly with L near the origin.Even if γ L scales superlinearly with L near the origin, quadratic scaling may be approached, but not attained.
Let us first single out the divergence near the origin by defining Γ (ex) L = I L + J L , with and J L the remainder of the integral.Note ξ (ex) (t 0 , ω) now contains the thermal occupation density explicitly.
The integral J L is O(L) because, on the domain of integration, ξ (ex) is bounded [also see Eq. ( 26)].We now turn to I L .Given an ε, there exists a constant C 2 such that ξ (ex) ≤ C 2 ω d−1 on (0, ε].Thus, Since ω d−1 is monotonically decreasing, the largest possible value of I L occurs when γ L (ω) is peaked at low ω.Herein it is constrained by γ L (ω) ≤ L 2 [see Eq. ( 24)] and (25)].Under these constraints I L is largest when γ L (ω) is a bump function, where the bump height is L 2 , the left of the bump coincides with the origin, and the width of the bump is 2π/(V L).Therefore, Thus, quadratic scaling of Γ (ex) L (t 0 ), and thereby quadratic scaling of Γ L (t 0 ), cannot be obtained in subohmic continuous thermal reservoirs.
To approach superlinear scaling, it is essential that a superlinear peak of γ L (ω) must be able to approach the origin arbitrarily closely as a function of L. In fact, if, on the contrary, there is a δ > 0 such that γ L (ω) = O(L) for all 0 ≤ x ≤ δ, then I L = O(L).This is shown as follows.Assume there is a δ > 0 such that 0 < δ < ε and γ L (ω) = O(L) for all 0 ≤ ω ≤ δ.Then because ω d−1 is finite on [δ, ε], and because there is a constant C 3 such that γ L (ω) ≤ C 3 L for all [0, δ), we have An example in which this occurs is the array model, in the specific case that the dephasing susceptibility is given by γ GHZ (ω) [Eq.(21)].To show this, let δ = π/(2a).Then cos 2 (aω/2) ≥ cos 2 (π/4) ≥ 1/2 for all 0 ≤ ω ≤ δ, and thus γ GHZ ≤ 2 sin 2 (aωL/2) ≤ 2 for all 0 ≤ ω ≤ δ.This means that the off-diagonal matrix element of GHZ does not suffer from superdecoherence in subohmic thermal reservoirs, despite the fact that ξ (ex) is unbounded at the origin.
The result Eq. ( 30) is an upper bound, so the question remains if it may be attained.This is not clear a priori because a dephasing susceptibility cannot attain the form of a bump function as in the proof.This is because it is the spectral density of a vector with a finite number of elements [Eq.(23)].We now show by explicit construction that the upper bound may also be attained.This construction is in the subohmic version of the array model (Sec.3), with dephasing susceptibility γ GHZ [Eq.( 21)].Roughly speaking, our strategy is to show that γ GHZ (ω) is a close enough approximation of the bump function.There are two main steps.The first is to show that for all Consider the expansion of γ GHZ (ω) in aω around aω = 0. Using the original definition of the dephasing susceptibility [Eq.( 8)], we have γ GHZ (ω) = j=2,4,... c j (aω) j , with The radius of convergence of the expansion is infinite.
Using the coefficients, we have The last step can be checked most easily by plotting both functions.The last inequality holds specifically for aω ≤ 1/(2L 2 ).After substitution, we have, therefore, that γ GHZ (ω) ≤ 1/L 2 for all 0 ≤ ω ≤ 1/(2aL 2 ).The second step is to show that Eq. ( 31) enables us to approach quadratic scaling of I L with L arbitrary closely.First, note that, from Eq. ( 28), for 1/(2aL 2 ) < ε.There exists an L 0 and a constant C 4 such that for all L > L 0 , ξ (ex) (t 0 , ω) ≥ C 4 ω d−1 on the entire domain of integration.Informally, this means that there is a C 4 such that, close enough to the origin, ξ (ex) (t 0 , ω) ≥ C 4 ω d−1 .Thus, for this C 4 , Now using Eq. ( 31), this leads to Here the meaning of Ω(x) is similar to that of O(x), but Ω(x) refers to a lower instead of an upper bound. 5hus, in the array model with a subohmic continuous thermal reservoirs, quadratic scaling of Γ (ex) L , and thereby Γ L , may be approached arbitrarily closely by the off-diagonal matrix element of |GHZ .

Infinite time limit
In our discussion of the system size scaling until now, we assumed the time t to be fixed.Here we consider the infinite time limit of the isotropic case [Eq.( 15)], with J(ω) as given in Eq. ( 16).In the following, we no longer assume d < 1 and Nω = N th ω as in the previous subsection.With ∂ t τ (t, ω) = sin(ωt)/ω [cf.Eq. ( 3)], Thus, the infinite time behavior of Γ L (t) depends only on the integrand at the origin, which is always nonnegative.If the limit on the right hand side of Eq. ( 32) is positive, Γ L (t) keeps growing indefinitely as a function of t.If, on the other hand, this limit is zero, Γ L (t) increases at most sublinearly with t as t goes to infinity.We call this a quasi-plateau, which naturally includes proper plateaus.These proper plateaus are also referred to as incomplete dephasing [40] or coherence trapping [14].In Appendix C.3 we compute the height of the proper plateaus of Γ (vac) L explicitly in the array model.
As an example, we can read off that for γ GHZ , which is O(ω 2 ) as ω goes to zero [see Eq. ( 21)], in a completely thermal reservoir [Eq.( 5)], a (quasi-)plateau is reached for all T ≥ 0 and d ≥ 0. From Eq. (32) alone we cannot defer anything about the height of the (quasi-)plateau.

Finite-size effects
In the previous section, we focused on the asymptotic system size scaling of the decoherence function.We saw that, in that case, a sharp delineation could be placed between cases of superlinear and linear scaling.For finite system sizes, the situation becomes less clear.This is because the decoherence function may scale quadratically up to some potentially large system size L 0 , and show linear scaling only for L > L 0 .
Even though the main goal of this manuscript is to investigate the asymptotic scaling of the decoherence function with the system size, we discuss some finitesize effects in this section.

Time
Assume, for simplicity, a linear, isotropic dispersion relation, ω k = ω = |k|, in units where the proportionality constant equals unity.Consider the temporal factor τ (t, ω k ) [Eq. ( 3)] as a function of k.The function is peaked at the origin, with height t 2 /2.Away from the origin, it drops to zero at |k| = 2π/t and remains small afterwards [O(1/|k| 2 )].Thus, for large t, τ gives large weight to wave vectors with a length below 2π/t, and ever smaller weight to wave vectors with a length above 2π/t.
In Sec. 3, we showed that, if γ L (k), as a function of k, has a peak of height L 2 , the support of that peak must scale as 1/L.This effect causes the absence of superdecoherence in continuous reservoirs.However, if this peak is located at the origin, but t is such that the peak of τ (t, ω k ) is much narrower than that of γ L (k), we have that γ L (k) is approximately constant on the interval where τ (t, ω k ) is non-negligible.Thus, the reducing bandwidth of γ L (k) is only guaranteed to have an effect if aL t. ( Therefore, the actual scaling of Γ L as a function of L may approach its asymptotic scaling only at times small compared to the system size.This seems to form an important caveat to our asymptotic results.However, it only applies in special cases.Firstly, γ L (ω) needs to scale superlinearly as a function of L near the origin, which is rarely the case [see App.B].Secondly, even if γ L (ω) scales superlinearly, the remaining factors ξ (ex) and ξ (vac)  may kill the entire integrand around the origin [see Eq. (11)], for example when ξ (vac) (t, ω) = O(ω d ) and ξ (ex) = O(ω d ) as ω → 0, with d ≥ 1.Then for every γ L (ω) that scales superlinearly at the origin and fixed L 0 , there is a continuous crossover from superlinear to linear behavior in L around L 0 as a function of d.See Fig. 4 for two concrete examples.

Peaked occupation density
A similar finite-size effect occurs if the occupation density has a peak that coincides with a superlinear peak of the dephasing susceptibility.To separate this effect from the one in the previous subsection, consider as an example the state GHZ , in the array model, with a Gaussian occupation density Nω .The Gaussian has mean π/a, variance σ and an integrated number of bosons Ntot := the 1/L bandwidth of the dephasing susceptibility has an effect only after the peak of the dephasing susceptibility becomes narrower than that of the occupation density.That is, we only expect linear scaling of the decoherence function for 2π aL < σ.
See Fig. 5 for plots of the leading order in time of the decoherence function Γ GHZ .Again, in many situations, the effect discussed in this subsection does not have significant effects.Firstly, note that the integral in Eq. ( 15) is over infinitely many periods of the dephasing susceptibility.In the example above, the peak of the occupation density occurs only at a single frequency.In this case, the effect described in this subsection will thus only occur at one of the periods of γ GHZ (ω).Secondly, the center of the peak of the occupation density has to coincide exactly with the peak of the dephasing susceptibility.
The latter situation, where there is a single peak in the occupation density that overlaps exactly with the peak in the dephasing susceptibility, occurred in the ion-trap experiment by Monz et al., where the model of single-reservoir dephasing is applicable [20].The state |GHZ was prepared in a semistatic, semi-uniform magnetic field, which was produced with a Helmholtz coil.Fluctuations of the field, caused by current fluctuations in the coil, excited long-wavelength modes (with k ≈ 0).The dephasing susceptibility of the off-diagonal matrix element of |GHZ scales as L 2 at k = 0 (see Fig. 3), and exactly the modes k ≈ 0 where heavily excited in the experiment.If modes were excited away from In units where a = 1, the occupation density Nω is taken to be a Gaussian, with mean ω0 = π, standard deviation σ and integrated occupation Ntot = 10.We see linear scaling with L is obtained for 1/L σ.In the limit σ → 0, the occupation density becomes unbounded, and consequently, it is only in the limit that the decoherence function scales as L 2 for all L. The plot uses an analytic solution of Eq. ( 15), with where Θ is the step function.This form of the spectral density is chosen to accentuate the finite-size effects; the function γ GHZ (ω) has peaks at ω = π, π + 2π, . .., whereas, in this example, Nω only has a peak at ω0 = π.Finite-size effects only occur at places where the two peaks overlap, and including frequencies higher than 2π into J(ω) means including more effects that scale with L rather than L 2 .
the origin, there would not have been superdecoherence.This follows from the explicit form of the dephasing susceptibility of the off-diagonal matrix element of |GHZ (also see Fig. 3).Furthermore, even if modes near the origin were excited, almost any matrix element other than the off-diagonal element of |GHZ would not have suffered superdecoherence.This statement is shown in more detail in Appendix B.

Conclusion and outlook
In this paper, we studied superdecoherence in the model of single-reservoir dephasing for asymptotic system sizes.We have shown that if the density of modes in k-space, D(k), and the occupation density Nk are bounded, superdecoherence is not possible.This is because if there is a k such that the dephasing susceptibility scales quadratically with the system size, γ L (k) ∝ L 2 , the support of this peak in the dephasing susceptibility always scales inversely with L.
Superdecoherence may thus only be obtained if D(k) or Nk is unbounded.The former happens if the reservoir supports only a discrete set of modes.The latter happens if the reservoir supports a continuum of modes, but only perfectly isolated modes are excited.In both cases, the unbounded point must coincide exactly with the mode for which the dephasing susceptibility γ L (k) scales quadratically.
For completely thermal continuous reservoirs, the occupation density N th ω diverges algebraically at the origin.Nevertheless, even in this case, the decoherence function scales at most linearly with the system size.There is one exception.This is the subohmic continuous thermal reservoir with nonzero temperature, where, furthermore, the dephasing susceptibility must scale superlinearly (which includes quadratic scaling) at low frequencies.In this case, the decoherence function may approach, but not attain, quadratic scaling with the system size.
All effects discussed int this manuscript can in principle be observed experimentally.One could compare the effects of narrow-band versus broadband noise at a frequency to which the system is highly susceptible.Or, the system could be placed in a high-Q cavity, in the vacuum state, that supports exactly the mode the system is highly susceptible to.This is to be compared to a situation where the cavity is slightly longer.
Other applications lie in quantum metrology, in which superdecoherence can be used as a means of enhancing sensitivity.
In this context, it is well-known that the GHZ state is highly susceptible to long-wavelength modes [41,42].The dephasing susceptibility, as defined in this paper, offers an effective way to extend metrology to other states and wavelengths; any state for which there is a ω 0 such that γ L (ω 0 ) ∝ L 2 is suitable for quantum metrology of the mode with wavelength ω 0 .An example is a linear array of qubits, with spacing a, in the state GHZ [Eq.( 17)].This system is highly susceptible to the staggered mode ω 0 = π/a.The dephasing susceptibility shows the added benefit that with increasing system size, the array becomes less sensitive to frequencies other than ω 0 .

A Spin-boson dephasing for arbitrary reservoir states
In this section, we generalize the decoherence function of single-reservoir dephasing, as it is found in Refs.[6][7][8], to more general reservoir states.We pay specific attention to Gaussian states, Gaussian product states, and a product of displaced thermal states.The latter solution is included in the main text as Eqs.( 2) and (7).
For a completely general reservoir state, it can be shown that [6][7][8] where χ is the characteristic function of the reservoir state, Equations ( 34) and ( 35) hold for the interaction as well as the Schrödinger picture density operators.The set D contains all wave vectors that are supported by the reservoir.The argument of the characteristic function, λ ∈ C |D| , depends on the matrix index (i, j) and the time t, but the notation of this dependence is suppressed.The kth entry of λ is given by [6][7][8] with d the Fourier transform of d = i−j [see Eq. ( 16)].
The exponent in Eq. (36) stems from the internal time evolution of the reservoir.Equations ( 34) and (35) give the most general form of the absolute value of the time evolved reduced density matrix in the singlereservoir dephasing model.For the class of Gaussian states [43,44], the absolute value of the characteristic function is given by This is the most general form of the decoherence function.Here, writing λ ki as λ i for short, The 2|D| × 2|D| matrix σ is the covariance matrix, The expectation value is with respect to the Gaussian initial reservoir state ρ B (0).The vector R is defined by RT = (q 1 , p1 , . . .q|D| , p|D| ), with {•, •} the anti-commutator.To avoid confusion about operators versus numbers, in this section we write operators (and vectors containing operators) with hats, as opposed to in the main text.The quadrature operators qm and pm , in turn, are defined by To obtain the decoherence function as a function of time and the density matrix index (i, j), the expression for λ [Eq.( 36)] has to be inserted into Eq.( 37).We may consider various simplifications of the decoherence function as it is given in Eq. (37).If the reservoir modes are unentangled, that is, ρ B (0) = k ρ B,k (0) with all ρ B,k (0) Gaussian, the covariance matrix is block-diagonal.Each block corresponds to a 2 × 2 single-mode covariance matrix, which we denote by σ k .In this case, we may write where, by Eq. ( 38), the entries of the single-mode covariance matrix read If the mode k is initially in the thermal state, with temperature T k , its density matrix reads ρ B,k (0) ∝ e −ω k â † k âk /T k .In this case, the single-mode covariance matrix is diagonal, with Nk the occupation number [Eq. ( 4)].A special thermal state is the vacuum, where Nk = 0.
If the reservoir modes are unentangled, and every mode is thermally excited with its own temperature, we have from combining Eqs. ( 39) and (40) that Inserting the equation for λ k [Eq.( 36)], we obtain Eq. (7).
In general, a single-mode Gaussian state can also be represented as a squeezed and displaced thermal state [43], Here r is the squeezing magnitude, and ϕ the squeezing angle.Note these expressions are invariant under displacement.Therefore, the decoherence function of a displaced thermal state is equal to Eq. ( 41), with N the regular Bose-Einstein distribution.Squeezing, on the other hand, does affect the covariance matrix, and would alter Eq. ( 41) straightforwardly.In the main text, we assume for simplicity that the reservoir modes are not squeezed.
Displaced vacuum states are precisely the coherent states.Thus, even if a reservoir mode is in a highly excited coherent state, this mode does not contribute more to the dephasing process than the same mode in the vacuum state would have done.A mixture of coherent states does lead to extra dephasing.However, the only mixture that can be described in the Gaussian state formalism is the thermal state.
To summarize, in single-reservoir dephasing, the decoherence process of the system is completely determined by the reservoir characteristic function; The argument of the characteristic function, λ, is a complex vector which depends on the matrix index (i, j) and time [Eq.(36)].For completely general reservoir states, the characteristic function is given by Eq. (35).For general Gaussian reservoir states, | χ(λ)| = e −Γ(λ) , with Γ(λ) = 1  2 Λ T σΛ the decoherence function, generalized to Gaussian states [Eq.(37)].In case the reservoir modes are unentangled, Γ(λ) may be written using a single sum over k [Eq.(39)].If, furthermore, each of these modes is a (possibly) displaced thermal state, the decoherence function simplifies further to Eq. (41).It is this form of the decoherence function that we use in the main text.Interestingly, displacing a reservoir state has no effect on the dephasing process.For example, this means that it does not matter for the dephasing process if either a mode is in a highly excited coherent state or the vacuum state.

B Typical values of the dephasing susceptibility
Here, we ask the question if there many i−j such that γ i−j (k) ≈ L 2 , given fixed values for L and k.We show this is not the case: as we go over all (i, j), the values of γ i−j (k) are distributed according to a Gaussian that has a standard deviation that is at most L/(2π).This means that, for a random (i, j), γ i−j (k) is typically on the order of L/(2π), or less.
To show this, fix L and k, and consider the function Consider the frequency distribution of this function.This is a table that, per possible value D 0 of D ij , shows the number of inputs (i, j) such that D ij = D 0 .To obtain this distribution, we see D as the distance from the origin of a random walker on the complex plane.The walker takes L steps, where the th step is given by d e ik•r , with d = i − j.For the th step, the walker has a probability 1/2 to make no step at all, a probability of 1/4 to take the step +e ik•r , and a probability of 1/4 to take the step −e ik•r .After L steps, the walker is a distance D ij away from the origin of the complex plane.
Naturally, the variance in the distances form the origin is largest if the walker is restricted to move on a single line, which happens if k = 0. Let us therefore put k = 0, keeping in mind that, at worst, we are overestimating the variance of D ij for other values of k.For a 1D random walker that can take the steps +1 and −1 with equal probability, it is well-known that, after L steps, the distribution of distances from the origin is well approximated by a Gaussian with standard deviation 2L/π.In our situation, half of the time the 1D walker does not take a step at all.Therefore, the distribution of D will be approximated by a Gaussian with variance L/(2π).Since γ i−j = (D ij ) 2 , the distribution of γ i−j over (i, j) is approximated by a Gaussian with standard deviation L/(2π).This means that for fixed L and k, and given a random (i, j), the decoherence function is, at most, typically on the order of L/(2π).Additionally, it means that, if we are given a random (i, j), where also the dimension L of i and j is random but equal, the probability that γ i−j ≥ κL 2 goes to zero as L goes to infinity, for all k and κ > 0.

C Explicit expressions for the vacuum contribution
Here, we derive the explicit solution of the vacuum part of the decoherence function, Γ (vac) d (t), in the array model of Sec. 3. The assumptions are that the qubits form a linear array with spacing a, coupling to a one-dimensional reservoir via the single-reservoir dephasing Hamiltonian.The spectral density of the reservoir is assumed to be given by Eq. ( 16).
In principle, in the array model, d = 1, but we will analytically extend our solutions to arbitrary d.After absorbing the integral over the solid angle, which in d = 1 dimensions gives a factor of 2, into α d , the vacuum decoherence function reads with In this section, we derive the full solution of Eq. (42).Additionally, we derive simplified approxi-mate solutions for the limits of infinitesimal and infinite time.For infinitesimal times, we find where Γ is the regular gamma function Γ(j + 1) = j!, not to be confused with the decoherence function.In the infinite time limit, Γ

C.1 General solution
We start by rewriting γ d (ω) as where Written this way, γ d (ω) is the cosine transform of f dr .
For later reference, we note that for the states |GHZ and GHZ , where ζ = 1 for |GHZ and ζ = −1 for GHZ .Going back to the general case, we have from Eqs. (42) and (43), and with This integral is solved using standard identities for Gaussian integrals.For d > 0, d = 1, Here c.c. stands for the complex conjugate of the preceding term, and with i the imaginary unit.For d = 1, We now have Γ

C.2 Infinitesimal time limit
The leading, second order in time of the integral in Eq. ( 46) equals Solving this integral, we obtain, for r = 0, For r > 0, we find These two solutions hold for all d > 0. Up to a factor α d d 2 , the first term of the decoherence function in Eq. ( 45) is given by Eq. (50).For the remaining terms, with r > 0, note that |f dr | ≤ 2f d0 = 2 d 2 .Thus, Note |I r (t)| is proportional to Thus, we obtain  The time interval between these extrema is equal to the time required by the mode to travel a distance a.There are no extrema after the mode has had the time to travel the distance aL, which is the total length of the array.For |GHZ , the extrema are alternating local maxima and minima.For |GHZ , there are only local minima at these points.After the series of extrema, Γ GHZ and Γ GHZ reach the same plateaus, the height of which is given by Eq. (58).Plots for higher L, and odd L, show the same behavior.
There are extra factors of L hiding in the O(t/a) 4 term.We can disregard this L dependence because, in this subsection, we are interested in the limit of infinitesimal time.Then for any L there is a t/a 1 such that the second term in Eq. ( 52) is negligible.
Thus, for small times and d > 0, the final result is where E contains both the error from the α d 2 I 0 term, and all of the remaining terms in Eq. ( 45), Given an L and ω c a 1, the relative error 4 , is negligible for t small compared to 1/ω c and a.

C.3 Infinite time limit
If d > 1 and j = 0, the function (Q rj ) 1−d vanishes in the limit that t goes to infinity.For j = 0, on the other hand, (Q rj ) 1−d is time-independent and nonzero.Thus, from Eq. ( 42), for d > 1.Therefore, lim t→∞ Γ (vac) d (t) exists for d > 1, and its value can be found by substituting Eq. (54) into Eq.(45).The existence of this limit means the vacuum decoherence function always reaches a proper plateau for d > 1 (cf.Sec.4.2).
We now show the height of this plateau scales linearly with L for d ≥ 2, and, for these d, simplify the exact expression for the height of the plateau.(This result need not imply superlinear scaling of the height of the plateau for d < 2.) Firstly, With f d0 = d 2 , Eq. (54), and Q 00 = 1/(aω c ), the first term (r = 0) equals The remaining terms in Eq. ( 55) can be neglected.This is because they are upper bounded by For r ≥ 1, |Q r0 |>1, and d ≥ 2, we have

Thus, with
L−1 r=1 1/r 2 < 2, we have for the sum in Eq. (56) that In conclusion, we have for d ≥ 2, with relative error The latter is negligible for aω c 1. Note that Γ

D Dynamical fidelity susceptibility of dephasing
In Ref. [34], we studied the leading order effect of an perturbation to the system-reservoir coupling on states in a decoherence-free subspace.We defined this leading order as the dynamical fidelity susceptibility of decoherence-free subspaces.In this previous work we did not assume a particular Hamiltonian, in contrast to the current work, where we focus on the single-reservoir dephasing Hamiltonian.So, on the one hand, the previous work applies more generally.On the other hand, in the previous work we assumed a pure initial reservoir state, whereas in the current work, we assume a product of displaced thermal states.These two assumptions on the initial reservoir state describe different situations.It is only for the vacuum state that both assumptions are simultaneously satisfied.
In this appendix, we compute the dynamical fidelity susceptibility of decoherence-free subspaces in the specific case that the Hamiltonian is given by the pure single-reservoir dephasing Hamiltonian.This is done in two ways: first by using the general result from our previous work, and then by a more direct computation that does not require the results of the previous work.The expressions we find are identical.However, there is one subtle difference: in the former method, the expression contains the expectation value of the number operator, N k ϕ , with respect to the pure initial reservoir state |ϕ (in agreement with the assumptions in the previous work).In the latter method, the same expression contains instead the expectation value of the number operator with respect to the thermal part of a product of displaced thermal states Nk (in agreement with the assumptions in the current work).Naturally, the two expressions agree on the subset of initial reservoir states that satisfy both the assumptions on the reservoir states in the current and the previous work.
We now briefly introduce the result of our previous work, and sequentially use this result to derive the dynamical fidelity susceptibility of the singlereservoir spin-boson dephasing Hamiltonian.Consider a quantum register that is coupled to a reservoir.In general, the overall Hamiltonian is of the Here H S is the system Hamiltonian, H B the reservoir Hamiltonian, and H SB the interaction term.Assume that the initial overall state is a product state between the system and the reservoir, |Ψ = |ψ ⊗ |ϕ , with |ψ the initial system state and |ϕ the initial reservoir state.A decoherence-free subspace (DFS) is a subspace of the system's Hilbert space that does not entangle with the reservoir as |Ψ is evolved under the Hamiltonian H 0 , despite the coupling H SB .States in a DFS only experience the unitary evolution due to H S .That is, by definition of a DFS, we may write the system state, in the Schrödinger picture, as where H ε contains an extra interaction term εV ; with ε 1.In general V may be written as To compare the state (59) to (60), we computed the dynamical fidelity, which is defined as the fidelity between ρ Sch 0 (t) and ρ Sch ε (t).Because the former state remains pure, this fidelity has the simple form (62) In the interaction picture, where H 0 is the bare Hamiltonian, and εV the interaction term, it is straightforward to show that Here |ψ is the initial system state, and ρ Int ε (t) the reduced system state in the interaction picture, evolved in time through the interaction-picture time-evolution operator, with ε = 0.
The leading order of F in both t and ε is the second.We defined this leading order as the dynamical fidelity susceptibility of decoherence-free subspaces, This quantifies the leading order effect of the added system-reservoir coupling on states in a DFS.We have shown that where the expectation values are with respect to the initial system and initial reservoir state.This equation is not specific to the single-reservoir dephasing model, but holds in general.Using this general result, we can compute the dynamical fidelity susceptibility of the spin-boson singlereservoir dephasing model.After making the substitution g k → εg k , the Hamiltonian of the singlereservoir dephasing model [Eq.( 6)] is of the form of Eq. ( 61), with H SB = 0 and Note that, because H SB = 0, actually the entire Hilbert space of the system is a DFS.(66) Note that we have now obtained χ without ever solving for the reduced time evolution of the system state.This illustrates that Eq. ( 65) can be used to study superdecoherence in models that have not been solved.The condition for superdecoherence would read χ ∝ L 2 instead of Γ L ∝ L 2 .One caveat is that χ is a leading order in time.With Eq. (66) we have already obtained a result that could not be obtained by using the full time evolution of the spinboson single-reservoir dephasing model, because the latter relies on the assumption that the initial reservoir state is a displaced thermal state.Equation (66) holds for general pure initial reservoirs states, some of which cannot be described as a displaced thermal state (with T = 0).Furthermore, we can take the continuum limit of Eq. (66), just as we have done in the main text for Γ L , and show, in exactly the same way, that χ = O(L) in continuous reservoirs.This extends the results in the main text to include arbitrary pure reservoir states, albeit only for the leading order in time.
For the single-reservoir dephasing model, we actually have the full solution of the reduced system density matrix at hand, which means the dynamical fidelity susceptibility may be obtained by other means.We may solve for the interaction-picture states ρ Int 0 (t) and ρ Int ε (t), and compute derivatives of Eq. (62).In the main text, are interested in the absolute value of the system density matrix only.Here, however, we need the full solution, which, in the interactionpicture, reads [7] ρ Int ij (t) = e i[Θ ij (t)−Λ ij (t)] e −Γ i−j (t) ρ Int ij (0), where and This is not in the continuum limit.
Examples include (i, j) = (i GHZ , j GHZ ) and (i, j) = (i GHZ , j GHZ ).Computing the derivatives of F [see Eq. ( 64)], we obtain Note the similarity between Eqs. (66) and (69).The assumptions involved in deriving these two equations, however, are differexnt.In the former, the assumption is that the initial reservoir state is pure.In the latter, it is assumed that the reservoir state is a product of unentangled displaced thermal states (see Sec. 2).The only reservoir state for which these two assumptions coincide is the vacuum state.In this case, both expressions agree, and read The similarity of Eqs.(66) and (69) indicates that our assumptions on the initial reservoir state may be relaxed in both situations.

Figure 3 :
Figure 3: The dephasing susceptibility of the off-diagonal matrix element of the state |GHZ , for system sizes L = 2, 4, 6, 8.The peaks have hight L 2 and width ∼ 1/L.For the off-diagonal matrix element of the state |GHZ , the entire graph is translated in such a way that the peaks lie above k = 0.

∞−∞Figure 4 :
Figure 4: The decoherence function in the array model (see Sec. 3), as a function of L, for the off-diagonal matrix element of the state |GHZ (left) and |GHZ (right).In both plots, we use units where a = 1, and set t = 20, J(ω) = α d ω d e −ω/ωc , with α d = 1 (for both d = 1 and d = 2), ωc = 20, and Nω = 0. We have used the analytical expressions for the decoherence function that are derived in appendix C. (Left) For d = 1 the decoherence function increases quadratically initially, after which it scales (sub)linearly.For d = 2 there is no quadratic scaling, even for aL t. (Right) No superlinear scaling for any t, L and d (including d other than d = 1, 2, which are not shown).The lines for d = 2 in the left and right plot are similar, but not exactly equal.

Figure 5 :
Figure 5: The leading order in time of the decoherence function of the off-diagonal matrix element of |GHZ , as a function of the system size L. The setting is that of the array model (Sec.3), with d = 1.In units where a = 1, the occupation density Nω is taken to be a Gaussian, with mean ω0 = π, standard deviation σ and integrated occupation Ntot = 10.We see linear scaling with L is obtained for 1/L σ.In the limit σ → 0, the occupation density becomes unbounded, and consequently, it is only in the limit that the decoherence function scales as L 2 for all L. The plot uses an analytic solution of Eq. (15), withJ(ω) = α 1 2 ω[1 − Θ(ω − 2π)],where Θ is the step function.This form of the spectral density is chosen to accentuate the finite-size effects; the function γ GHZ (ω) has peaks at ω = π, π + 2π, . .., whereas, in this example, Nω only has a peak at ω0 = π.Finite-size effects only occur at places where the two peaks overlap, and including frequencies higher than 2π into J(ω) means including more effects that scale with L rather than L 2 .

.
reaches a plateau for all d > 1.For d ≥ 2, we show that the height this plateau equalslim t→∞ Γ (vac) d (t) ≈ α d d 2 Γ(d − 1)ω d−1 cThis result extends that of Sec.4.2 for the current, specific setting.Note that, because ||d|| 2 ≤ L, the decoherence function scales at most linearly with L in the limits of infinitesimal and infinite time, in accordance with the results in the main text.
in closed form, except for the sum over a single index in Eq. (45).Using this analytic solution, Γ GHZ and Γ GHZ are plotted in Fig.6.

Figure 6 :
Figure 6: The decoherence function of the off-diagonal matrix element of the density matrix of |GHZ (top) and |GHZ (bottom) in d = 2 dimensions, with ωc = 10/a and T = 0. We observe L − 1 sharp extrema at t/a = 1, 2, . . ., L − 1.The time interval between these extrema is equal to the time required by the mode to travel a distance a.There are no extrema after the mode has had the time to travel the distance aL, which is the total length of the array.For |GHZ , the extrema are alternating local maxima and minima.For |GHZ , there are only local minima at these points.After the series of extrema, Γ GHZ and Γ GHZ reach the same plateaus, the height of which is given by Eq. (58).Plots for higher L, and odd L, show the same behavior.
= O(L) even if this condition does not hold.

= 1 4 mk d d m g k g * km ( 1 + 2 N= 1 4 k
For system states of the form |ψ = (|i + |j )/ √ 2, we findχ = 1 4 mkk d d m g k g * k m (δ kk + 2 a † k a k ϕ ) + g k g k m a † k a † k ϕ + c.c. ,with d = i − j, and where c.c. stands for the complex conjugate of the preceding term only.If the initial reservoir state |ϕ is a product of number states, thenχ k ϕ ).Using g k = g k e ik•r and writing γ d (k) as the spectral density of d (see Sec. 3), we obtain the resultχ |g k | 2 γ d (k)(1 + 2 N k ϕ ).