Gaussian time-dependent variational principle for the finite-temperature anharmonic lattice dynamics

The anharmonic lattice is a representative example of an interacting bosonic many-body system. The self-consistent harmonic approximation has proven versatile for the study of the equilibrium properties of anharmonic lattices. However, the study of dynamical properties therewithin resorts to an ansatz, whose validity has not yet been theoretically proven. Here, we apply the time-dependent variational principle, a recently emerging useful tool for studying the dynamic properties of interacting many-body systems, to the anharmonic lattice Hamiltonian at finite temperature using the Gaussian states as the variational manifold. We derive an analytic formula for the position-position correlation function and the phonon self-energy, proving the dynamical ansatz of the self-consistent harmonic approximation. We establish a fruitful connection between time-dependent variational principle and the anharmonic lattice Hamiltonian, providing insights in both fields. Our work expands the range of applicability of time-dependent variational principle to first-principles lattice Hamiltonians and lays the groundwork for the study of dynamical properties of the anharmonic lattice using a fully variational framework.

Introduction -Variational methods form the basis of our understanding of quantum mechanical many-body systems. In a variational method, the wavefunctions or density matrices of a system are parametrized by a set of parameters whose size is much smaller than the dimension of the Hilbert space. Static and time-dependent [1-3] variational methods are being actively used to study interacting many-body model Hamiltonians [4][5][6][7][8][9][10][11][12][13].
However, SCHA is limited in that one needs to resort to a specific ansatz to study the dynamical properties. It is known that the SCHA ansatz for the position-position Green function is correct in the static limit of zero frequency and the perturbative limit of weak anharmonicity [18]. However, the validity of the SCHA ansatz in the nonperturbative and dynamic regime [20,21,33,35], where the dynamical theory is most necessary, has not been theoretically justified.
In this Letter, we solve this important problem by applying the time-dependent variational principle (TDVP) with Gaussian variational states [7,13,36, 37] to the anharmonic lattice Hamiltonian at finite temperature.
Gaussian TDVP expands the static variational states of SCHA to states with nonzero momenta. We use the linearized time evolution to derive the position-position correlation function and prove the SCHA dynamical ansatz. We illustrate that the Gaussian TDVP is successful in describing the dynamics because it includes the 2-phonon states as true dynamical excitations. Our work connects the TDVP theory, whose application was mostly focused on model Hamiltonians for cold atoms, with anharmonic lattice dynamics and the SCHA method. Such connection gives fruitful results on both sides. In TDVP theory, the linearized time evolution and the projected Hamiltonian method [4,6,9] are two different ways to compute the excitation spectrum, whose superiority over the other varies across systems [4,7,13]. We use the anharmonic lattice model to show that only the linearized time evolution gives correct excitation energies in the perturbative limit. On the SCHA side, we illustrate ways to systematically expand the SCHA theory by leveraging recent developments of non-Gaussian TDVP [6,11,12].
Self-consistent harmonic approximation -We briefly review the key results of SCHA. Within the adiabatic Born-Oppenheimer approximation, the anharmonic lattice Hamiltonian iŝ H = N a=1ˆ p 2 a 2M a + V (ˆ r 1 , ...,ˆ r N ). (1) Here, a is the combined index for atoms and Cartesian directions, N = N atm ×d with N atm and d the numbers of the atoms and the spatial dimensions, respectively, M a the atomic mass,ˆ r a andˆ p a the position and momentum operators, and V the Born-Oppenheimer potential energy. We set = 1. In SCHA, the true thermal equilibrium state of the anharmonic Hamiltonian is approximated by that of a harmonic HamiltonianĤ (H) : (2) Since we study the dynamics around the SCHA equilibrium, we assume that the optimal harmonic potential V (H) is already found. The SCHA density matrix iŝ where β = 1/k B T is the inverse temperature. For later use, we define Â 0 ≡ Tr ρ 0Â .
Hereafter, we use the normal-mode representation, where the SCHA harmonic Hamiltonian becomeŝ with ω m the eigenvalue of the SCHA dynamical matrix, andr m andp m the normal-mode position and momentum operators. The anharmonic Hamiltonian [Eq.
(1)] can be written aŝ with V (r) = V (ˆ r) the potential energy in the normalmode representation.
Also, sinceρ 0 is a thermal state, we find p m 0 = 0, p mpn 0 = n m + 1 2 δ m,n , with n m = 1/(e βωm − 1) the occupation number. Gaussian time-dependent variational principle -Next, we discuss the general principles of Gaussian TDVP for a multimode bosonic system at finite temperature. We use the set of states obtained by applying a Gaussian unitary transformationÛ (x) to the SCHA density matrix as the variational manifold: Here, x is a real-valued vector that encodes all the variational parameters. We parametrize the Gaussian transformation asÛ whereD andŜ are the displacement and squeezing operators, respectively: where b mn ≡ 1/ 4(n m + n n + 1) if m = n 1/ 2(n m + n n + 1) if m = n , The variational parameters α m , β mn , and γ mn are complex numbers. The parameter β mn (γ mn ) is defined only for m ≤ n (m < n). Here, we assume for simplicity that ω m 's are nondegenerate and satisfy ω 1 < ω 2 < · · · < ω N . The total number of complex variational parameters is N 2 + N . In the linear response regime, degeneracy does not pose any theoretical difficulty: if modes m and n are degenerate, one just needs to exclude γ mn from the set of variational parameters. This exclusion is done because the infinitesimal transformation parametrized by γ mn does not changeρ 0 [38]. Each group of parameters describes a different type of excitation. Parameters α, β, and γ correspond to 1phonon excitations, 2-phonon excitations with two creations or two annihilations of phonons, and 2-phonon excitations with one creation and one annihilation, respectively.
The imaginary parts of the parameters generate dynamics. For example, Im α generates a finite atomic momentum through the displacement operator. The SCHA theory does not contain these imaginary parameters because the variational states are limited to the thermal state of a harmonic Hamiltonian. In contrast, Gaussian TDVP, which allows both the real and imaginary parts of the variational parameters to vary, naturally allows one to study the dynamics.
We define x, the vector of variational parameters as Sinceρ(x = 0) =ρ 0 is the variational solution that minimizes the SCHA free energy, x = 0 is a stationary point of the variational time evolution [38].
To apply TDVP to mixed states, we map the variational density matrices to wavefunctions by purification [11,39]. For each physical state in the number basis, we add an auxiliary state so that the purified wavefunction becomes where ⊗ denotes a tensor product, and |Φ + is a maximally entangled state [39] between the physical and the auxiliary modes [see Eq. (S15) and related discussions]. For the purified wavefunction, the expectation value of a physical operatorÂ 0 is whereÂ The variational time evolution is obtained by projecting the true dynamics of the wavefunction to the tangent space of the variational manifold. The tangent space is spanned by the tangent vectors, which at x = 0 are Using the variational linear response theory [13,38], one can show that the retarded correlation function G Here, the matrix Green function G(z) is defined as where K is the linearized time-evolution generator defined as with E(x) = Tr ρ(x)Ĥ . The symplectic form Ω is defined by By computing K and the corresponding matrix Green function G(z), one can find the physical correlation function G (R) AB (ω) using Eq. (19). Anharmonic lattice dynamics -Now, we study the dynamical properties of the anharmonic lattice Hamiltonian using Gaussian TDVP. First, the symplectic form is [38] with ⊕ the direct sum.
The three matrices correspond to the subspace spanned by the tangent vectors for the variation of α, β, and γ, respectively. In each matrix, the bases for the first (second) block of rows and columns are the tangent vectors for the real (imaginary) parts of the parameters.
For later use, we define P 1 , P 2+ , and P 2− as the projection operators to the bases of each of the three matrices. The subscripts 1, 2+, and 2− indicate the nature of the tangent vectors: 1-phonon excitations, 2-phonon excitations with two creations or two annihilations, and 2-phonon excitations with one creation and one annihilation. We also define the projection to the whole 2-phonon sector: P 2 = P 2+ + P 2− . Evaluating Eq. (21), we find that the time evolution generator K is the sum of the non-interacting part, 3-phonon interaction, and 4-phonon interaction (see Sec. S4 C of the Supplementary Material [38]): where Here, we defined the diagonal matrices: [ω ± ] mn,pq = (ω m ± ω n )δ mn,pq , B mn,pq =b mn (n m + n n + 1)δ mn,pq , C mn,pq = − c mn (n m − n n )δ mn,pq .
The implicit summation over a pair of mode indices m and n implies the constraint m ≤ n unless otherwise noted. We also defined the anharmonicity tensor The non-interacting part H (0) describes the free evolution of 1-and 2-phonon excitations in the SCHA Hamiltonian. The 3-phonon interaction V (3) couples the 1and 2-phonon excitations. The 4-phonon interaction V (4) couples the 2-phonon excitations to each other. Finally, we study the linear response of the anharmonic lattice and compute the position-position correlation function. First, we define the non-interacting Green function G (0) : From Eq. (25), one finds where and Next, we include the 4-phonon interaction V (4) . We define the partially interacting Green function G (4) (z): Since the 4-phonon interaction V (4) does not act on the 1-phonon sector, we find For the 2-phonon sector, we obtain the Dyson equation Finally, we study the fully interacting Green function G(z) by including the 3-phonon interaction V (3) . From the definitions of G and G (4) , we obtain the Dyson equation Here, we defined B + = B and B − = C. In Eq. (41), we omitted the direct sum of the zero matrix in the P 2 subspace for brevity.
From Eqs. (S1, S2), one finds that the matrix elements for the position operator is nonzero only for the variation of Re α: Then, from Eqs. (19,41), one can derive the Dyson equation for the interacting retarded position-position correlation function [38]: The self-energy is where W is a diagonal matrix defined as By recovering the mode indices and defining (47) In Eq. (47), the implicit summation over the mode indices is done without any constraints. Equation (47) and its derivation is the main result of this Letter. When transformed to the Cartesian representation, Eq. (47) becomes identical to the SCHA dynamical ansatz [Eq. (70) of Ref . 18]. We emphasize that we rigorously derived the phonon self-energy Π rr (z) using Gaussian TDVP. Our derivation theoretically proves the SCHA dynamical ansatz. The physical interpretation of the self-energy formula we obtained vary significantly from that of the SCHA dynamical ansatz. In Gaussian TDVP, the 2-phonon states are true dynamical excitations. However, in SCHA, the 2-phonon states do not have their own dynamics and appear only indirectly through the position dependence of the SCHA force constants. The presence of the dynamical 2-phonon excitations is the essential reason why Gaussian TDVP can describe dynamical properties while the SCHA theory cannot.
For example, the phonon lifetime is an important dynamical property of an anharmonic lattice. In Gaussian TDVP, the 1-phonon states acquire a finite lifetime by decaying to the continuum of 2-phonon states through the 3-phonon interaction. In contrast, in SCHA, there are no  continuum states to which the 1-phonon states can decay. Hence, in SCHA, the phonon lifetimes can only be described with a perturbative approximation [32] unless one resorts to an ansatz.
Discussion -A common alternative to the linearized time evolution is the projected Hamiltonian method [4,6,9]. There, the Hamiltonian is projected onto the tangent space of the variational manifold. Let us consider a single-mode anharmonic oscillator at T = 0, whose Hamiltonian iŝ (48) Here, λ is the perturbation strength. The SCHA variational Hamiltonian iŝ and the variational ground-state energy is ω 0 /2. In Table I we list the excitation energy, the difference of the ground-and first-excited state energy, computed using different methods [38]. Comparing the variational methods to the perturbation theory, we find that the linearized time evolution is correct in the perturbative limit λ → 0, while the projected Hamiltonian method is not. Since the SCHA dynamical ansatz is exact in the perturbative limit [18], this finding also holds for a general multimode anharmonic lattice at finite temperatures.
This difference occurs because the projected Hamiltonian method fails to describe the effect of virtual 3-and 4-phonon states. In Fig. 1, we show the two processes that appear in the time domain representation of the bubble diagram for the phonon self-energy. Figure 1(b) describes a process involving a 4-phonon state. Since the Gaussian projected Hamiltonian method completely neglects the 3-and 4-phonon excitations, it only includes the process described in Fig. 1(a), not that of Fig. 1(b). In contrast, in the linearized time evolution, the coupling of the 1-and 2-phonon states to virtual 3-and 4-phonon states is included by an additional term related to the derivative of the tangent vectors, which is neglected in the projected Hamiltonian method [13]. Thanks to this additional term, the linearized time evolution gives the correct perturbative limit, while the projected Hamiltonian cannot.
A promising future research direction based on our study is a rigorous, systematic expansion of the SCHA method to go beyond the harmonic approximation by using non-Gaussian variational transformations [6]. Also, the use of mixed fermionic and bosonic variational states [6,11,12] will allow the study of nontrivial electron-phonon correlation such as in phonon-mediated superconductivity or polarons in anharmonic lattices.
Recently, Monacelli and Mauri also reported a proof of the SCHA dynamical self-energy in an independent work [41]. While Ref.
[41] additionally presents a numerical algorithm to compute the correlation functions, our work focuses on the link between TDVP and SCHA. Also, while the proof for the finite-temperature case in Ref.
[41] is based on an analogy with the T =0 case, our proof uses purification to rigorously derive the finite-temperature equation of motion. The results of the two works are consistent when there is an overlap.
Conclusion -In summary, we developed a variational theory for the dynamical properties of anharmonic lattices using Gaussian TDVP, establishing a firm link between Gaussian TDVP and SCHA. We provided solid theoretical groundwork for the use of the SCHA dynamical ansatz in studying spectral properties. The presence of dynamical 2-phonon excitations in Gaussian TDVP was essential to obtain correct dynamics of the 1-phonon excitations. We compared the linearized time evolution and the projected Hamiltonian methods to find that only the former is correct in the perturbative limit. Our work establishes a useful connection between TDVP and SCHA, allowing further developments in both fields. In this section, we detail the physical meaning of the transformations and the tangent vectors by inspecting the infinitesimal transformation of the position and momentum operators. Using the definition of the operator transformation [Eq. (17) Hereafter, all derivatives with respect to the variational parameters are evaluated at x = 0 unless otherwise stated. The derivatives of the Gaussian transformation operator at x = 0 are We calculate how the position and momentum operators transform for an infinitesimal change of each variational parameter. For conciseness, we write the real and imaginary parts of the variational parameters as follows: First, for the displacement parameter α, the infinitesimal transformation of the position and momentum operators are For the real part of the squeezing parameters β and γ, the infinitesimal transformation ofr andp are = b mn (r m δ n,p +r n δ m,p ), (S13) = −c mn (r m δ n,p +r n δ m,p ).
From Eqs. (S9-S14), one can understand the role of each variational parameter. The real part of the displacement parameter, α r m , parametrizes the displacement of the position operator for mode m. These N degrees of freedom corresponds to the center position R in the SCHA harmonic Hamiltonian. The real parts of the squeezing parameters, β r mn and γ r mn , parametrize the change in the normal mode frequency and eigenvectors. Especially, γ r mn parametrizes the linear combination of the two eigenmodes m and n.
If modes m and n are nondegenerate, setting γ r mn = 0 mixes two modes with different frequencies, inducing a nontrivial transformation of the thermal density matrix. In contrast, if modes m and n are degenerate (i.e. ω m = ω n ), the linear combination parametrized by γ r mn is a gauge transformation that does not change the density matrix. Hence, it is justified to exclude γ r mn from the variational parameters when modes m and n are degenerate, as mentioned in the main text. From a theoretical point of view, including γ mn in the set of variational parameters for degenerate modes m and n makes the symplectic form [Eq. (22) of the main text] noninvertible and thus should be avoided [S1].
The imaginary parts of the Gaussian parameters generate dynamics of the variational states. The displacement parameter α i m parametrizes the generation of finite atomic momentum. The squeezing parameters β i mn and γ i mn parametrize the linear combination of the position coordinates with the momentum coordinates and vice versa.

S2. LINEAR RESPONSE FORMULATION OF TDVP AT FINITE TEMPERATURES
In this section, we derive and summarize the key results of the linear response formulation of TDVP at finite temperatures, following Ref. [S1].
In Eq. (15) of the main text, we mapped the variational density matrices to pure state wavefunctions by purification. The maximally entangled state |Φ + is defined as Φ + ∼ n1,··· ,n N |n 1 , · · · , n N ⊗ |n 1 , · · · , n N . (S15) Thanks to the unitarity ofÛ , the variational wavefunction |Ψ(x) is always normalized to unity. The original density matrix is recovered by taking a partial trace of the auxiliary system: Note that although |Ψ(0) is a purification of the thermal stateρ 0 of the harmonic HamiltonianĤ (H) , it is not a stationary state of the time evolution withĤ (H) : for any choice of the phase φ(t). In the second equality of Eq. (S17), we used the fact thatĤ (H) is diagonal in the eigenmode basis. Still, the corresponding density matrix that is obtained by taking the partial trace of the auxiliary system is time-independent. Hence, the time evolution of the purified wavefunction is not a true dynamics in the physical system. It is an auxiliary dynamics that occurs due to the non-uniqueness of the purification up to a unitary transformation at the auxiliary system. This artificial dynamics does not occur in our variational approach because we do not allow any variational degree of freedom to the auxiliary system. The time evolution of the variational wavefunction is obtained by projecting the change in the wavefunction to the tangent space of the variational manifold. The tangent space is spanned by the tangent vectors, which are the derivatives of the variational wavefunction orthogonalized to the original wavefunction. Formally, the tangent vectors are defined as (S18) whereQ(x) a projection operator:Q According to TDVP, the dynamics of the variational parameters can be described by a classical Hamilton equation of motion. To determine the equation of motion, we need the symplectic form and the derivatives of the energy expectation value E(x) = Ψ(x)|Ĥ|Ψ(x) [S1].
The symplectic form Ω µν (x) is the inverse of ω µν (x), which is twice the imaginary part of the inner product of the tangent vectors: We use Greek indices to denote the components of the real-valued vector x defined in Eq. (14) of the main text. We use Einstein's summation convention for repeated indices. According to the Lagrangian action principle, the equation of motion of the variational parameters is [S1, S2] We note that since the Gaussian variational manifold is a Kähler mainfold, the Lagrangian, McLachlan, and Dirac-Frenkel TDVP equations are all equivalent [S1]. Now, we illustrate how to compute dynamical and spectral properties using the linear response formulation of TDVP. As we are interested only in small changes of the wavefunction around the stationary state, we linearize the equation of motion Eq. (S22) around x = 0 to find [S1-S3] where the linearized time-evolution generator K is as shown in Eq. (21) of the main text. Here, we used (∂E/∂x ρ )| x=0 = 0 which is true because x = 0 is a stationary point. From now on, we denote ∂/∂x µ by ∂ µ . Also, we use Ω µρ to refer to Ω µρ (x = 0) unless otherwise noted. The solution of the linearized equation of motion is where dΦ(t) is the linearized free evolution flow defined as Let us consider a standard linear response setting, where an infinitesimal time-dependent perturbation is added to the Hamiltonian:Ĥ (t) =Ĥ + ϕ(t)Â. (S27) Here,Â is an arbitrary Hermitian operator in the Hilbert space of purified wavefunctions, ϕ(t) is a real-valued function, and is a real variable parametrizing the strength of the perturbation. We write the solution of the corresponding variational time evolution as |Ψ (t) ≡ |Ψ(x (t)) . The linear response of the variational parameter is defined as .

(S28)
According to Proposition 8 of Ref. [S1], δ A x µ (t) is given as where The linear response of the expectation value of an operatorB at time t is [S1] Now, we use the spectral decomposition of K to compute dΦ(t). One can decompose iK with eigenvalues λ l , eigenvectors E µ (λ l ) and dual eigenvectors E ν (λ l ) [S1]: The dual eigenvectors satisfy Then, the linearized free evolution flow becomes Using Eq. (S31) and Eq. (S34), we find By taking the Fourier transform of Eq. (S35), we find The retarded correlation function G From Eqs. (S37) and (S36), we find Then, using the definition of the matrix Green function [Eq. (20)], we find Eq. (19) of the main text.

S3. DERIVATION OF THE SYMPLECTIC FORM
In this section, we calculate the overlap of the tangent vectors to calculate the metric and the symplectic form. From the definition of the tangent vectors [Eq. (18)], one finds To evaluate Eq. (S39), we use the derivatives of the variational transformation, Eqs. (S2-S7). Now, we compute the overlap. First, since the thermal expectation value of an operator containing uneven numbers of creation and annihilation operators is zero, one finds and Next, we calculate the nonzero inner products. First, for two displacement parameters α m and α n , we find and Next, we consider the tangent vectors of the squeezing parameters β r mn and β r pq . Note that δ m,p δ n,q + δ m,q δ n,p = δ mn,pq if m = n 2δ mn,pq if m = n (S45) holds since m ≤ n and p ≤ q. Then, using Eqs. (S4, S5), we find =b mn b pq (δ m,p δ n,q + δ m,q δ n,p )(2n m n n + n m + n n + 1) = 2n m n n + n m + n n + 1 2(n m + n n + 1) δ mn,pq , and Finally, for the tangent vectors of the squeezing parameters γ r mn and γ r pq , m < n and p < q holds by definition. Thus, we find =c mn c pq δ m,p δ n,q (2n m n n + n m + n n ) = 2n m n n + n m + n n 2(n m − n n ) δ mn,pq , and The only inner products with nonzero imaginary parts are those in Eqs. (S43, S47, S49) and their complex conjugates. Using this result and the definition of the symplectic form [Eq. (22)], one obtains Eq. (23) of the main text.

S4. DERIVATIVES OF THE ENERGY
In this section, we calculate the first and second derivatives of the energy expectation value with respect to the variational parameters. In this section, all derivatives are evaluated at x = 0 unless otherwise noted. (S52) Here, ρ 0 (r) is the diagonal part ofρ 0 in the normal mode position basis [S4]: In the third equality of Eq. (S52), we used a partial integration with respect to r m [see also Eqs. (C1-C3) of Ref. [S5]]. In addition, using

B. First derivatives
Now, we compute the first derivatives of the energy expectation value and show that the SCHA solution is also the stationary state of the Gaussian TDVP. By settingÔ =Ĥ in Eq. (S1) and taking the equilibrium expectation value, the first-order derivative of the energy expectation value becomes So, the first-order derivatives can be computed using the derivatives of the Gaussian transformation operator, Eqs. (S2-S7).
Using the identities [Eqs. (S50-S59)] as well as the properties of the SCHA density matrix [Eqs. (6, 7)], the firstorder derivatives of energy at x = 0 can be computed as follows. We find that all first-order derivatives of the energy are zero. For the variational parameters included in the SCHA theory, the centroid position and the force constants, the stationarity ofρ 0 is expected sinceρ 0 is the variational solution that minimizes the SCHA free energy.ρ 0 is also stationary with respect to the variation of other parameters such as the atomic momentum parameter α i m because it is a thermal density matrix whose momentum expectation value is zero.

C. Second derivatives
Next, we calculate the second derivatives of energy. The result of this subsection can be summarized in a matrix form: The remaining part of this subsection is the derivation of Eq. (S68). By taking derivative of Eq. (17) withÔ =Ĥ, the second derivative of energy at x = 0 is given by When the two derivatives are for the same parameter type (displacement or squeezing), the second derivative of the transformation matrix becomes In this case, using ∂ µÛ † = −∂ µÛ , the second derivative of the energy can be written as For mixed second derivatives in which the derivatives are with respect to one displacement and one squeezing parameter, one finds and the same for γ instead of β. In this case, the second derivative of energy becomes and the same for γ instead of β.

S6. DERIVATION OF THE SCHA ANSATZ FROM THE TDVP SELF-ENERGY
In this section, we derive the SCHA ansatz Eq. (47) from the self-energy formula Eq. (44) which is derived from TDVP. In this section, the constraint m ≤ n in the summation over mode indices m and n is not implied. The constraint is made explicit whenever necessary by using smaller matrices which are defined only on the constrained indices: and W m n ,r s = W m n ,r s .
Here and in the remaining part of this section, we denote the constrained indices with primes: the index m n implies the constraint m ≤ n . Using these smaller matrices, Eq. (44) can be written as Next, we define a rectangular matrix R with matrix elements R m n ,rs = 1 if (r, s) = (m , n ) or (r, s) = (n , m ) 0 otherwise . (S132) By multiplying R to the smaller matrices, one can recover the full matrix: and R Φ (4) R = Φ (4) . (S134) These identities hold because Φ (S136) Equations (S135) and (S136) imply Using Eqs. (S133), (S134), and (S137), we can write Eq. (S131) as Equation (S138) is identical to Eq. (47) of the main text.

S7. DEGENERATE AND NEAR-DEGENERATE MODES
In this section, we detail the treatment of degenerate and near-degenerate modes. First, let us assume that states m and n are almost but not exactly degenerate. Although it is true that c mn [Eq. (13)] becomes large, what is important in the dynamics of the variational parameter γ mn in the linear response regime is the linearized time-evolution generator K. When computing K, the large c mn factors are counteracted by the small (n m − n n ) factors that appear when taking derivatives of E(x). Equations (S49, S90, S94) are examples of such counteraction. The fact that the equation of motion does not suffer any problems can be seen from Eqs. (26,27,31). One finds that c mn enters the time-evolution generator K only through C mn,pq [Eq. (31)]. There, the c mn factor is multiplied by (n m − n n ), so that C mn,pq is proportional to √ n m − n n . Hence, the matrix element of K involving near-degenerate states converges to zero in the limit of degeneracy. Next, let us consider exactly degenerate states. In the main text, we explained that if modes m and n degenerate, one needs to exclude γ mn from the set of variational parameters when studying the linear-response regime. The reason is that the infinitesimal transformation parametrized by γ mn does not change the variational equilibrium density matrixρ 0 . One should exclude parameters that do not change the variational density matrixρ(x) at that given x. The parameter γ mn for degenerate m and n is such a parameter for the equilibrium, x = 0. We note that such exclusion is necessary and justified only in the linear response regime. When considering large deviations from x = 0, one must include γ mn in the set of variational parameters.
Even if the energy of degenerate states m and n are slightly perturbed due to numerical inaccuracies, no problem will occur. As explained in the near-degenerate case, the corresponding equation of motion will be suppressed by a factor of √ n m − n n , so that γ mn stays at is equilibrium value, γ mn = 0.

S8. ZERO TEMPERATURE CASE
In the main text, we have focused only on the finite temperature case. At zero temperature, one should apply TDVP directly to the Gaussian wavefunctions without purification. The main difference with the finite temperature case is that the squeezing transformation parametrized by γ becomes a do-nothing operation at T = 0. This difference can be noticed by calculating the tangent vector by applying ∂Û /∂γ [Eqs. (S6, S7)] to the stationary state wavefunction. At T > 0, the purified stationary state wavefunction in the number basis has nonzero coefficients for states with nonzero phonon populations; hence, the tangent vectors do not vanish. On the contrary, at T = 0, the stationary state wavefunction is a vacuum state of the SCHA harmonic Hamiltonian. Hence, the rightmost annihilation operators in Eqs. (S6, S7) nullify the wavefunction and the corresponding tangent vectors become null vectors. So, at zero temperature, only α and β should be used as the variational parameters.
One can follow the same steps as in the finite temperature case to calculate the linearized time evolution generator and the position-position correlation function at zero temperature. The final form of the phonon self-energy is identical to the finite-temperature result, Eq. (47). The only difference is that the second term in the definition of χ [Eq. (46)] that originates from the variation of the γ parameter vanishes. Still, the equations need not be modified because the second term of Eq. (46) is already zero at T = 0 since n m = n n = 0.

S9. SINGLE-MODE ANHARMONIC HAMILTONIAN
In this section, we compute the excitation energy of the single-mode anharmonic Hamiltonian [Eq. (48)] using three different methods: perturbation theory, linearized time evolution, and projected Hamiltonian.