Matsubara Series - Part IV
绝热与非绝热 Matsubara 实时动力学推导
The scalar Matsubara formula combines a classical-looking nuclear generator with a complex quantum thermal weight. What changes when the potential becomes an electronic matrix? The nuclear Fourier transform and the origin of the phase remain the same, but the electronic trace must retain its order, and the force becomes a matrix anticommutator. This article derives the two position-correlation formulas together, keeping the electronic space discrete throughout.
1. The two formulas and the question they answer
Part I introduced the Kubo correlation and its phase-space representation. Part II introduced smooth imaginary-time modes and the Matsubara phase, while Part III examined a scalar estimator. Those ingredients do not yet explain where an electronic matrix belongs in the correlation function. We now reconstruct the necessary operator and thermal-kernel steps, using the same nuclear variables on both sides.
The target is the normalized nuclear position autocorrelation. For a single potential-energy surface, the Matsubara construction gives
Keeping the electronic space as matrices gives the corresponding expression
Here \(\beta=1/(k_BT)\), \(N\) is the imaginary-time bead number, and \(M\) is the odd number of retained nuclear modes. Their coordinates and momenta are \(Q_M,P_M\); \(Q_0\) is the nuclear centroid. The Hamiltonians, phase \(\theta_M\), cyclic electronic weight \(G_{e,N}\), and generators will be derived below. The electronic identity in the propagated observable is displayed explicitly in (2).
Equations (1) and (2) are Matsubara approximations, not alternative exact real-time quantum formulas. In (2), the electronic space still has dimension \(K^N\), where \(K\) is the number of physical electronic states. The notation suppresses an auxiliary \(N\) dependence: it denotes a family of expressions interpreted through the smooth-path construction, not an established finite-dimensional electronic limit. The derivation separates exact identities from the mode truncation and the limiting assumptions required to reach this compact form.
The scalar construction follows Hele, Willatt, Muolo, and Althorpe [1]; the nonadiabatic mode and thermal-kernel construction is compared with Chowdhury and Huo [2]. The explicit cyclic electronic matrix organization used here is a reconstruction of those steps. No additional classical electronic variables or trajectory closure are introduced.
2. Fix the model before changing its representation
To compare the derivations, the kinetic operator must be identical. We use one Cartesian nuclear coordinate \(R\), its conjugate momentum \(\hat p\), and nuclear mass \(m\). The two Hamiltonians are
The nonadiabatic potential \(\mathbf V(R)\) is a Hermitian \(K\times K\) matrix in a coordinate-independent electronic basis; \(V_0\) is the state-independent potential. A coordinate-dependent adiabatic basis would also introduce derivative couplings in the kinetic operator. The present derivation does not silently omit those terms: it fixes the coordinate-independent basis from the outset. “Adiabatic” here means a given scalar surface, not a derivation of the Born–Oppenheimer approximation.
Write either total potential as \(\mathsf V(R)\), a scalar or a matrix as appropriate. With \(Z=\operatorname{Tr}e^{-\beta\hat H}\), the common starting point is
The superscript \(K\) on the correlation means Kubo transformation; the separate symbol \(K\) for electronic dimension appears only in descriptions of the electronic space. The trace in (4) includes nuclei and, when present, electrons. This expression fixes both the initial ensemble and the real-time generator. Deriving a partition function alone would not fix the latter.
3. Turn the Kubo insertion into a cyclic bead average
The centroid factors in (1) and (2) must come from the Kubo insertion, rather than an assumed centroid estimator. Put \(\beta_N=\beta/N\) and \(E=e^{-\beta_N\hat H}\). A left-endpoint discretization of (4) is
At finite \(N\), this is a quadrature of the Kubo integral. To retain every insertion position without carrying its index through the derivation, replicate the complete Hilbert space \(N\) times. Define the left cyclic permutation \(\mathsf S\) by the trace identity
In a product basis, the left-hand side sums \((X_1)_{a_1a_2}(X_2)_{a_2a_3}\cdots(X_N)_{a_Na_1}\). Thus the permutation preserves a single ordered trace; it does not produce a product of independent traces.
Let \(\hat H^{(l)}\) act on copy \(l\), and define
Equation (6) gives \(\operatorname{Tr}_{\otimes}\Omega_N=Z\). Expanding \(A_NB_N(t)\) gives \(N^2\) insertion pairs. Each relative bead separation occurs \(N\) times, so their cyclic contraction reproduces (5):
Why use the trapezoidal rule here?
The next operation is a Wigner transform of the thermal position insertion. A one-sided operator product would leave a derivative of the thermal symbol. To obtain the simple centroid multiplier used in the final formula, we choose a symmetric insertion. This choice corresponds to a particular quadrature of the Kubo integral; the trapezoidal rule is convenient for this reason, rather than being an additional law of quantum dynamics.
Define the unnormalized Kubo integrand \(f(\lambda,t)=\operatorname{Tr}[e^{-(\beta-\lambda)\hat H}\hat R e^{-\lambda\hat H}\hat R(t)]\). The left sum is (5). The right sum samples \(k=1,\ldots,N\), so averaging the two gives
In the replicated trace, moving the first insertion across \(\Omega_N\) advances its imaginary-time position by one segment. Thus the two sums are represented by \(\Omega_N A_N\) and \(A_N\Omega_N\), respectively. Their average is
For a sufficiently smooth \(f\), this quadrature has an \(O(N^{-2})\) integration error at fixed \(\beta\), and converges to (4). We do not assume \(f(0,t)=f(\beta,t)\): cyclicity alone does not identify these differently ordered real-time insertions. The closed thermal path and the equality of those two integrand values are distinct questions. The replication is an exact identity for the chosen finite quadrature; it does not create \(N\) physical electrons or nuclei.
4. Transform the nuclei, leaving the electronic matrices intact
The replicated representation retains the exact time evolution but still contains nuclear operators. A nuclear Wigner transform makes the mode truncation possible without replacing electronic multiplication by scalar multiplication. For any replicated operator \(O\), use
The variables \(p_l\) are nuclear bead momenta and \(\Delta_l\) are nuclear endpoint differences. In the nonadiabatic case, \(O_W\) remains a matrix on the replicated electronic space. The associated nuclear star product is
For nonanalytic symbols the operator product defines the star product; an unrestricted derivative series is not assumed to converge. Heisenberg evolution then gives the exact phase-space generator
The scalar and matrix derivations have used the same transformation. Their difference is the electronic matrix multiplication inside (10) and (11), not a different nuclear Fourier transform.
Define \(Q_0=N^{-1}\sum_lR_l\). Both position insertions have Wigner symbol \(Q_0\mathbf I\). Because it is linear and electronic-scalar, the first derivatives in its two star products cancel, and all higher derivatives vanish:
In particular, before symmetrization the product \(\Omega_{N,W}\star_n Q_0\) contains \(-i\hbar(2N)^{-1}\sum_l\partial_{p_l}\Omega_{N,W}\). Reversing the product changes this sign. The cancellation in (12) explains the practical purpose of the trapezoidal choice: the first insertion becomes precisely \(Q_0\Omega_{N,W}\), without silently discarding a derivative term. It does not approximate the real-time propagator.
Consequently the finite trapezoidal correlation is exactly
For the scalar problem, remove the electronic trace. Equation (13) establishes why the final formulas contain two centroid insertions. It also shows why a general electronic observable cannot be substituted into them without rederiving its thermal insertion.
5. Resolve the thermal kernel before simplifying its electronic trace
To find the weight in (13), we need the endpoints of each imaginary-time segment. Define the exact one-segment electronic matrix kernel
The full cyclic permutation factors into nuclear and electronic permutations. With the direction fixed by (6), the nuclear Wigner transform is
This is a finite-\(N\) identity. In particular, the electronic kernel is not yet the local matrix \(e^{-\beta_N\mathbf V(R_l)}\).
For an explicit thermal reduction, choose the first-order product \(e^{-\beta_N\hat H}\simeq e^{-\beta_N\mathsf V(\hat R)}e^{-\beta_N\hat p^2/(2m)}\). It converges under the usual heat-kernel assumptions as \(N\) increases. This ordering makes the subsequent electronic placement explicit; symmetric splitting is also possible, but initially places half-potentials on both segment endpoints. The free nuclear kernel is
The nuclear displacement in segment \(l\) is therefore
Both problems have the same Gaussian \(\exp[-m\sum_la_l^2/(2\beta_N\hbar^2)]\). The scalar potential multiplies it as a number; the electronic potentials remain an ordered tensor product followed by \(\mathsf S_e\). This identifies the common nuclear calculation that will generate the phase.
6. Scale the normal modes and identify the dynamical approximation
6.1 Why a ring suggests Fourier modes
The Gaussian alone does not justify deleting nuclear modes from real-time evolution. We first need a basis that distinguishes slow and rapid variation around the imaginary-time ring. The trace closes the nuclear path, \(R_{l+N}=R_l\). A cyclic shift \((\mathsf T R)_l=R_{l+1}\) has the Fourier vectors \(v_{ln}=e^{2\pi i nl/N}\) as eigenvectors:
The second operator is the discrete ring Laplacian: its quadratic form is \(\sum_l(R_l-R_{l+1})^2\), exactly the free-particle thermal spring term. Fourier modes therefore both respect the cyclic boundary and diagonalize the universal nuclear Gaussian. A potential can still couple these modes; Fourier transformation does not diagonalize an arbitrary interacting Hamiltonian.
With imaginary time \(\tau_l=\beta\hbar l/N\), the same vector is \(e^{i\widetilde\omega_n\tau_l}\), where \(\widetilde\omega_n=2\pi n/(\beta\hbar)\). Small \(|n|\) means a path varying slowly along imaginary time. These are imaginary-time frequencies, not a selection of low real-time vibrational or electronic excitation frequencies. We transform the periodic nuclear path even in the matrix problem; no electronic antiperiodic boundary condition is being imposed on it.
For odd \(N\), take the real sine/cosine combinations of these eigenvectors:
Its orthogonality is \(N^{-1}\sum_lu_{ln}u_{lm}=\delta_{nm}\). Define the scaled modes, including endpoint-difference modes \(D_n\), by
The inverse relations for \(p_l\) and \(\Delta_l\) have the same form. Orthogonality gives the three factors that must remain consistent:
Here \(\Lambda\) denotes the antisymmetric bidifferential operator in (10). Thus the nuclear star parameter becomes \(\hbar/N\), while the summed Hamiltonian carries \(N\). The physical Planck constant has not been sent to zero.
6.2 Split all three integration variables before removing any
Equation (13) displays position and momentum integrals, but its thermal symbol (15) contains a third integral over endpoint differences. We must expose this integral to see what is eliminated. Define \(\mathcal I_N=\{-(N-1)/2,\ldots,(N-1)/2\}\), \(\mathcal I_M=\{-(M-1)/2,\ldots,(M-1)/2\}\), and \(\mathcal I_\perp=\mathcal I_N\setminus\mathcal I_M\). Use the same partition for each transformed variable:
Here \(Q\) describes the midpoint path, \(D\) the separation of the two endpoints used in the Wigner transform, and \(P\) is Fourier-conjugate to that separation. The normal-mode transformation has a constant Jacobian absorbed below into \(c_N\). Splitting a measure is an identity; it does not yet factor the interacting integrand into independent modes.
After the specified Trotter splitting, put \(x_l=R_l+\Delta_l/2\). The potential insertion on the replicated electronic space is
Writing \(\operatorname{tr}\) for the replicated electronic trace (or simply the identity operation in the scalar problem), define a common integral functional
The approximation sign here is the thermal Trotter step; no modes have yet been discarded. Equation (20c) makes the three integrations in the numerator and denominator explicit. The observable has no \(D\) argument because those variables belong to the initial thermal transform. It generally depends on all \(Q\) and \(P\), so the high momenta cannot yet be integrated as a bare Fourier factor.
6.3 The dynamical truncation permits the first integration
Retain the index set \(\mathcal I_M=\{-(M-1)/2,\ldots,(M-1)/2\}\). Split the exact generator by definition:
In \(\mathcal L_{M;N}\), nuclear derivatives and nuclear streaming are restricted to \(\mathcal I_M\); the difference includes high-mode and mixed derivative contributions. Equation (21) is the real-time Matsubara approximation. Its coefficients still depend on all nuclear coordinates. In particular, deleting derivatives is not the same as substituting a low-mode path into the potential.
Starting from \(Q_0\mathbf I\), this truncated generator does not create dependence on high-mode momenta, although it generally retains dependence on high-mode positions. This specific fact permits the first exact integration in the thermal reduction.
The semicolon emphasizes that \(Q_\perp\) still enters as a parameter in the force and electronic potential. Setting its derivative in the generator to zero freezes its dynamical role; it does not set its value to zero in the thermal ensemble.
7. Integrate the difference coordinates and derive the phase
7.1 High-mode momenta eliminate high-mode endpoint differences
Since the propagated observable no longer depends on \(P_n\) outside \(\mathcal I_M\), each such momentum occurs only in the thermal Fourier factor. Its integral is
This sets high-mode endpoint differences to zero. It does not set high-mode positions \(Q_n\) to zero. The distinction determines which remaining terms can be integrated next.
For all omitted modes together, the two successive integrations are
Here \(\Phi\) denotes the remaining thermal integrand, including its potential matrices. The Fourier integral is a distributional identity, understood under the remaining kernel integration. Thus \(P_\perp\) is integrated over, while \(D_\perp\) is constrained to zero. We have not made the replacement \(P_\perp=0\), nor have we integrated the physical midpoint coordinates \(Q_\perp\) yet.
7.2 The cyclic Gaussian fixes the phase sign
To evaluate the remaining low-mode integrals, put \(\xi_n=\pi n/N\) and transform (17). In a sine/cosine pair, shifting the bead index rotates the pair by \(2\xi_n\). The mean \((\Delta_l+\Delta_{l+1})/2\) supplies a factor \(\cos\xi_n\), while \(R_l-R_{l+1}\) supplies \(2\sin\xi_n\). Pairing the coefficients gives
Both signs of each mode index are included in the sums. For \(n=0\), the first term is simply \(ND_0^2\). Equation (23) is a finite discrete Fourier identity after (22); it is common to the scalar and matrix kernels.
The low-mode Gaussian has width of order \(N^{-1}\) and center \(-2\tan\xi_n Q_{-n}\), also of order \(N^{-1}\) at fixed retained index. Under regularity and domination assumptions, the slowly varying potential factors can therefore be evaluated at zero endpoint difference in the fixed-\(M\), large-\(N\) limit. The rapidly oscillatory Fourier factor must remain inside the integral. This is a thermal limiting step, not a finite-\(N\) identity.
After that step, set \(a_n=mN^2\cos^2\xi_n/(2\beta\hbar^2)\) and use the elementary Fourier transform of a Gaussian:
For example, shifting the integration variable to \(y=D_n+2\tan\xi_nQ_{-n}\) separates the phase from \(\int dy\,e^{-a_ny^2-iNP_ny/\hbar}\). This immediately produces (24), including the positive phase sign. No extra low-mode spring factor survives this completed Gaussian integral.
Now take fixed retained modes as \(N\to\infty\). Since \(\cos\xi_n\to1\) and \(2N\tan\xi_n/(\beta\hbar)\to\widetilde\omega_n=2\pi n/(\beta\hbar)\), the low-mode factor becomes
The common phase in (1) and (2) has now been derived. It is the remnant of the cyclic nuclear endpoint geometry, not an electronic approximation or an independently chosen correction to the distribution.
All difference coordinates are now gone, but for two different reasons: the high differences were eliminated by delta functions, and the low differences were integrated as shifted Gaussians. Setting \(D_M=0\) at the start would remove their Fourier integral and lose both the momentum Maxwell factor and the phase. Evaluating only the slowly varying potential factors at zero difference must not be confused with deleting the difference integral itself.
7.3 High-mode positions require a separate smooth-path limit
Equation (23) leaves a high-mode position weight
Potential matrices and the evolved observable still depend on these positions. An exact finite-\(M\) integration would average them with that dependence present; it would not evaluate them at \(Q_{\perp}=0\). We therefore keep this last integral visible before invoking the smooth-path reduction.
To locate the remaining dependence directly in \(C(t)\), write \(\mathcal T_N(Q)=\mathcal T_N(Q,D=0)\) using (20b), and define the finite-bead factors obtained from (24):
These finite-bead factors keep track of the Gaussian integration; the preceding evaluation of potential factors at \(D=0\) is still a limiting approximation. At this stage the correlation is
This is the missing bridge between the original three groups of integrals and the final two. It applies to both cases through the definition of \(\mathcal T_N\). The same Gaussian occurs in numerator and denominator, but it cannot cancel yet: it averages different functions of \(Q_\perp\). The nuclear potential, the electronic thermal matrices, and the propagated observable all retain that argument. Exact marginalization would keep their joint average, and does not in general replace the average of a propagator by the propagator of an averaged generator.
The smooth-path replacement instead removes \(Q_\perp\) from those arguments. It evaluates the thermal potential on the path below and constructs the reduced generator on that same path, rather than asserting that a conditional average has already done so:
The free high-mode covariance motivates this step: after taking the continuum bead limit, its omitted mean-square amplitude scales as \(\sum_{|n|>M/2}(\beta m\widetilde\omega_n^2)^{-1}=O(M^{-1})\). Turning that scale estimate into a replacement inside a thermal integral requires bounds on the potential and the observable. Doing so for a time-evolved, growing-dimensional matrix function requires further uniform control. Neither follows from the covariance estimate alone.
At fixed \(M\), the first omitted modes have finite continuum frequencies and finite Gaussian variances. Merely increasing \(N\) does not squeeze every omitted \(Q_n\) to zero. The decreasing tail estimate concerns increasing \(M\), with \(M\ll N\), and must be applied to the coupled integrand under suitable bounds.
We adopt the smooth-path reduction used in the Matsubara construction [1], [2], retaining this limitation explicitly. Once the non-Gaussian factors no longer depend on \(Q_\perp\), the remaining finite-\(N\) integral is
This common factor cancels before taking the continuum limit. The high positions have been integrated out, not set to zero as integration variables; only their occurrence in the other factors was replaced by the smooth path. This is the conditional step that permits the compact formulas, not a proof of universal real-time convergence.
7.4 Which variables survive?
The final measure contains \(Q_M,P_M\) because each of the other groups has now been accounted for. The operations and their logical status can be compared directly:
| Variables | Operation | Result and condition |
|---|---|---|
| \(P_\perp\) | Fourier integration | Produces \(\delta(D_\perp)\), after the dynamical truncation removes high-momentum dependence from the observable. |
| \(D_\perp\) | Delta-function integration | Sets high endpoint differences to zero; leaves midpoint positions intact. |
| \(D_M\) | Shifted Gaussian integration | Produces kinetic weight and phase, after the stated low-difference thermal limit. |
| \(Q_\perp\) | Gaussian integration after smooth-path replacement | Produces \(\mathcal N_\perp\), which cancels. The replacement is not exact finite-\(M\) marginalization. |
| \(Q_M,P_M\) | No elimination | Remain in the initial ensemble and the real-time propagated observable. |
The order matters: truncating the propagator enables the high-momentum integration; that integration constrains the high differences; the low-difference Gaussian creates the phase; and only a separate smooth-path argument permits the final high-position cancellation.
8. Assemble the scalar and matrix thermal weights
The nuclear integrations have produced kinetic energy and phase but have not removed the potential factors. Their remaining form distinguishes the two correlations. For the scalar surface, define the continuous imaginary-time average
The product of scalar potential factors combines with (25) to give \(W_M^{\mathrm a}=e^{-\beta(H_M^{\mathrm a}-i\theta_M)}\).
For the matrix problem, only the state-independent part can be collected into the scalar Hamiltonian. Write
The scalar average tends to its continuous-path integral as \(N\) increases. The electronic operator left by (15) and our specified Trotter ordering is
Tracing (30) by the index contraction in (6) yields
The right-hand side does not in general equal \(\prod_l\operatorname{Tr}_e E_l\); for noncommuting potentials, it also cannot in general be replaced by \(\operatorname{Tr}_e\exp[-\beta_N\sum_l\mathbf V_l]\). The cyclic permutation is therefore essential even though the real-time electronic Hamiltonian is a sum over copies.
For a general evolved matrix observable \(F\), one needs \(\operatorname{Tr}(G_{e,N}F)\), not merely \(\operatorname{Tr}G_{e,N}\). Replacing the operator by its scalar trace before propagation would lose the electronic information needed by the correlation. Also, \(G_{e,N}\) need not be Hermitian or positive: it is a cyclic contraction weight. The relevant pairing is \(\int\operatorname{Tr}(WF)\), without inserting a Hermitian conjugate or interpreting \(W\) as a positive sampling density.
9. Derive the two retained-mode generators
The thermal weights are now identified. To finish (1) and (2), we must specify what acts on the second \(Q_0\). In normal modes define the averaged potential on the replicated electronic space, with \(\mathbf V^{(l)}\) acting only on factor \(l\):
Insert (32) and the star parameter \(\hbar/N\) into (11), restricting nuclear derivatives as in (21). To display the surviving and discarded structures, expand in these retained derivatives. With subscripts denoting differentiation, this gives
For instance, the first star correction has coefficient \((i/\hbar)N(i\hbar/2N)=-1/2\); adding the two product orders gives the force anticommutator. The second has coefficient \((i/\hbar)N(-\hbar^2/8N^2)=-i\hbar/(8N)\); subtracting the product orders gives the displayed commutator. These coefficients distinguish a matrix potential from a scalar one.
For fixed \(M\) and function families with controlled derivatives and remainders, the higher nuclear terms vanish as \(N\) grows. In the matrix case, the leading displayed correction can be \(O(\hbar/N)\); in the scalar case that commutator vanishes, and the first remaining correction is \(O(\hbar^2/N^2)\). These are local generator estimates, not a uniform convergence proof for the propagator on a \(K^N\)-dimensional electronic space.
Combining this nuclear low-order limit with the already specified smooth-path replacement gives the matrix generator
In the scalar problem, all commutators vanish and the force anticommutator becomes twice ordinary multiplication. The same derivation gives
Equation (35) is a scalar Hamiltonian flow. Equation (34) is a matrix-valued partial differential operator. It has not been closed into one nuclear trajectory and one electronic wavefunction per bead. The absence of a nuclear spring force in both generators is consistent with the low-mode Gaussian integration: the statistical phase has been retained, rather than converted into an approximate real ring-polymer propagation.
10. Insert the weight and generator into the correlation
We can now return to (13). Its two insertions are fixed by (12), its reduced thermal weight by (25), (28), and (30), and its propagation by (34) or (35). All coordinate Jacobians and observable-independent Gaussian constants cancel in the normalized ratio. Write \(dQ_M\,dP_M=\prod_{n\in\mathcal I_M}dQ_n\,dP_n\).
In the scalar case the substitution gives
In the nonadiabatic case, keep the trace until after applying the matrix propagator:
Here \(H_{0,M}\) abbreviates the scalar factor in (29), with its bead quadrature understood in the stated limit. Square brackets specify the propagator's operand: it acts on the second position insertion, not on \(G_{e,N}\) or the initial \(Q_0\). Thus (36) and (37) are precisely the target formulas (1) and (2), with the electronic identity and the auxiliary bead dependence made explicit.
At \(t=0\), their numerators contain \(Q_0^2\). They therefore give the centroid-square estimator for the Kubo-transformed position correlation, not the equal-time ordinary position variance. A finite-mode weight is still a truncation; recovery of the full quantum static quantity requires the corresponding static convergence.
10.1 What changes in the nonadiabatic correlation?
The two final ratios look similar, so we must identify which changes carry electronic physics. Both retain exactly the same nuclear integration variables, linear centroid insertions, and Matsubara phase. Their difference enters in both the thermal operator and the propagated observable:
| Ingredient in \(C(t)\) | Scalar surface | Electronic matrix potential |
|---|---|---|
| Nuclear variables and phase | \(Q_M,P_M,\theta_M\) | The same nuclear variables and phase |
| Potential in the thermal weight | One scalar factor \(e^{-\beta U_M^{\mathrm a}}\) | \(e^{-\beta U_{0,M}}G_{e,N}\), retaining the cyclic electronic contraction |
| Time-evolved position insertion | The scalar function \(e^{t\mathcal L_M^{\mathrm a}}Q_0\) | The matrix function \(F_t=e^{t\mathcal L_M^{\mathrm{mat}}}(Q_0\mathbf I)\) |
| Nuclear force term | Ordinary multiplication by a potential gradient | Anticommutator with a matrix gradient |
| Electronic evolution | No electronic commutator | \(i[\mathsf V_\Sigma,F_t]/\hbar\), with \(\mathsf V_\Sigma=\sum_l\mathbf V^{(l)}\) |
| Final contraction | Multiply scalar weight and observable | Take \(\operatorname{Tr}(G_{e,N}F_t)\) after propagation |
In particular, starting with a nuclear-only observable does not make the propagated observable electronic-scalar. Apply (34) twice to its initial value:
The first derivative is scalar because a position operator initially streams with its momentum. The second derivative contains the electronic-state-dependent force. It is generally a matrix; subsequent evolution can also act on it through electronic commutators. In the scalar problem the corresponding force is just \(-\partial_{Q_0}U_M^{\mathrm a}\), and the observable remains a scalar along a Hamiltonian flow.
This distinction already appears in the initial curvature of the correlation. With \(w_0=e^{-\beta(H_{0,M}-i\theta_M)}\), the matrix expression gives, when differentiation under the integral is valid,
The force is contracted with the electronic thermal operator, not introduced as an independently chosen scalar surface. Equation (37b) is an analytic consequence of the proposed generator and ensemble; it does not by itself establish agreement with the exact quantum curvature.
10.2 Why tracing the thermal weight first does not finish the dynamics
In the denominator of (37), \(\operatorname{Tr}G_{e,N}\) is sufficient because the inserted observable is the identity. At \(t=0\) the numerator also has a scalar insertion, \(Q_0^2\mathbf I\). At later times (37a) shows why that simplification generally stops working:
Where \(\operatorname{Tr}G_{e,N}\ne0\), one could define a scalar ratio \(f_t:=\operatorname{Tr}(G_{e,N}F_t)/\operatorname{Tr}G_{e,N}\). That is a rewriting of the result, not a closed scalar equation of motion: computing it still requires the matrix evolution, and no autonomous Hamiltonian flow for it has been established. A replacement by one electronic state, a mean force, or independent bead electronic trajectories would require further assumptions. The ordered thermal weight and the matrix propagation together constitute the nonadiabatic difference in \(C(t)\).
11. Analytic checks of the result
11.1 One electronic state recovers the scalar formula
A matrix extension must reduce to the scalar expression when its electronic space has one state. For \(K=1\), the cyclic electronic permutation is trivial and
This factor combines with \(e^{-\beta H_{0,M}}\) into the scalar weight for the total surface \(V_0+V\). The commutator in (34) vanishes, and its anticommutator becomes the force in (35). Both numerator and denominator of (37) therefore reduce to (36). The thermal and dynamical reductions agree, rather than merely producing the same partition function.
11.2 Decoupled electrons cancel from a nuclear correlation
The next check keeps multiple electronic states but removes their nuclear coupling: \(\mathbf V(R)=\mathbf V_*\), a constant matrix. Then the electronic contribution to the nuclear force vanishes. Starting from \(Q_0\mathbf I\), the observable remains nuclear-scalar, so its electronic commutator also remains zero. Meanwhile, (31) gives
This constant cancels between numerator and denominator. The position correlation is exactly the scalar Matsubara expression for \(V_0\), as required by electronic–nuclear separation. No spurious electronic degeneracy factor remains.
11.3 The scalar harmonic correlation fixes the normalization
To test the distinction between Kubo and ordinary variance analytically, take \(V_{\mathrm a}(R)=m\Omega^2R^2/2\), where \(\Omega\) is the physical oscillator frequency. Fourier orthogonality gives \(U_M^{\mathrm a}=m\Omega^2\sum_nQ_n^2/2\). The phase couples only internal modes because \(\widetilde\omega_0=0\); their integrals cancel in the normalized centroid correlation. The centroid satisfies
Equation (36) consequently yields
This is the exact harmonic Kubo position correlation. Its value at zero time differs from the ordinary quantum variance \(\hbar\coth(\beta\hbar\Omega/2)/(2m\Omega)\). The check confirms the observable definition and normalization; it does not establish accuracy for anharmonic coupled systems.
12. What has been derived, and what does not follow
The two formulas share their nuclear Fourier geometry, phase, and linear position insertion. The nonadiabatic expression adds an ordered electronic thermal contraction and a matrix real-time generator. Keeping both is essential: replacing the electronic trace by independent bead populations, or the matrix force by an expectation value, would introduce another approximation beyond this article.
The logical status of the steps is now explicit. Cyclic replication, nuclear Wigner transformation, linear-insertion cancellation, and the Gaussian integral (24) are identities at their stated stages. The heat-kernel splitting has a Trotter limit. Deleting high-mode derivatives is a dynamical approximation. Replacing the remaining high-coordinate dependence by the smooth path requires the separate limiting assumptions in section 7.3. The nuclear star expansion then supplies the retained generator under controlled derivative assumptions. These operations do not prove that increasing \(M\) recovers arbitrary exact real-time dynamics.
Finally, a thermal weight appearing in a correlation formula is not by itself a proof that the approximate propagation preserves it. For the scalar generator, the continuous imaginary-time translation symmetry gives \(\mathcal L_M^{\mathrm a}\theta_M=0\), and energy conservation gives \(\mathcal L_M^{\mathrm a}H_M^{\mathrm a}=0\). The matrix case requires the entire cyclic weight. With vanishing integration boundary terms, the formal dual of (34) under \(\int\operatorname{Tr}(WF)\) is
A general conservation claim would require proving \((\mathcal L_M^{\mathrm{mat}})^\vee W_{M;N}^{\mathrm{na}}=0\) in the appropriate limit. Neither the shared nuclear phase nor the mere existence of (37) proves that identity. The completed result here is the conditional matrix correlation construction and its scalar, decoupled, and harmonic consistency checks.
References and derivation scope
- T. J. H. Hele, M. J. Willatt, A. Muolo, and S. C. Althorpe, “Boltzmann-conserving classical dynamics in quantum time-correlation functions: ‘Matsubara dynamics’,” Journal of Chemical Physics 142, 134103 (2015). University of Cambridge. DOI; author-hosted article. Relevant material: the normal-mode construction, Eqs. (69)–(82), and Appendix C on the thermal reduction.
- S. N. Chowdhury and P. Huo, “Non-adiabatic Matsubara dynamics and non-adiabatic ring-polymer molecular dynamics,” Journal of Chemical Physics 154, 124124 (2021). University of Rochester. DOI; author-hosted article. Relevant material: the nuclear-linear-observable restriction preceding Eq. (33), Eqs. (52)–(68), and Appendix B. The body text points to Appendix C for high-mode elimination, but the actual derivation is in Appendix B; Appendix C treats an electronic–nuclear decoupling case.
The cyclic replicated-space formulation, explicit matrix star expansion, and the Gaussian evaluation in this note are presented as a direct derivation with stated conventions. They should not be read as quotations of a separately established general matrix convergence theorem. This article contains analytical results only; it makes no numerical reproduction claim.