Practical recipes for the model order reduction, dynamical simulation, and compressive sampling of large-scale open quantum systems

This article presents numerical recipes for simulating high-temperature and non-equilibrium quantum spin systems that are continuously measured and controlled. The notion of a spin system is broadly conceived, in order to encompass macroscopic test masses as the limiting case of large-j spins. The simulation technique has three stages: first the deliberate introduction of noise into the simulation, then the conversion of that noise into an equivalent continuous measurement and control process, and finally, projection of the trajectory onto a state-space manifold having reduced dimensionality and possessing a Kahler potential of multi-linear form. The resulting simulation formalism is used to construct a positive P-representation for the thermal density matrix. Single-spin detection by magnetic resonance force microscopy (MRFM) is simulated, and the data statistics are shown to be those of a random telegraph signal with additive white noise. Larger-scale spin-dust models are simulated, having no spatial symmetry and no spatial ordering; the high-fidelity projection of numerically computed quantum trajectories onto low-dimensionality Kahler state-space manifolds is demonstrated. The reconstruction of quantum trajectories from sparse random projections is demonstrated, the onset of Donoho-Stodden breakdown at the Candes-Tao sparsity limit is observed, a deterministic construction for sampling matrices is given, and methods for quantum state optimization by Dantzig selection are given.

This article presents practical numerical recipes for simulating high-temperature and nonequilibrium quantum spin systems that are continuously measured and controlled. The notion of a "spin system" is broadly conceived, in order to encompass macroscopic test masses as the limiting case of large-j spins. The simulation technique has three stages: first the deliberate introduction of noise into the simulation, then the conversion of that noise into an informatically equivalent continuous measurement and control process, and finally, projection of the trajectory onto a Kählerian state-space manifold having reduced dimensionality and possessing a Kähler potential of multilinear (i.e., product-sum) functional form. These state-spaces can be regarded as ruled algebraic varieties upon which a projective quantum model order reduction (QMOR) is performed. The Riemannian sectional curvature of ruled Kählerian varieties is analyzed, and proved to be non-positive upon all sections that contain a rule. It is further shown that the class of ruled Kählerian state-spaces includes the Slater determinant wave-functions of quantum chemistry as a special case, and that these Slater determinant manifolds have a Fubini-Study metric that is Kähler-Einstein; hence they are solitons under Ricci flow. It is suggested that these negative sectional curvature properties geometrically account for the fidelity, efficiency, and robustness of projective trajectory simulation on ruled Kählerian state-spaces. Some implications of trajectory compression for geometric quantum mechanics are discussed. The resulting simulation formalism is used to construct a positive P -representation for the thermal density matrix and to derive a quantum limit for force noise and measurement noise in monitoring both macroscopic and microscopic test-masses; this quantum noise limit is shown to be consistent with wellestablished quantum noise limits for linear amplifiers and for monitoring linear dynamical systems. Single-spin detection by magnetic resonance force microscopy (MRFM) is then simulated, and the data statistics are shown to be those of a random telegraph signal with additive white noise, to all orders, in excellent agreement with experimental results. Then a larger-scale spin-dust model is simulated, having no spatial symmetry and no spatial ordering; the high-fidelity projection of numerically computed quantum trajectories onto low-dimensionality Kähler state-space manifolds is demonstrated. Finally, the high-fidelity reconstruction of quantum trajectories from sparse random projections is demonstrated, the onset of Donoho-Stodden breakdown at the Candès-Tao sparsity limit is observed, and methods for quantum state optimization by Dantzig selection are given.

Introduction
This article describes practical recipes for the simulation of large-scale open quantum spin systems. Our overall objective is to enable the reader to design and implement practical quantum simulations, guided by an appreciation of the geometric, informatic, and algebraic principles that govern simulation accuracy, robustness, and efficiency.

How does the Stern-Gerlach effect really work?
This article had its origin in a question that Dan Rugar asked of us about five years ago: "How does the Stern-Gerlach Effect really work?" The word "really" is noteworthy because hundreds of articles and books on the Stern-Gerlach effect have been written since the original experiments in 1921 [88,89,90] . . . including articles by the authors [178,182] and by Dan Rugar [179] himself. Yet we were unable to find, within this large literature, an answer that was satisfactory in the context in which the question was asked, that circumstance being the (ultimately successful) endeavor by Rugar's IBM research group to detect the magnetic moment of a single electron spin by magnetic resonance force microscopy (MRFM) [170].
1.1.1. Constraints upon the analysis Quantum theory has a reputation for mystery. But as Peter Shor has remarked, "Interpretations of quantum mechanics, unlike Gods, are not jealous, and thus it is safe to believe in more than one at the same time." In particular, it is well known-and we will review the literature in this article-that advances in quantum information theory have provided Shor's principle with rigorous foundations. We will build upon these informatic foundations in answering Dan Rugar's question in accord with the following constraints: our analysis will be orthodox in its respect for established principles of quantum physics. It will be operational in the sense that all its predictions are traceable to explicitly hardware-based measurement processes. The analysis will be scalable to accommodate large-dimension quantum systems (such as the spins in protein molecules that are the ultimate targets of MRFM microscopy). The analysis will be reductive in the sense that the analysis will yield simple design rules that are in reasonable quantitative accord with the predictions of more accurate-but more complicated-largescale numerical simulations. The analysis will be synoptic in the sense that when we are required to choose between equivalent analysis formalisms, a rationale for these choices will be provided, and the consequences of alternative choices noted. And finally, the analysis will be extensible-at least in principle-to the analysis and simulation of general quantum systems (such as spintronic devices, nanomechanical devices, and biomolecules).
There are of course strong practical motivations for seeking to analyze quantum systems by methods that are orthodox, operational, scalable, reductive, synoptic, and extensible: these same attributes are essential to practical methods for analyzing large-scale classical systems [20,113].

The feasibility of generic large-scale quantum simulation
We did not begin our investigations with the idea that the numerical simulation of large-scale quantum spin systems was feasible. Indeed, we were under the opposite impression, based upon the no-simulation arguments of Feynman [67] in the early 1980s. These arguments have been widely-and usually uncritically-repeated in textbooks [146, sec. 4.7]. But Feynman's arguments do not formally apply to noisy systems, and in the course of our analysis, it became apparent that this provides a loophole for developing efficient simulation algorithms. Furthermore, it became apparent that the class of noisy systems encompasses as a special case the low-temperature and strongly correlated systems that are studied in quantum chemistry and condensed matter physics. In the concluding section of this article (Section 4.6), we will develop the point of view that any quantum state that has been in contact with a thermal reservoir is an algorithmically compressible object.
This loophole helped us understand why-from an empirical point of view-simulation capabilities in quantum chemistry and condensed matter physics have been improving exponentially in recent decades [80,149,151]. The analysis and simulation methods that we will present in this article broadly define a geometric and quantum informatic program for sustaining this progress.
1.2.1. The geometry of reduced-order state-spaces This article's mathematical methods are novel mainly in their focus upon the geometry of reduced-order quantum state-spaces. We will show that the quantum state-spaces that are most useful for large-scale simulation purposes generically have an algebraic structure that can be geometrically interpreted as a network of geodesic curves (rules) having nonpositive Gaussian curvature for all sections that contain a rule (Section 2). We will see that these curvature properties are essential to the efficiency and robustness of model order reduction.

The central role of covert measurements
A technique that is central to our simulation recipes is to simulate all noise processes (including thermal baths) as equivalent covert measurement and control processes (Section 3). From a quantum informatic point of view, covert quantum measurement processes act to quench high-order quantum correlations that otherwise would be infeasibly costly to compute and store (Section 4). Thus the presence of noise can allow quantum simulations to evade the no-simulation arguments of Feynman [67].

1.2.3.
Background assumed by the presentation No reader will be expert in all of the disciplines that our analysis and simulation recipes embody, which are (chiefly) quantum mechanics in both its physical and informatic aspects, the engineering theory of model order reduction (MOR) and dynamical control, and the mathematical tools and theorems of algebraic and differential geometry. Indeed, the writing this article has made the authors acutely aware of their own considerable deficiencies in all of these areas. Recognizing this, we will describe all aspects of our recipes at a level that is intended to be broadly comprehensible to nonspecialists.

Overview of the analysis and simulation recipes
We begin by surveying our chain of reasoning in its entirety. Figures 1-4 concisely summarize the simulation recipes and their geometric basis. In a nutshell, the recipes embody the orthodox quantum formalism of Fig. 1, as translated into the practical numerical algorithm of Fig. 2, which is based upon the algebraic structures of Fig. 3, whose functionality depends upon the fundamental geometric concepts of Fig. 4.
The following overview summarizes those aspects of the simulation recipes that are fundamental, multidisciplinary, or novel, and it also seeks to describe the embedding of these recipes within the larger literature.

Overview of the formal simulation algorithm
Formally, our simulation algorithms will be of the general quantum information-theoretic form that is summarized in Fig. 1. Steps A.1-2 of the algorithm are adopted, without essential change, from the axioms of Nielsen and Chaung [146], and our discussion will assume a background knowledge of quantum information theory at the levels of Chapters 2 and 8 of their text.
Step A.3 of the simulation algorithm-projection of the quantum trajectory onto a statespace of reduced dimensionality-will be familiar to system engineers as projective model order reduction (as we will review in Sec. 1.4). We will also establish that projective MOR Steps A.1-2 summarize the formal theory of the simulation of quantum systems (see, e.g., Nielsen and Chuang [146, chs. 2 and 8]). Step A.3 is a model order reduction of the Hilbert states |ψ n by projection onto a reduced-dimension Kähler manifold K (see e.g. Rewieński [164]). Equivalently, Step A.3 may be viewed as a variety of Dirac-Frenkel variational projection (see, e.g., [133,161]).
is formally identical to a method that is familiar to physicists and chemists as a variational order reduction of Dirac-Frenkel-McLachlan type [59,76,136] (see also the recent references [98,133,161]).
For purposes of exposition, we define quantum model order reduction (QMOR) to be simply classical MOR extended to the complex state-space of quantum simulations.
1.3.1. The operational approach to quantum simulation Our simulation recipes will adopt a strictly operational approach to measurement and control, in the sense that we will require that the only information stream used for purposes of communication and control is the stream of binary stochastic outcomes of the measurement operations of Step A.2 of Fig. 1. Although it is not mathematically necessary, we will associate these binary outcomes with the classical "clicks" of physical measurement apparatuses, and we will develop a calibrated physical model of these clicks that will guide both our physical intuition and our simulation design.
1.3.2. The embrace of quantum orthodoxy Because the binary "clicks" of measurement outcomes are all that we seek to simulate, our analysis will regard the state trajectories {ψ 1 , ψ 2 , . . . , ψ n , . . . } as wholly inaccessible for all purposes associated with measurement and control, which is to say, as inaccessible for engineering purposes. We will analyze quantum state trajectories only with the goal of tuning the simulation algorithms to compress the trajectories onto low-dimension manifolds. In practice, this will mean that we mainly care about the geometric properties of quantum trajectories; this will be the organizing theme of our analysis.
In the course of our analysis we will confirm-mainly to check our algebraic manipulations-that several of the traditional quantum measurement short-cuts that deal directly with wave-functions (e.g., uncertainty principles, wave function collapse, quantum Zeno effects) yield the same results as our "clicks-only" reductive formalism. But our simulations will not use these short-cuts, and in particular, we will never simulate quantum measurement processes in terms of von Neumann-style projection operators.
The resulting simulation formalism is wholly operational, and can be informally described as "ultra-orthodox." The operational approach will require some extra mathematical work-mainly in the area of stochastic analysis-but it will also yield some novel mathematical results, including a closed-form positive P -representation [153] of the thermal density matrix. We will derive this P -representation by methods that provably simulate finite-temperature baths. Thus the gain in practical simulation power will be worth the effort of the extra mathematical analysis.

The unitary invariance of quantum operations
Our analysis will focus considerable attention upon the sole mathematical invariance of the simulation algorithm of Fig. 1, which is a unitary invariance associated in the choice of the quantum operations M in Step A.2. Our main mathematical discussion of this invariance will be in Section 3.2.1, our main discussion of its causal aspects will be in Section 3.3.6, our main review of the literature will be in Section 3.3.7, and it will be central to the discussion of all the simulations that we present in Section 4.
We will see that the short answer to the question "What is this unitary invariance all about?" is that (1) it ensures that measured quantities respect physical causality, and (2) it allows quantum simulations to be tuned for improved efficiency and fidelity.
It is prudent for students to browse among these works to find congenial points of view, because two of the above references are alike in the significance that they ascribe to the unitary invariance of quantum operations. This diversity arises because the invariance can be understood in multiple ways, including physically, algebraically, informatically, and geometrically. Our analysis will touch upon all these aspects, but much more than any of the above references, our approach will be geometric.
1.3.4. Naming and applying the Theorema Dilectum It is vexing that no short name for the unitary invariance associated with quantum operations has been generally adopted. For example, this theorem is indexed by Nielsen and Chuang under the unwieldy phrase "theorem: unitary freedom in the operator-sum representation" [146, thm. 8.2, sec. 8.2].
Because we require a short descriptive name, we will call this invariance the Theorema Dilectum, which means "the theorem of choosing, picking out, or selecting" (from the Latin deligo). As our discussions will demonstrate, this name is appropriate in both its literal sense and in its evocation of Gauss' Theorema Egregium.
In this article we will develop a geometric point of view in which the Theorema Dilectum is mainly a theorem about trajectories in state-space, and that the central practical role of the theorem in quantum simulations is to enable noisy quantum trajectories to be algorithmically compressed, such that efficient large-scale quantum simulation is feasible. To the best of our knowledge, no existing articles or textbooks have assigned to the Theorema Dilectum the central geometric role that this article focusses upon.
The Theorema Dilectum is the first of two main technical terms that we will introduce in this review. To anticipate, the other is gabion, which is the name that we will give to state-space manifolds that support a certain kind of affine algebraic structure (see Figs. 3-4 and Section 1.5). When gabion manifolds are endowed with a Kähler metric, we will call the result a gabion-Kähler manifold (GK manifold 1 ).
GK manifolds are the state-spaces onto which we will projectively compress the quantum trajectories of our simulations by exploiting the Theorema Dilectum. When we further impose an antisymmetry condition upon the state-space the result is a Grassmannian gabion-Kähler manifold (GGK manifold), and we will identify these manifolds as being both the well-known Slater determinants of quantum chemistry, and the equally well-known Grassmannian varieties of algebraic geometry.
1.3.5. Relation to geometric quantum mechanics Our recipes will embrace the strictly orthodox point of view that linear quantum mechanics is "the truth" to which our reducedorder Kählerian state-spaces are merely a useful low-order approximation. However, at several points our results will be relevant to a logically conjugate point of view, known as geometric quantum mechanics, which is described by Ashtekar and Schilling as follows [9] (see also [13,172]): [In geometric quantum mechanics] the linear structure which is at the forefront in text-book treatments of quantum mechanics is, primarily, only a technical convenience and the essential ingredients-the manifold of states, the symplectic structure and the Riemannian metric-do not share this linearity.
Thus in geometric quantum mechanics, Kählerian geometry is regarded as a fundamental aspect of nature, while in our quantum MOR discussion, this same geometry is a matter of deliberate design, whose objective is optimizing simulation capability. Because our main focus is upon quantum MOR, we will comment only in passing upon those results that are relevant to geometric quantum mechanics (e.g., see the discussion in Section 2.12).

Overview of the numerical simulation algorithm
The numerical simulation algorithm of Fig. 2 is simply the formal algorithm of Fig. 1 expressed in a form suitable for efficient computation. Note that Fig. 2 adopts the MATLABstyle engineering nomenclature of model order reduction, as contrasted with the physicsstyle bra-ket notation of Fig. 1.
The algorithm of Fig. 2 is a fairly typical example of what engineers call model order reduction (MOR) [7,81,143,147,148]. Rewieński's thesis is particularly recommended as a review of modern nonlinear MOR ( [164], see also [165]).

The main ideas of projective model order reduction
We will now briefly summarize the main ideas of projective MOR in a form that well-adapted to quantum simulation purposes. We consider a generic MOR problem defined by the linear equation δψ = Gψ.
Here ψ is a state vector, δψ is a state vector increment, and G is a (square) matrix. For the present it is not relevant whether ψ is real or complex. It commonly happens that ψ includes many degrees of freedom that are irrelevant to the practical interests that motivate a simulation. The central physical idea of MOR is to adopt a reduced order representation ψ(c), where c is a vector of model coordinates, having dim c dim ψ. The central mathematical problem of MOR is to describe the large-dimension increment δψ by a reduced-order increment δc. It is convenient to organize the partial derivatives of ψ(c) as a non-square matrix A(c) whose elements are [A(c)] ij ≡ ∂ψ i /∂c j . The reduced-order increment having least mean-square error is obtained by the following sequence of matrix manipulations: Here " P " is the Moore-Penrose pseudo-inverse that is ubiquitous in data-fitting and model order reduction problems [160], " † " is Hermitian conjugation, and the final step relies upon the pseudo-inverse identity X P = (X † X) P X † , which is exact for any matrix X [132]. This is the key step by which the master simulation equation is obtained that appears as Step B.3 at the bottom of Fig. 2. The great virtue of (1) for purposes of large-scale simulation is that (A † A) is a lowdimension matrix and (A † G ψ) is a low-dimension vector. Provided that both (A † A) and (A † G ψ) can be evaluated efficiently, and provided also that ψ(c) represents the "true" ψ with acceptable fidelity, substantial economies in simulation resources can be achieved. We will see that the required objectives of efficiency and fidelity both can be attained. For example, the low-dimension matrix∂⊗∂κ ≡ 1 2 A † A is obviously Hermitian (whether ψ is real or complex). Of what manifold is it the Hermitian metric tensor? How does this manifold's geometry influence the simulation's efficiency, fidelity, and robustness? Steps B.1-3 are a numerical recipe that implements the simulation algorithm of Fig. 1. The expressions (∂⊗∂κ) and (∂φ) that are introduced in Step B.3 serve solely as variable names for the stored partial derivatives of the Kähler potential κ(c, c) ≡ 1 2 ψ(c)|ψ(c) and the dynamic potential φ(c, c) ≡ 1 2 ψ(c)|δG|ψ(c) ; it is evident that these partial derivatives wholly determine the simulation's geometry and dynamics.
To answer these questions, we will show that κ is the Kähler potential of differential geometry, that the metric tensor∂⊗∂κ determines the Riemannian curvature of our reduced order state-space, and that the choice of an appropriate curvature for this state-space is vital to the simulation's efficiency, fidelity, and robustness.
In preparing this article, our search of the literature did not find a previous analysis of MOR state-space geometry from this Riemannian/Kählerian point of view. We did, however, find recent work in communication theory by Cavalcante [42] and coworkers [43,44] that adopts a similarly geometric point of view in the design of digital signal codes. Like us, these authors are unaware of previous similarly geometric work [43] "To the best of our knowledge this [geometric] approach was not considered previously in the context of designing signal sets for digital communication systems." Like us, they recognize that "[These state-spaces] have rich algebraic structures and geometric properties so far not fully explored." Also similar to us, they find [44] "The performance of a digital communication system depends on the sectional curvature of the manifold . . . the best performance is achieved when the sectional curvature is constant and negative." Our analysis will reach similar conclusions regarding the desirable properties of nonpositive sectional curvature in the context of quantum MOR. Because the mathematical basis of this apparent convergence of geometric ideas between MOR theory and coding theory is not presently understood (by us at least), we will not comment further upon it.
It is entirely possible that related work exists of which we are not aware. Model order reduction is ubiquitously practiced by essentially every discipline of mathematics, science, engineering, and business: the resulting literature is so vast, and the nomenclature so varied, that a comprehensive review is infeasible.
It is fair to say, however, that the central role of Riemannian and Kählerian geometry in model order reduction is not widely appreciated. A major goal of our article, therefore, is to analyze the Riemannian/Kählerian aspects of MOR, and especially, to link the Kählerian geometry of quantum MOR to the fundamental quantum informatic invariance of the Theorema Dilectum.

Preparing for a Kählerian geometric analysis
To prepare the way for our geometric analysis, at the bottom of Fig. 2 the pseudo-code defines storage variables named "(∂⊗∂κ)" and "(∂φ)." For coding purposes these names are of course purely conventional (an arbitrary string of characters would suffice), but these particular names are deliberately suggestive of partial derivatives of two scalar functions: κ and φ.
To anticipate, κ will turn out to be the Kähler potential of complex differential geometry, which determines the differential geometry of the complex state-space, and φ will turn out to be a stochastic dynamical potential, which determines the drift and diffusion of quantum trajectories on the Kählerian state-space.
The link between geometry and simulation efficiency thus arises naturally because both the geometry and the physics of our quantum trajectory simulations are determined by the same two scalar functions.

Overview of the unifying geometric ideas
The main algebraic and geometric features of our state-space are summarized in Figs. 3-4.
1.5.1. The algebraic structure of the reduced-order state space The state-space of all our simulations will have the algebraic structure shown in Fig. 3. We will regard this algebraic structure as a geometric object that is embedded in a larger Hilbert space, and we will seek to understand its geometric properties, including especially its curvature, in relation to our central topic of quantum model order reduction. The algebraic definition of a gabion-Kähler (GK) state-space (top) expressed equivalently as a matrix product state (MPS, bottom). By definition, the order of |ψ is the number of elements (spins) in each row's outer product, the rank of |ψ is the number of rows. The matrices A [l]m that appear at bottom are, by definition, r × r matrices-hence rank r-having diagonal elements (A [l]m ) kk ≡ l c m k and vanishing off-diagonal elements. Note that the matrix products are Abelian, such that the geometric properties of the state-space are invariant under permutation of the spins. Note also that when the above algebraic structure is antisymmetrized with respect to interchange of spins (equivalent to interchange of columns), the state becomes a sum of Slater determinants, or equivalently a join of Grasssmanian manifolds (a GGK manifold).
In the language of algebraic geometry [53,100], the geometric objects we will study are the algebraic manifolds that are associated with the projective algebraic varieties defined by the product-sums of Fig. 3. Although the literature on algebraic varieties is vast (and it includes many engineering applications [53]) and the literature on Riemannian sectional curvature is similarly vast, the intersection of these two subjects has apparently been little studied from an engineering point of view. This intersection, and especially its practical implications for quantum model order reduction, will be the main focus of our geometric investigations.
The general algebraic structure of Fig. 3 is known by various names in various disciplines. As noted in the caption to Fig. 3, these structures are known to physicists as a matrix product states (often abbreviated MPS) which are widely used in condensed matter physics and ab initio quantum chemistry [55,119,156,157,173,174,191]; these references provide entry to a rapidly growing body of MPS-related literature.
Quantum chemists have known the algebraic structures of Fig. 3 as Hartree product states [102] since 1928. Upon antisymmetrizing the outer products, we obtain the Slater determinants [184] that are the fundamental building-blocks of modern quantum chemistry; upon summing Slater determinants, and (optionally) imposing linear constraints upon these sums, we obtain post-Hartree-Fock quantum states [54]. All of the theorems we derive will apply to Slater determinants and post-Hartree-Fock states as special cases (see Section 2.9). We will comment later in this section, too, upon the intimate relation of these ideas to density functional theory (DFT). Nuclear physicists embrace these same ideas under the name of wave function factorization [150]. Beylkin and Mohlenkamp [16] note that statisticians call essentially the same mathematical objects canonical decompositions and also parallel factors.
As Leggett and co-authors have remarked with regard to the similarly immense literature on two-state quantum systems: "The topic of [this] paper is of course formally a problem in applied mathematics. . . . Ideas well known in one context have been discovered afresh in another, often in a language sufficiently different that it is not altogether trivial to make the connection. . . . [In such circumstances] the primary purposes of citations are to help the reader understand the paper, and the references in the text are chosen with this in mind" [129]. These same considerations will guide our discussion.
The general utility of affine algebraic structures for modeling purposes first came to our attention in a highly readable Mathematical Intelligencer article by Mohlenkamp and Monzón [140]; two further articles by Beylkin and Mohlenkamp [15,16] are particularly recommended also. Beylkin and Mohlenkamp call the algebraic structure of Fig. 3 a separated representation, and they have this to say about it [15]: When an algorithm in dimension one is extended to dimension d, in nearly every case its computational cost is taken to the power d. This fundamental difficulty is the single greatest impediment to solving many important problems and has been dubbed the curse of dimensionality. For numerical analysis in dimension d, we propose to use a representation for vectors and matrices that generalizes separation of variables while allowing controlled accuracy. . . . The contribution of this paper is twofold. First, we present a computational paradigm. With hindsight it is very natural, but this perspective was the most difficult part to achieve, and it has far-reaching consequences. Second, we start the development of a theory that demonstrates that separation ranks are low for many problems of interest.
In a subsequent article Beylkin and Mohlenkamp go on to say [16] "The representation seems rather simple and familiar, but it actually has a surprisingly rich structure and is not well understood." These remarks are remarkably similar in spirit to the coding theory observations of Cavalcante et al. that were reviewed in Section 1.4. For us, Kählerian algebraic geometry will provide a shared foundation for understanding the accelerating progress that all of the above large-scale computational disciplines have witnessed in recent decades.
1.5.2. The medieval idea of a gabion, and its mathematical parallels Deciding what to call the geometric state-space of quantum model order reduction is a vexing problem. We have seen that various plausible names include "Hartree products," "Slater determinants," "separated representations," "matrix product states," "wave function factorizations," "productsum states," "canonical decompositions," and "parallel factors." A shared disadvantage of the above names is that there is no precedent for associating them with the geometric properties that are the main focus of our investigations. We therefore seek an encompassing name for these state-spaces viewed as geometric entities.
Finding no precedent in the literature, and desiring a short name having a long-established etymology, we will call them by the medieval name of gabions [192], or more formally, gabion manifolds (this name arose spontaneously in the course of a seminar).
Most readers will have seen gabions numerous times, perhaps without recognizing that they have a well-established name. "Gabion" is the generic engineering name for a mesh basket that is filled with a weighty but irregularly-shaped material such as rocks or lumber, then stacked for purposes of reinforcement, erosion control, and fortification. In medieval times gabions were made of wicker or reed; Fig. 4(A) shows a typical medieval gabion. We will see that the defining geometry property of gabion manifolds is that they are possessed of a web of geodesic lines that constrain the curvature of the manifold, rather as the wicker reeds of a physical gabion constrain the curved rocks and boulders held inside. Like physical gabions, gabion manifolds come in a wide variety of sizes and shapes that are suitable for numerous practical purposes.
We will postpone giving a formal-and necessarily rather abstract-mathematical definition of a gabion until Section 2.5. For the present our main objective is to informally describe the geometric properties that will motivate this formal definition.

The geometric properties of gabion-Kähler (GK) manifolds
The main geometric properties of GK manifolds that are relevant to quantum simulation are depicted in Fig. 4(B-H). We will now survey these properties, and in doing so, we will introduce some of the nomenclature of Kählerian geometry.
We begin our geometric overview by remarking that even though Hilbert space is a complex state-space, a common viewpoint among mathematicians is that a complex manifold is a real manifold that is endowed with an extra symmetry, called its complex structure (see Sections 2.2 and 2.6 for details). For purposes of our geometric analysis, we will simply ignore this complex structure until we are ready to apply the quantum Theorema Dilectum (Section 2.6). Until then we will treat gabion manifolds as real manifolds.
In particular, the state-space of quantum mechanics has a natural real-valued measure of length. Specifically, along a time-dependent quantum trajectory |ψ(t) it is natural to define a real-valued velocity v(t) whose formal expression can be written equivalently in several notations: Hereψ(t) ≡ ∂ψ(t)/∂t, and we have used first the abstract notation of differential geometry (in which g(. . .) is a metric function), then the Dirac bra-ket notation of physics, then the matrix-vector notation of engineering and numerical computation, and finally the cumbersome but universal notation of components and sums over indices. We will assume an entry-level familiarity with all four notations, since this is a prerequisite for reading the literature.
As a token of considerations to come, the factor of dim H/2 in the index limit of (2) above arises because we will regard a complex manifold like C n as being a real manifold of dimension 2n. Thus we will regard the complex plane C as a two-dimensional (real) manifold, and an spin-1/2 quantum state as a point in a Hilbert space H having dim H = 4 (real) dimensions. This viewpoint leads to an ensemble of conventions that we will review in detail in Section 2.6. For now, we note that the arc length s along a trajectory is s = v(t) dt, so that geometric lengths in quantum state-spaces are dimensionless. An equivalent differential definition is to assign a length increment ds to a state increment |dψ via (ds) 2 = dψ|dψ . See Section 1.5 for a discussion of these principles.
Since we can now compute the real-valued length of an arbitrary curve on the gabion manifold, all of the usual techniques of differential geometry can be applied, without special regard for the fact that the state-space is complex.
1.5.4. GK manifolds are endowed with rule fields Beginning a pictoral summary of our geometric results, we first note that state-space gabions resemble physical gabions in that they are naturally endowed with a geometric mesh, which is comprised of a network of lines called rules, as depicted in Fig. 4(B). More formally, they are equipped with vector fields having certain mathematical properties (see (14) and Definition 2.1) such that the integral curves of the rule fields have the depicted properties.
Postponing a more rigorous and general definition of gabion manifolds until later (see Section 2.5), we can informally define a gabion rule to be the quantum trajectory associated with the variation of a single coordinate n c k r in the algebraic structure of Fig. 3, holding all the other coordinates fixed to some arbitrary set of initial values. We see that gabion rules are rays (straight lines) in the embedding Hilbert space, and hence, the gabion rules are geodesics (shortest paths) on the gabion manifold itself. As depicted in Fig. 4 (C), the set of all gabion points that belong to a rule is (trivially) the set of all gabion points itself, which is the defining characteristic of a gabion being ruled. Furthermore, we will show that at any given point, the vectors tangent to the rules that pass through that point are a basis set.
1.5.5. GK geometry has singularities Are gabion manifolds geometrically smooth, or do they have singularities? As depicted in Fig. 4(D), we will show that gabion manifolds have pinch-like geometric singularities. Algebraically these singularities appear whenever two or more rows of the product-sum in Fig. 3 are equal. Geometrically, we will show that the Riemann curvature diverges in the neighborhood of these singular points. However, it will turn out that the continuity of the geodesic rules is respected even at the singular points, so that gabion manifolds are geodesically complete. Pragmatically, this means that our numerical simulations will not become "stuck" at geometric singularities.
1.5.6. GK projection yields compressed representations As depicted in Fig. 4(E) model order reduction is achieved by the high-fidelity projection of an "exact" state |ψ in the large-dimension Hilbert space onto a nearby point |ψ K of the small-dimension gabion. It can be helpful to view this projection as a data-compression process. By analogy, the state |ψ is like an image in TIFF format; this format can store an arbitrary image with perfect fidelity, but consumes an inconveniently large amount of storage space. The projected state |ψ K on the gabion is like an image in JPEG format; lesser fidelity, but good enough for many practical purposes, and small enough for convenient storage and manipulation. We thus appreciate that data compression can be regarded as a kind of model order reduction; the two processes are fundamentally the same.
1.5.7. GK manifolds have negative sectional curvature Practical quantum simulations require that the computation of order-reducing projections be efficient and robust, just as we require image compression programs to be efficient and robust. As depicted in Fig. 4(F), order-reduction projection becomes ill-conditioned when the state-space manifold is "bumpy", in which case a numerical search for a high-fidelity projection can become stuck at local minima that yield poor fidelity. We will prove that the presence of a ruled net guarantees that gabion manifolds are always smooth rather than bumpy.
Resorting to slightly technical language to say exactly what we mean when we assert that gabion manifolds are not bumpy, in our Theorem 2.1 we will prove that a gabion has nonpositive sectional curvature for all sections on its geodesic net. This means that gabion manifolds can be envisioned as a net of surfaces that have the special property of being saddle-shaped everywhere (as contrasted with generic surfaces having dome-shaped "bumps"). As depicted in Fig. 4(G), the saddle-shaped curvature helps ensure that orderreducing projection onto gabion manifolds is a numerically well-conditioned operation.
1.5.8. GK manifolds have an efflorescing global geometry As depicted in Fig. 4(H), when the number of state-space dimensions becomes very large, it becomes helpful to envision nonpositively curved manifolds as flower-shaped objects composed of a large number of locally Euclidean "petals." This physical picture has been vividly conveyed by recent collaborative work between mathematicians and fabric artists [12]; the work of Taimina and Henderson on hyperbolic manifolds is particularly recommended [109].
When working in large-dimension spaces we will heed also Dantzig's remark that "one's intuition in higher dimensional space is not worth a damn!" [5]. For purposes of quantitative analysis we will rely upon Gauss' Theorema Egregium [86] to analyze the Riemannian and Kählerian geometric properties of gabion manifolds. We will prove that "number of petals" becomes exponentially large, relative to dim K, such that the petals loosely fill the embedding Hilbert space. In this respect our geometric analysis will parallel the informatic analysis of Nielsen and Chuang [146]; their Fig. 4.18 is broadly equivalent to our Fig. 4(H).
Our analysis will therefore establish two geometric properties of gabion manifolds: they are strongly curved, and they are richly endowed with straight-line rules. We will show that these gabion properties are essential to the efficiency, robustness, and fidelity of large-scale MOR. Later on, in Sections 4.6.4-4.6.6, we will establish a relation between these properties and compressive sampling (CS) theory.
1.5.9. GK basis vectors are over-complete We will nowhere assume that the basis vectors of the underlying algebraic structure of Fig. 3 are orthonormal; they might refer for example to the non-orthonormal gaussian basis states of quantum chemistry. A major geometric theme of our analysis, therefore, is that the negative sectional curvature of gabion manifolds helps generically account for the observed efficiency, fidelity, and robustness of gabion-based modeling techniques in many branches of science, engineering, and mathematics.
1.5.10. GK manifolds allow efficient algebraic computations Upon restricting our attention to the special case of gabion-Kähler (GK) manifolds, we will show that the existence of a ruled geodesic net allows the sectional curvature and the Riemann curvature tensors of GK manifolds to be calculated easily and efficiently. To anticipate, we will present data from Riemann curvature tensors having dimension up to 188, which we believe are the largestdimension curvature tensors yet numerically computed. We will see that it is Kraus' "long list of miracles" that makes large-scale numerical curvature computations feasible, and that these same miracles are equally essential to large-scale quantum dynamical calculations.
1.5.11. GK manifolds support the Theorema Dilectum One geometric idea remains that is key to our simulation recipes. For gabion manifolds to represent quantum trajectories with good fidelity, some physical mechanism must be invoked to compress quantum trajectories onto the petals of the gabion state-space. That key mechanism is, of course, the Theorema Dilectum that was mentioned in Section 1.3.4.
From a geometric perspective, the Theorema Dilectum guarantees that noise can always be modeled as a measurement process that acts to compress trajectories onto the GK petals. As depicted in Fig. 4(I), quantum simulation can be envisioned geometrically as a process in which compression toward the GK petals, induced by measurement processes, competes with expansion away from the petals, induced by quantum dynamical processes. The balance of these two competing mechanisms determines the MOR dimensionality that is required for good fidelity-the "petalthickness." Algebraically this petal thickness increases in proportion to the rank of the product-sum algebraic structure of Fig. 3.
Thus for us, trajectory compression is not a mathematical "trick," but rather is a reasonably well-understood and well-validated quantum physical mechanism, originating in the Theorema Dilectum, that compresses quantum trajectories to within an exponentially small fraction of the Hilbert phase space.
This noise-induced trajectory compression is the loophole by which QMOR simulations evade the no-simulation arguments of Feynman [67], as reviewed by Nielsen and Chuang [146, see their Section 4.7].
1.5.12. GK manifolds support thermal equilibria We will see that this covert-measurement approach encompasses numerical searches for ground states. Specifically, by explicit construction, we will show that contact with a zero-temperature thermal reservoir can be modeled as an equivalent process of covert measurement and control, in which the role of "temperature" is played by the control gain, such that zero temperature is associated with optimal control.
From this QMOR point of view, the calculation of a ground-state quantum wave function is a special kind of noisy quantum simulation, in which noise is present but masked by optimal control. This is how QMOR reconciles the strong arguments for the general infeasibility of ab initio condensed-matter calculations (as reviewed by, e.g., Kohn [121]) with the widespread experience that numerically computing the ground states of condensed matter systems is often, in practice, reasonably tractable [80].
1.5.13. GK manifolds support fermionic states Readers familiar with ab initio quantum chemistry, and in particular with density functional theory (DFT) [28,114,121, see Cappele [39] for an introduction] will by now recognize that QMOR and DFT are conceptually parallel in numerous fundamental respects: the central role of the low-dimension Kähler manifold of QMOR parallels the central role of the low-dimension density functional of DFT; the closed-loop measurement and control processes of QMOR parallel the iterative calculation of the DFT ground state; QMOR's fundamental limitation of being formally applicable only to noisy quantum systems parallels DFT's fundamental limitation of being formally applicable only to ground states; QMOR and DFT share a favorable computational scaling with system size.
Yet to the best of our knowledge-and surprisingly-the geometric techniques that that this article will deploy in service of QMOR have not yet been applied to DFT and related techniques of quantum chemistry and condensed matter physics [80]. A plausible starting point is to impose an antisymmetrizing Slater determinant-type structure upon the algebraic outer products of (3).
Some analytic results that we have obtained regarding the Kählerian geometry of Slater determinants are summarized in Section 2.9. With further work along these lines, we believe that there are reasonable prospects of establishing a geometric/informatic interpretation, via the Theorema Egregium and the Theorema Dilectum, of the celebrated Hohenberg-Kohn and Kohn-Sham Theorems of DFT [114] and their time-dependent generalization the Runge-Gross Theorem [28].
A physical motivation for this line of research is that the Theorema Dilectum of QMOR and the Hohenberg-Kohn Theorem of DFT embody essentially the same physical insight: the details of exponentially complicated details of quantum wave functions are only marginally relevant to the practical simulation of both noisy systems (QMOR) and systems near their ground-state (DFT).
At present, the two formalisms differ mainly in their domain of application: QMOR is well-suited to simulating spatially localized systems at high temperature (e.g., spin systems) while DFT is particularly well-suited to simulating spatially delocalized systems (e.g., molecules and conduction bands) at low temperature. In the future, as QMOR is extended to delocalized systems while the methods of quantum chemistry are increasingly extended to dynamical systems [28,80], opportunities will in our view arise for cross-fertilization of these two fields, both in terms of fundamental mathematics and in terms of practical applications.

Overview of contrasts between quantum and classical simulation
In aggregate, the formal, numerical, algebraic, and geometric concepts summarized in the preceding sections and in Figs. 1-4 are in many respects strikingly parallel to similar concepts in the computational fluid dynamics (CFD), solid mechanics, combustion theory, and many other engineering disciplines that entail large-scale simulation using MOR. However, it is evident that quantum MOR is distinguished from real-valued (classical) MOR by at least four major differences, which we will now summarize.
1.6.1. The Theorema Dilectum is fundamental and universal The first difference is that the Theorema Dilectum describes an invariance of quantum dynamics that is fundamental and universal. Its physical meaning, as we will see, is that it enforces causality. Nonlinear classical system do not possess any similarly universal invariance, which is in our view a major contributing reason that "developing effective and efficient MOR strategies for nonlinear systems remains a challenging and relatively open problem" [164, p. 20].
Our results, both analytical and numerical, will suggest that noisy quantum systems are fundamentally no harder to simulate than nonlinear classical systems, provided that the Theorema Dilectum is exploited to allow high-fidelity dynamical projection of quantum trajectories onto a reduced-order state-space.
1.6.2. Quantum state-spaces are veiled The second difference is a consequence of the first. As discussed in Sec. 1.3, to fully exploit the power of the Theorema Dilectum we are required to embrace the ultra-orthodox principle of never looking at the quantum state space. Furthermore, when we examine classical state-spaces more closely, we find that they too are encumbered with ontological ambiguities that precisely mirror the "spooky mysteries" of quantum state-spaces. As discussed in Sections 3.2.7, this modern recognition of spooky mysteries in classical physics echoes work in the 1940s by Wheeler and Feynman [196,197].
1.6.3. Noise makes quantum simulation easier The third difference is that higher noise levels are beneficial to QMOR simulations, because they ensure stronger compression onto the GK petals, which allows lower-rank, faster-running GK state-spaces to be adopted. Later we will discuss the interesting question of whether this principle, together with the concomitant principle "never look directly at the quantum state-space," have classical analogs. We will tentatively conclude that the Theorema Dilectum does have classical analogs, but that the power of this theorem is much greater in quantum simulations than in classical ones.
1.6.4. Kählerian manifolds are geometrically special Broadly speaking, Kählerian geometry is to Riemannian geometry what analytic functions are to ordinary functions. This additional structure is one of the reasons why the mathematician Shing-Tung Yau has expressed the view [200, p. 46] "The most interesting geometric structure is the Kähler structure." From this point of view, the geometry of real-valued MOR state-spaces is mathematically interesting, and the analytic extension of this geometry to Kählerian MOR state-spaces is even more interesting.
Let us state explicitly some of the analogies between analytic functions and Kähler manifolds. We recall that generically speaking, analytic functions have cuts and poles. These cuts and poles are of course exceedingly useful to scientists and engineers, since they can be intimately linked to physical properties of modeled systems. Similarly, the GK manifolds that concern us have singularities, as depicted in Fig. 4(D). Physically speaking, they are associated with regions of quantum state-space that locally are more nearly "classical" than the surrounding regions, in the sense that the local tangent vectors that generate highorder quantum correlations become degenerate. It is fair to say, however, that the deeper geometric significance of Kählerian MOR singularities remains to be elucidated.
Just as contour integrals of analytic functions can be geometrically adjusted to make practical reckoning easier, we will see the Theorema Egregium allows the trajectories arising from the drift and diffusion of noise and measurement models to be geometrically (and informatically) adjusted to match state-space geometry, and thereby improve simulation fidelity, efficiency, and robustness.
More broadly, Yau notes [200, p. 21]: "While we see great accomplishments for Kähler manifolds with positive curvature, very little is known for Kähler manifolds [having] strongly negative curvature." It is precisely these negatively-curved Kähler manifolds that will concern us in this article, and we believe that their negative curvature is intimately linked to the presence of the singularities mentioned in the preceding paragraph. We hope that further mathematical research will help us understand these connections better.

The sectional curvature of gabion-Kähler (GK) state-spaces
We will now proceed with a detailed derivation and analysis of our quantum simulation recipes. Our analysis will "unwind" the preceding overview: first we analyze the geometry of Fig. 4, as embodying the algebraic structure of Fig. 3, using the numerical techniques of Fig. 2. Only at the very end will we calibrate our recipes in physical terms, via the quantum physics of Fig. 1.

Quantum MOR state-spaces viewed as manifolds
To construct our initial example of a gabion state-space, we will consider the following algebraic function ψ(c), whose domain is a four-dimensional manifold of complex coordinates c = {c 1 , c 2 , c 3 , c 4 } and whose range in a four-dimensional Hilbert space is the set of points that can be algebraically represented as {c 1 c 3 , c 1 c 4 , c 2 c 3 , c 2 c 4 }. In the notation of Fig. 2 this function is where "⊗" is the outer product. The superscripts on the c i variables are indices rather than powers, as will be true throughout this section. From an algebraic geometry point of view, (3) defines a projective algebraic variety [53] (also called a homogeneous algebraic variety) over variables {ψ i : i ∈ 1, 4} that is specified above in parametric form in terms of parameters {c i : i ∈ 1, 4}. By definition, our example of a gabion state-space manifold is the solution set of this algebraic variety, and thus our state space is an algebraic manifold. Physically speaking, ψ(c) is the most general (unnormalized) quantum state of two spin 1/2 particles sharing no quantum entanglement.
We will now show that this state-space is a Kählerian manifold that has negative sectional curvature (under circumstances that we will describe) and that this property is ben-eficial for simulation purposes (for reasons that we will describe).
Practical computational considerations: We begin by remarking that the basic algebraic construct "arg1 ⊗ arg2" that appears in (3) can be readily implemented by the built-in functions of most scientific programming languages and libraries; for example in MATLAB by the construct "reshape((arg1*arg2')',[],1)" and in Mathematica by the construct "Outer[Times,arg1,arg2]//Flatten". Similar idioms exist for the efficient evaluation of more complex product-sum structures. Although we will not describe our computational codes in detail, they are implemented in MATLAB and Mathematica in accord with the general ideas and principles for efficient addition, inner products, and matrix-vector multiplication that are described by Beylkin and Mohlenkamp [16].
The abstract geometric point of view: From an abstract point of view, the algebraic structure (3) can be regarded as a sequence of maps where C is the manifold of complex variables {c 1 , c 2 , c 3 , c 4 }, the gabion manifold K is the range of ψ in H, and H is the larger Hilbert space within which K is embedded. To appreciate the surjective and injective nature of these maps, we notice in (3) More generally, it is clear that one coordinate can be set to any fixed nonzero value, without altering ψ(c), by an appropriate rescaling of the other three variables. In our example (3), the dimensions of the three manifolds C, K, and H are therefore where the factors of two arise because these are complex manifolds. We see that the map C → K is surjective (because dim K < dim C), while K → H is injective (because K is immersed in H).
2.1.1. Defining gabion pseudo-coordinates We will call the variables {c 1 , c 2 , c 3 , c 4 } pseudocoordinates. They are not ordinary coordinates because C → K is surjective rather than bijective, or to say it another way, open sets on C are not charts on K. Whenever we require an explicit coordinate basis, we can simply designate any one c k to be some arbitrary fixed (nonzero) value, and take the remaining {c i : i = k} to be coordinate functions.
In practical numerical calculations-where these algebraic structures are called "separated representations," "matrix product states," or "Slater determinants"-pseudocoordinate representations are adopted almost universally. Therefore, we will sometimes simply call the c's "coordinates"; this will make it easier to link the numerical algorithm of Fig. 2 to the geometric properties of K.

Regarding gabion manifolds as real manifolds
Now we will begin analyzing in detail the curvature of the gabion manifold K. For geometric purposes it is convenient to regard H not as a complex vector space, but as a Euclidean space, such that ψ is a vector of real numbers that in our simple example has the eight components {ψ m } = { (ψ 1 ), . . . , (ψ 4 ), (ψ 1 ), . . . , (ψ 4 )}. Similarly, we specify real coordinates on C via c k = x k + iy k , and with a long-term view toward interfacing with the Kähler geometry literature we agree to specify these real coordinates in the conventional order {x 1 , . . . , x 4 , y 1 , . . . , y 4 } ≡ {r 1 , . . . , r 8 } = {r a }. Thus {∂/∂r a } is a complete set of vectors on C.
2.2.1. Constructing the metric tensor Then the map ψ : C → H induces a metric tensor g upon C via the Euclidean metric of H. The components of g evidently are This also suffices to define g as a metric tensor on K, provided we restrict our attention-as for MOR purposes we always will-to functions on K having the functional form f (ψ(c(r))), such that the tangent vectors {∂/∂r a } always act either directly or indirectly (via the chain rule) upon ψ(c(r)). Then knowledge of g allows us to compute via (2) the velocities and path lengths of arbitrary trajectories on K, as is required of a metric on K.

2.2.2.
Raising and lowering the indices of a pseudo-coordinate basis Considered as a covariant matrix, the indices of g ab range over an over-complete basis set, and therefore g ab is singular. It follows that we cannot construct a contravariant matrix g ab in the usual manner, by taking a matrix inverse of g ab . To evade this difficulty we define the contravariant metric tensor to have components g ab ≡ (g ab ) P , where " P " is the same Moore-Penrose matrix pseudoinverse that appears in Step B.3 of Fig. 2, and that was discussed following Eq. (1). It is easy to verify that g ab and g ab act to raise, lower, and contract tensor indices in the usual manner, with a single important difference: the operation of raising followed by lowering is no longer the identity operator, but rather is a projection operator, in consequence of the general pseudoinverse identity (X P X) 2 = X P X. Physically this projection annihilates tangent vectors on K whose length is zero.

2.2.3.
Constructing projection operators in the tangent space From these identities it follows that at a specified point ψ K of K, the local operator P K (ψ K ) that projects vectors in H onto the tangent space of K at |ψ is In the interest of compactness, we will often write P K rather than P K (ψ K ). We readily verify that the projective property P K P K = P K follows from the definition of g ab given in (6) and the general pseudo-inverse identity X P XX P = X P . The ability to construct the projection P K solely from tangent vectors and the local metric tensor g will play a central role in our geometric analysis of K.

"Push-button" strategies for curvature analysis
At this point we can analyze K's intrinsic geometry by either of two strategies. The first strategy, which can be wholly automated, is to fix in our example problem (say) r 7 = 1 and r 8 = 0 so that the remaining {r 1 , r 2 , . . . , r 6 } can be regarded as conventional coordinate functions on the six-dimensional gabion K. The now-restricted set of tangent vectors associated with {r 1 , r 2 , . . . , r 6 } constitutes a coordinate basis, such that (6) specifies the metric tensor for this basis. By construction, this metric has no null vectors, and hence g ab is invertible.
The intrinsic geometric properties of K can then be automatically computed by any of the many symbolic manipulation packages that are available for research in general relativity. This automated approach allows us to "push the button" and discover that for our simple example the scalar Riemann curvature R of K is given by the remarkably simple expression R = −8/(ψ·ψ).

The deficiencies of push-button curvature analysis
What is unsatisfying about an automated coordinate-based analysis, however, is that this simple integer result is obtained as the result of seemingly miraculous cancellations of high-order polynomials. This method produces no insight as to why such a simple integer result is obtained, or why the sign of the curvature is negative, or whether this simplicity is linked to K's ruled structure.
Another objection to coordinate-based analysis is that it forces us to "break the symmetry" of the coordinate manifold C by designating arbitrary fixed values for arbitrarily selected gabion coordinates. This is undesirable because our quantum simulation algorithms in Figs. 1-2 do not break this symmetry. To do so in our geometric analysis would unnecessarily obstruct our goal of linking quantum simulation physics to the Kählerian geometry of K.
We will therefore develop a Riemannian/Kählerian curvature analysis of the gabion manifold K that fully respects the algebraic symmetries, not of K, but of C. For this purpose the sectional curvature proves to be an ideal mathematical tool.
2.4. The sectional curvature of gabion state-spaces Because K has a natural embedding in the Euclidean manifold H, our analysis of sectional curvature is able to follow quite closely the embedded geometric reasoning of Gauss' original derivation of the Theorema Egregium [86]. This approach has the advantage of yielding immediate physical insight. Equally important, this approach can be readily adapted to accommodate the pseudo-coordinate tangent basis that is most natural for analyzing gabion geometry.
As depicted in Fig. 4(E), we choose an arbitrary point on K and define tangent vectors U and V on K to be directional derivatives Because the map C → K is surjective, the representation of U and V as a sum over components u a and v a is nonunique, and we will take care to establish that our sectional curvature calculations are not thereby affected.

Remarks on gabion normal vectors
Conjugate to the tangent space of K at a given point is the space of vectors in H that are normal to the tangent space. We specifyn to be an (arbitrarily chosen) unit vector in that normal space, i.e., to be a vector satisfyinĝ n ·n = 1 andn · ∂ψ(r) ∂r a ≡n · ψ ,a = 0 , where we have adopted the usual notation that a comma preceding a subscript(s) indicates partial differentiation with respect to the indexed variable(s). The sign ofn will not be relevant. We remark thatn is not unique because the codimension of K (by definition codim K ≡ dim H − dim K) is in general greater than unity. Looking ahead, in some of our large-scale numerical examples codim K will be very large indeed, of order 2×2 18 512, 000.

2.4.2.
Computing the directed sectional curvature With reference to the vectors U , V , and n depicted in Fig. 4(E), we define a scalar function S(U, V,n), which we will call the directed sectional curvature, to be Here |. . .| denotes the determinant.
If we recall thatn · ψ ,a = 0, per (9), then it is straightforward to verify that S is a scalar under coordinate transformations. It further satisfies the identity for arbitrary real α, β, γ, and δ. Thus S(U, V,n) is a real-valued geometric invariant of the two-dimensional tangent subspace spanned by U and V .
In the preceding paragraph we emphasize that α, β, γ, and δ are real-valued because later on, when we admit a complex structure, it will not be true that phase-shifting U and/or V leaves the sectional curvature invariant (see Section 2.7).
The denominator of (10) has a simple physical interpretation as the geometric area of the section defined by U and V ; this quantity is often written as |U ∧ V | 2 . In terms of the 2.4.3. Physical interpretation of the directed sectional curvature For a two-dimensional surface embedded in a three dimensional space-the case considered by Gauss-the above expression reduces to the familiar expression S(U, V,n) = 1/(R 1 R 2 ), where R 1 and R 2 are the principal radii of curvature of the surface. In higher dimensions S(U, V,n) describes the Gaussian curvature of a two-dimensional section of K-a two-dimensional submanifold that is locally tangent to U and V -that has been projected onto the three-space spanned by {U, V,n}.
For MOR purposes, this means that whenever S(U, V,n) is negative we are guaranteed local concavity of the state-space K as viewed alongn (i.e, as viewed "from above") along at least one curve that is locally tangent to some linear combination of U and V . The resulting physical picture is shown in Fig. 4(G). For MOR purposes our goal will be, therefore, to choose state-space manifolds such that S(U, V,n) is negative, in the expectation that the associated local concavity of K will improve the robustness of projective model order reduction.

Definition of the intrinsic sectional curvature
Mathematicians usually prefer to describe the curvature of K in intrinsic terms. To accomplish this it is convenient to sum over a complete orthonormal set {n i } of vectors tangent to K. We use the identity in ⊗n =P K = I − P K , where I is the identity operator and P K is the projection matrix given in (7), to obtain Now all reference to unit normals has disappeared, because we have already established in (7) that P K can be described in intrinsic terms. However, the above expression still refers to the embedding Hilbert space via ψ. It will not be until later on (specifically, following (26)) that we show that S(U, V ) is determined solely by the metric tensor and its derivatives.

The formal definition of a gabion manifold
For general tangent vectors U and V , the directed sectional curvature S(U, V,n) and the intrinsic sectional curvature S(U, V ) can be either positive or negative. We will now derive a condition on U and V which is sufficient for S(U, V,n) and S(U, V ) to be nonpositive, and we will use this condition to motivate a formal definition of a gabion manifold. We recall that K is a reduced-dimension state-space manifold that is embedded in a larger-dimension Euclidean manifold H. Thus each individual component of the statevector ψ of H defines a scalar function on K. If we wish, we may regard the metric g of (6) and the normal vectorn of (9) as intrinsically defined in terms of the scalar functions ψ; this eliminates any formal reference to the embedding manifold H. We define a rule vector field, or simply rule field, to be any vector field V on K satisfying The motivation for this definition is simply that the above equation is both intrinsic and geometrically covariant, and furthermore, it is manifestly satisfied by the vector field that is associated with each gabion pseudo-coordinate; these coordinates thus are canonical examples of rule fields. We define rule lines, or simply rules, to be the integral curves of a rule field. It is straightforward to show that rule lines are geodesics; this formally justifies our earlier depiction of rules as "straight lines" in Fig. 4 and in Section 1.5. We define rule tangent vectors, or simply rule vectors, to be vectors that are locally tangent to a rule line. We now formally define a gabion manifold as follows: Definition 2.1. A gabion manifold is a manifold endowed with rule fields whose rule tangent vectors constitute a local basis at every point of the manifold.
This definition is more restrictive than the usual definition of a ruled manifold [71], in which there is no requirement that the rule vectors provide a local basis. Roughly speaking, therefore, a gabion manifold is a ruled manifold that is exceptionally rich in rule structure. Associated with a rule vector V , we define local rule coordinates such that ∇ V ∇ V ψ = 0 takes the component form a,b v a v b ψ ,ab = 0. Thus gabion pseudo-coordinates are local rule coordinates. Evaluating (13) in local rule coordinates, we see that whenever either U or V is a rule vector, the first numerator term vanishes, and the remaining numerator term is nonpositive. This proves the following theorem: Theorem 2.1. Let U be a rule vector at an arbitrary point on a manifold K, let V be an arbitrary tangent vector at that same point, and letn be an arbitrary unit vector normal to the tangent space. Then the directed sectional curvature satisfies S(U, V,n) ≤ 0 and therefore, the intrinsic sectional curvature satisfies S(U, V ) ≤ 0.
In physical terms, any two-dimensional section of K that includes a rule vector has negative sectional curvature. Since for gabion manifolds, the local rule tangents form an over-complete local basis at each point in K, we see that negative sectional curvature is ubiquitously present in our gabion state-spaces.

Recipes for constructing rules and rule fields
In the context of MOR analysis, rule lines are easy to construct, and they have a clear algorithmic significance.
Rules can be readily constructed via the product-sum algebraic structure of Fig. 3, by varying any one { l c m n } while holding the others fixed. The tangents to the rule lines then constitute an overcomplete local basis, as depicted in Fig. 4 (C).
More generally, rule fields can be constructed by selecting an arbitrary order (column) in the product-sum algebraic structure of Fig. 3, selecting arbitrary basis vectors for the Hilbert subspace associated with that order (equivalent to imposing an arbitrary rotation on the basis vectors of that column's subspace), selecting an arbitrary rank (row), selecting an arbitrary element of the substate of that order and rank, choosing a coordinate system (in the strict sense of Section 2.3) such that the selected element is one of the coordinate functions, and identifying a rule field on K with the partial derivative with respect to that coordinate.
From an algorithmic point of view, a rule vector is a direction in the state-space along which trajectories can move with great algorithmic efficiency, since only one state-space coordinate need be updated.

The set of gabion rules is geodesically complete
Whenever any two rows of the product-sum in Fig. 3 are algebraically degenerate, the tangent space of K will be geometrically degenerate, yet according to the construction of the geometric rules, and in particular because the underlying algebraic structure is polynomial, the rules pass through curvature singularities without disruption of their geodesic properties, as depicted in Fig. 4

(D).
Furthermore, it is evident that the sectional curvature (13) diverges in the neighborhood of a rule singularity, in consequence of the divergence of the pseudo-inverse metric g ab that appears in the projection operator P K as given in (7). Numerical experiments confirm this expectation. These phenomena suggest that gabion manifolds might fruitfully be analyzed in terms of affine algebraic varieties. The authors have not pursued this line of analysis.

Gabion-Kähler (GK) manifolds
We now specialize (13) to complex manifolds having a Kählerian metric, which we will call gabion-Kähler (GK) manifolds. In so doing, we will adopt certain "tricks" of indexing that Kählerian geometers use. We begin by writing the numerator and denominator of (13) in matrix notation, with all quantities still real: . Now we reason as follows. S(U, V ) is a real number that is independent of coordinate system. We are therefore free to analytically continue our coordinates, transforming (for example) the coordinate pair {x 1 , Since the sectional curvature is a geometric invariant, the (real) value of S(U, V ) will not be altered thereby, even though the coordinates themselves are now complex.

Kählerian indexing and coordinate conventions
It is evident that analytic continuation to complex coordinates treats c 1 andc 1 as independent coordinates for symbolic manipulation purposes (such as partial differentiation), just as x 1 and y 1 are independent coordinates. It is only at the very end of a calculation, when we assign (complex) numerical values to c 1 andc 1 , that they can no longer be varied independently. It is helpful too to replace Latin indices with unbarred and barred Greek indices, in a convention that associates barred indices with barred coordinates (see [74, p. 8] or [135]). Then the vector V has components {v 1 , v 2 , . . . ,v1,v2, . . .}, for example, and is represented in terms of partial derivatives by In this convention a barred indexk ≡ dim C/2 + k, such that the index1 is a shorthand for the integer dim C/2 + 1, and (v α ) =vᾱ. With regard to the embedding Hilbert space H, we will adopt the physics convention that the "ket" vector ψ(c) ≡ |ψ(c) is a complex vector of dimension dim H/2, with the "bra" vectorψ(c) ≡ ψ (c)| being the conjugate vector.
Thus ψ(c) is a holomorphic function (also known as an "analytic function") of complex pseudocoordinates c, andψ(c) is similarly a holomorphic function ofc. Defining the biholomorphic Kähler potential function κ(c, c) to be the components of the metric tensor g are given in terms of κ(c, c) by It immediately follows that in any holomorphic coordinate system g αβ = gβ α and g αβ = gᾱβ = 0 .
The simplest explicit example of this convention is the complex plane regarded as a two-dimensional Kähler manifold. Indexing its two coordinates by {1,1} yields coordinates {c 1 ,c1}, which in analytic function theory are conventionally called {z,z}. The real-valued Kähler potential is κ = c 1c1 /2 = zz/2, the components of the metric tensor are g ab = [ 0 1 1 0 ]/2 and g ab = 2 [ 0 1 1 0 ], and the length element ds 2 is which is the usual normalization. Note the ubiquitous factors of 2 and 1/2, which require careful attention in practical calculations. A general scalar function on this manifold is of the form f (z, z), and functions of the special form f (z) are the holomorphic (or analytic) functions.
In programming calculations on Kähler manifolds of larger dimension, it is helpful that (18) and (19) take the block-matrix forms where "[. . .]" is a square matrix. Environments like MATLAB and Mathematica provide built-in functions for block-matrix constructs of this type. We further see that in consequence of g ab = g ba and g αβ =ḡᾱ β , which follow from (18), the individual block matrices in (21) are Hermitian and semipositive; for this reason the submatrix g αβ is sometimes called the Hermitian metric of the complex manifold. In practical calculations it is considerably more efficient to work solely with the Hermitian metric and its pseudo-inverse, than with the larger matrix g ab . Numerically-minded readers who are new to the literature of Kähler manifolds will appreciate that the above indexing conventions elegantly resolve a contradiction of our intuitions. On the one hand, we expect that a coordinate transformation cannot change the range of an index. On the other hand, we expect on physical grounds that a manifold described by complex coordinates will need only half the number of coordinate variables as the same manifold described by real coordinates. The resolution of this dilemma is in the block structure and symmetry properties of g, which ensure that in practical geometric calculations, only half the index range need be summed over.
2.6.2. GK sectional curvature in physics bra-ket notation It is then a straightforward exercise to write (15a-15b) compactly, in the bra-ket notation of physics: where the partial derivatives and the projection operatorP K ≡ I − P K are given in terms of components by The above expressions show explicitly that bra-ket notation allows the sectional curvature to be computed by summing over half-ranges of coordinate indices. This notational compactness constitutes (from a Kählerian geometry point of view) the main practical rationale for the bra-ket notation of the physics literature.
2.6.3. Defining the Riemann curvature tensor Now we define the Riemann curvature tensor to be that scalar function R(A, B, C, D), defined at each point on the gabion manifold K, with A, B, C, D being arbitrary vectors, such that the local sectional curvature and the local Riemann curvature are related by It is known (see [82,Theorem 3.8] or [135,Theorem 7.51]) that the Riemann curvature so defined is unique, for both real and Kähler manifolds, provided that the following index symmetries are imposed: which are the conventional antisymmetries of the Riemann tensor, and provided in addition the following identity is satisfied which is called the first Bianchi identity.
On real manifolds, the implicit definition (23) of the Riemann curvature in terms of the sectional curvature is difficult to work with, in the sense that the general expression for R abcd as a function of g turns out to be too complicated to readily derive by inspection or manipulation of (13).
Fortunately, on Kähler manifolds we have the simpler definition of S(U, V ) given in (22a-22d), from which the Riemann curvature can be read off as where κ(c, c) is the biholomorphic Kähler potential introduced in (17) and g µν is the pseudoinverse of g µν = κ ,µν introduced in (21). We further specify that all components of R not fixed by (26) vanish, save those required by the symmetries (24). Since the following components of R cannot be obtained from (26) by symmetry, we take them to be zero and we see that the resulting Kählerian R abcd has a block structure similar to that of the Kählerian metric g ab in (21). In consequence of this block structure, it is straightforward to verify that the Bianchi identity (24) is equivalent to the following Kählerian Bianchi index symmetries (which we have not found explicitly given in the literature): and which (26) respects. Thus the definition, symmetries, and identities (23)(24)(25) all are satisfied, and we conclude that the Kählerian sectional curvature (22a-22d) uniquely determines the Riemann curvature to be (26). We remark that the elegant functional simplicity of the Kähler-Riemann curvature tensor (26), which in many textbooks mysteriously appears only at the end of a long algebraic derivation, emerges quite simply and naturally in our Gauss-style, immersive, bra-ket derivation.
2.6.4. The Theorema Egregium on GK manifolds Because S(U, V ) and R(A, B, C, D) depend solely on the intrinsic metric g, we have thus derived-solely by analysis of sectional curvature-the celebrated Gauss/Riemann Theorema Egregium as it applies to Kähler manifolds. As mentioned above, we have not found in the Kählerian geometry literature any similar derivation of the Riemann tensor by sectional curvature analysis. However, the Kähler geometry literature is so vast, that we can say with confidence only that the above derivation of (26) is quite different from the usual (intrinsic) derivations in the literature (see [135,Theorem 12.5.6] and [142,Ch. 6]). We remark also that minor misprints are commonplace in the mathematical literature (e.g., [120, eq. 7.11 of ch. 1] is identical to (26) save for an incorrect index-pair contraction). That is why we will derive some closed-form results in Section 2.8 against which large numerical codes can be compared.
To better serve our MOR purposes, we have generalized R αβγδ by allowing its indices to range over dim C rather than dim K, and defined the contravariant tensor g µν in terms of the matrix pseudoinverse. These extensions do not alter the functional form of (23-28), and they make practical numerical calculations very much easier to program, as we will see in Section 2.9.
It is apparent that the generalized metric tensor g µν that appears in the Riemann curvature tensor (26) is identical to the matrix (∂⊗∂κ) P that appears in the simulation algorithm of Step B.3 of Fig. 2. This is the first of many links that we will establish between Kählerian geometry and quantum simulation.

Readings in Kählerian geometry
Having derived the Kählerian Riemann curvature tensor, we will attempt a brief survey of the Kähler geometry literature. The literature on Kähler geometry is comparably vast to the literature on quantum measurement and information, and so our review necessarily will be exceedingly sparse and subjective.
Our development and indexing conventions in this section have paralleled the Riemannian conventions of Weinberg [195], as extended to Kählerian geometry by Flaherty [74], as further extended to abstract notation by Martin [135] and Moroianu [142]. However, our analysis has been centered upon sectional curvature, rather than Riemann curvature as in the preceding texts. Other suitable texts include Kobayashi [120], Frenkel [77], Gallot [82], Flanders [75], Jost [115], Hou and Ho [111], and a lengthy review by Yau [200]. Many more textbooks and review articles on Kähler geometry exist, and it is largely a matter of individual taste to chose among them.

Remarks upon holomorphic bisectional curvature
We have seen that the sectional curvature of both real and Kählerian gabion manifolds is constrained by Theorem 2.1. We now review additional (known) sectional curvature theorems that relate particularly to gabion-Kähler manifolds. These theorems concern a geometric measure called the holomorphic bisectional curvature [91].
We begin by remarking that physicists in particular are accustomed to thinking of ψ and iψ as being physically the same vector. In Kählerian notation, the notion of "multiplying by i" is associated with an almost complex structure J, which is defined to be a linear map J, defined on every point in K, satisfying J 2 V = −V for V an arbitrary tangent vector, and differentially smooth, that leaves the metric tensor invariant, i.e., g(U, V ) = g(JU, JV ). The effect of J upon the components {v α ,vᾱ} of an arbitrary vector V is simple: The functional form of (22a-22b) then implies the identities S(U, V ) = S(JU, JV ) and S(JU, V ) = S(U, JV ). We now consider the sum S(U, V ) + S(U, JV ) = S(U, V ) + S(JU, V ), which by definition is called the holomorphic bisectional curvature [91]. The fundamental inequality of holomorphic bisectional curvature follows immediately from (22a-22b), because the "multiply by i" rule implies that the first term in the denominator cancels in the sum. We will call (30) the holomorphic bisectional curvature nonpositivity theorem (HBCN Theorem). The HBCN Theorem is a long-known result [91] that applies in general to any Kähler manifold that is a complex submanifold of a Euclidean space. Physically speaking, if a given Kählerian section {U, V } has positive curvature, then the "rotated" section {U, JV } will have negative curvature. It is readily shown [91] that the HBCN Theorem implies the nonpositivity of the eigenvalues of the Ricci tensor, and therefore, the nonpositivity of the scalar curvature.
As a technical point, the HBCN Theorem applies only to Kähler manifolds that have (or can be given) a complex embedding in a larger Euclidean manifold, which is the case of greatest interest in quantum MOR applications. The HBCN Theorem does not apply to Kähler manifolds that have no complex Euclidean embedding. For example, Kähler manifolds having a Fubini-Study metric can (and sometimes do) exhibit positive scalar curvature. This is incompatible with the HBCN Theorem, and these manifolds therefore have no complex Euclidean embedding.
The functional form of (22a) allows us to immediately extend the HBCN Theorem to encompass the directed sectional curvature S(U, V,n) as follows: Theorem 2.2 (directed extension of the HBCN Theorem). At an arbitrary point on a Kähler manifold K that has a complex embedding within a Euclidean space H, let U and V be arbitrary tangent vectors and letn be an arbitrary unit vector normal to K. Then the directed sectional curvature satisfies S(U, V,n) + S(U, JV,n) ≤ 0.  (14) that if W is a rule field, then so is JW, hence if V is a rule vector, then so is JV . It then follows from Theorem 2.1 that S(U, V ) and S(U, JV ) are both nonpositive, which is a strictly stronger condition than the nonpositivity of their sum that is implied by the HBCN Theorem.
It further follows immediately from (15a-15b) that if V is a rule vector, then S(U, V ) = S(U, JV ). More generally, if U and V are both rule vectors, then S(U, V ) = S(U, JV ) = S(JU, V ) = S(JU, JV ), and all of these sectional curvatures are nonpositive. From this "gabionic" point of view, we see that on gabion-Kähler manifolds the rule vectors U and JU are effectively the same vector, i.e. rule vector phases are irrelevant to sectional curvature properties, as is natural in quantum mechanical analysis.
The ubiquity of nonpositive sectional curvature on GK manifolds therefore can be viewed as arising from the confluence of rule structure and complex structure, both of which are associated with strong theorems that imply negative sectional curvature.

Practical implications of sectional curvature theorems
From an MOR point of view, the nonpositive sectional curvature implied by Theorem 2.1 can be regarded as helping to ensure the robustness of MOR on both real and complex manifolds, while the HBCN Theorem helps makes MOR even more robust on complex manifolds. This is one of two fundamental reasons why quantum MOR can be regarded as intrinsically easier than classical MOR (the other reason being the algorithmic resources that are provided by the Theorema Egregium, as we will discuss in the following section).
We remark that not all sectional curvatures S(U, V ) are negative on GK manifolds. In small-dimension GK manifolds we have constructed analytical examples of positivecurvature sections-with some trouble because neither U nor V can be rule vectors-and in large-dimension gabions numerical searches find them. According to our present (limited) understanding, these positive-curvature Kählerian sections seem rather artificial, and so we will not discuss them further. It is possible that their significance has eluded us.

Analytic gold standards for GK curvature calculations
The Riemann tensor specifies the Gaussian curvature of every section of K, and on largedimension MOR manifolds the Riemann tensor therefore carries a vast amount of information. To compress this information (among other purposes) it is conventional to condense the four-index Riemann tensor R abcd into the two-index Ricci tensor R ab by which in Kählerian index notation can be written Further compression is achieved by the scalar curvature R defined by Caveat: various authors entertain diverse conventions regarding which indices should be contracted to obtain the Ricci tensor, and the factors of -1 and 2 that appear above are commonly associated with minor errors and imprecisions. There is one sectional curvature convention, however, that is universal: the directed sectional curvature on a unit hypersphere S n (the ordinary sphere embedded in R 3 being S 2 ) satisfies S(U, V,n) = 1 for all linearly independent U and V . In a locally orthonormal basis there are n(n − 1) such pairs; it follows that on a unit hypersphere the eigenvalues of the Ricci tensor are all n − 1, and the scalar curvature itself is therefore R = n(n − 1). This result is a useful "gold standard" for testing symbolic and numerical calculations on real manifolds.
To construct a similar gold standard for testing symbolic and numerical curvature calculations on Kählerian manifolds, it is convenient to consider the manifold of rank-one, order-n gabion states (see Fig. 3). We allow the spin quantum numbers {j i : i ∈ 1, n} associated with successive product subspaces to vary independently. Then by a straightforward (but not short) calculation, it can be shown that the scalar curvature of a general rank one, order-n product state is where κ is the Kähler potential function of (17). See the following section for a summary of the assorted algebraic techniques used to obtain this result.
To the best of our knowledge, this general analytic result has not previously appeared in the literature. For the simple example manifold of (3), we have order n = 2 and j 1 = j 2 = 1/2, so the above yields R = −4/κ = −8/ ψ|ψ = −8/(ψ·ψ), in agreement with the automated "push-button" analysis of Section 2.3.

The Riemann-Kähler curvature of Slater determinants
We now have all the tools we need to compute the scalar Riemann curvature of Slater determinant states, which are the main state-space of quantum chemistry [54,80,184]. The strategy of the calculation is straightforward, and although the details are lengthy, the final result is simple.
We begin by considering, as the simplest example having nontrivial curvature, the quantum states that are obtained by antisymmetrizing the outer products of a rank-one order-two product-sum state of spin j = 3/2 (see Fig. 3) In the language of quantum chemistry, we can equivalently regard ψ(c) as the variational state-space of the (unnormalized) Slater determinant states of two electrons "a" and "b" occupying linear combinations of four orbitals. With ψ(c) given, we can compute its scalar curvature R by either numerical or analytic means. To compute R numerically, we can simply set the eight pseudo-coordinates to any desired value (thereby choosing the point in state-space at which R is to be evaluated) then compute first the Kähler potential κ from (17), then the Riemann curvature R from (26), and (31)(32)(33), evaluating the pseudo-inverse metric g µν of (26) numerically.
Empirically we find that for our simple example (35) these numerical calculations invariably yield a Riemann scalar curvature of R = −8/κ (to machine precision) for all values of the pseudo-coordinates {c i a , c i b }. The simplicity of the numerical result motivates and guides the following analytic evaluation of R in closed form. We designate by ψ 0 the point in the state-space ψ(c) at which the Riemann curvature is to be evaluated. We write ψ 0 in the following standard form by an appropriate choice of basis vectors An arbitrary state in a neighborhood of ψ 0 can be written as an explicit holomorphic function of five complex coordinates That the above variety is in fact equivalent to (35) can be demonstrated by the "pushbutton" method of computing their respective Gröebner bases [53] and verifying that these bases generate the same homogeneous ideal in the variables{ψ i }. 2 In the transition from (35) to (37), four pseudo-coordinates {c 1 a , c 2 a , c 1 b , c 2 b } have disappeared and been replaced by the "0" and "1" entries. These are physically accounted as follows: c 1 a and c 2 b have been merged into the overall (complex) normalization c 0 by a simple rescaling; c 1 b and c 2 1 have tangent vectors that vanish at ψ 0 , such that they do not induce coordinate charts on ψ(c) at ψ 0 and can safely be dropped. In numerical calculations, their tangent vectors at ψ 0 are in the null-space of g αβ and g αβ .
The remainder of the analytic calculation is straightforward. The metric tensor components g µν are given from (37) by (17)(18), and when evaluated at ψ 0 , yield a 5×5 matrix that is diagonal and nonsingular. Computing the inverse gν µ is therefore trivial. The Riemann components R αβγδ are given by (26) and the Ricci components R µν are given by (32).
Finally, the scalar Riemann curvature given by (33) is found to be R = −8/κ, in accord with the numerical result.
This reasoning is readily generalized. Upon antisymmetrizing a general rank-one, ordern product-sum state (see Fig. 3) under the exchange of all pairs of spins-thus converting it to a single n-particle Slater determinant-and evaluating the metric and Riemann tensors in the diagonal basis of (36)(37), the dimensionality of the Kähler manifold is evidently dim K = 2(1 + n(n orb − n)), and the scalar curvature is found in closed analytic form to be Here n orb = 2j + 1 is the dimension of the individual states in the Slater product-sum, with the mnemonic "orb" referring to "orbitals" in token of the practical use of these states in quantum chemistry. Since the scalar curvature is a geometric invariant, this result holds in any basis. The curvature of a Slater determinant manifold having a Fubini-Study metric, i.e., having a Kähler potential κ = 1/2 log ψ |ψ instead of κ = 1/2 ψ |ψ , can similarly be calculated from (37) and its multidimensional generalizations. The result is Physically speaking, the Fubini-Study metric describes a manifold of normalized states in which a phase rotation of δφ radians has a path length of zero rather than δφ. Furthermore, it is a straightforward (but lengthy) algebraic exercise to verify during the course of the above calculation that g αβ ∝ R αβ , i.e., Slater determinant manifolds having a Fubini-Study metric are Ricci-Einstein manifolds [14]. An immediate consequence is that Slater determinant manifolds are solitons under Ricci flow [49]. We note that Ricci soliton manifolds are of central importance to the mathematical community, because they represent (in sense that can be made precise) the unique "smoothest" manifolds of a given topological class. In recent years this idea has been central to the Ricci flow [50] proofs of several long-standing topological conjectures.
To the best of our knowledge, the above algebraic/geometric properties of Slater determinants have never been noted in the literature. This gap is noteworthy, in view of the central role that Slater determinants play in chemistry and condensed-matter physics, and the similarly central role of Kählerian algebraic geometry in numerous branches of mathematics. It reflects what the National Research Council has called [149] "[. . . the anomaly that] although theoretical chemists understand sophisticated mathematics and make heavy use of the mathematical literature, they have typically not involved mathematicians directly in either the development of models or algorithms or the derivation of formal properties of equations and solutions. In fact, theoretical chemists have become accustomed to self-reliance in mathematics." A central objective of this article's geometric approach to quantum simulation is to more closely link the formal mathematical tools of algebraic geometry to practical problems in quantum simulation, and thereby, to help ensure that the emerging discipline of quantum system engineering is not needlessly "self-reliant in mathematics."

Slater determinants are Grassmannian GK (GGK) manifolds
We thank our colleague Joshua Kantor for directing our attention to the facts that the Slater determinants associated with n particles distributed among n orb orbitals have a natural isomorphism to the Grassmannian manifolds classified as G(n, n orb ) [100], that Grassmannian manifolds have a known presentation as projective algebraic varieties via the Plücker embedding [71], and that the Plücker embedding is possessed of isometries having a Lie group structure that is known to be compatible with a Kähler-Einstein metric [14].
From this algebraic geometry point of view, the example of a Slater determinant that we present in (35), having n = 2 and n orb = 4, can be identified with what Harris' classic textbook Algebraic Geometry calls [100, p. 65] "the first nontrivial Grassmannian-the first that is not a projective space-[namely] G(2, 4)." Further investigation of this convergence of algebraic geometry with quantum chemistry is in progress.

Practical curvature calculations for QMOR on GK manifolds
Now our appetite is whetted to ask even more difficult questions: what are the curvature properties of higher-rank GK manifolds? And what are the implications of these curvature properties for quantum simulation? We will turn to numerical experiments to learn more about these issues.
As we begin numerical curvature calculations on Kähler state-spaces of increasingly large dimension, we first consider what we may expect to learn from these calculations that would have practical implications for MOR.
It is known [78,144] that at each point on K we can construct locally Euclidean Riemann normal coordinates such that the local curvature tensor has the expansion and the local volume element dV ∝ | det g | 1/2 therefore has the expansion We see that the square roots of the eigenvalues of the Ricci tensor determine the length scale over which curvature effects dominate the volume of the state-space. Physically speaking, this sets the length scale over which gabion petals are locally Euclidean. We further see that negative curvature is associated with an exponential "flowering" of the state-space volume, in accord with the beautiful knitted representations of hyperbolic spaces by Taimina and Henderson [109]. We note that since the Ricci tensor eigenvalues are geometric invariants (meaning specifically, the eigenvalues of the mixed Ricci tensor R a b = c g ac R cb are geometric invariants), we can compute them in any coordinate system we please. We further note that our adoption of a over-complete Kählerian basis creates no anomalies in the Ricci eigenvalue distribution, because the extra eigenvalues of R α β = γ g αγ Rγ β vanish identically, due to the projective definition (21) of g αγ .
Our first goal, therefore, will be to compute the Ricci tensor eigenvalues for high-rank GK manifolds, expecting thereby to gain quantitative insight into the "flowering" gabion geometry that is depicted in Fig. 4(H).
As our first test case, we will consider the rank-6, order-9, spin-1/2 GK manifold. This manifold has (real) dimension dim K = 2 × 9 × (6 + 1) = 126 and it is embedded in a Points on the gabion-Kähler manifold are randomly selected by first randomly generating an independent normalized product state for each gabion rank (i.e., each row of Fig. 3), then summing the ranks, then normalizing the final state. The sparsely dotted line contains typical Ricci eigenvalues of the spin-1/2, rank-9, order-6 gabion, the densely dotted line contains the eigenvalues of the spin-1/2, rank-9, order-10 gabion. The dashed line is an empirical "rank×order" estimate of the mean Ricci eigenvalue of gabion-Kähler manifolds having large-codimension.
Choosing a random point in K and computing the Ricci eigenvalues numerically yields the typical eigenvalues shown in Fig. 5. We take the square root of the largest eigenvalue to be (roughly) the linear extent of a gabion petal; these petals evidently have an extent of order 0.01.
We emphasize that K is large enough to contain exponentially many such petals. See Nielsen and Chuang for a quantitative analysis [146,Section 4.54], noting that the "patches" of Nielsen and Chuang are broadly equivalent to our petals and therefore, their Fig. 4.18 is broadly equivalent to our Fig. 4(H). This mathematically justifies our physical picture of a gabion as a geometric "flower" having exponentially many petals. Later on, in Sections 4.6.4-4.6.6, we will rigorously count the number of petals using ideas from coding theory.
2.12. Numerical results for projective QMOR onto GK manifolds Suppose we generate a random (normalized) quantum ψ 0 , and numerically search for a highest-fidelity projection of ψ 0 onto K. We will call this high-fidelity image point ψ K .
Based on our geometric analysis so far, we have two strong expectations relating to quantum MOR by projection onto K. We expect first, that projections exist for which |ψ K − ψ 0 | is no greater than ∼ 0.01, since at greater separations the negative curvature of K is great enough to ensure a better solution, via the mechanism of Fig. 4(G). Second, we expect that the numerical search for ψ K will be well-conditioned. That is, it will robustly converge to a high-fidelity representation, despite the exponentially convoluted geometry of K, without becoming trapped in local minima (because every local section that contains a rule is negatively curved).
In numerical trials both expectations are fulfilled. The achieved fidelity is in good accord with the geometric expectation of ∼0.01: the median value of |ψ K − ψ 0 | in a trial of 100 projective reductions was 0.005, and the maximal value was 0.025.
We computed the reduced-order GK representation ψ K by a simple gradient search. Specifically, we integrated to convergence the dynamical equation of Step B.3 of Fig. 2, with the potential φ replaced by Upon conversion to Kähler index notation, the resulting reduction equation is manifestly geometrically covariant: From an algorithmic point of view, this reduction equation can be regarded as a downhill gradient search for a reduced-order representation ψ K = lim t→∞ ψ(c(t)) of ψ 0 . The search evidently converges when ψ K is directly "underneath" ψ 0 , i.e, when P K (ψ K ) ψ 0 −ψ K = 0, where we recall that P K (ψ K ) projects vectors in H onto the tangent space of K at ψ K . As an alternative to gradient search, we note that Beylkin and Mohlenkamp [16, sec. 3.1] describe an alternating least-squares algorithm that in their hands gives excellent results, but we have not ourselves tried this method. In Section 4.6.1 we derive the above equation (43) using the language (and nomenclature) of compressive sampling (CS), and we describe the numerical calculations in greater detail.
No attempt was made to optimize the efficiency of the gradient search; instead we simply reused the existing dynamical code implementing Step B.3 of of Fig. 2 (this code exploits the gabion algebraic structure of ψ(c) to evaluate (43) efficiently). This duplication of internal algorithmic structure again illustrates the intimate relation between trajectory calculations and geometric calculations.
No gross failures of convergence, such as might be expected from trapping of the trajectory induced by (43) on a distant gabion "petal," were observed in this (or any) of our numerical trials. Our geometric analysis explains this robustness as originating in the negative sectional curvature of the ruled net of the gabion state-space K, as depicted in Fig. 4(G).
We now increase the order of the gabion to 10, keeping the rank at 9. This increases the dimensionality of the Hilbert space to 2×2 10 = 2048, and the dimensionality of the Kählerian gabion to dim K = 2 × 9 × (10 + 1) = 188. With codim K = 1048 − 188 = 1860, our model order reduction is now discarding ∼ 90% of the dimensions of the larger Hilbert space. Thus, in deliberate contrast to the previous example, the MOR is now fairly aggressive. Typical resulting Ricci eigenvalues are shown in Fig. 5.
A pronounced flattening of the eigenvalue distribution is evident in this large-codim example. The plotted straight line is simply rank × order, which seems empirically to describe the average Ricci eigenvalue in this particular case, and also in many other trials that we have run. This rule-of-thumb is simply the analytic result for rank-1 curvature (34) multiplied by the rank.
There is at present no analytic theory that justifies this rule-of-thumb, or explains the observed flattening of the eigenvalue distribution. Obviously, such a theory would be welcome. We remark that to the best of our knowledge, the Ricci eigenvalues of Fig. 5, having dim K = 188, are the largest-dimension curvature eigenvalues ever numerically computed.
We then calculated projective MOR approximations by integrating (43) as before. We generated our targets ψ 0 by randomly selecting target points on K, and moving a distance of 0.05 along a random vectorn perpendicular to K.
Again, robust convergence to high-fidelity reduced-order representations was observed. Of 100 trials, in 98 cases the separation distance was precisely 0.050, representing convergence to the correct "petal." In the remaining two cases the separation distances were 0.123 and 0.120, representing sporadic convergence to a wrong-but-nearby "petal." Thus even wrong-petal convergence yielded a high-fidelity representation. This robustness can again be ascribed to the negative sectional curvature of the state-space manifold's ruled net, according to the geometric mechanism depicted in Fig. 4(G).

Avenues for research in geometric quantum mechanics
To begin, the results of the preceding section are understood only qualitatively, and rigorous bounds on projective fidelity and robustness would be very welcome. How much of the preceding section can be understood solely as a consequence of the known sectional curvature properties of the GK manifolds?
Although our presentation focuses on practical applications of quantum MOR, a broad class of fundamental physics questions can be given a geometric interpretation. We temporarily adopt the point of view of geometric quantum mechanics in which-as reviewed in Section 1.3-the manifold K is regarded as the "real" arena on which physics takes place.
We begin by noting that that the rule-field equation (14) specifies the dimensionality of Hilbert space to be the number of (linearly independent) scalar rule-fields that the (postulated) "real" GK manifold of geometric quantum mechanics supports. The question then arises, what determines the number of these embedding state-space fields? The present article suggests that this question is best investigated from a blended informatic-algebraicgeometric point of view.
As another example, suppose p and q are arbitrary linear Hermitian operators in the embedding Hilbert space H. At a given point of the Kahler manifold K, specified by a state |ψ(c) , we can construct tangent vectors V q and V p whose components are and similarly for V p . With this normalization we have g(V q , V q ) = ψ |qP K q|ψ . Physically speaking, V q defines a Kählerian velocity field along which the projected dynamical equation ∂|ψ(t) /∂t = P K q|ψ(t) moves state trajectories. The sectional curvature S(V p , V q ) then can be regarded as a fundamental property of the "true" physical manifold K. The Theorema Egregium guarantees that the sectional curvature is an intrinsic property of K, and our physical intuition suggests that it should therefore be measurable.
A host of fundamental questions then arise quite naturally, that intimately unite physics and mathematics in the context of quantum simulation. If the quantum sectional curvature is physically measurable-whether in reality or within projective simulations-then by what kinds of experiment? Are these experiments practical in our real world? Is there any experimental evidence already at-hand to indicate that quantum sectional curvature vanishes in the physical world? If so, to what precision has this been verified?
It is apparent that quantum chemists can apply a reverse strategy to create a geometric context for assessing the fidelity of chemical quantum simulations. Set p to be the kinetic energy of the electrons of a molecule, and q to be the potential energy (including perturbations due to applied potentials), and set the gabion state-space to be a sum of Slater determinants. Then the vanishing of the sectional curvature S(V p , V q ) indicates that the state-space section generated by {p, q} is Euclidean, as is (presumably?) desirable in chemical simulations, most particularly, density functional theory (DFT) simulations.
And finally, to anticipate, in Section 3.2.1 we will discuss how the Theorema Dilectum manifests itself upon GK state-spaces. This will turn out to raise thorny issues of how causality works in geometric quantum mechanics. For the present we will say no more about these difficult fundamental questions, instead referring the reader to Leggett's recent review [126], and we return instead to our central topic of practical quantum spin simulations.

Summary of the geometric analysis
We have seen that the ruled net of a GK manifold, with its associated nonpositive sectional curvature, plays a geometric role in quantum simulation that is broadly similar to the role of polytope convexity in linear programming, namely, the role of providing geometric foundations for developing efficient and robust algorithms. However, understanding polytope convexity has been a slow process, and this research is still being pursued along multiple mathematical lines that include algebra, geometry, information theory, and their numerous hybridizations. We foresee that the role of algebraic geometry in classical and quantum simulations will be elucidated along similarly multidisciplinary lines, and we are conscious that the analysis presented here has barely begun this enterprise.
Even in its present early stage, geometric analysis provides grounds for confidence that quantum model order reduction and trajectory simulation can be achieved with high efficiency, fidelity, and robustness, provided that some physical mechanism can be found to compress quantum trajectories onto the petals of gabion-type state-spaces. That needed physical mechanism is, of course, the Theorema Dilectum, which will be the topic of the section that follows.

Designing and Implementing Large-Scale Quantum Simulations
Our goal in this section is to join the geometric ideas and theorems of the preceding section with the well-known ideas and established design principles of linear quantum mechanics.
We have reason to worry that perhaps very few principles of linear quantum mechanics will survive the transition to reduced-order quantum MOR mechanics. After all, we have seen that gabion-Kähler (GK) state-spaces are strongly curved, so that when we project quantum trajectories onto them, we (seemingly) dispense with all of the mathematical properties of Hilbert space that depend upon its linearity. Furthermore, depending upon the degree of the model order reduction imposed, the GK projection of QMOR mechanics may even discard all but an exponentially small fraction of the dimensions of the embedding Hilbert space.
We seek, therefore, to establish that the principles of linear quantum mechanics hold true in QMOR simulations, and in particular, to establish precisely the mechanisms by which they can hold true.

Organization and nomenclature of the presentation
For convenience, whenever we derive a result that is novel or expressed in a new form, we state it as a formal (numbered) design rule that is accompanied (as needed) by formal (numbered) definitions. On the other hand, whenever a result is unsurprising, or can be found in the literature, and is obtained from previously given equations by standard manipulations and reductions, we outline the derivation but omit details.
We define the subject of our analysis to be QMOR mechanics: We seek to construct recipes by which QMOR mechanics simulates linear quantum mechanics as closely as feasible. To recapitulate the objectives of Section 1.1.1, our analysis seeks to be orthodox in its respect for linear quantum mechanics, operational in the traceability of its predictions to measurement processes, and reductive in the sense that its principles are summarized by closed-form analytic design rules. Our analysis strives also to be synoptic in the sense that whenever we choose between equivalent analysis formalisms, we state a rationale for our choice, and we note the practical consequences of alternative choices.

QMOR respects the principles of quantum mechanics
We begin by establishing that the mathematical structure of QMOR mechanics is sufficiently rich to respect the following principles of quantum and classical mechanics: • the causal invariance of the Theorema Dilectum is respected (Sec. 3.2.1), • the entropy of systems in thermodynamic equilibrium is respected (Sec. 3.2.2), • the principles of classical linear control theory are respected (Sec. 3.2.3).
• the quantum limits to measurement noise and back-action are respected (Sec. 3.2.6). If QMOR mechanics did not respect these principles, it would scarcely be useful for simulating real-world quantum systems. Conversely, by stating these principles in a quantitative, we gain a fairly clear picture of the path that our QMOR analysis needs to take.

QMOR respects the Theorema Dilectum
We now state the Theorema Dilectum in an algebraic form that is well-adapted to the formal quantum MOR algorithm of Fig. 1. Our mathematical analysis parallels that of Nielsen and Chuang [146, see their Theorem 8.2 in Section 8.2], whose analysis in turn derives from a 1975 theorem of Choi [48] (see Section 1.3.4).
We suppose that at the end of step n the QMOR simulation algorithm of Fig. 1 has computed a wave-function |ψ n . For simplicity, we further suppose that precisely one measurement operator pair {M (+) , M (−) } acts during the subsequent time-step. These operators satisfy the normalization condition where I is the identity operator. In the absence of GK projection the post-timestep density matrix ρ n+1 is readily shown to be where ρ n ≡ |ψ n ψ n |. This expression implicitly defines the well-known linear superoperator L to be that linear operation on (Hermitian) matrices that takes ρ n → ρ n+1 . The existence and strict positivity of this linear map is one of the main defining characteristics of linear quantum mechanics (as reviewed in Section 1.3.4). Now we ask "What mathematical operations upon {M (+) , M (−) } leave L invariant?" The Theorema Dilectum in its algebraic form due to Choi [48] states that all such invariance operations are of the general form where U is an arbitrary 2 × 2 unitary matrix of complex numbers (i.e. a matrix acting on the linear space of measurement operators, not the Hilbert space of |ψ n , such that the matrix elements of U are c-numbers). This is the sole general mathematical invariance of the first two steps of the simulation algorithm of Fig. 1, and so (from the QMOR point of view) it is the most fundamental invariance of linear quantum mechanics. In Section 3.3.6, we will establish that U -transform invariance enforces physical causality. We are thereby motivated to ask "How is the U -transform invariance of the Theorema Dilectum affected by GK projection?" According to the algorithm of Fig. 1, GK projection modifies (46) to where we recall that (P K ) n ≡ P K (|ψ n ) was defined in (7) as the operator that projects onto the local tangent space of the GK manifold at |ψ n . It is easy to check that tr[ρ n+1 ] = 1 (i.e., probability is conserved), and that above expression reduces to the linear result (46) whenever the commutators [(P K ) n , M (+) ] and [(P K ) n , M (−) ] both vanish.
and so we can regard all operator products involving δM (+) and δM (−) to be O(1).
In aggregate, our definitions and normalizations ensure that ρ n+1 as given by (48) has a well-defined power series expansion in ; it is therefore a straightforward (but not short) calculation to verify that this expansion can be written in the form: whereP K = I − P K . We note that the leading terms are determined solely by L and P K , and therefore are invariant under the U -transform (47) of the Theorema Dilectum. Furthermore, it can be shown that the value of the small parameter is itself invariant under the U -transform, and so is the sum |c| 2 + |s| 2 . And finally, we have calculated the (exact) c-and s-dependence of the O( 3 ) terms. We will establish in Section 3.3 that in the continuum limit of infinitesimally small time step intervals δt, physical quantities (for example, relaxation rates) are O( 2 /δt). Physically this means that the O( 3 ) terms in (54) are negligible in the continuum limit, provided the technical conditions |c| > 0 and |s| > 0 are satisfied (these technical conditions provide the rationale for calculating the c-and s-dependence of the O( 3 ) terms in (54) ). These results motivate us to adopt from Carlton Caves [46] the following definition: The result (54) then expresses the following fundamental design rule of QMOR mechanics: Design Rule 3.1. In the continuum limit, quantum trajectory simulations of first-class measurement processes, as projected onto state-spaces of gabion-Kähler type (QMOR simulations for short), respect the unitary invariance of the Theorema Dilectum.
In physical terms, first-class measurements are characterized in the continuum limit by stochastic drift and diffusion processes on GK manifolds, rather than quantum jumps. To ensure that the Theorema Dilectum is respected, QMOR mechanics must therefore simulate all quantum systems wholly in terms of drift and diffusion processes upon strongly-curved, low-dimension state-space manifolds of gabion-Kähler type.
This low-dimension curved-geometry stochastic description of QMOR mechanics thus contrasts sharply with the high-dimension linear-geometry "jump-oriented" description of quantum mechanics that is commonly given in textbooks; that is why many of our derivation methods and design rules are novel.
Turning our attention briefly to the broader context of geometric quantum mechanics (see Sections 1.3.5 and 2.13), we emphasize that we do not regard issues of causality in geometric quantum mechanics as settled by the projective generalization of the Theorema Dilectum given in (54). We will not discuss this topic further, partly for reasons of space, but mainly because we regard it as being an exceptionally difficult and subtle topic that is intimately bound-up with the question of the existence (or nonexistence) of quantum field theory in geometric quantum mechanics.

QMOR respects thermal equilibrium
Our recipes for simulating contact with thermal reservoirs in QMOR mechanics yield an algebraic result (not previously known) that holds exactly even in linear quantum mechanics, and that can be verified without reference to drift and diffusion equations. We state and prove this result now, so that it can provide a well-defined mathematical target for our subsequent QMOR analysis.
We start by identifying a spin=j state having z-axis quantum number m with the ketvector |j, m . Then a coherent state |x associated with a unit-vector spin directionx is by definition |x = D(φ, θ, 0)|j, j , where the rotation operator D that carriest = (0, 0, 1) intô The rotation operators are well understood. In particular, an identity due to Wigner [93,169] gives j, m|x in closed form as It follows that x|x = 1 and x|s|x = jx. It is well known [85,153,154,162] (and not hard to show from (56) ) that a resolution of the identity operator I is The Q-representation and P -representation of a Hermitian operator ρ are then defined (following Perelomov's conventions [153]) as Given an arbitrary Hermitian operator ρ, it is known that in general a P -representation P (x|ρ) can always be constructed. In brief, the construction is as follows: from the ansatz P (x|ρ) = ∞ l=0 j m=−l a l,m Y l m (θ, φ), with Y l m (θ, φ) a spherical harmonic, a set of linear equations for the coefficients a l,m is obtained by substituting (59) into (58) and expanding both sides in spherical harmonics. Solving the resulting linear equations always yields a valid P -representation. However, such P -representation constructions are non-unique in consequence of the identity [153,154] which can be proved directly from (56) as a consequence of the addition law for angular momenta. This result shows explicitly that the set of coherent states |x is over-complete (as is well-known).
With the above as background, we now slow the pace of presentation. We consider the problem of finding a positive P -representation for a given operator ρ, that is to say, a representation for which P (x|ρ) ≥ 0 for allx.
Positive P -representations have the useful property (for simulation purposes) of allowing us to interpret P (x|ρ) as a probability distribution over spin directionsx. But from a mathematical point of view, distressingly little is known about positive P -representations, in the sense that there is no known general method for constructing them, or even for determining whether they exist in a given case.
We will now consider an operator that often appears in practical QMOR simulations: the thermal operator ρ th j defined by where {s 1 , s 2 , s 3 } are the usual spin-j operators satisfying [s 1 , s 2 ] = is 3 (and cyclic permutations), andt is a unit axis along which a spin is thermally polarized with inverse-temperature β. By inserting a complete set of states into (58) and then substituting Wigner's expression (56), a well-known [153,162] closed-form expression for the Q-representation of ρ th j can be obtained: We now exhibit a closed-form analytic expression for a positive P -representation of the thermal operator as a distribution over coherent states, which is exact for all j and β: Design Rule 3.2. A positive P -representation for the spin-j thermal operator ρ th j is given in terms of the Q-representation by To our knowledge this is the first such P -representation given, other than the j → ∞ limit of quantum optics in which P and Q are both simple Gaussians. If we regard the above result solely as a mathematical theorem to be proved by the most expedient means, we can do so by treating it as an ansatz. The resulting proof is short. Taking matrix elements of the defining relation (59) between states j, m| and |j, m , and without loss of generality settingt = (0, 0, 1), Design Rule 3.2 is equivalent to the definite integral e −βm δ mm = 2j + 1 4π Substituting the Wigner representation (56) and the Q-representation (62) into (63), we can check by numerical integration that (63) is correct for randomly chosen values of j, m, m , and β. Thereby encouraged, we soon discover an integration strategy by which (63) yields to analytic evaluation in the general case. In brief, the φ-angle integration yields the requisite δ mm factor; the θ-angle integration can be transformed into an integral over rational functions in z = cos θ via identities like cos 2j 1 2 θ = (1 + z) j /2 j ; and the resulting integral is recognizably a representation of the Gauss hypergeometric series [1, eqs. 15 (63).
Design Rule 3.2 is simultaneously frustrating, reassuring, and intriguing. It is frustrating because our QMOR analysis will answer the natural question "Where did the ansatz come from?" by "It is the solution to a Fokker-Planck equation that describes spin-systems in thermal equilibrium." But this answer provides no satisfying rationale for why the P -representation exists, or for its simple analytic form. Furthermore, QMOR analysis will provide no answer at all to the natural follow-on question, "Is there a reason why the thermal operator's positive P -representation is simply given in terms of its Q-representation?" Design Rule 3.2 is reassuring because it tells us that QMOR simulations can simulate thermal equilibrium (at least in simple cases). From a physical point of view this reassures us that the dimensional reduction associated with QMOR mechanics preserves at least some crucial thermodynamic physics. From a practical point of view it provides a reassuring consistency check that the (rather lengthy) chain of theorems and stochastic analysis that leads to the P -representation is free of algebraic errors.
And finally, Design Rule 3.2 is intriguing because it suggests that a partial answer to the open question "Which operators have positive P -representations?" might be "Only thoses operators that are proportional to a density matrix that is associated with a stationary measurement-and-control process that drives input states to coherent states." This conjecture is intriguing because its negation is logically equivalent to "Density matrices exist that cannot be the result of any stationary measurement-and-control process that drives input states to coherent states." Such density matrices would be remarkable entities from both a mathematical and physical point of view, and in particular their existence (or nonexistence) is directly relevant to the further development of practical design rules for QMOR simulations. The conjectured answer is therefore intriguing whether it is true or false.
3.2.3. QMOR respects classical linear control theory QMOR simulations of macroscopic objects (like MRFM cantilevers) regard them as spin-j quantum objects having very large j. We will see that the resulting dynamics typically are linear. Engineers have their own idioms for describing linear dynamical systems, which are summarized in block diagrams like this one: (64) To show that QMOR simulations can accurately model systems like the above, we will transform the above diagram-using strictly classical methods-into an operationally equivalent form that more naturally maps onto QMOR algorithms and geometric structures.
For convenience, these classical equivalences are summarized as design rules in Figs. 6-7. Experienced researchers will recognize these equivalences as being elementary, but to the best of our knowledge, they have not previously been recognized in the literature of engineering or physics.
As with the Feynman diagrams of physics, the block diagrams of engineering depict systems of equations. Our diagram conventions are standard, as briefly follows. The block diagram (64) corresponds to a set of linear relations between force noise f n (t), measurement noise q n (t), input external force f ext (t), and output measurement q m (t), which we take to be classical real-valued functions. In particular we specify that f n (t) and q n (t) are white noise processes having correlation functions C satisfying so that S q n and S f n are (one-sided) white-noise spectral densities. A circle is a node whose inputs are added and subtracted, a cross is a node whose inputs are multiplied, a triangle indicates a positive real scalar gain γ, and a square box indicates convolution  with a general real-valued stationary kernel K such that Alternatively, convolution blocks can be specified in the Fourier domain. Our Fourier transform convention is thatã(ω) is defined to bẽ and similarly forb(ω),K(ω), etc. Therefore a frequency-domain description of (66) is We build our QMOR block diagrams from three classes of linear classical kernels: dynamical kernels ( ), feedback kernels ( ), and backaction kernels ( ), whose defining mathematical properties are specified in Fig. 6.
In brief, dynamical kernels by definition are time-reversal invariant, feedback kernels by definition are causal, and the backaction kernel is the Hilbert transform that is well-known (and much-used) by signal processing engineers. We remark that the Hilbert transform is formally non-causal, but in practical narrow-bandwidth applications (like radio transmitters, acoustic processors, MRFM cantilever controllers, etc.) its effects can be closely approximated by a causal derivative transform.
For causal kernels, analytic continuation from the Fourier variable ω to the Laplace variable s = iω is well-defined. Partly because causal kernels are of central importance in control engineering, and partly by tradition, Laplace variables are more commonly adopted in the engineering literature than Fourier variables (although both are used). However, the Laplace analytic continuation of the non-causal Hilbert transform kernelH(ω) = i sgn(ω) is not well-defined. For this reason our analysis will focus exclusively upon upon time-domain and frequency-domain (Fourier) kernel representations.
Central to QMOR analysis and simulation is the purely mathematical fact (which also appears in Fig. 7 as Design Rule 3.4) that the following classical systems are operationally equivalent: By "operationally equivalent" we mean that the dynamical relation between the applied forces f ext (t) and the measured position q m (t) is identical for the two sets of equations, as are the stochastic properties of q m (t). The internal state of the system is of course very different for the two cases, but so long as we adhere to the strict operational principle "never observe the internal state of a system (even a classical system)" this difference is immaterial.
Developing design rules for QMOR mechanics is considerably easier if we habituate ourselves to the Hilbert transform that appears in all three classical design rules of Fig. 7. We begin by remarking that from an abstract point of view, the Hilbert transform defines a complex structure on the space R of real-valued functions r, that is, a map H : R → R that satisfies H 2 = −I. This corresponds to two diagrammatic identities = and (70) that make it obvious (since the right-hand sides are causal) that the non-causal aspects of vanish when it is applied an even number of times. Because the even-numbered moments wholly determine the statistical properties of (zero-mean) Gaussian noise, we begin to see how it is that noncausal Hilbert transforms can appear in the equivalences of our classical design rules without inducing observable causality violations. In the broader context of quantum simulations, causality is assured by the Theorema Dilectum, and we will establish in a later section that the Theorema Dilectum is directly responsible for the appearance of Hilbert transforms in the classical limit.
We saw already (in Section 1.5.3) that complex structures are a defining geometric characteristic of Kähler manifolds, so the appearance of Hilbert transforms in real-valued classical dynamics is a mathematical hint that noisy classical dynamical trajectories support a complex structure that projects naturally onto gabion-Kähler manifolds. This complex structure is manifest not in the (real-valued) state variables of classical systems, but in the causal properties of the response of classical systems to noise.

Remarks upon the spooky mysteries of classical physics
For teaching purposes, it is helpful (and amusing) to pretend that we live in a world in which linear control theory is taught according to a nonstandard ontology in which the Hilbert transform has a central role. This ontology was conceived as a philosophical provocation (a herausforderung [175]), but it has subsequently proved to be useful in teaching and a fertile source of technical inspiration. We will call it the classical Hilbert ontology, or sometimes just the Hilbert ontology. Our motivation for emphasizing the mysterious properties of classical measurement theory is similar to the motivation of Wheeler and Feynman [196,197] in proposing their non-causal classical electrodynamic ontology proposed in the 1940s.
The tenets of Hilbert ontology are taken to be: Figure 7: Three design rules that reflect "spooky mysteries of classical physics." The above design rules follow from elementary classical considerations as (briefly) follows. Design Rule 3.3 follows from the time-reversal invariance of the dynamical kernel G 0 (see Fig. 6) plus the statistical properties of Gaussian noise. Design Rule 3.4 is then obtained by adjoining a causal feedback kernel H to the diagram of Rule 3.3. Design Rule 3.5 then follows as a practical application of Rule 3.4, in which a causal kernel H approximates a non-causal Hilbert backaction kernel H within a device's finite operational bandwidth.
• By Occam's Razor, ontologies that invoke one noise source are preferred over ontologies that invoke two. Therefore the Hilbert ontology regards backaction as being a physically real phenomenon that is universally present in all noisy dynamical systems (its mathematical expression being Design Rule 3.3 of Fig. 7).
• It follows (by Design Rule 3.4) that measurement noise always back-acts upon system dynamics in such a way as to "drag" the state of the system into agreement with the measurement. This state-dragging Hilbert backaction has a central ontological role: it is the fundamental mechanism of nature that causes measurement processes to agree with reality.
• As a measurement process approaches the zero-noise limit, the increasingly strong state-dragging Hilbert backaction from that measurement process dynamically "collapses" the variable being measured, so as to force exact agreement with the measurement.
• As a corollary, zero-noise measurements are unphysical, because they imply infinitely strong backaction. That is the explanation in Hilbert ontology of why all real-world measurement processes are noisy.
• In narrow-band systems, it is possible to cancel the Hilbert backaction noise via causal feedback control. This means that zero-temperature narrow-band systems can be simulated, or even realized in practice, provided that all noise processes are accessible for purposes of feedback (as mathematically expressed by Design Rule 3.5 of Fig. 7).
• Although the mathematical kernel associated with Hilbert feedback is non-causal, this noncausality cannot be exploited for purposes of communication, for the physical reason that the backaction kernel transmits only noise. The mathematical proof is simply that Design Rules 3.3-5 can be expressed as equivalent systems having all kernels causal (as shown in Fig. 7). But in Hilbert ontology this causality proof is viewed as being a purely mathematical artifice, because it postulates "spooky" uncorrelated forces that "obviously" have no physical reality (according to our Occam's Razor tenet) since no mechanism is given for them.
For purposes of this article, we designate the above Hilbert ontology to be the "true" classical reality of the world, and we seek to provide a microscopic justification of it from orthodox quantum mechanics (recognizing of course that the Hilbert ontology has been constructed specifically to ensure that our overall herausforderung yields practical and interesting results). From a teaching point of view, the Hilbert ontology helps students appreciate that the mysteries and ambiguities that traditionally are taught as belonging exclusively to quantum mechanics-like "wave functions collapsing to agree with measurement" and "noncausal correlations"-are manifest too in (at least one) wholly classical ontology.
It is traditional in both the popular and the scientific literature to call these ontological mysteries and ambiguities "spooky." The popularity of "spooky" in the scientific literature can be traced to an influential article by Mermin [138], who adopted it from an idiomatic phrase "spukhafte Fernwirkungen" that Einstein uses in the Einstein-Bohr correspondence; this phrase is generally translated as "spooky action at a distance" [99,105,141]. Given the central importance of the spookiness of quantum mechanics, it would be astonishing if this spookiness was wholly invisible at the classical level.
Our Hilbert ontology therefore serves both to guide our calculations and to help us celebrate the "spooky mysteries of classical physics."

Experimental protocols for measuring the Hilbert parameters
A crucial test of an ontology is whether it motivates us to proceed to practical calculations that yield useful and/or surprising results; we now do so. We consider a laboratory course in which students are guided by the Hilbert ontology to explore the fundamental and practical limits of lownoise sensing and amplification. For definiteness, we assume that the students work with nanomechanical oscillators as force detectors (as in MRFM technology), but a similar course could feasibly be organized around radio-frequency (RF) sensing and amplification, optical sensing and amplification, or acoustic sensing and amplification.
We consider the experimental problem of measuring the two Hilbert parameters of a nanomechanical oscillator, namely the measurement noise S q and the Hilbert gain γ. We neglect the intrinsic damping of the oscillator, such that the dynamical kernel isG 0 (ω) = 1/[m(ω 2 0 − ω 2 )], where m is the mass of the oscillator and ω 0 is the resonant frequency. The measurement is readily accomplished by the following protocol. Derivative feedback control is applied having a kernelH(ω) = iΓω, and the value of the controller gain Γ is adjusted until the spectral density S q m of the measured cantilever displacement q m (t) is observed to be flat in the neighborhood of the cantilever frequency ω 0 . The control gain adjustment is straightforward: if the spectrum of q m (t) has a peak, then Γ is too small; if the spectrum has a hole, then Γ is too large.
The preceding protocol is perfectly feasible in practice (i.e., it is not a Gedankenexperiment). The protocol fails only when the required dashpot gain is impractically large, such that Γω 0 k, where k = mω 2 0 is the spring constant of the cantilever. Such failures typically indicate that a measurement process has a noise level that is too large to be of practical interest.
In practical cases the required control gain typically satisfies Γ 0 Γ k, where Γ 0 = k/(ω 0 Q) is the intrinsic damping of a cantilever having quality Q. In such cases the intrinsic damping Γ 0 has negligible effect on the time-scale of the controlled response, and so we can regard the cantilever's dynamical kernel as having the time-reversal-invariant dynamical form G 0 described in Fig. 6. The design rules of our Hilbert classical ontology therefore can be applied without modification. Specifically, Design Rule 3.5 determines the Hilbert backaction to be γ = Γω 0 , and determines the measurement noise to be S q n = S q m (ω 0 ).
Of course, Design Rule 3.5 mathematically assures us that the adherents of traditional classical ontology can-without operational contradiction-ascribe these same measurements to a fictitious force noise f n (t) having spectral density S f n = γ 2 S q n , so that the vexing question of whether this force noise is "fictitious" or "real" is operationally immaterial.

QMOR simulations respect the fundamental quantum limits
We now apply the results of the preceding section to establish criteria that QMOR simulations must satisfy in order to respect the known fundamental quantum limits to amplifier noise and force measurement.
We suppose that a force signal f (t) = f 0 cos(ω 0 t + φ 0 ) is applied to the oscillator. The carrier frequency ω 0 is tuned to match the oscillator resonance frequency, and the unknown magnitude f 0 and unknown phase φ 0 of the signal are to be determined from measurement. This is a common task in practice.
When the oscillator is configured according to the measurement protocol of the preceding section, then according to the right-hand block diagram of Design Rule 3.5 of Fig. 7, the mean power p 0 absorbed by the oscillator's (wholly classical) feedback controller during the measurement process is p 0 = f 2 0 /(2γ). The absorbed power inferred from the (wholly classical) measurement record q m (t) has an equivalent (one-sided) noise PSD S output p 0 whose expression in terms of the Hilbert parameters is readily shown to be S output Now we connect this result with a known fundamental quantum limit on noise in power amplifiers. We adopt the definition of Caves [45]: An amplifier is any device that takes an input signal, carried by a collection of bosonic modes, and processes the output to produce an output signal, also carried by a (possibly different) collection of bosonic modes. A linear amplifier is an amplifier whose output signal is linearly related to its input signal.
Following a line of reasoning put forth in the 1960s by Heffner [106] and by Haus [104], which Caves article develops in detail [45], we suppose that the input power is supplied by a bosonic mode-type device (like a resonant circuit, or an RF wave-guide, or a single-mode optical fiber) whose power level has been independently measured with a shot-noise-limited photon counting device. According to orthodox quantum mechanics such counting processes have Poisson statistics (we will establish in Section 3.4 that QMOR simulations respect this rule), and therefore the input power has a quantum-limited noise PSD given by The measured power and phase suffice to create an arbitrary-gain replica of the input signal, and the noise-figure (NF) of this effectively infinite-gain power amplifier is simply given from (71) and (72) by In Caves' nomenclature, we are regarding our continuously measured oscillator as an equivalent "phase-insensitive linear amplifier" having infinite gain. The analyses of Heffner, Haus, and Caves [45,104,106] establish that what Caves calls the "fundamental theorem for phase-insensitive power amplifiers" is simply NF ≥ 2 (in the infinite-gain limit), which in decibels is 10 log 10 2 3 dB. We remark that the 3 dB quantum limit to power amplifier noise has been experimentally observed [108] and theoretically analyzed [176] since the early days of maser amplifiers in the 1950s; our review focuses upon the work of Heffner, Haus, and Caves solely because their analyses are notably rigorous, general, clearly stated, and (importantly) their predicted quantum limits are mutually consistent and consonant with subsequent experiments.
Our result (73) then establishes that the Heffner-Haus-Caves noise-figure limit finds its expression in QMOR analysis in the following three equivalent ways: The middle expression we recognize as the continuous-measurement version of the standard quantum limit to force and position measurement, in precisely the quantitative form derived by Braginsky and Khalili [22], which in turn derives from earlier seminal work by Braginsky, Vorontsov, and Thorne [23]. Although the Braginsky-Khalili limit was derived by very different methods from the Heffner-Haus-Caves' limit, we see that the two quantum limits are equivalent. To the best of our knowledge, this equivalence has not previously been stated in the above quantitative form. What we have chosen to call "the Hilbert limit" on the right-hand side of (74) has (to our knowledge) not previously been recognized anywhere in the literature, and yet from the viewpoint of Hilbert ontology it is the most fundamental of the three. We express the above three-fold equivalence as a fundamental QMOR design rule: Design Rule 3.6. QMOR simulations respect the quantum measurement limit in all its equivalent forms: the noise-figure (NF) limit in power amplifiers (NF ≥ 2), the standard quantum limit to the measurement of canonically conjugate variables (S f n S q n ≥ 2 ), and the Hilbert limit that measurement noise and state-dragging Hilbert backaction cannot both be small (γS q n ≥ ).
We present a quantum derivation of these limits in Section 3.4, in the context of a more general analysis that encompasses nonlinear quantum systems.
3.2.7. Teaching the ontological ambiguity of classical measurement What should we teach students about the internal state of the system during the preceding (strictly classical) protocols for measuring the Hilbert parameters S q n and γ? From a strictly logical point of view, this question need not be answered, since our measurement protocols and design rules are careful to make no reference to the internal state (even at the classical level). But in practice, some answer must be given, for the pragmatic reason that students require some coherent story about what the systems they are measuring are "really" doing. This coherence is especially necessary to those science and engineering students (a substantial portion, in our view) whose careers will require that they extend their professional expertise from the classical to the quantum domain.
It is therefore desirable that common-sense questions from students receive commonsense answers from teachers, and it is equally desirable that these answers prepare for a classical-to-quantum educational transition that is as nearly seamless as is practicable.
Adherents of traditional classical ontology can argue cogently as follows: "It may be operationally correct to ascribe experimental results wholly to measurement noise and Hilbert backaction, but it is physically wrong, because we have strong physical reasons to believe that the excitation of the oscillator is being driven by force noise from a thermal reservoir whose internal dynamics we do not observe." Adherents of Hilbert classical ontology can offer an similarly cogent counter-argument, which however requires the assertion of a definite mathematical result (in italics): "We too believe that the observed excitation of the oscillator is in fact being driven an unobserved thermal reservoir, but the action of a unobserved thermal reservoir can be ascribed to a covert process of measurement, Hilbert backaction, and control." The practical consequence of the above reasoning is that this article's embrace of Hilbert ontology has practical utility only if we can develop well-posed QMOR algorithms by which the action of thermal reservoirs upon a dynamical system is simulated by unobserved/covert processes of quantum measurement, Hilbert backaction, and classical control. We develop the necessary algorithms in Section 3.4, using as our main mathematical tool the Theorema Dilectum that was already given as Design Rule 3.1. A key mathematical result will be the positive P -representation that was already given as Design Rule 3.2; this helps us foresee the analysis path by which the QMOR analysis program will succeed. And of course Design Rules 3.3-5 will emerge naturally in the course of our analysis.
In summary, for students and teachers alike, a sufficient justification for embracing the Hilbert ontology is that it leads to QMOR simulation algorithms that are computationally efficient, operationally orthodox, and mathematically novel. A provocative side-effect is that we learn to perceive the "spooky mysteries of quantum physics" as being manifest in all noisy dynamical systems, even classical ones.

Physical aspects of QMOR
We now turn our attention to the physical aspects of measurement, our goal being to establish connections between the concrete description of measurement in terms of hardware and experimental protocols on the one hand, and the measurement operator formalism of our QMOR algorithms on the other hand.

Measurement modeled as scattering
We adopt as the fundamental building block of our simulations a single particle of spin j, described by a wave function |ψ , and numerically encoded as a complex vector with dim |ψ = 2j + 1. We will simulate all noise and all measurement processes by scattering photons off the spin, one photon at a time, and we associate each scattering event with a single time step in Fig. 1.
We envision photon scattering as the sole mechanism by which noise is injected into our simulations, and interferometry as the sole means of measurement. We describe photon scattering as a unitary transformation acting on the spin state (before scattering) |ψ n → exp(i2θs op ) |ψ n (after scattering).
A purely conventional factor of 2 is inserted in the above to simplify our calibration rules. We will call s op a measurement generator. In general s op can be any Hermitian matrix, but in our simulations it will suffice to confine our attention to s op ∈ {s 1 , s 2 , s 3 }, where {s 1 , s 2 , s 3 } are rotation matrices satisfying the commutation relation [s j , s k ] = i jkl s l . These matrices generate the rotation group, and our discussion will assume a basic knowledge of their algebraic properties, which are discussed at length in many textbooks (e.g., [93,139,169,198]).
We adopt the near-universal convention of working exclusively with spin operator matrix representations that are irreducible and sparse, having dimension 2j + 1 for j the spin quantum number, with s 3 a real diagonal matrix, s 1 real and bidiagonal, and s 2 imaginary and bidiagonal. The scattering strength is set by the real number θ, which in our simulations will always satisfy θ 1. Figure 8 provides three idealized diagrams for thinking about the physical, geometric, and algebraic aspects of photon interferometry, similar to the idealized diagrams in the Feynman Lectures [69, chs. 5-6] that depict the Stern-Gerlach effect. Figure 8(a) depicts physical fiber-optic interferometers that confine photons within lowloss optical fibers [26,171]. The region of overlap between the fibers allows photons to transfer from one fiber to another with an amplitude that is subject to engineering control. This overlap region plays the role of the semi-silvered mirror in a traditional Michelson interferometer. Figure 8(b) depicts optical couplers geometrically, as general-purpose devices for linking incoming and outgoing optical amplitudes. Couplers thus can be braided into optical networks of essentially arbitrary topology. Figure 8(c) depicts optical couplers algebraically, as linear maps between incoming complex optical amplitudes a in = {a in tl , a in tr , a in bl , a in br } and outgoing optical amplitudes a out = {a out tl , a out tr , a out bl , a out br } such that a out = S · a in , with S the optical scattering matrix. With regard to Fig. 8(c), the index "tl" is the "top left" port, "br" is the "bottom right" port, etc.. Our normalization convention is such that the probability of detection of an outgoing photon is |a out | 2 . Optical losses in real-world couplers are small, such that |a out | 2 = |a in | 2 , i.e. the optical scattering matrix S is unitary.

Survey of interferometric measurement methods
The stochastic measurement and noise processes in our simulations can be conceived as interferometric measurements performed on each scattered photon, and this physical picture will prove very useful in designing MOR techniques. Such measurements require that each incoming photon be interferometrically split before it scatters from the spin, to allow subsequent interferometric recombination and measurement. It is convenient to conceive of this initial splitting as performed by a 2×2 single-mode optical fiber coupler, as illustrated in Figs. 9(a-c). The devices of this figure may be regarded as physical embodiments of the simulation algorithm of Fig. 1.
This physical picture embodies the idealizing assumptions that optical couplers are exactly unitary, that photon emission into the fiber takes place at equally spaced intervals δt, and that photon detectors register a single classical "click" upon detection of each photon, which occurs with detection probability |a out | 2 . The unitarity of photon scattering and interferometric propagation then ensures that each incoming photon results in precisely one detection click.
With respect to the algorithm of Fig. 1, the stochastic selection of operator M k (+) versus M k (−) is physically identified with these clicks, such that the sole data set resulting from a simulation is a set of classical binary data streams, with each stream comprising the recorded clicks for a measurement operator pair. Such binary streams closely resemble realworld MRFM experimental records, in which signals corresponding to photoelectrons from an optically-monitored cantilever are low-pass filtered and recorded.
We remark that the interferometers of Fig. 9(a-c) can be regarded with equal validity either as idealized abstractions or as schematic descriptions of real-world experiments. For example, in our simulations we will regard MRFM cantilevers as spins of large quantum number j, in which case the interferometer geometry of Fig. 9(c) is identical (in the topology of its light path) to the real-world fiber-optic interferometers used in typical MRFM experiments [26,170,171].
Trapped-ion experiments are examples of continuously observed small-j quantum systems, since the quantum mechanics of a two-state atom is identical to the quantum mechanics of a spin-1 2 particle. Such experiments presently are conducted with photons detected directly as in Fig. 9(b) [18,130,158]. We remark that it would be theoretically interesting, and experimentally feasible, to conduct trapped ion experiments with weakly interacting photons detected interferometrically, as in Fig. 9(c).
In trapped ion experiments, the observed transitions in the photon detection rates are observed to have random telegraph statistics, which are typically attributed to "quantum jumps" [18,130] or to "instantaneous transitions between energy levels" [158]. As was recognized by Kraus more than twenty years ago [123, p. 98], "such reasoning is unfounded," and we will see explicitly that observation of telegraph statistics does not imply discontinuous evolution of |ψ . Thus, readers having a classical MOR background need not regard discussions of quantum jumps in the physics literature as being literally true. This physical embodiment of the formal QMOR simulation algorithm of Fig. 1, using 2 × 2 fiber-optic couplers to scatter photons off a spin state. In (a) the photons are detected with equal probability, such that no information about the state is obtained. In (b) the phase is detected via homodyne (self-interfering) interferometry, such that a small amount of information about about the state is obtained. To an outside observer, the fate of the downstream photons is immaterial, hence (a) and (b) must embody physically indistinguishable noise processes, despite their differing quantum description. This is the physical content of the Theorema Dilectum. From a mathematical point of view, the free choice of an arbitrary downstream unitary transform upon the photon paths is manifested in the U -transform invariance of (47). (On-line version only: photons in red are causally downstream of the scattering, and thus can be measured in arbitrary fashion without altering the physically observable properties of the noise being simulated.)

Physical calibration of scattering amplitudes
The calibration of our simulations will rely upon a physical principle that is well-established, but somewhat counterintuitive. The principle is this: measurements of the scattering phase (75) suffice to provide detailed information about the Hamiltonian that is responsible for the scattering. Gottfried's discussion of Fredholm theory in atomic scattering [93,Section 49.2] provides a good introduction to this topic.
As a calibration example, we consider here the measurement process associated with the simple interferometer of Fig. 9(a). Using the S-matrix of Fig. 8(c) to propagate the input photon in Fig. 9(a) The above equation follows the convention that the optical amplitudes of the top (bottom) fiber path in Fig. 9(a) are listed as the top (bottom) element of the above column arrays, such that successive interferometric couplings are described by successive 2 × 2 unitary matrix operations. This convention is common in optical engineering because it unites the geometric and algebraic descriptions of Fig. 8(b-c). In the context of our quantum simulations, each element of the above arrays is formally an operator on the wave function |ψ , but for c-number (complex number) array elements like "0", "1" and "i" an implicit identity operator is omitted for compactness of notation. These identity operators physically correspond to events that are not dynamically coupled to the spin state, like photon emission, propagation through interferometers, and subsequent detection. The resulting overall measurement operators {M

Noise-induced Stark shifts and renormalization
We now consider an optical scattering effect known as the AC Stark shift, as induced by the scattering process of Fig. 9(a). Applying the simulation algorithm, we find that during a time ∆t δt the final state of the simulation accumulates a state-dependent phase shift such that where n (+) is the number of photons detected on the (+)-channel. It is easy to show that n (+) has mean µ (+) = ∆t/(2δt) and standard deviation σ (+) = µ (+) /2. We see that the average effect of the photon scattering is equivalent to a dynamical Hamiltonian H op , such that which we identify as the effective Hamiltonian of a Stark shift. The Stark shift fluctuates due to statistical fluctuations in the number of photons detected on the (+)-channel, such that in the continuum limit the photon detection rate r(t) is a stochastic process having mean µ r = 1/(2δt) and white-noise spectral density S r = 1/(4δt) (in a two-sided spectral density convention). This implies that the Stark shift fluctuations have an operator-valued power spectral density (PSD) This result calibrates the externally-observable Stark shift parameters {H op , S H op } in terms of the internal simulation parameters {θs op , δt} and vice versa.
The preceding two equations (78)(79) reflect the well-known phenomenon in physics that interaction with a measurement process or thermal reservoir renormalizes the physical properties of a system. But as presented here, these same equations exhibit a known pathology in the limit δt → 0: if we take θ ∝ √ δt such that the Stark noise (79) has a finite limit, then the magnitude of the Stark shift Hamiltonian (78) diverges. The origin of this (unphysical) infinite energy shift is our (equally unphysical) assumption that measurement "clicks" occur at infinitely short intervals, such that the PSD of the measurement noise is white; such white noise processes inherently are associated with infinite energy densities.
In the present article we will simply repair this white-noise pathology by simply adding a counter-term to the measurement process, such that the Stark shift is zero in all our measurement processes. In the language of renormalization theory, we redefine all of our measurement operators so that they refer to the "dressed" states of the system. See (e.g.) [180] and [128] for further discussion of this delicate point, whose detailed analysis is beyond the scope of the present article.
3.3.6. Causality and the Theorema Dilectum In Section 3.2.1 we discussed the Theorema Dilectum from a mathematical point of view that defined it terms of the (unitary) Utransformation of (47). In Fig. 9(a-c) we associate the U -transformation with the physical choice of what to do "downstream" with photons that have scattered off the spin. If we regard this downstream choice as being made made by Alice, and we assume that Bob is independently monitoring the same spin as Alice, then the physical content of the Theorema Dilectum is that Alice's downstream measurement choices can have no observable consequences for Bob, and in particular, Alice cannot establish a communication channel with Bob via her measurement choices.
We therefore physically identify U in (47) with the downstream optical couplers in Fig. 9(a-c), such that the algebraic freedom to specify an arbitrary unitary transform U is identified with the physical freedom to "tune" the interferometer by adjusting its coupling ratio and fiber lengths (which adjust the phase of the output amplitudes relative to the input amplitudes). Because experimentalists are familiar with this process and with the phrase "interferometer tuning" to describe it, we adopt the word "tuning" to describe the process of adapting U to optimize {M (+) , M (−) } for our simulation purposes.

The
Theorema Dilectum in the literature The invariance associated with U has received most of its attention from the physics community only recently; it is not mentioned in most older quantum physics texts. Both Preskill's class notes [159] and Nielsen and Chuang's text [146] give physics-oriented proofs of the Theorema Dilectum. The crux of all such proofs is that the choice of U does not alter the density matrix associated with the ensemble average of all possible trajectories, such that the choice has no observable consequences.
Carmichael [40, p. 122] seems to have been among the first to use the now-widely-used term unraveling in describing quantum trajectory simulations, and he explicitly recognized that this unraveling is not unique: We will refer to the quantum trajectories as an unravelling of the source dynamics since it is an unraveling of the many tangled paths that the master equation evolves forward in time as a single package. It is clear [. . . ] that the unraveling is not unique.
(emphasis in the original). Diverse points of view regarding the ambiguity of trajectory unraveling can be found in the physics literature. Rigo and Gisin [166] argue that it is central to our understanding of the emergence of the classical world, and they make their case by presenting four different unravelings of a single physical process; we will adopt a similar multiple-unraveling approach in analyzing the IBM single-spin experiment (see Section 4.2). Percival's text adopts the equally valid but sharply contrasting view that [152, p. 46]: "the ambiguity [. . . ] is a nuisance, so it is helpful to adopt a convention which reduces this choice." Breuer and Petruccione [25] simply state a result equivalent to the Theorema Dilectum, without further comment or attribution.
Preskill's course notes and Nielsen and Chuang's text both take a middle point of view. They briefly describe the Theorema Dilectum as a "surprising ambiguity" [159, sec. 3.3] that is "surprisingly useful" [146, sec. 8.2.4]. The word "surprising" invites readers to think for themselves about unravelling, and the word "ambiguity" suggests (correctly) that the implications of this invariance are not fully understood. Nielsen and Chuang's text notes that up to the present time, the main practical application of the Theorema Dilectum has come in the theory of quantum error correction, where the freedom to choose unravelings that facilitate the design of error correction algorithms has been "crucial to a good understanding of quantum error correction" [146, sec. 8.2.4].
In the context of quantum computing theory, Buhrman et al. [27] have exploited the informatic invariance of the Theorema Dilectum to show that certain logical gates that are essential to universal quantum computing, when made noisy, can be indistinguishably replaced by randomly-selected gates from a restricted set of gates that can be simulated classically. Their proofs build upon earlier work by many authors [3,24,101] and in particular upon the Gottesman-Knill theorem [92]. To our knowledge, this is the first formal quantum informatic proof that additional noise makes systems easier to simulate.

Designs for spinometers
We now present some basic designs for measurement operators constructed from the fundamental set of spin operators {s 1 , s 2 , s 3 } for particles of arbitrary spin j.
We call this family of measurement operators spinometers, and we will characterize their properties in such sufficient detail that in Section 4.2 we will be able to simulate the single-spin MRFM experiment both numerically and by closed-form analysis. Most of these results have not appeared previously in the literature.  (76) and physically illustrated in Fig. 9(a), ergodic spinometers are constructed by applying the following tuning Note that immediately following photon detection, and independent of the channel on which the photon is detected, a compensating Hamiltonian θs op is applied to cancel the mean Stark shift that was noted in Section 3.3.4. The fluctuating portion of the Stark shift is not thereby cancelled, and it follows that ergodic spinometers are well-suited for simulating physically realistic noise, e.g., random magnetic fields that decohere spin states. We construct batrachian spinometers from ergodic spinometers by adding a downstream coupler, as shown in Fig. 9(b). Our phase and tuning conventions are: The upper-right output port is tuned to be as dark as possible, such that detection clicks occur only sporadically; in experimental interferometery this is called dark port tuning. Each click that a dark port records is accompanied by a discrete jump in the wave function, hence the name "batrachian" for these measurement operators. We will see that batrachian tuning is well-suited to the analysis of data statistics. We construct synoptic spinometers similarly, but with a different phase tuning, as shown in Fig. 9(c). Algebraically our tuning convention is We will see that synoptic spinometers do provide information about the quantum statehence the name "synoptic"-and also that they compress quantum trajectories.
We can now discern the general strategy of QMOR analysis: we model physical noise in terms of ergodic operators, we predict data statistics by the analysis of batrachian operators, and we compress simulated trajectories by applying synoptic operators.

Spinometers as agents of trajectory compression
The following derivations assume a knowledge of basic quantum mechanics at the level of Chapters 2 and 8 of Nielsen and Chuang [146] (an alternative text is Griffiths [95]), knowledge of coherent spin states at the level of Perelomov [153, eqs. 4.3.21-45] (alternatively see Klauder [118] or del Castillo [58]), and knowledge of stochastic differential equations at the level of Gardiner [84, sec. 4.3] (alternatively see Rogers [168]).
Generally speaking, the design rules of this section were first heuristically suggested by the Hilbert ontology of Section 3.2, then confirmed by numerical experiments, and finally proved by analysis. Some of the lengthier proofs would have been difficult to discover otherwise; this shows the utility of the Hilbert ontology backed-up by numerical exploration.
Nowhere in the derivations of this section will we make any assumption about the dimensionality of the Hilbert space in which the trajectory {|ψ n } resides. Therefore we are free to regard |ψ n as describing a spin-j particle that is embedded in a larger multi-spin Hilbert space. Thus all the theorems and calibrations that we will derive will be applicable both to the single-spin MRFM Hilbert space of Section 4.2 and to the large-dimension "spin-dust" spaces that we will discuss in Section 4.5.

Spinometers that einselect eigenstates
We define a uniaxial spinometer to be a measurement process associated with a single pair of measurement operators having generator s op . We can regard s op as an arbitrary Hermitian matrix, since in a uniaxial measurement there are no other operators for it to commute with. We consider ergodic, synoptic, and batrachian tunings as defined in (80)(81). Without loss of generality we assume tr(s op ) 2 = dim |ψ , i.e., the mean-square eigenvalues of s op are unity, which sets the scale of the coupling θ.
For a general state |ψ n and general Hermitian operator s op , we define the operator variance ∆ n (s op ) to be ∆ n (s op ) = ψ n | (s op − ψ n |s op |ψ n ) 2 |ψ n (83) remarking that in a finite-j Hilbert space Physically speaking, the smaller the variance ∆ n (s op ), the smaller the quantum fluctuations in the expectation value ψ n |s op |ψ n . We will now calculate the rate at which measurement operators act to minimize this variance.
Considering an ensemble of simulation trajectories, we define the ensemble-averaged variance at the n-th simulation step to be E[∆ n (s op )]. The algorithm of Fig. 1 evolves this mean variance according to For compactness we write the increment of the variance as δ∆ n (s op ) ≡ ∆ n+1 (s op )]−∆ n (s op ). Then for ergodic, synoptic, and bratrachian tunings the mean increment is batrachian tuning.
Each term in the sequence {E[∆ 1 ], E[∆ 2 ], . . .} is nonnegative by (83), and yet for synoptic and batrachian tuning the successive terms in the sequence are non-increasing (because in (86) the quantities ∆ 2 n (s op ) and F n (s op ) are nonnegative and there is an overall minus sign acting on them); the sequence therefore has a limit. For synoptic tuning the limiting states are evidently such that ∆ n (s op ) → 0, while for batrachian tunings F n (s op ) → 0, which in both cases implies that the limiting states are eigenstates of s op . This proves Design Rule 3.7. Uniaxial spinometers with synoptic or batrachian tunings, but not ergodic tunings, asymptotically einselect eigenstates of the measurement generator.

Convergence bounds for the einselection of eigenstates
We now prove a bound on the convergence rate of Design Rule 3.7. For QMOR purposes, this bound provides an important practical assurance that an ensemble of uniaxially observed spins never becomes trapped in a "dead zone" of state space.
To prove the convergence bound, we notice that in the continuum limit θ 1 the increment (86) can be regarded as a differential equation in simulation time t ≡ n δt.
This implies that the large-n variance is O(n −1 ). This proof nowhere assumes that the initial ensemble is randomly chosen; therefore the above bound applies to all ensembles, even those whose initial quantum states are chosen to exhibit the slowest possible einselection. We conclude that for synoptic tuning the approach to the zero-variance limit is never pathologically slow. We have not been able to prove a similar bound for batrachian tuning, but numerical experiments suggest that both tunings require a time t ∼ δt/θ 2 to achieve einselection. Proofs of stronger bounds would be valuable for the design of large-scale QMOR simulations.

Triaxial spinometers
We now consider triaxial spinometers, in which three pairs of synoptic measurement operators (82) are applied, having as generators the spin operators {s x , s y , s z }, applied with couplings {θ x , θ y , θ z }.
3.4.6. The Bloch equations for general triaxial spinometers In the general case we take θ 1 = θ 2 = θ 3 . We define x n = {x n , y n , z n } = j ψ n |s|ψ n to be the polarization vector at the n'th simulation step. This vector is normalized such that |x n | ≤ 1, with |x n | = 1 if and only if |ψ n is a coherent state. We further define δx n = x n − x n−1 . Taking as before E[. . .] to be an ensemble average over simulations, such that the density matrix of the ensemble is ρ n = E[|ψ n ψ n |], and therefore E[x n ] = j tr sρ n , we readily calculate that the Bloch equation that describes the average polarization of the ensemble of simulations is Since it depends only linearly upon ρ n , the above expression is invariant under the Utransform of the Theorema Dilectum. We are free, therefore, to regard our spinometers as being ergodically tuned (80), such that the simulation can be equivalently regarded, not as three competing axial measurement processes, but as independent random rotations being applied along the x-axis, y-axis, and z-axis. The above Bloch equation therefore has the functional form that we expect upon purely classical grounds.
3.4.7. The einselection of coherent states Now we confine our attention to balanced triaxial spinometers, i.e., those having with θ 1 = θ 1 = θ 1 ≡ θ, such that no one axis dominates the measurement process. Numerical simulations suggest that for synoptically tuned measurement processes, in the absence of entangling Hamiltonian interactions, simulated quantum trajectories swiftly converge to coherent state trajectories, regardless of the starting quantum state. We adopt Zurek's (exceedingly useful) concept of einselection [201] to describe this process. We now prove that synoptic spinometric observation processes always induces einselection by calculating a rigorous lower bound upon the rate at which einselection occurs.
Given an arbitrary state |ψ , we define a spin covariance matrix Λ n to be the following 3 × 3 Hermitian matrix (of c-numbers): (Λ n ) kl ≡ ψ n |s k s l |ψ n − ψ n |s k |ψ n ψ n |s l |ψ n .
This matrix covariance is a natural generalization of the scalar variance ∆ n (s op ) (83), and in particular it satisfies a trace relation that is similar to (84) tr Λ n = j if |ψ n is a coherent spin state, > j otherwise.
Here a coherent spin state is any spin-j state |x , conventionally labeled by a unit vectorx, such that x|s|x = jx (see, e.g., Perelomov [153, eq. 4.3.35]). The algorithm of Fig. 1 evolves the mean spin covariance according to For compactness we define the Λ-increment δΛ n ≡ Λ n+1 − Λ n . Then by a series expansion of (92) similar to that which led from (85) to (86)-but with more indices-we find that for ergodic, synoptic, and batrachian tuning the mean increment is batrachian tuning.
The "see text" for batrachian tuning indicates that we have found no closed-form expression simpler than several dozen terms; numerical experiments show that for this tuning the covariance exhibits random jump-type fluctuations that seemingly have no simple limiting behavior. In contrast, synoptic tuning's increment has a strikingly simple analytic form, which was guessed as an ansatz and subsequently verified by machine algebra. Proceeding as in the proof of Theorem 1, and temporarily omitting the subscript n for compactness, we now prove that for Λ computed from |ψ by (90), the scalar quantity tr Λ·Λ is non-negative for all |ψ , and vanishes if and only if |ψ is a coherent state. We remark that this proof is nontrivial because tr Λ·Λ is by no means a positive-definite quantity for general Hermitian Λ; an example is Λ = 0 +i −i 0 , for which tr Λ·Λ = tr −1 0 0 −1 = −2. Because Λ is a Hermitian 3 × 3 matrix, it can be decomposed uniquely into a real symmetric matrixΛ and a real vector v = ψ|s|ψ /j by Λ ik =Λ ik +i/2 3 l=1 ikl v l . Because the increment (93) is a scalar under rotations, without loss of generality we can choose a reference frame having basis vectors {x,ŷ,ẑ} such that {v·x, v·ŷ, v·ẑ} = {0, 0, z} and z ≥ 0. In this reference frame, the following decomposition is valid for any Hermitian matrix Λ (i.e., it holds forΛ an arbitrary symmetric matrix and z an arbitrary real number): where the residual terms p a , p b , . . . , p f are We now prove that each term in this decomposition is non-negative. The terms p a , p b , p c , and p d are non-negative prima facie. The term p e vanishes for arbitrary |ψ in consequence of the spin operator identityΛ 11 +Λ 22 +Λ 33 + j 2 z 2 = ψ|s 2 1 + s 2 2 + s 2 3 |ψ = j(j + 1). That the remaining terms are non-negative in general follows from the spin operator inequalities −j ≤ ψ|s 3 |ψ ≤ j and 0 ≤ ψ|s 2 3 |ψ ≤ j 2 , which together with our reference-frame convention imply the inequalities −1 ≤ z ≤ 1 andΛ 33 < j 2 (1 − z 2 ).
Next, we show that the sum of terms (94) vanishes if and only if |ψ is coherent, i.e., if and only if |ψ = |ẑ . It is a straightforward exercise in spin operator algebra to show that |ψ = |ẑ if and only if all of the following are true: z = 1,Λ 33 = j 2 ,Λ 11 =Λ 22 and Λ 12 =Λ 23 =Λ 13 = 0; it follows that (94) vanishes if and only if |ψ = |ẑ . By reasoning similar to Theorem 1, we conclude: Design Rule 3.8. Triaxial spinometers with synoptic tunings asymptotically einselect coherent spin states.
3.4.8. Convergence bounds for the einselection of coherent states We now exploit the identity (94) to prove a bound on the convergence of Design Rule 3.8. Our strategy is similar to our previous proof of the bound for Design Rule 3.7. Substituting the identity j 2 (1 − z 2 ) = tr Λ n − j in the first two terms of (94), and taking into account the nonnegativity of the remaining terms p a , p b , p c , p d , and p e , we obtain the following quadratic inequality in (tr Λ n − j): As an aside, our starting identity (94) was devised so as to imply a general inequality having the above quadratic functional form, in service of the proof that follows, but we have not been able to prove that the above coefficients { 1 2 , 1 4 } are the largest possible. Upon taking an ensemble average of the above inequality, followed by substitution in (93), followed by a continuum-limit integration, we obtain the following convergence bound for the ensemble-averaged trace covariance: It is instructive to restate this bound in terms of simulation time t = n δt. Taking the continuum limit Λ n → Λ(t), noting that the timescale T 1 = δt/θ 2 is the conventional T 1 that appears in the Bloch equations (89), defining for compactness of notation the initial trace covariance to be κ 0 ≡ tr E[Λ(0)] − j, and assuming for the sake of discussion that κ 0 1 (i.e, we assume that the initial ensemble is far-from-classical) the functional form of the above bound exhibits three asymptotic intervals, whose t-dependence is respectively O(1), O(1/t), and O(exp(−2t/T 1 )): The O(1) and O(1/t) behavior is functionally similar to the convergence bound established in (88) for the eigenvalue variance of Theorem 1, namely, an initial linear decrease, followed by an O(1/t) fall-off. Unique to triaxial spinometry (as far as the authors know) is the final exponentially rapid convergence to a coherent state.
We note that convergence is complete within a time ∼ T 1 that is independent of both the spin quantum number j and the overall dimensionality of the Hilbert space in which the spin is embedded. As with Theorem 1, this is a worst-case bound that applies to all ensembles, including (for example) exotic ensembles initialized with "Schröedinger's cat" states. More particularly, it applies to large-dimensional ensembles in which each of n spins in a Hilbert space of overall dimension (2j + 1) n is synoptically observed.

Implications of einselection bounds for quantum simulations
We now begin to have a quantitative appreciation of the geometric assertion of Fig. 4(i), that quantum simulations can be regarded as theaters in which the trajectory compression of synoptic observation opposes the creation of entanglement by Hamiltonian dynamics, with the balance between compression and expansion determining the dimensionality of the QMOR state space required for accurate simulation.
Even stronger convergence bounds than those we have proved would be valuable in designing QMOR simulations. Especially useful would be more tunings in which noise is realized as an entangled measurement. Physically speaking, an entangled measurement is performed by interferometrically splitting a photon along n paths, scattering the photon from a different spin along each path, then recombining and measuring the photon by freely choosing among any of an exponentially large set of braidings and interferometric couplings of the downstream optical fibers.
The analysis of such noise-equivalent tunings would require mathematical methods considerably more sophisticated than those we have deployed in this article. A known consequence of the Holevo-Schumacher-Westmoreland (HSW) theorem (which is the quantum analog of the Shannon channel capacity theorem) is that entangled measurements are necessary to maximize the information capacity of quantum channels [146, sec. 12.3.2].
If we hypothesize that quantum trajectory compression is in some sense proportional to information extracted by measurement, then the HSW theorem tells us that entangled measures will be more effective for QMOR purposes than the single-spin measures that we consider in this article. It is likely, therefore, that the search for more efficient QMOR techniques will benefit considerably from continued progress in quantum information theory.

3.4.10.
Positive P -representations of the thermal density matrix Now we focus upon control and thermodynamics. Fort the thermal axis defined in (61), we modify the synoptic spinometer matrices such that where α is the control gain. We will call this a closed-loop triaxial spinometer with unitary feedback, because (as we will see) the unitary operators exp(±iαθt × s) act cumulatively to align the spin axis witht. Closing the control loop does not alter the coherent einselection because the sole effect of a post hoc unitary operator on σ n is a spatial rotation. Since tr σ n ·σ n is a rotational scalar, (93) still holds. Thus we have Design Rule 3.9. Closed-loop triaxial spinometers with unitary feedback asymptotically einselect coherent states.
The density matrix ρ of an ensemble of closed-loop triaxial spinometer simulations is described by sequence {ρ 1 , ρ 2 , . . .} whose increment is By a straightforward (but not short) calculation we find that δρ n vanishes for ρ n = ρ th if and only if the closed-loop gain α satisfies The following two trigonometric identities hold for either choice, and will be used in Section 4.1 to establish that the choice is immaterial in practical numerical simulations.
1/α + α = −2 coth 1 2 β and 1/α − α = −2 csch 1 2 β Defining as usual the dimensionless temperature T = 1/β, we see that an optimal control gain α → ±1 establishes a temperature T → ∓0, while a control gain |α| = 1 establishes a finite temperature. We will establish later on that ρ th solves δρ n = 0 uniquely, because the Fokker-Planck equation for ρ has a unique stationary solution (thus the approach of the density matrix ρ to thermodynamic equilibrium never "stalls" or becomes trapped at false solutions). These results prove Design Rule 3.10. The density matrix of an ensemble of closed-loop triaxial spinometer simulations is asymptotically thermal.

The spin-1/2 thermal equilibrium Bloch equations
The special case of spin-1/2 particles in thermal equilibrium often arises in practice. Setting the polarization axist =ẑ, and allowing independent spinometric couplings {θ x , θ y , θ z } as in (89), we find that the finite-temperature synoptic measurement operators (99) imply the following asymmetric Bloch equations (valid for j = 1/2 only): As expected on thermodynamic grounds, we see that the equilibrium polarization is E[z] = − tanh 1 2 β. These equations are a generalization of the usual Bloch equations, in the sense that the relaxation rates along the x-, y-, and z-axes can differ independently. We remark that for j > 1/2 the thermal Bloch equations do not have a closed analytic form; that is why this more general case is not considered here.
3.4.12. The spinometric Itô and Fokker-Planck equations Now we focus on Itô and Fokker-Planck equations, aiming by our analysis to obtain both the already-validated positive Prepresentation of Design Rule 3.2 and (in the large-j limit) both the linear Design Rules 3.3-5 and the fundamental quantum limits of Design Rule 3.6.
We define a binary data three-vector d n = (d 1 n , d 2 n , d 3 n ) by is the spinometer gain. We define a zero-mean stochastic variable W n by such that (to leading order in θ) W n has the second-order stochastic properties of a discrete Wiener increment: Then via an identity valid for |ψ n a coherent state, the spinometer increments (99-b) are equivalent to an Itô increment For the drift vector a and diffusion matrix b we find Design Rule 3.9, asserts that the Itô increment (111) confines the trajectory of x n to the unit sphere. The mean increment of the m'th radial moment |x| m must therefore vanish when |x| = 1. We check this by direct calculation, finding which indeed vanishes for the a and b of (112a-b). By well-known methods [84], the Itô increment (111) immediately yields a Fokker-Planck equation for the PDF p j (x). Setting z =x·t we obtain the stationary state equation which (when properly normalized) has a unique solution in which we see that the symmetry α → 1/α is indeed respected. Substituting α + 1/α = −2 coth 1 2 β per (103), and adjusting the normalization to match the P -representation convention (59) yields Design Rule 3.2.
3.4.13. The standard quantum limits to linear measurement To connect these results to Design Rules 3.2-3.6, we first write the Itô equation (111) in Langevin form by substituting where r = 1/δt is the rate at which increments occur, and x N (t) is white noise with crosscorrelation Then (111) becomes the integral of the Langevin equatioṅ where x M (t) = x(t) + x N (t) is the measured spin axis. We see that x(t) is dynamically attracted toward the measured axis x M (t). Even openloop spinometers exhibit this attraction, since for α = 0 we finḋ We remark that in uniaxial spinometry we saw that a similar einselection-by-attraction generates the "collapse" of |ψ n to an eigenstate, as described by Design Rule 3.7. This attraction is of course a fundamental tenet of the Hilbert ontology of Section 3.2.7. We now transform (121b) to the second-order Newtonian equation of an oscillator. To do this, we introduce a spring k and frequency ω 0 by defining the operators We confine our attention to those coherent states that have z −1, which with our sign conventions means systems having positive inverse temperature β, negative Hilbert feedback gain α, and oppositely directed spinx and polarizationt, such thatx ·t −1. For these states the canonical commutator [q op , p op ] = −i s 3 /j i holds in the large-j limit. Defining the coherent oscillator coordinate q(t) to be we find that (121b) takes the linearized Newtonian form Here the spring k, mass m = k/ω 2 0 , and coordinate q are to be understood in a generalized sense in which the system energy is 1 2 mq 2 + 1 2 kq 2 . The measurement noise q n (t) is given from (123) in terms of the spinometer noises x N (t) and y N (t) of (120) by and we find from (120) and (123) that the measurement noise q n (t) has a PSD S q n of The force noise f n (t) is then determined from (121b) to be f n (t) = γH[q n (t)], where H is the Hilbert transform and the Hilbert gain γ is found to be 3.4.14. Multiple expressions of the quantum noise limit It follows from the preceding results that the PSD of the spinometric force noise f n (t) can be expressed in multiple equivalent forms: = γ force noise ∝ Hilbert backaction gain (128c) = 4k rjθ 2 /ω 0 raw spinometer parameters (128d) Each of the above relations has a plausible claim to expressing the "most natural" or "most fundamental" relation between measurement noise and force noise . . . despite the fact that no two physical interpretations are the same, and even though the interpretations given (128a) and (128b) seem contradictory. We further see by Design Rule 3.6 that these spinometric relations saturate the Hilbert noise limit (γS q n = ), the Braginsky-Khalili limit (S q n S f n = 2 ), and the Hefner-Haus-Caves limit (NF = 2); thus in some sense all of these fundamental quantum limits are embodied in the above family of spinometric relations. Acknowledging the self-consistency of this diversity, and appreciating its mathematical origin in the diversity of equivalent noise models that are supported by the Theorema Dilectum, helps us appreciate how the quantum noise literature can be so immensely large, and support so many different notations, physical arguments, and conclusions, and yet maintain its internal consistency.
In a teaching environment, it is not practical to sustain a dispassionately anarchical equality among physical interpretations (128a-d). This article's Hilbert ontology (Section 3.2.4) designates (128a) to be the fundamental relation, because it embodies the central Hilbert tenet that "measurement noise always back-acts upon system dynamics in such a way as to bring the state of the system into agreement with the measurement." This choice is justified solely because it yields useful guiding principles for efficient quantum simulations.

Summary of the design rules
In summary, we have established by Design Rules 3.7, 3.8 and 3.9 the quantum mechanism by which synoptic noise processes compress simulated quantum trajectories onto lowerdimension GK manifolds (as was promised in Section 1.5.11). We have established by Design Rule 3.10 that the effects of thermal reservoirs can be modeled as equivalent processes of covert measurement and control (as was promised in Section 1.5.12). And we have established by Design Rule 3.6 that the Hefner-Haus-Caves, Braginsky-Khalili, and Hilbert quantum noise limits are respected by QMOR simulations (as was promised in Section 3.2.6).
The focus of the remainder of this article is to show, by explicit examples, that these design rules are sufficient to "enable the reader to proceed to the design and implementation of practical quantum simulations, guided by an appreciation of the geometric and informatic principles that are responsible for the overall simulation accuracy, robustness, and efficiency" (as was promised in the Introduction).

Examples of quantum simulation
Now we turn our attention toward applying the preceding results in implementing practical quantum simulations.

Calibrating practical simulations
Our simulations provide data via the binary stream of defined in (106), which is low-pass filtered to produce a classical data record. We now work through, in detail, the process of computing and calibrating this data stream, and ensuring that it is numerically wellconditioned.
We begin by considering the problem of determining, from physical system parameters, the measurement operation parameters {θ, α} in (99) and the clock rate r = 1/δt. In essence this calibration process requires that we invert systems of equations that include the Bloch equations (89) and (105), the Langevin equation (121b), and the mapping of spinometer parameters onto oscillator parameters (refeq: q definition).
Whenever our simulations include a projective step, we must also ensure that theparameter of the projected Theorema Dilectum (54) satisfies 1. Physically speaking, imposing the small-condition ensures that the simulated trajectories evolve by drift and diffusion, rather than by "quantum jumps" that may be projectively ill-conditioned. Equally importantly, it ensures that they respect the informatic causality that is guaranteed by the Theorema Dilectum.
4.1.1. Calibrating the Bloch equations A common system to be simulated is a spin j = 1 2 in contact with a thermal reservoir. We desire that the three thermal relaxation rates along the x-, y-, and z-axes be {Γ x , Γ y , Γ z } = {1/T x , 1/T y , 1/T z } and that the equilibrium thermal density matrix be ρ th ∝ exp(−βs z ), such that the equilibrium spin polarization p 0 is − tanh 1 2 β as in (105). We thus have four physical parameters {Γ x , Γ y , Γ z , β} with which to determine five raw spinometer parameters {θ x , θ y , θ z , α, r}. Needing one more physical parameter, we note that the -parameters of the GK-projected Theorema Dilectum (54) is given for j = 1/2 spinometers from (105) by 2 = θ 2 (1+α 2 )/4, and so we impose as our fifth condition 0.1 for all three spinometers (the cut-off 0.1 yielding in our experience well-converged numerical results). The equations below then follow from (101), (103), and (105): Bloch-parameter calibration proceeds as follows. We first determine the spinometer gain α from (129a), and we will see that the choice "gain too big" versus "gain too small" is immaterial. The value of the spinometer click-rate r is then set from (129b-d) by requiring that min{ x , y , z } 0.1 for the reasons noted above. The values of the three spinometer phases {θ x , θ y , θ z } are then determined from (129e). We remark upon three features. First, we see that insofar as simulation efficiency and numerical conditioning are concerned, the choice between the two options for the feedback gain α in (129a) is immaterial, since according to the above construction the spinometer click-rate r and the -parameters are unaffected. Second, for the special case T z = T 1 and T x = T y = T 2 the above results reduce to the usual Bloch equations. Third, the positivity of the -parameters in (129b-d) requires that the Bloch relaxation rates satisfy the inequality The authors suspect that the above Bloch inequality is tight, in the sense that no spin-1/2 Lindblad-form master equation can violate it, but we have not proved this. We remark that the above triaxially asymmetric Bloch equations and their associated relaxation rate inequality have (to the best of our knowledge) not appeared in the literature before.

Calibrating test-mass dynamics in practical simulations
Now we consider test-masses (e.g., MRFM cantilevers) in contact with thermal reservoirs. Calibration proceeds by a line of reasoning similar to the above. We take the test-mass to be described by two physical parameters: the (dimensionless) temperature β = ω 0 /(k B T ) and the (dimensionless) quality Q of the ring-down wave-form q(t) ∝ cos(ω 0 t) exp(ω 0 t/(2Q)). The four spinometer parameters to be determined are {j, θ, α, r}. The spin number j we take to satisfy j 1, and we will find that the precise value chosen for j is immaterial. For coherent states with z −1 as discussed in Section 3.4.12, the -parameter is found to be jθ 2 (1 + α 2 ). Calibration proceeds as follows: We first determine the spinometer gain α from (131a). The value of the spinometer clickrate r is then set from (131b) by requiring that 0.1 as in the Bloch equation case, and the spinometer phase θ is determined from (131c). We remark again that insofar as simulation efficiency is concerned, the choice between the two options for the feedback gain α in (131a) is immaterial, since the spinometer click-rate r is unaffected. We see also that for fixed quality Q, the simulation rate r is O(coth 1 2 β), which physically means that hot cantilevers are computationally more expensive to simulate than cold ones, as is reasonable.

Calibrating purely observation processes
It can happen that we wish to directly observe a spin-1/2 particle along a single axis, nominally the z-axis, in an observation process in which the measured z-axis polarization z m (t) = z(t) + z n (t) has a specified (onesided) noise PSD S z n . No thermodynamical feedback is applied. Then by reasoning similar to the preceding cases, we find from (120) that the calibration relations are and as before, the spinometer click rate r is determined by requiring 0.1. Similarly, if we wish to simulate the continuous observation of a test-mass coordinate q(t), with no thermodynamical feedback, the required calibration equations are given from (120) and (123) in terms of the measurement noise PSD S q n by where again the spinometer click rate r is determined by requiring 0.1. We remark that in all of the above cases the raw binary data stream (106) of the simulation must be low-pass filtered in order to obtain a (noisy) data record that will have the above statistical properties within the filter passband having. This filtering closely models the way that real experimental signals (for example, a continuously measured cantilever motion) are displayed upon oscilloscopes.

Three single-spin MRFM simulations
With reference to Fig. 10, we now turn our attention to the simulation of the IBM singlespin MRFM experiment [170] . We will initially present the simplest possible class of simulations that reproduce the data of that experiment, postponing a discussion of more detailed simulations until Section 4.4. Our goal is to illuminate the central role of the Theorema Dilectum in answering the question that was raised in Section 1.1: "How does the Stern-Gerlach experiment work?" All three columns of Fig. 10 show a simulated thirteen-hour experiment (the length of the IBM experiment). The time-spacing δt = 1/r between spinometric clicks is set to 7.1 ms; thus approximately 6.6 × 10 6 time-steps were simulated. In each column, the experimental data are simulated as arising from three competing spinometric processes. Spin relaxation was simulated by x-axis and y-axis spinometers having θ x = θ y = 0.093. The consequent spin relaxation time from (105) is T z = 2/ r(θ 2 x + θ 2 y ) 0.76 s as was observed in the IBM experiment. Measurement effects were simulated by a z-axis spinometer having θ z = 0.026, and the consequent measurement noise PSD from (120) is S z n = 2/(rθ 2 z ), which numerically corresponds to a noise level of 11.5 Bohr magnetons of noise in one root-Hertz of bandwidth, as was observed in the IBM experiment.
For visualization purposes only, all time-domain data streams shown in Fig. 10 were low-pass filtered with a time constant τ = T z = 0.76 s.
The following discussion is insensitive to the above experimental details, and applies to all experiments of Stern-Gerlach type in which the signal-to-noise ratio (SNR) is low, such that continuous monitoring over extended periods of time is required to observe the effect.
In the next three sections, we simulate the spin relaxation of the IBM single-spin MRFM experiment by three different unravelings: batrachian, ergodic, and synoptic. We will see that the three unravelings lead to three very different classes of quantum trajectories, and hence, three very different answers to the question "How does the Stern-Gerlach experiment work?" Nonetheless, as guaranteed by the Theorema Dilectum, we will find that the simulated experimental data are identical for all three unravelings. We now work through the mathematical and physical details of how this comes about.

A batrachian single-spin unraveling
The left-hand column of Fig. 10 shows a simulation in which thermal noise is unravelled as a batrachian process, whose measurement operations are given algebraically in (81) and which are depicted in hardware-equivalent form in Figure 9(b). This is by far the easiest simulation to analyze in closed form: the spin polarization jumps randomly between ±1, driven by the batrachian jumps of the thermal reservoir, while being continuously measured by the (noisy) cantilever.
The simulated data stream is therefore a random telegraph signal with added white noise, such that the mean-square quantum spin polarization inferred from the data is unity. We conclude that from the batrachian point of view, the Stern-Gerlach effect (meaning, that the mean-square spin polarization is measured to be unity) comes about because noise is a quantized jump process, such that the mean-square spin polarization always is unity.

An ergodic single-spin unravelling
The middle column of Fig. 10(b) shows a simulation in which thermal noise is unravelled as an ergodic process, whose measurement operations are given algebraically in (80) and which are physically depicted in Figure 9(a). Physically speaking, the spin polarization is driven by random magnetic fields, such that the mean-square quantum polarization is 1/3. Now a subtle effect comes into play. The z-axis measurement process back-acts upon the spin state, such that whenever an "up" fluctuation in the data is observed, the spin state is "dragged" toward a positive polarization. This effect is evident in the simulated data. The consequence of state-dragging back-action is that the measured mean-square polarization is larger than the mean-square polarization of the underlying quantum state.
It would be quite a complicated task to calculate the resulting data statistics from (e.g.) the appropriate Itô, Langevin, and Fokker-Planck equations. Fortunately, the Theorema Dilectum does this mathematical work for us: the data statistics are guaranteed to be exactly the same random telegraph statistics as in the Batrachian case.
We conclude that from the ergodic point of view, the Stern-Gerlach effect (meaning, that the mean-square spin polarization is measured to be unity) comes about because measurement is a Hilbert process (meaning, it accords with the state-dragging Hilbert back-action ontology of Section 3.2.4).

A synoptic single-spin unravelling
The right-hand column of Fig. 10(c) shows a simulation in which thermal noise is unravelled as an synoptic process, whose measurement operations are given algebraically in (82) and which are physically depicted in Figure 9(c).
In synoptic unravelling, all processes are measurement processes, and each process seeks to align the spin polarization along its own axis. In our simulation, the x-axis and y-axis measurement processes are considerably stronger than the z-axis process. In consequence of the Hilbert state-dragging effect, the spin polarization now points predominantly in the equatorial direction, such that the mean-square quantum polarization is only ∼ 0.05.
Again it would be quite a complicated task to calculate the resulting data statistic from Itô equations, etc., and again the Theorema Dilectum does this mathematical work for us: as in the preceding two cases, the data statistics are random telegraph statistics with added white noise. We conclude that from the synoptic point of view, the Stern-Gerlach effect is not associated with "wave function collapse," but rather comes about (as in the ergodic case) because measurement is a Hilbert process.

So how does the Stern-Gerlach effect really work?
We are now in a position to answer more completely the question "How does the Stern-Gerlach effect really work?" We answer as follows: "Nothing definite can be said about the internal state of noisy systems, either at the classical or at the quantum level. It is best to pick an ontology that facilitates rapid calculations and suggests interesting mathematics. For purposes of large-scale quantum simulation, a particularly useful ontology is one in which all noise processes are conceived as equivalent covert measurement processes. In this ontology, the Stern-Gerlach effect works because competing measurement processes exert a Hilbert back-action mechanism that 'drags' quantum states into agreement with measurement. In consequence of these competing Hilbert measurements, experimental data having the statistics of random telegraph signals are obtained even when no quantum jumps are present." Of course, we saw in the three simulations of Fig. 10a-c that other explanations are perfectly reasonable, and that is why in Section 1.1.1 we embraced Peter Shor's maxim: "Interpretations of quantum mechanics, unlike Gods, are not jealous, and thus it is safe to believe in more than one at the same time," to which we now append the caveat "provided that all interpretations respect the fundamental mathematical and physical invariance of the Theorema Dilectum."

Was the IBM cantilever a macroscopic quantum object?
The least realistic element of the proceeding simulations is the modeling of the cantilever as a single z-axis spinometer having quantum number j = 1/2. A more realistic model would have treated the cantilever as a large-j quantum object subject both to thermal noise processes and to experimental measurement processes. However, we can appeal to the Theorema Dilectum to show that these refinements would not change the simulated data at all. The reason is that both the cantilever thermal reservoir and the experimental (interferometric) cantilever measurement process can be modeled as synoptic processes that compress the cantilever's quantum state to a coherent state. Then modeling the spin relaxation as a batrachian process, the output of the resulting (effectively semi-classical) batrachian simulation will be precisely the random telegraph signal that was obtained in the simpler batrachian simulation of Section 4.2.1 above.
It follows by the Theorema Dilectum that all quantum simulations of the cantilever, even elaborate large-j simulations in which a non-coherent quantum cantilever state is entangled with the quantum spin state, will simulate the same random telegraph data statistics as the simpler simulations already given, and in particular, will yield an observed mean-square polarization of unity.
This leads to an interesting question: what was the "real" quantum state of the IBM cantilever? We have seen that this question has a well-posed answer only insofar as there is agreement upon the "real" noise and observation processes acting upon the cantilever, such that the tuning ambiguity of the Theorema Dilectum does not come into play.
If we stipulate that the "real" cantilever thermal noise and the "real" spin relaxation are due to ergodic physical processes, then the IBM experiment can only be "really" described in terms of a spin-state that is quantum-entangled with the cantilever state, in which the observed mean-square polarization of unity is "really" due to the state-dragging Hilbert back-action associated with the cantilever measurement process. In other words, the IBM experiment "really" observed the cantilever to be a macroscopic quantum object.
As quantum objects go, the IBM cantilever was exceptionally large [170]: its resonant frequency was ω 0 /(2π) = 5.5 kHz, its spring constant was k = 0.011 mN/m, and its motional mass was m = k/ω 2 0 = 9.1 pg. The preceding paragraphs are an argument for regarding this cantilever to be among the stiffest, slowest, most massive dynamical systems whose quantum nature has been experimentally confirmed. Such measurements are significant from a fundamental physics point of view, in probing the limits of quantum descriptions of macroscopic objects, as reviewed by Leggett [125,126,127,128].
Any line of reasoning that is as brief as the preceding one, about a subject that is as subtle as macroscopic quantum mechanics, is sure to have loopholes in it. A major loophole is our modeling of decoherent noise as a Markovian process. As reviewed by Leggett et al. [129], spin decoherence in real experiments is (of course) due to non-Markovian quantumentangling interactions. We now turn our attention to the algorithmic and numerical challenges of simulating such systems.

The fidelity of projective QMOR in spin-dust simulations
As test cases, we computed what we will call spin-dust simulations. Spin-dusts are quantum systems that are deliberately constructed so as to have no symmetries or spatial ordering. Their sole purpose is to provide a well-defined test-bed for numerical and analytic studies of the fidelity of projective quantum model order reduction.
Spin-dusts couple pairs of spin-1/2 particles {j, k} via a dipole-dipole interaction Hamiltonian H jk that is given by The unit vectors n jk are chosen randomly and independently for each {j, k}, and we note that self-coupling is allowed. Physically we can think of spin-dusts as broadly analogous to-but less structured than-systems such as the interacting spins in a protein molecule. In our simulations each spin is randomly coupled to four other spins, in addition to its self-interaction. Then it is easy to show that tr H = 0 and tr H 2 / dim H = n spin , which is to say, the per-spin energy of our spin-dusts has zero mean and unit variance. The time-scale of the spin dynamics of the system is therefore unity. We further stipulate that each spin is subject to a triaxial spinometric observation process having relaxation time T x = T y = T z = 10. Thus the time-scale of decoherent observation is ten-fold longer than the dynamical time-scale.
Simulations were computed with a time-step δt = 0.1 and spinometric couplings θ x = θ y = θ z = 0.1, using the sparse matrix routines of Mathematica. The numerical result was an "exact" (meaning, full Hilbert space) quantum trajectory |ψ 0 (t) . These trajectories were then projected on GK manifolds of various order and rank by the numerical methods of Section 2.12. The main focus of our numerical investigations was the fidelity of the projected states |ψ K (t) relative to the exact states |ψ 0 (t) .
4.5.1. The fidelity of quantum state projection onto GK manifolds With reference to Fig. 11(a), simulations were conducted with numbers of spins n ∈ 1, 18, having random dipole coupling links as depicted. The median quantum fidelity was then computed, as a function of n, for GK rank r ∈ {1, 2, 5, 10, 20, 30} (see Fig. 3 for the definition of GK rank). Both synoptic and ergodic unravelings were simulated. Typically |ψ 0 (t) was projected at thirty different time-points along each simulated trajectory, always at times t > 100 to ensure that memory of the (randomly chosen) initial state was lost. We remark that numbers of spins n > 18 could not feasibly be simulated on our modest computer (an Apple G5). The quantum fidelity of a projected state |ψ K was defined to be [146, Section 9.2.2] As shown in Fig. 11(b), for ergodic unravellings large-n quantum fidelity fell-off exponentially, while for synoptic unravelings large-n fidelity remained high. No mathematical explanation for the observed exponential fall-off in ergodic unravelling fidelity is known. The asymptotic large-n behavior of the synoptic fidelity also is unknown. In particular, for systems of hundreds or thousands os spins, would the empirical rule-ofthumb "GK rank fifty yields high fidelity for spin-dust systems" still hold true? These are important topics for further investigation.
The achieved high-fidelity algorithmic compression was large: an 18-spin exact quantum state |ψ 0 (t) is described by 2 18 independent complex numbers, while an order-18 rank 30 GK state |ψ K (t) -as seen at lower right in Fig. 11(b)-is described by 30 × (18 + 1) = 570 independent complex variables. The dimensional reduction is therefore 460-to-1. cosines ψ|s ·m|ψ /j for randomly spins, randomly chosen trajectory points, and randomly chosen unit vectorsm. One hundred randomly chosen data points are shown. We observe that the rank-one GK manifold does an excellent job of representing the spin direction cosines, which can be regarded as (essentially) classical quantities.

4.5.
3. The fidelity of operator covariance in projective QMOR As a measure of pair-wise quantum correlation, we examined the spin operator covariance. With reference to (90), this quantity is given by which vanishes for rank-1 (product states). The second row of Fig. 12 plots Σ kl for one hundred randomly chosen trajectory points, and randomly chosen spins, having randomly chosen indices j and k. We observe that GK ranks in the range 20-50 are necessary for projection to preserve pair-wise quantum correlation with good accuracy. 4.5.4. The fidelity of quantum concurrence in projective QMOR As a measure of pairwise quantum entanglement, we examined Wooters' quantum concurrence [199]. The concur-rence is computed as follows. Let ρ AB be the reduced density matrix associated with spins A and B. Let λ i be eigenvalues of the non-Hermitian matrix ρ ABρAB , in decreasing order, wherẽ ρ AB = (σ y ⊗ σ y )ρ AB (σ y ⊗ σ y ). Then the concurrence c is defined to be It can be shown that the concurrence vanishes for product states, and that spins are pairwise entangled if and only if c > 0.
In the third row of Fig. 12, we observe that coupled spin-pairs are far more likely to be quantum-entangled than non-coupled spin-pairs, as expected on physical grounds. We further observe that GK ranks in the range 20-50 are necessary for projection to preserve concurrence with good accuracy. 4.5.5. The fidelity of mutual information in projective QMOR As a measure of pair-wise quantum information, we examined von Neumann's mutual information [146,Section 11.3], which is computed as follows. For a general density matrix ρ we define the von Neummann entropy S(ρ) = − tr ρ log 2 ρ. Then for spins A and B the mutual information is given by It is known that the mutual information vanishes for product states, and that this quantity is otherwise positive in general.
In the fourth row of Fig. 12, we observe that coupled spin-pairs share more mutual quantum information than non-coupled pairs, as expected on physical grounds. We further observe that GK ranks in the range 20-50 are necessary for projection to preserve mutual quantum information with good accuracy.
As a quantitative summary of this observation, the 15-spin simulations of Fig. 12 predict 15 single-spin density matrices ρ A and 105 pairwise reduced density matrices ρ AB (in addition to higher-order correlations). Each of the single-spin density matrices has 3 (real) degrees of freedom, and each pairwise density matrix introduces 9 more (real) independent degrees of freedom, for a total of 990 independent degrees of freedom associated with the one-spin and two-spin reduced density matrices. In comparison, the rank-50 GK manifold onto which the quantum states are projected is described by Kählerian coordinates having (it can be shown) 1600 locally independent coordinates. Using 1600 state-space coordinates to encode 990 physical degrees of freedom represents a level of MOR fidelity that (obviously) cannot be improved by more than another factor of two or so. The mathematical origin of this empirical algorithmic efficiency is not known.

Quantum state reconstruction from sparse random projections
We will conclude our survey of spin-dust simulations with some concrete calculations that are motivated by recent advances in the theory of compressive sampling (CS) and sparse reconstruction. It will become apparent that synoptic simulations of quantum trajectories mesh very naturally with CS methods and ideas. To the best of our knowledge, this is the first description of CS methods applied to quantum state-spaces.
Our analysis will mainly draw upon the ideas and methods of Donoho [61] and of Candès and Tao [37], and our discussion will assume a basic familiarity with these and similar CS articles [30,31,32,34,36], especially a recent series of articles and commentaries on the Dantzig selector [17,29,38,64,79,137,167]. Our analysis can alternatively be viewed as an extension to the quantum domain of the approach of Baraniuk, Hegde, and Wakin [11,47,194] to manifold learning [188] from sparse random projections.
Our objectives in this section are: • establish that synoptically simulated wave functions ψ 0 are compressible objects in the sense of Candès and Tao [37], • establish that high-fidelity quantum state reconstruction from sparse random projections is algorithmically tractable, • describe how nonlinear GK projection can be described as an embedding within a larger linear state-space of a convex optimization problem, and thereby • specify algorithms for optimization over quantum states in terms of the Dantzig selector (a linear convex optimization algorithm) of Candès and Tao [37].
At the time of writing, the general field of compressive sensing, sampling and simulation is evolving rapidly-"Nowadays, novel exciting results seem to come out at a furious pace, and this testifies to the vitality and intensity of the field" [38]-and our overall goal is to provide mathematical recipes by which researchers in quantum sensing, sampling and simulation can participate in this enterprise.  [37] call the design matrix to be an n × p matrix X. The projected state φ 0 = Xψ 0 is given, and our reconstruction task is to estimate the ψ 0 (the "model") from the φ 0 (the "sensor data"). In general n ≤ p and we particularly focus upon the case n p. We initialize the elements of the design matrix X with i.i.d. zero-mean unit-norm (complex) Gaussian random variables. Then the rows and columns of X are approximately pairwise orthogonal, such that X satisfies the approximate orthogonality relation XX † p I and therefore satisfies the approximate projective relation (X † X) 2 ∼ pX † X. As a remark, if we adjust X to make these orthogonality and projective relations exact instead of approximate-for example by setting all n nonzero singular values of X to unity-our sparse reconstructions are qualitatively unaltered.
In CS language, we have specified random design matrices X that satisfy the uniform uncertainty principle (UUP) [33,34,37], meaning (loosely) that the columns of X are approximately pairwise orthogonal. See [37] for a definition of UUP design matrices that is more rigorous and general.
From a geometric point of view, this means we can regard X † X-which will turn out to be the mathematical object of interest-as a projection operator from our (large) pdimensional quantum state-space onto a (much smaller) n-dimensional subspace.
We have already seen in Sections 2.12 and 4.5 that the following minimization problem can be tractably solved by steepest-descent methods: where we have adopted the CS literature's practice of specifying the l 2 norm explicitly. Here ψ κ (c) is a vector of multilinear gabion-Kähler (GK) polynomials as defined in Section 2 and depicted in Fig. 3. Inspired by the CS literature, we investigate the following CS generalization of (139): Now we are minimizing not on the full Hilbert space, but on the n-dimensional subspace projected onto by X. We recognize the right-hand expression as a nonlinear Kählerian generalization of a standard minimization problem (it is discussed e.g. by Donoho and Stodden [62, eq. 3] and by Candès and Tao [37, eq. 1.15]). To make this parallelism more readily apparent, we can write the above minimization problem in the form for some choice of c, where we have substituted φ 0 → y and introduced β as an auxiliary variable. Comparing the above to the well-known LASSO minimization problem [62,79] for some t, we see that the sole change is that the LASSO problem's l 1 sparsity constraint β l 1 ≤ t has been replaced with the GK representability constraint β = ψ κ (c). We remark upon the parallelism that both constraints are highly nonlinear in β.
But this parallelism in itself does not give us much reason to expect that the minimization (141) is tractable, since we saw in Section 2 that the space of feasible solutions ψ κ (c) is (floridly) nonconvex. Consequently, unless some "GK magic" of comparable algorithmic power to the well-known "l 1 magic" of CS theory [35] should come our rescue, there seems to be little prospect of computing the minimum (141) in practice.
Persisting nonetheless, we compute successive approximations {c 1 , c 2 , . . . , c i } by a projective generalization of the same steepest-descent method that produced the results of Figs. 11 and 12. Specifically, we expand the GK coordinates via c i+1 = c i + δc i and iterate the resulting linearized equations in δc i Empirically, good minima are obtained from O(dim K) iterations of this equation from randomly-chosen starting-points. This benign behavior is surprising, given that our objective function (140) is a polynomial in O(dim K) variables having O (dim H) 2 independent terms, because generically speaking, finding minima of large polynomials is computationally infeasible.
According to the geometric analysis of Section 2.5, the existence of feasibly computed minima is explained by the rule structure of GK state-space, which ensures that almost all state-space points at which the increment (143) vanishes are saddle points rather than local minima, in consequence of the nonpositive directed sectional curvature that is guaranteed by Theorem 2.1.
We now discuss GK geometry from the alternative viewpoint of CS theory, further developing the idea that GK rule structure provides the underlying geometric reason why CS "works" on GK state-spaces. 4.6.2. Randomly projected GK manifolds are GK manifolds With reference to the algorithm of Fig. 2, we immediately identify (A † X † XA) P in (143) as the Kählerian metric of a GK manifold having an algebraic Kähler potential (see (17)) that is simply Since X is constant, we see that the projected Kähler potential is a biholomorphic polynomial in the same variables and of the same order as the original Kähler potential. It follows that a projected GK manifold is itself a GK manifold, and in particular the GK rule structure is (of course) preserved under projection, and this means that all of the sectional curvature theorems of Section 2 apply immediately to QMOR-CS on GK state-spaces.
This GK inheritance property is mathematically reminiscent of the inheritance properties of convex sets and convex functions, and it suggests that a calculus of GK polynomials and manifolds might be developed along lines broadly similar in both logical structure and practical motivation to the calculus of convex sets and convex functions that is presented in the standard textbooks of CS [21] (we discuss this further in Section 4.6.9). 4.6.3. Donoho-Stoddard breakdown at the Candès-Tao bound Putting these ideas to numerical test, using the same spin-dust model as in previous sections, we find that random compressive sampling does allow high-fidelity quantum state reconstruction, provided that the state trajectory to be reconstructed has been synoptically unraveled (Fig. 13).
These numerical results vividly illustrate what Donoho and Stodden [62] have called "the breakdown point of model selection" and we note that Candès and Tao have described similar breakdown effects in the context of error-correcting codes [33]. Surprisingly, it does not appear to have been recognized that a similar breakdown occurs in quantum modeling whenever too many wave function coefficients are reconstructed from too few projections.
For discussion, we direct our attention to the rank-30 block of Fig. 13(f), where (as labeled) we reconstruct the p = dim H = 2 n spin = 2048 (complex) components of ψ 0 from n random projections onto a GK manifold having S = dim K = (GK rank) × (1 + log 2 (dim H)) = 30 × 12 = 360 (complex) dimensions. All six blocks (a-f) of the figure are similarly labeled with p and S.
Given our example with p = 2048 and S = 360, what does CS theory predict for the minimum number n of random projections required for accurate reconstruction? According to Candès and Tao [37] With overwhelming probability, the condition [for sparse reconstruction] holds for S = O n/ log(p/n) . In other words, this setup only requires O log(p/n) observations per nonzero parameter value; for example, when n is a nonnegligible fraction of p, one only needs a handful of observations per nonzero coefficient. In practice, this number is quite small, as few as 5 or 6 observations per unknown generally suffice (over a large range of the ratio p/n).
We naively adapt the above Candès-Tao big-O sampling bound to the case at hand by recalling that S = n/ log(p/n) implies n S log(p/S) for p n [52]. We therefore expect to observe Donoho-Stodden breakdown at a (complex) projective dimension n sb (p, S) (which we will call the sampling bound ) that to leading order in S/p is Here we recall that dim H is the (complex) dimension of the Hilbert space within which the efflorescent GK state-space manifold of (complex) dimension dim K is embedded (see Sections 2.6.4 and 2.11, and also (5), for discussion of how to calculate dim K).
The above sparsity bound accords remarkably well with the numerical results of Fig. 13. This empirical agreement suggests that QMOR and CS may be intimately related, but on the other hand, there are the following countervailing reasons to regard the agreement as being possibly fortuitous: 1. the Candès-Tao bound applies to state-spaces that are globally linear, whereas we are minimizing on a GK state-space that is only locally linear, and 2. the onset of Donoho-Stodden breakdown in Fig. 13 is (experimentally) accompanied by the onset of multiple local minima of (140), which are not present in the convex objective function of Candès and Tao, and Figure 13: Quantum state reconstruction from sparse random projections. A (typical) state from an 11-spin trajectory was reconstructed from sparse projections onto random subspaces (horizontal axis), and the resulting quantum fidelity was evaluated (vertical axis). Each point represents a single minimization of (143), by iteration of (143) with a conjugate gradient correction, from a random starting point chosen independently for each minimization. Convergence to "false" local minima was sporadically encountered for low-rank GK projections (graphs a-b, GK ranks 1 and 2) but not for higher-rank GK projections (graphs c-f, GK ranks 5, 10, 20, and 30). The onset of Donoho-Stodden breakdown was observed to occur near the Candes-Tao bound (145)-plotted as a dotted vertical line-for all GK ranks tested. The ergodic spin-dust simulation yielded states whose reconstruction properties were indistinguishable from random states, as expected.
3. high-accuracy numerical agreement with a "big-O" estimate is fortuitous; the agreement seen in Fig. 13 is better than we have a reason to expect.
So although the Candès-Tao bound seems empirically to be the right answer to "when does Donoho-Stodden breakdown occur in quantum model order reduction?" the numerical calculations do not explain why it is the right answer. We now present some partial results-which however are rigorous and deterministic insofar as they go-that begin to provide a nontrivial explanation of why the Candès-Tao bound applies in the sparse reconstruction of quantum states. The basic idea is to embed the nonlinear minimization (139) within a larger-dimension problem that is formally convex. We will show that this larger-dimension optimization problem can be written explicitly as a Dantzig selection.
The main mathematical tool that we will need to develop is sampling matrices X whose (small) row dimension is n = dim H, and whose (large) column dimension p is a power of dim H. These matrices are too large to be evaluated explicitly-they are what Cai and Lv call "ultrahigh-dimensional" [29]. A novel aspect of our analysis is that we construct these matrices deterministically, such that their analytic form allows the efficient evaluation of matrix products. 4.6.4. Wedge products are Hamming metrics on GK manifolds Let us consider how the efflorescent GK geometry that we described in Sections 1.5.8 and 2.11 can be made the basis of a deterministic algorithm for constructing good sampling matrices.
Our basic approach is to construct a deterministic lattice of points on GK manifolds, together with a labeling of the lattice for which the Hamming distance between two labels is a monotonic function solely of the wedge product between that pair of points, such that the larger the Hamming distance between two points, the closer they approach to mutual orthogonality. The problem of constructing good sampling matrices then becomes equivalent to the problem of constructing good error correction codes.
We now construct the desired GK lattice. We consider sampling matrices X whose columns are not random vectors, but rather are constrained to satisfy ψ = ψ κ (c) for some GK state-space ψ κ (c). We wish these vectors to be approximately orthogonal. To construct these vectors (and simultaneously assign each vector a unique code-word), we specify an alphabet of four characters {a, b, d, e}, we identify the four characters with the four vertices of a tetrahedron having unit vectors {n a ,n b ,n d ,n e }, and we identify the unit vectors with the four spin-j coherent states {|n a , |n b , |n d , |n e }, such that for s = {s x , s y , s z } the usual spin operators, the tetrahedron vertices aren a = n a |s|n a /j, etc. Soon it will become apparent that the vertices of any polytope, not only a tetrahedron, suffice for this construction, and that the vertices of Platonic solids are a particularly good choice.
We recall from our study of GK geometry that a wedge product (12) can be associated to each letter-pair {a, b} as follows |a ∧ b| 2 ≡ n a |n a n b |n b − n a |n b n b |n a From Wigner's identity (56) we have | n a |n b | 2 = |D j jj (0, θ ab , 0)| 2 = cos(θ ab /2) 4j where cos(θ ab ) =n a ·n b , so the spin-j wedge product is easily evaluated in closed form as which for our tetrahedral alphabet is simply Here and henceforth we have added a subscript j to all wedge products for which an analytic form is given that depends explicitly on the total spin j. Now we specify a dictionary to be a set of n-character words {w k } with each word associated with an ordered set of tetrahedral characters w k = {c k 1 , c k 2 , . . . , c k n }. We further associate with each word a petal-vector |w k Given two words {w i , w k }, their mutual Hamming distance h(w i , w k ) is defined to be the number of symbols that differ between w i and w k . We also can associate with any two petal-vectors their mutual wedge product |w i ∧ w k | defined by The Hamming distance h(w i , w k ) is a Hamming metric on our codeword dictionary, and it is easy to show that this code metric is related to the petal-vector wedge product |w i ∧ w k | by the simple expressions or equivalently for a tetrahedral petal-vector dictionary The main result of this section is the above monotonic functional relation between a Hamming distance and a wedge product. We are not aware of previous CS work establishing such a relation. The practical consequence is that given a dictionary of n-character words {w k } having mutually large Hamming distances (i.e., a good error correcting code), this construction deterministically specifies a set of nearly-orthogonal petal-vectors {|w k } (i.e., good vectors for sparse random sampling). Conversely, the general problem of constructing a deterministic set of nearly-orthogonal petal-vectors is seen to be precisely as difficult as deterministically constructing a good error correcting code. 4.6.5. The n and p dimensions of deterministic sampling matrices The column dimension p of the petal-vector sampling matrices thus constructed is given in Table 1 for Hamming distances 1-4 as a function of the row dimension n, the number of polytope vertices n ver , and the dimensionality of the polytope space dim V 0 (e.g., for our tetrahedral construction n ver = 4 and dim V 0 = 2). We see that for fixed row dimension n, larger Hamming distances are associated with smaller column dimension p, as is intuitively reasonable: the more stringent the pairwise orthogonality constraint, the smaller the maximal dictionary of sampling vectors that meet this constraint.
The construction has a further dimensional constraint as follows: it is straightforward for values of n that are powers of two (because the tetrahedral construction can be used), more complicated when n has factors other than 2 (because larger-dimension polytope vertices must be specified), and infeasible when n is large and prime. These constraints are reminiscent of similar constraints that act upon the fast Fourier transform, and arise for basically the same number-theoretic reason. identity code Hamming code and parity check b " « Summary of symbol definitions p = number of sampling matrix columns wC = length of codewords, thus n = (dim V0) w C n = number of sampling matrix rows wD = length of data words, thus p = (nver) w D nver = number of polytope vertices m = (integer) index of Hamming code dim V0 = dimensionality of the polytope space a Parity characters are calculated mod nver rather than the (more common in the literature) mod 2. b Hamming-and-parity codes are sometimes called SECDED codes (single-error correct, double-error detect). c The given functions p(n, nver, dim V0) are exact, and are valid for both real and complex vector spaces. Table 1: Recipes for deterministically constructing sampling matrices by the methods of Section 4.6.4. The primary design variables are taken to be the dimensionality of polytope space dim V 0 , the number of polytope vertices n ver , and the desired number of sampling matrix rows n. All expressions are exact, and the results apply to both real and complex vector spaces. The expressions are organized so as to make manifest that increased minimal Hamming distance is associated with decreased column-length p; this is the central design trade-off. 4.6.6. Petal-counting in GK geometry via coding theory These sampling theory results have a direct quantitative relation to the efflorescent GK geometry that we discussed in Sections 1.5.8 and 2.11. Specifically, we are now able to construct a petal-vector description of GK manifolds, and verify that they indeed have exponentially many petals.
We consider a spin-1 2 rank-1 GK manifold having n spin spins. The preceding tetrahedral GK construction deterministically generates a dictionary of petal-words {w k : k ∈ 1, 4 n spin } in one-to-one correspondence with petal-vector states {|w k : k ∈ 1, 4 n spin }. This dictionary of states of course exponentially over-complete, since its number of words is 4 n spin = 2 n spin dim H. Yet we also know that random pairs of petal-vectors in our tetrahedral dictionary are pairwise orthogonal to an excellent approximation, because their median Hamming distance is 3n spin /4, and consequently from (152) their median pairwise wedge product is As a concrete exercise in petal-counting, we consider a system of n spin = 16 spin-1 2 particles. The tetrahedral construction generates a dictionary of 4 16 = 2 32 petal-vectors for this system, each word of which labels a petal whose state-vector has a wedge separation of |w i ∧ w k | 2 ≥ 2/3 from the state-vector of all other petals. A subset of that petal dictionary having minimal Hamming distance 4 is specified by the SECDED code of Table 1. This SECDED subset has Hamming parameter m = 4, and hence m + 1 = 5 characters out of 16 in each word are devoted to error-correcting. The resulting (smaller) error-corrected dictionary has 4 16−5 = 4 11 = 2 22 petal-vectors, and the sampling matrix whose columns are  With only a bit of exaggeration, we can say that, if you formulate a practical problem as a convex optimization problem, then you have solved the original problem.
But this assertion must be regarded with caution when it comes to convex optimization over quantum state-spaces, because the matrices and vectors involved are of enormously larger dimension than is usually the case in convex optimization. Figueiredo, Nowak, and Wright [70] and also Cai and Lv [29] discuss this domain, and it is clear that Cai and Lv's conclusion "Clearly, there is much work ahead of us" applies especially to compressive quantum sensing, sampling, and simulation.
4.6.8. RIP properties of deterministic versus random sampling matrices It is clear from the preceding discussion that over-complete dictionaries of word-states {|w i } having the approximate orthogonality property w i |w j = (XX † ) ij δ ij are desirable both for simulation purposes and for sampling purposes. Stimulated by the work of Candes and Tao [37], an extensive and rapidly growing body of work characterizes such matrices in terms of the restricted isometry property (RIP). We now briefly discuss the RIP of tetrahedral sampling matrices, mainly following the notation and discussion of Baraniuk et al. [10]. We regard the word indices i and j in w i |w j as the rows index and column index of a Hermitian matrix. We specify a subset T of word indices, and we define the sparsity S of that subset to be S = #T . Then w i |w j T ≡ w i |w j : i, j ∈ T is an S × S Hermitian matrix, which we take to have minimal (maximal) eigenvalues λ min (λ max ). Then the isometry constant δ S of Candes and Tao is by definition Our word-state dictionary is said to have the restricted isometry property for order S iff δ S ∈ (0, 1) for all subsets T having sparsity S. Physically speaking, a dictionary of p wordstates having the RIP property in order S has the property that any set of S words is (approximately) mutually orthogonal. Testing for the RIP property is computationally inefficient, since (at present) no known algorithm is significantly faster than directly evaluating λ min and λ max for all p S distinct subsets T . Referring to Table 1, we see that a spin- 1 2 tetrahedral dictionary of three-letter words, one of which is a parity-check character, such that the minimal Hamming distance is two, yields a sampling matrix having p = 4 2 = 16 columns and with n = 2 3 = 8 rows. To calculated the RIP properties of this dictionary, the maximum sparsity we need to investigate is S = n = 8, for which p n = 16 8 = 12870 subsets must be evaluated, which is a feasible number. As summarized in Table 2, the tetrahedral construction yielded sampling matrices having the RIP property for sparsity S = 1, 2, 3, while for higher values of S the fraction of subsets having δ S ∈ (0, 1) dropped sharply.
For purposes of comparison, we computed also the median RIP properties of 8 × 16 random Gaussian matrices. We found for all values of sparsity, the RIP properties of Gaussian random matrices were inferior to those of the deterministic tetrahedral construction. We are not aware of any previous such random-versus-deterministic comparisons in the literature. Since is is known that the Gaussian random matrices are RIP in the large=p limit, we were surprised to find that their RIP properties are unimpressive for moderate values of p.
In preliminary studies of larger matrices, we found that known asymptotic expressions for the extremal singular values of Gaussian random sampling matrices-due to Marčenko and Pastur [134], Geman [87], and Silverstein [183], as summarized for CS purposes by Candès and Tao [33, see their Sec. III]-were empirically accurate for tetrahedral sampling matrices too, for all values of the row dimension n ≤ 256 and all values of the sampling parameter S ≤ n.
We emphasize however that although the average-case performance of these petal-vector sampling matrices is empirically comparable to Gaussian sampling matrices, their worstcase performance is presently unknown, and in particular such key parameters as their worst-case isometry constants are not known.
As Baraniuk et al. [10] note: "the question now before us is how can we construct matrices that satisfy the RIP for the largest possible range of S." It is clear that answering this question, in the context of the deterministic geometric construction given here, comprises a challenging problem in coding theory, packing theory, and spectral theory, involving sophisticated trade-offs among the competing goals of determinate construction, large (and adjustable) p/n ratio in the sampling matrix, and small isometry constants for all values of the sparsity parameter S ≤ n. 4.6.9. Why do CS principles work in QMOR simulations? Guided by the preceding analysis, we now try to appreciate more broadly why CS principles "work' in QMOR simulations by systematically noting mathematical parallels between the two disciplines. We will see that these parallels amount to an outline for extending the mathematical foundations of CS to provide foundations for QMOR-CS.
As our first parallel, we remark that what Candès and Tao call [37] compressible objects are ubiquitous in both the classical and quantum worlds. This ubiquity is not easily explained classically, and so almost always it is simply accepted as a fact of nature; for example almost all visual fields of interest to human beings are compressible images. In contrast, we have seen that the ubiquity of quantum compressible objects has a reasonably simple explanation: most real-world quantum systems are noisy, and noisy systems can be modeled as synoptic measurement processes that compress state trajectories; working through the mathematical details of this synoptic compression was of course our main concern in Section 3. From this quantum informatic point of view, it is a fundamental law of nature that any quantum system that has been in contact with a thermal reservoir (or equivalently, a measurement-and-control system) is a compressible object. The second parallel is the availability of what the CS field calls dictionaries [190] of the natural elements onto which both classical and quantum compressible objects are projected. For example, wavelet dictionaries are well-suited to image reconstruction. In the QMOR formalism of this article, the parallel quantum dictionary is (of course) the class of multilinear biholomorphic GK polynomials that define the Kählerian geometry of QMOR state-spaces (Section 2). This is not a linear dictionary of the type generally discussed in the CS literature, but rather is an algebraic generalization of such dictionaries. In the language of Donoho [61], open quantum systems exhibit a generalized transform sparsity whose working definition is the existence of high-fidelity projections onto GK manifolds.
The third parallel is the existence of robust, numerically efficient methods for projection and reconstruction. It is here that the mathematical challenges of aligning QMOR with CS are greatest. In our own research we have tried non-CS/non-QMOR optimization techniques-like regarding ψ 0 − ψ(c) = 0 as the definition of an algebraic variety, and decomposing it into a Groëbner basis-but in our hands these methods perform poorly. Turning this observation around, it is possible that the efficient methods of QMOR-CS might find application in the calculation of (specialized algebraic forms of) Groëbner bases.
Although a substantial body of literature exists [189] for minimizing functions that are convex along geodesic paths on Riemannian state-spaces-which generalizes the notion of convexity on Euclidean spaces-there does not seem to be any similar body of literature on the convexity properties of holomorphic functions on Kählerian state-spaces.
We have previously quoted Shing Tung Yau's remark [200, p. 21]: "While we see great accomplishments for Kähler manifolds with positive curvature, very little is known for Kähler manifolds [having] strongly negative curvature." By the preceding construction, we now appreciate that (negatively curved) GK manifolds have embedded within them lattices that display all the intricate mathematical structure of coding theory-so that it is not surprising that the geometric properties of these manifolds resists easy analysis. It seems that the ultrahigh-dimensional model selection of Cai and Lv [29] can be describedwith more-or-less equal mathematical justification-in terms of the differential geometry of ruled manifolds, or alternatively in terms of coding theory, or alternatively in terms of optimization theory.
There is also the as-yet unexplored practical issue of whether quantum optimization of l 2type functions over GK polynomials like (139) is more efficient, less efficient, or comparably efficient to CS-type optimization over petal-words of l 1 -type functions like (156a-c). This question is analogous to the long-standing issue in CS of whether interior-point methods are superior to edge-and-vertex polytope methods ... the answer after five decades of CS research being "yes, sometimes." Efficient numerical means for evaluating the Penrose pseudo-inverse of (143) are needed, as this inversion is the most computationally costly step of our sparse reconstruction codes as they are presently implemented. Preconditioned conjugate gradient techniques are one attractive possibility [51,94,96], because these techniques lend themselves well to the largescale parallel processing. The algebraic structure of the GK metric tensor creates additional algorithmic challenges and opportunities that (so far as the authors are aware) have not been addressed in the computing literature.
Finally, suites of test problems and open-source software tools have contributed greatly to the rapid development of CS theory and practice [190], and it would be valuable to have a similar suite of problems and tools for the simulation of open quantum systems.

Conclusions
As Terence Tao has remarked [187] The field of [partial differential equations] has proven to be the type of mathematics where progress generally starts in the concrete and then flows to the abstract, rather than vice versa.
The recipes in this article have resulted from a flow in the opposite direction, from abstract to concrete, in which abstract ideas from quantum information theory, algebraic geometry, quantum physics, and compressive sensing have found find concrete embodiment in practical recipes for large-scale quantum simulation.
Key abstract ideas from quantum information theory (QIT) Our recipes have adopted from quantum information theory the key idea that noise processes can be modeled as covert measurement processes. This leads naturally to the idea that quantum states that have been in contact with a thermal reservoir (or equivalently, a measurement and control process) are compressible objects.
Key ideas from algebraic geometry Our recipes have adopted from algebraic geometry the key idea that reduced-order quantum state-spaces can be described as geometric objects, using the language and methods of algebraic and differential geometry. In particular, quantum trajectories can be described in terms of drift and diffusion processes upon state-space manifolds, just as in classical modeling and simulation theory.
Key ideas from quantum physics theory Our recipes have adopted from theoretical quantum physics the key idea that quantum states have multiple unravelings, and that the efficiency of a calculation can be optimized by choosing an appropriate unravelling.
Key ideas from quantum physics experiments Our recipes have adopted from experimental quantum physics the key idea that mathematical ingredients of quantum simulation map one-to-one onto familiar physical systems such as measuring devices, and also the operational principle that the main deliverable of a quantum simulation is accurate prediction of the results of physical measurements.
Key ideas from compressive sensing, sampling, and simulation (CS) Our recipes have adopted from CS the idea that optimization problems involving compressible objects (like quantum states) can often be transformed into convex optimization problems. This can lead to both faster algorithms and improved physical insight.

Concrete applications of large-scale quantum simulation
By combining the preceding abstract ideas, the objective that began this article . . . to enable the reader to design and implement practical quantum simulations, guided by an appreciation of the geometric, informatic, and algebraic principles that govern simulation accuracy, robustness, and efficiency. has now been achieved in a preliminary sense, albeit there is much further work to be done. Now we consider some practical applications.

The goal of atomic-resolution biomicroscopy
The goal of atomic-resolution microscopy was the main motivation for developing the simulation algorithms described in this article. This goal was proposed as early as 1946 by Linus Pauling who envisioned "If it were possible to make visible the individual molecules of the serum proteins and other proteins of similar molecular weight, all the uncertainty which now exists regarding the shapes of these molecules would be dispelled" [151]. Later that same year, John von Neumann (possibly speaking as a reviewer of Pauling's proposal [116]) wrote a letter to Norbert Wiener [145] that embraced and extended Pauling's vision. The letter expressed a strikingly modern vision of atomic-level structural and systems biology: There is no telling what really advanced electron-microscopic techniques will do. . . . A "true" understanding of [viral-scale] organisms may be the first step forward and possibly the greatest step that may at all be required. I would, however, put on "true" understanding the most stringent interpretation possible: That is, understanding the organism in the exacting sense in which one may want to understand a detailed drawing of a machine, i.e. finding out where every individual nut and bolt is located.
It was not until 1959 that Richard Feynman-who spent a sabbatical year working as a biochemist [63]-issued his famous challenge: "Is there no way to make the electron microscope more powerful? . . . Make the microscope one hundred times more powerful, and many problems of biology would be made very much easier" [68].
Unfortunately, the problem of electron-beam radiation damage to fragile biological molecules proved intractable [110], and so the Pauling-von Neumann-Feynman challenge of achieving atomic-resolution biomicroscopy remained unanswered for several decades.
New ideas were needed, and three key ideas that emerged in ensuing decades were magnetic resonance imaging, nanotechnology, and quantum measurement theory. The early stages of development of each of these new fields was slow, not because fundamentally new concepts of mathematics or physics were required-the key concepts were reasonably familiar to Pauling, von Neumann, and Feynman's generation-but because each field sought to push familiar concepts to extreme limits. Magnetic resonance was a familiar concept; exploiting the tiny magnetic resonance signals for 3D imaging purposes was novel. Making devices smaller was a familiar concept; fabricating micron-scale and nanometer-scale devices was novel. Quantum measurement was a familiar concept; studying in detail the quantum evolution of small, continuously observed systems was novel.
These research fields are united in magnetic resonance force microscopy (MRFM), which was conceived explicitly as a means of meeting the Pauling-von Neumann-Feynman challenge of atomic-resolution biomicroscopy [177,178,179,181].

The acceleration of classical and quantum simulation capability
Simulation technologies, both classical and quantum, began an immense surge of progress during the Pauling-von Neumann-Feynman era, and this surge has continued to the present day. This is true especially in the classical domain, where simulation tools have become essential to system-level engineering [20,113]. In an era in which a new aircraft, a new processor chip, or a new drug can readily incur system development costs in excess of one billion dollars, engineering processes that maximize confidence, reliability, and economy have become a practical necessity.
Feynman, in his seminal 1982 article Simulating physics with computers [67], argued that generic quantum physics problems are exponentially hard to simulate on a classical computer. But Feynman's analysis was vague about what constitutes generic quantum physics. The viewpoint we have developed in this article is that any quantum system that has been in contact with a thermal reservoir is-in principle at least-a compressible object (in the language of CS theory) and thus is amenable to simulation with classical resources.
The present status of quantum simulation algorithms has, in our experience, striking parallels to the status of linear programming algorithms during 1947-2004 [56]. It was clear for many decades that linear programming methods were exceedingly useful for solving practical problems, but their mathematical foundations in convex set theory were established only very slowly (see Spielman and Teng's 2004 article [185] for a review and a reasonably definitive solution).
This slow-but-steady progress illustrates Dantzig's principle [5]: In brief, one's intuition in higher dimensional space is not worth a damn! Only now, almost forty years after the time the simplex method was first proposed, are people beginning to get some insight into why it works as well as it does. and also Feynman's words [66] I think the problem is not to find the best or most efficient method to proceed to a discovery, but to find any method at all. . . . [That is why] it is useful to have a wide range of physical viewpoints and mathematical expressions of the same theory.

The practical realities of quantum system engineering in MRFM
The practical experience of operating an MRFM devices is very much like operating a small satellite that is distant from the experimenter not in space, but in scale. In particular, the state-space of the MRFM device includes the quantum state-space of spins in the sample. The QMOR analysis methods of this article were conceived specifically to allow the efficient a priori modeling of this quantum state-space, by methods that have the potential to substantially extend present capabilities in modeling large-scale spin-systems [4]. To the extent that future progress in QMOR analysis allows this goal to be achieved, then the optimization of MRFM technology can proceed partially in silico, which will help retire technical risk, speed development, and build alliances among sponsors, researchers, enterprises, and customers.

Future roles for large-scale quantum simulation
In our view, the single most important role for quantum system engineering (QSE) will be to help sustain the exponentially cumulative technological progress that characterized the 20st century. In the context of computing this exponentiation is popularly known as "Moore's Law," but it is fair to say that similar exponentially cumulative progress is evident in fields such as nanotechnology, information theory, and (especially) biology. A recent theme issue of IBM Journal of Research and Development describes largescale simulation codes running on "Blue Gene" hardware that is approaching petaflop-scale computation speeds [112]. Both classical [65,72,73,124,163] and quantum [19,97,193] simulations are reviewed, and it is fair to say that the boundary between these two kinds of simulations is becoming indistinct, in particular when it comes to computing inter-atomic potentials that are both numerically efficient (classical) and accurate (quantum). This continues a sixty-year record of mutually supportive progress in hardware, software, and algorithm development [60].
From a geometric point of view, modern multi-processor computer architectures are exceptionally well-suited to the efficient computation of fundamental geometric objects such as polytopes, metrics, drift vectors, gradients, and diffusion tensors. These fundamental geometric objects are the raw building blocks-both conceptually and as software libraries-for broad classes of system simulations. In particular, the recipes of this article demonstrate that both classical and quantum systems can be simulated using a shared set of abstract mathematical ideas and concrete software tools that are well-matched to distributed computing architectures.
In summary, twenty-first century technologies seek to maximize the pace, coordination, and reliability of technology development, to create products that press against the quantum quantum and thermodynamic limits of device speed, sensitivity, size and power. For all such technologies, the quantum simulation recipes of this article promise to be useful. contract #W911NF-05-1-0403, and by the National Institutes of Health NIH/NIBIB under grant #EB02018. Prior support was provided by the ARO under contact #DAAD19-02-1-0344, by IBM under subcontract #A0550287 through ARO contract #W911NF-04-C-0052 from DARPA, and by the NSF/ENG/ECS under grant #0097544.