Introduction. Obtaining low-energy effective models (LEMs) that accurately encode the low-energy physical properties of the original models is one of the central tasks in condensed matter physics. The LEMs typically have an exponentially smaller Hilbert space, which provides significant convenience for both analytical and numerical studies. Moreover, emergent interactions exhibited in LEMs can offer a clearer and more straightforward understanding of low-energy physical properties. For example, the pioneering derivation of antiferromagnetic superexchange [1,2] provides an intuitive understanding of antiferromagnetism in Mott insulators. Today in condensed matter physics, perhaps the most commonly used LEMs are various spin models reduced from electronic degrees of freedom [3], which have achieved great success in describing low-energy excitations in materials [4,5].

Typically, the derivation of LEMs utilizes a unitary transformation to decouple the designated low-energy sector from the remainder of the Hilbert space, often called a canonical transformation [6]. An anti-Hermitian generator S transforms the Hamiltonian via H eff = e -S He S , which block diagonalizes the Hamiltonian, with the low-energy sector forming the LEM. When there is an energy gap separating the two sectors, which is large compared to the couplings between them, the transformation can be carried out perturbatively. One celebrated example is the Schrieffer-Wolff transformation of the Anderson model into the Kondo model [7]. Another wellknown case is the Hubbard model [8][9][10]:

with ij denoting nearest-neighboring sites, which can be reduced to a spin-1/2o rt-J model in the large U/t limit by integrating out double occupancy [11,12].

In contrast to various types of canonical transformations [6,[13][14][15][16][17][18] that directly transform Hamiltonians, another frequently used but less systematic approach, facilitated by the recent development of numerical simulation, derives LEMs by matching their low-energy wave functions with the original models [19][20][21][22]. This can be useful when the canonical transformation is difficult to carry out or where it yields a large number of significant terms, provided that one only aims to match a few low-lying states instead of the whole lowenergy sector. For example, in one approach [19], an effective Hamiltonian of simpler form is found that reproduces just the ground state rather than matches a spectrum. One first obtains the ground state of the original Hamiltonian, then designs an effective Hamiltonian ansatz H ( α) parametrized by α, and uses multiple simulations to find the optimal α that maximizes the overlap of the new and original ground state. However, the optimization cost can be high if the ansatz contains many parameters and a large number of simulations are required. In some cases, the original Hamiltonian simulation may identify degrees of freedom as being unoccupied, allowing a very simple truncation. For example, in a recent work on the three-band Hubbard models, density-matrix renormalization group (DMRG) simulations revealed that many one-particle orbitals were almost completely unoccupied, allowing them to be truncated without the need for two-particle unitary transformations [23]. The LEM was then formed by a simple single-particle rotation to the Wannier functions spanning the occupied orbitals. However, in the general case, when the higher bands have finite occupancy or when the high-energy degrees of freedom exhibit many-body effects such as significant double occupancy, a many-body unitary transformation is required to derive LEMs.

In this Letter, we propose a systematic and efficient method to obtain LEMs based on a unitary transformation of the ground state. The algorithm is based on tensor networks, which find a proper unitary circuit G made of local gates that transforms the ground state of the original model to a projected wave function with only low-energy degrees of freedom. The effective model can then be generated via H eff = GHG † , whose terms are local and few-body by the design of the unitary circuit, with their coefficients insensitive to system sizes. It is highly efficient since no more simulations other than obtaining the ground state of the original model are required. It is also nonperturbative and derives more accurate LEMs compared to the standard perturbative approach for intermediate coupling. The validity of the generated LEMs can be verified by simulation of their low-lying states.

Algorithm. The algorithm begins by using DMRG [24][25][26] to simulate the original Hamiltonian H 0 to obtain its ground state |ψ 0 . We then project out certain high-energy degrees of freedom and get |φ P = P|ψ 0 (with normalization). In the Hubbard model, which is our main focus, we take P as the Gutzwiller projection [27], which projects out double occupancy. In a multiband model, one could choose a projection operator that removes higher bands of single-particle states. The goal of the projection is to identify the target wave function with only the effective low-energy degrees of freedom. Consider a two-site half-filled Hubbard model as an example. Its ground state is |ψ 0 =c 1 (|↑, ↓ -| ↓, ↑ ) + c 2 (|0, ↑↓ +|↑↓, 0 ). We would like the ground state of the LEM |φ P to omit the doubly occupied states and retain the same singlet form as the original |ψ 0 , i.e., |φ P ∼c 1 (| ↑, ↓ -| ↓, ↑ ). Therefore, |φ P is just the Gutzwiller projection of |ψ 0 , where the spins are the effective degrees of freedom. We note that the primary use of the projected state is to find an appropriate unitary transformation G satisfying G|ψ 0 =|φ P . In particular, we are not interested in the energy of the projected state using the original Hamiltonian. The LEM is generated via H eff = GHG † , which guarantees its ground state to be |φ P .

Since |ψ 0 and |φ P are only two states in the exponentially large Hilbert space, there are many unitary transformations Ĝ that satisfy Ĝ|ψ 0 =|φ P while transforming other states differently. To obtain a Ĝ that can produce a useful H eff with local and few-body terms, it needs to be constructed with local and few-body unitary transformations as well, Ĝ = i ĝi . Therefore, a natural ansatz for Ĝ is a unitary circuit (also known as a quantum circuit). In this paper, we have considered two different structures of unitary circuits, the sequential circuit and the brick wall circuit, shown in Figs. 1(a) and 1(b), respectively. The sequential circuit turns out to be more useful as it can achieve a higher fidelity between the transformed state and the target state (with a similar number of gates), potentially due to its capacity to generate long-ranged correlations [28].

After the structure of the circuit is chosen, we optimize the unitaries to maximize φ P | Ĝ|ψ 0 . Here, we optimize one unitary at a time with the other unitaries fixed, and sweep through the whole circuit back and forth, as shown in Fig. 1(c). This "DMRG-like" local optimization has also been studied in various works with good performance [28][29][30].

At each step, the optimal unitary is obtained using the Evenbly-Vidal algorithm [30], but with additional steps to avoid encoding unnecessary transformations in the unitaries. As shown in Figs. 1(c) and 1(d), we start by contracting the whole circuit with the two states except the desired gate, to obtain the environmental tensor E . We perform a singular value decomposition (SVD) on E :

FIG. 1. Two structures of the unitary circuits considered in this work: (a) a sequential circuit and (b) a brick wall circuit. (c) To find the local unitaries that transform |ψ original to |φ projected , we optimize one unitary at a time with the others fixed and sweep through the whole circuit back and forth. The optimal individual unitary is constructed based on the environmental tensor E , which is obtained by leaving out the corresponding unitary in the contraction of the circuit with wave functions. (d) To optimize each unitary, one performs an SVD of the tensor E , adjusts the tensor V † into Ṽ † , and recombines with the tensor U to form the optimal unitary gate.

The unitary ĝ maximizing the overlap is constructed by removing the singular values S i :

Note that for ĝ to be a unitary, both {|U i } and {|V i } need to be a complete orthonormal basis set. However, not only are there vanishing individual singular values, but there are entire blocks that vanish. For example, if |φ P has double occupancy projected out, then the relevant blocks in E will be zero. Any vectors {| Ûj } and {| Vj } with zero singular values are not uniquely determined by the SVD in Eq. ( 2). While these different unitaries can achieve the same overlap (which is solely determined by the unique part in SVD with finite S i ), the derived LEMs will be different. To maximally preserve the locality and few-body nature of the original model, we need to minimize the rotation encoded in this flexible part of the unitary. To achieve this, the {| Ûj } and {| Vj } need to be rotated to match as much as possible. Inspired by the procedure in generating Wannier functions [23], for each |U j in the zero-singular value space, we make the linear combination of the

then symmetrically orthonormalize [31] {W j } to obtain {| Ṽj }. This subset updates the complete basis resulting proper unitary is

In practice, we set a cutoff of 10 -5 for S i below which we implement the above scheme for the corresponding | Ûi and | Vi . This has a minimal effect on the overlap but can avoid some unnecessary transformations being encoded in the unitaries.

We find that a modest number of sweeps [5-10, see Fig. 2(b)] are sufficient to make the transformed state converge to the target state. Then, the LEM can be generated with H eff MPO = GHG † . As a tensor network, this is a matrix product operator with additional unitary gates on both sides, which can be simplified to a single-matrix product operator.

We seek to convert the matrix product operator (MPO) H eff MPO to an analytic form as a sum of local terms, which we call H eff ( h) with h denoting the magnitude of the terms. There are subtleties associated with this process. For example, let us consider extracting a spin exchange between sites 1 and 2 in a fermionic model. One way is to contract the MPO tensors from sites 3 to N with an identity operator, which projects out the terms involving other sites. However, the above procedure does not work if the MPO has a term S z 1 S z 2 n 3 , which is essentially the same as S z 1 S z 2 in the low-energy spin space ( n 3 =1) of the half-filled Hubbard model we will be dealing with later. Therefore, it makes sense to extract the Hamiltonian terms under a basis that resides in the low-energy sector, i.e., to set a spin background, which can be done in an unbiased way using perfect sampling [32] of the projected ground state |φ P to obtain a product state of spins.

With the appropriate background, we can proceed to extract the terms. The effects of different terms can be mixed together. For example, both J 1 and J c [defined in Eq. ( 6)] can cause a pair of spin flips. Therefore, we evaluate H eff MPO on a cluster instead of a single bond, such that multiple terms are supported and considered simultaneously. The detailed procedure is as follows: (1) obtain a product state |α by sampling the projected ground state |φ P , (2) choose a cluster (e.g., a 2 × 2 plaquette) and construct a basis {α i } by listing all the spin configurations of cluster with the background spins outside unchanged as |α , (3) compute all the matrix elements of α i ′ |H eff MPO |α i , and (4) find the analytic H eff ( h) that reproduces those matrix elements:

. This procedure is repeated with multiple samplings, as well as varying the size and position of the selected cluster. If all of them yield the same H eff ( h), then the conversion from H eff MPO to an analytic form H eff ( h) can be confirmed.

An example of converting H eff MPO to H eff ( h) is presented in Fig. 3. Here, the sampled product state |α is the Néel state, based on which the basis {α i } is constructed by changing the spin configuration of the central 2 × 2 plaquette. Figure 3(a) shows that all the matrix elements of H eff MPO under this basis match well with a Hamiltonian containing four types of spin exchanges H eff (J h , J v , J 2 , J c ) [defined in Eq. ( 6), J h and J v are J 1 interaction on horizontal and vertical bonds, which are approximately the same), with an average difference of only 0.0059. Also, the extracted h = (J h , J v , J 2 , J c ) only varies slightly between different samples as shown in Fig. 3(b).W e further checked different cluster shapes and found no significant longer-ranged term. The extracted exchanges also do not depend on the position of the selected cluster. These evidences confirm that the extracted analytic H eff (J h , J v , J 2 , J c ) is a good approximation of the original H eff MPO from the unitary circuit. Results. We first test our algorithm on a half-filled 12-site Hubbard chain at an intermediate coupling U = 6( t = 1i s set as the energy unit throughout the paper). The results are presented in Fig. 2. The simple sequential-1 circuit made of two-site gates only already reduces the infidelity 1 -F := 1 -| φ P |G seq1 |ψ 0 | 2 to O(10 -3 ) within two sweeps. To improve upon it, we have tried to expand the circuit with two-site gates, either by replacing each gate with three sequential gates or by adding another stack. However, these structures generally do not offer a substantial improvement, suggesting the limitation of local two-site gates. Therefore, a further improved circuit, namely, the sequential-2 circuit, adds a stack of four-site gates on top of the sequential-1 circuit to deal with longer-ranged and four-body interactions [34]. When optimizing the four-site gates, the lower stack of two-site gates

3. Converting the numerical H eff MPO to an analytic H eff ( h)f o rt h e2× 12 Hubbard ladder at U = 6i nF i g .4 (t = 1 as the energy unit). (a) The 70 nonzero matrix elements of α i ′ |H eff MPO |α i and their corresponding α i ′ |H eff (J h , J v , J 2 , J c )|α i , sorted in ascending order. The close match indicates H eff MPO can be well described by these four spin exchanges with almost nothing left out. All the diagonal matrix elements α i |H eff MPO |α i and α i |H eff (J h , J v , J 2 , J c )|α i are shifted by a rough total energy estimate ∼-N (J h + J v )/2 to fit in the plot. {|α i } is a basis set constructed from one sampled product state |α (see text). (b) The extracted (J h , J v , J 2 , J c ) for each sample and the average values ( Jh , Jv , J2 , Jc ) over all samples.

is fixed to avoid absorbing the two-site transformation into the four-site ones. On the other hand, we have also tried the brick-1 circuit, which has one layer of two-site gates on odd bonds followed by another layer on even bonds. Its infidelity is one magnitude bigger, possibly because of its inability to generate long-range correlation, similar to the observation in Ref. [28]. Furthermore, the H eff it yields differs between even and odd bonds. For these reasons, we consider the sequential circuit to be a better choice.

For the resulting H eff ( h) for the one-dimensional (1D) chain, the only significant interaction is the nearest-neighbor (NN) Heisenberg exchanges J 1 , while the next-NN and fourbody exchange coefficients are 0.01. This is verified by a simulation of the Heisenberg chain, whose ground state is almost identical to that of the original Hubbard ground state projected into spin space, with an infidelity of only 1.6 × 10 -4 . Alternatively, H eff MPO as an MPO has a bond dimension of 16, since it also contains interaction involving doubly occupied or empty states. After projecting out those states and applying a cutoff of 1×10 -4 , its bond dimension reduces to 5, the same as a Heisenberg chain. In Fig. 2(c), we show the extracted J 1 from different circuits, which converges with the increase of the fidelity. Matching the spin gap between the Hubbard chain and the Heisenberg model also yields J 1 ≈ 0.60. This is expected since all the energy levels are preserved in a unitary transformation of the Hamiltonian.

Next, we test our algorithm on a two-leg Hubbard ladder. The sequential circuit is shown in Fig. 4(a). Its structure resembles a staircase: starting from the left end of the system, a gate is acted on the rung, followed by two gates acting on the legs, and this pattern continues to the right. The infidelity is further reduced by an order of magnitude by adding another stack of four-site gates on top of the two-site gates, as shown in Fig. 4(b).

For the two-leg ladder, there are four types of spinexchange terms (J 1 , J 2 , J 3 , J c ) according to the fourth-order t/U perturbative expansion: [11,12]. We first extract those acting within a two-by-two plaquette as illustrated in Fig. 4(c). We also checked J 3 and other long-ranged terms, which are found to be less than 0.02. The form of the approximate effective Hamiltonian H eff (J 1 , J 2 , J c )is

where ij 1(2) denotes the first(second) nearest-neighbor pair of sites and ijkl denotes a 2 × 2 plaquette. The process of extracting spin exchanges is shown in Fig. 3 as mentioned previously. Averaging the NN horizontal exchange J h and vertical exchange J v to J 1 , the resulting (J 1 , J 2 , J c ) at three different couplings U are listed in Fig. 4(e), and compared with those obtained from the fourth-order t/U expansion. As expected, the parameters match well at strong coupling U = 12. However, at intermediate couplings U = 8 and U = 6, while the J 1 remains similar for both approaches, the four-spin cyclic exchange term J c is much smaller in H eff . To test the validity of the obtained H eff (J 1 , J 2 , J c ), we simulate it using DMRG and compute the overlap of its ground state with the original projected Hubbard wave function. As shown in Fig. 4(e), H eff (J 1 , J 2 , J c ) produces essentially the same ground state (infidelity <0.01) for all three cases, while the t/U expansion becomes much less accurate as U decreases.

We further tested our method in two more challenging systems. The first one is a 4 × 8 cylinder at half-filling with U = 8. In this case, with a fidelity of 0.977 between the transformed wave function and the targeted one, the J c extracted from H eff is only around half of that from the U/t expansion, similar to the two-leg ladder case. We note that the spin exchanges obtained (J 1 = 0.437, J 2 = -0.004, J c = 0.083) on a four-leg cylinder is similar to that in the two-leg ladder case, indicating the validity and generality of the effective model derived on relatively small cluster. Secondly, we apply our algorithm to a doped two-leg ladder with U = 7.5 and doping = 0.125 and obtain an effective t-J model. In addition to the various spin exchanges, we extract a three-site correlated hopping term [35] with t 3-site /t = 0.11, comparable to the t 3-site /t = 0.13 from the second-order t/U expansion. However, in the doped case H eff is not as well approximated by these standard terms as the undoped case. Although these cases are more challenging, the results support the algorithm's applicability to both two-dimensional and doped models.

Discussions and summary. In the derivation of LEMs, ideally, one wishes to find an effective Hamiltonian that fully captures a certain low-energy sector of the original model in terms of both energy levels and the corresponding states. The unitary transformation preserves the energy levels. Although in our algorithm there is no guarantee that the excited states will also match, in our test on the two-leg half-filled Hubbard ladder at U = 6 (Fig. 4), the first excited state also matches with nearly the same overlap as the ground state. One future improvement is to require the unitary circuit to transform not only the ground states but also several low-lying excited states, potentially utilizing state-averaged matrix product state.

Regarding the structure of the unitary circuit, we have explored removing the two-site gates and using only four-site gates for the circuit in Fig. 4. While a similar fidelity between the transformed wave function and the targeted one can be reached, the H eff extracted from the circuit cannot be as well approximated by a simple H eff (J 1 , J 2 , J c ), with the average difference in the Hamiltonian elements around six times larger than that from the sequential-2 circuit in Fig. 4(a). This confirms that the circuit should be constructed primarily using local and few-body gates, with larger gates employed only to transform the remaining parts that are beyond the reach of the smaller ones. This ensures that the resulting H eff remains local and few-body.

We note that an alternative approach that utilizes ground states to derive LEMs is the parental Hamiltonian methods [36][37][38][39]. In this case, one directly optimizes H eff such that given target state becomes its In these approaches, one needs to take an extra step to determine the scale of H eff , which can be done by matching the excitation gap.

In summary, we propose a tensor-network-based method to derive LEMs based on a unitary transformation of the ground state. The algorithm finds a proper unitary circuit G made of local gates that integrates out the high-energy degrees of freedom in the ground state, which is then used to derive the effective model via H eff = GHG + . The resulting effective model is general and largely insensitive to system size. The algorithm is highly efficient and only requires a single simulation to obtain the ground state of the original model. Its effectiveness is tested on the square-lattice Hubbard model at half-filling and shown to generate more accurate LEMs compared to the perturbative approach at intermediate coupling.