Phase transition in SSH model

ITensor C++ Questions
1

Hi everyone,
I am trying to observe a phase transition in the non interacting SSH model, described by the following Hamiltonian( for reference see https://arxiv.org/pdf/1810.13286.pdf):

H=\sum_{i\in A ,j\in B} -[\hat{S^+_i}\hat{S^-_j}+\hat{S^+_j}\hat{S^-_i}]-J[\hat{S^+_j}\hat{S^-_i}+\hat{S^+_i}\hat{S^-_j}]

The SSH model is formulated on a one-dimensional lattice, divided into two sub-lattices: A = {1, 3, . . . , N − 1}, involving odd lattice sites, and B = {2, 4, . . . , N}, with an even number of sites N and staggered hopping of particles between A and B ( here I considered one hopping strength equal to 1 and the other equal to J). Depending on the value of J <1 or >1 , I should observe a phase transition, from a trivial state with to a topological state, characterized by the presence of edge state (see the following image for context).

Capture
Capture768×760 132 KB

I calculated the density operator for each site ( since I mapped the model to a spin half chain I used the following density operator)

\hat{n}=\frac{1}{2}Id - \hat{S_z}

and printed the result and I always get 0,5 for every site, no matter the value of J and I don’t see why the code it is not working as it should. I attach the full code:

#include "itensor/all.h" 
#include "myclass.h"
#include <complex>  
#include <cmath>
#include<iostream>
#include<stdlib.h>
#include<time.h>
#include <vector>

using namespace itensor;
using std::vector;


int main()
    {
    int N = 100;
    auto sites = SpinHalf(N,{"ConserveQNs=",false}); 
    auto state = InitState(sites);
    auto psi0 = randomMPS(state); 
    auto ampo = AutoMPO(sites); 
    std ::vector<double> J{0.1, 0.2, 0.3, 0.5, 0.7, 0.8, 1, 2, 3, 4, 5, 10};
    for(int l=0; l < J.size(); l++)
    {
        for(int j = 1; j <= N-1; ++j)
        {
            if (j%2 == 1)
            {
                ampo += -1,"S+",j,"S-",j+1;
                ampo += -1,"S+",j+1,"S-",j;
            }
            else (j%2 ==0);
            {
                ampo += -J[l],"S+",j+1,"S-",j;
                ampo += -J[l],"S+",j,"S-",j+1;
            }
        }
        auto H = toMPO(ampo);
        auto sweeps = Sweeps(5);
        sweeps.maxdim() = 10,40,50,60,100;
        sweeps.cutoff() = 1E-10;
        auto [energy,psi] = dmrg(H,psi0,sweeps,{"Quiet",true});
        auto psidag = dag(prime(psi,"Site"));

        for(int i = 1; i <= N; ++i)
        {
            auto k = sites(i);
            auto n = ITensor(k,prime(k)); 
            n.set(k(2),prime(k)(2), 1);
            psi.position(i);
            auto phi = psi(i); 
            auto bra_phi = dag(prime(phi,"Site"));
            auto exp_n = elt(bra_phi * n * phi);
            print(" ",exp_n);
        }
    }

    return 0;
    }

I always get the following:

Capture 2
Capture 21580×126 31.3 KB

Thanks a lot for your time and patience to go through all of this !
Francesco

2

Hi Francesco,
Thanks for the question. So as best I can tell, thinking of the model as a spin model which only involves S^+ and S^- operators, there is no reason for the spins to have an expected value along the z direction. That is, I would expect that \langle S^z_j \rangle=0 for all sites j. And then that would give

\langle n_j \rangle = \langle \frac{1}{2} I_j \rangle - \langle S^z_j \rangle = \frac{1}{2} \langle I_j \rangle = \frac{1}{2}

which is what you observe.

Is there a reason to expect that \langle S^z_j \rangle \neq 0 for this system?