Fisher zeros and persistent temporal oscillations in non-unitary quantum circuits

We present a quantum circuit with measurements and post-selection that exhibits a panoply of space- and/or time-ordered phases, from ferromagnetic order to spin-density waves to time crystals. Unlike the time crystals that have been found in unitary models, those that occur here are \emph{incommensurate} with the drive frequency. The period of the incommensurate time-crystal phase may be tuned by adjusting the circuit parameters. We demonstrate that the phases of our quantum circuit, including the inherently non-equilibrium dynamical ones, correspond to complex-temperature equilibrium phases of the exactly solvable square-lattice anisotropic Ising model.


I. INTRODUCTION
For a many-body quantum system with Hamiltonian oper-atorĤ, there is an evident formal similarity between the unitary time-evolution operator, e −iĤt/ , and the density operator for a thermal equilibrium state, e −βĤ . Since the 1950s this has led to very fruitful cross-fertilization between the theory of quantum dynamics and the equilibrium statistical mechanics of quantum systems. Perhaps the most influential of these is the Matsubara formalism 1 , where the thermal density operator is regarded as an evolution operator in imaginary time: this allows many of the tools of diagrammatic perturbation theory to be copied more or less directly from the dynamical to the statistical case.
In recent decades, the development of the theory of open quantum systems has led to a broadening of interest on the dynamical side of the dynamical-statistical correspondence, since interactions between the quantum system of interest and its environment generically induce (effectively) non-unitary evolution. The quantum circuits in which we shall be mainly interested here exhibit many-body mixed dynamics, with unitary evolution interrupted by projection operations meant to model measurements by a classical environment. Crucially, the many-body system is allowed to continue evolving after such measurements, and displays a host of novel phenomena due to the tunable interplay of non-unitary measurements and the intrinsic unitary dynamics 2- 20 .
There have been parallel broadenings of interest on the statistical side. Starting from the early 1950s, Lee, Yang, Fisher and others [21][22][23][24] pioneered the extension of conventional statistical mechanics to the case where the coupling constants in the Hamiltonian, or even the inverse temperature β itself, are considered to be complex quantities. This opens up the possibility of points in the complex β-plane where the partition function vanishes, something which is not possible for real temperature. For simple models, such as the isotropic zerofield Ising model on the square lattice, these 'Fisher zeros' occur on contours in the complex β-plane which cut the real β-axis at positions corresponding to the critical temperatures of phase transitions in the model. If the density of zeros vanishes as the real β-axis is approached, the transition is continuous; if the density remains finite, the transition is first-order.
Complex-coupling approaches to statistical mechanics have tended to be seen as an essentially formal tool. In light of recent progress in understanding the rich phase structure of non-unitary circuits such as those under continuous measurement, it is timely to revisit canonical statistical mechanics models at complex temperature considered as descriptions of non-unitary evolution. This is the perspective we adopt here.
Classical measurements generically introduce randomness in discrete space-time. By contrast, the canonical Ising model is disorder-free. The correspondence that we shall demonstrate therefore requires post-selection of measurement outcomes [25][26][27][28][29] , i.e. only certain outcomes are allowed to continue evolving and contribute to the eventual disorder-free ensemble of trajectories. Some of the features of the statistical side of our correspondence were anticipated in prior work by some of us 30 , demonstrating the existence of long-range incommensurately modulated correlations underlying Fisher zeros in the thermodynamics of Ising ladders.
We shall show that the dynamics of the corresponding M -qubit circuits exhibit long-range correlations that coalesce, in the M → ∞ limit, into extended ordered regions (see Fig. 1). These can be interpreted as phases of the anisotropic 2D Ising model at complex temperature. They include relatively conventional short-range ordered 'paramagnetic' phases and long-range-ordered ferromagnetic and antiferromagnetic phases, but also somewhat peculiar incommensurate critical phases. These latter phases exhibit spatially and/or temporally modulated correlators with a dynamically determined modulation period untethered from the underlying lattice. At least one of these latter phases bears a phenomenological resemblance to the time crystals recently discussed in the context of unitary dynamics of isolated many-body localized systems 31 . However, it does not fit into the classification presented in that work, since the circuits we consider are nonunitary, and the no-go theorems 32 forbidding time-crystalline order consequently do not apply.
Before turning to our results, we briefly discuss some connections to superficially similar questions discussed in previous literature. Temporally modulated phases in open quantum systems (i.e. limit cycles) have been shown to exist in more than two spatial dimensions 33,34 . The non-unitary quantum circuits we consider can be regarded as Trotterized non-Hermitian Hamiltonians, which have been extensively explored [25][26][27]35 . Unlike these works, we keep the Trotter "timestep" finite, so the models we consider are two-dimensional statistical mechanics models with a transfer matrix that may be contracted either sideways or from top to bottom. In addition, many-body entanglement properties have been computed by contracting transfer matrices sideways in a series of recent works [36][37][38][39][40][41][42][43][44][45] , but primarily in contexts where the dynamics is unitary along one or both directions. In particular, Refs. 38, 43, and 46 have proposed experimental protocols to study the non-Hermitian dynamics of large systems using spacetime duality. On the statistical mechanics side, Ref. 47 used tensor-network methods similar to those we use here 48 to characterize the thermodynamics of the Yang-Lee model. The present work applies tools from the complex-temperature statistical mechanics literature 47,48 to discuss the unexplored physics of spatio-temporal correlations in non-unitary quantum circuits. So far these circuits have primarily been studied for their entanglement properties; we demonstrate here that even their conventional correlation functions can exhibit striking phenomena that would be forbidden by unitarity (in closed systems) or by dimensionality (in open systems 33,34 ).
The remainder of the paper is organized as follows. In Section II we define our quantum circuit by explicitly constructing the local gates necessary to reproduce the complextemperature statistical mechanics of the Ising model. We also discuss observables of interest, both 'thermodynamic' quantities and two-point correlation functions, and construct the transfer matrices that govern the complex-temperature statistical mechanics of interest. In Section III we study the anisotropic complex-temperature Ising model analytically using fermionization and also numerically using a tensornetwork coarse-graining scheme. The fermionization treatment is restricted to the case without an externally applied magnetic field, but the tensor-network coarse-graining approach allows us to identify additional transitions as a function of field strength. We conclude with a summary and a discussion of open problems.

A. Quantum circuit formalism
A quantum circuit consists of a set of qubits, which we may label {1, . . . , M }, successively subjected to operations in the form of quantum gates. In the scenario investigated here, the qubits are Ising spins and the gates are of two types: two-qubit gates W m,m+1 = exp(JZ m Z m+1 ) and single qubit gates the single-step evolution operator T for each temporal slice of the circuit is T = VW, and the evolution operator for the full circuit is T L , where L is the total number of temporal slices. Note that our qubits are arranged in a finite chain vis-à-vis the operator W, rather than a ring with periodic boundary conditions. This is a common context for models of quantum circuits; what is new here is that the parameters J, γ, and A are all allowed to be complex, hence the individual quantum gates are not in general unitary. Acting on a density matrix ρ j , the result of each temporal slice of the circuit is the update ρ j+1 = T ρ j T † . Requiring that Tr ρ j be preserved then imposes a relationship among the complex parameters {J, γ, A}.

B. Single-qubit transfer matrix
At each site, our single-qubit transfer matrix V is a product of a unitary U = exp(iωd·σ) and a POVM (positive operatorvalued measure), where bothd andn are unit vectors on S 2 , and where α = ±1 is the specified measurement outcome. The fact that α P † (n, φ | α) P (n, φ | α) = 1 is what makes P (n, φ | α) a POVM; the fact that there is one such operator for each measurement outcome α means that the measurement is 'efficient'. We choosed =n =x, and we write (3) where A = cos(2φ)/2 and We define the operator Below we shall post-select α m = 1 for all sites, hence The single site transfer matrix V = A exp(γX) is identical to that of the one-dimensional classical Ising model with tanh γ = e −2βJx and A = √ coth γ − tanh γ . Throughout the remainder of this paper we shall set J x ≡ 1. Assuming periodicity in this (temporal) direction, the classical partition function is Z L (β) = Tr V L = cosh L β + sinh L β and the condition Z L (β) = 0 requires tanh β = e (2 +1)πi/L , occurring at L equally spaced values ∈ {0, . . . , L − 1} around the circle | tanh β | = 1 in the complex tanh β plane, as noted by Beichert et al. 30 Interleaved with these Fisher zeros are L points tanh β = e 2πi /L where the correlation function is long ranged, with C(r, β) = cos(2π r/L) (restricting 0 r L), corresponding to a wavevector q = arg tanh β = 2π /L. In the thermodynamic limit L → ∞, the Fisher zeros coalesce into a branch cut along the unit circle, with the free energy exhibiting a simple first-order-like cusp nonanalyticity across the cut. Along the contour of Fisher zeros tanh β = e iθ , with θ real. It follows that e −2β = tanh γ = −i tan(θ/2). Thus the entire unit circle in the complex tanh β plane corresponds to simple unitary stroboscopic precession, i.e. coherent spinflipping. The value θ = 0 corresponds to a perfectly static spin with no flipping and perfect persistence, analogous to the ferromagnetic ground state of the corresponding statistical mechanics problem. The value θ = π corresponds to a complete spin flip with no identity component and no persistence, analogous to the antiferromagnetic ground state.
Values of tanh β that do not lie on the unit circle are associated with circuits that include measurement. Such circuits in the single-qubit case generically exhibit exponentially decaying temporal correlations, just as the analogous complextemperature Ising models exhibit exponentially decaying spatial correlations. However, as we shall see, the multi-qubit case is richer: when M > 1 we shall find 'decoherencefree subspaces' along lines in the complex-temperature plane, though fully unitary evolution occurs only at isolated points.

C. Two-qubit transfer matrix
As noted above, our single step evolution operator is a product of single qubit and two qubit terms. For M > 1 we introduce the two-qubit transfer matrix W m,m+1 , which can be expressed as conditioning a symplectic operation exp(ηJZ m ) on the POVM P m+1 (ẑ, π 4 | η), where η = ±1, for each qubit m ∈ {1, . . . , L − 1}. Explicitly, we have and W ≡ M −1 m=1 W m,m+1 as given in eqn. (1).

D. Mapping to complex-temperature statistical mechanics
The correspondence between statistical mechanics in d + 1 dimensions and quantum mechanics in d dimensions is well known 49 . Traditionally it requires fine-tuning to a critical point of some sort (the so-called τ -continuum limit) to enable passing to continuum time where the correspondence is most powerful. More generally, however, the correspondence is to a discrete-time 'kicked' quantum evolution, of the type exhibited by our quantum circuit.
By construction, the one-step evolution operator of our quantum circuit resembles the transfer matrix of a statistical mechanical system: specifically, an anisotropic 2D Ising model. That model is characterized by its couplings in the xand y-directions, J x and J y respectively, and by its inverse temperature, β. These are functions of the quantum circuit parameters {J, γ, A}. In what follows, we shall consider the subspace of circuit parameters for which the inverse temperature β is complex, while the coupling constants J x and J y are real.
One important question is which direction in our twodimensional statistical model will be identified as the time direction in the quantum circuit. In isotropic models one typically chooses a diagonal direction; here, by contrast, we choose the x-axis of our 2D anisotropic Ising model -the one with the stronger coupling (J x > J y ) -to correspond to the time direction of our quantum circuit.

E. The M -qubit transfer matrix
For complex {J, γ, A}, the transfer matrix T = VW, which is of dimension 2 M × 2 M , is in general not normal, i.e. it does not in general commute with its Hermitian conjugate. Nevertheless, any non-normal complex matrix V can be brought to Jordan canonical form by a similarity transformation T = R −1 T R, where R is invertible. If we assume there are no Jordan blocks, then T may be decomposed in terms of its eigenvalues and its left and right eigenvectors, viz.
where L a || R b = δ ab , and where there is no complex conjugation implied in the bra vector L a || with a doubled bracket. The eigenvalues {λ a } are in general complex. If we order the eigenvalues such that |λ a | > |λ a+1 | for all a, then assuming the largest eigenvalue λ 0 is nondegenerate, after a sufficiently large number of iterations s we have It is convenient to here and henceforth implement a similarity transformation and redefine T ≡ V 1/2 W V 1/2 , which is manifestly symmetric: T = T . The corresponding right and left eigenvectors of T are then mutual transposes, with no complex conjugation, which we write as || Ψ a and Ψ a || , respectively. We consider two natural correlation functions which may be used to characterize the properties of the circuit. The first is the quantum two time correlator, With ρ 0 = 1, we have The second is the statistical correlator, and thus the s-independent terms in the above two correlators both vanish. We then have that both C(s; i, j) and C L (s; i, j) decay exponentially in the time direction with a correlation time τ = 1/ ln |λ 0 /λ 1 | and a frequency ω = arg(λ 1 /λ 0 ) which is generally incommensurate (i.e. irrational). When the spectral gap collapses, both correlation functions become long-ranged.
At short times we do not expect them to agree, e.g., for unitary circuits, quantum correlators obey rigid Lieb-Robinson bounds with strictly vanishing correlators outside the light cone, while statistical correlators are small but finite for spacelike separations at short times.
It will also be useful to extend some of our expressions from real-temperature thermodynamics to complex inverse temperature, β. We examine the modulus the partition function and define the free energy density accordingly f ≡ log |Z|/M Lβ, where M × L is the total number of spins, followed by the internal energy density and the specific heat capacity We typically plot our results not in the complex β-plane, but rather in the complex tanh β plane, which we shall refer to simply as the 'complex temperature plane'.

III. LARGE-M LIMIT: THE ANISOTROPIC 2D ISING MODEL
In this Section we consider our quantum circuit in the limit of a large number of qubits, M 1. In the M → ∞ limit, the dynamical correlation functions of the circuit may be written in terms of the statistical correlations of an anisotropic 2D Ising model at complex temperature. We analyze this model using three complementary methods: analytic continuation of the Onsager solution; numerical evaluation of a tensor-network representation of the partition and correlation functions; and exact fermionization of the zero-field problem using the Jordan-Wigner transformation.

A. Thermodynamics from the Onsager solution
In the isotropic Ising model, the zeros of the partition function lie on linear contours in the complex-temperature plane. In such a case, provided that the linear density of zeros reaches a finite value in the thermodynamic limit, we expect a simple slope discontinuity in the free energy, as already noted by Fisher 50 . The anisotropic Ising model was examined similarly; however, in that case we observe a far more complicated situation with patterns of zeros that appear to occupy extended regions in the complex-temperature plane. This makes the expected behavior of the free energy less clear.
We note at this point that Fisher's observation and the majority of others that have followed it are in fact based on a portion of Onsager's result; as demonstrated in appendix A, it is manifestly incorrect for the case of finite-width Ising ladders (corresponding to circuits with a finite number of qubits). Nevertheless, one might anticipate that the approximation remains asymptotically exact for 2D bulk (intensive) quantities. The Onsager expression for the real part of the (dimensionless) free energy per spin of the anisotropic model is where j x,y ≡ 2βJ x,y , where J x ≡ 1 and J y ≡ J. For the case J y = 0.1 and complex β we have evaluated this for the infinite system numerically. Fig. 1 was obtained by taking numerous cuts through the complex temperature plane. One particular cut that is especially revealing is a radial cut away from the real-temperature axis (Fig. 3) that clearly displays the continuous nature of the PM-NFM1 transition.

B. Correlations from tensor-network renormalization
We would like to characterize the different phases that appear in Fig. 1, especially the NFM1 and NFM2 phases that are not simple continuations of real-temperature phases. For this we need to know the spin-spin correlation functions in both the x (time) and y (qubit array) directions. Our most general method for determining these, which has the additional advantage of allowing the inclusion of a longitudinal magnetic field, is via a renormalization group algorithm based on tensor networks -in particular, the Tensor Renormalization Group (TRG) 48 , a method that involves representing classical partition functions as tensor networks and coarse-graining these tensor networks numerically.
Most of our results are obtained with bond dimensions up to 50, and we use relatively modest convergence goals which we check throughout, e.g. that the free energy density is converged to ∼ 0.001. In the J y → 0 limit, our system is a set of uncoupled Ising chains; we know that, in this limit, the entire complex-temperature plane is paramagnetic with the exception of the unit circle | tanh β| = 1. We therefore expect that,  Fig. 1). Top right: The specific heat capacity along the same contour (see text for the precise definition of 'specific heat capacity' at complex temperature). Bottom panels: The same specific heat graph, but multiplied by the specified factor, demonstrating that the two singularities are one-sided square-root singularities.
for J y 1, correlated non-paramagnetic phases will be concentrated near the unit circle.
We have used finite-field TRG to establish the disappearance of uniform ferromagnetic order as we traverse the unit circle | tanh β| = 1 counterclockwise from the realtemperature line, for several values of the coupling J y -see Fig. 4. It is clear that the FM phase is progressively suppressed as the interchain coupling is reduced. We interpret this as the gradual reversion to the incommensurately modulated order seen on the unit circle in the case of decoupled chains. This perspective already suggests that the NFM1 phase exhibits some form of long-range incommensurate order; we show below that that is essentially true.
Additional work is required to use TRG to compute correlation functions. As the method is multiplicative, i.e. distances are reduced by a factor of 2 per iteration, evaluating the correlation function at separations 2 n is relatively easy. These are useful in cases for which we expect simple power-law decays. Here, however, we are interested in modulated correlators, which requires careful renormalization at short distances.
Our TRG-computed correlation functions in the NFM1 phase are shown in the left-hand panels of Fig. 5. We observe modulated correlations along the direction of strong coupling (x in the statistical mechanics setting; time in the quantum circuit picture). The correlations along the weak direction (y in the statmech picture; inter-qubit in the quantum circuit picture) are non-oscillatory and apparently power-law decaying. In the next section we shall provide an interpretation of these results in terms of a fermionized version of the model.
Our TRG results for the correlation functions in the NFM2 phase are shown in the middle panels of Fig. 5. Surprisingly, we discover that the directions of the modulated and simply-decaying correlations are swapped in NFM2 relative to NFM1. It is worth noting that for the NFM1 phase the J y → 0 limit is solvable and contains modulated correlations already, whereas the NFM2 phase does not exist in the decoupled limit, owing its existence to interchain interactions.
The upper right-hand panel of Fig. 5 shows the correlation function in the x (time) direction at a different point in the NFM1 phase, where the period of the temporal oscillations is shorter. It also shows the effect on these correlations of the application of a longitudinal magnetic field. We see that they survive almost entirely unaltered, i.e. that this phenomenon is insensitive to the breaking of integrability.
The lower right-hand panel of Fig. 5 shows the magnetization as a function of longitudinal field in the modulated phases. Up to some critical field (which depends on the extent of anisotropy), these modulated phases are stable. At this critical field a metamagnetic transition takes place and the system exits the modulated phase.
The TRG method makes it straightforward to compute correlations even for large systems, for particular exponentially spaced sets of distances. Evaluating the spatial decay of the correlator in the NFM1 phase, we see clear signs of algebraic decay (Fig. 6), exactly as one would expect in a onedimensional critical phase.

C. Fermionization and the origin of oscillations
Finally, we present an approach to the zero-field problem that uses the Jordan-Wigner transformation to represent the spins/qubits in terms of fermionic degrees of freedom. We shall show that the occurrence of correlations is due to a certain type of resonance between two eigenvalues of these fermionic operators. We may use this picture to predict both the temporal period of the oscillations in the NFM1 phase and the spatial period of the oscillations in the NFM2 phase; in both cases, we find good agreement with our TRG results pre-sented above.
In the Jordan-Wigner representation, we write the Pauli matrices on site j as follows: where the operator c j annihilates a (spinless) fermion on site j of the qubit register. It follows that We now Fourier transform the fermion operators with respect to the index j, i.e. we move to a plane-wave basis in our qubit register: where the prime on the sum restricts k to the first Brillouin zone, i.e. k ∈ [−π, π). In terms of these plane-wave operators, the components of our single-step time-evolution operator become with T = V 1/2 W V 1/2 . We may streamline our notation using Anderson pseudospin operators τ α k , defined as follows: where α ∈ {0, 1, 2, 3}. In terms of these operators, the gates can now be re-written as follows: The partition function Z = Tr T L may be expressed as the product For each wavevector k > 0, the operator Θ k has two eigenvalues, λ k,± (see appendix C), where h k and δ k are given by h k = cosh(2β) cosh(2J) + sinh(2J) cos(k) (28) and δ k = sinh 2 (2β) sinh 2 (2J) sin 2 k (29) where we have used a mixed notation, trading coupling constants γ and A for β -recall the relations tanh γ = e −2β and A = √ coth γ − tanh γ. To analyze the late-time properties of the evolution, we find the largest-amplitude eigenvalue of Θ k (eq.26) for each wavenumber k. We denote this as λ 0 (k) and the corresponding right eigenvector as || ψ 0 (k) . These eigenvalues determine the decay rate and precession of a typical initial condition at late times (t → ∞): We now discuss the behavior of this late-time state in terms of properties of the fermionic spectrum. Outside the NFM1 phase, either the + or − branch is consistently largeramplitude throughout the Brillouin zone, and there is generically a unique steady state || 0 . In the NFM1 phase, however, the branches "invert" as a function of k, which is to say |λ + (k)| > |λ − (k)| for k ≈ 0 but the opposite inequality holds for k ≈ ±π. At special momenta ±k * , the two eigenvalues are degenerate. Therefore, in the subspace of ±k * , the system never reaches a unique steady state, and instead one has  I. Comparison of numerically estimated periods of order parameter oscillations (TTRG) against exact Jordan-Wigner fermion periods, TJW ≡ 4π/ω * . Note the addtional factor of 2 due to a twosite unit-cell in the time-direction implicit in the definition of the one step evolution operator. First three rows corresponds to points inside NFM1, where modulations are along the x-axis, while the last row is in NFM2. persistent oscillations at the frequency The band inversion point k * sweeps across the Brillouin zone as one progresses through the NFM1 phase, leading to incommensurate temporal modulations of varying frequency. By comparing with the numerical solutions we see that this phase exhibits temporal oscillations but no apparent spatial oscillations. While the momentum 2k * appears to be special in some sense, from the above argument, there is no simple relation between spectral degeneracies of the sort described above and spatial oscillations. To capture modulated correlations in the NFM2 phase, it is convenient instead to fermionize the model sideways, by performing the Jordan-Wigner transformation along the x axis (which hitherto we took to be the temporal direction). The "band inversion" described above now happens in the NFM2 phase, leading to oscillations in the spatial direction (i.e., along y).

IV. DISCUSSION
In the present work we have exploited the correspondence between non-unitary quantum circuits and complex temperature statistical mechanics to construct a simple quantum circuit that has a surprisingly rich phase diagram, including a phase with incommensurate temporal order. Such incommensurate time crystals do not seem to occur in closed systems; nor do they occur in one-dimensional open quantum systems with short-range interactions, for entropic reasons. Our results suggest that an important class of quantum circuits that exhibit incommensurate time-crystalline order are spacetime duals of circuits that realize incommensurate density-wave phases. In this simple one-dimensional case, such phases occur only for complex couplings, but in more general settings it might be possible to write down quantum circuits that cool the system into a ground state with density-wave order 51 . These would also be spacetime dual to temporally modulated phases.
In practice, post-selection is an expensive operation requiring effort that scales exponentially in the area of the quantum circuit. Thus, practical realizations of the physics discussed here will be restricted to circuits that are either very shallow or involve only a small number of qubits evolved for a long time. These map onto Ising ladders at complex temperature, which can be solved using the methods discussed above (App. A). We find that signatures of the modulated phases are present even for systems with modest numbers of spins (M = 5), which should be realistic to explore in a variety of presentday experiments. to be the only type of correlated 'phase' that occurs in the isotropic case. Anisotropic lattices, however, appear to support another type of 'gapless' correlated phase, which exhibits multiple long length-scales, and which shows precursor signatures in the behavior of the correlation lengths in the finite-M case.
We begin by reviewing our prior results 30 , in which the limiting behavior of Fisher zeros in ladders was computed -see Fig. 7. It can be shown, via a low-temperature expansion, that the number of contours (branch cuts) emerging from the two zero-temperature points tanh β = ±1 is equal to M . It is less clear how to compute the total number of contours, although contours that do not emerge from tanh β = ±1 do not appear to proliferate, and may be strongly dependent on boundary conditions.
Fisher's original proposal overlooked this behavior entirely. The exact solutions of the Ising model by Onsager and several others in the years that followed 52 usually consist of several contributions, only one of which dominates in the thermodynamic limit. Fisher's original argument for generalizing Yang-Lee results was based on a seemingly incorrect procedure whereby he analytically continued only the portion of the result that was important at real temperature. As explicitly demonstrated by Beichert et al., this produces entirely wrong patterns of zeros in ladders. Remarkably, however, Fisher's approximate solution is accurately reproduced by the unbiased TRG computational scheme applied to the 2D lattice.
Next we examine the growth of correlation lengths in the x (time) direction as we increase M . Each correlation length is controlled by the ratio between one of the subdominant eigenvalues of the transfer matrix, λ j (j > 0), and the dominant one, λ 0 . In the left-hand panels of Fig. 8 areas mark the regions of the complex-temperature plane in which the longest correlation length, i.e. the one controlled by λ 1 /λ 0 , is greater than ten lattice spacings. This is shown for the two-qubit case (top left) and for the five-qubit case (bottom left). It is clear that our arbitrarily determined threshold of 10 lattice sites is already exceeded for M = 5 in the entire crescent region outlined by Fisher's original 2D proposal. The right-hand panels show the same information but for the second-longest correlation length, i.e. the one controlled by λ 2 /λ 0 . We note that there is no other long correlation length present; as we can see from these plots, the eigenvalue λ 2 does not approach λ 0 except near the unitary point tanh β = i, i.e. the ferromagnetic phase is 'gapped'. We next turn to the case of anisotropic couplings, J y < 1. It was observed by van Saarloos and Kurtze 53 that, in this case, Fisher's approximation to the partition function produces highly complex patterns of zeros. For simple integer fractions J y = 1/n it is possible to compute and plot the contours 54 and observe an erratic pattern that does not show simple convergence to the limit of decoupled qubits. Numerically exact computations for finite-M transfer matrices, however, produce a nicely regular and convergent pattern, from which Fig. 9 was obtained. The ordered phase is reduced in its extent as we have reduced J y tenfold. As with the isotropic model we observe a growth of the longest correlation length as we increase M . However, in this anisotropic case we also see additional regions away from the unitary point in which the second-longest correlation length also becomes long: one tracking the unit circle, and another small patch on the real tanh β axes at values exceeding 1. These locations are suggestively similar in shape and location to the NFM1 and NFM2 phases in Fig. 1. Appendix B: TRG results for the magnetization : isotropic case In this appendix we present our TRG results for the behavior of the isotropic Ising model in a longitudinal field.
Based on the patterns of Fisher zeros, it seems that the isotropic Ising model in zero applied magnetic field exhibits, for general complex temperature, a first-order transition between a paramagnetic (PM) and a ferromagnetic (FM) phase. We explicitly verify this through a direct computation of the magnetization process in the vicinity of that phase boundary.
The phase boundary between the PM and FM phases forms the famous double crescent in the tanh β plane, but becomes a yet simpler unit circle in the complex sinh 2β plane: this appears to be a more natural variable for the isotropic case. A trajectory with a fixed modulus of sinh 2β slightly above/below 1 allows us to study the evolution of the character of the transition as we move away from the real-temperature critical point.
In Fig. 10, we show our results for the magnetization as a function of field for a variety of angles along two such circular contours in the complex sinh 2β plane. In the upper panel, we observe that the FM phase indeed exhibits spontaneous magnetization, the magnitude of which drops to zero as the point arg sinh 2β = π/2 is approached. The PM phase shows initially linear response, followed by a metamagnetic jump, as expected near conventional first-order transitions 55 . Interestingly, there is a large domain of diamagnetic response, which would be a rather peculiar state of affairs for a more conventional statistical problem.

Appendix C: Fermionization
In this appendix, we provide the details of the various steps of our fermionization approach to the zero-field case.
In §III C it was shown that the partition function Z = Tr V 1/2 W V 1/2 L is given by Z = k>0 Tr Θ L k , where Θ k ≡ A 2 e γτ z k e 2J(τ z k cos k+τ y k sin k) e γτ z k .
Steady state -We now obtain an expression for the steady state, where each (k, −k) mode pair is in an eigenstate of Θ k .
In the t → ∞ limit, and at each wavevector k ∈ (0, π), one of these states is selected -the one corresponding to the greater value of λ ± (k) . The surviving state's wavefunction is given by the appropriately normalized right eigenvector, and the asymptotic state is thus of the BCS form, where C k = 1/ 2 |Ω k,+ | + |Ω k,− | , and where | 0 is the Fock space vacuum, equivalent to the state | ↓↓ · · · ↓ for all the k-space Anderson pseudospins.