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.