{'fraction_non_alphanumeric': 0.029284108817526546, 'fraction_numerical': 0.009520241816725656, 'mean_word_length': 4.474506753412064, 'pattern_counts': {'":': 0, '<': 1, '<?xml version=': 0, '>': 10, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 7, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Enhanced Perturbative Continuous Unitary Transformations

Unitary transformations are an essential tool for the theoretical understanding of many systems by mapping them to simpler effective models. A systematically controlled variant to perform such a mapping is a perturbative continuous unitary transformation (pCUT) among others. So far, this approach required an equidistant unperturbed spectrum. Here, we pursue two goals: First, we extend its applicability to non-equidistant spectra with the particular focus on an efficient derivation of the differential flow equations, which define the enhanced perturbative continuous unitary transformation (epCUT). Second, we show that the numerical integration of the flow equations yields a robust scheme to extract data from the epCUT. The method is illustrated by the perturbation of the harmonic oscillator with a quartic term and of the two-leg spin ladders in the strong-rung-coupling limit for uniform and alternating rung couplings. The latter case provides an example of perturbation around a non-equidistant spectrum.


I. INTRODUCTION
Quantum many-body systems with correlations are notoriously difficult to describe theoretically. Many analytical and numerical tools have been developed to tackle such problems. A tool which is employed ubiquitously are unitary transformations. Famous applications are the fermionic Bogoliubov transformations in the mean-field theory of superconductivity by Bardeen, Cooper and Schrieffer 1 or the bosonic Bogoliubov transformations arising in linear spin wave theory of quantum antiferromagnets 2 . These are exact transformations which use the algebraic properties of fermions and bosons, respectively. They yield diagonal Hamiltonians if they are applied to bilinear initial Hamiltonians.
Another class of unitary transformations are those which are not exact but approximate because they rely on an expansion in a small parameter. A well-known example is the antiferromagnetic Heisenberg exchange coupling J as it is derived from a half-filled Hubbard model with hopping t and local repulsion U implying J = 4t 2 /U , see for instance Ref. 3. Obviously, higher contributions O t 3 /U 2 are neglected. But they can also be computed systematically [4][5][6][7][8][9][10] .
Moreover, the Hubbard model is not diagonalized by the transformation but mapped to an effective spin model. This mapping implies a simplification because the relevant part of the Hilbert space (here: spin degrees of freedom) has been separated from the remainder (charge degrees of freedom). The remainder does not need to be considered anymore. It is said that it has been eliminated or integrated out.
Another famous example in the same line is the Fröhlich transformation 11 which eliminates phononic degrees of freedom from an electron-phonon system in leading order of the coupling to derive an electron-electron interaction from an electron-phonon coupling. Starting from the Hamiltonian this transformation generates an attractive interaction in the BCS-channel with the matrix element where q := k − k, ∆ε = ε k −ε k . This explains the formation of Cooper pairs and conventional superconductivity. It is interesting to note that in standard treatments the interaction is usually approximated by a constant, leaving out any discussion of the resonance singularity in (3). It is, however, possible to achieve the elimination of the phonon degrees of freedom by a different, continuous unitary transformation (CUT). This approach relies on a continously parametrized antihermitean generator η( ) = −η † ( ) of the differential unitary transformation The eye-catching fact in V W k, k is that it does not have a resonant energy denominator. Hence this result is much smoother. In particular it implies an attractive interaction for all parameters.
The standard BCS interaction is a constant up to some phononic cutoff energy ω Debye . This result can be derived rigorously by a modification of the generator. In an eigen basis of H D the matrix elements of η MKU are chosen [17][18][19] to be η MKU ij := sgn(E i − E j )H ij . Then we find where Θ(x) is the Heaviside step function. Again, there is only attractive interaction. In addition, the interaction is only active in a restricted energy interval and zero outside.
It is very remarkable that all three approaches 3, 5, and 6 are different in their outcome although they do the same: Eliminating the linear electron-phonon coupling. We stress that this is not a spurious result, but relies on the fact that the unitary transformations are indeed different even in leading order. They express virtual processes in a different way. But the energy conserving processes at ∆ε = 0 are the same in all three results. This has to be so because such scattering processes can in principle be measured which implies that they have to be independent from the chosen basis.
We include the above results to demonstrate the power of the freedom to choose an appropriate basis. The derivation of the BCS electron-electron interaction illustrates the versatility of continuous unitary transformations. Moreover, it exemplifies that even a seemingly simple task, the derivation of an effective interaction in leading order, provides unexpected surprises. We advocate to do it in a less singular way with effective interactions which are restricted in energy range.
The main goals of the present article are twofold. First, we show how CUTs can be used to perturbatively derive effective Hamilton operators in real space. This goal has been realized for an unperturbed Hamiltonian with equidistant spectrum by perturbative CUT (pCUT) [18][19][20] . The gist of the pCUT is recalled below. In the present work we enhance the applicability of such an approach to unperturbed non-equidistant spectra by formulating the CUT directly in second quantization. The resulting transformation will be called enhanced perturbative CUT (epCUT) for distinction. The approach is exemplified for a uniform and for an alternating spin ladder. The latter has a non-equidistant spectrum if only the rung couplings are considered.
The second main goal is to establish a robust extrapolation of the perturbative results of the epCUT. We will show that a direct evaluation of the perturbatively established flow equation provides a very robust and reliable way to extrapolate the perturbative results. This approach will be called directly evaluated enhanced perturbative CUT (deepCUT).
The article is set up as follows. In the remainder of this section we briefly introduce the perturbative CUT and the self-similar CUT (sCUT) as predecessors of the epCUT and the deepCUT. In Sec. II, we introduce our paradigm model, spin ladders, for which we illustrate the general approaches. In Sec. III, we derive the epCUT and develop the deepCUT from it. Many technical aspects are discussed; a focus are simplifying rules which allow us to compute high orders efficiently. In Sec. IV, results of the epCUT and the deepCUT are presented for the uniform antiferromagnetic spin ladder with S = 1/2. Results for the alternating spin ladder, which does not have an equidistant unperturbed spectrum, are shown in Sec. V. The article terminates by the conclusions in Sec. VI.

A. Perturbative Continuous Unitary Transformation
We draw the readers' attention to the fact that we are dealing from now on with CUTs with a unique reference state. This means that the ground state is mapped by the CUT to the vacuum of excitations. In the Introduction, the mapping to effective models such as the Heisenberg exchange model or the BCS-model still left a many-body problem to be solved.
High order series expansions have long been used to compute reliable ground state energies and dispersions in strongly-correlated systems 21 . These quantities can be computed straightforwardly because the states are uniquely determined by their quantum number, for instance the momentum, even before the perturbation is switched on. This means it is sufficient to perform the perturbation for a one-dimensional subspace of the Hilbert space. States of two and more particles are more subtle because their subspaces are extensively large. For instance, binding energies cannot be computed as series unless the binding occurs already in linear order. Generally, unitary or orthogonal transformations must be introduced to define the perturbative approach on large subspaces 18,19,[21][22][23] .
The perturbative CUT was the first approach to achieve this goal 18,19 . Its starting point is a Hamiltonian which can be written in the form where H 0 is the unperturbed Hamiltonian with an equidistant spectrum. For simplicity we set its energy spacing to unity. Each energy quantum can be seen as an elementary excitation, a quasi-particle, so that H = 0 counts the number of quasi-particle up to an irrelevant offset. The expansion parameter is x and the terms in the perturbation are split according to their effect on the quasi-particle number H 0 : The terms in T m increase the number of energy quanta by m. Obviously, T −m = T † m holds. Generically, there is an upper bound N ≥ |m| to the change of energy quanta. Choosing η MKU = N m=−N sgn(m)T m a flow is defined via (4) which generates a hierarchy of differential equations in powers of x. In each finite order, the differential equations are closed and can be solved by computer aided analytics. Eventually, one obtains the general expansion for an effective Hamiltonian which conserves the number of elementary excitations The components of the vector m take the values −N, −N + 1, . . . , N − 1, N ; the vector has the dimension dim( m) = k. The coefficients C( m) are simple real fractions and T ( m) stands for the product The conservation of the number of quasi-particles is implied by the cross sum (8) is very general; to put it to practical use its irreducible effect on zero, one, two, and more quasiparticles is computed. In this way, the effective Hamiltonian is obtained in second quantized form 20 . Remarkable achievements of this approach are a quantitative understanding of inelastic scattering in spin ladders [24][25][26][27][28][29] , of spectral densities in spin chains 30 , of excitations in the Kitaev model 31,32 , of excitations in the toric code [33][34][35] , and of the ionic Hubbard model 36 to name a few extended systems where the ground state is described as a vacuum of excitations.
Conceptually, the most significant achievement of pCUT is that whole subspaces are treated perturbatively. The generality of the pCUT result (8) is surely one of its advantages. The fact that it can only deal with equidistant spectra is a certain caveat. Another caveat is that the approach does not allow for modifications of the generator.

B. Self-Similar Continuous Unitary Transformations
One way to circumvent the above mentioned restrictions concerning the unperturbed spectrum and the choice of the generator is to pass from a perturbative evaluation to a self-similar one. The approach follows a straightforward strategy. One chooses a set of operators Figure 1. (Color online) Schematic representation of the symmetric (alternating) S = 1 /2 Heisenberg ladder in the thermodynamic limit.
(DES) in the coefficients of the basis operators. A more detailed description follows in Sec. III below. We stress that in the latter step a certain truncation is required 37 . Unless the set of operators is closed under commutation, the commutator [η, H] comprises terms which cannot be expanded exactly in the operator basis. Thus this step generically requires an approximation. The Hamiltonian is kept in a self-similar form defined by the selected operator basis. Depending on the system the truncating approximation can be controlled by a small parameter 15,38,39 or by the spatial locality of the selected set of operators 7,40 . We stress that in the sCUT approach the choice of the operator basis and of the generator uniquely defines the DES of the flow equation. Clearly, the advantage of the sCUT over the pCUT is its larger versatility. Yet it is less general in the sense that the flow equation has to be solved for each model and each operator basis anew.
To derive a systematic perturbative expansion by sCUT is not an obvious step. It is one of our two main goals to show how this can be done and how it can be done efficiently. Thus the derivations and considerations in Sec. III are based on the sCUT approach and combine it with a perturbative expansion in order to reach the enhanced perturbative CUT.

II. MODEL
To illustrate the performance of the (de)epCUT we consider the S= 1 /2 antiferromagnetic two-leg Heisenberg ladder (uniform spin ladder) and an extension with an alternating rung coupling (alternating spin ladder) as testing ground, see Fig. 1

. The Hamiltonian reads
where r ∈ Z. The rung number is denoted by r and the legs by L and R. We define the ratio between the leg coupling J and the even rung coupling J e ⊥ as x := J /J e ⊥ and the ratio between the odd rung coupling J o ⊥ and the even rung coupling J e ⊥ as y := J o ⊥/J e ⊥ . In the limit of J e ⊥ = J o ⊥ , i.e., y = 1, the Hamiltonian describes the uniform spin ladder. This model has been subject of intensive studies, see Refs. 27, 41, and 42 and references therein. Thus it constitutes a suitable reference model to test the epCUT. It has been investigated by several different methods, such as density matrix renormalization 43,44 , exact diagonalization 45 , continuum field theory 46,47 , quantum Monte Carlo 48 , high order series expansions 49,50 , including methods based on CUTs, such as sCUT 40,51,52 and pCUT 20,26,27,53 . If the results of the epCUT agree with this data, the efficiency of the epCUT for the expansion around an unperturbed equidistant spectrum is verified.
To illustrate that the epCUT represents an advancement compared to pCUT we will show results for the alternating spin ladder as well. This system does not have an equidistant spectrum because the rung couplings are not equal J e ⊥ = J o ⊥ . Hence it cannot be dealt with by pCUT. Without loss of generality we consider J o ⊥ > J e ⊥ implying y > 1.
For the alternating spin ladder we expect a lowering of the ground state energy upon rising y because the expectation value of S L r · S R r is negative. The unit cell includes two rungs which implies two triplon branches in the Brillouin zone (BZ). For y = 1 the branches meet at the BZ boundary (k = ± π /2). For y > 1 a band gap of the order of |y − 1| opens at k = ± π /2 separating the two band.
To define a starting point for the CUT the bond operator representation 54,55 is used. A possible eigen basis of the local operators S L r , S R r is given by the singlet state and the three triplet states For x = 0, the ground state of the system is given by This vacuum of triplets serves as our reference state. The local operators t † x,r , t † y,r and t † z,r (t x,r , t y,r and t z,r ) create (annihilate) an excitation on rung r. They satisfy the hardcore-boson commutation relation The elementary magnetic excitations (S=1), known as triplons 27,56 can be continuously linked to the local triplets.
Represented in second quantization in terms of the triplon creation and annihilation operators the Hamiltonian reads where where a, r ∈ Z. This form of the Hamiltonian enters all the calculation described in the following.

A. Flow equation in second quantization
Similar to the implementation of previous CUT methods 7,12,20,40 , we formulate the flow equation (4) for the coefficients of the monomials {A i } of operators in second quantization. The Hamiltonian is parametrized by with the -dependent coefficients h i ( ). The generator reads withη being a superoperator denoting the application of a particular generator scheme such as those discussed in Ref. 40. Expanded in the operator basis {A i } the flow equation (4) reads Comparing the coefficients of different monomials, the flow equation (4) becomes equivalent to a set of ordinary differential equations for the coefficients h i ( ) The commutator relations between the basis operators are encoded in the coefficients D ijk of the bilinear differential equation system (DES). These coefficients D ijk are in general complex numbers. For the spin ladders under study they are given by integers or fractional numbers. We call a single D ijk a 'contribution' of the DES. The contributions are obtained from by comparing the coefficients of the expansion of the commutator monomial by monomial.
In this way, the problem of solving the flow equation is transformed into the algebraic problem of calculating the coefficients of the DES (21b) and of the subsequent numerical solution of Eq. (21a).

B. Perturbative expansion of the flow equation
Here we consider the perturbative solution of the flow equation which yields the resulting effective Hamiltonian in the form of a perturbative series. Hence this solution generalizes the established pCUT approach 18,19 . To this end, we decompose the initial Hamiltonian into an unperturbed part H 0 and a perturbation V . In contrast to pCUT 19 , we do not require the unperturbed part to have an equidistant spectrum. We do not require other restrictions either although we will see that the method works best for a (block-)diagonal H 0 . We aim at the perturbation series up to and including order n in x. Thus we expand the flowing Hamiltonian into terms of order x m up to m ≤ n. Expanding the H (m) in the operator basis {A i } we perform the expansion in powers of x by expanding the coefficient Figure 2. Sketch of the epCUT algorithm to calculate the DES for the iterative calculation of ∂ H (4) . Due to the commutators [η (1)  At l = 0, the initial values f For the prefactors of x m this implies We stress that the contributions D ijk do not depend on the order m of the coefficients, but only on the algebraic relations between the corresponding monomials. Hence they need to be calculated only once. Moreover, Eq. (26) defines a hierarchy between the coefficients because f (m) i ( ) is influenced only by coefficients of the same order m or lower, but not by coefficients of higher orders.

C. Algorithm
A key task in the implementation of epCUT is the design of an efficient algorithm to identify the monomials and to calculate exactly the commutators which are relevant for the transformed Hamiltonian in the order of interest n. Henceforth, we call the order we are aiming at the 'targeted' order.
Based on Eq. (26) we can calculate each order m based on the results of lower orders. Order zero is trivially given by the representation (22) if η[H 0 ] = 0 which means that H 0 is block-diagonal. The calculation of the commutators η (1) , H (m−1) , . . . , η (m−1) , H (1) can be carried out independently, see Fig. 2. According to Eq. (21b), the commutator [η[A j ], A k ] can be written as linear combination of monomials A i of which the prefactors define the contributions D ijk of the DES. For those monomials not yet present in the Hamiltonian, a new monomial has to be included in the operator basis with a unique index. We call the order in which a monomial occurs for the first time its minimum order O min (A i ).
We stress that in the evaluation of η The calculation of the commutators for η (m) , H (0) is special because its result may include new monomials of the same minimum order m which were not considered so far. Since these monomials also enter the commutator via η (m) , the block η (m) , H (0) has to be iterated until no new monomials occur: Then self-consistency is reached. This should be done once the inner blocks η (p>0) , H (q>0) are finished.
If the unperturbed Hamiltonian H 0 is local, the commutation of monomials from η (m) and H 0 lead to monomials acting on the same local cluster or smaller subclusters. Furthermore, if the local Hilbert space of the cluster is finite, the number of new monomials which can be generated by iterative commutations with H 0 is bounded by the finite number of linearly independent matrices on this finite-dimensional Hilbert space. Then the iterative loop is guaranteed to terminate after a finite number of cycles.
In the symmetric ladder model (see Eq. (9d)) the local Hilbert spaces are finite so that a finite number of cycles is certainly enough. Even better, the commutation of the monomials in terms of triplon creation and annihilation operators with H 0 = J e ⊥ H e ⊥ + J o ⊥ H o ⊥ does not generate any new monomials so that no iterations are needed in the calculation of η (p>0) , H (q>0) .
If the unperturbed Hamiltonian H 0 has also non-(block-)diagonal terms, the generator includes terms of order zero. Therefore, the blocks η (0) , H (m) have to be evaluated self-consistently as well. Since any term of the Hamiltonian may also appear in the generator, the blocks η (m) , H (0) and η (0) , H (m) have to be calculated simultaneously within a joint self-consistency loop.
For the sake of completeness, we note that in the special case H 0 := H, i.e., considering the total Hamiltonian as the unperturbed one, the whole algorithm constructing the DES reduces to the calculation of the block η (0) , H (0) . This has to be done self-consistently with respect to both the generator and the Hamiltonian. This approach is the one employed in the self-similar CUT (sCUT) previously 39,40 . Since for H 0 = H "the unperturbed" part H 0 includes non-local and non-(block-)diagonal terms and perhaps refers even to an infinite local Hilbert space for any but the simplest models, the iteration of commutators will not terminate. Thus additional truncation criteria are needed whose validity needs to be justified.

D. Perturbative evaluation
To evaluate the perturbation series for the ground state energy or the dispersion relation of a concrete system, the first step is to write the Hamiltonian in second quantization and to identify the relevant monomials. We illustrate the approach by the second order calculation for the symmetric spin ladder in Eq. (9d) using the quasiparticle conserving generator 19,40 .
The operator basis {A i } is given in Tab. I with A 0 and A 1 for the terms in H 0 (O min = 0) and A 2 to A 5 for the terms in V (O min = 1). We combined certain monomials whose prefactors must be the same due to symmetry and/or hermiticity into one element of the operator basis A i (cf. Sec. III E 2). The advantage is that less operators need to be tracked. The algorithm is not affected by this step except that the comparison of coefficients is a bit more complex.
Following the algorithm described above, the commutators of the block η (1) , H (0) are calculated to complete the first order. The contributions to the DES obtained by comparison of coefficients are given in Tab. II. Then, the contributions in second order are evaluated in the blocks η (1) , H (1) and η (2) , H (0) leading to the new basis operators A 6−17 with O min = 2.
Next, the perturbative flow equation (26) has to be solved. We do this numerically using a standard fourth order Runge-Kutta method 57 . The initial values for the coefficients in different orders of x are read off the initial Hamiltonian. They are zero for all basis operators and all orders which are not present in the initial Hamiltonian. We use a basis of only normal-ordered operators except for A 0 = r 1 so that the series expansion of the ground state energy per rung E 0 is obtained in the limit of → ∞ from the prefactor of A 0 Note that this result requires only three equations in the DES. Likewise, the dispersion relation is determined from the renormalized coefficients of the hopping terms Table I. Basis operators Ai (simple combinations of monomials obeying symmetry and/or hermiticity) occuring in a second order epCUT calculation for the uniform spin ladder using the particle conserving generator scheme ηpc. The third column contains the initial coefficients hi( = 0), the fourth the final renormalized coefficients hi( = ∞). The minimum order O min is given in which the corresponding operator occurs for the first time; O 0qp max is the highest relevant order of the coefficient for computing the ground state energy and O 1qp max is the highest relevant order for computing the dispersion. The terms marked in light gray are irrelevant for the computation of the ground state energy in second order; the terms in dark gray are irrelevant if the dispersion is computed.  Table II. Non-vanishing contributions D ijk to the differential equation system (DES) of the uniform spin ladder using the particle conserving generator scheme. Operators and contributions marked in light gray are irrelevant for the computation of the ground state energy in second order; those in dark gray are irrelevant for the dispersion.
which require five equations, only two more than the ground state energy.

E. Optimizations
The epCUT method presented so far can be applied to a wide range of models in order to calculate a perturbative expansion of decoupled quasi-particle spaces. With increasing order the number of representatives in the effective Hamiltonian, the runtime and the memory consumption rise exponentially, see Fig. 3. One is interested  in increasing the order of the calculation as high as possible because this generically enhances the accuracy of the calculation: More and more orders kept imply that more and more physical processes with an increasing range are taken into account.
To increase the order, more efficient generator schemes and the symmetries known from sCUT can be exploited. Focussing on selected quantities of interest, the perturbative foundation of epCUT allows us to optimize the algorithm even further. Generic performance data possible with full optimizations are given in Tab. III.
In practice, every optimization is carefully checked by comparing the results of the optimized faster program to the results from the slower program before optimization. In this way, one can be sure that no errors are introduced by incorrect assumptions.

Generator scheme
The (quasi-)particle-conserving generator scheme η pc used in our example decouples all subspaces of differ-ing numbers of excitations, i.e., quasi-particles, and sorts them in ascending order of their energy [17][18][19]40 . In most applications, however, only the ground state and the lowlying excitations are of interest. Consequently, the computational effort can be reduced by choosing a more efficient generator scheme which targets the quantities of interest only. In 2010, Fischer, Duffe and Uhrig 40 proposed a family of generator schemes based on modifications of η pc where only the first q quasi-particle spaces are decoupled from the Hilbert space. The corresponding generator reads In this notation, H i j comprises all monomials of the Hamiltonian creating i and annihilating j quasi-particles. For instance, the ground state generator incorporates monomials which consist purely of either creation or annihilation operators. Compared to the full quasi-particle conserving generator, the effort to compute the corresponding DES is reduced significantly. In analogy to η pc , the decoupled quasi-particle spaces are sorted according to energy. Thus the ground state energy is given by the vacuum energy of the effective Hamiltonian H(∞). If additionally the dispersion is calculated, the one-quasi-particle subspace has to be decoupled using η 1 . For decoupling higher quasi-particle spaces, analogous generator schemes can be used. But the increase in efficiency compared to the full quasi-particle conserving generator becomes less and less significant because the generator schemes η q do not conserve the block-banddiagonal structure of the Hamiltonian in contrast to η pc .

Symmetries
For models defined on infinite lattices, it is necessary to use the translation symmetry in order to be able to work directly in the thermodynamic limit. In addition, the presence of other symmetries leads to linear dependences of coefficients of monomials which are linked by the symmetry transformations of the Hamiltonian. As in sCUT 7,58 , this redundance can be significantly reduced by passing from simple monomials to symmetric linear combinations of them. Each of these polynomials is invariant under symmetry transformations of the Hamiltonian and requires only one prefactor where the single monomials would need many more. In our example (Tab. I) and in the following calculations the size of the operator basis is reduced by a factor of almost 24 exploiting self-adjointness, reflection and spin symmetry.

Reduction of the differential equation system (DES)
Targeting only certain quantities up to order n, such as the ground state energy or the one-particle dispersion, the DES can be reduced. Here we discuss how this can be done in practice.
Besides the minimum order O min , a maximum order O max can be a assigned to each monomial and its coefficient h i . The maximum order is the highest order of the series of h i which still has an influence on the targeted coefficients up to order n. For instance, complicated processes involving many quasi-particles do not influence the ground state energy directly, but only via other process. Then their O max is much lower than the targeted order n. Technically, this is due to the hierarchy of the DES (26) which implies where the equality holds if we consider only a single contribution D ijk = 0. The inequality takes into account that there may be many pairs (i, k) for a given j. Thus, If A j is targeted, for instance the ground state energy per rung A 0 , its O max is n by definition. For illustration, we consider the DES for the uniform spin ladder in second order, see Tab. II. If we only target the ground state energy h 0 up to order 2, the maximum order of monomial A 4 is given by where we deal with equalities because there is only one contribution for ∂ h 0 ( ) in the DES and the O max (A 0 ) is known. In this case the maximum order O max of A 4 is lower than the targeted order 2.
The O max of all coefficients can be calculated on the basis of the entire DES and of the minimum orders. Note that Eq. (32) defines O max implicitly, i.e., one has to find the correct self-consistent solution. This is done by starting from A monomial A i is targeted if we want to compute its coefficient h i in the given order n. Starting from the initial choice (34) Thus they do not matter for the relevant quantities up to order n and can be discarded completely. Moreover, all contributions to the DES which use these terms can be neglected. In addition, all contribution can discarded for which holds. These considerations allow us to reduce the DES significantly. In Tab. II the reduction of the DES for the uniform spin ladder in second order is marked for the ground state energy (light gray) and for the dispersion (dark gray), respectively.
We stress that one has to know the entire DES to apply the O max concept as described above.

Simplification rules
The reduction of the DES discards a large number of monomials and of the contributions D ijk , see for instance Fig. 3a, which is essential for an efficient evaluation. But it would be even more advantageous if one avoided the calculation of the omitted terms before they are tediously computed. The minimum orders O min are known at each step of the iterative setup of the DES so that they can be used on the fly. But due to their implicit definition the maximum orders O max are not known during the set-up of the DES.
Fortunately, estimates help. An upper bound for the maximum order is enough to accelerate the algorithm setting up the relevant part of the DES. Concomitantly, the memory consumption is reduced significantly 10 . Henceforth, we call such estimates 'simplification rules'. Their concrete form depends on the structure of the perturbed and the unperturbed Hamiltonian, for instance the blockdiagonality of the latter. We emphasize that the simplification rules constitute the part of the epCUT method which depends on the model.
In the following, we aim at a quantitative description up to order n of the effective Hamiltonian's blocks pertaining to at most q quasi-particles. For instance, q = 0 provides the correct perturbative expansion of the ground state energy and q = 1 allows us to calculate the dispersion relation up to order n.
A monomial creating c triplons and annihilating a triplons is targeted if both c ≤ q and a ≤ q hold. Its maximum order is the targeted order O max = n. If it is not targeted it can influence the targeted terms by affecting terms consisting of fewer creation and annihilation operators via the DES. For the Heisenberg ladder, the unperturbed Hamiltonian (Eq. (16)) is block-diagonal. Hence no commutation of generator terms with H 0 changes the number of created and annihilated triplons. The leading order of the generator is 1, i.e., η = O (x).
In the commutation of a monomial with a generator term some of the local creation and annihilation operators may cancel due to normal-ordering. In order to yield a term affecting the first subspaces with q quasi-particles local creation operators and local annihilation operators have to cancel. First, we consider commutations with lowest-order generator terms stemming from the initial Hamiltonian. In the spin ladder, these terms have order 1 and create or annihilate ∆QP = 2 quasi-particles on adjacent rungs. Because each commutation with η (1) increases the order of the affected coefficients by one, the maximum order is bounded by where the tilde on the left side means that one is dealing with an upper bound and y stands for the smallest integer that is still larger or equal to y. If in the calculation of ∂ H (m) the estimate O max of a monomial is lower than m, this contribution is irrelevant and can be omitted. This reduces the size of both the DES and of the Hamiltonian to be tracked. Moreover, discarding irrelevant monomials avoids the calculation of unnecessary commutators in the following iterations of the algorithm.
Clearly, the number of created and annihilated quasiparticles can be reduced by a number ∆QP larger than 2 by means of commutations with generator terms involving more quasi-particles which may have developed during the flow from the basic terms. But the generator terms involving more quasi-particles have a higher minimum order O min so that a single commutation with them affects coefficients only in a higher order m + O min . In fact, for the used generator schemes, the ratio between ∆QP and O min for new terms developed during the flow can not exceed the corresponding ratio for generator terms present in the initial Hamiltonian. Therefore, it is sufficient to consider only commutations with the initial terms in our simplification rules.
The above generic simplification rule can be easily adapted to other models as long as the unperturbed Hamiltonian H 0 is block-diagonal. Otherwise, H 0 will lead to generator terms of order zero which means that terms with high quasi-particle number can influence the coefficients of terms with low quasi-particle number in the same order. This is why it is desirable to set up the perturbation in such a way that H 0 is block-diagonal in the number of quasi-particles.
Applying the simplification rule reduces the number of representatives considerably, see Fig. 3, leading to a significant improvement of runtime and memory consumption. This basic simplification rule can be improved further by taking more model-specific information into account. A possibility to exploit the real-space structure of the monomials to lower the upper bound O max is described in Appendix A.
The computationally most costly part in the calculation of the DES is the evaluation of commutators. Because the simplification rules sketched above can only be applied after the commutation, we refer to them as aposteriori rules. For the sake of efficiency, it is highly desirable to extend them to a-priori rules estimating whether a commutator has to be evaluated at all prior to its computation. We describe such a-priori simplification rules in Appendices B and C.
Because these rules are necessarily less strict than there a-posteriori analogues, one should use the combination of both kinds in practice. The additional use of a-priori rules does not reduce the number of representatives or the memory consumption. But it boosts the speed of the calculation significantly because the vast majority of commutators can be discarded, see Fig. 3b, and the apriori rules help to avoid the laborious computation of these unnecessary commutators.

F. Directly Evaluated epCUT
In addition to the perturbative evaluation, the reduced DES computed by epCUT in a given order n can be evaluated non-pertubatively. After the reduction step described in Sec. III E 3, the DES consists exclusively of contributions which are relevant to the targeted quantities in the desired order n. This reduced DES in Eq. (21a) can be numerically integrated for any given value of x to obtain the coefficients of the Hamiltonian h i ( ) directly without passing by an expansion in x. In such a calculation all coefficients influence one another to infinite order. The numerical solution depends on the expansion parameter in an intricate manner and can no longer be understood as finite partial sum of an infinite series. In this sense, the perturbative reduced DES in order n is extrapolated by the direct evaluation in a non-perturbative way. To stress the difference to perturbation series computed by epCUT, we call this technique directly evaluated epCUT (deepCUT).
We emphasize that the reduction of the DES before the numeric integration is essential. It enhances the performance of the integration because the reduced DES is much smaller. But the crucial observation is that the reduction renders the integration much more robust. Numerical integrations of the full DES diverge for high orders and high values of x. We conclude that the reduced DES represents the relevant physical processes in a more consistent way. The integration of the full DES generates spurious higher order contributions which overestimate certain effects. In an exact solution the spurious higher order contributions would be compensated by other processes which are captured only in a higher order calculation.
Analogous observations are known from diagrammatic perturbation theory where the inclusion of subsets of diagrams in infinite order does not guarantee improved results. Improved results can only be expected from systematically controlled calculations. The inclusion of infinite orders is indicated if this achieves conserving selfconsistent approximations. For instance, the shift of poles in a propagator is not captured by any finite pertubation series in the propagator, but it follows easily from a perturbation of the self-energy 59 . Figure 4 compares the ground state energy per rung E 0 of the symmetric spin ladder as function of the relative leg coupling x obtained from the plain perturbative series in order 17, from various Padé extrapolations, and from the direct evaluation. The plain series shoots up at about x ≈ 0.7 while the Padé extrapolations start to scatter strongly beyond x ≈ 1. The direct evaluation lies between the two stiffest Padé extrapolations and remains stable up to even very large values of the expansion parameter x ≈ 3. Comparing various orders, see Fig. 5, the results of deepCUT converge rapidly and display only minor corrections for large values of x indicating a high reliability. In Sec. IV, a comparison of the deepCUT result with those of other methods will be presented.
This deepCUT bears similarities to the sCUT approach 7,37,39,40,52,58 . In sCUT, a set of basis operators is selected by a truncation scheme and for this set the full DES is computed. It comprises all commutation relations between the selected basis operators. In deepCUT, the order of the expansion parameter takes over the role of the truncation scheme. But we stress that deepCUT is not self-similar: In sCUT all commutators between the selected monomials are considered. In epCUT and thus in deepCUT only the commutators between specific subblocks based on the minimum orders O min are considered, see Sec. III C. Moreover, targeting certain subspaces with q quasi-particles and the concomitant reduction of the DES does not only discard irrelevant monomials. Also contributions linking relevant monomials are cancelled if their effect is of too high order. Therefore, the 'truncation' taking place in (de)epCUT, controlled by the expansion parameter, is a truncation of the DES rather than a truncation of operators as it is done in the sCUT approach.
One practical advantage of the deepCUT over the sCUT is that only one parameter, the maximum order of the expansion parameter, needs to be fixed in order to define the approximation. In the sCUT, generically many parameters define the truncation scheme 7,37,40,52,58 which leaves some ambiguity about how to systematically improve the approximation.
Another comparison of approaches is in order. Recently, Yang and Schmidt proposed a CUT approach based on graph theory (gCUT) 60 . To compute a certain quantity such as the ground state energy the irreducible contributions of subgraphs, i.e., of linked clusters, of the lattice are summed. The size of the largest subgraph considered determines the approximation. The larger it is the better the system is described because physical processes with a larger range are kept. Thus the fundamental idea of the approach is similar to the one of deepCUT: Truncation in the range of processes, but local processes are kept to infinite order.
The main difference is that the actual CUT is done on clusters. So gCUT has advantages and disadvantages. An advantage is that it is sufficient to deal with finite dimensional Hilbert spaces and the transformations can be performed on matrices. A disadvantage is that momentum conservation cannot be exploited on the level of the clusters because they are finite which restricts the choice of applicable generators 60 . The deepCUT is based on second quantization 22 and can take advantage of all symmetries of the problem under study. A detailed comparison of both approaches is left to future studies.

G. Transformation of Observables
In order to calculate spectral densities for instance, the coefficients of the corresponding observable must be known with respect to the same basis as the effective Hamiltonian. Thus the observables must be transformed as well. This can be realized by integration of the flow equation for observables introduced by Kehrein and Mielke 61,62 .
In analogy to the transformation of the Hamiltonian discussed in Sec. III A, we introduce an operator basis B i for the observable shifting the dependence on from where the second equation stands for the perturbative expansion of these coefficients. Hence the flow equation for observables (39) leads to a DES for their coefficients The contributions D obs ikj are obtained by calculating the commutators between the monomials of the generator and the monomials of the observable followed by a comparison of the coefficients The differential equations (41a) imply a hierarchical DES for the perturbative series (40b) for the coefficients ∂ f The algorithm for the calculation of the DES in Sec. III C can be easily adapted for the transformation of observables. Each order of the differential ∂ O (m) is calculated recursively, cf. Fig. 6. Since the generator η is defined solely by the Hamiltonian, it is not influenced by the outcome of the transformation of observables. For this reason the evaluation of η (m) , O (0) does not need to be carried out self-consistently. After the calculation of the commutators η (1) , O (m−1) . . . η (m−1) , O (1) , only the block η (0) , O (m) has to be treated self-consistently. But recall that η (0) only occurs if the unperturbed Hamiltonian H 0 is not (block-)diagonal. Because both differential equations (4) and (39) are coupled by the generator, their integrations have to be done simultaneously.
For the transformation of the Hamiltonian we extensively discussed that only certain contributions really matter. We introduced the concept of a maximum order in which the coefficient of a physical process needs to be known in order to influence the targeted quantities. This concept allowed us to reach significantly higher orders. Thus we want to extend the concept of a maximum order also to the transformation of observables. It turns out that this extension is more subtle than one may think.
Before, in the flow of the Hamiltonian, the maximum order of a generator coefficient O η,H max (A i ) is the maximum order of the same monomial O H max (A i ) in the Hamiltonian. Now we also target certain blocks of the observable and they are influenced by the monomials in the generator. This leads to maximum orders for both the observable term O O max (B i ) and the generator terms O η,O max (A i ). The latter does not need to coincide with the maximum order O η,H max (A i ) resulting from the consideration of the Hamilonian flow alone. Thus one has to find a unique and unambiguous way to fix O η max (A i ). We discuss three alternatives: (A) The maximum order of the generator terms is chosen in such a way that the targeted quantities in both the Hamiltonian and the observable(s) can be computed up to the targeted order 63 Then the iterative calculation of the O max must be realized within a single self-consistent loop. The perturbative evaluation yields a perturbative series for the coefficients of the observables under the transformation with the full generator up to order n. It may happen that in this way some generator terms are assigned a higher O η max ≥ O η,H max than in the transformation of the Hamiltonian alone so that the DES of the Hamiltonian comprises additional contributions. By construction, this does not affect the perturbative evaluation of the epCUT. But it will affect its direct evaluation (deepCUT) although it should be absolutely minor in a parameter regime of good convergence of the flow. Monomials which are discarded due to the reduction of the Hamiltonian will not be considered for the DES of the observables even though this may affect the targeted coefficients of the observable. Hence the transformation of the observables in perturbative evaluation is not realized with respect to the complete generator. We stress that this does not violate the unitarity of the transformation up to the order of calculation since the generator is still antihermitian and it is essentially the same as for the transformation of the Hamiltonian. No signficant deviations are expected in the regime of good convergence of the flow. Note also that any generator whose coefficients differ only by orders larger than O H max (A i ) leads to the same perturbative series for the relevant quantities in the Hamiltonian.
(C) A third alternative consists in taking over the O η,H max (A i ) for the reduction of the DES for the observables. Then only the values O O max (B i ) are computed selfconsistently.
For deepCUT, alternatives (B) and (C) ensure that the DES for the Hamiltonian is independent of the considered observables. Generally, we expect that the precision in the derivation of effective Hamiltonians is more important than the precision of matrix elements. Also in experiment, energies are generically known to much higher accuracy than matrix elements.
In order to keep the effective Hamiltonian in direct evaluation independent of the observables, we decide to use alternative (B) for deepCUT. For the perturbative evaluation, however, we favor alternative (A) because it makes the rigorous determination of the perturbation series of matrix elements possible.
The computational performance can again be increased decisively by applying simplification rules. They can be used directly for observables if both the Hamiltonian and the observables meet their requirements. This can often be achieved by appropriate definitions. For instance, the observable is needed for the calculation of the dynamic structure factor relevant for inelastic neutron scattering. But this observable includes non-block-diagonal terms in order zero.
To circumvent this problem, we consider the observable x · S L,z 0 instead. In this way, the non-block-diagonal monomials in the observable are shifted to order 1 so that they behave like the non-block-diagonal pertubation in the Hamiltonian. One loses an order of accuracy for a given fixed order n of the calculation. But all the simplification rules relying on block-diagonality in order zero can be used as before.

A. Ground State Energy
The ground state energy per rung is calculated up to and including order 17. The results of the direct (black solid line) and of the perturbative evaluation (dashed black line) are displayed in Fig. 7. The coefficients of the perturbative series, see Tab. IV, agree perfectly with the fractions from pCUT (up to order 14) 20 and with the decimal numbers (up to order 23) given by Zheng et al. 49 . This agreement shows that the epCUT works for a system with equidistant spectrum in H 0 . The plain series is trustworthy up to x ≈ 0.7. For larger x extrapolations are needed. The dotted black line shows the "best", i.e., stiffest, Padé-extrapolant of order [11,6] for this series. Other Padé-extrapolants are shown in Fig. 4. The result of the deepCUT are depicted as solid black line in Fig. 7. It fits perfectly to the perturbative result for weak leg couplings. For larger values of x, it serves as an excellent extrapolation of the perturbative results. In order to support that the directly evaluated results are quantitatively reliable even for larger x the deepCUT data is compared with results from sCUT 40,51 and from DMRG [64][65][66] . The sCUT results, represented by the solid light gray (orange) line, are calculated with the ground state generator and the truncation scheme d= (12,10,10,6,6,5,5,4,4).  The results of these methods agree perfectly. The deviations between the DMRG results and the results of the direct evaluation are shown in the lower inset of Fig. 7. They increase with rising x, but they remain still small. For x = 1.5, the deviation is less than 10 −3 J ⊥ and for x = 3 it is still less than 10 −2 J ⊥ .

B. Dispersion
The one-triplon dispersion is calculated up to order 15. The dispersion is obtained by a Fourier transform of the one-triplon sector of H eff . The dispersion reads where t d is a hopping element over the distance d. internal parameter 67 p(x) where ∆(x) denotes the gap. In order to use the above mapping x → p a reliable extrapolation of the gap ∆(x) is needed. This was achieved by dlog-Padé extrapolations which work very robustly for the gap. The dashed dark gray (blue) lines in Fig. 8 represent the plain series in this internal parameter without any further extrapolation. The solid black lines are the results of the direct evaluation. For small x the perturbative and the direct evaluation agree very well. The direct results are again very robust for larger x as well.
To corroborate that the deepCUT results are reliable even for relatively large values of x its spin gap is compared to the one obtained in DMRG [64][65][66] in Fig. 9. The solid black lines represent the deepCUT results. The solid dark gray (blue) line depicts the results of a 1 /L 2 finite-size scaling of the DMRG results based on the ansatz The deepCUT data shows that for larger x a higher order is needed to compute the gap accurately. This is understood on the basis of the correlation length of the system. A larger value of x implies a larger correlation length ξ which is given by v /∆ where v is the spin wave velocity if no gap were present. By construction of the (de)epCUT, the order n defines the range of processes which are still captured. Hence one can expect a reliable result as long as where the lattice constant is set to unity. The velocity v can be estimated by fitting v sin(k) to the maximum of the dispersion. We find indeed that (49) is satisfied up to x ≈ 3 for n = 15. The deepCUT curve for n = 15 agrees well with the DMRG-results. The deviations shown in the inset are rather small. For x = 2, it is below 10 −2 J ⊥ . Furthermore, the dispersions shown in Fig. 8 agree with exact diagonalization results 45 . We conclude that the reliable results for the uniform spin ladders beyond x = 1 illustrate the efficiency of the deepCUT approach for a model with an equidistant spectrum.

C. Spectral weights
Here we use the transformation of self-adjoint observables O (cf. Sec. III G) to address the issue of spectral weights. We denote the subspace spanned by the states with q quasi-particles by QP q . As in previous work 26,27,30,56 we split the total spectral weight at zero temperature into its contributions from the different subspaces QP q where O q p stands for the sum over all terms of the transformed observable consisting of q creation operators and p annihilation operators in normal-ordering. The state |0 denotes the vacuum state of the effective model, i.e., the ground state of the Hamiltonian.
If the subspaces QP q have been separated by the CUT, i.e., the effective Hamiltonian does no longer mix them, the spectral weights defined by (50) coincide with the ones defined previously 26,27 . The spectral weights correspond to the integral over momentum and and frequency of the corresponding dynamic structure factor S q (k, ω) where the subscript q denotes the contribution of the subspace QP q . Thus, separate sum rules exist for each QP q . Such a split-up is only possible because the dynamics does not mix the subspaces according to the above assumption. We recall that the dynamic structure factors encode the response of various inelastic scattering experiments.
If the subspaces QP q are not or not all separated the equal-time definition (50) is still well-defined. But I q can no longer be interpreted as the sum rule of S q (k, ω) because the subspaces mix in course of the dynamics induced by the Hamiltonian. Nevertheless, the values I q provide a plausible measure of the importance of the subspaces of different number of excitations. A large spectral weight for low numbers of quasi-particles indicates  that results of scattering experiments can be understood from the spectral densities involving only low numbers of quasi-particles.
In the present work, we concentrate on the spectral weights for the observables to illustrate the transformation of observables. The observable O I induces a local spin S = 1 excitation which can be studied experimentally by inelastic neutron scattering 28,29,44 . The observables O II and O III induce S = 0 excitations which can be studied by optical probes, e.g., Raman scattering 25 or infrared absorption 24,43 , in polarizations parallel and perpendicular to the ladder, respectively. Because triplons are S = 1 states, both observables O II and O III induce no contributions in the one-triplon channel: I 1 = 0. The calculation of the corresponding spectral densities S q (k, ω) is left to future studies.
Since the description in terms of triplons on rungs is obviously best for low values of x we expect that more and more triplons need to be addressed upon increasing x. To assess the relative importance of different triplon channels we introduce the relative weights I q,rel = Iq /Itot. They can be calculated using the sum rule for the total spectral weight I tot . For the observables defined in (51), the total weights are given by with the variables where we use the ground state energy per rung E 0 . We focus on the spectral weights in the first fourtriplon (four-quasi-particle) channels I 1 , I 2 , I 3 , and I 4 up to large values of the relative leg coupling x = 3 using deepCUT. In this region, a complete decoupling of all subspaces using η pc or η 4 is no longer possible because divergences occur in the numerical integration of the flow. This problem is well-known from sCUT; it stems from the overlap of continua of different quasiparticle number 40,51,52 . In this situation, the sorting of quasi-particle spaces ascending by energy is no longer possible. In the perturbative evaluation of epCUT, no divergencies appear if the energies in H 0 are separated and indeed ordered according to ascending number of quasiparticles because the hierarchy in Eq. (26) precludes any feedback of high order coeffients to low order coefficients. In the alternating spin ladder the epCUT based on the quasiparticle number must be modified for y ≥ 3. But we stress that this reflects a more sophisticated physics which must be considered in the choice of the generator. It does not represent a conceptual problem of epCUT.
To avoid convergence problems due to overlapping continua in deepCUT, we aim only to decouple subspaces with at most 2 quasi-particles using η 2 while keeping monomials linking subspaces with higher number of quasi-particles, for instance QP 3 ↔ QP 5 . As a consequence, the observables are transformed to a quasiparticle basis where 3 and 4 quasi-particle states still couple to other subspaces.
The spectral weights for the S = 1 observable O I depicted in panel (a) of Fig. 10 agree well with pCUT results 20,27,53 for small values of x. Note that only the one-and the two-quasi-particle channel can be compared quantitatively because the pCUT separated also the three-and four-quasi-particle subspaces, but the present calculation does not.
For O I , most weight is concentrated in the first two quasi-particle channels. Even at x = 3, the one-triplon channel still contains 57.9% of the total weight. The relative weight of the two-and three-triplon channel rises up to 35.2% and 11.3%, respectively. The four-triplon weight remains negligible. The sum rule is slightly violated because the accumulated relative weights exceed 100%. This inaccuracy is related to the finite order of calculation. The degree of the violation of the sum rule can be used as a measure for the reliability of the results. Even at x = 3, the excess weight is only 5.5%.
Panel (b) of Fig. 10 shows the spectral weights I q for the S = 0 observable O II . Since triplons have spin S = 1, there cannot be any weight in the one-triplon channel. Instead, most weight is concentrated in the two-triplon channel which agrees well with pCUT results 20,26,27,53 . Compared to O I , the three-triplon channel is much more pronounced displaying a relative weight of 44.4% at x = 3. At this value, the sum rule is fullfilled within 6.7%.
The observable O III is symmetric with respect to the ladder's centerline. This implies that this observable does not change the parity of a state 26,27 . A single triplon is an odd excitation with respect to the ground state. Hence O III can create or annihilate triplons only in pairs. There is no weight in odd channels in Fig. 10(c). As a consequence, the spectral weight is distributed over the twoand four-triplon channels only. Our results do not indi- cate a sizeable contribution from six and more triplons, but this has not been studied quantitatively. At x = 3, the sum rule is violated by 7.5%. Indeed, this violation sets in at about x = 0.6 when the four-triplon weight becomes significant. Thus we presume that the latter is a bit overestimated, but we could not identify the mechanism for this effect. The perturbative results for the weights fulfill the sum rule to the required order.

A. Ground State Energy
For J o ⊥ = J e ⊥ , the ground state energy is calculated up to order 16. Due to the doubled unit cell, only a slightly lower order can be reached than for the uniform spin ladder. Roughly, we need double the number of coefficients for the alternating spin ladder. The ground state energy per rung is given by E 0 = h0 /2.
The perturbative results from epCUT are shown in Fig.  11. As expected the ground state energy decreases upon rising y. The black lines represent the results of the plain series for various y. The coefficients are given in Tab. IV. The light gray (orange) lines correspond to various Padéextrapolants. The plain series is reliable up to x ≈ 0.75 for y = 1 and up to x ≈ 0.85 for y = 1.5. So the x up to which the series is reliable depend on the value of y. Since a larger value of y supports the dominance of the unperturbed Hamiltonian H 0 it is clear that an increasing y supports the validity of the perturbation. The dark gray (blue) lines represent the results of the deepCUT. These results again represent a very robust extrapolation of the perturbative results up to larger x.
To show the efficiency of the epCUT the results for the ground state energy per rung are compared to the results of an sCUT calculation and a DMRG 64-66 calculation. The sCUT was performed with the ground state generator and a d= (12,10,10,6,6,5,5,4,4) truncation. The DMRG results (L = 20, 40, . . . , 100, m = 100) are extrapolated according to Eq. (45). Figure 12 compares the results of the various approaches. They agree very well with one other. The deviations between the results of the direct evaluation and the DMRG calculation are small, see upper inset, are small. For x = 1.5 the deviation is less than 10 −3 J ⊥ .

B. Dispersion
The dispersion is calculated up to order 13 for the alternating ladder. An important step is the Fourier transform of the hopping in the one-triplon sector of H eff . But the doubled unit cell has to be taken into account which halves the Brillouin zone. In return, the dispersion acquires two branches reading where M ij stands for the Fourier transform of the hopping processes from a rung of parity i to a rung of parity j. Figure 13(a) displays the dispersions for x = 0.5 and various values of y. The solid lines represent y = 1, the dashed ones y = 1.2 and the dotted ones y = 1.4. The dark gray (blue) lines stand for directly evaluated results and the black lines for the perturbatively evaluated ones. For x = 0.5 the plain series is used. Both results agree very well. For y = 1 we retrieve the uniform ladder and the two branches meet at k = π /2. As expected the branches split at k = π /2 once y > 1 holds due to the reduced translational symmetry.
To show the efficiency of the epCUT, the dispersion relations for y = 1.2 is compared to the dispersion from an sCUT calculation in Fig. 14. The sCUT was performed with the generator η 1 and the truncation d= (12,10,10,6,6,5,5,4,4). The dispersions match perfectly. The deviation between sCUT and the perturbative evaluation is less than 10 −3 J ⊥ . The differences in the upper branch are larger than the ones in the lower branch.
Furthermore, the gap for y = 1.2 is compared to the gap obtained by a DMRG calculation [64][65][66] . The finitesize scaling is carried out again based on Eq. (48). In Fig. 15 the solid black line shows the perturbative result and the dashed black line the result of the deepCUT. The result of the DMRG calculation is depicted as solid dark gray (blue) line. The deviations between the deepCUT and the DMRG calculation are shown in the inset. Again the results agree very well, e.g., the deviation is less than 10 −2 J ⊥ even at x = 2 .
The dispersions at x = 1 of the direct and of the perturbative evaluation are plotted in Fig. 13(b). The perturbative result are rendered using the plain series in an internal parameter 67 . The parameter, however, defined in Eq. (47) does not work because at y = 1 it behaves like p(x) ∝ x 2 and not linearly in x anymore. Thus the series in x cannot be re-expressed in a series in p. So we modify the internal parameter where all matrix elements are taken at vanishing wave vector k = 0. We choose this parameter because it reproduces the previous definition (47) for y = 1 and for x → ∞. In addition, it fulfills p a ∝ x for x → 0 for all values of y. Otherwise, the extrapolation can be performed as before 67 . The Fourier transformed matrix elements M ij (x) are obtained by robust dlog-Padé extrapolations. The results of deepCUT and of the series in this internal parameter agree very well.
The epCUT results for the alternating ladder exemplify the efficiency of this CUT for a system with a nonequidistant spectrum in H 0 . Thereby, the range of applicability of perturbation by CUTs is crucially enhanced because the previous pCUT 18,19 is restricted to equidistant unperturbed spectra.

VI. CONCLUSIONS
In this article, we presented a methodical development and illustrated it for a well understood model. We extended the previously known perturbative continuous unitary transformation (pCUT) in two ways.
First, we formulated the perturbative realization of the CUTs directly in second quantization. Thereby, the un- perturbed part is no longer restricted to an equidistant spectrum of energy eigen values. The direct expansion of all coefficients in the effective Hamiltonian is not efficient enough. But tracking the powers in the expansion parameter x of all the physical processes it is possible to identify the relevant ones for the low-energy effective model: Ground state energy, single quasi-particle dispersion, and two-quasi-particle interactions. We could show that this leads to an efficient and competitive approach to obtain effective models. Their parameters are computed as series in the expansion parameter. For distinction, we baptized the enhanced approach enhanced perturbative CUT (epCUT). Second, we found that the system of differential flow equations, which has been reduced to provide the perturbative series representation of the effective model, can also be directly evaluated. It appears that this directly evaluated perturbative CUT (deepCUT) yields a very robust and reliable way to exploit the information in the perturbative differential flow equations. In some sense, one can think of it as a robust extrapolation though we stress that it is not an algorithm applied to a series. The deepCUT provides the parameters of the effective models for given initial Hamiltonian. Each set of initial parameters requires a numerical integration of the flow equations which is a moderate numerical task. The essential effort lies in deriving the system of differential flow equations which is the same as for the epCUT.
The two abstract key results were illustrated by calculations for antiferromagnetic S = 1/2 spin ladders with two legs. Two types of spin ladders were studied. The expansion parameter is the leg coupling relative to the (smallest) rung coupling. The uniform spin ladder with the same value of the rung coupling is the standard model which is very well studied. The alternating spin ladder with alternating rung couplings has not yet been studied to our knowledge. For the present purposes, it constitutes a model with a non-equidistant unperturbed spectrum if the perturbation is set up around the rung Hamiltonian.
For the uniform spin ladder, the known series coefficients could be retrieved by epCUT. The corresponding results for the alternating spin ladder have not been published elsewhere. They show that general unperturbed spectra can be treated by epCUT.
The data obtained by deepCUT illustrate that this approach yields surprisingly robust results. The uniform spin ladder could be treated up into the strong-leg limit with values of the relative leg coupling x = J /J ⊥ of up to x = 3. This is a parameter regime which was not accessible by CUTs before 27,28 .
The limit of the applicability of deepCUT can be understood in terms of the correlation length. A deepCUT calculation in order n in a perturbation linking adjacent sites allows us to capture processes up to the range n · a where a is the lattice constant. Hence reliable results can be expected if the correlation length ξ = v/∆ is lower than n · a. The deepCUT results for the alternating ladder are also very robust, though a little less than for the uniform ladder.
The idea behind the deepCUT approach bears some similarity to the graph-based CUT (gCUT) recently introduced by Yang and Schmidt 60 . In the gCUT the irreducible contributions of subclusters of increasing number M of sites of the lattice are summed. The number M plays a similar role as the number n in deepCUT. It restricts the maximum range of the tracked physical processes. These processes are treated in infinite order in the ratios of the coupling so that the evaluation is non-perturbative. Noticeable differences between deepCUT and gCUT are the following.
The gCUT approach performs calculation on finite clusters and can be done on the level of states in the finite-dimensional Hilbert space, assuming that the Hilbert space of a single site is finite. Only the Hilbert space of the clusters with up to M sites needs to be considered. Translational invariance is lost on the level of the single cluster which limits the choice of applicable generators.
The deepCUT approach works on the level of monomials of creation and annihilation operators, i.e., in second quantization. Thus essentially all symmetries of the lattice problem under study can be preserved by construction. A large variety of the generators can be realized. A detailed comparison of deepCUT and gCUT is an interesting issue left for future research.
In a nutshell, we advocate two approaches to derive effective models in a systematically controlled way in this article. They have been illustrated for spin ladders and we expect that applications to many other models will soon be possible.

ACKNOWLEDGMENTS
We gratefully acknowledge many useful discussions with Sebastian Duffe, Tim Fischer, Mohsen Hafez, Carsten Raas and Kai P. Schmidt. One of us (NAD) is funded by the NRW-Forschungsschule "Forschung mit Synchrotronstrahlung in den Nano-und Biowissenschaften".

Appendix A: Extended A-Posteriori Simplification Rule
The upper bound O max for the maximum order can be reduced by considering the real-space structure of the monomial. For clarity, we restrict ourselves to onedimensional models. As for the basic simplification rule, we discuss the effect of commutations with first-order terms present in the initial Hamiltonian. This is sufficient because any more complicated monomials in the generator have been induced by commutations of a number of first-order terms. Hence their gain in number of involved quasi-particles is paid for by a correspondingly higher order in x. Thus one may safely restrict the consideration to the basic building blocks present in the initial Hamiltonian.
For the spin ladder in terms of triplon operators (Eq. (9d)), a commutation with the generator η (1) cancels at most two local creation or annihilation operators on adjacent sites. Therefore, sparse and extended monomials require more commutations in order to reduce their local operators compared to monomials with the same numbers of operators which are more localized in real space. Figure 16. Decomposition of the sites with creation operators (or the annihilation operators, respectively) of a monomial into linked subclusters ki and its covering with first-order generator terms. Each circle stands for a rung of the spin ladder. Filled circles represent rungs where the local action of the monomial differs from identity. At most two adjacent local operators can be cancelled by a single commutation with η (1) ; this is represented by ellipses.
At first, we study the ground state energy per rung, i.e., the coefficient of the identity operators summed over all rungs, in highest order. The clusters of the creation and of the annihilation operators are treated separately. Both are decomposed into linked subclusters of size k c i and k a i , see Fig. 16. To cancel all local operators, each subcluster needs to be covered by ki 2 first-order generator terms. In conclusion, commutations with η (1) are needed for the clusters of creation or annihilation operators to be reduced to the coefficient of the identity operator. This argument leads to the extended upper bound for the maximum order For a single linked cluster, this formula resembles the result obtained for the basic simplification rule (38). The formula (A2) can be generalized to for the subspace QP q of states with q quasi-particles leading to modified cluster sums K q . Let q be the number of the targeted subspace with the highest number of quasiparticles. This means that q ist the maximum number of local creation and annihilation operators allowed in a monomial targeted up to order n. Terms which affect more quasi-particles have to be reduced to affecting at most q quasi-particles by commutations with η (1) until at most q local creation and annihilation operators are left. To obtain an upper bound O max , one has to choose q positions for local operators to be kept in the cluster in such a way that the other creation and/or annihilation operators can be cancelled by a minimum number of commutations. To this end, one also has to consider that the commutations with hopping terms stemming from H (1) may also shift creation and/or annihilation operators so that they form adjacent pairs which can be cancelled by pair creation or annihilation. But it turns out that this mechanism can reduce the cluster sum K 0 at most by unity while the elimination of a pair of adjacent local operators always reduces the cluster sum by unity. Hence the latter process dominates and provides the correct upper bound O max .
For the hopping in the symmetric ladder model, the above approach means to select sites at the edges of odd subclusters first. This saves one commutation for each local operator kept. Let α be the number of odd clusters. The cluster sum K 0 is reduced in this way by d 1 = min(α, q). (A4) If more local operators remain, i.e., α < q, the most efficient way to place them is in pairs on even subclusters. This reduces the cluster sum additionally by where y is the largest integer which is still smaller or equal to y.
In conclusion, the cluster sums are reduced when one is aiming at higher quasi-particle subspaces according to To avoid unreasonable negative results, this expression has to be checked against zero to obtain the final result K q = max(K q , 0).
We remark that the extended simplification rule can be easily adapted to other models with monomials of first order in the generator create or annihilate an arbitrary number ∆QP of quasi-particles on adjacent sites. For further refinements of O max , one may the triplon polarizations x, y, z as well. But the derivation and application of an appropriate polarization-sensitive simplification rule is beyond the scope of this article which aims primarily at the proof-of-principle demonstration of ep-CUT and deepCUT.
separately. Here we discuss T D explicitly; DT is treated in the same way. For an analogue of the basic simplification rule (Sec. III E 4), we estimate the minimum numbers of creation and annihilation operators c T D and a T D which can appear in the monomials of the normalordering of T D. We use the numbers c T , c D , a T , and a D from each factor as input. At most s T D = min(a T , c D ) (B2) pairs of local operators can cancel in the process of normal-ordering. Hence it follows Using these estimates in Eq. (38), one obtains an upper bound with q being the number of the targeted quasi-particle subspace. Considering also the inverse product DT , the commutator [T, D] does not need to be calculated while holds.
As an example, we consider the second order calculation (n = 2) given in Tabs. I and II for the ground state energy (q = 0). Calculating ∂ H (2) , the commutator of the monomials But in a calculation of order n > 2 or aiming at a higher quasi-particle subspace q > 0, Eq. (B4) yields higher upper bounds for the maximum order and thus the commutator must be evaluated explicitly. Note that this basic a-priori rule is only sensitive to changes of the quasi-particle numbers. It can not anticipate that the commutator in this example actually vanishes due to other properties of the hardcore algebra of the triplons.

Appendix C: Extended A-Priori Simplification Rule
The real-space structure of the commutator arguments T and D allows us to extend the above a-priori rule in analogy to the extended a-posteriori rule in App. A. Let C T and C D be the clusters of the creation operators in T and in D, respectively. Analogously, A T , A D are the clusters of their respective annihilation operators. Normalordering the product T D can cancel local operators only on the intersection Due to the locality of the triplon algebra, the commutator vanishes if none of the clusters overlap Thus the normal-ordered product T D definitely has local creation operators on the union cluster and local annihilation operators on the union cluster There may be additional creation or annihilation operators, but no general statements can be made on their existence. In this sense, the right hand sides of Eqs. Although it is less strict, the basic a-priori has the advantage to be much more lightweight in comparison to the extended a-priori rule because it requires mere counting of operators. Furthermore, it can be used very efficiently in the context of translation symmetry. Because it does not rely on the real-space structure of a term, it can be applied to all terms in the translation group in contrast to the extended rule. Therefore, for best performance it turns out to be most efficient to combine both rules in practice.