VUMPS in a two leg ladder model

ITensor Julia Questions
1

Dear all,

I am currently trying to use the implementation of the vumps algorithm from the ITensorInfiniteMPS.jl
package.

The model I want to study consists of two wires of Josephson Junctions described by
bosonic modes b_{\alpha,j}, \alpha = a,b. The modes are defined w.r.t. the mean filling of the
superconducting islands, i.e. n = diag(\dots, -2, -1, 0 , 1 ,2, \dots).
The Hamiltonian looks similar to:

H = \sum_{j,\alpha} b_{j,\alpha}^\dagger b_{j,\alpha} + \text{h.c.} + E_c\, n_{j,\alpha}^2 + \sum_j + \Omega b_{j,a}^\dagger b_{j,b} + \text{h.c.} + V n_{j, a} n_{j, b}

In previous finite size DMRG studies, I “enrolled” the unitcell to only have a local hilbertspace dimension of 2n_b + 1 with a chain of length 2L. I.e. I defined new operators:

\tilde{b}_j = \begin{cases} b_{(j + 1)/2, a}\ ,\ \text{for } j\ \text{odd}\\ b_{j/2, b}\ ,\ \text{for } j\ \text{even}\\ \end{cases}

The above Hamiltonian gets then transformed into a model with next nearest neighbor hoppings,
and interactions + hoppings only between even and odd sites 2j - 1 and 2j, but not between
2j and 2j+1.

I tried to define a new model with unit_cell terms as:

function ITensorInfiniteMPS.unit_cell_terms(::Model"JosphsonLadder";  Ej = 1, Ec = 1, V = 0, Omega = 0)
    opsum = OpSum()
    opsum += -Ej/2, "Cr", 1, "Anh", 3
    opsum += -Ej/2, "Cr", 3, "Anh", 1
    opsum += -Ej/2, "Cr", 2, "Anh", 4
    opsum += -Ej/2, "Cr", 4, "Anh", 2
    opsum +=  Ec, "N * N", 1

    opsum += V, "N", 1, "N", 2
    
    opsum += Omega, "Cr", 1, "Anh", 2
    opsum += Omega, "Cr", 2, "Anh", 1
    
    return opsum
end

But this is too naive, since it will not recognize that the \Omega term needs to be repeated only every
second site. What would be the correct way to define such a Hamiltonian?
Moreover, for \Omega = 0, the subspace expansion will fail with the error message:

Impossible to do a subspace expansion, probably due to conservation constraints

It seems, that there needs to be a nearest neighbor term for the subspace expansion.

Thanks for any help in advance!

2

Ok, I think I defined the model in a correct way now:

function ITensorInfiniteMPS.unit_cell_terms(::Model"JosphsonLadder";  Ej = 1, Ec = 1, V = 0, Omega = 0)
    opsum = OpSum()
    opsum += -Ej/2, "Cr", 1, "Anh", 3
    opsum += -Ej/2, "Cr", 3, "Anh", 1
    opsum += -Ej/2, "Cr", 2, "Anh", 4
    opsum += -Ej/2, "Cr", 4, "Anh", 2
    
    opsum +=  Ec, "N * N", 1   
    opsum +=  Ec, "N * N", 2

    opsum += V, "N", 1, "N", 2
    
    opsum += Omega, "Cr", 1, "Anh", 2
    opsum += Omega, "Cr", 2, "Anh", 1
    
    return opsum
end

At least this definition gives the same results in a finite system by using MPO(Model("JosephsonLadder"), ...). I think I also found a bug. Namely if there is only one term in the
OpSum acting on one position in the unit cell, the creation will fail. This is because in Line 136 of
models.jl the line

opsum_cell_dict = groupreduce(minimum ∘ ITensors.sites, +, opsum)

will not convert such a single term into a OpSum object, but it will simply stay a Scaled{ComplexF64, Prod{Op}}. I think this will be solved by using

opsum_cell_dict = groupreduce(minimum ∘ ITensors.sites, identity, +, opsum; init = OpSum())

instead.

And is there already an implementation for correlation functions?

Thanks!

3

Ah and I notice that the energy reported by vumps is not equal for both positions in the unit-cell, is this bad or something which can happen?

4

Could you make a PR for that fix for single terms in the OpSum? Pull requests · ITensor/ITensorInfiniteMPS.jl · GitHub

We’re adding an example for computing correlation matrices here: Added example file with Infinite MPS to finite MPS function by nbaldelli · Pull Request #72 · ITensor/ITensorInfiniteMPS.jl · GitHub (the goal is to just use finite MPS tools to do that, to minimize the amount of code I have to support in that library, since I don’t have much time to develop it and it will be superseded by next-generation infinite tensor network tools in GitHub - mtfishman/ITensorNetworks.jl: A package with general tools for working with higher-dimensional tensor networks based on ITensor. once those become available).

Unequal energies would depend on the physics, for example if there is spontaneous symmetry breaking (or the Hamiltonian already breaks the symmetry). Otherwise it could just be due to convergence, it of course depends on details of the calculation and the magnitude of the difference (and how the difference changes as you change the bond dimension, i.e. does it trend towards zero in the infinite bond dimension limit?).

5

Hey thanks!
yes I can do that later.

Ok that approach sounds good, will look into that. Are these tools already on Github in some branch?
If yes I would rather switch to that since it will be supported in the future.

Mhm ok, I guess its because of the additional interaction term between site 1 and 2 not present between
site 2 and 3. On the other hand, all local observables are homogeneous on the unitcell.

Thanks!