Open quantum dynamics with variational non-Gaussian states and the truncated Wigner approximation

Bond, Liam J.; Gerritsen, Bas; Minář, Jiří; Young, Jeremy T.; Schachenmayer, Johannes; Safavi-Naini, Arghavan

doi:10.1063/5.0226268

We present a framework for simulating the open dynamics of spin–boson systems by combining variational non-Gaussian states with a quantum trajectories approach. We apply this method to a generic spin–boson Hamiltonian that has both Tavis–Cummings and Holstein type couplings and which has broad applications to a variety of quantum simulation platforms, polaritonic physics, and quantum chemistry. Additionally, we discuss how the recently developed truncated Wigner approximation for open quantum systems can be applied to the same Hamiltonian. We benchmark the performance of both methods and identify the regimes where each method is best suited. Finally, we discuss strategies to improve each technique.

I. INTRODUCTION

Advances in the control of quantum systems over the past decade have led to the development of a wide variety of different platforms for investigating almost coherent quantum dynamics. However, in the absence of robust fault-tolerant operations, studies of these platforms must contend with various sources of noise induced by couplings to an environment. Moreover, these noisy intermediate-scale quantum (NISQ¹) devices can often be arbitrarily controlled, opening possibilities of creating states far out of equilibrium. In light of this, understanding the out-of-equilibrium dynamics of open quantum many-body systems has become of great interest, e.g., for understanding and realizing a quantum advantage in the NISQ era.

A particularly challenging class of problems arises for systems containing both spin and bosonic degrees of freedom with an (even locally) unbounded Hilbert space, applicable to quantum simulation and computation platforms ranging from superconducting circuits to trapped ions,^2–16 to paradigmatic problems in impurity physics,^17,18 quantum chemistry,¹⁹ or polaritonic chemistry.²⁰ The ubiquity and complexity of spin–boson Hamiltonians have led to the development of various techniques for their study and characterization. These include methods for the bosonic space, namely path integral techniques,^21–24 effective Hamiltonian formulation,^25–28 or lightcone conformal truncation used predominantly in high-energy physics,^29–31 as well as methods such as non-equilibrium Monte Carlo or tensor networks,^32–34 which have allowed for the simulation of (open) out-of-equilibrium dynamics of quantum many-body systems,^35–48 also in large-system scenarios. However, these methods are often constrained to one-dimensional setups, closed systems, a mesoscopic number of particles, or a combination thereof.

The number and breadth of the aforementioned approaches illustrate that no single method can tackle the range of systems described by spin–boson type Hamiltonians or even the full parameter space of a specific model. As such, it is important to pinpoint the strengths and shortcomings of each method. In this work, we undertake a comparative study between two methods: (i) the time-dependent variational ansatz using non-Gaussian states (NGS)^49,50 and (ii) the truncated Wigner approximation (TWA) combined with its generalization to discrete spaces, discrete truncated Wigner approximation (DTWA). Figures 1(a) and 1(b) show a schematic of each method. These methods may allow us to circumvent some of the aforementioned limitations^49–53 and have recently been generalized to open quantum systems.^54–61 Analyzing fundamental quantum effects in macroscopic limits can thus be enabled with NGS and TWA approaches.

FIG. 1.

(a) Schematic overview of NGS. The ansatz |ψz⃗(t)〉 is a point in a manifold embedded in Hilbert space. For different variational parameters, we can describe a variety of quantum states (some examples of accessible Wigner functions are shown). The action of an operator O, which can be a Hamiltonian H, non-Hermitian Hamiltonian Heff, or quantum jump c, can cause the state to leave the variational manifold, so the state is projected back to the manifold P|ψ〉. (b) Schematic depiction of TWA. (i) The initial quantum state of the spin and bosonic degrees of freedom can be represented by the spin Wigner functions Wi(θ0, ϕ0) and bosonic Wigner functions Wj(A0,A0*), respectively. (ii) Classical phase space points are individually sampled from the initial Wigner functions for each degree of freedom and (iii) evolve according to the classical (stochastic) equations of motion. (iv) Expectation values of observables are evaluated in phase space by computing the average of the associated Weyl symbols over ntraj phase space trajectories. (c) Both NGS and TWA apply to generic spin–boson systems, but we depict here the system studied in this work: Ns spins-1/2 interact with a common mode a with strength g and particle loss at rate κ. Each spin also interacts with a mode b at strength λ. The spins also undergo collective loss at rate Γ. (d) Schematic summary of our results illustrating the regions where each method tends to perform well together with reference to figures studying the dynamics in the respective parameter regimes (see Fig. 5 and Sec. V for more details).

View large Download slide

(a) Schematic overview of NGS. The ansatz $| ψ (\vec{z} (t)) 〉$ is a point in a manifold embedded in Hilbert space. For different variational parameters, we can describe a variety of quantum states (some examples of accessible Wigner functions are shown). The action of an operator $O$ ⁠, which can be a Hamiltonian H, non-Hermitian Hamiltonian H_eff, or quantum jump c, can cause the state to leave the variational manifold, so the state is projected back to the manifold $P | ψ 〉$ ⁠. (b) Schematic depiction of TWA. (i) The initial quantum state of the spin and bosonic degrees of freedom can be represented by the spin Wigner functions W_i(θ₀, ϕ₀) and bosonic Wigner functions $W_{j} (A_{0}, A_{0}^{*})$ ⁠, respectively. (ii) Classical phase space points are individually sampled from the initial Wigner functions for each degree of freedom and (iii) evolve according to the classical (stochastic) equations of motion. (iv) Expectation values of observables are evaluated in phase space by computing the average of the associated Weyl symbols over n_traj phase space trajectories. (c) Both NGS and TWA apply to generic spin–boson systems, but we depict here the system studied in this work: N_s spins-1/2 interact with a common mode a with strength g and particle loss at rate κ. Each spin also interacts with a mode b at strength λ. The spins also undergo collective loss at rate Γ. (d) Schematic summary of our results illustrating the regions where each method tends to perform well together with reference to figures studying the dynamics in the respective parameter regimes (see Fig. 5 and Sec. V for more details).

In NGS, one exploits the continuous variable structure of bosonic states to build a time-dependent variational wave function ansatz of non-Gaussian states. Here, we specifically use a superposition of squeezed displaced bosonic states, which converges to the true wave function due to the over-completeness of the set of coherent states. Since each state in the superposition is Gaussian, much of the previously developed machinery for Gaussian states can be reutilized. This method has been successfully applied to studies of systems ranging from the Kondo impurity problem,⁶² central spin,⁶³ spin-Holstein models,⁶⁴ Bose and Fermi polarons,^65–67 and (sub/super) Ohmic spin–boson model.⁶⁸ We also note that the closely related Davydov state ansatz has been applied in the studies of molecular crystals and polaritonic physics.^69–75

TWA is a semi-classical approach that factorizes the phase space functions (Weyl symbols) that describe a quantum observable in the phase space representation. As such, TWA is reminiscent of a product-state mean-field ansatz on Hilbert space. Like the latter, TWA allows one to treat systems with very large sizes [ $O (1 0^{4})$ particles] while still capturing some essential quantum features such as spin-squeezing or entanglement.^51,76–80 TWA can be easily adapted to systems with both bosonic and discrete degrees of freedom, combining sampling strategies from continuous⁸¹ and discrete Wigner functions.^51–53

Here, we consider the open and closed dynamics of a spin–boson Hamiltonian featuring multiple spins coupled to a discrete set of bosonic modes via Holstein and Tavis–Cummings couplings, thus ensuring the broad applicability of our results. We begin by introducing the two methods: the variational NGS method using the formulation introduced in Ref. 68 is discussed in Sec. II. We extend the method to open quantum systems using the quantum trajectories method in Sec. III. We discuss TWA with its discrete variant DTWA for closed and open systems in Sec. IV. In Sec. V, we introduce the Holstein–Tavis–Cummings spin–boson Hamiltonian, which we use to compare both methods over a range of parameters. Finally, in Sec. VI, we summarize our findings and discuss how each method can be improved to increase its accuracy and/or applicability both in terms of systems to which they can be applied and also in terms of the observables that can be accessed.

II. NON-GAUSSIAN ANSATZ FOR A CLOSED SYSTEM

We begin by introducing the non-Gaussian ansatz before discussing how to compute the equations of motion for the variational parameters. We consider a system of N_s spins-1/2 and N_b bosonic modes governed by some Hamiltonian H(t). Our wave function,

| ψ (\vec{z}) 〉

⁠, is a variational ansatz in the form of a non-Gaussian state parameterized by a set of real numbers

\vec{z}

⁠,

| ψ (\vec{z}) 〉 = \sum_{σ = 1}^{2^{N_{s}}} \sum_{p = 1}^{N_{p}} U_{p}^{(σ)} | σ, 0 〉,

(1)

where the summation over σ is over all

2^{N_{s}}

spin basis states.⁸² The summation over p produces a superposition of N_p bosonic states

U_{p}^{(σ)} | 0 〉

for each spin degree of freedom. We choose the operator

U_{p}^{(σ)}

to be of the form of a Gaussian unitary,

U_{p}^{(σ)} = e^{κ_{p}^{(σ)} + i θ_{p}^{(σ)}} D ({\vec{α}}_{p}^{(σ)}) S ({\vec{ζ}}_{p}^{(σ)}) .

(2)

The parameters κ_p and θ_p determine the weight and phase factors, while the many-mode displacement

D (\vec{α})

and squeezing

S (\vec{ζ})

operators are defined as

D (\vec{α}) = \prod_{k = 1}^{N_{b}} D (α_{k}) = \prod_{k = 1}^{N_{b}} \exp [α_{k} a_{k}^{†} - α_{k}^{*} a_{k}],

(3a)

S (\vec{ζ}) = \prod_{k = 1}^{N_{b}} S (ζ_{k}) = \prod_{k = 1}^{N_{b}} \exp [\frac{1}{2} (ζ_{k}^{*} a_{k}^{2} - ζ_{k} {a_{k}^{†}}^{2})],

(3b)

where we dropped the σ and p indices of

\vec{α}, \vec{ζ}

for simplicity. In this article, we always follow the convention that the many-mode operators appear with vector parameters

\vec{α}

and

\vec{ζ}

⁠, i.e.,

D (\vec{α})

⁠,

S (\vec{ζ})

⁠. Here, a_k

(a_{k}^{†})

is the annihilation (creation) operator of the kth mode satisfying the commutation relations

[a_{j}, a_{k}^{†}] = δ_{j k}

⁠. The complex numbers α_k = x_k + iy_k and

ζ_{k} = r_{k} e^{i ϕ_{k}}

describe the displacement and squeezing amplitudes of the kth mode, respectively. The parameters x_k, y_k, r_k, and ϕ_k are real, and we collect them, along with κ_p and θ_p for each Gaussian state, into a set {κ, θ, x, y, r, ϕ} indexed by ℓ. The total set of variational parameters

\vec{z}

is indexed by four indices:

σ \in {1, \dots, 2^{N_{s}}}

⁠, p ∈ {1, …, N_p}, k ∈ {1, …, N_b}, and ℓ ∈ {1, …, 6}. The total number of variational parameters is then

M = 2^{N_{s}} N_{p} (2 + 4 N_{b})

⁠.

The set of all coherent states forms an over-complete basis for the bosonic Hilbert space. Therefore, in the limit N_p → ∞, the NGS ansatz |ψ⟩, even without squeezing $({\vec{ζ}}_{p}^{(σ)} = 0 \forall p, σ)$ ⁠, approaches the true wave function |Ψ⟩. The inclusion of diagonal squeezing in our formalism, Eq. (2), is intended to enhance the descriptive power of the ansatz at finite N_p. However, in Fig. 2 and Sec. V, we show that for the systems studied in this work, a superposition of coherent states (i.e., $\vec{ζ}$ = 0) for reasonable N_p is typically sufficient to describe the relevant physics. One could further generalize the Gaussian unitary to $U \propto D (\vec{α}) \exp (- i A^{T} M A)$ instead of Eq. (2), where $A = {(a_{1}, \dots, a_{N_{b}}, a_{1}^{†}, \dots, a_{N_{b}}^{†})}^{T}$ and $M$ is a symmetric matrix.^49,50 However, we find that this significantly complicates all subsequent manipulations while not substantially reducing the N_p required to capture the relevant features of the systems studied in this work.

FIG. 2.

View large Download slide

Comparison between the NGS ansatz Eq. (1) with a superposition of displaced squeezed states vs a superposition of coherent states for evolution under the anharmonic oscillator, both with N_p = 4. The initial state is |ψ(0)⟩ = |α = 1⟩. The inclusion of squeezing (blue dashed) results in closer agreement with exact numerics (solid black) compared to the absence of squeezing in the ansatz (orange dashed). However, both capture the dynamics with small infidelities $1 - F < O (1 0^{- 2})$ ⁠. Here, $F (t) = | ⟨ Ψ (t) | ψ (t) ⟩ |$ ⁠, where |Ψ(t)⟩ is the quantum state obtained with exact numerics.

We start our analysis by adopting the Dirac–Frenkel variational principle.⁸³ In this framework, for a given H(t), one can derive equations of motion (EOMs) for the variational parameters

\vec{z} (t)

describing either real- or imaginary-time evolution of the wave function,⁵⁰

Real - time  ev.: {\dot{z}}^{ν} = - {(ω_{μ ν})}^{- 1} \partial_{μ} E (\vec{z}, t),

(4a)

Imag - time  ev.: {\dot{z}}^{ν} = - {(g_{μ ν})}^{- 1} \partial_{μ} ϵ (\vec{z}, t),

(4b)

where μ, ν = (σ, p, k, ℓ) index the variational parameters

\vec{z}

and ∂_μ = ∂/∂z^μ. Here,

E (\vec{z}, t) = 〈 ψ (\vec{z}) | H (t) | ψ (\vec{z}) 〉

is the energy and

ϵ (\vec{z}, t) = E (\vec{z}, t) / 〈 ψ (\vec{z}) | ψ (\vec{z}) 〉

is the normalized energy. We introduce the tangent vectors of the variational manifold at point

\vec{z}

⁠,

| v_{μ} 〉 = \partial_{μ} | ψ (\vec{z}) 〉 .

(5)

In terms of the tangent vectors, the symplectic form ω and the metric of the tangent space g are defined as

ω_{μ ν} = 2 Im 〈 v_{μ} | v_{ν} 〉,

(6a)

g_{μ ν} = 2 Re 〈 v_{μ} | v_{ν} 〉,

(6b)

with their inverses denoted

Ω^{μ ν} \equiv {(ω_{μ ν})}^{- 1}

and

G^{μ ν} \equiv {(g_{μ ν})}^{- 1}

⁠.

Next, we discuss a subtle property of the employed variational principle following the discussion in Ref. 50. We assume that the variational parameters are real,

\vec{z} \in R

⁠. Then, the tangent space

T_{ψ}

of the variational manifold at each point

| ψ (\vec{z}) 〉

is a real vector space spanned by the tangent vectors |v_μ⟩ embedded in complex Hilbert space. Therefore, for each basis vector |v_μ⟩, i|v_μ⟩ is not guaranteed to lie in the tangent space and has to be projected onto

T_{ψ}

⁠. As the tangent space is not a complex linear space, this projection takes the form

P_{ψ} i | v_{ν} 〉 = 2 | v_{μ} 〉 G^{μ σ} R e 〈 v_{σ} | i | v_{ν} 〉 = {J^{μ}}_{ν} | v_{μ} 〉,

(7)

where the complex structure, i.e., the representation of the projection of the imaginary unit, is introduced as

{J^{μ}}_{ν} = - G^{μ σ} ω_{σ ν}

⁠. If J² ≠ − 1, the projection is non-trivial, that is, i|v_μ⟩ does not lie in the tangent space. On the other hand, when J² = −1 on every tangent space, the tangent space of the variational manifold is called a Kähler space. Then,

i | v_{ν} 〉 = {J^{μ}}_{ν} | v_{μ} 〉

⁠. In this case, J specifies the decomposition of i|v_ν⟩ on the tangent space vectors |v_μ⟩. Having a tangent space

T_{ψ}

to be a Kähler space ensures that the three variational principles (Lagrangian, McLachlan, and Dirac–Frenkel) coincide and prevents incorrect couplings in the equations of motion, as shown in Ref. 50. Finally, we note that satisfying J² = −1 in the case of real variational parameters is equivalent to requiring a holomorphic parametrization in the case of complex parameters.

In the supplementary material, we prove that J² = −1 for the NGS ansatz, both in the case of a superposition of displaced squeezed states and a superposition of coherent states (supplementary material).

A. Analytic expressions for energies and energy gradients

Having introduced the NGS ansatz in Eq. (1) and its equations of motion in Eq. (4), we now show how to obtain analytic expressions for the two crucial ingredients to the equations of motion: the energy gradients ∂_μE and the geometric structures g, ω.

Consider a generic spin–boson Hamiltonian. Any such Hamiltonian can be cast in the form H = ∑(H_s ⊗ H_b), where H_s describes the spin degrees of freedom and H_b describes the bosonic degrees of freedom and can be decomposed into monomials of bosonic operators as $H_{b} = \otimes_{k = 1}^{N_{b}} H_{b}^{(k)}$ ⁠.

For a single spin, N_s = 1, the Pauli matrices {σ₀, σ_x, σ_y, σ_z} form a complete basis for H_s. For H_s = σ_z and H_s = σ_x, we obtain, respectively,

〈 σ_{z} {H_{b} 〉}_{ψ} = \sum_{p, p^{'}}^{N_{p}} 〈 {U_{p^{'}}^{(↑)}}^{†} H_{b} U_{p}^{(↑)} 〉 - 〈 {U_{p^{'}}^{(↓)}}^{†} H_{b} U_{p}^{(↓)} 〉,

(8a)

〈 σ_{x} {H_{b} 〉}_{ψ} = \sum_{p, p^{'}}^{N_{p}} 〈 {U_{p^{'}}^{(↑)}}^{†} H_{b} U_{p}^{(↓)} 〉 + 〈 {U_{p^{'}}^{(↓)}}^{†} H_{b} U_{p}^{(↑)} 〉,

(8b)

where the expectation values on the right-hand side are evaluated with respect to the bosonic vacuum |0⟩. Equivalent expressions for H_s = σ₀ and H_s = σ_y can be obtained using the same procedure. Therefore, the computation reduces to evaluating the many-mode overlap

⟨ 0 | {U_{p^{'}}^{(σ^{'})}}^{†} H_{b} U_{p}^{(σ)} | 0 ⟩

for each p, p′, σ, σ′. Because both the multi-mode displacement

D (\vec{α})

and squeezing

S (\vec{ζ})

operators in our ansatz are diagonal in mode operators, we can write the many-mode overlap as a product of single-mode overlaps,

\begin{align} ⟨ {U_{p^{'}}^{(σ^{'})}}^{†} H_{b} U_{p}^{(σ)} ⟩ \\ = \prod_{k = 1}^{N} 〈 0 | S^{†} (ζ_{k}^{(p^{'})}) D^{†} (α_{k}^{(p^{'})}) H_{b}^{(k)} D (α_{k}^{(p)}) S (ζ_{k}^{(p)}) | 0 〉, \end{align}

(9)

where

H_{b}^{(k)}

is the bosonic operators of H_b that act on mode k, as defined earlier. To evaluate the single mode overlaps in Eq. without specifying

H_{b}^{(k)}

⁠, we use the following identity:

〈 Θ | a^{†^{m}} a^{n} | {Φ 〉}_{s} = {{(\frac{\partial}{\partial β})}^{m} {(- \frac{\partial}{\partial β^{*}})}^{n} χ_{g}^{s} (β)|}_{β = 0},

(10)

where |Θ⟩, |Φ⟩ are any two quantum states⁸⁴ and

χ_{g}^{s} (β) = 〈 Θ | D (β) | Φ 〉 e^{(s / 2) β β^{*}},

(11)

is a generalized s-ordered characteristic function with s = 1 (s = −1) denoting (anti-)normal ordering of the bosonic operators a, a^† that appear in the left-hand side of Eq. . We note that in the calculation of the single mode overlaps in Eq. using Eq. , |Θ⟩ and |Φ⟩ are single mode Gaussian states. The overlap between any two single mode Gaussian states can be evaluated analytically,⁸⁵ which gives us an analytic expression for

χ_{g}^{s} (β)

⁠. Therefore, we can evaluate the partial derivatives in Eq. with respect to β, β* for any m, n and find analytic expressions for the energy of any Hamiltonian that is polynomial in a, a^† operators. From this, we can also obtain analytic expressions for the energy gradients, ∂_μE, as needed in the EOMs.

Example: Energy gradient of harmonic oscillator To demonstrate the above-mentioned machinery, we compute the energy gradient of the harmonic oscillator H = a^†a with respect to the NGS ansatz with coherent states only, i.e.,

\vec{ζ} = 0

⁠. We set N_p = 2, so

| ψ 〉 = e^{κ_{1} + i θ_{1}} | α_{1} 〉 + e^{κ_{2} + i θ_{2}} | α_{2} 〉

⁠, with α_i = x_i + iy_i. The energy is

\begin{array}{l} E = ⟨ H ⟩ & = e^{2 κ_{1}} (x_{1}^{2} + y_{1}^{2}) + e^{2 κ_{2}} (x_{2}^{2} + y_{2}^{2}) \\ + (e^{κ_{1} + κ_{2} + i (θ_{1} - θ_{2})} (x_{1} + i y_{1}) (x_{2} - i y_{2}) + h.c.), \end{array}

and its partial derivatives with respect to x₁ and κ₁ are given by

\partial_{x_{1}} E = 2 e^{2 κ_{1}} x_{1} + [e^{κ_{1} + κ_{2} + i (θ_{1} - θ_{2})} (x_{2} - i y_{2}) + h.c.]

and

\partial_{κ_{1}} E = 2 e^{2 κ_{1}} (x_{1}^{2} + y_{1}^{2}) + [e^{κ_{1} + κ_{2} + i (θ_{1} - θ_{2})} (x_{1} + i y_{1}) (x_{2} - i y_{2})]

⁠, respectively. The partial derivatives

\partial_{y_{1}} E

⁠,

\partial_{x_{2}} E

⁠,

\partial_{y_{2}} E

⁠,

\partial_{κ_{2}} E

⁠,

\partial_{θ_{1}} E

⁠, and

\partial_{θ_{2}} E

can be evaluated using the same procedure. We note that these expressions can be extended to any N_p and to any bosonic Hamiltonian that is polynomial in a, a^† using Eq. .

B. Tangent vectors and the overlap matrix

Next, we explain how to compute the tangent vectors of the ansatz,

| v_{μ} 〉 = \partial_{μ} | ψ (\vec{z}) 〉 = \partial_{σ, p, k, ℓ} | ψ (\vec{z}) 〉

⁠. Plugging Eq. into Eq. , we find

\begin{align} | v_{μ} 〉 & = \frac{\partial}{\partial z_{σ, p, k, ℓ}} U_{p}^{(σ)} | σ, 0 〉 = \prod_{k^{'} \neq k}^{N_{b}} [D (α_{k^{'}}) S (ζ_{k^{'}})] \\ \times \frac{\partial}{\partial z_{σ, p, k, ℓ}} [e^{κ + i θ} D (α_{k}) S (ζ_{k})] | σ, 0 〉 . \end{align}

(12)

The calculation involving many modes and a superposition of squeezed-displaced states now reduces to calculating the tangent vector of a single-mode Gaussian state for a single Gaussian, i.e., $\frac{\partial}{\partial z_{σ, p, k, ℓ}} e^{κ + i θ} D (α_{k}) S (ζ_{k}) | 0 〉$ ⁠. We proceed as follows: (i) using the disentangled forms of the displacement and squeezing operators, we normal order $D (α_{k}) S (ζ_{k}) | 0 〉$ such that only a^† operators remain, which (ii) enables us to obtain concise analytic expressions for the tangent vectors of $D (α_{k}) S (ζ_{k}) | 0 〉$ ⁠, and (iii) obtain expressions for the general NGS ansatz and outline the computation of the overlap matrix ω_μν.

The disentangled forms of the single-mode displacement and squeezing operators are given by

D (α) = e^{α a^{†} - α^{*} a} = e^{- \frac{| α |^{2}}{2}} e^{α a^{†}} e^{- α^{*} a},

(13)

S (ζ) = e^{\frac{1}{2} (ζ {a^{†}}^{2} - ζ^{*} a^{2})} = e^{\bar{r} {a^{†}}^{2}} e^{- \overset{̌}{r} (a^{†} a + \frac{1}{2})} e^{- \bar{r} a^{2}},

(14)

where ζ = re^iϕ,

\bar{r} = e^{i ϕ} \tanh (r) / 2

⁠, and

\overset{̌}{r} = \ln (\cosh (r))

⁠. When applied to

| ψ (\vec{z}) 〉

⁠,

S (ζ)

always acts directly on the bosonic vacuum, which simplifies its action to

S (ζ) | 0 〉 = \frac{1}{\sqrt{\cosh r}} e^{\bar{r} {a^{†}}^{2}} | 0 〉 .

(15)

After normal-ordering

D (α) S (ζ) | 0 〉

⁠, we arrive at

D (α) S (ζ) | 0 〉 = \frac{e^{- \frac{| α |^{2}}{2}} e^{\bar{r} {α^{*}}^{2}}}{\sqrt{\cosh (r)}} e^{α a^{†}} e^{\bar{r} ({a^{†}}^{2} - 2 α^{*} a^{†})} | 0 〉,

(16)

where we used the relation e^xaf(a,a^†)e^−xa = f(a,a^† + x) [Ref. 86, §3.3 Theorem 2] and the fact that

e^{- α^{*} a} | 0 〉 = | 0 〉

⁠. Note that this expression contains only a^† and, therefore, all terms commute, enabling us to take the derivative with respect to the variational parameters to obtain the tangent vectors.

Computing tangent vectors. For each p, σ, and k, there are six variational parameters: {κ, θ, x, y, r, ϕ}. The tangent vectors for the norm κ and phase θ are simply

\partial_{κ} U_{p}^{(σ)} | σ, 0 〉 = U_{p}^{(σ)} | σ, 0 〉,

(17a)

\partial_{θ} U_{p}^{(σ)} | σ, 0 〉 = i U_{p}^{(σ)} | σ, 0 〉,

(17b)

and are independent of the mode number k. We use the normal-ordered form of

D (α) S (ζ) | 0 〉

⁠, Eq. , and define

| α, ζ 〉 \equiv D (α) S (ζ) | 0 〉

to find the tangent vectors for x, y, r, and ϕ,

\partial_{x} | α, ζ 〉 = (f_{x} + g_{x} a^{†}) | α, ζ 〉,

(18a)

\partial_{y} | α, ζ 〉 = (f_{y} + g_{y} a^{†}) | α, ζ 〉,

(18b)

\partial_{r} | α, ζ 〉 = (f_{r} + g_{r} a^{†} + h_{r} {a^{†}}^{2}) | α, ζ 〉,

(18c)

\partial_{ϕ} | α, ζ 〉 = (f_{ϕ} + g_{ϕ} a^{†} + h_{ϕ} {a^{†}}^{2}) | α, ζ 〉,

(18d)

where we have introduced the following c-number functions f, g:

f_{x} = e^{i ϕ} \tanh (r) (x - i y) - x,

(19a)

g_{x} = 1 - e^{i ϕ} \tanh (r),

(19b)

f_{y} = - e^{i ϕ} \tanh (r) (y + i x) - y,

(19c)

g_{y} = i (1 + e^{i ϕ} \tanh (r)),

(19d)

f_{r} = (1 / 2) [- \tanh (r) + e^{i ϕ} s e c h {(r)}^{2} {(x - i y)}^{2}],

(19e)

g_{r} = (1 / 2) e^{i ϕ} s e c h {(r)}^{2} (2 i y - 2 x),

(19f)

h_{r} = (1 / 2) e^{i ϕ} s e c h {(r)}^{2},

(19g)

f_{ϕ} = (i / 2) e^{i ϕ} \tanh (r) {(x - i y)}^{2},

(19h)

g_{ϕ} = (i / 2) e^{i ϕ} \tanh (r) (2 i y - 2 x),

(19i)

h_{ϕ} = (i / 2) e^{i ϕ} \tanh (r) .

(19j)

Note that further simplification of the terms with creation operators in Eq. (18) in the form of $D (α) S (ζ) | 1 〉$ and $D (α) S (ζ) | 2 〉$ is not possible as we imposed that $S (ζ)$ must act directly on the vacuum to obtain the simplified version of $S (ζ) | 0 〉$ in Eq. (15).

Equations (18a)–(18d) are simple analytic expressions for the single mode tangent vectors of $D (α) S (ζ) | 0 〉$ ⁠. When combined with Eq. (12), we can, therefore, construct all the tangent vectors for any σ, p, and k.

Computing tangent vector overlaps to construct geometric structures. Having obtained the tangent vectors, we finally briefly comment on the calculation of the overlap matrix ⟨v_μ|v_ν⟩, which is used to construct the metric g_μν and symplectic form ω_μν. We note that due to Eq. (12), the overlap between two tangent vectors ⟨v_μ|v_ν⟩ can be written as a product of single-mode overlaps. Each single-mode overlap can be evaluated using the expressions in Eq. (18) and applying Eq. (10) to evaluate expectation values of the type $〈 α^{'}, ζ^{'} | {(a^{†})}^{m} a^{n} | α, ζ 〉$ ⁠.

Example: Quantum anharmonic oscillator As an example of the results of the above-mentioned machinery and to illustrate the role of squeezing in the ansatz, we use Eq. (4) and solve for the dynamics of a simple anharmonic oscillator, $H = ω a^{†} a + μ {(a^{†} a)}^{2}$ ⁠.^87–90 Results for the strong coupling regime setting ω = μ = 1 are shown in Fig. 2 with N_p = 4. We observe that the ansatz containing squeezing (dashed blue line) has a better agreement with the actual state |Ψ⟩ (solid black line, obtained from the exact numerical solution of the Schrödinger equation) than the ansatz without squeezing (dashed orange line). While this observation always depends on the specific system being studied, we observe that for many applications in this work, the improvement in accuracy from including squeezing comes at the expense of additional computational resources. Motivated by this, for the remainder of this article our results utilize the NGS ansatz of a superposition of (many-mode) coherent states.

III. OPEN DYNAMICS: COMBINING NGS WITH QUANTUM TRAJECTORIES

To extend the NGS machinery to open dynamics, there are several options available. For instance, in Ref. 55, which in turn builds on the developments in Ref. 54, an ansatz for the density matrix was developed. This was then used to formulate the master equation and find the corresponding equations of motion governed by the Lindbladian.

Here, we discuss how the ansatz introduced previously, $| ψ (\vec{z}) 〉$ in Eq. (1), can be used within the quantum trajectories framework. In Sec. III A, we recall the basic formulation of the quantum trajectories approach before formulating equations of motion for the effective non-Hermitian Hamiltonian in Sec. III B. We give a recipe for how to formulate the action of quantum jumps within the NGS ansatz in Sec. III C, and we finally discuss the relevant issues related to the implementation of the equations of motion in Sec. III D. To elucidate the formalism, we provide specific examples throughout this section. Our examples are motivated by typical sources of decoherence in spin–boson quantum simulation platforms governed by Hamiltonians of the form that will be introduced in Sec. V.

A. Quantum trajectories

Many problems of the dynamics of open quantum systems are amenable to the standard time-local master equation in the Lindblad form^91–93

\dot{ρ} = - i [H, ρ] + \sum_{m} c_{m} ρ c_{m}^{†} - \frac{1}{2} \{c_{m}^{†} c_{m}, ρ\},

(20)

where c_m are the jump operators and {⋅, ⋅} the anticommutator.

An alternative approach is to use the quantum trajectories method.^47,94–96 Here, rather than evolving the full density matrix described by

2^{2 N_{s}}

elements (for N_s spin-1/2 particles), one stochastically evolves a pure quantum state |ψ⟩ described by

2^{N_{s}}

elements. This approach consists of two steps: (i) continuous evolution under the Schrödinger equation with an effective non-Hermitian Hamiltonian,

H_{eff} = H - \frac{i}{2} \sum_{m} c_{m}^{†} c_{m} = H - i K,

(21)

and (ii) discrete evolution under the action of a quantum jump operator c_m. Such stochastic evolution of the wave function |ψ⟩ constitutes a quantum trajectory. Observables of interest are then evaluated by averaging over n_traj such quantum trajectories. In scenarios where it is possible to obtain the observables of interest with sufficient accuracy using

n_{traj} < 2^{N_{s}}

⁠, quantum trajectories can be more efficient as compared to solving the master equation.^47,48,95

Our key contribution is to perform both (i) and (ii) using NGS. The details of the quantum trajectories method and a step-by-step description of its implementation can be found in Ref. 47.

B. Equations of motion for non-Hermitian Hamiltonians

In between the quantum jumps, the wave function evolves continuously under the effective non-Hermitian Hamiltonian Eq. . Here, we follow the procedure in Ref. 97 to derive equations of motion for evolution under H_eff. Note that the same procedure is used to derive the real-time and imaginary-time equations of motion given in Eq. describing evolution under a Hermitian Hamiltonian H. The NGS wave function evolves according to the Schrödinger equation

\partial_{t} | ψ 〉 = - i H_{eff} | ψ 〉 .

(22)

McLachlan’s variational principle requires the variation of the norm resulting from the Schrödinger equation to vanish,

δ ‖ (d / d t + i H_{eff}) | ψ 〉 ‖ = 0 .

(23)

Since it is more convenient to work with the square norm, we rewrite Eq. as δ‖(d/dt + iH_eff)|ψ⟩‖² = 0, obtaining

\begin{align} ‖ (d / d t + i H_{eff}) | ψ 〉 ‖^{2} & = \sum_{μ, ν} 〈 v_{μ} | v_{ν} 〉 {\dot{z}}_{μ} {\dot{z}}_{ν} + (i 〈 v_{μ} | H_{eff} | ψ 〉 {\dot{z}}_{μ} + h.c.) \\ + 〈 ψ | {H_{eff}}^{†} H_{eff} | ψ 〉 . \end{align}

(24)

After making the following substitutions:

A_{μ ν} = 〈 v_{μ} | v_{ν} 〉,

(25a)

C_{μ} = 〈 ψ | H | v_{μ} 〉,

(25b)

D_{μ} = 〈 ψ | K | v_{μ} 〉,

(25c)

the variation of the square of the norm is

\begin{align} δ ‖ (d / d t + i H_{eff}) | ψ 〉 ‖^{2} & = \sum_{μ} \sum_{ν} (A_{μ ν} + h.c.) {\dot{z}}_{ν} δ z_{μ} \\ + [i (C_{μ}^{†} - i D_{μ}^{†}) + h.c.] δ z_{μ}, \end{align}

(26)

which yields the equations of motion,

\sum_{ν} R e [A_{μ ν}] {\dot{z}}_{ν} = R e [D_{μ}] + I m [C_{μ}] .

(27)

The object 2Re[A_μν] is precisely the metric g of the tangent space introduced in Eq. , which we showed how to compute in Sec. II B. Re[D_μ] can be related to ∂_μ⟨K⟩ via

2 R e [〈 v_{μ} | K | ψ 〉] = 〈 v_{μ} | K | ψ 〉 + h.c. = \partial_{μ} ⟨ K ⟩,

(28)

where we substituted the definition of the tangent vector into the definition of D in Eq. and used the fact that K = K^†. We showed how to analytically compute the gradient of expectation value such as ∂_μ⟨K⟩ in Sec. II A.

Although Im[C_μ] can be computed by directly calculating the overlaps in the definition in Eq. , if the tangent space is a Kähler manifold (i.e., J² = −1, as discussed in Sec. II), then from the relation

{J^{μ}}_{ν} | v_{μ} 〉 = i | v_{ν} 〉

we can relate Im[C_μ] to the complex structure

{J^{μ}}_{ν}

and the gradient ∂_μ⟨H⟩ using

\begin{align} 2 I m [C_{ν}] & = i 〈 v_{ν} | H | ψ 〉 - i 〈 ψ | H | v_{ν} 〉 \\ = - 〈 v_{μ} | {({J^{μ}}_{ν})}^{T} H | ψ 〉 - 〈 ψ | H {J^{μ}}_{ν} | v_{μ} 〉 \\ = - \sum_{μ} {J^{μ}}_{ν} \partial_{μ} E . \end{align}

(29)

As such, the non-Hermitian equations of motion in Eq. can be computed from g, ∂_μ⟨H⟩, and ∂_μ⟨K⟩. The total number of elements to compute scales as M² + 2M, compared to M² + M for purely real- or imaginary-time evolution, where M is the number of variational parameters.

Example: Computing Im[C_μ] for the coherent state ansatz. Here, we relate Im[C_μ] to the energy gradients ∂_μ⟨H⟩ for the coherent state ansatz with explicit normalization and phase factors,

| ψ 〉 = e^{κ + i θ} e^{- (x^{2} + y^{2}) / 2} e^{(x + i y) a^{†}} | 0 〉

⁠, with

\vec{z} = (κ, θ, x, y) \in R

⁠. The tangent vectors corresponding to μ ∈ {1, 2, 3, 4} are

| v_{1} 〉 = | ψ 〉,

(30a)

| v_{2} 〉 = i | ψ 〉,

(30b)

| v_{3} 〉 = (a^{†} - x) | ψ 〉,

(30c)

| v_{4} 〉 = (i a^{†} - y) | ψ 〉 .

(30d)

We begin by relating Im[C₁] = Im[⟨ψ|H|v₁⟩] to

\partial_{z_{2}} E

⁠,

\begin{align} I m [C_{1}] & = - \frac{1}{2} [i 〈 ψ | H | v_{1} 〉 - i 〈 v_{1} | H | ψ 〉] \\ = - \frac{1}{2} (〈 ψ | H | v_{2} 〉 + 〈 v_{2} | H | ψ 〉) = - \frac{1}{2} \partial_{z_{2}} E . \end{align}

(31)

Similarly for Im[C₂] = Im[⟨ψ|H|v₂⟩], we have

I m [C_{2}] = \frac{1}{2} \partial_{z_{1}} E .

(32)

To relate Im[C₃] = Im[⟨ψ|H|v₃⟩] to the gradients of E, we notice that we can write i|v₃⟩ as a real span of the tangent vectors {v_μ}. That is, i|v₃⟩ = |v₄⟩ + y|v₁⟩ − x|v₂⟩. Using this, we obtain

I m [C_{3}] = - \frac{1}{2} (\partial_{z_{4}} E + y \partial_{z_{1}} E - x \partial_{z_{2}} E),

(33)

and finally for Im[⟨v₄|H|ψ⟩], using i|v₄⟩ = −|v₃⟩ − x|v₁⟩ − y|v₂⟩, we obtain

I m [〈 v_{4} | H | ψ 〉] = \frac{1}{2} (\partial_{z_{3}} E + x \partial_{z_{1}} E + y \partial_{z_{2}} E) .

(34)

Therefore, we have related the computation of Im[C_μ] to the energy gradients of E = ⟨H⟩. In this example, we explicitly noticed that i|v_μ⟩ can be written as a real span of tangent vectors {|v_μ⟩}. In the supplementary material, we provide the constructions of g_μν, ω_μν, and

{J^{μ}}_{ν}

for this single-mode coherent state ansatz, showing that

| v_{ν} 〉 = {J^{μ}}_{ν} | μ 〉

⁠, as derived in Eq. . We provide

{J^{μ}}_{ν}

for the NGS ansatz and prove that J² = −1.

Finally, we provide a simple example of non-Hermitian dynamics using NGS and the equations of motion in Eq. (27) by computing the evolution of the coherent state ansatz for N_b = 1 mode under H_eff = ξ(a + a^†) − (i/2)κa^†a. Since the Hamiltonian is a Gaussian operator, if the initial state of interest can be accurately described by the NGS ansatz with N_p coherent states, NGS is exact in that it captures precisely the dynamics at all times also with N_p coherent states. The results for an N_p = 2 initial state and a comparison against exact numerics are shown in Fig. 3, depicting perfect agreement as expected.

FIG. 3.

Dynamics of the effective Hamiltonian Heff = ξ(a + a†) − (i/2)κa†a with ξ = 0.5, κ = 1.0 using the NGS ansatz without squeezing (i.e., ζ⃗=0) with Np = 2. The non-Hermitian term corresponds to a cavity loss jump operator. The initial state is a random superposition of two coherent states. The NGS results (dashed colors) agree perfectly with exact numerics (solid black), including the decay of the norm ⟨ψ|ψ⟩ due to the non-hermitian term.

View large Download slide

Dynamics of the effective Hamiltonian H_eff = ξ(a + a^†) − (i/2)κa^†a with ξ = 0.5, κ = 1.0 using the NGS ansatz without squeezing (i.e., $\vec{ζ} = 0$ ⁠) with N_p = 2. The non-Hermitian term corresponds to a cavity loss jump operator. The initial state is a random superposition of two coherent states. The NGS results (dashed colors) agree perfectly with exact numerics (solid black), including the decay of the norm ⟨ψ|ψ⟩ due to the non-hermitian term.

C. Quantum jumps

We will now incorporate the action of a quantum jump c_m in the NGS formalism. After each quantum jump, the wave function may (i) remain in the variational manifold or (ii) leave it. We discuss these two possibilities using the examples of single particle loss and gain, respectively. We have chosen these two processes as they are often the dominant sources of single particle decoherence in a variety of systems.

1. Jumps inside the manifold

To demonstrate the effect of a jump operator that produces a state that remains within the variational manifold, we consider single particle loss at rate κ⁽¹⁾. The jump operator is

c = \sqrt{κ^{(1)}} a

⁠. The action of c on the single-mode NGS ansatz in Eq. without squeezing, i.e.,

| \vec{ζ} | = 0 \forall p

⁠, is given by

\begin{align} \sqrt{κ^{(1)}} a | ψ 〉 & = \sqrt{κ^{(1)}} a \sum_{p} e^{κ_{p} + i θ_{p}} D (α_{p}) | 0 〉 \\ = \sum_{p} \sqrt{κ^{(1)}} e^{κ_{p} + i θ_{p}} α_{p} D (α_{p}) | 0 〉 \\ = \sum_{p} e^{κ_{p}^{'} + i θ_{p}^{'}} D (α_{p}) | 0 〉, \end{align}

(35)

where we used a|α⟩ = α|α⟩ and defined the updated norm and phase factors as

κ_{p}^{'} = \log (\sqrt{κ^{(1)}} e^{κ_{p}} | α_{p} |),

(36a)

θ_{p}^{'} = \arg (| κ^{(1)} | e^{i θ_{p} + κ_{p}} α_{p}) .

(36b)

We can easily extend the above-mentioned analysis to the two-photon loss case with the jump operator

c = \sqrt{κ^{(2)}} a^{2}

⁠, where the updated norm and phase factors are now defined as

κ_{p}^{'} = \log (\sqrt{κ^{(2)}} e^{κ_{p}} | α_{p} |^{2}),

(37a)

θ_{p}^{'} = \arg (\sqrt{κ^{(2)}} e^{i θ_{p} + κ_{p}} α_{p}^{2}) .

(37b)

It is important to note that our results rely on the coherent state being the eigenstate of the jump operator considered earlier. Therefore, for any other state, e.g., the squeezed coherent state |α, ζ⟩, the single-particle loss jump operator will take the state out of the variational manifold. This will be a more generic scenario for most jump operators. In Sec. III C 2, we describe how to deal with such situations.

2. Jumps outside the manifold

If the state after the jump is not within the variational manifold, we project it back to the manifold by maximizing its fidelity with the variational ansatz. To do so, we use gradient descent (GD) as an efficient numerical procedure. We note that while simulated annealing (SA) finds a global extremum of a function in the asymptotic limit of infinitely slow cooling rate, we find that in all cases studied in this work, the performance of GD is comparable to that of SA (the infidelity difference between the post-jump variational state found by each of the two methods is ≲ 10⁻³), with the advantage that GD is typically significantly faster.

Again, we consider a single mode NGS ansatz without squeezing (i.e.,

\vec{ζ} = 0

⁠), which we write as

| ψ 〉 = \sum_{i = 1}^{N_{p}} c_{i} | α_{i} 〉,

(38)

with complex coefficients c_i and complex amplitudes α_i. For single photon gain with jump operator c = a^†, the state after applying the jump operator a^†|ψ⟩ is projected back onto a generic variational state

| \tilde{ψ} 〉

by optimizing its variational parameters

{\tilde{c}}_{i}, {\tilde{α}}_{i}

to maximize the normalized fidelity given by

F = \frac{| ⟨ \tilde{ψ} | a^{†} | ψ ⟩ |^{2}}{| ⟨ \tilde{ψ} | \tilde{ψ} ⟩ | | ⟨ ψ | a a^{†} | ψ ⟩ |},

(39)

where the unnormalized overlap is

⟨ \tilde{ψ} | a^{†} | ψ ⟩ = \sum_{i, j}^{N_{p}} {\tilde{c}}_{i}^{*} c_{j} {\tilde{α}}_{i}^{*} e^{- \frac{1}{2} (| α |^{2} + | \tilde{α} |^{2}) + {\tilde{α}}_{i}^{*} α_{j}},

(40)

and the normalization factors are

〈 ψ | ψ 〉 = \sum_{i, j}^{N_{p}} c_{i}^{*} c_{j} e^{- \frac{1}{2} (| α_{i} |^{2} + | α_{j} |^{2} + α_{i}^{*} α_{j})},

(41a)

〈 ψ | a a^{†} | ψ 〉 = \sum_{i, j}^{N_{p}} c_{i}^{*} c_{j} (1 + α_{i}^{*} α_{j}) e^{- \frac{1}{2} (| α_{i} |^{2} + | α_{j} |^{2}) + α_{i}^{*} α_{j}} .

(41b)

In the case of N_p = 1, with $α = | α | e^{i φ_{α}}$ ⁠, $\tilde{α} = | \tilde{α} | e^{i φ_{\tilde{α}}}$ ⁠, Eq. (40) is maximized for $φ_{\tilde{α}} = φ_{α}$ and for the amplitude of $\tilde{α}$ ⁠, which is the solution to $| \tilde{α} |^{2} - | α | | \tilde{α} | - 1 = 0$ ⁠. For instance, if α = 0, $| \tilde{α} | = 1, \forall φ_{\tilde{α}}$ ⁠, cf. Fig. 4(a), and similarly for α ≠ 0, cf. Fig. 4(b). The starting point of the GD search is denoted by the red cross and the maximum of the numerically found maximum of the overlap, Eq. (40), by the white cross. We remark that the achievable fidelity after projecting back to the manifold is strongly dependent on the jump operator and the number of coherent states. For instance, for the case shown in Fig. 4(a), the state after the jump is a Fock state |1⟩. As such, after projecting it back to a single coherent state, its fidelity is given by $e^{- | \tilde{α} |^{2} / 2} | α |^{n} / \sqrt{n!} \to e^{- 1 / 2} \approx 0.61$ with n = 1 and $| \tilde{α} | = 1$ ⁠.

FIG. 4.

(a) and (b) Fidelity of a†|α⟩ projected back onto the variational manifold here formed by the set of all coherent states |α̃=x̃+iỹ〉. The red cross denotes the initial state upon which acts the jump operator a†, and the white cross denotes the maximum of Eq. (40) found via GD. The used initial states are (a) |α = 0⟩ (with a ring of global maxima due to the symmetry of a†|0⟩) and (b) |α = 0.3 + i0.8⟩. (c) and (d) Difference in optimized fidelities, cf. Eq. (39), found by GD vs bounded SA. The jump operators are (c) a† and (d) x = a + a†.

View large Download slide

(a) and (b) Fidelity of a^†|α⟩ projected back onto the variational manifold here formed by the set of all coherent states $| \tilde{α} = \tilde{x} + i \tilde{y} 〉$ ⁠. The red cross denotes the initial state upon which acts the jump operator a^†, and the white cross denotes the maximum of Eq. (40) found via GD. The used initial states are (a) |α = 0⟩ (with a ring of global maxima due to the symmetry of a^†|0⟩) and (b) |α = 0.3 + i0.8⟩. (c) and (d) Difference in optimized fidelities, cf. Eq. (39), found by GD vs bounded SA. The jump operators are (c) a^† and (d) x = a + a^†.

In the case of the NGS ansatz with N_p > 1, the optimization landscape becomes more complex. In Figs. 4(c) and 4(d), for N_p = 1 (blue), N_p = 2 (orange), and N_p = 3 (green), we use Eq. (39) and plot the difference in optimized fidelities $| F_{GD} - F_{SA} |$ ⁠, with the optimization performed by GD with backtracking $(F_{GD})$ and bounded SA $(F_{SA})$ ⁠. We consider two quantum jumps, (c) single-particle gain c = a^† and (d) momentum kicks c = x = a + a^†. Our choice of single-particle gain is motivated by its relevance to many spin–boson systems (see Secs. V and VI), while the momentum kick jump plays a crucial role in laser cooling large ion crystals.⁹⁸ We generate 10³ initial (pre-jump) random states |ψ⟩. As the generic variational starting state $| \tilde{ψ} 〉$ to be optimized, for GD we use the pre-jump state |ψ⟩, while for the SA we seed a random starting state $| \tilde{ψ} 〉$ whose coefficients [see Eq. (38)] are drawn from a uniform unit distribution ${c_{i}, α_{i}, {\tilde{c}}_{i}, {\tilde{α}}_{i}} \in [0, 1]$ ⁠. We set a sufficiently slow SA cooling rate such that the algorithm converges to the same local maxima irrespective of the randomly chosen starting point. For all the studied cases, the fidelities of the states obtained by the two numerical optimizers agree within ≲10⁻³.

D. Time evolution

We are now equipped to implement the quantum trajectories program for the NGS ansatz outlined in Sec. III A. In principle, one could evolve the wave function |ψ⟩ between the jumps according to the equations of motion in Eq. while tracking the decay of the norm to identify the time t_j of a jump. In practice, we find it convenient to Trotterize the time evolution between jumps governed by H_eff, Eq. , in the usual way as

e^{- i (H - i K) δ t} \approx e^{- i H δ t} e^{- K δ t},

(42)

for sufficiently small δt. The norm during the unitary dynamics under e^−iHδt is preserved due to the use of the norm factors κ as variational parameters (we note that the opposite case, namely the absence of a global norm factor as a variational parameter, can result in unphysical couplings, see Ref. 50). During the imaginary-time evolution e^−Kδt, we track the decay of the norm to determine the time of the jump t_j, at which point we apply the corresponding jump operator c_m.

IV. TRUNCATED WIGNER APPROXIMATION FOR SPINS AND BOSONS

The phase space picture of quantum mechanics provides alternative means to simulate and analyze the quantum many-body dynamics in systems with mixed spin and bosonic degrees of freedom, in particular in a semi-classical framework known as the truncated Wigner approximation (TWA).⁸¹ The Wigner–Weyl transform maps Hilbert space operators O of a quantum system to functions of classical phase space variables, known as Weyl symbols O_W. The Weyl symbol corresponding to the density matrix is known as the Wigner function and provides a full ensemble description of arbitrary quantum states in terms of a (potentially negative) quasi-probability distribution.

A general Wigner–Weyl transform can be defined using the framework of phase point operators.⁹⁹ For example, considering particles in 1D with positions x and momenta p, operators $A (x, p)$ for each point in phase space define the Wigner–Weyl transformation via $O_{W} (x, p) = t r [A (x, p) O] / \sqrt{2 π}$ and vice versa $O = \int d x d p O_{W} (x, p) A (x, p)$ ⁠. Given a proper orthonormal definition of phase point operators,⁹⁹ for any state ρ and any observable Q, expectation values can be evaluated from the Wigner function $W (x, p) = t r [A (x, p) ρ] / \sqrt{2 π}$ via ⟨Q⟩ = ∫dxdp Q_W(x, p)W(x, p). Equivalent constructions can be made for spin phase spaces, either using spin–boson mappings (suitable for large spin scenarios),⁸¹ using spherical coordinate representations of spins A(θ, ϕ),⁵⁹ or for phase spaces using only a discrete set of points.⁹⁹

Closed-system time-evolution equations of motion can be obtained by Wigner–Weyl transforming the Heisenberg equations of motion, which leads to the exact quantum dynamics for Weyl symbols being governed by ${\dot{Q}}_{W} (x, p) = {Q_{W} (x, p), H_{W} (x, p)}_{MB}$ ⁠, using the Weyl symbol of the Hamiltonian H_W, and the Moyal bracket defined as {Q_W, H_W} = 2Q_W sin(Λ/2)H_W, with Λ the symplectic operator (with ℏ ≡ 1) $Λ = \sum_{i} \frac{\overset{⃖}{\partial}}{\partial x_{i}} \frac{\vec{\partial}}{\partial p_{i}} - \frac{\overset{⃖}{\partial}}{\partial p_{i}} \frac{\vec{\partial}}{\partial x_{i}}$ ⁠. Expanding the sine function in the Moyal bracket at the lowest order is known as TWA and leads to a classical evolution of Weyl symbols ${\dot{Q}}_{W} (x, p) \approx {Q_{W} (x, p), H_{W} (x, p)}_{P}$ ⁠, where {⋅,⋅}_P now denotes the classical Poisson bracket. The Poisson bracket ensures that the Weyl symbols for any complex observable will always factorize into phase space variables and, therefore, in TWA it suffices to only compute the classical evolution of the phase space variables.⁵³ This makes TWA a very practical and efficient numerical method for the case of a positive initial Wigner function: Random positions and momenta can be sampled from the Wigner function and evolved in parallel using classical equations of motion, giving x_η(t) and p_η(t) for trajectory η. Expectation values in TWA are then statistically approximated by $⟨ Q ⟩ \approx \frac{1}{n_{traj}} \sum_{η}^{n_{traj}} Q_{W} (x_{η} (t), p_{η} (t))$ ⁠, using n_traj trajectories.

Importantly, for small-spin systems, and in particular for spin-1/2 models as considered here, TWA can be drastically improved when using a sampling of the initial Wigner function using only a discrete set of initial phase points.⁵¹ Considering a system consisting of a single spin-1/2 described by the Pauli operators $σ = (σ_{x}, σ_{y}, σ_{z})$ ⁠, we define the corresponding phase space variables as $S = (S_{x}, S_{y}, S_{z})$ ⁠. One can then define discrete Wigner functions, which are only defined for the eight different discrete points $S_{0} = (\pm 1, \pm 1, \pm 1)$ ⁠. For example, taking a state of N_s spin-1/2 particles of the form $| ψ 〉 = \otimes_{i = 1}^{N_{s}} | {↓ 〉}_{i}$ ⁠, it is straightforward to show that any possible observable can be exactly described by sampling each spin from a discrete Wigner function with only non-zero values of $W_{i}^{↓} (S_{x} = \pm 1, S_{y} = \pm 1, S_{z} = - 1) = 1 / 4$ ⁠. Correspondingly, the state |↑⟩_i is exactly described by the discrete distribution with non-zero elements $W_{i}^{↑} (S_{x} = \pm 1, S_{y} = \pm 1, S_{z} = + 1) = 1 / 4$ ⁠.

Furthermore, it can be shown in general that equivalent discrete sampling strategies can lead to exact quantum state descriptions for general discrete D-level systems and for eigenstates of general spin-S operators, in the sense that the measurement statistics for any observable can be exactly reproduced from sampling the Wigner function.⁵³ Discrete sampling in combination with classical evolution is known as (generalized) discrete truncated Wigner approximation, (G)DTWA.^51,53 Classical equations of motion for the spin-variables can be derived by Wigner–Weyl transforming the Heisenberg equations of motion while factorizing the Weyl symbols into the phase space variables. (G)DTWA has been shown to capture quantum features in spin-model dynamics in several theory settings^{77,79,100–102} and in comparison with experiments.^{76,78,103,104}

Below we will consider a system consisting not only of spins but also of bosonic

a (a^{†})

degrees of freedom. For the bosonic part, we will consider the complex numbers a → A and a^† → A* as the classical phase space. We note that for additional bosonic degrees of freedom with operators denoted as

b (b^{†})

⁠, one can introduce a corresponding classical phase space with

b (b^{†}) \to B (B^{*})

(see, for example, Sec. V). Then, considering a system with N_s spin-1/2 particles coupled to N_a/b bosonic modes a/b, computing expectation values of an observable Q(σ_i, a_j, b_k) with TWA at time t corresponds to numerically evaluating

\begin{align} ⟨ O (t; {σ_{i}, a_{j}, b_{k}}) ⟩ & \approx \int \prod_{i = 1}^{N_{s}} \prod_{j = 1}^{N_{a}} \prod_{k = 1}^{N_{b}} d S_{0}^{i} d^{2} A_{0}^{j} d^{2} B_{0}^{k} \\ \times W_{i} (S_{0}^{i}) W_{j} (A_{0}^{j}) W_{k} (B_{0}^{k}) \\ \times O_{W} ({S_{cl}^{i} (t), A_{cl}^{j} (t), B_{cl}^{k} (t)}), \end{align}

(43)

where

d S_{0} = d S_{0}^{x} d S_{0}^{y} d S_{0}^{z}

⁠, d²A₀ = dReA₀dImA₀/π, and the subscript 0 indicates the initial values at t = 0. The classical variables for spin i and bosons j and k are sampled from the initial Wigner functions

W_{i} (S_{0}^{i})

⁠,

W_{j} (A_{0}^{j})

⁠, and

W_{k} (B_{0}^{k})

⁠, respectively. Note that we always assume an initial product state between all degrees of freedom such that the Wigner functions factorize. For the spins we will use the discrete distributions

W_{i}^{↑ / ↓}

defined earlier, while for the bosonic modes we use standard continuous Wigner functions, in particular

W_{j} (A_{0}^{j}) = \frac{1}{2 π w_{j}^{2}} \exp \{- | A_{j} - {\bar{α}}_{j} |^{2} / (2 w_{j}^{2})\},

(44)

where

w_{j}^{2} = ({\bar{n}}_{j} + 1 / 2) / 2

and α center of the Wigner function. For the vacuum state

{\bar{n}}_{j} = {\bar{α}}_{j} = 0

⁠, while for a coherent state |α⟩,

{\bar{n}}_{j} = 0

and

{\bar{α}}_{j} = α

(see Refs. 52 and 81 for more details).

O_{W} ({S_{cl}^{i} (t), A_{cl}^{j} (t), B_{cl}^{k} (t)})

is the Weyl symbol corresponding to the observable of interest. We use the subscript cl on the time-dependent variables

S_{cl}^{i} (t)

⁠,

A_{cl}^{j} (t)

⁠, and

B_{cl}^{k} (t)

to indicate that they obey the classical equations of motion for spin i, boson j, and boson k, respectively. In Sec. V, we will provide the full classical equations of motion for our problem of interest for examples of both closed and open system dynamics. (G)DTWA methods have been recently developed further to also include open-system dynamics under Lindblad master equations.^56–60 For our simulations, we follow the procedure in Ref. 59 and use a spherical coordinate parametrization of the phase space for spin i with phase point operators defined as

A_{i} (θ_{i}, ϕ_{i}) = \frac{1}{2} [1_{i} + s (θ_{i}, ϕ_{i}) \cdot σ_{i}],

(45)

where we use the vector on the surface of a sphere with radius

\sqrt{3}

⁠,

s_{i} = \sqrt{3} {(\sin θ_{i} \cos ϕ_{i}, - \sin θ_{i} \sin ϕ_{i}, - \cos θ_{i})}^{⊺}

⁠.

In Ref. 59, it was shown that for open spin-1/2 models, it is convenient to work with flattened Wigner functions of the form $χ_{i} (θ_{i}, ϕ_{i}) \equiv W_{i} (θ_{i}, ϕ_{i}) \frac{\sin θ_{i}}{2 π}$ ⁠. The equations of motion for an open system can then be found by deriving correspondence rules (reminiscent of Bopp representations for bosonic systems⁸¹), i.e., rules for mapping terms such as X_iρ_i, ρ_iX_i, or $X_{i} ρ_{i} X_{i}^{†}$ with Pauli operators $X = σ_{i}^{x, y, z}$ to the phase space, which leads to terms incorporating the four linearly independent differential operations $1, \frac{d}{d θ_{i}}, \frac{d}{d ϕ_{i}}, \frac{d^{2}}{d ϕ_{i}^{2}}$ acting on χ(θ_i, ϕ_i). It can be shown that the resulting EOMs lead to standard Fokker–Planck equations. This is only valid, without further approximations, for systems of non-interacting spins or if the initial state is a large coherent spin state. In these scenarios, the dynamics are given by the solution to the Fokker–Planck equations. Rather than solving these equations directly, we solve the corresponding Itô stochastic differential equations and average the relevant expectation values over many trajectories.⁹³

In our discrete sampling, we select the initial angles θ_i, ϕ_i according to the parametrization given in Eq. . However, rather than sampling from the discrete

W_{i}^{↑ / ↓}

Wigner distribution, we sample from a slightly modified flattened Wigner function

χ^{↑ / ↓} (θ, ϕ) = \frac{1}{2 π} δ (θ \pm \arccos \frac{1}{\sqrt{3}}),

(46)

which is generated by rotating the discrete Wigner function

W_{i}^{↑ / ↓}

around the z-axis. This initial Wigner function is uniform in ϕ and, therefore, guarantees that

\frac{\partial^{2}}{\partial ϕ_{m} \partial (θ_{n} - θ_{n + 1})}

cross diffusion terms, which are dropped in TWA, vanish at t = 0.⁵⁹ For more details, we refer the reader to Ref. 59 and Sec. V, where the detailed equations for our problem of interest are introduced.

V. RESULTS

To evaluate the performance of our two numerical methods, we consider the disordered Holstein–Tavis–Cummings model describing a system of N_s spins coupled to a single bosonic mode (a), representing a coupling to a cavity mode, and N_s local vibrational modes (b_i) associated with each spin. The system is described by the Hamiltonian,³³

H = H_{TC} + H_{vib} + H_{H} + H_{dis},

(47)

where

H_{TC} = \frac{Δ}{2} \sum_{j = 1}^{N_{s}} (σ_{j}^{z} + 1) + \frac{g}{\sqrt{N_{s}}} \sum_{j = 1}^{N_{s}} a σ_{j}^{+} + a^{†} σ_{j}^{-},

(48a)

H_{vib} = ν \sum_{j = 1}^{N_{s}} b_{j}^{†} b_{j},

(48b)

H_{H} = - \frac{λ ν}{2} \sum_{j = 1}^{N_{s}} (b_{j} + b_{j}^{†}) (σ_{j}^{z} + 1),

(48c)

H_{dis} = - \frac{1}{2} \sum_{j = 1}^{N_{s}} ϵ_{j} (σ_{j}^{z} + 1),

(48d)

where Δ describes the detuning of the spin transition frequency relative to the cavity mode,

g / \sqrt{N_{s}}

the single-spin coupling to the cavity, ν the frequency of the vibrational modes, λ the relative strength of the Holstein coupling, and ϵ_j the disorder in the transition frequency for spin j.

The dynamics of the Holstein–Tavis–Cummings model, in particular in the presence of disorder, have importance, e.g., in the field of polaritonic chemistry.¹⁰⁵ It has been previously studied using a matrix product state method³³ and also using a similar non-Gaussian state framework to the one discussed here.⁷¹ By tuning the relative strength of the various terms, this Hamiltonian can be reduced to spin–boson Hamiltonians applicable, e.g., to trapped ion quantum simulators, impurity models, and quantum chemistry. The relative strength of g, λ, and ν allows us to go from the weak coupling regime between the spin and bosonic degrees of freedom to a model governed by a Tavis–Cummings type interaction to a Holstein coupling or a combination thereof. We schematically illustrate the performance of the two numerical methods in each of these regimes in Fig. 5(a).

FIG. 5.

View large Download slide

Schematic of the performance of TWA and NGS for (a) closed and (b) open dynamics. In the present study, NGS is limited to at most N_p = 16 multi-mode coherent states and no squeezing. (a) The performance in g − λ parameter space for closed dynamics. When g, λ ≲ 1, TWA performs well, while when g ≲ 1, λ ≳ 1, NGS is the better choice. When g ≳ 1, TWA captures short-time dynamics but can produce incorrect mid-to late-time results. In comparison, NGS typically does not capture quantitative details beyond the first spin relaxation but does provide qualitative insights into the dynamics by correctly capturing the magnitude of the persistent spin–cavity dynamics. (b) Open dynamics. We operate in the g ≲ 1, λ ≳ 1 region, so NGS outperforms TWA at small κ. At large κ, TWA performs well again, with both NGS and TWA agreeing. When evolving under only H_TC with the spins decaying collectively at rate Γ, NGS and TWA agree when using the large N_s Holstein–Primakoff (HP) transformation at small Γ. At large Γ the collective spin quickly decays to the ground state, which is inaccessible as the first order expansion of the root in HP is no longer sufficient near the ground state since it was expanded around the excited state (small $a_{s}^{†} a_{s}$ ⁠), so neither method is able to capture intermediate-to-late-time dynamics.

Furthermore, we investigate how both methods perform in the presence of sources of decoherence. In this case, we are interested in the evolution of the density matrix ρ as described by the Lindblad master equation given in Eq. (20). We study open dynamics with the following types of Lindblad jump operators:

Cavity decay (rate κ): $c_{m} = \sqrt{κ} a$ ⁠.
Single spin decay (rate γ): $c_{m} = \sqrt{γ} σ_{i}^{-}$ ⁠.
Collective spin decay (rate Γ): $c_{m} = \sqrt{Γ} \sum_{i} σ_{i}^{-}$ ⁠.

We summarize the performance of the two methods in the presence of these decoherence sources schematically in Fig. 5(b).

A. Details of NGS and TWA simulations

In principle, the NGS ansatz can be used directly to treat the spin degrees of freedom. However, to avoid the explicit exponential scaling with N_s, we use a Holstein–Primakoff (HP) transformation to map each spin-1/2 to a bosonic mode. We use the following form of the HP transformation:¹⁰⁶

S_{+} = (1 - a_{s}^{†} a_{s}) a_{s},

(49a)

S_{z} = 1 / 2 - a_{s}^{†} a_{s},

(49b)

which is exact for spin-1/2, and where we use the s subscript to denote a bosonic operator acting on the spin degree of freedom. As such, the dynamics are restricted to the |0⟩ ≡|e⟩ and |1⟩ ≡|g⟩ subspace, with

| g 〉 〈 g | = a_{s}^{†} a_{s}

and

| e 〉 〈 e | = 1 - a_{s}^{†} a_{s}

⁠. The vacuum is a Gaussian state and can, therefore, be described using only N_p = 1. The Fock state |1⟩ can be described using N_p = 2, by

| 1 〉 = 1 / N \lim_{α \to 0} (| α 〉 - | - α 〉)

⁠, where

N

is a normalization factor,¹⁰⁷ with α ∼ 0.001 sufficient to describe the state with high fidelity. As such, each spin can be captured using N_p = 3 coherent states.

For the NGS simulations, we employ the NGS ansatz with N_b = 2N_s + 1 modes. We include cavity decay at strength γ using the simple parameter update prescription outlined in Sec. V. Note that, as a consequence of the HP mapping, the ansatz is not well-suited to some scenarios. For example, there is no spin–spin coupling when evolving under only H_dis, so an unentangled initial state will evolve as a tensor product of single spins, each requiring N_p = 3. The number of coherent states therefore scales exponentially in the number of spins, $3^{N_{s}}$ ⁠. This limitation could potentially be mitigated by modifying the ansatz to be a superposition of squeezed displaced Fock states.⁵⁵

For the TWA simulations, to more accurately capture the dynamics of the cavity and vibrational modes, we extend the set of TWA equations by including the classical equations of motion for the mode excitation numbers

a^{†} a \to N_{A}, b_{k}^{†} b_{k} \to N_{B_{k}}

⁠. Including all three sources of decoherence described earlier, the stochastic equations of motion for the spin degrees of freedom are

\begin{align} {\dot{θ}}_{i, c l} & = (Γ + γ) (\cot θ_{i, c l} - \frac{\csc θ_{i, c l}}{\sqrt{3}}) + \frac{2 g}{\sqrt{N_{s}}} I m [e^{- i ϕ_{i, c l}} A_{cl}] \\ - \frac{Γ \sqrt{3}}{2} \sum_{j} \cos (ϕ_{i, c l} - ϕ_{j, c l}) \sin θ_{j, c l}, \end{align}

(50a)

\begin{align} d ϕ_{i, c l} & = (- Δ + ϵ_{n} + 2 λ ν R e [B_{i, c l}] - \frac{2 g}{\sqrt{N_{s}}} \cot θ_{i, c l} R e [A_{cl} e^{- i ϕ_{i, c l}}] \\ - \frac{Γ \sqrt{3}}{2} \sum_{j} \cot θ_{i, c l} \sin θ_{j, c l} \sin (ϕ_{i, c l} - ϕ_{j, c l})) d t \\ + \sqrt{(Γ + γ) f (θ_{i, c l})} d W_{ϕ_{i}}, \end{align}

(50b)

where we introduced the function

f (θ_{i, c l}) = 1 + 2 \cot {θ_{i, c l}}^{2} - 2 \cot θ_{i, c l} \csc θ_{i, c l} / \sqrt{3}

⁠. The stochastic behavior of the equations of motion is generated by the Wiener increments

d W_{ϕ_{i}}

acting on ϕ_i. For each timestep dt, the Wiener increment

d W_{ϕ_{i}}

for each angle is independently drawn from a normal distribution with zero mean and a variance of dt.

The equations of motion for the bosonic degrees of freedom are given by

\dot{A_{cl}} = - \frac{κ}{2} A_{cl} - i \frac{\sqrt{3} g}{2 \sqrt{N_{s}}} \sum_{i} e^{i ϕ_{i, c l}} \sin θ_{i, c l},

(51a)

{\dot{N}}_{A, c l} = - κ N_{A, c l} - \frac{\sqrt{3} g}{\sqrt{N_{s}}} \sum_{i} \sin θ_{i, c l} I m [A_{cl} e^{- i ϕ_{i, c l}}],

(51b)

{\dot{B}}_{k, c l} = - i ν B_{k, c l} + i \frac{λ ν}{2} (1 - \sqrt{3} \cos θ_{k, c l}),

(51c)

{\dot{N}}_{B_{k, c l}} = λ ν (1 - \sqrt{3} \cos θ_{k, c l}) I m [B_{k, c l}] .

(51d)

The initial state used for all NGS simulations includes a small amount of randomness for each variational parameter to break the degeneracy of the Gaussian states. We draw random values from a uniform distribution between (0, 10⁻⁴). We find that this is sufficient to ensure each Gaussian state in the superposition evolves independently while preserving extremely small infidelity with the true initial state, $1 - F (t = 0) < O (1 0^{- 6})$ ⁠, as seen in Fig. S1 of the supplementary material.

B. Closed system dynamics

In this section, we compare the performance of the two numerical methods. We consider a system with N_s = 3 spins, the corresponding three phonon modes, and one cavity mode. For this system size and suitable initial states and Hamiltonian parameters, choosing a moderate Fock state truncation $\sim 10$ allows us to compare our results to exact numerics. Our initial state is spins polarized up, the vibrational modes in the vacuum, and the cavity in a coherent state, $| ψ (0) 〉 = | {↑ 〉}^{\otimes N_{s}} | α = {1 〉}_{a} | {0 〉}_{b}^{\otimes N_{s}}$ ⁠. Note that, in contrast to Ref. 71, our initial state is a superposition of several excitation manifolds, precluding further simplifications to the ansatz.

For all closed dynamics simulations, we set Δ = 0. For evolution under the Holstein–Tavis–Cummings model in Eq. (47), we find that while both methods are generally able to capture the short time dynamics, at later times they outperform one another in different parameter regimes. We summarize our findings in Fig. 5(a), where we qualitatively depict the performance for different parameter regimes in the absence of decoherence. More detailed dynamics for each regime are plotted in Figs. 6–9, with the first and second rows of each figure showing dynamics without and with disorder ϵ, respectively.

FIG. 6.

View large Download slide

Closed dynamics: g = λ = 0.1, Δ = 0. Top row: Without disorder. Both NGS and TWA capture the initial spin relaxation, but NGS incorrectly predicts a slower revival. Bottom row: Disordered, $\vec{ϵ} = [2 g, 3 g, 4 g]$ ⁠. Here, TWA captures the dynamics more accurately compared to NGS. NGS is with N_p = 4, TWA is with n_traj = 10⁴ with standard error shaded.

FIG. 7.

View large Download slide

Closed dynamics: g = 0.1, λ = 1, Δ = 0. Top row: Without disorder. Neither TWA nor NGS completely capture the spin–cavity observables, but NGS more closely tracks the dynamics. Bottom row: Disordered, $\vec{ϵ} = [2 g, 3 g, 4 g]$ ⁠. NGS performs well, capturing all dynamics with small error. Interestingly, here TWA overestimates the magnitude of changes in the spin–cavity observables. Here NGS uses N_p = 12, TWA is with n_traj = 10⁴ with standard error shaded.

FIG. 8.

View large Download slide

Closed dynamics: g = 1, λ = 0.1, Δ = 0. Top row: Without disorder. Beyond initial spin relaxation, both TWA and NGS perform relatively poorly. TWA incorrectly predicts equilibration of the spin–cavity dynamics. NGS does continue to produce spin–cavity dynamics but overestimates the magnitude of the oscillations. Bottom row: Disordered, $\vec{ϵ} = [2.6 g, 3.2 g, 4.2 g]$ ⁠. TWA incorrectly predicts equilibration of both spin and cavity observables after the first oscillation, while NGS produces qualitatively correct dynamics, even with only N_p = 4 coherent states. Here NGS uses N_p = 2 for the first row and N_p = 4 for the second row, and TWA is with n_traj = 10⁴ with standard error shaded.

FIG. 9.

View large Download slide

Closed dynamics: g = 1, λ = 1, Δ = 0. Top row: Without disorder. TWA accurately captures the dynamics at early times as compared to the NGS. Bottom row: Disordered, $\vec{ϵ} = [2.6 g, 3.2 g, 4.2 g]$ ⁠. TWA correctly captures key oscillations, even at late times. NGS also captures some of these features, but less accurately. Both methods perform similarly for the vibrational dynamics. Here NGS uses N_p = 4, TWA is with n_traj = 10⁴ with standard error shaded.

We begin in the weak coupling regime g = λ = 0.1, shown in Fig. 6. Here, the dynamics is slow on the considered time-scale. In both the disorder-free (first row) and disordered (second row) settings, NGS and TWA accurately capture the small fluctuations of the vibrational modes at all times. Both methods also capture the initial spin relaxation; however, NGS misses the revival time. TWA’s ability to capture the first oscillation appears universally in all of the parameter regimes considered in this work. This behavior can be understood as follows: due to our choice of a factorisable initial state with a corresponding positive semi-definite Wigner function, the TWA sampling is able to reproduce the initial state and the mean-field and low order correlations that are generated during the short-time dynamics.

A generic feature of TWA is that when extending into the medium-to-long-time dynamics, the potential buildup of higher order correlations is not captured by the method. While this is not visible in Fig. 6, it can be seen clearly in the figures corresponding to the regimes discussed below.

The second regime we consider is the strong spin–vibrational coupling λ ≫ g, shown in Fig. 7. Here, we expect NGS to perform well, as NGS is exact with any coherent state number N_p for H_H. In the disorder-free regime (top row), although NGS with N_p = 12 does not fully capture the dynamics at late times, it does outperform TWA, which majorly underestimates the spin decay. In this case, the lack of disorder is challenging for our NGS ansatz: each spin evolves identically, requiring the superposition of coherent states to be factored into a product. This symmetry is broken when introducing disorder (second row), and we see that NGS captures accurately the dynamics for all considered times, including the initial decay and then revival of ⟨S_z⟩. For TWA, although the numerics match better for the vibrational dynamics in the presence of disorder, this is primarily a consequence of the fact that the disorder causes the Holstein interactions to dominate, and the spin and cavity dynamics continue to disagree with the exact solution.

Third, we move to the strong spin–cavity coupling regime, g ≫ λ, shown in Fig. 8. After accurately capturing the first oscillation, the TWA spin dynamics equilibrate about ⟨S_z⟩∼0, unlike the exact dynamics which, although they do oscillate about ⟨S_z⟩∼0, exhibit persistent oscillations with magnitude ⟨S_z⟩∼1/2. Similarly, NGS with N_p = 4 coherent states fails to capture any of the spin, vibration, or cavity dynamics. This is unsurprising because, in this regime where the Tavis–Cummings term dominates, we expect the number of coherent states required to scale as $\sim 3^{N_{s}}$ ⁠, as each spin must be described using N_p = 3 coherent states. Although increasing N_p would eventually improve the accuracy, we found that increasing it up to N_p ≲ 16 did not provide a substantial increase in accuracy while increasing the computational cost. In principle, TWA does not suffer from the same problem. However, if one needs to access higher order correlations with TWA, the introduction of higher order cumulants and the BBGKY hierarchy may become necessary. This poses an analogous problem to the number of Gaussian states in the ansatz: an exponentially increasing number of equations and potential numerical instabilities.¹⁰⁸

Fourth, we consider the strong spin–cavity and spin–vibration regime λ = g = 1, shown in Fig. 9. Without disorder, both methods struggle to capture the dynamics at late times, although TWA in particular is able to capture qualitative features with reasonable accuracy. Introducing disorder breaks the collective nature of the spins, enabling both methods to more accurately track the dynamics. TWA is able to qualitatively reproduce the periodic peaks in the spin and cavity dynamics at even later times than NGS. Both methods correctly obtain the vibrational dynamics.

Finally, we note that an advantage of NGS is the accessibility of the wave function. This means that any desired quantity, including entanglement entropy, can be computed. Furthermore, for small systems, strict performance measures such as fidelity can be easily computed. These are shown in the supplementary material for the four different coupling regimes (g, λ) considered in Figs. 6–9. The analogous plots for TWA cannot be generated.

C. Open system dynamics

1. Holstein Tavis–Cummings

Next, we introduce decoherence to our simulations. Figures 10 and 11 show the spin and bosonic dynamics for N_s = 3 in the presence of cavity loss a at strengths κ = g (Fig. 10) and κ = 10g (Fig. 11), and in the regime λ ≫ g, where NGS provides more accurate predictions in the closed setting (cf. Fig. 7). In the top row we set the disorder ϵ = 0, while the bottom row shows the dynamics in the presence of disorder, $\vec{ϵ} = [2 g, 3 g, 4 g]$ ⁠. Although exact dynamics were accessible for a closed system of this size, obtaining exact numerics in the open dynamics setting is challenging. In the supplementary material, we compare the two methods against exact numerics for a smaller system of N_s = 1. (supplementary material)

FIG. 10.

Open dynamics. g = 0.1, λ = 1, Δ = 0, and cavity decay κ = g. Top row: Without disorder. Bottom row: With disorder ϵ⃗=[2g,3g,4g]. NGS and TWA capture the same short-time dynamics but disagree beyond νt ∼ 5. Based on the corresponding closed dynamics results in Fig. 7, we expect NGS to be more reliable in this regime. Exact numerics is challenging for open dynamics of systems of this size (see Fig. S2) for benchmarking of smaller systems. Here, NGS uses Np = 8 and ntraj = 40 with standard error indicated by the error bars; TWA uses ntraj = 104 with standard error shaded.

View large Download slide

Open dynamics. g = 0.1, λ = 1, Δ = 0, and cavity decay κ = g. Top row: Without disorder. Bottom row: With disorder $\vec{ϵ} = [2 g, 3 g, 4 g]$ ⁠. NGS and TWA capture the same short-time dynamics but disagree beyond νt ∼ 5. Based on the corresponding closed dynamics results in Fig. 7, we expect NGS to be more reliable in this regime. Exact numerics is challenging for open dynamics of systems of this size (see Fig. S2) for benchmarking of smaller systems. Here, NGS uses N_p = 8 and n_traj = 40 with standard error indicated by the error bars; TWA uses n_traj = 10⁴ with standard error shaded.

FIG. 11.

Open dynamics. g = 0.1, λ = 1, Δ = 0, and cavity decay κ = 10g. Top row: Without disorder. Bottom row: With disorder ϵ⃗=[2g,3g,4g]. We observe compatible results for the spins and a remarkable agreement between NGS and TWA for the vibrational and cavity dynamics, including the small amplitude oscillations in ⟨ncav⟩ at late times. Here NGS uses Np = 8 and ntraj = 40 with standard error indicated by the error bars; TWA uses ntraj = 104 with standard error shaded.

View large Download slide

Open dynamics. g = 0.1, λ = 1, Δ = 0, and cavity decay κ = 10g. Top row: Without disorder. Bottom row: With disorder $\vec{ϵ} = [2 g, 3 g, 4 g]$ ⁠. We observe compatible results for the spins and a remarkable agreement between NGS and TWA for the vibrational and cavity dynamics, including the small amplitude oscillations in ⟨n_cav⟩ at late times. Here NGS uses N_p = 8 and n_traj = 40 with standard error indicated by the error bars; TWA uses n_traj = 10⁴ with standard error shaded.

For these parameters, we find that, perhaps unsurprisingly, NGS continues to perform well in both κ = g and κ = 10g decoherence regimes. In the weaker decay limit shown in Fig. 10, NGS and TWA agree only at short times. The under and over-estimation of spin–cavity dynamics by TWA in the non-disordered and disordered systems, respectively, is consistent with the behavior of TWA in the closed system (see Fig. 7). In the large decay limit shown in Fig. 11, NGS and TWA are in reasonable agreement with one another for the spin dynamics and in near total agreement for the vibrational and cavity dynamics. Both show fast decay of the cavity to the vacuum and remarkably both capture small oscillatory dynamics at late times with excellent agreement. Physically, both methods demonstrate that strong cavity decay stabilizes the spin dynamics, which we attribute to the reduction of the effective Tavis–Cummings coupling strength and the prevention of the buildup of correlations in the system between the spins, as well as spin–boson correlations, due to the loss of cavity excitations.

2. Tavis–Cummings

Next, we consider the effect of spin decay. We set ϵ, λ = 0 in Hamiltonian (47) such that evolution is only under H_TC and include collective spin decay at rate Γ. We set Δ = 1 and g = 0.1. Within the TWA formalism, we treat the dynamics using two methods. First, we continue to treat the system as a collection of individual spins, as described in Sec. IV. In Fig. 12, we plot the resulting dynamics, where TWA (CD) refers to implementing this sampling in the presence of collective spin decay at $Γ = 0.1 g / \sqrt{N_{s}}$ and $g / \sqrt{N_{s}}$ ⁠. TWA (SSD) uses the same sampling, but the decay mechanism is single spin at the corresponding rate. The agreement between the two and the disagreement with exact numerics (EXA) indicate that treating the spins individually with either of these two methods is inadequate to simulate collective spin decay.

FIG. 12.

Collective spin dynamics. ν = ϵ = λ = 0, Δ = 1, g = 0.1. Top row: Γ=0.1g/Ns. NGS (HP) and TWA (HP) use the Holstein–Primakoff representation of the collective spin and agree excellently with exact numerics (EXA). On the other hand, when TWA treats the spins individually, the collective decay dynamics TWA (CD) agrees with the single spin decay dynamics TWA (SSD), but both deviate from the exact evolution due to TWA’s failure to capture the correlations. Bottom row: Γ=g/Ns. NGS (HP) and TWA (HP) agree with exact numerics (EXA) until tν ∼5, when the first order Taylor series expansion of the HP mapping breaks down. Here, NGS uses Np = 8 and ntraj = 40 with standard error indicated by the error bars; TWA uses ntraj = 104 with standard error shaded.

View large Download slide

Collective spin dynamics. ν = ϵ = λ = 0, Δ = 1, g = 0.1. Top row: $Γ = 0.1 g / \sqrt{N_{s}}$ ⁠. NGS (HP) and TWA (HP) use the Holstein–Primakoff representation of the collective spin and agree excellently with exact numerics (EXA). On the other hand, when TWA treats the spins individually, the collective decay dynamics TWA (CD) agrees with the single spin decay dynamics TWA (SSD), but both deviate from the exact evolution due to TWA’s failure to capture the correlations. Bottom row: $Γ = g / \sqrt{N_{s}}$ ⁠. NGS (HP) and TWA (HP) agree with exact numerics (EXA) until tν ∼5, when the first order Taylor series expansion of the HP mapping breaks down. Here, NGS uses N_p = 8 and n_traj = 40 with standard error indicated by the error bars; TWA uses n_traj = 10⁴ with standard error shaded.

This motivates our second strategy. Because the Hamiltonian is H_TC, the dynamics takes place in the collective Tavis–Cummings manifold as H_TC conserves the total excitation number

N_{ex} = a^{†} a + \sum_{i} σ_{i}^{z}

⁠. Then, we can use the Holstein–Primakoff (HP) transformation to map the collective spin S to a single bosonic mode,

S_{+} = \sqrt{2 s - a_{s}^{†} a_{s}} a_{s},

(52a)

S_{z} = s - a_{s}^{†} a_{s},

(52b)

where s ≡ N_s/2. We use a Taylor series to expand the square root in powers of 1/s to first order. The NGS simulation then proceeds using the NGS ansatz with N_b = 2 bosonic modes. We can include collective spin decay

\sum_{i} σ_{i}^{-}

at strength Γ using the manifold projection technique described in Sec. III C.

For the TWA simulations, the equations of motion for the cavity mode A and the large spin HP mode B are

{\dot{A}}_{cl} = - \frac{κ}{2} A_{cl} - i \frac{g}{2 \sqrt{2 s}} B_{cl}^{*} (4 s - N_{B, c l}),

(53a)

{\dot{N}}_{A, c l} = - κ N_{A, c l} - \frac{g}{\sqrt{2 s}} I m [A_{cl} B_{cl}] (4 s - N_{B, c l}),

(53b)

\begin{align} {\dot{B}}_{cl} & = Γ B_{cl} [s - \frac{1}{8 s} + N_{B, c l} (\frac{| B_{cl} |^{2}}{16 s} - \frac{1}{2})] \\ + \frac{i g}{2 \sqrt{2 s}} [A_{cl} B_{cl}^{2} - 2 A_{cl}^{*} (2 s - N_{B, c l})] + i Δ B_{cl}, \end{align}

(53c)

\begin{align} {\dot{N}}_{B, c l} & = Γ [2 S + N_{B, c l} (\frac{1}{4 s} + 2 (s - 1) + F (B_{cl}))] \\ - \frac{g}{\sqrt{2 s}} I m [A_{cl} B_{cl}] (4 s - N_{B, c l}), \end{align}

(53d)

where we introduce the function

F (B_{cl}) = \frac{1}{8 s} | B_{cl} |^{2} (4 - 8 s + | B_{cl}^{2} |)

⁠. We sample the cavity (mode A) assuming a coherent state |α = 1⟩ and sample the large spin (mode B) assuming the spins are polarized pointing up, which corresponds to sampling the vacuum for mode B.

Our results using this approach, including collective spin decay, are shown in Fig. 12 in the lines labeled as TWA (HP) and NGS (HP). In the weak spin decay limit $Γ = 0.1 g / \sqrt{N_{s}}$ ⁠, both NGS (HP) and TWA (HP) are in excellent agreement with the exact numerics (EXA). NGS in particular captures the cavity dynamics with little error, while the error bars on the spin dynamics are still somewhat large due to our use of relatively few trajectories, n_traj = 40. In the large decay limit $Γ = g / \sqrt{N_{s}}$ ⁠, both NGS and TWA correctly capture the rapid spin decay until νt ∼5. Beyond this point, the first order Taylor series expansion of the HP mapping breaks down as highly excited Fock states are populated. One can potentially circumvent this issue by simulating collective spin decay, as was performed in Ref. 80 with DTWA. There, the spins were collectively coupled to a single cavity whose cavity loss was much stronger than the collective spin–cavity coupling, resulting in effective collective spin decay and without the utilization of the HP mapping. Comparing the performance of these approaches in different parameter regimes and for different models, e.g., for more complex forms of collective spin decay, represents an interesting future direction.

VI. CONCLUSIONS AND OUTLOOK

Summary and conclusions: In this work, we presented a non-Gaussian variational ansatz approach to studying the dynamics of open quantum systems composed of spin and bosonic degrees of freedom. While several other works in recent years have utilized NGS to study the time evolution of open quantum systems, previous efforts have focused on developing an equivalent ansatz for the density matrix and simulating the Lindblad equations. Here, we utilized the quantum trajectories method, allowing us to take advantage of the previously developed machinery and analytic expressions obtained for real- and imaginary-time dynamics.

In addition to providing a comprehensive overview of this method, we performed extensive numerical simulations over a broad range of parameters of a spin–boson Hamiltonian [Eq. (47)] with Tavis–Cummings (TC) and Holstein couplings, which is applicable to a broad range of quantum simulation platforms as well as problems of interest in quantum chemistry, atomic physics, and condensed matter theory. We compared the performance of NGS with a method using the truncated Wigner approximation for systems with mixed bosonic and spin degrees of freedom, extended to open quantum systems following the approach in Ref. 59.

In the absence of decoherence, our findings are as follows: for strong TC coupling, TWA is the more accurate method, while for strong Holstein couplings, NGS is the better choice. When neither term dominates, for both weak and strong coupling regimes, TWA captures the short-time dynamics, while NGS generally displays the correct qualitative behavior, even at late times. After introducing spin disorder, the performance of NGS typically improves, while for TWA the vibrational dynamics match the exact dynamics better, which is attributable to the fact that disorder causes the Holstein interaction to dominate.

For open quantum dynamics, we focused on the regime where the Holstein term dominates and considered the effect of cavity loss. At weak decay rates, the NGS continues to perform well. TWA improves as the cavity decay rate increases due to the loss of quantum correlations, with both NGS and TWA showing excellent agreement. In the presence of collective spin decay, we considered the TC model only, finding that in the limit of small collective decay, both methods perform well when using a Holstein–Primakoff transformation for the large spin. In the limit of large collective decay, NGS and TWA both only capture the short time dynamics, as the Holstein–Primakoff transformation is no longer accurate at later times. Using TWA, we were able to also treat each spin individually. However, TWA does not capture the collective nature of the decay, with the results closely matching the effect of single spin decay.

Further considerations should also be made when deciding between the two methods. Although the NGS ansatz can be made less computationally demanding by reducing the number N_p of coherent states, in general, TWA methods are easier to implement and require fewer computational resources. The resource requirement and the complexity of NGS are offset by the advantages that it is a controlled approximation and gives access to the wave function, allowing one to access any observable, including higher-order correlations. The TWA framework needs to be amended if one hopes to capture these correlations accurately. A potential strategy to remedy this can be to use a cluster TWA approach.¹⁰⁹ There, several sub-system parts are grouped into a single (discrete or continuous) large sub-system that, provided a proper sampling strategy, follows the fully exact quantum evolution, while correlations between clusters are approximated in TWA. Such approaches allow to use the cluster size as a controllable parameter to enhance the simulation toward the exact one.

Outlook: The extension of NGS to open quantum systems using the quantum trajectories formulation that we presented in this work is perhaps the most natural pathway. For Hermitian jump operators, an alternative approach and a simple extension of this work would be to instead solve the stochastic partial differential equations that result from the unraveling of the master equation. Furthermore, one could explore the impact of the chosen unraveling on the performance of the NGS method at a fixed N_p, as it has been shown that this choice can have a large effect on the entanglement buildup in the trajectory.^110,111

Another limitation of the present formulation is that the spin states in Eq. (1) are exact and span the full spin Hilbert space of dimension $2^{N_{s}}$ ⁠, thus limiting the use of the ansatz to a handful of spins unless approximations such as the large N_s expansion in the Holstein–Primakoff mapping used in Sec. V C 2 are invoked. On the one hand, studies using a fermionic Gaussian state representation of spins have been performed,¹¹² and these might be combined in principle in a straightforward way with the NGS ansatz for bosons.⁴⁹ On the other hand, it would be highly interesting to combine the NGS ansatz with other variational techniques highly suitable for the spins such as tensor network based approaches.^39–41 Furthermore, similar to the present comparison between NGS and TWA, it would be beneficial to apply the here-presented non-Gaussian ansatz to the study of other systems, which might be challenging to simulate otherwise. These include, for instance, purely bosonic models, such as the Bose–Hubbard model, with disorder and on non-regular lattices.¹¹³ Such extensive studies will allow for a comparison between our approach and the corresponding master equation approach based on extending the ansatz of Eq. (1) to density matrices.^54,55

Finally, we note that a particular promising application field of the methodologies introduced here could be in the emerging field of polaritonic chemistry.^114–117 Recent experiments have demonstrated that large collective strong cavity couplings (e.g., $g_{c} = g \sqrt{N}$ ⁠) can be functionalized to modify chemical reactivity. A theoretical understanding for such modifications is currently centered around the question of how a delocalized polaritonic state can play a role in changing chemistry on the single-molecule level as local amplitudes of collective polaritons vanish in the thermodynamic limit. In spin–boson approximations to the problems, in particular for the disordered Holstein–Tavis–Cummings¹⁰⁵ model that we studied here, it was recently discovered that the interplay of disorder and collective cavity couplings can give rise to robust local quantum effects in the large-N limit in the form of non-Gaussian distributions of the nuclear coordinate.³³ Using matrix product state methods, it was possible to push simulations to systems with 160 effective molecules, but in particular the TWA approach discussed here would allow access to much larger systems. Typical parameter regimes discussed here are well covered by the TWA approach (e.g., the typical parameters from Ref. 33 correspond to λ ∼ 0.1ν and g ≪ ν) and, therefore, hint to a general applicability of the method in the relevant regime. Furthermore, the TWA approach can be straightforwardly adapted to simulate nuclear dynamics not only on harmonic but also arbitrary potential energy surfaces, while the NGS ansatz using the machinery presented here can also model anharmonic Hamiltonian terms. In the future, this may allow for the analysis of quantum effects in realistic chemical reaction models, even in a macroscopic limit.

SUPPLEMENTARY MATERIAL

In the supplementary material, we show that the tangent space of the variational manifold for the NGS ansatz is a Kähler manifold and benchmark NGS and TWA against exact numerics for a single spin Ns = 1.

ACKNOWLEDGMENTS

A.S.N., B.G., L.J.B., and J.M. are supported by the Dutch Research Council (NWO/OCW) as part of the Quantum Software Consortium program (Project No. 024.003.037), Quantum Delta NL (Project No. NGF.1582.22.030), and ENW-XL grant (Project No. OCENW.XL21.XL21.122). J.M. was partly supported by a Proof-of-Concept grant through the Innovation Exchange Amsterdam (IXA) and the NWO Take-off Phase 1 grant (Project No. 20593). J.T.Y. is supported by the NWO Talent Program (Project No. VI.Veni.222.312), which is (partly) financed by the Dutch Research Council (NWO). J.S. is supported by the Interdisciplinary Thematic Institute QMat, as part of the ITI 2021–2028 program of the University of Strasbourg, CNRS and Inserm, and was supported by IdEx Unistra (Grant No. ANR-10-IDEX-0002), SFRI STRAT’US project (Grant No. ANR-20-SFRI-0012), and EUR QMAT Grant No. ANR-17-EURE-0024 under the framework of the French Investments for the Future Program. The work was supported by the CNRS through the EMERGENCE@INC2024 project DINOPARC and by the French National Research Agency under the Investments of the Future Program Project No. ANR-21-ESRE-0032 (aQCess).

AUTHOR DECLARATIONS

Conflict of Interest

The authors have no conflicts to disclose.

Author Contributions

Liam J. Bond: Conceptualization (equal); Data curation (equal); Investigation (lead); Methodology (equal); Software (equal); Writing – original draft (equal); Writing – review & editing (equal). Bas Gerritsen: Conceptualization (equal); Data curation (equal); Investigation (equal); Methodology (equal); Software (equal); Writing – original draft (equal); Writing – review & editing (equal). Jiří Minář: Conceptualization (equal); Methodology (equal); Supervision (lead); Writing – original draft (equal); Writing – review & editing (equal). Jeremy T. Young: Conceptualization (equal); Methodology (equal); Supervision (supporting); Writing – original draft (equal); Writing – review & editing (equal). Johannes Schachenmayer: Conceptualization (equal); Methodology (equal); Supervision (supporting); Writing – review & editing (equal). Arghavan Safavi-Naini: Conceptualization (equal); Methodology (equal); Supervision (lead); Writing – original draft (equal); Writing – review & editing (equal).

DATA AVAILABILITY

The data that support the findings of this study are available from the corresponding author upon reasonable request.

REFERENCES

1.

J.

Preskill

,

Quantum

2

,

79

(

2018

); arXiv:1801.00862.

https://doi.org/10.22331/q-2018-08-06-79

2.

B.

Peropadre

,

D.

Zueco

,

D.

Porras

, and

J. J.

García-Ripoll

,

Phys. Rev. Lett.

111

,

243602

(

2013

).

https://doi.org/10.1103/physrevlett.111.243602

Google Scholar

Crossref

PubMed

3.

F.

Yoshihara

,

T.

Fuse

,

S.

Ashhab

,

K.

Kakuyanagi

,

S.

Saito

, and

K.

Semba

,

Nat. Phys.

13

,

44

(

2017

).

https://doi.org/10.1038/nphys3906

Google Scholar

Crossref

4.

P.

Forn-Díaz

,

J. J.

García-Ripoll

,

B.

Peropadre

,

J.-L.

Orgiazzi

,

M.

Yurtalan

,

R.

Belyansky

,

C. M.

Wilson

, and

A.

Lupascu

,

Nat. Phys.

13

,

39

(

2017

).

https://doi.org/10.1038/nphys3905

Google Scholar

Crossref

5.

L.

Magazzù

,

P.

Forn-Díaz

,

R.

Belyansky

,

J.-L.

Orgiazzi

,

M.

Yurtalan

,

M. R.

Otto

,

A.

Lupascu

,

C.

Wilson

, and

M.

Grifoni

,

Nat. Commun.

9

,

1403

(

2018

).

https://doi.org/10.1038/s41467-018-03626-w

Google Scholar

Crossref

PubMed

6.

M.

Marcuzzi

,

J. c. v.

Minář

,

D.

Barredo

,

S.

de Léséleuc

,

H.

Labuhn

,

T.

Lahaye

,

A.

Browaeys

,

E.

Levi

, and

I.

Lesanovsky

,

Phys. Rev. Lett.

118

,

063606

(

2017

).

https://doi.org/10.1103/physrevlett.118.063606

Google Scholar

Crossref

PubMed

7.

F. M.

Gambetta

,

W.

Li

,

F.

Schmidt-Kaler

, and

I.

Lesanovsky

,

Phys. Rev. Lett.

124

,

043402

(

2020

).

https://doi.org/10.1103/physrevlett.124.043402

Google Scholar

Crossref

PubMed

8.

H.

Tamura

,

T.

Yamakoshi

, and

K.

Nakagawa

,

Phys. Rev. A

101

,

043421

(

2020

).

https://doi.org/10.1103/physreva.101.043421

Google Scholar

Crossref

9.

R. V.

Skannrup

,

R.

Gerritsma

, and

S.

Kokkelmans

, arXiv:2008.13622 (

2020

).

10.

P.

Méhaignerie

,

C.

Sayrin

,

J.-M.

Raimond

,

M.

Brune

, and

G.

Roux

,

Phys. Rev. A

107

,

063106

(

2023

).

https://doi.org/10.1103/PhysRevA.107.063106

Google Scholar

Crossref

11.

D. F. V.

James

,

Appl. Phys. B

66

,

181

–

190

(

1998

).

https://doi.org/10.1007/s003400050373

Google Scholar

Crossref

12.

D.

Porras

and

J. I.

Cirac

,

Phys. Rev. Lett.

93

,

263602

(

2004

).

https://doi.org/10.1103/physrevlett.93.263602

Google Scholar

Crossref

PubMed

13.

C.

Schneider

,

D.

Porras

, and

T.

Schaetz

,

Rep. Prog. Phys.

75

,

024401

(

2012

).

https://doi.org/10.1088/0034-4885/75/2/024401

Google Scholar

Crossref

PubMed

14.

D.

Kienzler

,

H.-Y.

Lo

,

B.

Keitch

,

L.

De Clercq

,

F.

Leupold

,

F.

Lindenfelser

,

M.

Marinelli

,

V.

Negnevitsky

, and

J.

Home

,

Science

347

,

53

(

2015

).

https://doi.org/10.1126/science.1261033

Google Scholar

Crossref

PubMed

15.

H.-Y.

Lo

,

D.

Kienzler

,

L.

de Clercq

,

M.

Marinelli

,

V.

Negnevitsky

,

B. C.

Keitch

, and

J. P.

Home

,

Nature

521

,

336

(

2015

).

https://doi.org/10.1038/nature14458

Google Scholar

Crossref

PubMed

16.

D.

Kienzler

,

H.-Y.

Lo

,

V.

Negnevitsky

,

C.

Flühmann

,

M.

Marinelli

, and

J. P.

Home

,

Phys. Rev. Lett.

119

,

033602

(

2017

).

https://doi.org/10.1103/physrevlett.119.033602

Google Scholar

Crossref

PubMed

17.

J.

Pérez-Ríos

,

Mol. Phys.

119

,

e1881637

(

2021

).

https://doi.org/10.1080/00268976.2021.1881637

Google Scholar

Crossref

18.

R. S.

Lous

and

R.

Gerritsma

, in

Advances in Atomic, Molecular, and Optical Physics

,

Advances In Atomic, Molecular, and Optical Physics

, edited by

L. F.

DiMauro

,

H.

Perrin

and

S. F.

Yelin

(

Academic Press

,

2022

), Vol.

71

, pp.

65

–

133

.

Google Scholar

Crossref

19.

C. H.

Valahu

,

V. C.

Olaya-Agudelo

,

R. J.

MacDonell

,

T.

Navickas

,

A. D.

Rao

,

M. J.

Millican

,

J. B.

Pérez-Sánchez

,

J.

Yuen-Zhou

,

M. J.

Biercuk

,

C.

Hempel

,

T. R.

Tan

, and

I.

Kassal

,

Nat. Chem.

15

,

1503

(

2023

).

https://doi.org/10.1038/s41557-023-01300-3

Google Scholar

Crossref

PubMed

20.

F. J.

Garcia-Vidal

,

C.

Ciuti

, and

T. W.

Ebbesen

,

Science

373

,

eabd0336

(

2021

).

https://doi.org/10.1126/science.abd0336

Google Scholar

Crossref

PubMed

21.

P.

Nalbach

and

M.

Thorwart

,

Phys. Rev. B

81

,

054308

(

2010

).

https://doi.org/10.1103/physrevb.81.054308

Google Scholar

Crossref

22.

D.

Kast

and

J.

Ankerhold

,

Phys. Rev. Lett.

110

,

010402

(

2013

).

https://doi.org/10.1103/physrevlett.110.010402

Google Scholar

Crossref

PubMed

23.

P.

Nalbach

and

M.

Thorwart

,

Phys. Rev. B

87

,

014116

(

2013

).

https://doi.org/10.1103/physrevb.87.014116

Google Scholar

Crossref

24.

F.

Otterpohl

,

P.

Nalbach

, and

M.

Thorwart

,

Phys. Rev. Lett.

129

,

120406

(

2022

).

https://doi.org/10.1103/physrevlett.129.120406

Google Scholar

Crossref

PubMed

25.

D.

Lee

,

N.

Salwen

, and

D.

Lee

,

Phys. Lett. B

503

,

223

(

2001

).

https://doi.org/10.1016/s0370-2693(01)00197-6

Google Scholar

Crossref

26.

S.

Rychkov

and

L. G.

Vitale

,

Phys. Rev. D

91

,

085011

(

2015

).

https://doi.org/10.1103/physrevd.91.085011

Google Scholar

Crossref

27.

D.

Szász-Schagrin

and

G.

Takács

,

Phys. Rev. D

106

,

025008

(

2022

).

https://doi.org/10.1103/PhysRevD.106.025008

Google Scholar

Crossref

28.

T.

Rakovszky

,

M.

Mestyán

,

M.

Collura

,

M.

Kormos

, and

G.

Takács

,

Nucl. Phys. B

911

,

805

(

2016

).

https://doi.org/10.1016/j.nuclphysb.2016.08.024

Google Scholar

Crossref

29.

N.

Anand

,

A. L.

Fitzpatrick

,

E.

Katz

,

Z. U.

Khandker

,

M. T.

Walters

, and

Y.

Xin

, arXiv:2005.13544 (

2020

).

30.

H.

Chen

,

A. L.

Fitzpatrick

, and

D.

Karateev

,

J. High Energy Phys.

2022

,

146

.

https://doi.org/10.1007/JHEP02(2022)146

Crossref

31.

L. V.

Delacrétaz

,

A. L.

Fitzpatrick

,

E.

Katz

, and

M. T.

Walters

,

J. High Energy Phys.

2023

,

45

.

https://doi.org/10.1007/JHEP02(2023)045

Crossref

32.

J.

del Pino

,

F. A. Y. N.

Schröder

,

A. W.

Chin

,

J.

Feist

, and

F. J.

Garcia-Vidal

,

Phys. Rev. Lett.

121

,

227401

(

2018

).

https://doi.org/10.1103/physrevlett.121.227401

Google Scholar

Crossref

PubMed

33.

D.

Wellnitz

,

G.

Pupillo

, and

J.

Schachenmayer

,

Commun. Phys.

5

,

120

(

2022

).

https://doi.org/10.1038/s42005-022-00892-5

Google Scholar

Crossref

34.

M. L.

Wall

,

A.

Safavi-Naini

, and

A. M.

Rey

,

Phys. Rev. A

94

,

053637

(

2016

).

https://doi.org/10.1103/physreva.94.053637

Google Scholar

Crossref

35.

N.

Makri

and

D. E.

Makarov

,

J. Chem. Phys.

102

,

4611

(

1995

).

https://doi.org/10.1063/1.469509

Google Scholar

Crossref

36.

M.

Thorwart

,

P.

Reimann

,

P.

Jung

, and

R.

Fox

,

Chem. Phys.

235

,

61

(

1998

).

https://doi.org/10.1016/s0301-0104(98)00128-1

Google Scholar

Crossref

37.

M.

Thorwart

,

P.

Reimann

, and

P.

Hänggi

,

Phys. Rev. E

62

,

5808

(

2000

).

https://doi.org/10.1103/physreve.62.5808

Google Scholar

Crossref

38.

T. L.

Schmidt

,

P.

Werner

,

L.

Mühlbacher

, and

A.

Komnik

,

Phys. Rev. B

78

,

235110

(

2008

).

https://doi.org/10.1103/physrevb.78.235110

Google Scholar

Crossref

39.

U.

Schollwöck

,

Ann. Phys.

326

,

96

(

2011

).

https://doi.org/10.1016/j.aop.2010.09.012

Google Scholar

Crossref

40.

R.

Orús

,

Ann. Phys.

349

,

117

(

2014

).

https://doi.org/10.1016/j.aop.2014.06.013

Google Scholar

Crossref

41.

S.

Montangero

,

Introduction to Tensor Network Methods

(

Springer

,

2018

).

Google Scholar

Crossref

42.

S. R.

White

and

A. E.

Feiguin

,

Phys. Rev. Lett.

93

,

076401

(

2004

).

https://doi.org/10.1103/physrevlett.93.076401

Google Scholar

Crossref

PubMed

43.

P.

Schmitteckert

,

Phys. Rev. B

70

,

121302

(

2004

).

https://doi.org/10.1103/physrevb.70.121302

Google Scholar

Crossref

44.

M.

Nuss

,

M.

Ganahl

,

E.

Arrigoni

,

W.

von der Linden

, and

H. G.

Evertz

,

Phys. Rev. B

91

,

085127

(

2015

).

https://doi.org/10.1103/physrevb.91.085127

Google Scholar

Crossref

45.

B.

Dóra

,

M. A.

Werner

, and

C. P.

Moca

,

Phys. Rev. B

96

,

155116

(

2017

).

https://doi.org/10.1103/physrevb.96.155116

Google Scholar

Crossref

46.

M.

Zwolak

and

G.

Vidal

,

Phys. Rev. Lett.

93

,

207205

(

2004

).

https://doi.org/10.1103/physrevlett.93.207205

Google Scholar

Crossref

PubMed

47.

A. J.

Daley

,

Adv. Phys.

63

,

77

(

2014

).

https://doi.org/10.1080/00018732.2014.933502

Google Scholar

Crossref

48.

G.

Preisser

,

D.

Wellnitz

,

T.

Botzung

, and

J.

Schachenmayer

,

Phys. Rev. A

108

,

012616

(

2023

).

https://doi.org/10.1103/physreva.108.012616

Google Scholar

Crossref

49.

T.

Shi

,

E.

Demler

, and

J.

Ignacio Cirac

,

Ann. Phys.

390

,

245

(

2018

).

https://doi.org/10.1016/j.aop.2017.11.014

Google Scholar

Crossref

50.

L.

Hackl

,

T.

Guaita

,

T.

Shi

,

J.

Haegeman

,

E.

Demler

, and

I.

Cirac

,

SciPost Phys.

9

,

048

(

2020

).

https://doi.org/10.21468/scipostphys.9.4.048

Google Scholar

Crossref

51.

J.

Schachenmayer

,

A.

Pikovski

, and

A. M.

Rey

,

Phys. Rev. X

5

,

011022

(

2015

); arXiv:1408.4441.

https://doi.org/10.1103/physrevx.5.011022

52.

A.

Piñeiro Orioli

,

A.

Safavi-Naini

,

M. L.

Wall

, and

A. M.

Rey

,

Phys. Rev. A

96

,

033607

(

2017

).

https://doi.org/10.1103/physreva.96.033607

Google Scholar

Crossref

53.

B.

Zhu

,

A. M.

Rey

, and

J.

Schachenmayer

,

New J. Phys.

21

,

082001

(

2019

).

https://doi.org/10.1088/1367-2630/ab354d

Google Scholar

Crossref

54.

L.

Joubert-Doriol

and

A. F.

Izmaylov

,

J. Chem. Phys.

142

,

134107

(

2015

).

https://doi.org/10.1063/1.4916384

Google Scholar

Crossref

55.

D. S.

Schlegel

,

F.

Minganti

, and

V.

Savona

, arXiv:2306.13708 (

2023

).

56.

J.

Huber

,

P.

Kirton

, and

P.

Rabl

,

SciPost Phys.

10

,

045

(

2021

); arXiv:2011.10049.

https://doi.org/10.21468/scipostphys.10.2.045

57.

J.

Huber

,

A. M.

Rey

, and

P.

Rabl

,

Phys. Rev. A

105

,

013716

(

2022

); arXiv:2105.00004.

https://doi.org/10.1103/physreva.105.013716

58.

V. P.

Singh

and

H.

Weimer

,

Phys. Rev. Lett.

128

,

200602

(

2022

); arXiv:2108.07273.

https://doi.org/10.1103/physrevlett.128.200602

59.

C. D.

Mink

,

D.

Petrosyan

, and

M.

Fleischhauer

,

Phys. Rev. Res.

4

,

043136

(

2022

).

https://doi.org/10.1103/physrevresearch.4.043136

Google Scholar

Crossref

60.

C. D.

Mink

and

M.

Fleischhauer

,

SciPost Phys.

15

,

233

(

2023

).

https://doi.org/10.21468/SciPostPhys.15.6.233

Google Scholar

Crossref

61.

K.

Nagao

,

I.

Danshita

, and

S.

Yunoki

, “

Semiclassical descriptions of dissipative dynamics of strongly interacting Bose gases in optical lattices

,” arXiv:2307.16170 [cond-mat] (

2024

).

Google Scholar

62.

Y.

Ashida

,

T.

Shi

,

R.

Schmidt

,

H. R.

Sadeghpour

,

J. I.

Cirac

, and

E.

Demler

,

Phys. Rev. A

100

,

043618

(

2019

).

https://doi.org/10.1103/physreva.100.043618

Google Scholar

Crossref

63.

Y.

Ashida

,

T.

Shi

,

R.

Schmidt

,

H. R.

Sadeghpour

,

J. I.

Cirac

, and

E.

Demler

,

Phys. Rev. Lett.

123

,

183001

(

2019

).

https://doi.org/10.1103/physrevlett.123.183001

Google Scholar

Crossref

PubMed

64.

J.

Knörzer

,

T.

Shi

,

E.

Demler

, and

J. I.

Cirac

,

Phys. Rev. Lett.

128

,

120404

(

2022

).

https://doi.org/10.1103/physrevlett.128.120404

Google Scholar

Crossref

PubMed

65.

A.

Christianen

,

J. I.

Cirac

, and

R.

Schmidt

,

Phys. Rev. Lett.

128

,

183401

(

2022

).

https://doi.org/10.1103/physrevlett.128.183401

Google Scholar

Crossref

PubMed

66.

A.

Christianen

,

J. I.

Cirac

, and

R.

Schmidt

,

Phys. Rev. A

105

,

053302

(

2022

).

https://doi.org/10.1103/physreva.105.053302

Google Scholar

Crossref

67.

P. E.

Dolgirev

,

Y.-F.

Qu

,

M. B.

Zvonarev

,

T.

Shi

, and

E.

Demler

,

Phys. Rev. X

11

,

041015

(

2021

).

https://doi.org/10.1103/physrevx.11.041015

Google Scholar

Crossref

68.

L. J.

Bond

,

A.

Safavi-Naini

, and

J. c. v.

Minář

,

Phys. Rev. Lett.

132

,

170401

(

2024

).

https://doi.org/10.1103/physrevlett.132.170401

Google Scholar

Crossref

PubMed

69.

L.

Chen

,

Y.

Yan

,

M. F.

Gelin

, and

Z.

Lü

,

J. Chem. Phys.

158

,

104109

(

2023

).

https://doi.org/10.1063/5.0138399

Google Scholar

Crossref

70.

Y.

Zhao

,

J. Chem. Phys.

158

,

080901

(

2023

).

https://doi.org/10.1063/5.0140002

Google Scholar

Crossref

71.

K.

Sun

,

C.

Dou

,

M. F.

Gelin

, and

Y.

Zhao

,

J. Chem. Phys.

156

,

024102

(

2022

).

https://doi.org/10.1063/5.0076485

Google Scholar

Crossref

72.

N.

Zhou

,

L.

Chen

,

D.

Xu

,

V.

Chernyak

, and

Y.

Zhao

,

Phys. Rev. B

91

,

195129

(

2015

).

https://doi.org/10.1103/physrevb.91.195129

Google Scholar

Crossref

73.

N.

Zhou

,

L.

Chen

,

Y.

Zhao

,

D.

Mozyrsky

,

V.

Chernyak

, and

Y.

Zhao

,

Phys. Rev. B

90

,

155135

(

2014

).

https://doi.org/10.1103/physrevb.90.155135

Google Scholar

Crossref

74.

N.

Wu

,

L.

Duan

,

X.

Li

, and

Y.

Zhao

,

J. Chem. Phys.

138

,

084111

(

2013

).

https://doi.org/10.1063/1.4792502

Google Scholar

Crossref

75.

S.

Bera

,

A.

Nazir

,

A. W.

Chin

,

H. U.

Baranger

, and

S.

Florens

,

Phys. Rev. B

90

,

075110

(

2014

).

https://doi.org/10.1103/physrevb.90.075110

Google Scholar

Crossref

76.

S.

Lepoutre

,

J.

Schachenmayer

,

L.

Gabardos

,

B.

Zhu

,

B.

Naylor

,

E.

Maréchal

,

O.

Gorceix

,

A. M.

Rey

,

L.

Vernac

, and

B.

Laburthe-Tolra

,

Nat. Commun.

10

,

1714

(

2019

).

https://doi.org/10.1038/s41467-019-09699-5

Google Scholar

Crossref

PubMed

77.

M. A.

Perlin

,

C.

Qu

, and

A. M.

Rey

,

Phys. Rev. Lett.

125

,

223401

(

2020

); arXiv:2006.00723.

https://doi.org/10.1103/physrevlett.125.223401

78.

J.

Franke

,

S. R.

Muleady

,

R.

Kaubruegger

,

F.

Kranzl

,

R.

Blatt

,

A. M.

Rey

,

M. K.

Joshi

, and

C. F.

Roos

,

Nature

621

,

740

(

2023

); arXiv:2303.10688.

https://doi.org/10.1038/s41586-023-06472-z

79.

S. R.

Muleady

,

M.

Yang

,

S. R.

White

, and

A. M.

Rey

,

Phys. Rev. Lett.

131

,

150401

(

2023

); arXiv:2305.17242.

https://doi.org/10.1103/physrevlett.131.150401

80.

J. T.

Young

,

E.

Chaparro

,

A. P.

Orioli

,

J. K.

Thompson

, and

A. M.

Rey

, arXiv:2401.06222 (

2024

).

81.

A.

Polkovnikov

,

Ann. Phys.

325

,

1790

(

2010

); arXiv:0905.3384.

https://doi.org/10.1016/j.aop.2010.02.006

82.

The number of spin configurations in the present ansatz scales exponentially with the number of spins and is the bottleneck in the use of the ansatz in Eq. (1) for many spins. Addressing this issue, such as combining the non-Gaussian ansatz for the bosonic modes with tensor-network techniques for spins, is a matter of future work. We discuss an alternative approach, namely using the Holstein–Primakoff representation of the spins, in Sec. V C 2.

83.

A remark is that instead of the Dirac–Frenkel variational principle, one could employ the Lagrangian or McLachlan ones. Importantly, these coincide in that they yield the same equations of motion if the tangent space $T_{ψ}$ of the variational manifold is a Kähler space.

84.

Note a similar identity holds for open quantum states, obtained by using ρ and the trace.

85.

K. B.

Mo/ller

,

T. G.

Jo/rgensen

, and

J. P.

Dahl

,

Phys. Rev. A

54

,

5378

(

1996

).

https://doi.org/10.1103/physreva.54.5378

Google Scholar

Crossref

PubMed

86.

W.

Louisell

,

Quantum Statistical Properties of Radiation, A

(

Wiley-Interscience Publication

,

1973

).

Google Scholar

87.

G. J.

Milburn

,

Phys. Rev. A

33

,

674

(

1986

).

https://doi.org/10.1103/physreva.33.674

Google Scholar

Crossref

88.

F. R.

Halpern

,

J. Math. Phys.

14

,

219

(

1973

).

https://doi.org/10.1063/1.1666299

Google Scholar

Crossref

89.

C. M.

Bender

and

T. T.

Wu

,

Phys. Rev. D

7

,

1620

(

1973

).

https://doi.org/10.1103/physrevd.7.1620

Google Scholar

Crossref

90.

C. M.

Bender

and

T. T.

Wu

,

Phys. Rev.

184

,

1231

(

1969

).

https://doi.org/10.1103/physrev.184.1231

Google Scholar

Crossref

91.

H.-P.

Breuer

and

F.

Petruccione

,

The Theory of Open Quantum Systems

(

Oxford University Press

,

2002

).

Google Scholar

92.

U.

Weiss

,

Quantum Dissipative Systems

(

World Scientific

,

2012

).

Google Scholar

Crossref

93.

C.

Gardiner

and

P.

Zoller

,

Quantum Noise

(

Springer

,

2004

).

Google Scholar

94.

J.

Dalibard

,

Y.

Castin

, and

K.

Mølmer

,

Phys. Rev. Lett.

68

,

580

(

1992

).

https://doi.org/10.1103/physrevlett.68.580

Google Scholar

Crossref

PubMed

95.

K.

Mølmer

,

Y.

Castin

, and

J.

Dalibard

,

J. Opt. Soc. Am. B

10

,

524

(

1993

).

https://doi.org/10.1364/josab.10.000524

Google Scholar

Crossref

96.

M. B.

Plenio

and

P. L.

Knight

,

Rev. Mod. Phys.

70

,

101

(

1998

).

https://doi.org/10.1103/revmodphys.70.101

Google Scholar

Crossref

97.

X.

Yuan

,

S.

Endo

,

Q.

Zhao

,

Y.

Li

, and

S.

Benjamin

,

Quantum

3

,

191

(

2019

); arxiv:1812.08767 [quant-ph].

https://doi.org/10.22331/q-2019-10-07-191

98.

B.

Gerritsen

,

L.

Bond

,

J.

Minar

,

R.

Spreeuw

,

R.

Gerritsma

, and

A.

Safavi-Naini

, “

Sub-Doppler cooling of trapped ions using optical tweezers

” (unpublished).

99.

W. K.

Wootters

,

Ann. Phys.

176

,

1

(

1987

).

https://doi.org/10.1016/0003-4916(87)90176-x

Google Scholar

Crossref

100.

J.

Schachenmayer

,

A.

Pikovski

, and

A. M.

Rey

,

New J. Phys.

17

,

065009

(

2015

); arXiv:1501.06593.

https://doi.org/10.1088/1367-2630/17/6/065009

101.

O. L.

Acevedo

,

A.

Safavi-Naini

,

J.

Schachenmayer

,

M. L.

Wall

,

R.

Nandkishore

, and

A. M.

Rey

,

Phys. Rev. A

96

,

033604

(

2017

).

https://doi.org/10.1103/physreva.96.033604

Google Scholar

Crossref

102.

M.

Kunimi

,

K.

Nagao

,

S.

Goto

, and

I.

Danshita

,

Phys. Rev. Res.

3

,

013060

(

2021

).

https://doi.org/10.1103/physrevresearch.3.013060

Google Scholar

Crossref

103.

A. P.

Orioli

,

A.

Signoles

,

H.

Wildhagen

,

G.

Günter

,

J.

Berges

,

S.

Whitlock

, and

M.

Weidemüller

,

Phys. Rev. Lett.

120

,

063601

(

2018

).

https://doi.org/10.1103/physrevlett.120.063601

Google Scholar

Crossref

PubMed

104.

Y. A.

Alaoui

,

S. R.

Muleady

,

E.

Chaparro

,

Y.

Trifa

,

A. M.

Rey

,

T.

Roscilde

,

B.

Laburthe-Tolra

, and

L.

Vernac

, arXiv:2404.10531,

1

(

2024

).

105.

F.

Herrera

and

F. C.

Spano

,

Phys. Rev. Lett.

116

,

238301

(

2016

).

https://doi.org/10.1103/physrevlett.116.238301

Google Scholar

Crossref

PubMed

106.

M.

Vogl

,

P.

Laurell

,

H.

Zhang

,

S.

Okamoto

, and

G. A.

Fiete

,

Phys. Rev. Res.

2

,

043243

(

2020

).

https://doi.org/10.1103/physrevresearch.2.043243

Google Scholar

Crossref

107.

J.

Marshall

and

N.

Anand

,

Opt. Quantum

1

,

78

(

2023

).

https://doi.org/10.1364/opticaq.504311

Google Scholar

Crossref

108.

In the numerical implementation of NGS, we observe that for certain initial states the equations of motion become numerically unstable, which we attribute to the singularity of the metric g and the associated need to perform a pseudoinverse instead of the inverse in the evaluation of G = g⁻¹. Dealing with this issue is a matter of current and future work.

109.

J.

Wurtz

,

A.

Polkovnikov

, and

D.

Sels

,

Ann. Phys.

395

,

341

(

2018

).

https://doi.org/10.1016/j.aop.2018.06.001

Google Scholar

Crossref

110.

T.

Vovk

and

H.

Pichler

,

Phys. Rev. Lett.

128

,

243601

(

2022

).

https://doi.org/10.1103/physrevlett.128.243601

Google Scholar

Crossref

PubMed

111.

T.

Vovk

and

H.

Pichler

, “

Quantum trajectory entanglement in various unravelings of Markovian dynamics

,”

Phys. Rev. A

110

,

012207

(

2024

); arXiv:2404.12167 [quant-ph].

https://doi.org/10.1103/physreva.110.012207

Google Scholar

112.

M. P.

Kaicher

,

D.

Vodola

, and

S. B.

Jäger

,

Phys. Rev. B

107

,

165144

(

2023

).

https://doi.org/10.1103/physrevb.107.165144

Google Scholar

Crossref

113.

E.

Gottlob

and

U.

Schneider

,

Phys. Rev. B

107

,

144202

(

2023

).

https://doi.org/10.1103/physrevb.107.144202

Google Scholar

Crossref

114.

J. A.

Hutchison

,

T.

Schwartz

,

C.

Genet

,

E.

Devaux

, and

T. W.

Ebbesen

,

Angew. Chem., Int. Ed.

51

,

1592

(

2012

).

https://doi.org/10.1002/anie.201107033

Google Scholar

Crossref

115.

B. S.

Simpkins

,

A. D.

Dunkelberger

, and

J. C.

Owrutsky

,

J. Phys. Chem. C

125

,

19081

(

2021

).

https://doi.org/10.1021/acs.jpcc.1c05362

Google Scholar

Crossref

116.

D. S.

Wang

and

S. F.

Yelin

,

ACS Photonics

8

,

2818

(

2021

).

https://doi.org/10.1021/acsphotonics.1c01028

Google Scholar

Crossref

117.

J. A.

Campos-Gonzalez-Angulo

,

Y. R.

Poh

,

M.

Du

, and

J.

Yuen-Zhou

,

J. Chem. Phys.

158

,

230901

(

2023

).

https://doi.org/10.1063/5.0143253

Google Scholar

Crossref

© 2024 Author(s). All article content, except where otherwise noted, is licensed under a Creative Commons Attribution-NonCommercial-NoDerivs 4.0 International (CC BY-NC-ND) license (https://creativecommons.org/licenses/by-nc-nd/4.0/).

Open quantum dynamics with variational non-Gaussian states and the truncated Wigner approximation

I. INTRODUCTION

II. NON-GAUSSIAN ANSATZ FOR A CLOSED SYSTEM

A. Analytic expressions for energies and energy gradients

B. Tangent vectors and the overlap matrix

III. OPEN DYNAMICS: COMBINING NGS WITH QUANTUM TRAJECTORIES

A. Quantum trajectories

B. Equations of motion for non-Hermitian Hamiltonians

C. Quantum jumps

1. Jumps inside the manifold

2. Jumps outside the manifold

D. Time evolution

IV. TRUNCATED WIGNER APPROXIMATION FOR SPINS AND BOSONS

V. RESULTS

A. Details of NGS and TWA simulations

B. Closed system dynamics

C. Open system dynamics

1. Holstein Tavis–Cummings

2. Tavis–Cummings

VI. CONCLUSIONS AND OUTLOOK

SUPPLEMENTARY MATERIAL

ACKNOWLEDGMENTS

AUTHOR DECLARATIONS

Conflict of Interest

Author Contributions

DATA AVAILABILITY

REFERENCES

Citing articles via

Submit your article

Sign up for alerts

Resources

Explore

pubs.aip.org

Connect with AIP Publishing

Open quantum dynamics with variational non-Gaussian states and the truncated Wigner approximation

I. INTRODUCTION

II. NON-GAUSSIAN ANSATZ FOR A CLOSED SYSTEM

A. Analytic expressions for energies and energy gradients

B. Tangent vectors and the overlap matrix

III. OPEN DYNAMICS: COMBINING NGS WITH QUANTUM TRAJECTORIES

A. Quantum trajectories

B. Equations of motion for non-Hermitian Hamiltonians

C. Quantum jumps

1. Jumps inside the manifold

2. Jumps outside the manifold

D. Time evolution

IV. TRUNCATED WIGNER APPROXIMATION FOR SPINS AND BOSONS

V. RESULTS

A. Details of NGS and TWA simulations

B. Closed system dynamics

C. Open system dynamics

1. Holstein Tavis–Cummings

2. Tavis–Cummings

VI. CONCLUSIONS AND OUTLOOK

SUPPLEMENTARY MATERIAL

ACKNOWLEDGMENTS

AUTHOR DECLARATIONS

Conflict of Interest

Author Contributions

DATA AVAILABILITY

REFERENCES

Citing articles via

Submit your article

Sign up for alerts

Related Content

Resources

Explore

pubs.aip.org

Connect with AIP Publishing

This Feature Is Available To Subscribers Only