Support biquadratic interactions for SU(N) coherent states, and beyond

[kbarros]

Currently we support scalar biquadratic interactions in :dipole mode. Eventually we should allow this also for :SUN mode. The latter would traditionally require thinking about a basis for SU(N) generators that explicitly defines spin operators $\hat S$ as well as operators $\hat Q$ quadratic in spin.

Rather than hand-coding these generators, it might be better to work entirely with polynomials of the spin operators, and their tensor products, e.g. $S_j = S \otimes I$ and $S_k = I \otimes S$ to represent spins for sites $j$ and $k$ in the context of the bond $(j,k)$. This approach would generalize to any order exchange interaction, as well as to entangled units.

[kbarros]

Here is some code to construct the spin operators as tensor products, and to use SVD to compress operators.

using Sunny, LinearAlgebra

# Produces a matrix representation of a tensor product of operators, C = A⊗B.
# Like built-in `kron` but with permutation. Returns C_{acbd} = A_{ab} B_{cd}.
function kron_operator(A::AbstractMatrix, B::AbstractMatrix)
    TS = promote_type(eltype(A), eltype(B))
    C = zeros(TS, size(A,1), size(B,1), size(A,2), size(B,2))
    for ci in CartesianIndices(C)
        a, c, b, d = Tuple(ci)
        C[ci] = A[a,b] * B[c,d]
    end
    return reshape(C, size(A,1)*size(B,1), size(A,2)*size(B,2))
end

# Use SVD to find the decomposition D = ∑ₖ Aₖ ⊗ Bₖ, where Aₖ is N₁×N₁ and Bₖ is
# N₂×N₂. Returns the list of matrices [(A₁, B₁), (A₂, B₂), ...].
function svd_tensor_expansion(D::Matrix{T}, N1, N2) where T
    @assert size(D, 1) == size(D, 2) == N1*N2
    D = permutedims(reshape(D, N1, N2, N1, N2), (1,3,2,4))
    (; S, U, V) = svd(reshape(D, N1*N1, N2*N2))
    ret = []
    for (k, σ) in enumerate(S)
        if abs(σ) > 1e-12
            u = reshape(U[:, k], N1, N1)
            v = reshape(V[:, k], N2, N2)
            push!(ret, (σ*u, conj(v)))
        end
    end
    return ret
end

# Returns the spin operators for two sites, with Hilbert space dimensions N₁ and
# N₂, respectively.
function spin_pair(N1, N2)
    S1 = Sunny.spin_matrices(N1)
    S2 = Sunny.spin_matrices(N2)
    I1 = Ref(Matrix(I, N1, N1))
    I2 = Ref(Matrix(I, N2, N2))
    return (kron_operator.(S1, I2), kron_operator.(I1, S2))
end

N1 = 3
N2 = 3

# Check Kronecker product of operators
A1 = randn(N1,N1)
A2 = randn(N1,N1)
B1 = randn(N2,N2)
B2 = randn(N2,N2)
@assert kron_operator(A1, B1) * kron_operator(A2, B2) ≈ kron_operator(A1*A2, B1*B2)

# Check SVD decomposition
S1, S2 = spin_pair(N1, N2)
B = (S1' * S2)^2                                # biquadratic interaction
D = svd_tensor_expansion(B, N1, N2)             # a sum of 9 tensor products
@assert sum(kron_operator(d...) for d in D) ≈ B # consistency check

[kbarros]

Update: User interface design has been moved to #177. This comment is out of date.

Here is a proposal for the user-facing interface. The command below will simultaneously set a spin-spin exchange matrix J1 as well as a biquadratic interaction with general matrix J2.

set_exchange!(sys, (S1, S2) -> S1'*J1*S2 + (S1'*J2*S2)^2, bond)

Note that we would now pass an anonymous function, which allows to construct a totally general interaction. The parameters S1 and S2 each contain three spin operators for their respective sites. For example, S1[3] is the z-component of spin for the first site, in a matrix representation.

For convenience, we may consider defining

set_exchange!(sys, J::Number, bond) = set_exchange!(sys, (S1, S2) -> J * S1'*S2 , bond)
set_exchange!(sys, J::AbstractMatrix, bond) = set_exchange!(sys, (S1, S2) -> S1'*J*S2 , bond)

which is a bit shorter, and hopefully avoids breaking most existing scripts.

We will need to get rid of the current set_biquadratic_exchange! and force people to use the more general set_exchange! above. In this new interface, each call to set_exchange! will override any existing exchange for that bond (bilinear and biquadratic treated on equal footing, and must be specified together).

Anisotropy. For consistency, I would propose to adopt an analogous notation for setting single-ion anisotropy,

set_anisotropy!(sys, S -> -D*S[3]^2, i) # set an easy axis
set_anisotropy_as_stevens!(sys, O -> -(D/3)*O[2,0]^2, i) # equivalent

This allows some internal simplifications, as it is no longer needed to create symbolic representations of polynomials (everything becomes just a matrix).

Interesting observation: because set_exchange! is so general, it now allows a different path for constructing anisotropies. For example, this command will set an easy axis anisotropy for the atoms at bond.i and bond.j

set_exchange!(sys, (S1, S2) -> -D*(S1[3]^2 + S2[3]^2), bond)

I would not recommend this style, but it is interesting to see that it is possible.

Inhomogeneous systems. For homogeneous systems,set_exchange! will propagate interactions by symmetry. For inhomogeneous systems, one will use instead,

sys2 = to_inhomogeneous!(sys)
set_exchange_at!(sys2, (S1, S2) -> S1'*J*S2 , site1, site2)

Looking forward. All of this is designed to work in :SUN mode. There will be many restrictions on the allowed exchange interactions in :dipole mode. Probably we will restrict to the interactions allowed by SpinW: 3x3 exchange matrices, scalar biquadratic, and arbitrary single-ion anisotropy.

[kbarros]

Initial support landed in #176.

Questions about user interface have been pushed to #177.