About static spin model

ITensor Julia Questions
,
1

Dear there,
I’ve a question about using dmrg for electron-spin model. So I have a model about interacting electrons and pseudospins. The Hamiltonian is:

H=-\sum_{j,\sigma} c^\dag_{j,\sigma}c_{j+1,\sigma}+h.c+\sum_j aS^z_j(c^\dag_{j,\sigma}c_{j+1,\sigma}+h.c)+bS^x_j

Here S^z_j is a pseudospin which is same as a regular spin but just not the spin of the electrons in the system, so the picture is I have a chain of pseudospins and a chain of electrons and some coupling between chains (I’ve asked about setting this up before, and thank you so much for your help for that!)

And my question is, when I set b=0, so the pseudospin configuration is static with no dynamics. Then if I give the system an initial pseudospin configuration for example (-1/2)^j, i.e. a antiferromagnatic pattern (I do this by using a staggered initial MPS as below). I wonder will Sweeping in dmrg be able to find the true group state pattern for pseudospins that makes E_{tot} minimum or I will be stuck in initial MPS pattern because I don’t have pseudospin dynamics (S^x term)? Thank you!

for i in 1: N
    state_spin[i] = (isodd(i) ? "Dn" : "Up")
end
psi0 = productMPS(sites, state, 10)
2

It’s a good question… interesting. I would say to be really cautious here and not assume that DMRG will be able to find the right pattern of spins for you if they don’t have any dynamics (off-diagonal matrix elements in the Hamiltonian).

I think this situation calls for a bit of creativity. Personally one thing I would try is to extrapolate your results from the b > 0 limit to the b \rightarrow 0 limit. Then for each finite value of b you could hopefully get nicely converged results then reliably extrapolate properties to the limit of b=0.

Another option is if you have some kind of theoretical argument for what the minimum-energy pattern of \langle S^z_j \rangle should be, then you could just not even set the S^z_j as actual operators but just replace them with the numbers \langle S^z_j \rangle and solve that Hamiltonian instead. If you imagined doing that for some unknown \langle S^z_j \rangle values and take derivatives with respect to them and set these derivatives to zero, you might be able to use calculus arguments combined with numerics to work out what the optimal (extremal) values of \langle S^z_j \rangle are.

1 Like
4

Thanks for your suggestions! I got some time to try both ways during holiday, and find something interesting (but also confusing…). Basically, I combined two options. I first used mean-field theory approximation, assuming a minimum-energy pattern of both \langle S^z_j\rangle and bond operator \langle c^\dag_{j,\sigma}c_{j+1,\sigma}+h.c\rangle, then I got a free mean-field H​. By minimizing the energy, I got a range for both coupling constants a, b for where there is an ordered phase. By ordered I mean both |\langle S^z_j\rangle| and |\langle c^\dag_{j,\sigma}c_{j+1,\sigma}+h.c\rangle| nonzero. And you know, as in usual mean-field calculatioin, these two order parameters are related to each other, so if one is finite/ordered, another should also be ordered.

Then I solved my original H using DMRG with a pair of (a,b) in that ordered parameter range. I choose a spot in that range where b is 0.5, so that hopefully my pseudospin has enough dynamics to make dmrg converge and it’s also not too large to make my mean-field result unreliable. Then I calculated both |\langle S^z_j\rangle| and |\langle c^\dag_{j,\sigma}c_{j+1,\sigma}+h.c\rangle| with ground-state MPS, if they are indeed nonzero then numerics agree with the mean-field results.

But my problem is, I got an indeed antiferromagnetic ordered \langle S^z_j\rangle with a finite amplitude~0.1. But at same time, \langle c^\dag_{j,\sigma}c_{j+1,\sigma}+h.c\rangle looks random over sites with amplitude~10^{-5}, so enssentially zero. I feel this problem is not related to MF part, because I still use original H in dmrg, it just seems like the bond order can’t be detected… So I was wondering do you happen to know what can caused such situation. My system size is only 20, and I use bonddim=1500, 50 sweeps, I’m guessing maybe it’s not really converged? And my code for calculating bond expectation is below. Maybe this is wrong…?

Any suggestion would be really appreciated!! And I really want to thank you for providing such a creative way to treat this static spin situation!

for j in 1:N-1
      orthogonalize!(psi1, j)
      psi1_j = psi1[j]*psi1[j+1]
      psidag1_j = dag(prime(psi1[j]*psi1[j+1], "Site"))
      ampo_b= -b*(op(sites,"Cdagup",j)*op(sites,"Cup",j+1)+op(sites,"Cdagup",j+1)*op(sites,"Cup",j)+op(sites,"Cdagdn",j)*op(sites,"Cdn",j+1)+op(sites,"Cdagdn",j+1)*op(sites,"Cdn",j))
      bond1[j] = b*scalar(psidag1_j * ampo_b * psi1_j) 
end
5

Glad you are making some progress. I want to be cautious though and not say that I could really suggest any approach that I could guarantee would work here. There could be all sorts of approaches one could try, and various reasons why the results could come out in a not ideal way when you have non-dynamic degrees of freedom in your Hamiltonian.

One suggestion, though, is just to not make your spins part of your Hilbert space at all (meaning don’t even have sites for them). Instead you could just think of them as classical parameters, similar to a or b or some external fields applying to your system. Then you would just program in a certain choice of these spins when you make the OpSum to define your Hamiltonian. That way you could be sure that DMRG isn’t doing anything to change your S^z_j values becaus they would not even be part of the “quantum mechanics” that is accessible to the DMRG algorithm (i.e. they wouldn’t be part of your Hilbert space at that level of the formalism you would be using).

As far as your code to compute the bond operator, while it looks ok at a brief glance, I would say always be careful with fermions and the best thing to do would be to test it out on some example states where you know the right answer and/or have another reliable way to compute it. A nice thing to do here would be to make some random MPS and compute those bond expectation values using both your code and using the OpSum system to make some of the same operators as MPOs. Then you can use the inner function to evaluate the expected values of those MPOs with the same random MPS and confirm that everything matches. Finally you can go back to using your code, which will be much faster than using MPOs, once it is tested and you feel you can trust it.

1 Like