SU(4) antiferromagnetic Heisenberg modle on honeycomb lattice

Magnon bands and inelastic neutron spectra by SU(4) linear spin wave theory

The following codes could compute the flavor wave dispersions of the SU(4) antiferromagnetic Heisenberg model, i.e. $H=J∑\\_{⟨ij⟩}[∑\\_{1<=a<b<=5} Γᵢᵃᵇ Γⱼᵃᵇ - ∑ₐΓᵢᵃΓⱼᵃ]$ :

using SUNSpinWaveTheory
using Plots: plot, plot!
using QuantumLattices: Lattice, Point, Hilbert, Algorithm, ReciprocalPath, atol,ReciprocalZone
using TightBindingApproximation: EnergyBands
using QuantumLattices: Fock, @σ_str, Onsite, expand, bonds, MatrixCoupling, FID, reciprocals, @hexagon_str
using StaticArrays: @SVector
using QuantumLattices: @indexes, Index, Coupling, Coulomb

lattice = Lattice(
    [0.0, 0.0], [1/(sqrt(3)),0.0];
    vectors = [[sqrt(3)/2, -0.50], [0.0, 1.0]],
    )
cell = lattice

#construct the Hamiltonian
σx = [0 1; 1 0]
σy = [0 -im; im 0]
σz = [1 0; 0 -1]
σ0 = [1 0; 0  1]
Γ¹ = kron(σz, σy); Γ² = kron(σz, σx); Γ³ = kron(σy, σ0);Γ⁴ = kron(σx, σ0);Γ⁵ = kron(σz, σz);

#H=J∑\_{⟨ij⟩}[∑\_{1<=a<b<=5}Γᵢᵃᵇ Γⱼᵃᵇ-∑ₐΓᵢᵃΓⱼᵃ]
Γ(x,y) = (x*y-y*x)/2im
Γ¹²=(Γ¹*Γ²-Γ²*Γ¹)/2im; Γ¹³=Γ(Γ¹,Γ³);Γ¹⁴=Γ(Γ¹,Γ⁴);Γ¹⁵=Γ(Γ¹,Γ⁵)
Γ²³=Γ(Γ²,Γ³);Γ²⁴=Γ(Γ²,Γ⁴); Γ²⁵=Γ(Γ²,Γ⁵)
Γ³⁴=Γ(Γ³,Γ⁴);Γ³⁵=Γ(Γ³,Γ⁵); Γ⁴⁵=Γ(Γ⁴,Γ⁵)
Γab = [Γ¹²,Γ¹³,Γ¹⁴,Γ¹⁵,Γ²³,Γ²⁴,Γ²⁵,Γ³⁴,Γ³⁵,Γ⁴⁵];
Γa = [Γ¹, Γ² , Γ³ ,Γ⁴ ,Γ⁵ ]

j, h = 1.0, 0.001
Jmat = zeros(ComplexF64,16,16)

for m in Γab
    Jmat[:,:] += kron(m, m)
end
for m in Γa
    Jmat[:,:] -= kron(m, m)
end
hmat = Γ¹²
f1(bond) = iseven(bond[1].site) ? 1 : -1

J = SUNTerm(:J, Jmat, 4, 4, 1; value=j)
J1 = SUNTerm(:B, hmat, 4, 4, 0; value=h, amplitude=f1)

hilbert = Hilbert(pid=>Fock{:b}(4, 1) for pid in 1:length(cell))

magneticstructure = MagneticStructure(
    cell,
    Dict(1 => [1.5707960439947803, 1.5707963267641, 1.570796326798081, 3.1432634733226568, 3.1430008310022868, 4.2168412327339055],
         2 => [2.8285525202602516e-7, 1.5708000153983595, 1.5707999995656055, 3.1429367971485704, 3.1430659136725807, 1.075295154664448],
    )
)
eng = SUNLSWT(lattice, hilbert, (J, J1), magneticstructure)

#start optimize:
op = optimorder(eng; numrand = 1);

#the classical energy
println(expand(op[1].Ω).contents[()].value)

antiferromagnet = Algorithm(:SquareAFM, op[1] )

# order parameters
px = Dict(pid => Γ¹² for pid in 1:length(cell))
py = Dict(pid => Γ³⁴ for pid in 1:length(cell))
pz = Dict(pid => Γ⁵  for pid in 1:length(cell))
orderpara = [localorder(antiferromagnet.frontend, px),localorder(antiferromagnet.frontend, py),localorder(antiferromagnet.frontend, pz)]
println(orderpara)
-9.002 - 6.660796046664699e-16im
Dict{Int64, ComplexF64}[Dict(2 => -0.999999999999999 + 0.0im, 1 => 0.9999999999999988 + 0.0im), Dict(2 => -0.999999999999999 + 0.0im, 1 => 0.9999999999999988 + 0.0im), Dict(2 => 1.0000000000000002 + 0.0im, 1 => 0.9999999999999999 + 0.0im)]

Band structure

sx = [0       sqrt(3)/2      0      0;
     sqrt(3)/2   0            1.0    0;
       0        1.0           0     sqrt(3)/2;
       0         0           sqrt(3)/2  0
]
mx = MatrixCoupling(:, FID, sx, :, :)
sy = [0       -im*sqrt(3)/2      0      0;
     im*sqrt(3)/2   0            -1im    0;
       0        1im           0     -im*sqrt(3)/2;
       0         0           im*sqrt(3)/2  0
]
my = MatrixCoupling(:, FID, sy, :, :)
sz = [3/2   0       0      0;
       0    1/2     0      0;
       0    0       -1/2   0;
       0    0       0     -3/2
]
mz = MatrixCoupling(:, FID, sz, :, :)

path = ReciprocalPath(reciprocals(lattice), hexagon"Γ-K-M-Γ,60°", length=100)
sx = expand(Onsite(:mu, 1.0+0.0im, mx,amplitude=nothing, modulate=true),bonds(cell,0),hilbert,half=false)
sy = expand(Onsite(:mu, 1.0+0.0im, my,amplitude=nothing, modulate=true),bonds(cell,0),hilbert,half=false)
sz = expand(Onsite(:mu, 1.0+0.0im, mz,amplitude=nothing, modulate=true),bonds(cell,0),hilbert,half=false)
ss = @SVector [sx,sy,sz]

ins = antiferromagnet(:INS,
    Spectra{InelasticNeutron}(path, range(0.0, 20, length=501), (ss,ss); fwhm=0.2, scale=identity, gauss=false, atol=1e-9)
    )
energybands = antiferromagnet(:EB, EnergyBands(path; atol=1e-9))

plt = plot()
plot!(plt, ins)
plot!(plt, energybands, color=:white, linestyle=:dash, linealpha=0.1)
#display(plt)
Example block output

Spectra of multipole operators

#multipole
quadrupolar = [multipoleOperator(3//2,3//2,2,m) for m=-2:2]
Mcouplings = [MatrixCoupling(:, FID, i, :, :) for i in quadrupolar]
M₋₂ = expand(Onsite(:mu1, 1.0+0.0im, Mcouplings[1],amplitude=nothing, modulate=true), bonds(cell,0), hilbert, half=false)
M₋₁ = expand(Onsite(:mu2, 1.0+0.0im, Mcouplings[2],amplitude=nothing, modulate=true), bonds(cell,0), hilbert, half=false)
M₀ = expand(Onsite(:mu3,  1.0+0.0im, Mcouplings[3],amplitude=nothing, modulate=true), bonds(cell,0), hilbert, half=false)
M₁ = expand(Onsite(:mu4,  1.0+0.0im, Mcouplings[4],amplitude=nothing, modulate=true), bonds(cell,0), hilbert, half=false)
M₂ = expand(Onsite(:mu5,  1.0+0.0im, Mcouplings[5],amplitude=nothing, modulate=true), bonds(cell,0), hilbert, half=false)
Mss = [M₋₂, M₋₁, M₀, M₁, M₂]
Mssd = [M₋₂', M₋₁', M₀', M₁', M₂']

insmultipole = antiferromagnet(:Multipole,
    Spectra{Multipole}(path, range(0.0, 20, length=501), (Mss, Mssd); fwhm=0.2, gauss=true, scale=x->x, atol=1e-9)
    )
plt1 = plot()
plot!(plt1, insmultipole, clims=(0, 100))
plot!(plt1, energybands, color=:white, linestyle=:dash, linealpha=0.1)
#display(plt1)
Example block output

Spectra of multipole operators with E=1.0

#Ecut
nx, ny = 32, 32
b₁, b₂ = reciprocals(lattice)
zone = ReciprocalZone([b₁, b₂-b₁/2],
    0=>1, 0=>1; length=(nx, ny),
    ends=((true, true), (true, true))
    )
inszone = antiferromagnet(:MultipoleZone,
    Spectra{Multipole}(zone, range(0.0, 20, length=501), (Mss, Mssd); fwhm=0.2, gauss=true, scale=identity, atol=1e-9)
    )
ecut, dE = 3.0, 0.4
# data = spectraEcut(inszone[2], ecut, dE, nx, ny)
plt2 = plot()
plot!(plt2, inszone, ecut, dE)
#display(plt2)
Example block output