Geometry of variational methods: dynamics of closed quantum systems

We present a systematic geometric framework to study closed quantum systems based on suitably chosen variational families. For the purpose of (A) real time evolution, (B) excitation spectra, (C) spectral functions and (D) imaginary time evolution, we show how the geometric approach highlights the necessity to distinguish between two classes of manifolds: K\"ahler and non-K\"ahler. Traditional variational methods typically require the variational family to be a K\"ahler manifold, where multiplication by the imaginary unit preserves the tangent spaces. This covers the vast majority of cases studied in the literature. However, recently proposed classes of generalized Gaussian states make it necessary to also include the non-K\"ahler case, which has already been encountered occasionally. We illustrate our approach in detail with a range of concrete examples where the geometric structures of the considered manifolds are particularly relevant. These go from Gaussian states and group theoretic coherent states to generalized Gaussian states.


I. INTRODUCTION
Variational methods are of utmost importance in quantum physics. They have played a crucial role in the discovery and characterization of paradigmatic phenomena in many-body system, like Bose-Einstein condensation [1,2], superconductivity [3], superfluidity [4], the fractional quantum hall [5] or the Kondo effect [6]. They are the basis of Hartree-Fock methods [7], Bardeen-Cooper-Schrieffer theory [3], Gutzwiller [8] or Laughlin wavefunctions [9], the Gross-Pitaevsky equation [1,2], and the density matrix renormalization group [10], which are nowadays part of standard textbooks in quantum physics. Variational methods are particularly well suited to describe complex systems where exact or perturbative techniques cannot be applied. This is typically the case in many-body problems: on the one hand, the exponential growth of the Hilbert space dimension with the system size restricts exact computational methods to relatively small systems; on the other hand, as perturbations are generally extensive, they cannot be considered as small as the system size grows. Furthermore, variational methods are becoming especially relevant in recent times due to the continuous growth of computational power, as this enables to enlarge the number of variational parameters, for instance, to scale polynomially with the system size. Their power and scope can be further extended in combination with other methods, like Monte-Carlo, or even in the context of quantum computing.
A variational method parametrizes a family of states |ψ(x) or, in case of mixed states, ρ(x), in terms of so-called variational parameters x = (x 1 , . . . , x n ). The choice of the family of states is crucial as it has to encompass the physical behavior we want to describe, as well as to be amenable of efficient computations, circumventing the exponential growth in computational resources that appears in exact computations. A variety of variational principles can then be used, depending on the problem at hand. At thermal equilibrium, one can rely on the fact that the state minimizes the free energy, which reduces to the energy at zero temperature. In that case, for instance, one can just compute the expectation value E(x) of the Hamiltonian in the state |ψ(x) , and find the x 0 that minimizes that quantity, yielding a state, |ψ(x 0 ) that should resemble the real ground state of the sys-tem. For time-dependent problems, one can use Dirac's variational principle. There, one computes the action S[x(t),ẋ(t)] for the state |ψ(x(t)) and extracts a set of differential equations for x(t) requiring it to be stationary. Thus, the computational problem is reduced to solving this set of equations, which can usually be done even for very large systems. While the use of time-dependent variational methods is not so widespread as those for thermal equilibrium, the first have experienced a renewed interest thanks to the recent experimental progress in taming and studying the dynamics of many-body quantum systems in diverse setups. They include cold atoms in bulk or in optical lattices [11], trapped ions [12,13], boson-fermion mixtures [14], quantum impurity problems [15] and pump and probe experiments in condensed matter systems [16][17][18]. Recently, such methods have been used in the context of matrix product states to analyze a variety of phenomena, or with Gaussian states in the study of impurity problems, Holstein models, or Rydberg states in cold atomic systems.
Time-dependent variational methods can also be formulated in geometric terms. Here, the family of states is seen as a manifold in Hilbert space, and the differential equations for the variational parameters are derived by projecting the infinitesimal change of the state onto the tangent space of the manifold. This approach offers a very intuitive understanding of the variational methods through geometry. The translation between the different formulations is straightforward in the case of complex parametrizations: that is, where the x µ are complex variables, in which case the corresponding manifold is, from the geometric point of view, usually referred to as a Kähler manifold. If this is not the case, the geometric formalism has the advantage of highlighting several subtleties that appear and that have to be treated with care.
The main aim of the present paper is two-fold. First, to give a complete formulation of the geometric variational principle in the more general terms, not restricting ourselves to the case of complex parametrizations: that is, when the x µ are taken to be real parameters 1 . For all the variational methods that will be introduced, we will provide a detailed analysis of the differences that emerge in the non-complex case, most importantly the existence of inequivalent time dependent variational principles. The motivation to address the case of real parametrization stems from the fact that, in some situations, one has to impose that some of the variational parameters are real since otherwise the variational problem becomes intractable. This occurs, for instance, when one deals with a family of the form |ψ(x) = U(x)|φ , where φ is some fiducial state and U(x) is unitary. By taking x complex, U(x) ceases to be unitary, and thus even the computation 1 Note that a complex parametrization (in terms of z µ ∈ C) can always be expressed in terms of a real parametrization just by replacing z µ = x µ 1 + ix µ 2 , with x µ 1,2 ∈ R.
of the normalization of the state may require unreasonable computational resources. Second, even though there exists a vast literature on geometrical methods [19][20][21][22][23], it is mostly addressed to mathematicians and it may be hard to practitioners to extract from it readily applicable methods. Here, we present a comprehensive, but at the same time rigorous illustration of geometric methods that is accessible to readers ranging from mathematical physicists to condensed matter physicists. For this, we first give a simple and compact formulation, and then present the mathematical subtleties together with simple examples and illustrations. We will address some of the issues which are most important when it comes to the practical application of these methods: the conservation of physical quantities, the computation of excitations above the ground state, and the evaluation of spectral functions as suggested by the geometrical approach. For each of them, we will provide a motivation and derivation from physical considerations and, where we find inequivalent feasible approaches, give a detailed discussion of the differences and subtleties. Moreover, we discuss how the geometrical method can be naturally extended to imaginary time evolution, providing us with a very practical tool for analyzing systems at zero temperature.
To illustrate our results and to connect to applications, we discuss representative examples of variational classes, for which the presented methods are suitable. In particular, we will recast the prominent families of bosonic and fermionic Gaussian states in a geometric language, which makes their variational properties transparent. We will further show how the geometric structures discussed in this paper emerge in a natural way in the context of Gilmore-Perelomov group theoretic coherent states [24,25], of which traditional coherent and Gaussian states are examples. Finally, we will discuss possible generalizations going beyond ansätze of this type.
The paper is structured as follows: In section II, we motivate our geometric approach and present its key ingredients without requiring any background in differential geometry. In section III, we give a pedagogical introduction to the differential geometry of Kähler manifolds and fix conventions for the following sections. In section IV, we define our formalism in geometric terms and discuss various subtleties for the most important applications ranging from (A) real time evolution, (B) excitation spectra, (C) spectral functions to (D) imaginary time evolution. In section V, we apply our formalism to the variational families of Gaussian states, group theoretic coherent states (Gilmore-Perelomov) and certain classes of non-Gaussian states. In section VI, we summarize and discuss our results.

II. VARIATIONAL PRINCIPLE AND ITS GEOMETRY
In this section we give a brief summary and motivation of the main content of the paper, namely a geometric approach to time-dependent variational methods. We begin with Dirac's time-dependent variational principle, which is widely used in physics, highlighting subtleties and lessknown issues that can arise in the general case. Then, we give a geometric interpretation of the corresponding equations of motion for the variational parameters in terms of projections of infinitesimal time-evolutions onto the tangent space of the variational manifold. We show how this perspective allows a clearer understanding of the more intricate details and is useful to maximize the results one can extract from the formalism of the timedependent variational principle.
As we will see, the greatest benefits of this geometrical perspective arise when one considers the most general variational families, i.e., classes of states parametrized by arbitrary real parameters. The fact that the variational parameters are in general real has to be correctly taken into account when projecting time evolution, as the tangent space becomes a real vector space, and thus is different from the complex Hilbert space, in which the problem is originally formulated. This has consequences not only for the approximate time evolution, but also for determining which quantities are conserved, for computing an approximate excitation spectrum and other measurable quantities, such as spectral functions. Furthermore, the geometric approach will enable us to generalize the time-dependent variational principle to imaginary-time evolution to find approximate ground states.
All this motivates the geometric approach to study time-dependent variational methods. While in this section, we give a brief and simple overview of these methods, focusing on its most important aspects, the next sections will provide the reader with a deeper, mathematically rigorous and comprehensive account of the geometric approach and its consequences.
We consider closed quantum systems with Hamiltonian H acting on some Hilbert space H. Here, we have a many-body quantum system in mind, where H is a tensor product of local Hilbert spaces, although we will not use this fact in the general description. We would like to find the evolution of an initial state, |ψ(0) , according to the real-time Schrödinger equation and also according to the imaginary-time evolution The latter will converge to a ground state ofĤ as long as |ψ(0) possesses a non-vanishing overlap with the corresponding ground subspace.

A. Variational families
We will study variational families of states described as |ψ(x) ∈ H, where x µ ∈ R N is a set of real parameters. The goal is to approximate the time evolution within this class of states, i.e., to find a set of differential equations for x µ (t) so that, provided the variational family accounts well for the physically relevant properties of the states, we can approximate the exact evolution by |ψ[x(t)] . In the case of imaginary-time evolution, the goal will be to find x 0 such that |ψ(x 0 ) minimizes the energy, i.e., the expectation value ofĤ, within the variational family.
In principle, one could restrict oneself to variational families that admit a complex parametrization, i.e., where |ψ(z) ∈ H is a holomorphic in z µ ∈ C M and thus independent of z * . As we will see, this leads to enormous simplifications, as in the geometric language we are dealing with so-called Kähler manifolds, which have very friendly properties. However, in general, we want to use real parametrizations, which cover the complex case (taking the real and imaginary part of z as independent real parameters), but apply to more general situations.
While in certain situations, it is easy to extend or map a real parametrization to a complex one, this is not always the case. This applies, in particular, to parametrizations of the form where |φ is a suitably chosen reference state and U(x) is a unitary operator that depends on x µ ∈ R N . Such parametrizations are often used to describe various many-body phenomena appearing in impurity models [26][27][28][29] electron-phonon interactions [30,31] and lattice gauge theory [32], and the fact that U(x) is unitary is crucial to compute physical properties efficiently. However, extending x µ analytically to complexify our parametrizations, will break unitarity of U and often make computations inefficient, thereby limiting the applicability of the variational class. We review such examples in section V, including bosonic and fermionic Gaussian states and certain non-Gaussian generalizations.
The following example shows an important issue about different possibilities for parametrizations: Example 1. For a single bosonic degree of freedom (with creation operatorâ † and annihilation operatorâ), we define normalized coherent states as where x = (Rez, Imz). This parametrization is complex but not holomorphic. We can define the extended family whose parametrization is holomorphic in (z 1 , z 2 ). The latter parametrization differs from the former as it allows the total phase and normalization of the state to vary freely.
Given a family with a generic parametrization |ψ(x) we can always include two other parameters, (κ, ϕ), to allow for a variation of normalization and complex phase, so that the new family is |ψ (x ) = e κ+iϕ |ψ(x) , (6) where now the total set of variational parameters is x = (κ, ϕ, x). While the global phase does not have a physical meaning on its own, if we want to study the evolution of a superposition of one or several variational states, or quantities like spectral functions, the phase will be relevant and thus should be included in the computation. This extension of the variational parameteres can always be done at little extra computational cost and the variational principle can be formulated most simply in terms of the extended variables x . For this reason, in the rest of this section (except for subsection II B 1 where we add some extra observations on this issue) we will assume that this extension has been done and we will drop the primes, for the sake of an easier notation.

B. Time-dependent variational principle
One can get Schrödinger's equation (1) from the action as the Euler-Lagrange equation ensuring stationarity of S. This immediately yields a variational principle for the real-time evolution, the so-called Dirac principle. For this, we compute the Euler-Lagrange equations 2 as where E(x) = ψ(x)|Ĥ|ψ(x) is the expectation value of H on the unnormalized state |ψ(x) , ω µν = 2 Im v µ |v ν and |v µ = ∂ µ |ψ(x) . We exploited the antisymmetry of ω, resulting from the antilinearity of the Hermitian inner product. Here and in the following, we use Einstein's convention of summing over repeated indices and we omit to indicate the explicit dependence on x of some quantities, such as ω. Furthermore, we refer to the time derivative d dt by a dot, and to the partial derivative with respect to x ν by ∂ ν .
If ω is not invertible, this means that the evolution equations for x µ are underdetermined. The reason may be an overparametrization, in which case one can simply drop some of the parameters. However, when we discuss the geometric approach, we will see there can be other reasons related to the fact that the parameters x µ are real, in which case one has to proceed in a different way.
Let us also remark that if we have a complex representation of the state, i.e., |ψ(z) is holomorphic in z ∈ C M , we can get the equations directly for z, namely 3 where (Ω −1 ) µν = −i v µ |v ν with |v µ = d dz µ |ψ and E(z, z * ) = ψ(z * )|Ĥ|ψ(z) . Notice that in this case, Ω is invertible unless there is some redundancy in the parametrization. This is a consequence of the fact that such variational families are, from the geometric point of view, what is known as a Kähler manifold (see definition Section III).
In what follows, we will see that many desirable properties are naturally satisfied when dealing with Kähler manifolds, while we also point out various subtleties that arise otherwise.

Dynamics of phase and normalization
Let us now briefly consider some more details related to the inclusion of the normalization and phase (κ, ϕ) as variational parameters. For this, we will temporarily reintroduce the distinction between |ψ(x) and |ψ (x ) as in equation (6). If we consider the Euler-Lagrange equations corresponding specifically to each of the parameters (κ, ϕ, x), we have ϕ = − E(x ) +ẋ µ Im ψ (x )|v µ e 2κ ψ(x)|ψ(x) (12) and equations for x that do not depend on ϕ and are proportional to e 2κ , so that one can replace the solution of (11) in those equations and solve them independently. If one is interested in the evolution of the phase, then one just has to plug the solutions for x in (12) and integrate that differential equation separately. It is important to note that using the Lagrangian L from (7) without having introduced the extra parameters (κ, ϕ) or, more precisely, without ensuring that both 3 We have L = i 2 (ż µ ψ|vµ −ż * µ vµ|ψ ) − E(z, z * ). The Euler-Lagrange equations d dt ∂L ∂ż * µ = ∂L ∂z * µ are therefore given by the expression d dt ∂L ∂ż * µ = − i 2 (ż ν vµ|vν +ż * ν ∂ * ν vµ|ψ ) and also ∂L ∂z * µ = i 2 (ż ν vµ|vν −ż * ν ∂ * µ vν |ψ ) − ∂ * µ E.
phase and normalisation can be freely varied, can lead to unexpected results. More specifically, it can produce equations of motion which leave some parameters undetermined or where the unwanted coupling between phases and physical degrees of freedom leads to artificial dynamics. Nonetheless, one can also equivalently derive the equations for the x directly from a Lagrangian formulation, without introducing the extra parameters (κ, ϕ). It is sufficient to use the alternative Lagrangian which is invariant under |ψ(x) → c(x) |ψ (up to a total derivative) and thus differs from (7). We discuss more in detail how these two definitions are related in Section IV A 3.

Conserved quantities
An important feature of the time-dependent variational principle is that energy expectation value is conserved if the Hamiltonian is time-independent. This can be readily seen because from (11) we know that the states remain normalized during the evolution, so for an initially normalized state, the energy E will always coincide with the function E, for which we find as Ω is antisymmetric. However, in general, other observablesÂ =Â † that commute with the Hamiltonian may not be conserved by the time-dependent variational principle. Indeed, for every variational family, one can find symmetry generatorŝ A with [Â,Ĥ] = 0 which will not be preserved. The question is if those quantities are conserved, that are relevant for describing the physics of the problem at hand. In the special case where we have a complex parametrization, we will show in Section IV A 2 what further conditionsÂ has to satisfy for it to be conserved. More specifically, it must fulfil a compatibility requirement with the chosen variational family. Importantly, we will also discuss how this simple picture is no longer true in the case that no complex parametrization is available.
If the observables of interest happen not to be conserved, it may be wise to consider an alternative variational family, but one can also enforce conservation by hand at the expense of effectively reducing the number of parameters. There are indeed several possibilities to enforce the conservation of observables other than the energy.
For instance, one may think of including timedependent Lagrange multipliers in the Lagrangian action to ensure that property [33]. However, this can only work in a restricted number of cases, as can be already seen if one wants to conserve just a single observableÂ. Denoting by A(x) = ψ(x)|Â|ψ(x) , and adding to the Lagrangian L the term λ(t)A, it is easy to see that both

E[x(t)] and A[x(t)] remain constant iḟ
for all times. The function λ(t) can be chosen such thatÄ(t) = 0 for all times, namely taking λ(t) = On top of that, one has to choose an initial state and a parametrization such that at the initial timeȦ(0) = 0. Furthermore, the denominator in the definition of λ(t) must not vanish and since the addition of a Lagrange multiplier modifies the Schrödinger equation, one has to compensate for that.
In particular, at the final time T , one has to apply the operator exp(i T 0 λ(t)dtÂ) to the final state, which may be difficult in practice. This severely limits the range of applicability of the Lagrange multiplier method.
Another possibility is to solve A(x) = A 0 for one of the variables, e.g. leading to x N = f (x 1 , . . . , x N −1 ), and choose a new reduced variational family with parametersx = (x 1 , . . . , x N −1 ) as |ψ(x) = |ψ(x, f (x)) . On this reduced family, A will have the constant value A 0 by construction. However, this requires to find the function f which may be difficult in practice. In Section IV A 2 we will discuss how, thanks to the geometric understanding, this condition can be easily enforced locally without having to explicitly solve for f . In the same section we will also discuss how to deal with the fact that reducing the variational family by an odd number of real degrees of freedom, as proposed here, will inevitably make ω degenerate and thus non-invertible.

Excitation spectra
A standard approach for computing the energy of elementary excitations is to linearize the equations of motion (9) around the approximate ground state x 0 to find The spectrum of K comes in conjugate imaginary pairs ±iω . The underlying idea is that if we slightly perturb the state within the variational manifold and solve the linearized equations of motion, we can approximate the excitation energies as the resulting oscillation frequencies ω of the normal mode perturbations around the approximate ground state.
We will see how our geometric perspective provides us with another possibility to compute the excitation spectrum. Both methods have advantages and drawbacks which we carefully explain in IV B.

Spectral functions
In the literature one can find several approaches [34][35][36] for estimating spectral functions by relying on a variational family. In section IV C, we argue that the approach that at the same time is most in line with the spirit of variational principles and better adapts to being used with generic ansätze consists in performing linear response theory directly on the variational manifold. Furthermore, this approach leads to a simple closed formula for the spectral function based only on the generator K µ ν of the linearized equations of motion introduced in (16). Let us decompose K µ ν in terms of eigenvectors as 4 where E ν (λ) refers to the dual basis of left eigenvectors, chosen such that E µ (λ) E µ (λ ) = δ λλ . Further, we use the normalization E µ (iω ) * ω µν E ν (iω ) = i sgn(ω ), where we apply complex conjugation component-wise. Then, the spectral function associated to a perturbationV is calculated as

C. Geometric approach
Let us now make the connection between the timedependent variational method reviewed above with a differential geometry description. The basic idea is to consider the states |ψ(x) as constituting a manifold M embedded in Hilbert space, and define a tangent space at each point. Then the evolution can be viewed as a projection on that tangent space after each infinitesimal time step. The main issue here is that, if our parametrization is real, the tangent space is not a complex vector space. Therefore, we cannot utilize projection operators in Hilbert space, but rather need to define them on the real tangent spaces. Before entering the general case, let us briefly analyze the one of complex parametrization.
where |v µ = ∂ µ |ψ(z) . Thus, it lies in the tangent space, which is spanned by the |v µ . The right hand side of the equation, however, does not necessarily do so, asĤ |ψ(z) will have components outside that span. If we evolve infinitesimally and we want to remain in the manifold, we will have to project the change of |ψ(z) onto the tangent space. In fact, in this way we get the optimal approximation to the real evolution within our manifold. In practice, this amounts to projecting the right hand side of (1) on that tangent space. This can be achieved by just taking the scalar product on both sides of the equation with |v µ which leads exactly to (10).
If we do not have a complex parametrization, this procedure needs to be modified. In the rest of this section, we will explain how this is done.

Tangent space and Kähler structures
The tangent space T ψ to the manifold M at its point |ψ is the space of all possible linear variations on the manifold around |ψ(x) . We can write them aṡ x µ ∂ µ |ψ(x) and thus the tangent space can be defined as the span of the tangent vectors |v µ = ∂ µ |ψ(x) . However, as our parameters x are taken to be real to maintain generality, this span should only allow real coefficients. The tangent space should therefore be understood as a real linear space embedded in the Hilbert space H. In particular, this implies that for |v ∈ T ψ , the direction i |v should be seen as linearly independent of |v and therefore may itself not belong to the tangent space.
Note that if, on the other hand, M has a complex holomorphic parametrization then both |v µ and i |v µ naturally belong to the tangent space as they correspond to ∂ ∂Rez µ |ψ and ∂ ∂Imz µ |ψ , respectively. In this case, T ψ is clearly a complex subspace of the Hilbert space.
From the Hilbert inner product we can derive the two real-valued bilinear forms g µν = 2 Re v µ |v ν and ω µν = 2 Im v µ |v ν . (21) We define their inverses respectively as G µν and Ω µν . As mentioned before, ω may not necessarily admit a regular inverse, in which case we can still define a meaningful pseudo-inverse as discussed in section III C. Given any Hilbert space vector |φ , we define its projection on the tangent space T ψ as the vector P ψ |φ ∈ T ψ that minimizes the distance from |φ in state norm. As we are not dealing with a complex linear space, this will not be given by the standard Hermitian projection operator in Hilbert space. Rather, it takes the form Finally, let us introduce J µ ν = −G µσ ω σν , which represents the projection of the imaginary unit, as seen from where we used (22) and (21). As highlighted previously, i |v ν may not lie in the tangent space, in which case the projection is non-trivial and we have J 2 = −1 in contrast to i 2 = −1. We will explain that J satisfying J 2 = −1 on every tangent space is equivalent to having a manifold that admits a complex holomorphic parametrization, in such case we will speak of a Kähler manifold. If, on the other hand, it is somewhere not satisfied, we speak of a non-Kähler manifold and in this case there exist tangent vectors |v µ , for which i |v µ will not belong to the tangent space. Moreover, the projection P ψ will not commute with multiplication by the imaginary unit. Example 2. Following example 1, normalized coherent states have tangent vectors For z = 0, i |v µ will not be a tangent vector, i.e., i |v µ / ∈ span R (|v 1 , |v 2 ). This changes if we allowed for a variation of phase and normalization (from the complex holomorphic parametrization), such that we had the additional basis vectors |v 3 = |Ψ(x) and |v 4 = i |Ψ(x) .
Note that in this section we have chosen to use a simplified notation to be able to convey the most fundamental aspects in a simple and clear way. In the detailed discussion, that can be found in section III, we will see that for a general parametrization |ψ(x) that does not necessarily contain (r, ϕ), we need to define the tangent space in a slightly different way by projecting out the directions |Ψ(x) that correspond to changing phase or normalization of the state. To avoid confusion between these different definitions we use the symbols ω, g, J, P ψ and |v µ to indicate the quantities introduced in this section, while later we will use ω, g, J , P ψ and |V µ .

Real time evolution
We already mentioned how the time dependent variational principle is equivalent to projecting infinitesimal time evolution steps onto the tangent space. In the general case of non-Kähler manifolds, there exist two inequivalent projections of Schrödinger's equation given by which are obviously equivalent on a complex vector space, as the two forms only differ by a factor of i. However, the defining property of a non-Kähler manifold is precisely that its tangent space is not a complex, but merely a real vector space and multiplication by i will not commute with the projection P ψ . The first choice of projecting Schrödinger's equation in (25) can be shown to be equivalent to the formulation in terms of a Lagrangian L introduced earlier. It consequently leads to the equations of motion (9). The second choice of (25), often referred to as the McLachlan variational principle, corresponds to minimizing the local error d dt |ψ − (−iĤ) |ψ made at every step of the approximation of the evolution and leads to the equationṡ In section IV, we will argue that in most cases the Lagrangian action principle presents the more desirable properties. In particular, it leads to simple equations of motion that only depend on the gradient ∂ µ E and whose dynamics necessarily preserve the energy itself. However, for some aspects, the McLachlan evolution still has some advantages, such as the conservation of observables that commute with the Hamiltonian and are compatible with the variational family, in the sense defined in Section IV A 2. We will further explain, how one can construct a restricted evolution that maintains the desirable properties of both projections in (25) for non-Kähler families, but at the expense of locally reducing the number of free parameters.
Finally, our geometric formalism provides a simple notation to understand and describe the methods reviewed so far.

D. Imaginary time evolution
So far, our discussion was purely focused on real-time dynamics. In the context of excitations and spectral functions, we referred to an approximate ground state |ψ 0 in our variational family, that minimizes the energy function E(x). While there are many numerical methods to finding minima, our geometric perspective leads to a natural approach based on approximating imaginary time evolution, which we defined in (2) for the full Hilbert space. We would like to approximate this evolution as it is known to converge to a true ground state of the Hamiltonian, provided one starts from a state with a nonvanishing overlap with such ground state. However, as this evolution does not derive from an action principle, one cannot naively generalise for it Dirac's time dependent variational principle. On the other hand, the tangent space projection can be straightforwardly applied to equation (2), leading, as we prove in Section IV D, to the time evolution where E(x) = ψ(x)|Ĥ|ψ(x) / ψ(x)|ψ(x) is the energy expectation value function. The evolution defined in this way always decreases the energy E of the state, as can be seen from [37] where we used that G is positive definite.
Indeed, the dynamics defined by (27) can be simply recognised as a gradient descent on the manifold with respect to the energy function and the natural notion of distance given by the metric g. Consequently, this evolution will converge to a (possibly only local) minimum of the energy. In conclusion, we recognize imaginary time evolution projected onto the variational manifold as a natural method to find the approximate ground the state |ψ 0 = |ψ(x 0 ) .

III. GEOMETRY OF VARIATIONAL FAMILIES
In this section, we review the mathematical structures of variational families. To be as general as possible, we only assume that such a family M consists of pure quantum states ψ(x) which are differentiable functions of some real parameters x µ . Consequently, our family M is a real differentiable manifold and the main task of this section is to reformulate quantum mechanics in terms of real vector spaces and real differential geometry.
First, we explain how a complex Hilbert space can be described as real vector space equipped with so called Kähler structures. Second, we describe the manifold of all pure quantum states as projective Hilbert space, which is a real differentiable manifold whose tangent spaces can be embedded as complex subspaces in Hilbert space and thereby inherit Kähler structures themselves. Third, we introduce general variational families as real submanifolds, whose tangent spaces may lose the Kähler property. Fourth, we study this potential violation and possible cures.

A. Hilbert space as Kähler space
Given a separable Hilbert space H with inner product ·|· , we can always describe vectors by a set of complex number ψ n with respect to a basis {|n }, i.e., |ψ = n ψ n |n .
We will see that the tangent space of a general variational manifold is a real subspace of Hilbert space. Given a set of vectors {|n }, we distinguish the real and complex span On a real vector space, |ψ = 0 and i |ψ are linearly independent vectors, because one cannot be expressed as linear combination with real coefficients of the other. A real basis {|V µ } of H has therefore twice as many elements as the complex basis {|n }, such as Given any real basis {|V µ } of vectors, we can express every vector |X as real linear combination where we use Einstein's summation convention 5 . A general real linear mapÂ : H → H will satisfŷ A(α |X ) = αÂ |X only for real α. If it also holds for complex α, we refer toÂ as complex-linear. The imaginary unit i becomes itself a linear map, which only commutes with complex-linear maps.
The Hermitian inner product ·|· can be decomposed into its real and imaginary parts given by with g µν = 2 N Re V µ |V ν , ω µν = 2 N Im V µ |V ν and N being a normalization which we will fix in (51). This gives rise to the following set of structures. Definition 1. A real vector space is called Kähler space if it is equipped with the following two bilinear forms • Metric 6 g µν being symmetric and positive-definite with inverse G µν , so that G µσ g σµ = δ µ ν , • Symplectic form ω µν being antisymmetric and non-degenerate 7 with Ω µν , so that Ω µσ ω σν = δ µ ν , and such that the linear map J µ ν := −G µσ ω σν is a The last condition is also called compatibility between g and ω. We refer to (g, ω, J ) as Kähler structures.
Clearly, our g and ω are a metric and symplectic form, respectively. Furthermore, we will see that they are indeed compatible and define a complex structure J . For this, it is useful to introduce the real dual vectors Re X| and Im X| that act on a vector |Y via as one may expect. The identity 1 = n |n n| is then Similarly, the matrix representation of an operatorÂ is 5 We will be careful to only write equations with indices that are truly independent of the choice of basis, such that the symbol X µ may very well stand for the vector |X itself. This notation is known as abstract index notation (see appendix A 2). 6 Here, "metric" refers to a metric tensor, i.e., an inner product on a vector space. It should not be confused with the notion of metric spaces in analysis and topology. 7 A bilinear form bµν is called non-degenerate, if it is invertible.
For this, we can check det(b) = 0 in any basis of our choice. In particular, we compute the matrix representation of the imaginary unit i to be given by as anticipated in our definition. From i 2 = −1, we conclude that the so defined J indeed satisfies J 2 = −1 and is thus a complex structures. Therefore, g and ω as defined in (33) are compatible.
Example 3. A qubit is described by the Hilbert space H = C 2 with complex basis {|0 , |1 } and real basis With respect to this real basis g µν , ω µν and J µ ν are where 1 is the 2 × 2 identity matrix. We can represent a complex-linear mapÂ = n,m a nm |n m|, i.e., with [A, J] = 0, as the matrix where A = Re(a) and B = Im(a) in above basis.
In summary, every Hilbert space is a real Kähler space with metric, symplectic form and complex structure.

B. Projective Hilbert space
Multiplying a Hilbert space vector |ψ with a non-zero complex number does not change the quantum state it represents. Therefore, the manifold representing all physical states is given by the projective Hilbert space P(H), which we will define and analyze in this section. Variational families, which we will discuss in the following section, should then naturally be understood as submanifolds M of projective Hilbert space P(H).
The projective Hilbert space of H is given by the equivalence classes of non-zero Hilbert space vectors with respect to the equivalence relation Thus, a state ψ ∈ P(H) is a ray in Hilbert space consisting of all non-zero vectors that are related by multiplication with a non-zero complex number c.
The tangent space T ψ P(H) represents the space of changes δψ around an element ψ ⊂ P(H). Changing a representative |ψ in the direction of itself, i.e., |δψ ∝ |ψ , corresponds to changing |ψ by a complex factor and thus does not change the underlying state ψ. Two Hilbert space vectors |X , |X ∈ H therefore represent the same change |δψ of the state |ψ ∈ ψ, if they only differ by some α |ψ . We thus define tangent space as where we introduced the equivalence relation leading to a regular (not projective) vector space. We can pick a unique representative |X of the class [|δψ ] at the state |ψ by requiring ψ|X = 0. Viceversa, two vectors |X = |X both satisfying ψ|X = ψ|X = 0 belong to different equivalence classes. We thus identify T ψ P(H) with Given a general representative |δψ ∈ [|δψ ], we compute the unique representative mentioned above as |X = Q ψ |δψ with There is a further subtlety: representing a change δψ of a state ψ as vector |δψ will always be with respect to a representative |ψ . If we choose a different representative |ψ = c |ψ ∈ ψ, the same change δψ would be represented by a different Hilbert space vector |δψ = c |δψ . It therefore does not suffice to specify a Hilbert space vector |δψ , but we always need to say with respect to which representative |ψ it was chosen. This could be avoided when moving to density operators 8 .
The fact we can identify the tangent space at each point with a Hilbert space H ⊥ ψ enables us, given a local real basis {|V µ } at ψ, such that H ⊥ ψ = span R {|V µ }, to induce a canonical metric g µν , symplectic form ω µν and J µ ν onto the tangent space, which thus is a Kähler space, as discussed previously. We see at this point that on the tangent space T ψ P(H), it is convenient to choose N = ψ|ψ as normalization for the Kähler structures. The rescaled metric 1 2 g µν is well-known as the Fubini-Study metric [38,39], while the symplectic form gives projective Hilbert space a natural phase space structure.
Manifolds such as P(H), whose tangent spaces are equipped with differentiable Kähler structures, are called almost-Hermitian manifolds. In appendix C, we show that P(H) satisfies even stronger conditions, which make it a so called Kähler manifold.
Example 4. The projective Hilbert space of a Qubit is P(C 2 ) = S 2 , equivalent to the Bloch sphere. Using spherical coordinates x µ ≡ (θ, φ) and the complex Hilbert space basis {|0 , |1 }, we can parametrize the set of states as The elements of P(H) are the equivalence classes ψ(x) = c |ψ(x) c ∈ C, c = 0 . Consequently, the tangent space Using (33), we can compute the matrix representations We recognize g µν dx µ dx ν = 1 2 (dθ 2 + sin 2 (θ)dφ 2 ) to be the standard metric of a sphere with radius 1/ √ 2. Similarly, we recognize ω µν dx µ dx ν = 1 2 sin θdθ ∧dφ to be the standard volume form on this sphere. Finally, it is easy to verify that J 2 = −1 everywhere.
In summary, a given pure state can be represented by the equivalence class ψ ∈ P(H) of all states related by multiplication with a non-zero complex number. Similarly, a tangent vector [|X ] ∈ T ψ P(H) at a state ψ is initially defined as the affine space [|X ] of all vectors |X differing by a complex multiple of |ψ . A unique representative |X can be chosen requiring ψ|X = 0. This leads to the identification T ψ P(H) H ⊥ ψ , such that the Hilbert space inner product ·|· induces local Kähler structures onto T ψ P(H). 8 We can equivalently define projective Hilbert space as the set of pure density operators, i.e., Hermitian, positive operators ρ with Trρ = Trρ 2 = 1. The state ψ is then given by the density operator ρ ψ = |ψ ψ| ψ|ψ and its change δψ by δρ ψ = |δψ ψ|+|ψ δψ| ψ|ψ FIG. 2. Tangent vectors. We sketch the basis vectors |Vµ of tangent space T ψ M along some coordinate lines x µ .

C. General variational manifold
The most general variational family is a real differentiable submanifold M ⊂ P(H). At every point ψ ∈ M, we have the tangent space T ψ M of tangent vectors |X ψ . T ψ M can be embedded into Hilbert space by defining the local frame |V µ ψ ∈ H, such that |X = X µ |V µ , as before. Note, however, that in general the tangent space T ψ M = span R {|V µ } is only a real, but not necessarily a complex subspace of H. Thus, we will encounter families, for which |X is a tangent vector, but not i|X .
In practice, we often parametrize ψ(x) ∈ M by choosing a representative |ψ(x) ∈ H. This allows us to construct the local basis at the state |ψ(x) , where Q ψ was defined 9 in (46). To simplify notation, we will usually drop the reference to ψ(x) or x and only write |V µ , whenever it is clear at which state we are. Similar to projective Hilbert space, we define restricted Kähler structures on tangent space T ψ M ⊂ T ψ P(H) as There are two important differences to the respective definition (33) in full Hilbert pace. First, with a slight abuse of notation, |V µ here does not span Hilbert space, but rather the typically much smaller tangent space. Second, we set N = ψ|ψ just like for P(H), such that 9 The projector Q ψ is important to ensure that |Vµ can be identified with an element of H ⊥ ψ T ψ P(H) as discussed in Section III B, i.e., ψ|Vµ = 0. For derivations, it can be useful to go into a local coordinate system of x µ , in which |Vµ = ∂µ |ψ , i.e., the action of Q ψ can be ignored. This can always be achieved locally at a point and any invariant expressions derived this way, will be valid in any coordinate system. This has the important consequence that the restricted Kähler structures are invariant under the change of representative |ψ of the physical state. Namely, under the transformation |ψ → c |ψ with |V µ → c |V µ , our Kähler structures will not change. This ensures that equations involving restricted Kähler structures are manifestly defined on projective Hilbert space and thus independent of the representative |ψ(x) ∈ H, we use to represent the abstract state ψ(x) ∈ M ⊂ P(H).
We have T ψ M ⊂ H and for |X , |Y ∈ H, we have the real inner product Re X|Y on H inducing the norm |X = Re X|X = X|X , which is nothing more than the regular Hilbert space norm. We then define the orthogonal projector P ψ from H onto T ψ M with respect to Re ·|· , i.e., for each vector |X ∈ H we find the vector P ψ |X ∈ T ψ M minimizing the norm |X − P ψ |X .
We can write this orthogonal projector in two ways: such that we have P ψ = |V µ P µ ψ . The difference lies in the co-domain: While P ψ : H → H maps back onto Hilbert space, e.g., to compute P 2 ψ = P ψ , we have that P µ ψ : H → T ψ M is a map from Hilbert space into tangent space. Due to T ψ M ⊂ T ψ P(H), we have which follows from Q ψ |V µ = |V µ and Q † ψ = Q ψ . In contrast to Q ψ , the projector P ψ is in general not Hermitian.
Provided that there are no redundancies or gauge directions (only changing phase or normalization) in our choice of parameters, g µν will still be positive-definite and invertible with inverse G µν . We find that is the projection of the multiplication by the imaginary unit (as real-linear map) onto T ψ M. It will not square to minus identity if multiplication by i in full Hilbert space does not preserve the tangent space. If g µν is not invertible, it means that there exists a set of coefficients X µ such that X µ g µν X ν = 0, that is X µ |V µ = 0 and therefore X µ |V µ = 0. In other words, not all vectors |V µ are linearly independent and thus also not all parameters are independent. If this is the case, it is not a real problem as the formalism introduced can still be used with little modifications. More precisely, the projectors (53), as well as all other objects we will introduce, are meaningfully defined if we indicate with G µν the Moore-Penrose pseudo-inverse of g µν , i.e., we invert g µν only on the orthogonal complement to its kernel (orthogonal with respect of the flat metric δ µν in our coordinates 10 ). Indeed, all directions in the kernel correspond to a vanishing vector in the tangent space and therefore do not matter. In this case, also Ω µν , should be defined as the inverse of ω µν on the orthogonal complement to the kernel of g µν . 11 However, it is still possible that ω and J are not invertible even on this reduced subspace.
In this case, in order to define Ω one has to reduce oneself to working on an even smaller subspace, that is one that does not contain the kernel of ω and J . Here, however, the way in which we reduce these extra dimensions is not equivalent, as these directions are not anymore just redundant gauge choices of our parametrization. The reduction here effectively corresponds to working on a physically smaller manifold, as we will explain better in the next section. For what follows we will always suppose that Ω is defined by inverting ω on the tangent subspace orthogonal, with respect to the metric g µν , to the kernel of J . That is, Ω is the Moore-Penrose pseudo-inverse of ω with respect to g, i.e., the pseudo-inverse is evaluated in an orthonormal basis.
In conclusion, we see that we are able to define the restricted structures (g, ω, J ) which, however, do not necessarily satisfy the Kähler property. This is due to the fact that the tangent space, as we have pointed out, is a real, but not necessarily complex subspace of H. Note that these objects are locally defined in each tangent space T ψ M for ψ ∈ M.
Example 5. For the Hilbert space H = (C 2 ) ⊗2 of two Qubits, we can choose the variational manifold M of symmetric product states represented by is a single Qubit state as parametrized in (47). The tangent space is spanned by where |V µ was defined in (48). With this, we find and ω µν ≡ 0 sin θ − sin θ 0 (58) 11 Note that the kernel of ωµν itself does not necessarily correspond to redundant directions of the parametrization as X µ ωµν = 0 does not imply X µ |Vµ = 0.
leading to J 2 = −1 everywhere. We therefore conclude that the tangent space T ψ M satisfies the Kähler property everywhere.
Example 6. For the single Qubit Hilbert space H = C 2 , we can choose the equator of the Bloch sphere as our variational manifold M. This amounts to fixing θ = π/2 in (47) leading to the representatives with a single variational parameter φ. We have the single tangent vector |V = |V 1 as defined in (48). From the inner product V |V = 1 4 , we find g = 1 2 and ω = 0 implying J = 0. Consequently and not surprising due to the odd dimension, the tangent spaces of our variational manifold M are not Kähler spaces. Moreover, neither ω nor J are invertible.
In summary, we introduced general variational manifolds as real differentiable submanifolds M of projective Hilbert space P(H). By embedding the tangent spaces T ψ M into Hilbert space, the Hilbert space inner product defines restricted Kähler structures on the tangent spaces, whose properties we will explore next.

D. Kähler and non-Kähler manifolds
We categorize variational manifolds depending on whether their tangent spaces are Kähler spaces or not. We will see in the following sections that this distinction has some important consequences for the application of variational methods on the given family. • Non-Kähler, if it is not Kähler. If ω is degenerate, we define Ω as the pseudo-inverse.
This classification refers to the manifold as a whole.
Many well-known variational families, such as Gaussian states [41], coherent states [24,25,42], matrix product states [43] and projected entangled pair states [44], are Kähler. On the other hand, one naturally encounters non-Kähler manifolds when one parametrizes states through a family of general unitaries U(x) applied to a reference state |φ , i.e., For example, this issue arises for the classes of generalized Gaussian states introduced in [30], for the Multiscale Entanglement Renormalisation Ansatz states [45] or if one applies Gaussian transformations U(x) to general non-Gaussian states. 12 A general manifold M, whose tangent spaces are equipped with compatible Kähler structures, is known as an almost Hermitian manifold. However, if an almost Hermitian manifold is the submanifold of a Kähler manifold (as defined in appendix C), then it is also a Kähler manifold itself. Thus, due to the fact that P(H) is a Kähler manifold, all almost Hermitian submanifolds M ⊂ P(H) are also Kähler manifolds, which is why we use the term Kähler in this context. Example 8. We already reviewed examples for these three cases in the previous section. More precisely, example 5 is Kähler, example 6 is non-Kähler with degenerate ω and example 7 is non-Kähler with non-degenerate ω.
Given a submanifold M ⊂ P(H), we can use the embedding in the manifold P(H) to constrain the form that the restricted complex structure J can take on M.
Proposition 1. On a tangent space T ψ M ⊂ H of a submanifold M ⊂ P(H) we can always find an orthonormal basis {|V µ }, such that g µν ≡ 1 and the restricted complex structure is represented by the block matrix with 0 < c i < 1. This standard form induces the decomposition of T ψ M into the three orthogonal parts where T ψ M is the largest Kähler subspace and T ψ M is the largest space on which J and ω are invertible.
Proof. We present a proof in appendix B.
Proposition 1 is also known in the context of classifying real subspaces of complex Hilbert spaces. Interestingly, it is linked to the entanglement structure of fermionic Gaussian states, as made explicit in [46].
The manifold M is Kähler if there is only the first block in (66) everywhere. The symplectic form ω is nondegenerate if we only have the first two diagonal blocks. The next corollary provides some futher intuition for the non-Kähler case.
Proposition 2. The Kähler property is equivalent to requiring that T ψ M is not just a real, but also a complex subspace, i.e., for all |X ∈ T ψ M, we also have i|X ∈ T ψ M. Therefore, the multiplication by i commutes with the projector P ψ , i.e., P ψ i = iP ψ and P ψ is complex-linear.
Proof. We present a proof in appendix B.
If a manifold admits a complex holomorphic parametrization, i.e., a parametrization that depends on the complex parameters z µ , but not on z * µ , then the manifold will be Kähler. Indeed, taking Rez µ and Imz µ as real parameters gives the tangent vectors It is actually also possible to show that, viceversa, a Kähler manifold is also a complex manifold, that is it admits, at least locally, a complex holomorphic parametrization.
As mentioned before, in order to define the inverse of ω it is necessary to restrict ourselves to work only in a subspace of T ψ M. We now see that the definition we gave previously of always defining Ω as the pseudoinverse with respect to g coincides with always choosing to consider only the tangent directions in In order to apply variational methods as explained in the following sections, it may be necessary to at least locally restore the Kähler property. We can achieve this by locally further restricting ourselves to Using the bases {|V µ } and {|V µ }, we can define the restricted Kähler structures (g, ω, J ), which are compatible, and (g, ω, J ), where ω and J are non-degenerate. Our assumption on the definition of Ω can be understood as taking Ω to be zero on the subspace D ψ M, where ω is not invertible, and equal to the inverse of ω on T ψ M. Note that this definition is only possible because the tangent space is also equipped with a metric g, which makes the orthogonal decomposition In summary, a general variational family M ⊂ P(H) is not necessarily a Kähler manifold. We can check locally, if the restricted Kähler structures fail to satisfy the Kähler condition. If this happens, we can always choose local subspaces on which the restricted Kähler structures satisfy the Kähler properties or are at least invertible. Defining Ω as the pseudo-inverse with respect to g is equivalent to inverting ω only on T ψ M. In what follows, we therefore do not need to distinguish between the non-Kähler cases with degenerate or non-degenerate structures, as we will always be able to apply the same variational techniques based on Ω.

E. Observables and Poisson bracket
Any Hermitian operatorÂ defines a real-valued function Â on the manifold M and in fact on the whole projective Hilbert space. The function is given by the expectation value It is invariant under rescalings of |ψ by complex factors and is thus a well-defined map on projective Hilbert space P(H). We will use the notation Â and A(x) interchangeably. For the function deriving from the Hamiltonian operatorĤ, we use the symbol E = Ĥ and call it the energy. Given a Hermitian operatorÂ and the representative |ψ(x) , we have the important relation which is invariant under the change of representative |ψ → c |ψ and |V µ → c |V µ . It follows from where we used product rule and (50).
The following definition will play an important role in the context of Poisson brackets and conserved quantities. Every operatorÂ defines a vector field given by Q ψÂ |ψ . If this vector field is tangent to M for all ψ ∈ M, the following definition applies.
The symplectic structure of the manifold naturally induces a Poisson bracket on the space of differentiable functions, which is given by In some special cases this can be related to the commutator of the related operators.
Proposition 3. Given two Hermitian operatorsÂ and B of which one preserves the Kähler manifold M, i.e., the Poisson bracket is related to the commutator via Proof. We compute where we used (74) and J = −J −1 = Ωg for a Kähler manifold.
For M = P(H), above conditions are clearly met for any Hermitian operatorsÂ andB. For a general Kähler submanifold M ⊂ P(H), however, the validity of (78) depends on the choice of operators considered. On a submanifold which is not Kähler the statement is in general no longer valid.

IV. VARIATIONAL METHODS
Having introduced the mathematical background in the previous section, we can now study how variational methods allow us to describe closed quantum systems approximately. Given a system H governed by a Hamilto-nianĤ, we assume that a choice of a variational manifold M ⊂ P(H), as defined in the previous section, has been made and show how to (A) describe real time dynamics, (B) approximate excitation energies, (C) compute spectral functions, (D) search for approximate ground states. In doing so, we will emphasize the differences that arise between the cases where the chosen variational manifold is or is not of the Kähler type.
In this section we will consider the variational manifold M as a submanifold of projective Hilbert space, in the sense described in Section III. All information about the global phase and normalization of states is therefore considered as pure gauge and projected out of all tangent space quantities. For this reason, the objects we use in this section, which have been introduced in Section III, are slightly different from the ones defined in Section II, where we assumed for simplicity that global phase and normalization were included as independent parameters of the manifold. To avoid confusion we use here the symbols ω, g, J , P and |V µ , while in Section II we have used ω, g, J, P and |v µ .
The advantage of this formalism is that it allows to describe the most general parametrizations, including ones where global phase and normalization are fixed or vary either independently or as functions of the other parameters.

A. Real time evolution
For what concerns real time evolution, we would like to approximate the Schrödinger equation (1) on our variational manifold M. There are different principles, used extensively in the literature, according to which this ap- proximation can be performed. We will see that only in the case of Kähler manifolds they are all equivalent.

Variational principles
Following the literature, we can define the following variational principles for |ψ := |ψ(t) .
Lagrangian action principle. The most commonly used variational principle relies on the Lagrangian action, already introduced in (13), whose stationary solution satisfies for all times and all allowed variations |δψ(t) with Q ψ |δψ = |δψ − ψ|δψ ψ|ψ |ψ from (46). This is equivalent to Schrödinger's equation on projective Hilbert space 13 . On a variational manifold M ⊂ P(H), i.e., where we require Q ψ |δψ(t) ∈ T ψ(t) M in (82), we instead have This gives rise to the equations of motion (9) anticipated in Section II, which we derive in Proposition 4. For a time-independent Hamiltonian, they always preserve the energy expectation value.
McLachlan minimal error principle. Alternatively, we can try to minimize the error between the approximate trajectory and the true solution. As we do not 13 The fact that the projector Q ψ onto projective tangent space H ⊥ ψ appears, shows that the resulting dynamics is defined on projective Hilbert space, while global phase and normalization are left undetermined. We will explain how to recover the dynamics of phase and normalization in Section IV A 3. II. Action principles. We review the different action principles and how they relate to the respective manifolds.

Lagrangian
McLachlan Dirac-Frenkel know the latter, we cannot compute the total error, but at least we can quantify the local error in state norm due to imposing that d dt |ψ(x) represents a variation tangent to the manifold, i.e., Q ψ d dt |ψ(x) ∈ T ψ M. It is minimized by the projection This gives rise to the equations of motion (26) anticipated in Section II, which we derive in Proposition 5.
The resulting equations of motion only agree with the Lagrangian action if M is a Kähler manifold. Otherwise, they may not preserve the energy expectation value. Dirac-Frenkel variational principle. Another variational principle requires for all allowed variations |δψ(t) . It is easy to see that the real and imaginary parts of (86) are equivalent to (83) and (85) respectively. Therefore, this principle is welldefined (and equivalent to the other two) only in the cases in which they are equivalent between themselves, that is, as we will see, if and only if M is a Kähler manifold. Otherwise, the resulting equations will be overdetermined. Expressing equations (83) and (85) in coordinates leads to flow equations for the manifold parameters x(t). We can then define a real time evolution vector field X µ everywhere on M, such that Integrating such equations allows us to define the flow map Φ t that maps an initial set of coordinates x µ (0) to the values x µ (t) that they assume after evolving for a time t.
In the case of the Lagrangian action principle, the vector field X takes the following form.
where E(x) is the energy function, defined in the context of equation (72). Such evolution always conserves the energy expectation value.
Proof. From the definition (50), we have We substitute this in (83) and then expand the projectors using the relations (53), (55) and P ψ i |ψ = 0 to obtain We further simplify the expression by using (74) and (55). This leads to To obtain the variation of the energy expectation value E we compute directly where we used the antisymmetry of Ω µν .
The most important lesson of (88) is that projected time evolution on a Kähler manifold is equivalent to Hamiltonian evolution with respect to energy function E(x). As was pointed out in [19], already the time evolution in full Hilbert space, i.e., M = P(H), follows the classical Hamilton equations if we use the natural symplectic form Ω µν . Let us point out that the sign in equation (88) depends on the convention chosen for the symplectic form, which in classical mechanics differs from the one adopted here. One further consequence of equation (88) is that the real time evolution vector field X (x) vanishes in stationary points of the energy, that is points x 0 such that ∂ µ E(x 0 ) = 0. These points will therefore also be stationary points of the evolution governed by X .
Let us here recall that, if ω µν is not invertible, Ω µν refers to the pseudo-inverse, as discussed in sections III C and III D. This convention means that in practice the Lagrangian evolution will always take place in the submanifold of M on which ω is invertible. There may be pathological cases where ω vanishes on the whole tangent space and therefore the Lagrangian principle does not lead to any evolution.
In the case of the McLachlan minimal error principle, the evolution equations take the following form which cannot be simplified further. It is also in general not true that this evolution conserves the energy or that has a stationary point in energy minima.  (85) is Proof. By substituting (50) in (85) we havė from which the proposition follows by expanding the projector according to (53).
To perform real time evolution in practice, either based on (88) for Lagrangian evolution or based on (93) for McLachlan evolution, we will typically employ a numerical integration scheme [47,48] to evolve individual steps.
Let us now relate the different variational principles. Proof. To prove the statement, it is sufficient to see that equations (83) and (85) can be written simply as applying the tangent space projector P ψ to two different forms of the Schrödinger equation, i.e.,

Lagrangian:
These two forms only differ by a factor of i. However, as we discussed in proposition 2, one equivalent way to define the Kähler property of our manifold is that multiplication by i commutes with the projector P ψ . Therefore, if the manifold is Kähler, an imaginary unit can be factored out of equations (95) and (96) making them coincide. If, on the other hand, the manifold is non-Kähler, this operation is forbidden and they are in general not equivalent.
As discussed in Section III D, if the chosen manifold does not respect the Kähler condition, we always have the choice to locally restrict ourselves to consider only a subset of tangent directions with respect to which the manifold is again Kähler, i.e., T ψ M. Then both principles will again give the same equation of motion, which will have the same form as (88) where we just replace Ω µν with Ω µν , which will conserve the energy and have stationary points in the minima of the energy. We will refer to this procedure as Kählerization.
We can compute explicitly how the vector fields of the Lagrangian and McLachlan variational principles differ. For this, we only consider the subspaces, in which the complex structure fails to be Kähler, i.e., where we have as discussed in (66). On the enlarged tangent space including all vectors i |V µ , the enlarged complex structurě associated to the Lagrangian and the McLachlan principle, respectively. We see explicitly that they agree for c i = 1, but also when α i = β i = 0.
Example 9. We consider the variational family from example 7 for a system with two bosonic degrees of freedom. We choose the Hamiltonian where 1 and 2 are the excitation energies with ± = 1 ± 2 , while φ is a coupling constant, such thatĤ = 1n1 + 2n2 for φ = 0. Figure 4 shows the time evolution of the expectation valuesz α ≡ (q,p) for the two operatorŝ whereb was defined in (60). For r = 0, the variational family is Kähler and the two variational principles give rise to the same evolution. For r = 0, the two principles generally disagree. The explicit formulas can be efficiently derived using the framework of Gaussian states reviewed in section V A, where we reconsider the present scenario in example 20.
Kähler vs. non-Kähler. On a Kähler manifold all three variational principles are well-defined and equivalent. They all give the same energy conserving equations of motion (88). On a non-Kähler manifold, only the Lagrangian and McLachlan variational principles are welldefined, but they give in general inequivalent equations of motion given by (88) and (93). Only the Lagrangian ones will manifestly conserve the energy and have stationary points in the minima of the energy. In table II, we review advantages and drawbacks discussed in the following. While in most cases, the Lagrangian principle appears to be a natural choice, the McLachlan principle is often preferable if ω is highly degenerate-in particular, if ω = 0 its pseudo-inverse is Ω = 0 and the evolution would vanish everywhere independent ofĤ, such that the McLachlan principle appears to be the better choice.

Conserved quantities
Given the generatorÂ of a symmetry of the Hamilto-nianĤ, i.e., [Ĥ,Â] = 0, the expectation value A(t) = ψ t |Â |ψ t is necessarily preserved by unitary time evolution on the full Hilbert space We now consider if this continues to be true for projected time evolution on a manifold. For a time-independent HamiltonianĤ, we have seen that the energy expectation value E is always conserved by Lagrangian projected real time evolution. However, projected time evolution will not in general preserve expectation values of an operatorÂ with [Ĥ,Â] = 0. To guarantee this, one has to further require thatÂ preserves the manifold.
the expectation value A(x(t)), defined as in equation (72), is preserved under real time evolution projected according to the McLachlan variational principle. It is also true for Lagrangian variational principle, if the two principles agree, i.e., if the manifold is Kähler.
where in the first line we used relation (74), the definition of McLachlan evolution (85) and that P µ ψ Â |ψ = 0, in the second step we used that, thanks to the condition (104), the restricted bilinear form g in the first line coincides with the full Hilbert space one in the second line of (105).
This result only applies the McLachlan projected real time evolution, for which (94) holds. In the Lagrangian case, we would havė which is in general not equal to i ψ|[Ĥ,Â]|ψ ψ|ψ −1 on a non-Kähler manifold 14 and thus not necessarily zero.
We see here the main advantage of the Kählerization procedure described in the previous section. Indeed, through Kählerization we are able to define, even for general non-Kähler manifolds, a projected real time evolution that shares the desirable properties of both, the Lagrangian and the McLachlan projections, i.e., it is a symplectic, energy preserving evolution with stationary points in the energy minima and at the same time preserves the expectation value of symmetry generators satisfying (104). Note that Kählerization may spoil the conservation laws of observablesÂ, for which Q ψÂ |ψ does not lie in the Kähler subspace T ψ M, in which case we will need to enforce conservation by hand, discussed next.
For operatorsÂ where (104) is not satisfied, we have two options to correct for this: (a) Enlarge the variational manifold M, such that condition (104) is satisfied.
(b) Enforce conservation by hand, for which we modify the projected time evolution vector field X µ .
While option (a) is typically more desirable, it requires creativity to find a suitable extension of a given family M. Of course, ifÂ is an important physical observable that is relevant to the problem, a manifold that does not preserve it may not be a good choice to approximate the system's behavior. In practice, however, it may still be worthwhile to check the predictions of an approximated time-evolution adopting (b). This is done by adding a further projection of the real time evolution flow onto the subspace of the tangent space orthogonal (with respect to g) to the direction P µ ψÂ |ψ = G µν ∂ ν A. This is equivalent to restricting ourselves to the submanifold where A 0 is the initial value ψ(0)|Â|ψ(0) ψ(0)|ψ(0) −1 . Note that this modified evolution may spoil other conservation laws (e.g., energy) that were previously intact.
To preserve several quantitiesÂ I , we can project onto the subspace orthogonal to the span of X I = P ψÂI |ψ . If we also want to preserve the Kähler property, we should choose X I = (P ψÂI |ψ , iP ψÂI |ψ ). We can then define g IJ = X µ I g µν X ν J to define the projector where G IJ is the inverse (or pseudo inverse, if not all vectors X I are linearly independent) of g IJ . The modified Lagrangian evolution vector field X µ is while for the McLachlan evolution, we find where X µ represents the unmodified evolution vector field in the McLachlan case. It will conserve all expectation values A I (t) by construction. In the Lagrangian case also the energy will continue to be conserved by construction, which would need to be enforced by hand for the McLachlan case. Kähler vs. non-Kähler. On a non-Kähler manifold, where we have two inequivalent definitions of the evolution, only the one coming from the McLachlan principle preserves the expectation value of symmetry generators satisfying (104). Thus a key reason to Kählerize a non-Kähler manifold is to conserve these expectation values also in the Lagrangian evolution.

Dynamics of global phase
Up to now we have always considered our variational manifolds M as submanifolds of projective Hilbert space and thus the tangent space T ψ M as a subspace of H ⊥ ψ . This means all states are only defined up to a complex factor. In practice, our family ψ(x) ∈ P(H) will be described by a choice |ψ(x) ∈H , i.e., for every set of parameters x µ , we will have a Hilbert space vector |ψ(x) representing the quantum state ψ(x) ∈ P(H).
If the parametrization x µ happens to contain the global phase or normalization of the state as an independent parameter, we are overparametrizing our family and the evolution equations (88) or (93) will keep the evolution of some parameters undetermined leading to some gauge redundancy. This is due to the fact that normalization and phase do not change the quantum state |ψ(x) represents and our equations of motion only determine the physical evolution of the quantum state and not of its Hilbert space representative.
We can include the time evolution of the global phase and the state normalization by extending our parametrization by defining where κ and ϕ are two additional real parameters. If phase or normalization were already contained in x µ this will lead to an overparametrization, but we have already explained how to take care of this in Section III C. Then, on top of the real time evolution equations (83) or (85), we obtain equations for these extra parameters by projecting Schrödinger's equation on the corresponding tangent space directions, i.e., |V κ = |Ψ and |V ϕ = i |Ψ , to find the two equations Equivalently, as anticipated in Section II, we can use the Lagrangian action principle to find the same equations by extremizing the alternative action for the full set of parameters (x, κ, ϕ) rather than S from (81) for only x.
In both cases, the time evolution of x µ (t) is unchanged, but we find the additional equationṡ relating the evolution of phase and normalization with |ψ(x(t)) . Interestingly, the time evolution of κ will ensure that |Ψ(x, κ, ϕ) does not change normalization.
The procedure can be understood as follows. Global phase and normalization are conjugate parameters when considering Hilbert space as Kähler space, as can be seen from |V ϕ = i |V κ . When considering a variational manifold M ⊂ P(H), we have the following options: 1. When we are only interested in the time evolution of physical states ψ, we must project out the information about global phase and normalization using P µ ψ . Consequently, our evolution equations will not determine how to change global phase or normalization as this information is pure gauge. We followed this philosophy until the current section.
2. When we are also interested in the time evolution of global phase and normalization, we can always extend M to include both phase and normalization as independent parameters. Given a generic parametrization |ψ(x) , we can extend it to |Ψ(x, κ, ϕ) to ensure that it satisfies the Kähler property in the phase/normalization subspace. We can then find evolution equations for ϕ and κ. This is what we explained in the current subsection.
Finally, let us emphasize that using equations (112) or extremizing action (113) without first ensuring both phase and normalization are included as independent parameters may lead to unphysical results.

B. Excitation spectra
We would like to use a variational family M to approximate the excitation energies E i of some eigenstates |E i of the Hamiltonian. Typically, we are interested in low energy eigenstates, that is eigenstates close to the groundstate of the Hamiltonian. Suppose then that on M we are able to find an approximate ground state |ψ 0 , that is the state with energy ω 0 that represents the global energy minimum on M (we will describe a method for finding such state in Section IV D). Then there are two distinct approaches of deriving a spectrum: the projection of the Hamiltonian and the linearization of the equations of motion.

Projected Hamiltonian
Given a tangent space T ψ0 M at a state ψ 0 ∈ M, we can approximate the excitation spectrum of the Hamil-tonianĤ from its projection onto T ψ0 M.
Based on the two action principles (Lagrangian vs. McLachlan), we define two different projections given by On a Kähler manifold M, we will have R µ ν = −J µ σ H σ ν and [J , H] = [R, H] = 0. In this case, H represents a Hermitian operator on tangent space (which is complex sub Hilbert space) and R is anti-Hermitian. In this case, the eigenvalues of H are real and come in pairs (ω , ω ), while the ones of R come are purely imaginary and come in conjugate pairs (iω , −iω ). The two associated eigenvectors of R are related by multiplication of J and also span the respective eigenspace of H. On a non-Kähler manifold M, the relation between H and R as well as their respective spectra is non-trivial. The eigenvalues ω of H µ ν will still be real, while the ones of R µ ν continue to appear in conjugate pairs. The latter also implies that for an odd-dimensional manifold R µ ν must have a vanishing eigenvalue, which is a pure artefact of the projection.
The projected Hamiltonian H µ ν represents the full Hamiltonian restricted to the tangent space.
The Courant-Fischer-Weyl min-max principle states that the eigenvalues E ofĤ and the eigenvalues ω of H µ ν satisfy with N = dim R H and n = dim R T ψ0 M, where we assume that all eigenvalues are sorted and appear with their multiplicity. Therefore, every approximate eigenvalue ω bounds a corresponding true eigenvalue E i from above. How good this approximation is will highly depend on the choice of variational manifold. Note that the energy differences ω − ω 0 instead do not necessarily bound E − E 0 , because the ground state energy ω 0 might not be exact, i.e., ω 0 > E 0 . Furthermore, the eigenvalues ω are variational in the sense that if X µ is an eigenvector of H µ ν such that H µ ν X ν = ω X µ , then the corresponding Hilbert space vector |X = X µ |V µ is satisfies Kähler vs. non-Kähler. On a Kähler manifold, H µ ν and R µ ν are related via R = −J H and they will be the representations of a complex Hermitian and anti-Hermitian operators, respectively. Real eigenvalue pairs (ω , ω ) of H will be related to imaginary eigenvalue pairs (iω , −iω ) of R. On a non-Kähler manifold, the eigenvalues ω of H could all be different and unrelated to the ones R, which are still imaginary appearing in conjugate pairs.

Linearized equations of motion
A common alternative is to linearize the equations of motion around a fixed point x 0 such that X (x 0 ) = 0 Here, δx µ represents a small perturbation around the approximate ground state. The frequencies ω appearing in conjugate pairs ±iω in the spectrum of K µ ν represent the frequencies with which such perturbations oscillate around the ground state and thus provide an approximation to the excitation energies E − E 0 of the Hamiltonian.
As pointed out in Section IV A 1, the fixed point x 0 only coincides with the approximate ground state ψ 0 if the real time evolution is defined in terms of the Lagrangian action principle. We thus assume the equations of motion (88) based on Lagrangian action principle. In this case, we find where everything is evaluated at x 0 after taking the derivatives. We used that ∂ ρ E = 0 at the fixed point 15 . By construction, K µ ν is a symplectic generator, because it satisfies KΩ + ΩK = 0, which implies that M = e K preserves the symplectic form, i.e., M ΩM = Ω. Provided that ψ 0 is an energy minimum, the bilinear form h µν = ∂ ν ∂ µ E is positive definite. By Williamson's theorem [49], K µ ν is diagonalizable and the resulting eigenvalues appear in conjugate pairs (iω , −iω ).
From a geometric point of view, δx µ represents a tangent vector |δψ = δx µ |V µ living in the tangent space T ψ0 M at the approximate ground state |ψ 0 . The time evolution of a tangent vector at a fixed point |ψ 0 is described by the linearized evolution flow 16 K is the generator of the flow dΦ t leading to the important relation 15 Mathematically speaking, K µ ν is the derivative of the vector field X µ . It is well-known in differential geometry that such a derivative is coordinate dependent, unless one uses the so-called covariant derivative ∇ν X µ = ∂ν X µ + Γ µ νρ X ρ , where Γ µ νρ need to be chosen. However, in our case X µ vanishes at the fixed point, which implies that the covariant derivative at |ψ 0 is actually independent of the choice of Γ µ νρ and K µ ν is canonically defined. 16 Mathematically, the linearized flow dΦt is defined as the differential (also known as push-forward) of the flow map Φt defined after equation (87). In general it is a map from the tangent space In the special case of |ψ 0 being a fixed point of the time evolution, it reduces to a linear map from T ψ 0 M onto itself. One can then show that this map is generated by the linearization K µ ν of the vector field X µ that defines the evolution flow.
which shows that dΦ t is symplectic. Unitary evolution on Hilbert space leads to a flow on projective Hilbert space that preserves all three Kähler structures. However, when we project this flow onto a variational manifold to find X µ , we will project out the part of the vector field orthogonal to tangent space. When using the Lagrangian action principle, the projected flow will continue to be symplectic, i.e., preserve Ω, but none of the other two Kähler structures 17 . Geometrically, this implies that the trajectories of states near the fixed point ψ 0 will be elliptic rather than circular, when measured with respect to G.
Therefore, even if M is a Kähler manifold, K µ ν will in general neither commute with J nor be antisymmetric with respect to G, i.e., satisfy KG = −GK . This has the following consequences: , but need to be computed independently. This is important when computing spectral functions in section IV C.
• There does in general not exist a Hilbert space op-eratorK, such that K µ ν is its restriction in the sense of K µ ν = P µ ψ0K |V ν or K µ ν = P µ ψ0 iK |V ν . Thus, K is not a restriction of a Hamiltonian.
Kähler vs. non-Kähler. On a non-Kähler manifold, where we have two inequivalent definitions of the equations of motion, it only makes sense to linearize the ones coming from the Lagrangian action principle, as their fixed point coincides with the approximate ground state. The resulting generator K µ ν will in generally not commute with J , even for Kähler manifolds, which has important consequences for its eigenvectors relevant for spectral functions.

Comparison: projection vs. linearization
In the following, we will compare the previously introduced approaches of approximating excitation energy. This comparison is particularly illuminating in the case of Kähler manifold.
At a stationary point, i.e., ∂ µ E = 2Re V µ |Ĥ|ψ = 0, we consider the symplectic generator K µ ν defined as where we only have J −1 = −J for Kähler manifolds. Evaluating the derivative in (126) gives the two pieces where we evaluate everything at ψ 0 after computing the derivatives. We recognize H σ ν = P σ ψĤ |V ν and define In summary, we see that the linearization K µ ν consists of the two pieces: First, the projected Hamiltonian H and second, the derivative of the projector. These terms are multiplied with the inverse complex structure J −1 .
In the case of a Kähler manifold there is a further way to understand these two term that make up K µ ν . In this case, we can use J 2 = −1 to decompose any linear operator K on T ψ0 M as K = K + + K − with We will see that this decomposition coincides exactly with the one of K µ ν in (128). To do this, we use the fact that a Kähler manifold of dimension 2n always admits 18 a parametrization x µ = (x 1 , · · · , x 2n ) satisfying i.e., the coordinate x j is conjugate to x n+j . In this basis, J and Ω are where η jk = V j |V k . Then the structure of matrices that commute or anti-commute with J is We can evaluate K µ ν to find exactly this form where its two pieces are explicitly given by where h jk = V j |Ĥ|V k and f jk = ∂ j V k |Ĥ|ψ 0 . This clearly shows that the two pieces are given by K − = J H and K + = J F as defined before (128).
In conclusion, from the decomposition (133) we immediately see again that K µ ν has two contributions. One is related to the projected Hamiltonian H µ ν and commutes with J . The other is related to the overlap ofĤ |ψ 0 with the double tangent vectors |∂ α V β , which coincides with the one we previously described in terms of the derivative of the projector and anti-commutes with J . Thus K − is a complex linear map, while K + is a contribution that makes K µ ν non-complex linear. Finally, if we complexify tangent space, i.e., treat complex linear combinations of |V µ as linearly independent, there exists a basis transformation that makes J diagonal and brings K + and K − respectively, into block diagonal and block off-diagonal form, given by i.e., for Kähler manifolds the terms in (128) decouple nicely. For a non-Kähler manifold, neither H nor F may commute with J , but even the decomposition (129) will not work for J 2 = −1.
In the next section, we will see that the term K − can be a blessing and a curse: on the one hand, it can ensure that in systems with spontaneously broken symmetry the eigenvalues of K µ ν contain a Goldstone mode. On the other hand, for unfortunate choices of the variational family M we may encounter such massless modes even if there is no spontaneously broken symmetry (spurious Goldstone mode).
Kähler vs. non-Kähler. We can relate the linearization K µ ν with the projected Hamiltonian H µ ν via (128). For Kähler manifolds, this decomposition becomes particularly geometric, as the two pieces correspond to its complex linear and complex anti-linear part.

Spurious Goldstone mode
The spectrum of K µ ν is not variational. In contrast to a variational approximation of an eigenstate, our eigenvector |E µ (λ) of K µ ν with and |E(λ) = E µ (λ) |V µ , does not satisfy The expectation value of the full Hamiltonian with respect to |E(λ) is in general not easily related to λ, as it would be for a variational state. It is also not true that for every eigenvalue pair ±iω , there exists a true eigenstate |E ofĤ with excitation energy E − E 0 ≤ ω . In fact, there are situations, where the true ground state |E 0 is non-degenerate, but K µ ν still has a zero eigenvalue associated to a massless Goldstone mode. This typically occurs if we have a conserved quantityÂ with [Â,Ĥ], such that −iÂ |ψ ∈ T ψ M everywhere as discussed in the context of (104). At this point, the question is if the global energy minimum |ψ 0 on M is invariant under e −iÂ or not. Whenever the global minimum on |ψ 0 on M is not invariant, i.e., there is a whole family |ψ 0 (ϕ) = e −iÂϕ |ψ 0 of approximate ground states, the generator K µ ν will have a massless Goldstone mode Whenever the true ground state |E 0 of the system is invariant under e −iÂ , this Goldstone mode is spurious and merely an artefact of a spontaneous symmetry breaking on M, but not on full P(H). This was pointed out in [50] as an important problem of approximating the spectrum via linearized equations of motion rather than using the projected Hamiltonian. However, in [51] we found that this can also be desirable to capture features of the thermodynamic limit for finite system size. In particular, the gapless Bogoliubov excitation spectrum of the Bose-Hubbard model can be shown to result from the diagonalization of (123) on the manifold of coherent states. The Bogoliubov spectrum is gapless independent from the system size, even though the true ground state |E 0 becomes only degenerate in the thermodynamic limit. This also extends from Bogoliubov theory to the full Gaussian time dependent variational principle, as discussed in [51].
We illustrate this issue in figure 5, where the Hamiltonian is spontaneously broken only on the variational manifold M, but not in the full projective Hilbert space, where it has a unique ground state |E 0 .

C. Spectral functions
Next, we would like to use the variational manifold M to estimate the spectral function of a system with respect to the perturbation operatorV .
Given a Hermitian operatorV , the spectral function is where G R is the retarded Green's function with Θ(t) being the step function andV (t) the Heisenberg evolved operator under the system HamiltonianĤ.
The definition in terms of the retarded Green's function stems from linear response theory. Indeed, let us suppose that a small external perturbing probe field ϕ(t) couples to our system through the operatorV . That is, the system state |ψ (t) evolves under the perturbed HamiltonianĤ + ϕ(t)V . Then, let us measure the response of the system through the expectation value of the same observableV . As the perturbation is ideally infinitesimally small, we only consider such response up to linear order in . Consequently, we define the timedomain linear response as Then in frequency domain we have that That is, G R is exactly the so-called linear susceptibility of the system, which is an experimentally accessible quantity. Given a variational manifold there are two possible paths to trying to approximate A(ω).
1. We can calculate the quantity (142) after having projected the evolution of |ψ (t) on the manifold. In other words, we perform linear response theory directly on the variational manifold. This leads us to express A(ω) in terms of the eigendecomposition of the generator of linearized real time evolution K µ ν introduced in (123).
2. Alternatively, one can try to approximate on the manifold the quantity that appears in equation (141). In this case one should note that in generalV |ψ 0 does not belong to the variational manifold, so one has to perform some truncation even before applying the time evolution operator e −iĤt . The other subtlety here is that one must make sure that the quantity (144) is calculated with the correct global phase, as we explained in Section IV A 3.
It seems to us that method 2 captures less the spirit of variational manifolds. Indeed one has that the quan-tityV |ψ 0 would morally represent a small perturbation around the groundstate |ψ 0 and would thus naturally 6. Linear response theory. We consider an approximate ground state ψ0 ∈ M. While |ψ0 does not evolve in time, certain trajectories of nearby states are approximately elliptic. A finite perturbation at time t changes the state to |ψ (t ) = e i V |ψ0 . This state will then evolve according to the equations of motion. We linearize by taking the limit → 0, where we find that the tangent vector |δψ(t) = d d |ψ (t) | =0 can be decomposed into eigenvectors of K µ ν that rotate on ellipses.
live in the tangent space to the manifold of states at |ψ 0 . RepresentingV |ψ 0 as a vector of M therefore is only meaningful if the manifold itself is a good representation of its own tangent space. But this is not true for general manifolds and indeed there is no uniquely defined method for representingV |ψ 0 on M. The first method, on the other hand, can alternatively be thought of precisely as representing the perturbations generated byV on T ψ0 M. Furthermore, we will show that method 1 leads to a closed expression for the spectral function from which it is immediate to see that sgnA(ω) = sgn ω (as it is in the full Hilbert space), while this cannot be shown in general for method 2. For these reasons in the next subsections we will focus on the details of the first method, giving a final expression for the spectral function estimated in this way in Proposition 10.

Linear response theory
As mentioned, a possible way of calculating spectral functions is to perform linear response theory directly on the variational manifold. In this subsection we will then briefly explain how this can be done in a slightly more general setting.
Let us consider a possibly time-dependent perturba-tionÂ(t) of our unperturbed HamiltonianĤ 0 , such that H (t) =Ĥ 0 + Â (t), and an observableB, whose response we are interested in. For spectral functions, we will be interested in the particular case whereÂ(t) = ϕ(t)V for arbitrary functions ϕ(t) andB =V , but we will for the moment keep our treatment general.
Our perturbed Hamiltonian gives rise to the time dependent real time evolution vector field X µ (t), which is where X 0 and X A (t) are the evolution vector fields associated to the HamiltoniansĤ 0 andÂ(t) respectively. The solution of this perturbed evolution is |ψ (t) satisfying For the following Proposition 8 it would not be important whether the evolution vector field is defined according to the Lagrangian or McLachlan variational principles, as long as it has the form (145). However, later we will be interested in the case in which the perturbed evolution happens around the approximate ground state |ψ 0 and it will be important that this state is also a fixed point of the time evolution. So, as was the case in Section IV B 2, from now on we will suppose that the evolution vector fields are defined according to the Lagrangian evolution (88).
We are interested in the response in expectation value of the observableB at linear order in , that is where we defined the propagated perturbation which can be evaluated as follows.
Proposition 8. Given a variational manifold M we define (according to the Lagrangian action principle) the free projected real time evolution |ψ(t) as governed by the free HamiltonianĤ 0 and the perturbed projected real time evolution |ψ (t) as governed by the perturbed Hamilto-nianĤ (t) =Ĥ 0 + Â (t), both with the same initial state |ψ(0) . Then, the propagated perturbation, defined according to (148), is given by where dΦ t is the linearized free evolution flow 19 .
Proof. This can be shown in a standard way by using the interaction representation. We sketch a proof in Appendix B.
Put simply, δx µ is the superposition of all propagated perturbations, i.e., a perturbation at time t is evolved with the linearized free evolution dΦ t−t to time t where it contributes towards δx µ (t).
If we now take as initial state |ψ(0) the approximate ground state |ψ 0 , that is a fixed point of the projected evolution, we have that the free evolution is trivial |ψ(t) = |ψ 0 . It also follows that dΦ t is a linear map from T ψ0 M onto itself given by where K µ ν is the generator of the linearized flow introduced in (123). The map dΦ t can therefore be evaluated in terms of the spectral decomposition of K µ ν . Proposition 9. The linear response to a perturbation A(t), measured in terms of the observableB, for a system initially in the state ψ 0 ∈ M is given by where all derivatives are evaluated at |ψ 0 and E µ (λ ) is an eigenvector of K µ ν such that and normalised so that E µ (λ ) ω µν E ν (λ ) * = i sgn(iλ ).
Proof. We can always decompose K µ ν in terms eigenvectors E µ (λ) with eigenvalues λ and dual eigenvectors 20 E µ (λ), such that The eigenvalues λ will come in conjugate pairs ±iω , which implies that the associated eigenvectors and dual eigenvectors are complex and mathematically speaking lie the complexified tangent space. However, as K µ ν is a real map, we must have E µ (iω) = E µ (−iω) * .

Spectral response
To calculate spectral functions we now just need to evaluate the result (152) forÂ(t) = ϕ(t)V andB =V and then take the Fourier transform.
Proposition 10. The spectral function with respect to the perturbation operatorV , estimated by performing linear response theory on the variational manifold M, is where E i (iω ) are the eigenvectors of K µ ν , normalized such that E µ (iω ) * ω µν E ν (iω ) = i sgn(ω ), and the sum runs over all possible values of ω (appearing in pairs of opposite signs).
Proof. Evaluating the Fourier transform of (152) and comparing with (143) leads us to the estimate for the retarded Green's function where the Sokhotski-Plemelj formula has been used. The imaginary part of this expression can be then be inserted into equation (140), leading to the result (156).
Spectral functions calculated in this way have the desirable property sgnA(ω) = sgn ω.
Kähler vs. non-Kähler. On a non-Kähler manifold, where we have two inequivalent definitions of the equations of motion, it only makes sense to perform linear response theory with the ones coming from the Lagrangian action principle, as their fixed point coincides with the approximate ground state.

D. Imaginary time evolution
In the previous sections we have assumed we knew the state |ψ 0 that minimizes the energy on the variational manifold M. Solving this optimization problem is often 7. Imaginary time evolution. We illustrate the imaginary time evolution vector field F i on the variational family.
non-trivial and different methods may be more appropriate in different situations. However, we would like here to present a method, known as projected imaginary time evolution, that makes use of the same geometric notions introduced in Section IV A for real time evolution. On full Hilbert space, imaginary time evolution is which can be integrated to the solution This will converge in the limit τ → ∞ to a true ground state if and only if the initial state |ψ(0) had some nonzero overlap with the ground state space. Given a variational manifold M, we can approximate imaginary time evolution on it and hope that it will also converge to the approximate ground state |ψ 0 . This can be done by projecting (158) onto tangent space. Contrary to real time evolution, there does not exist a formulation of imaginary time evolution in terms of an action principle, so the projection can only be done according to the McLachlan minimal error principle.
We would like to minimize the local projection error imposing that d dτ |ψ(τ ) ∈ T ψ M, which leads to where we used P ψ |ψ = 0. This leads to the projected evolution equation from which we can define the imaginary time evolution vector field F i everywhere on M, such that This vector field can be understood as follows.
Proposition 11. Given a manifold M, the projected imaginary time evolution is given by where E(x) is the energy function, defined in the context of equation (72). Its solution x(τ ) monotonically decreases the energy.
Proof. We apply the projector P µ ψ in (163) to find We simplify this by using (74). Plugging this back into the previous equations, we arrive at (164). To show that the energy monotonically decreases, we find which follows from the positivity of G µν .
We recognize projected imaginary time evolution (164) as gradient descent of the energy function E(x) with respect to the natural geometry encoded in the metric G µν on the manifold M. This converges typically better than a naive gradient descent without taking the natural metric G µν into account. Hence, when replacing the fixed time step by a line search, imaginary time evolution becomes equivalent to Riemannian gradient descent. Indeed, the literature on Riemannian optimization [52][53][54][55] describes how the Riemannian geometry (i.e., the metric) of a manifold can be into account in each of the standard optimization algorithms such as the gradient descent method, Newton's method, the conjugate gradient method, and quasi-Newton methods such as the (limited memory) Broyden-Fletcher-Goldfarb-Shanno scheme [56][57][58][59].
Kähler vs. non-Kähler. The results discussed in this section do not rely on the manifold M being a Kähler manifold. The McLachlan projection principle is the only one that can be resonably defined for imaginary time evolution and leads to the desirable gradient descent result for any real differentiable manifold, independently of the Kähler property.

Conserved quantities
In many situations, one would like to further constrain our variational manifold by requiring that certain opera-torsÂ I have fixed expectation values A I . Geometrically, this amounts to restricting the search to the submanifold For example, for Hamiltonians commuting with the total particle number operatorN , one often wants to find lowest energy state within an eigenspace ofN witĥ N |ψ = N |ψ . To approximate such a state on a variational manifold M, we can search for minimal energy state on the submanifold of states with N = N .
In general, this manifold M will not satisfy the Kähler property anymore. In particular, if we only fix a single expectation value, we will generically reduce the dimension of a Kähler M to an odd dimension, which cannot be again a Kähler manifold. However, we have seen that for the purpose of finding the state of minimal energy, we can apply formula (158) on the reduced manifold, regardless of the Kähler property.
Instead of finding a new parametrization of the reduced variational manifold, as long as we choose an initial state for the imaginary time evolution that satisfies the desired constraints, we can then just implement them locally. We can indeed modify the imaginary time evolution vector field F by further projecting it onto the restricted tangent space T M. In this way the respective expectation values are preserved by construction.
If there are several quantitiesÂ I that we wish to fix, T ψ M is given by the sub tangent space orthogonal to the span of X µ I = P µ ψÂ I |ψ . To project onto it, we define which gives rise to the projector where G IJ is the inverse of g IJ (or pseudo inverse, if not all constraints are independent). The modified imaginary time evolution vector field is then which will conserve all the expectation values A I (τ ). In analogy to (109), this is equivalent to If we want to fix the expectation value of the number operatorN , we have the scalar function N (x) = N with X µ = P µ ψN |ψ = G µν ∂ ν N , such that which clearly satisfies dN dτ = (∂ µ N ) F µ = 0.

V. APPLICATIONS
In this section, we apply the formalism developed in the previous section to the important class of states of the form |ψ(x) = U(x) |φ where U(x) is subset of unitary transformations and |φ an appropriately chosen reference state. Here, we will consider the case where the unitary transformations form a (projective) representation U(M ) where M ∈ G for some finite-dimensional Lie group G. We will be particularly interested in the case where the reference state |φ is a highest weight vector of the group representation U(M ). Such a state is annihilated by the maximal number of certain raising op-eratorsÊ α , just like the maximal spin state |l, m = l is annihilated by the spin raising operatorĴ + =Ĵ x + iĴ y for G = SU(2). We can distinguish the following cases: where the reference state |φ(x) depends on additional parameters. As |φ(x) can include highest weight states, the full manifold will in general consist of individual sheets labelled by x, which may be partially Kähler and partially non-Kähler (making the full family non-Kähler). Such families were recently [30] used to construct non-Gaussian states with correlations between bosonic and fermionic degrees of freedom, as we discuss in section V C.

A. Gaussian states
We consider pure bosonic and fermionic Gaussian states, which are also known as quasi-free states. Bosonic Gaussian states (squeezed coherent states) form a prominent variational family in the study of bosonic systems, such as Bose-Einstein condensates [60], cold atoms in optical lattices [51] and photonic systems [61]. Fermionic Gaussian states (generalized Hartree-Fock states, including Slater determinants) are equally important for the study of fermionic systems, including Bardeen-Cooper-Schrieffer theory [3] and the Hartree-Fock method [7]. Other applications range from field theory [62], continuous variable quantum information [63], relativistic quantum information [64] and quantum fields in curved spacetime [65]. Most importantly, they are completely determined by their one-and two-point function and they satisfy Wick's theorem. Interestingly, pure Gaussian states are in one-to-one correspondence to compatible Kähler structures (g ab , ω ab , J a b ) on the classical phase space. These Kähler structures are distinct from those (g µν , ω µν , J µ ν ) on tangent space, but we will derive relations between them.
While the relationship between pure Gaussian states and Kähler structures has been pointed out before in the context of quantum fields in curved spacetime [65][66][67][68], the linear complex structure J was only recently used as a convenient parametrization of states. This allowed for a unified formulation of bosonic and fermionic Gaussian states [69,70], which enabled new results in the context of their entanglement dynamics [71][72][73], the typicality of energy eigenstate entanglement [74][75][76] and in the study of circuit complexity in quantum field theory [77,78].

Definition
The theory of N bosonic or fermionic degrees of freedom can be constructed on the classical phase space V R 2N . We denote a phase space vector by ξ a ∈ V and a dual vector by w a ∈ V * , where we use the Latin indices a, b, c, d, e for phase space to distinguish from tangent space indices µ, ν, σ, λ. We will see in a moment how this relates to standard bosonic and fermionic creation and annihilation operators.
The classical bosonic phase space V is always equipped with symplectic form Ω ab and its inverse ω ab satisfying Ω ac ω cb = δ a b , which define the classical Poisson brackets and then give rise to the canonical commutation relations (CCR). Similarly, one can define classical fermionic phase space V by equipping it with a positive definite metric G ab and its inverse g ab satisfying G ac g cb = δ a b giving rise to the canonical anti-commutation relations (CAR).
We can always choose a basis of classical linear observ-ables ξ a ≡ (q 1 , · · · , q N , p 1 , · · · , p N ), which can be bosonic or fermionic, such that (ω, Ω) or (g, G), respectively, take their standard form given by Quantization promotes these linear observables ξ a to the quantum operatorsξ a ≡ (q 1 , · · · ,q N ,p 1 , · · · ,p N ) satisfying the commutation or anti-commutation relations We refer toξ a as quadratures for bosons and as Majorana modes for fermions. These relations ensure that the bosonic or fermionic creation and annihilation operatorŝ 2 will satisfy their standard commutation or anti-commutation relations, respectively.
Given a normalized quantum state |ψ , we define its one and two point correlation functions as The fact thatξ a is Hermitian implies that z a is real. This does not apply to C ab 2 , which we therefore decompose into its real and imaginary part given by It is easy to verify that G ab is real, symmetric and positive definite with inverse g ab and Ω ab is anti-symmetric. For bosons, Ω ab is fixed by the commutation relations, while for fermions, G ab is fixed by the anticommutation relations. With respect to our basisξ a ≡ (q 1 , · · · ,q N ,p 1 , · · · ,p N ), they are given by (173). We define the covariance matrix which is the only part of C ab 2 that depends on |ψ . Later on, we will use Γ ab as parameters of our variational manifold.
We define the linear complex structure J of |ψ as as before. Note, however, that for fermions Ω ab may not be invertible, in which case we define ω ab as its unique pseudoinverse with respect to G ab , i.e., we invert Ω ab on the orthogonal complement of its kernel. We define Put differently, pure bosonic and fermionic Gaussian states are those states, for which (g ab , ω ab , J a b ) are compatible Kähler structures as defined in section III. For fermions, this also implies 23 z a = 0, while for bosons z a is a free parameter, which we call the displacement. We denote a pure Gaussian state with displacement z a and complex structure J by |J, z , which is uniquely determined (up to normalization and phase) by requiring Given a Gaussian state |J, z , we define the operatorŝ which satisfyξ a =ξ a + +ξ a − + z a andξ − |J, z = 0. Put differently, we project onto the eigenspaces of J with eigenvalues ±i and displaced by z. These spaces, i.e., complex linear combinations ofξ a + orξ a − , correspond to the spaces of creation and annihilation operators,respectively. Most importantly,ξ a ± satisfy the commutation and anti-commutation relations given by Using these relations together withξ − |J, z = 0, one can show the equivalence of (176) and (181). Moreover, one can show that any higher order correlation function with n > 2 can be computed purely from C ab 2 via Wick's theorem, which states the following: where the permutations σ satisfy σ(2i − 1) < σ(2i) and |σ| = 1 for bosons and |σ| = sgn(σ), called parity, for fermions.
This enables us to compute arbitrary expectation values of observables O written as polynomials inξ a efficiently.
Example 11. The harmonic oscillatorĤ = 1 2 (p 2 +ω 2q2 ) with standard position and momentum operators [q,p] = i 23 There exist fermionic coherent states [79], where z a is a Grassmann number, which we do not consider here.
has ground state wave function ψ(q) = ω 1/4 e −ωq 2 /2 from which we find the one-and two point functions This leads to the linear complex structure We can computeξ a − = 1 whereâ = 1 2 √ ω (ωq + ip). Example 12. The fermionic oscillator of a single degree of freedom is given byĤ = ω(n − 1 2 ), wheren =â †â withâ † andâ being fermionic creation and annihilation operators. We can go to Majorana modes q = 1 There are only two distinct fermionic Gaussian states given by |0 and |1 . The associated complex structures with respect to the basisξ a ≡ (q,p) are given by respectively. We have the annihilation operatorsξ a − ≡ . In summary, we showed how Kähler structures can be used as unifying framework to treat bosonic and fermionic Gaussian states on an equal footing. In particular, we can label pure Gaussian states |J, z by their associated complex structure J, rather than using their covariance matrix. In praxis, we can quickly switch between the different Kähler structures using their relations as reviewed in appendix A 4.

Gaussian transformations
There is a subgroup of unitary transformations that map pure Gaussian states onto themselves. We refer to them as Gaussian transformations and we will see that they are given by the exponential of linear and quadratic operators in terms ofξ a .
In this context, we will encounter the symplectic and orthogonal Lie groups and algebras. We identify them with linear maps M a b and K a b on the classical phase space V . Using the symplectic form Ω ab (bosons) and the metric G ab (fermions), we define the symplectic and orthogonal group G as We can represent Lie algebra elements K faithfully as anti-Hermitian quadratic operators K given by Using the canonical commutation, we can verify i.e., our quadratic operators K form a representation of the Lie algebra g. The respective exponentials of K give rise to a projective representation of the associated Lie group G, i.e., we have the identification between Lie group elements M and unitary operators S(M ), which we refer to as squeezing transformations. For bosons, we also have displacement transformations i.e., we identify a phase space vector z ∈ V with the respective unitary operator D(z). For fermions, products of M = e K for K ∈ so(2N, R) will only generate the subgroup SO(2N, R), whose group elements satisfy det M = 1. To generate other group elements M ∈ O(2N, R) with det M = −1, we can take any dual vector v a ∈ V * satisfying v a G ab v b = 2 to define where " ∼ =" refers to equality up to a complex phase. Using Baker-Campbell-Hausdorff, it is not difficult to show Therefore, it is a projective representation with where (Mz) a = M a bz b . The action of U(M, z) onto a Gaussian state is particularly simple and given by Given a linear complex structure J 0 , there is a unique group of transformations M preserving J 0 , i.e., This group is indeed the complex linear group GL (N, C), because the complex structure J 0 can turn V R 2N into a complex vector space V C N with scalar multiplication : With this in mind, a transformation M ∈ GL(N, C) can be seen as a linear map on V C N that commute with the representation J 0 of the imaginary unit, i.e., Consequently, the subgroup of G which also preserves J 0 is the intersection This turns out to be the group of transformations that preserve all Kähler structures (G, Ω, J 0 ) on V , which are just the unitary transformations preserving the Kähler induced Hermitian inner product on V given by Given a Gaussian state |J 0 , 0 , we can reach another 24 Gaussian state |J, z by applying the transformation where we introduced the relative covariance matrix Above transformation follows from J = e K J 0 e −K .
Example 13. The squeezing transformations of a single bosonic modeξ a ≡ (q,p) are G = Sp(2, R), generated by with their associated quadratic operators The Gaussian state |J 0 , 0 is preserved by u ∈ U(1) with where U(1) = Sp(2, R) ∩ GL(1, C) and u = e ϕZ with [Z, J 0 ] = 0. The most general Gaussian state of one bosonic mode has complex structure for which we can verify J 2 = −1. From the relative complex structure ∆ = −JJ 0 , we compute the generator such that e K |J 0 , 0 ∼ = |J, 0 .
Example 14. The squeezing transformations of a single fermionic modeξ a ≡ (q,p) are G = O(2, R), generated by and our choice of the additional group element with v a ≡ ( √ 2, 0). There are exactly two distinct complex structures given by J + and J − from (190) which are related by J + = M v J − M −1 v and which both satisfy uJu −1 = J for u = e X . This is because u is the generator SO(2, R) which is identical to the subgroup U(1) preserving J ± .
Example 15. The squeezing transformations of two fermionic modesξ a ≡ (q 1 ,p 1 ,q 2 ,p 2 ) are G = O(4, R) with six linearly independent generators K i satisfying KG + GK = 0. Given the linear complex structure there is a 4-dimensional subspace of these generators also satisfying [K, J 0 ], which generates U(2) ⊂ O(4, R). The most general fermionic complex structure is (215) for ∆ = J + J 0 . From the explicit form of J ± , we see that (θ, φ) behave as spherical coordinates. This agrees with the fact that the manifold of fermionic Gaussian states with two modes consists of two disjoint spheres.
In summary, Gaussian transformations consist of squeezing (bosons and fermions) and displacements (only bosons). The former changes both the complex structure (equivalently: covariance matrices) and displacements, while the latter only displaces the state. Whenever we choose a Gaussian state |J 0 , z 0 , it defines with the respective background structure (symplectic form for bosons or metric for fermions) a subgroup U(N ) of transformations in G that preserve J 0 .

Geometry of variational family
We will now shift gears. Rather than looking at the Kähler structures (g ab , ω ab , J a b ) associated to an individual Gaussian state |J, z , we now consider the full variational family M of pure bosonic or fermionic Gaussian states to study the Kähler structures (g µν , ω µν , J µ ν ) on tangent space T ψ M. For this, we will need to compute the tangent basis vectors |V µ , which we will split into two types |V a and |V ab .
We introduced Gaussian states |J, z as being uniquely (up to a complex phase) characterized by their complex structure J and their displacement vector z. Consequently, a general tangent vector is characterized by the pair (δJ a b , δz a ) describing the respective changes of J a b and z a . However, the resulting expressions simplify if we express everything in terms of the (bosonic or fermionic) covariance matrix Γ ab as defined in (178). We will therefore label states by |Γ, z and tangent vectors by (δΓ, δz) with |δΓ, δz = δΓ ab |V ab + δz a |V a ∈ H ⊥ |Γ,z .
The tangent space T (Γ,z) M thus decomposes as where S (Γ,z) corresponds to the changes of J and D (Γ,z) refers to the changes of z (only for bosons). We will see that these spaces can be naturally identified with the space of one-and two-particle excitations in Hilbert space (constructed as Fock space). Squeezing tangent space S (Γ,z) . The tangent space of squeezings is described by variations δΓ ab of the covariance matrix. Note that such changes are constrained to preserve the complex structure property J 2 = −1 and reflect the symmetry/antisymmetry of Γ for bosons/fermions, respectively. We therefore have . (220) The set of tangent vectors |V ab is overcomplete, but this does not cause any problems as for any change δΓ ab , the associated tangent vector is given by |∆Γ = δΓ ab |V ab . Using Wick's theorem, we can compute the inner product δΓ|δΓ = 1 16 (g ac g bd + iω ac g bd )δΓ ab δΓ cd .
We thus see that the inner product between squeezing tangent vectors can be computed from the Kähler structure (g, ω, J) on the classical phase space V . Displacement tangent space D (Γ,z) . The tangent space of displacements is described by the variations δz a of z a , which can be changed freely. Therefore, we can identify the tangent space T |Γ,z with the classical phase space V , i.e., δz a ∈ V . The local frame |V a is We find the inner product which implies that the tangent space is isomorphic to V embedded with metric g ab and symplectic form −ω ab . 26 The spaces S (Γ,z) and D (Γ,z) are orthogonal, because V a |V bc = 0 follows from C abc 3 = 0 in Wick's theorem. Kähler structures on tangent space. The tangent space of bosonic Gaussian states can be decomposed into the direct sum of displacements and squeezings, which are orthogonal due to V ab |V c = 0. The tangent space of fermionic Gaussian states, only consists of the squeezings. Evaluating the respective inner products gives δΓ, δz|δΓ, δz = 1 16 (g ce g df − iω ce g df )δΓ cd δΓ ef 25 The most general tangent vector of the squeezing manifold is given by |V K = Q (Γ,z) K |Γ, z , where K is a general quadratic operator from (193). We can then relate |V K to the change δΓ = 2Re Γ, z|(ξ aξb +ξ bξa − 2z a z b )|V K for bosons and δΓ = 2Im Γ, 0|(ξ aξb −ξ bξa )|V K for fermions to find (220). For bosonic states |Γ, z with z a = 0, |V K also has a component |δz = δz a |Va with δz a = K a b z b in the displacement tangent space D (Γ,z) , discussed next. 26 The sign difference is due to our chosen conventions and related to the issue discussed in footnote 27. giving rise to the Kähler structures Here, we do not have V ab |V cd = 1 16 (g ac g bd −iω ac g bd ), but only δΓ|δΓ = 1 16 (g ac g bd − iω ac g bd )δΓ ab δΓ cd , where δΓ ab needs to satisfy (219). Inverting them on this subspace leads to the dual Kähler structures They are defined on the dual subspace spanned by (w a , W ab ) with W ab satisfying dual to (219). Given a general W ab , we can define its projection W onto this subspace as with W (ab) = 1 2 (W ab + W ba ) and W [ab] = 1 2 (W ab − W ba ) referring to symmetrization and anti-symmetrization, respectively. The resulting projection W ab satisfies W ab δΓ ab = W ab δΓ ab for δΓ ab satisfying (228). We find the action of i onto |δΓ, δz to be Example 16. We consider the Gaussian state |J 0 , 0 of a single bosonic modeξ a ≡ (q,p) with Kähler structures A change δΓ of the covariance matrix Γ = G is constrained to satisfy (219), such that we have The associated tangent vectors are given by which leads to a general tangent vector Example 17. We saw in example 14 that fermionic Gaussian states for N = 1 consists of two points and thus is zero dimensional. Instead, we consider for N = 2 the state |Γ, 0 and allowed change δΓ with The non-zero tangent vectors are then In summary, the Kähler structure (g ab , ω ab , J a b ) on the classical phase space V are intimately linked to the ones (g µν , ω µν , J µ ν ) on tangent space. For bosons, we saw that the displacement tangent space D |J,z can be naturally identified with phase space V as Kähler space.

Variational methods
After having clarified the Kähler geometry of bosonic and fermionic Gaussian states, we can use them to illustrate the variational methods discussed in section IV.

Applications include (A) real time evolution, (B) excitation spectra, (D) spectral functions and (D) imaginary time evolution.
Real time evolution. Based on (88), we have the Lagrangian evolution equation dx µ dt = −Ω µν ∂E ∂x ν . As Gaussian states are parametrized by x µ = (Γ, z), we need to compute ∂E ∂Γ ab and ∂E ∂z a . As Ω µν is only defined on the subspace satisfying (228), we need to project ∂E ∂Γ onto this subspace to find d dt Γ ab = 8Ω ac G bd ∂E ∂Γ cd and thus 27 in agreement with the respective equations 28 of [30]. We introduced the evolution vector field X µ = (X a , X ab ). We can also define the instantaneous Hamiltonian 27 By construction the energy will only depend on the symmetric or antisymmetric part of Γ ab for bosons and fermions, respectively, such that ∂E ∂Γ will be automatically symmetric or antisymmetric. Note the sign difference inż = Ω∂E compared toẊ = −Ω∂E, which is due to the chosen conventions mentioned in 26. 28 See equations (31) of [30], where z = ∆ R / √ 2, Ω = σ and Γ = G = Γ b for bosons and G = 1 and Γ = Ω = −Γm for fermions. with l = ∂E ∂z and k = ±2 ∂E ∂Γ , where (+) applies to bosons and (−) to fermions. This is the quadratic Hamiltonian whose time evolution on the Gaussian family agrees with the projection X µ , i.e., −iĤ |Γ, z = X µ |V µ . We further define the instantaneous Lie algebra element which allows us to write the time evolution equation for the linear complex structure and displacement simply aṡ Note that these equations are in general non-linear, as K and l depend on the state and thus on J and z. Imaginary time evolution. We recall that on a general manifold, we have the the imaginary time evolution vector field dx µ dτ = −G µν ∂E ∂x ν . This translates to d dτ Γ ab = 8Ω ac G bd ∂E ∂Γ cd and thus to 29 which agrees with the respective equations 30 in [30]. Excitation spectrum. We recall that we can approximate the excitation spectrum either by linearizing the equations of motion or by projecting the Hamiltonian onto tangent space. The latter can be straightforwardly done in the number basis by expressingĤ in terms of creation and annihilation operators associated to the approximate ground state |ψ 0 . We therefore focus here on the former case, where the spectrum is encoded in the eigenvalues of K µ ν = −Ω µσ (∂ ν ∂ σ E). We evaluate K for Gaussian states at a stationary point, i.e., at x 0 with (∂ µ E)(x 0 ) = 0. We use the real time evolution vector field X µ = (X a , X ab ) from (239) to compute

Displacements (bosons) Squeezings (bosons) Squeezings (fermions)
Change δz a δΓ ab δΓ ab Constraints none Γ ab = Γ ba , ∆ΓJ = JδΓ which can be explicitly computed as 31 Note that for fermions, we only have the last block ∂X ∂Γ , as the others are related to the displacement z a , which vanishes for fermions. For bosons, the spectrum of the first block ∂X a ∂z c can be related to Bogoliubov mean field theory, as discussed in [51], i.e., it captures the oneparticle spectrum.
Spectral functions. We can then evaluate spectral functions based on formula (156) for any operatorV . In practice, all eigenvectors E µ (λ ) will be represented as such that we can compute the dual vectors Note that we will need to remove unphysical eigenvectors and dual eigenvectors, as discussed in the next paragraph. The spectral function is then directly computed from (156). This was done explicitly in [51] to study the excitation spectrum and spectral functions of the paradigmatic Bose-Hubbard model, whereV was chosen to either describe density variations or lattice modulations.
As explained in table III, the manifold of bosonic states is N (N + 3)-dimensional, while the manifold of fermionic states is N (N − 1)-dimensional. In contrast, the matrix representation (246) has a much larger dimension, due to the fact that we do not implement the constraints on Γ directly in the basis of K, but rather by applying the projection (229) in the definition of K. By construction, any 31 We make the assumptions discussed in footnote 27. forbidden change δΓ violating (219) will be projected out and thus contributes the eigenvalue zero to the spectrum of K. In practice, we therefore have two options: (a) We can just compute the spectrum of K as written in (246) and drop N (N ∓ 1) vanishing eigenvalues for bosons and fermions, respectively. If there are more zero eigenvalues, the spectrum will have one or more zero modes. To identify the corresponding (physical) eigenvector (z c , δΓ cd ), one would need identify those eigenvectors which do not violate (219). All eigenvectors with non-vanishing eigenvalues are physical and necessarily satisfy the constraint.
(b) We can also reduce the dimension of K by constructing an orthonormal basis δΓ cd i explicitly. This allows us restrict/project K onto this subspace. This is particulary important when There are large classes of examples where variational methods have been successfully applied to the family of Gaussian states to describe physical systems. Given a single fixed unitary transformation U 0 , we can define the variational family of transformed Gaussian states M = {U 0 |Γ, z }. Using this variational family for a given HamiltonianĤ is equivalent to applying our methods directly to the transformed HamiltonianĤ = U † 0Ĥ U 0 . This approach has been successfully applied to various systems, such as the Kondo problem [27].
Example 18. We consider a free bosonic system with one degree of freedom. Let its quadratic Hamiltonian bê The manifold of Gaussian states was explicitly parametrized in (210), such that we find the The change of z a and Γ ab under time evolution are dz dt ≡ ω leading toρ = 0 andφ = 2ω. Similarly, imaginary time evolution is given by dΓ dτ ≡ −ω 2 sinh 2 ρ + sin φ sinh 2ρ cos φ sinh 2ρ cos φ sinh 2ρ leading to ρ = −2 sinh ρ and φ = 0. Finally, we find from which we can read off the 1-and 2-particle excitation energies of the free Hamiltonian. For the matrix representation of K, we chose the orthonormal basis (z a , Γ i ) Example 19. We consider a free fermionic system with two degrees of freedom. Let its quadratic Hamiltonian bê This manifold of Gaussian states was explicitly parametrized in (215) with two coordinates x µ ≡ (θ, φ) leading to the energy of the state with J ± given by The resulting equations for real time evolution are leading to the equationsθ = 0 andφ = (ω 2 ± ω 1 ). Similarly, imaginary time evolution is described by leading to θ = (ω 2 ± ω 1 ) sin θ and φ = 0. We can compute the matrix representation of K as whose eigenvalues correspond to the 2-particle spectrum of the free HamiltonianĤ. For the matrix representation of K, we chose the orthonormal basis of variations Example 20. We now derive the time evolution for example 9 of a free bosonic system with two degrees of freedom. We fix a basisξ a ≡ (q 1 ,p 1 ,q 1 ,p 1 ), with respect to which symplectic form and covariance matrix of |Γ, z are We now choose the 2-dimensional variational families of coherent states |Γ, z , where z a ≡ (q 1 , p 1 , q 2 , p 2 ) depends on the two variational parametersz α = (q,p) via which agrees with |ψ(α) from example 7. Thez α are the expectation values of the canonically conjugated operatorŝ p =p 1 cosh r −p 2 sinh r .
We represent the orthogonal projector P (Γ,z) as 2-by-4 where P acts on the tangent space of all displacements δz from (222) and orthogonally projects onto the tangent space of our family described by δz. The Hamiltonian (101) can be rewritten asĤ where c 1 = 1 cos 2 φ + 2 sin 2 φ, c 1 = 1 sin 2 φ + 2 cos 2 φ and c 3 = 1− 2 2 sin 2φ. We consider the time evolution underĤ of the expectation valuesz ≡ (q,p) for the following scenarios. True evolution. The time evolution equationż a = X a follows from (239) and is given by where we introduced the constants given by We can compare the time evolution of the variational parametersz α ≡ (q,p) for the two variational principles with the exact evolution of the expectation values of (q,q) for the same initial state |Γ, z , as shown in figure 4.

Approximating expectation values
In this section, we compare how the expectation value of observables changes depending on if we evolve in the full Hilbert space with −iĤ or on our variational manifold with P ψ (−iĤ). Clearly, any observable O that can be expanded in powers ofξ a can, in principle, be computed exactly from the covariance matrix Γ ab and displacement vector z a (for bosons) using Wick's theorem.
As it turns out, the derivative d dt Â for linearquadratic observablesÂ agrees in the two cases on the Gaussian manifold, i.e., linear-quadratic observables are insensitive to non-Gaussianities to linear order on the Gaussian manifold. Put differently, the variation δ Â is insensitive if we perturb the Gaussian state |Γ, z either into the direction |δψ , which will generally leave the class of Gaussian states, or into the projected direction P |Γ,z |ψ , which is projected onto the Gaussian manifold. Here, we have δ Â = d dt Â if we choose |δψ = d dt |ψ = −iĤ |ψ . Proposition 12. Given a Gaussian state |Γ, z and an arbitrary tangent vector |φ with Γ, z|φ = 0, the change of linear-quadratic observablesÂ = f aξ a + 1 2 h abξ aξb is Proof. The proof is rather simple and goes in two steps: First, we note that a linear-quadratic operatorÂ acting on |Γ, z allows for the decomposition Second, the inner product betweenÂ |Γ, z and the tangent vector |φ only happen in this subspace. Consequently, equation (277) follows.
Corollary 1. The time derivatives of displacement vectorż a and covariance matrixΓ ab at a Gaussian state |Γ, z are the same for the true time evolution with some interacting HamiltonianĤ and for its projection onto the Gaussian manifold. In formulas, this meanṡ where Γ = G for bosons and Γ = Ω for fermions.
In summary, projecting onto the Gaussian manifold is equivalent to truncating the equations of motion of the n-point functions at second order, i.e., we integrate equations (279) to (281) and ignore non-Gaussian evolution of higher n-point functions C n for n > 2.

B. Group theoretic coherent states
Standard coherent states of the form |α = e αâ † −α * â |0 of a single bosonic degree of freedom were introduced by Glauber [42]. From the perspective of Gaussian states, as presented in the previous section, coherent states correspond to a submanifold of bosonic Gaussian states with fixed complex structure J a b (or covariance matrix Γ ab ), but variable displacement vector z a . Given a fixed Gaussian state |J, 0 , we can reach the set of coherent states |J, z with all possible z a ∈ V by applying the displacement transformation |J, z = D(z) |J, 0 . Here, D(z) forms a projective representation of the group of vector addition V .
In this section, we generalize this construction to semisimple Lie groups to construct group theoretic coherent states. They were independently introduced by Gilmore [24,80] and Perelomov [25,81]. While we follow the excellent review [82], we will particularly emphasize the geometric structure of the resulting variational families and their advantages. In particular, we will show that group theoretic coherent states constructed from highest weight vectors always give rise to Kähler manifolds. Furthermore, we will connect to the previous section to show explicitly how the full families of bosonic and fermionic Gaussian states can be naturally understood as group theoretic coherent states with respect to groups G = Sp (2N, R) and G = O(2N, R) for bosons and fermions, respectively.

Definition
We consider a separable Hilbert space H, a real Lie group G with real Lie algebra g and a possibly projective unitary representation U of G on H, i.e., we have indicates equality up to a complex phase. We may later impose conditions, such as requiring that the Lie group is compact or semi-simple. Given a basis Ξ i of the Lie algebra g, we have the Lie brackets 32 where c k ij are called structure constants. Our group representation induces a representation of Ξ i as operatorŝ which are anti-Hermitian 33 due to U † = U −1 . A general Lie algebra element A ∈ g is then represented as anti-Hermitian operatorÂ = A iΞ i . 32 In physics, some authors [83] choose the basis X i = −iΞ i , such that A = A i Ξ i = iA i X i and [X i , X j ] = ic k ij Ξ k . 33 A basis X i = −iΞ i would lead to Hermitian operatorsX i = −iΞ i . We can always represent group and Lie algebra through their action on the Lie algebra itself. This is called the adjoint representation. Here, a Lie group element M is represented as the linear map Ad M : g → g with which reduces for matrices to Ad Similarly, the adjoint representation of a Lie algebra element Ξ i is given by the linear map ad i : g → g with which implies that the adjoint representation of Ξ i always takes the matrix form (ad i ) k j = c k ij with respect to Ξ j . The adjoint representation defines the Killing form which is non-degenerate for semi-simple Lie groups.
Example 21. We consider the Lie group SU(2) consisting of complex, unitary 2-by-2 matrices with unit determinant. Its algebra is well known to be spanned by the matrices Ξ i ≡ i 2 σ i , where σ i are the Pauli matrices, with structure constants c k ij ≡ k ij given by the totally antisymmetric tensor with 1 11 = 1. The Killing form is thus K ij = k il l jk = −2δ ij . We consider the following two representations: Spin-1 2 . This is the fundamental representation and thus coincides with the definition. We represent the group on the Hilbert space H = C 2 and it is generated by the anti-Hermitian operatorsΞ i ≡ − i 2 σ i witĥ which coincides with the definition Ξ i . Spin-1. This representation is also the adjoint representation and the matrices can be chosen to be real. It is then given by real 3-by-3 rotation matrices acting on the Hilbert space H = C 3 . It is generated bŷ which are infinitesimal rotations in the three planes. As matrices, they correspond to (Ξ i ) k j = k ij , which confirms that this is the adjoint representation.
We will now explain how the choice of a reference state φ together with a projective representation U(M ) on some Hilbert space H defines a variational family M φ ⊂ P(H).
When applying variational methods to M φ , it is useful to parametrize states by group elements, i.e., |ψ(M ) = U(M ) |φ . Rather than taking derivatives with respect to some artificial coordinates, we define local coordinates around every state |ψ(M ) with tangent vectors at x = 0 This allows us to introduce a map between the Lie algebra g and the tangent spaces of the manifold M φ as follows.
Definition 5. For any M ∈ G, we associate to every Lie algebra element A = A i Ξ i the tangent vector This map is only an isomorphism if h φ = 0, i.e., there are no Lie algebra elements A that only generate a change of complex phase and are thus mapped to A i |V i = 0. Such Lie algebra elements define a subalgebra h φ generating the stabilizer group H φ , defined as As a manifold, we thus have If the Lie algebra is semi-simple, we can use the then non-degenerate Killing form to choose a uniquely represent g/h φ as We can now proceed to calculate the restricted Kähler strucures for the manifold M φ . We can use the Lie algebra induced tangent space vectors |V i , introduced in (292), to derive simple expressions of the restricted Kähler structures independent of M . Proposition 13. The restricted Kähler structures are which are independent of M and thus everywhere the same.
Proof. We can straightforwardly compute Proposition 13 implies a dramatic simplification of the variational manifold M φ , because it suffices to choose a single basis Ξ i of generators to bring the Kähler structures into a standard form which extends to all tangent spaces via the map of definition 5. This is particularly important for numerical implementations as discussed in section V B 5.
Example 22. We will reconsider the Lie group SU(2) from example 21. As we are interested in the possible families M φ , we will need to understand which φ are inequalivant, i.e., do not give rise to the same family M φ . We can compute the tangent vectors |V i based on definition 5 for the following representations. Spin-1 2 . Up to a multiplication with a complex number c, we can use our fundamental representation U(M ) to transform any non-zero vector into any other complex vector. As our definition of M φ ignores complex rescalings, we therefore must have M φ = P(C 2 ), i.e., for any non-zero state φ, the resulting family is the full projective Hilbert space as discussed in example 4. To compute the tangent space vectors, we take |φ ≡ (1, 0) and find which clearly satisfy J 2 = −1. Spin-1. Our representation U(M ) consists of standard 3-by-3 rotation matrices acting on C 3 . In the basis of these matrices, we can thus represent any vector as column |φ ≡ a + i b, with a, b ∈ R 3 . We choose this vector to be normalized, such that a 2 + b 2 = 1. Multiplying with a complex phase e iϕ corresponds to a rotation a → cos(ϕ) a + sin(ϕ) b and b → cos (ϕ) b − sin(ϕ) a, which we can use to ensure a · b = 0 and a 2 ≥ b 2 , such that we can always choose 0 ≤ θ ≤ π 4 with | a| = cos θ and | b| = sin θ. We can then apply the rotation matrices U(M ), such that a = (cos θ, 0, 0) and b = (0, sin θ, 0). Furthermore, we find the tangent vectors For this choice, we can compute the Kähler structures leading to the linear complex structure For 0 < θ < π 4 , we have h φ = span(|V 1 , |V 2 , |V 3 ) and M φ is degenerate non-Kähler. For θ = 0, we have h φ = span(|V 2 , |V 3 ), on which ω vanishes and M φ is again degenerate non-Kähler. Only for θ = π 4 , we have h φ = span(|V 1 , |V 2 ), on which J 2 = −1 holds, and M φ is Kähler. The families M φ ⊂ P(C 3 ) are 3-dimensional copies of SO(3, R) for 0 < θ < π 4 and spheres S 2 , otherwise. Together these orbits foliate the projective Hilbert space P(C 3 ) just like a sphere can be foliated in circles of latitudes with single points at the poles.

Compact Lie groups
We are interested in geometry of the variational family M φ ⊂ P(H), namely its restricted Kähler structures (297) and (298). In particular, we would like to find a simple criterion when such a family of group theoretic coherent states is a Kähler manifold. In this section, we will focus on the representation theory of compact semisimple Lie algebras, as discussed in [83]. Later, we will also consider Lie algebras that are not compact or not semi-simple, as discussed in [84].
For readers familiar with the theory of Lie algebras, let us emphasize that we need to carefully distinguish between the theory of real and complex Lie algebras. In our case, we have a real Lie algebra g, because we started from a real Lie group G and its unitary representation U. We will be lead to also consider the complexification but only real elements A ∈ g ⊂ g C will be represented by anti-Hermitian operatorsÂ. For a general element A = A i Ξ i ∈ g C , we define its conjugation as We consider a compact semi-simple Lie group G with real Lie algebra g, i.e., the Killing form K is negative definite, such that K(A, B) < 0 for all non-zero A, B ∈ g. Our construction relies on first choosing a Cartan subalgebra h ⊂ g (not to be confused with h φ ), characterized by the property that if [X, Y ] ∈ h for all X ∈ h also Y ∈ h. For semi-simple Lie algebras, this implies that h is Abelian, i.e., [X, Y ] = 0 for all X, Y ∈ h. We choose a basis of h as H I = H i I Ξ i , which is smaller than Ξ i that spans g. While the choice of h is not unique, they are all isomorphic for compact semi-simple Lie algebras and they give rise to the following structures 34 : • The adjoint representation ad I = H i I ad i of the Cartan basis H I has joint eigenspaces v α ⊂ g C with 34 The same structures arise for non-compact semi-simple Lie algebras, but they are not necessarily isomorphic anymore.
where the set of eigenvalues α I is called a root 35 , which always come in pairs (α, −α).
• The root system ∆ is the set of all non-zero roots α. It can be split into the two disjoint sets of positive roots ∆ + and negative roots ∆ − , such that for every α ∈ ∆ + , we have −α ∈ ∆ − .
• We have the root space decomposition given by where all v α are complex and one-dimensional.
• For every root α ∈ ∆, we have the generator 36 • We can choose for each eigenspace v α an eigenvector E α , such that K(E α , E −α ) > 0 and such that where N αβ are real and satisfy N α,β = −N −α,−β , while all other brackets [E α , E β ] vanish.
One can use the property of g being compact to show that above choices of E α satisfy E * α = −E −α and that all α I are imaginary. Thus, we have the real basis where α ∈ ∆ + and we introduced the real generators When we represent Ξ i by anti-Hermitian operatorsΞ i on a Hilbert space H, allĤ I = H i IΞ i commute with each other. A common eigenvector |µ ∈ H of allĤ I witĥ is called a weight vector and the joint eigenvalues µ I are called its weight 37 . A weight vector |µ is called highest weight vector 38 if it is annihilated by all positive root operatorsÊ α with α ∈ ∆ + , i.e., |µ highest weight ⇔Ê α |µ = 0 ∀ α ∈ ∆ + . (313) 35 Each root is a linear map α : h → C with α(A I H I ) = A I α I . 36 There is a slight abuse of notation: H I refers to a real basis of h and is distinct from Hα ∈ h C defined in (308). 37 Weights are linear maps µ : h → C with µ(A I H I ) = A I µ I . 38 Technically, we could also call it lowest weight vector, if we switched the roles of positive and negative roots, which is why some authors use the term 'extremal weight'. Different choices of h and different assignments of which roots are positive will lead to different highest weight vectors. We refer to all vectors |µ as highest weight vectors, for which there exists a choice of Cartan subalgebra h and an assignment of positive roots, such that |µ is a highest weight vector in the above sense. For compact Lie groups, all such highest weight vectors are related by applying U(M ) for some M ∈ G, i.e., the family M φ for φ being a highest weight vector is actually the set of all highest weight vectors with respect to all possible choices of Cartan subalgebras and positivity of roots. The following proposition shows that such M φ is Kähler.

Proposition 14.
If φ is a highest weight vector and G a semi-simple compact group, the manifold M φ is Kähler.
Proof. As the group is compact, because of the discussion in section V B 2, a basis of the corresponding algebra is given by (310). Let |φ = |µ the highest weight vector with respect to the Cartan subalgebra h C spanned by H α . We split our set ∆ + of positive root into thoseα ∈ ∆ + withÊ −α |µ = 0 and those α ∈ ∆ + withÊ −α |µ = 0. Note thatÊ α |µ = 0 for all α ∈ ∆ + due to |µ being highest weight. With this, we can split our basis Ξ i further into the three parts Ξ i ≡ (H I , Qα, Pα, Q α , P α ). We can construct the induced tangent space basis |V i from their definition (5) to find Here, (|V 1 α , |V 2 α ) forms a basis of tangent space. Indeed E −α |µ is an eigenvectors ofĤ I with eigenvalue µ I − α I , such thatĤ IÊ−α |µ = (µ I − α I )Ê −α |µ and as such are orthonormal. Clearly, we have the pairs |V 1 α and |V 2 α = i |V 1 α ensuring that tangent space satisfies the Kähler property.
Assuming |µ is a weight vector, i.e., an eigenvector of the Cartan subalgebra operators, we can explicitly compute the matrix representations with respect to (|V 1 α , |V 2 α ). We see immediately that such structures then take the canonical Kähler form (J 2 = −1) only if |µ is the highest weight, that is µ|Ê αÊ−α |µ = 0. In this case, they only depend on the following factor which we evaluate explicitly using (308): This can be explicitly checked in the following example.
Example 23. We reconsider example 22 of SU(2) for the spin-1 2 and spin-1 representation. SU(2) is a compact group. We choose as real Cartan subalgebra h = span R (Ξ 3 ), which gives rise to E ± = i 2 (Ξ 1 ± iΞ 2 ), where we E + corresponds to the only positive root. Spin- 1 2 . There is only one M φ , which is the full space P(C 2 ). Therefore every state in the representation is a highest weight state with respect to some choice of h C and selection of positive roots. For our choice, the highest weight state is |φ = (1, 0), for which we already computed the tangent vector in example 22 and verified that M φ is Kähler. Spin-1. Not every state is a highest weight state anymore. With respect to our choice of positive root vector E + , the highest weight state is |φ = ( 1 √ 2 , i √ 2 , 0), which corresponds to the boundary case θ = π/4 from our previous considerations. This agrees with our finding that M φ is only Kähler for θ = π 4 . Note that we can always include more generators to construct a larger group (here: SU(3)), so that also states for π 4 = θ = 0 are highest weight states (with respect to the larger state) and M φ will be Kähler (here: full P(H)). For large Hilbert spaces, we need to strike a balance between the properties of M φ and its dimension.

General Lie groups
Proposition 14 shows that M φ is a Kähler manifold if G is a compact semisimple Lie grup and φ a highest weight vector. How much of this analysis can be carried out for non-compact or non-semi-simple Lie groups/algebras? For a semi-simple, non-compact real Lie algebra g, i.e., K is not negative definite anymore (but still nondegenerate), we can still choose a Cartan subalgebra h ⊂ g and use the same root space decomposition (307). Not all choices of h are isomorphic anymore, as K restricted to h may have different signatures. This may lead to additional requirements for a highest weight vector |µ to give rise to Kähler manifolds. For any choice of Cartan subalgebra h and a positive root system ∆ + , we consider representations 39 with unique highest weight vector |µ annihilated by all positive root operatorsÊ α . We can distinguish the following cases: Compact Cartan subalgebra. We refer to a chosen Cartan subalgebra h as compact if the Killing form K restricted to h is negative definite. In this case, all roots α ∈ ∆ are imaginary, i.e., all α(H I ) ∈ iR, and the positive roots ∆ + can be divided into two sets 40 : the set ∆ k of compact roots and the set ∆ p of non-compact roots. The associated root operators E α satisfy E * α = ∓E −α with (−) for compact roots and (+) for non-compact roots. A real basis of the real algebra g is then where α k ∈ ∆ k , α p ∈ ∆ p . At this stage, the proof of proposition 14 can be directly applied to above basis Ξ i and M φ is Kähler. Non-compact Cartan subalgebra. Whenever the Killing form restricted to our chosen Cartan subalgebra is not negative definite, some roots α ∈ ∆ may be nonimaginary, in which case M µ constructed from the associated highest weight vector |µ may not be Kähler. Only if all the non-imaginary root operatorsÊ α annihilate the highest weight state |µ , i.e., E α |µ = E −α |µ = 0 for all non-imaginary roots α, we can once again apply the proof of proposition 14. Indeed, in this case the basis of the real Lie algebra g will be given by (323) plus some real basis vectors constructed from the non-imaginary root spaces. Those additional basis vectors, however, annihilate the reference state and thus do not play a role in the properties of the tangent space and M φ will again be Kähler. In summary, a highest weight vector |µ will give rise to a Kähler manifold M µ if it is not only annihilated by the positive imaginary root operatorsÊ α , but also by all non-imaginary root operatorsÊ ±α .
Much less is known for general Lie groups that are not semi-simple, because their Lie algebras are still not classified. Instead one can attempt to apply the presented analysis case by case. Most importantly, this analysis works for the prominent family of regular coherent states, constructed from the Heisenberg algebra that is not semisimple, as we will explain in example 25.
weights FIG. 8. Root and weight diagrams for sp(4, R) and so(4, R). We consider two bosonic (left) and two fermionic (right) modes with representation S(M ) of squeezing transformations G. We chooseĤI ≡ (in1 ± i 2 , in2 ± i 2 ) with eight bosonic and four fermionic root vectors. The root vectors are purely imaginary, allowing us to represent Im(αI ) by arrows (red = positive, green = negative) in two dimensions. The weight vectors (black = even sector, blue = odd sector) are the number eigenstates |n1, n2 with highest weight vector |0, 0 .
even and odd part of the Hilbert space, i.e., H = H + ⊕H − where H ± are the eigenspace of the parity operator P = e iπN withN being a total number operator. For fermions, the group representation is irreducible due S(M v ) defined in (197), which mixes the two sectors, while the algebra representation still splits over the even and odd sector, as the bosonic case. The algebra is N (2N + 1) dimensional for bosons and N (2N − 1) dimensional for fermions. We can always select creation and annihilation operatorŝ a † i andâ i with number operatorsn i =â † iâ i . This allows us to choose the N -dimensional Cartan subalgebra, whose basis is represented as H I ≡ (in 1 ± i 2 , . . . , in N ± i 2 ), where (+) applies to bosons and (−) to fermions. We divide the associated N (2N − 1 ± 1) roots into positive and negative roots and pick the corresponding basis vectors aŝ With this choice of real Cartan subalgebra all roots are imaginary. Only the roots associated toâ iâj andâ † iâ † j for bosons are non-compact, while all others are compact. We verifyÊ † α = ±Ê −α , with (+) for compact and (−) for non-compact roots. Vice versa E * α = ∓E −α , with (−) for compact and (+) for non-compact roots. With this choice, the weight vectors of the representation on H + are the states |n 1 , · · · , n N with fixed excitation numbers n i for alln i and i n i even. The highest weight vector is |0, . . . , 0 . The representation of the representation on H − is the same, except that we require i n i to be odd, such that the highest weight vector is |1, 0, . . . , 0 . For fermions, the highest weight families M φ of the two sectors are actually a single family (consisting of two disconnected components), because we also represent group elements M with det M = −1 in our representation. For bosons, only the family constructed from the highest weight, |0, . . . , 0 , in the even representation is called Gaussian. Let us highlight that the family M φ for |φ = |1, 0, . . . , 0 is an interesting Kähler family whose properties are not fully explored. For N = 2, the Cartan subalgebra is 2-dimensional, which allows us to plot roots and weights in the plane, as done in figure 8.
Example 25. Group theoretic coherent states are a generalization of the well-known bosonic coherent states used in physics. These are generated from the displacement transformations D(z) introduced in (196). While D(z) is a projective representation of the Abelian group of phase space translations (V, +) with D(z)D(z) = e iz a ω abz b D(z + z), this is actually not true for its Lie algebra. Instead, one can consider e iϕ D(z) as a representation of the Heisenberg group with (2N + 1)-dimensional Lie algebra represented by the anti-Hermitian operatorsΞ i ≡ (iq 1 , ip 1 , . . . , iq N , ip N , i1). This Lie algebra is not semisimple and the standard approach of computing roots fails. However, if we extend the Lie algebra even further to also include the N elements represented as number operators in i = iâ † iâ i = i 2 (q 2 i +p 2 i − 1), we can construct a root diagram. The Cartan subalgebra is then represented byĤ I ≡ (in 1 , . . . , in N , i1) with root vectorsâ i and a † i . Consequently, the roots form an orthonormal basis, the eigenstates |n 1 , . . . , n N of the number operatorsn are the weight states and |0, . . . , 0 is the highest weight state, as for bosonic Gaussian states.
To summarize we have the following cases: • Semi-simple, compact algebra. Any compact group G has a compact Lie algebra g, for which the family M φ is Kähler if |φ is the highest weight state of the representation. See example 24 of fermionic Gaussian states.
• Semi-simple, non-compact algebra. If the group is non-compact, not all highest weight vec-tors give rise to Kähler manifolds. For every highest weight vectors associated to the real Cartan subalgebra h ⊂ g, we need to split the corresponding roots into imaginary and non-imaginary roots.
Only highest weight vectors that are also annihilated by all (not just the positive) non-imaginary root vectors E α will give rise to Kähler manifolds. This includes in particular highest weight vectors, whose Cartan subalgebra is compact and thus has only imaginary roots. See example 24 of bosonic Gaussian states.
• Not semi-simple algebra. The Kähler properties of the manifold have to be checked case by case, as there is no general classification theory. In the physical important case of the Heisenberg algebra (regular coherent states), the presented construction still works. See example 25 of coherent states.
Our proofs and discussions only showed one direction, namely that φ being a highest weight vector (and annihilated by non-imaginary roots, if there are any) implies that M φ is Kähler. According to [82,85], the opposite is also true, i.e., group theoretic coherent states are Kähler if and only if they constructed from such a state φ. According to [82], variational families M φ that are Kähler also satisfy the conditions of so called symmetric spaces [86].

Co-adjoint orbits
In the context of group theoretic coherent states, one often transforms the problem of describing the geometry of M φ to a real submanifold O φ of the dual Lie algebra g * , which is called co-adjoint orbit. Using this language is particularly useful when we can establish an isomorphism M φ O φ , i.e., when they are equivalent manifolds. Definition 6. Given a non-zero state |φ ∈ H and M ∈ G, we define the dual Lie algebra element β M ∈ g * as This gives rise to its coadjoint orbit as the submanifold where we can compute (β M ) i = (β 1 ) j (Ad M ) j i . The motivation for introducing the coadjoint orbit is that O φ and M φ will turn out to coincide under certain conditions, including the important case where M φ is of Kähler type. Analogous to H φ and h φ , we define which are the stabilizer subgroup of the dual vector β 1 = i φ|Ξ i |φ / φ|φ and the associated Lie algebra. In other words, S φ behaves to the orbit O φ just like H φ to M φ .
Proposition 15. We always have h φ ⊂ s φ and H φ ⊂ S φ . If they are equal, O φ and M φ are isomorphic, i.e., Consequently, we also have h φ ⊂ s φ .
Only if H φ = S φ , we will have the equality while M φ is in general larger than O φ , i.e., there may be distinct states U(M ) |ψ ∈ M φ for which β M are the same. The following proposition allows the efficient computation of ω ij from β 1 .
It is non-degenerate if and only if h φ = s φ or, equiva- Proof. We compute explicitly Degeneracy of ω ij on tangent space means that there is a non-zero A i |V i = 0 with A i ω ij = 0. Recalling implies such A i cannot exist if and only if h φ = s φ .
The conclusion of this proposition is that ω ij is only non-degenerate if M φ = O φ and vice versa. It is wellknown in the theory of Lie groups [87] that any coadjoint orbit comes naturally equipped with a non-degenerate symplectic form and here we see that it agrees with the one on Example 26. Let us consider the example of SU(2) for the spin-1 2 and spin-1 representation a final time. The co-adjoint representation is isomorphic to the spin-1 representation (just like the adjoint) and can be understood as real 3-by-3 rotation matrices acting on real dual vectors (β 1 ) k . Therefore, all co-adjoint orbits for β 1 = 0 are 2-spheres, while the orbit β 1 = 0 is the single point {0}. Spin-1 2 . We recall that there is only a single M φ = P(C 2 ), so that we can just compute β 1 for the representative |φ ≡ (1, 0). This gives the dual vector β 1 ≡ (0, 0, 1 2 ). Unsurprisingly, we have the orbit O φ M φ S 2 .

Numerical implementation
The goal of this section is to explain how any class of coherent states, be it Kähler or non-Kähler, allows for an efficient implementation of real and imaginary time evolution. Here, we use the Lagrangian action for real time dynamics. The key advantage of coherent states lies in the fact that we can use the Lie algebra to identify the different tangent spaces, such that symplectic form Ω ij and metric G ij does not need to be evaluated at every step. This provides a significant computational advantage compared to a naive implementation without taking the natural group structure into account.
In practice, this reduces the calculation of real and imaginary time evolution tremendously. Instead of needing to compute metric and symplectic form at every step with respect to given coordinates, we can parametrize states by matrices M and use the identification from definition 5 to evaluate G ij or Ω ij once based on proposition 13. This gives the following dynamics.

C. Generalized Gaussian states
In the previous section we have considered manifolds made up of states of the form U(M ) |φ , where, as in definition 4, U is a unitary representation of the Lie group G, M is an element of G and |φ is a chosen reference state. We have seen that the geometric properties of the manifold defined in this way crucially depend on the choice of |φ . Indeed, if |φ is chosen as a highest weight vector of the representation then the Kähler property of the manifold follows naturally from the group structure. On the other hand if |φ is not a highest weight vector, then the Kähler property is not guaranteed.
We have also seen in example 24 that the group of Gaussian unitaries, that is unitaries that can be written as the exponential of anti-Hermitian operators at most quadratic in the linear operatorsξ a , gives a representation of the Lie groups ISp(2N, R) for bosons and O(2N, R) for fermions.
We will now consider in more detail the case of this representation, while we analyse different possible choices for the reference state |φ . As warm up, we will take the simple case of a single bosonic mode. Here we will already see all the different cases that can emerge. Then we will move to the case of an arbitrary number of bosonic and fermionic modes, where we will define the generalised Gaussian states first introduced in [30] highlighting how they fit in the present setting.

Warm up examples
Let us consider a single bosonic mode defined by the creation and annihilation operatorsb † andb or the quadraturesq = 1 √ 2 (b † +b) andp = i √ 2 (b † −b). As discussed in example 24, in this case the Gaussian algebra is spanned by the Cartan subalgebra elementĤ = i(1 + 2b †b ) and by the elements E + =bb and E − =b †b † , corresponding to the one positive and one negative root of the algebra. Then, the real basis (310) of the algebra is represented by the operatorŝ where we followed the conventions of example 13. Indeed, the most general Gaussian squeezing operator in such system can be then written as We then consider the manifold of group theoretic coherent states of the form |ψ(x, y, z) = U(x, y, z) |φ .
As discussed in the context of defintion 5, for states of this form there exists a natural isomorphism between the tangent spaces at different points, which are all equivalent to a subspace of the associated Lie algebra. It is thus sufficient to consider the tangent space at the point |ψ(0, 0, 0) = |φ . At this point, tangent space is spanned by the vectors which also coincide with the ones introduced in (292). The form of these states, after applying the projector Q, depends on the specific choice of |φ . We will now review some possible choices, showing that if we choose a highest weight state of the representation (first case) we will obtain a Kähler manifold, while if we do not (subsequent cases) we can construct non-Kähler manifolds of the different types discussed in the previous sections.
Kähler (Gaussian) case |φ = |0 . As already discussed extensively, choosing |φ as the Fock vacuum leads to the manifold of one-mode bosonic squeezed states. We have seen that these form a Kähler manifold, as can be expected in the light of what discussed in Section V B. Indeed, |0 is the highest weight state of the representation of Sp(2, R) we are using, so the resulting manifold is Kähler by Proposition 14. More concretely, we see that the last vector in (348) will be proportional to |φ and thus will vanish once the projector Q φ is applied. We are then left with a two dimensional tangent space, spanned by the first two which form a Kähler pair. Overall we indeed have Note that, as we are only considering squeezing and not displacements, a similar behaviour will appear also for |φ = |1 . This is indeed the the highest weight state of the odd sector of the representation, as discussed in example 24.
Non-Kähler non-degenerate case |φ = |2 . If we choose |φ as a Fock state |n with n ≥ 2 then we will similarly have that the the last vector in (348) vanishes after applying Q φ . The remaining two vectors will however not form a conjugate pair. For n = 2, we find Through simple calculations one can see that in this case the symplectic form ω will be invertible, but the complex structure J will, in the basis {|V 1 , |V 2 }, have the form which clearly does not square to −1.
Non-Kähler degenerate case |φ = 1 √ 2 (|0 + |2 ). If we choose a |φ that is not an eigenstate ofn, for example a superposition of two different Fock states, then it will not be a weight state of the representation. In this case none of the vectors in (348) will vanish after applying Q φ . Therefore, we will have a three dimensional tangent space. Being odd-dimensional, it cannot admit an invertible symplectic form ω. In particular, for |φ = 1 √ 2 (|0 + |2 ) we have This leads to the complex structure from which we have that J 2 has eigenvalues 0, −1 and −1. We see therefore that the manifold is not naturally Kähler, however if we invert ω only on the twodimensional subspace where this is possible (i.e., with the pseudo-inverse), as discussed in section III D, then the resulting manifold is Kähler.
Variable reference state |φ(α, Γ, z) = e iαn 2 |Γ, z . Finally, we can also consider the case where |φ is not a fixed state, but rather depends itself on some parameters that will then be part of the total parameter set of the manifold. In this case, different instances of the possibilities discussed above may occur at different points of the manifold. In particular, let us consider for |φ the states obtained by applying the unitary e iαn 2 , depending on the single real parameter α, to the family of Gaussian states, parametrized by |Γ, z according to the conventions of section V A. Here, at a generic point of the manifold, varying the parameter α will lead to a single independent unpaired tangent vector Q φ ∂ α e iαn 2 |Γ, z , which will necessarily make the tangent space non-Kähler. However at the special points where J = 1 and z = 0, that is |φ = |0 for any α, we have that the tangent space will be isomorphic to the one of Gaussian states, that is Kähler.

Definition
As seen in the previous section, choosing variable reference state |φ(x) can lead to manifolds with rather elaborate geometric structures. In this section, we will introduce a set of states that can be considered as a generalization of such example to the case of an arbitrary number of bosonic and fermionic modes. These states were first introduced in [30] as a possible generalizations of Gaussian states as a variational ansatz in many-body physics. The understanding of the properties of such nontrivial state manifolds is one of the main motivations of the present work.
We consider a Hilbert space containing both bosonic and fermionic degrees of freedom, i.e., H = H b ⊗ H f where H b is the bosonic Fock space of N b bosonic modes and H f is the fermionic Fock space of N f fermionic modes. We refer to the classical phase spaces V b R 2N b and V f R 2N f . On this space, we fix a basis of bosonic and fermionic linear operatorŝ which also determine the number operatorŝ for all i = 1, . . . , N b/f . Gaussian states of systems with bosonic and fermionic degrees of freedom are defined as the tensor product and are thus unable to capture correlations between bosons and fermions. This is an important drawback when studying correlated boson-fermion mixtures, which generalized Gaussian states are able to overcome. Gaussian transformations on the mixed Hilbert space H = H b ⊗ H f are defined as the representation of the group , such that U b and U f are the representations of respectively bosonic and fermionic Gaussian unitaries defined in section V A. We now consider states of the form U(M b , M f ) |φ(x) for variable choices of the state |φ(x) . More precisely, we consider the non-Gaussian reference states Here, U NG is a non-Gaussian unitary, acting on the full Hilbert space H = H b ⊗ H f , parametrized by the real parameters α and |Γ b , z b , Γ f is Gaussian. There are two important choices for the unitary U NG , depending on the type of correlation between the bosonic and fermionic sector one wishes to introduce, where the indices i, j run over the bosonic modes and i,j run over the fermionic ones. The resulting family of non-Gaussian states is then given by where we have M = (M b , z b , M f ), α = (α b , α f , α bf ) and Γ = (Γ b , z b , Γ f ). Note that this parametrizations certainly contains redundancies and careful analysis of the resulting family is desirable. This choice of non-Gaussian unitaries is motivated by the fact that, while going beyond the space of Gaussian transformation, they still allow for efficient computations. They satisfy the property U † NGξ a U NG = (U G ) a bξ b , where U G is a Gaussian unitary combined with a linear transformation of theξ a . Moreover, these states can be understood as true generalizations of Gaussian states in the sense of a generalized Wick's theorem. As discussed in [30], the evaluation of n-function follows from a quadratic generating function, just like in Wick's theorem. However, what makes them truly different from Gaussian states is that this generating functional is different for different n, such that more interesting correlation structures (such as bosonfermion correlations, but even within one sector) can be captured. For this reason any n-point function can be efficiently computed for the states introduced in this section, generalizing the property that in the case of Gaussian states follows from Wick's theorem. For this property to hold, the unitarity of U NG is essential. With this in mind, the parameters α should be considered real and not extended to the complex plane.
Therefore, the previous considerations about Kähler manifolds apply directly. Clearly, the manifold decomposes into equivalence classes (sheets) of states U |φ(x) related by Gaussian transformations. In particular, the sheets constructed this way will only be Kähler where the reference state |φ(x) is the highest weight state of the representation, i.e., the vacuum state |Γ b , z b , Γ f = |0 b ⊗|0 f , while the overall manifold (collection of sheets) will not be Kähler. Consequently, these states provide a natural playground for the concepts and methods introduced in this paper and vice versa a rigorous understanding of variational methods for non-Kähler manifolds is required to use generalized Gaussian states in practice.

VI. SUMMARY AND DISCUSSION
We have presented a systematic geometric framework to use variational methods for the study of closed quantum systems. Our results build on extensive previous work ranging from the geometric formulation of the time dependent variational principle [20] by Saraceno and Kramer to the formulation of group theoretic coherent states [24,25] due to Gilmore and Perelomov. Our approach places a particular emphasis on the restricted Kähler structures of a given variational family, which enabled us to classify families based on the spectrum of the restricted linear complex structure J . Only for J 2 = −1, the variational family is a so called Kähler manifold, which leads to many desirable properties. Most importantly, the different variational principles (Lagrangian, McLachlan, Dirac-Frenkel) all agree.
Geometric formulations of quantum theory have a long tradition [19,[21][22][23] in the mathematical physics community, but their usage for everyday applications in quantum many body physics and quantum optics has been limited for several reasons. Studying real and imaginary time evolution on high dimensional variational families requires extensive computational resources, which were not available thirty years ago. While it was often sufficient in the past to do a first order calculation based on standard families, such as coherent states, more complex models and more complicated physical questions (often inspired by experiments) often require new approaches, such as larger variational families. An important example are strongly correlated boson-fermion mixtures [31], whose properties have been successfully explored by combining generalized Gaussian states with exact diagonalization.
In the introduction, we formulated two main aims of the present paper, namely (a) to present a complete formulation of variational methods (covering real and imaginary time evolution together with approximating energy spectra and spectral functions) which is (b) accessible to practitioners and mathematical physicists alike. We addressed these goals by presenting an intuitive formalism based on defining tangent vectors |V µ , such that the vector fields X µ and F µ , the bilinear forms g µ and ω µν , and linear maps J µ ν and K µ ν can be covariantly expressed as tensors with indices, where we carefully distinguish between upper and lower indices. This formalism is based on Roger Penrose's abstract index notation, which provides an intuitive, yet rigorous approach of writing tensors. We believe that this formalism is particularly suitable for practicioners who need to seamlessly switch between abstract equations (without reference to a concrete basis) and concrete calculations (with respect to a fixed basis), whose numerical evaluation is based on storing tensor components in arrays. Finally, we provided a large number of examples to illustrate the discussed methods in familiar scenarios.
Kähler manifolds play a key role in our treatment, which is well motivated by the large number of desirable properties that they exhibit. On the other hand, there are good reasons to go beyond Kähler manifolds to construct new variational families, such as generalized Gaussian states or group theoretic coherent states (with |φ not being of highest weight). Our in-depth discussion highlights the most important differences between Kähler and non-Kähler manifolds with regards to variational methods. Moreover, we demonstrate how the orthogonal pseudo-inverse Ω provides a natural inverse ω to implement real time evolution, even when the restricted symplectic form ω is degenerate. We argue that Kähler and non-Kähler manifolds behave exactly the same with regards to imaginary time evolution and largely the same if we apply the Lagrangian action principle to compute real time evolution.
Our applications focus on coherent and Gaussian states as well as their generalizations (group theoretic coherent states, generalized Gaussian states). The reason for this is reviewed in section V B 5, where we explain that group theoretic coherent states of the form |ψ(M ) = U(M ) |φ with [81,82], our treatment emphasizes their geometric structures with respect to Kähler geometry and their utilization as variational families. In particular, we present a simple proof showing under which conditions group theoretic coherent states form Kähler manifolds.
Another motivation of the present paper is to lay the theoretical foundations for the exploration of generalized Gaussian states as proposed in [30] and reviewed in section V C. They present a new approach for the variational study of various systems, ranging from boson-fermion mixtures in Holstein models to Su-Schrieer-Heeger models and Kondo models. While these states have already been used in several studies [30,31], their geometric and mathematical structures have been largely unexplored. They form manifest non-Kähler manifolds and are constructed by applying certain non-Gaussian transformation U NG to bosonic and fermionic Gaussian states. This makes our formalism particularly suited for the avenue of exploring and classifying the mathematical properties of generalized Gaussian states, which may prove vital for their utilization as variational families.
We expect that the systematic application of variational methods to existing and newly proposed models in many body quantum problems will contribute to a better understanding of the involved physics with the potential of making new predictions relevant for experimental studies. Even revisiting known models with enlarged variational families can reveal new properties. For example, moving from coherent to the larger family of bosonic Gaussian states revealed recently that the ground state of trapped Bose-Einstein condensates with attractive s-wave interaction exhibits features of a squeezed state [89,90].
The present paper focused exclusively on variational families of pure quantum states, used for the study of closed quantum systems, i.e., at zero temperature. A natural extension of our formalism would be to also incorporate open quantum systems by allowing mixed states within our variational family. The approximate ground states ψ 0 minimizing the energy will be replaced by density operator ρ 0 minimizing the free energy. Nonequilibrium phenomena are often treated in their Markovian approximation, where time evolution is governed by master equations of the Gorini-Kossakowski-Sudarshan-Lindblad form [91,92]. Using this to derive meaningful variational equations is an important challenge, which has been only partially accomplished in the context of specific variational families [93,94]. We believe that the geometric perspective laid out in the current manuscript may also be helpful for developing variational methods to study the dynamics of open quantum systems. X µ real time evolutation vector field F µ imaginary time evolution vector field P µ ν projector onto conservation laws subspace X µ , F µ conservation laws preserving vector fields M manifolds of constant conserved quantities |Ψ = e κ+iϕ |ψ family with variable phase/normalization K µ ν linearized time evolution flow λ = ±iω eigenvalues of K µ ν E µ (λ) right-eigenvectors of K µ ν Eµ(λ) left-eigenvectors of K µ ν A(ω) spectral function V classical bosonic/fermionic phasespace a, b, c, d classical phase space indices g ab , ω ab , J a b Gaussian Kähler structures z a , C ab 2 one-point and two-point function Γ ab bosonic/fermionic covariance matrix C a 1 ...an n n-point function K representation of generator K S(M ) Bogoliubov transformation D(z) displacement transformation U(M, z) Gaussian transformation ∆ = −JJ0 relative covariance matrix |Va ∈ D (Γ,z) displacement tangent vector |V ab ∈ S (Γ,z) squeezing tangent vector G, g Lie group and Lie algebra i, j, k, l Lie algebra indices Ξi Lie algebra basiŝ Ξi operator representation of Lie algebra Kij Killing form |Vi Lie algebra induced tangent space basis φ, |φ reference state and state vector M φ ⊂ P(H) group theoretic coherent states H φ , h φ stabilizer group and algebra of φ h real Cartan subalgebra I, J, K, L Cartan subalgebra indices HI Cartan subalgebra basis vα weight spaces with weights α (βM ) k ∈ g * expectation value ofΞ k O φ co-adjoint orbit of φ S φ , s φ stabilizer group and algebra of β 1

Abstract index notation
Throughout this paper, all equations containing indices follow the conventions of abstract index notation. This formalism is commonly used in the research field of general relativity and gravity, where differential geometry plays an important role, but we believe that it is also of great benefit when studying the geometry of variational manifolds.
The formalism is suitable to conveniently keep track of tensors built on a vector space. Given a finite dimensional real vector space V with dual V * , a (r, s)-tensor T is a linear map T : V * × · · · × V * × V × · · · × V → R . (A1) In particular, a (1, 0)-tensor is a vector, a (0, 1)-tensor is a dual vector and a (1, 1)-tensor is a linear map. To keep track of the type of tensor, abstract index notation refers to the (r, s)-tensor T as T i1···ir j1···js , i.e., we assign r upper indices and s lower indices. Typically, we choose the indices from some alphabet to indicate which vector space, we are referring to. For tangent space T ψ M of some variational manifold M, we use Greek letters µ, ν, γ, δ, while we use Latin letters a, b, c, d to refer to the classical phase space V of a bosonic or fermionic system. In the context of group theoretic states, we use i, j, k, l to refer to the Lie algebra and I, J, K, L to refer to its Cartan subalgebra.
The key advantage of abstract index notation in the context of variational manifolds is that it helps us to keep track of what types of tensors, we are dealing with and which contractions are allowed. Apart from vectors X µ and dual vectors w µ , we are mostly dealing with tensors that have two indices, namely linear maps J µ ν , bilinear forms g µν and dual bilinear forms Ω µν In the present paper, we often deal with linear maps and bilinear form, i.e., tensors that have two indices. They are naturally represented as matrices, in particular, for numerical evaluation. For convenience, we will also use the notation, where tensors with suppressed indices are implied to be contracted, just as standard matrix multiplication works. Obviously, this means that only such expressions are allowed where the adjacent indices are given by one upper and one lower index.

Special tensors and tensor operations
In many area of physics, it is not necessary to distinguish much between vectors and covectors.
Identity. Every vector space V comes with the canonical identity map δ µ ν satisfying δ µ ν X µ = X µ . Note that the notation 1 µ ν would also be consistent, but we stayed with the commonly used Kronecker delta. There does not exist a canonical analogue as bilinear form, e.g., a form δ µν or δ µν which only make sense with respect to a specific basis and are therefore not canonical, but rather a specific choice, such as a metric g µν .
Transformation rules. An invertible linear map M a b : V → V of the vector space V acts on a general (r, s)-tensor T i1···ir j1···js and transforms it to M a1 c1 · · · M ar cr (M −1 ) d1 b1 · · · (M −1 ) ds bs T c1···cr d1···ds (A2) In particular, a vector X a transforms as The determinant of a bilinear form s ab or S ab is ill defined, unless we a reference object, such as a metric g ab or G ab . Then, we can compute the determinant of the matrix of the linear maps S ac g cb or G ac s cb .
Trace. The trace tr(M ) = M a a is only defined for a linear map, not for bilinear forms S ab or s ab , unless we again have a reference object, such as a metric.
Eigenvalues. Without additional structures, we can only defined eigenvalues for a linear map M a b , where an eigenvalue λ associated to an eigenvector X a satisfies This is well-known from linear algebra. A bilinear form X ab does not have intrinsic eigenvalues, but we can compute its eigenvalues relative to another bilinear form. Given a bilinear form s ab and a metric G ab or symplectic form Ω ab , we can define the metric or symplectic eigenvalues as the regular eigenvalues of the linear map G ac s cb or Ω ac s cb , respectively. Transpose. The transpose of a linear map M a b : V → V is the map (M ) a b : V * → V * . We have the relation (M ) a b = M b a , which means that the two represent the same tensor and typically one does not distinguish between the two in abstract index notation. However, for our shorthand notation, it is important to keep the order of indices in right order. From the perspective of abstract index notation, there is not really much point to use the transpose operation, but we will still write the respective expressions for convenience, so that they can be easily converted to matrix expressions, e.g., for numerical implementations.
Gradient. Given a function f (x) on some manifold M, its gradient (df ) µ = ∂ µ f is field of covectors, also known as 1-form. This means that the gradient alone does not define a tangent space direction, e.g., to move in the direction of steepest ascent. Indeed, only if we have a metric G µν , we can define typical gradient vector field F µ = G µν (∂ ν f ). The reason is that the gradient df as dual vector encodes the linearized change df (X) = df µ X µ of the function f , when performing a step in the direction X µ . Clearly, by increasing our step size, we can make this change arbitrarily large, so there is no "steepest" direction. Only if we have an absolute measure of our step size, e.g., a norm X = X µ g µν X ν induced by an inner product g µν , we can find a unique direction, for which a step of fixed size maximizes the change.
plex structure is represented by the block matrix with 0 < c i < 1. This standard form induces the decomposition of T ψ M into the three orthogonal parts where T ψ M is the largest Kähler subspace and T ψ M is the largest space on which J and ω are invertible.
Proof. We focus on a single tangent space T ψ M ⊂ H and refer to the Kähler structures on H, rather than the restricted ones on T ψ M, as (g, ω, J). To shorten notation, we define A := T ψ M and B as its orthogonal complement in H with respect to g, so that H = A ⊕ B. We will refer to the restricted Kähler structures on A or B by (g A , ω A , J A ) and (g B , ω B , J B ), respectively. The relation J = Gω = −Ωg implies g = −ωJ or, equivalently, g(v, w) = −ω(v, Jw), and also g(Jv, Jw) = g(v, w) for all v, w ∈ H. From here, g(Jv, w) = −g(v, Jw) follows and we can derive for a, a ∈ A g A (J A a, a ) = g(Ja, a ) = −g A (a, J A a ) , which implies that J A is antisymmetric with respect to g A and thus is diagonalizable, has either vanishing or purely imaginary eigenvalues with the latter appearing in pairs ±c i . Furthermore, we can always choose an orthonormal basis, such that g A = 1 and J A is represented by (B1). Next, we show that c i ∈ (0, 1]. We define the orthogonal projectors P We write J 2 − 1 = 0 in blocks to find We consider an eigenvector a ∈ A of J A with J A a = ica for non-zero c, which implies J 2 A a = −c 2 a. We compute g(a, a) = g(J A a + J BA a, J A a + J BA a) = g(J A a, J A a) + g(J BA a, J BA a) ≥ g(J A a, J A a) = −g(a, J 2 A a) = c 2 g(a, a) where we used that A and B are orthogonal which eliminates crossing terms. This implies the inequality c 2 ≤ 1 and thus, we can choose c i ∈ (0, 1] as in (66).
Proposition 3. The Kähler property is equivalent to requiring that T ψ M is not just a real, but also a complex subspace, i.e., for all |X ∈ T ψ M, we also have i|X ∈ T ψ M. Therefore, the multiplication by i commutes with the projector P ψ , i.e., P ψ i = iP ψ and P ψ is complex-linear.
Proof. We want to show that J 2 A = −1 A implies that for all a ∈ A, we also have ia = Ja ∈ A. Therefore, we need to show that Ja = J A a, which is equivalent to J BA = 0. For arbitrary a ∈ A, we compute g(J BA a, J BA a) = g(Ja, J BA a) = −g(a, JJ BA a) = −g(a, J AB J BA a) .
This expression vanishes if J 2 A = −1 A , because in that case J AB J BA = J 2 A − 1 A = 0 follows from the first block in (B5). Since g is non-degenerate, this implies J BA = 0. Similarly, we can use the last block in (B5) to conclude J 2 B = −1 B , which implies J BA = 0. With vanishing J AB and J BA , J is block diagonal and commutes with the projectors. In the language of complex vector spaces, this implies that P ψ i = iP ψ . Proposition 8. Given a variational manifold M we define (according to the Lagrangian action principle) the free projected real time evolution |ψ(t) as governed by the free HamiltonianĤ 0 and the perturbed projected real time evolution |ψ (t) as governed by the perturbed Hamilto-nianĤ (t) =Ĥ 0 + Â (t), both with the same initial state |ψ(0) . Then, the propagated perturbation, defined according to (148), is given by where dΦ t is the linearized free evolution flow.
Proof. Let us define the perturbed evolution flow Φ t as the map that sends the coordinates of an initial state x µ (0) to the coordinates x µ (t) of the state time evolved under the projected perturbed real time evolution. It is governed by where X 0 and X A are the evolution vector fields associated to the HamiltoniansĤ 0 andÂ respectively. We define the free evolution flow Φ 0 t analogously by just setting = 0 in the previous expressions. Let us now define the interaction picture flow Φ t = Φ 0 −t •Φ t . It has the useful property that its time evolution only depends on the perturbing vector field: • We restrictÑJ onM to M to find (ÑJ ) µ νσ = Jλ σ ∂λJ µ ν − Jλ ν ∂λJ µ σ + J µλ (∂ ν Jλ σ − ∂ σ Jλ ν ) (C4) which is not obviously equal to (N J ) µ νσ due to the contraction overλ, which takes the full manifold into account. However, our previous considerations showed thatJ = J ⊕ J . This implies that Jλ µ = J λ µ , which proves the equality. Consequently NJ = 0 implies N J = 0.
We therefore conclude that any submanifold M of a Kähler manifoldM with J 2 = −1 everywhere is again a Kähler manifold. Note that this implies in particular that M is also a complex and a symplectic manifold.
Appendix D: Normalized states as principal bundle We introduced projective Hilbert space P(H) as the space of distinguishable quantum states, where normalization and phases of Hilbert space vectors can be ignored. However, in practice we usually parametrizing a variational manifold M ⊂ P(H) by a set normalized states |ψ(x) which depend on some real parameters x i . It is therefore useful to introduce the manifold of normalized states N (H) = |ψ ∈ H ψ|ψ = 1 . (D1) This manifold will play an important when we are interested in relative phases between states. Mathematically, N (H) is a principal fiber bundle over the base manifold P(H). Given a quantum state ψ ∈ P(H), we a corresponding fiber e iϕ |ψ ∈ N (H) of normalized Hilbert state vector representing this state. We also have a natural U(1) group action onto such fibers given by multiplication with a complex phase, i.e., |ψ → e iϕ |ψ with e iϕ ∈ U(1).
We can now ask if we can use any structures of the Hilbert space to define a natural notion of parallel transport, i.e., how to choose the complex phases when changing the quantum state continuously. This question is well-studied in the context of gauge theories and amounts to choosing a natural notion of horizontal tangen spaces, which encode how to move naturally through the principal fiber bundle. Interestingly, N (H) is equipped with a natural notion of moving horizontally, namely by requiring that a horizontal curve |ψ(t) satisfies i.e., we require that the tangent vector d dt |ψ(t) is always orthogonal to the state |ψ(t) . The tangent space of normalized states is given by T |ψ N (H) = |ϕ ∈ H Re ψ|φ = 0 (D3) Locally, we can decompose this space into where the former the former is vertical space along the fiber and the latter is our natural choice of horizontal subspace. We can understand the local choice of complex phase e iϕ as pure gauge and U(1) as the corresponding gauge group. This implies that our description of normalized states N (H) is equivalent to electromagnetism on the base manifold P(H). In particular, we can compute the gauge field A µ and its field strength tensor F µν .
We can compute the gauge field A i as A µ = Im ψ(x)|∂ µ |ψ(x) ≡ (0, sin 2 θ 2 ) . (D6) As differential form, we therefore have A = sin 2 θ 2 dφ. We find the field strength as its differential We can compute the change of phase ∆ϕ, also called holonomy, if we move horizontally in N (H), such that our path describes a circle at constant latitude θ when projected onto the Bloch sphere. We find ∆ϕ = 2π 0 dφA = 2π sin 2 θ 2 . (D8) We can also use Hamiltonian evolution bŷ to compute the time evolution. The Hamiltonian vector field is given byẋ such that the solution of the equation of motions is just φ(t) = −ωt. The quantum state will return to its original state at t = π ω , where φ(t) = −2π. However, we will pick up a relative phase given by the integral where we used E = E 0 + ω cos θ.
In general, we conclude that the change of complex phase after returning to the same quantum state through time evolution is given by where E(x) = E(x 0 ) is constant for time-independent Hamiltonians. If we go to variational families M ⊂ P(H), we can still use (D12) to compute the holonomy associated to projected time evolution. This is important if we want to compute spectral functions from real time evolution on M rather than via the linearization around a stationary point.
In practice, we therefore see that the only required ad-ditional structure on a variational family M is the computation of the gauge field A µ . Once this is computed, we can derive the change of relative phases from integration over A µ and the energy expectation value, which is constant for time-independent Hamiltonians. Once, we have taken the relative phase in the time evolution into account, we can also use it to compute inner products between different state vectors, while ensuring the correct complex phase.