# Phase transition in SSH model

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).

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:

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

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?