Reaction-Coordinate Mapping

Reaction-coordinate mapping is a way to redraw the boundary between a quantum system and its environment. A collective environmental coordinate that is strongly coupled to the system, slowly relaxing, or spectrally structured is promoted into the system. The remaining residual bath is then treated as weaker, broader, and more nearly Markovian.

The guiding picture is:

original system + structured bath
    = enlarged system + residual bath

The original reduced system may show memory. The enlarged system can often be described by a time-local master equation. The memory has not disappeared; it has been carried by the reaction coordinate, which is now explicit.

For a closely related pole-based construction, see Pseudomode Methods. Pseudomodes are often tied to Lorentzian response poles or leaky resonances. Reaction-coordinate mappings emphasize collective coordinates of a bosonic environment and are especially useful for strong-coupling and thermodynamic questions.

A system coupled to a bosonic bath

Take a dimensionless Hermitian system operator $A$ and real mode couplings $g_k$ with units of energy. A common bosonic system–bath model is

$$H = H_S + A\otimes B + H_B,$$

with

$$B = \sum_k g_k \left( b_k+b_k^\dagger \right), \qquad H_B = \sum_k \hbar\omega_k b_k^\dagger b_k.$$

Use the spectral-density convention

$$J(\omega) = \sum_k g_k^2\delta(\omega-\omega_k), \qquad \omega>0.$$

Thus $J$ has units of energy squared times time, because $\omega$ is an angular frequency. If $J(\omega)$ is broad and smooth near the system transition frequencies, a weak-coupling Markovian master equation may be adequate. If $J(\omega)$ has a narrow peak, a strong low-frequency component, or an underdamped vibrational structure, the system may retain correlations with a collective environmental degree of freedom. The reaction-coordinate mapping makes that degree of freedom explicit.

Mapped Hamiltonian

After the mapping, one writes an enlarged Hamiltonian of the form

$$H = H_{S+\mathrm{RC}} + H_{\mathrm{res}} + H_{\mathrm{RC-res}},$$

where

$$H_{S+\mathrm{RC}} = H_S + \hbar\Omega a^\dagger a + \lambda A \left( a+a^\dagger \right) + H_{\mathrm{ct}}.$$

Here $a$ is the reaction-coordinate mode, $\Omega$ is its frequency, $\lambda$ is its coupling to the system operator $A$, and $H_{\mathrm{ct}}$ denotes any counterterm or renormalization term required by the convention.

The residual bath has modes $c_q$:

$$H_{\mathrm{res}} = \sum_q \hbar\nu_q c_q^\dagger c_q,$$

and the reaction coordinate couples to it through

$$H_{\mathrm{RC-res}} = \left( a+a^\dagger \right) \sum_q h_q \left( c_q+c_q^\dagger \right).$$

The mapping is chosen so that tracing out the reaction coordinate and residual bath reproduces the original bath influence on $S$. The residual spectral density is not arbitrary; it is determined by the original $J(\omega)$ and the mapping convention. The initial bath state must be transformed as well: an original thermal bath generally becomes a correlated reaction-coordinate–residual-bath state. Replacing that state by a product is an additional approximation, even if the Hamiltonian mapping is exact.

An exact two-mode coordinate mapping

A bath with two oscillators makes the change of variables explicit. Introduce mass-weighted coordinates

$$x_k=\sqrt{\frac{\hbar}{2\omega_k}}(b_k+b_k^\dagger), \qquad p_k=-i\sqrt{\frac{\hbar\omega_k}{2}}(b_k-b_k^\dagger), \qquad k=1,2.$$

Then $[x_j,p_k]=i\hbar\delta_{jk}$ and $B=f_1x_1+f_2x_2$, where $f_k=g_k\sqrt{2\omega_k/\hbar}$. For $F=\sqrt{f_1^2+f_2^2}>0$, rotate both coordinates and momenta by the same orthogonal matrix:

$$\begin{pmatrix}Q\\R\end{pmatrix} =\frac1F\begin{pmatrix}f_1&f_2\\-f_2&f_1\end{pmatrix} \begin{pmatrix}x_1\\x_2\end{pmatrix}, \qquad \begin{pmatrix}\Pi_Q\\\Pi_R\end{pmatrix} =\frac1F\begin{pmatrix}f_1&f_2\\-f_2&f_1\end{pmatrix} \begin{pmatrix}p_1\\p_2\end{pmatrix}.$$

This preserves the canonical commutators and gives $B=FQ$: only $Q$ couples directly to $S$. Direct substitution into the bath Hamiltonian gives

$$H_B=\frac{\Pi_Q^2+\Omega^2Q^2+\Pi_R^2+\nu^2R^2}{2} +KQR-\frac{\hbar(\omega_1+\omega_2)}2,$$

where

$$\Omega^2=\frac{f_1^2\omega_1^2+f_2^2\omega_2^2}{F^2}, \qquad \nu^2=\frac{f_2^2\omega_1^2+f_1^2\omega_2^2}{F^2}, \qquad K=\frac{f_1f_2(\omega_2^2-\omega_1^2)}{F^2}.$$

The scalar term retains the original convention of omitting zero-point energies. Writing $Q=\sqrt{\hbar/(2\Omega)}(a+a^\dagger)$ and $R=\sqrt{\hbar/(2\nu)}(c+c^\dagger)$ recovers the position-coupling form above, with $\lambda=F\sqrt{\hbar/(2\Omega)}$ and residual coupling $h=\hbar K/(2\sqrt{\Omega\nu})$. Existing system counterterms stay part of the same Hamiltonian; this rotation does not add a physical interaction.

The two normal frequencies remain $\omega_1,\omega_2$, as checked by $\Omega^2\nu^2-K^2=\omega_1^2\omega_2^2>0$. In the original spectral convention these coefficients also obey

$$\Omega^2=\frac{\sum_{k=1}^2g_k^2\omega_k^3}{\sum_{k=1}^2g_k^2\omega_k}, \qquad \lambda^2=\frac{\sum_{k=1}^2g_k^2\omega_k}{\Omega}.$$

These are the finite-mode counterparts of the phonon-mapping moments in Nazir and Schaller (2018), Sec. II.B.1, with $\hbar$ restored. They require the corresponding spectral moments to exist when extended to a continuum. No rotating-wave or Markov approximation entered this example. In particular, the remaining oscillator is a finite memory carrier, not automatically a reservoir.

Treating the Strong Interaction Explicitly

The coupling $\lambda A(a+a^\dagger)$ is part of the enlarged Hamiltonian and can be treated nonperturbatively. Only the remaining interaction with the residual bath is expanded when deriving a master equation. The residual bath may have a shorter correlation time and weaker coupling, but both properties must be checked; an exact change of variables alone does not establish a Markov approximation.

Brownian-Oscillator Example

For a coordinate damped by an Ohmic residual bath, with its frequency counterterm included and damping constant $\gamma$ in the oscillator equation $\ddot x+\gamma\dot x+\Omega^2x=0$, the corresponding Brownian-oscillator spectrum is

$$J_{\mathrm{BO}}(\omega) = \frac{ 2\lambda^2\gamma\Omega\omega/\pi } { (\Omega^2-\omega^2)^2+\gamma^2\omega^2 }, \qquad \omega>0.$$

The parameters suggest the mapping directly:

• $\Omega$ is the reaction-coordinate frequency;
• $\lambda$ sets the system-reaction-coordinate coupling strength;
• $\gamma$ sets the damping of the reaction coordinate by the residual bath.

In the mapped model, the system couples to a single oscillator of frequency $\Omega$, and that oscillator is damped by a smoother residual bath. If $\gamma$ is small, the reaction coordinate stores excitation for a long time and memory is prominent. Eliminating the coordinate also requires its actual correlation time to be short on system scales. In particular, very strong overdamping can produce a slow relaxation rate of order $\Omega^2/\gamma$; increasing $\gamma$ indefinitely does not guarantee rapid forgetting.

The numerator has three powers of angular frequency and the denominator has four, so the result has the required units $\lambda^2\times\text{time}$. In the narrow-width limit it tends to $\lambda^2\delta(\omega-\Omega)$ on positive frequencies. These checks fix the normalization for the convention above. The Ohmic expression assumes the appropriate large residual cutoff; finite cutoffs can add frequency-dependent shifts. It is the form derived by Iles-Smith, Lambert, and Nazir, with their dimensionless residual coupling converted to the damping rate $\gamma$.

Master Equation on the Enlarged System

Once the reaction coordinate is included in the system, the density matrix $\varrho$ for $S+\mathrm{RC}$ obeys an equation schematically like

$$\frac{d\varrho}{dt} = - \frac{i}{\hbar} [H_{S+\mathrm{RC}},\varrho] + \mathcal L_{\mathrm{res}}(\varrho).$$

The residual dissipator $\mathcal L_{\mathrm{res}}$ should be derived from the coupling between the reaction coordinate and the residual bath. In a careful weak-coupling treatment, one diagonalizes $H_{S+\mathrm{RC}}$ and constructs jump operators at the Bohr frequencies of the enlarged Hamiltonian.

This is important. A local dissipator such as

$$\kappa\mathcal D[a]\varrho$$

may be a useful approximation in some regimes, but it is not automatically thermodynamically consistent when $a$ is strongly hybridized with $S$. The residual bath sees the dressed eigenstates of the enlarged system, not the bare reaction coordinate in isolation.

Reduced Memory

The original system state is obtained by tracing out the reaction coordinate:

$$\rho_S(t) = \operatorname{Tr}_{\mathrm{RC}}\varrho(t).$$

Even if $\varrho(t)$ follows a Markovian master equation, $\rho_S(t)$ can show memory because information can move from $S$ into the reaction coordinate and later return. The reaction coordinate stores correlations that would otherwise be absent from the system state.

This is the same structural lesson as Pseudomode Methods: Markovianity depends on where the system boundary is drawn. The enlarged description may be CP-divisible while the smaller reduced description is not.

Strong-Coupling Equilibrium

Reaction-coordinate mappings are especially useful for equilibrium and thermodynamic questions at strong coupling.

If the enlarged system thermalizes with a residual bath at inverse temperature $\beta$, its equilibrium state may be close to

$$\varrho_{\mathrm{eq}} = \frac{ e^{-\beta H_{S+\mathrm{RC}}} } { \operatorname{Tr}_{S,\mathrm{RC}} e^{-\beta H_{S+\mathrm{RC}}} }.$$

The reduced system equilibrium is then

$$\rho_S^{\mathrm{eq}} = \operatorname{Tr}_{\mathrm{RC}} \varrho_{\mathrm{eq}}.$$

This is generally not the bare Gibbs state

$$\frac{e^{-\beta H_S}}{\operatorname{Tr}_S e^{-\beta H_S}}.$$

It is instead related to a Hamiltonian of mean force. This distinction explains why demanding a bare-Gibbs steady state can be wrong outside the weak-coupling limit. See Thermal Master Equations for the weak-coupling thermal case.

Relation to Chain Mappings

Reaction-coordinate mapping can be viewed as the first step of a star-to-chain transformation. The original bath modes form a star around the system. Orthogonal-polynomial or Lanczos methods transform the bath into a chain:

system - site 0 - site 1 - site 2 - ...

The first site is the reaction coordinate. Keeping only that site in the enlarged system and treating the rest as a residual bath gives the basic reaction-coordinate picture. Keeping more sites gives a longer explicit memory chain and can improve accuracy for strongly structured environments.

This perspective is common in tensor-network simulations and numerically exact open-system methods. The tradeoff is transparent: more explicit bath sites mean less memory left in the residual bath but a larger system Hilbert space.

Practical Workflow

A reaction-coordinate calculation usually follows this sequence.

1. Specify the original $H_S$, coupling operator $A$, bath state, and spectral density $J(\omega)$.
2. Choose a reaction coordinate that captures the dominant structured part of the bath.
3. Derive or fit $\Omega$, $\lambda$, counterterms, and the residual spectral density.
4. Build $H_{S+\mathrm{RC}}$ and choose a reaction-coordinate Hilbert-space cutoff.
5. Derive $\mathcal L_{\mathrm{res}}$ for the residual bath, preferably in the dressed eigenbasis when coupling is strong.
6. Evolve $\varrho(t)$ or compute its steady state.
7. Trace out the reaction coordinate to obtain $\rho_S(t)$ and system observables.
8. Test convergence against cutoff, mapping depth, residual-bath approximation, and known limiting cases.

The method is constructive, but not automatic. The residual bath assumptions and numerical truncations are part of the model.

When to Use It

Reaction-coordinate mapping is useful when:

• one environmental feature is strongly coupled to the system;
• the bath spectrum has a pronounced peak or underdamped mode;
• weak-coupling rates built from the original $J(\omega)$ fail;
• equilibrium at strong coupling matters;
• the environment is better viewed as a slow coordinate plus a faster background;
• a finite explicit memory mode makes the model easier to simulate or interpret.

It is less natural when the environment has broad algebraic memory, many comparable peaks, a strongly non-Gaussian spin structure, or no clear coordinate to isolate. In those cases, hierarchical methods, tensor-network chain simulations, collision models, or direct microscopic modeling may be more appropriate.

For Gaussian baths whose correlation functions can be expanded into exponentials, Hierarchical Equations of Motion provide a systematic auxiliary-density-operator alternative to choosing a small number of explicit coordinates.

Exercises

Mapping the initial bath state

In the two-mode example, let the original bath be the product thermal state at inverse temperature $\beta$. Derive $\langle QR\rangle$. Then take its zero-temperature limit with $f_1=f_2\ne0$, $\omega_1=\omega_0$, and $\omega_2=2\omega_0$. Find the change in bath energy if the mapped state is replaced by the product of its exact $Q$-mode and $R$-mode marginals. Does that replacement merely change coordinates?

Solution

The original means and cross correlations vanish, while $\langle x_k^2\rangle=\hbar\coth(\beta\hbar\omega_k/2)/(2\omega_k)$. Rotating gives

$$\langle QR\rangle =\frac{f_1f_2}{F^2} \left[ \frac{\hbar}{2\omega_2}\coth\!\left(\frac{\beta\hbar\omega_2}{2}\right) -\frac{\hbar}{2\omega_1}\coth\!\left(\frac{\beta\hbar\omega_1}{2}\right) \right].$$

For the specified vacuum, $\Omega^2=\nu^2=5\omega_0^2/2$ and $K=3\omega_0^2/2$. Its mapped moments are

$$\langle Q^2\rangle=\langle R^2\rangle=\frac{3\hbar}{8\omega_0}, \qquad \langle QR\rangle=-\frac{\hbar}{8\omega_0},$$ $$\langle\Pi_Q^2\rangle=\langle\Pi_R^2\rangle=\frac{3\hbar\omega_0}{4}, \qquad \langle\Pi_Q\Pi_R\rangle=\frac{\hbar\omega_0}{4}.$$

Thus even the original product vacuum is correlated across this new division of modes. The product of the mapped marginals keeps all local moments but sets $\langle QR\rangle$ to zero. Since the only coupling term in $H_B$ is $KQR$, it raises the bath energy by

$$\Delta E_B=-K\langle QR\rangle=\frac{3\hbar\omega_0}{16}.$$

This is a change of preparation. It is not part of the exact coordinate transformation, and it need not reproduce the original bath influence. The energy difference is independent of the scalar zero-point convention.

Memory from an Explicit Coordinate

Explain how a Markovian master equation for $S+\mathrm{RC}$ can lead to non-Markovian reduced dynamics for $S$.

Solution

The enlarged state $\varrho(t)$ contains correlations between $S$ and the reaction coordinate. At an intermediate time $s$, the reduced state

$$\rho_S(s)=\operatorname{Tr}_{\mathrm{RC}}\varrho(s)$$

does not determine the full $\varrho(s)$. The future of $S$ can depend on excitation or phase information stored in the reaction coordinate. Therefore there may be no completely positive intermediate map acting only on $\rho_S(s)$, even though $S+\mathrm{RC}$ evolves Markovianly.

Local Versus Dressed Dissipation

Why can a local dissipator $\kappa\mathcal D[a]$ be unreliable when $\lambda A(a+a^\dagger)$ is large?

Solution

The residual bath couples to the reaction-coordinate coordinate, but the energy eigenstates of the enlarged system are eigenstates of

$$H_{S+\mathrm{RC}} = H_S+\hbar\Omega a^\dagger a+\lambda A(a+a^\dagger).$$

When $\lambda$ is large, these eigenstates are dressed mixtures of system and reaction-coordinate states. Thermal transitions induced by the residual bath should occur at Bohr frequencies of $H_{S+\mathrm{RC}}$. A bare local dissipator may impose rates at the wrong frequencies and can give an incorrect steady state or heat current.

One More Chain Site

In a chain mapping, what is the expected effect of keeping two explicit bath sites instead of only the first reaction coordinate?

Solution

Keeping a second site moves more environmental memory into the enlarged system. The residual bath seen by the second site is often smoother and shorter-correlated than the residual bath seen by the first site. The price is a larger Hilbert space and more numerical cost. If results change significantly when the second site is added, the one-coordinate mapping plus residual-bath approximation was not yet converged.

References

• A. W. Chin, A. Rivas, S. F. Huelga, and M. B. Plenio, “Exact mapping between system-reservoir quantum models and semi-infinite discrete chains using orthogonal polynomials,” Journal of Mathematical Physics 51, 092109 (2010).
• I. de Vega and D. Alonso, “Dynamics of non-Markovian open quantum systems,” Reviews of Modern Physics 89, 015001 (2017).
• A. Garg, J. N. Onuchic, and V. Ambegaokar, “Effect of friction on electron transfer in biomolecules,” Journal of Chemical Physics 83, 4491-4503 (1985).
• J. Iles-Smith, N. Lambert, and A. Nazir, “Environmental dynamics, correlations, and the emergence of noncanonical equilibrium states in open quantum systems,” Physical Review A 90, 032114 (2014), doi:10.1103/PhysRevA.90.032114.
• R. Martinazzo, B. Vacchini, K. H. Hughes, and I. Burghardt, “Communication: Universal Markovian reduction of Brownian particle dynamics,” Journal of Chemical Physics 134, 011101 (2011).
• A. Nazir and G. Schaller, “The reaction coordinate mapping in quantum thermodynamics,” in Thermodynamics in the Quantum Regime: Fundamental Aspects and New Directions, edited by F. Binder, L. A. Correa, C. Gogolin, J. Anders, and G. Adesso, Springer (2018), pp. 551–577, Sec. II.B.1. DOI; author manuscript.
• P. Strasberg, G. Schaller, N. Lambert, and T. Brandes, “Nonequilibrium thermodynamics in the strong coupling and non-Markovian regime based on a reaction coordinate mapping,” New Journal of Physics 18, 073007 (2016).