# Ion-trap electrode preparation with Ne+ bombardment

K. S. McKay  D. A. Hite  Y. Colombe  Present Address: Institut für
Experimentalphysik, Universität Innsbruck, Technikerstr. 25, A-6020 Innsbruck,
Austria  R. Jördens  A. C. Wilson  D. H. Slichter  D. T. C. Allcock  D.
Leibfried  D. J. Wineland  D. P. Pappas  National Institute of Standards and
Technology, 325 Broadway, Boulder, Colorado 80305

###### Abstract

We describe an ex-situ surface-cleaning procedure that is shown to reduce
motional heating from ion-trap electrodes. This precleaning treatment, to be
implemented immediately before the final assembly and vacuum processing of ion
traps, removes surface contaminants remaining after the electrode-fabrication
process. We incorporate a multi-angle ion-bombardment treatment intended to
clean the electrode surfaces and interelectrode gaps of microfabricated traps.
This procedure helps to minimize redeposition in the gaps between electrodes
that can cause electrical shorts. We report heating rates in a stylus-type ion
trap prepared in this way that are lower by one order of magnitude compared to
a similar untreated stylus-type trap using the same experimental setup.

Electric-field noise from ion-trap surfaces can cause motional heating and
limit the coherence time of quantum operations to less than that required for
fault-tolerant quantum information processing faulttolerant . Heating rates
tend to decrease with increasing ion-electrode distance d, however smaller ion
traps are desired for scalability and faster quantum gates. The motional
heating rate $\dot{\overline{n}}$, defined as the time rate of change of the
average motional-state occupation number, is related to the spectral density
of the electric-field noise $S_{E}$ by

$S_{E}(\omega)={4m\hbar\omega\over q^{2}}\dot{\overline{n}}(\omega),$ (1)

where $\omega/2\pi$ is the motional frequency of the ion in the trap, m is the
ion’s mass, q is its charge, and $\hbar$ is Planck’s constant divided by
2$\pi$ Turchette00 . The electric-field noise that causes motional heating
originates from thermally driven processes Deslauriers06 ; Labaziewicz08a ,
and can be significantly reduced by operating at cryogenic temperatures. For
example, relative to room temperature, reductions in electric-field noise of
22 dB Brown11 ; Hite12 , 19 dB Chiaverini14 , and 34 dB Labaziewicz08b  have
been reported, corresponding to cryogenic-electrode heating rates as low as 70
quanta/s (9Be+, $\omega/2\pi$ = 2.3 MHz), 6 quanta/s (88Sr+, 1.3 MHz), and 2
quanta/s (88Sr+, 1.0 MHz), at ion-electrode distances of 40, 50, and 75
$\mu$m, respectively. In addition, electric-field noise from room-temperature
surface-electrode traps has been reduced by approximately 22 dB following in-
situ treatments of ion bombardment, corresponding to heating rates as low as
43 quanta/s (9Be+, 3.6 MHz, d = 40 $\mu$m) Hite12 , and 3.8 quanta/s (40Ca+,
1.0 MHz, d = 100 $\mu$m) Daniilidis13 . For comparison, heating rate data and
other experimental parameters for some traps are summarized in Table 1, and
the corresponding normalized electric-field noise spectral densities $\omega
S_{E}$ are plotted in Fig. 1 as a function of d Hite12 .

Table 1: Heating rates $\dot{\overline{n}}$ from electric-field noise in some microfabricated ion traps of various sizes. Electrode treatments and experimental parameters are listed to compare the “frequency-normalized” electric-field noise spectral densities $\omega S_{E}$ for various ion-electrode distances d and electrode temperatures. The data are grouped by electrode temperature and the following electrode treatments: untreated (UT), apart from typical microfabrication cleaning with solvents, in-situ treated (IT), and pretreated (PT), described in this work. |  | $\dot{\overline{n}}$ | d | $\omega/2\pi$ |  | Electrode | $\omega S_{E}$ |
---|---|---|---|---|---|---|---|---
| Treatment | (s-1) | ($\mu$m) | (MHz) | Ion | Material | (10-6 V2/m2) | Reference
1 | UT | 600 | 40 | 4.6 |  25Mg+ | Au | 340 | Epstein07
2 | (300 K) | 500 | 50 | 1.3 |  88Sr+ | Au | 83 | Chiaverini14
3 |  | 410 | 63 | 1.4 |  40Ca+ | Al | 33 | Doret12
4 | UT | 70 | 40 | 2.3 |  9Be+ | Au | 3.6 | Brown11
5 | (5 K) | 6 | 50 | 1.3 |  88Sr+ | Au | 1.0 | Chiaverini14
6 |  | 2.1 | 100 | 0.8 |  88Sr+ | Nb | 0.025 | Wang10
7 | IT | 43 | 40 | 3.6 |  9Be+ | Au | 5.4 | Hite12
8 | (300 K) | 3.8 | 100 | 1.0 |  40Ca+ | Cu/Al | 0.15 | Daniilidis13
9 |  UT (300 K) | 16,000 | 40 | 3.6 |  9Be+ | Au | 2100 | Hite13
10 |  PT (300 K) | 930 | 40 | 3.6 |  9Be+ | Au | 120 | Hite13
11 |  UT (300 K) | 387 | 62 | 4.0 |  25Mg+ | Au | 170 | Arrington13
12 |  PT (300 K) | 30 | 60 | 4.4 |  25Mg+ | Au | 14 | This Work
Figure 1: Normalized electric-field noise spectral densities in selected low-
noise microfabricated traps of various sizes and electrode materials. Since
$S_{E}$ exhibits an approximate $1/\omega$ dependence in many ion-trap
experiments, we normalize data taken at different frequencies $\omega/2\pi$
onto a common curve by plotting $\omega S_{E}$. Dotted lines indicating a
1/$\textit{d}^{4}$ distance scaling are plotted as guides to the eye. Open
symbols indicate data from cryogenic traps with electrodes cooled to liquid
helium temperatures. Electrode materials are indicated in the legend.

Despite this progress in reducing motional heating rates through cryogenic
operation or in-situ treatments, electric-field noise in these traps typically
remains well above estimates for thermal voltage noise (Johnson noise) from
the electrode circuitry, and therefore is commonly referred to as “anomalous
heating”. While this anomalous heating is thought to be related to surface
contamination on the trap electrodes, a deeper understanding of its physical
origin will help to solve this problem. In an experiment that tested the
effects of contamination from air exposure and vacuum baking on a surface-
electrode trap that was previously treated with in-situ Ar+ bombardment Hite12
, heating rates were measured to increase from 50 to 930 quanta/s (9Be+, 3.6
MHz, d = 40 $\mu$m), after a 48-h exposure to air followed by vacuum baking at
180 ∘C for approximately 100 h (see Fig. 5(c) in Ref. Hite13 ). This
corresponds to a 12 dB reduction in the electric-field noise spectral density
as compared to that measured before the Ar+-cleaning treatments (Table I, # 9
and 10). These data indicate a potential benefit resulting from a precleaning
treatment prior to trap assembly. While it is believed that in-situ treatments
are required for full removal of the dominant source of this anomalous heating
Hite13 , it would be advantageous to implement an ex-situ precleaning
treatment to reduce heating rates, thereby enabling new experiments, operating
at either room or cryogenic temperatures, without the need for in-situ
cleaning capabilities. It may be that in-situ cleaning treatments can be
combined with cryogenic cooling to reduce heating rates further; however,
implementing in-situ ion-bombardment treatments in a cryogenic ion-trap setup
is challenging, given the necessity to provide optical access, and to
incorporate cryogenic, vacuum, and gas-handling components.

Ion bombardment tends to redeposit sputtered metal into the interelectrode
gaps. This redeposited metal can cause electrical shorts in regions of the
gaps that are shielded by the electrodes from ion bombardment. Since the etch
rate from sputtering is observed to be greater than the rate of redeposition,
designing the trap-electrode and ion-beam geometry to eliminate the shielding
of the gaps can mitigate this problem. One could also optimize the dose of the
in-situ treatment required to achieve the desired heating rate by employing
treatment parameters (i.e., beam energy, mass, and incident angle) that
minimize the sputter yield of the electrode material and concurrently maximize
that of the surface contaminants Eckstein07 . In addition, one could perform
an ex-situ precleaning treatment just prior to final trap assembly, allowing
ion bombardment from multiple angles to remove any metal redeposited in the
gaps. There is evidence from Ref. Hite12  that higher sputtering energies may
be necessary to achieve very low heating rates at room temperature. This
suggests a requirement for relatively high-energy surface modifications, e.g.,
removal of constituents with relatively low sputter yield footnote1 . However,
higher sputtering energies likely exacerbate the redeposition problem. In
addition to the potential benefit to traps without in-situ cleaning
capabilities, ex-situ precleaning treatments may also help to reduce this
redeposition problem by allowing a gentler in-situ treatment.

In this work, we used a multi-angle precleaning procedure using high sputter
doses to remove contaminants while simultaneously removing redeposited metal
from the gaps between electrodes. We compare heating rates in two similar (but
not identical) stylus-type ion traps, one with and one without the precleaning
treatment (Table I, # 12 and 11, resp.), and find a heating rate in the
precleaned trap lower by an order of magnitude. Making use of in-situ Auger
electron spectroscopy, we correlate this reduction in electric-field noise
with the reduction in surface contaminants due to the precleaning procedure.

Using a microfabricated stylus trap, described in detail in Ref. Arrington13
(Table I, # 11), we reported an untreated trap-electrode heating rate of 387
$\pm$ 15 quanta/s for $d=$ 62 $\mu$m above an 80-$\mu$m diameter Au stylus
post (Fig. 1, hexagon symbol). The corresponding electric-field noise, which
is well above estimates for thermal voltage (Johnson) noise, is thought to
arise from dynamical processes of contaminants on the surface that cause
contact-potential fluctuations. These contaminants have many potential
sources, including residues from the fabrication processing, vacuum baking for
ultra-high vacuum (UHV) operation, or other adventitious contamination from
air exposure.

In order to test the dependence of ion motional heating on the amount of
surface contamination, and the benefits of a precleaning treatment, in this
work we prepared two nominally identical traps that were fabricated on the
same wafer. Both traps were precleaned as described next, exposed to air for
$\sim$ 24 h, and then vacuum-baked together in the ion-trap chamber. The
precleaning consisted of a series of treatments of Ne+ bombardment at three
different angles: 10 min at + 30∘ incident angle (with respect to the normal
of the surface, see Fig. 2), 10 min at - 30∘, and 20 min at 0∘, all treatments
using 2 kV and 30 $\mu$A/cm2 with a $\sim$ 3-mm focus diameter of the incident
ion beam, at a Ne pressure of $\sim$ 6 $\times$ 10-3 Pa footnote2 . The
initial angled treatments are intended to remove most of the contaminants away
from the surface and minimize any embedded contaminants. The final normal-
incidence treatment cleans shadowed redeposition of metal (and other
deleterious impurities) from the gaps in the electrodes. One of the traps,
shown in Fig. 2, was used to trap ions and compare heating rates with the
untreated trap of Ref. Arrington13 . The duplicate trap was mounted to a
sample transfer stage in the ion-trap chamber and subsequently transferred, in
UHV, to a surface-science chamber for the post-bake analysis.

Figure 2: Top-view optical micrograph of the as-fabricated stylus trap used in
this work before the precleaning treatment. The electroplated Au electrodes
are formed in two layers, a 10-$\mu$m thick base layer (grainy gold color) and
the 200-$\mu$m tall stylus posts (appearing white) protruding from the base
layer Arrington13 . The interelectrode gaps are 10 $\mu$m wide. The
orientation of the Ne+ sputter beam with respect to the electrodes is
indicated. The precleaning treatment angles are $\pm$ 30∘ and 0∘ with respect
to the normal of the surface (out of the page). 25Mg+ ions are trapped at a
height $d=$ 60 $\mu$m above the 60-$\mu$m diameter center post. Other than the
geometric differences on the chip, this trap is electrically the same as the
trap in Ref. Arrington13 , and uses the same experimental setup for the
heating rate measurements.

Prior to the precleaning treatment or vacuum baking, Auger electron
spectroscopy (AES) of the untreated Au electrodes on the as-fabricated
duplicate trap chip indicates the presence of carbonaceous contamination (1 -
2 carbon monolayers (ML) thick), previously correlated with anomalous heating
Hite12 , along with small amounts of other contaminants (S, O, Fe) from the
electroplating bath used for gold deposition (Fig. 3-a). Immediately after
applying the Ne+ precleaning treatment to the duplicate trap-electrode
surfaces, Auger spectra collected in situ signify a clean Au surface free of
these contaminants (Fig. 3-b). Before vacuum baking, air exposure of these
precleaned surfaces for $\sim$ 24 h typically results in a reacquiring of less
than 0.5 ML of oxygen-free adventitious carbon. The level of contamination on
the actual ion-trap electrodes was also monitored using AES during the steps
of the precleaning treatment to confirm comparable results. Upon vacuum baking
for UHV processing of the assembled ion-trap chamber, the thickness of the
carbonaceous layer grows slightly, presumably from hydrocarbons desorbing from
the vacuum-chamber walls and other surfaces. Nevertheless, vacuum baking of
these precleaned and air-exposed surfaces typically results in less surface
contamination than the untreated, as-fabricated electrodes. The approximate
carbon coverage on the precleaned, air-exposed, then vacuum-baked duplicate
trap was determined from in-situ AES to be $\sim$ 0.6 ML (Fig. 3-c). Since the
actual trap chip and the duplicate trap chip were exposed to air and vacuum
baked together in the assembled ion-trap chamber, we assume that the actual
trap chip has the same level of contamination as the duplicate trap.

Figure 3: Surface analysis of the electroplated-Au electrodes on the duplicate
trap chip, (a) before and (b) immediately after the precleaning treatment.
This companion sample was then exposed to identical environments as the actual
trap chip used to measure heating rates, namely air exposure for 24 h and
vacuum baking at 180 ∘C for approximately 100 h (together in the same
chamber). Spectrum (c) was acquired from the duplicate trap chip in situ,
after operating the Mg source used in the ion-trap chamber and then UHV
transfer to the surface analysis chamber. No Mg was detected on the duplicate
trap chip. We note the presence of a signal from thallium, which is driven to
the surface from the bulk material upon baking of the electroplated, soft-gold
electrodes. Traces are offset vertically for clarity.

The lowest motional heating rate of the precleaned ion trap was measured to be
30 $\pm$ 2 quanta/s (Table I, # 12), for $d=$ 60 $\mu$m above a 60-$\mu$m
diameter Au stylus post (Fig. 1, star symbol). This corresponds to an 11 dB
improvement in electric-field noise compared to the untreated stylus trap of
Ref. Arrington13  (Table I, # 11). Moreover, this is consistent with the level
of improvement (12 dB) observed in the air-exposure/vacuum-baking test of Ref.
Hite13  (Table I, # 9 and 10), discussed above. Therefore, this suggests that
the significantly lower heating rate reported here is due to surface
modifications made during the precleaning treatment, and that the electric-
field noise from surfaces sensed by trapped ions is directly related to the
coverage of adsorbates on the surface.

There are slight differences between the untreated trap of Ref. Arrington13
and the precleaned trap, whose heating rates are being compared in this study.
Most of these differences are geometric design modifications made to the
stylus posts, primarily to minimize shadowing of the gaps near the posts
during sputter treatments that can cause electrical shorting. Beside
increasing the height of the stylus posts to 200 $\mu$m above the chip
(compared to the 150-$\mu$m posts used in Ref. Arrington13 ), the diameter of
the center post was reduced to 60 $\mu$m (Fig. 2), from the Ref. Arrington13
design of 80 $\mu$m. Due to these geometric changes, simulations of the trap
indicate that the ion is trapped approximately 60 $\mu$m above the center
post, as compared to 62 $\mu$m in the trap in Ref. Arrington13 . The different
geometries of the untreated and precleaned traps should correspond to slightly
different electric-field noise spectral densities at the ion for a given
voltage noise on the electrode surfaces. Assuming either a 1/$\textit{d}^{4}$
or 1/$\textit{d}^{2}$ distance scaling for the noise, we estimate that the
precleaned trap should have 10$\%$ or 15$\%$ less noise, respectively, due to
geometrical considerations. Therefore, we conclude that these geometric
differences cannot account for the observed reduction in heating rate.

We checked that the heating rate in the precleaned trap was not limited by
extraneous technical noise. This can arise from voltage noise from the
computer-controlled digital-to-analog converters and amplifiers that apply
static potentials to the trap. Therefore, to exclude this possibility, we
applied static (dc) potentials generated by batteries. We also checked for ion
state de-pumping during the wait period in the heating rate measurements by
using a fast shutter to eliminate possible stray light from the acousto-
optical switches. During these tests, we observed no measurable differences in
the heating rate.

Preliminary tests of the cleaning procedure on other electrode materials
indicate that, in addition to gold, Cu and Al electrodes may also benefit from
such a precleaning treatment. However, for Nb electrodes, our Auger analysis
indicates the presence of a very low-sputter-yield refractory-metal carbide
following the precleaning treatment. Niobium carbide forms when the native
oxide of Nb is covered with adventitious carbon and undergoes ion bombardment
Grundner80 . This may preclude the use of ion-bombardment precleaning
treatments on Nb electrodes.

This research was funded by the Office of the Director of National
Intelligence (ODNI), Intelligence Advanced Research Projects Activity (IARPA).
All statements of fact, opinion, or conclusions contained herein are those of
the authors and should not be construed as representing the official views or
policies of IARPA or the ODNI. We also acknowledge the support of ONR and the
NIST Quantum Information Program. We thank the Sandia National Laboratories’
Metal Micromachining Program and CAMD for fabricating the stylus traps used in
this work. This article is a contribution of the U.S. Government, and is not
subject to U.S. copyright.

## References

  * (1) J. Preskill, Proc. R. Soc. London, Ser. A 454, 385 (1998); A. S$\o$rensen and K. M$\o$lmer, Phys. Rev. A 62, 022311 (2000); E. Knill, Nature 463, 441 (2010).
  * (2) Q. A. Turchette, et al., Phys. Rev. A 61, 063418 (2000).
  * (3) L. Deslauriers, S. Olmschenk, D. Stick, W. K. Hensinger, J. Sterk, and C. Monroe, Phys. Rev. Lett. 97, 103007 (2006).
  * (4) J. Labaziewicz, Y. Ge, P. Antohi, D. Leibrandt, K. R. Brown, and I. L. Chuang, Phys. Rev. Lett. 100, 013001 (2008).
  * (5) K. R. Brown, C. Ospelkaus, Y. Colombe, A. C. Wilson, D. Leibfried and D. J. Wineland, Nature 471, 196 (2011).
  * (6) D. A. Hite, et al., Phys. Rev. Lett. 109, 103001 (2012).
  * (7) J. Chiaverini and J. M. Sage, Phys. Rev. A 89, 012318 (2014).
  * (8) J. Labaziewicz, Y. Ge, D. R. Leibrandt, S. X. Wang, R. Shewmon, and I. L. Chuang, Phys. Rev. Lett. 101, 180602 (2008).
  * (9) N. Daniilidis, S. Gerber, G. Bolloten, M. Ramm, A. Ransford, E. Ulin-Avila, I. Talukdar, and H. Häffner, ArXiv 1307.7194v1 (2013).
  * (10) D. A. Hite, Y. Colombe, A. C. Wilson,D. T. C. Allcock, D. Leibfried, D. J. Wineland, and D. P. Pappas, MRS Bull. 38, 826 (2013).
  * (11) W. Eckstein, in Sputtering by Particle Bombardment, Topics in Applied Physics Vol. 110, edited by R. Behrisch and W. Eckstein (Springer, Berlin, 2007), p.33.
  * (12) Generally, use of Ar+ results in higher sputter yields, however the lower sputter yields for Au using Ne+ are desired to mitigate redeposition. See Ref. Eckstein07 .
  * (13) C. L. Arrington, et al., Rev. Sci. Instrum. 84, 085001 (2013).
  * (14) Base pressures should be low enough that the adsorption rate of residual gas species on the surface is less than the sputter rate of the impinging ion beam, e.g., to avoid oxide or nitride formation Eckstein07 .
  * (15) M. Grundner and J. Halbritter, J. Appl. Phys. 51, 397 (1980).
  * (16) R. J. Epstein, et al., Phys. Rev. A 76, 033411 (2007).
  * (17) S. C. Doret et al., New J. Phys. 14, 073012 (2012).
  * (18) S. X. Wang, Y. F. Ge, J. Labaziewicz, E. Dauler, K. Berggren, I. L. Chuang, Appl. Phys. Lett. 97, 244102 (2010).


# Impact of Secondary User Interference
on Primary Network in Cognitive Radio Systems

Amit Kachroo and Sabit Ekin  School of Electrical and Computer Engineering
Oklahoma State University, Stillwater, OK, USA
Email: amit.kachroo<EMAIL_ADDRESS>
###### Abstract

Most of the research in cognitive radio field is primarily focused on finding
and improving secondary user (SU) performance parameters such as bit error
rate, outage probability and capacity etc. Less attention is being paid
towards the other side of the network that is the primary network which is
under interference from SU. Also, it is the primary user (PU) that decides
upon the interference temperature constraint for power adaptation to maintain
a certain level of quality of service while providing access to SUs. However,
given the random nature of wireless communication, interference temperature
can be regulated dynamically to overcome the bottlenecks in entire network
performance. In order to do so, we need to analyze the primary network
carefully. This study tries to fill this gap by analytically finding the
closed form theoretical expressions for signal to interference and noise ratio
(SINR), mean SINR, instantaneous capacity, mean capacity and outage
probability of PU, while taking peak transmit power adaptation at SU into
picture. Furthermore, the expressions generated are validated with the
simulation results and it is found that our theoretical derivations are in
perfect accord with the simulation outcomes.

###### Index Terms:

Cognitive Radio Network, Interference Temperature, Mean Capacity, SINR, Outage
Probability

##  I Introduction

In cognitive radio network, a secondary user (SU) is allowed to access the
primary user (PU) spectrum completely if the available spectrum is not used by
PU (interweave) or concurrently (underlay) with PU. The concurrent
transmission is allowed if and only if the SU maintains a certain power
threshold constraint known as interference temperature [1, 2]. Most of the
studies [1, 2, 3, 4, 5, 6] that involve this interference temperature model
utilize either peak or average transmit power adaptation for the purpose of
analyzing or improving the performance of secondary network. However, the
impact of SU interference on primary network is let off completely. Therefore,
to asses the performance and other quality of service (QoS) parameters, closed
form expression need to be derived and validated. This study focuses on the
mathematical foundation to derive these necessary performance expressions.
This theoretical analysis is done by first considering interference from a
single SU and then extended to the case of interference from multiple SUs on
primary network. As per our knowledge, this is the first paper to analytically
analyze the effect of SU interference on primary network. The contribution of
this paper can be summarized as follows:

  * •

Probability distribution function (PDF) and cumulative distribution function
(CDF) expressions for noise plus interference, signal to noise and
interference ratio (SINR) are derived for both cases of interference from a
single SU and multiple SUs on PU network.

  * •

Closed form mean SINR expression, mean capacity and outage probability
expressions are derived.

  * •

The expression generated above are validated with simulation results to show
the accuracy of the theoretical expressions.

The rest of the paper is organized as follows. In Section II, system model
with underlying assumptions is described in detail. Section III presents the
extensive theoretical analysis and comparisons with simulation results and
finally, conclusions are given in Section IV.

##  II System Model

The cognitive radio network that is under consideration in this work is shown
in Fig. 1. It consists of $n$ SUs (SU network) and $k$ PUs with corresponding
SBS and PBS. Since we are analyzing the effect of SU on PU network, we don’t
need to consider the channel gain between PU and SBS, and also we don’t need
to consider the interference among PUs because of orthogonal resource
allocation between them.

Figure 1: Underlay cognitive network with $n$-SUs sharing the spectrum with PU
network of $k$-PUs. The channel power gain between any SU-i ($i^{th}$ user)
and PBS is denoted by $\alpha_{i}$, between any SU-i ($i^{th}$ user) and SBS
by $\beta_{i}$ and between $i^{th}$-PU and PBS by $\gamma_{i}$.

In addition, since the channel fading is assumed to be Rayleigh distributed,
the channel power gains follow an exponential distribution. Considering the
peak power adaptation [6, 4] for the case of interference from single SU on
primary network, the secondary transmit power is given by:

$P^{tx}=\text{min}\left\\{p,\frac{q}{\alpha}\right\\}$ (1)

In the following sections, for theoretical analysis purpose, we assume that
there are $n$ SUs that form the underlay cognitive network with primary user,
where $n=\\{1,2,3\dots,n\\}$. Furthermore, the thermal additive white Gaussian
noise (AWGN) in the network is assumed to have circularly symmetric complex
Gaussian distribution with zero mean and variance as $\sigma^{2}$, i.e.,
$\mathcal{CN}(0,\sigma^{2})$.

##  III Theoretical Analysis

In this section, theoretical expressions for PU performance parameters with
interference from a single SU will be derived first and then extended to the
case of interference from multiple SUs. The interference observed at primary
receiver because of a single SU and multiple SUs with peak power adaptation
will be given as:

$\begin{split}I_{\text{single}}&=\alpha P_{sec}=\text{min}\\{\alpha p,q\\},\\\
I_{\text{multi}}&=\text{min}\left\\{\sum_{i=1}^{n}\alpha_{i}p,q\right\\}.\end{split}$
(2)

Eq. (2) represents a minimum of a random variable and a constant111 For
illustration purposes, the value’s of peak power $p$ and interference
temperature $q$ in this paper are chosen to be in linear scale. However, the
expressions derived in this paper hold for any value of $p$ and $q$ for any
scale.. From the theory of mixed random variables[7, 8, 9], a constant $c$ can
be modelled as a random variable with PDF equal to $\delta(x-c)$ and CDF equal
to $H(x-c)$, where $H(x)$ is a Heaviside function and $\delta(x)$ is a Dirac
Delta function. So with interference from a single SU, the CDF and PDF of
minimum of two independent random variables is then given by:

${}\begin{split}F_{I}(x)&=F_{\alpha p_{0}}(x)+F_{q}(x)-F_{\alpha
p_{0}}(x)F_{q}(x),\\\ F_{I}(x)&=1-e^{-\frac{\lambda
x}{p}}(1-H(x-q)).\end{split}$ (3)

On including noise $\mathcal{CN}(0,\sigma^{2})$, the CDF of noise plus
interference will be then,

$F_{IN}(x)=1-e^{-\frac{\lambda(x-\sigma^{2})}{p}}(1-H(x-\sigma^{2}-q)).$ (4)

Correspondingly, the PDF of interference and noise is given by differentiating
the CDF with respect to noise and interference variable $x$, i.e.

$\begin{split}f_{IN}(x)&=\frac{\lambda}{p}e^{-\frac{\lambda(x-\sigma^{2})}{p}}\Bigl{(}1-H(x-\sigma^{2}-q)\\\
&~{}+\frac{p}{\lambda}\delta(x-\sigma^{2}-q)\Bigr{)},\ \ \forall\
\sigma^{2}\leq x\leq\infty.\end{split}$ (5)

For the case of interference from multiple SUs, the distribution of of
interference given in Eq. (2) will follow Gamma distribution222 The
distribution of sum of independent exponential random variables with the same
rate parameters follows Gamma distribution. Also, to distinguish between the
channel $\gamma$ between PU and PBS, the Gamma distribution is denoted as
$\bar{\gamma}$ in this study., $f_{\bar{\gamma}}(x,\kappa,\theta)$, where
$\kappa$ and $\theta$ represent the shape and rate parameter. The PDF and CDF
of Gamma distribution is given as

$\begin{split}f_{\bar{\gamma}}(x)&=\sum_{i=1}^{n}\alpha_{i}p=\bar{\gamma}\left(x,n,\frac{\lambda_{2}}{p}\right)=\bar{\gamma}\left(x,n,\bar{\lambda}\right),\\\
&=\frac{\bar{\lambda}^{n}x^{n-1}}{\Gamma(n,0)}e^{-\bar{\lambda}x},\ \forall\
\left\\{x\geq 0,n>0,\bar{\lambda}>0\right\\},\\\
F_{\bar{\gamma}}(x)&=1-\frac{\Gamma\left(n,\bar{\lambda}x\right)}{\Gamma(n,0)},\
\forall\ \left\\{x\geq 0,n>0,\bar{\lambda}>0\right\\},\end{split}$ (6)

where $n$ is the total number of SUs in the underlay network,
$\bar{\lambda}=\lambda_{2}/p$ is the scaled rate parameter between SU and PBS
and $\Gamma(a,x)$ is an incomplete gamma function defined as:

$\Gamma(a,x)=\int_{a}^{\infty}t^{a-1}e^{-t}dt,\ \forall\ a>0,x\geq 0.$

Following the same mathematical approach that was used in single SU case, the
distribution of noise plus interference in multiple SUs case is then derived
as:

$\begin{split}F^{m}_{NI}(x)&=1-\frac{\Gamma\left(n,\bar{\lambda}(x-\sigma^{2})\right)}{\Gamma(n,0)}+H(x-\sigma^{2}-q)\\\
&~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}\times\frac{\Gamma\left(n,\bar{\lambda}(x-\sigma^{2})\right)}{\Gamma(n,0)},\\\
f_{NI}^{m}(x)&=\frac{\Gamma\left(n,\bar{\lambda}(x-\sigma^{2})\right)}{\Gamma(n,0)}\delta(x-\sigma^{2}-q)+\bar{\lambda}^{n}\\\
&\times\frac{(x-\sigma^{2})^{n-1}}{\Gamma(n,0)}e^{-\bar{\lambda}(x-\sigma^{2})}[1-H(x-\sigma^{2}-q)],\end{split}$
(7)

where $\sigma^{2}$ is the $\mathcal{CN}(0,\sigma^{2})$. Fig. 2 and Fig. 3
plots the CDF and PDF for the theoretical expression (Eq. (7) ) with the
simulation result of $p>q$ and $q<p$ for different SU densities of $n=1,2,3$.

Figure 2: PDF and CDF of noise and interference for different number of SUs
($n=1,2,3$), when $p>q$, where $p=4$, $q=2$ and $\sigma^{2}=1$ with support
region from $\sigma^{2}\leq x\leq\infty$. Figure 3: PDF and CDF of noise and
interference for different number of SUs ($n=1,2,3$) when $p<q$, where $p=2$
$q=4$ and $\sigma^{2}=1$ with support region from $\sigma^{2}\leq
x\leq\infty$.

###  III-A Instantaneous SINR

The instantaneous SINR at PBS considering the system model (Fig. 1) is given
by:

$\text{SINR}=\frac{\gamma p}{\sigma^{2}+I},$ (8)

where $I$ is the interference from SUs given by Eq. (2). The distribution of
numerator is a scaled exponential distribution and the distribution of
denominator is already derived in the previous section (Eq. (5) and Eq. (7)).
Therefore, the PDF of ratio of two independent random variables [10] i.e.
$z=x/y$, where $x=\gamma p$ and $y=\sigma^{2}+I$ will be given as

$\begin{split}f_{z}(z)&=\int_{\sigma^{2}}^{\infty}y\cdot
f_{x,y}(yz,y)dy=\int_{\sigma^{2}}^{\infty}y\cdot f_{x}(yz)f_{y}(y)dy\\\
&~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}\forall\
y\geq 0.\end{split}$ (9)

For the interference from a single SU user, the SINR distribution will be as
follows:

$\begin{split}f_{z}(z)&=\int_{\sigma^{2}}^{\infty}y\cdot\frac{\lambda_{1}e^{\frac{-\lambda_{1}yz}{p}}}{p}\frac{\lambda_{2}}{p}e^{-\frac{\lambda_{2}(y-\sigma^{2})}{p}},\\\
&\times\Bigl{\\{}1-H(y-\sigma^{2}-q)+\frac{p}{\lambda_{2}}\delta(y-\sigma^{2}-q)\Bigr{\\}}dy\\\
&=\frac{\lambda_{1}\lambda_{2}}{p^{2}}e^{\frac{\lambda_{2}\sigma^{2}}{p}}\left\\{\int_{\sigma^{2}}^{\infty}y\cdot
e^{\frac{-y(\lambda_{1}z+\lambda_{2})}{p}}dy-\int_{\sigma^{2}+q}^{\infty}y\right.\\\
&\times
e^{\frac{-y(\lambda_{1}z+\lambda_{2})}{p}}dy\left.+\frac{p}{\lambda_{2}}(\sigma^{2}+q)e^{\frac{-(\sigma^{2}+q)(\lambda_{1}z+\lambda_{2})}{p}}\right\\}.\end{split}$

By using integration by parts and on further simplification, the PDF is
reduced to:

$\begin{split}f_{z}(z)&=\frac{\lambda_{1}\lambda_{2}}{\Lambda
p}\left\\{e^{\frac{-\sigma^{2}\lambda_{1}z}{p}}\left(\sigma^{2}+\frac{p}{\Lambda}\right)\right.\\\
&\left.+\
e^{-\frac{\sigma^{2}\lambda_{1}z+q\Lambda}{p}}\left(\frac{(\sigma^{2}+q)\lambda_{1}z}{\lambda_{2}}-\frac{p}{\Lambda}\right)\right\\},\end{split}$

where $\Lambda$ is the scaled and shifted random variable version
333$\lambda_{1}$ is the channel rate parameter between PU and PBS, whereas
$\lambda_{2}$ is the channel rate parameter between SU and PBS. of $z$ given
by $\Lambda=\lambda_{1}+\lambda_{2}z$ . Under the scenario of
$\lambda_{1}=\lambda_{2}=1$, with AWGN as $\mathcal{CN}(0,\sigma^{2}=1)$, the
PDF can be further simplified to

$\begin{split}f_{z}(z)&=\frac{1}{p(z+1)}\left\\{e^{\frac{-z}{p}}\left(1+\frac{p}{z+1}\right)+e^{-\frac{z+q(z+1)}{p}}\right.\\\
&~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}\left.\times\left((1+q)z-\frac{p}{z+1}\right)\right\\}.\end{split}$
(10)

Following the same analytical framework used for the single SU case, the
distribution with interference from multiple SUs will be:

$\begin{split}f_{z}^{m}(z)&=\int_{\sigma^{2}}^{\infty}y\bar{\lambda_{1}}e^{-\bar{\lambda_{1}}yz}\Big{\\{}f_{\bar{\gamma}}(y-\sigma^{2},n,\bar{\lambda})+\delta(y-\sigma^{2}-q)\\\
&~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}\times(1-F_{\bar{\gamma}}(y-\sigma^{2},n,\bar{\lambda}))-f_{\bar{\gamma}}(y-\sigma^{2},n,\bar{\lambda})\\\
&~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}\times
H(y-\sigma^{2}-q)\Big{\\}}dy.\\\ \end{split}$

which on further evaluation and simplification reduces to

$\begin{split}f_{z}^{m}(z)&=\bar{\lambda_{1}}\bar{\lambda}^{n}e^{-\sigma^{2}\bar{\lambda_{1}}z}\Theta^{-1-n}\Bigg{[}n+\bigg{(}\sigma^{2}\Theta\\\
&\times\bigg{[}1-\frac{\Gamma(n,q\Theta)}{\Gamma(n,0)}\bigg{]}\bigg{)}-\frac{\Gamma(n+1,q\Theta)}{\Gamma(n,0)}\Bigg{]}\\\
&~{}~{}~{}~{}~{}~{}~{}~{}+\bar{\lambda_{1}}(\sigma^{2}+q)\frac{\Gamma(n,q\bar{\lambda})}{\Gamma(n,0)}e^{-\bar{\lambda_{1}}(\sigma^{2}+q)z},\end{split}$
(11)

where $\Theta$ is the scaled and shifted random variable version 444 Here
$\bar{\lambda}=\lambda_{2}/p$ is the scaled rate parameter of SUs and
$\bar{\lambda_{1}}=\lambda_{1}/p$ is the scaled rate parameter for PU. of $z$
given by $\Theta=\bar{\lambda}+\bar{\lambda_{1}}z$. Fig. 4 shows the plot of
the derived theoretical expression with simulation data for the two cases of
$p<q$ and $p>q$ with different SU densities ($n=1,2,3$).

Figure 4: PDF of SINR for two cases of $p<q$ and $p>q$ for different number of
SUs ($n=1,2,3$).

In the following sections, we will look into the crucial performance metrics
(outage probability and capacity) of underlay cognitive network. The mentioned
approach can be extended to the case of interference from multiple SUs given
that the important SINR expression Eq. (11) for multiple SUs is already been
derived. However, given the space limitations, the derivations considering
multiple SUs are not detailed herein in coming sections. Nonetheless, the
fundamental case of interference from single SU case has been presented in
detail.

###  III-B Mean SINR

The mean SINR is given as $\mu=\int_{0}^{\infty}zf_{z}(z)dz$ where the PDF of
SINR $f_{z}(z)$ was derived in Eq. (10). Thus,

$\begin{split}\mu&=\int_{0}^{\infty}z\Bigl{\\{}\frac{1}{p(z+1)}\Bigl{\\{}e^{\frac{-z}{p}}\left(1+\frac{p}{z+1}\right)\\\
&~{}~{}~{}+e^{-\left(\frac{z+q(z+1)}{p}\right)}\left((1+q)z-\frac{p}{z+1}\right)\Bigr{\\}}dz,\\\
&=\int_{0}^{\infty}\left(\frac{ze^{\frac{-z}{p}}}{p(z+1)}\right)dz+\int_{0}^{\infty}\left(\frac{ze^{\frac{-z}{p}}}{(z+1)^{2}}\right)\
dz\\\
&~{}~{}~{}~{}~{}~{}~{}~{}+\int_{0}^{\infty}\left(\frac{z^{2}(1+q)e^{-\left(\frac{z+q(z+1)}{p}\right)}}{p(z+1)}\right)dz\\\
&~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}-\int_{0}^{\infty}\left(\frac{ze^{-\left(\frac{z+q(z+1)}{p}\right)}}{(z+1)^{2}}\right)dz.\\\
\end{split}$

which on further simplification reduces to

$\begin{split}\mu&=e^{\frac{1}{p}}\left\\{\Gamma{\left(0,\frac{1}{p}\right)}-\Gamma{\left(0,\frac{1+q}{p}\right)}\right\\}+\frac{pe^{\frac{-q}{p}}}{1+q}.\end{split}$
(12)

Figure 5: Mean SINR vs Interference Temperature, $q$, for $p=2$ and $p=4$.

Fig. 5 shows that the change in mean SINR while varying the interference
temperature $q$ for a constant peak transmit power $p$. The higher transmit
power (for both PU and SU) with lower interference temperature gives better
mean SINR than lower transmit power (for both PU and SU) with high IT
constraint.

###  III-C Outage Probability of Primary Network

The outage probability is defined as the probability when the instantaneous
SINR drops below a given threshold. Mathematically, this is given as:
$Pr(\gamma\leq\psi)=F_{z}(\psi)$, which is nothing but the CDF of SINR.
Therefore,

$\begin{split}F_{z}(\psi)&=\int_{0}^{\psi}f_{z}(z)dz,\\\
&=\int_{0}^{\psi}\left(\frac{e^{\frac{-z}{p}}}{p(z+1)}\right)dz+\int_{0}^{\psi}\left(\frac{e^{\frac{-z}{p}}}{(z+1)^{2}}\right)dz\\\
&~{}~{}~{}~{}+\int_{0}^{\psi}\left(\frac{z(1+q)e^{-\left(\frac{z+q(z+1)}{p}\right)}}{p(z+1)}\right)dz\\\
&~{}~{}~{}~{}~{}~{}~{}~{}~{}-\int_{0}^{\psi}\left(\frac{e^{-\left(\frac{z+q(z+1)}{p}\right)}}{(z+1)^{2}}\right)dz.\end{split}$

which on further integration and simplification reduces to

$\begin{split}F_{z}(\psi)&=1-\frac{e^{\frac{-\psi}{p}}}{\psi+1}\left(1+\psi
e^{\frac{-q(\psi+1)}{p}}\right)\end{split}$

Figure 6: Outage probability of PU for $p<q$, where $p=2$ and $q=4$ and for
$p>q$, where $p=4$ and $q=2$.

It can be directly inferred from Fig. 6 that if $q>p$, the outage probability
is higher than in the case of $p>q$. In addition to this inference, it can be
also observed that the theoretical expressions derived are in sync with the
simulation results, i.e., increase the spectral efficiency of the network.

###  III-D Instantaneous Capacity of Primary Network

The PDF of instantaneous capacity can be readily found from the PDF of
instantaneous SINR by using transformation of random variables method [7, 8].
This can be obtained by using:

$\displaystyle f_{x}(x)=f_{z}(z)\Bigl{|}\frac{dz}{dx}\Bigr{|}_{z=e^{x}-1},$

where $f_{z}(z)$ is derived in Eq. (10) for the case of interference from
single SU on primary network.

Figure 7: Instantaneous PDF of Capacity for $p<q$, where $p=2$ and $q=4$ and
for $p>q$, where $p=4$ and $q=2$.

It can be seen from Fig. 7 that there is a point where the instantaneous
capacity for $p<q$ goes below $p>q$. It proves the point that the interference
temperature should not be kept constant rather should be dynamic in nature to
exploit full potential of the network.

###  III-E Mean Capacity

The average capacity from the PDF of instantaneous SINR ($f_{z}(z)$) is given
as:

$\displaystyle\bar{C}=\int_{0}^{\infty}\text{log}(1+z)f_{z}(z)dz.$

Substituting Eq. (10) in the above expression and on further evaluation.

$\begin{split}\bar{C}&=\frac{e^{\frac{1}{p}}}{p}\left[\Gamma\left(0,\frac{z+1}{p}\right)-\Gamma\left(0,\frac{(q+1)(z+1)}{p}\right)\right.\\\
&~{}~{}~{}~{}~{}~{}~{}~{}\times(p+q+1)\bigg{]}-\frac{e^{\frac{1-(q+1)(z+1)}{p}}}{z+1}\left\\{e^{\frac{q(z+1)}{p}}\right.\\\
&~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}\times\Big{(}1+\text{log}(z+1)\Big{)}+z\text{log}(z+1)-1\Big{\\}}\Bigg{|}_{0}^{\infty}.\end{split}$

At $z=\infty$, $\Gamma(0,z)\to 0$ and also,
$\frac{e^{\frac{1-(q+1)(z+1)}{p}}}{z+1}\to 0$.
Therefore, the final mean capacity expression will be evaluated at $z=0$:

$\begin{split}\bar{C}&=1-e^{-\frac{q}{p}}+\frac{e^{\frac{1}{p}}}{p}\Bigg{[}(p+q+1)\Gamma\left(0,\frac{q+1}{p}\right)-\Gamma\left(0,\frac{1}{p}\right)\Bigg{]}\end{split}$

Fig. 8 shows the plot of this theoretical expression with simulation results
for two cases of transmit power $p=2$ and $p=4$. Intuitively, high transmit
power $p=4$ will result in high capacity for the network than the low transmit
power of $p=2$ but when the interference temperature is relaxed, the
interference caused due to secondary user on primary will also increase that
in turn will reduce the overall capacity.

Figure 8: Theoretical and simulation result plots for capacity at $p=2$ and
$p=4$ with varying interference temperature: $q$.

##  IV Conclusion

In this paper, the performance of primary network is studied considering
interference from the SU network. The analysis is done under peak power
adaptation method at secondary transmitter. Given the importance of dynamic
interference temperature for network performance, closed form expressions for
the PDF and CDF of interference and noise, SINR for interference from single
and multiple SUs are derived. Furthermore, instantaneous capacity with
theoretical expressions for mean SINR, mean capacity and outage probability
are deduced for simplistic network consisting of interference from single SU.
Finally, the theoretical expressions are validated with the simulation
results.

## Acknowledgment

This work was supported by NASA Oklahoma Space Grant Consortium (EPSCoR)
Research Initiation Grant. Also, the authors would like to thank the TPC and
reviewers for their valuable feedback and suggestions.

## References

  * [1] S. Ekin, M. M. Abdallah, K. A. Qaraqe, and E. Serpedin, “Random subcarrier allocation in ofdm-based cognitive radio networks,” _IEEE Transactions on Signal Processing_ , vol. 60, no. 9, pp. 4758–4774, 2012.
  * [2] L. Musavian and S. Aïssa, “Capacity and power allocation for spectrum-sharing communications in fading channels,” _IEEE Transactions on Wireless Communications_ , vol. 8, no. 1, pp. 148–156, 2009.
  * [3] H. Tran, M. A. Hagos, M. Mohamed, and H.-J. Zepernick, “Impact of primary networks on the performance of secondary networks,” in _Computing, Management and Telecommunications (ComManTel), 2013 International Conference on_. IEEE, 2013, pp. 43–48.
  * [4] T. W. Ban, W. Choi, B. C. Jung, and D. K. Sung, “Multi-user diversity in a spectrum sharing system,” _IEEE Transactions on Wireless Communications_ , vol. 8, no. 1, pp. 102–106, 2009.
  * [5] S. Ekin, F. Yilmaz, H. Celebi, K. A. Qaraqe, M.-S. Alouini, and E. Serpedin, “Capacity limits of spectrum-sharing systems over hyper-fading channels,” _Wiley Wireless Communications and Mobile Computing_ , vol. 12, no. 16, pp. 1471–1480, 2012.
  * [6] R. Zhang, “On peak versus average interference power constraints for protecting primary users in cognitive radio networks,” _IEEE Transactions on Wireless Communications_ , vol. 8, no. 4, 2009.
  * [7] A. Papoulis and S. U. Pillai, _Probability, random variables, and stochastic processes_. Tata McGraw-Hill Education, 2002.
  * [8] S. Miller and D. Childers, _Probability and random processes: With applications to signal processing and communications_. Academic Press, 2012.
  * [9] H. Pishro-Nik, _Introduction to probability, statistics, and random processes_. Kappa Research LLC, 2014.
  * [10] J. H. Curtiss, “On the distribution of the quotient of two chance variables,” _The Annals of Mathematical Statistics_ , no. 4, pp. 409–421, 1941.


# Electron transfer at electrode interfaces via a straightforward
quasiclassical fermionic mapping approach

Kenneth A. Jung  Joseph Kelly  Thomas E. Markland<EMAIL_ADDRESS>Department of Chemistry, Stanford University, Stanford, California, 94305, USA

(August 28, 2024)

###### Abstract

Electron transfer at electrode interfaces to molecules in solution or at the
electrode surface plays a vital role in numerous technological processes.
However, treating these processes requires a unified and accurate treatment of
the fermionic states of the electrode and their coupling to the molecule being
oxidized or reduced in the electrochemical processes and, in turn, the way the
molecular energy levels are modulated by the bosonic nuclear modes of the
molecule and solvent. Here we present a physically transparent quasiclassical
scheme to treat these electrochemical electron transfer processes in the
presence of molecular vibrations by using an appropriately chosen mapping of
the fermionic variables. We demonstrate that this approach, which is exact in
the limit of non-interacting fermions, is able to accurately capture the
electron transfer dynamics from the electrode even when the process is coupled
to vibrational motions in regimes of weak coupling. This approach thus
provides a scalable strategy to explicitly treat electron transfer from
electrode interfaces in condensed-phase molecular systems.

##  I Introduction

Reactions at electrode interfaces are ubiquitous in industrial processes and
omnipresent in chemical, medical and energy research. Accurately capturing
these processes requires treating the quantum dynamics of the electrode states
and their coupling to the electronic and nuclear dynamics of the molecules in
the electrolyte. The vastly different timescales of these dynamical processes
as well as the combined size of the electronic and nuclear Hilbert spaces
makes accurate simulations of these systems with atomisitic detail highly
challenging.

Semiclassical and quasiclassical methods, which use classical trajectories to
approximate exact quantum dynamics, provide a potentially appealing route to
treat the dynamics of condensed phase systems with lead(electrode)-molecule
interactions due to their efficiency and low scaling with dimensionality. The
two approaches that have generated the most methods for treating these systems
semiclassically are those based on surface hopping and those that arise from
mapping Hamiltonians.

Recently introduced surface hopping based methods include the independent
electron surface hoppingShenvi, Roy, and Tully (2009a, b); Gardner  _et al._
(2023) (IESH) and the classical master equation evaluated with surface
hoppingDou, Nitzan, and Subotnik (2015a, b, 2016) (CME-SH) approaches. The
CME-SH approach has been combined with a diffusive description of the
surrounding condensed phase environment to describe a range of electrochemical
properties.Coffman and Subotnik (2018); Coffman  _et al._ (2019, 2020) The CME
also provides a link to methods which incorporate the effects of the electrons
in the lead using electronic frictionHead-Gordon and Tully (1995); Dou, Miao,
and Subotnik (2017); Dou (2018). However, the derivation of the CME assumes
weak lead-molecule coupling and requires introducing an ad hoc broadening
procedure to properly treat the lead-molecule hybridization for stronger
coupling regimes.

Semiclassical mapping methods, which map the discrete, single-particle
creation and annihilation operators to continuous, classical degrees of
freedom, provide an appealing approach to obtain an efficient classical-like
description of the dynamics of lead-molecule systems. The Li-Miller mappingLi
and Miller (2012) (LMM) has been shown to give accurate results when applied
to describing electronic transport through molecular junctionsLi  _et al._
(2013, 2014); Levy  _et al._ (2019) (the Anderson impurity model). The
apparent success of LMM for these systems occurs despite its mapping of
fermionic creation annihilation operators, which obey fermionic symmetries,
onto Cartesian degrees of freedom that do not obey the fermionic anti-
commutation identities. The more recent complete quasiclassical mappingLevy
_et al._ (2019) (CQM) has extended the LMM to capture more general classes of
fermionic observables. The success of the LMM and CQM approaches in accurately
reproducing the dynamics in benchmark systems of nanoscopic transport despite
not satisfying the fundamental anticommutivity of the individual fermionic
operators initially seems quite remarkable and that a deeper reason exists for
their success. Recent work has provided general insights into the reasons that
mapping the fermionic dynamics of the lead onto bosons in many cases can give
accurate or even exact resultsSun, Sasmal, and Vendrell (2021); Montoya-
Castillo and Markland (2023). In particular, it has been shown that, in the
case of non-interacting electrons, one can exchange the fermionic and bosonic
dynamics, while retaining the fermionic statistics, and yet still obtain the
exact quantum dynamicsSun, Sasmal, and Vendrell (2021). More generally the
conditions (i.e. combination of the form of the Hamiltonian and observables)
under which individual fermionic creation and annihilation operators can be
rigorously replaced with their bosonic counterparts has been derivedMontoya-
Castillo and Markland (2023).

Given the realization that under particular conditions one can obtain the
exact dynamics by mapping the fermionic lead to bosons thus allows for the use
of the Meyer-Miller (MM) mappingMeyer and Miller (1980), which has long been
appreciated as an exact mapping in the case of bosonsStock and Thoss (1997),
to be applied to develop a classical-like description of fermionic processes.
While Meyer-Miller mapping of the fermionic operators is no longer exact in
the case of interacting fermions, recent work has shown that it can still
provide an accurate description of purely fermionic systems.Sun, Sasmal, and
Vendrell (2021) Here we investigate and analyze the accuracy of MM mapping in
treating systems involving a mixture of fermionic (lead and molecule
electronic occupations) and bosonic (nuclear motion) degrees of freedom. In
particular, we investigate an Anderson-Holstein model consisting of a lead
(electrode) that can transfer electrons to a molecule in the electrolyte with
the transfer mediated by a vibrational mode. This Hamiltonian has formed the
basis of numerous previous studies to benchmark techniques for use in
electrochemical simulations Navrotskaya and Hammes-Schiffer (2009); Dou,
Nitzan, and Subotnik (2015a); Coffman  _et al._ (2019, 2020). From this we are
able to contrast the benefits of MM mapping to surface hopping based
approaches. We show that the mapping approach has advantages over surface
hopping approaches in that the states of the lead are straightforward to
discretize and it is accurate for strong lead-molecule coupling. However, the
mapping approach suffers from detailed balance issues that become increasingly
apparent in the case of strong molecule-vibration coupling.

##  II Theory

We focus on the dynamics arising from a Hamiltonian of the Anderson-Holstein
form,

$\hat{H}=\underbrace{\hat{H}_{\textrm{elec}}+\hat{H}^{(U)}_{\textrm{vib}}+\hat{H}_{\textrm{elec-
vib}}}_{\hat{H}_{M}}+\hat{H}_{L}+\hat{H}_{\textrm{elec-}L}.$ (1)

In particular, here we consider the case of a molecule $\hat{H}_{M}$ coupled
to a single lead (electrode) with Hamiltonian $\hat{H}_{L}$ by
$\hat{H}_{\textrm{elec-}L}$. The molecular Hamiltonian is comprised of the
electronic states, vibrational states, and the coupling between them described
by $\hat{H}_{\textrm{elec}}$, $\hat{H}^{(U)}_{\textrm{vib}}$, and
$\hat{H}_{\textrm{elec-vib}}$ respectively. The molecule receives charge
through its interaction with the lead which can alter the charge from its
initial “uncharged” state to a charged state. This system is shown
schematically in Fig. 1.

The particular form of the Anderson-Holstein Hamiltonian in Eq. 1 (also
sometimes referred to as the Anderson-Newns modelNEWNS (1969)) has been
previously been used to describe electron and proton transfer in
electrochemical reactionsNavrotskaya and Hammes-Schiffer (2009) as well as
adsorption/desorption from metal interfacesErpenbeck  _et al._ (2018). In the
most general form, the molecular Hamiltonian can be written as a sum of the
charged (C) and uncharged (U) states as

$\displaystyle\hat{H}_{M}$ $\displaystyle=$ $\displaystyle
H^{(C)}_{\textrm{vib}}(\hat{\bm{Q}},\hat{\bm{P}})\hat{d}^{\dagger}\hat{d}+H^{(U)}_{\textrm{vib}}(\hat{\bm{Q}},\hat{\bm{P}})[1-\hat{d}^{\dagger}\hat{d}],$
(2)

where

$H^{(U)}_{\textrm{vib}}(\hat{\bm{Q}},\hat{\bm{P}})=\frac{1}{2}\hat{\bm{P}}^{T}\hat{\bm{M}}^{-1}\hat{\bm{P}}+V^{(U)}(\hat{\bm{Q}}),$
(3)

and

$H^{(C)}_{\textrm{vib}}(\hat{\bm{Q}},\hat{\bm{P}})=\frac{1}{2}\hat{\bm{P}}^{T}\hat{\bm{M}}^{-1}\hat{\bm{P}}+V^{(C)}(\hat{\bm{Q}})+\epsilon_{M}.$
(4)

The set of operators $\bm{\hat{Q}}$ and $\bm{\hat{P}}$ describe the position
and momentum of the nuclear degrees of freedom and $\epsilon_{M}$ is the
energy difference between the charged and uncharged states. The diagonal mass
matrix is denoted by $\bm{M}$ and the fermionic raising and lowering operators
corresponding to the molecule are denoted by $\hat{d}$ and $\hat{d}^{\dagger}$
which satisfy the anticommutation rule $[\hat{d},\hat{d}^{\dagger}]_{+}=1$.
Eq. (2) can be rewritten as

$\hat{H}_{M}=\epsilon_{M}d^{\dagger}d+H^{(U)}_{\textrm{vib}}(\hat{\bm{Q}},\hat{\bm{P}})+\Delta
V(\hat{\bm{Q}})d^{\dagger}d,$ (5)

where we have defined $\Delta
V(\bm{\hat{Q}})=V^{(C)}(\hat{\bm{Q}})-V^{(U)}(\hat{\bm{Q}})$. By rewriting Eq.
(2) in this form and comparing with Eq. (1) one can see that

$\hat{H}_{\textrm{elec}}=\epsilon_{M}\hat{d}^{\dagger}\hat{d},$ (6)

and

$\hat{H}_{\textrm{elec-vib}}=\Delta V(\hat{\bm{Q}})\hat{d}^{\dagger}\hat{d}.$
(7)

The forms of $V^{(U)}(\hat{\bm{Q}})$ and $V^{(C)}(\hat{\bm{Q}})$ are kept
general for now but will be given specified forms in Sec. III. Although in the
present work we will use a simple form of the molecular Hamiltonian to allow
for comparison to exact results we emphasize that the mapping approach to
treating the fermionic dynamics described in Sec. II.1 is fully compatible
with atomistic treatments of molecular vibrations using diabatization
schemes.Cave and Newton (1996); Senthilkumar  _et al._ (2003); Kondov  _et
al._ (2007); Difley and Voorhis (2011); Oberhofer and Blumberger (2012); Mao,
Montoya-Castillo, and Markland (2019, 2020)

Figure 1: Schematic representation of the type of Anderson-Holstein models
studied in this work. The lead is described by $\hat{H}_{L}$, the initially
uncharged molecular system by
$\hat{H}_{M}=\hat{H}_{\textrm{elec}}+\hat{H}^{(U)}_{\textrm{vib}}+\hat{H}_{\textrm{elec-
vib}}$, and the coupling between the two by $\hat{H}_{\textrm{elec-}L}$.

Here the lead is assumed to be comprised of free electrons giving the specific
form of the Hamiltonian

$\hat{H}_{L}=\sum^{G}_{k}\epsilon_{k}\hat{c}^{\dagger}_{k}\hat{c}_{k},$ (8)

and the lead-molecule coupling is taken to be via an exchange mechanism of the
form

$\hat{H}_{C}=\sum^{G}_{k}t_{k}(\hat{d}^{\dagger}\hat{c}_{k}+\hat{c}^{\dagger}_{k}\hat{d}).$
(9)

The annihilation and creation operators $\hat{c}_{k}$ and
$\hat{c}^{\dagger}_{k}$ destroy and create electrons of energy $\epsilon_{k}$
in the lead and satisfy
$[\hat{c}_{k},\hat{c}_{k^{\prime}}^{\dagger}]_{+}=\delta_{kk^{\prime}}$ while
$\\{t_{k}\\}$ are the transfer coefficients between the lead states and the
electronic states of the molecule. The distribution of fermionic states in the
lead is determined by its spectral function $J_{L}(\epsilon)$ which allows the
flexibility to treat leads ranging from metals to semiconductors.

###  II.1 Classical mapping of a many-body bosonic and fermionic system

While it has been long established how to classically describe bosonic degrees
of freedom via the formal relationship between ladder operators and Cartesian
variablesMeyer and Miller (1980); Stock and Thoss (1997) the case for fermions
has faced considerably more difficulties. Previous work in this area has made
use of quaternions to map fermionic operators into classical variables.Li and
Miller (2012); Li  _et al._ (2013, 2014); Levy  _et al._ (2019) Recently it
was demonstrated that for non-interacting fermionic systems it is possible to
exactly replace the fermionic operators by classical-like phase space
variables using the MM mappingSun, Sasmal, and Vendrell (2021); Montoya-
Castillo and Markland (2023) giving the relations

$\hat{c}\to\frac{1}{\sqrt{2}}(q+ip),$ (10)
$\hat{c}^{\dagger}\to\frac{1}{\sqrt{2}}(q-ip).$ (11)

This mapping and the resulting classical equations of motion for $q$ and $p$
exactly reproduce the evolution of the one-body density matrix for
noninteracting fermionic systems while outside of this limit, it is an
approximation. It is important to emphasize that the classical oscillators of
the MM mapping are used here to count the occupation of the single-particle
fermionic states whose sum adds up to the total number of particles which
differs from its traditional use in counting the total electronic population
which always sums to one. Applying this mapping to the Hamiltonian of the
Anderson-Holstein form in Eq. (1) gives

$\displaystyle
H(\bm{q},\bm{p},\bm{Q},\bm{P})=\frac{1}{2}U(\bm{Q})(q^{2}_{d}+p^{2}_{d}-2\gamma_{M})+H^{(U)}_{\textrm{vib}}(\bm{Q},\bm{P})$
$\displaystyle+\sum^{G}_{k}\frac{\epsilon_{k}}{2}(q^{2}_{k}+p^{2}_{k}-2\gamma_{L})+\sum^{G}_{k}t_{k}(q_{d}q_{k}+p_{d}p_{k}),$
(12)

with

$U(\bm{Q})=\Delta V(\bm{Q})+\epsilon_{M}.$ (13)

Here $\bm{q},\bm{p}$ are the mapped Cartesian coordinates that describe the
fermionic variables of both the molecule ($q_{d},p_{d}$) and lead
($\\{q_{k},p_{k}\\}$). In Eq. (II.1) we have also replaced the nuclear
coordinates with classical variables $\bm{Q},\bm{P}$ which is a separate
approximation.

Each of these degrees of freedom evolves classically according to Hamilton’s
equations of motion

$\dot{q}_{d}=\frac{dH}{dp_{d}}=Up_{d}+\sum^{G}_{k}t_{k}p_{k},$ (14)
$\dot{p}_{d}=-\frac{dH}{dq_{d}}=-Uq_{d}-\sum^{G}_{k}t_{k}q_{k},$ (15)
$\dot{Q}_{j}=\frac{dH}{dP_{j}}=\frac{dH^{(U)}_{\textrm{vib}}}{dP_{j}},$ (16)
$\displaystyle\dot{P}_{j}=-\frac{dH}{dQ_{j}}=-\frac{1}{2}\frac{dU}{dQ_{j}}(q^{2}_{d}+p^{2}_{d}-2\gamma_{M})-\frac{dH^{(U)}_{\textrm{vib}}}{dQ_{j}}$
$\displaystyle-\sum^{G}_{k}\frac{dt_{k}}{dQ_{j}}(q_{d}q_{k}+p_{d}p_{k}),$ (17)
$\dot{q}_{k}=\frac{dH}{dp_{k}}=\epsilon_{k}p_{k}+t_{k}p_{d},$ (18)
$\dot{p}_{k}=-\frac{dH}{dq_{k}}=-\epsilon_{k}q_{k}-t_{k}q_{d},$ (19)

where the dot refers to the derivative with respect to time. The second line
of Eq. (II.1) accounts for the fact that the transfer coefficients may depend
on the nuclear coordinates, which can easily be incorporated within the
mapping framework but is not studied in this work.

The factors $\gamma_{M}$ and $\gamma_{L}$ parameterize the zero-point energy
of the molecular and lead fermionic states respectively. Only for the choice
of $\gamma_{M}=\gamma_{L}=\gamma$ does the MM mapping exactly describe the
Anderson-Holstein model in the absence of nuclear modes, thus we will limit
our attention to this case. The value of $\gamma$ has traditionally been
restricted to lie between $0$ and $1$,Stock and Thoss (1997); Golosov and
Reichman (2001); Cotton and Miller (2013); Sun, Sasmal, and Vendrell (2021)
although recent work has discussed the use of negative values.He  _et al._
(2021) For noninteracting systems the classical equations of motion are
independent of $\gamma$ and it only affects the initial conditions (see Sec.
III). However, for interacting systems such as the case studied here the
choice of this $\gamma$ affects the resulting dynamics as can be seen in Eq.
(II.1). For the particular value of $\gamma$, we investigated values of
$\gamma=(\sqrt{3}-1)/2$Cotton and Miller (2013); Runeson and Richardson
(2019), which gives the best result for 2-state systems such as the spin-boson
problem, and $\gamma=0$, which has been shown to be the best choice as the
number of states becomes large.Runeson and Richardson (2019) Overall we find
that the latter of these choices yields slightly improved results over a wider
range of Hamiltonian parameters and thus are shown in Sec. IV while the
results using the former are provided in SI Sec. I.

Upon integrating the equations of motion given in Eqs. (14)-(19) one can use
the trajectories to compute time-dependent quantities such as the electronic
population of the molecule

$n_{M}(t)=\langle\hat{d}^{\dagger}\hat{d}\rangle\to\frac{1}{2}\langle
q^{2}_{d}+p^{2}_{d}-2\gamma\rangle,$ (20)

and the current

$I(t)=-\frac{d}{dt}\sum_{k}\langle\hat{c}^{\dagger}_{k}\hat{c}_{k}\rangle\to-\sum_{k}t_{k}\langle
p_{d}q_{k}-q_{d}p_{k}\rangle.$ (21)

where the mapping from the fermionic operators to the classical variables are
obtained by using Eqs. (10) and (11). Here the angular brackets denote
averaging with respect to the initial density matrix. Appendix A shows that
the MM mapping guarantees the total charge and energy are conserved just as is
the case with the exact quantum dynamics.

As discussed in the introduction, the MM mapping is not the only choice of a
classical mapping that is available for fermionic operators. Indeed, the LMMLi
and Miller (2012) and its generalized form known as the complete
quasiclassical map (CQM)Levy  _et al._ (2019) are also exact in the same
limits as the MM mapping. It was also shown in Ref. 17 that each of these
mappings perform similarly outside of the non-interacting fermion case. The
reason for this is likely due to the fact that despite each of these mappings
trying to capture fermionic properties they can actually be expressed in a
unified mapping framework, along with the MM mapping, which makes no
assumptions of fermionic canonical commutation relations Liu (2016, 2017). In
SI Sec. III we compute the results obtained using the CQM approach for the
Anderson-Holstein model we investigate here and demonstrate that it gives
almost numerically indistinguishable results to those obtained using MM
mapping. Since MM mapping only requires half the number of degrees of freedom
to describe the same system as the LMM and CQM in Sec. IV, we show the results
obtained using MM mapping.

##  III Model system and generation of initial conditions

To examine the ability of the MM mapping to describe the many-body physics of
the Anderson-Holstein model we specify the nuclear vibrational Hamiltonians in
Eq. (1) to take the forms,

$H^{(U)}_{\textrm{vib}}=\frac{1}{2}\hbar\omega(P^{2}+Q^{2}),$ (22)

and

$\displaystyle H^{(C)}_{\textrm{vib}}$ $\displaystyle=$
$\displaystyle\frac{1}{2}\hbar\omega(P^{2}+Q^{2})+\sqrt{2}gQ+\epsilon_{M},$
(23) $\displaystyle=$
$\displaystyle\frac{1}{2}\hbar\omega\left[P^{2}+\left(Q+\frac{\sqrt{2}g}{\hbar\omega}\right)^{2}\right]+\tilde{\epsilon}_{M}.$

In our case, the nuclear Hamiltonian thus consists of one harmonic vibrational
degree of freedom with frequency $\omega$ coupled to the two electronic states
of the molecule where $g$ mediates the electron-vibration coupling,
$\epsilon_{M}$ is the bias of the charged state and
$\tilde{\epsilon}_{M}=\epsilon_{M}-g^{2}/\hbar\omega$ is the renormalized
molecular energy in the presence of the vibrational mode. This model was
previously investigated using the CME-SH approach that we compare to here.Dou,
Nitzan, and Subotnik (2015a) The initial density matrix is

$\rho=\rho_{M}\otimes\rho^{\textrm{eq}}_{\textrm{vib}}\otimes\rho^{\textrm{eq}}_{L},$
(24)

where $\rho_{M}$ represents the initial density of the fermionic degrees of
freedom of the molecule where the molecule is initialized such that only
$H^{(U)}_{\textrm{vib}}$ is populated i.e. $n_{M}(0)=0$ (see Eq. 20). The
nuclear vibrational degrees of freedom are thermalized in this state according
to the Wigner distribution given by,

$\rho^{\textrm{eq}}_{\textrm{vib}}=\frac{\alpha}{\pi}e^{-\alpha(Q^{2}+P^{2})},$
(25)

where $\alpha=\tanh(\beta\hbar\omega/2)$ and $\beta=1/(k_{B}T)$ is the inverse
temperature. In the limit of $\beta\hbar\to 0$ this reduces to the classical
Boltzmann distribution.

$\rho^{\textrm{eq}}_{L}$ is the thermal equilibrium density matrix of the lead
states. Since the electrons in the lead are in thermal equilibrium, the
population of a state with energy $\epsilon$ is given by the Fermi-Dirac
distribution

$n_{f}(\epsilon)=\frac{1}{1+e^{\beta(\epsilon-\mu)}},$ (26)

where $\mu$ is the chemical potential. Here, we assume a lead that is metallic
in nature such that the fermionic states of the lead form a continuum and are
described by the continuous spectral function

$J_{L}(\epsilon)=\frac{\Gamma}{(1+e^{A(\epsilon-B/2)})(1+e^{-A(\epsilon+B/2)})},$
(27)

where $A$ and $B$ define the cutoff and width of the lead distribution and
$\Gamma$ controls the coupling strength of the lead to the molecule. The
transfer coefficients, $t_{k}$ (see Eq. 9), are modeled in the wide band limit
following the procedure of Refs. 15 and 16 giving

$t_{k}=\sqrt{\frac{J_{L}(\epsilon_{k})\Delta\epsilon}{2\pi}},$ (28)

where $\Delta\epsilon$ is the spacing between lead states. The initial
conditions for the molecular fermionic degrees of freedom are

$q_{d}(0)=\sqrt{2\gamma+2n_{M}(0)}\cos(\theta),$ (29)

and

$p_{d}(0)=-\sqrt{2\gamma+2n_{M}(0)}\sin(\theta),$ (30)

where $\theta$ is a uniform random number between $0$ and $2\pi$ and
$n_{M}(0)$ is initial charge state of the molecule. The lead degrees of
freedom are initialized similarly as

$q_{k}(0)=\sqrt{2\gamma+2n_{k}}\cos(\theta),$ (31)

and

$p_{k}(0)=-\sqrt{2\gamma+2n_{k}}\sin(\theta),$ (32)

where the $n_{k}$ are sampled using the procedure of Ref. Levy  _et al._ ,
2019 where each $n_{k}$ is either zero or one such that the Fermi-Dirac
distribution is reproduced upon ensemble averaging. In what follows all units
are reported in atomic units unless otherwise stated. Additional Computational
details can be found in the SI Sec. IV.

##  IV Results

To investigate the performance of the MM mapping in describing the Anderson-
Holstein model we first look at cases without the presence of the vibrational
mode, in which the mapping should yield the exact result regardless of the
strength of the lead-molecule couplingSun, Sasmal, and Vendrell (2021);
Montoya-Castillo and Markland (2023). We then consider the effect of coupling
this system to a vibrational mode of the lead and how this alters the ability
to capture the properties of the system.

###  IV.1 Lead-molecule system

When $g=0$ in Eq. (23) the nuclear vibrational mode is decoupled from the
system and it reduces to an Anderson model consisting of just the lead and
molecule (often referred to as a dot). In the limit of small lead-molecule
coupling, $\Gamma$, the equilibrium molecular population in the presence of a
wide band lead is given byBruus and Flensberg (2004); Stefanucci and van
Leeuwen (2013)

$n^{eq}_{M}=\int\frac{d\epsilon}{2\pi}\frac{\Gamma}{(\epsilon-\tilde{\epsilon}_{M})^{2}+(\Gamma/2)^{2}}n_{f}(\epsilon).$
(33)

In the case of $g\to 0$ (no electron-vibration coupling) Eq. (33) is exact for
all values of $\Gamma$. Using Eq. (33) we can test how well the classical
dynamics of the MM mapping recovers equilibrium properties both in the absence
and presence of nuclear motion as a function of the lead-molecule coupling
strength. Since, in the absence of nuclear modes the MM mapping gives a
quadratic Hamiltonian, which ensures the dynamics are exact, it should recover
the correct equilibrium population regardless of the lead-molecule coupling
strength.

Figure 2 shows that in this case ($g=0$) the MM mapping is able to capture the
correct equilibrium molecular population as the energy bias $\epsilon_{M}$ is
scanned revealing that the hybridization between the molecular level and the
lead is being properly described. This is not the case for the CME-SH
methodology,Dou, Nitzan, and Subotnik (2015a, b) which only recovers correct
equilibrium populations after being “broadened” to account for the
hybridization between the lead and molecule. While the broadening schemes that
have been used in conjunction with CME-SH are effective in correcting the
equilibrium behavior in the absence of nuclear motion it is unclear how to
extend them beyond cases where analytic results provide guidance. In contrast,
the MM mapping naturally incorporates the hybridization through the explicit
dynamics of the lead degrees of freedom and does not need to be modified to
capture the equilibrium behavior in the absence of nuclear motions. In
addition, the current (Eq. (21)) is correctly zero when the system is in
equilibrium.

Figure 2: Left: Transient dynamics of the molecule population with
$\Gamma=0.006$, $\mu=0.0$ and $\beta=100$. Dotted lines are the transient
populations computed using the MM mapping, dashed lines are the equilibrium
value predicted by the CME-SH method without broadening and solid lines are
the exact values obtained from Eq. (33). Right: Lead distribution after 750
fs. Colored shapes are distributions computed from the MM mapping and the
solid black line is the Fermi-Dirac distribution. 10000 trajectories were used
to obtain these results.

###  IV.2 Lead-molecule-vibration system

We now consider the case where the lead-molecule system is now coupled to a
nuclear vibrational mode (see Fig. 1) which can dissipate energy. The strength
of the coupling of the molecule’s electronic population to the vibrational is
controlled by $g$ (Eq. 23). Due to the presence of the nuclear motion coupled
to the molecular electronic degrees of freedom the MM representation is no
longer exact.

Figure 3 shows the equilibrium molecule population (obtained at 1.5 ps) as a
function of the energy bias for a series of electron-vibration coupling ($g$)
values. The parameter regimes shown correspond to those in Ref. 4 such that
comparisons can be made to the CME-SH approach. Despite the MM mapping not
being exact we see that it is able to accurately predict the correct
population for the cases when $g=0.001$, $g=0.0025$, and $g=0.005$. However,
for the strongest coupling to the vibrational mode $g=0.0075$ we find that the
MM mapping underestimates the long-time population. We also observe that MM
mapped dynamics more accurately capture the equilibrium molecule population at
large positive bias ($\epsilon_{M}>0$) than negative bias. This is likely due
to the fact that at high positive bias the molecule’s energy is shifted high
into the Fermi distribution where the population of lead states is much
smaller which effectively reduces the amount of coupling between the lead and
the molecule.

Figure 3: Equilibrium molecule population as a function of $\epsilon_{M}$.
$\Gamma=0.003$, $\mu=0.0$, $\hbar\omega=0.003$ and $\beta=100$. The solid line
are the exact NRG results and dots are from the explicit MM mapping dynamics
using 10000 trajectories. Here $\gamma=0$. $1200$ lead states were used
ranging in energy from $-0.1$ to $0.1$ except for the case where $g=0.0075$ in
which $2400$ states were used ranging from $-0.2$ to $0.2$.

To further test MM mapping’s ability to capture the correct long-time
molecular population in Fig. 4 we plot the molecular population as a function
of the renormalized bias, $\tilde{\epsilon}_{M}$, for a case with strong
vibrational coupling and weak lead coupling (top panel) and a case with strong
lead coupling and weak vibrational coupling (bottom panel). Our MM mapping
results are compared with the numerically exact numerical renormalization
group (NRG) and CME-SH calculations without broadening in Fig. 4. As expected
from the observations in Fig. 3 the MM mapping is not able to reproduce the
correct equilibrium molecule population in the case of strong vibrational
coupling but is very accurate for the case of weak vibrational coupling even
when the lead is strongly coupled to the molecule. The latter result follows
from the fact that the MM mapping describes the bare fermionic system exactly
regardless of the value of $\Gamma$ (see Sec. IV.1). A weakly coupled nuclear
mode does not significantly affect this description and demonstrates that
there are regimes in which a mapping description that employs fully linearized
trajectories may be more accurate than the CME-SH since there is no ambiguity
in how to broaden the results. However, in cases with strong coupling to the
vibrational mode, even without broadening, the CME-SH is able to more
accurately capture the exact result. The failure of mean-field dynamics and
related methods such as classical mappings to capture detailed balance for
electronic systems coupled to bosonic modes has been well documented in the
literature and again presents a challenge in the cases of stronger coupling to
the vibrational mode.Bellonzi, Jain, and Subotnik (2016); Dou, Nitzan, and
Subotnik (2016); Jain and Subotnik (2017)

Figure 4: Equilibrium molecule population as a function of
$\tilde{\epsilon}_{M}$. $\mu=0.0$, $\hbar\omega=0.003$ and $\beta=100$. Top
panel: $\Gamma=0.003$ and $g=0.075$ Bottom panel: $\Gamma=0.01$ and $g=0.025$
Solid lines are the NRG results, dots are from the MM mapping using 10000
trajectories and dashed lines are unbroadened CME-SH results taken from Ref.
4. Broadened CME-SH results are indistinguishable from the NRG results and are
not shown.

To investigate the failure of the MM mapping in capturing detailed balance, we
performed a calculation in which the molecule with no bias between the charged
and uncharged states ($\tilde{\epsilon}=0$) and the lead are initialized in
thermal equilibrium. Due to the absence of a bias between the charged states
and since everything is at the same temperature there should be no net
exchange of energy between the lead and molecular degrees of freedom and hence
the kinetic energy should not evolve with time. The left panel of Fig. 5 shows
the time evolution of the kinetic energy of the vibrational mode where we see
that the MM mapping shows a drift in the average kinetic energy until it
reaches a steady state value of $\sim 75\%$ of its initial value. The right
panel of Fig. 5 shows the initial (black) and final (blue) populations of the
lead states where, after magnifying the difference between them (green), we
observe that those near the value of $\tilde{\epsilon}=0$ have lost population
which has been redistributed to the surrounding states. From the green line we
see that more of the population has been pushed into states of higher energy.
This increase of energy in the lead is responsible for the decrease in the
average kinetic energy of the molecule. Similar results are found when the
number of states in the lead is increased verifying that this effect is not an
artifact of improper lead discretization (see SI. Sec. II).

Figure 5: The molecule with zero bias and the lead are both initialized in
thermal equilibrium. Left: the average kinetic energy from the MM mapping
dynamics (red line) and exact value (black dashed line). Right: The initial
(black line) and final (blue line) lead distribution. The green is the
difference between the two and has been magnified by a factor of 4 for
clarity. $\Gamma=0.01$, $\mu=0.0$, $\hbar\omega=0.003$, $g=0.025$,
$\beta=100$.

##  V Conclusion

We have presented a simple and practical method for simulating electrochemical
systems using classical trajectories. By utilizing the MM mapping to describe
fermionic degrees of freedom we have shown it is possible to capture the
fundamental physics of lead-molecule and electron-vibration coupling in the
Anderson-Holstein model. Most notable is the mapping approach’s ability to
correctly capture the lead-molecule hybridization. The mapping approach thus
provides advantages over other methods that could be applied to an atomistic
description of the molecular degrees of freedom such as the CME-SH
approachDou, Nitzan, and Subotnik (2015a, b) since it naturally captures the
correct electronic population of the bare fermionic system without the need
for an ad hoc broadening procedure. However, while surface hopping methods
accurately capture detailed balance,Dou, Nitzan, and Subotnik (2015a); Jain
and Subotnik (2017); Schmidt, Parandekar, and Tully (2008) we have shown that,
consistent with other MM mapped systems, when strong coupling to the nuclear
motion is present detailed balance issues arise leading to incorrect molecular
populations. Given that a number of strategies have been introduced to improve
the detailed balance failures of mapping approaches to semiclassical dynamics
these provide a potential avenue for improvements in future work.

Ultimately we believe this straightforward procedure should be usable in a
variety of contexts where strong molecule-lead and weak electron-vibration
couplings are present making it complementary to other trajectory-based
approaches as well as perturbative master equations. The flexibility that this
method affords should make it possible to extend it to describe a range of
electrochemically relevant properties such as adsorption, charge transfer in
atomistic environments as well as bond breaking/formation at charged surfaces
when combined with _ab initio_ descriptions of the forces.

## Acknowledgments

This work was supported by National Science Foundation Grant No. CHE-2154291.

## Data Availability

The data that supports the findings of this study are available within the
article and its supplementary material.

##  Appendix A Conservation of charge and energy

The total charge is given by the sum of the fermionic population in the lead
and the molecule

$\hat{T}=\hat{d}^{\dagger}\hat{d}+\sum^{G}_{k}\hat{c}^{\dagger}_{k}\hat{c}_{k}.$
(34)

Using the Heisenberg equation of motion and Eq. (1) it is straightforward to
show that

$\frac{d}{dt}\hat{T}=i[\hat{H},\hat{T}]=0,$ (35)

meaning that charge is conserved in this system when treated quantum
mechanically. Likewise, it can also be shown that when the system is mapped
classically using the MM mapping that charge is also conserved. Within the
mapping the total charge is given by

$T=\frac{1}{2}(q^{2}_{d}+p^{2}_{d}-2\gamma)+\frac{1}{2}\sum^{G}_{k}(q^{2}_{k}+p^{2}_{k}-2\gamma).$
(36)

From this definition and using the classical equations of motion given in Eqs.
(14)-(19) it follows that

$\frac{d}{dt}T=q_{d}\dot{q}_{d}+p_{d}\dot{p}_{d}+\sum^{G}_{k}(q_{k}\dot{q}_{k}+p_{k}\dot{p}_{k})=0.$
(37)

Similarly the total energy is conserved, which trivially follows from the fact
that quantum Hamiltonian commutes with itself and Hamilton’s equations of
motion (Eqs. (14)-(19)) conserve the total energy of the system. These results
reveal that regardless of the value of $\gamma$ the charge and energy of the
total system is conserved on the trajectory level. While it is true that
charge and energy are conserved at the trajectory level this does not imply
that individual trajectories are physically meaningful. One must still average
a sufficient number of configurations to properly approximate the quantum
dynamics of the system.

## References

  * Shenvi, Roy, and Tully (2009a) N. Shenvi, S. Roy, and J. C. Tully, The Journal of Chemical Physics 130, 174107 (2009a).
  * Shenvi, Roy, and Tully (2009b) N. Shenvi, S. Roy, and J. C. Tully, Science 326, 829 (2009b).
  * Gardner  _et al._ (2023) J. Gardner, D. Corken, S. M. Janke, S. Habershon, and R. J. Maurer, The Journal of Chemical Physics 158, 064101 (2023).
  * Dou, Nitzan, and Subotnik (2015a) W. Dou, A. Nitzan, and J. E. Subotnik, The Journal of Chemical Physics 142, 084110 (2015a).
  * Dou, Nitzan, and Subotnik (2015b) W. Dou, A. Nitzan, and J. E. Subotnik, The Journal of Chemical Physics 142, 234106 (2015b).
  * Dou, Nitzan, and Subotnik (2016) W. Dou, A. Nitzan, and J. E. Subotnik, The Journal of Chemical Physics 144, 074109 (2016).
  * Coffman and Subotnik (2018) A. J. Coffman and J. E. Subotnik, Physical Chemistry Chemical Physics 20, 9847 (2018).
  * Coffman  _et al._ (2019) A. J. Coffman, A. K. Harshan, S. Hammes-Schiffer, and J. E. Subotnik, The Journal of Physical Chemistry C 123, 13304 (2019).
  * Coffman  _et al._ (2020) A. J. Coffman, W. Dou, S. Hammes-Schiffer, and J. E. Subotnik, The Journal of Chemical Physics 152, 234108 (2020).
  * Head-Gordon and Tully (1995) M. Head-Gordon and J. C. Tully, The Journal of Chemical Physics 103, 10137 (1995).
  * Dou, Miao, and Subotnik (2017) W. Dou, G. Miao, and J. E. Subotnik, Physical Review Letters 119 (2017), 10.1103/physrevlett.119.046001.
  * Dou (2018) W. Dou,  _Modeling Nonadiabatic Dynamics At Molecule-Metal Interfaces_ , Ph.D. thesis, University of Pennsylvania (2018).
  * Li and Miller (2012) B. Li and W. H. Miller, The Journal of Chemical Physics 137, 154107 (2012).
  * Li  _et al._ (2013) B. Li, T. J. Levy, D. W. H. Swenson, E. Rabani, and W. H. Miller, The Journal of Chemical Physics 138, 104110 (2013).
  * Li  _et al._ (2014) B. Li, E. Y. Wilner, M. Thoss, E. Rabani, and W. H. Miller, The Journal of Chemical Physics 140, 104110 (2014).
  * Levy  _et al._ (2019) A. Levy, W. Dou, E. Rabani, and D. T. Limmer, The Journal of Chemical Physics 150, 234112 (2019).
  * Sun, Sasmal, and Vendrell (2021) J. Sun, S. Sasmal, and O. Vendrell, The Journal of Chemical Physics 155, 134110 (2021).
  * Montoya-Castillo and Markland (2023) A. Montoya-Castillo and T. E. Markland, The Journal of Chemical Physics 158, 094112 (2023).
  * Meyer and Miller (1980) H.-D. Meyer and W. H. Miller, The Journal of Chemical Physics 72, 2272 (1980).
  * Stock and Thoss (1997) G. Stock and M. Thoss, Physical Review Letters 78, 578 (1997).
  * Navrotskaya and Hammes-Schiffer (2009) I. Navrotskaya and S. Hammes-Schiffer, The Journal of Chemical Physics 131, 024112 (2009).
  * NEWNS (1969) D. M. NEWNS, Physical Review 178, 1123 (1969).
  * Erpenbeck  _et al._ (2018) A. Erpenbeck, C. Schinabeck, U. Peskin, and M. Thoss, Physical Review B 97 (2018), 10.1103/physrevb.97.235452.
  * Cave and Newton (1996) R. J. Cave and M. D. Newton, Chemical Physics Letters 249, 15 (1996).
  * Senthilkumar  _et al._ (2003) K. Senthilkumar, F. C. Grozema, F. M. Bickelhaupt, and L. D. A. Siebbeles, The Journal of Chemical Physics 119, 9809 (2003).
  * Kondov  _et al._ (2007) I. Kondov, M. Čížek, C. Benesch, H. Wang, and M. Thoss, The Journal of Physical Chemistry C 111, 11970 (2007).
  * Difley and Voorhis (2011) S. Difley and T. V. Voorhis, Journal of Chemical Theory and Computation 7, 594 (2011).
  * Oberhofer and Blumberger (2012) H. Oberhofer and J. Blumberger, Physical Chemistry Chemical Physics 14, 13846 (2012).
  * Mao, Montoya-Castillo, and Markland (2019) Y. Mao, A. Montoya-Castillo, and T. E. Markland, The Journal of Chemical Physics 151, 164114 (2019).
  * Mao, Montoya-Castillo, and Markland (2020) Y. Mao, A. Montoya-Castillo, and T. E. Markland, The Journal of Chemical Physics 153, 244111 (2020).
  * Golosov and Reichman (2001) A. A. Golosov and D. R. Reichman, The Journal of Chemical Physics 114, 1065 (2001).
  * Cotton and Miller (2013) S. J. Cotton and W. H. Miller, The Journal of Chemical Physics 139, 234112 (2013).
  * He  _et al._ (2021) X. He, Z. Gong, B. Wu, and J. Liu, The Journal of Physical Chemistry Letters 12, 2496 (2021).
  * Runeson and Richardson (2019) J. E. Runeson and J. O. Richardson, The Journal of Chemical Physics 151, 044119 (2019).
  * Liu (2016) J. Liu, The Journal of Chemical Physics 145, 204105 (2016).
  * Liu (2017) J. Liu, The Journal of Chemical Physics 146, 024110 (2017).
  * Bruus and Flensberg (2004) H. Bruus and K. Flensberg,  _Many-Body Quantum Theory in Condensed Matter Physics: An Introduction (Oxford Graduate Texts)_, hardcover ed. (Oxford University Press, 2004) p. 435.
  * Stefanucci and van Leeuwen (2013) G. Stefanucci and R. van Leeuwen,  _Nonequilibrium many-body theory of quantum systems_ (Cambridge University Press, Cambridge, England, 2013).
  * Bellonzi, Jain, and Subotnik (2016) N. Bellonzi, A. Jain, and J. E. Subotnik, The Journal of Chemical Physics 144, 154110 (2016).
  * Jain and Subotnik (2017) A. Jain and J. E. Subotnik, The Journal of Physical Chemistry A 122, 16 (2017).
  * Schmidt, Parandekar, and Tully (2008) J. R. Schmidt, P. V. Parandekar, and J. C. Tully, The Journal of Chemical Physics 129, 044104 (2008).