Local entanglement and String Order Parameter in dimerized models

TL;DR

This paper applies String Order Parameter to fermionic systems, specifically a dimerized fermionic chain, to identify topological phases and analyze phase transitions using analytical and entanglement measures.

Contribution

It introduces a novel application of SOP to fermionic models and analytically characterizes topological phase transitions via Toeplitz determinants and entanglement measures.

Findings

SOP successfully identifies topological phases in fermionic chains.

Phase transition characterized by non-analytic behavior of concurrence derivative.

Analytical expressions for SOP and entanglement measures in the model.

Abstract

In this letter, we propose an application of String Order Parameter (SOP), commonly used in quantum spin systems, to identify symmetry-protected topological phase (SPT) in fermionic systems in the example of the dimerized fermionic chain. As a generalized form of dimerized model, we consider a one-dimensional spin-1/2 XX model with alternating spin couplings. We employ Jordan-Wigner fermionization to map this model to the spinless Su-Schrieffer-Heeger fermionic model (SSH) with generalized hopping signs. We demonstrate a phase transition between a trivial insulating phase and the Haldane phase by the exact analytical evaluation of reconstructed SOPs which are represented as determinants of Toeplitz matrices with the given generating functions. To get more insight into the topological quantum phase transition (tQPT) and microscopic correlations, we study the pairwise concurrence as a local entanglement measure of the model. We show that the first derivative of the concurrence has a non-analytic behaviour in the vicinity of the tQPT, like in the second order trivial QPTs.

Equations

$$H_1=J\sum _{i=1}^N(S_{2i-1}^y{S}_{2i}^y+S_{2i-1}^x{S}_{2i}^x)$$

$$H_2=J^{\prime }\sum _{i=1}^N(S_{2i}^y{S}_{2i+1}^y+S_{2i}^x{S}_{2i+1}^x)$$

$$H=H_1+H_2$$

$$H_1=\cos (\theta )\sum _{i=1}^{\frac{N}{2}}(S_i^{+(a)}{S}_i^{-(b)}+S_i^{-(a)}{S}_i^{+(b)})$$

$$H_2=\sin (\theta )\sum _{i=1}^{\frac{N}{2}-1}(S_i^{+(b)}{S}_{i+1}^{-(a)}+S_i^{-(b)}{S}_{i+1}^{+(a)})$$

$$S_i^{+(a)}=a_i^{†}{e}^{i\pi {\sum }_{m<i}(a_m^{†}{a}_m+b_m^{†}{b}_m)}$$

$$S_i^{+(b)}=b_i^{†}{e}^{i\pi {\sum }_{k<i}(a_k^{†}{a}_k+b_k^{†}{b}_k+a_i^{†}{a}_i)}$$

$$S_i^{z(a)}=a_i^{†}{a}_i-1/2$$

$$\mathcal{H}=\sum _{i=1}^{N_c}\cos (\theta )(a_i^{†}{b}_i+b_i^{†}{a}_i)+\sin (\theta )(b_i^{†}{a}_{i+1}+a_{i+1}^{†}{b}_i)$$

$$a_j=\frac{1}{\sqrt{N_c}}\sum _{k\in BZ}{e}^{-ikj}{a}_k,$$

$$b_j=\frac{1}{\sqrt{N_c}}\sum _{k\in BZ}{e}^{-ikj}{b}_k,$$

$$\mathcal{H}$$

$${\hat{F}}_k=\left[\begin{array}{cc}0& e^{-i\frac{k}{2}}\sin (\theta )+e^{i\frac{k}{2}}\cos (\theta )\\ \sin (\theta )e^{i\frac{k}{2}}+e^{-i\frac{k}{2}}\cos (\theta )& 0\end{array}\right]$$

$$\mathcal{H}$$

$$E^{\pm }(k,\theta )=\pm \sqrt{1+\sin (2\theta )\cos (k)}$$

$${\Psi }_k={\hat{U}}_k{\Gamma }_k$$

$$\begin{gathered} v^{\pm }=\frac{1}{\sqrt{2}}\text{matrixquantity}(\pm \frac{E^{+}(k,\theta )}{e^{i\frac{k}{2}}\sin (\theta )+e^{-i\frac{k}{2}}\cos (\theta )} \\ 1) \end{gathered}$$

$${\Delta }_{G1}(k=\pm \pi )=2\sqrt{1-\sin (2\theta )}.$$

$${\Delta }_{G2}(k=0)=2\sqrt{1+\sin (2\theta )}$$

$$O^{S(\gamma )}=\underset{r\to \infty }{\lim }{O}^{\gamma }(r)$$

$$O_{l,m}^{\gamma }(r)=-4\text{expectationvalue}{S}_l^{\gamma (b)}{e}^{i\pi (S_{l+1}^{\gamma (a)}+S_{l+2}^{\gamma (b)}+\cdots +S_{m-1}^{\gamma (b)})}{S}_m^{\gamma (a)}$$

$$O_{l,m}^z=\text{expectationvalue}{e}^{i\pi (S_l^{z(b)}+S_{l+1}^{z(a)}+S_{l+2}^{z(b)}+\cdots +S_{m-1}^{z(b)}+S_m^{z(a)})}$$

$$O_{l,m}^{x,y}=\text{expectationvalue}(S_l^{+(b)}-S_l^{-(b)})(S_{l+1}^{+(a)}-S_{l+1}^{-(a)})\dots (S_m^{+(a)}-S_m^{-(a)})$$

$$O_{l,m}^z=\text{expectationvalue}{e}^{i\pi (b_l^{†}{b}_l+a_{l+1}^{†}{a}_{l+1}+\cdots +a_m^{†}{a}_m)}$$

$$O_{l,m}^{x,y}=\text{expectationvalue}(b_l^{†}+b_l)(a_{l+1}^{†}-a_{l+1})\dots (a_m^{†}-a_m)$$

$$A_m=(a_m^{†}-a_m)$$

$$B_m=(a_m^{†}+a_m)$$

$$C_m=(b_m^{†}-b_m)$$

$$D_m=(b_m^{†}+b_m)$$

$$\text{anticommutator}{A}_m{A}_n=\text{anticommutator}{C}_m{C}_n=-2{\delta }_{m,n}$$

$$\text{anticommutator}{B}_m{B}_n=\text{anticommutator}{D}_m{D}_n=2{\delta }_{m,n}$$

Full text

I Introduction

A detailed study of the phases of matter and also transition between them have long been an actual problem of condensed matter physics. While Landau-Ginzburg theory can provide such understanding and explain trivial quantum phase transitions in terms of ’spontaneous symmetry-breaking’, topological phases of matter and their transitions can not be characterized in the frame of the theory. Particularly, this is due to the lack of a local order parameter that can identify a topological phase. Consequently, exploration of these non-trivial quantum phases and their transitions has been an actively studied topic in the field of condensed matter physics.

The oldest example of a symmetry-protected topological (SPT) phase of 1D quantum spin systems is the Haldane phase, which is the ground state phase of the standard antiferromagnetic Heisenberg model with $S=1$. The gaped nature of the ground state phase and the existence of the edge states were predicted a long time ago Haldane (1983a, b); Affleck et al. (1988), while no local order was identified to characterize the phase. Later on, Den Nijs and Rommelse den Nijs and Rommelse (1989) have shown that the Haldane phase has a hidden Neel order, which can be identified by the non-local string order parameter (SOP). Kennedy and Tasaki Kennedy and Tasaki (1992) explained the origin of the gap in terms of the hidden $\mathbb{Z}_{2}\times\mathbb{Z}_{2}$ symmetry breaking using non-local unitary transformation. This transformation usually converts SOPs to a simple ferromagnetic order parameter shedding light on the symmetry structure which is protecting the Haldane phase. As a result, for a long time the only condition for existence of the Haldane phase in spin systems was a non-zero value of SOP. Employing this, in the early 1990’s the Haldane phase was detected in modified Heisenberg spin chains with $S=1$ and Heisenberg ladders with $S=\frac{1}{2}$ both numerically and analytically Hida (1991); Kennedy (1990). In particular, the Haldane phase of the bond-alternating Heisenberg model was also studied numerically using SOP via exact diagonalization methods Hida (1992a, b).

While the Haldane phase was identified via SOPs in spin systems, in fermionic and bosonic systems rigorous and generalized tools such as the Berry phase have been used. As a result, the question of generalization of fermionic order parameters for spin systems has arisen. Hastugai et al. Hirano et al. (2008) introduced a local spin twist as a generic parameter of Berry phase for gapped spin systems, which has been used extensively Maruyama et al. (2009). Along this line of thought, Rosch and Anfuso Anfuso and Rosch (2007a) considered the reverse process, i.e using SOP as an SPT order parameter for fermionic systems. They reconstructed SOP for identification of the SPT phase in band insulators and studied the robustness of the SOP against perturbation terms Anfuso and Rosch (2007b). More recently, SOP was used to identify the Haldane phase in a 1D bosonic lattice, a topological Kondo insulator and a Kitaev Ladder Sorensen et al. (2019); Xu et al. (2018); Lisandrini et al. (2017), demonstrating their solid application for non-spin systems. However, these works were performed numerically using DMRG White (1992); Schollwöck (2011) and SOPs were not evaluated analytically. As the first result of the present paper, we analytically evaluate SOP to identify the SPT phase in the example of a dimerized fermionic model. The current well-studied model is considered due to it’s simplicity and it can demonstrate clearly the physical picture of SOPs. This result, particularly, expands the list of exactly calculable order parameters of SPT phases.

Another question which we address in this letter is the behaviour of local two-site entanglement in the vicinity of topological quantum phase transitions (tQPT). There have been numerous studies carried out on the characterization of topological phases and their transitions in which entanglement entropy and entanglement spectrum are introduced to identify tQPT Jiang et al. (2012); Li and Haldane (2008); Pollmann et al. (2010) and the corresponding phases. However, only a few studies were carried on the behaviour of local pairwise entanglement in tQPTs Nehra et al. (2018). A good measure of local bipartite entanglement is concurrence, defined by Wootters Wootters (1998), and it has been used extensively to study trivial QPT occuring in 1D quantum spin chains Glaser et al. (2003); Dillenschneider (2008); Arnesen et al. (2001); Osborne and Nielsen (2002); Vidal et al. (2003); Osterloh et al. (2002). In particular, in Wu et al. (2004) it was shown, that in trivial QPT non-analytic behaviour of the concurrence or its derivative determine the order of the quantum transition. Hence, an interesting question to investigate is the behaviour of microscopic correlations and the concurrence in the vicinity of tQPT, which we discussed in the second part of the paper.

The paper is organized as follows. In Sec. II we introduce the considered spin model and describe basic steps of fermionization and exact analytical diagonalization of it. Fermionized versions of the SOPs, their exact evaluation and the final phase diagram are presented in Sec. III. We discuss the behaviour of local bipartite entanglement in the vicinity of topological QPT in Sec. IV. Finally, conclusive notes and future outlooks are presented in the last section Sec. V.

II Model definition

We consider a generalized spin-1/2 $XX$ chain of length $N$ with bond-alternation defined as:

[TABLE]

[TABLE]

where $J$,$J^{\prime}$ are spin couplings. The sum of these two terms defines the Hamiltionian of our model:

[TABLE]

We note that $J,J^{\prime}$ can be ferromagnetic (FM) and antiferromagnetic (AFM) . Since the relative absolute value of couplings define ground state phases of (3), we parameterize $J=2\cos(\theta)$ and $J^{\prime}=2\sin(\theta)$, where $\theta\in[0,2\pi)$.

Due to broken translational symmetry, one can relabel spins within the effective unit cell and rewrite Eqs.(1-2) in the following form:

[TABLE]

[TABLE]

where (a) and (b) refer to spins located within the $i$-th unit cell. For $\theta=\frac{n\pi}{2}$ , $n\in\mathbb{Z}$ one has localized spin dimers along the chain . Hereafter, we assume periodic boundary conditions (PBC).

To fermionize (4-5) we use the Jordan - Wigner transformation Jordan and Wigner (1928) which we define as:

[TABLE]

[TABLE]

[TABLE]

where operators $a_{i},b_{i}$ are spinless fermionic operators with standard anticommutation rules. This transformation maps (3) to the following fermionic model:

[TABLE]

where $N_{c}$ is the number of effective unit cells. In derivation of (9) we neglected boundary terms, assuming that the chain is long enough, so that the boundary terms’ contribution to the physical quantities $\sim O(\frac{1}{N})$ is negligible.

The Hamiltonian (9) is a generalized non-interacting Su-Schrieffer-Heeger (SSH) model. It should be noted, that in contrast to the well-studied standard SSH model Su et al. (1980); Asbóth et al. (2016), here one may have different signs of the fermionic hopping parameters (due to the freedom of the signs for spin couplings). Furthermore, if the isotropic Dzyaloshinskiy-Moriya interaction along the chain $H_{DM}=D^{z}\sum_{i=1}^{L}(S^{x}_{i}S^{y}_{i+1}-S^{y}_{i}S^{x}_{i+1})$ is introduced to the model (3), one may have complex values of the hoping couplings $J,J^{\prime}$. This representation freedom of spin chains in terms of fermions enriches the physical phenomena and may correspond to fictitious fermionic models.

To express (9) in $k$-space, we use the discrete Fourier transformation of the form:

[TABLE]

Then, the Hamiltonian (9) in the momentum space can be written as:

[TABLE]

where

[TABLE]

and $\Gamma_{k}^{\dagger}=(a_{k}^{\dagger},b_{k}^{\dagger})$.

One may get spectrum of the Bogolyubov quasiparticles by diagonalization of the kernel (12),

[TABLE]

with $\Psi^{\dagger}_{k}=(\alpha_{k}^{\dagger},\beta^{\dagger}_{k})$. The matrix elements of the new diagonal kernel $\Lambda_{k}$ define the spectrum of the quasiparticles:

[TABLE]

The vector $\Psi_{k}$ of qusiparticle operators is connected to the vector $\Gamma_{k}$ via,

[TABLE]

where $\hat{U}_{k}$ is the transformation matrix formed out by the eigenvectors,

[TABLE]

In Fig.1 we plot the spectrum of the quasiparticles for different values of intracell and intercell couplings. For $\theta=0$ the intercell coupling vanishes, thus one has localised fermions within the unit cell. This flattens the quisiparticle bands, making them to have infinite mass. This character in fact is preserved for any $\theta=\frac{n\pi}{2}$ , $n\in\mathbb{Z}$.

For other values of the couplings, the bands get dispersed with a finite density of states (DOS). Due to the freedom on couplings’ signs, there are two different modes for the system: couplings with the same signs (FM-FM and AFM-AFM, red sectors in Fig.1 (b)) and couplings with alternating signs (AFM-FM and FM-AFM, blue sectors in Fig.1(b)).

The first mode (red sectors of the unitary circle) has a particle - hole symmetric excitation spectrum with the minimum gap at $k=\pm\pi$, as shown in Fig.1(a) (for $\theta=\frac{\pi}{6}$). The value of the gap at $k=\pm\pi$ is

[TABLE]

These red sectors correspond to the SSH limit, since both couplings have the same signs. In the second mode, which corresponds to the blue sectors of the unitary circle, the intercell (intracell) coupling is ferromagnetic, while intracell (intercell) is antiferromagnetic. Then minimum value of the gap is reached at $k=0$ and can be determined as,

[TABLE]

From the last expressions (18-19) for the gap value, on can see that there are gap closures which occure for $\theta=\frac{(2n+1)\pi}{4}$, $n\in\mathbb{Z}$. These gap closures hint for a quantum phase transition, which happens to have topological nature.

III Fermionized String Order parameters

Previously, the tQPT of the model in the Heisenberg limit was studied via DMRG using SOP and symmetry fractionalization technique Wang et al. (2013); Haghshenas et al. (2014). In this section, we show an exact calculation of SOP employing the fermionization approach for our $XX$ model.

For the generalized model (3) the SOP can be defined as Hida (1991),

[TABLE]

[TABLE]

where $r=\absolutevalue{m-l}$ and $\gamma\in[x,y,z]$. For convenience, we further call (21) Long-Range String (LRS). From Eq.(20) it is clear that a long-distant limit of LRS defines SOP. Here, the prefactor $-4$ is used to normalize, so that in the dimerized topological limit $\theta=\frac{(2n+1)\pi}{2}$ it yields $O^{\gamma}_{l,m}=1$.

For longitudinal LRS, using $e^{i\pi S^{z}}=\frac{S^{z}}{2i}$, Eq.(20) can be written as:

[TABLE]

Since the model is isotropic, $O^{x}_{l,m}=O^{y}_{l,m}$, one needs to evaluate only one of them. Using $e^{i\pi S^{y}_{i}}=S^{+}_{i}-S^{-}_{i}$, LRS for the y-component can be expressed as:

[TABLE]

By the use of JWT, we convert LRS for spins to fermionic LRS :

[TABLE]

which is obtained using the fact that the number of operators in the exponent is always even. For the transverse LRS we obtain,

[TABLE]

We introduce the following operators for a more convenient notation:

[TABLE]

[TABLE]

[TABLE]

[TABLE]

It is straightforward to derive anticommutation relations of the pairs for introduced operators:

[TABLE]

[TABLE]

where $\delta_{m,n}$ is Kronecker delta function. While for any other pair $N$ and $M$ from the set above, we have:

[TABLE]

Then, fermionized LRSs (24-25) can be expressed in terms of the introduced $A-B-C-D$ operators,

[TABLE]

[TABLE]

Using Wick theorem, including only non-zero contractions, we get the following result:

[TABLE]

[TABLE]

where by site $l=1$ we mean any chosen $b$ site as reference and $m$ any chosen site $a$ which follows $b$. Matrix elements $K_{l,m}$ and $Q_{l,m}$ are defined as,

[TABLE]

[TABLE]

Matrices $K$ and $Q$ are $(m-1)\crossproduct(m-1)$ Toeplitz matrices, with their corresponding generating functions.

Finally, $O^{S(x,y,z)}$ can be evaluated as the determinant of the Toeplitz matrices in the long range limit,

[TABLE]

[TABLE]

To evaluate SOP $O^{S(\gamma)}$ analytically, one may use Szegö theorem Szegö (1915) in the Haldane phase region of the phase diagram. However, on the other insulating limit this theorem is not applicable, due to the emerging zeros in the corresponding generating functions.

It can be easily seen that SOPs are unity when the wavefunction is the tensor product of maximally entangled pairs between the unit cells ($b-a$ pair), i.e

[TABLE]

This corresponds to the Haldane phase edge limit, i.e $\cot(\theta)=0$. In this limit, free edge spins can be observed in a chain with open boundary conditions as shown in Fig.2, similar to the edge states of the $S=1$ Heisenberg model. However, when the value of the intracell hopping is increased, maximally entangled spin dimers become less entangled, which results in a decreased value of $O^{S(\gamma)}$. Indeed, using Szegö theorem one can show that the value of $O^{S(z)}$ for infinitesimal intracell hopping,

[TABLE]

This result shows that for non-zero intracell hopping, $O^{S(z)}$ decreases from unity. The same behaviour can be found for $O^{S(x,y)}$. Thus, $O^{S(\gamma)}$ can be connected to the entanglement measure of $b-a$ spin pairs .

In Fig.3 we plot the longitudinal LRS $O^{z}(r)$ for $\theta\in[\frac{\pi}{6},\frac{\pi}{4},\frac{\pi}{3}]$ and corresponding SOPs $O^{S(x,z)}$ within the unit circle using Eqs.(35-40). Fig.3(a) shows that $O^{z}(r)$ exponentially converges to the finite value $O^{S(z)}\approx 0.8$ in the Haldane phase at $\theta=\frac{\pi}{3}$, while exponential convergence to $O^{S(z)}=0$ can be observed in the trivial phase, i.e at $\theta=\frac{\pi}{6}$. In the vicinity of quantum critical points (QCP), $O^{S(z)}(r)$ decreases accroding to a power law, which leads to the divergence of correlation length $\zeta$ at the QCP. Thus, one needs to evaluate $O^{z}(r)$ in the very large distance at these points to extract the value of $O^{S(z)}$. A similar behaviour can be observed on the corresponding phases for $O^{x,y}(r)$.

In Fig.3(b) we show the phase diagram of the model based on the $O^{S}$. It demonstrates that for the values of $\absolutevalue{\cot(\theta)}<1$, the system is in the Haldane phase, which is characterized by the finite value of SOPs. The maximum values of $O^{S}$ at $\cot(\theta)=0$. In the QCPs $O^{S}$ vanish and for the remained circle sectors, i.e $\absolutevalue{\tan(\theta)}<1$ the system is in the trivial phase. In analogy to (21), one may consider the following relation for LRS,

[TABLE]

In contrast to the introduced SOPs, new SOPs derived from Eq.(43), detect trivial phase of the model, i.e they are unity when the maximally entangled pairs are $a-b$ spin pairs.

It should be noted that SOPs are not robust against any symmetry-breaking terms. As an example, an applied magnetic field in $z$ direction which breaks $SU(2)$ symmetry, makes $O^{x,y}_{l,m}$ converge to a finite value even in the trivial insulating phase. An alternative example is the interacting limit, i.e bond-alternating $XXZ$ chain, for which SOP is a valid order parameter only for the $SU(2)$ symmetric case, i.e in the $XXX$ limit.

IV Entanglement and Correlations

More information about the tQPT discussed in the previous section, can be extracted from a measure of local entanglement. To gain an insight about the entanglement structure of the model, we calculate the two-site density matrix and use two-site concurrence $C^{a,b}_{i,j}$ as an entanglement measure Wootters (1998). In trivial QPTs, non-analytic divergent character of the concurrence usually corresponds to first order quantum phase transitions (1QPT), while this behaviour in the derivative of the concurrence is peculiar to the second order transitions (2QPT). Thus, it is interesting to refer to which class of the trivial QPTs the tQPT corresponds. Here, we will exactly evaluate $C^{a,b}_{i,j}$ in terms of two-point correlation functions and investigate the behaviour of it and $\frac{dC^{a,b}_{i,j}}{d\theta}$ in QCPs.

The functional dependence of the two-site concurrence on standard correlation functions in generalized quantum spin chains was derived earlier Glaser et al. (2003) from the reduced density matrix for any two given sites of the chain. Here, we simplify this derivation using symmetries of the Hamiltonian (3).

The reduced density matrix for any given two sites $i$ and $j$ can be expressed in terms of two-point correlation functions in the following manner (in $\ket{\uparrow\uparrow}$, $\ket{\uparrow\downarrow}$,$\ket{\downarrow\uparrow}$,$\ket{\downarrow\downarrow}$ basis):

[TABLE]

where $u^{(a,b)}_{i,j}=\frac{1}{4}+G^{zz(a,b)}_{i,j}$, $w^{(a,b)}_{i,j}=\frac{1}{4}-G^{zz(a,b)}_{i,j}$ and $t^{(a,b)}_{i,j}=2G^{xx(a,b)}_{i,j}$. Derivation of two-point correlation functions $G^{\gamma\gamma(a,b)}_{i,j}$ for all spin components can be found in Appendix.

To derive (44) we exploited the spin-flip symmetry (which leads to $\expectationvalue{S^{z}}=0$) and also the symmetry of Hamiltonian (3) with respect to the rotation in the $x-y$ plane. Due to broken inversion symmetry, one has two inequivalent $\rho^{a,b}_{i,j}$.

The concurrence of two sites $C_{i,j}$ may be computed from the density matrix $\rho_{i,j}$ as $C_{i,j}=[\lambda_{1}-\lambda_{2}-\lambda_{3}-\lambda_{4}]$, where the $\lambda_{i}$ are the eigenvalues in decreasing order of the matrix $R=\sqrt{\sqrt{\rho}\tilde{\rho}\sqrt{\rho}}$. The matrix $\tilde{\rho}_{i,j}$ can be expressed as $\tilde{\rho}_{i,j}=(\sigma^{y}\otimes\sigma^{y})\rho^{*}_{i,j}(\sigma^{y}\otimes\sigma^{y})$. Then, the two site concurrence for any given sites $i$ and $j$ $C^{a,b}_{i,j}$ can be expressed in terms of two-point correlation functions as:

[TABLE]

In Fig.4 we present concurrence $C^{a}_{i,i+1}$ and $C^{b}_{i,i+1}$ for nearest-neighbours calculated from Eq.(45). For vanishing intracell couplings (i.e at $\theta=\frac{(2n+1)\pi}{2},n\in\mathbb{Z}$), $C^{b}_{i,i+1}=1$, showing maximal entanglement of the $b-a$ spin pairs. This result is consistent to our SOP result corresponding to maximally entangled $b-a$ spin pairs. Non-zero intracell spin couplings, perturb this entanglement structure of pairs, inducing a decrease in $C^{b}_{i,i+1}$. This corresponds to the breaking of entanglement monogamy of $b-a$ spin pairs. Further decrease of $C^{b}_{i,i+1}$ continues when the intracell coupling is incremented. However, unlike the SOPs, it doesn’t vanish at QCPs. In these gapless Luttinger liquid points every site is entangled with the nearest neighbour with the same concurrence $C^{a}_{i,i+1}=C^{b}_{i,i+1}\approx 0.35$ (marked with a dashed line). Further increment of $J$ results in the sudden death of entanglement which happens at $\theta=\frac{(2n+1)\pi}{2}+\frac{\pi}{3},n\in\mathbb{Z}$.

In the opposite way $C^{a}_{i,i+1}$ behaves, which measures an entanglement of $a-b$ spin pairs. Obviously, it is unity for vanishing intercell hopping values $J$.

As pointed out, the derivative of the concurrence should diverge for QCPs, at least this is what one expects in trivial QPTs. Indeed, from Fig.4(b) it is clear that for $\theta=\frac{(2n+1)\pi}{4},n\in\mathbb{Z}$, one has sharp peaks in $\frac{dC^{a,b}}{d\theta}$ signalling about the tQPT. This character is inherent for the trivial 2QPTs. Similarly, one can show that the second derivative of the ground state energy $\frac{d^{2}E}{d\theta^{2}}$ (exact form of $E$ is presented in Appendix) is also singular at QCPs, which is also a peculiarity of the 2QPTs. Except for these peaks, step like drops of the function can be noticed. These does not effect the phase diagram of the model, however, it is related to the sudden death (or birth) of the entanglement.

The concurrences for next neighbours $C^{a,b}_{i,i+r}=0$ for all values of $\theta$. To show this and to clarify the origin of the singular peaks, in Fig.5 we plot the absolute value of the two-point correlation functions $|G^{xx}(r,\theta)|=|\expectationvalue{S^{x(a)}_{i}S^{x(a,b)}_{i+r}}|$ and $|G^{xx}(r,\theta)|=|\expectationvalue{S^{x(a)}_{i}S^{x(a,b)}_{i+r}}|$ for $a-x$ pair of spins. To begin with, the derivative of nearest-neighbour correlators have singular peaks at QCP (not shown), resulting in similar peaks of the concurrence. However, derivatives of these correlators for other distances vanish at QCPs, since these points are stationary points. This can be seen from Fig.5(a-b) where the absolute value of correlators for other distances have a maximum value at QCPs. These values decrease with distance, which is known to have a power law behaviour. One may notice that for $r>1$, the value of $2\absolutevalue{G^{xx}}-\absolutevalue{1/4+G^{zz}}<0$. Thus, from Eq.(45) one concludes that $C^{a,b}_{i,i+r}$ vanish for $(r>1)$, i.e starting from the next nearest neighbours.

In contrast to derived SOPs, two-site concurrences are universal: they identify diverse type of quantum phase transitions even in the presence of symmetry-breaking terms. In this respect, SOPs are limited on application, one should modify them in a necessary way to determine the phase transition. For this, of course, one needs to have a priliminary knowledge of the corresponding phases.

V Summary and conclusions

In this paper, we discussed the implementation algorithm of SOP as an exactly calculable order parameter for the SPT phases and studied the local pairwise entanglement behaviour at tQPT. We mapped the bond-alternating spin chain via Jordan - Wigner fermionization to the generalized version of the SSH model and demonstrated the correspondence of the topological insulating phase of the standard SSH model and the Haldane phase of the bond-alternating spin model. We derived exact analytical expressions for two-point transverse and longitudinal correlation functions. By the fermionization of the longitudinal and transverse SOPs we showed the exact expression of them in terms of the determinant of the Toeplitz matrices with the given generating functions. We demonstrated that the derived expressions of SOPs take their maximal value of unity in the fully dimerized SPT Haldane phase and vanish in the trivial insulating phase. To get more insight about the phase transition, we studied pairwise local entanglement (namely, two-site concurrence) in the vicinity of tQPTs. We found that the two-site concurrence vanishes for all other distances except the nearest neighbour, like in the standard Heisenberg model. The derivative of the concurrence has a singularity at the tQPT, like in trivial second order quantum phase transitions. In contrast to the SOPs, we found the two-site concurrence universal and robust against symmetry-breaking terms.

For future outlook, we note that our result of fermionized SOPs can be extended for other 1D spin chains and ladder models. However, the more important topic of fragility of SOPs needs to be investigated in a more detailed way, considering symmetries of the models.

Acknowledgements.

M.S.B thanks Prof.Aikaterini Mandilara for useful discussions.

Appendix A Calculation of zero-temperature correlation functions

To evaluate the ground states energy per spin, one needs to integrate the band below Fermi energy, which leads to:

[TABLE]

where $E$ is the elliptic integral of the second kind. To evaluate static two-point correlator as the function of spin couplings, we rewrite creation and annihilation operators in terms of quasiparticle operators:

[TABLE]

[TABLE]

From the definition of quasiparticle vacuum one can get following relationships:

[TABLE]

[TABLE]

where

[TABLE]

One can show that,

[TABLE]

Using anticommutation relations one can get a remained matrix elements. Obviously, for even $r=n-m$ correlation functions $\expectationvalue{a^{\dagger}_{m}a_{n}}=\expectationvalue{b^{\dagger}_{m}b_{n}}$ vanish :

[TABLE]

Here, indexes do not correspond to the effective unit cell, but to a site, thus we will work in a reduced Brillouin zone.

Next, we evaluate non-vanishing correlation functions $W_{n-m}=\expectationvalue{a^{\dagger}_{m}b_{n}}$ and $Z_{n-m}=\expectationvalue{b^{\dagger}_{m}a_{n}}$; Further, we assume $n>m$.

[TABLE]

Another non-vanishing correlator is:

[TABLE]

It is can be noticed that $W_{r}=Z_{-r}$. Integrals above for nearest neighbours (i.e $\absolutevalue{r}=1$) are the linear composition of $K(\sin(\theta),\cos(\theta))$ and $E(\sin(\theta),\cos(\theta))$ which are elliptic integrals of the first and the second kinds. For remaining distances, the analytical solution is involved and the integration consists other terms.

Transverse correlation functions

We consider correlation function of x - components of spins, i.e $G^{xx(a)}_{m,n}=\expectationvalue{S_{m}^{(a)x}S_{n}^{(b)x}}$

[TABLE]

Index $a$ in $G^{xx(a)}_{m,n}$ shows that we look x-component correlation of any $a-b$ spin pair. In fermionic representation it has the following form:

[TABLE]

Thus, we have for every site one has exponential string operator:

[TABLE]

where $p=a,b$ is fermionic operator. Equation (A.12) can be written in terms of $A-B-C-D$ operators (26-29) as:

[TABLE]

One can evaluate this correlation function using Wick theorem, since operators anticommute. Thus, we firstly evaluate expectation values of the form $\expectationvalue{N_{m}M_{n}}$ where $N$ and $M$ belong to a set of operators defined in (26-29).

Using (53) it can be shown that for different sites $\expectationvalue{N_{m}N_{n}}=0$ for any $N$. Furthermore, $\expectationvalue{A_{m}B_{n}}=\expectationvalue{C_{m}D_{n}}=0$.

The only non-zero elements are:

[TABLE]

[TABLE]

At this point, this information is enough to evaluate expression (56) using Wick theorem which leads to the following results:

[TABLE]

where matrix elements of $G^{xx}$ are defined by $W_{n-m}$ and $Z_{n-m}$. While for the standard XY model Lieb et al. (1961) McCoy and Wu (1967) the matrix $G^{xx}$ is a well studied Toeplitz matrix, in our case the matrix is a block Toeplitz matrix, where diagonal elements oscillate between two values.

We note that using (62) $G^{xx(a)}$ can be evaluated for any $a-x$ pair, where $x\in[a,b]$. Obviously, $G^{xx(b)}_{m,n}\neq G^{xx(a)}_{m,n}$. It is straightforward task to show that an equivalent matrix for $G^{xx(b)}_{m,n}$ can be obtained by changing $W_{n-m}\rightarrow Z_{n-m}$.

Finally, we note that due to the isotropy of the model in $x-y$ plane, the the correlation function for $y$ - component behaves in the same way.

Longitudinal correlation function

Calculation of $G^{zz}_{m,n}$ is not as complex job as transverse correlation functions. We firstly consider $\expectationvalue{S^{(a)z}_{m}S^{b(z)}_{n}}$ case, longitudinal correlation of $a-b$ spin pairs:

[TABLE]

Fermionized expression is,

[TABLE]

In terms of $A-B-C-D$ operators,

[TABLE]

The only non-zero contraction is:

[TABLE]

Thus,

[TABLE]

Similarly, an expression for a longitudinal correlation function of $b-a$ spin pairs is,

[TABLE]

For even distances ( $a-a$ spin pair or $b-b$ spin pairs), $C^{zz(a,b)}_{m,n}$ vanishes, due to the fact that all contractions of the same type of fermions (i.e consisting only $a$ ’s or $b$’s) vanish in any permutations.