Modeling additional system coupled to the bulk

DMRG and Numerical Methods
1

Hi,

I’m wondering if it’s possible to model an additional system coupled to the bulk Hamiltonian in ITensor. More concretely, let’s say the Hilbert space of our system takes the form

|\alpha \rangle \otimes |s_1s_2\cdots s_N\rangle
In the simplest case, s_i are fermion sites with QN conservation, and |\alpha\rangle is a two-level system coupled by some interaction hamiltonian
\langle \alpha s_1\cdots s_N| H_{int} | \alpha' s'_1\cdots s'_N \rangle \sim ( 1 - \delta_{\alpha, \alpha'}) f( s_1, s_1'\cdots )
And there is some diagonal bulk Hamiltonian only acting on the s_i sites. In ED this is easily written as
\begin{pmatrix} H_{bulk} & H_{int} \\ H_{int} & H_{bulk} \end{pmatrix}

My first implementation is to simply define the sites as

  sites = Vector{Index}(undef, L + 1)

  #two level system
  sites[1] = siteind("Fermion"; addtags="n=1", conserve_qns =false)

  for mid = 2: L + 1
    sites[mid] =  siteind("Fermion"; addtags="n=$mid", conserve_qns =true)
  end

It seems even with the simplest bulk/interaction Hamiltonian, it could not be properly transformed into a MPO by calling H = MPO(ampo, sites), as I get the following errors

ERROR: MethodError: no method matching contraction_output(::NDTensors.DenseTensor{Float64, 2, Tuple{Index{Int64}, Index{Int64}}, NDTensors.Dense{Float64, Vector{Float64}}}, ::NDTensors.BlockSparseTensor{Float64, 2, Tuple{Index{Vector{Pair{QN, Int64}}}, Index{Vector{Pair{QN, Int64}}}}, NDTensors.BlockSparse{Float64, Vector{Float64}, 2}}, ::Tuple{Index{Int64}, Index{Int64}, Index{Vector{Pair{QN, Int64}}}, Index{Vector{Pair{QN, Int64}}}})
Closest candidates are:
  contraction_output(::NDTensors.Tensor, ::Any, ::NDTensors.Tensor, ::Any, ::Any) at ~/.julia/packages/NDTensors/c2BpJ/src/dense.jl:602
  contraction_output(::NDTensors.Tensor, ::NDTensors.DiagBlockSparseTensor, ::Any) at ~/.julia/packages/NDTensors/c2BpJ/src/blocksparse/diagblocksparse.jl:290
  contraction_output(::NDTensors.DiagBlockSparseTensor, ::NDTensors.Tensor, ::Any) at ~/.julia/packages/NDTensors/c2BpJ/src/blocksparse/diagblocksparse.jl:287
  ...
Stacktrace:
  [1] contraction_output(T1::NDTensors.DenseTensor{Float64, 2, Tuple{Index{Int64}, Index{Int64}}, NDTensors.Dense{Float64, Vector{Float64}}}, labelsT1::Tuple{Int64, Int64}, T2::NDTensors.BlockSparseTensor{Float64, 2, Tuple{Index{Vector{Pair{QN, Int64}}}, Index{Vector{Pair{QN, Int64}}}}, NDTensors.BlockSparse{Float64, Vector{Float64}, 2}}, labelsT2::Tuple{Int64, Int64}, labelsR::NTuple{4, Int64})
    @ NDTensors ~/.julia/packages/NDTensors/c2BpJ/src/dense.jl:604
  [2] contract(T1::NDTensors.DenseTensor{Float64, 2, Tuple{Index{Int64}, Index{Int64}}, NDTensors.Dense{Float64, Vector{Float64}}}, labelsT1::Tuple{Int64, Int64}, T2::NDTensors.BlockSparseTensor{Float64, 2, Tuple{Index{Vector{Pair{QN, Int64}}}, Index{Vector{Pair{QN, Int64}}}}, NDTensors.BlockSparse{Float64, Vector{Float64}, 2}}, labelsT2::Tuple{Int64, Int64}, labelsR::NTuple{4, Int64})
    @ NDTensors ~/.julia/packages/NDTensors/c2BpJ/src/dense.jl:620
  [3] contract(T1::NDTensors.DenseTensor{Float64, 2, Tuple{Index{Int64}, Index{Int64}}, NDTensors.Dense{Float64, Vector{Float64}}}, labelsT1::Tuple{Int64, Int64}, T2::NDTensors.BlockSparseTensor{Float64, 2, Tuple{Index{Vector{Pair{QN, Int64}}}, Index{Vector{Pair{QN, Int64}}}}, NDTensors.BlockSparse{Float64, Vector{Float64}, 2}}, labelsT2::Tuple{Int64, Int64})
    @ NDTensors ~/.julia/packages/NDTensors/c2BpJ/src/dense.jl:611
  [4] _contract(A::NDTensors.DenseTensor{Float64, 2, Tuple{Index{Int64}, Index{Int64}}, NDTensors.Dense{Float64, Vector{Float64}}}, B::NDTensors.BlockSparseTensor{Float64, 2, Tuple{Index{Vector{Pair{QN, Int64}}}, Index{Vector{Pair{QN, Int64}}}}, NDTensors.BlockSparse{Float64, Vector{Float64}, 2}})
    @ ITensors ~/.julia/packages/ITensors/OjQuG/src/itensor.jl:1864
  [5] _contract(A::ITensor, B::ITensor)
    @ ITensors ~/.julia/packages/ITensors/OjQuG/src/itensor.jl:1870
  [6] contract(A::ITensor, B::ITensor)
    @ ITensors ~/.julia/packages/ITensors/OjQuG/src/itensor.jl:1974
  [7] *
    @ ~/.julia/packages/ITensors/OjQuG/src/itensor.jl:1961 [inlined]
  [8] svdMPO(os::Sum{Scaled{ComplexF64, Prod{Op}}}, sites::Vector{Index}; kwargs::Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}})
    @ ITensors ~/.julia/packages/ITensors/OjQuG/src/physics/autompo/opsum_to_mpo.jl:119
  [9] svdMPO
    @ ~/.julia/packages/ITensors/OjQuG/src/physics/autompo/opsum_to_mpo.jl:2 [inlined]
 [10] MPO(os::Sum{Scaled{ComplexF64, Prod{Op}}}, sites::Vector{Index}; kwargs::Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}})
    @ ITensors ~/.julia/packages/ITensors/OjQuG/src/physics/autompo/opsum_to_mpo_generic.jl:226
 [11] MPO(os::Sum{Scaled{ComplexF64, Prod{Op}}}, sites::Vector{Index})
    @ ITensors ~/.julia/packages/ITensors/OjQuG/src/physics/autompo/opsum_to_mpo_generic.jl:217

As far as I know, QN conservation is a U(1) symmetry that enforces certain ‘flux’ structure on state tensors, so I’m wondering if this has to do with the specifics of the Hamiltonian or it’s just some general incompatibility.

2

Thanks for the question. The error message could definitely be clearer – that is on us. But it is stemming from the mixture of QN-conserving and non-QN-conserving sites you are attempting here, which is not currently something you can do in ITensor.

So what you’ll need to do is to use all QN-conserving sites.

To know how to help you more with that, we’d need to talk more specifically about what your Hamiltonian is, because yes QN conservation is exactly about what symmetries your Hamiltonian has (globally, across the whole system). If even one of the terms does not respect the U(1) symmetry then you won’t be able to use it at all for your system. On the other hand if your whole Hamiltonian does respect U(1), then you should also set conserve_qns=true for your first site.

3

Hi Miiles,

Thank you for the reply.

The system we have in mind is as follows, the ‘bulk’ system is a regular fermionic chain. For simplicity we can just assume there’s only NN hopping H = \sum_i c_i^{\dagger}c_{i + 1} + h.c. .

The two-level system is coupled to the bulk, and the transition between levels induces an effective (diagonal) dipole potential

\langle \alpha \bold{s_{bulk}}| H_{int}| \alpha'\bold{s'_{bulk}}\rangle \sim \ \delta_{\bold{s} \bold{s' }}(1 - \delta_{\alpha\alpha'})\frac{pd }{r ^3},\ \text{where}\ d\sim \sum_ii\hat{n}_i|\bold{s}\rangle

that is felt by all the sites in the bulk, i.e., the induced Hamiltonian on each bulk site i has the form

H_i \sim (c_{\alpha} + c^{\dagger}_{\alpha}) \frac{1}{r_i} \sum_j j \hat{n_i}\hat{n_j}

Added to the Hamiltonian. The problem is that obviously QN is conserved in the bulk, but there is transition in the two-level system.

4

I gave it a bit more thought, and I can see that the ( c_{\alpha} + c^{\dagger}_{\alpha}) here breaks the symmetry (or that PN is not a good quantum number).
So if I use the algorithm without QN conservation, that the algorithm now explores the ‘whole’ space of potential particle numbers, would it achieve the same effect (albeit less efficient) if I define some kind of ‘chemical potential’, that is an additional term in the Hamiltonian of the form
\Lambda |\sum_i \hat{n}_i - N|
(or if I can define this operator at all in the MPS language), where \Lambda is a very large number compared to the GS and lowly excited state, and N is the target number of particles in my system? So that any states with the ‘wrong’ particle number will be high up in the spectrum and not be found by the GS algorithms. Will this incur any numerical issue?