11institutetext: Sano Centre for Computational Medicine, Cracow, Poland
22institutetext: Warsaw School of Economics, Warsaw, Poland 33institutetext:
Massachusetts General Hospital, Harvard Medical School, Boston, MA, USA

# Estimation of the Impact of COVID-19 Pandemic Lockdowns on Breast Cancer
Deaths and Costs in Poland using Markovian Monte Carlo Simulation

Magdalena Dul 1122 Michal K. Grzeszczyk 11 0000-0002-5304-1020 Ewelina
Nojszewska 22 0000-0003-3176-0240 Arkadiusz Sitek 33 0000-0002-0677-4002

###### Abstract

This study examines the effect of COVID-19 pandemic and associated lockdowns
on access to crucial diagnostic procedures for breast cancer patients,
including screenings and treatments. To quantify the impact of the lockdowns
on patient outcomes and cost, the study employs a mathematical model of breast
cancer progression. The model includes ten different states that represent
various stages of health and disease, along with the four different stages of
cancer that can be diagnosed or undiagnosed. The study employs a natural
history stochastic model to simulate the progression of breast cancer in
patients. The model includes transition probabilities between states,
estimated using both literature and empirical data. The study utilized a
Markov Chain Monte Carlo simulation to model the natural history of each
simulated patient over a seven-year period from 2019 to 2025. The simulation
was repeated 100 times to estimate the variance in outcome variables. The
study found that the COVID-19 pandemic and associated lockdowns caused a
significant increase in breast cancer costs, with an average rise of 172.5
million PLN (95% CI [82.4, 262.6]) and an additional 1005 breast cancer deaths
(95% CI [426, 1584]) in Poland during the simulated period. While these
results are preliminary, they highlight the potential harmful impact of
lockdowns on breast cancer treatment outcomes and costs.

###### Keywords:

Breast Cancer Costs Markov Model Covid Lockdowns.

††M. Dul and M. K. Grzeszczyk – Authors contributed equally.

##  1 Introduction

The COVID-19 pandemic impacted the lives of people around the world. To slow
down the spread of the disease, many countries introduced lockdown
restrictions in form of banning gatherings, limiting outdoor trips and
canceling public events [23]. While lockdowns positively influenced the
pandemic progression (decreased doubling time) [25] or even environment [6],
the negative impact on mental health [1], physical fitness [49], dietary
habits [8] and other important aspects of our lives are evident. In this work
we analyze the effect of pandemic lockdowns on breast cancer care in Poland.

Breast cancer is the most frequent cause of cancer deaths among women [48] and
is a high burden to public finance. There is an estimated 2.3 million women
diagnosed with breast cancer and 685,000 deaths globally in 2020 [46]. The
direct cause of breast cancer is unknown, but there exist a number of risk
factors like obesity, late menopause or alcohol use [22]. Since there are few
to no symptoms at the early stage of breast cancer, many countries introduced
screening programs in the form of free mammography procedures to support the
early detection of the disease [47]. COVID lockdowns resulted in restricted
access to healthcare [19] which consequently reduced the number of diagnosed
and treated breast cancer patients.

In this paper, we present a Markov Model-based approach to the Monte Carlo
simulation of breast cancer disease progression. The nodes of the model are
different states or cancer stages that the subject can be in at a specific
point in time. The probabilities of transitions between states are computed
based on the existing literature and empirical experiments. We use this method
to conduct 100 repetitions of seven-year-long simulations on 1% of the total
women population in Poland. In the simulation, we consider the direct costs
(medicines, surgeries, chemotherapy), indirect costs (premature death,
absenteeism, disability) of breast cancer and statistics of the number of
subjects in all states. We conduct two types of experiments. First, we perform
the simulation taking into consideration the impact of COVID lockdowns on the
accessibility of public healthcare, screening programs and treatment.
Secondly, we conduct the simulation as if there was no pandemic. We
extrapolate results to the population of the entire country.

The main contributions of this paper are as follows:

  1. 1.

We present a Markov Model-based simulation of the progression of breast
cancer.

  2. 2.

We analyze the impact of COVID lockdowns on mortality and healthcare costs
using a comparison of simulations conducted on the population of Polish women
with and without the simulated effect of pandemic.

  3. 3.

We provide a publicly available code to simulate the progression of breast
cancer: https://github.com/SanoScience/BC-MM.

The rest of the paper is structured as follows. In Section 2, we describe the
existing methods for the simulation of disease progression and present our
approach to breast cancer modeling based on Markov Models in Section 3. We
show the results of simulations with and without the effects of pandemic and
discuss how COVID-19 impacted breast cancer patients and the costs of the
disease in Section 4 and conclude in Section 5.

##  2 Related work

In this section, we describe works presented in the literature related to the
investigation of the impact of the COVID-19 pandemic on healthcare, modeling
the progression of diseases and analysis of disease costs in public finance.

###  2.1 The impact of COVID-19 pandemic on healthcare

Since the beginning of the pandemic, researchers have been concerned about the
possible, negative side effects of lockdowns [41]. Paltrinieri et al. [33]
reported that 35.1% lifestyle survey participants encountered worsening of
physical activity during lockdowns. Similar concerns were presented by Tsoukos
and Bogdanis [49] who described lower body fitness, poorer agility tests
results and increased body mass index in adolescent students as the effects of
a five-month lockdown. The negative impact of lockdowns does not end on the
deterioration of physical fitness. Mental health is one of the factors that
suffered the worst during the pandemic. Adams et al. [1] discussed a
significant decrease in self-reported mental health in the United States. The
self-harm incidents due to stress related to COVID-19 were reported in India
[40]. Cases of depression, anxiety and post-traumatic stress disorders were
expected to rise in Sub-Saharan Africa [43].

During the pandemic, access to healthcare, especially related to the treatment
of other diseases was limited [19]. Many global healthcare resources were
reallocated to prevent and treat coronavirus infections [14]. The expected
results of the depletion of healthcare resources were the increase of COVID-19
and all-cause mortality [38]. Additionally, more than 28 million surgeries
were expected to be canceled in the UK due to lockdowns [15]. In most cases,
those were operations for benign diseases, however, the effect cannot be
neglected. Jiang et. al [21] described examples of the co-epidemics of
COVID-19 and other infectious diseases and potential negative effects on the
treatment of non-communicable and chronic diseases.

Concerns regarding the impact of COVID-19 on the treatment of diseases give a
justified basis for the analysis the influence of lockdowns on breast cancer
prevalence and costs. As reported by Gathani et al. [18], there was a steep
decrease in the number of referrals for suspected breast cancer (28% lower)
and breast cancer diagnosis (16% lower) in the UK in 2020. Yin et al. [50]
describe the decline in the usage number of 3 services (breast imaging, breast
surgery and genetics consultation) in the early stages of the pandemic. In
[17], the pauses in screening programs that occurred in various countries were
described, and disease modeling was mentioned as one of the possible
approaches to analyse the repercussions of COVID-19 on breast cancer. In this
paper, we analyse the impact of those radical changes in breast cancer
diagnosis and treatment on the costs of breast cancer in public finance.

###  2.2 Modelling progression of the disease

There are multiple methods for developing disease models for the purpose of
conducting simulations [24]. One of the approaches to stochastic process
modeling (like the progression of the chronic disease) and economic impact
analysis is the utilization of Markov Modelling [11]. In such a graph model,
nodes represent stages of the disease and edges the probabilities of moving
from one state to another. For instance, Liu et al. [26] presented the
Continuous-Time Hidden Markov Model for Alzheimer’s disease progression. In
[37], a multi-state semi-Markov model was used to investigate the impact of
type 2 diabetes on the co-occurrence of cardiovascular diseases.

Markov Model can be successfully utilized to conduct an analysis of breast
cancer. Momenzadeh et al. [29] used hidden Markov Model to predict the
recurrence of breast cancer, while Pobiruchin et al. [35] presented a method
for Markov Model derivation out of real-world datasets (cancer registry’s
database). Breast cancer modeling was also used to investigate the decline in
screening, delays in diagnosis and delays in treatment during COVID-19
pandemic in the USA [3]. Alagoz et al. [3] developed three models representing
each issue. The conducted simulation exposed that there is a projected excess
of breast cancer deaths by 2030 due to the pandemic. In this paper, we present
the results of the experiments conducted with Monte Carlo simulation based on
the Markov Model of breast cancer progression in Poland.

###  2.3 Costs of breast cancer care

The analysis of disease costs for public finance is a difficult task as there
are different methods that could be used and various types of costs that have
to be taken into consideration [12]. Costs in pharmacoeconomics can be divided
into four categories: direct medical costs, direct non-medical costs, indirect
costs and intangible costs [34]. Direct medical costs are the easiest to
determine. They include the costs of medicines, diagnostic tests, hospital
visits etc. Direct non-medical costs are costs mainly related to the treatment
of the patient, but not having a medical basis. Another group of costs are
indirect costs. They are mainly related to the loss of productivity associated
with the patient’s illness or death. The intangible costs are the costs
associated with pain, suffering, fatigue and anxiety associated with the
disease, as well as side effects of treatment such as nausea. They are
difficult to estimate and measure, as they are mainly related to the patient’s
feelings [39]. In this paper, we take into consideration direct and indirect
costs only.

Depending on the methodology the calculated costs may highly vary (e.g.
$US20,000 to $US100,000 of the per-patient cost) [12]. Different studies
analyse different types of costs, making it difficult to compare them. In [7],
the mathematical model of continuous tumor growth and screening strategies was
applied for the evaluation of screening policies. In [16], cost-effectiveness
studies were used to estimate the costs of breast cancer screening per year of
life saved to be $13,200-$28,000. The total cost of breast cancer in the USA
was estimated to be $3.8 billion. Blumen et al. [10] conducted a retrospective
study to compare the treatment costs by tumor stage and type of service. The
analysis was undertaken on a population selected from the commercial claims
database which facilitated the study as costs-related data was directly
available.

##  3 Methodology

Figure 1: Breast cancer Markov Model with 10 states. Blue arrows indicate
transitions between different stages of breast cancer, yellow ones show the
diagnosis of breast cancer, green arrows indicate that the patient was healed
and red represent the death event related to breast cancer.

In this section, we describe the Markov Model used for the simulation of
breast cancer progression. We define the types and values of costs used and
provide details on the parameters used in simulations. For more details on the
actual algorithmic implementation refer to the source code available at
https://github.com/SanoScience/BC-MM.

####  3.0.1 Breast cancer Markov Model

We use the Monte Carlo simulation based on the Markov Model to describe the
course of breast cancer disease. The time horizon of the analysis is divided
into equal time increments, called Markov cycles. During each cycle, the
patient can transition from one state to another. Arrows connecting two
different states indicate allowed transitions. We applied values of those
transitions using clinical information, derived from previous studies or
estimated empirically. The arrows from the state to itself indicate that the
patient may remain in that state during cycles [44]. Transition probabilities
are used to estimate the proportion of patients who will transfer from one
state to another. The probability of an event occurring at a constant rate
($r$) during time ($t$) can be expressed by the equation:

$p=1-e^{-rt}$ (1)

In Figure 1, we present the Markov Model describing the progression of breast
cancer. There are ten states in our model: healthy, four states describing a
non-diagnosed person with breast cancer at four stages [5] of the disease
($stageN_{i}$, where $i\in\\{1,2,3,4\\}$), four states for diagnosed stages
($stageD_{i}$, where $i\in\\{1,2,3,4\\}$) and deceased. Deceased is a terminal
stage of our model and if a subject reaches this state its simulation is
terminated. We follow The American Joint Committee on Cancer which defines
four stages of breast cancer [5].

####  3.0.2 Simulation

We assume that each Markov cycle is equal to one week. Breast cancer is a
disease that develops for many years, however, the longer the simulated period
is, the less reliable are the results due to assumptions made about the
future. Therefore, we set the number of cycles to 364, corresponding to 7
years (assuming that every year has 52 weeks). This period allows us to
measure the long-term effects of the COVID-19 pandemic. We set the beginning
of the simulation to January 1st 2019 so that the simulation stabilizes
(reaches equilibrium) before the 2020 year with a pandemic. We conduct two
types of simulations, one taking into account COVID-19 lockdowns, and one
which assumes that there was no effect of COVID-19 on breast cancer treatment
and progression. We assume that lockdowns in Poland were lasting from March
2020 till the beginning of March 2021. We repeat 100 times each type and
average the collected results.

In the simulation, we take into account malignant breast cancer only (C50
ICD-10 code). Stage 0 of breast cancer (which has a different ICD-10 code)
often has a 100% 5-year survival rate. Thus, we omit this stage in the
analysis as it should not have a significant impact on costs and survival. We
simulate the breast cancer progression in women as this sex accounts for most
of the cases of the disease. We conduct computation on 1% of the
representative women population in Poland with the age distribution according
to Table 1 \- after the end of simulations, we multiply results by 100. The
minimum age of simulated patients is set to 25 because below this age the
occurrence of breast cancer is rare (Table 2). We increase the age of each
simulated person every 52 cycles.

To find the number of women diagnosed with breast cancer in 2019 we compute a
linear trend line based on the available data. There were 143,911, 151,831,
158,534, 166,031, and 174,005 patients with breast cancer in 2010, 2011, 2012,
2013, and 2014 respectively [31]. According to Agency for Health Technology
Assessment and Tariff System in Poland, those numbers rose to 227,784 and
242,838 in 2016 and 2017 [2]. The projected trend line indicated that in 2018
and 2019 there were 247,013 and 263,590 women with this disease in Poland.
Taking into consideration the distribution of cancer stages among women with
the diagnosed disease in the UK [13] and the distribution of age (Table 1) we
derive the number of diagnosed women in Poland by cancer stage and by age
(Table 3). In addition, we estimate the number of undiagnosed people. We
assume that breast cancer would only be detected using mammography and follow-
up diagnostic regimen, and around 71% of patients show up at this procedure
[42]. In 2019 the number of mammography tests was 1.041 million with 19620
cases detected (2%). The number of people who fell ill in 2019 was 263,590.
This is 2% of the 71% of people who appeared on mammograms. On this basis, the
remaining people who did not show up on the mammogram and would have a
positive result can be calculated. They are 2% of the remaining 29% that
should come for a mammogram. The estimated number of people is 108,427. Using
this number and information about the percentage of patients in a specific
stage, we calculate the number of undiagnosed patients in stages II, III and
IV (Table 3). The same strategy for stage I destabilizes the model. Thus, we
set the number of undiagnosed patients in the first stage to the same value as
for those diagnosed in the first stage in 2019. We make this assumption due to
the fact that people in stage I very often do not have symptoms yet.

Table 1: The age distribution of Polish women above 25 years in 2019 [45]. Age | Number of women | Percentage
---|---|---
25-29 | 1,233,777 | 8%
30-34 | 1,436,161 | 10%
35-39 | 1,596,757 | 11%
40-44 | 1,502,164 | 10%
45-49 | 1,294,636 | 9%
50-54 | 1,142,203 | 8%
55-59 | 1,234,478 | 8%
60-64 | 1,460,924 | 10%
65-69 | 1,359,815 | 9%
70-74 | 1,013,368 | 7%
75-79 | 640,118 | 4%
80-84 | 581,529 | 4%
85+ | 583,545 | 4%
Total | 15,079,475 | 100%
Table 2: The distribution of patients diagnosed with breast cancer in 2019 in Poland [36]. Age | 0-24 | 25-29 | 30-34 | 35-39 | 40-44 | 45-49 | 50-54 | 55-59 | 60-64 | 65-69 | 70-74 | 75-79 | 80-84 | 85+
---|---|---|---|---|---|---|---|---|---|---|---|---|---|---
# | 8 | 79 | 292 | 697 | 1,252 | 1,591 | 1,832 | 2,120 | 2,970 | 3,377 | 2,039 | 1,442 | 1,107 | 814
% | 0.0 | 0.4 | 1.5 | 3.6 | 6.4 | 8.1 | 9.3 | 10.8 | 15.1 | 17.2 | 10.4 | 7.3 | 5.6 | 4.1
Table 3: The projected distribution of diagnosed and undiagnosed women with breast cancer in 2019 in Poland. | Diagnosed | Undiagnosed
---|---|---
Age \Stage | I | II | III | IV | I | II | III | IV
25-29 | 461 | 441 | 102 | 57 | 461 | 181 | 42 | 24
30-34 | 1,705 | 1,630 | 377 | 211 | 1,705 | 670 | 155 | 87
35-39 | 4,069 | 3,891 | 900 | 504 | 4,069 | 1,600 | 370 | 208
40-44 | 7,308 | 6,989 | 1,617 | 906 | 7,308 | 2,875 | 665 | 373
45-49 | 9,287 | 8,881 | 2,055 | 1,151 | 9,287 | 3,653 | 845 | 474
50-54 | 10,694 | 10,227 | 2,366 | 1,326 | 10,694 | 4,207 | 973 | 545
55-59 | 12,375 | 11,834 | 2,738 | 1,534 | 12,375 | 4,868 | 1,126 | 631
60-64 | 17,337 | 16,579 | 3,836 | 2,150 | 17,337 | 6,820 | 1,578 | 884
65-69 | 19,713 | 18,851 | 4,361 | 2,444 | 19,713 | 7,754 | 1,794 | 1,005
70-74 | 11,902 | 11,382 | 2,633 | 1,476 | 11,902 | 4,682 | 1,083 | 607
75-79 | 8,417 | 8,050 | 1,862 | 1,044 | 8,417 | 3,311 | 766 | 429
80-84 | 6,462 | 6,179 | 1,430 | 801 | 6,462 | 2,542 | 588 | 330
85+ | 4,752 | 4,544 | 1,051 | 589 | 4,752 | 1,869 | 432 | 242

####  3.0.3 State transition probabilities

We derive the following state transition probabilities:

  1. 1.

$P(healthy\rightarrow stageN_{1})$ \- the probability of developing breast
cancer,

  2. 2.

$P(healthy\rightarrow deceased)$ \- the probability of non-cancer related
death,

  3. 3.

$P(stageN_{i}\rightarrow stageN_{i+1})$ \- the probability of cancer stage
increase,

  4. 4.

$P(stageN_{i}\rightarrow stageD_{i})$ \- the probability of cancer diagnosis,

  5. 5.

$P(stageN_{i}\rightarrow deceased)$ \- the probability of breast cancer death,

  6. 6.

$P(stageD_{i}\rightarrow stageD_{i+1})$ \- the probability of cancer stage
increase,

  7. 7.

$P(stageD_{i}\rightarrow healthy)$ \- the probability of healing,

  8. 8.

$P(stageD_{i}\rightarrow deceased)$ \- the probability of breast cancer death.

To simulate the effects of covid lockdowns we modify three transition
probabilities. The probability of cancer diagnosis is decreased because of
lockdowns and restricted access to healthcare. The probability of breast
cancer-related death is increased due to a lack of proper healthcare
assistance, and the probability of healing is decreased due to poorer
treatment during the COVID-19 pandemic. Numerically we implement the models as
follows.

We assume probability of developing breast cancer is only dependent on women’s
age. We set this probability in the following manner - age (probability;
probability in one cycle): 20 (0.1%; 0.0002%), 30 (0.5%; 0.001%), 40 (1.5%;
0.0029%), 50 (2.4%; 0.0046%), 60 (3.5%; 0.0067%), 70 (4.1%; 0.0079%), 80
(3.0%; 0.0058%) [4]. We define the probability of non-cancer related death
according to the life tables from 2019 [45]. There are multiple resources
defining the Progression-Free Survival parameter which is a time during which
the patient lives with the disease but their state is not worsening. For
example, Haba-Rodriguez et al. [20] state that the median of PFS varies
between 4 to 18 months. Thus, we empirically set the probability of cancer
stage increase for diagnosed patient to $p=k(1-e^{-\lambda t})$ where $k$ is
0.0009, $t$ is the number of cycles and $\lambda$ is 10, 15, 20 or 25
depending on the stage of the disease. It is difficult and highly uncertain to
assess the progression of breast cancer in an undiagnosed patient as no data
exist that describe those transitions. Therefore, we define the probability of
cancer stage increase for undiagnosed women in the same manner as in the case
of diagnosed cases and set $\lambda$ to 20, 25, 30 or 35 depending on the
stage which is a reasonable approximation.

We define the probability of healing based on the 5-year Disease Free Survival
(DFS) parameter which changes depending on the cancer stage [32]. The 5-year
DFS in 2019 was 0.987 (stage I), 0.873 (stage II), 0.52 (stage III), and 0.037
(stage IV). We decrease those values during lockdowns by 19% (due to a 19%
decrease in hospitalizations [30]) to 0.801 (stage I), 0.708 (stage II), 0.422
(stage III), and 0.03 (stage IV) and set the probability of healing to
$p=k(1-e^{-\lambda t})$ where $k$ is set to $\frac{1}{3}$ and $\lambda$ is
computed based on the 5-year DFS. The probability of death for diagnosed
patient is computed from the 5-year survival rate which indirectly provides
the probability of death within 5 years. Taking into consideration the 5-year
survival in stages I, II, III and IV of 0.975, 0.856, 0.44, 0.23 [32], we
compute $\lambda$ parameter of the probability of death in cycle $\leq t$
($p(x\leq t)=1-e^{-\lambda t}$) to be 0.0061, 0.0302, 0.1642 and 0.2939
respectively. For covid simulation according to predictions that the mortality
rate might increase by 9.6% [27] the $\lambda$ parameter is set to 0.0056,
0.0344, 0.1903, 0.3715 for every stage. The 3-month delay in cancer diagnosis
may increase the chances of premature death by 12% [9]. We, therefore, find
the probability of death for undiagnosed patient by increasing the 5-year
death probability for diagnosed patients and compute $\lambda$ for those
probabilities equal to 0.0061, 0.0349, 0.1989, 0.3932 depending on the stage.

In 2019, approximately 7% of all women aged 25 and over had a mammogram. The
situation changed in 2020 when the number of mammograms decreased by over
305,000, which was a decrease of about 29%. We assume that malignant breast
cancer can only be detected by mammography and diagnostic follow-ups. The
newly detected cases in 2019 (19,620) accounted for 2% of all mammograms
performed. In 2020, the number of detected cases decreased due to the COVID-19
pandemic. The average annual growth rate of new cases of breast cancer is 3%.
This means that around 20,209 new cases of cancer should have been diagnosed
in 2020. We assume that the percentage of positive mammograms did not change,
which will bring the number of detected cases to about 13,873. This is a
difference of about 6,000 cases compared to what should have been detected.
About 12.5% of mammograms are thought to have a false-negative result and data
shows that only 71% of all women show up for the examination. In 2020, the
number of these women probably decreased even more (by 29%). Therefore, we
define the probability of diagnosis in one year as
$p=\frac{l_{pos}}{l_{pos}+l_{fneg}+l_{nm}}$ where $l_{pos}$ is the number of
positive mammography cases, $l_{fneg}$ is the number of false negative
mammography and $l_{nm}$ is the number of women that should have taken part in
the mammography and would have had positive results. Thus, this probability
for 2019 is 12.4% and 11.6% during lockdowns.

####  3.0.4 Costs of breast cancer in Poland

We collect two types of costs during simulation: direct costs and indirect
costs. We divide the latter into indirect costs related to premature death and
other indirect costs related to absenteeism, presenteeism, absenteeism of
informal caregivers, presenteeism of informal carers, disability etc.

We derive direct per-person, per-stage costs in the following way. We estimate
the total direct costs for 2019 based on the estimated number of breast cancer
patients and direct costs in 2010-2014 years [31] to be 846,653 thousand PLN.
We follow the distribution of the stage-specific costs in [28]. We compute
that direct per-stage yearly costs, based on the number of patients in every
stage, in the simulation are: stage I (25% of total costs, 207,560 thousand
PLN, 1881 PLN per person), stage II (38%, 325,108 thousand PLN, 3,185 PLN per
person), stage III (25%, 215,076 thousand PLN, 5,573 PLN per person), stage IV
(12%, 98,909 thousand PLN, 7,869 PLN per person). We add those costs (divided
by the number of cycles in one year) for every diagnosed woman in each cycle
of the simulation.

We compute indirect death costs by averaging the per-person death costs
related to breast cancer patients in 2010-2014 years [31]. The average value
of premature death cost is 123,564 PLN. We add this value every time a
simulated person enters deceased state from one of the breast cancer stages.
We estimate other indirect costs in the same way. The average value of
indirect per-patient cost in 2010-2014 years is 13,159 PLN. We add this value
(divided by the number of cycles in year) for every patient in the
$stageD_{i}$ state in every cycle.

####  3.0.5 Experimental setup

We develop the Monte Carlo simulation with Python 3.10 programming language.
We conduct all simulations on a 1.80 GHz Intel Core i7 CPU and 16 GB RAM. The
execution of all simulations took 3 hours to complete.

##  4 Results and discussion

####  4.0.1 Costs of breast cancer

In Table 4, we present changes in average direct and indirect costs over the
7-years period. In all cases, the total costs incurred in the absence of
lockdowns are smaller than the ones that resulted from the simulation with the
COVID-19 pandemic. However, the only statistically significant difference
(p-value < 0.001) is in the case of indirect costs related to premature death.
This is reasonable because breast cancer is a long-lasting disease and the
costs of treatment or patient care did not change drastically due to
lockdowns. On the other hand, delayed diagnoses and surgeries resulted in more
premature deaths. The impact of the pandemic is also reflected in the total
costs of breast cancer (Table 5). the pandemic resulted in an increase in
breast cancer costs of 172.5 million PLN (average total costs with covid -
average total costs without covid) with 95% confidence interval (CI) of [82.4,
262.6]. The difference between total costs with and without lockdowns is
statistically significant (p-value < 0.001). The positive influence of
lockdowns on the progression of the pandemic should not be neglected. However,
the presented results suggest that lockdowns had a negative impact on overall
disease treatment, both socially and economically.

Table 4: Direct and indirect simulated costs (in million PLN) of breast cancer in Poland with and without COVID-19 pandemic. | DIRECT | INDIRECT_DEATH | INDIRECT_OTHER
---|---|---|---
year | NO COVID | COVID | NO COVID | COVID | NO COVID | COVID
2019 | $763\pm 4$ | $764\pm 4$ | $625\pm 79$ | $625\pm 93$ | $3297\pm 14$ | $3302\pm 13$
2020 | $770\pm 7$ | $771\pm 7$ | $827\pm 93$ | $875\pm 106$ | $3317\pm 26$ | $3321\pm 28$
2021 | $775\pm 10$ | $775\pm 9$ | $990\pm 107$ | $996\pm 113$ | $3337\pm 36$ | $3338\pm 34$
2022 | $777\pm 11$ | $777\pm 12$ | $1087\pm 118$ | $1119\pm 120$ | $3352\pm 40$ | $3355\pm 44$
2023 | $777\pm 11$ | $779\pm 13$ | $1186\pm 116$ | $1197\pm 121$ | $3367\pm 43$ | $3372\pm 50$
2024 | $775\pm 12$ | $779\pm 14$ | $1275\pm 120$ | $1270\pm 119$ | $3377\pm 49$ | $3387\pm 51$
2025 | $774\pm 13$ | $777\pm 15$ | $1306\pm 150$ | $1340\pm 30$ | $3389\pm 53$ | $3399\pm 56$
TOTAL | $5411\pm 55$ | $5422\pm 64$ | $7296\pm 246$ | $7420\pm 267$ | $23437\pm 213$ | $23474\pm 236$
Table 5: Total simulated costs (in million PLN) of breast cancer in Poland. | TOTAL
---|---
year | NO COVID | COVID
2019 | $4686\pm 77$ | $4691\pm 91$
2020 | $4915\pm 91$ | $4967\pm 107$
2021 | $5101\pm 116$ | $5109\pm 118$
2022 | $5216\pm 133$ | $5251\pm 122$
2023 | $5329\pm 122$ | $5348\pm 135$
2024 | $5428\pm 130$ | $5435\pm 133$
2025 | $5469\pm 170$ | $5516\pm 140$
TOTAL | $36144\pm 316$ | $36317\pm 333$

####  4.0.2 Breast cancer with and without COVID-19 pandemic

Simulations also showed that there was a significant difference in the number
of women deaths due to COVID-19. On average, during 7 years, 60,052 women died
taking into consideration lockdowns. This number would be smaller by 1005
deaths if the pandemic did not occur (95% CI [426, 1584]). Year-by-year
visualization of deaths is presented in Figure 2. It can be noticed that
delayed diagnoses and poorer treatment resulted in an overall increase in the
number of deaths. The long-term effects will be visible in years to come.
Figure 2 depicts also the average number of cases of diagnosed breast cancer
in Poland. There is a sharp decline in the number of diagnoses between covid
and no covid simulations in 2020. The delayed diagnoses resulted in an
increased probability of complications and death. In the following years, the
trend was reversed and more of the delayed cases were diagnosed in covid
simulation. However, the inefficient healthcare system is not capable of
making up for the lost diagnostic time.

Figure 2: The average number of breast cancer-related deaths (top) and the
average number of breast cancer diagnoses (bottom) with and without lockdowns.

####  4.0.3 Limitations

Our study is subject to limitations. First, the model and simulations
presented are only approximations of the real-world situation, and therefore,
the results should be interpreted with caution. Second, the impact of the
COVID-19 pandemic on the costs associated with breast cancer is complex and
rapidly evolving, and our study provides only a snapshot of the situation at
the time of the analysis. Third, in order to build the model, we had to make
several assumptions and rely on estimates or information from countries other
than Poland due to a lack of national data. Access to healthcare and treatment
may vary across different countries, and this may have resulted in
overestimated or underestimated data in our model. Therefore, our findings
should be interpreted in the context of these limitations and further research
is needed to validate and expand our results. To account for the uncertainty
around the course of the tumor, empirical fitting of transition probabilities
was necessary. This is because, upon diagnosis, patients are immediately
referred for treatment, leaving no research data on the disease’s development.
Furthermore, the study assumes that people with cancer did not directly die
from coronavirus infection, but those at an advanced stage of the disease may
have had a higher risk of succumbing faster after being infected with the
pathogen. It is also worth noting that the model does not consider potential
changes in healthcare policies or treatment protocols during the pandemic,
which could have affected breast cancer care costs and patient outcomes.
Despite these limitations, the study provides valuable insights into the
potential impact of the pandemic on breast cancer care costs, and its findings
could be beneficial to healthcare policymakers, clinicians, and researchers.
Nevertheless, more research is necessary to confirm and expand the results
presented in this study.

##  5 Conclusion

In this study, we have used a Monte Carlo simulation approach and a Markov
Model to analyze the effects of COVID-19 lockdowns on the costs and mortality
of breast cancer in Poland. Our findings indicate a significant negative
impact on breast cancer treatment, resulting in increased costs and higher
mortality rates. Although these findings are preliminary, they offer important
insights for future discussions on strategies that could be employed during
future pandemics. As part of our ongoing research, we plan to conduct a
sensitivity analysis of model parameters and expand our analysis to estimate
the impacts of lockdowns on other diseases.

## Acknowledgements

This publication is partly supported by the European Union’s Horizon 2020
research and innovation programme under grant agreement Sano No. 857533 and
the International Research Agendas programme of the Foundation for Polish
Science, co-financed by the European Union under the European Regional
Development Fund. We would like to thank Tadeusz Satława, Katarzyna Tabor,
Karolina Tkaczuk, and Maja Więckiewicz from Sano for help and their initial
contribution in the development of the model.

## References

  * [1] Adams-Prassl, A., Boneva, T., Golin, M., Rauh, C.: The impact of the coronavirus lockdown on mental health: evidence from the united states. Economic Policy 37(109), 139–155 (2022)
  * [2] Agencja Oceny Technologii Medycznych i Taryfikacji: Kompleksowa diagnostyka i leczenie nowotworów piersi (breast unit). https://bipold.aotm.gov.pl/assets/files/zlecenia_mz/2017/033/RPT/WS.434.3.2017_bcu.pdf (2018), [Online; accessed 25.01.2023]
  * [3] Alagoz, O., Lowry, K.P., Kurian, A.W., Mandelblatt, J.S., Ergun, M.A., Huang, H., Lee, S.J., Schechter, C.B., Tosteson, A.N., Miglioretti, D.L., et al.: Impact of the covid-19 pandemic on breast cancer mortality in the us: estimates from collaborative simulation modeling. JNCI: Journal of the National Cancer Institute 113(11), 1484–1494 (2021)
  * [4] American Cancer Society: Breast Cancer Facts & Figures (2019)
  * [5] Amin, M.B., Edge, S.B., Greene, F.L., Byrd, D.R., Brookland, R.K., Washington, M.K., Gershenwald, J.E., Compton, C.C., Hess, K.R., Sullivan, D.C., et al.: AJCC cancer staging manual, vol. 1024. Springer (2017)
  * [6] Arora, S., Bhaukhandi, K.D., Mishra, P.K.: Coronavirus lockdown helped the environment to bounce back. Science of the Total Environment 742, 140573 (2020)
  * [7] Baker, R.D.: Use of a mathematical model to evaluate breast cancer screening policy. Health Care Management Science 1(2), 103–113 (1998)
  * [8] Bennett, G., Young, E., Butler, I., Coe, S.: The impact of lockdown during the covid-19 outbreak on dietary habits in various population groups: a scoping review. Frontiers in nutrition 8, 626432 (2021)
  * [9] Bish, A., Ramirez, A., Burgess, C., Hunter, M.: Understanding why women delay in seeking help for breast cancer symptoms. Journal of psychosomatic research 58(4), 321–326 (2005)
  * [10] Blumen, H., Fitch, K., Polkus, V.: Comparison of treatment costs for breast cancer, by tumor stage and type of service. American health & drug benefits 9(1), 23 (2016)
  * [11] Briggs, A., Sculpher, M.: An introduction to markov modelling for economic evaluation. Pharmacoeconomics 13(4), 397–409 (1998)
  * [12] Campbell, J.D., Ramsey, S.D.: The costs of treating breast cancer in the us. Pharmacoeconomics 27(3), 199–209 (2009)
  * [13] Cancer Research UK: Breast cancer survival statistics. https://www.cancerresearchuk.org/health-professional/cancer-statistics/statistics-by-cancer-type/breast-cancer/survival#heading-Zero, [Online; accessed 25.01.2023]
  * [14] Chudasama, Y.V., Gillies, C.L., Zaccardi, F., Coles, B., Davies, M.J., Seidu, S., Khunti, K.: Impact of covid-19 on routine care for chronic diseases: a global survey of views from healthcare professionals. Diabetes & Metabolic Syndrome: Clinical Research & Reviews 14(5), 965–967 (2020)
  * [15] Collaborative, C.: Elective surgery cancellations due to the COVID-19 pandemic: global predictive modelling to inform surgical recovery plans. British Journal of Surgery 107(11), 1440–1449 (06 2020). https://doi.org/10.1002/bjs.11746, https://doi.org/10.1002/bjs.11746
  * [16] Elixhauser, A.: Costs of breast cancer and the cost-effectiveness of breast cancer screening. International journal of technology assessment in health care 7(4), 604–615 (1991)
  * [17] Figueroa, J.D., et al.: The impact of the covid-19 pandemic on breast cancer early detection and screening. Preventive Medicine 151, 106585 (2021)
  * [18] Gathani, T., Clayton, G., MacInnes, E., Horgan, K.: The covid-19 pandemic and impact on breast cancer diagnoses: what happened in england in the first half of 2020. British journal of cancer 124(4), 710–712 (2021)
  * [19] Goyal, D.K., et al.: Restricted access to the nhs during the covid-19 pandemic: Is it time to move away from the rationed clinical response? The Lancet Regional Health–Europe 8 (2021)
  * [20] de la Haba-Rodríguez, J., et al.: Time-to-progression in breast cancer: A stratification model for clinical trials. The Breast 17(3), 239–244 (2008)
  * [21] Jiang, P., Klemeš, J.J., Fan, Y.V., Fu, X., Bee, Y.M.: More is not enough: A deeper understanding of the covid-19 impacts on healthcare, energy and environment is crucial. International Journal of Environmental Research and Public Health 18(2), 684 (2021)
  * [22] Key, T.J., Verkasalo, P.K., Banks, E.: Epidemiology of breast cancer. The lancet oncology 2(3), 133–140 (2001)
  * [23] Koh, D.: Covid-19 lockdowns throughout the world. Occupational Medicine 70(5), 322–322 (2020)
  * [24] Kopec, J.A., et al.: Validation of population-based disease simulation models: a review of concepts and methods. BMC public health 10(1), 1–13 (2010)
  * [25] Lau, H., et al.: The positive impact of lockdown in wuhan on containing the covid-19 outbreak in china. Journal of travel medicine (2020)
  * [26] Liu, Y.Y., Li, S., Li, F., Song, L., Rehg, J.M.: Efficient learning of continuous-time hidden markov models for disease progression. Advances in neural information processing systems 28 (2015)
  * [27] Maringe, C., et al.: The impact of the covid-19 pandemic on cancer deaths due to delays in diagnosis in england, uk: a national, population-based, modelling study. The lancet oncology 21(8), 1023–1034 (2020)
  * [28] Mittmann, N., et al.: Health system costs for stage-specific breast cancer: a population-based approach. Current Oncology 21(6), 281–293 (2014)
  * [29] Momenzadeh, M., Sehhati, M., Rabbani, H.: Using hidden markov model to predict recurrence of breast cancer based on sequential patterns in gene expression profiles. Journal of Biomedical Informatics 111, 103570 (2020)
  * [30] Narodowy Fundusz Zdrowia: https://statystyki.nfz.gov.pl/Benefits/1a, [Online; accessed 15.08.2022]
  * [31] Nojszewska, E., Bodnar, L., Łyszczarz, B., Sznurkowski, J., Śliwczyński, A.: Ocena strat ekonomicznych i kosztów leczenia nowotworów piersi, szyjki macicy i jajnika w Polsce. Fundacja Instytut Innowacyjna Gospodarka (2016)
  * [32] Nowikiewicz, T., et al.: Overall survival and disease-free survival in breast cancer patients treated at the oncology centre in bydgoszcz–analysis of more than six years of follow-up. Contemporary Oncology/Współczesna Onkologia 19(4), 284–289 (2015)
  * [33] Paltrinieri, S., et al.: Beyond lockdown: the potential side effects of the sars-cov-2 pandemic on public health. Nutrients 13(5), 1600 (2021)
  * [34] Petryszyn, P.W., Kempiński, R., Michałowicz, J., Poniewierka, E.: Non-medical costs of colonoscopy. Gastroenterology Review/Przegląd Gastroenterologiczny 9(5), 270–274 (2014)
  * [35] Pobiruchin, M., Bochum, S., Martens, U.M., Kieser, M., Schramm, W.: A method for using real world data in breast cancer modeling. Journal of Biomedical Informatics 60, 385–394 (2016)
  * [36] Polish National Cancer Registry: https://onkologia.org.pl/pl/raporty, [Online; accessed 15.08.2022]
  * [37] Ramezankhani, A., Azizi, F., Hadaegh, F., Momenan, A.A.: Diabetes and number of years of life lost with and without cardiovascular disease: a multi-state homogeneous semi-markov model. Acta Diabetologica 55(3), 253–262 (2018)
  * [38] Randolph, H.E., Barreiro, L.B.: Herd immunity: understanding covid-19. Immunity 52(5), 737–741 (2020)
  * [39] Rascati, K.: Essentials of pharmacoeconomics. Lippincott Williams & Wilkins (2013)
  * [40] Sahoo, S., Rani, S., Parveen, S., Singh, A.P., Mehra, A., Chakrabarti, S., Grover, S., Tandup, C.: Self-harm and covid-19 pandemic: An emerging concern–a report of 2 cases from india. Asian journal of psychiatry 51, 102104 (2020)
  * [41] Schippers, M.C.: For the greater good? the devastating ripple effects of the covid-19 crisis. Frontiers in Psychology 11, 577740 (2020)
  * [42] Screening & Immunisations Team, NHS Digital: Breast Screening Programme, England 2017-18 (2018)
  * [43] Semo, B.w., Frissa, S.M.: The mental health impact of the covid-19 pandemic: implications for sub-saharan africa. Psychology research and behavior management 13, 713 (2020)
  * [44] Sonnenberg, F.A., Beck, J.R.: Markov models in medical decision making: a practical guide. Medical decision making 13(4), 322–338 (1993)
  * [45] Statistics Poland: https://bdl.stat.gov.pl/bdl/start, [Online; accessed 15.08.2022]
  * [46] Sung, H., Ferlay, J., Siegel, R.L., Laversanne, M., Soerjomataram, I., Jemal, A., Bray, F.: Global cancer statistics 2020: Globocan estimates of incidence and mortality worldwide for 36 cancers in 185 countries. CA: a cancer journal for clinicians 71(3), 209–249 (2021)
  * [47] Tabar, L., et al.: Mammography service screening and mortality in breast cancer patients: 20-year follow-up before and after introduction of screening. The Lancet 361(9367), 1405–1410 (2003)
  * [48] Torre, L.A., et al.: Global cancer statistics, 2012. CA: a cancer journal for clinicians 65(2), 87–108 (2015)
  * [49] Tsoukos, A., Bogdanis, G.C.: The effects of a five-month lockdown due to covid-19 on physical fitness parameters in adolescent students: A comparison between cohorts. International Journal of Environmental Research and Public Health 19(1), 326 (2022)
  * [50] Yin, K., Singh, P., Drohan, B., Hughes, K.S.: Breast imaging, breast surgery, and cancer genetics in the age of covid-19. Cancer 126(20), 4466–4472 (2020)


# Hysteresis-Free High Mobility Graphene Encapsulated in Tungsten Disulfide

Karuppasamy Pandian Soundarapandian1 Domenico De Fazio1,2
<EMAIL_ADDRESS>Sefaattin Tongay3 Frank H. L. Koppens1
<EMAIL_ADDRESS>1ICFO-Institut de Ciencies Fotoniques, The Barcelona
Institute of Science and Technology, 08860 Castelldefels (Barcelona), Spain
2Department of Molecular Sciences and Nanosystems, Ca’ Foscari University of
Venice, 30172 Venezia (Veneto), Italy  3School for Engineering of Matter,
Transport and Energy, Arizona State University, 85287 Tempe (Arizona), USA

###### Abstract

High mobility is a crucial requirement for a large variety of electronic
device applications. The state-of-the-art for high quality graphene devices is
based on heterostructures made with graphene encapsulated in $>80\,$nm-thick
flakes of hexagonal boron nitride (hBN). Unfortunately, scaling up multilayer
hBN while precisely controlling the number of layers remains an elusive
challenge, resulting in a rough material unable to enhance the mobility of
graphene. This leads to the pursuit of alternative, scalable materials, which
can be simultaneously used as substrate and encapsulant for graphene. Tungsten
disulfide (WS2) is a transition metal dichalcogenide, which was successfully
grown in large ($\sim$mm-size) multi-layers by chemical vapour deposition.
However, the resistance vs gate voltage characteristics when gating graphene
through WS2 exhibit largely hysteretic shifts of the charge neutrality point
(CNP) in the order of $\Delta n\sim$2.6$\cdot$1011 cm-2, hindering the use of
WS2 as a reliable encapsulant. The hysteresis originates due to the charge
traps from sulfur vacancies present in WS2. In this work, we report for the
first time the use of WS2 as a substrate and the overcoming of hysteresis
issues by chemically treating WS2 with a super-acid, which passivates these
vacancies and strips the surface from contaminants. The hysteresis is
significantly reduced below the noise level by at least a factor five (to
$\Delta n<$5$\cdot$1010 cm-2) and, simultaneously, the room-temperature
mobility of WS2-encapsulated graphene is as high as $\sim$6.2$\cdot$104
cm-2V-1s-1 at a carrier density $n$ $\sim$1$\cdot 10^{12}$ cm-2. Our results
promote WS2 to a valid alternative to hBN as encapsulant for high-performance
graphene devices.

Graphene, transition metal dichalcogenides, encapsulation, hysteresis,
mobility

High mobility is a crucial requirement for multiple electronic and
optoelectronic applications such as field effect transistors[1],
modulators[2], photodetectors[3, 4] and sensors[5]. Owing to its ultra-high
($>$105 cm-2V-1s-1) room-temperature mobility[6], complemented by broadband
absorption[7], scalability[8, 9] and compatibility to the complementary metal-
oxide-semiconductor (CMOS) platform[10], graphene is quickly rising for
consideration in advanced multi-purpose technologies[11, 2]. The room-
temperature mobility in graphene is limited by the scattering of carriers with
acoustic phonons[12] and it is typically inversely proportional to the
residual charge carrier density ($n$*)[13], which arises from local strain
fluctuations in graphene[13, 14]. In this context, the choice of the perfect
substrate, dictated by the need of atomic flatness and absence of charge
traps, plays a vital role for preserving the extra-ordinary properties of
graphene[15].

Suspended graphene exhibits extremely high mobility ($\sim$2$\cdot$105
cm-2V-1s-1)[16] but the integration of such devices is impractical[17].
Several substrates such as aluminium dioxide (Al2O3)[18], aluminium nitride
(AlN)[19], fused silica[20], strontium titanate (SrTiO3)[21] and calcium
fluoride (CaF2)[22] were tested to preserve the properties of suspended
graphene. Nevertheless, hexagonal boron nitride (hBN), an atomically flat,
layered dielectric material is so far the unrivalled choice for encapsulating
graphene[23]. State-of-the-art (SOTA) high-quality graphene (Gr) devices are
encapsulated with hBN flakes and assembled through a novel stamping and
cleaning technique[24]. The highest values of mobility measured in such
hBN/Gr/hBN heterostructure at room temperature exceeds $\sim$105 cm2V-1s-1 at
a carrier density of $\sim$1012 cm-2 with a $n$* of $2\cdot$109 cm-2 [24, 25].
In such experiments atomically-flat multilayer flakes of hBN $>80$ nm-thick
were used both as a substrate and as a capping layer, in order to completely
screen charge traps and roughness from the substrate underneath, but also to
protect graphene from exposure to air-contaminants[24]. However, growing
multilayer hBN films with equal flatness to the exfoliated counterpart, and
being able to precisely control the number of layers, remains an open
challenge[26, 27]. Additionally, hBN grown by chemical vapour deposition (CVD)
tends to possess a high concentration of point defects, wrinkles and grain
boundaries, resulting in a significant increase of the overall roughness[27].
Consequently, multilayer CVD-grown hBN does not represent at the moment an
ideal candidate for scalable encapsulation of high quality graphene based
devices[28, 17]

Ultra-high room temperature mobility ($>$2.5$\cdot$105 cm2V-1s-1) at a carrier
density of ($\sim$1012 cm-2) was reported when using TMD flakes as capping
layers while still using hBN as a substrate[29]. Ref.29 claimed that this
increase in mobility with respect to fully-hBN-encapsulated Gr might be due to
a modification of the acoustic phonon bands in Gr, although the origin of this
mechanism has yet to be understood. Ref.30 performed Raman characterization of
hBN/Gr heterostructures placed on a variety of substrates. The 2D peak in the
Raman spectrum of graphene originates from a double-resonant process and it is
the most intense measurable feature[31, 32]. The full-width-at-half-maximum
(FWHM) of the 2D peak (FWHM(2D)) has been shown to be related to the amount of
nanometre-scale strain variations in the sample[14]. hBN/Gr heterostructures
with TMDs as a substrate measured in Ref.30 exhibited the smallest values of
FWHM(2D), with values similar to the ones obtained with full hBN
encapsulation[30]. Refs.33, 34, 35 reported growing uniform multilayer
transition metal dichalcogenides (TMDs) by CVD, as well as demonstrating the
possibility to integrate these materials in the CMOS back-end-of-line
(BEOL)[35]. Thus, TMDs are potentially promising substrates for Gr.

One potential disadvantage of using TMDs as a substrate, and possibly as a
gate dielectric in Gr-based devices, is their proneness to hysteretic
behavior[36, 37]. Hysteresis can be calculated by extracting the difference in
charge carrier density ($\Delta n$) at the charge neutrality point (CNP)
between the forward and reverse sweeps[38], with the CNP being the transition
between electron and hole doping in the graphene channel[39]. The origin of
the hysteresis in TMDs can be explained by the presence of trap states in the
TMD layer[40]. Residues from the fabrication, active absorption of molecules
at atomic vacancies and other defect sites are all possible forms of trap
states[41]. A high density of defect sites $\sim~{}10^{13}$ cm-2[42] not only
degrades the electrical performance of the material but also represents a
detrimental factor in the photoluminescence brightness of semiconducting TMDs
by creating channels for defect-mediated non-radiative recombination
processes[43]. Ref.43 showed that a chemical passivation of the defects in
sulfur-based monolayer TMDs, obtained by using the super-acid
bis(trifluoromethane) sulfonimide (TFSI), leads to a factor of 190-fold
improvement in the material photoluminescence intensity. Although the exact
passivating mechanism is not fully understood, the authors hypothesize a
combination of two processes: on the one hand the removal of contaminants,
driven by the protonating nature of the acid and, on the other hand, an
energy-favourable reconstruction of the sulfur vacancies promoted by
hydrogenation and the rich presence of sulfur atoms in the super-acid[43].
Many optimization steps have been performed to improve the passivation of
monolayer TMD defects/trap-states by chemical treatment[44, 45], although the
approach of Ref.43 appears to be more straightforward. Furthermore, the impact
of this chemical treatment on multilayer TMDs has yet to be tested.

In this work, we report for the first time, a hysteresis-free ($\Delta
n<10^{10}$ cm-2) electrical behaviour of a high quality heterostructure formed
by graphene encapsulated in TFSI-treated multilayer tungsten disulfide
(labeled here as Tr-WS2). Specifically, we found that the hysteresis in the
Tr-WS2/Gr/Tr-WS2 heterostructure can reach values of $\Delta n<5\cdot 10^{10}$
cm-2, which is the measured noise level, and at least a factor five lower than
$\Delta n\sim 2.6\cdot 10^{11}$ cm-2 measured for a similar device but without
the chemical treatment. We attribute the large hysteresis in this latter
untreated heterostructure mostly to the charge traps present in the WS2 used
as a gate dielectric. We confirmed the validity of the TFSI treatment on the
electrical behaviour of the Tr-WS2/Gr/Tr-WS2 heterostructure by also
performing sweep rate and sweep range measurements, which display no change in
$\Delta n$. Notably, the hysteresis is also unaffected by ageing after 120
days, suggesting the resiliency and persistence of the acid treatment in three
devices. Moreover, the room temperature mobility of the Tr-WS2/Gr/Tr-WS2
heterostructure at carrier density of $\sim$1012 cm-2 is as high as
$\sim$6.2$\cdot$104 cm-2V-1s-1 with a $n$* of $\sim 10^{11}$ cm-2. Mobility is
stable and almost flat even at higher carrier density $n\sim$3.8$\cdot$1012
cm-2. The above results suggest that graphene encapsulated in chemically Tr-
WS2 could be an interesting candidate for integrated applications[1, 2, 3, 4],
especially in a context where the heterostructure is replicated with all-CVD
grown materials.

##  I Results

Figure 1: (a) Exfoliation of WS2 on SiO2/Si, (b) chemical treatment of the
exfoliated flake using TFSI for $\sim$30 minutes and (c) dried WS2 shown on
SiO2/Si, stored for later use. Processes (a-c) were performed inside a
nitrogen glovebox. (d) A treated WS2 flake on SiO2/Si is picked up by a PC
film and it is used as the top layer of the heterostructure, then (e)
exfoliated graphene is picked up with the top WS2. (f) The WS2/Gr stack is
then dropped on top of another treated WS2 flake to form the final treated
heterostructure (Tr-WS2/Gr/Tr-WS2). (g) Schematic front view of the
measurement configuration of the Hall bar shaped Tr-heterostructure. (h)
Resistance of graphene measured as a function of the back gate voltage
V${}_{\text{GS}}$ under a forward and backward sweep, comparing data between a
treated and an untreated heterostructure.

We start our work with fabricating an heterostructure made of Tr-WS2/Gr/Tr-
WS2. The fabrication procedure starts with the micro-mechanical cleavage[46]
of $\sim$6-25 nm-thick WS2 flakes and single layer graphene on 285 nm SiO2+Si
chips. Optical microscopy[47] and Raman spectroscopy[31] were used to select
the required flakes before proceeding with further processing, as depicted in
Fig.1. Next, we prepared the acid solution by dissolving TFSI in
1,2-dichlorethane (DCE) at a concentration of 5 mg/ml. The exfoliated WS2
flakes on SiO2/Si substrates were immersed in the acid for about 30 minutes
and then blow dried using a nitrogen gun as shown in Fig.1(a-c). All these
processes were performed inside a nitrogen glovebox as suggested by Ref.45. In
parallel, we prepared a polycarbonate (PC) stamp and placed it over a
polydimethylsiloxane (PDMS) support on a glass slide, in order to be able to
dry-transfer the exfoliated flakes, as reported by Ref.23. We then started
assembling the heterostructure, by first picking up a Tr-WS2 flake from a
SiO2/Si substrate at 40°C, as shown in Fig.1(d). Next, a single layer graphene
is picked up using the stamp+Tr-WS2, again at 40°C. Finally, the whole stack
is dropped on the bottom Tr-WS2 (Fig.1(c)), followed by release of the PC film
at 120°C 1(e-f) and dissolution of PC in chloroform lasting $\sim$10 minutes.
The above heterostructure is then shaped in the form of an Hall bar (see
MethodsIII). In order to disentangle the role of the top and bottom WS2 layers
in the device hysteresis, we also prepared an Hall bar where the WS2 layers
were not treated with TFSI (WS2/Gr/WS2, untreated heterostucture).

Fig.1(g) is a schematic of our four probe measurement setup: a bias current is
driven between the source and drain external contacts, whilst the voltage
difference is read between the two inner probes. A back gate voltage is
applied through the bottom (Tr-WS2 or WS2)+SiO2 insulating layers. Fig.1(h)
shows the four probe resistance measurements performed on the treated and
untreated heterostructures by sweeping the back gate voltage
(V${}_{\text{GS}}$). In order to evaluate the hysteresis, we performed forward
(red) and reverse (blue) sweeps of the untreated (solid line) and treated
(line with circles) heterostructures. The untreated heterostructure exhibits a
considerable shift in the CNP, $\Delta$V${}_{\text{CNP}}$, of $\sim$5V). This
shift could be explained as follows: when we sweep the back gate voltage,
depending on the polarity, electrons or holes transfer from the graphene to
trap sites on the substrate and become trapped. On sweeping back the gate
voltage, the trapped charges electrostatically dope graphene, which manifests
in shifting its CNP[38]. In our heterostructure, hysteresis is caused by the
charge trapping of carriers from graphene within the WS2 vacancy sites[48, 49]
and also at the interface with the SiO2 surface[38], with the former being
orders of magnitude more dominant with respect to the latter[42, 50, 38]. The
heterostructure sample shows less than an order of magnitude hysteresis
compared with the untreated heterostructure, with a $\Delta$V${}_{\text{CNP}}$
of $\sim$0.2V. The charge trap density $\Delta n$ is then also a measure of
the sample hysteresis and can be calculated from the
$\Delta$V${}_{\text{CNP}}$ as[38]:

$\Delta n=\frac{C_{g}\cdot\Delta V_{\text{CNP}}}{2e}.$ (1)

Figure 2:  (a) Graphical representation of the as-exfoliated WS2 on SiO2/Si
with vacancies and substitutional oxygen occupancies (the red dashed line
shows an area with two substitutional O2 atoms) and (b) change of the surface
after the treatment with TFSI (here the red dashed line represents how sulfur
atoms have regained their place in the lattice instead of O2 substitutional
atoms). (c) AFM histogram comparison of the treated (red) and untreated
(black) WS2 on SiO2/Si. (d) Raman spectrum of the treated (red) and untreated
(black) WS2 on SiO2/Si.

where Cg is the gate capacitance calculated from:

$C_{g}=\frac{\epsilon_{0}\cdot\epsilon_{1}\cdot\epsilon_{2}}{\epsilon_{2}\cdot
d_{1}+\epsilon_{1}\cdot d_{2}}.$ (2)

As we gate graphene through two dielectrics, here $\epsilon_{1}$ is the
relative permittivity of SiO2 (3.9), $\epsilon_{2}$ is the relative
permittivity of WS2 (6.2, from Ref.51), $\epsilon_{0}$ is the vacuum
permittivity, while $d_{1}$ and $d_{2}$ are the thickness values of SiO2 and
WS2, respectively. We measured a trap density $\Delta n$ of $\sim 2.6\cdot
10^{11}$ cm-2 for the untreated heterostructure and $\sim 7\cdot 10^{9}$ cm-2
for the Tr-WS2/Gr/Tr-WS2 (treated heterostructure) by taking the absolute
difference between the maximum values in resistance in forward and backward
sweeps. However, we note that the resistance of the sample fluctuates,
potentially masking the correct position of the CNP. A difference of 0.7 V
(corresponding to $\Delta n\sim 5\cdot 10^{10}$ cm-2) is indeed observed
between the position of the CNP taken from the raw data as the maximum value
of resistance and the same value taken while smoothing the data, instead. We
then infer that a more conservative value of $\Delta n\sim 5\cdot 10^{10}$
cm-2 can be claimed as upper bound of hysteresis, which is certainly not
overcome by the treated heterostructure. From now on we refer to this value as
the noise level. This clearly indicates that the TFSI treatment played a
crucial role in suppressing the defect sites, drastically decreasing the
device hysteresis by almost one order of magnitude compared to the untreated
WS2/Gr/WS2 heterostructure.

In order to be able to address the exhibited behaviour of the treated
heterostructure in Fig.1(h) one has to carefully analyse the existing problem
of hysteresis in TMD-based devices and the origin of such effect. The
predominant intrinsic defect types in as-exfoliated and as-grown TMDs are
chalcogen vacancies, depicted as hollow spheres in Fig.2[52]. At ambient
conditions, these vacancies are often filled with substitutional atoms, mainly
oxygen[53, 54, 55], represented as white spheres in Fig.2. Depending on their
concentration, defects can considerably alter the electronic and optical
properties of TMD-based devices[56, 57]. Vacancies represent the main source
of trapping/detrapping of charges[48], which in turn results in hysteresis
effects[36, 37], especially when a TMD is used as a dielectric. Many methods
have been proposed to eliminate the substitutional atoms and/or “passivate”
the vacancies with the original chalcogen atom species[43, 44, 45, 58].

Figure 3: (a) Device resistivity ($\rho$) comparison of the TFSI-treated
samples as a function of gate voltage V${}_{\text{GS}}$. (b) Conductivity of
the treated samples as a function of $n$ plotted in log scale; the inset plots
$\mu_{D}$ as a function of $n$*. (c) Drude mobility ($\mu_{D}$) comparison of
the TFSI-treated samples as a function of carrier density $n$. The grey shaded
area with the black dotted borders indicates the electron-phonon limited
mobility.

In this work, we chose to perform a non invasive chemical approach as proposed
in Ref.43. As schematized by the red dashed circle in Fig.2(a) and (b), the
treatment with TFSI removes the substitutional oxygen atoms or any other
absorbed contaminant and passivates the vacancy site[43, 45]. Performing a
chemical treatment could theoretically deteriorate the morphological
properties of the TMD. For this reason we scanned by atomic force microscopy
(AFM) an area of 15 $\mu$m$\times$15 $\mu$m of the sample before and after the
treatment (2(c) and insets) and we plotted an histogram of the height
profiles. Black and red circles represent height profiles of the samples
acquired before and after the chemical treatment, respectively. On fitting
these height distributions with a Gaussian, we obtain a standard deviation
which points towards a roughness around $\sim$155 pm and $\sim$145 pm for the
sample before and after the treatment, respectively. These AFM results suggest
that the chemical treatment had produced no significant morphological change
to the material surface. Raman characterisation was also performed using a
Renishaw Invia Spectrometer with 532 nm laser excitation and a 100$\times$
objective lens. Raman spectra of the WS2 flake on SiO2/Si before and after
treatment are represented as the black and red lines in Fig.2(d),
respectively. On carefully examining the E${}^{1}_{2g}$ mode at $\sim$350 cm-1
and the A${}^{1}_{g}$ mode at $\sim$420 cm-1, we observe a slight blue-shift
after the treatment, which has been associated with a reduction of the doping
in the sample[58, 60, 61]. This effect is in agreement with the removal of the
substitutional oxygen and passivation of the vacancy sites by the TFSI
treatment explained earlier. Results from AFM and Raman demonstrate the
robustness of the TFSI treatment, and therefore open a path for the
exploration of more complicated heterostructures.

Figure 4: Systematic study of $\Delta n$ as a function of (a) sweep rate, (b)
sweep range and (c) ageing over time for Sample 1 (blue), Sample 2 (red) and
Sample 3 (black). The dotted lines in the plots are the noise levels of each
sample.

We then fabricated three different Hall bars out of Tr-WS2/Gr/Tr-WS2
heterostructures (which we name here as Sample 1, 2, 3, respectively), with
thickness of the top and bottom WS2 layers in the range $\sim$6-25 nm. The
length $L$ and width $W$ of the Hall bars are $L$=2 $\mu$m-$W$=7 $\mu$m for
Sample 1, $L$=4 $\mu$m-$W$=10 $\mu$m for Sample 2, and $L$=1.5 $\mu$m-$W$=2.5
$\mu$m for Sample 3. We focused on the transport properties of such Hall bars
measured at room temperature. Fig.3(a) is the resistivity ($\rho$) of the
treated Samples 1-3 measured while sweeping the backgate voltage
V${}_{\text{GS}}$. Here we use resistivity $\rho$ instead of resistance to
have a more direct comparison among the three samples without accounting for
geometrical differences among the Hall bars. The CNP in the three samples is
positioned within ± 2.7 V from zero gate bias, which correspond to variations
in the Fermi level of graphene (at zero gate bias) in the three samples within
± 74 meV[59], indicating a minimal doping. The $\rho$ values in Samples 1-3
are comparable, varying from $\sim$640-990 $\Omega$ at CNP to $\sim$100-120
$\Omega$ at high doping. Fig.3(b) shows the measured conductivity of graphene
at room temperature, plotted on a double log scale as a function of $n$. The
carrier density $n$ is calculated from $n$ =
$C_{g}$(V${}_{\text{GS}}-$V${}_{\text{CNP}}$)/$e$[23]. $n^{*}$ can be
extracted from this plot by noting the intersection point between the two
quasi-linear sections of the conductivity plot as highlighted by the green
dashed line. The $n^{*}$ for Sample 1, 2, and 3 are $\sim$1.7, $\sim$1.8 and
$\sim$2.2$\cdot 10^{11}$ cm-2, respectively. The Drude mobility $\mu_{D}$ is
calculated from $\mu_{D}$ = $\sigma$/($ne$), where $\sigma$ = 1/$\rho$ is the
conductivity of graphene. At technologically relevant charge carrier densities
$n\sim$ 1012 cm-2, Sample 1, 2 and 3, exhibit a $\mu_{D}$ of $\sim$6.2,
$\sim$5.0 and $\sim$2.5$\cdot$104 cm-2V-1s-1, respectively, as shown in the
inset of Fig.3(b). The slight difference in the mobility among the samples,
can be correlated to the $n^{*}$ and agrees well with the inversely linear
relationship between $\mu_{D}$ and $n^{*}$[13].

In general, at higher $n$ ($>10^{12}$ cm-2) carriers couple to acoustic
phonons and the mobility degrades[62, 63, 12]. Fig.3(c) is the plot of the
$\mu_{D}$ as a function of the charge carrier density $n$. The gray shaded
region with the black dotted borders illustrates the upper and lower limit of
the electron-phonon limited mobility model from Ref.12, which are calculated
assuming deformation potential (D) values of 18-20 eV[12]. We then decided to
probe the sample with intermediate values of mobility up to a
$n\sim$-3.8$\cdot 10^{12}$ cm-2 and we recorded a factor $\sim$30% drop in
mobility compared to the maximum acceptable value measured at $n\sim$-1$\cdot
10^{12}$ cm-2. This is a factor two lower compared to the mobility drop
recorded in hBN-encapsulated graphene devices, which typically undergo a
$\sim$70% drop in mobility in the same doping range[24, 6]. This behaviour is
one crucial advantage of TMD-encapsulated graphene devices over hBN-
encapsulated graphene devices, especially in the context of applications such
as photodetectors or modulators operating at high doping. A comparison between
the record-high mobility obtained with hBN-encapsulated graphene samples and
the treated-TMD-encapsulated graphene samples used here would not be fair, as
the reported values of the former refer to using hBN thick flakes of $\sim>$80
nm [23, 24]. In this work the thickness of both top and bottom WS2 layers in
all samples were kept below $\sim$25 nm. We strived to contain the WS2
thickness below $\sim$25 nm so that the used materials would be as similar as
possible, in terms of thickness, to those successfully grown in the research
community by CVD[33, 34, 35]. While comparing our samples with those using hBN
flakes of the same thickness, we obtain similar values of $\mu_{D}$[23].

In order to probe the reliability of the treatment and its durability in time,
we monitored the $\Delta n$ in Samples 1-3 while changing the sweep rate,
range and re-testing the samples over time. The voltage sweeping rate is very
important as reducing the sweep rate gives charges sufficient time to become
trapped in the defect site, resulting in higher hysteresis. Therefore, one
should expect an increase in hysteresis, while reducing the sweep rate if
traps are present in the Tr-WS2. However, Fig.4(a) shows that by reducing the
sweep rates from 0.1 Vs-1 to 0.005 Vs-1 no significant change in the
hysteresis values $\Delta n$ occurs: remarkably, all the samples show very low
variations of $\Delta n$ below the noise level for any of the utilised sweep
rates, confirming the high efficiency of the treatment in passivating the
defect sites and hindering charge trap mechanisms. The $\Delta n$ of Samples
1-3 are lower than the noise level measured for each of the samples, indicated
as color-matched dashed lines in Fig.4(a). Next, we monitored the $\Delta n$
with respect to the V${}_{\text{GS}}$ sweep range, as we should expect an
increase in the device hysteresis with the increase in the voltage sweep
range[38]. This is due to the fact that, by increasing the V${}_{\text{GS}}$
sweep range, more and more carriers become available for transport, but they
are also more likely to become trapped in the defect sites. In Fig.4(b) we
plot $\Delta n$ against various gate voltage ranges, where, for example, a
V${}_{\text{GS}}$ sweep from -10 V to 10 V corresponds to a voltage range of
10 V, considering the CNP being close to 0 V. The measurement range of each
sample is limited by the WS2 and SiO2 breakdown voltages, assumed to be
$\sim$0.1 MVcm-1[64] and $\sim$15 MVcm-1[65], respectively. For these
measurements the sweep rate is kept constant at 0.0125 Vs-1. Samples show
almost no dependence on the sweep range with $\Delta n$ kept consistently
below $\sim 5\cdot 10^{10}$ cm-2. However, when the V${}_{\text{GS}}$ sweep
range was increased to 20 V for Sample 1, having 6 nm bottom Tr-WS2 thickness,
we start to notice an increase in $\Delta n$ to $\sim 7\cdot 10^{10}$ cm-2.
This particular behaviour could be attributed to the onset of the dielectric
breakdown of the the bottom Tr-WS2 flake[64].

Next, the $\Delta n$ of the treated heterostructures is monitored at regular
time intervals to check the integrity of the TFSI treatment over time. For
these measurements we also adopted a sweep rate of 0.0125 Vs-1. Samples were
measured until 120 days from the fabrication date and data are shown in
Fig.4(c). All the three samples show again no dependency towards the ageing,
confirming the resilience of the chemical treatment. It is worth mentioning
that the samples were stored in ambient lab atmosphere for the entire duration
of this experiment. In any of the measurement conditions $\Delta n$ has never
overcome $\sim 5\cdot 10^{10}$ cm-2. To highlight the significance of the
hysteresis treatment when measuring a Tr-WS2/Gr/Tr-WS2 heterostructure, we
outline the noise level of the Samples 1-3 in each of the graphs of Fig.4, and
all the data-points are below such values of noise, proving that the
hysteresis is indeed minimal.

##  II Conclusion

In conclusion, we have demonstrated that TFSI-treated WS2 can be an excellent
choice for encapsulating graphene and also performs well as gate dielectric.
The TFSI treatment is highly stable, resilient and it drastically reduces the
hysteresis in high mobility WS2/Gr/WS2 heterostructures. This could be due to
the removal and replacement of substitutional atoms responsible for charge
traps, as suggested by a Raman analysis and transport measurements. The
hysteresis remains well under $\Delta n\sim 5\cdot 10^{10}$ cm-2, even when
the sample is subjected to large sweep rates, range and ageing tests. This
hysteresis value is at least five times smaller than that measured in an
untreated heterostructure. In one treated heterostructure we achieved a
mobility as high as $\sim$6.0$\cdot$104 cm-2V-1s-1. Due to the weak influence
of surface acoustic phonons, the mobility of the treated heterostructures are
almost flat up to $n\sim 3.8\cdot 10^{12}$ cm-2. These results suggest that
TMD/Gr/TMD heterostructures could be adopted in advanced optoelectronic
applications requiring low hysteresis and high mobility at high carrier
concentrations.

##  III Methods

Materials: WS2 highly crystalline van der Waals (vdW) crystals were
synthesized through the self-flux technique. In the first step, the
stoichiometric ratio of elemental puratronic grade (99.9999% pure) purity
tungsten powder and sublimation purified sulfur powder (99.9999% pure) were
mixed and sealed under 10-6 Torr inside a quartz ampoule. The powder was
heated at 950°C for 5 days to create vdW powders. In the second step, these
pre-reacted powders were resealed in a quartz ampoule with extra sulfur to
reach 1:2 stoichiometry. In the second step, all the ingredients were kept at
a high-temperature zone (1040°C) with a temperature drop of 50°C to initiate
the vapor transport and realize vdW WS2 crystals after 7 weeks. The resultant
crystals measured 2-6 mm in size. Overall, the crystal appearance and
spectroscopic signatures were superior to chemical vapor transported crystals.
Single layer graphene was exfoliated from highly-oriented pyrolytic graphite
(HOPG) commercially bought from HQ Graphene labelled as natural graphite. TFSI
(product number 15220) and 1,2-dichlorethane (DCE) (product number 284505)
were bought from Sigma-Aldrich.

Fabrication: Lithography was performed by using an Inspect F50 RAITH Elphy
Plus system and a poly(methyl methacrylate) (PMMA) mask. The heterostructure
was then dry-etched (reactive ion etching) using a trifluoromethane/oxygen
(CHF3/O2) mixture of gases with a flow rate of 40/4 sccm in a Plasmalab System
100 from Oxford Instruments. Cr was deposited thermally with thickness of 5 nm
using a Lesker LAB18 evaporator, while Au was deposited thermally with the
same instrument and with a 50 nm thickness. The evaporation was followed by an
overnight lift-off in acetone, dipping in isopropanol and nitrogen blow-
drying.

Electrical Transport Measurement: transport measurements were performed using
a Lakeshore cryogenic probe station at a pressure of $\sim 10^{-3}$ mbar with
a lock-in amplifier from Stanford Research Systems SR810. A bias of 100 nA was
applied between the source and drain electrodes with a lock-in modulation
frequency of 373 Hz. The gate voltages were applied using a Keithley 2460
sourcemeter.

##  IV Acknowledgements

The authors acknowledge Matteo Ceccanti for graphical support and Francisco
Bernal Texca for his help towards setting up the fabrication steps inside the
glovebox. F.H.L.K. acknowledges financial support from the Government of
Catalonia trough the SGR grant and from the Spanish Ministry of Economy and
Competitiveness through the “Severo Ochoa” Programme for Centres of Excellence
in R&D (CEX2019-000910-S), support by Fundacio Cellex Barcelona, Generalitat
de Catalunya, through the CERCA program, and the Mineco Grants Plan Nacional
(FIS2016-81044-P) and the Agency for Management of University and Research
Grants (AGAUR) 2017 SGR 1656. Furthermore, the research leading to these
results has received funding from the European Commission in the Horizon 2020
Framework Programme under Grant Agreements Nos. 785219 (Core2) and 881603
(Core3) of the Graphene Flagship. This work was also financially supported by
the German Science Foundation (DFG) within the priority program FFlexCom
Project “GLECS” (Contract No. NE1633/3). D.D.F. acknowledges funding from the
“Severo Ochoa” Excellence Programme. S.T acknowledges support from DOE-
SC0020653 (materials synthesis), NSF DMR-2206987 (magnetic characterization),
NSF ECCS 2052527 (electrical optimization), and NSF DMR 2111812 (optical
optimization).

## References

  * [1] S. M. Sze and K. K. Ng, Physics of Semiconductor Devices (John Wiley and Sons, Inc., 2006).
  * [2] M. Romagnoli, V. Sorianello, M. Midrio, F. H. L. Koppens, C. Huyghebaert, D. Neumaier, P. Galli, W. Templ, A. D’Errico, and A. C. Ferrari, Graphene-Based Integrated Photonics for next-Generation Datacom and Telecom, Nat. Rev. Mater. 3, 392 (2018).
  * [3] Y. Kang, H. D. Liu, M. Morse, M. J. Paniccia, M. Zadka, S. Litski, G. Sarid, A. Pauchard, Y. H. Kuo, H. W. Chen, W. S. Zaoui, J. E. Bowers, A. Beling, D. C. McIntosh, X. Zheng, and J. C. Campbell, Monolithic Germanium/Silicon Avalanche Photodiodes with 340 GHz Gain–Bandwidth Product, Nat. Photonics 3, 59 (2008).
  * [4] J. E. Muench, A. Ruocco, M. A. Giambra, V. Miseikis, D. Zhang, J. Wang, H. F. Y. Watson, G. C. Park, S. Akhavan, V. Sorianello, M. Midrio, A. Tomadin, C. Coletti, M. Romagnoli, A. C. Ferrari, and I. Goykhman, Waveguide-Integrated, Plasmonic Enhanced Graphene Photodetectors, Nano Lett. 19, 7632 (2019).
  * [5] Z. Wang, L. Banszerus, M. Otto, K. Watanabe, T. Taniguchi, C. Stampfer, and D. Neumaier, Encapsulated Graphene-Based Hall Sensors on Foil with Increased Sensitivity, Phys. Status Solidi 253, 2316 (2016).
  * [6] L. Wang, I. Meric, P. Y. Huang, Q. Gao, Y. Gao, H. Tran, T. Taniguchi, K. Watanabe, L. M. Campos, D. A. Muller, J. Guo, P. Kim, J. Hone, K. L. Shepard, and C. R. Dean, One-Dimensional Electrical Contact to a Two-Dimensional Material, Science 342, 614 (2013).
  * [7] R. R. Nair, P. Blake, A. N. Grigorenko, K. S. Novoselov, T. J. Booth, T. Stauber, N. M. R. Peres, and A. K. Geim, Fine Structure Constant Defines Visual Transparency of Graphene, Science 320, 1308 (2008).
  * [8] S. Bae, H. Kim, Y. Lee, X. Xu, J. S. Park, Y. Zheng, J. Balakrishnan, T. Lei, H. Ri Kim, Y. Il Song, Y. J. Kim, K. S. Kim, B. Özyilmaz, J. H. Ahn, B. H. Hong, and S. Iijima, Roll-to-Roll Production of 30-Inch Graphene Films for Transparent Electrodes, Nat. Nanotechnol. 5, 574 (2010).
  * [9] E. S. Polsen, D. Q. McNerny, B. Viswanath, S. W. Pattinson, and A. John Hart, High-Speed Roll-to-Roll Manufacturing of Graphene Using a Concentric Tube CVD Reactor, Sci. Rep. 5, 1 (2015).
  * [10] S. Goossens, G. Navickaite, C. Monasterio, S. Gupta, J. J. Piqueras, R. Pérez, G. Burwell, I. Nikitskiy, T. Lasanta, T. Galán, E. Puma, A. Centeno, A. Pesquera, A. Zurutuza, G. Konstantatos, and F. Koppens, Broadband Image Sensor Array Based on Graphene–CMOS Integration, Nat. Photonics 11, 366 (2017).
  * [11] D. Akinwande, C. Huyghebaert, C. H. Wang, M. I. Serna, S. Goossens, L. J. Li, H. S. P. Wong, and F. H. L. Koppens, Graphene and Two-Dimensional Materials for Silicon Technology, Nature 573, 507 (2019).
  * [12] E. H. Hwang and S. Das Sarma, Acoustic Phonon Scattering Limited Carrier Mobility in Two-Dimensional Extrinsic Graphene, Phys. Rev. B 77, 115449 (2008).
  * [13] N. J. G. Couto, D. Costanzo, S. Engels, D. K. Ki, K. Watanabe, T. Taniguchi, C. Stampfer, F. Guinea, and A. F. Morpurgo, Random Strain Fluctuations as Dominant Disorder Source for High-Quality on-Substrate Graphene Devices, Phys. Rev. X 4, 041019 (2014).
  * [14] C. Neumann, S. Reichardt, P. Venezuela, M. Drögeler, L. Banszerus, M. Schmitz, K. Watanabe, T. Taniguchi, F. Mauri, B. Beschoten, S. V. Rotkin, and C. Stampfer, Raman Spectroscopy as Probe of Nanometre-Scale Strain Variations in Graphene, Nat. Commun. 6, 1 (2015).
  * [15] C. R. Dean, A. F. Young, I. Meric, C. Lee, L. Wang, S. Sorgenfrei, K. Watanabe, T. Taniguchi, P. Kim, K. L. Shepard, and J. Hone, Boron Nitride Substrates for High-Quality Graphene Electronics, Nat. Nanotechnol. 5, 722 (2010).
  * [16] K. I. Bolotin, K. J. Sikes, Z. Jiang, M. Klima, G. Fudenberg, J. Hone, P. Kim, and H. L. Stormer, Ultrahigh Electron Mobility in Suspended Graphene, Solid State Commun. 146, 351 (2008).
  * [17] M. A. Giambra, V. Mišeikis, S. Pezzini, S. Marconi, A. Montanaro, F. Fabbri, V. Sorianello, A. C. Ferrari, C. Coletti, and M. Romagnoli, Wafer-Scale Integration of Graphene-Based Photonic Devices, ACS Nano 15, 3171 (2021).
  * [18] S. Kim, J. Nah, I. Jo, D. Shahrjerdi, L. Colombo, Z. Yao, E. Tutuc, and S. K. Banerjee, Realization of a High Mobility Dual-Gated Graphene Field-Effect Transistor with Al2O3 Dielectric, Appl. Phys. Lett. 94, 062107 (2009).
  * [19] J. Gun Oh, S. Ki Hong, C. K. Kim, J. Hoon Bong, J. Shin, S. Y. Choi, and B. Jin Cho, High Performance Graphene Field Effect Transistors on an Aluminum Nitride Substrate with High Surface Phonon Energy, Appl. Phys. Lett. 104, 193112 (2014).
  * [20] M. E. Ramón, K. N. Parrish, S. F. Chowdhury, C. W. Magnuson, H. C. P. Movva, R. S. Ruoff, S. K. Banerjee and D. Akinwande, IEEE Trans. Nanotechnol 11, 877 (2012).
  * [21] N. J. G. Couto, B. Sacépé, and A. F. Morpurgo, Transport through Graphene on SrTiO3, Phys. Rev. Lett. 107, 225501 (2011).
  * [22] P. Klar, E. Lidorikis, A. Eckmann, I. A. Verzhbitskiy, A. C. Ferrari, and C. Casiraghi, Raman Scattering Efficiency of Graphene, Phys. Rev. B 87, 205435 (2013).
  * [23] D. G. Purdie, N. M. Pugno, T. Taniguchi, K. Watanabe, A. C. Ferrari, and A. Lombardo, Cleaning Interfaces in Layered Materials Heterostructures, Nat. Commun. 9, 1 (2018).
  * [24] Z. Huang, A. Alharbi, W. Mayer, E. Cuniberto, T. Taniguchi, K. Watanabe, J. Shabani, and D. Shahrjerdi, Versatile Construction of van Der Waals Heterostructures Using a Dual-Function Polymeric Film, Nat. Commun. 11, 1 (2020).
  * [25] S. Pezzini, V. Mišeikis, S. Pace, F. Rossella, K. Watanabe, T. Taniguchi, and C. Coletti, High-Quality Electrical Transport Using Scalable CVD Graphene, 2D Mater. 7, 041003 (2020).
  * [26] S. M. Kim, A. Hsu, M. H. Park, S. H. Chae, S. J. Yun, J. S. Lee, D. H. Cho, W. Fang, C. Lee, T. Palacios, M. Dresselhaus, K. K. Kim, Y. H. Lee, and J. Kong, Synthesis of Large-Area Multilayer Hexagonal Boron Nitride for High Material Performance, Nat. Commun. 6, 1 (2015).
  * [27] Y. Shen, W. Zheng, K. Zhu, Y. Xiao, C. Wen, Y. Liu, X. Jing, M. Lanza, Y. Shen, W. Zheng, Y. Xiao, C. Wen, Y. Liu, X. Jing, K. Zhu, and M. Lanza, Variability and Yield in H-BN-Based Memristive Circuits: The Role of Each Type of Defect, Adv. Mater. 33, 2103656 (2021).
  * [28] V. Shautsova, A. M. Gilbertson, N. C. G. Black, S. A. Maier, and L. F. Cohen, Hexagonal Boron Nitride Assisted Transfer and Encapsulation of Large Area CVD Graphene, Sci. Rep. 6, 1 (2016).
  * [29] L. Banszerus, T. Sohier, A. Epping, F. Winkler, F. Libisch, F. Haupt, K. Watanabe, T. Taniguchi, K. Müller-Caspary, N. Marzari, F. Mauri, B. Beschoten, and C. Stampfer, Extraordinary High Room-Temperature Carrier Mobility in Graphene-WSe2 Heterostructures, arXiv:1909.09523 (2019).
  * [30] L. Banszerus, H. Janssen, M. Otto, A. Epping, T. Taniguchi, K. Watanabe, B. Beschoten, D. Neumaier, and C. Stampfer, Identifying Suitable Substrates for High-Quality Graphene-Based Heterostructures, 2D Mater. 4, 025030 (2017).
  * [31] A. C. Ferrari, J. C. Meyer, V. Scardaci, C. Casiraghi, M. Lazzeri, F. Mauri, S. Piscanec, D. Jiang, K. S. Novoselov, S. Roth, and A. K. Geim, Raman Spectrum of Graphene and Graphene Layers, Phys. Rev. Lett. 97, 187401 (2006).
  * [32] A. C. Ferrari and D. M. Basko, Raman Spectroscopy as a Versatile Tool for Studying the Properties of Graphene, Nat. Nanotechnol. 8, 235 (2013).
  * [33] C. Jung, S. M. Kim, H. Moon, G. Han, J. Kwon, Y. K. Hong, I. Omkaram, Y. Yoon, S. Kim, and J. Park, Highly Crystalline CVD-Grown Multilayer MoSe2 Thin Film Transistor for Fast Photodetector, Sci. Rep. 5, 1 (2015).
  * [34] P. D. Cunningham, K. M. McCreary, A. T. Hanbicki, M. Currie, B. T. Jonker, and L. M. Hayden, Charge Trapping and Exciton Dynamics in Large-Area CVD Grown MoS2, J. Phys. Chem. C 120, 5819 (2016).
  * [35] T. Schram, Q. Smets, B. Groven, M. H. Heyne, E. Kunnen, A. Thiam, K. Devriendt, A. Delabie, D. Lin, M. Lux, D. Chiappe, I. Asselberghs, S. Brus, C. Huyghebaert, S. Sayan, A. Juncker, M. Caymax, and I. P. Radu, WS2 Transistors on 300 Mm Wafers with BEOL Compatibility, Eur. Solid-State Device Res. Conf. 212 (2017).
  * [36] H. Wang, D. Ren, C. Lu, and X. Yan, Investigation of Multilayer WS2 Flakes as Charge Trapping Stack Layers in Non-Volatile Memories, Appl. Phys. Lett. 112, 231903 (2018).
  * [37] H. An, Y. H. Lee, J. H. Lee, C. Wu, B. M. Koo, and T. W. Kim, Highly Stable and Flexible Memristive Devices Based on Polyvinylpyrrolidone: WS2 Quantum Dots, Sci. Rep. 10, 1 (2020).
  * [38] H. Wang, Y. Wu, C. Cong, J. Shang, and T. Yu, Hysteresis of Electronic Transport in Graphene Transistors, ACS Nano 4, 7221 (2010).
  * [39] K. S. Novoselov, A. K. Geim, S. V. Morozov, D. Jiang, Y. Zhang, S. V. Dubonos, I. V. Grigorieva, and A. A. Firsov, Electric Field in Atomically Thin Carbon Films, Science 306, 666 (2004).
  * [40] J. Roh, J. H. Lee, S. H. Jin, and C. Lee, Negligible Hysteresis of Molybdenum Disulfide Field-Effect Transistors through Thermal Annealing, J. Inf. Disp. 17, 103 (2016).
  * [41] C. Lan, X. Kang, Y. Meng, R. Wei, X. Bu, S. Yip, and J. C. Ho, The Origin of Gate Bias Stress Instability and Hysteresis in Monolayer WS2 Transistors, Nano Res. 13, 3278 (2020).
  * [42] H. Qiu, T. Xu, Z. Wang, W. Ren, H. Nan, Z. Ni, Q. Chen, S. Yuan, F. Miao, F. Song, G. Long, Y. Shi, L. Sun, J. Wang, and X. Wang, Hopping Transport through Defect-Induced Localized States in Molybdenum Disulphide, Nat. Commun. 4, 1 (2013).
  * [43] M. Amani, D. H. Lien, D. Kiriya, J. Xiao, A. Azcatl, J. Noh, S. R. Madhvapathy, R. Addou, K. C. Santosh, M. Dubey, K. Cho, R. M. Wallace, S. C. Lee, J. H. He, J. W. Ager, X. Zhang, E. Yablonovitch, and A. Javey, Near-Unity Photoluminescence Quantum Yield in MoS2, Science 350, 1065 (2015).
  * [44] H. Kim, D. H. Lien, M. Amani, J. W. Ager, and A. Javey, Highly Stable Near-Unity Photoluminescence Yield in Monolayer MoS2 by Fluoropolymer Encapsulation and Superacid Treatment, ACS Nano 11, 5179 (2017).
  * [45] H. Bretscher, Z. Li, J. Xiao, D. Y. Qiu, S. Refaely-Abramson, J. A. Alexander-Webber, A. Tanoh, Y. Fan, G. Delport, C. A. Williams, S. D. Stranks, S. Hofmann, J. B. Neaton, S. G. Louie, and A. Rao, Rational Passivation of Sulfur Vacancy Defects in Two-Dimensional Transition Metal Dichalcogenides, ACS Nano 15, 8780 (2021).
  * [46] K. S. Novoselov, D. Jiang, F. Schedin, T. J. Booth, V. V. Khotkevich, S. V. Morozov, and A. K. Geim, Two-Dimensional Atomic Crystals, Proc. Natl. Acad. Sci. 102, 10451 (2005).
  * [47] C. Casiraghi, A. Hartschuh, E. Lidorikis, H. Qian, H. Harutyunyan, T. Gokus, K. S. Novoselov, and A. C. Ferrari, Rayleigh Imaging of Graphene and Graphene Layers, Nano Lett. 7, 2711 (2007).
  * [48] H. Nan, Z. Wang, W. Wang, Z. Liang, Y. Lu, Q. Chen, D. He, P. Tan, F. Miao, X. Wang, J. Wang, and Z. Ni, Strong Photoluminescence Enhancement of MoS2 through Defect Engineering and Oxygen Bonding, ACS Nano 8, 5738 (2014).
  * [49] L. Li and E. A. Carter, Defect-Mediated Charge-Carrier Trapping and Nonradiative Recombination in WSe2 Monolayers, J. Am. Chem. Soc. 141, 10451 (2019).
  * [50] T. Ando, A. B. Fowler, and F. Stern, Electronic Properties of Two-Dimensional Systems, Rev. Mod. Phys. 54, 437 (1982).
  * [51] A. Laturia, M. L. Van de Put, and W. G. Vandenberghe, Dielectric Properties of Hexagonal Boron Nitride and Transition Metal Dichalcogenides: From Monolayer to Bulk, Npj 2D Mater. Appl. 2, 1 (2018).
  * [52] Z. Lin, B. R. Carvalho, E. Kahn, R. Lv, R. Rao, H. Terrones, M. A. Pimenta, and M. Terrones, Defect Engineering of Two-Dimensional Transition Metal Dichalcogenides, 2D Mater. 3, 022002 (2016).
  * [53] Y. Liu, P. Stradins, S.-H. Wei, Y. Liu, P. Stradins, and S.-H. Wei, Air Passivation of Chalcogen Vacancies in Two-Dimensional Semiconductors, Angew. Chemie Int. Ed. 55, 965 (2016).
  * [54] J. Pető, T. Ollár, P. Vancsó, Z. I. Popov, G. Z. Magda, G. Dobrik, C. Hwang, P. B. Sorokin, and L. Tapasztó, Spontaneous Doping of the Basal Plane of MoS2 Single Layers through Oxygen Substitution under Ambient Conditions, Nat. Chem. 10, 1246 (2018).
  * [55] S. Barja, S. Refaely-Abramson, B. Schuler, D. Y. Qiu, A. Pulkin, S. Wickenburg, H. Ryu, M. M. Ugeda, C. Kastl, C. Chen, C. Hwang, A. Schwartzberg, S. Aloni, S. K. Mo, D. Frank Ogletree, M. F. Crommie, O. V. Yazyev, S. G. Louie, J. B. Neaton, and A. Weber-Bargioni, Identifying Substitutional Oxygen as a Prolific Point Defect in Monolayer Transition Metal Dichalcogenides, Nat. Commun. 10, 1 (2019).
  * [56] W. Wang, L. Bai, C. Yang, K. Fan, Y. Xie, and M. Li, The Electronic Properties of O-Doped Pure and Sulfur Vacancy-Defect Monolayer WS2: A First-Principles Study, Mater. 11, 218 (2018).
  * [57] Z. He, X. Wang, W. Xu, Y. Zhou, Y. Sheng, Y. Rong, J. M. Smith, and J. H. Warner, Revealing Defect-State Photoluminescence in Monolayer WS2 by Cryogenic Laser Processing, ACS Nano 10, 5847 (2016).
  * [58] Y. Liu, H. Liu, J. Wang, and D. Liu, Defect-Type-Dependent Carrier Lifetimes in Monolayer WS2 Films, J. Phys. Chem. C 126, 4929 (2022).
  * [59] A. Das, S. Pisana, B. Chakraborty, S. Piscanec, S. K. Saha, U. V. Waghmare, K. S. Novoselov, H. R. Krishnamurthy, A. K. Geim, A. C. Ferrari, and A. K. Sood, Monitoring Dopants by Raman Scattering in an Electrochemically Top-Gated Graphene Transistor, Nat. Nanotechnol. 3, 210 (2008).
  * [60] J. Yang, J. U. Lee, and H. Cheong, Excitation Energy Dependence of Raman Spectra of Few-Layer WS2, FlatChem 3, 64 (2017).
  * [61] W. Zhao, Z. Ghorannevis, K. K. Amara, J. R. Pang, M. Toh, X. Zhang, C. Kloc, P. H. Tan, and G. Eda, Lattice Dynamics in Mono- and Few-Layer Sheets of WS2 and WSe2, Nanoscale 5, 9677 (2013).
  * [62] A. Principi, M. Carrega, M. B. Lundeberg, A. Woessner, F. H. L. Koppens, G. Vignale, and M. Polini, Plasmon Losses Due to Electron-Phonon Scattering: The Case of Graphene Encapsulated in Hexagonal Boron Nitride, Phys. Rev. B 90, 165408 (2014).
  * [63] M. A. Yamoah, W. Yang, E. Pop, and D. Goldhaber-Gordon, High-Velocity Saturation in Graphene Encapsulated by Hexagonal Boron Nitride, ACS Nano 11, 9914 (2017).
  * [64] Y. Fan, A. W. Robertson, Y. Zhou, Q. Chen, X. Zhang, N. D. Browning, H. Zheng, M. H. Rümmeli, and J. H. Warner, Electrical Breakdown of Suspended Mono- and Few-Layer Tungsten Disulfide via Sulfur Depletion Identified by in Situ Atomic Imaging, ACS Nano 11, 9435 (2017).
  * [65] J. W. McPherson, J. Kim, A. Shanware, H. Mogul, and J. Rodriguez, Trends in the Ultimate Breakdown Strength of High Dielectric-Constant Materials, IEEE Trans. Electron Devices 50, 1771 (2003).