Local entanglement and string order parameter in dimerized models

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.


Introduction
A detailed study of the phases of matter and also the trans ition 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 symmetrybreaking', topological phases of matter and their trans itions 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, explora tion of these nontrivial 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 antifer romagnetic 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 [1][2][3], while no local order was identified to characterize the phase. Later on, Den Nijs and Rommelse [4] have shown that the Haldane phase has a hidden Neel order, which can be identified by the nonlocal string order parameter (SOP). Kennedy and Tasaki [5] explained the origin of the gap in terms of the hidden Z 2 × Z 2 symmetry breaking using nonlocal unitary transformation. This trans formation 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 nonzero value of SOP. Employing this, in the early 1990s the Haldane phase was detected in modified Heisenberg spin chains with S = 1 and Heisenberg ladders with S = 1 2 both numerically and analytically [6,7]. In par ticular, the Haldane phase of the bondalternating Heisenberg model was also studied numerically using SOP via exact diag onalization methods [8,9].
While the Haldane phase was identified via SOPs in spin systems, in fermionic and bosonic systems rigorous and gen eralized 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 [10] introduced a local spin twist as a generic parameter of Berry phase for gapped spin systems, which has been used exten sively [11]. Along this line of thought, Rosch and Anfuso [12] 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 [13]. More recently, SOP was used to identify the Haldane phase in a 1D bosonic lattice, a topological Kondo insulator and a Kitaev Ladder [14][15][16], demonstrating their solid application for nonspin systems. However, these works were performed numerically using DMRG [17,18] and SOPs were not evalu ated 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 twosite 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 [19][20][21] and the corresponding phases. However, only a few studies were carried on the behaviour of local pairwise entanglement in tQPTs [22]. A good measure of local bipartite entanglement is concurrence, defined by Wootters [23], and it has been used extensively to study trivial QPT occuring in 1D quantum spin chains [24][25][26][27][28][29]. In particular, in [30] it was shown, that in trivial QPT nonanalytic behaviour of the con currence 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 section 2, we intro duce 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 section 3. We discuss the behaviour of local bipartite entanglement in the vicinity of topological QPT in section 4. Finally, conclusive notes and future outlooks are presented in the last section section 5.

Model definition
We consider a generalized spin1/2 XX chain of length N with bondalternation defined as: where J, J are spin couplings. The sum of these two terms defines the Hamiltionian of our model: We note that J, J can be ferromagnetic (FM) and anti ferromagnetic (AFM) . Since the relative absolute value of couplings define ground state phases of (3), we parameterize J = 2 cos(θ) and J = 2 sin(θ), where θ ∈ [0, 2π).
Due to broken translational symmetry, one can relabel spins within the effective unit cell and rewrite equations (1) and (2) in the following form: where (a) and (b) refer to spins located within the ith unit cell. For θ = nπ 2 , n ∈ Z one has localized spin dimers along the chain . Hereafter, we assume periodic boundary conditions (PBC).
To fermionize (4) and (5), we use the Jordan-Wigner trans formation [31] which we define as: where operators a i , b i are spinless fermionic operators with standard anticommutation rules. This transformation maps (3) to the following fermionic model: 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 ∼ O( 1 N ) is negligible.
The Hamiltonian (9) is a generalized noninteracting Su-Schrieffer-Heeger (SSH) model. It should be noted, that in contrast to the wellstudied standard SSH model [32,33], here one may have different signs of the fermionic hopping param eters (due to the freedom of the signs for spin couplings). Furthermore, if the isotropic Dzyaloshinskiy-Moriya inter action along the chain is introduced to the model (3), one may have complex values of the hoping couplings J, J . This representation freedom of spin chains in terms of fermions enriches the physical phe nomena and may correspond to fictitious fermionic models.
To express (9) in kspace, we use the discrete Fourier trans formation of the form: Then, the Hamiltonian (9) in the momentum space can be written as: wherê One may get spectrum of the Bogolyubov quasiparticles by diagonalization of the kernel (12), The matrix elements of the new diagonal kernel Λ k define the spectrum of the quasiparticles: The vector Ψ k of qusiparticle operators is connected to the vector Γ k via, where Û k is the transformation matrix formed out by the eigenvectors, In figure 1, we plot the spectrum of the quasiparticles for dif ferent values of intracell and intercell couplings. For θ = 0 the intercell coupling vanishes, thus one has localised fermions within the unit cell. This flattens the quisiparticle bands, making them have infinite mass. This character in fact is pre served for any θ = nπ 2 , n ∈ 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 (FMFM and AFMAFM, red sectors in figure 1(b)) and couplings with alternating signs (AFMFM and FMAFM, blue sectors in figure 1(b)).
The first mode (red sectors of the unitary circle) has a par ticle-hole symmetric excitation spectrum with the minimum gap at k = ±π, as shown in figure 1(a) (for θ = π 6 ). The value of the gap at k = ±π is These red sectors correspond to the SSH limit, since both cou plings have the same signs. In the second mode, which cor responds to the blue sectors of the unitary circle, the intercell (intracell) coupling is ferromagnetic, while intracell (inter cell) is antiferromagnetic. Then minimum value of the gap is reached at k = 0 and can be determined as, From the last expressions (18) and (19) for the gap value, on can see that there are gap closures which occure for θ = (2n+1)π 4 , n ∈ Z. These gap closures hint for a quantum phase transition, which happens to have topological nature.

Fermionized string order parameters
Previously, the tQPT of the model in the Heisenberg limit was studied via DMRG using SOP and symmetry fractionalization technique [34,35]. In this section, we show an exact calcul ation of SOP employing the fermionization approach for our XX model. For the generalized model (3) the SOP can be defined as [8,9], where r = |m − l| and γ ∈ [x, y, z]. For convenience, we fur ther call (21) long-range string (LRS). From equation (20) it is clear that a longdistant limit of LRS defines SOP. Here, the prefactor −4 is used to normalize, so that in the dimerized (20) can be written as: Since the model is isotropic, O x l,m = O y l,m , one needs to eval uate only one of them. Using e iπS y i = S + i − S − i , LRS for the y component can be expressed as: ) .
(23) By the use of JWT, we convert LRS for spins to fermionic LRS : which is obtained using the fact that the number of operators in the exponent is always even. For the transverse LRS we obtain, We introduce the following operators for a more convenient notation: It is straightforward to derive anticommutation relations of the pairs for introduced operators: where δ m,n is Kronecker delta function. While for any other pair N and M from the set above, we have: Then, fermionized LRSs (24) and (25) can be expressed in Using Wick theorem, including only nonzero contrac tions, we get the following result: 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, 1 + sin(2θ) cos(k) dk. (38) Matrices K and Q are (m − 1) × (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, To evaluate SOP O S(γ) analytically, one may use Szegö the orem [36] 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 wave function is the tensor product of maximally entangled pairs between the unit cells (b − a pair), i.e.
This corresponds to the Haldane phase edge limit, i.e. cot(θ) = 0. In this limit, free edge spins can be observed in a chain with open boundary conditions as shown in figure 2, similar to the edge states of the S = 1 Heisenberg model. However, when the value of the intracell hopping is increased, This result shows that for nonzero intracell hopping, O S(z) decreases from unity. The same behaviour can be found for O S(x,y ) . Thus, O S(γ) can be connected to the entanglement mea sure of b − a spin pairs . In figure 3 we plot the longitudinal LRS O z (r) for θ ∈ [ π 6 , π 4 , π 3 ] and corresponding SOPs O S(x,z) within the unit circle using equations (35)-(40). Figure 3(a) shows that O z (r) exponentially converges to the finite value O S(z) ≈ 0.8 in the Haldane phase at θ = π 3 , while exponential conv ergence to O S(z) = 0 can be observed in the trivial phase, i.e. at θ = π 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 ζ 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). We found the corre sponding power law decrement of SOPs at the QCP to be O S(z) (r) ≈ r −0.3 . This is in a good agreement with bosoniza tion result presented in [9], while the result there is presented for a dimerized isotropic Heisenberg model. Furthermore, we note that the obtained results are consistent with a superior accuracy with the result of DMRG calculations presented in [39].
In figure 3(b) we show the phase diagram of the model based on the O S . It demonstrates that for the values of |cot(θ)| < 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(θ) = 0. In the QCPs O S vanish and for the remained circle sectors, i.e. |tan(θ)| < 1 the system is in the trivial phase. In analogy to (21), one may consider the fol lowing relation for LRS, In contrast to the introduced SOPs, new SOPs derived from equation (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 symmetrybreaking terms. As an example, an applied magn etic 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. bondalternating XXZ chain, for which SOP is a valid order parameter only for the SU(2) symmetric case, i.e. in the XXX limit.

Entanglement and correlations
More information about the tQPT discussed in the previous section, can be extracted from a measure of local entangle ment. To gain an insight about the entanglement structure of the model, we calculate the twosite density matrix and use twosite concurrence C a,b i,j as an entanglement measure [23]. In trivial QPTs, nonanalytic divergent character of the con currence 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 twopoint correlation functions and investigate the behaviour of it and dC a,b i,j dθ in QCPs. The functional dependence of the twosite concurrence on standard correlation functions in generalized quantum spin chains was derived earlier [24] 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 twopoint correlation functions in the following manner (in | ↑↑ , | ↑↓ ,| ↓↑ ,| ↓↓ basis): for all spin components can be found in appendix.  To derive (44) we exploited the spinflip symmetry (which leads to 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 ρ a,b i,j . The concurrence of two sites C i,j may be computed from the density matrix ρ i,j as where the λ i are the eigenvalues in decreasing order of the matrix R = √ ρρ √ ρ. The matrix ρ i,j can be expressed as Then, the two site concurrence for any given sites i and j C a,b i,j can be expressed in terms of twopoint correlation functions as: In figure 4 we present concurrence C a i,i+1 and C b i,i+1 for nearestneighbours calculated from equation (45). For vanishing intracell couplings (i.e. at θ = (2n+1)π 2 , n ∈ Z), 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. Nonzero intra cell 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 cou pling is incremented. However, unlike the SOPs, it does not 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 ≈ 0.35 (marked with a dashed line). Further increment of J results in the sudden death of entanglement which happens at θ ≈ (2n+1)π 2 + 50π 171 , n ∈ 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 figure 4

(b) it is clear that for
, n ∈ Z, one has sharp peaks in dC a,b dθ signal ling about the tQPT. This character is inherent for the trivial 2QPTs. Similarly, one can show that the second derivative of the ground state energy d 2 E dθ 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 dia gram of the model, however, it is related to the sudden death (or birth) of the entanglement.
In the vicinity of the QCP, θ → θ c , W 1 (θ) is continious 1+sin(2θ) )−2 π . Exploiting an asymptotic expansion in the leading order, K(1 − Similar to the critical behaviour of the XY model [29], the derivative of the concurrence has a logarithmic divergence, indicating the Ising universality class of the transition. The concurrences for next neighbours C a,b i,i+r = 0 for all values of θ. To show this and to clarify the origin of the sin gular peaks, in figure 5 we plot the absolute value of the two point correlation | for a − r pair of spins. To begin with, the derivative of nearestneighbour correlators have sin gular 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 figures 5(a) and (b) where the absolute value of correlators for other distances have a max imum 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|G xx | − |1/4 + G zz | < 0. Thus, from equation (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 the derived SOPs, the twosite concurrences are universal: they identify diverse types of quantum phase transitions even in the presence of symmetrybreaking terms. In this respect, SOPs are limited on application, one should modify them in a necessary way to determine the phase trans ition. For this, of course, one needs to have a priliminary knowledge of the corresponding phases. On the other hand, the twosite concurrence similar to the other quantum information tools to study QPTs, does not clarify the nature of transitions. Since it captures trivial and topological QPTs equivalently in the same manner, there is no way to distinguish phases of the transition. This again opens the question of classification of phases, which can be particularly performed by nonlocal order parameters such as SOPs.

Summary and conclusions
In this paper, we discussed the implementation algorithm of string order parameter (SOP) as an exactly calculable order parameter for the symmetryprotected topological phase (SPT) phases and studied the local pairwise entanglement behaviour at topological quantum phase transition (tQPT). We mapped the bondalternating spin chain via Jordan-Wigner fermionization to the generalized version of the SSH model and demonstrated the correspondence of the topological insu lating phase of the standard SSH model and the Haldane phase of the bondalternating spin model. We derived exact ana lytical expressions for twopoint transverse and longitudinal correlation functions. By the fermioniz ation of the longitu dinal 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, twosite concurrence) in the vicinity of tQPTs. We found that the twosite 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 robust against symmetrybreaking terms and its application as an additional order parameter to be practical.
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, consid ering symmetries of the models.

Acknowledgments
MSB thanks Prof Aikaterini Mandilara for useful discussions and acknowledges ORAU grant "Dissecting the collective dynamics of arrays of superconducting circuits and quantum metamaterials" (no. SST2017031) for the support of his short stay in Nazarbayev University.

Appendix. Calculation of zero-temperature correlation functions
To evaluate the ground state energy per spin, one needs to integrate the band below the Fermi energy, which leads to: where E(x) is the Elliptic integral of the Second kind. To eval uate static twopoint correlator as the function of spin cou plings, we rewrite creation and annihilation operators in terms of quasiparticle operators: 3) From the definition of quasiparticle vacuum one can get following relationships: One can show that, (A.7) Using anticommutation relations one can get a remained matrix elements. Obviously, for even r = n − m correlation functions a † m a n = b † m b n vanish: 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 nonvanishing correlation functions W n−m = a † m b n and Z n−m = b † m a n ; Further, we assume n > m.
(A.10) It is can be noticed that W r = Z −r. Integrals above for nearest neighbours (i.e. |r| = 1) are the linear composition of K(sin(θ), cos(θ)) and E(sin(θ), cos(θ)) which are elliptic integrals of the first and the second kinds. For remaining dis tances, the analytical solution is involved and the integration consists other terms. m,n shows that we look xcomponent correlation of any a − b spin pair. In fermionic representation it has the following form: (A.12) Thus, we have for every site one has exponential string operator: (A.14) One can evaluate this correlation function using Wick the orem, since operators anticommute. Thus, we firstly evaluate expectation values of the form N m M n where N and M belong to a set of operators defined in (26)- (29).
Using (A.8) it can be shown that for different sites N m N n = 0 for any N. Furthermore, A m B n = C m D n = 0.
The only nonzero elements are: 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.

A.2. Longitudinal correlation function
Calculation of G zz m,n is not as complex job as transverse cor relation functions. We firstly consider S Finally, we note that analytically obtained resuts of phase diagram, SOPs and ground state properties are consistent with numerically performed DMRG studies, which are presented in [39].