Generalized Phase-Space Techniques to Explore Quantum Phase Transitions in Critical Quantum Spin Systems

We apply the generalized Wigner function formalism to detect and characterize a range of quantum phase transitions in several cyclic, finite-length, spin-$\frac{1}{2}$ one-dimensional spin-chain models, viz., the Ising and anisotropic $XY$ models in a transverse field, and the $XXZ$ anisotropic Heisenberg model. We make use of the finite system size to provide an exhaustive exploration of each system's single-site, bipartite and multi-partite correlation functions. In turn, we are able to demonstrate the utility of phase-space techniques in witnessing and characterizing first-, second- and infinite-order quantum phase transitions, while also enabling an in-depth analysis of the correlations present within critical systems. We also highlight the method's ability to capture other features of spin systems such as ground-state factorization and critical system scaling. Finally, we demonstrate the generalized Wigner function's utility for state verification by determining the state of each system and their constituent sub-systems at points of interest across the quantum phase transitions, enabling interesting features of critical systems to be intuitively analyzed.


I. INTRODUCTION
Quantum many-body systems often present interesting and unexpected properties that have no classical counterparts, and which are ultimately due to the strongly correlated and highly entangled nature of the underlying many-body states. Quantum spin-lattice systems in particular are able to elucidate some of the underlying properties and characteristics of more general many-body systems by providing a simple yet rich and valuable framework that constitutes the basic models of many physical systems such as magnetic insulators [1][2][3] and coupled qubits [4]. The study of spin-lattice models of quantum magnetism has thus become an area of enormous theoretical study in recent years. Additional impetus has come from the growing realizations that quantum many-body systems themselves provide valuable resources for quantum computing (see, e.g., Refs. [4][5][6]), and that quantum information theory concepts can themselves provide valuable insights into the properties and behaviors of quantum many-body systems, typically via a detailed study of the intimate entanglement structure of many-body wave functions [7].
From a modern perspective quantum entanglement [8] has become a key element of various quantum technologies, due to the significance of entangled states as a re- * raymond.bishop@manchester.ac.uk † m.j.everitt@physics.org source in applications such as the development of protocols for teleportation, secure quantum key distribution schemes, etc. Means to produce various entangled states both systematically and efficiently are hence of growing importance. For example, Einstein-Podolsky-Rosen (EPR) pairs, which may be considered as qubit pairs in the form of a Bell state, furnish one of the simplest illustrations of maximum entanglement at the bipartite level.
The demonstration that such EPR pairs can be rapidly and robustly generated with spin chains [9] has also provided additional impetus for their study.
One especially interesting feature of spin system models is their ability to exhibit and transition between different quantum phases [10,11]. Such phase transitions occur in interacting many-body ensembles in the thermodynamic limit. The more well-known thermal phase transitions that occur at specific critical temperatures T c are due to thermal fluctuations that represent the competition between the entropy and the internal energy of the system, which themselves depend on such parameters as the pressure (P ), volume (V ) and temperature (T ). Such thermal phase transitions (in one-component systems) are typically studied by thermodynamic free energy state functions such as the Gibbs free energy, G(P, T ) or Helmholtz free energy, F (V, T ), which become doubly degenerate at the corresponding phase boundary between two thermal phases.
By contrast, the quantum phase transitions (QPTs) [12][13][14] of interest here occur at zero temperature (T = 0) between differing ground states of the system, and are due to quantum fluctuations that represent the competition between at least two non-commuting terms in the system Hamiltonian, each of which by itself promotes a different form of ground state. When acting together they therefore act to frustrate one another. A prototypical frustrated spin-lattice system is the spin- 1 2 J 1 -J 2 model on a one-dimensional (1D) chain [15,16] or a two-dimensional square lattice [17][18][19], in which the two competing terms in the Hamiltonian are isotropic Heisenberg spin-spin interactions between nearest-neighbor and next-nearest-neighbor pairs, with strength parameters J 1 ≡ λ 1 and J 2 ≡ λ 2 , respectively.
If the set of variables {λ 1 , λ 2 , · · · , λ n } represents the strength parameters of the various non-commuting parts in the model Hamiltonian, then the groundstate energy of the system, E 0 (λ 1 , λ 2 , · · · , λ n ), correspondingly becomes doubly degenerate at the boundary (λ c 1 , λ c 2 , · · · , λ c n ) between two quantum phases. Such QPTs are said to be a first-order QPT (1QPT) or a second-order QPT (2QPT, also known as a continuous QPT) depending, respectively, on whether the firstorder partial derivatives {∂E 0 /∂λ i , i = 1, 2, . . . , n} are discontinuous or continuous at the transition boundary. Whereas a 2QPT is typically characterized by the existence of an infinite correlation length in the system and a corresponding power-law decay in its correlations, which is often expressed via a divergence (or finite discontinuity) in one or more second-order derivatives of the ground-state energy, there also exist infinite-order QPTs (∞QPTs) at which all finite-order derivatives of the ground-state energy remain continuous. In the present work we will consider spin-lattice systems that display all three types of QPTs (i.e., 1QPT, 2QPT, ∞QPT).
The study of QPTs has itself acquired additional impetus from the increasingly widespread use of ultra-cold atoms trapped in optical lattices formed by a periodic potential, which itself has been created by standing waves formed from a suitable array of lasers, in order to simulate a wide variety of of systems of interest in condensed matter and many-body physics [20][21][22][23], especially spinlattice systems. In such simulations it is then often experimentally possible to vary the strength parameters {λ 1 , λ 2 , · · · , λ n } discussed above, and hence to map out a QPT and to find its boundary (λ c 1 , λ c 2 , · · · , λ c n ) in this parameter space. Examples include that from a superfluid to a Mott insulator [24,25], as well as many others in quantum magnetism, both for one-dimensional spin chains (see, e.g., Ref. [26]) and for periodic lattices in two or more dimensions (see, e.g., Refs. [20,27,28]). For example, a two-dimensional regular honeycomb lattice may be formed by interfering three coplanar laser beams propagating at relative angles of ±120 • , and Duan et al. have shown how to engineer specific spin Hamiltonians on such a lattice accordingly [20]. Similarly, concrete proposals to form optical lattices representing honeycomb-lattice bilayers in both AA stacking [27] and AB (or Bernal) stacking [28], both of which use five lasers, have also been discussed.
Exact solutions for the ground-state energy of quantum many-body systems in general [29], and spin-lattice systems in particular [2], are known only in a few special cases. Otherwise, one has to resort to various tools of quantum many-body theory to find approximate solutions. Furthermore, witnessing critical behavior in these systems, even when some exact results are known, can be quite difficult and cumbersome (see, e.g., Refs. [11,30]). However, the crucial role that quantum correlations play in critical many-body systems has highlighted the use of entanglement and correlation measures as possible methods for witnessing and characterizing QPTs [31,32]. The subsequent application of quantum information tools to critical quantum systems has by now firmly demonstrated their unparalleled ability to measure and characterize quantum critical behavior in a wide range of manybody systems [33][34][35][36][37][38][39][40][41][42][43].
Within the wider context of the phase-space formulations of quantum mechanics and the associated quantum distribution functions (see, e.g., Ref. [44]), the Wigner function has found a large number of applications to a wide variety of problems in physics. Examples include the calculation of quantum mechanical observables in quantum ballistic transport studies [45] and many other areas such as quantum optics, quantum electronics, quantum chemistry, signal processing, and quantum information theory. Reference [46] presents a nice review of some of these latter applications.
More recently, it has been realized that the Wigner function in particular can also provide a useful measure of quantum correlations [47][48][49][50]. Consequently, it is natural to explore the application of this formalism to the study of QPTs and quantum critical behavior. Extensions of the original Wigner function formalism for continuous systems in the usual position-momentum phase space [45,51,52] to discrete, finite-dimensional systems, such as ensembles of spins, has been a challenging endeavor.
There have been many attempts to create a finite analog to the Wigner function [53][54][55][56], each with their own formulation. Here we follow the formulation of a generalized Wigner function (GWF) introduced in Ref. [57], that is able accurately to construct a complete, continuous Wigner function for any arbitrary quantum system, as well as its corresponding Weyl function [58]. The latter work in particular has completed the phase-space formulation of quantum mechanics by extending the original work of Wigner [51], Weyl [59], Moyal [60], Groenewold [61], and others to any quantum system.
Research into the application of Wigner function formalisms to the study of critical quantum spin systems was recently carried out by Mzaouali et al. [62]. They successfully demonstrated the utility of phase-space techniques for witnessing, characterizing and distinguishing first-, second-and infinite-order QPTs. In addition to this, they demonstrated the ability of the method to detect more nuanced features of these systems such as ground-state factorization. In this work we aim both to verify these results and to build upon them by applying the GWF formalism to several finite-sized (N = 6) spin-1 2 one-dimensional (1D) models, viz., the Ising ferromagnetic chain in a transverse magnetic field, the anisotropic XY chain in a transverse magnetic field, and the XXZ chain. We will explore their individual properties and critical behavior by providing an exhaustive exploration of the correlations present within these systems. We will further demonstrate the ability of the GWF fomalism to witness and characterize first-, second-and infinite-order QPTs, as well as ground-state factorization, through single, bipartite and multi-partite correlation functions. In addition, we demonstrate the formalism's ability to capture scaling of finite spin chain systems. We also employ the GWF's ability for intuitive state analysis [63] by visualizing the state of the system and its constituent sub-systems at points of interest across QPTs through spin Wigner function plot visualizations. In turn, we demonstrate the utility of this visualization technique specifically to witness features of infinite-order QPTs and ground-state factorization.
The remainder of this paper is organized as follows. In Sec. II we outline the GWF formalism and discuss its application to spin systems. We then apply this formalism in Sec. III to the three spin-1 2 chain models of interest here that we have cited above. For each model we produce a phase line plot to witness critical behavior in the system and then determine the state of the system at points of interest across the phase line through the use of spin Wigner function visualizations and reference state plots. Finally, in Sec. IV we conclude and summarize our findings.

II. THE WIGNER FUNCTION
We make use of a continuous, informationally complete, spin Wigner function [46,57,64] with similar properties to the one originally formulated by Wigner in 1932 [51]. For a two-level quantum system it takes the form of a quasi-probability distribution defined over the same Euler angles as the Bloch sphere. The spin-Wigner function can be expressed in terms of the expectation value of an analog to the usual displaced parity-operator, given by∆ where the Euler angles θ ∈ [0, π] and φ ∈ [0, 2π) are set by the rotation operatorR(θ, φ, Φ) ≡ e −iσ z φ/2 e −iσ y θ/2 e −iσ z Φ/2 and parametrize the phase space, andΠ ≡ 1 . For a composite system of qubits the corresponding displaced parity operator is simply the tensor product, where N is the number of spins. The Wigner function, whereρ is the total density matrix of the spin system, is a function on the product of N spheres. For the purposes of visualization (see, for example, Fig. 4) it often suffices to restrict to the equal-angle slice; In this representation, high-order spherical harmonics are then a witness to many-spin entanglement (see, for example, Fig. 7 of Ref. [63]). Following Ref. [62] we will also use, as a witness of phase transitions, just one specific point of the Wigner function equal-angle slice, W EA (0, 0) (see, for example, Fig. 3). The GWF is our tool to visualize the density matrixρ (orρ tot ) of the ground state of the full system. We will also visualize correlations through the corresponding reduced density matricesρ I , where I ⊆ tot is a subset of the spin labels tot = {i : 1 ≤ i ≤ N }. The reduced density matrix is the partial trace of the density matrix over the remaining qubits tot \ I. The corresponding reduced Wigner function is the marginal obtained by integrating the full Wigner function over the angular variables of the remaining qubits. The correspondence follows from repeated use of Our choices of reduced Wigner functions for the 6-spin systems studied hare are shown in Fig. 1, where we label ρ tot ≡ W EA (θ, φ). We first calculate a reduced Wigner function W I (θ, φ) by integrating out the non-selected spin degrees of freedom from the full Wigner function in Eq. and then take the equal-angle slice to plot it as the corresponding correlation function ρ I . For instance, for a system with N = 6 spins, ρ 12345 is the equal-angle slice of the reduced Wigner function such that Likewise, all of the correlation functions shown in Fig. 1 refer to the equal-angle slices of the corresponding reduced Wigner functions.

III. APPLICATIONS TO SPIN SYSTEMS
In this Section we apply the generalized Wigner function formalism to a range of physical models: the spin-1 2 Ising ferromagnetic cyclic chain in a transverse magnetic field, the spin-1 2 anisotropic XY cyclic chain in a transverse magnetic field, and the spin-1 2 XXZ anisotropic cyclic Heisenberg chain. These models have been chosen due to the wide range of quantum phase transitions that they exhibit [2], which in turn enables us to demonstrate the utility of phase-space techniques in exploring quantum critical behavior. Each model is considered as a 1D chain of identical spin-1 2 particles, with periodic boundary conditions and a finite number of spins (N = 6). The finite number of spins enables exact solutions of the ground state to be calculated through direct diagonalization of the respective Hamiltonian of each model [2]. The finite size also keeps the dimensionality of the system relatively low to enable more intuitive analysis of the results. In addition to this, the low number of spins and the GWFs ability to visualize high dimensional systems enables an exhaustive exploration of single, bipartite and multipartite correlations within the models, which are known to provide valuable insights into critical behavior [41,66]. Figure 1 shows all correlations explored in these models. Previous demonstrations of critical behavior in small finite-sized systems have demonstrated their utility in exploring multi-body quantum systems and QPTs [67,68]. Explicit discontinuities in the phase lines at critical points and more exotic phases are not always expected however, as our models are not implemented in the thermodynamic limit (N → ∞) [12,13]. Due to this limitation we have chosen to include spin Wigner functions plotted on the Bloch sphere for each correlation function of the system to better witness and characterize critical behavior by enabling a deeper appreciation of the state of the system at points of interest. This also enables us to infer the state of the system and its subsystems at these points. As we show in detail later, we can thereby infer the presence of critical behavior in an infinite system. Figure 2 shows a collection of reference spin Wigner functions that will be used to interpret the results presented below.
A. Spin-1/2 Transverse-Field Ising Spin Chain Model The Hamiltonian for a spin-1 2 1D transverse-field Ising model with periodic boundary conditions is given bŷ where N is the number of spins, λ is the strength parameter of the nearest-neighbor Ising exchange interactions (which are ferromagnetic when λ > 0, as henceforth assumed here),σ α i , α = x, y, z are the usual Pauli matrices, and h is the external magnetic field strength [2,34,69,70]. Periodic boundary conditions are assumed, so thatσ α N +1 =σ α 1 , α = x, y, z. We are at complete liberty to set the overall energy scale, and for simplicity we will henceforth set h = 1.
Let us first consider the model at the classical level, in which case σ i (= 2s i ) is an arbitrary vector of unit length. In the presence of the external magnetic field the classical spins all cant at an angle α to the x axis, and it is easy to see that the classical ground-state energy, E cl 0 , is minimized for λ ≥ 1 2 by the choice α = sin −1 ( 1 2λ ). A that corresponds to the direction on the Bloch sphere. In (b)-(l) we now need to display the functions in the the equal-angle slice. The first of these states is the two-qubit product state |↑↑ in (b). Since this is a doubled product state, each point in the equal-angle slice is equal to the squared value of the equivalent point in the single-qubit phase space, resulting in an equal-angle slice that is non-negative. If instead we want to take a product of the state |↑ with the orthogonal |↓ , we get the state in (c). The single-qubit Wigner function of |↓ is the same as in (a) with the maximum of the distribution pointing at the south pole, resulting in a Wigner function that is flipped around the equator. The equal-angle slice of the product of these two states is then the point-wise multiplication of the individual distributions, resulting in (c). Alternatively one may be interested in two-qubit entangled states that start to manifest correlations not possible in product states. These can be seen in (d) and (e), that give us an insight into the behavior of entangled states in phase space. A comparable four-qubit state to that of (e) is shown in (f). More iconic is the appearance of the GHZ states, GHZ ± z in (g) and (h), that manifest in phase space as two coherent states, one on each of the poles, with oscillating positive and negative values around the equator; note that the number of oscillations matches the number of qubits that make up the GHZ state. Alternative entangled states are shown in (i) and (j) where entanglement isn't between aligned states; such states will be found later. We will also need to consider mixed states as reduced Wigner functions will also be important in our analysis. When taking the reduced Wigner function for one qubit of a multi-qubit maximally entangled state, the result is the fully mixed state in (k) that presents a uniform value over all phase space. Further, when we remove the quantum correlations of a maximally entangled GHZ state, the result is a statistical mixture of two spin coherent states where the quantum correlations around the equator are no longer present. An example of this can be seen in (l).
classical phase transition then ensues at λ = λ cl c = 1 2 such that for λ < λ cl c the classical spins all align along the z direction of the external field, with α = π 2 . The groundstate energy per spin is thus given at the classical level by Clearly, both the energy and its first derivative with respect to λ are continuous at the classical transition point. Similarly, the classical values of the components of the magnetization (viz., the particle spin) vector in both the x (viz., the Ising) and the z (viz., the transverse field) directions are given by and The quantum spin-1 2 version of the model is also exactly solvable for both finite values of N and in the thermodynamic limit (N → ∞) through the following sequence of transformations: (i) a Jordan-Wigner transformation [71], which transforms the spin operators into fermionic operators; (ii) a Fourier transformation from lattice position to lattice momentum; and (iii) a Bogoliubov transformation [2,70]. It is important to note too that the model possesses a Z 2 symmetry group due to its invariance under the unitary operation of flipping all spins in the x direction. It is precisely the breaking of this symmetry that causes the sole QPT that this model possesses in the thermodynamic limit (N → ∞) [70]. The exact ground-state energy per spin for the quantum spin-1 2 model of Eq. (8) in the thermodynamic limit (N → ∞) is given by [70] which is an elliptic integral of the second kind. Although the expressions of Eqs. (9) and (12) are very dissimilar, they agree numerically with one another to a few percent or less. For example, at λ = 1, we have E cl 0 (λ = 1)/N = −1.25 and E 0 (λ = 1)/N = − 4 π ≈ −1.273. It is also trivial to confirm that in the two limiting cases λ → ∞ and λ → 0, the expression of Eq. (12) reduces respectively to the values E 0 (λ → ∞)/N → −λ and E 0 (λ = 0) = −1, which are exactly also as given by the classical expression of Eq. (9). The agreement of the classical and quantum results for the ground-state energy in these two limiting cases is exactly as expected, since in both of these two extremes the spins are fully aligned, and such fully aligned states are also eigenstates of the corresponding quantum Hamiltonian.
Despite the above levels of numerical agreement, the classical and quantum results for the ground-state energy differ in one profound aspect. Thus, the integral expression in Eq. (12) is nonanalytic at the point λ = λ c = 1, which is now the quantum phase transition point, which is considerably shifted from its classical counterpart at the quite different value λ = λ cl c = 1 2 . One may readily confirm from Eq. (12) that both E 0 and ∂E 0 /∂λ are continuous and finite at λ = λ c = 1, whereas all higherorder derivatives ∂ n E 0 /∂λ n with n > 1 diverge at the same point, so that the transition there is of 2QPT type.
The exact result for the x component (viz., in the Ising direction) of the magnetization (i.e., the groundstate expectation value of the spin vector, s i = 1 2 σ i , at a given site) for the system in the thermodynamic limit (N → ∞) is given by [70] which exhibits the quantum phase transition at λ = λ c = 1 much more clearly than does Eq. (12) for the groundstate energy. Finally, the corresponding exact result for the z component (viz., in the transverse field direction) of the magnetization for the system in the thermodynamic limit (N → ∞) is given by [70] M z = 1 2π an elliptic integral of the first kind. Once again, Eq. (14) is nonanalytic at the quantum phase transition point, λ c = 1, where it may also readily be evaluated, M z (λ = 1) = 1 π . One may readily confirm from Eq. (14) that as λ increases from zero to infinity M z decreases smoothly and monotonically from 1 2 to zero, the same extreme values as in the corresponding classical result of Eq. (11).
Thus, in the limit λ → ∞, the model has perfect ferromagnetic order with all the spins aligned along either the positive or negative x direction. This two-fold degeneracy of the ground state is maintained for the system in the thermodynamic limit as λ is decreased (and, accordingly, the magnetization in the spin-space x direction decreases from its maximal value at λ → ∞) until at the critical point λ c = 1 the magnetization in the x direction disappears. For all values λ > λ c the ground state of the system is a (ferromagnetic) ordered state in which the spin-spin interactions are able to prevail over the interaction with the external magnetic field, and the system thus has a non-vanishing value for the magnetization in the x direction (which now plays the role of order parameter). By contrast, for 0 < λ < λ c the system is in a disordered state, which is a paramagnet (i.e., with vanishing magnetization in the x direction), and which preserves the Z 2 symmetry of the Hamiltonian, and is hence non-degenerate. For example, in the limiting case λ = 0 the ground state has all of the spins aligned in the z direction (of the external field). Both the ferromagnetic and paramagnetic phases are gapped. It is only at the precise 2QPT point for the infinite (N → ∞) chain, λ = λ c = 1, that the ground state is gapless.
We note that the presence of a ground-state degeneracy in general always implies the need to choose among the possible ground states. For the sake of consistency between regions with degeneracy (e.g., as in the ordered ferromagnetic phase of the current Ising model) and without degeneracy (e.g., as in the disordered paramagnetic phase of the current Ising model), we shall always choose the ground state to preserve all of the symmetries of the Hamiltonian. Thus, for the case of double degeneracy, as here, this choice simply implies a ground state, in the limiting case λ → ∞ of the pure Ising model, for example, of the Greenberger-Horne-Zeilinger (GHZ) type [72], both in the thermodynamic limit and for finite values of N , which we discuss in more detail below.
Consider the case of N finite first. Thus for zero external field (i.e., λ → ∞) the ground state is doubly degenerate, with the two fully aligned states having values M x = ± 1 2 of the magnetization. This degeneracy is then lifted by the addition of the external magnetic field, as in Eq. (8) with λ finite, so that the new ground state is symmetric, and is unchanged under the transformation σ x i → −σ x i , such that M x = 0, while the first excited state is antisymmetric. It has been shown [70] that these two states have a gap of order λ −N for λ > 1, which hence disappears rapidly as N → ∞. Thus, in the thermodynamic limit (N → ∞), the ground and first excited states become degenerate, allowing a non-zero or- der parameter. Conversely, for λ < 1, the ground state remains non-degenerate in the thermodynamic limit, and no order appears (i.e., M x = 0).
As stated previously, our main aim here is to demonstrate the utility of our generalized Wigner function formalism to detect, characterize and distinguish QPTs of various kinds, including the 2QPT that is exhibited by the present model, as discussed above. The finite size (N = 6) that we choose for our calculations of the system below enables exact solutions to be calculated through direct diagonalization of its Hamiltonian and the subsequent evaluation of the density matrix and the complete set of reduced density matrices (as shown in Fig. 1). Figure 3 shows the behavior of the N = 6 transverse Ising model as the coupling parameter λ is varied in terms of equal-angle slices of the GWF for the correlation functions shown in Fig. 1. We choose to take the equal-angle slice of the GWF ( as this enables appreciation of high-dimensional systems while also remaining invariant under particle exchange [55]. Clear discontinuities are not seen in the phase line plot of Fig. 3 for the transverse Ising model due both to the finite system size and the nature of 2QPTs. A continuous gradual change in the value of the GWF for each correlation function is seen across the range with higher-order correlations experiencing a larger shift. It has been shown elsewhere that exploration of the first derivative of the GWF can witness and characterize 2QPT critical points [62]. Accordingly, we plot the first derivatives with respect to λ of ρ tot and ρ 12345 in the insert to Fig. 3. Both derivatives experience a minimum at λ ≈ 0.9, highlighting the presence of the 2QPT pseudocritical point for the finite system. We note that all of the other correlation functions shown in Fig. 3 experience a very similar minimum in their first derivatives. The slight deviation from the minimum away from the infinite-chain critical point λ c = 1 is almost certainly due to the finite size of the system [33,73]. A similar finite-size scaling behavior is also witnessed in the transverse-field XY model discussed below in Sec. III B. Interestingly however, the pseudo-critical point of the current transverse-field Ising model seems to occur at decreasing values as the size of the system decreases, rather than at the increasing values seen in the corresponding transverse-field XY model. This scaling behavior is worthy of further exploration. Figure 4 shows the spin Wigner plots at points of interest in the transverse-field Ising model. Figure 4(a) shows the state of the system at the initial value λ = 0. All correlation functions in the system clearly correspond to a state with all spins aligned along the z axis [and see Figs. 2(a) and (b)] at this point. Figure 4(b) shows the state of the system at the critical point λ c = 1 of the 2QPT in the infinite chain, at which point the finite chain should be showing signs of the disordered paramagnetic phase in the infinite chain, according to the results shown in the insert to Fig. 3, and as discussed above. Thus, the GWF plot for ρ tot in Fig. 4(b) clearly witnesses an x-axis aligned GHZ type state [see Fig. 2(g) for the corresponding GHZ state in the z direction], with all of the lower-order correlations also presenting features of the GHZ state.
Finally, Fig. 4(c) shows the state of the system at λ = 10. The GWF plot for ρ tot now clearly shows a more uniform x-axis aligned GHZ state, with the southpole region at a higher amplitude than for the previous case λ = 1. ρ 12345 continues to present features of the GHZ state now with lower amplitude nodes of negativity. Lower order correlations no longer present clear features of the GHZ state with many tending towards a mixed state.
The Hamiltonian of Eq. (8) is stoquastic (i.e., all offdiagonal matrix elements are non-positive) in the basis {|→ = (|↑ + |↓ )/ √ 2, |← = (|↑ − |↓ )/ √ 2} for λ ≥ 0; therefore the ground state can be written as a positive linear combination in this basis. For finite N the ground state evolves from the z-axis-aligned state at λ = 0, to the GHZ state that is a positive superposition of the two x-axis Ising-ordered states, in the limit λ → ∞. So long as we take the limit λ → ∞ before the thermodynamic limit N → ∞, an alternative argument is that as the coupling is increased adiabatically from the value λ = 0, where it contains the This results in a ground state that is fully aligned along the z axis, where the state is fully separable and each spin is in the same state. (b) shows the ground state for λ = 1; in the infinite chain case this is a critical point of the 2QPT. In the finite-chain case here, we see signs of the disordered paramagnetic phase; in phase space, we see this as a gradual shift towards a GHZ state aligned along the x axis. This trend is then continued in (c), and as λ increases further the ground state finally reaches the state GHZ + x at λ → ∞. equally weighted admixture of the 2 N terms indicated in Eq. (15), all but the two fully aligned terms of Eq. (16) will smoothly disappear as λ → ∞ and, by symmetry, their phase relationship will be retained. The entanglement therefore increases with λ; in the GHZ limit given by Eq. (16) all reduced density matrices are mixtures 1 2 (|→ · · · → → · · · →| + |← · · · ← ← · · · ←|). Accordingly, only the full Wigner function will show interference terms, as seen in Fig 4(c).

B. Spin-1/2 Transverse-Field XY Spin Chain Model
The spin-1 2 transverse-field anisotropic XY model on a 1D chain with periodic boundary conditions is a generalization of its Ising counterpart that we have discussed above in Sec. III A. Its Hamiltonian is given bŷ where γ is now the spin anisotropy parameter, and the remaining parameters are exactly as in the Hamiltonian of Eq. (8) for its corresponding Ising limit counterpart, and to which Eq. (17) reduces for the special case γ = 1. Once again, periodic boundary conditions are assumed, such thatσ α N +1 =σ α 1 , α = x, y, z. Just as in Sec. III A for the special case γ = 1 of the transverse-field Ising model, we again set the external field strength to h = 1 in order to set the overall energy scale and to simplify the parameter set. We also restrict ourselves to the case where the nearest-neighbor spin interactions are wholly ferro-magnetic, such that λ > 0 and −1 ≤ γ ≤ 1. (Clearly, it is actually sufficient to restrict ourselves to the range 0 ≤ γ ≤ 1, since the Hamiltonian is invariant under the replacementsσ x i ↔σ y i , i = 1, . . . , N ; γ → −γ). The special case γ = 0 is simply the transverse-field isotropic XY (i.e., the transverse-field XX) model.
The Hamiltonian of Eq. (17) again has a Z 2 symmetry associated with the fact that it is invariant under spin rotations of π about the global spin-space z axis, under whichσ µ i → −σ µ i , µ = x, y ; i = 1, . . . , N . It thus commutes with the spin parity operator defined aŝ Hence, all non-degenerate eigenstates of the Hamiltonian of Eq. (17), including the ground state, are also eigenstates ofP z with eigenvalues ±1. At the XX isotropic limiting point (γ = 0) of the model, the Hamiltonian has a much larger symmetry group, since it is now invariant under rotations by an arbitrary angle about the global spin-space z axis. The transverse-field spin-1 2 anisotropic XY model on the 1D chain with periodic boundary conditions described by Eq. (17) plays an important archetypal role in quantum many-body theory and quantum statistical mechanics for two important reasons. First, it is one of the relatively rare models for which an exact analytic solution is known, as described below. Second, the model provides a good approximation to a variety of real physical systems, which can hence be used to simulate it. Examples include ultra-cold neutral atoms loaded onto an optical lattice [20], and a quantum circuit that processes log N qubits, which simulates the model with N spins (or qubits) [74].
We note that the model, like its special-case (γ = 1) Ising counterpart, is also exactly solvable both for all finite values of N and in the thermodynamic limit (N → ∞) by exactly the same sequence of transformations as discussed in Sec. III A. In the case of zero external field the solution was first obtained by Lieb, Schultz, and Mattis [75]. This method of solution is easily generalized to the case of nonzero external field for any finite value of N (see, e.g., Ref. [76]), and the exact solution in the thermodynamic limit (N → ∞) has also been discussed separately [77,78]. Of course, for the small-chain (N = 6) results presented here, direct diagonalization of the Hamiltonian is readily performed, and the groundstate density matrix and subsequent reduced matrices of the system can then be calculated, In the λ-γ plane, in our region of interest (viz., when h = 1, λ > 0, |γ| < 1), the phase diagram of the system in the thermodynamic limit is separated into three regions by the lines λ = 1 and {γ = 0, λ ≥ 1}. The line segment {γ = 0, λ ≥ 1} simply demarcates the anisotropic phase transition from a ferromagnetic phase with magnetic ordering in the spin-space x direction (for γ > 0) to one ordered in the spin-space y direction (for γ < 0).
In the latter paramagnetic phase (for λ < λ c ) the ground state is non-degenerate, and the energy spectrum shows a gap to the first excited state for all values of γ, By contrast, in the ordered phases (for λ > λ c ) the system has a doubly degenerate ground state with a gapped energy spectrum for all nonzero values of γ, while precisely at the isotropic point γ = 0 the ground state is non-degenerate and the energy spectrum is gapless. The spectrum is again gapless precisely on the critical line λ = λ c . Clearly, the point (λ = 1, γ = 0) is a multicritical point in the quantum phase diagram of the model. We note that it has also been shown [33,73] that finite-sized chains display a pseudo-critical point λ c (N ), which deviates from the true critical point λ c = λ c (∞), such that λ c (N ) > λ c for all finite values of N , and λ c (N ) → λ c as N → ∞. With no loss of generality, as explained above, we henceforth restrict ourselves to values of the anisotropy parameter in the range 0 ≤ γ ≤ 1.
Once again, as in Sec. III A, we note that, in order to be consistent between regions with degeneracy (e.g., as in the ordered ferromagnetic phases of the current anisotropic XY model for γ = 0) and without degeneracy (e.g., as in the disordered paramagnetic phase of the current anisotropic XY model), for all our calculations below we always choose the ground state to preserve all of the symmetries of the Hamiltonian.
The Hamiltonian of Eq. (17) also possesses another remarkable property. Namely, for any value of the anisotropy parameter in the range 0 < γ ≤ 1 and for all values of the system size N , there exists a ground-state factorization point λ f = λ f (γ), given by for which the system contains two degenerate ground states, both of which are fully factorized in the sense that they are simple products of single-site states [68,76,[80][81][82]. Clearly, for these states themselves entanglement vanishes in principle. However, one must take note that these factorized states themselves are not eigenstates of the spin parity operatorP z , as we now discuss more fully. For the case 0 < γ ≤ 1, one can rather readily show [76] that one of the two factorized states, which we denote by |ϑ , has all of the spins fully aligned in the x-z spin-space plane and at at an angle ϑ to the spin z axis, given by The other factorized state, degenerate in energy with |ϑ , is simply |−ϑ =P z |ϑ . The two states |±ϑ are nonorthogonal for ϑ = π 2 (i.e., for γ = 1). In that case the correct orthonormal basis that conserves spin parity is spanned by the two entangled states |ϑ ± , defined by which satisfyP z |ϑ ± = ± |ϑ ± . We note that an equalangle slice of the reduced GWF, ρ I (θ, φ) ≡ W EA I (θ, φ), for a similar wave function to that of Eq. (21), viz., one with ϑ = π 4 , such that the two aligned states have a relative angle π 2 between their alignment directions, is shown in Fig. 6(f) of Ref. [63].
The states |ϑ ± , as defined in Eq. (21), are the actual eigenstates of the HamiltonianĤ of Eq. (17) in each spin parity subspace at the factorization point λ = λ f , and are just the corresponding limits of the exact definiteparity ground eigenstates |Ψ ± (λ) as λ → λ f . Thus, the factorization point λ f represents a crossing point of the two lowest opposite-parity levels [76], and is hence a parity transition point. We note that for the limiting value γ = 0 (i.e., when the Hamiltonian becomes isotropic), for which λ f = 1 the alignment angle ϑ = 0 and the groundstate degeneracy vanishes, in accordance with our discussion above. Indeed, the state |0 is a trivial eigenstate of the Hamiltonian of Eq. (17) in this isotropic limit for all values of the parameter λ. In the opposite limiting case γ = 1, which is just the Ising model in a transverse field, the alignment angle ϑ = π 2 , and λ f becomes infinite. In this case the overlap − π 2 π 2 = 0, and, as expected, the states π 2 ± simply become x-axis aligned GHZ states, GHZ ± x , one of which has been explicitly displayed in Eq. (16) above for the case of the positive superposition of the two x-axis Ising-ordered states. These states are just the counterparts of the z-axis aligned states, GHZ ± z , the equal-angle slices of the reduced GWFs, ρ I (θ, φ) ≡ W EA I (θ, φ), of which have been shown in Figs. 2(g) and (h).
Such GHZ states, although globally entangled, display no two-spin entanglement (for N > 2). As an aside here we note that the entanglement properties of the N = 2 version of the model have also been studied in great detail [35]). By contrast, for 0 < ϑ < π 2 , the two-spin entanglement in the states |ϑ ± depends critically on the nonzero value of the overlap −ϑ|ϑ . Indeed, the two-spin entanglement has been shown to undergo a change as λ is varied across λ f , being zero exactly at the factorization point. For this reason this point has also become known as the entanglement transition point [68], although its real origin lies in the fact that it is a point of accidental spin-parity symmetry breaking where two branches of eigenstates with different spin parities cross one another.
We note finally that, for the system in the thermodynamic limit, at all points in the ordered phase (for which λ > λ c = 1), which includes the factorization point, the ground state is exactly doubly degenerate, and hence the factorization point cannot be expected to be a prominent feature. Similarly, the (e.g., entanglement) effects are small in large anisotropic chains, for which the ground states in each of the two spin-parity sectors are nearly degenerate in energy. However, as we shall see explicitly below, the factorization point is quite prominently visible for small finite cyclic chains (here with N = 6). Indeed, they have also been shown to remain appreciable for increasingly larger values of N as the size of the anisotropy parameter γ is accordingly decreased [76]. Figure 5 shows the behavior of the XY model as the coupling strength λ is varied for a fixed value of the anisotropy parameter γ = 0.5, in terms of equal-angle slices of the GWFs. Phase lines for both the full system correlation function and the complete set of sub-system correlation functions (shown in Fig. 1) are shown.
Abrupt changes in the ground state of the system are clearly seen from Fig. 5 at both λ = λ f (γ = 0.5) ≈ 1.155, which corresponds to the factorization point, and λ ≈ 1.545, which represents the pseudo-critical point in the N = 6 system corresponding to the critical point of the 2QPT at λ = λ c = 1 in the infinite chain. Discontinuities at the factorization point λ f are seen in all correlation functions of the system, with higher-order correlations better capturing the transition. This highlights the fact that the factorization point coincides with an energy level crossing of the ground state and first excited state [68,83,84]. The mean of the two values of W EA (θ, φ) just before and just after the factorization point is precisely that of an equal mixture of the two (positive and negative) parity states |ϑ ± . Thus, the interference terms precisely cancel there, giving the corresponding value of the state |ϑ , viz., W EA (0, 0) = 2 −n (1 +  Fig. 4, where here (a) is also an eigenstate ofσz; resulting in a pure, separable state. This can be seen by considering the reduced states, especially the single-spin state. (b)-(d) show a steady loss in the purity (negative values) of the single-spin state as λ increases. This, in addition to the full state being a pure state, is an indicator that the overall state is entangled. In the 2-, 3-, 4-, and 5-spin functions, we can see the spin-coherent state from λ = 0 state split in two and move closer to the equator at orthogonal points on the Bloch sphere. At λ = 4 the two coherent states are distinct and almost orthogonal -it is clear that this approaches a GHZ-type state, which is the case as λ → ∞. Note that (a) and (b) follow similarly to the transverse-field Ising model in Fig. 4. After this point there is a parity change; similar to the difference between the GHZ states in Figs. 2(g) and 2(h), the shift from the state in (b) to (c) shows a rotation of the interference terms around the x axis. This results in a larger negative value at ρtot(0, 0), which is reflected in the phase transition in Fig. 5. Note that the correlations in ρ12345 indicate that the state is not yet a perfect GHZ state, for which all reduced states are equal statistical mixtures of |↑ ⊗n and |↓ ⊗n . A second discontinuity is seen in many of the system correlation functions at λ ≈ 1.545, which coincides with the 2QPT pseudo-critical point. Higher-order correlations such as ρ tot and ρ 12345 experience abrupt jumps in their phase lines. Nearest-neighbor two-site and threesite correlation functions, ρ 12 and ρ 123 respectively, also present clear discontinuities in their phase lines, highlighting their utility for witnessing 2QPTs in larger systems [33,34,69]. Despite this, the 2QPT is less clearly seen in the ρ 1234 and ρ 135 correlation functions. The fact that ρ 1234 is unable to witness the 2QPT clearly is somewhat surprising as higher-order correlations such as the five-and six-site correlation functions present the clearest indications of the 2QPT. Figure 6 shows the spin Wigner plots at points of in-terest across the phase plot. Figure 6(a) shows the initial state of the system at λ = 0. Here the system is in a state with all spins aligned in the z direction of the external field. This is clearly reflected in all the spin Wigner plots as they match the aligned spin reference plots shown in Figs. 2(a) and (b).
As λ is increased, the total state of the system ρ tot gradually transitions towards a state of the GHZ + x type. The lower-order correlations, while only able to capture some features of the GHZ state, display how the spins migrate from a coherent state oriented at the north pole and split into two coherent states aligned in the z-x plane, moving in opposite directions. For the λ = 4 plot in Fig. 6(d), we can see the two coherent states more clearly as they approach the x axis, anti-aligned on opposite sides of the equator. The higher-order correlations, such as ρ 12345 and ρ tot , clearly start to show oscillating positive and negative values -characteristic features of a GHZ state. Note that the existence of correlations be-tween the two coherent states in ρ 12345 shows that the state is not a perfect GHZ state, which only happens in the limit λ → ∞.
As with the transverse-field Ising model, it is clear that this is a gradual change to the GHZ + x state; the difference here is at the factorization point λ f that causes an abrupt change in the symmetry of the system, which can be seen in Fig. 6(c). The state of the system then transitions abruptly from a positive-parity state to a negativeparity state. These positive-and negative-parity states are simply the superposition of the two coherent states that are evident in the lower-order correlations, with either a positive phase or negative phase, which at the limit where the two states aligned along the x axis is similar to the difference between a GHZ + x -type state and a GHZ − x -type state. This behavior is expected and has been documented in previous research [68,83,84].
Finally, at λ = λ c we witness the system transition back to a positive-parity state as a critical point of the 2QPT is passed. The actual critical point of the system cannot be determined with any certainty as the factorization point obscures this area of the phase line. Further exploration of this phenomenon and its scaling properties as a function of system size could be fruitful.

C. Spin-1/2 XXZ Spin Chain Model
The Hamiltonian for a cyclic, spin-1 2 anisotropic 1D XXZ Heisenberg model with nearest-neighbor interactions is given bŷ where J is the coupling strength and ∆ is the anisotropy parameter [2,69,85]. As before,σ α i , α = x, y, z are the usual Pauli matrices, and periodic boundary conditions are assumed, so thatσ α N +1 =σ α 1 , α = x, y, z. Although the XXZ Hamiltonian of Eq. (22) contains two parameters, the essential physics is captured by just one, as we now show.
For simplicity we restrict attention to the case when N is even and, by analogy with Eq. (18), we consider the unitary operatorÛ z , defined as follows, Considered as a similarity transformation, its mode of action is to perform a rotation by π about the global spin-space z axis of the spins on the even sites, under whichσ µ 2m → −σ µ 2m , µ = x, y ; m = 1, . . . , N 2 , with all other spin components unchanged. Clearly, it leaves the underlying SU (2) algebra unchanged. However, its mode of action on the XXZ Hamiltonian of Eq. (22) is thus From Eq. (24) we immediately conclude that it is only the relative sign between the parameters J and ∆ that is relevant, rather than their separate signs. Thus, with no real loss of generality we may take J > 0 so long as we also consider ∆ over the full range −∞ < ∆ < ∞. Thus, henceforth we put J = +1 to set the overall energy scale.
Only for the isotropic case when ∆ = 1 does the XXZ Hamiltonian of Eq. (22) preserve the full SU (2) symmetry of arbitrary rotations in spin space. However, it is easy to show that the Hamiltonian commutes with the z component of the total spin operator, which then leads to the smaller U (1) symmetry under arbitrary rotations about the global spin-space z axis, which are generated by the rotation operator, The mode of action ofR z (φ) on the basic spin operators is as follows, and it thus leaves the XXZ Hamiltonian of Eq.
To set the scene for the quantum spin-1 2 case, let us first consider the classical version of the XXZ model Hamiltonian of Eq. (22) with J > 0, as used here, and for the limiting case N → ∞. For ∆ < −1 the system is a ferromagnet aligned along either of the ±z spin-space directions. At a first critical point, ∆ c1 = −1, the system undergoes a first-order transition to a Néel antiferromagnet with the axis of alignment in an arbitrary direction in the spin-space x-y plane. This so-called easy-plane antiferromagnetic phase survives over the region −1 < ∆ < 1 until, at the second critical point, ∆ c2 = 1, the system undergoes a further first-order transition to a so-called easy-axis Néel antiferromagnet with the axis of alignment along the spin-space z direction for ∆ > 1. The classical ground-state is thus doubly degenerate for |∆| > 1 and infinitely degenerate for |∆| ≤ 1.
By writing the first two terms in the Hamiltonian of Eq. (22) are the usual Pauli spin raising and lowering operators, it is easy to see that the two spin-space z-aligned ferromagnetic states are always exact eigenstates of this Hamiltonian. By contrast, we also see at once that the corresponding z-aligned Néel antiferromagnetic states are not eigenstates. For all values ∆ < −1 these two ferromagnetic states form the doublydegenerate ground state, just as in the classical case. The order parameter is the absolute value of the z component of the magnetization, which assumes its maximally saturated value |M z | = 1 2 throughout the region ∆ < −1, over which the energy spectrum is gapped.
Precisely at the value ∆ = −1, the order parameter M z drops to zero discontinuously, and we have a 1QPT. Throughout the ensuing region −1 < ∆ < 1 the spectral threshold is also gapless. The quantum fluctuations completely destroy the classical easy-plane Néel order that exists there. Indeed, there is no long-range order present in the spin-1 2 model in this region, and accordingly the ground state also becomes non-degenerate there. Such ground-state phases that are completely devoid of all long-range order are called critical phases. We note that this critical phase also contains the very special point ∆ = 0, which can be exactly mapped by a Jordan-Wigner transformation onto free (lattice) fermions.
In this planar-regime critical phase with −1 < ∆ < 1 all equal-time two-spin correlation functions decay algebraically with separation distance, with an exponent that depends on the value of ∆. It is only at the exact limiting value ∆ = −1 + (with a positive infinitesimal) that the ground state again breaks the symmetry of the Hamiltonian. It is a doubly-degenerate state with a fully saturated value of the staggered magnetization (i.e., as defined on either of the two sublattices of even or odd sites) pointing in a direction perpendicular to the spin-space z axis, and with no correlated fluctuations. A priori , we do not expect to see features of such critical phases in small finite-sized systems far removed from the thermodynamic limit (N → ∞) [12,13].
At ∆ c2 = 1, the first-order classical transition now becomes a ∞QPT of the Berezinskii-Kosterlitz-Thouless type, and the ground state reverts to being doubly degenerate, such that for ∆ > 1 antiferromagnetic long-range order gradually develops as ∆ increases. The (Néel) order parameter for the region ∆ > 1 is the staggered magnetization, m z , in the spin-space z direction, which was first calculated by Baxter [92]. He showed explicitly that m z develops smoothly and monotonically as a function of ∆ from a value equal to zero for ∆ = 1, only reaching its fully saturated (classical) value in the limit ∆ → ∞. Similarly, the whole region ∆ > 1 has an energy threshold that is gapped, with an energy gap that also increases smoothly and monotonically from zero at ∆ = 1 as ∆ is increased further.
As previously stated, the finite size (N = 6) of our system enables exact solutions of the ground state to be calculated through direct diagonalization. In turn, the system's density matrix and complete set of reduced density matrices (see Fig. 1) can then be calculated. Figure 7 shows the behavior of the XXZ model as the anisotropy parameter ∆ is varied. Once again the equal-angle slice of the GWF has been calculated for the full system correlation function and the constituent subsystem correlation functions. Spin Wigner functions are plotted at points of interest across the phase diagram in Fig. 8.
A clear discontinuity is seen in the phase lines at the 1QPT critical point ∆ c1 , with all of the system's correlation functions exhibiting an abrupt jump. We note parenthetically that if the symmetry had been preserved by using the state GHZ + z with an equal mixture of ferromagnetic up and down states, the discontinuity in ρ 1 at ∆ c1 would presumably disappear. In general, higherorder correlations are better able to capture the behavior of the system at this critical point. Although past research was unable to witness 1QPT critical behavior in the single site GWF [62], in our findings a clear discontinuity at the 1QPT critical point is seen in ρ 1 . An explicit maximum or minimum at the ∞QPT critical point, ∆ c2 , is not seen however.
Despite this, interesting behavior is captured in many of the system's correlation functions across the −1 < ∆ < 10 range. This type of behavior has been captured As with the other spin chains, the first ground state is the pure, separable state where all spins are aligned along the z axis, which can be seen in (a). This model differs in that the ground state, ignoring the degeneracy, stays as this state until the phase transition at ∆ = −1, this is also evident from Fig. 7, where all the correlation functions are constant until the phase transition. (b) shows the GWFs for the state just after the phase transition; it demonstrates how the ground state after the phase transition becomes highly entangled. This is clear from noting that the single-qubit state is a fully mixed state while the full six-spin function is a pure state. Note that all subsequent states also have a fully mixed state in the single-spin Wigner function. They then all differ in how the entanglement is distributed throughout the state, which can be seen from the 2-, 3-, and 4-spin reduced Wigner functions. For instance, in (c) we have a completely uniform distribution in all equal-angle slices. For ∆ = 1 all equal-angle GWF are isotropic as the Hamiltonian is invariant under spin rotation and has a singlet ground state. The two-spin cases ρ 1k are Werner states. In (b) and (d) there are more features evident in the reduced states. This is because they also have contributions from the triplet state |Ψ+ shown in Fig. 2(d). As ∆ increases, in (d) we can see in some of the reduced states, such as ρ13 and ρ135, that this is similar to a GHZ-state state in that there are two antipodal coherent states on the Bloch sphere. We can further see that that the full state tends towards that in Fig. 2(j).
in past research [62]. ρ 1 , ρ 12345 and ρ tot all take constant values across the range. ρ 12 experiences a decrease in the value of the GWF across the range. The first derivatives of ρ 13 and ρ 135 experience a minimum at ∆ ≈ 1, while ρ 124 and ρ 1235 both experience a maximum in their first derivative at this point. This highlights how exploration of the first derivatives of the GWF can be used to witness the critical point of a ∞QPT. ρ 14 and ρ 124 both experience a maximum at ∆ ≈ −0.3, while ρ 1245 experiences a minimum at ∆ ≈ −0.3. The first derivatives of ρ 123 and ρ 1234 decrease across the range. It has been shown that through extremization procedures certain correlation measures can better witness ∞QPTs [32,68,69,95]. Past research has implemented these for the GWF with much success [62], and such further exploration of these techniques could be useful.
At ∆ c1 we have a 1QPT between a non-degenerate ground state and a two-dimensional ground-state space. Unlike in the case of the transverse Ising model, there is no adiabatic continuity criterion to select a vector in this space, so, for example, neither the state |↑↑ · · · with all spins aligned in the upwards z direction nor the state GHZ + z ≡ 1 √ 2 (|↑↑ · · · + |↓↓ · · · ) has any special claim. The unbiased choice for ∆ < −1 is presumably then a mixtureρ = 1 2 (|↑↑ · · · ↑↑ · · ·| + |↓↓ · · · ↓↓ · · ·|). This will have the same reduced density matrices as the state GHZ + z . Nevertheless, for simplicity Fig. 8 shows the state |↑↑ · · · . Thus, Fig. 8(a) clearly shows the state of the system at ∆ = −2 to be in an aligned ferromagnetic state; this can be inferred due to the resemblance of the state's spin Wigner plots to the reference plots of Fig. 2(l) for the ferromagnetic mixed state. Figure 8(b) shows the state of the system at the 1QPT critical point ∆ = −1 + . The state of ρ tot and ρ 1235 show similarities to the entangled, antiparallel aligned state seen in Fig. 2(j), while the state of ρ 1234 and ρ 1245 present similarities to the anti-aligned state in Fig. 2(c). On the other hand, ρ 123 and ρ 124 are both in a low amplitude singlet-type state [see Fig. 2(e)], while ρ 135 and ρ 13 both show similarities to the Bell state of Fig. 2(d). The states of ρ 12 and ρ 14 possess similarities to the antialigned entangled state in Fig. 2(f). The abrupt emergence of these states at the 1QPT critical point enable us to witness the 1QPT and infer the presence of an antiferromagnetic type phase. Figure 8(c) shows the state of the system at the ∞QPT critical point ∆ = 1. Here all equal-angle Wigner functions are constant over the Bloch sphere as the Hamiltonian is an isotropic antiferromagnet. The signs of the two-spin reduced Wigner functions show that the nearest-neighbor correlation (represented by ρ 12 ) is negative.
The states that have a constant value over the Bloch sphere for the equal-angle slice in Fig. 2 where |Ψ − = 1 √ 2 (|↑↓ − |↓↑ ) and x ∈ [0, 1]; note that ρ(x) is only entangled when x > 1 3 [96]. The equal-angle slice for the Wigner function of a Werner state is and is only an entangled state when the equal-angle slice has a negative value. We see in Fig. 8(c) that the value of the Wigner function changes as the distance between the spins changes, with explicit values of ≈ 0.623, is seen here. For ρ 13 and ρ 135 , we have odd-number indices for the spins, picking out every other spin from the full state. This, in effect, gives us an equal statistical mixture of |↑ and |↓ states. We can then see in the Wigner functions two coherent states, one on the north pole and the other on the south pole; the two-and threespin equivalents to Fig. 2(l).
In ρ 12 and ρ 14 we instead pick out anti-aligned features from the full state. We therefore tend towards an equal statistical mixture of |↑↓ and |↓↑ in both cases. In the equal-angle slice, both |↑↓ and |↓↑ have the same Wigner function representation. In addition, the Wigner function of the statistical mixture of two states is the same as the addition of the two Wigner functions. This statistical mixture then produces a Wigner function reminiscent of Fig. 2(c) for ρ 12 and ρ 14 . Similar results are seen in the three-spin functions, ρ 123 and ρ 124 where we have a statistical mixture of a state with two |↑ and one |↓ with a state that has one |↑ and two |↓ spins.
As we move onto the four-spin states, the same effect can be seen, where we have the equal statistical mixture of |↑↓↑↓ and |↓↑↓↑ in ρ 1234 and ρ 1245 . Conversely, ρ 1235 differs as it doesn't have the same number of |↑ and |↓ spins in each side, which can been seen in the spin indices in the full state.
This all then paints a picture that we have a GHZtype state that has every other spin flipped around the Bloch sphere, further indicating that the system is in an antiferromagnetic phase at ∆ = 10. This behavior has been recorded in past research [85,93]. In contrast to phase line techniques, spin Wigner function visualizations were able to witness the ∞QPT exactly at the critical point without the use of extremization procedures. This highlights the method's utility in witnessing features of a system that are less easily discerned through other techniques.

IV. CONCLUSIONS
We have demonstrated the utility of the generalized Wigner function formalism in exploring quantum critical behavior in a range of spin chain models. The generalized Wigner function's ability to explore the ground state of many-body systems, while also capturing the quantum correlations present within them, has been shown to enable first-, second-and infinite-order quantum phase transitions to be witnessed and characterized in a range of systems. In addition to this, its generalization to discrete, finite-sized spin systems has been shown to allow for exhaustive exploration of each model's correlation functions, which in turn has enabled a much deeper appreciation of the critical behavior of these spin systems. We also demonstrated the GWFs ability to determine the state of the system through the use of spin Wigner function visualization techniques and qubit state reference plots. These reference plots enable us reliably to infer the state of the system through observation alone, making this method much more intuitive and accessible than other leading techniques. In addition to this, the spin Wigner plots highlight important features of secondand infinite-order QPTs that fixed angle phase lines are either unable or less successful at capturing.
To further explore the utility of this method for witnessing second-and infinite-order QPTs, extremization of the GWF through careful choice of angles (θ, φ) should be explored, as this has shown to be useful in characterizing higher order QPTs [32,64,68,69,95]. The properties of ground-state factorization were also explored in the XY model with the phenomena being captured in a wider range of correlation functions than past research in this area [62]. We were able to capture both the entanglement transition and the characteristic change in symmetry of the system through the GWF alone.
One interesting area of inquiry left to future research is the obscuring of the 2QPT critical point by the XY model's ground-state factorization. The length over which the state of the system remains in the factorized state and its scaling with system size should be studied. This work also highlights the continued utility of finite spin systems in exploring quantum critical behavior as, despite their simplicity, they are able to produce a rich and insightful look into critical spin chains and manybody systems as a whole [67-69, 73, 80]. In addition to this, scaling of finite critical quantum systems is well captured by our method, thereby allowing for the properties of this phenomenon to be explored further. The ability of phase-space techniques to provide in-depth and insightful analyses of many-body systems highlights their utility for better understanding more exotic states of matter present in this field of study [97][98][99]. Our method also provides a foundation upon which future research can apply phase-space techniques to discern the factors that lead to the emergence of different states and phenomena in these systems. The fact that the Wigner function can be reconstructed from experimental data highlights this method for implementation in the experimental investigation of quantum phase transitions [46].
On a final note, is is interesting to compare the current GWF formalism with other widely used quantum many-body theory techniques that also have a similar phase-space underpinning. It is widely acknowledged that the coupled cluster method (CCM) [100,101] nowadays provides one of the most pervasive, most powerful, and most accurate at attainable levels of implementation of all fully microscopic formulations of quantum manybody theory. Thus, the CCM (in both its so-called normal and extended versions) provides an extremely convenient parametrization of the Hilbert space, which also exactly maps it onto a classical Hamiltonian quantum many-body/field theory in a complex symplectic phase space [102][103][104][105].
In CCM formalisms one constructs explicitly an energy functional that variationally determines both the ground-state wave function and the dynamical equations of motion for non-stationary states. The equations of motion thus take the familiar classical canonical Hamiltonian form in some well-defined many-body configuration space. The expectation-value functional of an arbitrary operator is also manifestly constructed in a way that, very importantly, automatically preserves the Hellmann-Feynman theory exactly at all natural levels of truncation in the configuration space [101,106]. Another key feature of the CCM techniques is that they also automatically satisfy the Goldstone linked-cluster theorem at all levels of implementation [100,101], so that the thermodynamic limit (N → ∞) can be taken from the very outset, thereby obviating the need for any finite-size scaling of the results, as is needed in most other many-body calculations using alternative techniques.
In the present context we note that the CCM has been very extensively applied with great success to a plethora of systems in quantum magnetism, at high levels of implementation in a systematic hierarchy of approximations that are tailor-made for systems confined to a regular spatial lattice in any number of dimensions [107]. Among many others, it has been applied in particular to the transverse Ising model [108], the XY model [109], and the XXZ model [110,111], both for the spin-1 2 chains considered here as well as for higher spin values and on specific lattices in two or more dimensions. It will surely be of considerable future interest to explore possible synergies between the CCM and the current GWF techniques in these and other applications to spin-lattice problems.
One particular area of such mutual interest is to extend the discussion to the dynamical phase transitions [112] that can occur when a quantum many-body system is driven out of equilibrium, e.g., by a quantum quench. After such a quench the long-time steady state can display a symmetry-broken phase, with singular features at the transition to the disordered phase. If the energy of the system is thereby shifted across a symmetry-restoration threshold and the system thermalizes, such a transition may be thought of as occurring in the microcanonical ensemble. Alternatively, non-ergodic systems can display long-time steady states that are unable to be described via any of the conventional ensembles of statistical mechanics. In such cases one can generate phases and phase transitions that are unable to be found in any equilibrium situation [113].
Spin-1 2 transverse-field Ising models provide a paradigmatic example of such systems [112], the dynamical phase transitions displayed by which have recently been studied [114] by means of a discrete truncated Wigner approximation [54], which shares some features with our own GWF formalism utilized here. It would hence surely be of considerable interest to investigate such dynamical phase transitions using the techniques discussed here, possibly in conjunction with additional input inspired by the CCM.