A single wavepacket is useful for visualizing dynamics. Thermal rate and correlation calculations need an ensemble. This note reorganizes the density-matrix workflow: thermal initialization, momentum filtering, multi-state propagation, and absorber-based population accounting.

Hilbert space and density matrix

The full space is a product of the nuclear DVR space and a finite electronic space:

\[ \mathcal{H}=\mathcal{H}_{\mathrm{nuc}}\otimes\mathcal{H}_{\mathrm{el}}. \]

The propagated object is the density matrix \(\rho(t)\), not a single wavefunction. This naturally supports thermal initial states and mixed incoming ensembles.

DVR and FFT operators

In the ensemble source note, the coordinate grid is uniform and the momentum operator is built through a discrete Fourier transform:

\[ F_{jk}=\frac{1}{\sqrt{N_x}}\exp\!\left(-\frac{2\pi i}{N_x}jk\right), \qquad p_\alpha=\hbar\,2\pi\,\mathrm{fftfreq}(N_x,d=\Delta x). \]
\[ P=F^\dagger\mathrm{diag}(p_\alpha)F,\qquad P^2=F^\dagger\mathrm{diag}(p_\alpha^2)F. \]

The FFT form is convenient for constructing thermal states and momentum filters on a uniform grid.

Thermal initialization

A harmonic reference Hamiltonian can define the incoming nuclear thermal state:

\[ H_{\mathrm{ho}}=\frac{P^2}{2m} +\frac{1}{2}m\omega_{\mathrm{init}}^2(X-x_{\mathrm{init}})^2, \qquad \rho_{\mathrm{nuc}}=\frac{e^{-\beta H_{\mathrm{ho}}}}{\mathrm{Tr}(e^{-\beta H_{\mathrm{ho}}})}. \]

Diagonalizing \(H_{\mathrm{ho}}\) gives a stable way to evaluate the matrix exponential.

Positive-momentum filtering

To represent an incoming ensemble from left to right, the thermal density can be filtered in momentum space. The sharp version is

\[ f(p)= \begin{cases} 1,& p\gt 0,\\ 1/2,& p=0,\\ 0,& p\lt 0, \end{cases} \qquad P_f=F^\dagger\mathrm{diag}(f(p_\alpha))F. \]

A smoother alternative is \(f(p)=\frac{1}{2}[1+\tanh(p/dp)]\). The filtered density is

\[ \rho_{\mathrm{nuc}}^{(+)} = \frac{P_f\rho_{\mathrm{nuc}}P_f} {\mathrm{Tr}(P_f\rho_{\mathrm{nuc}}P_f)}. \]

Multi-state lift

The diabatic Hamiltonian is

\[ \hat H= \frac{\hat P^2}{2m}\otimes I + \sum_{a,b}|a\rangle V_{ab}(\hat X)\langle b|. \]

On the product basis this is again

\[ H_{(i,a),(j,b)}=\delta_{ab}T_{ij}+\delta_{ij}V_{ab}(x_i). \]

The source note also tracks local adiabatic populations by diagonalizing \(V(x_i)\) at each grid point and using the resulting local eigenvectors as projectors.

Time propagation

For a Hermitian Hamiltonian, the density matrix evolves as

\[ \rho(t+\Delta t)= e^{-iH\Delta t/\hbar}\rho(t)e^{+iH\Delta t/\hbar}, \]

which is the finite-dimensional form of the Liouville-von Neumann equation \(i\hbar\,d\rho/dt=[H,\rho]\). If no absorber is applied, this can also be evaluated directly at any time from \(H=C\,\mathrm{diag}(E_\alpha)C^\dagger\): the energy-basis density matrix gains phase factors \(e^{-i(E_\alpha-E_\beta)t/\hbar}\). The step notation appears here because the workflow interleaves unitary propagation with counting and mask operations.

For the published SAC benchmark, the filtered density is propagated in the equivalent low-rank form

\[ \rho^{(+)}\simeq \sum_{\alpha=1}^{r} w_\alpha |\psi_\alpha\rangle\langle\psi_\alpha|, \qquad r=20, \]

so each \(|\psi_\alpha\rangle\) is propagated and masked while observables are accumulated with weight \(w_\alpha\). This avoids a large dense density-matrix multiply at every step while preserving the density-matrix ensemble interpretation for the linear population diagnostics used here.

Absorber bookkeeping

Instead of treating the final population as only what remains inside the numerical box, the source workflow records a ledger. A left keep-mask can be written as

\[ M_L(x)= \begin{cases} 0,& x\le x_1,\\ \sin^2\!\left[\frac{\pi}{2}\frac{x-x_1}{x_2-x_1}\right],& x_1\lt x\lt x_2,\\ 1,& x\ge x_2, \end{cases} \]

and a right keep-mask as

\[ M_R(x)= \begin{cases} 1,& x\le x_1,\\ \cos^2\!\left[\frac{\pi}{2}\frac{x-x_1}{x_2-x_1}\right],& x_1\lt x\lt x_2,\\ 0,& x\ge x_2. \end{cases} \]

The useful reported quantity is therefore

\[ P_k^{\mathrm{total}}(t)= P_k^{\mathrm{remain}}(t) +P_{k,\mathrm{left\ out}}(t) +P_{k,\mathrm{right\ out}}(t), \]

so transmitted, reflected, and still-inside probability are counted consistently.

Four core formulas

  • Thermal state: \(\rho_{\mathrm{nuc}}=e^{-\beta H_{\mathrm{ho}}}/Z\).
  • Full Hamiltonian: \(\hat H=\hat P^2/(2m)\otimes I+\sum_{ab}|a\rangle V_{ab}(\hat X)\langle b|\).
  • Density propagation: \(\rho(t+\Delta t)=U\rho(t)U^\dagger\).
  • Population ledger: remain plus left-absorbed plus right-absorbed probability.

Benchmarks

The first numerical check is a free-particle density-matrix propagation. After thermal initialization and positive-momentum filtering, the ensemble moves right while \(\langle p\rangle\) stays constant under the free Hamiltonian.

Positive momentum density matrix free particle benchmark
Free density-matrix benchmark generated by the executable script for this note. Left: the blue curve is the ensemble mean position \(\langle x\rangle\) as a function of time. Middle: the orange curve is the conserved mean momentum \(\langle p\rangle\), confirming that the filtered thermal ensemble is right-moving under free evolution. Right: the position density \(\rho(x,x;t)\) at \(t=0\) in blue, \(t=20\) in orange, and \(t=50\) in green; spreading and translation are both visible.

The multi-state benchmark uses the Tully simple avoided crossing model. A preliminary right-transmission-only diagnostic is not a final scattering result: it can keep slow and reflected components inside the box, so the transmitted curves need not have reached an asymptote. The corrected calculation uses a smooth right-moving momentum filter centered at \(p=0.45\) and width \(0.08\). The filtered density matrix is diagonalized into 20 weighted pure-state components, retaining \(0.999999999385\) of the density weight, and both left and right absorbers are included in the population ledger.

Convergence checked density matrix ensemble benchmark for the Tully simple avoided crossing model
Tully simple avoided crossing density-matrix ensemble benchmark generated by the executable convergence script. Top left: remaining position density \(n(x,t)\), with blue for \(t=0\), orange for \(t=150000\), and green for \(t=300000\); the gray band is the left absorber \([-34,-27]\), and the pale blue band is the right absorber \([9,16]\). Top middle: residual in-box probability \(\mathrm{Tr}\,\rho_{\mathrm{remain}}\) on a log scale, falling to \(9.16\times10^{-4}\). Top right: remaining local adiabatic populations, blue for channel 0 and orange for channel 1. Bottom left: cumulative absorber ledger, with solid curves for left outflow, dashed curves for right outflow, blue for channel 0, and orange for channel 1; the final ledgers are \(P_L=(0.0344,0.0089)\) and \(P_R=(0.7360,0.2198)\). Bottom middle: remaining plus absorbed population by channel, which has plateaued to \(P_0=0.7713\) and \(P_1=0.2287\). Bottom right: conservation error \(P_{\mathrm{remain}}+P_{\mathrm{absorbed}}-1\), which stays below \(4\times10^{-13}\). Over the last 20% of the trajectory, each outflow channel drifts by less than \(1.8\times10^{-4}\), and the residual norm changes by \(5.0\times10^{-4}\).

The next notes reuse this operator bookkeeping for equilibrium Kubo correlation functions.

Code used in this note

  • dvr_ensemble_demo.py is a complete, executable density-matrix benchmark that regenerates the free-ensemble figure above.
  • dvr_sac_ensemble_convergence.py regenerates the corrected Tully SAC convergence figure and prints the residual norm, outflow ledgers, conservation error, and late-window drift.
  • dvr_kubo_minimal.py supplies the shared DVR Hamiltonian, Gaussian packet, and propagation utilities used in the surrounding notes.

DVR series map

  1. Part I: Wavepacket Dynamics Methods - discretization choices and the harmonic-oscillator benchmark.
  2. Part II: Operator Matrices in the DVR Basis - matrix structure behind kinetic, side, and flux operators.
  3. Part III: Wavepacket Dynamics with DVR - exact finite-dimensional propagation of a prepared wavepacket.
  4. Part IV: Density-Matrix Ensemble Evolution - thermal ensembles, filtering, and absorber bookkeeping.
  5. Part V: Kubo-Transformed Position Correlation - spectral thermal correlation functions.
  6. Part VI: Flux-Side Kubo Correlation - rate-like observables from side and flux operators.
  7. Part VII: Excited-Ground Kubo Population Correlation - local adiabatic projectors and population-resolved transfer diagnostics.