Eigenvalues of spin S-1 chain

[HamidArianZad]

Dear EDKit developer(s),

I am writing to seek assistance regarding the computation of eigenvalues for a spin-1 diamond-like chain featuring two distinct exchange couplings in the presence of an external magnetic field along the z-axis, as depicted in the figure below.

I tried following code:

using EDKit
using LinearAlgebra

L = 12
for Bz=0:0.1:8
       mat = spin((1, "xx"), (1, "yy"), (1, "zz"), (1, "z"), D=3)
       mats = fill(mat, L)
       inds = [
                     [[i, mod(i+1, L)] for i=1:3:L];    # J_1 exchange coupling
                     [[i, mod(i+1, L)] for i=1:3:L];    # J_1 exchange coupling
                     [[i, mod(i+1, L)] for i=1:3:L];    # J_1 exchange coupling
                     [[i, mod(i+2, L)] for i=1:3:L];    # J_1 exchange coupling
                     [[i, mod(i+2, L)] for i=1:3:L];    # J_1 exchange coupling
                     [[i, mod(i+2, L)] for i=1:3:L];    # J_1 exchange coupling
                     [[i, mod(i+1, L)] for i=2:3:L];    # J_2 exchange coupling
                     [[i, mod(i+1, L)] for i=2:3:L];    # J_2 exchange coupling
                     [[i, mod(i+1, L)] for i=2:3:L];    # J_2 exchange coupling
                     [[i, mod(i+2, L)] for i=2:3:L];    # J_1 exchange coupling
                     [[i, mod(i+2, L)] for i=2:3:L];    # J_1 exchange coupling
                     [[i, mod(i+2, L)] for i=2:3:L];    # J_1 exchange coupling
                     [[i, mod(i+1, L)] for i=3:3:L];    # J_1 exchange coupling
                     [[i, mod(i+1, L)] for i=3:3:L];    # J_1 exchange coupling
                     [[i, mod(i+1, L)] for i=3:3:L];    # J_1 exchange coupling
                     [i for i=1:L]                        # Bz magnetic field
                   ]
       B = TranslationalBasis(L=L, N=L, base=3)
       H = operator(mats, inds, B)
      vals, vecs = Array(H) |> Hermitian |> eigen
      println(vals)

end

In my calculations, I have considered the conservation ofSzT (N=L). However, I have encountered several challenges in proceeding with the calculations. Despite consulting the examples provided in your documentation, I have been unable to resolve these issues. I would greatly appreciate it if you could address the following concerns:

1. How can I incorporate an onsite magnetic field along the z-direction? I attempted methods such as (1, "z") and [spin("Z") for i=1:L], but encountered errors.
2. What approach should I take to introduce distinct exchange interactions J_1 and J_2 into the model?
3. I am interested in obtaining the eigenenergies of the Hamiltonian for this model at various values of the magnetic field within the interval 0<=B_z<=8 (specifically, for B_z=0:0.1:8). Could you please provide guidance on modifying the code to achieve this?

Your assistance in resolving these issues would be greatly appreciated.

[jayren3996]

Hi Hamid,

Thank you for your message. I have created a Jupyter Notebook demonstrating the steps, and you can find it in the example folder. I hope you find it helpful.

I would like to share a few comments:

• Utilizing the trans_inv_operator is recommended for translational invariant Hamiltonians, as it simplifies the process. Additionally, leveraging translational symmetry can further reduce the dimension of the many-body Hamiltonian. In your example, specifying the momentum sector looks crucial, otherwise, the dimension of the matrix seems too large.
• Given that the translation involves a 3-site shift, the current approach involves merging the spins in a unit cell into a local quantum dot with increased dimensionality. It requires defining local terms in this enlarged system. In the future, I am considering an upgrade to the package for a more transparent procedure.

If you encounter any issues or notice mistakes in the provided example, please don't hesitate to bring them to my attention. As the package is not extensively documented and is primarily used within my group, I am pleased to assist others who find it useful.

Best regards, Jie

[HamidArianZad]

Hi Jie,

Thank you very much for your efficient code. I made some modifications and executed the results. Upon comparison, I noticed a discrepancy between the exact diagonalization (ED) results obtained from QuSpin and those from EDKit for a small system with L=2 (total number of spins N=6). In QuSpin, I did not consider any symmetry, thus obtaining the full ED results. Please see below figure for reference:

Blue broken line is EDKit results and black solid like is QuSpin results.

Below is the EDKit code I used:

using EDKit, LinearAlgebra
using DelimitedFiles

# Hamiltonian parameters:
J1 = 1.0
J2 = J1
J3 = 2.7
J4 = J1
J5 = J1
L = 2
k = 0
N = 3L
sz_min = 0
sz_max = N
B_min = 0
B_max = 8
uB = 0.1
# -------------------------------------------------------
all_energies = []
for (m_B, Bz) in enumerate(B_min:uB:B_max)
    # current_energies = []
    for i in sz_min:sz_max
    # ---------------------------------
        ha = Bz
        hb = ha
        hc = ha

        # onsite Hamiltonian matrix
        O1 = spin((J1, "xx1"), (J1, "yy1"), (J1, "zz1"), D=3)
        O2 = spin((J2, "x1x"), (J2, "y1y"), (J2, "z1z"), D=3)
        O3 = spin((J3, "1xx"), (J3, "1yy"), (J3, "1zz"), D=3)

        hz = spin((ha, "z11"), (hb, "1z1"), (hc, "11z"), D=3)

        h1 = O1 + O2 + O3 + hz
        C1 = spin((J4, "1x1x11"), (J4, "1y1y11"), (J4, "1y1y11"), D=3)
        C2 = spin((J5, "11xx11"), (J5, "11yy11"), (J5, "11yy11"), D=3)

        h2 = C1 + C2
        Nj = diag( spin((1,"z11"), (1,"1z1"), (1,"11z"), D=3) + 3*I )

        # since Julia vector start from 1, hence the `sj+1`
        total_z(s) = sum(Nj[sj+1] for sj in s)
        basis = TranslationalBasis(L=L, f = s -> total_z(s) == i, k=k, base=3^3)
        H = trans_inv_operator(h1, 1:1, basis) + trans_inv_operator(h2, 1:2, basis)
        energy = eigvals(Hermitian(H))
        for sublist in energy
            for element in sublist
                push!(all_energies, element)
            end
        end
    end 
    # push!(all_energies, current_energies)
    # ----------------------------------------------------------------------
    println("*************************************************************")
    println("-------------- Finished job: B = ", Bz, "-------------------")
    println("*************************************************************")
end

open("Energies_DiamondS1_cSzT_EDKit_J$(J3)_L$(L)_Bmin$(B_min)_Bmax$(B_max).txt", "w") do io
    for energies in all_energies
        for En in energies
            writedlm(io, En, '\n')
        end
    end
end

In above code I tried to keep fix total_z(s) in each step of magnetic field such that all possible fixed values of magnetization are taken into account.

I really like your package, especially for its concise code structure. I'm optimistic that it will assist me in achieving my objectives effectively.

Your insights and comments on this matter would be greatly appreciated.

[jayren3996]

Hi Hamid,

One thing I noticed in your code: if you are interested in obtaining the whole spectrum, you should do a for loop over k = 0,1,...,L; also, for the spin-1 system, the maximal magnon filling is 2N (sz_max = 2N).

If it still does not solve the discrepancy, could you provide the complete code that generates the results, both in EDKit and QuSpin so that I can check it myself?

[HamidArianZad]

Dear Jie,

I found the problem. Please edit the code that you have added to examples. Namely, the part

C1 = spin((J4, "1x1x11"), (J4, "1y1y11"), (J4, "1y1y11"), D=3)
C2 = spin((J5, "11xx11"), (J5, "11yy11"), (J5, "11yy11"), D=3)

should be

C1 = spin((J4, "1x1x11"), (J4, "1y1y11"), (J4, "1z1z11"), D=3)
C2 = spin((J5, "11xx11"), (J5, "11yy11"), (J5, "11zz11"), D=3)

Thanks! Now the code works well. For the case k=0, I can obtain the results that I want.

[jayren3996]

Hi Hamid,

Thanks for finding out the mistake I made. I am happy that it works for you.

[HamidArianZad]

Dear Jie,

I still see mistake in the line: C1 = spin((J4, "1x1x11"), (J4, "1y1y11"), (J4, "1y1y11"), D=3). It should be like: C1 = spin((J4, "1x1x11"), (J4, "1y1y11"), (J4, "1z1z11"), D=3). Namely, (J4, "1z1z11") is correct for correlation along z-axis.

[jayren3996]

Hi Hamid,

Thanks, I think I fixed it this time.