Exponential of single-site operators with QN

ITensor Julia Questions
1

Hi, I have a user-defined site-type as follows:

    # Reference: https://docs.itensor.org/ITensorMPS/stable/examples/Physics.html
    # User-defined spin type with Z2xZ2 charge
    function ITensors.space(::SiteType"SpinZ2xZ2=1"; conserve_qns=false)
        if conserve_qns
            return [
                QN(("X",1,2),("Z",1,2))=>1, # |+>
                QN(("X",1,2),("Z",0,2))=>1, # |0>
                QN(("X",0,2),("Z",1,2))=>1, # |->
            ]
        end
        return 3
    end

    # User-defined site name 
    ITensors.state(::StateName"+", ::SiteType"SpinZ2xZ2=1") = [1 0 0];
    ITensors.state(::StateName"0", ::SiteType"SpinZ2xZ2=1") = [0 1 0];
    ITensors.state(::StateName"-", ::SiteType"SpinZ2xZ2=1") = [0 0 1];

    # User-defined charge operators 
    E = [1 0 0; 0 1 0; 0 0 1];

    Sx = [0 1 0; 1 0 0; 0 0 0];
    Sy = [0 0 0; 0 0 im; 0 -im 0];
    Sz = [0 0 1; 0 0 0; 1 0 0];

    # Assign them to ITensors operators
    ITensors.op(::OpName"Id", ::SiteType"SpinZ2xZ2=1") = E
    ITensors.op(::OpName"Sx", ::SiteType"SpinZ2xZ2=1") = Sx
    ITensors.op(::OpName"Sy", ::SiteType"SpinZ2xZ2=1") = Sy
    ITensors.op(::OpName"Sz", ::SiteType"SpinZ2xZ2=1") = Sz

which I believe is implementing the user-defined site type correctly, and I was able to run a few correct DMRG simulations with these sites. However, when I try to compute operators like \exp(i\pi S^\alpha), via something like

ii = siteind("SpinZ2xZ2=1";conserve_qns=true)
U = op(ii,"Sz")
U = exp(im*pi*U)

I ran into the following error:

ERROR: exp currently supports only block-diagonal matrices

The version with conserve_qns=false works correctly, so I believe this is a not a bug of my own. Is there a way to fix this?

I need operators of this type because I wanted to compute string order parameters across phase transition points, and being able to utilize the underlying Z2xZ2 symmetries will be hugely beneficial.

2

Ok, I figure out a solution that works for the specific situation

Since the operators Sx, Sy, Sz have definite flux, the exponential of these operators must also carry definite flux. Therefore, a simple work around would be to define the exponential of these operators directly,

    # Exponential of operators 
    ExpIπSx = [-1 0 0; 0 -1 0; 0 0 1];
    ExpIπSy = [1 0 0; 0 -1 0; 0 0 -1];
    ExpIπSz = [-1 0 0; 0 1 0; 0 0 -1];
    ITensors.op(::OpName"ExpIπSx",::SiteType"SpinZ2xZ2=1") = ExpIπSx
    ITensors.op(::OpName"ExpIπSy",::SiteType"SpinZ2xZ2=1") = ExpIπSy
    ITensors.op(::OpName"ExpIπSz",::SiteType"SpinZ2xZ2=1") = ExpIπSz

However, I would love to know if there are better ways to do exp() for general operators within the current framework.

3

I wouldn’t expect there to be a better way in general. This is because \exp(i \theta S_x) only has definite flux for particular values of \theta. In fact, numerically, it only has definite flux for \theta = 0.

using LinearAlgebra

Sx = [0 1 0; 1 0 0; 0 0 0]
print("e^(i π Sx) = ")
display(exp(im * pi * Sx))

print("e^(i π/4 Sx) = ")
display(exp(im * pi / 4 * Sx))

Output:

e^(i π Sx) = 3×3 Matrix{ComplexF64}:
 -1.0+0.0im           0.0+1.52695e-16im  0.0+0.0im
  0.0+1.52695e-16im  -1.0+0.0im          0.0+0.0im
  0.0+0.0im           0.0+0.0im          1.0+0.0im
e^(i π/4 Sx) = 3×3 Matrix{ComplexF64}:
 0.707107-0.0im            0.0+0.707107im  0.0+0.0im
      0.0+0.707107im  0.707107-0.0im       0.0-0.0im
      0.0+0.0im            0.0+0.0im       1.0+0.0im

For something slightly more general, although also kinda janky you could do the following

function ITensors.op(::OpName"ExpIΘSx",::SiteType"SpinZ2xZ2=1"; θ::Real)::Matrix{ComplexF64}
  return ComplexF64[abs(aij) < 1e-10 ? 0 : aij for aij in exp(im * θ * Sx)]
end

@show flux(op(ii, "ExpIΘSx"; θ=π))