Observables¶

This section describes the theoretical framework for computing observables in cavity-coupled molecular dynamics, including ensemble averages, time-correlation functions, and material time for aging experiments.

Overview ¶

In molecular dynamics simulations, we compute ensemble averages and time-correlation functions to characterize system behavior. For non-equilibrium cavity-coupled systems, special care is needed to properly average over stochastic trajectories and handle time-dependent distributions.

Key Concepts:

Phase space: \(\Gamma = (\mathbf{R}, \mathbf{P}, \mathbf{q}, \mathbf{p})\) for positions, momenta, cavity coordinates, and cavity momenta
Distribution: \(\rho^{(\lambda)}(\Gamma, t)\) evolves according to Fokker-Planck equation
Ensemble average: Integral over phase space weighted by \(\rho\)
Time-correlation: Relates observable at two different times

Ensemble Averages ¶

Schrodinger vs Heisenberg Picture ¶

Schrodinger Picture:

Time dependence is in the phase-space distribution. For an observable \(A(\Gamma)\), the ensemble average at time \(t\) is:

\[\langle A(t) \rangle_\lambda = \int d\Gamma\,A(\Gamma)\,\rho^{(\lambda)}(\Gamma, t)\]

Heisenberg Picture:

Time dependence is in the observable. Given initial condition \(\Gamma\), evolve the observable averaged over all realizations of noise:

\[A_t(\Gamma; \lambda) = \mathbb{E}\left[ A(\Gamma(t; \lambda)) \mid \Gamma(0) = \Gamma \right]\]

The ensemble average becomes:

\[\langle A(t) \rangle_\lambda = \int d\Gamma\,A_t(\Gamma; \lambda)\,\rho_\mathrm{eq}^{(\lambda=0)}(\Gamma)\]

Both pictures are equivalent but offer different computational advantages.

Fokker-Planck Equation ¶

The phase-space distribution \(\rho^{(\lambda)}(\Gamma, t)\) evolves according to:

\[\frac{\partial \rho^{(\lambda)}}{\partial t} = \mathcal{L}_\lambda \rho^{(\lambda)}\]

where the Liouvillian operator is:

\[\mathcal{L}_\lambda = \mathcal{L}_0 + \mathcal{L}_\mathrm{int}(\lambda) + \mathcal{L}_\mathrm{b}\]

Components:

Bare Liouvillian \(\mathcal{L}_0\): Hamiltonian dynamics at zero coupling

\[\mathcal{L}_0 \rho = -\sum_i \frac{\mathbf{P}_i}{M_i} \cdot \nabla_{\mathbf{R}_i} \rho + \sum_i \nabla_{\mathbf{R}_i} V \cdot \nabla_{\mathbf{P}_i} \rho - \sum_\alpha p_\alpha \frac{\partial \rho}{\partial q_\alpha} + \sum_\alpha \omega_\mathrm{c}^2 q_\alpha \frac{\partial \rho}{\partial p_\alpha}\]

Interaction Liouvillian \(\mathcal{L}_\mathrm{int}(\lambda)\): Light-matter coupling

\[\mathcal{L}_\mathrm{int}(\lambda)\rho = \sum_i \sum_\alpha \lambda z_i e \left(\omega_\mathrm{c} q_\alpha + \lambda \mu_\alpha\right) \hat{\mathbf{e}}_\alpha \cdot \nabla_{\mathbf{P}_i} \rho - \sum_\alpha \lambda \omega_\mathrm{c} \mu_\alpha \frac{\partial \rho}{\partial p_\alpha}\]

Bath Liouvillian \(\mathcal{L}_\mathrm{b}\): Noise and dissipation from molecular and cavity baths

Initial Condition:

For nonthermal aging protocols starting at \(t = t_0\):

\[\rho^{(\lambda)}(\Gamma, t_0) = \rho_\mathrm{eq}^{(\lambda=0)}(\Gamma)\]

where the equilibrium distribution at zero coupling is:

\[\rho_\mathrm{eq}^{(\lambda)}(\Gamma) = \frac{1}{Z(T, \lambda)} e^{-\beta H(\Gamma; \lambda)}\]

Practical Computation ¶

Sample Mean Over Trajectories:

In simulations, we approximate ensemble averages by averaging over \(N_\mathrm{T}\) independent trajectories:

\[\langle A(t) \rangle_\lambda \approx \frac{1}{N_\mathrm{T}} \sum_{n=1}^{N_\mathrm{T}} A^{(n)}(t; \lambda)\]

where \(A^{(n)}(t; \lambda) = A(\Gamma^{(n)}(t; \lambda))\) is the observable evaluated on the \(n\)-th trajectory.

Equilibrium Averages:

At equilibrium, we can use both ensemble and time averaging (ergodic hypothesis):

\[\langle A \rangle_\lambda = \lim_{T \to \infty} \frac{1}{T} \int_0^T dt\,A(\Gamma(t; \lambda))\]

In practice:

\[\langle A \rangle_\lambda \approx \frac{1}{MN_\mathrm{T}} \sum_{n=1}^{N_\mathrm{T}} \sum_{p=0}^{M-1} A(\Gamma^{(n)}(t_p; \lambda))\]

where we subsample at times \(t_p = p\Delta t\).

Static Structure Factor Example ¶

For density-mode fluctuations at wavevector \(\mathbf{k}\):

\[\rho_{\mathbf{k}}^{(n)}(t; \lambda) = \sum_{i=1}^N e^{i\mathbf{k} \cdot \mathbf{R}_i^{(n)}(t; \lambda)}\]

The static structure factor is:

\[S_{\mathbf{k}}^{(\lambda)}(t_\mathrm{w}) = \left\langle \rho_{\mathbf{k}}(t_\mathrm{w}) \rho_{\mathbf{k}}^*(t_\mathrm{w}) \right\rangle_\lambda \approx \frac{1}{N_\mathrm{T}} \sum_{n=1}^{N_\mathrm{T}} \rho_{\mathbf{k}}^{(n)}(t_\mathrm{w}; \lambda) \rho_{\mathbf{k}}^{*(n)}(t_\mathrm{w}; \lambda)\]

At equilibrium, this becomes time-independent:

\[S_{\mathbf{k}}^{(\lambda)} \approx \frac{1}{N_\mathrm{T}M} \sum_{n=1}^{N_\mathrm{T}} \sum_{p=0}^{M-1} \rho_{\mathbf{k}}^{(n)}(t_p; \lambda) \rho_{\mathbf{k}}^{*(n)}(t_p; \lambda)\]

Time-Correlation Functions ¶

General Non-Equilibrium Case ¶

For an observable \(A(t)\), the time-correlation function in a non-equilibrium aging experiment is:

\[C_{AA}^{(\lambda)}(t; t_\mathrm{w}) = \langle A(t + t_\mathrm{w}) A(t_\mathrm{w}) \rangle_\lambda\]

where \(t_\mathrm{w}\) is the waiting time after initiating the protocol and \(t\) is the correlation time.

Interpretation:

Depends on both \(t\) (time difference) and \(t_\mathrm{w}\) (absolute time)
Reflects non-equilibrium, aging dynamics
Converges to equilibrium correlation as \(t_\mathrm{w} \to \infty\)

Ensemble Average Form:

In terms of the conditionally averaged observable and phase-space distribution:

\[C_{AA}^{(\lambda)}(t; t_\mathrm{w}) = \int d\Gamma\,A_{t+t_\mathrm{w}}(\Gamma; \lambda)\,A(\Gamma)\,\rho^{(\lambda)}(\Gamma, t_\mathrm{w})\]

Trajectory Sampling:

We approximate by generating \(N_\mathrm{T}\) trajectories from equilibrium at zero coupling, simulating directly to \(t + t_\mathrm{w}\):

\[C_{AA}^{(\lambda)}(t; t_\mathrm{w}) \approx \frac{1}{N_\mathrm{T}} \sum_{n=1}^{N_\mathrm{T}} A^{(n)}(t + t_\mathrm{w}; \lambda)\,A^{(n)}(t_\mathrm{w}; \lambda)\]

Intermediate Scattering Function ¶

The time-correlation function of density modes is the intermediate scattering function (ISF):

\[F_{\mathbf{k}}^{(\lambda)}(t; t_\mathrm{w}) = \left\langle \rho_{\mathbf{k}}(t + t_\mathrm{w}; \lambda) \rho_{\mathbf{k}}^*(t_\mathrm{w}; \lambda) \right\rangle_\lambda\]

Practical Implementation:

Average over \(N_{\mathbf{k}}\) wavevectors of equal magnitude distributed on a sphere:

\[F_{k}^{(\lambda)}(t; t_\mathrm{w}) \approx \frac{1}{N_\mathrm{T}N_{\mathbf{k}}} \sum_{n=1}^{N_\mathrm{T}} \sum_{s=1}^{N_{\mathbf{k}}} \rho_{\mathbf{k}_s}^{(n)}(t + t_\mathrm{w}; \lambda)\,\rho_{\mathbf{k}_s}^{*(n)}(t_\mathrm{w}; \lambda)\]

Normalized ISF:

\[\phi_k(t; t_\mathrm{w}) = \frac{F_k(t; t_\mathrm{w})}{S_k(t_\mathrm{w})}\]

Structural Relaxation Time:

Define \(\tau_\mathrm{s}(t_\mathrm{w}, \lambda)\) as the time when ISF decays to threshold \(a\) (typically 0.1):

\[\phi_k^{(\lambda)}(\tau_\mathrm{s}; t_\mathrm{w}) = a\]

Located by linear interpolation between adjacent data points.

Equilibrium Case ¶

In equilibrium (large \(t_\mathrm{w}\)), time-correlation functions obey time-translational invariance:

\[C_{AA}^{(\lambda)}(t; t_\mathrm{w}) = C_{AA,\mathrm{eq}}^{(\lambda)}(t)\]

depending only on time difference \(t\).

Proof:

At equilibrium, \(\rho^{(\lambda)}(\Gamma, t_\mathrm{w}) = \rho_\mathrm{eq}^{(\lambda)}(\Gamma)\) which satisfies \(\mathcal{L}_\lambda \rho_\mathrm{eq}^{(\lambda)} = 0\). Since the dynamics are time-translationally invariant, \(A_{t+t_\mathrm{w}}(\Gamma) = A_t(\Gamma)\) with initial condition at \(t=0\).

Ergodic Hypothesis:

At equilibrium, ensemble average equals long-time average:

\[C_{AA,\mathrm{eq}}^{(\lambda)}(t) = \lim_{T \to \infty} \frac{1}{T} \int_0^T dt'\,A(\Gamma(t + t'; \lambda))\,A(\Gamma(t'; \lambda))\]

Discrete Estimator:

Sample at times \(t_p = p\Delta t\), compute autocorrelation at lag \(\tau_k = k\Delta t\):

\[\bar{C}_{AA}^{(\lambda)}(\tau_k) \approx \frac{1}{N_\mathrm{T}(M-k)} \sum_{n=1}^{N_\mathrm{T}} \sum_{p=0}^{M-k-1} A^{(n)}(t_p; \lambda)\,A^{(n)}(t_{p+k}; \lambda)\]

This is the standard “sliding time-origin” implementation.

Dipole Autocorrelation and IR Spectrum ¶

Centered DACF:

For dipole moment, subtract mean to remove bias:

\[\delta\mu_\alpha(t) = \mu_\alpha(t) - \bar{\mu}_\alpha\]

\[\bar{C}_{\mu\mu}^{(\lambda)}(\tau_k) \approx \frac{1}{N_\mathrm{T}(M-k)} \sum_{n=1}^{N_\mathrm{T}} \sum_{p=0}^{M-k-1} \sum_{\alpha=1,2} \delta\mu_\alpha(t_p; \lambda)\,\delta\mu_\alpha(t_{p+k}; \lambda)\]

IR Absorption Spectrum:

From the fluctuation-dissipation theorem:

\[n(\omega)\alpha(\omega) = \frac{\beta\omega^2}{4\epsilon_0 Vc} \int_{-\infty}^{\infty} dt\,e^{i\omega t}\,\bar{C}_{\mu\mu}^{(\lambda)}(t)\]

Using time-reversibility, \(\bar{C}_{\mu\mu}(t) = \bar{C}_{\mu\mu}(-t)\):

\[n(\omega)\alpha(\omega) = \frac{\beta\omega^2}{4\epsilon_0 Vc} \cdot 2\int_0^{\tau_{\max}} dt\,\cos(\omega t)\,\bar{C}_{\mu\mu}^{(\lambda)}(t)\]

Discrete Cosine Transform:

Sample \(C_n = C_\mu(n\Delta t)\) uniformly with \(N\) points. Define frequencies:

\[\omega_k = \frac{\pi k}{(N-1)\Delta t}\]

Compute DCT:

\[X_k = (-1)^k C_{N-1} + 2\sum_{n=1}^{N-2} C_n \cos\left(\frac{\pi nk}{N-1}\right)\]

Trapezoidal rule yields:

\[n(\omega_k)\alpha(\omega_k) \approx \frac{\beta\omega_k^2}{4\epsilon_0 Vc} \cdot \frac{\Delta t}{2} X_k\]

Material Time ¶

Material-Time Translational Invariance ¶

In aging systems, time-correlation functions across different protocols can be reparameterized using material time \(h_\lambda(t)\):

MTTI Hypothesis:

All correlation functions collapse to a universal curve when plotted against material-time differences:

\[C_{AA}^{(\lambda)}(t; t_\mathrm{w}) = \Phi\left(h_\lambda(t + t_\mathrm{w}) - h_\lambda(t_\mathrm{w})\right)\]

Physical Interpretation:

Material time measures the “internal age” of the system. Equal material-time differences correspond to equal structural evolution, regardless of absolute time or coupling strength.

Reconstruction from Relaxation Times ¶

Structural Relaxation Constraint:

For each waiting time \(t_\mathrm{w}^{(r)}\), measure \(\tau_\mathrm{s}^{(r)}\) where \(\phi_k(\tau_\mathrm{s}^{(r)}; t_\mathrm{w}^{(r)}) = 0.1\). MTTI implies:

\[h_\lambda(t_\mathrm{w}^{(r)} + \tau_\mathrm{s}^{(r)}) - h_\lambda(t_\mathrm{w}^{(r)}) = 1, \quad r = 1, \ldots, N_\mathrm{w}\]

Least-Squares Formulation:

Discretize time on a grid \(t_m = m\Delta t\) for \(m = 0, \ldots, M_\mathrm{grid}-1\). Parameterize:

\[\tilde{h}_\lambda(t) = \sum_{m=0}^{M_\mathrm{grid}-1} h_\lambda^m \theta_m(t)\]

where \(\theta_m(t)\) are piecewise linear basis functions:

\[\begin{split}\theta_m(t) = \begin{cases} \frac{t - t_{m-1}}{t_m - t_{m-1}}, & t \in [t_{m-1}, t_m] \\ \frac{t_{m+1} - t}{t_{m+1} - t_m}, & t \in [t_m, t_{m+1}] \\ 0, & \text{otherwise} \end{cases}\end{split}\]

Constraint Matrix:

Define \(N_\mathrm{w} \times M_\mathrm{grid}\) matrix:

\[A_{rm} = \theta_m(t_\mathrm{w}^{(r)} + \tau_\mathrm{s}^{(r)}) - \theta_m(t_\mathrm{w}^{(r)})\]

so the constraints become linear system \(\mathbf{A}\mathbf{h}_\lambda = \mathbf{b}\) where \(\mathbf{b} = [1, \ldots, 1]^\mathsf{T}\).

Regularized Fitting ¶

Smoothness Penalty:

To avoid overfitting, add curvature penalty. Define second-derivative matrix \(\mathbf{D}\) (finite differences):

\[\begin{split}D_{rm} = \begin{cases} 1, & m = r \\ -2, & m = r+1 \\ 1, & m = r+2 \\ 0, & \text{otherwise} \end{cases}\end{split}\]

Regularized Objective:

\[\mathbf{h}_\lambda^* = \arg\min_{\mathbf{h}_\lambda} \left\{ \|\mathbf{A}\mathbf{h}_\lambda - \mathbf{b}\|^2 + \alpha\|\mathbf{D}\mathbf{h}_\lambda\|^2 \right\}\]

where \(\alpha > 0\) controls smoothness-constraint tradeoff.

Solution:

Taking gradient and setting to zero:

\[(\mathbf{A}^\mathsf{T}\mathbf{A} + \alpha\mathbf{D}^\mathsf{T}\mathbf{D})\mathbf{h}_\lambda^* = \mathbf{A}^\mathsf{T}\mathbf{b}\]

Solve using sparse linear solvers (conjugate gradient or Cholesky decomposition).