Add exponentially decaying coupling $\sum_{j!=i} \alpha \beta^{|i-j|} O_i O_j$

[sajantanand]

Currently there is functionality to add an exponentially decaying term between all sites (possibly of some subset) of the lattice of the form $\sum_i\sum_j \alpha \beta^{|i-j|} O_i O_j$. I would like the option to only sum over one index rather than both. So the goal is to have an exponentially decaying coupling "centered" on site $i$ and coupling site $i$ to all other sites $j$. If one were to have a exponentially decaying Ising Hamiltonian $H = \sum_i \sum_j \alpha \beta^{i-j} Z_i Z_j$, the energy density on site $i$ would be an example of what I'm trying to implement. The use case of this is to do Heisenberg simulations of the energy density to calculate energy transport of long-range interacting models.

I have modified the add_to_graph function of ExponentiallyDecayingTerms (https://github.com/tenpy/tenpy/blob/main/tenpy/networks/terms.py#L1363) to realize this functionality.

[sajantanand]

Also, MPOs with exponentially decaying couplings are not strictly block upper-triangular since the exponentially decaying term is added "after" IdR (on line 1424 of https://github.com/tenpy/tenpy/blob/main/tenpy/networks/terms.py#L1363, you note that you assign it a high label so that it comes at the end). If we want to do backwards solving of MPO environment a la the Ghent group, don't we want a block triangular MPO?

[sajantanand]

As a simple test case, I have tried making a exponentially decaying XY model with $L=10$. The first five sites have coupling prefactor 1 while the last five have coupling prefactor -1.

L = 10
site = SpinSite(S=0.5, sort_charge=True, conserve='Sz')
lat = Chain(L, site)
operator1 = CouplingModel(lat)
operator2 = CouplingModel(lat)
operator3 = CouplingModel(lat)

for i in range(0, L//2):
    operator1.add_exponentially_decaying_coupling(1, 0.1, 'Sp', 'Sm', fixed_site=i, plus_hc=True)
for i in range(L//2, L):
    operator1.add_exponentially_decaying_coupling(-1, 0.1, 'Sp', 'Sm', fixed_site=i, plus_hc=True)
operator1_MPO = operator1.calc_H_MPO()

for i in range(0, L//2):
    operator2.add_exponentially_decaying_coupling(1, 0.1, 'Sp', 'Sm', subsites=np.arange(0, L//2), fixed_site=i, plus_hc=True)
for i in range(L//2, L):
    operator2.add_exponentially_decaying_coupling(-1, 0.1, 'Sm', 'Sp', subsites=np.arange(L//2, L), fixed_site=i, plus_hc=True)
operator2_MPO = operator2.calc_H_MPO()

operator3.add_exponentially_decaying_coupling(2, 0.1, 'Sp', 'Sm', subsites=np.arange(0, L//2), plus_hc=True)
operator3.add_exponentially_decaying_coupling(-2, 0.1, 'Sm', 'Sp', subsites=np.arange(L//2, L), plus_hc=True)
operator3_MPO = operator3.calc_H_MPO()

assert operator1_MPO.is_hermitian()
assert operator2_MPO.is_hermitian()
assert operator3_MPO.is_hermitian()

assert operator2_MPO.is_equal(operator1_MPO)
assert operator3_MPO.is_equal(operator1_MPO)

operator3_MPO has the smallest bond dimension (3) compared to 12 and 7 for the first and second version. This is because the FSM used to construct the MPO knows that all of the exponentially decaying terms can share the same state. However, this demonstrates that the functionality to add exponentially decaying terms without summing over all sites is working.

[sajantanand]

As far as I can tell, this code can be merged. I edited the documentation to describe the new changes (which don't change the default behavior) and fixed an issue with to_TermList().

[sajantanand]

I added a function to tenpy.networks.mpo to create an MPO $\alpha 1 + \beta O$ from the MPO for $O$. This is only implemented for finite MPOs. All one needs to do is modify the left-most tensor to include the amplitudes and a new "onsite operator" (i.e. the identity) into the D block.