State space distribution and dynamical flow for closed and open quantum systems

One of the foundational concepts in chemical physics is that the observable properties of macroscopic systems can be represented as an average over an ensemble of identical but statistically independent microscopic subsystems. The field of statistical mechanics provides the theoretical formalism for characterizing these subsystem statistics and relating them to macroscopic observables. In this formalism, ensemble statistics and dynamics are conveniently expressed in terms of time evolving probability distributions over the subsystem state space. In classical mechanics, these probability distributions are formulated in Hamiltonian phase space and evolve according to Liouville’s equation. Generalizing this formulation to quantum mechanics has been a longstanding problem due to the wave nature of quantum states and the influence of measurements on system dynamics, which does not allow for a well-defined probability distribution in phase space. In this letter, we show that this problem can be solved by expressing quantum ensembles as probability distributions over the state space of density operators, thereby formally representing the lack of knowledge of the initial preparation of a system and its bath.

The state of a quantum ensemble is typically described by a single density operator, $ρ^$⁠. This representation efficiently encodes the statistics of ensemble observables and is central to quantum theories for dynamics,^1–4 optics,^5,6 thermodynamics,^7–9 information,^10,11 and control.^12–14 However, a single density operator is an incomplete description of an ensemble because it contains no information about the state of individual systems. Rather, $ρ^$ describes only the uncertainty of the ensemble observables, combining classical contributions, due to an unknown initial state of the individual system or bath, and quantum contributions, due to the random outcome of measurements on superposition states. As a result, different ensembles can be represented by the same $ρ^$⁠, complicating their microscopic interpretation. Classical uncertainty in the initial preparation of the bath of each system can lead to heterogeneity in the dynamics of systems within the ensemble. The importance of this physical insight has been exemplified in the study of light-induced biomolecular dynamics,^15–20 environment-conditioned qubit dynamics,²¹ and emerging ultrafast single-molecule spectroscopies.^22–27 

A general treatment of quantum dynamics as stochastic processes on Liouville space has also been proposed by Davies.²⁸ This approach has been effective at describing the effects of measurement on the dynamics of quantum systems and for the formulation of stochastic methods for numerical quantum dynamics. The quantum filtering theory of Belavkin and co-workers describes the gain of information in continuously monitored quantum systems and the resulting change in quantum state.²⁹ The measurement at each instant probabilistically collapses the wavefunction into an a posteriori state conditioned on the measurement outcome, producing a stochastic process that randomly drives system dynamics. Averaging over all possible measurement outcomes then reproduces the typical ensemble density matrix, $ρ^$⁠. The stochastic wavefunction unraveling methods first proposed by Mølmer, Castin, and Dalibard³⁰ and Dum, Zoller, and Ritsch³¹ exploit this observation to numerically propagate a density matrix as the average of a stochastic process on the system Hilbert space. This has led to the quantum jump³² and quantum state diffusion³³ unraveling methods which differ in the nature of the stochastic process used to generate trajectory ensembles. While these stochastic trajectories were first used as numerical intermediates, Wiseman and Milburn³⁴ showed that they can be physically interpreted as the dynamics of a system conditioned on the outcome of a continuous strong measurement of its bath. Each unraveling then corresponds to a specific bath measurement, with quantum jump trajectories being associated with continuous photon counting measurements and quantum state diffusion with heterodyne measurements.

The source of stochasticity in both quantum filtering theory and wavefunction unraveling methods is the probabilistic collapse of the wavefunction conditioned on a measurement outcome. This differs significantly from the type of uncertainty present in classical systems where measurement has no impact on the state of the system. Stochasticity in classical system arises due to the lack of knowledge of the initial state of the system or its surrounding bath. As such, the underlying stochastic process in classical systems is a model for approximating the deterministic dynamics of a very high dimensional environment with an unspecified initial state. The physical origin of classical stochasticity is statistical dispersion in the initial state of the system and bath. The same classical uncertainty is also present in quantum systems but cannot be directly observed in ensemble experiments where its effect is combined with the quantum uncertainty in measurement outcome. However, some forms of this classical uncertainty may be measurable in emerging single-molecule ultrafast spectroscopy experiments. In these experiments, individual molecules are impulsively excited and allowed to propagate for a short waiting time t ∼ 10 fs–1 ps after which they are measured with a probe pulse. By repeating the experiment several times on the same molecule, the statistics of the observable can then be measured on a system by system basis. Moreover, the time taken for each measurement repetition may be much faster than some of the bath modes, in principle allowing multiple measurements to be conducted with fixed values of these slow environment coordinates. This procedure then measures the short time dynamics of a system conditioned on the slow environment modes, resolving the classical uncertainty due to these sluggish coordinates.

As a concrete example, consider an ensemble of excitonic systems composed of a dye heterodimer embedded in a DNA origami scaffold. These systems allow for strict control of the position and orientation of dye molecules and have therefore been proposed as a platform for programmable excitonic coupling.³⁵ In a single molecule experiment, an ensemble of individual DNA-dye systems may be attached to the surface of the slide with random orientations, as shown schematically in Fig. 1(a). Now, consider three different types of ultrafast experiments that can be done on this randomly oriented ensemble: (i) a single molecule experiment where only one DNA-dye system is in the optical spot of the experiment, (ii) a spatially resolved experiment with a small subensemble of assemblies in the optical spot, and (iii) an ensemble experiment where the whole ensemble is sampled. Each of these experimental measurements samples different ensembles of states. If the randomly oriented systems are excited by a polarized light pulse, each excitation will be initialized in a different superposition of donor and acceptor dyes. Single molecule measurements on any one system will then observe a pure state superposition of donor and acceptor dyes, producing a distribution of initial pure state excitations sketched in red in Fig. 1(c). A spatially resolved experiment that measures a few assemblies will then observe a mixed state that averages the initial excitation of the single molecules in the optical spot, as shown in blue in Fig. 1(c), while an ensemble experiment will average over all possible orientations, giving the state plotted in black. The dynamics measured with these three types of experiments will differ significantly. Variations between single molecule measurements contain information that can be used to resolve classical disorder in the ensemble of initial states, while the ensemble measurement averages over this disorder and combines it with the quantum measurement uncertainty.

In addition to this random orientation of excitonic circuits, individual DNA origami scaffolds can deform, modifying the relative orientation of donor and acceptor dyes shown in Fig. 1(b). This bath induced orientation change modifies the excitonic coupling between the dyes, thereby changing their dynamics. However, these fluctuations are often several orders of magnitude slower than the time scale of experimental observation, allowing multiple measurements to be performed on a single orientation. A sequence of single molecule measurements taken in a short time window will thus sample a single configuration of these slow bath modes, resolving the effect of disorder in those coordinates. In addition to this experimentally relevant situation, we can consider the dynamics conditioned on an initial configuration of all bath modes, analogously to an explicit solvent classical molecular dynamics trajectory. The resulting trajectory ensemble encodes the effect of bath disorder on the system dynamics in a similar manner to classical statistical ensembles, providing a potentially new theoretical bridge between quantum and classical ensembles.

Here, we present a theoretical formalism for treating quantum ensembles analogously to classical ensembles. This approach resolves the challenges that arise due to the wave properties of quantum states on phase space by working in a natural quantum state space. We outline this formalism in Sec. II by first providing a brief summary of classical ensembles on phase space in Sec. II A. Then, by considering stochastic processes on Liouville space, as first proposed by Davies,²⁸ we define a quantum state space probability density in Sec. II B. We then derive equations of motion for these distributions in Sec. II C. This formalism is then compared to previously studied wavefunction unraveling approaches in Sec. II D. Remarkably, these take the same form as the classical Liouville’s equation. This equivalence systematically generalizes the methods and intuition of classical statistical mechanics to quantum systems, which we illustrate by proving a novel quantum Liouville theorem in Sec. III. We provide an example applying our approach to the DNA-dye system of Fig. 1 in Sec. IV. Finally, we conclude by summarizing the results of this paper in Sec. V.

We begin by briefly outlining the classical theory of phase space ensembles. The state of a classical system is described by enumerating the positions and momenta of all particles. For a system with N position coordinates ${qi}i=1N$ and momenta ${pi}i=1N$⁠, this defines a 2N-dimensional vector x = [q, p] in phase space, $R2N$⁠. The evolution of a closed system is governed by Hamilton’s equations of motion and defines a vector field $ẋ=[∂H/∂q,−∂H/∂p]$ called the dynamical flow field. This can be modified to treat dissipative systems by adding damping terms. For example, linear dissipators, given by positive matrices ${Γ^α}$⁠, yield a dissipative flow field, $ẋ=ẋC−∑αΓ^αx$⁠, where $ẋC$ is the closed system flow field. A system initially in x₀ propagates along the flow field to x(t|x₀) at time t, tracing out a trajectory {x(t|x₀)|t ≥ 0}. Example flow fields are shown for closed and damped harmonic motion in Figs. 2(a) and 2(c).

An ensemble is comprised of a collection of systems, each found in a different state. The state of an ensemble is then given by a probability distribution P(x; t), where x is a coordinate identifying a point in phase space. Liouville’s equation describes the dynamics of this distribution in terms of the flow field, giving

where $κ(x)≔−div(ẋ)$ is the compressibility of the flow field and $a⋅b≔∑i=12Naibi$ is the inner product on $R2N$⁠. If the single system dynamics, x(t|x_o), are known, Eq. (1) can be analytically solved to give P(x; t) = ∫dx′P(x′; 0)δ(x(t|x′) − x), where the convection Green’s function δ(x(t|x′) − x) selects initial conditions x′ that evolve to x at time t. Distribution and trajectory dynamics are illustrated for two harmonic oscillator in Figs. 2(b) and 2(d), respectively.

The isosurfaces in Fig. 2 reveal an important property of classical systems. Closed system distributions remain the same size at all times, while damped system distributions compress toward a low energy steady state. Liouville’s theorem relates this notion of phase space volume to the reversibility and determinism of the underlying dynamics. In a closed system, the phase space volume, or equivalently the local density evaluated along an evolving trajectory, P(x(t|x₀); t), is constant and the flow field is incompressible everywhere [κ(x) = 0], reflecting the fact that reversible, deterministic trajectories cannot cross. This property does not hold for open systems where linear dissipative flow has positive compressibility. The distribution collapses onto a steady state point x_SS where the dynamics are irreversible since the initial state of a trajectory that has relaxed to x_SS cannot be determined.

This classical statistical structure is effective because it is constructed on a space where each point uniquely identifies a state. Motivated by this observation, we formulate quantum ensembles on a quantum state space rather than on classical phase space. The state of an open quantum system is defined by a density operator $σ^$ in Liouville space, $L(H)$⁠, a complex vector space of linear operators acting on the system Hilbert space, $H$⁠. Given a basis ${|i}i=1N$ for $H$⁠, any operator can be written in vector form $σ^→σ≔[σ11,σ12,…,σNN]$ by listing its matrix elements $σij≔i|σ^|j$⁠. Moreover, it is equipped with a trace inner product, $α⋅β≔Tr{α^†β^}=∑ijαij*βij$⁠.

Similarly to classical ensembles, quantum ensembles are comprised of a collection of systems, each in a different state. The state of the ensemble can then be defined by a probability distribution P(σ; t) on Liouville space, where σ is a coordinate that identifies a point in Liouville space. This formulation provides a more detailed account of the state of a quantum ensemble than the ensemble density matrix $ρ^$⁠. While $ρ^$ describes only the statistics of the ensemble observables, P(σ; t) encodes the state of subsystems in the ensemble. In fact, $ρ^$ is the average of this distribution, ρ(t) ≔ ∫dσP(σ; t) = ⟨σ⟩(t), which integrates over the classical uncertainty. As a result, P(σ; t) resolves the ambiguity between quantum and classical uncertainty by treating each quantum state as a well-defined point. The classical uncertainty in the preparation of system states is described by the probability distribution, P(σ; t). The quantum uncertainty in observable outcomes of each subsystem is separately contained in the quantum description of its state, σ.

We now consider the dynamics of an ensemble of N dimensional quantum systems. First, we define a quantum dynamical flow field on Liouville space. This describes the evolution of a quantum system at point σ. For a closed system with Hamiltonian H, this is given by the Liouville von-Neumann equation,

where $L0=−iℏ[H,⋅]$ is the Liouville superoperator.

This can be extended to open systems by adding dissipator terms. For simplicity, we restrict our attention in the main text to Markovian dynamics, providing the general non-Markovian Nakajima-Zwanzig theory^2,36,37 in  Appendix B. The general form of Markovian dynamics is given by the Gorini-Kossakowski-Sudarshan-Lindblad (GKSL) equation,^1,2,38–40

where the Lindblad operators ${Lα}α=1N2$ form an orthonormal basis of $L(H)$⁠, $D[Lα]σ=LασLα†−12{Lα†Lα,σ}$ is the associated dissipators with rates γ_α ≥ 0 and {A, B} = AB + BA is the anticommutator. By convention, $LN2=I$ is the identity operator and ${Lα}α=1N2−1$ have vanishing trace. This expression takes a form similar to classical linear dissipation described in Sec. II A.

A system initialized in state σ₀ evolves for a time t by moving along the flow field. This yields a state σ(t|σ₀) and traces out a trajectory {σ(t|σ₀)|t ≥ 0} through Liouville space that is continuous and defined at all times. Since trajectories are defined at all times, dynamics never creates or destroys trajectories. This property is known as trajectory conservation.

Using these properties, we derive a continuity equation for trajectories in Liouville space. First, we consider the probability P_Ω(t) of finding a system in a region Ω at time t. This quantity can change in one of the 3 ways: (1) A subsystem state can enter or leave Ω, that is, there exists a time t′ where some subsystem trajectory is outside (inside) of Ω immediately before t′ and inside (outside) of it immediately after. (2) A subsystem is created in Ω, that is, at a time t′, a new trajectory is created inside Ω that was not defined for 0 ≤ t < t′. (3) A subsystem is destroyed in Ω, that is, at a time t′, a trajectory that was inside Ω is no longer defined for t > t′.

The evolution of a subsystem state is defined for all time t ≥ 0, since a system can be propagated by specifying an initial state σ₀. This property of the dynamics prevents the creation or destruction of trajectories required for processes (2) and (3) above to occur since they require a trajectory to be undefined for some times t ≥ 0. Trajectories generated by these dynamics are said to be conserved, and all changes in P_Ω must therefore occur due to a subsystem entering or leaving Ω through process (1). Moreover, since all trajectories produced by these dynamics are continuous, a trajectory cannot enter or leave Ω without passing through its boundary.

Consequently, the change in P_Ω must be related to the probability flux, $j(σ;t)≔P(σ;t)σ̇(t)$⁠, at its boundary, giving the integral form of the continuity equation,

where the second integral is taken over the boundary, ∂Ω, of Ω and dA(σ) denotes the infinitesimal surface normal at σ in the surface-flux integral. The gradient of a function f(σ) on Liouville space is given by $∇$f(σ) ≔ [∂f(σ)/∂σ₁₁, …, ∂f(σ)/∂σ_NN], and the divergence of a vector field of g(σ) is defined as $div(g(σ))≔∑ij∂gij(σ)/∂σij*$⁠.⁴¹ 

Equation (4) can be written in differential form by applying Gauss’s divergence theorem to the surface-flux integral. This gives

which applies to any region Ω. By taking Ω to be an infinitesimal volume element surrounding a point σ,⁴² we then obtain the differential form

where we have applied the divergence product rule to j(σ; t) and $κ(σ)≔−div(σ̇)$ is the flow field compressibility. Conveniently, this derivation mirrors the classical derivation of Eq. (1) exactly, highlighting the formal simplicity of the state space distribution formalism.

Equation (6) can then be applied to closed and open dynamics by substituting the appropriate flow field [Eq. (2) or (3)] to give

which directly propagate state-space distributions.

Remarkably, Eq. (6) describing quantum ensemble dynamics takes an identical form to Eq. (1) for classical ensembles, thus revealing that uncertainty in the initial preparation of systems in an ensemble propagates similarly for quantum and classical systems. Consequently, the problematic behavior of quantum mechanics on phase space arises due to the incompatibility of quantum uncertainty that prohibits the simultaneous knowledge of positions and momenta, with a position momentum state space. This difficulty can be avoided by treating quantum dynamics on its natural state space where the statistical structure of a dynamical system is conserved.

The similarity between Eqs. (1) and (6) makes it straightforward to propagate the dynamics of quantum ensembles. Methods used to solve Eq. (1) can be directly applied to Eq. (6). If the dynamics of the microscopic system, σ(t|σ₀), are known, the distribution dynamics can be obtained using the Green’s function method, yielding

where the Green’s function δ(σ(t|σ′) − σ) selects initial conditions σ′ that evolve to state σ at time t.

While the mathematical objects used to represent these ensembles are analogous to those used in stochastic unraveling and quantum filtering theories, they represent very different forms of uncertainty. This leads to several key physical differences between the two theories. First, the trajectories generated in the two methods have very different physical meaning. The trajectories generated by stochastic unravelings represent the time-dependent state of a system conditioned on a sequence of observations on its surrounding bath with different measurement schemes. Typically, these theories assume a known, potentially mixed, initial state of the system and bath. As such, the stochasticity in this trajectory ensemble originates from the quantum uncertainty of observables of superposition states. This quantum uncertainty manifests as a classical-like stochastic process, with a corresponding probability distribution, through probabilistic measurement-induced wavefunction collapse. Consequently, this form of trajectory ensemble has no classical counterpart, as classical measurement does not influence the system state.

The trajectories underlying the distributions discussed in this manuscript do not assume any measurement. Instead, the heterogeneity in system dynamics originates purely from an unknown initial preparation of the system and bath. These different initial conditions will then evolve differently producing a heterogeneous ensemble of dynamical trajectories. This form of uncertainty is analogous to uncertainty in classical statistical mechanics that physically originates from an unknown initial state of the system and bath.

As a result, the trajectory ensemble and resulting time-dependent probability distribution model different types of experiments. The stochastic unraveling trajectory ensembles describe experiments where a single system and its environment are continuously measured. In this case, the stochastic process is physically realized through the measurement process. These experiments have been realized in cavity and trapped ion systems where complete observation is possible but are not possible in more complicated condensed phase systems. In contrast, the trajectories and distributions discussed in this paper describe experiments where only a subsystem of interest is experimentally accessible and measured systems are discarded after measurement. This is the case, for example, in stroboscopic single molecule experiments where the system is reinitialized in the same initial state after each measurement. The reinitialization step in these experiments removes the effect of measurement on the trajectory ensemble as no system is measured twice.

Despite these differences between wavefunction unraveling approaches and the approach we describe here, these two approaches are complimentary and can be formally related to each other. In particular, they can be combined by performing a stochastic unraveling to each of the initial condition trajectories considered in this study, or equivalently by stochastically sampling the initial condition of a stochastic unraveling calculation from the initial state distribution. The resulting initial state conditioned unraveled trajectory ensemble will then show the effect of both initial state dispersion and continuous measurement, showing both quantum and classical forms of uncertainty. In this light, the relationship between the two types of stochastic processes is clear. Traditional stochastic unraveling processes average over initial state disorder but resolve the effect of quantum uncertainty rendered classical through continuous measurement. As such, the different types of unraveling approaches reflect the information gain from different measurement schemes. In contrast, the ensembles considered in this study resolve the initial state disorder but average over the measurement outcomes. These averaging processes are not necessarily undesirable as they mathematically remove information that is not conditioned on in a given experiment.

The identical structure of quantum and classical ensembles can also be exploited to systematically generalize key classical results built on Liouville’s equation. To provide a road map for extending classical results, we show that the classical Liouville’s theorem can be trivially extended to quantum ensembles. Notably, previous phase space studies^43–45 have concluded that quantum dynamical flow in phase space cannot be written in an incompressible form. This shows that quantum dynamics on phase space differs fundamentally from classical mechanics. In contrast, the applicability of Liouville’s theorem shows that classical notions of reversibility and determinism can still be applied to quantum ensembles on $L(H)$⁠.

Following the classical derivation, we consider the probability density evaluated along a Liouville space trajectory P(σ(t|σ₀); t). The rate of change of this quantity is

after an application of the chain rule and a substitution of Eq. (6). By substituting the Liouville-von Neumann flow field for a closed system [Eq. (2)], we find that $div(σ̇)=0$ for all σ since div(Hσ) = div(⁠$σH$⁠). This indicates that the classical picture of deterministic, reversible dynamics in terms of nonintersecting state space trajectories applies directly to quantum trajectories. Moreover, the derivation of this result directly mimics the classical proof, highlighting the ease of applying classical derivations in this formalism. The same process can be repeated for GKSL dynamics. Substituting Eq. (3) into Eq. (9) gives a uniform compressibility $κ(σ)=N∑α=1N2−1γα≥0$⁠, taking a similar form to classical linear dissipation. A glossary summarizing the quantum-classical analogies in our formalism is provided in  Appendix A.

To demonstrate the application of the formulation presented above, we consider the dynamics of excitations on the ensemble of dye dimer systems shown in Fig. 1. We first present a two-level Frenkel model for the exciton dynamics and a GKSL model for the effect of a Markovian bath of intramolecular nuclear vibrations. Using this model, we then treat two types of classical uncertainties and provide examples of two numerical approaches to propagating the resulting distributions. First, we consider the effect of stochastic orientation of the DNA origami scaffold which leads to uncertainty in the initial preparation of the exciton. This is solved by analytically performing the Green’s function integral in Eq. (8) using the analytical solution for single system dynamics. The resulting distribution dynamics show how a specific experiment can be encoded in the distribution and are consistent with the quantum Liouville theorem proved in Sec. III. We then consider the effect of slowly varying deformations of the DNA origami scaffolding on exciton dynamics by generating a trajectory ensemble that performs a Monte-Carlo integration of Eq. (8). This example shows the manifestation of classical uncertainty in single-molecule experiments and demonstrates the resulting dynamical heterogeneity.

Excitonic systems can be conveniently modeled using a finite set of discrete quantum states that represent an excitation localized on a given dye. For the dimeric systems considered in this section, this yields a two state basis {|D⟩, |A⟩} where each dye has excitation frequencies ω_α and transition dipole moments μ_α, where α = A or α = D denote the acceptor or donor molecule, respectively. Excitations on the two dyes are coupled, with an interaction given by

in the dipole approximation, where ϵ₀ is the permittivity of free space and r is the displacement vector between the two dyes. This coupling leads to delocalized eigenstates, {|1⟩, |2⟩}, separated by a frequency $ω=2ω02+4VDA2/ℏ2$⁠, where ω₀ ≔ |ω_D − ω_A| is the difference in excitation frequency between the dyes.

Letting {|α⟩, |β⟩} be any orthonormal basis for the two level systems, density matrices can be conveniently expressed in the Pauli basis ${Si}i=0,x,y,z$ in $L(H)$⁠, where $S0≔I/2$ is the normalized identity operator, and $Sx≔(|αβ|+|βα|)/2$⁠, $Sy≔−i(|αβ|−|βα|)/2$⁠, and $Sz≔(|ββ|−|αα|)/2$⁠. Any density matrix $σ=(S0+xSx+ySy+zSz)/2$ can be written as a real valued 3D vector σ → [x, y, z] with $x2+y2+z2≤1$⁠, where a basis change can be enacted by a rotation of the 3D vectors.

The effect of a Markovian environment on the exciton dynamics can be modeled using the GKSL formulation. The closed system evolves under the Hamiltonian $H=ℏω0Sz/2+2VDASx$ in the site basis. For simplicity, we restrict our attention to dye homodimers with ω₀ = 0 and μ_D = μ_A, both real. In this case, the excitonic states are simply given by |1⟩ ∝ |D⟩ + |A⟩ and |2⟩ ∝ |D⟩ − |A⟩. The effect of the bath can then be conveniently modeled by two Lindblad operators.

The first, L₁ = |1⟩⟨2|, induces |2⟩ → |1⟩ transitions that represent the transfer of electronic energy to the phonon bath of nuclear vibrations leading to dissipation with Lindblad rate Γ, while L₂ = S_z fluctuates the energy of the two sites leading to dephasing with Lindblad rate γ_ϕ. In the site basis, Eq. (3) gives a flow field,

plotted in Fig. 3(c). The unitary evolution gives a flow field, $σ̇U=[0,−ωz,ωy]$ [Fig. 3(a)], which precesses the state of the system about the x axis. This represents the coherent oscillation of the excitation between the |D⟩ and |A⟩ sites. Dissipative flow, $σ̇Diss=[Γ(1−x),−Γ/2y,−Γ/2z]$⁠, pushes the population into the lower energy excitonic state located at the point (1, 0, 0) on the Bloch sphere. Dephasing flow, $σ̇Dep=[0,−γϕy,−γx]$⁠, does not drive any transitions and so acts perpendicular to x. However, the fluctuating energy gap destroys phase information between excitonic states, driving pure states on the surface of the Bloch sphere to mixed states on the x axis. When combined, this leads to an excitation that oscillates between the localized sites before slowly relaxing to the lower energy excitonic eigenstate.

Substituting Eq. (11) into Eq. (9), closed systems with γ_ϕ = 0 = Γ have a vanishing compressibility, validating the quantum Liouville theorem. For open systems, this yields the expected uniform compressibility κ(σ) = 2(γ_ϕ + Γ), indicating that GKSL dynamics extend linear dissipation to quantum systems. In particular, the purely quantum dephasing process describing the loss of quantum phase information is formulated equivalently to classical friction.

The dynamics of this Markovian model are analytically known for all initial conditions σ₀ ≔ [x₀, y₀, z₀]. This gives

The closed system trajectory [Fig. 3(b)] simply precesses about the x axis, while the open system [Fig. 3(d)] spirals toward the |1⟩ state.

Combining Eqs. (8) and (11), an initial distribution P(σ; 0) can be directly propagated. We first consider the case where this initial distribution represents uncertainty in the initial excitation generated by a linearly polarized pump pulse in a single molecule experiment. As shown in Fig. 1, this arises when the excitonic circuits attach with an unknown orientation and when more than one circuit are found in the excitation spot. The stochastic orientation of the circuits leads to the preparation of random pure states leading to dispersion on the surface of the Bloch sphere. In addition, multiple circuits with different orientations may be found within the excitation spot. An experiment conducted on this subensemble of molecules will then observe the average excitation, producing a mixed state. This partial averaging procedure leads to a distribution that contains both mixed and pure states, thereby encoding the information that is available in a given experiment. The initial distribution can be directly calculated from the orientation statistics of the ensemble and the statistics that govern the number of circuits in the excitation spot.

For simplicity, consider a truncated Gaussian distribution, $P(σ;0)∝N(σ;σ¯0,Δ^0)Θ(1−|σ|)$ with initial mean $σ¯0$ and covariance matrix $Δ^0$ and Θ is the Heaviside function. The resulting time dependent distribution,

remains a truncated Gaussian at all times. The time-dependent mean $σ¯(t|σ¯0)$ is obtained by propagating the initial mean using Eq. (11). The time dependent covariance matrix $Δt=S(ν(t))Rx(ωt)Δ0Rx(ωt)TS(ν(t))$⁠, where $Rx(θ)$ is the x rotation by angle θ and $S(ν(t))$ is the scaling matrix that scales the Cartesian axes by ν_z = exp(−(γ_ϕ + Γ/2)t) = ν_y and ν_x = exp(−Γt).

As an example, we consider an experiment that attempts to selectively excite site |D⟩ (e.g., by using an excitation pulse orthogonal to the expected direction of μ_A) in Fig. 3. This yields a Gaussian that rotates around the x axis while collapsing toward |1⟩ as shown in Fig. 3(d). The closed system distribution, shown in Fig. 3(b), simply precesses about the x axis. We provide the solution for a general initial distribution in  Appendix C.

In addition to this unknown initial excitation of the system, the phonon bath of scaffolding vibrations can modify the relative orientation of the dyes on a time scale much slower than the experimental measurement time. This leads to an uncertainty in the initial state of the bath that can be experimentally conditioned on by performing a series of measurements in quick succession before the bath degrees of freedom change appreciably. Since varying the relative orientation of the dyes can substantially modify the interdye coupling, different systems in a single molecule experiment can show very different dynamics. The resulting dynamical heterogeneity is analogous to classical dynamical heterogeneity that originates from an unknown initial state of the bath. In the formalism presented in this paper, these two forms of stochasticity are encoded in the same way, enabling them to be treated analogously.

To demonstrate this, we treat the slow modes in a manner similar to classical molecular dynamics simulations with explicit solvent. That is, we first randomly sample the initial state of the slow bath modes and then propagate the dynamics conditioned on this initial state. By histogramming these trajectories, we obtain a time dependent probability density through a Monte-Carlo integration of the Green’s function integral Eq. (8). In this example, the back action of the exciton on the slow modes can be neglected as these modes typically involve large scale motions of the large DNA assembly. As a result, the slow bath modes are well approximated by a constant initial value and do not need to be propagated, significantly reducing the required computational effort. Moreover, the fast intramolecular modes that cannot be conditioned on in this proposed experiment are averaged over and treated using the GKSL approach, allowing for a selective sampling that reflects the experimental details. The resulting algorithm is analogous to the frozen modes’ approach that has recently been introduced to facilitate the treatment of non-Markovian quantum dynamics.^46,47

In Fig. 4, we present the results for a coupled heterodimer with varying initial state of the slow bath mode. We assume that the transition dipole moment of both dyes is orthogonal to their displacement, leading to a vanishing second term in Eq. (10). Moreover, we consider bath deformations that rotate the relative orientation of the dyes in this plane, modifying the angle θ between them. This leads to the following simplified form of the bath dependent dipolar coupling:

In the resulting heterodimer, the ratio μ_Dμ_A/ℏω₀r³ then quantifies the influence of the coupling on the exciton energetics and therefore the magnitude of the effect of variations in θ. To highlight the effect of bath heterogeneity, we consider systems that are consistently prepared on the localized state |D⟩ in order to remove the effect of uncertainty in system preparation. We take the ratio μ_Dμ_A/ℏω₀r³ = 1 and take a normal distribution of angles θ with a mean of 0 and standard deviation of π/3. Dissipation is neglected in this case and dephasing due to energy fluctuations occurs with a rate γ = ω/4.

A set of sample trajectories that would be observed in the single molecule experiment described above are shown in Fig. 4(a) alongside the ensemble averaged trajectory shown in red. This trajectory ensemble shows that systems prepared in the same system initial state |D⟩ will diverge when projected onto the Bloch sphere. These spreading trajectories demonstrate how, even in this modest model for bath heterogeneity, uncertainty in the initial preparation of the bath can substantially change system dynamics. Moreover, we see that as the trajectories dephase, they approach the z axis, moving toward an incoherent mixture. Notably, the ensemble averaged trajectory dephases more rapidly than the individual systems as the variations in coupling lead to systems that oscillate at different frequencies, destroying the well-defined phase information in the ensemble. Consequently, the distributions shown in Fig. 4(b) first spread from their well-defined initial state before beginning to contract along the z axis as they approach their incoherent steady state. The heterogeneity in the trajectory ensemble and the resulting spread of the probability distributions encode information about the bath in the statistics of single molecule trajectories, indicating that single molecule experiments may be able to probe bath heterogeneity.

In conclusion, the Liouville space of density operators provides a natural state space for the study and interpretation of quantum ensembles. Individual quantum states are represented by discrete points in this state space, and their dynamics obey flow properties identical to those of classical systems in phase space. By exploiting the familiar classical structure of this quantum state-space, it is possible to directly apply tools from classical statistical mechanics to quantum systems and to exercise classical intuition when interpreting their statistical properties. The formalism we have presented here thus casts the challenging problem of quantum mechanical mixed states in a form that is mathematically similar to classical ensembles, potentially enabling a unified treatment of quantum and classical statistical mechanics.

The authors thank Aurelia Chenu and David Coker for useful discussions. This work was supported by the National Science Foundation under award number CHE-1839155. A.D. acknowledges support from the Natural Science and Engineering Research Council of Canada Postgraduate Scholarship Doctoral (NSERC PGS D) program, and A.W. acknowledges support from the Research Corporation for Scientific Advancement (Cottrell Scholar).

The Liouville state space formalism presented in the main text is constructed to mirror the structure of classical statistical mechanics on Hamiltonian phase space. As such, nearly all objects in the classical theory have a quantum analog in the Liouville space theory. Below, we provide a glossary that summarizes the key analogs between the two theories.

To derive an equation of motion for open quantum systems, we follow the Nakajima-Zwanzig approach.² In this approach, the Liouville-von Neumann dynamics of a composite system and bath are projected (via the use of projection superoperators) onto two different Hilbert spaces called the relevant and irrelevant Hilbert spaces. We define the relevant projection superoperator, $P$⁠, as

where Tr_B indicates a trace over the bath and $ρB∈L(HB)$ is a stationary bath reference state that is normalized so that $TrBρB=1$⁠. Defined in this way, $Pσ$ yields a projected density operator $σrel∈LHS⊗HB$ that is related to the more familiar reduced system density operator via $σS=TrBσrel$⁠. The irrelevant projection superoperator is given by

where $I$ is the identity superoperator.

With these projection superoperators, the dynamics of the irrelevant space can be formally solved and expressed in terms of its effect on the dynamics within the relevant Hilbert space. Specifically, consider a general system-bath Hamiltonian,

where H_S is the system Hamiltonian, H_B is the bath Hamiltonian, and V describes the interaction between the system and bath. The dynamics of such a system can be expressed in the interaction picture using the Nakajima-Zwanzig equation,

where

and

where $L(t)σ≔V(t),σ(t)$ is the Liouville superoperator, $T←$ is the time-ordering superoperator, and t₀ corresponds to the time at which the system is initialized.

The Nakajima-Zwanzig equation can be simplified with the appropriate choice of initial conditions. For example, the first term in Eq. (B4a) vanishes in the case where the bath reference state is selected so that $TrBV(t)ρB=0$⁠.² Likewise, the second term (describing contributions arising due to entangled initial conditions) vanishes when σ(t₀) = σ_S(t₀) ⊗ ρ_B. For simplicity, we will restrict out attention to systems whose initial conditions cause the first two terms in Eq. (B4a) to vanish in this way. This leaves only a homogeneous integrodifferential equation, as described by the final term in Eq. (B4a).

Dynamical flow under the Nakajima-Zwanzig equation can be determined from the non-Markovian flow field, $σ̇S=TrB{σ̇rel}$⁠. To compute the divergence of this flow field, we take the component-wise functional derivative of $σ̇S$ with respect to σ_S to get

where the bath is referred to with Greek indices and the system is referred to with Latin indices.

The divergence can be computed from this expression to yield

which represents the compressibility of Nakajima-Zwanzig flow on the state space of reduced density operators. This expression corresponds to general non-Markovian dynamics, including the Markovian limit where the memory kernel $K(t,t′)∝δ(t−t′)$⁠. Notably, a similar functional derivative approach is also used to derive the Euler-La Grange equation. Finally, combining Eq. (B6) and the fluid mechanical equation of motion in the main text leads to an equation of motion for open quantum systems of the form

which is similar to that of the generalized Langevin equation.⁴⁸ This formalism can thus be used to develop a quantum analog to the generalized non-Markovian Fokker Planck equation.

We now consider the distribution dynamics of a spin system coupled to a Bosonic bath. For convenience, we will perform the derivation in the energy eigenbasis, noting that the results can be converted to the site basis of a homodimer by a rotation of π/2 about the y axis. It is convenient to express the dynamics of the system in the Pauli basis ${Si}i∈{0,x,y,z}$⁠, where $S0≔I/2$ is the normalized identity matrix and $Sx≔(|12+|21)/2$⁠, $Sy≔−i(|12|−|21|)/2$⁠, and $Sz≔(|22|−|11|)/2$ are the normalized Pauli matrices. The $1/2$ normalization constant comes about due to normalization with respect to the trace inner product. In this notation, the system Hamiltonian is given by

and the density matrix can be conveniently written in Bloch vector notation as $σ≔(xSx+ySy+zSz+I)/2→[x,y,z]$⁠. The bath interactions are treated phenomenologically through the following Lindblad dissipators:

where L₁ ≔ |1⟩⟨2| is the dissipation operator that induces transitions from the excited state to the ground state and L₂ ≔ σ_z is the pure dephasing operator that fluctuates the spin energy splitting. The Lindblad rates for these processes are given by the dissipation rate Γ and the dephasing rate γ_ϕ, respectively.

Using the GKSL equation [Eq. (3) of the main text], the dynamical flow field for this model system can then be written as

The subsystem dynamics generated by Eq. (C3) are analytically solvable. The dynamics are given by

where the initial condition is written in Bloch vector notation as σ₀ ≔ [x₀, y₀, z₀].

These subsystem results can be generalized to obtain the dynamics of an arbitrary initial distribution P(σ; 0) using convective Green’s function equation given in the main text, Eq. (6), as

where the convective Green’s function δ(σ(t|$σ′$⁠) − σ) selects initial conditions $σ′$ = [x′, y′, z′] that evolve to σ at time t. Using Eq. (C4), the convective Green’s function can be simplified to give

which expresses the Green’s function as explicit functions of the integrating variable $σ′$⁠.

Combining Eqs. (C5) and (C6), the time-dependent distribution can be found by transforming the coordinates of the initial distribution P(σ; 0) giving

where $S(ν)$ is the scaling matrix and $Rz(θ)$ is the matrix for a rotation by angle θ about the z axis of the Bloch sphere. The normalization factor in Eq. (C7a) accounts for the change in normalization due to the scaling transformation and arises from the determinant of the scaling transformation. These transformations are explicitly written in matrix form as

In the main text, we consider a Gaussian initial distribution with mean $σ¯0$ and covariance matrix $Δ^0$ given by

where $A$ denotes the determinant of A and A⁻¹, its inverse. Applying Eq. (C7) to a Gaussian initial distribution yields a Gaussian at all times with time varying mean and covariance

where $σ¯(t|σ¯0)$ is obtained by propagating the initial mean using Eq. (C4) and the time dependent covariance is given by the transformation $Δ^t≔S(ν(t))Rz(ωt)Δ^0Rz(ωt)TS(ν(t))$ where the scaling vector is given by $ν(t)≔[e−(γϕ+12Γ)t,e−(γϕ+12Γ)t,e−Γt]$⁠. Notably, the scaling factor in Eq. (C7a) accounts for the renormalization of the Gaussian due to the time varying covariance matrix.