These authors contributed equally to this work.

These authors contributed equally to this work.

[2]Rae Woong Park

[1]Edward Choi

1]KAIST 2]Ajou University School of Medicine

# KorMedMCQA: Multi-Choice Question Answering Benchmark for Korean Healthcare
Professional Licensing Examinations

Sunjun Kweon<EMAIL_ADDRESS>Byungjin Choi<EMAIL_ADDRESS>Minkyu
Kim<EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS>[  [

###### Abstract

We introduce KorMedMCQA, the first Korean multiple-choice question answering
(MCQA) benchmark derived from Korean healthcare professional licensing
examinations, covering from the year 2012 to year 2023. This dataset consists
of a selection of questions from the license examinations for doctors, nurses,
and pharmacists, featuring a diverse array of subjects. We conduct baseline
experiments on various large language models, including proprietary/open-
source, multilingual/Korean-additional pretrained, and clinical context
pretrained models, highlighting the potential for further enhancements. We
make our data publicly available on
HuggingFace(https://huggingface.co/datasets/sean0042/KorMedMCQA) and provide a
evaluation script via LM-Harness, inviting further exploration and advancement
in Korean healthcare environments.

##  1 Introduction

Evaluating language models’ expertise in the medical field can be effectively
done by testing them with questions found in medical licensing exams [1, 2, 3,
4, 5, 6]. This approach not only reveals their knowledge of medical subjects
but also their potential utility in healthcare environments. Datasets like
MedQA [7] from the U.S. and China, HeadQA [8] from Spain, FrenchMCQA [9] from
France, and MedMCQA [10] from India, all based on each country’s medical
licensing exams, have been used for this purpose.

Figure 1: Three sample questions of KorMedMCQA dataset respectively for
doctor, nurse, and pharmacist license examination. Note that the translation
has been added only to this figure to aid the understanding of global readers,
and the actual data is not translated. The samples are from the 2021 year
exams: 16th question of the doctor’s third session, 1st question of the
nurse’s first session, and 7th question of the pharmacist’s first session.

Yet, there’s a clear need for more diverse datasets, covering a broader range
of languages and country-specific exams. This is crucial for three reasons.
First, a model’s effectiveness is closely tied to its understanding of the
language used. No matter how well it knows clinical procedures or solves
problems, a model’s performance in other languages remains uncertain without a
proper testing. Second, medical practices and treatment outcomes can vary
significantly across different racial and national backgrounds. For example,
treatment guidelines for acute pyelonephritis(APN) differ between the U.S. and
Korea [11, 12], affecting how medical questions might be answered. Lastly,
medical exams test not just clinical knowledge but also understanding of
relevant laws and regulations, which vary by the country. Thus, excelling in
one country’s exam does not guarantee success in another due to these
linguistic, racial, and regulatory differences.

In conclusion, ensuring the effectiveness and relevance of language models in
the medical field of a specific country necessitates the development of
datasets that cater to the distinctive requirements of that country’s
examinations. To address this need, our study introduces the KorMedMCQA, the
first Korean medical multi-choice question answering (MCQA) dataset derived
from licensing examinations for doctors, nurses, and pharmacists in Korea
(refer to Figure 1 for an example). This dataset, comprising 5,345 questions
from exams administered between 2012 and 2023, offers a comprehensive tool for
evaluating the capability of models to understand medical knowledge within the
Korean context. We utilized this dataset to evaluate various large language
models, including proprietary and open-source models. Notably, the highest
average scores for each case were achieved by GPT4-0125 with an average score
of 83.06, and LLaMA2-70b with 53.09. For further research, the dataset is
available on Huggingface [13], with the evaluation script accessible via LM-
Harness [14].

##  2 Data

###  2.1 Data Collection

Our study gathered an extensive dataset of Korean medical licensing
examinations spanning from 2012 to 2023. This included exams for doctors from
2012 to 2023, and for nurses and pharmacists from 2019 to 2023. Every question
and its corresponding answer were obtained directly from the official releases
by the Korea Health Personnel Licensing Examination
Institute111https://www.kuksiwon.or.kr/EngHome/index.do, which made the exam
publicly available strictly for academic purposes.

The process of data collection involved systematically crawling publicly
available PDFs of the examination materials. These questions and answers were
then structured into datasets, with a Korean medical doctor overseeing the
conversion process to ensure accuracy and minimize errors. To maintain
consistency, questions with multiple correct answers were excluded.
Additionally, questions based on images or videos were also left out to avoid
copyright issues. Each dataset entry includes a unique identifier (denoting
the year, session, and question number), the question text, five potential
answer choices, and the index of the single correct answer.

###  2.2 Data Statistics

The whole dataset comprises 5,345 multi-choice question entries, including
2,339 for physicians, 1,460 for nurses, and 1,546 for pharmacists. The
distribution of questions across different years is detailed in table 1, with
specific partitions for model development and evaluation: 2021 was set aside
as the development set, while 2022 and 2023 served as the test set. The rest
formed the training set.

The dataset covers a broad spectrum of medical disciplines, with physicians’
exams featuring three main subjects (e.g., General Medicine, Specialized
Medicine), nurses’ exams encompassing seven subjects (e.g., Adult Nursing,
Maternal Nursing, Pediatric), and pharmacists’ exams including 21 subjects
(e.g., Biochemistry, Molecular Biology, Microbiology), showcasing a diverse
set of medical exams in Korean. The specific subjects are detailed further in
table 2.

Year | Doctor | Nurse | Pharm
---|---|---|---
2012 | 309 | - | -
2013 | 257 | - | -
2014 | 242 | - | -
2015 | 226 | - | -
2016 | 203 | - | -
2017 | 179 | - | -
2018 | 164 | - | -
2019 | 163 | 291 | 318
2020 | 147 | 291 | 314
2021 | 164 | 291 | 300
2022 | 135 | 295 | 301
2023 | 150 | 292 | 313
Total | 2339 | 1460 | 1546
(a) The number of entries by year

Split | Doctor | Nurse | Pharm
---|---|---|---
Train | 1890 | 582 | 632
Dev | 164 | 291 | 300
Test | 285 | 587 | 614
(b) The number of entries by split

Table 1: Statistics of KorMedMCQA dataset Table 2: Subjects across doctor, nurse, and pharmacist license examination in Korea | Subjects
---|---
Doctor | Healthcare Laws and Regulations, General Medicine, Specialized Medicine
Nurse |  | Healthcare Laws and Regulations, Adult Nursing Maternal, Nursing Pediatric,
---
Nursing Community Health Nursing, Psychiatric Nursing, Nursing Management,
Fundamental Nursing
Pharm |  | Biochemistry, Molecular Biology, Microbiology, Immunology, Pharmacology,
---
Preventive Medicine, Pathophysiology, Physical Pharmacy, Synthetic Chemistry,
Medicinal Chemistry, Pharmaceutical Analysis, Pharmaceutics, Pharmacognosy,
Herbal Medicine, Pharmacotherapy, Pharmacy Practice, Drug Manufacturing,
Drug Quality Control, Pharmacy Administration and Management (Social
Pharmacy),
Pharmacy Law

##  3 Evaluation

We conducted a comprehensive evaluation of various Large Language Models
(LLMs) utilizing the KorMedMCQA dataset. Our assessment methodology, inspired
by the work of [15], involved a consistent application of greedy decoding
across all models in a 5-shot scenario, where each “shot” was derived from the
development set (see Figure 2 for an example). The models’ performance was
measured based on their ability to generate answers that matched the correct
answer index, a technique applied to both proprietary and open-source LLMs.

Figure 2: Prompt for evaluating large language models in a 5-shot setting. The assessment is conducted by comparing the generated output with a specific answer index. Table 3:  Average accuracy across the KorMedMCQA dataset, with the highest-scoring large language models highlighted in bold for both proprietary and open-source categories. Model | Doctor | Nurse | Pharm | Avg
---|---|---|---|---
Proprietary Models
gpt-4-0125 | 76.49 | 89.95 | 82.74 | 83.06
gpt-4-1106 | 70.53 | 89.79 | 75.57 | 78.63
gpt-4-0613 | 75.09 | 85.86 | 83.22 | 81.39
gpt-3.5-turbo-0125 | 44.56 | 62.18 | 60.26 | 55.67
gpt-3.5-turbo-1106 | 41.4 | 60.31 | 58.14 | 53.28
gpt-3.5-turbo-0613 | 41.75 | 62.18 | 56.35 | 53.43
Multilingual Pretrained Models
llama2-70b | 42.46 | 63.54 | 53.26 | 53.09
llama2-13b | 29.82 | 39.52 | 33.22 | 34.19
llama2-7b | 25.96 | 24.02 | 28.01 | 26.00
Mistral-7b-v0.1 | 29.82 | 42.08 | 43.49 | 38.46
Yi-34b | 40 | 55.54 | 52.77 | 49.44
Yi-6b | 27.37 | 37.31 | 40.39 | 35.02
SOLAR-10.7B-v1.0 | 37.19 | 55.54 | 54.07 | 48.93
Multilingual Instruction-Tuned Models
llama2-70b-chat | 29.47 | 47.09 | 37.46 | 38.01
llama2-13b-chat | 29.47 | 36.63 | 30.94 | 32.35
llama2-7b-chat | 19.3 | 22.32 | 24.76 | 22.13
Mistral-7b-Instruct-v0.1 | 28.42 | 33.73 | 32.57 | 31.57
Yi-34b-chat | 1.05 | 5.96 | 17.26 | 8.09
Yi-6b-chat | 28.42 | 35.76 | 39.9 | 34.69
SOLAR-10.7B-Instruct-v1.0 | 40.7 | 60.14 | 53.42 | 51.42
Medical Continual Pretrained Models
Meditron-70b | 43.51 | 57.58 | 48.86 | 49.98
Meditron-7b | 24.91 | 24.7 | 26.22 | 25.28
Korean Continual Pretrained Models
Yi-Ko-6B | 38.6 | 54.68 | 47.07 | 46.78
llama-2-koen-13b | 30.18 | 46.68 | 36.97 | 37.94

In the case of proprietary models, our analysis included all available
versions of GPT-4 and GPT-3.5, specifically versions 0613, 1106, and 0125
[16]. For open-source models, we evaluated a selection of multilingual
pretrained LLMs: LLaMA2 [17], Mistral [18], Yi [19], SOLAR [20], and their
instruction-tuned counterparts. Additionally, we examined Meditron [3], which
represents an adaptation of LLaMA2 to an English medical context. Our
investigation also included two models [21, 22] that had been continually
pretrained with Korean text on multilingual LLMs. The results of these
evaluations are presented in Table 3.

The analysis revealed a consistent trend: the models tended to score lower in
doctor exams, with slightly higher scores in pharmacist exams, and the highest
scores in nurse exams. Among the proprietary models, the most recent GPT-4
version (0125) achieved the highest average score of 83.06. In contrast, among
open-source models, LLaMA2-70b was the top performer with a score of 53.09.
Interestingly, models that received further instruction fine-tuning on
multilingual pretrained models generally scored lower than those without fine-
tuning (e.g., LLaMA2-70b with 53.09 versus LLaMA2-70b-chat with 38.01). This
suggests that additional instruction fine-tuning did not necessarily translate
into improved performance.

Our comparison between Meditron and LLaMA2 revealed that pretraining with an
English medical context did not significantly enhance performance, a finding
that contrasts with the expected improvements seen with English-based medical
exams like MedQA or MedMCQA [3, 5]. This indicates that the Korean language
medical dataset may not benefit substantially from English medical context
pretraining. However, models continually pretrained on Korean text, even
though it is not medical related, showed improved performance (e.g., Yi-6b
with 35.02 versus Yi-Ko-6b with 46.78, and LLaMA2-13b with 34.19 versus
LLaMA2-KOEN-13b with 37.94), highlighting the value of aligning model training
with the specific language context of the dataset.

##  4 Conclusion

This study presents KorMedMCQA dataset, a collection of medical licensing
examination questions from Korea. The dataset, encompassing a wide range of
medical subjects and reflecting the intricacies of Korea’s healthcare
standards, regulations, and practices, can serve as a critical tool for
benchmarking the capabilities of language models in Korean healthcare
environments. Our findings reveal that while the latest version of GPT-4
demonstrated exceptional performance with an average score of 83.06, there
remains significant scope for improvement among open-source models, with
LLaMA2-70b leading at a score of 53.09. The findings also underscore the
necessity for models to deeply understand the Korean language and indicates
that pretraining with an English medical context does not yield substantial
benefits for Korean medical datasets. By making this dataset publicly
available, we aim to facilitate extensive research into the development of
models capable of navigating the complexities of Korean healthcare
environments.

## References

  *   * Nori et al. [2023a] Nori, H., King, N., McKinney, S.M., Carignan, D., Horvitz, E.: Capabilities of gpt-4 on medical challenge problems. arXiv preprint arXiv:2303.13375 (2023)
  * Nori et al. [2023b] Nori, H., Lee, Y.T., Zhang, S., Carignan, D., Edgar, R., Fusi, N., King, N., Larson, J., Li, Y., Liu, W., et al.: Can generalist foundation models outcompete special-purpose tuning? case study in medicine. arXiv preprint arXiv:2311.16452 (2023)
  * Chen et al. [2023] Chen, Z., Cano, A.H., Romanou, A., Bonnet, A., Matoba, K., Salvi, F., Pagliardini, M., Fan, S., Köpf, A., Mohtashami, A., et al.: Meditron-70b: Scaling medical pretraining for large language models. arXiv preprint arXiv:2311.16079 (2023)
  * Toma et al. [2023] Toma, A., Lawler, P.R., Ba, J., Krishnan, R.G., Rubin, B.B., Wang, B.: Clinical camel: An open-source expert-level medical language model with dialogue-based knowledge encoding. arXiv preprint arXiv:2305.12031 (2023)
  * Han et al. [2023] Han, T., Adams, L.C., Papaioannou, J.-M., Grundmann, P., Oberhauser, T., Löser, A., Truhn, D., Bressem, K.K.: Medalpaca–an open-source collection of medical conversational ai models and training data. arXiv preprint arXiv:2304.08247 (2023)
  * Dhakal et al. [2024] Dhakal, U., Singh, A.K., Devkota, S., Sapkota, Y., Lamichhane, B., Paudyal, S., Dhakal, C.: Gpt-4’s assessment of its performance in a usmle-based case study. arXiv preprint arXiv:2402.09654 (2024)
  * Jin et al. [2021] Jin, D., Pan, E., Oufattole, N., Weng, W.-H., Fang, H., Szolovits, P.: What disease does this patient have? a large-scale open domain question answering dataset from medical exams. Applied Sciences 11(14), 6421 (2021)
  * Vilares and Gómez-Rodríguez [2019] Vilares, D., Gómez-Rodríguez, C.: Head-qa: A healthcare dataset for complex reasoning. arXiv preprint arXiv:1906.04701 (2019)
  * Labrak et al. [2023] Labrak, Y., Bazoge, A., Dufour, R., Rouvier, M., Morin, E., Daille, B., Gourraud, P.-A.: Frenchmedmcqa: A french multiple-choice question answering dataset for medical domain. arXiv preprint arXiv:2304.04280 (2023)
  * Pal et al. [2022] Pal, A., Umapathi, L.K., Sankarasubbu, M.: Medmcqa: A large-scale multi-subject multi-choice dataset for medical domain question answering. In: Conference on Health, Inference, and Learning, pp. 248–260 (2022). PMLR
  * of Urogenital Tract Infection and Inflammation [2018] Urogenital Tract Infection, K.A., Inflammation: Guidelines for the Antibiotic Use in Urinary Tract Infections, (2018). {}{}}{https://www.uti.or.kr/include/pdf/antibiotic-use-2018.pdf}{cmtt}
  * Colgan et al. [2011] Colgan, R., Williams, M., Johnson, J.R.: Diagnosis and treatment of acute pyelonephritis in women. American family physician 84(5), 519–526 (2011)
  * Wolf et al. [2019] Wolf, T., Debut, L., Sanh, V., Chaumond, J., Delangue, C., Moi, A., Cistac, P., Rault, T., Louf, R., Funtowicz, M., et al.: Huggingface’s transformers: State-of-the-art natural language processing. arXiv preprint arXiv:1910.03771 (2019)
  * Gao et al. [2023] Gao, L., Tow, J., Abbasi, B., Biderman, S., Black, S., DiPofi, A., Foster, C., Golding, L., Hsu, J., Le Noac’h, A., Li, H., McDonell, K., Muennighoff, N., Ociepa, C., Phang, J., Reynolds, L., Schoelkopf, H., Skowron, A., Sutawika, L., Tang, E., Thite, A., Wang, B., Wang, K., Zou, A.: A framework for few-shot language model evaluation. Zenodo (2023). https://doi.org/10.5281/zenodo.10256836 . https://zenodo.org/records/10256836
  * Son et al. [2024] Son, G., Lee, H., Kim, S., Kim, S., Muennighoff, N., Choi, T., Park, C., Yoo, K.M., Biderman, S.: Kmmlu: Measuring massive multitask language understanding in korean. arXiv preprint arXiv:2402.11548 (2024)
  * Achiam et al. [2023] Achiam, J., Adler, S., Agarwal, S., Ahmad, L., Akkaya, I., Aleman, F.L., Almeida, D., Altenschmidt, J., Altman, S., Anadkat, S., et al.: Gpt-4 technical report. arXiv preprint arXiv:2303.08774 (2023)
  * Touvron et al. [2023] Touvron, H., Martin, L., Stone, K., Albert, P., Almahairi, A., Babaei, Y., Bashlykov, N., Batra, S., Bhargava, P., Bhosale, S., et al.: Llama 2: Open foundation and fine-tuned chat models. arXiv preprint arXiv:2307.09288 (2023)
  * Jiang et al. [2023] Jiang, A.Q., Sablayrolles, A., Mensch, A., Bamford, C., Chaplot, D.S., Casas, D.d.l., Bressand, F., Lengyel, G., Lample, G., Saulnier, L., et al.: Mistral 7b. arXiv preprint arXiv:2310.06825 (2023)
  * 01.AI [2024] 01.AI: Building the Next Generation of Open-source and Bilingual Llms., (2024). {}{}}{https://huggingface.co/01-ai/Yi-34B}{cmtt}
  * Kim et al. [2023] Kim, D., Park, C., Kim, S., Lee, W., Song, W., Kim, Y., Kim, H., Kim, Y., Lee, H., Kim, J., et al.: Solar 10.7 b: Scaling large language models with simple yet effective depth up-scaling. arXiv preprint arXiv:2312.15166 (2023)
  * L. Junbum, Taekyoon Choi [2023] L. Junbum, Taekyoon Choi: llama-2-koen-13b. Hugging Face (2023). https://doi.org/10.57967/hf/1280 . https://huggingface.co/beomi/llama-2-koen-13b
  * Lee Junbum [2024] Lee Junbum: Yi-Ko-6B (Revision 205083a). Hugging Face (2024). https://doi.org/10.57967/hf/1708 . https://huggingface.co/beomi/Yi-Ko-6B


# First-passage time theory of activated rate chemical processes in electronic
molecular junctions

Riley J. Preston  College of Science and Engineering, James Cook University,
Townsville, QLD, 4811, Australia  Maxim F. Gelin  School of Sciences, Hangzhou
Dianzi University, 310018 Hangzhou, China  Daniel S. Kosov  College of Science
and Engineering, James Cook University, Townsville, QLD, 4811, Australia

###### Abstract

Confined nanoscale spaces, electric fields and tunneling currents make the
molecular electronic junction an experimental device for the discovery of new,
out-of-equilibrium chemical reactions. Reaction-rate theory for current-
activated chemical reactions is developed by combining a Keldysh
nonequilibrium Green’s functions treatment of electrons, Fokker-Planck
description of the reaction coordinate, and Kramers’ first-passage time
calculations. The NEGF provide an adiabatic potential as well as a diffusion
coefficient and temperature with local dependence on the reaction coordinate.
Van Kampen’s Fokker-Planck equation, which describes a Brownian particle
moving in an external potential in an inhomogeneous medium with a position-
dependent friction and diffusion coefficient, is used to obtain an analytic
expression for the first-passage time. The theory is applied to several
transport scenarios: a molecular junction with a single, reaction coordinate
dependent molecular orbital, and a model diatomic molecular junction. We
demonstrate the natural emergence of Landauer’s blowtorch effect as a result
of the interplay between the configuration dependent viscosity and diffusion
coefficients. The resultant localized heating in conjunction with the bond-
deformation due to current-induced forces are shown to be the determining
factors when considering chemical reaction rates; each of which result from
highly tunable parameters within the system.

##  I INTRODUCTION

A molecular junction is a single molecule confined in the nanoscale gap
between two macroscopic, conducting leads. An applied voltage allows for the
flow of electronic current across the system through the valence states of the
molecule. Large current densities and power dissipation give way to strong
current-induced forces and bond-selective heating which act to destabilize the
molecular configuration within the junction, resulting in molecular
conformational changes including telegraphic switching Weick  _et al._ (2010);
Lü  _et al._ (2012); Preston  _et al._ (2020a); Kershaw and Kosov (2020),
along with providing the necessary energy for total bond ruptureDzhioev and
Kosov (2011); Dzhioev  _et al._ (2013); Erpenbeck  _et al._ (2018); Härtle and
Thoss (2011); Li  _et al._ (2015, 2016); Lü  _et al._ (2010). This is
obviously an undesirable feature for promoting molecular electronics as a
possible avenue for moving beyond the traditional silicon semiconductor
technology into a regime of highly efficient and tailorable molecular-scale
devices, and so a thorough theoretical understanding is required for further
progress. Nonetheless, it creates an exciting opportunity to explore and
produce new chemical reactions by providing a device which traps a single
molecule in a confined space of a few nanometers where the electric field and
current are applied locally and selectively Huang  _et al._ (2013); Aragonés
_et al._ (2016); Borca  _et al._ (2017).

The adequate and well established theories have been developed for reaction
rate calculations in gas and condensed phases Zwanzig (2001); Coffey  _et al._
(2004); Hänggi  _et al._ (1990); Melnikov (1991); Voth and Hochstrasser
(1996); Miller (1998); Pollak and Talkner (2005); Schüller  _et al._ (2020).
However, the development of similar theories for molecules in an electronic
junction environment is no simple task and as such, the scope of theoretical
work is still very limited. Three types of approaches have been proposed to
model current-induced dissociation. The first is based on the rate equation
approach where a single harmonic vibration is pumped beyond the dissociation
threshold limit Koch  _et al._ (2006); Brisker and Peskin (2008). The second
is a numerically exact scheme, which uses the hierarchical quantum master
equation method in conjunction with a discrete variable representation for the
nuclear degrees of freedom to numerically study current-induced dissociation
Erpenbeck  _et al._ (2020). The third uses Keldysh nonequilibrium Green’s
functions to obtain a Fokker-Planck equation for the reaction coordinate which
is used to compute average escape times and the accompanying reaction rates
Dzhioev and Kosov (2011); Dzhioev  _et al._ (2013). The further development of
this approach is the subject of this paper.

Our approach utilizes the Born-Oppenheimer approximation, in which nuclei
within the system are assumed to be slow-moving, classical particles
interacting with a sea of fast, quantum electrons. This separation of time-
scales enables the use of a Langevin description for the motion of nuclei, in
which the forces due to the quantum electrons are conceived through a
frictional force which acts to drive the nuclei into a motionless state, a
fluctuating force which re-energizes the nuclei, and an adiabatic force which
can deform the structure of the reaction potential. This enables the
consideration of highly non-trivial behaviour on the nuclear dynamics at the
cost of a fully quantum description. Nevertheless, the method has proven
successful in a range of circumstances Dou and Subotnik (2018, 2017); Dou  _et
al._ (2017); Bode  _et al._ (2011); Dzhioev  _et al._ (2013); Pistolesi  _et
al._ (2008); Weick  _et al._ (2010); Bode  _et al._ (2012); Kershaw and Kosov
(2020); Preston  _et al._ (2020a); Lü  _et al._ (2010); Preston  _et al._
(2020b).

The method further lends itself to the study of current-induced chemical
reaction rates Dzhioev and Kosov (2011); Dzhioev  _et al._ (2013). The use of
Langevin dynamics to compute reaction rates was first explored by Kramers in
his seminal 1940 paper Kramers (1940), in which the mean first-escape time of
a particle trapped in an arbitrary potential well subject to Langevin forces
was approximated. The next significant step was made in the 1990s, when
Kramers’ theory was extended to account for position-dependent friction and
generalized Langevin equations describing finite-memory (non-Markovian)
fluctuation-dissipation processes Pollak and Berezhkovskii (1993); Haynes  _et
al._ (1994); Haynes and Voth (1995); Neria and Karplus (1996). The effect of a
velocity-dependent friction on Kramers’ rates was investigated in Ref. Gelin
and Kosov (2007). However, these studies were limited to the regime of
thermodynamic equilibrium. Beyond this regime, the fluctuation-dissipation
theorem no longer holds, allowing for the emergence of localized heating
effects in analogy to Landauer’s proposed blowtorch effect Landauer (1975,
1993), in which specific configurations of the reaction coordinate may
experience heightened temperatures, which may have a significant affect on the
evolution of the system. Such systems are not limited to the realms of
molecular electronics; the most common examples include numerous molecular
motors, ratchets, and heat engines Reimann (2002); Hänggi and Marchesoni
(2009); Erbas-Cakmak  _et al._ (2015) as well as various confined nanosystems
Das and Ray (2015); Devine and Jack (2017); Franzese  _et al._ (2017); Holubec
_et al._ (2017); Arango-Restrepo and Rubi (2020), notably of biological
significance Basak  _et al._ (2019); Rubi (2019). Several explicit simulations
of Landauer’s blowtorch effects in double-well potentials have also been
performed recently Bekele  _et al._ (1999); Das  _et al._ (2015); Das and
Kantz (2019).

One of the aims of this paper is to further shed light on this topic. A
comprehensive understanding of the effects of localized heating on the
stability of molecular geometries is required to ensure the productive
development of specific functionalities of molecular-scale electronic systems.
In this paper, we relax the requirement of thermodynamic equilibrium, allowing
for the self-consistent study of the mean first-passage times in model
molecular electronic junctions in both the underdamped and overdamped regimes.
This is calculated through a Fokker-Planck equation, which arises due to our
Langevin description of the reaction coordinate within the junction. To a
certain extent this work is the continuation of two papers of one of the
authors (DSK) Dzhioev and Kosov (2011); Dzhioev  _et al._ (2013), however
Ref.Dzhioev and Kosov (2011) considered the problem employing the fluctuation-
dissipation theorem and Ref.Dzhioev  _et al._ (2013) focused on the
underdamped case only.

The paper is structured as follows. In section II, we demonstrate our
calculations for the mean first-passage time in the limiting regimes. This
involves the calculation of the current-induced forces in the system, from
which a Fokker-Planck description then yields an equation for the mean first-
passage times. This is then applied to a simple model of the blowtorch effect
in section III.1. In section III.2, we calculate the reaction rates for a
single-level junction model, in which current-induced forces are calculated
self-consistently within the model. This is then further applied to a model
two-level molecule within the junction in section III.3. Atomic units are used
in all calculations such that $\hbar=e=1$. Boltzmann’s constant is set to one
$k_{B}=1$ in all derivations and calculations meaning that the temperature is
measured in units of energy.

##  II THEORY

###  II.1 Hamiltonian

We begin with the general Hamiltonian describing a molecular junction as given
by

$H(t)={H}_{M}(t)+{H}_{L}+{H}_{R}+{H}_{LM}+{H}_{MR}.$ (1)

The total system Hamiltonian is partitioned into the following components; the
molecular Hamiltonian ${H}_{M}(t)$, the left and right leads Hamiltonians
${H}_{L}$ and ${H}_{R}$, and the leads-molecule coupling Hamiltonians
${H}_{LM}$ and ${H}_{MR}$ which describe the coupling between the electronic
states on the central molecule and the left and right leads, respectively.

The molecular Hamiltonian takes the form:

${H}_{M}(t)=\sum_{ij}h_{ij}(x(t))d^{{\dagger}}_{i}d_{j}+\frac{p^{2}}{2m}+U(x),$
(2)

where the operators $d^{{\dagger}}_{i}$ and $d_{j}$ denote the creation and
annihilation operators for the molecular electronic states whose energies are
given by the Hamiltonian matrix elements $h_{ij}(x(t))$. Note the explicit
time dependence here, which arises as a result of the time evolution of the
classical reaction coordinate $x$. The last two terms in $H_{M}$ are not
quantum mechanical operators, $p^{2}/2m$ is the kinetic energy for the motion
of the reaction coordinate and $U(x)$ is the classical potential energy.

The leads Hamiltonian is taken in the standard form of non-interacting
electrons reservoirs:

${H}_{L}+{H}_{R}=\sum_{k\alpha}\epsilon_{k\alpha}d^{\dagger}_{k\alpha}d_{k\alpha},$
(3)

where the creation and annihilation operators are given by
$d^{\dagger}_{k\alpha}$ and $d_{k\alpha}$, and the subscript $k\alpha$ denotes
the operator acting on state $k$ in the $\alpha$ lead which has energy
$\epsilon_{k\alpha}$.

Finally, the system-lead coupling Hamiltonians $H_{LM}$ and $H_{MR}$ are given
by:

${H}_{LM}+{H}_{MR}=\sum_{k\alpha i}\Big{(}t_{k\alpha
i}d^{\dagger}_{k\alpha}d_{i}+\text{h.c.}\Big{)}.$ (4)

The matrix elements $t_{k\alpha i}$ (and their conjugates) describe the
tunnelling amplitudes between lead states $k\alpha$ and the molecular orbitals
$i$.

###  II.2 Green’s Functions and Self-Energies

The foundational components from which we build our model for the classical
motion within a molecular junction are the adiabatic Green’s functions. The
required Green’s function components (advanced, retarded, lesser and
greater)Stefanucci and van Leeuwen (2013) can be solved for from the Keldysh-
Kadanoff-Baym equations, expressed in the Wigner space byHaug and Jauho (2010)

$\Big{(}\omega+\frac{i}{2}\partial_{t}-e^{\frac{1}{2i}\partial_{\omega}^{G}\partial_{t}^{h}}h(t)\Big{)}G^{R/A}(t,\omega)=I\\\
+e^{-\frac{1}{2i}\partial_{\omega}^{\Sigma}\partial_{t}^{\mathcal{G}}}\Sigma^{R/A}(\omega)G^{R/A}(t,\omega),$
(5)

and

$\Big{(}\omega+\frac{i}{2}\partial_{t}-e^{\frac{1}{2i}\partial_{\omega}^{G}\partial_{t}^{h}}h(t)\Big{)}G^{</>}(t,\omega)=\\\
e^{-\frac{1}{2i}\partial_{\omega}^{\Sigma}\partial_{t}^{G}}\Big{(}\Sigma^{R}(\omega)G^{</>}(t,\omega)+\Sigma^{</>}(\omega)G^{A}(t,\omega)\Big{)},$
(6)

where we have shown the retarded/advanced and lesser/greater equations
collectively. Here, we have already assumed that our molecule-leads coupling
components $t_{k\alpha i}$ are time-independent, which makes the self-energies
only $\omega$-dependent.

In alignment with our previous workKershaw and Kosov (2017, 2018, 2019, 2020);
Preston  _et al._ (2020a, b), we assume that the classical motion along the
reaction coordinate within the system occurs over long time-scales relative to
the characteristic electron tunnelling time. This provides us with the
required small parameter to be able to perturbatively solve (5) and (6) up to
the first order in expansion of the exponents with derivatives. Resultantly,
the adiabatic Green’s functions in the Wigner space (Green’s functions which
instantaneously follow the changes in the reaction coordinate) are given by
matrices of the form:

$G^{R/A}(x,\omega)=(\omega I-h(x)-\Sigma^{R/A})^{-1},$ (7)
$G^{</>}(x,\omega)=G^{R}(x,\omega)\Sigma^{</>}(\omega)G^{A}(x,\omega),$ (8)

where $I$ represents the identity operator in the molecular space. The
corresponding self-energy components are given by

$\Sigma_{\alpha,ij}^{R}=-\frac{i}{2}\Gamma_{\alpha,ij},$ (9)

$\Sigma_{\alpha,ij}^{A}=\frac{i}{2}\Gamma_{\alpha,ij},$ (10)

$\Sigma_{\alpha,ij}^{<}(\omega)=if_{\alpha}(\omega)\Gamma_{\alpha,ij},$ (11)

and

$\Sigma_{\alpha,ij}^{>}(\omega)=-i(1-f_{\alpha}(\omega))\Gamma_{\alpha,ij},$
(12)

where we have applied the wide-band approximation to the leads, eliminating
any $\omega$ dependence from $\Gamma$. As a result, the elements of the level
broadening function are given by

$\Gamma_{\alpha,ij}=2\pi\rho_{\alpha}t^{*}_{\alpha i}t_{\alpha j},$ (13)

where $\rho_{\alpha}$ is the constant, energy-independent density of states
for lead $\alpha$ and $t_{i\alpha}$ is the tunneling amplitude $t_{ik\alpha}$
which no longer depends on $k$.

We additionally solve for the first non-adiabatic correction to the Green’s
functions, which account for the non-zero velocity of the classical
coordinate. The Green’s function corrections are given by Bode  _et al._
(2012)

$G_{(1)}^{R/A}=\frac{1}{2i}G^{R/A}\Big{[}G^{R/A},\partial_{t}h\Big{]}_{-}G^{R/A},$
(14)

and

$G_{(1)}^{<}=G^{R}\Sigma^{<}G_{(1)}^{A}\ +G_{(1)}^{R}\Sigma^{<}G^{A}\\\ \
+\frac{1}{2i}G^{R}\Big{(}\partial_{T}hG^{R}\partial_{\omega}\Sigma^{<}+G^{<}\partial_{T}h+h.c.\Big{)}G^{A}.$
(15)

###  II.3 The Langevin Equation

Under the overarching assumption that the reaction coordinate $x$ along with
its corresponding momentum $p$ are classical variables in our approach due to
the separation of time-scales within the system, the equation of motion of the
reaction coordinate can be expressed in the form of a quasi-classical Langevin
equation, in which the classical motion is dictated by quantum mechanical
forces. Our Langevin equation is given by Bode  _et al._ (2012); Preston  _et
al._ (2020a); Pistolesi  _et al._ (2008); Weick  _et al._ (2010); Dzhioev  _et
al._ (2013); Kershaw and Kosov (2020)

$\frac{dp}{dt}=-\partial_{x}U(x)+F(x)-\frac{\xi(x)}{m}p+\delta f(t),$ (16)

in which we have the external classical potential $U(x)$, an adiabatic force
$F(x)$ arising due to the occupation of electronic levels within the molecule,
the frictional force and its corresponding viscosity coefficient $\xi(x)$, and
the stochastic force $\delta f(t)$. In our model, the adiabatic force takes
the form

$F(x)=\frac{i}{2\pi}\int
d\omega\text{Tr}\Big{\\{}\partial_{x}hG^{<}\Big{\\}}.$ (17)

Similarly, the viscosity coefficient is defined according toBode  _et al._
(2012); Dzhioev  _et al._ (2013)

$\frac{\xi(x)}{m}p=-\frac{i}{2\pi}\int
d\omega\text{Tr}\Big{\\{}\partial_{x}hG_{(1)}^{<}\Big{\\}}.$ (18)

It is assumed that the fluctuating force can be described by a Gaussian
process, defined by

$\langle\delta f(t)\rangle=0,$ (19)

and

$\langle\delta f(t)\delta f(t^{\prime})\rangle=D(x)\delta(t-t^{\prime}),$ (20)

where $D(x)$ is the quantum-mechanically calculated diffusion coefficient
which determines the variance in the fluctuating force. In our model, $D(x)$
is calculated according toBode  _et al._ (2012); Dzhioev  _et al._ (2013)

$D(x)=\frac{1}{2\pi}\int
d\omega\text{Tr}\Big{\\{}\partial_{x}hG^{>}\partial_{x}hG^{<}\Big{\\}}.$ (21)

###  II.4 Effective potential energy surface

The key quantity for chemical reaction rate calculations is the effective
potential energy surface. Here the effective nonequilibrium potential energy
surface is defined via the integration of the nonequilibrium force in the
Langevin equation

$U_{\text{eff}}(x)=U(x)-\int_{x_{0}}^{x}dy\;F(y).$ (22)

We chose $x_{0}$ in such a way that $U_{\text{eff}}(x)$ at the bottom of the
potential well, $x_{a}$, is zero, $U_{\text{eff}}(x_{a})=0$. This enforces
that the energy $E$ in the formulae below can vary in the range from 0 to
$\infty$. Suppose that the effective potential has the ”standard” Kramers
problem form, in which case we have a minimum at $x=x_{a}$ (reactant state)
and an energy barrier at $x=x_{b}$, separating the reactant and product state.

Let’s call

$U_{a}=U_{\text{eff}}(x_{a})=0$ (23)

and

$U_{b}=U_{\text{eff}}(x_{b}).$ (24)

###  II.5 Calculations of chemical reaction rates

####  II.5.1 Fokker-Planck equation with reaction-coordinate dependent
viscosity and diffusion

The starting point for all derivations of the first-passage times is van
Kampen’s Fokker-Planck equation for the probability density function
$\rho(x,p,t)$ which describes a Brownian particle of mass $m$, position $x$
and momentum $p$ moving in the external potential $U_{\text{eff}}(x)$ in an
inhomogeneous medium with position-dependent friction $\xi(x)$ and diffusion
coefficient $D(x)$ van Kampen (1988a)

$\partial_{t}\rho(x,p,t)=$
$\left(-\partial_{x}\frac{p}{m}+\partial_{p}U_{\text{eff}}^{{}^{\prime}}(x)+\left\\{\partial_{p}\frac{p}{m}\xi(x)+\partial_{p}^{2}D(x)\right\\}\right)\rho(x,p,t).$
(25)

Van Kampen’s Fokker-Planck equation (25) is internally consistent (for
example, it yields energy- and mass-flow transport coefficients obeying
Onsager’s reciprocity relations van Kampen (1991)) and can unambiguously be
derived from the Langevin equation (16) with Gaussian stochastic force obeying
(19) and (20). The Itô-Stratonovich dilemma van Kampen (1981) does not arise
in the derivation of van Kampen’s Fokker-Planck equation (25), because the
viscosity $\xi(x)$ of Eq. (18) and the diffusion coefficient $D(x)$ of Eq.
(21) are momentum-independent.

####  II.5.2 Underdamped limit

In this section, we present our equations for the mean first-passage times
from a reaction potential in an inhomogeneous medium for both the underdamped
and overdamped limiting regimes. Beginning with van Kampen’s Fokker-Planck
equation (25) for our reaction coordinate $x$ and momentum $p$, we can obtain
a one-dimensional energy-diffusion equation which is valid in the underdamped
limit (see Appendix A for the derivation):

$\partial_{t}\rho(E,t)=\partial_{E}D(E)\rho^{\text{un}}_{\text{st}}(E)\partial_{E}\left(\rho^{\text{un}}_{\text{st}}(E)\right)^{-1}\rho(E,t).$
(26)

Here, we have introduced the energy-diffusion $D(E)$ as

$D(E)=\frac{\nu(E)}{\Omega(E)},$ (27)

defined as the ratio of the two energy-dependent functions,

$\nu(E)=\sqrt{\frac{2}{m}}\int
dx\xi(x)T_{\text{eff}}(x)\sqrt{E-U_{\text{eff}}(x)},$ (28)

and

$\Omega(E)=\sqrt{\frac{m}{2}}\int\frac{dx}{\sqrt{E-U_{\text{eff}}(x)}},$ (29)

where we have defined the effective temperature $T_{\text{eff}}$ as

$T_{\text{eff}}(x)=\frac{D(x)}{2\xi(x)}.$ (30)

In Eqs. (28) and (29) and in all other similar expressions, the integration
domain corresponds to $E>U_{\text{eff}}(x)$.

The stationary distribution is written as

$\rho^{\text{un}}_{\text{st}}(E)=Z_{\text{un}}^{-1}\Omega(E)\exp\left\\{-\int_{0}^{E}dE^{\prime}\frac{\mu(E^{\prime})}{\nu(E^{\prime})}\right\\},$
(31)

with normalization constant

$Z_{\text{un}}=\int
dE\,\Omega(E)\exp\left\\{-\int_{0}^{E}dE^{\prime}\frac{\mu(E^{\prime})}{\nu(E^{\prime})}\right\\},$
(32)

and

$\mu(E)=\sqrt{\frac{2}{m}}\int dx\xi(x)\sqrt{E-U_{\text{eff}}(x)}.$ (33)

The energy diffusion equation can be used to calculate the mean first-passage
time in the underdamped regime as per

$\tau=\int_{U_{a}}^{U_{b}}dE^{\prime}\frac{1}{D(E^{\prime})\rho^{\text{un}}_{\text{st}}(E^{\prime})}\int_{0}^{E^{\prime}}dE^{\prime\prime}\rho^{\text{un}}_{\text{st}}(E^{\prime\prime}).$
(34)

####  II.5.3 Overdamped limit

In the overdamped limit, van Kampen’s Fokker-Planck equation (25) reduces to
the diffusion equation van Kampen (1988a, b) which can be cast into a form
similar to Eq. (26):

$\partial_{t}\rho(x,t)=\partial_{x}\frac{T_{\text{eff}}(x)}{\xi(x)}\rho^{\text{ov}}_{\text{st}}(x)\partial_{x}(\rho^{\text{ov}}_{\text{st}}(x))^{-1}\rho(x,t).$
(35)

The stationary distribution is given by

$\rho^{\text{ov}}_{\text{st}}(x)=\frac{1}{Z_{\text{ov}}T_{\text{eff}}(x)}\exp\left\\{-\int_{0}^{x}dx^{\prime}\frac{\partial_{x}U(x^{\prime})-F(x^{\prime})}{T_{\text{eff}}(x^{\prime})}\right\\},$
(36)

where

$Z_{\text{ov}}=\int_{-\infty}^{+\infty}dx\frac{1}{T_{\text{eff}}(x)}\exp\left\\{-\int_{0}^{x}dx^{\prime}\frac{\partial_{x}U(x^{\prime})-F(x^{\prime})}{T_{\text{eff}}(x^{\prime})}\right\\}.$
(37)

The first-passage time in the overdamped regime can be thus written as

$\tau=\int_{x_{a}}^{x_{b}}dx^{\prime}\frac{\xi(x^{\prime})}{T_{\text{eff}}(x^{\prime})\rho^{\text{ov}}_{\text{st}}(x^{\prime})}\int_{-\infty}^{x^{\prime}}dx^{\prime\prime}\rho^{\text{ov}}_{\text{st}}(x^{\prime\prime}).$
(38)

To summarize the derivations in section II: the main result is the combined
use of the adiabatic force (17), the position-dependent viscosity (18) and the
position-dependent diffusion coefficient (21) - all obtained from
nonequilibrium Green’s function calculations - with exact (in underdamped and
overdamped limits) expressions (34,38) for the mean first-passage time $\tau$.
The rates for corresponding chemical reactions can obtained as inverse of the
first-passage time $1/\tau$.

##  III RESULTS

###  III.1 The blowtorch effect

In this section we investigate Landauer’s proposed blowtorch effectLandauer
(1993), in which a non-equilibrium system allows for coordinate-dependent
variations to the dissipative forces acting on a particle, which then has an
effect on the properties of the steady-state distribution. Landauer’s
blowtorch effect plays a critical role in chemical reactions in molecular
electronic junctions, therefore we first discuss its general features which
will be relevant for our subsequent discussion. For consistency with Kramers’
seminal paperKramers (1940) on chemical reaction rates, our analysis is
formulated in terms of the mean first escape time from the left minimum of a
bistable potential, which we model according to a quartic of the form

$U(x)=-ax^{2}+bx^{4},$ (39)

where $x$ is our reaction coordinate. The constants $a$ and $b$ are adjustable
parameters which determine the width and depth of the minimum. In all tests in
this section, we set $a=0.04$ and $b=0.008$, such that $U_{b}=0.05$. In
addition, the particle always begins its trajectory with zero velocity at the
minimum of the reaction potential. To begin with, the viscosity $\xi_{0}$ and
diffusion coefficients $D_{0}$ are set to a constant value over the range of
$x$, yielding a constant temperature as determined by the fluctuation-
dissipation theorem.

In order to introduce an inhomogeneity into the temperature, a Gaussian spike
is applied to the diffusion coefficient locally at position $x_{0}$,

$D(x)=D_{0}+D_{peak}(x),$ (40)

where

$D_{peak}(x)=D_{m}e^{-\frac{(x-x_{0})^{2}}{\sigma^{2}}},$ (41)

with adjustable width $\sigma$ and magnitude $D_{m}$ parameters. The effective
temperature profile is given by (30), the effective temperature at the peak is
then given by

$T_{max}=T_{0}(1+\frac{D_{m}}{D_{0}}),$ (42)

where

$T_{0}=\frac{D_{0}}{2\xi_{0}}.$ (43)

This represents Laundauer’s so-called ”blowtorch” which heats a small segment
of the reaction coordinate, as shown diagramatically in Fig.1.

Figure 1: An adjustable temperature spike is introduced which heats a chosen
part of our reaction coordinate potential.

Here, the intention is to study the effect of shifting the temperature spike
along the reaction coordinate on the mean first-passage time $\tau$. We
analyse the overdamped and underdamped regimes separately for the same
parameters except for the mass $m$ of a Brownian particle, which is chosen to
satisfy the desired regime.

(a)

(b)

Figure 2: The calculated mean first-passage time $\tau$ as a function of the
temperature peak’s position along the reaction coordinate, for the (a)
overdamped and (b) underdamped regime ($m=1000$a.u). The black line denotes
$\tau$ in the absence of an applied blowtorch. Parameters: $D_{0}=0.01$,
$\xi_{0}=1$, $\sigma=0.05$.

In Fig.2, we observe the behaviour of the mean first-passage time as the
position of the temperature peak is shifted along the reaction coordinate
(shown in blue), while the reaction potential is shown as a reference in
orange. In the underdamped regime, we observe $\tau$ to be minimized when the
heating is applied to the bottom of the potential. This enables the molecule
to heat up quickly at low energies, and repeatedly attain more energy as it
passes through this region in a near-harmonic manner.

The overdamped regime differs, in that $\tau$ is minimized when the heating is
applied approximately halfway up the potential, around the point of steepest
ascent. In the overdamped regime, the escaping particle will very quickly
equilibrate to any given temperature fluctuation to which it is exposed. As
such, the heated region causes a flattening of the probability distribution in
that region, nullifying the dependence of the distribution on the reaction
coordinate. This causes an effective reduction to the height of the energy
barrier $U_{b}$ as elucidated by LandauerLandauer (1975, 1993); a phenomenon
which is maximized when the heating is applied in the region of steepest
ascent. We note the counter-intuitive observation that if the heating is
applied to the bottom of the potential in the overdamped case, this causes
only a small reduction to $\tau$. This is because the particle will quickly
lose the obtained energy as it returns to the cooler regions when it attempts
to escape.

It is insightful for us to also study the effect of the strength of
interaction of a Brownian particle with the environment, while maintaining a
homogeneous temperature. This entails that any changes in the diffusion
coefficient as a function of $x$ will be counteracted by a corresponding
changes in the viscosity coefficient at the same $x$, enforcing a homogeneous
temperature as per the fluctuation-dissipation theorem. Here, we perform a
similar analysis as above, such that we have a moveable peak of increased
interaction (simultaneously locally increased diffusion and viscosity) while
the temperature is homogeneous. The results of this are displayed in Fig.3.

(a)

(b)

Figure 3: The calculated mean first-passage time as a function of the
interaction peak’s position along the reaction coordinate, for the (a)
overdamped and (b) underdamped regime ($m_{e}=1000$a.u). The black line
denotes $\tau$ in the absence of an applied blowtorch. Parameters:
$D_{0}=0.01$, $\xi_{0}=1$, $\sigma=0.05$.

In the underdamped regime, we observe that the largest reduction to $\tau$
occurs when the interaction peak is placed at the minimum. The decrease in
$\tau$ agrees with the homogeneous-case solution, with the distinction that
reaction coordinates at higher energies in the potential have diminishing
contributions to the decreasing $\tau$. In the overdamped regime, it is seen
that the interaction peak results in an increase to $\tau$ as also predicted
in the homogeneous case. However, we observe that this is dominated by the
increased interaction strength near to the maximum of the potential, while
changing the interaction strength in the rest of the potential has negligible
effect. This demonstrates that $\tau$ has little regard for the interaction
strength in the overdamped regime, except in the region approaching the
maximum.

This general analysis arms us with the required physical intution before
proceeding to the next section, in which we first observe how a Landauer
blowtorch emerges naturally from a simple molecular junction model, then
demonstrate the effect on hypothetical chemical reaction rates.

###  III.2 Application to a molecule with a single current-carrying molecular
orbital

In this section, we analyze the calculated mean first-passage time $\tau$ for
our model of a molecular junction. Contrary to section III.1, the viscosity
and diffusion coefficients will be computed using nonequilibrium Green’s
functions according to eqs. (18) and (21). We consider the case of a single
electronic level coupled to the left and right leads under some applied bias
voltage. The molecular Hamiltonian is simplified to

$H_{M}=(h_{0}+\lambda x)d^{{\dagger}}d,$ (44)

where subscript $d^{\dagger}$ and $d^{\dagger}$ denote the creation and
annihilation operators for an electron on the molecular orbital, while all
previous matrix quantities are simplified to scalars in this regime. Here, the
dependency of $H_{M}$ on the reaction coordinate acts to shift the electronic
level, as scaled by the tuneable parameter $\lambda$. The left and right lead
are each at room temperature (0.00095a.u) and are symmetrically coupled to the
central electronic state such that our level-width function is given by
$\Gamma_{L}=\Gamma_{R}=0.03$.

In the interest of consistency, we again utilize the same quartic to describe
our classical nuclear potential for the reaction coordinate, which is now
acted on by an additional adiabatic force term computed using nonequilibrium
Green’s functions according to equation 17. This has the effect of shifting
and shallowing/deepening the reaction potential depending on the parameters
chosen. We also allow for the manual shift of the external potential along the
reaction coordinate according to some parameter $x_{a}$:

$U(x)=-a(x-x_{a})^{2}+b(x-x_{a})^{4}.$ (45)

This means that when $x_{a}=0$, the potential minimum (ignoring the effects of
$F$) occurs at $x=0$, while a positive $x_{a}$ shifts the input potential to
the right. Any bias voltage is applied symmetrically, such that
$\mu_{L}=-\mu_{R}=V/2$, where $\mu_{L}$ and $\mu_{R}$ are the chemical
potentials of the left and right leads.

(a)

(b)

(c)

(d)

Figure 4: The effect of an applied gate voltage to (a) the viscosity
coefficient and (b) the effective temperature. The mean first-passage time
$\tau$ in the (c) overdamped and (d) underdamped regime is plotted against the
peak coordinate of the viscosity and temperature (as determined by the applied
gate voltage) for different $\lambda$. The coordinates of the minimum and
maximum of the reaction potential are denoted by the vertical black lines in
(c) and (d).

We study the effect of applying a gate voltage to the system, as modelled by a
shift in the $h_{0}$ value. This allows for a degree of controllability of the
reaction rates for a given system. Figs.4(a) and 4(b) demonstrate the
resultant viscosity coefficient and temperature respectively, as a function of
the reaction coordinate. Application of a gate voltage acts to shift the curve
along the reaction coordinate. This analysis is performed for a non-zero bias
voltage such that the temperature is now inhomogeneous in addition to the
viscosity. In the underdamped regime shown in Fig.4(d), $\tau$ is minimized
when the viscosity and temperature peaks are shifted near to the minimum of
the reaction potential (note however, that the minimum in $\tau$ does not
occur exactly when the peaks are shifted to the minimum of the potential due
to the slight asymmetry of the reaction potential). In contrast, the
overdamped regime displays highly non-trivial behaviour, arising as a result
of the interplay between the strength of the viscosity and the temperature. In
our analysis of the overdamped regime in the previous section, we noted that
the dependence of $\tau$ on the temperature is dominated by the region of
steepest ascent up towards the maximum. Here, we again observe this behaviour
as the large peak in Fig.4(c) corresponds to when the dip in the temperature
occurs in this region (when the temperature peak has been shifted to the
right). A corresponding but smaller peak also occurs due to a shift to the
left in the temperature such that the low temperature aligns with the steep
region of the potential. The difference in peak sizes arises as a result of
the inhomogeneous viscosity, which per the previous section, we know is
important in the region near the maximum of the reaction potential. The large
peak in $\tau$ occurs when the temperature is low in the steep region, while
the viscosity is high towards the maximum. The small peak has a low viscosity
near the maximum, explaining its comparatively smaller magnitude.

###  III.3 Model of a two-level molecule

In this section, we expand the model to consider a two-level system. The
second molecular orbital is not a simple addition of an extra level here,
since the Green’s functions, self-energies and molecular Hamiltonian become
$2\times 2$ matrices and some nontrivial terms such as the commutator in (14)
will no longer be zero. In our model, the molecular energy levels is taken to
correspond to the bonding and anti-bonding states of a free $H_{2}^{+}$
moleculeMcQuarrie (2007). As such, the molecular Hamiltonian now reads

$H_{M}(t)=\sum_{ij}h_{ij}(q(t))d^{{\dagger}}_{i}d_{j},$ (46)

where $d_{i}^{{\dagger}}$ and $d_{j}$ are now in the molecular orbital basis.
The electronic Hamiltonian elements are represented in the form of a $2\times
2$ matrix according to

$h=\begin{pmatrix}H_{b}(q)&0\\\ 0&H_{a}(q)\end{pmatrix},$ (47)

where we use $H_{b}(q)$ and $H_{a}(q)$ to denote the bonding and anti-bonding
molecular orbitals, respectively, while $q$ is the bond-length. The values for
$H_{b}(q)$ and $H_{a}(q)$ are calculated according to molecular orbital
theoryMcQuarrie (2007) and are given by

$H_{b}(q)=\frac{H_{AA}+H_{AB}}{1+S_{AB}},$ (48)

and

$H_{a}(q)=\frac{H_{AA}-H_{AB}}{1-S_{AB}},$ (49)

where $H_{AA}$ and $H_{AB}$ are the Hamiltonian elements in the atomic basis
and $S_{AB}$ is the overlap integral between atomic 1s Slater orbitals. The
constituent components are then given by

$H_{AA}=-\frac{1}{2}+e^{-2q}\Big{(}1+\frac{1}{q}\Big{)}-\frac{1}{q},$ (50)

(a)

(b)

Figure 5: The adiabatic potential as a function of the bond-length presented
for (a) varying bias voltages and (b) varying the magnitude of leads coupling
to $H_{a}$. Parameters: $V=0$, $\Gamma_{aa}/\Gamma_{bb}=1$, unless otherwise
specified. $\Gamma_{bb}=0.03$ in all calculations.

(a)

(b)

Figure 6:  The effect of varying the bias voltage is shown for the (a)
viscosity coefficient and (b) the effective temperature, as a function of the
bond length. Insets: Shows the same quantity at $V=0.02$ for
$\Gamma_{aa}/\Gamma_{bb}=0$ (dashed) and $\Gamma_{aa}/\Gamma_{bb}=0.5$
(solid). Parameters; $\Gamma_{aa}/\Gamma_{bb}=1$ in the main plots.
$\Gamma_{bb}=0.03$ in all calculations.

(a)

(b)

Figure 7: The mean first-passage time $\tau$ as a function of the bias
voltage, varying the coupling to $H_{a}$ in the (a) overdamped and (b)
underdamped cases.

$H_{AB}=-\frac{S_{AB}(q)}{2}-e^{-q}(1+q),$ (51)

and

$S_{AB}=e^{-q}(1+q+q^{2}/3).$ (52)

In the interest of simplicity, each of the molecular orbitals is symmetrically
coupled to the left and right leads; as controlled by the parameter $\Gamma$
which now takes the form of a matrix as per

$\Gamma_{\alpha}=\begin{pmatrix}\Gamma_{\alpha,bb}&\Gamma_{\alpha,ba}\\\
\Gamma_{\alpha,ab}&\Gamma_{\alpha,aa}\end{pmatrix}$ (53)

for the $\alpha$ lead, where the off-diagonal components can be defined
according to,

$\Gamma_{\alpha,ba}=\Gamma_{\alpha,ab}=\sqrt{\Gamma_{\alpha,aa}\Gamma_{\alpha,bb}}.$
(54)

In each test, we have $\mu_{L}=-0.7$ and $\mu_{R}=\mu_{L}-V$, while the lead
temperatures are again set to room temperature.

The external potential now represents the classical nuclear repulsion, which
in atomic units is given by

$U(q)=\frac{1}{q}.$ (55)

Inclusion of the electronic forces allows us to generate modified electronic
potentials for varied parameters in order to assess the molecular stability.
Examples of these potentials are shown in Fig.5; where in (a) an applied bias
voltage is shown to decrease the energy required for bond rupture, while (b)
shows the effect of the additional electronic level which when occupied, acts
to increase the bond stability.

Along with the shape of the effective potential, the bond stability is also
determined by the electronic viscosity and effective temperature, which are
demonstrated in Fig.6. In the viscosity coefficient, each curve shows a peak
at small $q$, which approximately corresponds to when $H_{a}$ crosses the
fermi-level of the left lead. Likewise, the peaks at large $q$ are a result of
$H_{b}$ crossing the fermi-level of the right, then left, leads (these split
peaks merge together when $V=0$). The inset plot demonstrates the effect of
allowing an additional transport channel through the excited state which not
only introduces the peak at small $q$, but also increases the magnitude of the
viscous forces overall. The effective temperature is equal to the leads
temperature for $V=0$, while non-zero voltages yield a complex array of
localized heating and cooling effects, which arise as the energy levels shift
in and out of the resonance region as the bond-length is increased.

These competing effects culminate in our calculation of the mean first-passage
time, which is demonstrated in Fig.7 as a function of the bias voltage, for
different coupling values to the excited electronic state. In both limiting
regimes, an increase to the bias voltage acts to destabilize the bond and
decrease $\tau$, both due to the increased effective temperatures and the
weakening of the bond due to the current-induced forces. In the overdamped
regime, allowing the leads to be coupled to an additional level in the central
region has a stabilizing effect for all voltages tested, increasing the
average amount of time for bond rupture. The underdamped case shows similar
behavior for very low voltages; however, at higher voltages the availability
of the additional transport channel through the excited state increases the
current-induced forces such that the energy required for a bond-rupture is
found more easily, decreasing $\tau$.

##  IV CONCLUSIONS

In this paper, we have demonstrated that the rates of chemical reactions for
molecules in electronic junctions depend on three crucial ingredients; the
potential energy surface which defines the energy required for a configuration
change or bond rupture, the rate of the energy removal from vibrational to
electronic degrees of freedom given by the electronic viscosity coefficient,
and lastly, the effective temperature dynamically established in the molecule.
While the magnitude of these quantities is of high importance, the local
distribution of the viscosity and effective temperature along the potential
energy surface (Landauer’s blowtorch effect) also proves to be critical.

The addition of localized heating and cooling effects as a result of
inhomogeneity with respect to the molecular configuration has been shown to
induce significant variations in the mean first-passage time, as calculated
according to a Fokker-Planck description obtained by separating slow (reaction
coordinate) and fast (electronic) time-scales in the Keldysh-Kadanoff-Baym
equations for the nonequilibrium Green’s functions. This has been demonstrated
for a single-level molecular junction model, as well as a two-level model
inspired by $H_{2}^{+}$ molecular orbitals with the bond length considered as
the reaction coordinate. This interplay between the amount of energy required
for bond rupture and the energy supplied due to tunneling electrons has been
shown to be strongly dependent on the choice of experimentally tuneable
parameters for the system. This enables the possibility of a high degree of
controllability for molecular junction systems, with promises of controlled
initiation of chemical reactions or conversely, enforcing the stability of
specific configurations within the system.

DATA AVAILABILITY

The data that supports the findings of this study are available within the
article.

##  Appendix A Derivation of the energy diffusion equation in the underdamped
limit

Following Zwanzig Zwanzig (2001), we can introduce the projection operator

$\widehat{O}=\Omega^{-1}(E)\int dxdp\,\,\delta(E-H)$ (56)

where

$H=\frac{p^{2}}{2m}+U_{\text{eff}}(x)$ (57)

is the Hamiltonian of the bath-free Brownian particle and the microcanonical
partition function is defined as

$\Omega(E)=\int dxdp\,\,\delta(E-H)$ (58)

(equivalently, $\Omega(E)$ is determined via Eq. (29)). Note that
$\widehat{O}$ has an extra factor of $\Omega^{-1}(E)$ in comparison with the
operator introduced by Zwanzig. With our definition,
$\widehat{O}^{2}=\widehat{O}$. Applying the projection operator of Eq. (56) to
the Fokker-Planck equation (25), exploiting the identity

$\widehat{O}\left(-\partial_{x}\frac{p}{m}+\partial_{p}U_{\text{eff}}^{{}^{\prime}}(x)\right)=0$

and taking into account that $\xi(x)\ll 1$ (underdamped limit) we obtain

$\partial_{t}\widehat{O}\rho(x,p,t)=\widehat{O}\zeta(x)\left\\{\partial_{p}\frac{p}{m}+\partial_{p}^{2}T(x)\right\\}\widehat{O}\rho(x,p,t).$
(59)

Changing the differentiation variables ($\partial_{p}=m^{-1}\partial_{E}p$)
and using the identity
$\partial_{p}^{2}\widehat{O}=m^{-1}\partial_{p}p\partial_{E}\widehat{O}$, we
obtain the following equation for $\rho(E,t)\equiv\Omega(E)\widehat{O}\rho(t)$
(cf. Zwanzig (2001); Gelin and Kosov (2007)):

$\partial_{t}\rho(E,t)=\partial_{E}\left\\{\mu(E)+\nu(E)\partial_{E}\right\\}\Omega^{-1}(E)\rho(E,t).$
(60)

Here $\mu(E)$ and $\nu(E)$ are defined through Eqs. (33) and (28), and
$\rho(E,t)$ is normalized according to $\int dE\rho(E,t)=1$.

## References

  * Weick  _et al._ (2010) G. Weick, F. Pistolesi, E. Mariani, and F. von Oppen, Phys. Rev. B 81, 121409 (2010).
  * Lü  _et al._ (2012) J.-T. Lü, M. Brandbyge, P. Hedegård, T. N. Todorov, and D. Dundas, Phys. Rev. B 85, 245444 (2012).
  * Preston  _et al._ (2020a) R. J. Preston, V. F. Kershaw, and D. S. Kosov, Phys. Rev. B 101, 155415 (2020a).
  * Kershaw and Kosov (2020) V. F. Kershaw and D. S. Kosov, The Journal of Chemical Physics 153, 154101 (2020), https://doi.org/10.1063/5.0023275 .
  * Dzhioev and Kosov (2011) A. A. Dzhioev and D. S. Kosov, The Journal of Chemical Physics 135, 074701 (2011), https://doi.org/10.1063/1.3626521 .
  * Dzhioev  _et al._ (2013) A. A. Dzhioev, D. S. Kosov, and F. von Oppen, The Journal of Chemical Physics 138, 134103 (2013), https://doi.org/10.1063/1.4797495 .
  * Erpenbeck  _et al._ (2018) A. Erpenbeck, C. Schinabeck, U. Peskin, and M. Thoss, Phys. Rev. B 97, 235452 (2018).
  * Härtle and Thoss (2011) R. Härtle and M. Thoss, Phys. Rev. B 83, 125419 (2011).
  * Li  _et al._ (2015) H. Li, T. A. Su, V. Zhang, M. L. Steigerwald, C. Nuckolls, and L. Venkataraman, Journal of the American Chemical Society 137, 5028 (2015), pMID: 25675085, https://doi.org/10.1021/ja512523r .
  * Li  _et al._ (2016) H. Li, N. T. Kim, T. A. Su, M. L. Steigerwald, C. Nuckolls, P. Darancet, J. L. Leighton, and L. Venkataraman, Journal of the American Chemical Society 138, 16159 (2016), pMID: 27960303, https://doi.org/10.1021/jacs.6b10700 .
  * Lü  _et al._ (2010) J.-T. Lü, M. Brandbyge, and P. Hedegård, Nano Letters 10, 1657 (2010), https://doi.org/10.1021/nl904233u .
  * Huang  _et al._ (2013) K. Huang, L. Leung, T. Lim, Z. Ning, and J. C. Polanyi, Journal of the American Chemical Society 135, 6220 (2013).
  * Aragonés  _et al._ (2016) A. C. Aragonés, N. L. Haworth, N. Darwish, S. Ciampi, N. J. Bloomfield, G. G. Wallace, I. Diez-Perez, and M. L. Coote, Nature 531, 88 (2016).
  * Borca  _et al._ (2017) B. Borca, T. Michnowicz, R. Pétuya, M. Pristl, V. Schendel, I. Pentegov, U. Kraft, H. Klauk, P. Wahl, R. Gutzler, A. Arnau, U. Schlickum, and K. Kern, ACS Nano 11, 4703 (2017), pMID: 28437066, https://doi.org/10.1021/acsnano.7b00612 .
  * Zwanzig (2001) R. Zwanzig,  _Nonequilibrium statistical mechanics_ (Oxford University Press, Oxford, 2001).
  * Coffey  _et al._ (2004) W. T. Coffey, Y. P. Kalmykov, and J. T. Waldron,  _The Langevin Equation_ (World Scientific, Singapore, 2004).
  * Hänggi  _et al._ (1990) P. Hänggi, P. Talkner, and M. Borcovec, Reviews of Modern Physics 62, 251 (1990), https://doi.org/10.1063/1.527963 .
  * Melnikov (1991) V. I. Melnikov, Phys. Rep. 209, 1 (1991), https://doi.org/10.1063/1.527963 .
  * Voth and Hochstrasser (1996) G. A. Voth and R. M. Hochstrasser, J. Phys. Chem. 100, 13034 (1996), https://doi.org/10.1063/1.527963 .
  * Miller (1998) W. H. Miller, J. Phys. Chem. A 102, 793 (1998), https://doi.org/10.1063/1.527963 .
  * Pollak and Talkner (2005) E. Pollak and P. Talkner, Chaos 15, 026116 (2005), https://doi.org/10.1063/1.527963 .
  * Schüller  _et al._ (2020) B. Schüller, A. Meistrenko, H. Hees, Z. Xu, and C. Greiner, Annals of Physics 412, 168045 (2020), https://doi.org/10.1063/1.527963 .
  * Koch  _et al._ (2006) J. Koch, M. Semmelhack, F. von Oppen, and A. Nitzan, Phys. Rev. B 73, 155306 (2006).
  * Brisker and Peskin (2008) D. Brisker and U. Peskin, The Journal of Chemical Physics 129, 244709 (2008), https://doi.org/10.1063/1.3021288 .
  * Erpenbeck  _et al._ (2020) A. Erpenbeck, Y. Ke, U. Peskin, and M. Thoss, Phys. Rev. B 102, 195421 (2020).
  * Dou and Subotnik (2018) W. Dou and J. E. Subotnik, Phys. Rev. B 97, 064303 (2018).
  * Dou and Subotnik (2017) W. Dou and J. E. Subotnik, Phys. Rev. B 96, 104305 (2017).
  * Dou  _et al._ (2017) W. Dou, G. Miao, and J. E. Subotnik, Phys. Rev. Lett. 119, 046001 (2017).
  * Bode  _et al._ (2011) N. Bode, S. V. Kusminskiy, R. Egger, and F. von Oppen, Phys. Rev. Lett. 107, 036804 (2011).
  * Pistolesi  _et al._ (2008) F. Pistolesi, Y. M. Blanter, and I. Martin, Phys. Rev. B 78, 085127 (2008).
  * Bode  _et al._ (2012) N. Bode, S. V. Kusminskiy, R. Egger, and F. von Oppen, Beilstein J. Nanotechnol 3, 144 (2012).
  * Preston  _et al._ (2020b) R. J. Preston, T. D. Honeychurch, and D. S. Kosov, The Journal of Chemical Physics 153, 121102 (2020b), https://doi.org/10.1063/5.0019178 .
  * Kramers (1940) H. Kramers, Physica 7, 284 (1940).
  * Pollak and Berezhkovskii (1993) E. Pollak and A. M. Berezhkovskii, J. Chem. Phys. 99, 1344 (1993), https://doi.org/10.1063/1.527963 .
  * Haynes  _et al._ (1994) G. R. Haynes, E. Pollak, and G. A. Voth, J. Chem. Phys. 101, 7811 (1994), https://doi.org/10.1063/1.527963 .
  * Haynes and Voth (1995) G. R. Haynes and G. A. Voth, J. Chem. Phys. 103, 10176 (1995), https://doi.org/10.1063/1.527963 .
  * Neria and Karplus (1996) E. Neria and M. Karplus, J. Chem. Phys. 105, 10812 (1996), https://doi.org/10.1063/1.527963 .
  * Gelin and Kosov (2007) M. F. Gelin and D. K. Kosov, J. Chem. Phys. 126, 244501 (2007), https://doi.org/10.1063/1.527963 .
  * Landauer (1975) R. Landauer, Phys. Rev. A 12, 636 (1975).
  * Landauer (1993) R. Landauer, Physica A: Statistical Mechanics and its Applications 194, 551 (1993).
  * Reimann (2002) P. Reimann, Phys. Rep. 361, 57 (2002), https://doi.org/10.1063/1.527963 .
  * Hänggi and Marchesoni (2009) P. Hänggi and F. Marchesoni, Rev. Mod. Phys. 81, 387 (2009), https://doi.org/10.1063/1.527963 .
  * Erbas-Cakmak  _et al._ (2015) S. Erbas-Cakmak, D. A. Leigh, C. T. McTernan, and A. L. Nussbaumer, Chem. Rev. 115, 10081 (2015), https://doi.org/10.1063/1.527963 .
  * Das and Ray (2015) M. Das and D. S. Ray, Phys. Rev. E 92, 052133 (2015), https://doi.org/10.1063/1.527963 .
  * Devine and Jack (2017) J. Devine and M. W. Jack, Phys. Rev. E 96, 062130 (2017), https://doi.org/10.1063/1.527963 .
  * Franzese  _et al._ (2017) G. Franzese, I. Latella, and J. M. Rubi, Entropy 19, 507 (2017), https://doi.org/10.1063/1.527963 .
  * Holubec  _et al._ (2017) V. Holubec, A. Ryabov, M. H. Yaghoubi, M. Varga, A. Khodaee, M. E. Foulaadvand, and P. Chvosta, Entropy 19, 119 (2017), https://doi.org/10.1063/1.527963 .
  * Arango-Restrepo and Rubi (2020) A. Arango-Restrepo and J. M. Rubi, J. Chem. Phys. 153, 034108 (2020), https://doi.org/10.1063/1.527963 .
  * Basak  _et al._ (2019) S. Basak, S. Sengupta, and K. Chattopadhyay, Bioph. Rev. 11, 851 (2019), https://doi.org/10.1063/1.527963 .
  * Rubi (2019) J. M. Rubi, EPL 127, 10001 (2019), https://doi.org/10.1063/1.527963 .
  * Bekele  _et al._ (1999) M. Bekele, S. Rajesh, G. Ananthakrishna, and N. Kumar, Phys. Rev. E 59, 143 (1999), https://doi.org/10.1063/1.527963 .
  * Das  _et al._ (2015) M. Das, D. Das, D. Barik, and D. S. Ray, Phys. Rev. E 19, 052102 (2015), https://doi.org/10.1063/1.527963 .
  * Das and Kantz (2019) M. Das and H. Kantz, Phys. Rev. E 100, 032108 (2019), https://doi.org/10.1063/1.527963 .
  * Stefanucci and van Leeuwen (2013) G. Stefanucci and R. van Leeuwen,  _Nonequilibrium Many-Body Theory of Quantum Systems: A Modern Introduction_ (Cambridge University Press, 2013).
  * Haug and Jauho (2010) H. Haug and A. Jauho,  _Quantum Kinetics in Transport and Optics of Semiconductors_ (Springer, Berlin/Heidelberg, 2010).
  * Kershaw and Kosov (2017) V. F. Kershaw and D. S. Kosov, J. Chem. Phys. 147, 224109 (2017), https://doi.org/10.1063/1.5007071 .
  * Kershaw and Kosov (2018) V. F. Kershaw and D. S. Kosov, J. Chem. Phys. 149, 044121 (2018), https://doi.org/10.1063/1.5028333 .
  * Kershaw and Kosov (2019) V. F. Kershaw and D. S. Kosov, Journal of Chemical Physics 150, 074101 (2019), 1809.07140 .
  * van Kampen (1988a) N. G. van Kampen, IBM J. Res. Develop. 32, 107 (1988a), https://doi.org/10.1063/1.527963 .
  * van Kampen (1991) N. G. van Kampen, Journal of Statistical Physics 63, 1019 (1991), https://doi.org/10.1063/1.527963 .
  * van Kampen (1981) N. G. van Kampen, Journal of Statistical Physics 24, 175 (1981), https://doi.org/10.1063/1.527963 .
  * van Kampen (1988b) N. G. van Kampen, Journal of Mathematical Physics 29, 1220 (1988b), https://doi.org/10.1063/1.527963 .
  * McQuarrie (2007) D. A. McQuarrie,  _Quantum Chemistry_ (University Science Books, Sausalito, California, 2007).