# Nonequilibrium reaction rate theory: Formulation and implementation within
the hierarchical equations of motion approach

Yaling Ke<EMAIL_ADDRESS>Institute of Physics, Albert-
Ludwig University Freiburg, Hermann-Herder-Strasse 3, 79104 Freiburg, Germany
Christoph Kaspar  Institute of Physics, Albert-Ludwig University Freiburg,
Hermann-Herder-Strasse 3, 79104 Freiburg, Germany  André Erpenbeck  Department
of Physics, University of Michigan, Ann Arbor, Michigan 48109, USA  Uri Peskin
Schulich Faculty of Chemistry, Technion-Israel Institute of Technology, Haifa
32000, Israel  Michael Thoss<EMAIL_ADDRESS>Institute
of Physics, Albert-Ludwig University Freiburg, Hermann-Herder-Strasse 3, 79104
Freiburg, Germany  EUCOR Centre for Quantum Science and Quantum Computing,
Albert-Ludwig University Freiburg, Hermann-Herder-Strasse 3, 79104 Freiburg,
Germany

###### Abstract

The study of chemical reactions in environments under nonequilibrium
conditions has been of interest recently in a variety of contexts, including
current-induced reactions in molecular junctions and scanning tunneling
microscopy experiments. In this work, we outline a fully quantum mechanical,
numerically exact approach to describe chemical reaction rates in such
nonequilibrium situations. The approach is based on an extension of the flux
correlation function formalism to nonequilibrium conditions and uses a mixed
real and imaginary time hierarchical equations of motion approach for the
calculation of rate constants. As a specific example, we investigate current-
induced intramolecular proton transfer reactions in a molecular junction for
different applied bias voltages and molecule-lead coupling strengths.

##  I Introduction

Chemical reactions are often characterized by thermal rate constants. The
theory of chemical reaction rates is well developed for reactions in the gas
phase or in condensed phases, including a variety of approaches that range
from transition state theories or approaches that use classical dynamics to
numerically exact fully quantum mechanical methods.wyatt1996dynamics ;
RevModPhys.62.251 ; nitzan2006chemical ; miller1975semiclassical ;
miller1975path ; chandler1978statistical ; truhlar1983current ;
Topaler_1994_J.Chem.Phys._p7500 ; geva2001quantum ; pollak2005reaction ;
manthe1995new ; manthe2008state ; wang2000forward ;
Richardson_2014_J.Chem.Phys._p74106 A particular efficient approach to
calculating thermal reaction rates is the flux correlation function
formalism.RevModPhys.62.251 ; Miller_J.Chem.Phys._1983_p4889–4898 ;
Thompson_J.Chem.Phys._1997_p142–150 ; manthe2002reaction ;
Wang_J.Chem.Phys._2006_p174502 ; Craig_J.Chem.Phys._2007_p144503 A basic
assumption of these approaches is that the reaction occurs in equilibrium,
i.e. the reactants are in equilibrium at a given temperature.

This assumption is not fulfilled for chemical reactions under nonequilibrium
conditions, which have been of interest in different contexts in recent years.
Examples include chemical reactions induced in scanning tunneling microscopy
(STM) experiments at molecule-metal interfaces PhysRevLett.100.166101 ;
PhysRevB.79.035423 ; liljeroth2007current ; PhysRevB.81.045402 ; peller2020sub
; ladenthin2015hot ; kumagai2014controlling ; boeckmann2016direct as well as
current-induced chemical reactions in single-molecule
junctions.koch2006current ; Aradhya_2013_Nat.Nanotechnol._p399 ;
Zang_2020_NanoLett._p673 ; huang2019electric ; PhysRevLett.123.246601 ;
avriller2021photon ; chen2021electron To study chemical reactions in such
scenarios, we have recently proposed a fully quantum mechanical simulation
scheme, based on the hierarchical equations of motion (HEOM) approach.
Erpenbeck_2019_J.Chem.Phys._p191101 ; Erpenbeck_2020_Phys.Rev.B_p195421 ;
ke2021unraveling Applications to current-induced bond rupture at metal
surfaces, a process which is closely related to the breakdown of molecular
functionalities,Sabater_2015_BeilsteinJ.Nanotechnol._p2338 ;
Li_2015_J.Am.Chem.Soc._p5028 ; Li_2016_J.Am.Chem.Soc._p16159 ;
Capozzi_NanoLett._2016_p3949–3954 showed that the HEOM method provides an
accurate and efficient description of this process at high bias
voltages.Erpenbeck_2020_Phys.Rev.B_p195421 ; ke2021unraveling However, in the
low bias voltage regime, reactions occur on a much longer timescale. Instead
of tracking the dynamics of the slow transformation process from the reactants
to the products, in such cases, the formulation of a reaction rate theory and
the development of an efficient approach to calculate such rates are highly
desirable.

To address this issue, in this work, we extend the flux correlation function
formalism for the calculation of reaction rates to nonequilibrium scenarios.
Within this extended formalism, the evaluation of an extended “Boltzmannized”
flux operator is required, which involves not only the molecular degrees of
freedom (DoFs) but also the manifold of all electronic states of the metal
surface/electrodes. To obtain this operator and calculate reaction rate in an
efficient and exact way, we further present a mixed real and imaginary time
HEOM approach for open quantum systems coupled to fermionic reservoirs. It is
worth mentioning that the imaginary time HEOM approach can also be used as an
effective tool for evaluating thermal equilibrium properties as well as
preparing a correlated initial state for the real time HEOM method.

As a specific example to illustrate the validity and applicability of the
proposed nonequilibrium flux correlation function formalism using the mixed
real and imaginary time HEOM method, we consider intramolecular proton
transfer reactions in a single-molecule junction. Proton transfer has been
found to be involved in a variety of chemical and biological processes,
Hammes-Schiffer_J.Phys.Chem.B_2008_p14108 and in this case quantum effects
play an important role, which limits the applicability of approaches that are
based on a classical description of the reactive nuclear
DoFs.Tully_J.Chem.Phys._1971_p562 ; Shenvi_J.Chem.Phys._2009_p174107 ;
Fischer_J.Chem.Phys._2011_p144102 ; Head_J.Chem.Phys._1995_p10137 ;
Dou_J.Chem.Phys._2016_p24116 ; Galperin_J.Phys.Chem.Lett_2015_p4898 ;
Lu_Prog.Surf.Sci._2019_p21 ; Lu_Phys.Rev.B_2012_p245444 ;
Dzhioev_2011_J.Chem.Phys._p74701 ; Dzhioev_2013_J.Chem.Phys._p134103 ;
Pozner_2014_NanoLett._p6244 ; Preston_Phys.Rev.B_2020_p155415 ;
Erpenbeck_2018_Phys.Rev.B_p235452 In the context of molecular electronics, it
has been shown that current-induced proton transfer processes can be used to
realize molecular nanoswitches and molecular
transistors.Liljeroth_Science_2007_p1203–1206 ;
Sobolewski_Phys.Chem.Chem.Phys._2008_p1243–1247 ;
Benesch_J.Phys.Chem.C_2009_p10315–10318 ; Hofmeister_J.Mol.Model._2014_p1–5 ;
Weckbecker_J.Phys.Chem.Lett._2020_p413–417 To control and design such devices
requires a thorough theoretical study of current-induced proton transfer
processes at molecule-metal interfaces, which was so far limited to
perturbative master equations.Hofmeister_2017_J.Chem.Phys._p92317

The remainder of this work is arranged as follows: In Sec. II, we introduce
the nonequilibrium flux correlation function formalism for the specific
example of current-induced reactions in molecular junctions. Furthermore, the
description of the quantum dynamics with a correlated initial state and the
real and imaginary time HEOM approach are outlined. Applications of the method
to simulate current-induced proton transfer in molecular junctions are
presented in Sec. III. In particular, we investigate the influence of the
applied bias voltage and the molecule-lead coupling on the proton transfer
rate. Sec. IV concludes with a summary.

##  II Model and Method

###  II.1 Model

As a specific example of a nonequilibrium scenario, we consider a molecular
junction, where a single molecule is bounded to macroscopic electrodes. The
full system is modeled by the Hamiltonian,

$H=H_{s}+H_{b}+H_{c},$ (1)

where $H_{s}$ corresponds to the molecular system, $H_{b}$ to electrodes, and
$H_{c}$ describes their coupling.

The molecule comprises electrons and nuclei and the corresponding Hamiltonian
is given by

$H_{s}=T_{n}+H_{el},$ (2)

where $T_{n}=\sum_{i}p_{i}^{2}/2m_{i}$ represents the kinetic energy operator
of the nuclei. The variables $p_{i}$, $x_{i}$, and $m_{i}$ denote the
momentum, position, and the effective mass of the $i$-th nuclear DoF,
respectively. $H_{el}$ is the electronic Hamiltonian which includes electronic
kinetic energy and the interactions among electrons and nuclei.

The electrodes serve as electron reservoirs and are represented as a continuum
of non-interacting electrons,

$H_{b}=\sum_{\alpha}\sum_{k}\epsilon_{\alpha k}c_{\alpha k}^{\dagger}c_{\alpha
k},$ (3)

where $c_{\alpha k}^{\dagger}$ and $c_{\alpha k}$ are the creation and
annihilation operators for an electron in the $k$-th state of electrode
$\alpha$ with the corresponding energy

$\epsilon_{\alpha k}=\epsilon_{\alpha k}^{(0)}+\mu_{\alpha}.$ (4)

Here, $\mu_{\alpha}$ denotes the electrochemical potential of electrode
$\alpha$, which depends on the applied bias voltage $\Phi$, and
$\epsilon_{\alpha k}^{(0)}$ is the corresponding energy at zero bias voltage.

The coupling between the molecule and electrodes, which generally depends on
the nuclear position $\mathbf{x}$, leads to electron transport. We assume that
the interaction Hamiltonian takes the general form

$H_{c}=\sum_{\alpha k}v_{\alpha k}s(\mathbf{x})c_{\alpha
k}^{\dagger}d+v^{*}_{\alpha k}s(\mathbf{x})d^{\dagger}c_{\alpha k},$ (5)

where $v_{\alpha k}$ specifies the maximal coupling strength between the
molecular electronic level and the $k$-th state in electrode $\alpha$. Here, a
dimensionless real function $s(\mathbf{x})$ is introduced to model the
dependence of the coupling on the nuclear position. $d^{\dagger}$ $(d)$
denotes the operator for creating (annihilating) an electron in a specific
molecular electronic level. Note that for notational simplicity, we consider a
system with a single electronic level, but the extension to a multilevel
system is straightforward.

###  II.2 Nonequilibrium flux correlation function formalism

The connection between the equilibrium reaction rate and the flux correlation
function is well established.Miller_J.Chem.Phys._1983_p4889–4898 ;
Thompson_J.Chem.Phys._1997_p142–150 ; Wang_J.Chem.Phys._2006_p174502 ;
Craig_J.Chem.Phys._2007_p144503 However, reaction dynamics for molecular
systems in contact with macroscopic electrodes is generally subject to
nonequilibrium situations. For instance, the electrodes are usually held at
different chemical potentials. As a consequence, the molecular system is
driven out of equilibrium by the electron current. Our goal is to extend the
equilibrium flux correlation function formalism and formulate a nonequilibrium
quantum reaction rate theory. In the present work, we focus on a non-
equilibrium scenario where the molecule is bound to two electrodes under an
applied bias voltage $e\Phi=\mu_{L}-\mu_{R}$, and both leads have the same
temperature $T$ (or the inverse temperature $\beta=1/k_{\rm B}T$, with $k_{\rm
B}$ being the Boltzmann constant).

First, we briefly review the basic ideas of the flux correlation function
formalism pertaining to its application in this work. In the study of a
unimolecular chemical reaction, the system is usually divided into a reactant
and a product region, which are separated by a dividing surface located at
$x_{\rm ds}$. The quantum-mechanical projection operator is given by

$h=\theta(x-x_{\rm ds}),$ (6)

where $\theta(x)$ denotes the Heaviside function, and its Heisenberg time-
derivative yields the flux operator

$F=\frac{i}{\hbar}[H,h].$ (7)

We assume that the system is initially located in the reactant region and
represented by a non-stationary density operator $\rho_{r}$. Thereafter, the
time evolution of the reactant and product population are given by

$P_{r}(t)=\textrm{tr}\\{(1-h)\rho(t)\\}=\textrm{tr}\\{(1-h)e^{-\frac{i}{\hbar}Ht}\rho_{r}e^{\frac{i}{\hbar}Ht}\\}$
(8)

and

$P_{p}(t)=1-P_{r}(t)=\textrm{tr}\\{h\rho(t)\\}=\textrm{tr}\\{he^{-\frac{i}{\hbar}Ht}\rho_{r}e^{\frac{i}{\hbar}Ht}\\},$
(9)

where $\rho(t)=e^{-\frac{i}{\hbar}Ht}\rho_{r}e^{\frac{i}{\hbar}Ht}$ denotes
the density operator at time $t$.

The definition of a rate constant for the transition from reactant to product
(and vice versa) is based on the assumption that the population transfer is
governed by first order kinetics at sufficiently long times,

$\frac{d}{dt}P_{r}(t)=-\frac{d}{dt}P_{p}(t)=-kP_{r}(t)+k_{\rm b}P_{p}(t),$
(10)

in which $k$ and $k_{\rm b}$ are the forward and backward reaction rate
constant, respectively. The forward rate constant $k$ can be formally
expressed asCraig_J.Chem.Phys._2007_p144503

$\displaystyle k$ $\displaystyle\equiv$
$\displaystyle\lim_{t\rightarrow\infty}k(t)$ $\displaystyle=$
$\displaystyle\lim_{t\rightarrow\infty}\frac{C_{\rm
f}(t)}{P_{r}(0)+\chi[P_{r}(0)-1]-(1+\chi)\int_{0}^{t}C_{\rm f}(\tau)d\tau}.$

Here, $P_{r}(0)=\textrm{tr}\left\\{(1-h)\rho_{r}\right\\}$ denotes the initial
reactant population, and $\chi=k_{\rm b}/k$ is the ratio of the backward and
forward reaction rate, which depends on the temperature $\beta$ and applied
bias voltage $\Phi$. It can be obtained as the ratio of the steady state
solutions of $P_{r}(t)$ and $P_{p}(t)$,

$\chi(\beta,\Phi)=\frac{P^{\rm ss}_{r}}{P^{\rm ss}_{p}},$ (12)

since the time derivative of the reactant (or product) population, as given in
Eq. (10), is zero at steady state. Besides, we have introduced the reactive
flux correlation function

$C_{\rm
f}(t)=-\dot{P}_{r}(t)=\textrm{tr}\left\\{F_{\beta}e^{\frac{i}{\hbar}Ht}he^{-\frac{i}{\hbar}Ht}\right\\},$
(13)

and the “Boltzmannized” flux operatorMiller_J.Chem.Phys._1983_p4889–4898

$F_{\beta}=-\frac{i}{\hbar}[H,\rho_{r}].$ (14)

An alternative approach for calculating the rate constant is based on an
assumption that the population transfer between reactants and products is much
slower than any other dynamical processes in the system and its surroundings.
When this assumption holds, the forward rate, e.g., can be calculated directly
at “shorter times”,

$k\approx\lim_{t\rightarrow t_{\rm plateau}}C_{\rm f}(t),$ (15)

where $t_{\rm plateau}$ is the time at which $C_{\rm f}(t)$ reaches a plateau.
The above approximation is reliable when two conditions are met. The first
condition is that the initial reactant population $P_{r}(0)$ is equal or very
close to 1. Second, the reaction timescale is much longer than the timescales
of other transient dynamical processes. When these two conditions are
fulfilled, the loss of reactant population is negligible and the denominator
in Eq. (II.2) is essentially unity over the timescale considered. As such, one
can avoid the calculation of $\chi$ and obtain the reaction rate by computing
the flux correlation function using quantum dynamics approaches.

Notice that the selected form of $\rho_{r}$ determines the explicit expression
of the Boltzmannized flux operator in Eq. (14) and influences significantly
the short-time dynamics of $C_{\rm f}(t)$. Nevertheless, different choices of
$\rho_{r}$ which represent reliably the initial thermal ensemble (either
”reactants” or ”products”) are expected to yield similar values for the rate
constants. An optimal choice of $\rho_{r}$ should therefore be both a reliable
representation of the initial ensemble and facilitate rapid convergence of
$C_{\rm f}(t)$ to a plateau. As illustrated in the supporting material, in the
equilibrium cases ($\mu_{L}=\mu_{R}=0$), the flux correlation function shows
the least oscillating behavior with the choice

$\rho^{\rm eq}_{r}=\frac{1}{Q_{r}^{\rm eq}}e^{-\beta H/2}(1-h)e^{-\beta H/2},$
(16)

corresponding to the Boltzmannized flux operator

$F^{\rm eq}_{\beta}=\frac{1}{Q_{r}^{\rm eq}}e^{-\beta H/2}Fe^{-\beta H/2}.$
(17)

In nonequilibrium cases, a related choice is

$\rho_{r}^{\rm ss}=\frac{1}{{Q^{\rm
ss}_{r}}}e^{-\beta(H-Y)/2}(1-h)e^{-\beta(H-Y)/2}.$ (18)

Here, $Q_{r}^{\rm ss}=\mathrm{tr}\\{(1-h)e^{-\beta(H-Y)}\\}$, and the operator
$Y$ is given formally as

$Y=\sum_{\alpha=L/R}\Omega\mu_{\alpha}N_{\alpha}\Omega^{-1},$ (19)

where
$\Omega=\lim_{t\rightarrow\infty}e^{\frac{i}{\hbar}Ht}e^{-\frac{i}{\hbar}{H_{b}t}}$
is the Møller operator, and $N_{\alpha}=\sum_{k}c_{\alpha
k}^{\dagger}c_{\alpha k}$ represents the occupation number operator of the
electrons in lead $\alpha$. The above expressions are built upon the
nonequilibrium steady state (NESS) theory,Ness_Phys.Rev.E_2014_p_62119 ;
Ness_Entropy_2017_p158 which shows that the nonequilibrium steady state can be
described by a density matrix in a Gibbs-like form,

$\rho^{\rm ss}=e^{-\beta(H-Y)}/\mathrm{tr}\\{e^{-\beta(H-Y)}\\}.$ (20)

Using the fact that the operator $e^{-\beta(H-Y)}$ commutes with the
Hamiltonian $H$, the Boltzmannized flux operator is expressed as

$F^{\rm ss}_{\beta}=e^{-\beta(H-Y)/2}Fe^{-\beta(H-Y)/2}/Q_{r}^{\rm ss},$ (21)

which reduces for vanishing bias voltage ($\mu_{L}=\mu_{R}=0$) to the
equilibrium expression Eq. (17). Expression Eq. (21) for the flux operator,
together with Eqs. (II.2) - (13), generalizes the flux correlation function
formalism to nonequilibrium situations and represent the main theoretical
result of this work.

Although Eqs. (18) - (21) have an appealing structure, it is generally
difficult to obtain the exact expression of $Y$. The expression simplifies
significantly by replacing $H-Y$ by the zero-bias Hamiltonian
$H-\sum_{\alpha}\mu_{\alpha}N_{\alpha}$, which corresponds to

$\rho_{r}^{0}=\frac{1}{Q^{0}_{r}}e^{-\frac{\beta}{2}(H-\sum_{\alpha}\mu_{\alpha}N_{\alpha})}(1-h)e^{-\frac{\beta}{2}(H-\sum_{\alpha}\mu_{\alpha}N_{\alpha})}$
(22)

with the partition function of the reactant

$Q^{0}_{r}=\mathrm{tr}\left\\{(1-h)e^{-\beta(H-\sum_{\alpha}\mu_{\alpha}N_{\alpha})}\right\\}.$
(23)

The expression of $\rho_{r}^{0}$ in Eq. (22) reduces also to the equilibrium
expression, Eq. (16) for vanishing bias voltage. Although it does not lead to
an intuitive form of a Boltzmannized flux operator, it is significantly easier
to calculate than Eq. (18). Furthermore, from the physics perspective, notice
that this choice of $\rho_{r}$ implies that the rate constant, as defined by
Eq. (II.2) or (15), corresponds to the reaction induced by a sudden switch on
of the bias voltage. Therefore, this rate should be directly comparable to
experiments in which the reaction is induced by the applied voltage, and the
reaction rate is much slower than the rate of the voltage ramp. Moreover, this
approach applies also to voltage-induced dissociation rate
calculationsfung2019breaking ; ke2021unraveling by Eq. (II.2) or (15).

By performing the cyclic permutation within the trace expression, the reactive
flux correlation function in Eq. (13) is reformulated into the form

$C_{\rm
f}(t)=\mathrm{tr}\left\\{Fe^{-\frac{i}{\hbar}Ht}\rho_{r}e^{\frac{i}{\hbar}Ht}\right\\}=\mathrm{tr}_{s}\left\\{F\rho_{r}^{s}(t)\right\\},$
(24)

which is valid for both choices, Eq. (18) and (22). Here, we have used fact
that the flux operator $F$ acts only on the system subspace, because the
coupling Hamiltonian $H_{c}$ is only a function of $\mathbf{x}$ and commutes
with the projection operator $h$. The above equation indicates that the
reactive flux correlation function $C_{\rm f}(t)$ can be calculated as the
expectation value of the flux operator with respect to the reduced density
operator

$\rho_{r}^{s}(t)=\mathrm{tr}_{b}\\{e^{-\frac{i}{\hbar}Ht}\rho_{r}e^{\frac{i}{\hbar}Ht}\\}.$
(25)

Above, $\mathrm{tr}_{s}$ and $\mathrm{tr}_{b}$ denote tracing over system and
environmental DoFs, respectively.

The numerically exact calculations of $\rho_{r}^{s}(t)$ can be performed with
the HEOM approach. However, the method typically starts with a factorized
initial state between the system and environment parts, which is not the case
for $\rho_{r}$ in Eq. (22). Therefore, an additional imaginary time
propagation is required. Below we present a mixed real and imaginary time HEOM
approach that can be used to efficiently obtain nonequilibrium flux
correlation functions for chemical reaction dynamics in molecular systems that
are in contact with fermionic reservoirs.

###  II.3 Hierarchical equations of motion

In the seminal work by Tanimura and Kubo,Tanimura_1989_J.Phys.Soc.Jpn._p101
the HEOM approach was proposed to investigate relaxation dynamics. Over the
last three decades, considerable progress has been made to extend its
applicability.Tanimura_2020_J.Chem.Phys._p20901 ;
Yan_Chem.Phys.Lett._2004_p216–221 ; Xu_2005_J.Chem.Phys._p41103 ;
Ishizaki_2005_J.Phys.Soc.Jpn._p3131 ; Shi_2009_J.Chem.Phys._p84105 ;
Hu_2010_J.Chem.Phys._p101106 ; Kreisbeck_J.Chem.TheoryComput._2011_p2166–2174
; Suess_Phys.Rev.Lett._2014_p150403 ; Hou_2015_J.Chem.Phys._p104112 ;
Tang_2015_J.Chem.Phys._p224112 ;
Tsuchimoto_J.Chem.TheoryComput._2015_p3859–3865 ; Ke_J.Chem.Phys._2016_p024101
; Hsieh_2018_J.Chem.Phys._p14103 ; Nakamura_Phys.Rev.A_2018_p012109 ;
Shi_J.Chem.Phys._2018_p174102 ; Erpenbeck_2018_J.Chem.Phys._p64106 ;
Zhang_2020_J.Chem.Phys._p64107 ; Ikeda_J.Chem.Phys._2020_p204101 Tanimura
developed an imaginary time HEOM for open systems coupled to bosonic
environments to facilitate the calculations of thermal
quantities.Tanimura_2014_J.Chem.Phys._p44114 Yan and coworkers significantly
extend the method for studying charge transport in quantum impurity
models.Jin_2008_J.Chem.Phys._p234703 ; Zheng_2008_J.Chem.Phys._p184112 ;
Zheng_2013_Phys.Rev.Lett._p86601 ; Yan_2014_J.Chem.Phys._p54105 ;
Ye_2016_WIREsComputMolSci_p608 ; Zhang_2017_J.Chem.Phys._p44105 ;
Han_2018_J.Chem.Phys._p234108 Schinabeck et al. further extended it to
investigate charge transport in nanosystems with electronic-vibrational
couplings.Schinabeck_2016_Phys.Rev.B_p201407 ;
Schinabeck_2018_Phys.Rev.B_p235429 ;
Erpenbeck_2019_TheEuropeanPhysicalJournalSpecialTopics_p1981 ;
Schinabeck_Phys.Rev.B_2020_p075422 Erpenbeck et al. formulated and applied the
method for studying reaction dynamics at metal surfaces using real time
propagation.Erpenbeck_2019_J.Chem.Phys._p191101 Here, we present a mixed real
and imaginary time HEOM approach for open quantum systems which are initially
correlated with fermionic environments. The method is formulated within the
reduced system framework, where the reduced density operator $\rho_{r}^{s}(t)$
is obtained by tracing over all reservoir DoFs. For the Hamiltonian given in
Sec. II.1, this can be performed analytically in the Feynman path-integral
formalism,Feynman__2010_p ; Zinn-Justin__2010_p and the expression of
$\rho_{r}^{s}(t)$ reads

$\displaystyle\rho^{s}_{r}(\varphi^{+}_{t},\tilde{\varphi}^{-}_{t},t)$
$\displaystyle=$
$\displaystyle\langle\varphi_{t}|\mathrm{tr}_{b}\left\\{e^{-\frac{i}{\hbar}Ht}\hat{\rho}_{r}e^{\frac{i}{\hbar}Ht}\right\\}|\tilde{\varphi}_{t}\rangle$
$\displaystyle=$
$\displaystyle\frac{1}{Q_{r}}\int\mathrm{d}\varphi^{+}_{0}\mathrm{d}\varphi^{-}_{0}\int\mathrm{d}\tilde{\varphi}^{+}_{0}\mathrm{d}\tilde{\varphi}^{-}_{0}\int_{\varphi^{-}_{0}}^{\varphi^{+}_{t}}\mathcal{D}\varphi^{+}(t)\mathcal{D}\bm{\varphi}^{-}(t)\int_{\varphi^{+}_{0}}^{\tilde{\varphi}^{-}_{0}}\mathcal{D}\bm{\phi}^{+}(\tau)\mathcal{D}\bm{\phi}^{-}(\tau)\int_{\tilde{\varphi}^{+}_{0}}^{\tilde{\varphi}^{-}_{t}}\mathcal{D}\tilde{\bm{\varphi}}^{+}(t)\mathcal{D}\tilde{\bm{\varphi}}^{-}(t)$
$\displaystyle(1-h(\phi^{+}(\beta\hbar/2),\phi^{-}(\beta\hbar/2)))e^{\frac{i}{\hbar}S(\bm{\varphi}^{+}(t),\bm{\varphi}^{-}(t))-\frac{1}{\hbar}S_{\beta}(\bm{\phi}^{+}(\tau),\bm{\phi}^{-}(\tau))-\frac{i}{\hbar}S^{*}(\tilde{\bm{\varphi}}^{+}(t),\tilde{\bm{\varphi}}(t))}\mathcal{F}.$

Here, we have introduced the fermionic coherent state
$\left|\varphi\right\rangle$, which is the eigenstate of $d$ with eigenvalue
$\varphi^{-}$, i.e.,
$d\left|\varphi\right\rangle=\varphi^{-}\left|\varphi\right\rangle$. Its
conjugation yields
$\left\langle\varphi\right|d^{\dagger}=\left\langle\varphi\right|\varphi^{+}$
for the creation operator $d^{\dagger}$. $\varphi^{+}$ and $\varphi^{-}$ are
independent Grassmann variables and satisfy the anticommutation relations
$\\{\varphi^{\sigma},\varphi^{\bar{\sigma}}\\}=1$ and
$\\{\varphi^{\sigma},\varphi^{\sigma}\\}=0$. For conciseness, we introduce the
symbol $\sigma=\pm$ and $\bar{\sigma}=-\sigma$. $\bm{\varphi}^{\sigma}(t)$ and
$\tilde{\bm{\varphi}}^{\sigma}(t)$ are the forward and backward system paths
along the real time axis $t$. $\bm{\phi}^{\sigma}(\tau)$ denote the paths in
the imaginary time $\tau$ within the range of $0$ and $\beta\hbar$. $S$ and
$S_{\beta}$ denote the system action functional and their explicit expressions
can be found in the supplementary material.

All the statistical information about the electron reservoirs and their
influence on the system dynamics are encoded in the influence functional

$\begin{split}\mathcal{F}=&e^{-\frac{1}{\hbar}\int_{0}^{t}\mathrm{d}s\int_{0}^{s}\mathrm{d}s^{\prime}\sum_{\sigma=\pm}\left[\varphi^{\bar{\sigma}}(s)+\tilde{\varphi}^{\bar{\sigma}}(s)\right]\sum_{\alpha}\left[G_{\alpha}^{\sigma}(s-s^{\prime})\varphi^{\sigma}(s^{\prime})-G^{\bar{\sigma}*}_{\alpha}(s-s^{\prime})\tilde{\varphi}^{\sigma}(s^{\prime})\right]}\\\
&e^{\frac{i}{\hbar^{2}}\int_{0}^{t}\mathrm{d}s\sum_{\sigma=\pm}\left[\varphi^{\bar{\sigma}}(s)+\tilde{\varphi}^{\bar{\sigma}}(s)\right]\int_{0}^{\beta\hbar}\mathrm{d}\tau\sum_{\alpha}X_{\alpha}^{\bar{\sigma}}(s+i\tau)\phi^{\sigma}(\tau)\mathrm{d}\tau+\frac{1}{\hbar^{2}}\int_{0}^{\beta\hbar}\mathrm{d}\tau\int_{0}^{\tau}\mathrm{d}\tau^{\prime}\sum_{\alpha\sigma}\phi^{\bar{\sigma}}(\tau)Y^{\sigma}_{\alpha}(\tau-\tau^{\prime})\phi^{\sigma}(\tau^{\prime})}.\end{split}$
(27)

The thermal equilibrium correlation function of lead $\alpha$ is given by

$\begin{split}G_{\alpha}^{\sigma}(t)=\sum_{k}\frac{v_{\alpha k}^{*}v_{\alpha
k}}{Q_{\alpha}}\mathrm{tr}_{\alpha}\\{e^{\frac{i}{\hbar}H_{\alpha}t}c^{\sigma}_{\alpha
k}e^{-\frac{i}{\hbar}H_{\alpha}t}c^{\bar{\sigma}}_{\alpha
k}e^{-\beta(H_{\alpha}-\mu_{\alpha}N_{\alpha})}\\}=\frac{1}{2\pi}\int_{-\infty}^{\infty}e^{\frac{i\sigma}{\hbar}\epsilon
t}\Gamma_{\alpha}(\epsilon)f_{\alpha}^{\sigma}(\epsilon)\mathrm{d}\epsilon,\end{split}$
(28)

where
$Q_{\alpha}=\mathrm{tr}_{\alpha}\left\\{e^{-\beta(H_{\alpha}-\mu_{\alpha}N_{\alpha})}\right\\}$
is the partition function of lead $\alpha$. $\Gamma_{\alpha}(\epsilon)$ is the
level-width function, which is defined as

$\Gamma_{\alpha}(\epsilon)=2\pi\sum_{k}v^{*}_{\alpha k}v_{\alpha
k}\delta(\epsilon-\epsilon_{\alpha k}).$ (29)

In what follows, a Lorentzian form is adopted,

$\Gamma_{\alpha}(\epsilon)=\frac{\Gamma_{\alpha}\Omega_{\alpha}^{2}}{\Omega_{\alpha}^{2}+(\epsilon-\mu_{\alpha})^{2}},$
(30)

where $\Gamma_{\alpha}$ quantifies the coupling strength between the molecular
level and lead $\alpha$. $\Omega_{\alpha}$ is the bandwidth of the spectral
energy distribution. The electron distribution in lead $\alpha$ at thermal
equilibrium is given by the Fermi distribution function,

$f_{\alpha}^{\sigma}(\epsilon)=\frac{1}{1+e^{\sigma\beta(\epsilon-\mu_{\alpha})}}.$
(31)

Various sum-over-poles decomposition schemes for the Fermi distribution
function have been proposed, including the Matsubara spectrum decomposition
scheme, the Padé spectrum decomposition scheme, Hu_2010_J.Chem.Phys._p101106 ;
Hu_2011_J.Chem.Phys._p244106 and the Fano spectrum decomposition
scheme.Cui_2019_J.Chem.Phys._p24110 ; Zhang_2020_J.Chem.Phys._p64107 In this
work, we employ the Padé decomposition and expand the correlation functions in
Eq. (28) as a sum over exponential functions,

$G_{\alpha}^{\sigma}(t-s)\simeq\sum_{p=0}^{P}\eta_{\alpha
p}e^{-\gamma^{\sigma}_{\alpha p}(t-s)},$ (32)

The explicit expressions of the coefficients $\\{\eta_{\alpha p}\\}$ and
exponents $\\{\gamma^{\sigma}_{\alpha p}\\}$ are provided in the supporting
material.

The terms in the second line of Eq. (27) describe the influence of the initial
correlation between the system and reservoirs with the following correlation
functions,

$\displaystyle X_{\alpha}^{\bar{\sigma}}(t+i\tau)$ $\displaystyle=$
$\displaystyle\frac{e^{-\bar{\sigma}\mu_{\alpha}\frac{\tau}{\hbar}}}{2\pi}\int_{-\infty}^{\infty}\mathrm{d}\epsilon
e^{-\frac{i\bar{\sigma}}{\hbar}\epsilon(t+i\tau)}\Gamma_{\alpha}(\epsilon)f^{\bar{\sigma}}_{\alpha}(\epsilon),$
(33) $\displaystyle\simeq$ $\displaystyle\sum_{p=0}^{P}\eta_{\alpha
p}^{*}e^{-\gamma^{\sigma}_{\alpha p}t}e^{-i\gamma_{\alpha p}\tau},$

and

$\displaystyle Y^{\sigma}_{\alpha}(\tau-\tau^{\prime})$ $\displaystyle=$
$\displaystyle\frac{1}{2\pi}\int_{-\infty}^{\infty}\mathrm{d}\epsilon
e^{\sigma\frac{\epsilon-\mu_{\alpha}}{\hbar}(\tau-\tau^{\prime})}\Gamma_{\alpha}(\epsilon)f^{\sigma}_{\alpha}(\epsilon)$
(34) $\displaystyle\simeq$ $\displaystyle\tilde{\eta}_{\alpha
0}\cos\Omega_{\alpha}\left(\frac{\beta}{2}-\frac{1}{\hbar}(\tau-\tau^{\prime})\right)+\sum_{p=1}^{\infty}\tilde{\eta}_{\alpha
p}\left(\sin{\gamma_{\alpha p}(\tau-\tau^{\prime})}\right).$

with $\tilde{\eta}_{\alpha 0}=\eta_{\alpha 0}e^{i\beta\Omega_{\alpha}/2}$ and
$\tilde{\eta}_{\alpha p}=i\eta_{\alpha p}$.

Given the above expansions of the correlations functions, we introduce a
Grassmann variable,

$\displaystyle\mathcal{B}_{\bm{a}}$ $\displaystyle=$
$\displaystyle-\frac{i}{\hbar}\int_{0}^{t}e^{-\gamma^{\sigma}_{\alpha
p}(t-s)}\left[\eta_{\alpha p}\varphi^{\sigma}(s)-\eta_{\alpha
p}^{*}\varphi^{\prime\sigma}(s)\right]\mathrm{d}s-\frac{\eta_{\alpha
p}^{*}e^{-\gamma_{\alpha
p}^{\sigma}t}}{\hbar}\int_{0}^{\beta\hbar}\phi^{\sigma}(\tau)e^{-i\gamma_{\alpha
p}\tau}\mathrm{d}\tau,$ (35)

with a multi-index subscript $\bm{a}=(\alpha,\sigma,p)$. For a given
$\mathcal{B}_{\bm{a}}$, it can be understood as being related to the process
where an electron from lead $\alpha$ is transferred to/back from
($\sigma=+/-$) the molecular electronic level with the characteristic
correlation time $\gamma_{\alpha p}^{-1}$. To proceed, we introduce a set of
auxiliary influence functionals
$\mathcal{F}^{(n)}=\mathcal{B}_{\bm{a}_{n}}\cdots\mathcal{B}_{\bm{a}_{1}}\mathcal{F}$.
By substituting the influence functional $\mathcal{F}$ in Eq. (II.3) with
$\mathcal{F}^{(n)}$, one readily constitutes a group of auxiliary density
matrices (ADMs), $\\{\rho^{(n)}_{\bm{a}_{1}\cdots\bm{a}_{n}}(t)\\}$, and their
time-derivative yields the equation of motionJin_2008_J.Chem.Phys._p234703

$\begin{split}\frac{\partial\rho^{(n)}_{\bm{a}_{1},\cdots,\bm{a}_{n}}(t)}{\partial
t}=&-\frac{i}{\hbar}\left[H_{s},\rho^{(n)}_{\bm{a}_{1},\cdots,\bm{a}_{n}}(t)\right]-\left(\sum_{l=1}^{n}\gamma_{\bm{a}_{l}}\right)\rho^{(n)}_{\bm{a}_{1},\cdots,\bm{a}_{n}}(t)\\\
&-\frac{i}{\hbar}\sum_{\bm{a}}\left(d^{\bar{\sigma}}s(\mathbf{x})\rho^{(n+1)}_{\bm{a}_{1},\cdots,\bm{a}_{n},\bm{a}}(t)-(-1)^{n}\rho^{(n+1)}_{\bm{a}_{1},\cdots,\bm{a}_{n},\bm{a}}(t)d^{\bar{\sigma}}s(\mathbf{x})\right)\\\
&-\frac{i}{\hbar}\sum_{l=1}^{n}(-1)^{n-l}\left(\eta_{\alpha_{l}p_{l}}d^{\sigma_{l}}s(\mathbf{x})\rho^{(n-1)}_{\bm{a}_{1},\cdots,\bm{a}_{l-1},\bm{a}_{l+1},\cdots,\bm{a}_{n}}(t)+(-1)^{n}\eta_{\alpha_{l}p_{l}}^{*}\rho^{(n-1)}_{\bm{a}_{1},\cdots,\bm{a}_{l-1},\bm{a}_{l+1},\cdots,\bm{a}_{n}}(t)d^{\sigma}s(\mathbf{x})\right).\end{split}$
(36)

The above equation actually has the same form as the one obtained from a
factorized initial condition.Jin_2008_J.Chem.Phys._p234703 However, the
initial values of the ADMs are yet unkown, and they are given by definition as

$\begin{split}\rho^{(n)}_{\bm{a}_{1},\cdots,\bm{a}_{n}}(\psi^{+}_{0},\tilde{\psi}^{-}_{0},t=0)=&\frac{1}{Q_{r}}\int^{\psi^{+}_{0}}_{\tilde{\psi}^{-}_{0}}\mathcal{D}\bm{\phi}^{+}(\tau)\mathcal{D}\bm{\phi}^{-}(\tau)e^{-\frac{1}{\hbar}S_{\beta}(\bm{\phi}^{+}(\tau),\bm{\phi}^{-}(\tau))}\\\
&\exp\left\\{\frac{1}{\hbar^{2}}\int_{0}^{\beta\hbar}\mathrm{d}\tau\int_{0}^{\tau}\mathrm{d}\tau^{\prime}\sum_{\alpha\sigma}\phi^{\bar{\sigma}}(\tau)Y^{\sigma}_{\alpha}(\tau-\tau^{\prime})\phi^{\sigma}(\tau^{\prime})\right\\}\left(1-h(\phi^{+}(\beta\hbar/2),\phi^{-}(\beta\hbar/2))\right)\mathcal{G}_{\bf{a}_{n}}\cdots\mathcal{G}_{\bf{a}_{1}}\end{split}$
(37)

with $\mathcal{G}_{\bf{a}}=-\frac{\eta_{\alpha
p}^{*}}{\hbar}\int_{0}^{\beta\hbar}\phi^{\sigma}(\tau)e^{-i\gamma_{\alpha
p}\tau}\mathrm{d}\tau$ encoding the initial system-bath correlation
information.

Following the ideas in Ref. Tanimura_2014_J.Chem.Phys._p44114, , we introduce
a set of imaginary time ADMs,

$\begin{split}\tilde{\varrho}^{(n)}_{\bm{a}_{1}\kappa_{1},\cdots,\bm{a}_{n}\kappa_{n}}(\psi^{+}_{0},\tilde{\psi}^{-}_{0},\tau,\lambda)=&\frac{1}{\prod_{\alpha}Q_{\alpha}}\int^{\psi^{+}_{0}}_{\tilde{\psi}^{-}_{0}}\mathcal{D}\bm{\phi}^{+}(\tau)\mathcal{D}\bm{\phi}^{-}(\tau)e^{-\frac{1}{\hbar}S_{\tau}(\bm{\phi}^{+}(\tau),\bm{\phi}^{-}(\tau))}\left(1-\lambda\theta(\tau-\beta\hbar/2)h(\phi^{+}(\beta\hbar/2),\phi^{-}(\beta\hbar/2))\right)\\\
&\exp\left\\{\frac{1}{\hbar^{2}}\int_{0}^{\tau}\mathrm{d}\tau^{\prime}\int_{0}^{\tau^{\prime}}\mathrm{d}\tau^{\prime\prime}\sum_{\alpha\sigma}\phi_{j}^{\bar{\sigma}}(\tau^{\prime})\tilde{\eta}_{\alpha
0}\cos\gamma_{\alpha
0}\left(\frac{\beta\hbar}{2}-\tau^{\prime}+\tau^{\prime\prime}\right)\phi^{\sigma}(\tau^{\prime\prime})\right\\}\\\
&\exp\left\\{\frac{1}{\hbar^{2}}\int_{0}^{\tau}\mathrm{d}\tau^{\prime}\int_{0}^{\tau^{\prime}}\mathrm{d}\tau^{\prime\prime}\sum_{\alpha
p\sigma}\phi^{\bar{\sigma}}(\tau^{\prime})\tilde{\eta}_{\alpha
p}\sin{\gamma_{\alpha
p}(\tau^{\prime}-\tau^{\prime\prime})}\phi^{\sigma}(\tau^{\prime\prime})\right\\}\tilde{\mathcal{G}}_{\mathbf{a}_{n}\kappa_{n}}\cdots\tilde{\mathcal{G}}_{\mathbf{a}_{1}\kappa_{1}}.\end{split}$
(38)

with
$\tilde{\mathcal{G}}_{\mathbf{a}\kappa}=\frac{1}{\hbar}\int_{0}^{\tau}\phi^{\sigma}(\tau^{\prime})\left[\kappa\cos(\gamma_{\alpha
p}\tau^{\prime})+\bar{\kappa}\sin(\gamma_{\alpha
p}\tau^{\prime})\right]\mathrm{d}\tau^{\prime}$. The additional index $\kappa$
in the above definition takes only the value of either 0 or 1, and
$\bar{\kappa}=1-\kappa$. The imaginary time $\tau$ varies from 0 to
$\beta\hbar$. Different from real time ADMs, all imaginary time ADMs are real-
valued.

By comparing Eq. (37) and Eq. (38), one can verify that the real time ADMs at
$t=0$ are related to the imaginary time ADMs of the same tier and can be
obtained through the following transformation,

$\rho^{(n)}_{\bm{a}_{1},\cdots,\bm{a}_{n}}(t=0)=\frac{1}{Q_{r}^{s}}\sum_{\kappa_{1},\cdots,\kappa_{n}}\left(\prod_{l=1}^{n}\eta^{*}_{\alpha_{l}p_{l}}\left(-\kappa_{l}+i\bar{\kappa}_{l}\right)\right)\tilde{\varrho}^{(n)}_{\bm{a}_{1}\kappa_{1},\cdots,\bm{a}_{n}\kappa_{n}}(\beta\hbar,\lambda=1).$
(39)

The partition function of the reactant is given by

$Q_{r}^{s}=Q_{r}/\prod_{\alpha}Q_{\alpha}=\mathrm{tr}_{s}\left\\{\tilde{\varrho}^{(0)}(\beta\hbar,\lambda=1)\right\\}.$
(40)

Differentiating the imaginary time ADMs defined in Eq. (38) with respect to
$\tau$, we have

$\begin{split}\frac{\partial\tilde{\varrho}^{(n)}_{\bm{a}_{1}\kappa_{1},\cdots,\bm{a}_{n}\kappa_{n}}(\tau,\lambda)}{\partial\tau}=&-\frac{1}{\hbar}H_{s}\tilde{\varrho}^{(n)}_{\bm{a}_{1}\kappa_{1},\cdots,\bm{a}_{n}\kappa_{n}}(\tau,\lambda)\\\
&+\frac{1}{\hbar}\sum_{(\alpha\sigma\kappa)}d^{\bar{\sigma}}s(\mathbf{x})\tilde{\eta}_{\alpha
0}\left[\kappa\cos{\gamma_{\alpha
0}\left(\frac{\beta\hbar}{2}-\tau\right)}-\bar{\kappa}\sin{\gamma_{\alpha
0}\left(\frac{\beta\hbar}{2}-\tau\right)}\right]\tilde{\varrho}^{(n+1)}_{\bm{a}_{1}\kappa_{1},\cdots,\bm{a}_{n}\kappa_{n},(\alpha\sigma
0\kappa)}(\tau,\lambda)\\\
&+\frac{1}{\hbar}\sum_{p=1}^{P}\sum_{\alpha\sigma\kappa}d^{\bar{\sigma}}s(\mathbf{x})\tilde{\eta}_{\alpha
p}\left[\kappa\sin{\left(\gamma_{\alpha
p}\tau\right)}-\bar{\kappa}\cos{\left(\gamma_{\alpha
p}\tau\right)}\right]\tilde{\varrho}^{(n+1)}_{\bm{a}_{1}\kappa_{1},\cdots,\bm{a}_{n}\kappa_{n},(\alpha\sigma
p\kappa)}(\tau,\lambda)\\\
&+\frac{1}{\hbar}\sum_{l=1}^{n}(-1)^{n-l}d^{\sigma_{l}}s(\mathbf{x})\left[\kappa_{l}\cos(\gamma_{\alpha_{l}p_{l}}\tau)+\bar{\kappa}_{l}\sin(\gamma_{\alpha_{l}p_{l}}\tau)\right]\tilde{\varrho}^{(n-1)}_{\bm{a}_{1}\kappa_{1},\cdots,\bm{a}_{l-1}\kappa_{l-1},\bm{a}_{l+1}\kappa_{l+1},\cdots,\bm{a}_{n}\kappa_{n}}(\tau,\lambda),\end{split}$
(41)

which is numerically solvable with the explicit initial condition. The zeroth
tier imaginary time ADM is an identity matrix in the system subspace,
$\tilde{\varrho}^{(0)}(\tau=0)=\bm{I}$. All other ADMs beyond the zeroth tier
are zero, i.e., $\tilde{\varrho}^{(n>0)}(\tau=0)=\bm{0}$.

It is worth noting that the hierarchies in Eq. (36) and Eq. (41) both
automatically terminate at a certain tier. The Grassmann variables
$\mathcal{B}_{\bm{a}}$ and $\tilde{\mathcal{G}}_{\mathbf{a}\kappa}$ obey the
anti-commutation relation, which implies that all indices in
$\rho^{(n)}_{\bm{a}_{1}\cdots\bm{a}_{n}}(t)$ or
$\tilde{\varrho}^{(n)}_{\bm{a}_{1}\kappa_{1},\cdots,\bm{a}_{n}\kappa_{n}}$
should be different. The highest tier for nonzero real time ADMs is
$L=2N_{e}N_{\alpha}(P+1)$, with $N_{e}$, $N_{\alpha}$, and $P$ being the
number of molecular electronic levels, the number of reservoirs, and Padé
poles, respectively. For the nonzero imaginary time ADMs, the highest tier is
increased to $2L$. In general, $L$ is very large, and a truncation of the
hierarchy is necessary. To obtain the real time ADMs at $t=0$, as indicated by
Eq. (39), we should truncate the imaginary time hierarchy at least at the same
tier as the real time hierarchy. However, in our simulations, we found that
the imaginary time equation (Eq. (41)) generally converges faster than the
real time equation (Eq. (36)) with respect to the hierarchical truncation
tier.

With the zeroth tier real and imaginary time ADMs, one can obtain a variety of
thermal and dynamical properties for the system of interest. For instance, the
flux correlation function in Eq. (24) is evaluated as the trace expression

$C_{\rm f}(t)=\mathrm{tr}_{s}\left\\{F\rho^{(0)}(t)\right\\}/Q_{r}^{s},$ (42)

and the system population in joint equilibrium with its reservoirs can be
directly obtained as

$\rho_{s}^{\rm eq}(x)=\frac{\langle
x|\tilde{\varrho}^{(0)}(\beta\hbar,\lambda=0)|x\rangle}{\mathrm{tr}_{s}\left\\{\tilde{\varrho}^{(0)}(\beta\hbar,\lambda=0)\right\\}}.$
(43)

Bath-related properties are enciphered in the higher-tier ADMs. For example,
the electronic current flowing from lead $\alpha$ to the molecule can be
defined as the time derivative of the average number of electrons in the lead,
and is obtained through the first-tier ADMs,

$\begin{split}I_{\alpha}=&-e\frac{d\langle N_{\alpha}\rangle}{dt}\\\
=&\frac{ie}{\hbar}\sum_{p}\mathrm{tr}_{s}\left\\{ds(\mathbf{x})\rho^{(1)}_{(\alpha,+,p)}-d^{\dagger}s(\mathbf{x})\rho^{(1)}_{(\alpha,-,p)}\right\\},\end{split}$
(44)

where $e$ denotes the electron charge. The interaction energy at position
$\mathbf{x}$ is also related to the first-tier ADMs,

$\begin{split}\langle H_{c}(\mathbf{x})\rangle=&\sum_{\alpha
p}s(\mathbf{x})\langle\mathbf{x}|\mathrm{tr}_{s_{e}}\left\\{d\rho^{(1)}_{(\alpha,+,p)}+d^{\dagger}\rho^{(1)}_{(\alpha,-,p)}\right\\}|\mathbf{x}\rangle,\end{split}$
(45)

where $\mathrm{tr}_{s_{e}}$ denotes tracing over system electronic DoFs. The
derivation details are provided in the supporting material.

a)

b)

c)

Figure 1: (a) Schematic representation of intramolecular proton transfer in a
molecular junction. (b) Double well potential for the proton transfer reaction
model of the isolated molecule (see Eq. (47) with $l=0.5$ Å and $V_{b}=0.4$
eV). The black horizontal and colored wavy lines denote the energy and
wavefunction of the eigenstates, respectively. The first and the second state
are nearly degenerate. (c) Switch function defined in Eq. (48) with
$\sigma=0.1$.

##  III  Results

In this section, we use the method outlined above to investigate the reaction
rate of intramolecular proton transfer dynamics in a molecular junction. We
first introduce the model system used and provide the simulations details.
Then, results for a broad range of applied bias voltages and molecule-lead
couplings are analyzed and discussed.

###  III.1 Model

We consider a one-dimensional model for studying the proton transfer dynamics
in a molecular junction. The system Hamiltonian is given by

$H_{s}=\frac{p^{2}}{2m}+V_{0}(x)+\epsilon(x)d^{\dagger}d,$ (46)

As shown in Fig. 1 (a), the proton is translocated between the donor and
acceptor moieties within the molecule. We shall assume that the proton
transfer between donor and acceptor sites changes the molecule-lead couplings,
but without loss of generality, the local potential energy surfaces at donor
and acceptor sites are assumed to be identical. Therefore, the potential
energy of the neutral molecule along the reaction coordinate is assumed to
have the form of a symmetric double-well

$V_{0}(x)=\frac{V_{b}}{l^{4}}(x+l)^{2}(x-l)^{2},$ (47)

which is characterized by an energy barrier $V_{b}$ at $x=0$. Two local minima
are located at $x=\pm l$, corresponding to two stable conformations of an
isolated molecule.

The model was used in previous work to investigate the influence of proton
transfer on the conductance in molecular
junctions.Hofmeister_J.Chem.Phys._2017_p092317 Specifically, Fig. 1 (b)
depicts the potential with the parameters $l=0.5\textrm{ \AA}$ and $V_{b}=0.4$
eV. The energy levels and wavefunctions of the seven lowest eigenstates are
shown as well. The vibrational frequency is characterized as the curvature
near the local minima, i.e.,
$\hbar\omega_{0}=\left.\sqrt{\frac{1}{m}\frac{\partial^{2}V_{0}(x)}{\partial
x^{2}}}\right|_{x=\pm l}=0.2$ eV (the proton mass is 1 amu), which
characterizes the timescale of the nuclear motion. In general, the potential
energy surface for the charged molecule may be different from the neutral one.
Here, we assume for simplicity that the charging only leads to a constant
energy shift by the electron affinity $\epsilon_{0}$, i.e.,
$\epsilon(x)=\epsilon_{0}$. In the following calculations, we adopt a value of
$\epsilon_{0}=0.1$ eV. More general cases where $\epsilon(x)$ is a function of
$x$ will be the subject of future work.

Besides, we assume that the dependence of molecule-lead coupling on the proton
position (see Eq. (5)) is characterized by the function

$s(x)=\frac{1+\sigma}{2}-\frac{1-\sigma}{2}\tanh(x),$ (48)

which is visualized in Fig. 1 (c) with $\sigma=0.1$. As a result, in the
current model, the conductance is larger when the proton is located in the
left well and smaller in the right well. Furthermore, we point out that
because of the assumed uniform charging energy at any proton position, the
coupling of the nuclear motion to the electrons arises exclusively from the
coordinate-dependent coupling to the metal electrons.

###  III.2 Simulation details

To obtain the proton transfer rate, a convenient choice of the diving surface
is $x_{\rm ds}=0$, but we should stress that, in principle, the obtained rate
is invariant to the choice of the dividing surface. For all calculations
presented in this work, we assume that the molecule is coupled symmetrically
to the left and right leads, $\Gamma_{L/R}=\Gamma/2$ and a temperature of
$T=300$ K is adopted. The proton degree of freedom is represented using the
sine discrete variable representation (DVR).Colbert_1992_J.Chem.Phys._p1982 ;
Echave_1992_Chem.Phys.Lett._p225 ; Seideman_1992_J.Chem.Phys._p4412 A number
of 26 DVR grid points is adequate to give converged results. The bandwidth in
Eq. (30) is taken as $\Omega_{L/R}=1$ eV. At least 15 Padé poles are required.
The hierarchical equations (Eq. (36) and Eq. (41)) are truncated at the third
tier ($L=3$) to assure convergence unless otherwise specified. We use a simple
truncation scheme for both real and imaginary time equations, namely, setting
the ADMs beyond the $L$th tier to zero. Both real and imaginary time equations
are propagated using the fourth-order Runge-Kutta method. The ratio of
reactant and product steady state population $\chi$ in Eq. (II.2) is obtained
employing an efficient iterative steady state solver developed by Kaspar et
al. Kaspar_J.Phys.Chem.A_2021_p5190 . The numerical integration is performed
on NVIDIA Tesla V100 GPUs.

###  III.3 Reaction at different bias voltage regimes

First, we study the proton transfer dynamics in nonequilibrium situations,
where the molecule is connected to two electrodes at a nonzero applied bias
voltage. The bias voltage drops symmetrically on both electrodes, i.e.,
$\mu_{L}=-\mu_{R}=e\Phi/2$. The coupling to the two leads is fixed at
$\Gamma_{L/R}=0.05$ eV.

Fig. 2 (a) and (b) show the flux correlation function $C_{\rm f}(t)$ for two
different bias voltages in the off and near-resonant transport regime,
respectively. At low biases, for instance, $\Phi=0.1$ and $0.2$ V, the flux
correlation functions exhibit significant oscillations with a period of about
21 fs ($\sim 0.2$ eV), which can be assigned to the vibrational motion at the
bottom of the double well. After tens of picoseconds, the oscillations are
damped out and $C_{\rm f}(t)$ reaches a well-defined plateau. In these cases,
the approximate expression Eq. (15) can be used to obtain the nonequilibrium
rate constant for current-induced proton transfer process, as the time
separation condition is satisfied and the initial population of the reactant
is very close to one ($P_{r}(0)=0.99993$).

In the resonant transport regime, as illustrated in Fig. 2 (c) for $\Phi=4$ V,
the oscillations of the flux correlation function vanish after a few
picoseconds and the reaction takes place faster. As a consequence, the flux
correlation function $C_{\rm f}(t)$ decreases over time because the timescale
separation is no longer satisfied and the reactant population changes
significantly from its initial value. Nevertheless, the reaction obeys for
longer time the first order kinetic scheme, Eq. (10), with the rate constant
$k$ given by the generalized expression, Eq. (II.2), which exhibits a plateau,
as shown by the orange line.

Fig. 3 displays the forward reaction rate $k$ and the ratio $k/k_{b}$ in the
inset against the bias voltages $\Phi$ ranging from 0 to 4 V. The reaction
rate increases with bias voltage in a bimodal way with a steep rise at lower
bias voltages and a slow increase at higher bias voltages.

At zero bias voltage, there is no electrical current flowing through the
molecule. However, the proton transfer reaction rate is not zero due to
quantum tunneling and thermal activation effects. It should be noted that due
to the small mass of the proton, the vibrational frequency $\omega_{0}$ is
relatively high. As such, even at the room temperature, proton transfer due to
quantum tunneling (rate roughly estimated as
$\omega_{0}\exp\left(-V_{b}/(\hbar\omega_{0}/2\pi)\right)$) dominates over
thermally activation processes (rate
$\sim\omega_{0}\exp\left(-V_{b}/(k_{B}T)\right)$). Interestingly, although we
assume a symmetric molecular potential surface, the forward reaction rate is
smaller than the backward reaction rate ($k/k_{b}=0.167$). This is due to the
coupling of the molecule to the leads, which is a function of the proton
position $x$, where more details will be elaborated in Sec. III.4.

a) $\Phi=$ 0.1 V

b) $\Phi=$ 0.2 V

c) $\Phi=$ 4 V

Figure 2: Nonequilibrium flux correlation function $C_{\rm f}(t)$ (see Eq.
(24)) for three different bias voltages in the off-resonant, near resonant,
and deep resonant transport regime, respectively. In addition, the rate
expression $k(t)$ (see Eq. (II.2)) at a high bias voltage is shown in (c). In
(a) and (b) for $\Phi=0.1$ and 0.2 V, $k(t)$ is not shown as the results
overlap $C_{\rm f}(t)$.

Turning on the bias voltage, the electrons transfer into and out of the
molecule with a timescale characterized by $\hbar/\Gamma$, and this process
can be accompanied by vibrational excitations. Thus, the proton transfer is
additionally induced by current-induced vibrational
excitation,Elste_Appl.Phys.A_2008_p345–354 which is much faster than quantum
tunneling and thermal activation processes. At low bias voltages, more
vibrationally coupled electron transport channels are opened up with
increasing bias voltage. As a result, the forward reaction rate is
significantly increased. However, when the bias voltage is sufficiently large,
the energy transmitted into the molecule by a single electron transferring
from the left lead to the right lead can overcome the reaction barrier for
both forward and back reaction. In this case, by further increasing the bias
voltage, the reaction rate increases only slightly and the ratio $k/k_{b}$
approaches 1.

Figure 3: Forward reaction rate constant $k$ and the ratio $k/k_{b}$ as a
function of the applied bias voltage.

###  III.4 Influence of molecule-lead coupling

It is known that the coupling between molecule and leads can have a profound
impact on the molecular properties.Molen_J.Phys.:Condens.Matter_2010_p133001 ;
Chen_Chem.Rev._2020_p2879 In this subsection, we investigate the influence of
the molecule-lead coupling strength on the proton transfer rate and molecular
properties.

Fig. 4 displays the proton transfer rate as a function of the molecule-lead
coupling strength in the range of $\Gamma=0.01$ to 0.2 eV for both equilibrium
($\Phi=0$ V) and nonequilibrium ($\Phi=4$ V) cases. The ratios of the forward
and back reaction rates $k/k_{b}$ are shown in the insets. We use a smaller
barrier height $V_{b}=0.3$ eV, which facilitates the calculation of the
reaction rate as the oscillations are damped at a shorter timescale.

a) $\Phi=0$ V

b) $\Phi=4$ V

Figure 4: Forward reaction rate constant $k$ as a function of the molecule-
lead coupling strength $\Gamma$ for two different bias voltages. Shown in the
inset are the ratios of forward and backward reaction rate $k/k_{b}$ against
$\Gamma$. The barrier height is $V_{b}=0.3$ eV. The results for $\Phi=$ 0 V
are obtained with a fourth-tier hierarchical truncation.

At equilibrium with zero bias voltage, we observe that the forward reaction
rate $k$ and the ratio $k/k_{b}$ decrease monotonically with increasing
molecule-lead coupling strength $\Gamma$. For $\Gamma>0.2$ eV, the flux
correlation function still shows a strongly oscillatory structure on the scale
of hundreds of picoseconds (data not shown), and the propagation to a much
longer timescale is required to obtain the reaction rate.

At finite bias voltage in the resonant transport regime, as shown in Fig. 4
(b) for $\Phi=4$ V, the reaction rates are about two orders of magnitude
larger, but exhibit the same dependence on $\Gamma$. This suggests that
inherent molecular properties are the main reason for the behavior rather than
the external bias voltage.

a)

b)

c)

Figure 5: Interaction energy between the system and reservoirs (a) and
equilibrium populations (b) as well as the potential of mean force (c) at
$300$ K. The definitions of the above quantities are given in Eq. (45), Eq.
(43), and Eq. (49), respectively, and are obtained with a fourth-tier
hierarchical truncation. The lines with color varying from cyan to purple
correspond to different molecule-lead couplings ranging from 0 to 2 eV with a
spacing of 0.02 eV. Other parameters are: bias voltage $\Phi=0$ V, barrier
height $V_{b}=0.3$ eV.

To understand this behavior, we calculate the interaction energy $\langle
H_{c}(x)\rangle$ (Eq. (45)), the equilibrium population $\rho^{\rm eq}_{s}(x)$
(Eq. (43)), as well as the potential of mean force $V_{\rm PMF}(x)$ for the
reactive nuclear motion at equilibrium, which is given
byDou_J.Chem.Phys._2016_p74109

$\begin{split}V_{\rm PMF}(x)=&-\frac{1}{\beta}\ln\langle x|{\rm
tr}_{e}\\{e^{-\beta H_{el}}\\}|x\rangle\end{split}$ (49)

where $\mathrm{tr}_{e}$ denotes tracing over all electronic DOFs. ${\rm
tr}_{e}\\{e^{-\beta H_{el}}\\}$ can be obtained exactly through the imaginary
time HEOM approach, as presented in Eq. (41). The results are displayed in
Fig. 5 for various molecule-lead coupling strength ranging from 0 to 2 eV.

For vanishing molecule-lead coupling, the interaction energy is zero and the
potential of mean force reproduces the symmetric double-well potential of the
isolated neutral molecule. The equilibrium populations in the reactant and
product region are equal, thus we have $k/k_{b}=1$. Increasing the molecule-
lead coupling $\Gamma$, the interaction energy $\langle H_{c}(x)\rangle$ shows
a dip in the reactant region with a minimum at around $x=-0.5$ Å, which is due
to the stronger coupling of the molecule to the leads when the proton is
placed closer to the donor site. The negative interaction energy in the
reactant region leads to a stabilization effect. As a consequence, the
molecule is more likely to populate the reactant region for increasing
$\Gamma$, and $k/k_{b}$ decreases. In particular, for strong coupling
$\Gamma>0.5$ eV, the equilibrium population in the product region as well as
the ratio $k/k_{b}$ approach zero, as shown in Fig. 5 (b) and the inset of
Fig. 4 (a), respectively. The potential of mean force becomes asymmetric and
the barrier top shifts towards the positive $x$ direction. The barrier height
between the bottom of the left well and the barrier top is enhanced with
increasing $\Gamma$, which results in a decrease of the forward reaction rate
$k$.

The results shown here indicate the crucial role of a coordinate-dependent
molecule-lead coupling in reshaping the reactive potential surface and the
reaction rates. This not only deepens our understanding of chemical reactions
at molecule-metal interfaces. More importantly, it may be utilized to steer
the chemical reaction in the desired direction by engineering the molecule-
lead coupling strength.

##  IV Conclusions

We have presented a nonequilibrium reaction rate theory based on an extension
of the flux correlation function formalism to nonequilibrium conditions.
Moreover, to efficiently calculate flux correlation functions within this
approach, we have proposed a mixed real and imaginary time HEOM method for
open quantum systems that are coupled to fermionic reservoirs. The combined
method can be used to calculate rate constants of chemical reactions in
nonequilibrium situations in a fully quantum mechanical and numerically exact
way.

As an illustrative example, we have considered a model of intramolecular
proton transfer in a molecular junction and calculated reaction rates for
current-induced proton transfer for a range of different bias voltages and
molecule-lead coupling strengths. For small bias voltages, time-scale
separation allows to use directly the flux correlation function to determine
the rate constant. For higher bias voltages, a renormalized flux-correlation
function formulation is required to calculate the rate constant. The results
also show that a strong molecule-lead interaction that depends on the reaction
coordinate can reshape the effective potential energy surface and influence
the reaction dynamics profoundly.

In the present work, we have considered a simple model for proton transfer
with a single reaction coordinate. Using the recent extension of the HEOM
method, which combines it with the matrix product state
representation,yaling_heom_tt it should also be applicable to more general
models and different reactions involving several nuclear degrees of freedom
and multiple electronic states of the molecule.

##  V Acknowledgements

The authors thank Jakob Bätge helpful discussions. This work was supported by
the German Research Foundation (DFG). Y.K. was supported by the Alexander von
Humboldt Foundation. The work of A.E. was funded by the DFG (grant 453644843).
U.P. acknowledges support by the Israel Science foundation (1243/18) and by
the U.S. - Israel Binational Science foundation. Furthermore, the authors
acknowledge support by the state of Baden-Württemberg through bwHPC and the
German Research Foundation (DFG) through grant no INST 40/575-1 FUGG (JUSTUS 2
cluster).

## Supplementary Material

See the supplementary material for the details of (1) the explicit expressions
of the coefficients $\\{\eta_{\alpha p}\\}$ and exponents
$\\{\gamma^{\sigma}_{\alpha p}\\}$ in Eq. (32); (2) derivation details of the
interaction energy and current expression in Eq. (45) and Eq. (44),
respectively; (3) the flux correlation function $C_{\rm f}(t)$ calculated with
different initial density operators.

## Data Availability

The data that support the findings of this study are available from the
corresponding author upon reasonable request.

## References

  * (1) R. Wyatt,  Dynamics of molecules and chemical reactions,  CRC Press, 1996.
  * (2) P. Hänggi, P. Talkner, and M. Borkovec,  Rev. Mod. Phys. 62, 251 (1990).
  * (3) A. Nitzan,  Chemical dynamics in condensed phases: relaxation, transfer and reactions in condensed molecular systems,  Oxford university press, 2006.
  * (4) W. H. Miller,  J. Chem. Phys. 62, 1899 (1975).
  * (5) W. H. Miller,  J. Chem. Phys. 63, 1166 (1975).
  * (6) D. Chandler,  J. Chem. Phys. 68, 2959 (1978).
  * (7) D. G. Truhlar, W. L. Hase, and J. T. Hynes,  J. Phys. Chem. 87, 2664 (1983).
  * (8) M. Topaler and N. Makri,  J. Chem. Phys. 101, 7500 (1994).
  * (9) E. Geva, Q. Shi, and G. A. Voth,  J. Chem. Phys. 115, 9209 (2001).
  * (10) E. Pollak and P. Talkner,  Chaos 15, 026116 (2005).
  * (11) U. Manthe,  J. Chem. Phys. 102, 9205 (1995).
  * (12) U. Manthe,  J. Chem. Phys. 128, 064108 (2008).
  * (13) H. Wang, M. Thoss, and W. H. Miller,  J. Chem. Phys. 112, 47 (2000).
  * (14) J. O. Richardson and M. Thoss,  J. Chem. Phys. 141, 074106 (2014).
  * (15) W. H. Miller, S. D. Schwartz, and J. W. Tromp,  J. Chem. Phys. 79, 4889 (1983).
  * (16) W. H. Thompson and W. H. Miller,  J. Chem. Phys. 106, 142 (1997).
  * (17) U. Manthe,  J. Theor. Comput. Chem. 1, 153 (2002).
  * (18) H. Wang, D. E. Skinner, and M. Thoss,  J. Chem. Phys. 125, 174502 (2006).
  * (19) I. R. Craig, M. Thoss, and H. Wang,  J. Chem. Phys. 127, 144503 (2007).
  * (20) T. Kumagai et al.,  Phys. Rev. Lett. 100, 166101 (2008).
  * (21) T. Kumagai et al.,  Phys. Rev. B 79, 035423 (2009).
  * (22) P. Liljeroth, J. Repp, and G. Meyer,  Science 317, 1203 (2007).
  * (23) T. Kumagai et al.,  Phys. Rev. B 81, 045402 (2010).
  * (24) D. Peller et al.,  Nature 585, 58 (2020).
  * (25) J. N. Ladenthin et al.,  ACS nano 9, 7287 (2015).
  * (26) T. Kumagai et al.,  Nat. Chem. 6, 41 (2014).
  * (27) H. Böckmann et al.,  Nano Lett. 16, 1034 (2016).
  * (28) J. Koch, M. Semmelhack, F. Von Oppen, and A. Nitzan,  Phys. Rev. B 73, 155306 (2006).
  * (29) S. V. Aradhya and L. Venkataraman,  Nat. Nanotechnol. 8, 399 (2013).
  * (30) Y. Zang et al.,  Nano Lett. 21, 673 (2020).
  * (31) X. Huang et al.,  Sci. Adv. 5, eaaw3072 (2019).
  * (32) Q. Schaeverbeke, R. Avriller, T. Frederiksen, and F. Pistolesi,  Phys. Rev. Lett. 123, 246601 (2019).
  * (33) R. Avriller, Q. Schaeverbeke, T. Frederiksen, and F. Pistolesi,  Phys. Rev. B 104, L241403 (2021).
  * (34) H. Chen et al.,  J. Am. Chem. Soc. 143, 8476 (2021).
  * (35) A. Erpenbeck and M. Thoss,  J. Chem. Phys. 151, 191101 (2019).
  * (36) A. Erpenbeck, Y. Ke, U. Peskin, and M. Thoss,  Phys. Rev. B 102, 195421 (2020).
  * (37) Y. Ke, A. Erpenbeck, U. Peskin, and M. Thoss,  J. Chem. Phys. 154, 234702 (2021).
  * (38) C. Sabater, C. Untiedt, and J. M. van Ruitenbeek,  Beilstein J. Nanotechnol. 6, 2338 (2015).
  * (39) H. Li et al.,  J. Am. Chem. Soc. 137, 5028 (2015).
  * (40) H. Li et al.,  J. Am. Chem. Soc. 138, 16159 (2016).
  * (41) B. Capozzi et al.,  Nano Lett. 16, 3949 (2016).
  * (42) S. Hammes-Schiffer and A. V. Soudackov,  J. Phys. Chem. B 112, 14108 (2008).
  * (43) J. C. Tully and R. K. Preston,  J. Chem. Phys. 55, 562 (1971).
  * (44) N. Shenvi, S. Roy, and J. C. Tully,  J. Chem. Phys. 130, 174107 (2009).
  * (45) S. A. Fischer, C. T. Chapman, and X. Li,  J. Chem. Phys. 135, 144102 (2011).
  * (46) M. Head-Gordon and J. C. Tully,  J. Chem. Phys. 103, 10137 (1995).
  * (47) W. Dou and J. E. Subotnik,  J. Chem. Phys. 144, 024116 (2016).
  * (48) M. Galperin and A. Nitzan,  J. Phys. Chem. Lett 6, 4898 (2015).
  * (49) J.-T. Lü, B.-Z. Hu, P. Hedegård, and M. Brandbyge,  Prog. Surf. Sci. 94, 21 (2019).
  * (50) J.-T. Lü, M. Brandbyge, P. Hedegård, T. N. Todorov, and D. Dundas,  Phys. Rev. B 85, 245444 (2012).
  * (51) A. A. Dzhioev and D. Kosov,  J. Chem. Phys. 135, 074701 (2011).
  * (52) A. A. Dzhioev, D. S. Kosov, and F. Von Oppen,  J. Chem. Phys. 138, 134103 (2013).
  * (53) R. Pozner, E. Lifshitz, and U. Peskin,  Nano Lett. 14, 6244 (2014).
  * (54) R. J. Preston, V. F. Kershaw, and D. S. Kosov,  Phys. Rev. B 101, 155415 (2020).
  * (55) A. Erpenbeck, C. Schinabeck, U. Peskin, and M. Thoss,  Phys. Rev. B 97, 235452 (2018).
  * (56) P. Liljeroth, J. Repp, and G. Meyer,  Science 317, 1203 (2007).
  * (57) A. L. Sobolewski,  Phys. Chem. Chem. Phys. 10, 1243 (2008).
  * (58) C. Benesch et al.,  J. Phys. Chem. C 113, 10315 (2009).
  * (59) C. Hofmeister et al.,  J. Mol. Model. 20, 1 (2014).
  * (60) D. Weckbecker, P. Coto, and M. Thoss,  J. Phys. Chem. Lett. 12, 413 (2020).
  * (61) C. Hofmeister, P. B. Coto, and M. Thoss,  J. Chem. Phys. 146, 092317 (2017).
  * (62) H. Ness,  Phys. Rev. E 90, 062119 (2014).
  * (63) H. Ness,  Entropy 19, 158 (2017).
  * (64) E.-D. Fung et al.,  Nano Lett. 19, 2555 (2019).
  * (65) Y. Tanimura and R. Kubo,  J. Phys. Soc. Jpn. 58, 101 (1989).
  * (66) Y. Tanimura,  J. Chem. Phys. 153, 020901 (2020).
  * (67) Y.-a. Yan, F. Yang, Y. Liu, and J. Shao,  Chem. Phys. Lett. 395, 216 (2004).
  * (68) R.-X. Xu, P. Cui, X.-Q. Li, Y. Mo, and Y. Yan,  J. Chem. Phys. 122, 041103 (2005).
  * (69) A. Ishizaki and Y. Tanimura,  J. Phys. Soc. Jpn. 74, 3131 (2005).
  * (70) Q. Shi, L. Chen, G. Nan, R.-X. Xu, and Y. Yan,  J. Chem. Phys. 130, 084105 (2009).
  * (71) J. Hu, R.-X. Xu, and Y. Yan,  J. Chem. Phys. 133, 101106 (2010).
  * (72) C. Kreisbeck, T. Kramer, M. Rodriguez, and B. Hein,  J. Chem. Theory Comput. 7, 2166 (2011).
  * (73) D. Suess, A. Eisfeld, and W. Strunz,  Phys. Rev. Lett. 113, 150403 (2014).
  * (74) D. Hou et al.,  J. Chem. Phys. 142, 104112 (2015).
  * (75) Z. Tang, X. Ouyang, Z. Gong, H. Wang, and J. Wu,  J. Chem. Phys. 143, 224112 (2015).
  * (76) M. Tsuchimoto and Y. Tanimura,  J. Chem. Theory Comput. 11, 3859 (2015).
  * (77) Y. Ke and Y. Zhao,  J. Chem. Phys. 145, 024101 (2016).
  * (78) C.-Y. Hsieh and J. Cao,  J. Chem. Phys. 148, 014103 (2018).
  * (79) K. Nakamura and Y. Tanimura,  Phys. Rev. A 98, 012109 (2018).
  * (80) Q. Shi, Y. Xu, Y. Yan, and M. Xu,  J. Chem. Phys. 148, 174102 (2018).
  * (81) A. Erpenbeck, C. Hertlein, C. Schinabeck, and M. Thoss,  J. Chem. Phys. 149, 064106 (2018).
  * (82) H.-D. Zhang et al.,  J. Chem. Phys. 152, 064107 (2020).
  * (83) T. Ikeda and G. D. Scholes,  J. Chem. Phys. 152, 204101 (2020).
  * (84) Y. Tanimura,  J. Chem. Phys. 141, 044114 (2014).
  * (85) J. Jin, X. Zheng, and Y. Yan,  J. Chem. Phys. 128, 234703 (2008).
  * (86) X. Zheng, J. Jin, and Y. Yan,  J. Chem. Phys. 129, 184112 (2008).
  * (87) X. Zheng, Y. Yan, and M. Di Ventra,  Phys. Rev. Lett. 111, 086601 (2013).
  * (88) Y. Yan,  J. Chem. Phys. 140, 054105 (2014).
  * (89) L. Ye et al.,  WIREs Comput Mol Sci 6, 608 (2016).
  * (90) H.-D. Zhang, Q. Qiao, R.-X. Xu, X. Zheng, and Y. Yan,  J. Chem. Phys. 147, 044105 (2017).
  * (91) L. Han, H.-D. Zhang, X. Zheng, and Y. Yan,  J. Chem. Phys. 148, 234108 (2018).
  * (92) C. Schinabeck, A. Erpenbeck, R. Härtle, and M. Thoss,  Phys. Rev. B 94, 201407 (2016).
  * (93) C. Schinabeck, R. Härtle, and M. Thoss,  Phys. Rev. B 97, 235429 (2018).
  * (94) A. Erpenbeck, L. Götzendörfer, C. Schinabeck, and M. Thoss,  The European Physical Journal Special Topics 227, 1981 (2019).
  * (95) C. Schinabeck and M. Thoss,  Phys. Rev. B 101, 075422 (2020).
  * (96) R. P. Feynman, A. R. Hibbs, and D. F. Styer,  Quantum mechanics and path integrals,  Courier Corporation, 2010.
  * (97) J. Zinn-Justin,  Path integrals in quantum mechanics,  Oxford University Press, 2010.
  * (98) J. Hu, M. Luo, F. Jiang, R.-X. Xu, and Y. Yan,  J. Chem. Phys. 134, 244106 (2011).
  * (99) L. Cui, H.-D. Zhang, X. Zheng, R.-X. Xu, and Y. Yan,  J. Chem. Phys. 151, 024110 (2019).
  * (100) C. Hofmeister, P. B. Coto, and M. Thoss,  J. Chem. Phys. 146, 092317 (2017).
  * (101) D. T. Colbert and W. H. Miller,  J. Chem. Phys. 96, 1982 (1992).
  * (102) J. Echave and D. C. Clary,  Chem. Phys. Lett. 190, 225 (1992).
  * (103) T. Seideman and W. H. Miller,  J. Chem. Phys. 96, 4412 (1992).
  * (104) C. Kaspar and M. Thoss,  J. Phys. Chem. A 125, 5190 (2021).
  * (105) F. Elste, G. Weick, C. Timm, and F. von Oppen,  Appl. Phys. A 93, 345 (2008).
  * (106) S. J. van der Molen and P. Liljeroth,  J. Phys.: Condens. Matter 22, 133001 (2010).
  * (107) H. Chen, W. Zhang, M. Li, G. He, and X. Guo,  Chem. Rev. 120, 2879 (2020).
  * (108) W. Dou, A. Nitzan, and J. E. Subotnik,  J. Chem. Phys. 144, 074109 (2016).
  * (109) Y. Ke, R. Borrelli, and M. Thoss,  arXiv:2202.10273 (2022).


11institutetext: Department of Theoretical Physics,
Moscow State University, 119991 Moscow, Russia
11email<EMAIL_ADDRESS><EMAIL_ADDRESS>
# Wave packet treatment of neutrino flavour and spin oscillations in galactic
and extragalactic magnetic fields

Artem Popov  Alexander Studenikin

###### Abstract

We consider neutrino flavour and spin oscillations in a magnetic field using
formalism of wave packets. Decoherence effects due to neutrino wave packets
separation are studied. The considered effects are especially important for
describing astrophysical neutrino oscillations, since they propagate on
kiloparsec scale and bigger. The obtained results are of interest for neutrino
telescopes IceCube, Baikal-GVD and KM3NeT, and also can be applied for
description of supernovae neutrino oscillations effects would be detected by
JUNO and Hyper-Kamiokande.

###### Keywords:

Neutrino magnetic moment Neutrino oscillations Astrophysics.

##  1 Introduction

It is known that massive neutrinos possess nontrivial electromagnetic
properties, in particular nonzero anomalous magnetic moments. For a review of
theory and experiment of neutrino electromagnetic properties see [1] and [2].
Currently, the best upper bound on neutrino effective magnetic moment obtained
by the XENONnT experiment and is on the level $\mu_{\nu}\sim 6.4\times
10^{-12}\mu_{B}$ [3, 4]. Interaction of neutrino magnetic moments with a
magnetic field induces neutrino spin and spin-flavour precession, i.e. a
phenomena in which neutrinos are converted to right-handed states. In
particular, spin precession can be induced by neutrino interaction with cosmic
magnetic field, that is of order of $\mu$G in our Galaxy [5].

To study neutrinos propagation through an environment of cosmological scales
it is necessary to employ the wave packet approach that accounts for
decoherence of neutrino oscillations. Previously, the wave packet formalism
was developed for the case of vacuum neutrino oscillations within relativistic
quantum mechanics (see [6] for a review) and quantum field theory [7]
approaches, for neutrino oscillations in matter [9, 10] and collective
neutrino oscillations [11]. In this paper we extend the wave packet approach
to the case of neutrino propagation in a magnetic field.

##  2 Neutrino oscillations in a magnetic field in the wave packet formalism

Neutrino oscillations in a uniform magnetic field are described by the
following modified Dirac equation

$(i\gamma^{\mu}\partial_{\mu}-m_{i})\nu_{i}(x)-\mu_{i}\bm{\Sigma}\bm{B}\nu_{i}(x)=0,$
(1)

where $i=1,2,3$ enumerates neutrino massive states. In [12, 13], equation (1)
was solved in the plane wave approximation. To account for potentially
important decoherence effects in neutrino oscillations at long distances, we
solve (1) assuming that at initial moment $t=0$ neutrino wave function Fourier
transform is described by a Gaussian wave packet

$\nu_{i}(p,0)=f_{i}(p,p_{0})u^{-}_{i}(p),\;\;\;f_{i}(p,p_{0})=\frac{1}{(2\pi\sigma_{p}^{2})^{1/4}}\exp\left(-\frac{(p-p_{0})^{2}}{4\sigma_{p}^{2}}\right),$
(2)

where $p_{0}$ is the average wave packet momentum, $\sigma_{p}$ is wave packet
width in momentum space and $u^{-}_{i}$ is a left-handed solution of the Dirac
equation in vacuum. For clarity, in this paper we consider one dimensional
wave packets. Three-dimensional consideration leads to different wave packet
spreading times in longitudinal and transversal directions, but does not
significantly affect neutrino oscillations patterns [7, 8].

The dispersion relation for neutrinos interacting with a magnetic field is
given by the following expression [12]

$E_{i}^{s}(p)=\pm\sqrt{m_{i}^{2}+p^{2}+\mu_{i}^{2}B^{2}+2s\mu_{i}\sqrt{m_{i}^{2}B^{2}+p^{2}B_{\perp}^{2}}},$
(3)

where $s=\pm 1$. Here magnetic field $\bm{B}$ is decomposed into transverse
$\bm{B}_{\perp}$ and longitudinal $\bm{B}_{\parallel}$ components with respect
to the neutrino momentum.

For sufficiently small $\sigma_{p}$, (3) can be decomposed near average
momentum $p_{0}$

$E_{i}^{s}(p)=E(p_{0})+v_{i}^{s}(p_{0})(p-p_{0})+\mathcal{O}((p-p_{0})^{2}),$
(4)

where $\frac{\partial E_{i}^{s}(p)}{\partial p}\Big{|}_{p=p_{0}}$ are group
velocities. Using (3), we calculate neutrino wave packets group velocities
accounting for neutrino interaction with a magnetic field

$v_{i}^{s}(p_{0})=\frac{p_{0}}{E_{i}^{s}(p_{0})}\left(1+\frac{s\mu_{i}B^{2}_{\perp}}{\sqrt{m_{i}^{2}B^{2}+p^{2}_{0}B^{2}_{\perp}}}\right).$
(5)

Finally, for the probabilities of neutrino flavour and spin oscillations in a
magnetic field accounting for decoherence effects we get

$P_{\nu_{\alpha}^{h}\to\nu_{\beta}^{h^{\prime}}}(L)=\sum_{i,j}\sum_{s,s^{\prime}}U^{*}_{\beta
i}U_{\alpha i}U_{\beta j}U^{*}_{\alpha
j}C_{is}^{hh^{\prime}}C_{js^{\prime}}^{hh^{\prime}}\exp\Big{(}-i2\pi\frac{L}{L_{osc}^{ijss^{\prime}}}\Big{)}\exp\Big{(}-\frac{L^{2}}{(L_{coh}^{ijss^{\prime}})^{2}}\Big{)},$
(6)

where the corresponding oscillations and coherence lengths are introduced

$L_{osc}^{ijss^{\prime}}=\frac{\pi}{\omega_{ij}^{ss^{\prime}}},\;\;\;L_{coh}^{ijss^{\prime}}=\frac{2\sqrt{2}\sigma_{x}}{v_{i}^{s}-v_{j}^{s^{\prime}}}.$
(7)

The oscillations frequencies $\omega_{ij}^{s\sigma}$ are given by
$\omega_{ij}^{s\sigma}(p_{0})\approx\frac{\Delta
m^{2}_{ij}}{2p_{0}}+(\mu_{i}s-\mu_{j}\sigma)B_{\perp}$.

The coefficients $C^{hh^{\prime}}_{is}$ in (6) are given by

$C^{LL}_{is}\approx\frac{1}{2}+\mathcal{O}\Big{(}\frac{m_{i}^{2}}{p^{2}}\Big{)},\;\;\;C^{RL}_{is}\approx-\frac{s}{2}+\mathcal{O}\Big{(}\frac{m_{i}^{2}}{p^{2}}\Big{)}.$
(8)

The oscillations probabilities (6) generalize the expressions we obtained in
[13] and account for exponential damping of neutrino oscillations at large
distance due to wave packets separation. Using (5), we obtain the following
approximate expressions for the coherence lengths given that $p\gg
m_{i}\gg\mu_{i}B$:

$L_{coh}^{ijss}\approx\frac{4\sqrt{2}\sigma_{x}p^{2}}{\Delta
m_{ij}^{2}},\;\;\;L_{coh}^{ii-+}\sim\frac{\sigma_{x}p^{3}}{\mu_{i}Bm_{i}^{2}},\;\;\;L_{coh}^{ij-+}\approx
L_{coh}^{ijss}.$ (9)

Here the coherence lengths $L_{coh}^{ijss}$ and $L_{coh}^{ii-+}$ describe
damping of oscillations on vacuum lengths $L_{osc}^{ijss}=\frac{4\pi p}{\Delta
m^{2}_{ij}}$ and magnetic lengths
$L_{osc}^{ii-+}=\frac{\pi}{\mu_{i}B_{\perp}}$, correspondingly. Note that
unlike the coherence lengths of vacuum neutrino oscillations, the coherence
length $L_{coh}^{ii-+}$ is proportional to the cube of the neutrino momentum.

Interaction with a magnetic field can modify the flavour composition of
neutrino emanating from astrophysical objects. For example, consider high-
energy neutrinos originating from the Galactic centre. The flavour composition
at distance $L$ is given by

$r_{\alpha}(L)=\frac{\sum_{\alpha}r^{0}_{\alpha}P_{\alpha\beta}(L)}{\sum_{\alpha\beta}r^{0}_{\alpha}P_{\alpha\beta}(L)},$
(10)

where $r_{\beta}^{0}$ are the flavour ratios at the source and
$P_{\alpha\beta}=P_{\nu_{\alpha}^{L}\to\nu_{\beta}^{L}}$ are the probabilities
of flavour oscillations.

The mechanism of high-energy neutrinos production is presently unknown. We
study two cases considered in the literature: case of pion decay neutrinos and
case of muon decay neutrinos. The corresponding initial flavour compositions
are $r^{0}=(\frac{1}{3},\frac{2}{3},0)$ and $r^{0}=(0,1,0)$. We also consider
different values of neutrino magnetic moments from the interval
$(0.1,6.4)\times 10^{-12}\mu_{B}$, as well as different experimental values of
neutrino mixing parameters within their $3\sigma$ intervals. In Fig. 1 we show
the possible flavour ratios of neutrinos propagating from the Galactic centre
to a terrestrial neutrino telescope for different values of neutrino magnetic
moments $\mu_{1},\mu_{2},\mu_{3}$ and compare them to ones predicted by the
vacuum neutrino oscillations case. We assume that galactic magnetic field
strength is $B=0.2\;\mu$G, which corresponds to the average value of the
regular component of the Galactic magnetic field. Magnetic field of
approximately same strength is found in the Local Supercluster [14]. We find
that the neutrino interaction with the Galactic magnetic field can
significantly modify the observed flavour ratio. Thus, we conclude that
neutrino magnetic moments potentially can be probed by neutrino telescopes,
such as IceCube, Baikal-GVD and KM3NeT. The obtained results can be also
applied for description of supernovae neutrino oscillations effects would be
detected by JUNO and Hyper-Kamiokande.

This study was conducted within the Scientific Program of the National Center
for Physics and Mathematics, Section No. 8 (Stage 2023-2025).

Figure 1:  Flavour compositions of neutrinos coming from the Galactic centre
for the cases of pion decay and muon decay high-energy neutrinos production.

## References

  * [1] C. Giunti, A. Studenikin, Rev. Mod. Phys. 87, 531 (2015).
  * [2] A. Studenikin, PoS CORFU2021 (2022), 057.
  * [3] E. Aprile et al. [XENON Collaboration], Phys. Rev. Lett. 129 (2022) no.16, 161805.
  * [4] C. Giunti and C. A. Ternes, Phys. Rev. D 108 (2023) no.9, 095044.
  * [5] R. Beck, AIP Conf. Proc. 1085 (2009) no.1, 83-96.
  * [6] C. Giunti, Found. Phys. Lett. 17 (2004), 103-124.
  * [7] D. V. Naumov, V. A. Naumov, Phys. Part. Nucl. 51 (2020) no.1, 1-106.
  * [8] D. V. Naumov, Phys. Part. Nucl. Lett. 10 (2013), 642-650.
  * [9] J. T. Peltoniemi and V. Sipilainen, JHEP 06 (2000), 011.
  * [10] J. Kersten and A. Y. Smirnov, Eur. Phys. J. C 76 (2016) no.6, 339.
  * [11] E. Akhmedov, J. Kopp and M. Lindner, JCAP 09 (2017), 017.
  * [12] A. Popov and A. Studenikin, Eur. Phys. J. C 79 (2019) no.2, 144.
  * [13] A. Lichkunov, A. Popov and A. Studenikin, [arXiv:2207.12285 [hep-ph]].
  * [14] J.P. Vallée, The Astronomical Journal (2002), 124, 1322 - 1327.