{'fraction_non_alphanumeric': 0.044123982495521845, 'fraction_numerical': 0.016280071650454617, 'mean_word_length': 4.14392348961978, 'pattern_counts': {'":': 0, '<': 21, '<?xml version=': 0, '>': 19, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 10, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Effective mapping of spin-1 chains onto integrable fermionic models. A study of string and Neel correlation functions

We derive the dominant contribution to the large-distance decay of correlation functions for a spin chain model that exhibits both Haldane and Neel phases in its ground state phase diagram. The analytic results are obtained by means of an approximate mapping between a spin-1 anisotropic Hamiltonian onto a fermionic model of noninteracting Bogolioubov quasiparticles related in turn to the XY spin-1/2 chain in a transverse field. This approach allows us to express the spin-1 string operators in terms of fermionic operators so that the dominant contribution to the string correlators at large distances can be computed using the technique of Toeplitz determinants. As expected, we find long-range string order both in the longitudinal and in the transverse channel in the Haldane phase, while in the Neel phase only the longitudinal order survives. In this way, the long-range string order can be explicitly related to the components of the magnetization of the XY model. Moreover, apart from the critical line, where the decay is algebraic, we find that in the gapped phases the decay is governed by an exponential tail multiplied by algebraic factors. As regards the usual two points correlation functions, we show that the longitudinal one behaves in a 'dual' fashion with respect to the transverse string correlator, namely both the asymptotic values and the decay laws exchange when the transition line is crossed. For the transverse spin-spin correlator, we find a finite characteristic length which is an unexpected feature at the critical point. We also comment briefly the entanglement features of the original system versus those of the effective model. The goodness of the approximation and the analytical predictions are checked versus density-matrix renormalization group calculations.


Introduction
The Haldane phase [1], found in many low-dimensional spin systems, has attracted a great amount of attention in the last two decades both from the theoretical and from the experimental points of view. Its genuine quantum nature is signalled by two characteristic features. First, the excitation spectrum above the ground state (GS) displays a finite energy gap and, second, one can identify suitable long-ranged string correlation functions that measure a hidden topological order of the phase. The most intuitive idea to understand the physical features of the Haldane phase is probably the spin liquid picture [2]: In a spin-1 chain with Heisenberg interactions and quantization axis directed along z, let us assign the presence of an effective spin-1/2 particle with spin pointing up (down) if at the i-th lattice site S z i = +1 (−1) and no particles if S z i = 0. The Haldane phase is then interpreted as a liquid in which these effective particles carry no positional order along the chain but still retain antiferromagnetic (AFM) order in their effective spins. The positional disorder is associated with the absence of long-range order in the usual spin-1 correlation functions C α (R) ≡ (−1) R S α i S α i+R , α = x, y, z whereas the spin-1/2 magnetic order that we would get if all the sites with S z i = 0 were taken off from the chain is measured by the asymptotic value of the string correlators [3]: for R → ∞. Interestingly enough, the Haldane gap has been interpreted as the excitation energy associated with a "spinon" (or kink) with respect to the hidden order [4].
The nonvanishing values of the string-order parameters (SOP) can be understood as a spontaneous breaking of hidden (nonlocal) Z 2 symmetries of the λ − D Hamiltonian, as discussed thoroughly by Kennedy and Tasaki [1]. From a numerical inspection of the string correlation functions (1) computed on the first excited state with S z tot = 1, rather than on the GS, Elstner and Mikeska [4] argued that this excited wave function is characterized by a transition region with vanishing string correlations that connects two asymptotic limits with symmetry breaking and different values of the hidden order. In a field-theoretic approach to spin-1/2 Heisenberg chain [5] such a kink is described as an effective particle -a soliton -moving with relativistic dispersion relation. When the system is moved away from criticality, due to the action of a relevant field, the soliton acquires a nonvanishing mass or an energy gap, in the condensed matter language. In ref. [6] some of the authors have proposed a picture of the states that form the Haldane triplet at the isotropic point in terms of massive solitons and their bound states arising in the sine-Gordon formulation, valid in the neighbourhood of the critical line that marks the limit of the Haldane phase towards the large-D one (see below). The first solid numerical evidence of a nonzero Haldane gap has been provided by White and Huse [7] using the by now celebrated densitymatrix renormalization group (DMRG) method.
Actually, the Haldane phase is not restriced to spin-1 systems and can be found, for example, in spin-S Heisenberg chains for every integer value of S. According to ref. [8] the gap vanishes as the classical limit S → ∞ is approached as ∆ ∝ S −1 exp(−πS) while the behaviour of the string order is more subtle: in order to have a nonzero value one has to generalize the string correlation function of equation (1) using not π in the exponential but S-dependent optimal angles θ n = (2n + 1)π/S with n = 0, 1, . . . , S − 1. Again, when S → ∞ the resulting values of O α (θ n ) tend to zero.
It is interesting to examine also how the features of the Haldane phase are destroyed by varying the parameters of the Hamiltonian out of the isotropic spin-S Heisenberg model (S integer). In this paper we shall stick from now on to the case S = 1 and consider two types of anisotropies along z: Ising-like interactions (parametrized by λ) and single-ion terms (parametrized by D) The phase diagram of this model has been investigated in various papers with different approaches [9,3,1]. In order to fix the ideas we will refer to our recent determination [6], reported (in a simplified form) in figure 1. Fixing a nonnegative value for λ and varying D we encounter three gapped phases: the Large-D one in which O α = 0 ∀α indicating the absence of magnetic order in the effective spin-1/2 particles. Their positional degrees of freedom are also disordered. In the Haldane phase the spatial disorder persists but magnetic order emerges. As a consequence both longitudinal and transverse string order parameters (SOP) become nonzero: O α = 0. As pointed out in ref. [10] on the basis of an exact solution for an integrable variant of (2) with λ = 0 and in-plane anisotropy, the excitation gaps in the Large-D and in the Haldane phases have a rather different nature. Despite the fact that they are both found within the sector S z tot = 1, the former corresponds to a flip of a single spin out of the xy plane while the latter is related to the breaking of a two-site singlet composing the resonant-valence-bond GS similar to the one of the spin-1 chain exactly solved by Affleck, Kennedy, Lieb and Tasaki [11]. Finally, by decreasing further the value of D, we pass in the Néel phase where both positional and magnetic degrees of freedom orders are signalled by a nonvanishing (spontaneous) magnetization along z as order parameters to classify the three types of behaviour. The universality classes associated with the two transition lines will be frequently denoted using the language of conformal field theory (CFT -see, for instance, [12,9,6]), in particular by specifying the central charge c. We interpret the fully-disordered large-D phase with (O z = 0, O x = 0) and (O z = 0, M z = 0) as a spin gas. By crossing the c = 1 line we enter the Haldane phase where the effective spin-1/2 experience a first magnetic ordering: . Then, loosely speaking, at the c = 1/2 line the spin liquid crystallizes and the fully-ordered Néel phase can be interpreted as a spin solid with (O z = 0, O x = 0) and (O z = 0, M z = 0). Note the interchanged role of O x and M z (see below). In the Néel and Haldane phases G H (∞) refers to the positional order of nonzero spins [3], so that it vanishes in the Haldane phase but G H (∞) = 0 in the Néel one.
In order to determine the SOP numerically one has to extrapolate to the thermodynamic limit and to infinite distance the data computed on necessarily finite samples. However, apart from the qualitative statements made in ref. [3] about the exponential decay of the string correlation functions (as reported in table 1), the available literature contains scarce information about the spatial behaviour of such correlators and the extrapolation may become problematic, especially close to the transition lines where the bulk correlation length becomes very large. In a particular case, namely the transition from the Large-D to the Haldane phase, the low-energy physics is described by a compactified free boson field theory (c = 1 CFT). Once the compactification radius is known in some other independent way, one can read off the decay exponent of the string correlation functions from the set of scaling dimensions of the possible vertex operators. Interestingly, it turns out [13,14] that, even if the lattice model has periodic boundary conditions (PBC), the vertex operators to be associated with string correlators belong to the sector with antiperiodic boundary conditions.
The main purpose of this paper, instead, is to address the spatial behaviour of spinspin and string correlation functions in the Haldane and Néel phases, making use of a solvable theory of spinless fermions. Starting from well inside the Néel phase, where the density of sites with S z i = 0 is negligible, we approximate the problem by as- suming that the hidden magnetic order is frozen so that a given contribution to the GS wavefunction can be described by occupation numbers: no fermions if S z i = 0 and one fermion when |S z i | = 1, no matter the orientation, which is dictated by the underlying string order. The details of this approach will be presented in Sect. 2; actually it is very close to what done by Gómez-Santos in ref. [15]. The difference here is that we include also the single-ion anisotropy term and, in fact, the two formulations are related by a particle-hole transformation. The novelty is that we work out in detail the mapping of the spin-spin and string correlation functions (subsec. 2.1) onto fermionic correlators, so that we can derive in Sect. 3 the precise form of their asymptotic behaviour at large distances by exploiting the machinery of Toeplitz determinants. Sect. 4 reports a comparison with DMRG simulations of the system in equation (2) while Sect. 5 is devoted to conclusions.

Mapping onto spinless fermions
The basic idea underlying the approximation used in this work is the spin solid picture of the Néel state(s): which is, in fact, the GS of the Hamiltonian (2) for λ → ∞ at fixed D or D → −∞ and λ > 0. Actually the GS is doubly degenerate: for a given configuration of the type (3) with, say, S z i = 1 at the reference site i = 0, the energy is unchanged by the Z 2 transformation T = exp(iπ j S y j ) that performs a π-rotation about the y-axis (spinflip). We will refer to |N and |N = T |N as Néel and anti-Néel states, respectively. Now, in a perturbative fashion, when |D| and/or λ ≫ 1 the effect of the transverse terms in the Hamiltonian S x,y i S x,y i+1 is to: i) create pairs of adjacent sites with S z i = S z i+1 = 0 : | ↑↓ → |0 0 ; ii) move the zeroes in the AFM background, e.g.: | ↑ 0 → |0 ↑ ; iii) re-create a pair ↑↓ or ↓↑ in place of a pair of adjacent zeroes.
Notice that ii) preserves the AFM order, albeit not on nearest neighbours but mediated by string of zeroes (hidden order). Again, due to the AFM order (induced by λ > 0 and by the transverse terms), even if both states of iii) can be created in an "island" of zeroes, as far as the low-energy part of the spectrum is concerned, one of the two will be preferred according to the orientation of the surrounding spins, that is, by the hidden AFM order. Note also that |N and |N are connected through a large number of virtual processes, so that in the thermodynamic limit only one of the two will be selected by a spontaneous symmetry breaking mechanism induced by an infinitesimal staggered magnetic field. Alternatively, if the system under consideration is described by a thermal density matrix exp(−βH), when β → ∞ the GS reduces to a symmetric mixed state |N N | + |N N |.
Once the question of the GS is accounted for, from the scenario above one can see that the orientation of spins with nonzero component along z is determined by the hidden order and can be taken for granted. The validity of such an approximation is ultimately measured by the values of the longitudinal SOP: the closer is O z to unity the higher is the AFM order of nonzero spins. We then introduce the following fermionic picture: assign a spinless fermion |+ i ≡ c † i |− i at site i if S z i = 0 and no fermions |− i ≡ |0 i in the spin language if S z i = 0. (This notation for spinless fermions has a direct translation in the language of the XY model that will be introduced at the end of this section.) Process ii) is nothing but a hopping of spinless fermions, while processes i) and iii) represent annihilation and creation of pairs | + i + i+1 . The density of nonzero spins (S z i ) 2 is simply translated to the local fermion number n i = c † i c i , while, due to the underlying AFM order, the Ising-like term takes the form −λn i n i+1 that contributes with a negative energy when two fermions are present on adjacent sites.
Hence, under the hypothesis of hidden AFM order, the dynamics of equation (2) is reproduced by the following effective fermionic model (acting in a reduced Hilbert space i where H (2) denotes the local Hilbert space of a two-level system -as that of a spinless fermion or a spin-1/2 introduced below). It should be observed that equation (4) with D = 0 is essentially equivalent (apart from an additive constant) to equation (2) of ref. [15] once a particle-hole transformation n i → 1 − n i is performed at every site.
Following Gómez-Santos [15] we now proceed to a further approximation on the fermionic Hamiltonian that is not amenable to an exact treatment due to the λ-term. At the Hartree-Fock level this term can be approximated as: where the expectation values . . . now are taken with respect to the GS of the quadratic Hamiltonian where the parameters have to be determined self-consistently. The advantage of a Hamiltonian of the form (5) is that it can be diagonalized by means of a Bogoliubov transformation and θ k is given by Note that, as we are interested in the thermodynamic limit, we do not specify here the boundary conditions on the spin and fermionic Hamiltonians. The momenta are quantized as ∆k = 2π/L and their precise location within the first Brillouin zone depend on the conditions imposed on the end sites. However, for L → ∞ Apart from additive terms of O(L −1 ) the Hamiltonian in diagonal form is In the thermodynamic limit, the selfconsistent equations are The notation used in equation (6) is the one commonly used for the XY spin-1/2 model in a transverse field. In fact, by (inverse) Jordan-Wigner transform one gets [12] where the σ α j 's are Pauli matrices at site j. This model is known to be critical at h = ±1 for γ = 0, where it belongs to the c = 1/2 universality class (the same as the 2D classical Ising model) and at γ = 0 for h ∈ (−1, 1) where the universality class becomes that of the compactified free boson, c = 1.
From the numerical solutions of (7)-(9) it turns out that in the Haldane and Néel phases A < 0 and B < 0 so that γ > 1 (as long as λ|A| < 1 -some representative cases are listed in table 2), while most studies are limited to |γ| < 1. As a consequence the region γ 2 < 1 − h 2 corresponding to oscillations with wavenumber different from π [16,12] is not present in our case. However, having γ > 1 does not affect the critical condition we are interested in, that remains h = ±1. In these cases we have just c = 1/2, as reported before [9] for the Haldane-to-Néel transition. At λ = 0, D ∼ = −2 this transition line merges with the boundary towards the so-called XY phases corresponding to c = 1. Interestingly this change of universality class is captured also by our approximation since for λ = 0, D = −2 the self-consistent solution yields just γ = 1 and h = 1. From the data in table 2 one can also estimate, for example, the critical value of D at fixed λ = 1; the result is D c ∼ = −0.214, which is not in good quantitative agreement with the numerical value D c = −0.315 [9,14]. The perturbation of the isotropic Heisenberg Hamiltonian with λ > 1 and D = 0, instead, seems to be better described by the spinless fermions approach; already at this level of approximation the value λ c = 1.125 found in [15] is close to our best DMRG independent estimate λ c = 1.1856 [17]. Even if it is likely that the inclusion of configurations with nearest-neighbour parallel spins could improve the results, as discussed by Gómez-Santos [15], we do not insist along this line here because we are ultimately interested in the decay laws of correlation functions that are essentially dictated by the universality classes. In fact, it is important to stress that neither the extension to D = 0, nor the extension of the model as in equation (8) of ref. [15] modify the universality class of the transition, that remains of the c = 1/2 (or Ising) type for λ > 0. Although the location of the critical points and of the prefactors depend on the values of the parameters, the scaling dimensions of the operators in the Table 2: Self-consistent estimates of the three decoupling parameters n 0 , A and B of equations (7)-(9) for some choices of λ and D in the Haldane and Néel phases. It must be kept in mind that the continuum versions of the self-consistent equations neglect some O(L −1 ) terms coming from isolated contributions at wavenumber 0 or π. Last two column contain the corresponding parameters h and γ of the effective XY model according to equation (6). continuum field theory (i.e. the decay exponents of the correlation functions) do not change when we move along the c = 1/2 line. Nonetheless, due to the lack of an explicit mapping of the spin-1 strings onto the corresponding correlators in the Ising fermionic field theory, up to now the exponents appearing in the large-distance decay of string correlation functions were unknown. This is precisely the subject of subsecs. 3.1 and 3.3. Eventually, we note that alternative pictures of the Haldane gap in fermionic language can be derived by perturbation theory near the Babujian-Takhtajan integrable biquadratic spin-1 chain [18] or from two-leg ladders with ferromagnetic coupling on the rungs [19]. At this stage it is interesting to compare the entanglement properties of the original spin-1 model (eq. (2)) with those of the XY spin-1/2 chain resulting from the mapping. On the one hand, for the former it has been shown [20] that at the isotropic Heisenberg point λ = 1, D = 0 there is long-distance spin-1 (qutrit) entanglement in the thermodynamic limit for two sites arbitrarily far apart. It is reasonable to expect that this entanglement survives in a neighbourhood of the isotropic point. On the other hand, in ref. [21] it is stated that the qubit entanglement in the XY model with transverse field vanishes beyond a distance of order γ −1 . In our case γ > 1 and the degrees of freedom of the qubits represent the presence or the absence of an effective particle with |S z i | = 1. Therefore we are led to speculate that wherever there is full spin-1 entanglement in the vicinity of the Heisenberg point, this is due to the spin correlations between the sites with S z i = 0. Recalling the hypothesis of underlying string order and imagining to eliminate the sites with S z i = 0, the qualitative picture of the long-distance entangled states in the Haldane region if that of a Greenberger-Horne-Zeilinger state [22] with effective AFM order | . . . ↑↓↑↓ . . . + | . . . ↓↑↓↑ . . . .

Mapping for the spin-spin and string correlators
We shall exploit now the mapping from spin-1 to spinless fermions, based on the existence of an underlying string order, to translate the various spin-1 correlation functions onto expectation values of strings of fermionic operators that can be computed exactly when the Hamiltonian has the form (5). Let us start from the z-component of the spin is a Jordan-Wigner tail that accounts for the correct sign when S z j = 0 assuming, conventionally, that the first nonzero spin is pointing up. By inserting the expression S z j = 1+σ z j 2 i<j (−σ z j ) into the definition of the longitudinal spin-spin correlation function and using the properties of Pauli matrices one finds As far as the transverse correlation functions are concerned, it can be checked by direct inspection on a generic configuration with perfect string order that the identification has the correct action, since the only cases in which the l.h.s. does not break the string order are those with S z k = 0, that is σ z k = −1, on all sites between j and j + R. The product on the r.h.s. of (13) is exactly the expression involved in the so-called emptiness formation probability (see, for example, [23] and refs. therein).
Let us now study the spin-1 strings. Along the z-direction we have simply Again by using the relation S z j = 1+σ z j 2 i<j (−σ z j ) and plugging the string written above into eq. (1) one gets Note that in the language of the effective XY model, the Néel correlation function (12) involves a string of Pauli operators whereas the string correlation function (16) involves only one-and two-points correlators of the σ's. Thanks to equation (14) we easily obtain also the pure-string correlation function as: From equation (12) we see that, in this approach, G H (R) is nothing but the first term of the usual spin-spin correlation function C z (R) apart from the prefactor.
Along the x-direction, instead, we exploit the fact that that is, apart from an overall sign, the operator above performs a swap between S z j = 1 and S z j = −1 leaving the S z j = 0 component isolated. This swap is important because it can be checked by direct inspection that, for every possible combination of the spin at sites j and j + R that respects the hidden AFM order, both with an even or an odd number of nonzero spins in between, the action of O x (R) produces only one allowed configuration and some other forbidden ones. More precisely, using the spin-1/2 operator σ x that changes empty sites into occupied sites and viceversa we can write where the inner spin-1 transverse string contributes with the sign prefactors. Thanks to hidden order in our reduced Hilbert space, the spin-1/2 configurations generated by σ x represent the allowed spin-1 states and the forbidden ones are automatically filtered out. The coefficient 1/2 comes from the matrix elements of S x at sites j and j + R. Now, thanks to the fact that the Hamiltonian (5) is quadratic in the fermionic operators, all the correlation functions can be evaluated using Wick's theorem. Following the notation of the seminal paper by Lieb, Schultz and Mattis [24] we introduce the operators A j = c † j + c j and B j = c † j − c j that allow to express the basic two-point correlations as If we further assume translational invariance (i.e. PBC) along the chain we have where and O x (R) is given directly by σ z j σ z j+R and σ x j σ x j+R respectively. The ordinary correlators C x,z (R) require a step more since they involve strings of Pauli operators. For example, each of the terms in equations (12) and (17) has the form k B k A k . When R → ∞ all the four terms in equation (12) tend to coincide so that Exploiting Wick's theorem, Caianiello and Fubini [25] have shown that the expectation value above can be expressed as a Pfaffian Thanks to the fact that Q l =m = S l =m = 0 this Pfaffian reduces to a Toeplitz determinant [26] and we get So, the determinants of the matrices with entries G j becomes the central quantities of our analysis. The matter is more complicated for the transverse spin-spin correlator essentially because it eventually involves a Toeplitz determinant generated by a matrix-valued symbol that may also become singular. According to ref. [27] this case is not yet solved in the theory of Toeplitz determinants and in ref. [28] it has been suggested to extend directly the procedure valid in the nonsingular case. Fortunately in our case a workaround is possible: thanks to a suitable diagonalization, we are able to complete the calculation of the dominant contribution to C x (R) in terms of a product of Toeplitz determinants, each one computed using the Fisher-Hartwig conjecture (see, for instance, App. A.2 in [23]). The details of this procedure are reported in the Appendix.

Longitudinal string correlation O z (R)
The first object we will compute is the longitudinal string correlator. From equations (16) and (19) we get Following Barouch and McCoy [16] (in particular their eq. (6.12)) we express G R as follows with and .
Note that since γ > 1 the two roots of the numerator are always real; the behaviour for R → ∞ is controlled by λ 2 . From equations (6.17), (6.14) and (6.18) in ref. [16] we have, respectively: In every case the asymptotic value O z = 0 is simply interpreted as a non-saturated value of the magnetization along z in the XY model in transverse field where σ z j = G 0 (h, γ) can be computed using equation (21) at R = 0.

Longitudinal spin-spin correlation function C z (R) and pure string correlator G H (R)
The asymptotic behaviour of the Toeplitz determinant in equation (20) can be found using the same technique as in [29], since (apart from a sign) the generating symbol c(e ik ) is essentially the same used by Wu. Then we find: • Haldane phase h < 1 (λ 2 > 1) which corresponds to the known decay behaviour at the isotropic Heisenberg point, as predicted by the nonlinear σ-model approach (see, for example, [30]). Moreover, in refs. [31,32] it was argued that the same behaviour of the connected longitudinal correlation function persists also in presence of a staggered magnetic field; in this sense such a behaviour could be considered a signature of the Haldane phase, robust against anisotropic perturbations.
• Néel phase h > 1 (0 < λ 2 < 1) Apart from the nonzero asymptotic value for h > 1, that serves as an order parameter for the Néel phase (the ordered phase T < T c in Wu's paper [29]), it must be noticed that both the power of R in the denominator and the exponential constant are different on the two sides of the transition. The roots λ 1,2 and the bulk correlation length ξ are the same as in subsec. 3.1 (see eq. (22)).

Transverse string correlation function O x (R)
• Haldane phase h < 1 (λ 2 > 1). The nonzero asymptotic value O x comes from the long-range order lim R→∞ σ x 0 σ x R in the XY model with spontaneous breaking of the symmetry σ x → −σ x . The result can be borrowed directly from equation (4.1) of [16] O with ξ having the same meaning of subsec. 3.1.
• Critical line h = 1 (λ 2 = 1). There is no long-range-order in σ x 0 σ x R , that decays to zero as R −1/4 as expected from the scaling dimension 1/8 of the primary operator in the c = 1/2 CFT [12]. Using equation (4.7) in [16] • Néel phase h > 1 (0 < λ 2 < 1). Equation (4.25) in ref. [16] is translated to We should stress that the critical exponent in equation (25) differs from the one in equations (24) and (26); it is not possible to recover the decay behaviour at h = 1 from the functions found for h > 1 or h < 1 simply by letting R/ξ → 0 in the exponentials. Qualitatively, the reason is that the correlation functions should be described by a unique scaling function F (r) of the variable r = R/ξ, but the asymptotic expansions in the off-critical regime and in the critical regime are different. The former corresponds to r ≫ 1 while the latter to r → 0 for any large but finite value of R. A similar argument holds also for the longitudinal spin-spin correlation function C z (R) of the previous subsection. Although possible in principle (see, for instance, [33] for the 2D classical Ising model), the derivation and the usage of the whole scaling functions is beyond the scope of this paper. Finally, we notice from the equations above, as compared to those of subsec. (3.2), that the correlators O x (R) and C z (R) play a dual role above and below the transition line; when one order parameter is vanishing, the other is not. Here we do not have an explicitly duality relation between order and disorder lattice operators for the spin-1 model as in the Ising case (see, however, ref. [34] for the XY chain). Hence, what is a nontrivial fact is to see that also the decay laws interchange when the transition line is crossed.

Transverse spin-spin correlation C x (R)
From the analysis reported in the Appendix, we can prove a conjecture already put forth in ref. [15], namely that the transverse correlation function decays always exponentially even when one crosses the critical line. Here we can be more precise and derive also the power-law terms in front of the exponential where β and β ′ in the Haldane phase, along the critical line and in the Néel phase take, respectively, the form written in equations (37), (41), (38), (42), (39) and (43) in the Appendix. In particular we have checked that both β c and β ′ c for h = 1 and γ ≥ 1 are nonzero. Hence, despite the fact that the system is critical, the transverse correlation function exhibits a finite characteristic length Ξ. As far as the exponent η x is concerned: • Haldane phase h < 1 (λ 2 > 1) and Néel phaseh > 1 (λ 2 < 1): η x = 1/2; the algebraic prefactor is the same power-law that describes critical correlations in the quantum Ising model.
From the values of h and γ reported in the last two columns of table 2 we have computed Ξ[h(λ, D), γ(λ, D)]; for example when λ = 1 we find that Ξ decreases steadily as D is decreased towards larger negative values, passing from the Haldane to the Néel phase. This behaviour is consistent with the numerical best-fit estimates of Ξ made in the next section.

Comparison with DMRG results
The results of the previous section regarding the long-distance decay of ordinary and string correlation functions are summarized in table 3 where: Within the approximation of hidden order and for large R these asymptotic laws are to be considered exact and valid for the Haldane and Néel phases and associated transition line as specified in table 3. It should be noted that f 1 (R) and f 2 (R) agree with the general form argued for the d-dimensional Ising model (see, for example, [35] and refs. therein) although the derivation of the latter did not include the case of string correlation functions.
The two functional forms f 1,2 now can be used to extract, for example, the asymptotic value of string order correlation functions computed numerically; in this sense A 0 , A 1 and A 2 may be regarded as free fitting parameters. The goodness of the best-fit procedure can be assessed by computing the reduced χ 2 : squares of differences # of data points − # of fit parameters − 1 .
Clearly one could use many other different functions to extrapolate the correlators to R → ∞. However, as recalled in the introduction, the literature contains very few, empirical, information about the asymptotic approach to the limit values of the string correlation functions. Our study was motivated by this fact and so here we perform a comparison between f 0 , f 1 and f 2 by examining their capability to fit the spin-spin and string correlations evaluated numerically through the DMRG. Actually, following Table 3: Expected asymptotic behaviour of string (O) and usual (C) correlation functions in the Haldane and Néel phases of model (2) and along the critical transition line separating them. The fitting functions f 0,1,2 are defined in equations (28)- (30). Note the interchanged role of O x (R) and C z (R) above and below the transition line. Phase C.f. Decay law the idea of ref. [30], in order to take into account the PBC we employ the left-right symmetrized expressions of equations (28)-(30) at least for the correlation functions in the z-channel. As regards O x , while translational invariance implies that it depends on the difference between the sites at the ends of the string, it is not always guaranteed that it depends only on the distance on the ring. In other terms the expectation value may differ from the same expression with i and j interchanged. In fact, using the properties of the exponentials of spin-1 operators, exp(iπS x i ) = exp(−iπS x i ) and S x i exp(iπS x i ) = −S x i , it can be shown that the expression above can rewritten as The point is that in general the GS of an anisotropic spin chain is not invariant under the action of exp(iπS x tot ) and so a direct inspection is required case by case in order to decide if a symmetrized fitting function has to be used or not.
The asymptotic limits (i.e. the values of A 2 ) resulting from of a series of bestfits made on DMRG data obtained by fixing λ = 1 and letting D to vary across the Haldane-Néel transition from −0.125 to −0.875 are plotted in figure 2. It is seen that the nonvanishing values of O x and M z characterize, respectively, the Haldane and the Néel phase. It is reasonable to expect that the location of the critical point as the value of D at which the two order parameters vanish leads to two slightly different estimates. However, with more accurate methods the critical point was previously found to be D c = −0.315 [9,17].
In the simulations we have fixed the total length of the chain to be L = 100 sites and computed the GS properties by retaining from 243 to 324 DMRG states in the sector with S z tot = 0, which is the only good quantum number that we could use. All the functional forms derived above are asymptotic so we cannot expect them to be reliable for very short distances. Therefore, we have conventionally excluded the data with R ≤ 5 from the fitted points. In the Néel phase the GS tends to become doubly degenerate in the limit L → ∞; in order to take into account this difficulty we have built the reduced density matrix by targeting the two low-lying states rather than just the GS. Finally we have performed three finite-system sweeps to achieve a better accuracy. In the cases we have considered, the transverse string correlation O x (R) turned out to be symmetric with respect to the middle of the chain except for D = −0.75 and D = −0.875. For this reason we have repeated the fit using directly the functions of equations (28), (29) and (30) without symmetrization selecting only the points in the first half of the chain. The asymptotic values are essentially unaffected, with the exception of those referring to the critical point. In general when the results of the fit are such that A 1 ≫ L (typically close to criticality) we conclude that the exponential tail of the fitting function is essentially saturated to unity and an algebraic fit would produce the same result.
As far as the best-fitting functions for C z (R) and O x (R) are concerned, the passage from the type f 1 to the type f 2 going through f 0 at the critical point, as in table 3, actually takes place gradually in the interval D ∈ (−0.345, −0.315), the worst values ofχ 2 being of order 10 −5 . The best choice to fit the transverse spin-spin correlation function C x (R), instead, follows the prediction of table 3 (F 1 except at the critical line, where it becomes F 0 ) with a deviationχ 2 < 10 −8 . Finally, the longitudinal string correlator O z (R) is very well fitted by F 2 , in agreement with table 3, withχ 2 ∼ 10 −9 or better. It is also important to check quantitatively the goodness of the Hartree-Fock approximation. The decoupling parameter in the fermionic version are n 0 = n j=0 , A = c † 1 c 0 and B = c 1 c 0 where we have selected a reference site "0" invoking translational invariance. In the original spin-1 formalism it can be checked directly that The operator c 1 c 0 destroys a couple of fermions in adjacent sites; in the spin language they could be ↑↓ or ↓↑ depending on the surrounding sites in order to respect the AFM order. Let us express the GS in the form |GS = α| ⇑ + β| ⇓ , where | ⇑ denotes a linear combination of states in which the first nonzero spin along-z is directed upward and | ⇓ the same state with all the spin reversed. Only one of the terms in S + 0 S − 1 + S + 1 S − 0 will act on | ⇑ respecting the AFM order and the other term will thereby act on | ⇓ . When the scalar product with GS| is taken, the states from | ⇑ will not mix with those from | ⇓ . Therefore we try with the expression In table 4 we report the values of the decoupling parameters for a set of points in the Haldane and Néel phases, comparing the DMRG values with the numerical solution of the self-consistent equations using 100 iterations from different choices of initial conditions. Having the DMRG estimates for n 0 and A we may also produce a "hybrid" estimate of the critical point by setting h = 1 in equation (6) and then solving forD c (λ) = 2[λ(n 0,DMRG − A DMRG ) − 1]. With λ = 1 we find for ex-ampleD c = −0.254, that compares slightly better than the fully-self-consistent value (D c = −0.214) to the accepted numerical one D c ∼ = −0.315.
Apart from the value n 0 , which quantifies the number of spins with nonzero projection along z, we expect that the goodness of the mapping used in this work is higher when the hidden order is larger. Therefore, as a final check, we have repeated the passages of Sect. 2 (see [36] for details) including also a biquadratic term 1 3 i ( S i · S i+1 ) 2 in the spin-1 Hamiltonian (2). For λ = 1 and D = 0 the ground state of the model can be found exactly [11,1] using the valence-bond picture: each spin-1 is viewed as the triplet sector of a pair of spin-1/2 particles and the ground state is constructed by creating a sequence of singlets between adjacent sites. In this case the string correlation functions can be computed exactly and it turns out that O x,z (R) = −4/9 independent of R. At the isotropic point with biquadratic term the self-consistent equations are solved by n 0 = 2/3, A = B = −2/9 and the effective parameters of the XY model become h = 3/5 and γ = 4/5. From equation (22) we find λ 1 = λ 2 = 3 so that in Table 4: DMRG (L = 100) versus self-consistent (s-c) estimates of the three decoupling parameters n 0 , A and B of equations (32) and (33). It must be kept in mind that the continuum versions of the self-consistent equations neglect some O(L −1 ) terms coming from isolated contributions at wavenumber 0 or π.

Conclusions
In this paper we have reconsidered and extended the approach of ref. [15] to the GS properties of spin-1 anisotropic quantum chains. We have included a single-ion term in the Hamiltonian and, moreover, we have analyzed explicitly how the spin-1 correlation functions are written in the spinless fermions language and then in the framework of the XY model in a transverse field for effective spin-1/2 degrees of freedom. In particular, we have focused on the decay laws of the string correlators towards their asymptotic values which apparently were missing in the literature.
The decay laws of string and spin-spin correlation functions (in the longitudinal channel) are all related to the generating function c(e ik ) of equation (21) and the determinants of the Toeplitz matrices derived from it. The asymptotic behaviour of the transverse correlation function C x (R), instead, originates from a product of two Toeplitz determinants (see the Appendix, in particular eq. (35)). The leading terms in the regime R ≫ 1 for the various correlators are discussed in Sect. 3 and summarized in table 3. In brief the most interesting points unveiled by the approach used here are: • The nonvanishing string-order parameters of the spin-1 model (2) are simply interpreted as the magnetization along x and z in the XY chain with transverse field (eqs. (23) and (24)).
• There exists dual behaviour of O x (R) and C z (R) above and below the transition, both for the asymptotic order parameters and for the decay functional forms.
• The explicit calculation of C x (R) allows us to prove an unusual feature in statistical mechanics, already conjectured by Gómez-Santos [15]: the spin-spin transverse correlation function exhibits always a finite characteristic length Ξ (eq. (27)) even when the system becomes critical.
The analytical results are supported by comparison with a numerical (DMRG) study of the model, especially for the correlations C x and O z . A more detailed comparison between the analytical and the numerical estimates should take into account: i) finite-size effects due to a finite total length L while in Sect. 2 we passed readily to the thermodynamic limit; ii) corrections for finite distance R beyond the dominant ones. Although in principle they can be computed systematically, in this paper we have limited ourselves to the leading terms in order to derive analytical expressions with the smallest possible number of fitting parameters.

Appendix: Toeplitz formulation of C x (R)
The fermionic version of equation (13) reads By observing that we can write C x (R) = (−1) R S + j S − j+R + S − j S + j+R /4 (using translational and U(1) rotational invariance about z) as where the two terms come, respectively, from the Pfaffians and It is useful to note that F l−j = −F j−l and H l−j = H j−l . According to usual conventions, the Toeplitz matrix (34) is generated by matrix-valued symbol (analytically continued to the unit circle) Now, since the trace norm of M 0 is vanishing, when R → ∞ we get det M 1 = det M 2 so that so that C x (R) ≃ −2 √ det M 1 . In order to compute the determinant, we need to check the possible zeroes of det φ(z) = 1 2   1 + sign(z) z 2 + 1 − 2hz for instance, when z = ±1, det φ(±1) = 1 2 (1 + sign(±1)sign (1 ∓ h)) so that for h < 1 the symbol is singular at z = −1, while for h > 1 it is singular at z = +1. Unfortunately, as discussed also in ref. [27], known results for matrix-valued symbol do not cover the case of singular symbols with vanishing determinant. Hence, the strategy is to factorize the determinant of M 1 as a product of determinants of matrices generated by scalar-valued symbols. Fortunately, in this case this task is accomplished by transforming M 1 through the matrix Now, we first use a theorem by Widom and Silbermann (see, for instance, [37,38]) according to which F −1 is a Toeplitz matrix generated by f −1 (in the present case this result holds for even R) and then express the product H T F −1 H as another Toeplitz matrix generated by the symbol ih(z −1 )f −1 (z)h(z). The last identification can be done by using repeatedly a theorem by Brown and Halmos [39]: T (ϕ)T (ψ) is a Toeplitz operator iff either ϕ * (z) or ψ(z) are analytic functions; if the latter condition is satisfied then T (ϕ)T (ψ) = T (ϕψ) where T (ϕ) denotes the Toeplitz matrix generated by the function ϕ.
Let us start by computing det(−iF): • Haldane phase h < 1 (λ 2 > 1). From the analytic continuation to the unit circle we see that f vanishes at z = ±1 and is singular at z = λ −1 • Haldane phase h < 1 (λ 2 > 1). The Fisher-Hartwig conjecture now can be applied, thanks to the following factorization (40) with α 1,2 and β 1,2 as in table 6 and with satisfying Szegö's theorem. Consequently, the asymptotic behaviour is purely exponential: • Critical line h = 1 (λ 2 = 1). There are no singularities and a simple zero at z = −1, with exponents α and β as in the first row of table 6. Therefore, the net power of R in the algebraic prefactor vanishes and the decay is purely exponential with 2π 0 dk ln cos k − 1 + (cos k − 1) 2 + (γ sin k) 2 + ln γ 2 .