Quantum spin systems versus Schrödinger operators: A case study in spontaneous symmetry breaking

Spontaneous symmetry breaking (SSB) is mathematically tied to some limit, but must physically occur, approximately, before the limit. Approximate SSB has been independently understood for Schrödinger operators with double well potential in the classical limit [1, 2] and for quantum spin systems in the thermodynamic limit [3, 4]. We relate these to each other in the context of the Curie–Weiss model, establishing a remarkable relationship between this model (for finite N) and a discretized Schrödinger operator with double well potential. Copyright C. F. van de Ven et al. This work is licensed under the Creative Commons Attribution 4.0 International License. Published by the SciPost Foundation. Received 30-11-2018 Accepted 15-01-2020 Published 06-02-2020 Check for updates doi:10.21468/SciPostPhys.8.2.022

Here ω ħ h is defined as a normalized positive linear functional on the C*-algebra A ħ h = B 0 (L 2 ( )) of compact operators on L 2 ( ). The algebraic formalism is particularly useful for combining classical and quantum expressions. 6 For example, if δV is positive and is localized to the right, then the relative energy in the left-hand part of the double well is lowered, so that localization will be to the left. See §4 for details. This Hamiltonian has a 2 -symmetry (σ 1 , σ 2 , σ 3 ) → (σ 1 , −σ 2 , −σ 3 ), which at each site x is 79 implemented by u(x) = σ 1 (x). The ground state of this model is unique for any |Λ| < ∞ 80 and any B = 0, and yet, as for the double well potential, in the thermodynamic limit it has 81 two degenerate ground states, provided 0 < B < 1. As explained in Ref. [7,[29][30][31]], perhaps 82 unexpectedly this limit actually defines a classical theory, with phase space B 3 ⊂ 3 , i.e. the introducing the discretization techniques we use to non-specialists. 115 2 Reduction of the Curie-Weiss Hamiltonian 116 Since the spatial dimension is irrelevant, we may as well consider the Curie-Weiss Hamiltonian
where we we have used the bra-ket notation in the above expression. By linearity, we may separately compute this for the operators In order to compute (2.6), we need to know how σ 3 and σ 1 act on the vectors |n + 〉.
, where the minus sign appears only 143 if e n x = e n y . We conclude that the standard basis for the N -fold tensor product is a set of 144 eigenvectors for σ 3 ( y)σ 3 (x) with eigenvalues ±1. Thus we know that x, y σ 3 (x)σ 3 ( y) is a 145 diagonal matrix with respect to this standard basis. Note that |n + 〉 is a (normalized) sum of 146 permutations of such basis vectors, with n + times the vector e 1 and N − n + times the vector 147 e 2 . Since x, y σ 3 ( y)σ 3 (x) acts diagonally on any of these vectors, and it is also permutation 148 invariant, it follows that in the inner product any vector with itself yields the same contribution, 149 namely 150 〈e 1 ⊗ · · · ⊗ e 1 ⊗ e 2 · · · ⊗e 2 | x, y σ 3 ( y)σ 3 (x)|e 1 ⊗ · · · ⊗ e 1 ⊗ e 2 · · · ⊗e 2 〉, (2.9) where e 1 occurs n + times and e 2 occurs N − n + times. The above expression equals since there are 2n + (N − n + ) minus signs and hence N 2 − 2n + (N − n + ) = n 2 + + (N − n + ) 2 plus 151 signs, so that in total the correct value is indeed given by This shows that the contribution to the diagonal is given by (2.3 Both sets are clearly disjoint. Then we compute (2.14) We used the fact that the vectors β n + ,l are orthonormal, that with n + − 1 = n + , and that In the next section we will argue that for 0 < B < 1 the above (N + 1)-dimensional matrix, 170 denoted by J N +1 , can be linked to a Schrödinger operator with a symmetric double well 171 on L 2 ([0, 1]), for N sufficiently large. 9 Since for a sufficiently high and broad enough where T ±a is the translation operator over distance a (i.e., (T ±a ϕ 0 )(x) = ϕ 0 (x ± a)), 177 where ±a denotes the minima of the potential well. N (for which the claim is a theorem) but also the first excited state ψ (1) N lies in the symmetric 184 subspace Sym N ( 2 ). 10 For N ≤ 60 (a number that obviously depends on the machine precision 10 This second claim is not essential for our results themselves but is helpful for the following explanation thereof. The point is that the first excited state computed from the tridiagonal matrix J N +1 might not the same as the one from the original Hamiltonian h CW N , in which case Figure 2 would be misleading. Fortunately, we have shown numerically that up to N = 12 the first excited state of h CW N represented as a matrix on 2 N is the same as the one corresponding to the tridiagonal matrix J N +1 . For N > 12 this computation became unfeasible, as the dimension of the relevant subspace grows exponentially with N . 11 We checked this numerically, but omitted the plots. Moreover, we observed that for increasing N these two numerical degenerate states were randomly localized in (one of the) both wells.
(numerically) degenerate ground state eigenvectors are given by the functions where the functions ϕ n (x) now have to be understood as functions on a discrete grid. Once 199 again, for N = 60, the fact that ψ 206 with grid points x i = i∆ and grid size ∆ = X /N . The symbol ∆ is called the grid spacing. Note 207 this the grid spacing is chosen to be constant or uniform in this specific example. For the first 208 order derivatives we have (3.2) These derivatives are approximated with finite differences. There are basically three types of 210 such approximations: Since it is more accurate in our case, will focus on the central difference approximation method 212 and apply this to the second order differential operator d 2 /d x 2 .
By throwing away the error term O(∆ 2 ) in the above equation, it follows that we can 218 approximate the second derivative operator in matrix form This matrix is the standard discretization of the second order derivative on a uniform grid 220 consisting of N points of length ∆ · N , with uniform grid spacing ∆. In this specific case, we 221 have ∆ = 1/N . We denote this matrix also by 1

222
Suppose now that we are given a symmetric tridiagonal matrix A of dimension N with constant 223 off-and diagonal elements: We are going to extract a kinetic and a potential energy from this matrix. We write where the latter matrix is a diagonal matrix with the element b +2a on the diagonal. It follows for T = a[···1 −2 1···] N , and V = diag(b+2a). In view of the above, the matrix T corresponds 228 to a second order differential operator. This matrix plays the role of (3.5), but with uniform 229 grid spacing 1/ a on the grid of length N / a. Since the matrix V is diagonal, it can be seen 230 as a multiplication operator. Therefore, given such a symmetric tridiagonal matix A, we can 231 derive an operator that is the sum of a discretization of a second order differential operator 232 and a multiplication operator. The latter operator is identified with the potential energy of the 2 , is precisely an example of such a matrix.

238
The question we ask ourselves is if we can link such a matrix to a discretization of a Schrödinger 239 operator as well (see Appendix B for background).

240
Writing T = J N +1 , consider the ratios with non-uniform grid spacing h j and h j−1 . We divide the original tridiagonal matrix J N +1 by 242 N for scaling. Thus, we consider J N +1 /N . If we then compute the distances h j , we see that 243 they are almost all of O(1), except at the boundaries. We will see later that the corresponding 244 Schrödinger operator analog of the matrix J N +1 /N will be an operator on a domain of length First, we compute the ratios ρ j : Use the following approximations (3.11) and observe that for j >> 1 and N − j >> 1, we have (3.14) Moreover, we see that the ratio satisfies using the fact that that the big-O notation respects the product, that O( 1 In the next subsection, we will see from numerical simulations that to a good approximation  This is an approximation, since we neglect the (relatively small) function values of the Gaussian 262 that are more than O( N ) away from the central maximum. However, this approximation is 263 highly accurate, as the Gaussian decays to zero exponentially. This observation is extremely 264 important, as we will now see.

265
Let us first focus on the left-located Gaussian. For a point x j = j/N , clearly j ∈ O(N ).

266
Therefore, for N large enough, , and we find (3.17) We will now show that, in the present context where we work on a domain of order L (i.e. corresponding to the minimum x j 0 of the potential well by T j 0 , for the off-diagonal elements 275 within a range of order σ, we derive the next estimate: where we used the inequality (1 + 1/N ) σ ≤ 1 + C σ N as well as the fact that T j 0 is of order 1.

277
Here, C > 1 is a constant independent of N . Since the left peak of the Gaussian eigenfunction 278 is approximately non-zero corresponding to a matrix segment of length of order N , we 279 apply the above estimate to σ ≈ N . We see immediately that |T j 0 − T j 0 +σ | goes to zero.

280
Therefore, on matrix segment of length of order N centered around the left minimum x j 0 281 of the potential, the off-diagonal elements coincide in the limit N → ∞. This means that 282 the grid spacing becomes constant and hence that we have locally uniform discretization of 283 the domain. By symmetry, the same is true for the peak located on the right of the well. We 284 conclude that for large N the tridiagonal matrix locally behaves like a kinetic energy, and 285 therefore like a discretized Schrödinger operator. All this will be explained in more detail in 286 the next subsection.

288
Our aim is to show that for N large enough, the matrix J N +1 /N obtained from the Curie-Weiss

289
Hamiltonian by reduction (see §2) locally approximates a discretization matrix representing 290 a Schrödinger operator describing a particle moving in a symmetric double well. This means 291 that there exists a sub-block of J N +1 /N that has a form approximately given by the sum of (for a certain h) and a diagonal matrix playing the role of a potential. 293 We started with the symmetric tridiagonal matrix J N +1 /N with non-constant entries. In order 294 to link this matrix to a second derivative and a multiplication operator, we needed to apply 295 the non-uniform discretization procedure.

296
At first sight the off-diagonal matrix entries of J N +1 /N cannot immediately be identified 297 with a second order derivative operator, but we have seen that in the limit N → ∞ we do V N is a diagonal matrix (and hence a multiplication operator, as in the continuum) given by Note that d(x) corresponds to the non-uniform grid spacing. We rewrite: and hence the total length of the interval is given by 1 Moreover, the interval coordinate is then given by It follows that on each matrix segment of N entries (for N large), the matrix H N is where x is chosen appropriately and the segment describes an interval of length (3.27) in the total interval of length L. 14 We saw in the previous section that on length scales of order is approximately constant so that T N can be seen as a locally 317 uniform discretization of the second order derivative ∆ on the sub interval approximately Therefore, on these length scales the operator 319 ∆ indeed corresponds to a uniform part of We see in this section that to a very good approximation the spectral properties of both (a 321 priori different) matrices J N +1 /N andH N +1 , defined in (3.34) below, coincide, improving 322 with increasing N . Rescaling the interval [0, L] to unity yields a Schrödinger operator h, but 323 now defined on L 2 ([0, 1]). Hence h is given by where, via the potential V N defined in (3.22), the potentialṼ N is defined by Let us now consider the scaled tridiagonal matrix J N +1 /N . By definition it follows that the 326 elements on the diagonal, the lower diagonal, and the upper diagonal are given by respectively. The idea is to split the matrix J N +1 /N into two parts, one corresponding to the 328 kinetic energy and the other to the potential energy. However, since the off-diagonal elements 329 of J N +1 /N are non-constant and not an even function around x = 1/2, we cannot isolate these 330 elements and decompose J N +1 /N into two parts. Therefore, we approximate the off-diagonal 331 elements by the function (3.21), which is even in 1/2. This approximation makes perfect sense 332 in the semi-classical limit, and from this it is clear that where each segment with order N matrix entries now describes an interval of length (3.32) which indeed coincides with (3.28). The matrix (3.31) corresponds to a non-uniform 341 discretization of h, in such a way that it is locally uniform, namely on length scales of order 342 N . Finally, discretizing ∆ y on the unit interval with uniform grid spacing 1/N yields the 343 following discretization matrix: (3.33) Using this discretization of − 1 L 2 N 2 ∆ y , we show that for N → ∞ the spectrum of the matrix 345 J N +1 /N approximates the spectrum of the matrix Using the fact thatH N is a discretization of (3.28), this indeed establishes a link between the Theorem (Appendix A) that the ground state is always unique for any finite N . 16 We also the Curie-Weiss model will not display SSB, not even in the classical limit. For these values of 378 B, the well will be a single potential, as depicted in Figure 7 for B = 2.

379
In view of the corresponding Schrödinger operator, the ground state in the classical limit will 380 not break the symmetry for a single potential well, and is therefore also compatible with the corresponding to the sum of the uniform discretization K N of the second order derivative (viz.

384
(3.33)) andṼ N (viz. (3.29)). We will see that to a very good approximation the spectral 16 Due to this degeneracy, the computer picks or the symmetric combination χ + , or the anti-symmetric combination χ − , which depends on the algorithm. We changed N and observed that the location of the peak changed as well. This suggests that random superpositions of the two degenerate states are formed.  been compared to those of J N +1 /N (Table 1). We computed the first ten eigenvalues of the 388 matrix J N +1 /N , denoted by ε n , and those ofH N , denoted by λ n . In the left column the 389 eigenvalues ε n are displayed. In the right column the absolute difference |λ n −ε n | is displayed.  We see that the first ten eigenvalues for both matrices are the same up to at least six decimals. It 394 is also clear that these eigenvalues are doubly degenerate, at least up to six decimals. Moreover, 395 we plotted all the eigenvalues ε n and λ n corresponding to the bound states, i.e., the energy 396 levels within the well. This is depicted in Figure 8 below.  On the one hand, it is clear that this width will go to zero as N → ∞. On the other hand, also 413 the unit interval depends on N , as the latter has to be discretized with N +1 points. Therefore, 414 the total number of grid points in the ground state peak living on a subset of order N is (3. 35) In fact, due to the discretization of the grid we have a better approximation of the Gaussian 416 ground state when N increases. 417 We have computed the minimum of the potential, and subtracted this minimum from the 418 lowest eigenvalues. Then, we have set the potential minimum to zero. These shifted 419 eigenvalues then live in a positive potential well. For J N +1 /N with N = 1000, we now 420 consider its eigenvalues ε n . We have already seen above that the lowest eigenvalues of J This firstly shows that for large N the wells approximately decouple since tunneling is 434 suppressed in the limit, and secondly that each well of the double well potential is locally 435 quadratic and therefore approximately has the spectrum of a harmonic oscillator (the latter 436 approximation increasingly breaks down at higher excitation energies, however).

437
Finally, both tables have been computed for fixed N , but also different values of N need to be 438 considered. Table 3 shows the ground state eigenvalue ε N 0 of the matrix J N +1 /N .
where the potential V is C ∞ , non-negative, strictly positive at ∞, zero at two points m 1 and to one of these minima, in the following sense. We use the following notation: is the Agmon distance between the two minima of V ;

456
• d 1 = 2 min{d V (m 1 , supp δV ), d V (m 2 , supp δV )} is twice the Agmon distance between 457 the support of δV and the minimum that is closest to to this support;

458
• d 2 = 2 max{d V (m 1 , supp δV ), d V (m 2 , supp δV )} is twice the Agmon distance between 459 the support of δV and the minimum that is furthest away from this support.

460
Finally, and perhaps most crucially, as ħ h → 0 the perturbation should dominate the energy Though perhaps surprising at first sight, this is actually easy to understand, either from 469 energetic considerations or from a 2 × 2 matrix analogy [2,7]. First, the ground state tries 470 to minimize its energy according to the rules: The eigenvalues and eigenvectors of h (2) ħ h , respectively, are given by

and the (metaphorically) localized states would be
Now take δ > 0 and introduce a "flea" perturbation by changing h (2) ħ h to h The eigenvalues of h (2) ħ h + δ (2) (4.10) As long as lim ħ h→0 ∆(E) ħ h /δ = 0, we have so that under the influence of the flea perturbation the ground state localizes as ħ h → 0. There 483 is no need for a separate limit δ → 0, since it follows from the limit ħ h → 0.

484
Returning to the real thing, a (mathematically) very natural flea-like perturbation δV for the 485 Schrödinger operator h ħ h , and the one we shall mimic for the Curie-Weiss model, is where the parameters (b, c, d) represent the location of its center b, its width 2c and its height 487 d, respectively. Tuning these, the conditions above can be satisfied in many ways: for example,

491
The next step in our analysis, then, is to find an analogous perturbation to (4.12) but Here the suffix l in β k l labels the basis vector β k ∈ β within the same orbit O k . So for each orbit O k , we have N k vectors β k . Hence for each l = 1, ..., N k the image S N (β k l ) under S N is 508 always the same, namely the coordinate vector written with respect to β. 509 The perturbation we are going to define is very similar to the symmetriser S N . Of course, since 510 we have expressed our original Curie-Weiss Hamiltonian with respect to this |n + , n − 〉 basis, 511 we need to do the same for the perturbation. Since we have a partition of our 2 N -dimensional 512 basis β into N + 1 orbits, we define our perturbation ∆V N by where k ∈ {0, ..., N } is the unique number such that β ∈ O k , and δV b,c,d is defined by (4.12). 514 We will see that a specific choice of parameters results in localization of the ground state and that the ground state eigenfunction of the perturbed Hamiltonian h N + ∆V N is unique 520 and positive, 17 then we may conclude that the ground state lies in the subspace Sym N ( 2 ).

521
The reason for this is the same as for the unperturbed Curie-Weiss Hamiltonian: these 522 properties push this eigenvector into the subspace ran(S N ) = Sym N ( 2 ), so that we may 523 diagonalize this Hamiltonian represented as a matrix that can be written with respect to 524 the symmetric subspace, which will be a tridiagonal matrix of dimension N + 1 as well. never zero, the matrix can never be decomposed into two blocks, so that it remains irreducible.

537
Non-negativity is achieved when which is clearly satisfied for positive d. Therefore, taking d > 0 together with the fact that  (4.20) We now show numerically that our flea perturbation ∆V N forces the ground state to localize 546 for large N , leaving an analytic proof à la Simon (1985) to the future.

547
As in section 3, we extract the potential corresponding to the perturbed Hamiltonian h N +∆V N , 548 written with respect to the symmetric basis, scale this Hamiltonian by 1/N , and translate 549 the potential so that its minima are set to zero. We plot this perturbed potential on the unit 550 interval in Figure 9, where for convenience we scale the domain to the unit interval. Moreover, that the ground state is unique, and hence localization is not a result of numerical degeneracy 555 but is genuinely caused by the perturbation. A similar simulation for flea perturbation, but 556 now located on the left site of the barrier, is shown in Figure 11. As expected, we now see a 557 localization of the ground state to the right side of the barrier (Figure 12).

558
Our conclusion is that due to the flea perturbation, the ground state will localize in one of 559 the wells depending on where the flea is put. As in the continuous Schrödinger operator  and also in this case a collapse of the ground state takes place. This is also clear from the 2 × 2 570 matrix argument given in the previous section.  where we add that the same is true as ħ h becomes smaller at fixed system size, cf. [7] and 598 references therein, and that the precise form of the instability is that for small N or large ħ h the 599 perturbation plays almost no role (in either the spectral properties of the Hamiltonian or in its 600 eigenfunctions), whereas for large N or small ħ h it metaphorically 601 "can irritate the elephant enough so that it shifts its weight, i.e., we will see 602 that the ground state, instead of being asymptotically in both wells, may reside 603 asymptotically in only one well".

604
As explained in the Introduction, this accounts for the fact that real materials (which are 605 described by the quantum theory of finite systems) do display SSB, even though the theory 606 seems to forbid this. As we (and most condensed matter physicists) see it, these perturbations to the Curie-Weiss Hamiltonian (2.1) and arguing that the correct order of the limits in question is lim →0 lim N →∞ , which gives SSB by one of the two pure ground states on the limit algebra, the sign of determining the direction of symmetry breaking. In contrast, the 620 opposite order lim N →∞ lim →0 gives a symmetric but mixed and hence unstable or unphysical 621 ground state on the limit algebra. 20  combination of low-lying states) converges to some symmetry-breaking pure ground state on 637 the limit system (be it a classical system or an infinite quantum system). The ensuing limit is 638 then duly continuous in an appropriate meaning of the word (which it would not be without 639 the perturbation mechanism). 24

640
Of course, in order for this symmetry breaking to be spontaneous rather than explicit, the 641 perturbation should be small to begin with, and should disappear in the pertinent limit. 25

642
As we have seen, it is easy to endow the "flea" δV with these properties; for symmetry breaking 643 in the double well potential all we need is that ∆E → 0 more rapidly than δV → 0 as ħ h → 0, . It may seem more natural to just require that the ground state fails to be G-invariant, but in the C*-algebraic formalism we rely on ground states that are not necessarily pure, which leaves the possibility of forming G-invariant mixtures of non-invariant states that lose the purity properties one expects physical ground states to have. Similarly for equilibrium states, where 'pure' is replaced by 'primary', which is a mathematical property of a pure thermodynamical phase. Order parameters follow from this definition, cf. §10.3, loc. cit. 22 And similarly for the limit ħ h of the symmetric double well system, where, as pointed out in the Introduction, the ħ h = 0 limit of the ground state is the mixed state (1.4), as opposed to, for example, a Dirac delta-function located between the wells (i.e. at q = 0). See also Refs. [2,43]. This also explains the fact that the flea perturbations even work if their support is localized away from the bottoms of the two wells (or, equivalently, from the two peaks of the unperturbed ground state); see especially Ref.
[26] for a detailed explanation of this point. take x = δV and y = ∆E. In case of SSB, certain physical quantities m(x, y) (like the 652 magnetization) are discontinuous at (0, 0) and hence their value at (0, 0) depends on the 653 path towards the origin. Eq. (5.26) expresses this path dependence in a crude way, which 654 is captured by perturbations like (5.25), but the perturbations we consider follow specific

A digraph is called strongly connected if there is a directed path x to y between any two
698 vertices x, y. 699 We use the notion of the directed graph or digraph of a square N -dimensional matrix a, 700 denoted by G(a). We say that the digraph of a is the digraph with The following result links irreducibility of a non-negative matrix to strongly connectedness of 702 its corresponding digraph. The proof is easy and therefore omitted.  attains its minimum in n + = N /2, which is given by N 2 − 4 N 2 n + 4( N 2 ) 2 = 0. So indeed, there 731 are at least as many plus signs as minus signs, so that the corresponding diagonal term is 732 non-negative. The other term x∈Λ N σ 1 (x) does not contain any negative entries at all, so if 733 we apply this to any basis vector {e n 1 ⊗ ... ⊗ e n N }, we get a non-negative matrix. It follows that 734 both operators in (A.1) are non-negative in the basis under consideration.

735
Now we show that the matrix corresponding to the Curie-Weiss Hamiltonian is irreducible.

736
Note that irreducibility of a matrix does not depend on the basis in which the operator 737 is represented, since similar matrices define equivalent representations which preserve 738 irreducibility. We use Lemma A.2 to show that there is a direct path between any two vertices.

739
But this is obvious: the operator x σ 1 (x) flips the spins one by one, and therefore the 740 associated digraph is clearly strongly connected as we can find a directed path between any The approximations are obtained by neglecting the error terms indicated by the O-notation. The question is how small h has to be in order for the algebraic difference f (x+h)− f (x) h (i.e. in 761 this case the forward difference approximation) to be good approximation of the derivative.

762
It is clear from the above formulas that the error for the central difference formula is O(h 2 ).

763
Thus, central differences are significantly better than forward and backward differences. A forward difference approximation to f (x) is then and a centered difference approximation is for example Now we discretize the kinetic and potential energy operator. For simplicity, consider the The approximations are again obtained by neglecting the error terms.

775
Using this uniform grid with grid spacing h = 1/N , it follows that the second derivative 776 operator in one dimension is given by the tridiagonal matrix 1 h 2 [· · ·1 −2 1 · ··] N and the 777 potential which acts as multiplication, is given by a diagonal matrix. With the notation