Dynamics with Simultaneous Dissipations to Fermionic and Bosonic Reservoirs

Open systems evolve over time interacting with fermionic or bosonic reservoirs representing electrons in electrodes, lattice vibrations or electromagnetic fields. In electrochemical systems, the dynamics occurs under a moderate non-equilibrium condition influenced by a delicate balance of the multiple reservoirs. Particularly important property to characterize the systems is the rate of exchanging the energy and charge with reservoirs, which can be investigated in principle using a microscopic theory-like based on the Langevin equationLangevin (1908); Lemons and Gythiel (1997); Ford and Kac (1987); Ford et al. (1988) by incorporating the fluctuations and dissipation due to environments. The multiple reservoir effect, however, has not received much attention.

Over the years, the Langevin equation has been formally derivedFord et al. (1965); Schmid (1982) and re-derivedCaldeira and Leggett (1983a); Metiu and Schon (1984); Newns (1985); Sebastian (1987); Ford and Kac (1987); Ford et al. (1988) in various contexts, accounting for the dissipation of the impinging particle’s kinetic energy not only into substrate phonon excitations via phononic frictionNewns (1985), but also into the generation of coherent electron-hole (e-h) pairs via electronic frictionBrako and Newns (1980); Schonhammer and Gunnarsson (1980); Brandbyge et al. (1995); Head‐Gordon and Tully (1995); Plihal and Langreth (1998); Dou et al. (2017); Martinazzo and Burghardt (2022). The equation has subsequently been generalized for the formerNewns (1985); Sebastian (1987) and the latterHedegård and Caldeira (1987); Brandbyge et al. (1995), but their simultaneous operation has mostly been neglected, often justified by arguments based on the disparity of particle masses. This is, however, not always the case because the strength of the former is considerably dependent on the environments; it is reduced, for example, when the reaction center moves away from the electrode and correspondingly electrons must transfer farther away therefrom, as typically occurs with increasing pHKumeda et al. (2023); Li et al. (2022). The crossover of the dominant dissipation channels thus occurs in nature.

To address this issue, we reformulate the particle dynamics using the influence functional path integralFeynman and Vernon (2000) framework, which yields without phenomenological assumptions a quasiclassical Langevin equation that incorporates non-Markovian effects of both reservoirs. The dissipation kernel and fluctuations become Markovian in the limit of slow particle motion, allowing analytical evaluation of the electronic friction coefficient. We demonstrate this through numerical studies of (1) quantum vibrational relaxation of hydrogen on metal surfaces and (2) solvated proton discharge dynamics on metal electrodes. Considering the generality of this scheme, it should extend the frontier of the research on quasi-equilibrium systems.

The total Hamiltonian of a system comprising a particle coupled to two thermal reservoirs is expressed as $H=H_{p}+H_{\mathcal{R}}+H_{p\mathcal{R}}$, where $H_{p}$ is the particle Hamiltonian, $H_{\mathcal{R}}=H_{1}+H_{2}+H_{12}$ represents the composite reservoirs ($\mathcal{R}=\mathcal{R}_{1}\cup\mathcal{R}_{2}$) with $H_{12}$ accounting for inter-reservoir interactions, and $H_{p\mathcal{R}}$ describes particle-reservoir coupling. The formalism is initially developed using a bosonic representation, with the case of a fermionic reservoir recovered subsequently. We represent $\mathcal{R}$ as interacting sets of harmonic oscillators $\mathcal{R}_{1}=\left\{\Omega_{1j}\mid j\in J\right\}$ and $\mathcal{R}_{2}=\left\{\Omega_{2q}\mid q\in Q\right\}$ where $\Omega_{\alpha=1,2}$ is the frequency, and $J$ and $Q$ are sets of modes of $\mathcal{R}_{1}$ and $\mathcal{R}_{2}$, respectively. We describe the $\mathcal{R}_{1}$ and $\mathcal{R}_{2}$ interaction by the mode-resolved coupling $\beta_{jq}$. Eliminating $H_{12}$ (Sec. LABEL:sup:tot_Ham of Supplementary Materialsup ) results in

$$H_{\mathcal{R}}=\sum_{j}\Omega_{1j}B_{j}^{\dagger}B_{j}+\sum_{q}\Omega_{2q}b_{q}^{\dagger}b_{q}.$$

(1)

Here, $B_{j}(B_{j}^{\dagger})$ and $b_{q}(b_{q}^{\dagger})$ are the boson annihilation (creation) operators of $\mathcal{R}_{1}$ and $\mathcal{R}_{2}$, respectively. The eigenvalue equation, $H_{\mathcal{R}}|\varphi\rangle=\varepsilon_{\varphi}|\varphi\rangle$ forms the basis of our functional integral scheme. We factorize $|\varphi\rangle$ as $|\varphi\rangle=|n\rangle|m\rangle$, where $|n\rangle$ and $|m\rangle$ are the eigenstates of $\mathcal{R}_{1}$ and $\mathcal{R}_{2}$, respectively. Denoting the particle coupling with $\mathcal{R}_{1}$ and $\mathcal{R}_{2}$ as $c_{j}(x)$ and $f_{q}(x)$, respectively, with $x$ being the particle coordinate, the total Hamiltonian then reads

$$H=H_{\mathcal{R}}+\frac{p^{2}}{2m}+V(x)\ +\sum_{j}\left[c_{j}(x)B_{j}+c_{j}^{*}(x)B_{j}^{\dagger}\right]\ +\sum_{q}\left[f_{q}(x)b_{q}+f_{q}^{*}(x)b_{q}^{\dagger}\right],$$

(2)

where the second and third terms describe the particle dynamics. At the initial time $t=t_{i}$ we assume the reservoirs are in thermal equilibrium and the subsystems are noninteracting. This state is described by

$$\rho(t_{i})=\rho_{p}(t_{i})\otimes\frac{1}{Z}e^{-\beta H_{\mathcal{R}}}$$

where $\rho_{p}(t_{i})$ denotes the particle’s initial density matrix, and $Z=\sum_{\varphi}e^{-\beta\varepsilon_{\varphi}}$ is the partition function. Rapid subsystem mixing at $t>t_{i}$ justifies neglecting initial correlationsBrandbyge et al. (1995). The total density operator

$$\rho(t_{f},t_{i})=U(t_{f},t_{i})\rho(t_{i})U^{\dagger}(t_{f},t_{i})$$

governs the system’s evolution from $t_{i}$ to $t_{f}$, where $U(t_{f},t_{i})=\exp\left[-iH(t_{f}-t_{i})\right]$ and $U^{\dagger}(t_{f},t_{i})$ its corresponding adjoint. We express $\rho(t_{f})$ as a reduced density matrix defined as

$\displaystyle\rho_{red}(t_{f},t_{i})$

$\displaystyle=\ \int dx_{i}\int dy_{i}\rho_{p}(x_{i},y_{i},t_{i})\Big{<}U(x;t_{f},t_{i})U^{\dagger}(y;t_{f},t_{i})\Big{>},$

(3)

where $\langle\ldots\rangle=\frac{1}{Z}\mathrm{Tr}[e^{-\beta H_{\mathcal{R}}}\dots]$ is thermal average over the complete sets of bosonic states in Eq. (1), $\rho_{p}(x_{i},y_{i},t_{i})\equiv\langle x_{i}|\rho_{p}(t_{i})|y_{i}\rangle$, $U(x;t_{f},t_{i})\equiv\langle x_{f}|U(t_{f},t_{i})|x_{i}\rangle$ and $U^{\dagger}(y;t_{f},t_{i})\equiv\langle y_{f}|U^{\dagger}(t_{f},t_{i})|y_{i}\rangle$. $x(t)$ and $y(t)$ are the forward and backward moving paths in a contour that defines the time evolution of $U(x;t_{f},t_{i})U^{\dagger}(y;t_{f},t_{i})$. Expressing these operators in Feynman path integral representationFeynman and Hibbs (2005) gives $\rho_{red}(t_{f},t_{i})=\int dx_{i}\int dy_{i}\mathcal{J}(x_{f},y_{f},t_{f},x_{i},y_{i},t_{i})\rho_{p}(x_{i},y_{i},t_{i},)$ where the transition amplitude $\mathcal{J}(x_{f},y_{f},t_{f},x_{i},y_{i},t_{i})\equiv\mathcal{J}$ is given by

$$\mathcal{J}=\int\mathcal{D}x\int\mathcal{D}y\exp\left\{-i[S(x)-S(y)]\right\}\ \mathscr{F}_{1}[x,y]\mathscr{F}_{2}[x,y],$$

(4)

with the particle action $S(x)=\int dt\left[\frac{1}{2}m\dot{x}^{2}-V(x)\right]$. The influence functionals $\mathscr{F}_{\alpha}[x,y]$ ($\alpha=1,2$) are expressed (interaction picture) as

$\displaystyle\mathscr{F}_{\alpha}[x,y]$

$\displaystyle=\Bigg{<}\tilde{\mathcal{T}}\exp\left[i\int dtH_{I,\alpha}(x(t))\right]\mathcal{T}\exp\left[-i\int dsH_{I,\alpha}(y(s))\right]\Bigg{>},$

(5)

where, $\mathcal{T}(\tilde{\mathcal{T}})$ denotes the time (anti-time) ordering operator and

	$\displaystyle H_{I,1}(x(t))$	$\displaystyle=\sum_{j}\left[c_{j}(x(t))e^{-i\Omega_{1j}t}B_{j}+c_{j}^{*}(x(t))e^{i\Omega_{1j}t}B_{j}^{\dagger}\right]$
	$\displaystyle H_{I,2}(x(t))$	$\displaystyle=\sum_{q}\left[f_{q}(x(t))e^{-i\Omega_{2q}t}b_{q}+f_{q}^{*}(x(t))e^{i\Omega_{2q}t}b_{q}^{\dagger}\right].$

The bosonic trace is worked out in Sec. LABEL:sup:bos_trace of sup . We notice that the product of the influence functionals in Eq. (4) suggests that the reservoirs are essentially noninteractingFeynman and Hibbs (2005). We work with the influence phase $\Phi=-i\ln[\mathscr{F}]$ (Eq. (LABEL:eq:inf_phase)) instead of $\mathscr{F}$ and introduce a transformation in terms of the average $R=\frac{1}{2}(x+y)$ and relative $r=x-y$ coordinates. This is formally parallel to the generating functional methodSchmid (1982) where $R(t)$ and $r(t)$ are the $classical$ and $quantum$ fields obtained from Keldysh rotations of $x(t)$ and $y(t)$ residing on the forward and backward branches of the Keldysh contourKamenev (2011). Applying this transformation which has a Jacobian of unity modifies the position variables in Eq. (4) as $x\to R+r/2$ and $y\to R-r/2$, with their time-dependence implied. The potentials contained in $S[R,r]$ are not always conveniently defined in this transformationSchmid (1982); Newns (1985) requiring a functional Taylor expansion in $r(t)$ terminated at some order. The expansion gives

$$S(R,r)=-\int dt\left[m\ddot{R}(t)r(t)+V^{\prime}(R(t))r(t)\right].$$

(6)

where we neglected the terms $\geq\mathcal{O}(r^{3})$ that contain the quantum effects of the wavepacket spreadingNewns (1985). For harmonic $V(R)$, the higher derivatives are identically zero and $S(R,r)$ is undoubtedly quantum. For general $V(R)$, ignoring these terms constitutes the quasiclassical approximationSchmid (1982). We now split the influence phase as $\Phi_{\alpha}=\Pi_{\alpha}+i\Sigma_{\alpha}$, where $\Pi\equiv\mathrm{Re}[\Phi]$ and $\Sigma\equiv\mathrm{Im}[\Phi]$ whose functional expansions are given in Sec. LABEL:supp:fte of sup . Consequently, the transition amplitude may then be expressed as

$$\mathcal{J}=\int\mathcal{D}R\int\mathcal{D}r\exp\left\{-i\left[S(R,r)-\sum_{\alpha}\Pi_{\alpha}^{(1)}(R,r)\right]\right\}\ \exp\left[\sum_{\alpha}\Sigma^{(2)}_{\alpha}(R,r)\right],$$

(7)

where the superscript $(n)$ denotes the order of expansion. In Eq. (LABEL:eq:Sigma2) one can deduce that $\Sigma^{(2)}\propto\frac{\delta^{2}\Sigma[R,r]}{\delta r(t)\delta r(s)}$ provides the stochastic fluctuations to the dynamics via the random force autocorrelation function $K(t-s)$ defined as

$$K_{1}(t-s)=\frac{1}{2}\sum_{j}c^{\prime}_{j}(R(t))c_{j}^{\prime*}(R(s))\ \mathrm{coth}\left(\frac{\Omega_{1j}}{2k_{B}T}\right)\cos\Omega_{1j}(t-s).$$

(8)

Similarly, one can show [Eq. (LABEL:eq:Pi1)] that $\Pi_{1}^{(1)}\propto\frac{\delta\Pi[R,r]}{\delta r(t)}$ contains the dissipation through the nonlocal friction coefficient kernel

$$\Gamma_{1}(t-s)=\frac{1}{2}\sum_{j}\frac{c^{\prime}_{j}(R(t))c_{j}^{\prime*}(R(s))}{\Omega_{1j}}\ [\mathrm{sgn}(t-s)+1]\cos\Omega_{1j}(t-s).$$

(9)

and a potential shift $\Delta\Omega_{1}\equiv\sum_{j}\frac{|c_{j}(R)|^{2}}{2\Omega_{1j}}$ that renormalizes $V(R)$. Except for the presence of a sign function, the form of $\Gamma(t-s)$ above closely resembles the friction kernels previously obtained by Head-Gordon and TullyHead‐Gordon and Tully (1995). Expressions for $K_{2}(t-s)$ and $\Gamma_{2}(t-s)$ are retrieved by replacing in Eq. (9) and Eq. (8) their corresponding variables. The transition amplitude now takes the form

	$\displaystyle\mathcal{J}$	$\displaystyle=\int\mathcal{D}R\int\mathcal{D}r\exp\left\{i\int dt\left[m\ddot{R}(t)+\tilde{V}^{\prime}(R(t))\ +\int ds\sum_{\alpha}\Gamma_{\alpha}(t-s)\dot{R}(s)\right]r(t)\right.$		(10)
		$\displaystyle\left.-\frac{1}{2}\int dt\int dsr(t)\sum_{\alpha}K_{\alpha}(t-s)r(s)\right\},$

where $\tilde{V}(R)=V(R)-\sum_{\alpha}\Delta\Omega_{\alpha}$. The summations over $\alpha$ in Eq. (10) implies the existence of an effective random force $\xi(t)$ that accounts for the stochastic fluctuations coming from both reservoirs. $\xi(t)$ has the properties $\langle\xi(t)\rangle=0$ and $\langle\xi(t)\xi(s)\rangle=\sum_{\alpha}K_{\alpha}(t-s)$. After stochastic averaging of Eq. (10) over $\xi(t)$, the $r(t)$ path integral can be done explicitly (Sec.LABEL:supp:stoav of sup ). This yields a Dirac delta functional (Eq. (LABEL:eq:supp_J)) where we deduce the generalized quasiclassical Langevin equation

$\displaystyle m\ddot{R}(t)+\tilde{V}^{\prime}(R(t))\ +\int ds\sum_{\alpha}\Gamma_{\alpha}(t-s)\dot{R}(s)=\xi(t)$

(11)

describing the motion of the particle interacting simultaneously with two thermal reservoirs.

Let us now turn to the case where one of the reservoirs is fermionic. Specifically, we consider a system consisting of a particle interacting with a metal surface immersed in a solvent. The solvent modes and the electronic manifold of the metal represent the bosonic ($\mathcal{R}_{1}$) and fermionic ($\mathcal{R}_{2}$) reservoirs, respectively. The frequency and eigenstates are renormalized by the inter-reservoir interaction, but for simplicity, we assume that the simulation parameters have been already renormalized. We also assume that the particles couple linearly with the bosonic bath, thereby yielding the featureless bosonic friction coefficient $\gamma$ in the Caldeira-Leggett frameworkCaldeira and Leggett (1983b, a) (Eq. (LABEL:eq:fric_kerR1_2) of sup ). The dynamics of solvated particles on a metal is typically described by the time-dependent Newns-Anderson-Schmickler (TD-NAS) modelSchmickler (1986); Arguelles and Sugino (2024). In this model, the electronic coupling is treated in the wide-band limit in terms of $\Delta(t)=\pi\sum_{k}|V_{ak}(t)|^{2}\delta(\varepsilon-\varepsilon_{k})$, with $V_{ak}(t)$ being the time-dependent hopping integrals and $\varepsilon_{k}$ is the band energy. The particle level is renormalized as $\tilde{\varepsilon}_{a}(t)=\varepsilon_{a}(t)+(2Z-1)\lambda$ by the solvent reorganization energy $\lambda$ ($Z$ is the particle charge). We note that the TD-NAS model reduces to the usual time-dependent Anderson-Holstein (AH) HamiltonianAnderson (1961); Holstein (1959) for $Z=0$ when $\lambda$ is re-interpreted as an average electron-phonon coupling. Using the TD-NAS model, we obtained an effective HamiltonianBrako and Newns (1981, 1989); Arguelles and Sugino (2024) for the metal electron system perturbed by the slowly approaching particle. We then mappedTomonaga (1950); Brako and Newns (1989); Mahan (2000) it to $\mathcal{R}_{2}$ to obtain $f_{q}$ in Eq. (2). This effective Hamiltonian is given by Eq. (LABEL:eq:Heh1) in sup , where (after bosonization) we identified $f_{q}(R)=\left[\frac{2\omega_{q}}{\pi^{2}\rho(\varepsilon_{F})}\right]^{1/2}\int dt\dot{\delta}(\varepsilon_{F},R)\exp(i\omega_{q}t)$. Inserting $f_{q}(R(t))$ into the $\mathcal{R}_{2}$ version of Eq. (9), yields the electronic frictional force $\int ds\Gamma_{2}(t-s)\dot{R}(s)=\eta(R)\dot{R}(t)$ where (Eq. (LABEL:eq:fric_kerR2_2) of sup )

$$\eta(R)=\pi\left\{\frac{d\Delta(R)}{dR}\frac{\left[\varepsilon_{F}-\tilde{\varepsilon}_{a}(R)\right]}{\Delta(R)}\ +\frac{d\tilde{\varepsilon}_{a}(R)}{dR}\right\}^{2}\rho_{a}^{2}(\varepsilon_{F},R),$$

(12)

is the electronic friction coefficient. Here, $\rho_{a}(\varepsilon_{F},R)=\frac{1}{\pi}\frac{\Delta(R)}{[\varepsilon_{F}-\tilde{\varepsilon}_{a}(R)]^{2}+\Delta^{2}(R)}$ is the particle local density of states. Eq. (12) agrees with different forms of $\eta(R)$ that have been derived over the years with varying levels of generalizationsNourtier, A. (1977); Makoshi (1983); Brako and Newns (1981); Schonhammer and Gunnarsson (1980); Head‐Gordon and Tully (1995); Plihal and Langreth (1998); Trail et al. (2003); Mizielinski et al. (2005); Dou et al. (2017). Eq. (11) then reduces to

$\displaystyle m\ddot{R}(t)+\tilde{V}^{\prime}(R(t))\ +\gamma_{eff}(R)\dot{R}(t)=\xi(t),$

(13)

where $\gamma_{eff}(R)=\gamma+\eta(R)$ and $\xi(t)=\xi_{\gamma}(t)+\xi_{\eta}(t)$ is the effective random force with a corresponding autocorrelation function $K(t-s)=K_{\gamma}(t-s)+K_{\eta}(t-s)$, where

	$\displaystyle K_{\gamma}(t-s)$	$\displaystyle=\gamma\int_{0}^{\infty}d\Omega\Omega\ \coth\left(\frac{\Omega}{2k_{B}T}\right)\cos\Omega(t-s)$		(14)
	$\displaystyle K_{\eta}(t-s)$	$\displaystyle=\eta(R)\int_{0}^{\infty}d\omega\omega\ \coth\left(\frac{\omega}{2k_{B}T}\right)\cos\omega(t-s).$

Using the influence function methodFeynman and Vernon (2000), we have derived a non-phenomenological equation of motion for a particle accounting for simultaneous interactions with bosonic and fermionic thermal reservoirs. As an example of the fermionic reservoir, we take electrons in a metal electrode, which dissipate the particle energy into the creation of electron-hole pairs thereby causing electronic friction (EF). An explicit expression for the local dissipation kernel (Markovian kernel) is given in the limit of slow particle motion providing thereby a way to compute the rate of energy transfer via a stochastic simulation. For demonstration purposes, we applied the framework to prototypical electrochemical systems where a hydrogen (H) atom moves in contact with a metal electrode and a solvent and investigated the interplay of the reservoirs using two setups: (1) quantum vibrational relaxation of a hydrogen (H) confined on a metal surface and (2) solvated electrochemical proton discharge. In the former case we have seen that the electronic friction enhances the vibrational friction, and in addition, quantum fluctuation effect manifests as a spread of the wave function of H within the Gaussian white noise approximationFurutani and Salasnich (2023). In the latter case, we have seen that, contrary to the phononic friction, the electronic friction is active at a short time interval of electronic level crossing and contributes to delaying of the electron and proton transfers. On the other hand, the average energy dissipation is comparable between the two reservoirs despite the electronic friction being local near the level crossing. This work marks the first explicit consideration of friction from b oth solvent modes and electronic excitations in the metal during electrochemical proton discharge. Although initially motivated by electrochemistry, the presented framework is be broadly applicable to systems where interactions with multiple thermal baths play a significant role.

This research is supported by the New Energy and Industrial Technology Development Organization (NEDO) project, MEXT as “Program for Promoting Researches on the Supercomputer Fugaku” (Fugaku battery & Fuel Cell Project) (Grant No. JPMXP1020200301, Project No.: hp220177, hp210173, hp200131), Digital Transformation Initiative for Green Energy Materials (DX-GEM) and JSPS Grants-in-Aid for Scientific Research (Young Scientists) No. 19K15397. Some calculations were done using the supercomputing facilities of the Institute for Solid State Physics, The University of Tokyo.