The quantum Ising chain for beginners

We present here various techniques to work with clean and disordered quantum Ising chains, for the benefit of students and non-experts. Starting from the Jordan-Wigner transformation, which maps spin-1/2 systems into fermionic ones, we review some of the basic approaches to deal with the superconducting correlations that naturally emerge in this context. In particular, we analyse the form of the ground state and excitations of the model, relating them to the symmetry-breaking physics, and illustrate aspects connected to calculating dynamical quantities, thermal averages, correlation functions and entanglement entropy. A few problems provide simple applications of the techniques.


Jordan-Wigner transformation
The quantum Ising chain is an ideal playground for testing many of the ideas of statistical mechanics, including recent non-equilibrium aspects. As such, it is a standard test case in many recent papers. These notes are intended for students willing to begin working with quantum Ising chains. They can be also useful as a practical guide to researchers entering the field. We, unfortunately, do no justice to the immense literature where concepts and techniques were first introduced or derived, and even less so to the many papers where physical applications are presented: this would require too much effort for our limited goals. Therefore, our references -except for those useful for the problems provided -are reduced to the bare minimum, and we apologise for that with authors whose work is not duly cited. To remedy for this omission, most of the topics include detailed derivations which should make these notes reasonably self-standing.
The level of our presentation is roughly appropriate to graduate students, but also master students should be able to follow most of the developments, provided they acquire the necessary pre-requisites: second quantization [1] to deal with bosons and fermions, and basic knowledge of quantum mechanics of the spin-1/2 [2].
We start, in Sec. 1, from the Jordan-Wigner transformation. Next, we give in Sec. 2 a fermionic formulation of the transverse field Ising model in one dimension. Section 3 treats the ordered case with periodic boundary conditions, where a simple analytical reduction to an assembly of 2 × 2 problem is possible. In Sec. 4 we give the Nambu formalism for the general disordered case. Section 5 shows how to diagonalise the Hamiltonian in the general disordered case, while Sec. 6 deals with the dynamics for a time-dependent Hamiltonian. Next, in Sec. 7 we show how to calculate the overlap between different Fock states for two different Hamiltonians. In Sec. 8 we show how to calculate thermal averages. Sections 9 and 10, finally, contain the technicalities related to the calculation of correlation functions involving Jordan-Wigner string operators and the entanglement entropy.

Jordan-Wigner transformation
For systems of bosons and fermions, a large assembly of many-body techniques has been developed [1]. A difficulty in dealing with spin systems is that they are neither bosons nor fermions.
Consider, to start with, a single spin-1/2, and the three components of the spin-operators represented in terms of the usual Pauli matrices 1σα with α = x, y, z. The Hilbert space of a single spin is two-dimensional: for instance you can write a basis as {| ↑ , | ↓ }, in terms of the eigenstates ofσ z , withσ z | ↑ = | ↑ andσ z | ↓ = −| ↓ . Moreover, ifσ α j denote Pauli matrices at different lattice sites j, hence acting on "different" (distinguishable) two-dimensional Hilbert spaces, then σ α j ,σ α j = 0 for j = j .
However, if we decide to truncate such a Hilbert space to only two states, {|0 , |1 }, assuming (b † ) 2 |0 = 0, then the Hilbert space of a single spin-1/2 can be easily mimicked. Such a truncation, which can be thought as adding a large -ideally "infinite" -on-site repulsion term to the boson Hamiltonian, is known as hard-core boson. We transform the Pauli spin-1/2 operatorsσ α j (with α = x, y, z, and j a generic site index) into hard-core bosonsb † j , by identifying 2 at each site |0 ↔ | ↑ and |1 =b † |0 ↔ | ↓ . Recalling thatσ ± = (σ x ± iσ y )/2 act asσ + | ↓ = | ↑ , and σ − |↑ = |↓ , we must evidently have: These operatorsb † j commute at different sites -as the originalσ α j do -but are not ordinary bosonic operators. They anti-commute on the same site 3 b j ,b † j = 1 and they verify the hard-core constraint (b † j ) 2 |0 = 0, i.e., at most one boson is allowed on each site.
Info: The hard-core boson mapping might be viewed as a possibly useful way of rewriting spin-1/2 models in a rather general setting. For instance, if you have a Heisenberg model for spin-1/2 sitting on a lattice, whose sites are denoted by x and with nearest-neighbour pairs denoted by x, x we could write, definingn x =b † xbx : The second expression shows that we are dealing with hard-core bosons hopping on the lattice and repelling each other at nearest neighbours. Needless to say, this helps in no way in solving the problem. i The hard-core constraint seems to be ideally representable in terms of spinless fermionsĉ † j , where the absence of double occupancy is automatically enforced by the Pauli exclusion principle, and the anti-commutation on the same site comes for free. Unfortunately, whereas the mapping ofσ α j into hard-core bosonsb † j is true in any spatial dimension, writingb † j in terms of spinless fermionsĉ † j is straightforwardly useful only in one-dimension 2 This identification is not unique, as you can swap the two states. 3 Since on the same site σ + j ,σ − 1 Jordan-Wigner transformation (1D), where a natural ordering of sites is possible, 4 j = 1, 2, · · · , L. In other words, because fermion operators on different sites must anti-commute, the exact handling of the resulting minus signs -which are absent in the original spin problem -is very natural only in 1D.
The Jordan-Wigner transformation of hard-core bosons into fermions reads: where the non-local string operatorK j is simply a sign,K j = ±1, counting the parity of the number of fermions which sit before site j (i.e., between site 1 and site j − 1), multiplying the fermionic operatorĉ j . Notice thatK j =K † j =K −1 j , andK 2 j = 1. Heren j =ĉ † jĉj is the fermion number operator, but one can show that the phase-factorK j cancels out inn j , i.e.,n j =ĉ † jĉj =b † jbj . Let us prove that everything works fine. More precisely, we will now show that if theĉ j are taken to be standard fermionic operators, with canonical anti-commutation relations ĉ j ,ĉ † j = δ j,j and ĉ j ,ĉ j = ĉ † j ,ĉ † j = 0, then the following two properties of theb j follow:

P1
: which is a formal way of writing that theb j are hard-core bosons. Property P1 is straightforward because the string factorK j cancels completely: and, similarly,b jb † j =ĉ jĉ † j . In essence, on each site theb j inherit the anti-commutation property P1 from the fermionsĉ j . To prove P2, let us consider b j1 ,b † j2 , assuming j 2 > j 1 . Using Eq. (5), it is simple to show thatb j1b † j2 =ĉ j1 e −iπ j 2 −1 j=j 1n jĉ † j2 .
Similarly, you can show that where the change of sign in the second line is due to the fermionic anti-commutation, while the crucial final change of sign is due to the fact thatn j1 = 0 at the beginning of the second line, whilê n j1 = 1 in the final expression, because of the neighbouring action ofĉ j1 . Comparing Eq. (7) with Eq. (8), you immediately deduce that b j1 ,b † j2 = 0. All the other commutation relationships in P2 are proven similarly.
Here is a summary of a few useful expressions where the string operatorK j disappears exactly: 1 Jordan-Wigner transformation Notice the minus signs on the right-hand side, which should not be forgotten. Notice also that we have used (1 − 2n j ) = 1 − 2n j .
Jordan-Wigner transformation. Summarising, spins are mapped into fermions using: Armed with these expressions, it is simple to show that some spin operators transform in a simple way into local fermionic operators. Here is a summary: Unfortunately, a longitudinal field term involving a singleσ x j cannot be translated into a simple local fermionic operator.
i One important point to note concerns boundary conditions. One often assumes periodic boundary conditions (PBC) for the spin operators, which means that the model is defined on a ring geometry with L sites, j = 1, · · · , L, and the understanding thatσ α 0 ≡σ α L andσ α L+1 ≡σ α 1 . This immediately implies the same PBC conditions for the hard-core bosons. Hence, for instance: But observe what happens when we rewrite a term of this form using spinless where is the total number of particles (fermions or bosons, it does not matter), and the second equality follows because to the left ofĉ † L we certainly haven L = 1, and therefore e iπn L ≡ −1. Similarly, you can verify that:b Info: There is a whole class of one-dimensional spin systems where a fermionic re-formulation can be useful. Probably the most noteworthy is the XXZ Heisenberg chain, which would read The corresponding fermionic formulation would read: which shows that the fermions interact at nearest-neighbours, due to the J zz j -term.
i Let us now concentrate on a class of one-dimensional models where the resulting fermionic Hamiltonian can be exactly diagonalised, because it is quadratic in the fermions: such a class includes the XY model and the Ising model in a transverse field. After a rotation in spin-space, we can write the Hamiltonian (allowing for non-uniform, possibly random, couplings) as follows: whereσ α j are Pauli matrices. The couplings J x,y j and the transverse fields h j can be chosen, for instance, as independent random variables with uniform distribution. For a system of finite size L with open boundary condition (OBC), the first sum runs over j = 1, · · · , L − 1, or, equivalently, we would set J x,y L = 0. If periodic boundary conditions (PBC) are chosen, the sum runs over j = 1, · · · , L and one assumes thatσ α L+1 ≡σ α 1 . For J y j = 0 we have the Ising model in a transverse field, for J y j = J x j the isotropic XY model.
In terms of hard-core bosons, the Hamiltonian becomes: where we have introduced a shorthand notation 5 J ± j = J x j ± J y j .
Next, we switch to spinless fermions, since all terms appearing in the previous expression do not involve explicitly the string operator K j . In terms of fermions, the Hamiltonian is essentially identical. We would remark that in the fermionic context the pair creation and annihilation terms are characteristic of the BCS theory of superconductivity [3]. The only tricky point has to do with the boundary conditions. If one uses open boundary conditions, the first sum runs over j = 1, · · · , L − 1 and there is never a term involving site L + 1, hence we have: In the PBC-case, terms likeb † appear in the Hamiltonian, where N is the number of fermions operator. Therefore: 5 This notation should not generate confusion with the angular momentum ladder operators. Here there is no imaginary-unit i, and the couplings J ± j = J x j ± J y j are just real numbers.
Info: Notice that, although the number of fermions N is not conserved by Hamiltonian in Eq. (21), its parity (−1) N = e iπ N is a "constant of motion" with value 1 or −1. So, from the fermionic perspective, it is as if we apply anti-periodic boundary conditions (ABC), hencê c L+1 = −ĉ 1 , if there is an even number of fermions and periodic boundary condition (PBC), henceĉ L+1 =ĉ 1 , if there is an odd number of fermions. This symmetry can also be directly seen from the spin Hamiltonian in Eq. (18), where one should observe that the nearest-neighbour σ x jσ x j+1 andσ y jσ y j+1 can only flip pairs of spins, hence the parity of the overall magnetisation along the z direction is unchanged. Such a parity can be easily and equivalently expressed as: We remark that P flips all theσ x j andσ y j , i.e., Pσ x,y j P = −σ x,y j , in the Hamiltonian in Eq. (18), leaving it invariant. This parity symmetry is the Z 2 -symmetry which the system breaks in the ordered ferromagnetic phase, as we will better discuss later on. i Let us define the projectors on the subspaces with even and odd number of particles: With these projectors we can define two fermionic Hamiltonians acting on the 2 L−1 -dimensional even/odd parity subspaces of the full Hilbert space: in terms of which we might express the full fermionic Hamiltonian in block form as: Observe that if you write a fermionic Hamiltonian of the form: then you can regard H 1 as a legitimate PBC-fermionic Hamiltonian since with the interpretationĉ L+1 ≡ĉ 1 . Similarly, H 0 is a legitimate ABC-fermionic Hamiltonian where you should poseĉ L+1 ≡ −ĉ 1 . Neither of them, however, expresses the correct fermionic form of the PBC-spin Hamiltonian. However, they are useful in expressing the fermionic blocks: since H p=0,1 conserve the fermionic parity, hence they commute with P 0,1 .
It is customary to parameterise J x = J(1 + κ)/2 and J y = J(1 − κ)/2, so that J + = J and J − = κJ. The Hamiltonian is then: 6 for the OBC case, and: We assume from now on that the number of sites L is even: this is not a big restriction, and is useful.
In the spin-PBC case, if the number of fermions N takes an odd value, then we effectively havê c L+1 ≡ĉ 1 ; if, on the contrary, N takes an even value, then the L-th bond has an opposite sign to the remaining ones, which can also be reformulated asĉ L+1 ≡ −ĉ 1 . Since the Hamiltonian conserves the fermion parity, both the even and the odd particle subsectors of the fermionic Hilbert space have to be considered when diagonalising the model, precisely as in the general case of Eq. (25). Introducing the two fermionic Hamiltonians as in Eq. (26) we now have: where we recall that p = 0, 1 is associated with the fermionic parity -p = 0 for even and p = 1 for odd parity -and that this compact way of writing assumes that the boundary terms are treated with:ĉ 6 Notice that one can change the sign of the h-term by making a particle-hole transformationc j → (−1) jĉ † j , which transformsñ j → 1−n j , and 1−2ñ j → 2n j −1, while leaving the hopping term untouched (same sign of J). With the current choice of the h-term, the h → +∞ ground state in the spin representation |↑↑ · · · ↑ is mapped into the fermionic vacuum, which will be useful in discussing the ground state. (Notice that the phase factor (−1) j exchange the roles of k = 0 and k = π in the discussion of the ground state.) Similarly, the same particle-hole transformation but without phase factor (−1) j would also invert the sign of the J-term, from ferromagnetic to antiferromagnetic.
Let us now introduce the fermion operators in k-space,ĉ k andĉ † k . The direct and inverse transformations are defined as follows: where the overall phase e iφ does not affect the canonical anti-commutation relations, but might be useful to change the phase of the anomalous BCS pair-creation terms (see below) . Which values of k should be used in the previous transformation depends on p. For p = 1 we haveĉ L+1 ≡ĉ 1 , which in turn implies 7 that e ikL = 1, hence the standard PBC choice for the k's: For p = 0 we haveĉ L+1 ≡ −ĉ 1 , which implies that e ikL = −1, hence an anti-periodic boundary conditions (ABC) choice for the k's: In terms ofĉ k andĉ † k , with the appropriate choice of the k-vectors, H p becomes: 8 Notice the coupling of −k with k in the (anomalous) pair-creation term, with the exceptions, for the p = 1 (PBC) case, of k = 0 and k = π, which do not have a separate −k partner. It is useful to manipulate the (normal) number-conserving terms 9 to rewrite the Hamiltonian as: The two terms with k = 0 and k = π, present for p = 1 (PBC), taken together can be written as: H k=0,π = −2J (n 0 −n π ) + 2h (n 0 +n π − 2) .
The remaining p = 1 terms, and all terms for p = 0, come into pairs (k, −k). Let us define the positive k values as follows: (39) 7 From the expression forĉ j in terms ofĉ k 8 We use the standard fact that the sum over j introduces a Kronecker delta for the wave-vectors: 9 We use that where we used the anti-commutation relations and k cos k = 0, and Then we can write the Hamiltonians as: where we have grouped together terms with k and −k in a single Hamiltonian of the form: Interestingly, the Hamiltonians H k live in the 4-dimensional space generated by the states: where they have a 4 × 4 matrix of the form: Check of dimensions. Recall that both H p=0,1 have 2 L eigenvalues. Indeed, there are L 2 such terms for H 0 , hence a dimension 4 L 2 = 2 L . Notice that H k=0,π also works in a 4dimensional subspace, and there are L 2 − 1 wave-vectors in K + 1 , hence again a total dimension for H 1 of 4 L 2 −1 4 = 2 L . Recall, finally, that the correct eigenvalues are obtained from the block Hamiltonians H p=0,1 which have 2 L−1 eigenvalues each, those with even (p = 0) or odd (p = 1) fermion parity.
with anti-commutation relations (α = 1, 2 stands for the two components of Ψ) We can then rewrite each H k as: where we have highlighted a 2 × 2 Hermitean matrix H k which can be expressed in terms of new pseudo-spin Pauli matricesτ x,y,z as: Here we recognise an "effective magnetic field" R k given by: Info: Observe now the role of the arbitrary phase φ introduced in the transformation from real-space to momentum space, Eq. (33). For φ = 0 the effective magnetic field lives in the y −z plane in pseudo-spin space, while for φ = π with corresponding eigenvectors (u k± , v k± ) T which can be expressed in terms of spin eigenstates in the direction R k /|R k |. From now on we will fix φ = 0, so that the pseudo-spin effective magnetic field lives in the y − z plane. Define the shorthand R k = (0 , y k , z k ) T with z k = 2(h − J cos k) and y k = 2κJ sin k. For the positive energy eigenvector, we have: where we have introduced the shorthands u k = u k+ and v k = v k+ . Note, in passing, that u −k = u k , while v −k = −v k , since z k is even in k, while y k is odd. The negative-energy eigenvector (u k− , v k− ) T is related to the previous one by a simple transformation: 10 The unitary matrix U k having the two previous eigenvectors as columns: diagonalises H k : So, define new fermion two-component operators Φ k through where, in the second term, we have made use of the fact that u −k = u k and v −k = −v k . It is straightforward to verify thatγ k is indeed a fermion. 11 In terms of Φ k = (γ k ,γ † −k ) T and 10 Indeed, write the eigenvalue problem for k+ = + k : Now change sign to the first equation, take the complex-conjugate of both, and rewrite them in inverted order, to get: which is the eigenvalue equation for (u k− , v k− ) T .
the last equality following from the normalisation condition for the eigenvectors.
The form of the two bands ± k , as a function of k and for several values of h is noteworthy. Figures 2-3 show some plots that illustrate them.  Winding and topology. It is instructive to trace the behaviour of the "effective magnetic field" R k , of magnitude |R k | = k , that the system "sees" as the wave-vector k spans the Brillouin zone [−π, π). Fixing φ = 0 in Eq. (49), this effective magnetic field lies in the y − z plane, R k = (0, y k , z k ) T with y k = 2κJ sin k and z k = 2(h − J cos k), where it draws the ellipse of equation a y 2 as k spans the interval [−π, π). We show in Fig. 4 three examples of this ellipse (circles, for κ = 1), one for |h| < J, one for h > J, and that for h = J. For |h| < J we see that the vector R k turns around and comes back to its original position, making one complete revolution around the origin, as k varies in [−π, π). We term the number of revolutions as the index [4] (or winding number) of the vector, and here it equals 1. As we change h in the range −J < h < J, the index, for continuity reasons, keeps the constant value 1 (it can only assume discrete values). In the case h > J, the vector R k makes no revolution around the origin and its index is 0: it keeps this value for any h > J, for the same continuity argument as before. The transition of the index between the two values 1 and 0 occurs at h = J. At that point, the continuity of the index as a function of the curve is broken, because the index is not defined for h = J, as the curve passes through the origin for k = 0. The index is a topological quantity, invariant under continuous transformations. Because it takes different values for |h| < J and h > J we say that these two phases have different topology. We see that R k = 0 corresponds to a degeneracy point of the 2 × 2 Hamiltonian H k -realised for k = 0 and h = J (but also for k = π and h = −J) -and the discontinuity of the index corresponds to the closing of the gap in the single-quasiparticle spectrum shown in Fig. 3.
a The ellipse degenerates into a segment for κ = 0, corresponding to the XY model. Hence, our argument requires κ = 0. The expression for H k in Eq. (59) allows to immediately conclude that the ground state of the Hamiltonian must be the state |∅ γ which annihilates theγ k for all k, positive and negative, the so-called Bogoliubov vacuum:γ In principle, one can define two such states, one in the p = 0 (even, ABC) sector, and one in the p = 1 (odd, PBC) sector. However, one finds that the winner between the two, i.e., the actual global ground state, is the one in the p = 0 (even) sector, with an energy The ground state can be written explicitly as: where |0 is the vacuum for the original fermions,ĉ k |0 = 0. We get where we used that u −k = u k and v −k = −v k . By normalising the state, we arrive at the BCS form: where the second more standard [3] expression applies to the case where φ = 0, with the choice of phase in which u k is real and v k purely imaginary.
The PBC-sector ground state must contain an odd number of particles. Since a BCS-paired state is always fermion-even, the unpaired Hamiltonian terms H k=0,π must contribute with exactly one fermion in the ground state. It is simple to verify that, with our choice of the sign of h > 0, the ground state hasn k=0 → 1 andn k=π → 0, resulting in an extra term of the form The PBC-ground state is, therefore: where we definedγ 0 =ĉ † k=0 andγ π =ĉ k=π for the unpaired states. The corresponding energy is: And here comes an amusing subtlety of the thermodynamic limit L → ∞. You would naively expect that, when you consider the energy-per-site e 0 = E 0 /L, then the ground state energy should simply tend to an integral: It turns out that the whole subtlety is hidden in the way one treats the boundary points at 0 and π: notice in particular that Eq. (62) for E ABC 0 involves L/2 k-points in the interval (0, π), while Eq. (68) for E PBC 0 involves L/2 − 1 points in the interval (0, π) and an extra term −2J. If you refrain from being too cavalier with the L → ∞ limit, you discover that the energy splitting is, in the whole ferromagnetically ordered region −J < h < J, a quantity that goes to zero exponentially fast when L → ∞: in other words, the two sectors ABC and PBC provide the required double degeneracy of the ferromagnetic phase: you can see that easily for h = 0. Less trivial, but true, for all |h| < J. On the contrary, ∆E 0 is finite in the quantum disordered regions |h| > J, ∆E 0 = 2(|h| − J), 12 and goes to zero as a power-law, more precisely as π/(2L), at the critical points h c = ±J. In Fig. 5 we illustrate these facts by numerically evaluating ∆E 0 .
Regarding the excited states, let us start from the p = 0 (even, ABC) sector. Consider, as a warm-up, the stateγ † k1 |∅ γ ABC . A simple calculation shows that, regardless of the sign of k 1 : In essence, the application ofγ † k1 transforms the Cooper-pair at momentum (|k 1 |, −|k 1 |) into an unpaired fermion in the stateĉ † k1 |0 . This would cost an extra energy + k1 over the ground state: the gain − k1 obtained from pairing is indeed transformed into a no-gain (energy 0) for the unpaired stateĉ † k1 |0 , consistently with the 4 × 4 structure of Eq. (43) predicting two eigenvalues 0 for the unpaired states. There is a problem with parity, however: a single unpaired fermion changes the overall fermion parity of the state. Hence, the lowest allowed states must involve two creation The energy of such an excitation is E ABC 0 + k1 + k2 , because we loose two Cooper pairs. Quite amusingly, if you consider the special caseγ † k1γ † −k1 you find that: This means thatγ † k1γ † −k1 transforms the Cooper pair at momentum (|k 1 |, −|k 1 |) into the corresponding anti-bonding pair: This costs an energy 2 k1 , consistent with the previous expression Generalising, we can construct all the excited states in the even-fermion sector, by applying an even number ofγ † k to |∅ γ ABC , eachγ † k costing an energy k . In the occupation number (Fock) representation we have therefore: We see that there are 2 L−1 such states, as required.
Remark: An important remark and check is here in order. First: the counting of the excitation number is evidently correct, if the k in Eq. (73) are allowed to range among the L positive and negative wave-vectors allowed by ABC: 2 L Fock states if the parity check is not enforced, 2 L−1 if we enforce parity. Second: recall that we can transform a Cooper pair of energy − k into the corresponding anti-paired state, of energy + k . The state that realises that isγ † kγ † −k |∅ ABC . Its energy is 2 k above that of the ground state, consistently with the formula given in Eq. (73), since −k + k = 2 k . ! In the p = 1 (odd, PBC) sector, some care must be exercised. One should apply an even number ofγ † k to the ground state |∅ γ PBC , including in the choice the unpaired operatorsγ † 0 , amounting to removing the fermion from the k = 0 state, andγ † π , amounting to creating a fermion in the k = π state.

Relationship with the spin representation
It is instructive to comment on the relationship between the spectrum we have found in the fermionic representation and the corresponding physics in the original spin representation. In this section we will focus on the Ising case, fixing the anisotropy parameter to κ = 1.

Let us start with the classical Ising model (h = 0)
and consider the two degenerate ground states that you can easily construct in this case: where |± = 1 √ 2 (1, ±1) T denote the two eigenstates ofσ x with eigenvalues ±1. Recall that the parity operators reads, in terms of spins, as P = L j=1σ z j , and thatσ z |± = |∓ . Hence, you easily deduce that: This implies that the two eigenstates of the parity operator must be: These two opposite parity states must be represented by the two fermionic ground states belonging to the ABC and PBC sectors. They are exactly degenerate for h = 0. The states in this doublet are crucial for the symmetry-breaking in the thermodynamic limit.
Now consider the effect of a small h, taking for simplicity of argument the Ising Hamiltonian with OBC: , +, · · · , + with l = 1 · · · L − 1 .
For h = 0, these L − 1 lowest-energy domain-wall excitations are degenerate, and separated from the two ground states by a gap 2J. Therefore, we can study the effect of a small transverse-field term, for |h| J, using standard textbook degenerate perturbation theory [2]. The Hamiltonian restricted to the L − 1-dimensional subspace of the domain-wall excitations has the form This Hamiltonian is quite easy to diagonalise, as it resembles a standard tight-binding problem with open boundary conditions. As in the quantum mechanical example of an infinite square well [2], it is simple to verify that the appropriate sine combination of two opposite momenta plane-waves of momentum k satisfy the correct boundary conditions: for n = 1, . . . , L−1, and N k a normalization factor. These delocalised domain-walls have an energy: We notice that if we expand the excitation energies k in Eq. (50) up to lowest order in h/J we obtain k ≈ µ k + O((h/J) 2 ). So, we see that a state with a single quasiparticleγ † k has the same energy as the delocalised domain-wall state in Eq. (81). This justifies the picture that quasiparticles are indeed delocalised domain walls. 13 Let us continue our perturbative reasoning to see how we can estimate the separation between the two ground states originating from the h = 0 doublet discussed above, when h > 0. As mentioned above, the states |+, +, · · · , + and |−, −, · · · , − are degenerate for h = 0, and this doublet is separated from the other states by a gap ≥ 2J. The degenerate states |+ = |+, +, · · · , + and |− = |−, −, · · · , − are coupled only at order L in perturbation theory: we need to flip L spins, with theσ z j operators, to couple one to the other. Hence, we expect their splitting to be ∆E 0 ∼ (h/J) L , i.e., exponentially small in the system size L for small |h|: the resulting eigenstates |ψ ± (h) , even and odd under parity, approach the two eigenstates in Eq. (77) for h → 0. 14 This energy splitting is exactly the quantity ∆E 0 discussed in Sec. 3.1. So, in the thermodynamic limit, we break the Z 2 symmetry. At any finite size we have the symmetry preserving ground states |ψ ± (h) which tend to Eq. (77) for h → 0. These states can be regarded as superpositions of two macroscopically ordered 15 . So, in the subspace generated by |ψ ± (h) there can be an explicit symmetry-breaking 16 of Z 2 only in the thermodynamic limit, where the two states are degenerate and the slightest local perturbation selects one of the two macroscopically ordered superpositions |± h .

Nambu formalism for the general disordered case
As we have seen, in the ordered case the Hamiltonian can be diagonalised by a Fourier transformation, reducing the problem to a collection of 2 × 2 "pseudo-spin-1/2" problems, followed by a Bogoliubov transformation, as first shown by Lieb, Schultz and Mattis [5,6]. In the disordered case, we can proceed in similarly, but we cannot reduce ourselves to 2 × 2 problems in a simple way. 17 By using the Nambu formalism, we define a column vector Ψ and its Hermitean conjugate row vector Ψ † , each of length 2L, by 13 We performed our analysis for the case of OBC. The case with PBC is similar, the only difference being that the lowest-energy excitations for h = 0 have two domain walls. Nevertheless, with an analysis very similar to the one above, one can show that the excited states for |h| J have energy 2µ k and can be interpreted as states with two quasiparticles. 14 We notice that the gap of order J between the ground-state doublet and the higher excited states is essential to get this result. The excited states at h = 0 appear in massively degenerate multiplets and therefore this construction is impossible for them. The doublet structure appears only for the two quasi-degenerate ground states and the two quasi-degenerate highest-energy states. 15 "Macroscopically ordered" means that the longitudinal magnetization operatorM x = jσ x j has an expectation value which is extensive in L. 16 Nevertheless, all the states in this doublet, |φ = α|ψ + (h) + β|ψ − (h) with |α| 2 + |β| 2 = 1, show long-range correlations also at finite size. Indeed, the correlator φ|σ x jσ x j+l |φ is always finite (equal to 1 in the limit h → 0), and lim expressing the long-range order associated with symmetry-breaking. See the related problem in Sec. 9. 17 For the time-independent case, a theorem due to Bloch and Messiah guarantees that there is always an appropriate basis in which the problem reduces to 2×2 blocks, but this is not very useful if you are willing to tackle dynamical problems. See Sec. 7.
Warning: Notice that the Ψ satisfy quite standard fermionic anti-commutation relations for j, j = 1, ..., 2L, except that { Ψ j , Ψ j+L } = 1 for all j ≤ L, which brings about certain factors 2 in the Heisenberg's equations of motion (see later).

!
It is useful, for later purposes, to introduce the 2L × 2L swap matrix S so defined: 19 Consider now a general fermionic quadratic form 20 expressed in terms of Ψ as: For a general quadratic fermion Hamiltonian, the 2L × 2L matrix H should be Hermitean and its L × L blocks A and B should be, respectively, Hermitean (A = A † ) and anti-symmetric There is an intrinsic particle-hole symmetry in a fermionic Hamiltonian having this form. This symmetry is connected with the fact that: In the Ising case all couplings are real and we have two different fermionic Hamiltonians, one for each parity sector p = 0, 1, which we report here for convenience: 21 with the boundary condition:ĉ (90) 18 The notation forĉ might be a bit confusing and should be intended as a shorthand, rather than a column vector. We are not consistently assuming, for instance, thatĉ † is a row vector. The same shorthanded but imperfect notations will be assumed later on for the Bogoliubov rotated operatorsγ andγ † . 19 Interestingly, you can write: 20 Indeed one can show that the most general quadratic form in the fermion operators where A j j = A * jj (A = A † is Hermitean) and B jj = −B j j (B = −B T is anti-symmetric, becauseĉ jĉj is anti-symmetric under exchange of the two operators, and any symmetric part of B would not contribute) has exactly the form given in Eq. (87), plus a constant term Tr A. 21 For convenience we also re-write The corresponding 2L×2L matrices H p are now real and symmetric. Hence A is real and symmetric (A = A * = A T ), and B is real and anti-symmetric (B = B * = −B T ): The structure of the two blocks A and B is given, in the Ising case, by: where we have assumed, once again, that J x j = J j (1+κ)/2 and J y j = J j (1−κ)/2. In the PBC-spin case, we have additional matrix elements: and both depending on the fermion parity p. The OBC case is recovered by simply setting J L = 0, which makes H 1 = H 0 .
5 Diagonalisation of H: the time-independent case.
We start considering the eigenvalue problem for a general Hermitean 2L × 2L matrix showing that intrinsic particle-hole symmetry of the problem leads to the Bogoliubov-de Gennes (BdG) equations.

The Bogoliubov-de Gennes equations.
Let us consider the eigenvalue problem for a general Hermitean 2L where u, v are L-dimensional vectors and µ index refers to µ-th eigenvector. By explicitly writing the previous system, we find the so-called Bogoliubov-de Gennes equations: In the Ising case, A = A * and B = B * , and we can always take the solutions to be real. 23 coincides exactly with Eq. (96), after taking a complex conjugation, exchanging the two equations and reshuffling the terms. An alternative derivation uses the fact that and that HS = −SH * , see Eq. (88). 23 Since H is a real and symmetric matrix, it can be diagonalised by a real orthogonal matrix.
We can organise the eigenvectors in a unitary (orthogonal, if the solutions are real) 2L × 2L matrix U and V being L × L matrices (real, as we can choose to be, in the Ising case) with the u j and v j as columns. As a consequence: If we define new Nambu fermion 24 operators Φ and Φ † in such way that we can write H in diagonal form Similarly to Ψ, we can define new fermion operatorsγ such that More explicitly, we can write: 25 24 We have: The conditions for the transformation in Eq. (104) to be canonical are: since you realise that the block 22 is simply the * of block 11 and block 12 is the † of block 21. Interestingly, the condition V T U + U T V = 0 tells us that V T U is anti-symmetric, and the same happens for U T V. Similarly, one must have: which again implies the anti-symmetry of both UV † and V * U T .
which can be easily inverted, remembering that Ψ = U Φ, to express theĉ j operators in terms of theγ µ : Finally H in terms of theγ operators reads, assuming we have taken µ > 0: and the ground state is the state annihilated by allγ µ , which we denote by |∅ γ : The 2 L Fock states can be expressed as: Warning: The previous discussion applies to a generic quadratic fermion Hamiltonian H. Consequently, it also applies to the two different parity Hamiltonians H p relevant for the Ising case, which one could express as: This implies that there are two distinct Bogoliubov vacuum states |∅ p , one for each set of operatorsγ p,µ . Recall, however, that the block Hamiltonian H p = P p H p P p involves projectors on the appropriate sub-sectors, which must be handled appropriately. Moreover, the possible presence of zero-energy eigenvalues must be appropriately taken care of: see below. This is important in calculating thermal averages, as further discussed in Sec. 8.

!
Before ending, one note about zero-energy eigenvalues, which has a practical relevance when calculating thermal averages. If you calculate the eigenvalues { µ } by a numerical diagonalisation routine, the presence of zero-energy eigenvalues complicates the story. Indeed, the zero-energy eigenvalues, if present, must come in an even number. This is rather clear from the fact that total dimension is 2L and that every non-zero positive eigenvalue µ > 0 must have a negative partner − µ < 0. Unfortunately, the computer will produce eigenvectors associated with the degenerate zero-energy eigenvalues which do not have the structure alluded at in Eq. (97). To enforce such a structure you can exploit the swap matrix S.
Info: Let us consider the Ising case, where H is real and particle-hole symmetry reads reads Hence the zero-energy subspace -the so-called Ker H, whose even dimension we denote by N 0 -is invariant for S. Hence you can restrict S to such N 0 -dimensional subspace, and diagonalise it there. But S is such that S 2 = 1, hence its eigenvalues can be only ±1, the eigenstates of S being even or odd under swap of the first and last L components. Even more: you can show that S must have exactly as many +1 as −1 eigenvalues in that subspace. Now, if (u , u) T is a zero-energy even-swap eigenstate, and (v , −v) T a zero-energy odd-swap eigenstate -both normalised and orthogonal -then the two combinations: are both normalised, orthogonal, and have precisely the structure shown in Eq. (97). These states should be used to enforce the required structure of Eq. (97), crucial to fulfilling the correct anti-commutation rules.
i Show that the static Bogoliubov-de Gennes equations in Eq. (96) are equivalent, for general couplings, to the diagonalization of the following tight-binding problem for the two-component whereτ are pseudo-spin Pauli matrices acting on the two components of W jµ .
Next, consider the uniform case J j = J and h j = h. Use Fourier transforms W jµ = 1 L k e ikj W kµ (where the k-vectors used depend, as usual, from the boundary conditions) to show that: where H k = (2κJ sin k)τ y + 2(h − J cos k)τ z as in Eq. (47). This shows that the correct correspondence between the general BdG approach of Sec. 5 and the k-space approach of Sec. 3, is given by 2 µ ≡ k .

Problem 1 Tight-binding formulation of the BdG equations.
Consider now a uniform Ising chain with κ = 1 and a single impurity on site j = l. More in detail, take Show that the impurity induces two bound-state excitations, one above and one below the continuum [2|J − h|, 2|J + h|] of the extended excitations of the uniform chain. Answer: Hint: Assuming that 2 µ is outside of the spectrum of the unperturbed uniform chain, and using Fourier transforms, show that you arrive at the following equation determining the non-trivial solutions µ: Assuming himp h, show that, for L → ∞, 1 2himp U jµ V jµ are localised in space. This means that these eigenfunctions are uniformly bounded by a function exponentially decaying over a characteristic length-scale ξ loc , the so-called localisation length. More formally, fixing h max and J min , there exists a ξ loc such that where l µ depends on µ and C is a constant. This phenomenon can be clearly seen when the system size exceeds the localisation length, L > ξ loc . Study localisation also using the inverse participation ratio [7,8] Average IPR µ over µ and verify that it tends towards a constant value, for increasing L. a a Observe that plane-wave delocalised states have IPR = 1/L, while fully localised states have IPR = 1.

Problem 3 Anderson localisation of states for the disordered Ising chain.
The space localisation phenomenon discussed in the previous problem is an example of Anderson localisation, see Refs. [7,8]. In the fermionic representation, the non-interacting quasiparticle excitationsγ † µ are space localised. In the spin representation, this is equivalent to saying that the non-interacting dressed domain-wall excitations are localised. This can be seen using a heuristic argument which amounts to applying the Jordan-Wigner transformation, Eq. (11), to the definition ofγ † µ , Eq. (104). One obtainsγ † Let us consider the case h max J min , where the ground-state doublet is well approximated by the states in Eq. (77), see also discussion below. By applying the operator in Eq. (113) to the ground state,K j essentially acts only on the site l µ where U jµ and V jµ are localised. Its action consists of flipping all the spins for j < l µ thereby creating a domain wall.
The fact that excitations are localised domain walls has an important physical consequence: the excitations cannot destroy the long-range order of the ground state in the ferromagnetic phase. They can alter the state only locally, and long-range order is preserved. Indeed, in the thermodynamic limit the correlator lim l→∞ | σ x jσ x j+l | is always finite. 26 In order to see that for h max J min all the spectrum shows long-range order, we proceed by means of a perturbative argument similar to that of Sec where · · · denotes the average over the distribution of the random J j and h j , and the corrections of order O( √ L) result from the central limit theorem (assuming the h j and J j independently distributed). 27 So, in this disordered ferromagnetic case, the whole spectrum behaves as the ground state of the clean model discussed in Sec. 3.2. All the eigenstates are linear superpositions 28 of macroscopically ordered states, and break the Z 2 symmetry in the thermodynamic limit. In particular, any eigenstate shows long-range order, since | ψ ± |σ x jσ x j+l |ψ ± | = 1 for any value of j and l, as can be easily verified. The present analysis is valid in the limit h max J min , but long-range order of excited eigenstates was found also beyond this limit [10].

Open boundary conditions and Majorana fermions
The case of a chain with open boundary conditions is particularly interesting, because Majorana fermions, and the associated zero-energy modes, emerge quite naturally from the discussion [11]. In this section, we work out explicitly the case of a chain with open boundary conditions, introduce the Majorana fermions first as a formal device to perform the diagonalization and then discuss the physical role they have as boundary excitations at vanishing energy in the broken symmetry phase.
For illustration purposes, let us consider the case in which the spin chain has L = 4 sites, and couplings J 1 , J 2 , J 3 > 0, while J 4 = 0 (as dictated by OBC). 29 The 2L × 2L Hamiltonian matrix 26 See Sec. 9 for the evaluation of this correlator. 27 We can naively expect that the phase with all the spectrum breaking the Z 2 symmetry in the thermodynamic limit is destroyed when ∆ ∼ 1, that is to say when log h = log J . Remarkably, the result of this naive argument is confirmed (at least for the ground state) by a very accurate renormalization-group analysis [9]. 28 The so-called "cat states". 29 If you want to do numerical tests, we suggest you take different values for J j , for instance, J 1 = 1, J 2 = 2 and J 3 = 3, to avoid degeneracies, which might lead to a mixing of the corresponding eigenvectors.
(in this case, an 8 × 8 matrix) will have the form (we fix the anisotropy parameter κ = 1): Let us consider first the case with h = 0, corresponding to the classical Ising model with the given couplings. The corresponding eigenvalues/eigenvectors (disregarding the ordering of the non-zero eigenvalues) are found to be: where you should observe that the structure of Eq. (97) is correctly respected, except for the two zero eigenvalues, which the diagonalisation routine has decided to give us in this particular form. This form itself is particularly interesting. It suggests that the following fermionic combinations naturally emerge: Several things strike our attention. First: Φ 4 is Hermitean and Φ 8 is anti-Hermitean, and they are not Hermitean conjugate pairs, contrary to all other ( Φ j , Φ j+4 ) pairs. If you want to construct ordinary fermionic operators, then you should redefine: with an orthogonal transformation which leaves the subspace of two degenerate eigenvalues 0 invariant, precisely as alluded at in the last info-box of 5.1. Second: certain Hermitean combinations seem to play a peculiar role. In particular, let us define the Majorana fermions: 30 These operators are manifestly Hermitean. They allow to express the original fermions as: and satisfy the anti-commutation relations: Notice, in particular, that this implies that different Majorana anti-commute, but (č j,α ) 2 = 1.
In terms of these operators we have: There is something simple behind the previous story. If you rewrite the original Ising coupling in terms of fermions, you realise that, for instance: i.e., the Ising term couples in a precise way neighbouring Majorana operators. All we have done, to diagonalise it, is to introduce the appropriate combinationγ j = 1 2 (č j+1,1 − ič j,2 ) andγ † j = 1 2 (č j+1,1 + ič j,2 ) and re-express the coupling term as: which suggests that the ground state is the vacuum of thoseγ j operators. 31 There is a second simple case we can deal with. Take all h j > 0 and J j = 0, so that the Hamiltonian is now: This shows that the ground state, now the vacuum of theĉ j , still involves a "pairing" of Majorana fermions, but now on the same site j. The two different Majorana pairings are sketched in Fig. 6.
31 Interestingly, in the vacuum we gain an energy −J j from each bond. Breaking that, we would get a contribution +J j from the bond, hence an energy cost, referred to the vacuum, of 2J j : this explains, for instance, the factor 2 in front of µ in Eq. (106). Returning to the previous case with h j = 0, the ground states certainly verifŷ But there are two states satisfying such a condition, a degeneracy that is ultimately related to the presence of unpaired Majorana operators at the end of the chain, as emerging from Fig. 6 (left). Indeed, one such state is also the vacuum ofγ L=4 : On such state, we have (generalising now to arbitrary even L) The second possible ground state is |∅ 1 =γ † L |∅ 0 for which: These two ground states have opposite fermion parity, 32 because they differ by the application of γ † L .
32 Notice, incidentally, that the fermionic parity can be expressed as: Info: As discussed by Kitaev [11], the two zero-modes survive for 0 < h < J, with a splitting which is exponentially small in the length of the chain, as long the broken symmetry leads to two possible ground states.
The existence of these modes is deeply related to the topological considerations done when discussing Fig. 4. Indeed, a fermionic chain with |h| < J and OBC is equivalent to surrounding the chain with the fermionic vacuum -in turn equivalent to an Ising chain with h → ∞. But one cannot go continuously from a phase with winding index 1 to a phase with index 0. Therefore, at the border between two phases with different index, the gap must close to enforce this discontinuity (we saw when discussing Fig. 4, the deep connection between the discontinuity of the index and the closing of the gap). Hence, the gap must close at the boundary, and this effect appears as two zero-energy boundary modes which behave as Majorana excitations. As we saw above, there are only two ways of combining them into fermionic excitations, which are indeed very non-local objects. For any finite system size L, the two Majorana fermions have an overlap exponentially small in L. If we combine them into fermionic excitations, we find a gap between them which is exponentially small in the system size. This is the same gap we found in Secs. 3.1 and 3.2, that vanishes in the thermodynamic limit and leads to symmetry-breaking. Now we appreciate its intimate connection with topology.
To visualise these facts, we show in Fig. 7 the spectrum of eigenvalues µ ≥ 0 (evaluated numerically) of an ordered Ising chain with OBC. We mark in red one of the two zero-energy modes we have discussed above, surviving for all h ≤ h c = J. Actually, the mode is not exactly at zero, but, rather, exponentially small in L. Finite-size effects (here L = 256) lead to a visible rounding effect in the proximity of h c .

The BCS-form of the ground state.
The next problem we would like to solve is how to write the Bogoliubov vacuum |∅ γ in terms of theĉ † j in the general non-homogeneous case, in a way that generalizes the simple BCS form we have in k-space: For that purpose, let us make the Ansatz that |∅ γ can be written as a Gaussian state of the form: where Z will be our shorthand notation for the quadratic fermion form we exponentiate. Clearly, we can take the matrix Z to be antisymmetric (but complex, in general): any symmetric part of Z would give 0 contribution. The conditions that Z has to satisfy should be inferred from the fact that we pretend thatγ µ |∅ γ = 0, hence: Since Z is made of pairs ofĉ † s, it commutes withĉ † j , hence,ĉ † j e Z |0 = e Zĉ † j |0 . The first term, containingĉ j e Z |0 , is more problematic. We would like to commuteĉ j through e Z to bring it towards the |0 , where it annihilates. To do so, let us start calculating: where we have used the antisymmetry of Z. We see, therefore, that [ĉ j , Z], being a combination of c † j commutes with Z and with any function of Z. It takes then little algebra 33 to show that: The conditions in Eq. (134) therefore read: Noticing thatĉ j |0 = 0, substituting Eq. (135), and omitting irrelevant prefactors we therefore have: where we have exchanged the dummy indices j and j in the first term. Next, we collect the two terms by writing: 33 Simply expand the exponential in the usual way, realise that This is the condition that Z has to verify for the state |∅ γ to be annihilated by allγ µ . This is the so-called Thouless formula [12]. It takes very little algebra 34 to verify that, indeed, such a form of Z is antisymmetric.
According to a theorem of linear algebra, any antisymmetric matrix can always be reduced to a "standard canonical" form by applying a unitary matrix D as follows: [12] where in general the λ p are complex. If L is even, there are L 2 blocks 2×2 with some λ p , while if L is odd, Λ has an extra row/column of zeroes. The unitary matrix D allows us to define combinations of the fermions c † j which form natural "BCS-paired" orbitals, Labelling the consecutive columns of D as 1, 1, 2, 2, · · · , p, p, · · · , with p up to L/2, one can readily check that in terms of the d † s the Bogoliubov vacuum |∅ γ reads: It remains to evaluate the normalisation constant N. Now we calculate: 35 Summarising, we have derived the so-called Onishi formula [12], which states that: 34 Observe that: However, from block 12 in Eq. (102) we get: 35 In the derivation we use that:

32
6 Dynamics in the time-dependent case If we express the Bogoliubov vacuum in terms of the λ p we have: where we have defined u p = 1/ 1 + |λ p | 2 and v p = λ p / 1 + |λ p | 2 , which verify |u p | 2 + |v p | 2 = 1.

Dynamics in the time-dependent case
A time dependence can come from many different sources. The simplest case, which is used in the so-called quantum annealing approach, consists in assuming that the transverse fields are time-dependent h j (t), for instance they might be slowly annealed from a very large value towards zero, or changed in some time-periodic fashion. In this way, the diagonal elements of matrix A become time-dependent and consequently H → H(t). 36 We proceed now in general, assuming A(t) and B(t).
Start from Schrödinger's equation: Since the norm of |ψ(t) must be conserved, this implies the existence of a unitary evolution operator U (t, t 0 ) such that |ψ(t) = U (t, t 0 )|ψ(t 0 ) , which satisfies the same equation: Next, consider the expectation value of a time-dependent operatorÔ(t) in the Schrödinger's picture where we have introduced the Heisenberg's picture operator Therefore the equation of motion of an operator in Heisenberg's picture for the general case of a time-dependent Hamiltonian reads: 37 36 Alternatively, a time-dependence is found when doing time-dependent mean-field approaches. 37 Here we use: By calculating the commutator we see that we have a linear equation of motion and analogously for the operatorĉ † j . With a more compact notation, one can write the linear Heisenberg equations of motion for the elementary fermionic operators as: the factor 2 on the right-hand side originating from the off-diagonal { Ψ j , Ψ j+L } = 1 for j = 1 · · · L. The initial condition for these equations can be written as: whereγ are the Bogoliubov fermions that diagonalise H(t 0 ), and U 0 the corresponding rotation matrix.
We are not quite done: We have an explicit linear equation for Ψ H (t), but we need an explicit solution for this equation, obtained by some "simple enough" integration of a finite-dimensional linear problem. There are now at least two ways of getting the desired result.
First route. We make the Ansatz that |ψ(t) , the time-evolved state of the system, is a Bogoliubov vacuum annihilated by a set of time-dependent quasi-particle annihilation operatorsγ µ (t) This requirement immediately implies, by taking a total time-derivative, that: where we have added, in the last step, a termγ µ (t) |ψ(t) = 0. The last expression implies: 38 By considering the equation of motion of the Heisenberg operatorγ µH (t) we have where we have used Eq. (160) in the last step. So, sinceγ µH does not depend on t, it must coincide with its t = t 0 value. Let's call this valueγ µ =γ µH =γ µ (t = t 0 ).
Let us assume now, inspired by Eq. (105), that theĉ jH (t) are indeed expressed bŷ and let us see if this expression solves the required Heisenberg equations in Eq. (155) for an appropriate choice of the time-dependent coefficients U jµ (t) and V jµ (t). Substituting in Eq. (155) we get: By equating the coefficients ofγ µ andγ † µ we obtain the time-dependent Bogoliubov-De Gennes equations: or more compactly, collecting together µ = 1, · · · , L solutions in L × L blocks U and V: 39 T is also a solution, so we need to find only µ = 1, · · · , L solutions, as indeed alluded by the compact form (166), not 2L. Once we have the first L, it is automatically guaranteed that: 38 A mathematician would complain, here, that this is not a valid implication: an arbitrary linear combination of γ µ (t) could be added that, acting on |ψ(t) , gives 0. We are a bit swift here, but the result is correct. We will get to the same result by a second route in a short while. 39 In the time-independent case, the solution is equivalent to solving the time-independent Bogoliubov-De Gennes equations. Indeed in this case the time evolution of the solution is and, as you can easily verify, the same result can be obtained by using directly Eq. (166) with H(t) = H. i d dt Second route. It is reassuring to get to the same time-dependent Bogoliubov-de Gennes equations by a second, quicker, route. Let us recall the linear equation we want to solve, with its initial condition: whereγ are the Bogoliubov fermions that diagonalise H(t 0 ), and U 0 the corresponding 2L × 2L rotation matrix. Inspired by the form of the initial condition, let us search for a solution of the same form: with the same Φ used to diagonalise the initial t = t 0 problem. In order for this to be a solution, the time-dependent coefficients U(t) must satisfy the linear Bogoliubov-de Gennes time-dependent equations: with initial conditions U(t = t 0 ) = U 0 . The latter form is just a compact way of expressing Eq. (167).
It is easy to verify that this implies that the operatorsγ µ (t) in the Schrödinger picture are timedependent and annihilate the state |ψ(t) : this was indeed the starting point of the Bogoliubov Ansatz presented in the first route. Indeed, since we can immediately write, in the Schrödinger picture: If we go back to Sec. 5.3, we realise that the algebra carried out there is perfectly applicable here, and allows us to write the time-dependent state |ψ(t) in the explicit Gaussian form: with the anti-symmetric matrix Z(t) given by: It is not very hard to explicitly verify that such a state satisfies the Schrödinger equation: provided U(t) and V(t) satisfy the time-dependent BdG equations in Eq. (166). Indeed, the time derivative of the state |ψ(t) is simply: On the right-hand side, the Hamiltonian terms can be rewritten by using that, for instance: Rewriting similarly all the Hamiltonian terms we get: By explicitly calculating the derivative of Z(t) using the BdG equations one can check, after some lengthy algebra, that the two expressions indeed coincide. 40 Consider the uniform model with anisotropy κ = 1. Take J = 1 (constant in time, and taken as our unit of energy), and consider a time-dependent transverse magnetic field h(t). Show that, in analogy with the form of the ABC ground state in Eq. (65), the time-dependent state solves the time-dependent Schrödinger equation i |ψ(t) = H 0 (t)|ψ(t) in the fermionic ABC sector provided u k (t) and v k (t) satisfy, for all k, the following BdG equations: Problem 4 Time-dependent BdG equations for a uniform chain.
Consider now a slow annealing of the transverse field h(t) from the initial value h i J at time t = 0, to the final value h f = 0 at time t = τ , for instance linearly: h(t) = h i (1−t/τ ). Initialise the system in the ground state of H k (t = 0), and numerically study the BdG evolution for all k, for a sufficiently large L. Consider now the expectation value of the density of defects over the ferromagnetic ground state at the end of the non-equilibrium protocol Show that ρ def (τ ) ∼ τ −1/2 for sufficiently large τ , provided L is large enough. a This is the so-called Kibble-Zurek scaling, see Ref. [13] for further references on this topic.
a The power-law scaling of ρ def (τ ) holds for τ τ * L ∝ L 2 . For larger τ , the finite-size critical gap at hc/J = 1 becomes relevant, and the density of defects starts decaying faster.
Problem 5 Crossing the critical point with an out-of-equilibrium protocol.

Calculating time-dependent expectation values.
Once we have a solution to the time-dependent BdG equations, we can calculate time-dependent expectations of operators quite easily. Consider, for instance, the elementary one-body Green's function: We assume that the initial state |ψ(t 0 ) is the Bogoliubov vacuum of the operators γ which diagonalise H(t 0 ), i.e., |ψ(t 0 ) = |∅ γ . The algebra is most directly carried out by working with the 2L × 2L Nambu one-body Green's function matrix by using the fact that the corresponding transformed Green's function is simple, since |ψ(t 0 ) = |∅ γ : In matrix form, we immediately calculate: Summarising, the four L × L blocks of G read: Info: Notice that, quite generally, G is Hermitean, while F, as a consequence of the fermionic anti-commutations, see also Eq. (103), is anti-symmetric: and i Expectations of more complicated operators can be reduced to sums of products of Green's functions through the application of Wick's theorem [1]. Moreover, time-correlation functions with Heisenberg operators at different times can be calculated similarly.

Floquet time-dependent case.
A particular case of dynamics is that in which the Hamiltonian is periodic in time, i.e., a period τ exists such that H(t + τ ) = H(t). The Floquet theorem [14,15] guarantees the existence in the Hilbert space of a complete basis of solutions of the time-dependent Schrödinger equation which are periodic "up to a phase factor", i.e., such that: This way of writing is closely reminiscent of the time-independent case, except that the state |ψ Pα (t) , known as Floquet mode, is now periodic in time rather than a time-independent eigenstate of the Hamiltonian; the E α , which plays the role of the eigenenergy, is known as Floquet quasienergy. There are 2 L , as many as the dimension of the Hilbert space, Floquet solutions of this type, and these solutions can be used as a convenient time-dependent basis to expand states. Their usefulness consists in the fact that if we expand a general initial state as then the time-evolution can be written, for free, in a form that is reminiscent of the timeindependent case, i.e.: An explicit construction of the many-body Floquet states can be obtained through a Floquet analysis of the time-dependent Bogoljubov-de Gennes (BdG) equations, in a way similar to that used to construct the energy eigenstates from the solution of the static BdG equations (see Sec. 5). To do that, let us write the BdG equations (169) Since H(t + τ ) = H(t) is a periodic 2L × 2L matrix, the Floquet theorem guarantees the existence of a complete set of 2L solutions which are periodic up to a phase. L of them have the form: and the remaining L, by particle-hole symmetry, are automatically obtained as .
Collecting all the quasi-energies µ into a diagonal matrix = diag( µ ), and the various column vectors u Pµ (t) and v Pµ (t) into a L × L matrices U P (t) and V P (t), it is straightforward to show that the structure of the Floquet solutions of the BdG solutions is 41 Using these solutions, we can construct the Bogoljubov operatorsγ Fµ (t) which annihilate a vacuum Floquet state |∅ F (t) through the standard method employed in the general time-dependent case (see Eq. 104): or, more explicitly, for µ = 1, · · · , L: The Floquet vacuum state |∅ F (t) annihilated by all theγ Fµ (t) has the Gaussian form (see Eq. (172)): where, see Secs. 5 and 6, the Thouless and Onishi formulas hold: Let us show that the Floquet vacuum state is periodic, i.e., or, to put it differently, its many-body quasi-energy is E 0 = 0. To this aim, it suffices to show that Z F (t) and N F (t) are both periodic.
From these relationships, in turn, it follows immediately that the quasi-energy phase-factors cancel in Z F , i.e.: which is manifestly periodic in time, Z F (t + τ ) = Z F (t), because both U P and V P are periodic. The periodicity of N F (t) follows because det(U F (t)) = det(U P (t)) det(e −i t/ ) = det(U P (t)) e −i µ µt/ = det(U P (t)) , i.e., once again something manifestly periodic in time. At this point, we can easily, in principle, construct all the 2 L many-body Floquet states by simply applying any product ofγ † Fµ (t) to |∅ F (t) : where n µ = 0 or 1 is the occupation number of theγ † Fµ (t) operator. From Eq. (186) and the periodicity of the Floquet vacuum, it follows that the quasi-energy of |ψ F{nµ} (t) is given by:

Overlap between BCS states
Sometimes, for instance in the context of quantum quenches, where the Hamiltonian is abruptly changed, it is important to know how to calculate the overlap between BCS states belonging to two different Ising Hamiltonians H 0 and H 1 . Let us start considering the two BCS ground states of H 0 and H 1 . These two states are Bogoliubov vacua with respect to the fermionic operatorsγ 0µ andγ 1µ , and we denote them, for a more compact notation, as |∅ γ0 = |∅ 0 and |∅ γ1 = |∅ 1 We will first compute | ∅ 1 |∅ 0 | 2 , and then we will extend the result to the overlap of general excited states. The two sets of fermions can be written in terms of the original Jordan-Wigner fermions as: We can write the direct unitary transformation from one set to the other as follow: where: with: We will prove that, if |∅ 0 and |∅ 1 are not orthogonal, then: a relationship which is known as Onishi formula. Indeed, we have already given a proof of this relationship in Sec. 5.3, for the special case in which one of the two sets of fermions where the original Jordan-Wigner fermionsĉ j with associated vacuum state |0 . There we showed that, with the present notation: with: and | 0|∅ 0(1) | 2 = |N 0(1) | 2 = | det(U 0(1) )|. With exactly the same algebra, we could establish, for instance, that: with: and | ∅ 0 |∅ 1 | 2 = |N| 2 = | det(U)|. But there are several points where we were a bit swift, in this derivation: for instance, we assumed that U is invertible, which is not guaranteed. The issue of possible orthogonality was also not raised. Moreover, the case in which the resulting ground state is a pure Slater determinant without BCS-pairing, a case which is relevant for the XY model, is not discussed in this way.
We will now give an alternative proof which makes use of an interesting theorem due to Bloch and Messiah [12,16], and which clarifies all these issues. We will perform an intermediate canonical transformation which first allows us to write an explicit equation for |∅ 1 in terms of |∅ 0 , and then to compute easily ∅ 1 | ∅ 0 . The theorem that Bloch and Messiah proved [16] shows that matrices with the structure of U above can be decomposed into a product of three unitary transformations as follows: where D, C are L × L unitary matrices and U, V are L × L real matrices of the form: in which u p > 0, v p > 0 and u 2 p + v 2 p = 1. From these relations we have: The idea is the following. Since: we can think of the transformation as the product of (1) a first unitary transformation C which does not mix particles and holes for fermionsγ 1 , defined by followed by (2) a simple "canonical form" of a transformation leading to new fermions: The final transformation (3) leading to the fermionsγ 0 is again a unitary D which does not mix particles and holes: In essence, what the theorem guarantees is that one can always find a basis such that the transformed fermions,α 0 andα 1 , are coupled by a particularly simple matrix in which there are only three possibilities: (i) for some indices, which we denote by l, there is no transformation at all (the ones in the diagonal of U), i.e.,α 1l =α 0l ; (ii) for some other indices, which we denote by k, the transformation is a pure particle-holeα † 1k =α 0k : these indices correspond to the zeroes in the diagonal of U, and the ones in the diagonal of V; (iii) all other indices, denoted by (p, p), are BCS-paired in a simple way, and they form 2 × 2 blocks in the matrices U and V with coefficients u p and v p , such that: We must stress that the theorem does not tell us how many indices belong to the three categories above: in some cases, all the indices might be 2 × 2-paired, but it is also possible that the transformation is a pure particle-hole transformation without any pairing at all.
The construction of the relationship between |∅ 0 and |∅ 1 becomes particularly simple in terms for the fermionsα 0(1) . The key idea is theα 0(1) is related toγ 0(1) by a transformation which does not mix particles and holes, and therefore it is still true thatα 0n |∅ 0 = 0 andα 1n |∅ 1 = 0. Since |∅ 1 is the state which is annihilated by anyα 1n we can write it as: where N is a normalization constant. Notice that we included only BCS-paired indices and particlehole transformed k-indices but not l-indices, sinceα 1l =α 0l and the inclusion of such terms would give zero, sinceα 0l |∅ 0 = 0. Since, by hypothesis, the two states |∅ 0 and |∅ 1 are not orthogonal there should not be pure particles-holes k-indices either, and therefore: Finally, since U = D † UC † , and D, and C are unitary transformations: which is what we wanted to show.
The extension to the calculation of the overlap between |∅ 0 and any eigenstate |{n 1µ } = µ∈Iγ † 1µ |∅ 1 , where I is the set of occupied states (n 1µ = 1), is in principle straightforward. Here is a possible way to tackle the problem. This state can be thought as the vacuum of the following new set of fermions: in which we have performed a particle-hole transformation for the occupied modes. Now we can use the equation obtained for the scalar product between empty states, i.e., where the matrix U is: in which: A second approach to calculate these overlaps with excited states makes explicit use of the Gaussian nature of the states. The relevant algebra follows directly from that of Sec. 5.3. Let us start by considering the overlap bewteenγ † 0µ1γ † 0µ2 |∅ 0 and |∅ 1 = Ne Z |∅ 0 . This is given by: where in the second step we have made use of the commutation property: Notice that, in order for the overlap to be non-vanishing, we were forced to contractγ 0µ2 against γ † 0µ 1 in the final step. A similar calculation shows that, if we have an even number 2n of operator, the result is highly reminiscent of a Wick's theorem sum-of-products of contractions: while the overlap vanishes for an odd number ofγ 0µi . In the last expression, the Wick's sum is rewritten in terms of the so-called Pfaffian of the anti-symmetric matrix Z (or more properly, of the 2n × 2n elements of Z required by the indices µ 1 · · · µ 2n ): Notice that the Pfaffian is defined by a Wick's sum which contains n products of Z-matrix elements, and not 2n, as the familiar det (Z) 2n×2n . However, a remarkable identity exists [17] which links the two objects: Notice, however, that the link exists only if the dimension of the antisymmetric matrix we are considering is even: The determinant of an odd-dimension anti-symmetric matrix is simply zero, while the Pfaffian is not defined. Summarising, we have obtained the generalization of the Onishi formula in the form:

Thermal averages
Let us recall a few basic facts about the general structure of the Ising model Hamiltonian. The full Hamiltonian, when PBC are imposed to the spins reads: The two blocks of even and odd parity can be written as: where H p=0,1 both conserve the fermionic parity, hence they commute with P 0,1 , and are given by: with the boundary condition set by the requirementĉ L+1 ≡ (−1) p+1ĉ 1 .
Warning: Neither H p=0 nor H p=1 , alone, expresses the correct fermionic form of the PBCspin Hamiltonian. Indeed, after the BdG diagonalisation, we can write: hence we have, in general, two different vacuum states |∅ p such thatγ p,µ |∅ p = 0, and 2 L corresponding Fock states, with different eigenvalues, unless the spins have OBC, in which case J L = 0 and therefore H p=1 = H p=0 . But we should keep only 2 L−1 even and 2 L−1 odd eigenvalues! How to do correctly the sums involved in the thermal averages is the next issue we are going to consider.
This derivation uses standard properties of projectors and of the trace, and, in the final step, the assumption thatÔ commutes with the parity. The thing to remark is that now the fermionic Hamiltonians H p appear, accompanied by a single projector P p . Next we recall that Hence we arrive at: Tr Ô e −β H = 1 2 p=0,1 Tr Ô e −β Hp + (−1) p Tr Ô e iπ N e −β Hp .
The next thing to consider is how to deal with the term e iπ N , which we can always re-express as: Info: The meaning of such expression should be reasonably transparent. Parity is a good quantum number. Once you determine it on the Bogoliubov vacuum, calculating ∅ p |e iπ N |∅ p = ±1 -see Eq. (249) in Sec. 9 to learn how to calculate it -, then the parity of each Fock state simply amounts to counting the number ofγ † operators applied to the Bogoliubov vacuum. There is a slight ambiguity in the meaning of |∅ p that is good to clarify here. We have not defined |∅ p to be the ground state in the sub-sector with parity p -in which case you would directly anticipate that ∅ p |e iπ N |∅ p = (−1) p , but rather the Bogoliubov vacuum state associated with H p . So, depending on the couplings and boundary conditions, the parity of |∅ p might differ from (−1) p . There are cases where, for instance, there is a single H with a single associated |∅ , but also cases where such a single H can have two degenerate vacuum states |∅ p , as discussed in Sec. 5.2. Hence, defining: we finally get: In particular, the partition function can be expressed as: 43 Observe that to calculate the fermionic single-particle Green's functions we need the separate ingredients for the two sub-sectors p: where for every p-sector we defined (see Sec. 5) We have tested these formulas on a chain of length L = 10. In the left panel of Fig. 8 we show that for L = 10 the Jordan-Wigner results are equal to those obtained by Exact Diagonalisation (ED) of the problem with both open and periodic boundary conditions. In the right panel of Fig. 8 we compare, for L = 100, the thermal averages obtained with the Jordan-Wigner approach for OBC and PBC.

Spin-spin correlation functions
The question we address here is how to calculate spin-spin correlation functions for the quantum Ising chain, see Ref. [18] for a more in-depth and general discussion. Here we discuss in detail the equilibrium case.
Certain correlation functions, which do not involve the Jordan-Wigner non-local string, come automatically from using the Wick's theorem, once you have calculated the single-particle Green's functions, including the anomalous term, see Sec. 6.2.
But suppose you want to calculateσ x j1σ x j2 correlations 44 for j 2 > j 1 : Recall now that: where the last expression involves the definitions [19] A j = (ĉ † j +ĉ j ) ≡č j,1 and closely related to the Majorana fermions. Hence we can write: 44 Recall that, in our convention, x is the spin direction in which the symmetry-breaking occurs.
where we used that ( A j ) 2 = 1.
At this point we use Wick's theorem [1], since this is a product of fermion operators averaged on the ground state of a quadratic Hamiltonian (a Gaussian state, as we have seen), or on a thermal state of a quadratic Hamiltonian, depending on whether we are calculating correlations at T = 0 or finite T . There are 2(j 2 − j 1 ) elements in the product, which we should be contracted pairwise in all possible ways, with the appropriate permutation sign [1].

Recall:
The results of Sec. 6.2 concerning the elementary fermionic Green's functions tell us that: and Recall also that the anti-commutation of the fermionic destruction operators forces F to be anti-symmetric. Moreover, if the average is taken over the ground state, so that U and V can be taken to be both real a , then G = UU T is real and symmetric, and F = UV T is real and anti-symmetric.
a The same results can be easily obtained from thermal averages, see Sec. 8.

!
Let us start by observing that where we used that G = G T and (F * ) T = −F. Similarly, we have that: These findings simply eliminate contractions between operators of the same type. 45 We are therefore left with the contractions of the type: and as perhaps expected.
Warning: The previous considerations apply to expectations calculated on the ground state or to thermal expectations. When considering time-dependent expectations the condition of reality of the matrices G and F no longer applies, and appropriate modifications would be needed, see Ref. [18].
! 45 Observe that contractions of the type A j A j or B j B j never occur from the Wick's expansion.
We now have to account for all possible contractions of the B-A type, say, with the proper permutation sign. If you think for a while, you realise that you can organise the whole Wick's sum into the determinant of an appropriate matrix as follows: Here, to help the reader recognising the various contractions, we have used colours.
Info: A good way to understand the structure of the matrix determinant you see is to notice that the second (column) index is constant -going from j 1 + 1 to j 2 -and tells you which is the A operator in the contraction: the corresponding first (row) index tells you the B partner in the contraction, and as you see it grows from j 1 up to j 2 − 1, as appropriate for the B partners.
i As a simple consequence of the previous algebra, let us see how to calculate the ground-state fermionic parity operator we used in the evaluation of thermal averages (see Eqs. (228) and (229)). We find: Consider a uniform Ising chain with PBC. By working in the fermionic ABC sector, calculate numerically the spin-spin correlation function C xx 1,j for the three cases: a) h/J = 1/2, b) h/J = 2 and c) h/J = 1. Verify that the correlations tend to decrease with increasing j, until j ∼ L/2, and then increase back towards a value C xx 1L ≡ C xx 12 . Explain why this happens. Verify numerically that, when L → ∞ -practically, try increasing L in your calculation -, C xx 1L/2 tends towards a finite value for case a), it goes to zero exponentially fast in case b), and as a power-law in case c). Estimate the exponent of such a power-law. a a Compare with the analytical solution provided in Ref. [6].

Entanglement entropy
The question we address here is how to calculate the entanglement entropy for the quantum Ising chain [20]. To explain what it is, let us start from the concept of reduced density matrix.
For a system described by a quantum state |ψ(t) , we can equivalently adopt a density matrix formulation [21], usingρ(t) = |ψ(t) ψ(t)|. The reduced density matrix is a proper (generally nonpure) density matrix obtained by tracing out a part of the system. To be concrete, if our Ising chain has L sites, then we can write: whereρ {L} (t) denotes the pure-state density matrix, and Tr {L}\{l} indicates that we take a partial trace over the sites of the chain not belonging to the set we have indicated by {l}, for instance the first l sites of the chain, {l} = 1, 2, · · · l.
The reduced density matrixρ {L} (t) must be a positive Hermitean operator acting in the Hilbert space of the {l} spins, whose trace is 1. It has, for a given value of t, 2 l non-negative eigenvalues w i (t) which sum to 1. The only case in which it is itself a pure state is when only 1 of them is equal to 1, w 1 = 1, and all the other vanish. This, in turn, is only possible when the state of the two subchains {l} and {L} \ {l} is a product state, i.e., a state without entanglement. Hence, a good way to capture such entanglement is to calculate the so-called entanglement entropy which vanishes exactly when the state is a product state, and is positive otherwise.
The first question is how the coefficients C α1···α l (t) in this expansion can be written.
Info: Let us recall the analogous problem for the case of a single spin. If we writeρ = 1 2 α C ασ α , then while C 0 = 1 due to the unit trace ofρ, we can reconstruct the other coefficients C α=x,y,z from measurements of the different components of the spin, more precisely, from the spin expectation values σ α = Tr σ αρ = 1 2 where we used that Tr(σ ασα ) = 2δ α,α . The name tomography is often used in this context: you reconstruct full information on the state by appropriate measurements.
Remark: The C α1···α l (t) are 4 l complex coefficients which completely specify the reduced density matrix, an operator in a 2 l -dimensional space. How to get the eigenvalues of such a reduced density matrix, and extract the entanglement entropy, seems a highly non-trivial task, at this stage. Notice also that the previous considerations apply to any spin system, even those that cannot be solved. There is nothing special, so far, about the quantum Ising chain.
! Let us look more closely to an example of such a term. Consider, for a block of l = 8 sites the term C 000xz00y : First of all observe that such a term respects the parity invariance of the Hamiltonian, as it possesses an even number ofσ x andσ y operators: if there was an odd number of them, the coefficient would vanish due to the parity selection rule. Next, we map spins into fermions.
(259) Warning: In the notation of Sec. 9, see Eq. (240), such an expectation value would translate into: hence in equilibrium (i.e., for a ground state or thermal calculation) it would still vanish, simply because you cannot construct the correct number of non-vanishing Wick's contractions! In any case, you notice how the approach we have undertaken is essentially impossible to carry out. Even after calculating, through an appropriate application of Wick's theorem, all possible non-vanishing coefficients C α1···α l calculating the eigenvalues of the corresponding reduced density matrix seems an incredibly difficult task! ! 46 Use the fact that (1 − 2n j ) commutes with terms which do not involve fermions at site j, and that (1 − 2n j ) 2 = 1. This leads, in particular, to a cancellation of the two tails originating from the Jordan-Wigner string operatorŝ K 4 andK 8 . Moreover, recall that the square of a Majorana gives the identity: (č j,1 ) 2 = (č j,2 ) 2 = 1.
But there is something very special about a quantum Ising chain, encoded in the fact that its state |ψ(t) is a Gaussian state, and Wick's theorem allows calculating any expectation value of fermionic operators in terms of elementary one-particle Green's functions. It turns out that working with Majorana fermions is a good way of handling efficiently ordinary and anomalous fermionic Green's functions. To do that, let us be equipped with a matrix notation for the Majorana as well.
where we defined the 2L × 2L block matrix: One can define the Majorana 2L × 2L correlations matrix as: which in full matrix form is immediately related to the Nambu Green's function G(t): Upon substituting the block-expression for the Nambu Green's function G(t) in Eq. (179) we obtain, after simple block-matrix algebra: where the 2L × 2L matrix A(t) is real and anti-symmetric, 48 and both G(t) and F(t) are timedependent and (in general) complex.
Remark: Notice that, quite generally, F T = −F as a consequence of the fermionic anticommutation. If we are in equilibrium (ground state or thermal) then both G and F are time-independent and can be taken to be real. Moreover, G is symmetric, G = G T . This implies that in agreement with the Majorana equilibrium correlators seen in Sec. 9.
Step 3: For every fixed value of t, we can transform the matrix A {l} (t) to a canonical form, by a (real) orthogonal transformation R. The canonical form of a real anti-symmetric matrix is, by Schur's decomposition, made of l anti-symmetric 2 × 2 blocks along the diagonal: with λ q=1···l real. Using the rotation matrix R, we can define 2l new combinations of Majorana fermions:ď q (t) = 2l n=1 R nq (t)č n for q = 1, · · · , 2l .
Now we switch back to ordinary fermions, simply because we are much more trained and used to think in terms of them. We therefore define: d q (t) = 1 2 ď 2q−1 (t) + iď 2q (t) for q = 1, · · · , l .
By construction, given the simple correlations encoded by the matrix Λ, these fermions have averages which shows that λ q ∈ [−1, 1], in order for the average fermionic occupation to be 0 ≤ P q (t) ≤ 1, and that the different fermionsd q (t) are uncorrelated.
The evaluation of the correlation matrix and the entanglement entropy can be implemented numerically with rather standard techniques.
Calculate the ground-state entanglement entropy for a uniform Ising chain, with l = L/2, in the three cases: a) h/J = 1/2, b) h/J = 2, and c) h/J = 1. Shows that the half-chain entanglement entropy tends to a constant, for increasing L, in cases a) and b), while it grows logarithmically in the critical case c).
Problem 7 Ground state entanglement for a uniform Ising chain.