A Supplemental material: Bounds on energy absorption in quantum systems with long-range interactions

# Bounds on energy absorption, and prethermalization in quantum systems with long-range interactions

## Abstract

Long-range interacting systems such as nitrogen vacancy centers in diamond and trapped ions serve as useful experimental setups to probe a range of nonequilibrium many-body phenomena. In particular, via driving, various effective Hamiltonians with physics potentially quite distinct from short-range systems can be realized. In this Letter, we derive general rigorous bounds on the linear response energy absorption rates of periodically driven systems of spins or fermions with long-range interactions that are sign-changing and fall off as with . We show that the disorder averaged energy absorption rate at high temperatures decays exponentially with the driving frequency. This strongly suggests the presence of a prethermal plateau in which dynamics is governed by an effective, static Hamiltonian for long times, and we provide numerical evidence to support such a statement. Our results are relevant for understanding timescales of both heating and hence new dynamical regimes described by effective Hamiltonians in such long-range systems.

\LetLtxMacro\ORIGselectlanguage\ORIGselectlanguage

english

Introduction.—Quantum many-body physics far from equilibrium is an exciting frontier of condensed matter physics. Recent experimental advances in designing and controlling well-isolated many-body systems, such as ultracold atoms Bloch et al. (2008), trapped ions Blatt and Roos (2012) and nitrogen-vacancy (NV)-centers in diamond Choi et al. (2017a), have enabled a controlled study of a range of nonequilibrium phenomena such as thermalization and many-body localization in closed many-body systems Nandkishore and Huse (2015); Abanin and Papić (2017); Schreiber et al. (2015); Chien et al. (2015); Smith et al. (2016); Kucsko et al. (2016); Choi et al. (2017a). Additionally, periodic driving has emerged as a useful tool to engineer interactions in these systems and create various effective Hamiltonians Eckardt (2017); Goldman and Dalibard (2014); Aidelsburger et al. (2013); Jotzu et al. (2014); Aidelsburger et al. (2015); Potirniche et al. (2016); Yao et al. (2013), even enabling novel nonequilibrium phases of matter such as time crystals Khemani et al. (2016); Else et al. (2016); Zhang et al. (2017); Choi et al. (2017b).

These experimental platforms generically fall into two classes distinguished by the nature of interactions: short-ranged (e.g. cold atoms), and long-ranged, power-law decaying (e.g. NV-centers and trapped ions). Most theoretical work has focused on systems with short-range interactions; in contrast, systems with long-range interactions have not been studied as extensively. The difference between the nature of the interactions is important as long-range interactions can lead to physical phenomena distinct from the short-range case. For example, it was argued that depending on the power-law exponent, long-ranged interactions can either destroy localization Anderson (1958); Burin (2006); Levitov (1990); Gutman et al. (2016) or reinstate MBL nonperturbatively Nandkishore and Sondhi (2017). Furthermore, the existence of a new, critical regime of time crystals was recently uncovered in a driven dipolar spin system Choi et al. (2017b); Ho et al. (2017). Thus, studying long-range systems and effective regimes derivable from them opens up avenues to observe new and interesting physics.

Unbounded heating due to a drive can however invalidate proposals to engineer new non-equilibrium phases and dynamical regimes, a problem which plagues both classes of systems Lazarides et al. (2014); Ponte et al. (2015); D’Alessio and Rigol (2014). It is therefore important to understand the heating timescales in driven many-body systems. Known rigorous results along this front such as exponentially slow heating Abanin et al. (2015) and prethermalization Abanin et al. (2017a, b); Mori et al. (2016); Kuwahara et al. (2016) at high driving frequency fundamentally rely on the local nature of the interactions and thus only apply to systems with interactions that are sufficiently short-ranged. Given strong experimental interests in periodically driven long-range systems, it is important to understand whether general constraints on heating in long-range interacting systems exist.

In this Letter, we derive general rigorous bounds on the energy absorption rate at high temperatures for driven systems of spins (or fermions) in spatial dimensions with long-range interactions whose coupling strengths are random and sign-changing, and which decay as with . We note that the interaction couplings in long-range systems with  are inevitably sign-changing to ensure thermodynamic stability Sig (). For example, an ensemble of NV-centers represents an effective spin system with long-range dipolar interactions () that are sign changing in nature. Moreover, systems of trapped ions can have  Blatt and Roos (2012).

Specifically, we prove that the disorder-averaged linear response energy absorption rate is exponentially suppressed at high driving frequencies (relative to the typical interaction strength), for both local and global driving. In order to prove our rigorous results, we do not employ the techniques involved in the previous proofs of prethermalization, which as noted before, rely on the local nature of the interactions Abanin et al. (2017b); Mori et al. (2016); Abanin et al. (2017a); Kuwahara et al. (2016). Rather, we use that it is their random nature that accords a cancellation of many terms in the response function at high temperatures. These results strongly suggest the presence of a prethermal plateau in which dynamics is governed by an effective, static Hamiltonian for long times, and we also provide numerical evidence to support such a statement.

Set-up and results.— We consider a many-body system of spins (or fermions) in dimensions, placed either on a regular lattice or randomly distributed in space with a short distance cut-off. The spins interact via long-range disordered interactions, so that the Hamiltonian is

 H=∑μJμrαμOμ+κ∑i→hi⋅→σi. (1)

Here the sum in is over all links between two sites with the distance between them , which without loss of generality is measured relative to the lattice spacing or short distance cut-off so that . is a generic two-body interaction where are Pauli-matrices at site and fixed real coefficients so that its spectral norm . The exponent characterizing the decay of interactions and is taken to satisfy , although we are mostly interested in the ‘truly’ long-range case for which : for such , the mean field on a given site due to the interactions is not absolutely convergent not (). We have also allowed for a potentially random on-site field , and assume that for all .

We take the interaction strengths to be independent and identical bounded random variables with zero mean, so that the following conditions on the -th moments are automatically satisfied:

 J(1)=0,J(2)≡J2,|J(n)|≤(λJ)n, (2)

for some . Lastly, we note that the sum of the inverse squares of the interactions emanating from site converges: because and , and further assume is uniformly bounded from above:

We focus on the case of a harmonic drive of the system at frequency and strength :

 H(t)=H+gcos(ωt)V, (3)

where is a sum of local terms with for simplicity and we also take (without loss of generality) . Here acts on a single site , which describes, for example, driving the system by a magnetic field .

Assuming the system is initially at thermal equilibrium with inverse temperature , the energy absorption rate is related to , the dissipative part of the linear response function, via . For a quantum system in a finite volume with discrete spectrum, the Lehmann representation of in the high temperature limit is

 σ(ω)=∑nmπβωZ0⟨n|V|m⟩⟨m|V|n⟩δ(En−Em−ω), (4)

where are energy eigenstates of and is the dimension of the Hilbert space. A related quantity was studied in Mukerjee et al. (2006). Note that Eq. (4) is actually a distribution and not a bona fide function – to state precise results, we have to integrate over a finite frequency window. For us, the object of interest is the disordered averaged high-frequency spectral weight of the response function where denotes disorder averaging over (and possibly ). We derive a bound for , establishing the main result of this paper:

 σ([ω])≤Nπβωe−Bω, (5)

where is the total number of spins in the system, and some constant that depends only on the typical two-body interaction strength up to numerical factors. This indicates that the energy absorption of long-range systems at high temperatures is exponentially suppressed at high frequencies.

Local drive with no on-site field.— We first prove the local version of Eq. (5) for a local drive acting only on site and without the on-site field so that .

We can rewrite as follows, since the delta function enforces energy conservation :

 σ([ω])=⟨∫∞ωdω′∑nmπβω′Z0∣∣ ∣∣f(p)n,mω′p∣∣ ∣∣2δ(En−Em−ω′)⟩

where is the matrix element of a -nested commutator of with . The integral over picks out a subset of eigenstates in the double sum so that . Then, extending the sum over over all eigenstates, using the resolution of the identity and also the cyclicity of the trace allows us to upper bound as (see SOM () for more details)

 σ([ω]) Missing or unrecognized delimiter for \left (6)

which serves as the starting point of our analysis for tra (). Writing as in Eq. (1) with , we can further bound using the triangle inequality as

 σ([ω]) ≤πβωZ0ω2p∑→μ∣∣Tr(V0[[[V0,Oμ1],Oμ2],⋯,Oμ2p])∣∣ ×∣∣⟨Jμ1Jμ2⋯Jμ2p⟩∣∣rαμ1rαμ2⋯rαμ2p (7)

where .

We note that there are two ways that the infinitely many terms in Eq. (7) can be controlled. First, because the interaction strength is sign changing so that , under disorder averaging, is nonzero only when each link that appears in the sum is at least paired. This then means that in the denominator, distances come with a power of at least two, i.e. , . Anticipating that we will later sum over one of the sites in , we see that this higher power guarantees convergence of the sum. Note that we could not have straightforwardly bounded Eq. (7) without the disorder average due to the non-absolutely convergent mean field strength. Second, the -nested commutator enforces connectivity of a contributing term: the links must be such that the operators have overlapping support and are fully connected. Terms with disconnected parts therefore do not contribute SRc ().

Let us first make use of the fact that each in the sum must be at least paired. A natural way to account for this is to consider how many times each label appears: this makes us consider all (unordered) integer partitions of , which we denote by Here is the length of the integer partition . In doing so, we have assumed that are all distinct, that is, . However, specifying the number of times appear is not sufficient – one also has to consider all orderings of these links in the nested commutator, as different orderings give distinct terms. Thus, let us introduce for a given integer partition , the set denoting the partitions of the list into parts, each of which is a subset of the list. We take each to be ordered by the smallest element appearing in each part. Thus, the length of the -th part of is an integer part of , such that is just some permutation of the integer partition egg (). With this information, we can then define as the -nested commutator in which the indices are distributed according to : appears at the positions in the -nested commutator dictated by the first part of , appears at the positions dictated by the second part of and so on (see eg- () and SOM () for more details).

Eq. (7) can then be written exactly as

 σ([ω])≤ πβωZ0ω2p∑s∈S(2p)′∑μ1,⋯,μl(s)∑q∈Q(s)f(q)× |⟨Jn1(q)μ1⟩⟨Jn2(q)μ2⟩⋯⟨Jnl(s)(q)μl(s)⟩|rαn1(q)μ1rαn2(q)μ2⋯rαnl(s)(q)μl(s), (8)

where the sum over runs over links which are all distinct (denoted by the prime). We can further upper bound the above expression using Eq. (2) by restricting the integer partitions to be over in where is the set of integer partitions where each integer part is only or i.e. , while simultaneously allowing the sum over the links to run unrestricted, so that

 σ([ω])≤ ApZ0∑s∈S′(2p)∑μ1⋯μl∑q∈Q(s)f(q)rαn1(q)μ1⋯rαnl(q)μl, (9)

where . The sum of Eq. (9) counts every term of Eq. (8), and more.

Next, we represent each link as a pair of sites , where denotes a new site appearing in the above sum at the th commutator, while denotes some site from the set . This parameterization arises because of connectivity: must be site , the point where the perturbation acts, while must be sites or , and so on. We also use and to replace each by , so that

 σ([ω])≤ Ap∑s∈S′(2p)∑i1,⋯,il∑a1∈{0,i1}⋯∑al(s)∈{0,i1,⋯,il(s)−1}× 1r2α0,i1r2αa1,i2⋯r2αal(s)−1,il(s)∑q∈Q(s)f(q)Z0. (10)

If we let be the number of times appear in an integer partition i.e. , then it is straightforward to uniformly bound SOM (). We can also perform the sum over all sites to get a uniform upper bound of . Lastly we just have to count the number of terms of , but this grows at most linearly in , which we upper bound as for SOM (). Putting everything together, we obtain:

 σ([ω])<πβω(2√2C(2)λJeν/2ω)2p(2p)! (11)

Using and finding the optimal for which Eq. (11) is minimized yields SOM (), and

 σ([ω])<πβωe−Bω. (12)

for SOM () that depends on and some system size independent numerical factors.

Global drive.— We next consider the case where the drive is global, but for now, still without the on-site field (), so that . We have to therefore analyze the global version of Eq. (7), in which we replace one by and the other by , and where we further sum over all sites . The only difference between this expression and that of the previous case is the presence of additional off-diagonal terms. However, such terms can be accounted for rather simply: for a fixed , terms arising in the -nested commutator that might potentially contribute must have support that overlaps with site . This is because if the term in the -nested commutator emanating from site is disjoint from , then the trace over the system can be broken up into the two subregions and its exclusion , which evaluates to since . With this fact, we can proceed with the analysis almost verbatim, defining in analogy to , a new function where once again the placement of the indices depends on the partition in question.

The only difference in the subsequent analysis is Eq. (10), which becomes

 σ([ω]) ≤Ap∑x,y∑s∈S′(2p)∑i1,⋯,il∑a1∈{x,i1}⋯∑al(s)∈{x,i1,⋯,il(s)−1} ∑q∈Q(s)fxy(q)Z0δy,i1+δy,i2+⋯δy,il(s)rαn1(q)x,i1rαn2(q)a1,i2⋯rαnl(s)(q)al(s)−1,il(s), (13)

where the (at most number of) Kronecker-deltas enforce connectivity of the terms. Then, summing over and proceeding as before, we end up with the global version of Eq. (11) in which the bound picks up a factor of , the total number of sites in the system, reflecting the extensivity of the heating rate. At the optimal which minimizes the expression, we therefore obtain the bound Eq. (5) SOM ().

Lastly, the most general case when there is a static on-site field can be analyzed in almost identical fashion, as the on-site terms cannot ‘grow’ the terms in a -nested commutator. We hence once again arrive at the result Eq. (5) (see SOM () for details).

Prethermal effective Hamiltonian and numerics.— Let us now discuss the implications of our results. We have shown that heating (over one Floquet cycle) due to direct transitions between eigenstates of the undriven Hamiltonian separated by energy apart are exponentially suppressed in frequency. While this is a result obtained within linear response theory and at high temperatures, it strongly suggests that there should be a rotating frame of reference (effected by some time-dependent unitary ) in which stroboscopic dynamics is equivalently described by a new Hamiltonian

 H′(t)≡Q(t)†(H(t)−i∂t)Q(t)=Heff+geffVeff(t), (14)

such that is a static effective Hamiltonian, and is a remaining driving piece. Since in this frame is the term that drives direct transitions between states separated by energy apart, its effective coupling should be exponentially suppressed, , for some effective interaction strength which is related to the local energy scale , as per our result Eq. (5). Parameterizing the unitary as , this suggests that the generator , and hence also , are organized in a power series in : , where the local norm of . is then a ‘dressed’ version of the undriven Hamiltonian . Hence, dynamics for times before the effects of ‘kicks in’ should be well captured by the static effective Hamiltonian , i.e. a ‘prethermal’ regime , while for times , heating due to should result in an error in measurements of any local observable evolved by the exact Floquet dynamics and the effective Hamiltonian that grows linearly in time: .

In order to support this highly plausible and appealing picture, we turn to numerics. We consider a family of 1d long-range spin Hamiltonians

 H(t) =∑ijsijrαij(Jzzσziσzj+Jxxσxiσxj)+∑ihxσxi +g(1−2θ(t−T/2))∑i(σzi+σyi), (15)

where are random in with equal probability (see num () for more details) and with . Because the drive in this case is a step function, it is natural to utilize the Baker-Campbell-Hauserdoff (BCH) formula and construct a family of effective Hamiltonians labeled by , defined as the -th order truncation of the BCH expansion BCH (). We note that the zeroth order effective Hamiltonian is nothing but the time-averaged Hamiltonian which is the first line of Eq. (15). Initializing our quantum state to be the product state in the -basis with energy density closest to , we evolve it in time (via Krylov subspace methods) by both the exact Floquet dynamics and the effective Hamiltonians, and measure its energy density as with respect to .

Fig. 1 shows our results. From Fig. 1(a,b), we see that the exact Floquet dynamics shows a an initial lack of heating, before an eventual heating up (to infinite temperature), while the dynamics generated by the effective Hamiltonian never shows any heating of initial state. More importantly, characterizing this difference in errors of the measurement of the energy density between states evolved by exact Floquet dynamics and the higher-order effective Hamiltonians, we see from Fig. 1(c,d) that at sufficiently high frequencies (larger than but smaller than the many-body bandwidth), this error is small and constant for time , indicating the presence of a ‘prethermal’ plateau, while for there is an apparent linear increase in error . Extracting the dependence of the slope of this linear increase in Fig. 1(e,f) gives an excellent agreement with the statement that the heating rate is exponentially suppresed in frequency, , and allows for an extraction of the effective interaction strength . These numerical results are therefore consistent with, and support the presence of effective Hamiltonians in high-frequency driven disordered long-range interacting systems.

Summary and discussion.— We have shown that the heating rate of periodically driven, long-range systems is exponentially suppressed at high frequencies, and furthermore provided numerial evidence to indicate the presence of a prethermal, effective, static Hamiltonian governing well stroboscopic dynamics for exponentially long times. The bounds we have derived are relevant to understanding and constraining dynamics in experimentally accessible setups of long-range interacting degrees of freedom such as dipolar systems realized in ensembles of NV-centers in diamond or trapped ions, as they show that such systems should show a broad prethermalization window under a high-frequency drive. This thus opens up the possibility of realizing new prethermal phases and dynamical regimes, previously discussed only for short-ranged systems. On a mathematical level, we used the sign-changing nature of the interactions in long-range systems; in the future, it would be interesting to relax the assumption of high temperatures and prove non-perturbatively the existence of a prethermal effective Hamiltonian similar to Refs. Abanin et al. (2017b); Mori et al. (2016); Abanin et al. (2017a); Kuwahara et al. (2016). It would also be appealing to apply our techniques to obtain improved bounds on other dynamical properties in long-range systems, such as entanglement spreading.

###### Acknowledgements.
Acknowledgments.—We thank Curt von Keyserlingk, Vedika Khemani, Misha Lukin, Tomotaka Kuwahara and Rahul Nandkishore for useful discussions. We thank the Kavli Institute for Theoretical Physics, where this work was initiated, for hospitality during the program Synthetic Quantum Matter. W.W.H. thanks Joonhee Choi for help with the simulations. This research was supported by the Swiss National Science Foundation, and in part by the Russian Science Foundation under the grant No. 14-42-00044 (IP).

## Appendix A Supplemental material: Bounds on energy absorption in quantum systems with long-range interactions

In this supplemental material we provide additional details on the bounds derived in the main text, such as providing examples to illustrate the counting of terms using the restricted integer partitions and the partitions , and also more technical estimates for certain constants encountered.

In the main text, we analyzed three cases: (i) Local drive with no on-site field, (ii) Global drive with no on-site field, and (iii) Global drive with on-site field. In all three cases, the starting point for the bound on the high frequency part response function is

 σ([ω])≤πβωZ0ω2p∣∣Tr(V⟨[[[V,H],⋯],H](2p)⟩)∣∣. (16)

To get to this expression, we used the Lehmann representation of and performed the following manipulations:

 σ([ω]) :=⟨∫∞ωdω′∑nmπβω′Z0⟨n|V|m⟩⟨m|V|n⟩δ(En−Em−ω′)⟩ =⟨∫∞ωdω′∑nmπβω′Z0⟨n|[[[V,H],⋯],H](p)|m⟩⟨m|[[[V,H],⋯],H](p)|n⟩ω′2pδ(En−Em−ω′)⟩ ≤πβωZ0ω2p⟨∫∞ωdω′∑nm⟨n|[[[V,H],⋯],H](p)|m⟩⟨m|[[[V,H],⋯],H](p)|n⟩δ(En−Em−ω′)⟩ ≤πβωZ0ω2p⟨∑n∑m:Em∈(−∞,−ω+En)⟨n|[[[V,H],⋯],H](p)|m⟩⟨m|[[[V,H],⋯],H](p)|n⟩⟩ ≤πβωZ0ω2p⟨∑nm⟨n|[[[V,H],⋯],H](p)|m⟩⟨m|[[[V,H],⋯],H](p)|n⟩⟩ =πβωZ0ω2p⟨Tr([[[V,H],⋯],H](p)[[[V,H],⋯],H](p))⟩ Missing \left or extra \right =πβωZ0ω2p∣∣Tr(V⟨[[[V,H],⋯],H](2p)⟩)∣∣, (17)

where are eigenstates of (one disorder realization of) . In the second line, we introduced commutators of with together with the energy differency , which is an equality under the delta-function; in the third line we uniformly bounded the denominator of all terms by (since every term is positive); in the fourth line we performed the integral and in the fifth line we let the sum extend over all eigenstates. Then, in the second last line, we made use of the fact that repeatedly to transfer all the commutators of one term with -nested commutators to the other term, to end up with a -nested commutator, while in the last line we made use of the fact the disorder averaging commutes with the trace operation, which is possible only in the high temperature limit.

## Appendix B Local drive with no on-site field: additional details

In the case of local driving where which is assumed to act only on site , we have

 σ([ω]) ≤πβωZ0ω2p∑→μ∣∣Tr(V0[[[V0,Oμ1],Oμ2],⋯,Oμ2p])∣∣ ×∣∣⟨Jμ1Jμ2⋯Jμ2p⟩∣∣rαμ1rαμ2⋯rαμ2p (18)

where and each sum is over all links in the system. For a finite system of sites, there are links and each sum runs over some enumeration of the links.

### b.1 Counting using integer partitions S(2p),S′(2p) and partitions Q(s)

As explained in the main text, because each has be to be at least paired (owing to ), a natural way to organize the counting is through the number of times distinct () appear. This makes us consider the set of all integer partitions of the integer with each integer part , which we denote by , where is the length of the integer partition. For example, for , there are two integer partitions :

 2p=4={(4)(2,2) (19)

with lengths respectively. So Eq. (18) for the case of a four-nested commutator becomes organized as

 ∑→μ⋯ =∑μ1⋯+∑μ1∑μ2≠μ1⋯ =∑s∈S(2p=4)′∑μ1,⋯,μl(s)⋯, (20)

where the prime on the sum means distinct links.

However, as mentioned, we also need to consider the ways that links are distributed in the commutators. This means that for each integer partition , we have to consider its partitions of , which is the set of all partitionings of the list into parts, ordered by the smallest element appearing in each part, with the th part having elements where is an integer part of . For example, for the integer partition of , there are three partitions of :

 s=(2,2)→Q(s)=⎧⎪ ⎪ ⎪ ⎪⎨⎪ ⎪ ⎪ ⎪⎩(1,2),(3,4)(1,3),(2,4)(1,4)part,(2,3). (21)

As can be seen, each part of each partition is ordered (and the partition itself is ordered by the smallest element appearing within each part), and the number of elements in each part is this case is .

What gives us information about is the positionings of how the links are distributed: the first part of a partition tells us the positions in the -nested commutators where appears; the second part of that partition tells us the positions where appears, and so on. For the above example of , Eq. (18) reads

 ∑μ1|Tr(V0[[[[V0,Oμ1],Oμ1],Oμ1],Oμ1]||⟨J4μ1⟩|r4αμ1+∑μ1∑μ2≠μ1 (|Tr(V0[[[[V0,Oμ1],Oμ1],Oμ2],Oμ2]||⟨J2μ1⟩⟨J2μ2⟩|r2αμ1r2αμ2 |Tr(V0[[[[V0,Oμ1],Oμ2],Oμ1],Oμ2]||⟨J2μ1⟩⟨J2μ2⟩|r2αμ1r2αμ2 |Tr(V0[[[[V0,Oμ1],Oμ2],Oμ2],Oμ1]||⟨J2μ1⟩⟨J2μ2⟩|r2αμ1r2αμ2). (22)

More generally, this leads to the exact representation of Eq. (18)

 σ([ω])≤ πβωZ0ω2p∑s∈S(2p)′∑μ1,⋯,μl(s)∑q∈Q(s)f(q)× |⟨Jn1(q)μ1⟩⟨Jn2(q)μ2⟩⋯⟨Jnl(s)(q)μl(s)⟩|rαn1(q)μ1rαn2(q)μ2⋯rαnl(s)(q)μl(s), (23)

as stated in the main text.

Now, as stated in the main text, if we restrict the sum over integer partitions in to be over those that have integer parts being only or , i.e. , while simultaneously lifting the restriction of distinct links in the sum over , we get an upper bound

 σ([ω])≤ ApZ0∑s∈S′(2p)∑μ1⋯μl∑q∈Q(s)f(q)rαn1(q)μ1⋯rαnl(q)μl, (24)

where . This is because we can cover every integer partition of and its partitions (assuming links are distinct) non-uniquely by some integer partition of together with its partitions (assuming some links are distinct). Going back to the example of , if we allow to run over values of in Eq. (22), then we can write Eq. (24) as follows. The only integer partition in is , so the bound is

 ∑μ1∑μ2(|Tr(V0[[[[V0,Oμ1],Oμ1],Oμ2],Oμ2]||⟨J2μ1⟩⟨J2μ2⟩|r2αμ1r2αμ2|Tr(V0[[[[V0,Oμ1],Oμ2],Oμ1],Oμ2]||⟨J2μ1⟩⟨J2μ2⟩|r2αμ1r2αμ2 |Tr(V0[[[[V0,Oμ1],Oμ2],Oμ2],Oμ1]||⟨J2μ1⟩⟨J2μ2⟩|r2αμ1r2αμ2) =3∑μ1|Tr(V0[[[[V0,Oμ1],Oμ1],Oμ1],Oμ1]||⟨J4μ1⟩|r4αμ1+∑μ1∑μ2≠μ1(|Tr(V0[[[[V0,Oμ1],Oμ1],Oμ2],Oμ2]||⟨J2μ1⟩⟨J2μ2⟩|r2αμ1r2αμ2 |Tr(V0[[[[V0,Oμ1],Oμ2],Oμ1],Oμ2]||⟨J2μ1⟩⟨J2μ2⟩|r2αμ1r2αμ2+|Tr(V0[[[[V0,Oμ1],Oμ2],Oμ2],Oμ1]||⟨J2μ1⟩⟨J2μ2⟩|r2αμ1r2αμ2) (25)

which indeed over-estimates Eq. (22).