# Collapse transition in epidemic spreading subject to detection with limited
resources

Santiago Lamata-Otín  Department of Condensed Matter Physics, University of
Zaragoza, 50009 Zaragoza, Spain  GOTHAM lab, Institute of Biocomputation and
Physics of Complex Systems (BIFI), University of Zaragoza, 50018 Zaragoza,
Spain  Adriana Reyna-Lara  Instituto Tecnológico y de Estudios Superiores de
Monterrey, 64849 Monterrey, N.L., México  David Soriano-Paños  Institute
Gulbenkian of Science (IGC), 2780-156 Oeiras (Portugal)  GOTHAM lab, Institute
of Biocomputation and Physics of Complex Systems (BIFI), University of
Zaragoza, 50018 Zaragoza, Spain  Vito Latora  School of Mathematical Sciences,
Queen Mary University of London, London E1 4NS, United Kingdom  Dipartimento
di Fisica ed Astronomia, Università di Catania and INFN, Catania I-95123,
Italy  Complexity Science Hub Vienna, A-1080 Vienna, Austria  Jesús Gómez-
Gardeñes  Department of Condensed Matter Physics, University of Zaragoza,
50009 Zaragoza, Spain  GOTHAM lab, Institute of Biocomputation and Physics of
Complex Systems (BIFI), University of Zaragoza, 50018 Zaragoza, Spain

###### Abstract

Compartmental models are the most widely used framework for modeling
infectious diseases. These models have been continuously refined to
incorporate all the realistic mechanisms that can shape the course of an
epidemic outbreak. Building on a compartmental model that accounts for early
detection and isolation of infectious individuals through testing, in this
article we focus on the viability of detection processes under limited
availability of testing resources, and we study how the latter impacts on the
detection rate. Our results show that, in addition to the well-known epidemic
transition at ${\mathcal{R}}_{0}=1$, a second transition occurs at
${\mathcal{R}}^{*}_{0}>1$ pinpointing the collapse of the detection system
and, as a consequence, the switch from a regime of mitigation to a regime in
which the pathogen spreads freely. We characterize the epidemic phase diagram
of the model as a function of the relevant control parameters: the basic
reproduction number, the maximum detection capacity of the system, and the
fraction of individuals in shelter. Our analysis thus provides a valuable tool
for estimating the detection resources and the level of confinement needed to
face epidemic outbreaks.

###### pacs:

89.20.-a, 89.75.Hc, 89.75.Kd

##  I Introduction

The COVID-19 pandemic has affected the entire world, causing significant loss
of life, economic hardship, and widespread disruption of social and cultural
norms. As a result, it is essential to understand how communicable diseases
spread and what can be done to mitigate their impact. Mathematical models are
particularly useful in this regard Estrada (2020); Perra (2021), as they allow
researchers to gain valuable insights into the transmission dynamics of
infectious diseases Chinazzi et al. (2020); Pullano et al. (2020), inform
public health policies Arenas et al. (2020); Wong et al. (2020), and guide
efforts to control Meidan et al. (2021); Della Rossa et al. (2020) and prevent
future outbreaks Davis et al. (2021).

Compartmental models are the most widely used framework for modeling
infectious diseases Keeling and Rohani (2011); Anderson and May (1992). In
these models, a population is divided into compartments or states, with
transitions between them mediated by different parameters. While simple and
widely used, these basic models, such as the Susceptible-Infectious-Recovered
model Kermack and McKendrick (1927), have limitations in accounting for
complex demographic or social factors that may impact transmission dynamics
and the effectiveness of containment policies.

Over the past years, researchers have made significant efforts to overcome the
limitations of compartmental models used for modeling infectious diseases. One
approach has been to incorporate various elements that make the models more
realistic, thus broadening their range of applicability. These refinements
cover the use of complex networks to model interactions through which the
pathogen can spread Pastor-Satorras et al. (2015); Wang et al. (2017), thus
allowing to better capture the heterogeneity of the connection between
individuals and its impact on disease transmission dynamics. Another
refinement has been its combination with diffusion processes Barbosa et al.
(2018), which mimic mobility flows between densely populated areas Hazarie et
al. (2021), enabling the development of metapopulation frameworks to analyze
the role of travel and movement patterns in the spread of infectious diseases
Balcan et al. (2009); Tizzoni et al. (2014); Soriano-Paños et al. (2022).
Additionally, researchers have coupled the spreading dynamics of infectious
diseases with behavioral factors that impact the acceptance of interventions
Funk et al. (2010); Wang et al. (2016). This refinement acknowledges the
importance of human behavior in the success of disease control measures and
allows for the exploration of interventions more likely to be adopted by the
population Bauch (2005); Steinegger et al. (2020).

Following this line of research, recently new compartmental models have been
developed to study how early detection and isolation of infectious individuals
through testing and the subsequent activation of contact tracing strategies
can interrupt the advance of transmission chains Reyna-Lara et al. (2021);
Kojaku et al. (2021); Bianconi et al. (2021); Serafino et al. (2022). In this
study, we explore how limited availability of testing resources alter the
viability of detection processes and their impact on the ongoing epidemic
outbreak. We propose a minimal compartmental model that can simulate the
effects of different interventions, such as lockdowns and testing, and derive
the epidemic phase diagram analytically.

Our results show that, in addition to the well-known epidemic transition that
occurs when the basic reproduction number ${\mathcal{R}}_{0}$ is
${\mathcal{R}}_{0}^{c}=1$, a second transition takes place at
${\mathcal{R}}^{*}_{0}>1$, which depends on the maximum detection capacity of
the system. When ${\mathcal{R}}_{0}>{\mathcal{R}}^{*}_{0}$, the system moves
from a phase where detection can mitigate the epidemic outbreak to a phase
where the pathogen spreads freely. By characterizing this transition, we can
determine the precise value of ${\mathcal{R}}^{*}_{0}$ as a function of both
the detection capacity of the system and the fraction of individuals in
shelter. Our model provides a valuable tool for estimating the detection
resources and confinement needed to face epidemic outbreaks and can be adapted
for use in more elaborate models.

##  II Epidemic spreading dynamics

To simulate the time progression of infections during a single epidemic wave
without considering reinfections, we adopt a Susceptible-Infected-Recovered
($SIR$) framework. This compartmental modeling approach divides the population
into three epidemiological states: Susceptible ($S$), i.e. individuals who
lack prior exposure to the pathogen and, therefore, possess no immunity,
Infected ($I$), individuals infected with the pathogen and carrying a
sufficient viral load to infect others, and Recovered ($R$), individuals who
have overcome the infection and have developed immunity. In addition to the
three typical SIR categories, we introduce two further categories: Locked
($L$) and Detected ($D$), so that we can refer to our model as SLIDR. The $L$
category comprises epidemiologically $S$ individuals who are under strict
lockdown measures and thus cannot contract the disease. The $D$ group denotes
those individuals who were Infected but have been detected through testing,
and consequently, are no longer infectious due to their isolation.

As typical in compartmental models, each individual of a population can occupy
only one state at a given time, and the transition from one state to another
is governed by the flow diagram illustrated in Fig. 1.a. Initially, a fraction
$l_{0}$ of the population transitions from $S$ to $L$, while the remaining
individuals in $S$ are susceptible to infection at a rate $\beta$ per contact
with an agent in compartment $I$. Infectious individuals transition to either
recovered, $R$, at a rate $\mu$, or detected, $D$, at a rate $g(t)$. Note that
the detection rate, $g(t)$, is time-dependent, as it is contingent upon the
testing capacity, as we will discuss in details below. Finally, detected
individuals transition to the recovered state at a rate $\gamma$.

Figure 1: In panel (a) we show the flux diagram of the $SLIDR$ model. The
model has five compartments: susceptible ($S$), locked (L), infectious ($I$),
detected ($D$) and recovered ($R$). Arrows indicate the possible transitions
between different states. In panel (b) we sketch the dependence of the
detection rate on the fraction of infected individuals in the population,
given by Eq. (6).

The former transitions allow us to write a set of mean-field equations by
considering that each agent is involved in $\langle k\rangle$ contacts per
unit time in a population of size $N$. Considering the fractions of the
population in each compartment ($s=S/N$, $l=L/N$, $i=I/N$, $d=D/N$, and
$r=R/N$) fulfilling $s+l+i+d+r=1$, the differential equations that govern
their time evolution read as:

$\displaystyle\dot{s}$ $\displaystyle=$ $\displaystyle-\left(\langle
k\rangle-1\right)\beta si\;,$ (1) $\displaystyle\dot{l}$ $\displaystyle=$
$\displaystyle 0\;,$ (2) $\displaystyle\dot{i}$ $\displaystyle=$
$\displaystyle\left(\langle k\rangle-1\right)\beta
si-\left(\mu+g(t)\right)i\;,$ (3) $\displaystyle\dot{d}$ $\displaystyle=$
$\displaystyle-\gamma d+gi\;,$ (4) $\displaystyle\dot{r}$ $\displaystyle=$
$\displaystyle\mu i+\gamma d\;.$ (5)

Note that, as explained above, Eq. (2) implies that the fraction of population
initially set under lockdown remains constant during time ($l(t)=l_{0}$). As
mentioned above, the detection rate $g(t)$ is not constant in order to capture
the limited nature of testing resources. In particular we assume the time
dependence of $g(t)$ shown in Fig. 1.b, whose mathematical expression reads:

$g(t)=\begin{cases}g_{0}&\text{if }\quad i(t)<\theta\;,\\\ g_{0}e^{-\lambda
N(i(t)-\theta)}&\text{if }\quad i(t)>\theta\;.\end{cases}$ (6)

The previous functional form assumes that detection operates in a normal way,
i.e. with a constant rate $g_{0}$, provided that the fraction of infectious
agents remains below a specified capacity threshold $\theta$. Under normal
conditions, tests are readily available, and the entire detection process,
including identification of infectious agents, testing, and processing of
results, takes place optimally within an average time period of $1/g_{0}$.
However, when the number of infected individuals $i(t)$ exceeds the capacity
threshold $\theta$, we assume that a national health system begins to
experience delays, causing a reduction in detection efficiency. Ultimately,
when the detection system becomes too slow, infected individuals may recover
before being detected, so detection does not alter their infectious period. To
model the collapse of the detection system, we introduce an exponential decay
of the original rate, $g_{0}$, towards the detection compartment according to
the difference between the availability threshold $\theta$ and the demand
$i(t)$. The exponential decay is regulated by a tunable decay rate $\lambda$
times the size of the population $N$.

Figure 2: In panels (a)-(b) we show the temporal evolution of all the
compartments of the $SLIDR$ model in absence of lockdown ($l_{0}=0$), with
unlimited resources ($\theta=1$), an average connectivity $\langle
k\rangle=5$, a infectivity $\beta=0.137$, a recovery rate $\mu=1/7$, a
transition rate from $D$ to $R$ regulated by $\gamma=3/20$, and assuming a
baseline detection rate of $g_{0}=0.2$. In panels (c)-(d) we appreciate how
considering limited resources, $\theta=0.07$ (the rest of the parameters are
identical to panels (a)-(b)), gives rise to the acceleration of the infected
cases growth once the capacity threshold is reached. This acceleration depends
on the decay rate $\lambda=1$ and on the population size $N=700$. In panels
(e)-(f) we show how collapse can be avoided when a fraction of population is
under lockdown ($l_{0}=0.1$). $r^{\infty}$ reference curves display how
lockdown reduces the attack rate of the diseases even below the one
corresponding to unlimited resources (panel (a)). The rest of the parameters
are identical to panels (c)-(d) ($\theta=0.07$, $\langle k\rangle=5$,
$\beta=0.137$, $\mu=1/7$, $\gamma=3/20$, $g_{0}=0.2$). Note that in all panels
the basic reproduction number is fixed to $R_{0}=1.6$ and time is measured in
arbitrary units.

In Fig. 2, we study how the interplay between the implementation of lockdown
policies and the existence of limited testing resources affects epidemic
trajectories. To set a reference, we represent an epidemic trajectory in Fig.
2.a-b with a baseline detection rate $g_{0}=0.2$, unlimited resources, i.e.
$\theta=1$, and with no lockdown policies at a play, i.e. $l_{0}=0$. First, we
explore in Fig. 2.c-d the impact of limited resources by setting the capacity
threshold to $\theta=0.07$ and the decay rate $\lambda=1$. The selection of
the capacity threshold value draws inspiration from the typical testing ratios
observed during the most challenging weeks of the COVID-19 epidemic in Europe
COV (2023). The time evolution of each compartment clearly shows that, once
the infectious population reaches the value $\theta$ noticed by the horizontal
dashed red line in Fig. 2.d , the fraction of detected agents decays and, as a
consequence, the increase of the infectious population speeds up. As a
consequence, a much larger attack rate
$r^{\infty}=\lim_{t\rightarrow\infty}r(t)$ is observed. For longer times we
observe that detection has a second peak, pinpointing that after the epidemic
peak the infectious population falls back under the threshold $\theta$.
Finally, Fig. 2.e-f shows what happens when $\theta=0.07$, but a finite
fraction $l_{0}=0.1$ of the susceptible population is under lockdown. In this
case, the pool of susceptible individuals available to be infected is smaller
and, as a result, the capacity threshold $\theta$ is more effective in
stopping the disease spreading, leading to a better mitigation of the outbreak
and, consequently, to a smaller attack rate.

Figure 3: In panel (a) we show the $SLIDR$ attack rate $r^{\infty}$ for two
values of $\lambda$ in the absence of lockdown ($l_{0}=0$). The lower value of
$\lambda$ yields to a second order collapse transition, while the higher
displays a first order transition. The gray curve corresponds to a perfect
mitigation case and the red to the free propagation dynamics. In panel (b) the
complete phase diagram is shown. The horizontal dashed lines indicate the
values of $\lambda$ corresponding to the parameters used to represent curves
in panel (a). The vertical black dashed line indicates the critical value
$\mathcal{R}_{0}^{\star}$ computed according to Eq. (17). In both panels
simulations are performed for the range of infectivity values
$\beta\in[0,0.75]$, assuming a baseline detection rate $g_{0}=1$, a capacity
threshold $\theta=0.07$, a population size $N=700$ with average connectivity
$\langle k\rangle=5$, a recovery rate $\mu=1/7$ and a transition rate from $D$
to $R$ regulated by $\gamma=3/20$.

The overall effect of detection and its limited capacity can be analyzed by
computing the epidemic diagram, i.e. the impact of the contagion wave,
measured by the value of $r^{\infty}$, as a function of the basic reproduction
number of the pathogen, ${\cal R}_{0}$, whose expression for the $SLIDR$ model
considering an initially fully susceptible population is given by:

$\mathcal{R}_{0}=\frac{\beta(\langle k\rangle-1)}{g_{0}+\mu}\;.$ (7)

In Fig. 3.a, we show (thin curves) the epidemic diagrams,
$r^{\infty}({\mathcal{R}}_{0})$, in the case $l_{0}=0$ and a force of
detection characterized by a baseline rate $g_{0}=1$, a capacity threshold
$\theta=0.07$, and two different values of $\lambda$. From these diagrams, we
observe that there exist two transition points. First, the well-known epidemic
threshold at $\mathcal{R}_{0}=\mathcal{R}_{0}^{c}=1$, pinpointing that beyond
this point the infective power $\beta(\langle k\rangle-1)$ is larger than the
effective recovery rate $(g_{0}+\mu)$. In addition to the epidemic threshold,
a second transition point appears at ${\mathcal{R}}_{0}^{\star}>1$ which
corresponds to the collapse transition.

To better illustrate the collapse transition point, we also show the epidemic
diagrams (thick curves) corresponding to free propagation,
$r^{\infty}_{\text{FP}}(\mathcal{R}_{0})$, i.e. in the absence of the
detection policies, and that corresponding to perfect mitigation,
$r^{\infty}_{\text{PM}}(\mathcal{R}_{0})$, which is computed assuming the
availability of unlimited resources ($\theta=1$). With these two phase
diagrams as limiting cases, it is clear that the collapse point
${\mathcal{R}}_{0}^{\star}$ corresponds to the minimum reproduction number
that a pathogen needs to jeopardize the detection capacity and decrease the
mitigation effects of the early removal of infectious agents. Beyond this
point, the epidemic diagram separates from the one corresponding to unlimited
detection resources and eventually reaches that corresponding to null
detection at a value ${\mathcal{R}}_{0}>{\mathcal{R}}_{0}^{\star}$ for which
the mitigation power of detection is completely suppressed.

It is also important to discuss the role played by the decay rate $\lambda$.
Although $\lambda$ does not intervene in the precise value of
${\mathcal{R}}_{0}^{\star}$, it controls the transition between perfect-
mitigation and the free-propagation regimes. In particular, the larger the
values of $\lambda$, the more abrupt the collapse transition becomes. For
large values of $\lambda$, the explosive nature of the collapse transitions
implies that once the system reaches the collapse transition
$\mathcal{R}_{0}^{\star}$, the attack rate will be equivalent to the one
observed in an uncontrolled scenario. The increase of the sharpness of the
collapse transition is clear from Fig. 3.b where $r^{\infty}(\mathcal{R}_{0})$
is reported for a continuous range of $\lambda$ values.

To better characterize the collapse transition shown in Fig. 3 we have
introduced the mitigation efficiency $\eta$, which quantifies the global
impact of the detection collapse. The mitigation effectiveness $\eta$,
compares the actual epidemic diagram $r^{\infty}({\mathcal{R}}_{0})$ for a
given detection force ($g_{0}$, $\theta$, $\lambda$) with that obtained with
the same detection rate $g_{0}$ and unlimited resources ($\theta=1$) and is
defined as:

$\eta=\frac{\int_{\mathcal{R}_{0}}\left(r^{\infty}_{\text{FP}}(\mathcal{R}_{0})-r^{\infty}(\mathcal{R}_{0})\right)d\mathcal{R}_{0}}{\int_{\mathcal{R}_{0}}\left(r^{\infty}_{\text{FP}}(\mathcal{R}_{0})-r^{\infty}_{\text{PM}}(\mathcal{R}_{0})\right)d\mathcal{R}_{0}}\;.$
(8)

It takes values in the range $[0,1]$, with $\eta=0$ when the mitigation effect
is null, and $\eta=1$ when the mitigation attained is the best possible one.

Fig. 4 shows the mitigation effectiveness as a function of $\theta$ and
$\lambda$. The latter parameter becomes more relevant when little resources
are available, yielding a great variance in the parameter $\eta$ when
considering low $\theta$ values. However, as the availability of resources
increases, the transition is delayed and the nature of the transition is less
relevant. This is because for high basic reproduction number diseases, the
free-propagation and the perfect-mitigation curves are quite close, as shown
in Fig. 3.a.

Figure 4: The mitigation effectiveness $\eta$ (color code) is reported as a
function of the capacity threshold $\theta$ and of the decay rate $\lambda$.
Simulations are performed for a large range of infectivity values with no
lockdown polices implemented ($l_{0}=0$). We assume a baseline detection rate
$g_{0}=1$, a population size $N=700$ with average connectivity $\langle
k\rangle=5$, a recovery rate $\mu=1/7$ and a transition rate from $D$ to $R$
regulated by $\gamma=3/20$.

##  III The collapse threshold ${\cal R}_{0}^{\star}$

Once characterized numerically the existence of a collapse transition leading
to the failure of the epidemic containment through detection of infectious
agents, we now proceed to derive analytically the precise value of
${\cal{R}}_{0}^{\star}$. This analysis will shed light into the dependence of
${\cal{R}}_{0}^{\star}$ on the epidemiological parameters characterizing the
spread of the pathogen.

As mentioned earlier, the collapse threshold $\mathcal{R}_{0}^{\star}$ is the
minimum value of the basic reproduction number to keep the epidemic curve
always below the capacity threshold $\theta$. Therefore,
$\mathcal{R}_{0}^{\star}$ can be determined as the value of $\mathcal{R}_{0}$
that gives rise to $i_{\text{max}}=\theta$. Although our $SLIDR$ model
includes two additional compartments $L$ and $D$, we can effectively treat it
as an $SI\bar{R}$ model with a $\bar{R}$ compartment which aggregates the
populations $L$,$D$ and $R$ of the original $SLIDR$ formulation. Since we
assume that the fraction of population in $L$ is contained in the new
compartment $\bar{R}$, we have $\bar{r}_{0}=l_{0}$ and $s_{0}=1-l_{0}-i_{0}$.
Secondly, using our time-continuous approach, we obtain the transition rate
from the $I$ compartment to the new ${\bar{R}}$ compartment by adding the
detection rate $g_{0}$ and the recovery rate $\mu$. This way, we can transfer
the outgoing flow from $I$ to the effective compartment $\bar{R}$, neglecting
the internal dynamics $D\rightarrow R$ within the new effective compartment.
After this reformulation we have :

$\displaystyle\dot{s}$ $\displaystyle=$ $\displaystyle-\left(\langle
k\rangle-1\right)\beta si\;,$ (9) $\displaystyle\dot{i}$ $\displaystyle=$
$\displaystyle\left(\langle k\rangle-1\right)\beta
si-\left(\mu+g_{0}\right)i\;,$ (10) $\displaystyle\dot{\bar{r}}$
$\displaystyle=$ $\displaystyle(\mu+g_{0})i\;,$ (11)

with initial conditions: $s_{0}=1-l_{0}-i_{0}$, $i_{0}\ll 1$, and
$\bar{r}_{0}=l_{0}$. Note that in the formulation of this effective model we
have fixed the detection rate to $g_{0}$ as we are interested in the maximum
possible value of $\mathcal{R}_{0}^{\star}$ not triggering the decrease of the
detection rate during the epidemic evolution.

To obtain an analytical expression of $i_{\text{max}}$ we proceed as usual in
the SIR model by dividing Eq. (9) and Eq. (11), yielding:

$\frac{ds}{d\bar{r}}=-\frac{(\langle
k\rangle-1)\beta}{(g_{0}+\mu)}s=-\mathcal{R}_{0}s\;,$ (12)

which can be integrated with the initial conditions above to obtain the
evolution of infectious agents as a function of the fraction of susceptible
ones:

$i(t)=1-l_{0}-s(t)+\frac{1}{\mathcal{R}_{0}}\log\left(\frac{s(t)}{1-l_{0}}\right)\;,$
(13)

where we have set $\bar{r}_{0}=l_{0}$ and made the approximation $s_{0}\simeq
1-l_{0}$, considering that, at the beginning of an epidemic, the number of
infectious agents is small, $i_{0}\ll 1$.

The implicit expression $i(s)$ in Eq. (13) allows us to determine the density
of susceptibles when the outbreak reaches its peak. Setting $\frac{di}{ds}=0$
we obtain $s(i_{\text{max}})=\mathcal{R}_{0}^{-1}$ and, by inserting this
value into Eq. (13), we finally obtain:

$i_{\text{max}}=1-l_{0}-\frac{1}{\mathcal{R}_{0}}\left[1+\log\left({\mathcal{R}_{0}(1-l_{0})}\right)\right]\;,$
(14)

Figure 5: Phase diagram of the $SLIDR$ model in absence of lockdown
($l_{0}=0$). $\mathcal{R}_{0}^{c}=1$ is the epidemic threshold that separates
the absorbent phase from the active epidemic region. Within the active phase,
the model allows the computation of the minimum amount of resources
$\theta^{\star}$ that we need to avoid collapse, given by (18), as a function
of the basic reproduction number $R_{0}$. Thus, $\mathcal{R}_{0}^{\star}$ is
the critical point related to the transition between that Perfect-Mitigation
Active Phase and the Active Collapse Phase. Some $\mathcal{R}_{0}$ values have
been indicated as reference, with their associated $\theta^{\star}$.

Once derived the expression for $i_{\text{max}}$ we set in Eq. (14) the
condition fulfilled at $\mathcal{R}_{0}^{\star}$, i.e.
$\theta=i_{\text{max}}$, and after some algebra we obtain the implicit
relation for ${\mathcal{R}}_{0}^{*}$, given the values of $\theta$ and
$l_{0}$:

${\mathcal{R}_{0}^{\star}}(1-l_{0})e^{{\mathcal{R}}_{0}^{\star}(1-l_{0})\left[\frac{\theta}{(1-l_{0})}-1\right]}=\frac{1}{e}\;.$
(15)

To solve the former relation for ${\mathcal{R}}_{0}^{\star}$ we perform the
following change of variables $x=(\theta/(1-l_{0})-1)e^{-1}$ and
$y={\mathcal{R}}_{0}^{\star}(1-l_{0})\left[\theta/(1-l_{0})-1\right]$. With
this new variables the former expression reads:

$ye^{y}=x\;,$ (16)

which is a well-known transcendental equation whose solution is the Lambert
function, $y=W(x)$.

In order to assess the validity of the expression in Eq. (16) to derive
${\mathcal{R}}_{0}^{*}$ we recall some of the properties of the Lambert
function. First, we note that the Lambert $W(x)$ function is only real
valuated in the case $x\geq-\frac{1}{e}$. Considering the expression of $x$,
it is easy to derive that the former condition demands that $\theta\geq 0$,
which is automatically satisfied. Additionally, for $x<0$ ($\theta+l_{0}<1$),
the Lambert function possesses two branches, namely $y=W_{0}(x)$ and
$y=W_{-1}(x)$. This particular range of $x$ is of interest, as it corresponds
to situations in which both the capacity threshold and the fraction of locked
individuals are small enough for the collapse transition to show up.

After comparing the numerical solution of ${\mathcal{R}}_{0}^{\star}$ with the
analytical values obtained from the two branches of the Lambert function
predictions, we found that the correct behavior is captured by the
$y=W_{-1}(x)$ branch. Thus, the expression of the collapse threshold can be
finally written as:

${\mathcal{R}}_{0}^{\star}(\theta,l_{0})=\frac{1}{\theta-(1-l_{0})}W_{-1}\left(\frac{\theta-(1-l_{0})}{e(1-l_{0})}\right)\;.$
(17)

In Fig. 3.b we show (see vertical dashed line) that the analytical expression
in Eq. (17) works fairly well in reproducing the precise value of
$\mathcal{R}_{0}$ at which the numerical curve $r^{\infty}(\mathcal{R}_{0})$
detaches from the one corresponding to perfect mitigation,
$r_{\text{PM}}^{\infty}(\mathcal{R}_{0})$, and starts approaching that
corresponding to free propagation, $r_{\text{FP}}^{\infty}(\mathcal{R}_{0})$.

Figure 6:  In panel (a) we show the critical value
$\theta^{\star}(\mathcal{R}_{0},l_{0})$ according to Eq. (18). The blue line
separates the region of values where the collapse transition can occur from
the region where the transition never exists because there is always perfect
mitigation. In panel (b) the critical value
$l_{0}^{\star}(\mathcal{R}_{0},\theta)$ is drawn according to Eq. (19). The
blue line indicates the boundary $\theta+l_{0}=1$ and separates the region
where $R_{0}^{\star}$ is finite from the region where the transition never
exists. Note that the contour lines correspond to non equally spaced
$\mathcal{R}_{0}$ values.

We continue by addressing the question about what is the minimum amount of
resources $\theta$ needed to avoid collapse given that the spreading pathogen
is characterized by $\mathcal{R}_{0}$. This threshold value for $\theta$,
hereafter called $\theta^{*}$, can be straightforwardly calculated by just
imposing $\theta=i_{\text{max}}(R_{0})$ in Eq. (14) which yields:

$\theta^{\star}(\mathcal{R}_{0},l_{0})=(1-l_{0})\left\\{1-\frac{1}{\mathcal{R}_{0}(1-l_{0})}\left[1+\log\left({\mathcal{R}_{0}}(1-l_{0})\right)\right]\right\\}\;.$
(18)

Note that the expression (18) can be obtained as the inverse function of the
equation (17) in the parameter range ($\theta+l_{0}<1$) prescribed above.

To complete our analytical derivations, we will now focus on determining the
minimum lockdown fraction, denoted by $l_{0}$, required to maintain perfect
mitigation when facing a spreading pathogen characterized by the basic
reproduction number $\mathcal{R}_{0}$, assuming a fixed resource capacity
$\theta$.

We obtain the threshold value $l_{0}^{\star}$ by imposing
$\theta=i_{\text{max}}(R_{0})$ in Eq. (14). By constructing the Lambert
function, we obtain the following expression for the threshold value:

$l_{0}^{\star}(\mathcal{R}_{0},\theta)=1+\frac{1}{\mathcal{R}_{0}}W_{-1}\left(-e^{-\theta\mathcal{R}_{0}-1}\right)\;,$
(19)

where $W_{-1}$ denotes the branch of the Lambert function. It is worth noting
that the validity limits of the Lambert function are always satisfied, since
$\theta\mathcal{R}_{0}\geq 0$ by definition of the parameters. Additionally,
the branch $W_{-1}$ is well-defined, as the condition
$e^{-\theta\mathcal{R}_{0}-1}>0$ is always fulfilled.

##  IV Epidemic phase diagrams

The analytical results derived in section III allow us to analyze the epidemic
phase diagram of the system as a function of the relevant parameters: the
infectiousness of the pathogen, $\mathcal{R}_{0}$, the maximum detection
capacity of the system, $\theta$, and the fraction of the population out in
shelter, $l_{0}$. In the first case when $l_{0}=0$ so that no lockdown is
imposed, as shown in Fig. 5, the phase diagram can be plotted as a function of
$\mathcal{R}_{0}$ and $\theta$. The diagram exhibits three possible phases:
the disease-free phase when $\mathcal{R}_{0}<1$, the perfect-mitigation phase
when $\mathcal{R}_{0}>1$ and $\theta>\theta^{*}({\mathcal{R}_{0}})$, and the
collapse active phase provided $\mathcal{R}_{0}>1$ and
$\theta<\theta^{*}({\mathcal{R}_{0}})$. The two former phases are separated by
the curve provided by Eq. (17) (or alternatively Eq. (18)).

The phase diagram in Fig. 5 shows that for pathogens with $\mathcal{R}_{0}>3$
the maximum detection capacity should be $\theta^{*}>0.3$, pointing out that
the sole active detection demands an extraordinary amount of resources to
avoid the collapse phase. Thus, in the following we explore how combining
active detection and lockdown can suffice to achieve the mitigation of the
outbreak. From Eq. (18) it is clear that when lockdown enters into the game
($l_{0}>0$) two beneficial effects show up. First, the basic reproductive
number turns into an lower effective one
$\bar{\mathcal{R}}_{0}={\mathcal{R}}_{0}(1-l_{0})$ and, second, the maximum
capacity detection is effectively increased from $\theta^{*}$ to
$\bar{\theta}^{*}=\theta^{*}/(1-l_{0})$. These two effects combined allow to
increase the area of the perfect mitigation phase in the
$({\mathcal{R}}_{0},\theta)$-plane as shown in Fig. 6.a which reports the
curves $\theta^{*}(\mathcal{R}_{0})$ that separate the perfect mitigation and
the active collapse phases for different values of the fraction $l_{0}$.

The beneficial effects of combining partial lockdown with detection are also
illustrated in Fig. 6.b. Here we show for different values of
$\mathcal{R}_{0}$ the lockdown fraction, $l_{0}^{*}$, needed to remain in the
perfect mitigation phase considering a fixed maximum detection capacity
$\theta$.

##  V Conclusions

The implementation of reliable detection systems is key to ensure that
policies, such as contact tracing and isolation of infectious individuals have
the desired impact on outbreak control. In this study, we have introduced and
studied a compartmental model in which detection resources are limited. Using
a mean field approach, we have characterized the different dynamic regimes of
the system as a function of its parameters.

The most relevant result of this work is the observation of two transitions as
a function of the basic reproductive number: the epidemic
($\mathcal{R}_{0}=1$) and the collapse
($\mathcal{R}_{0}=\mathcal{R}_{0}^{\star}$) transitions. In the latter
transition the health system is unable to meet the demand for detection (for
$\mathcal{R}_{0}>\mathcal{R}_{0}^{\star}$) and we move from a controlled
regime, where detection drives the mitigation of the epidemic outbreak, to a
regime in which the pathogen spreads freely.

The existence of a collapse transition has motivated the analysis of a
combined implementation of detection and lockdowns Peak et al. (2020). We have
observed that the combination of the two strategies can help to avoid the
collapse point, specially for those pathogens with a large $\mathcal{R}_{0}$.
Besides, our results show how, for certain values of the decay constant
$\lambda$, the nature of the collapse transition turns out to be explosive.
This result is striking because it means that, above the collapse transition,
the attack rate can be the same as in the unmitigated dynamics. The way on how
this explosive behavior shows up is also remarkable since it arises in the
active phase, i.e. well beyond the epidemic threshold. Thus, it stands in
contrast to conventional forms of explosivity observed in contagion models
Gómez-Gardeñes (2016); Scarselli et al. (2021), in which deliberate delays in
epidemic onset lead to abrupt transitions deviating from the usual smooth ones
D’Souza et al. (2019).

Finally, it should be noted that the analytical results presented here are
based on a mean-field approach and, thus, some limitation are worth
mentioning. First our model does not account for all observed features of
social connections and pathogen performance. In this context, it would be
worthwhile to refine the mechanisms built into $SLIDR$ to be able to
incorporate real connectivity patterns given by networks of close contacts. In
this area, the inclusion of contact tracing strategies and not only
symptomatic detection could be of particular interest. This approach could
also be applied to reaction-diffusion processes that simultaneously
incorporate mobility flows and contact patterns Cota et al. (2021), paving the
way for the identification of optimal distributions of detection resources Zhu
et al. (2021); Zhang et al. (2022). In addition, the lockdown that complements
detection has been implemented in a stylized way, i.e. starting from the
beginning of the epidemic way rather than being applied in subsequent times.
For this particular scenario, we have checked that the reported results are
robust provided the time elapsed between the start of the epidemic wave and
the lockdown is small enough compared to the time associated to the epidemic
peak. Overall, our model is able to provide analytical insights as a benchmark
for more realistic models that can capture the full range of complexities
involved in infectious disease outbreaks.

## Acknowledgements

S.L.O and J.G.G. acknowledge financial support from the Departamento de
Industria e Innovación del Gobierno de Aragón y Fondo Social Europeo (FENOL
group grant E36-23R) and from grant PID2020-113582GB-I00 funded by
MCIN/AEI/10.13039/501100011033.

## References

  * Estrada (2020) E. Estrada, Physics Reports 869, 1 (2020).
  * Perra (2021) N. Perra, Physics Reports 913, 1 (2021).
  * Chinazzi et al. (2020) M. Chinazzi, J. T. Davis, M. Ajelli, C. Gioannini, M. Litvinova, S. Merler, A. P. y Piontti, K. Mu, L. Rossi, K. Sun, et al., Science 368, 395 (2020).
  * Pullano et al. (2020) G. Pullano, F. Pinotti, E. Valdano, P.-Y. Boëlle, C. Poletto, and V. Colizza, Eurosurveillance 25, 2000057 (2020).
  * Arenas et al. (2020) A. Arenas, W. Cota, J. Gómez-Gardeñes, S. Gómez, C. Granell, J. T. Matamalas, D. Soriano-Paños, and B. Steinegger, Phys. Rev. X 10, 041055 (2020).
  * Wong et al. (2020) G. N. Wong, Z. J. Weiner, A. V. Tkachenko, A. Elbanna, S. Maslov, and N. Goldenfeld, Physical Review X 10, 041033 (2020).
  * Meidan et al. (2021) D. Meidan, N. Schulmann, R. Cohen, S. Haber, E. Yaniv, R. Sarid, and B. Barzel, Nature Communications 12, 220 (2021).
  * Della Rossa et al. (2020) F. Della Rossa, D. Salzano, A. Di Meglio, F. De Lellis, M. Coraggio, C. Calabrese, A. Guarino, R. Cardona-Rivera, P. De Lellis, D. Liuzza, et al., Nature Communications 11, 5106 (2020).
  * Davis et al. (2021) J. T. Davis, M. Chinazzi, N. Perra, K. Mu, A. Pastore y Piontti, M. Ajelli, N. E. Dean, C. Gioannini, M. Litvinova, S. Merler, et al., Nature 600, 127 (2021).
  * Keeling and Rohani (2011) M. J. Keeling and P. Rohani, _Modeling infectious diseases in humans and animals_ (Princeton University Press, New Jersey, 2011).
  * Anderson and May (1992) R. M. Anderson and R. M. May, _Infectious Diseases of Humans: Dynamics and Control_ (Oxford University Presss, 1992).
  * Kermack and McKendrick (1927) W. O. Kermack and A. G. McKendrick, Proceedings of the royal society of london. Series A, Containing papers of a mathematical and physical character 115, 700 (1927).
  * Pastor-Satorras et al. (2015) R. Pastor-Satorras, C. Castellano, P. Van Mieghem, and A. Vespignani, Reviews of Modern Physics 87, 925 (2015).
  * Wang et al. (2017) W. Wang, M. Tang, H. E. Stanley, and L. A. Braunstein, Reports on Progress in Physics 80, 036603 (2017).
  * Barbosa et al. (2018) H. Barbosa, M. Barthelemy, G. Ghoshal, C. R. James, M. Lenormand, T. Louail, R. Menezes, J. J. Ramasco, F. Simini, and M. Tomasini, Physics Reports 734, 1 (2018).
  * Hazarie et al. (2021) S. Hazarie, D. Soriano-Paños, A. Arenas, J. Gómez-Gardeñes, and G. Ghoshal, Communications Physics 4, 1 (2021).
  * Balcan et al. (2009) D. Balcan, V. Colizza, B. Gonçalves, H. Hu, J. J. Ramasco, and A. Vespignani, Proceedings of the National Academy of Sciences 106, 21484 (2009).
  * Tizzoni et al. (2014) M. Tizzoni, P. Bajardi, A. Decuyper, G. Kon Kam King, C. M. Schneider, V. Blondel, Z. Smoreda, M. C. González, and V. Colizza, PLoS Computational Biology 10, e1003716 (2014).
  * Soriano-Paños et al. (2022) D. Soriano-Paños, W. Cota, S. C. Ferreira, G. Ghoshal, A. Arenas, and J. Gómez-Gardeñes, Annalen der Physik 534, 2100482 (2022).
  * Funk et al. (2010) S. Funk, M. Salathé, and V. A. A. Jansen, Journal of the Royal Society, Interface 7, 1247–1256 (2010).
  * Wang et al. (2016) Z. Wang, C. T. Bauch, S. Bhattacharyya, A. d’Onofrio, P. Manfredi, M. Perc, N. Perra, M. Salathé, and D. Zhao, Physics Reports 664, 1–113 (2016).
  * Bauch (2005) C. T. Bauch, Proceedings of the Royal Society B: Biological Sciences 272, 1669 (2005).
  * Steinegger et al. (2020) B. Steinegger, A. Arenas, J. Gómez-Gardeñes, and C. Granell, Physical Review Research 2, 023181 (2020).
  * Reyna-Lara et al. (2021) A. Reyna-Lara, D. Soriano-Paños, S. Gómez, C. Granell, J. T. Matamalas, B. Steinegger, A. Arenas, and J. Gómez-Gardeñes, Phys. Rev. Res. 3, 013163 (2021).
  * Kojaku et al. (2021) S. Kojaku, L. Hébert-Dufresne, E. Mones, S. Lehmann, and Y.-Y. Ahn, Nature Physics 17, 652 (2021).
  * Bianconi et al. (2021) G. Bianconi, H. Sun, G. Rapisardi, and A. Arenas, Physical Review Research 3, L012014 (2021).
  * Serafino et al. (2022) M. Serafino, H. S. Monteiro, S. Luo, S. D. Reis, C. Igual, A. S. Lima Neto, M. Travizano, J. S. Andrade Jr, and H. A. Makse, PLoS Computational Biology 18, e1009865 (2022).
  * COV (2023) European Centre for Disease Prevention and Control, European Union (2023), URL https://www.ecdc.europa.eu/en/publications-data/covid-19-testing.
  * Peak et al. (2020) C. M. Peak, R. Kahn, Y. H. Grad, L. M. Childs, R. Li, M. Lipsitch, and C. O. Buckee, The Lancet Infectious Diseases 20, 1025 (2020).
  * Gómez-Gardeñes (2016) L. L. T. S. e. a. Gómez-Gardeñes, J., Sci Rep 6, 19767 (2016).
  * Scarselli et al. (2021) D. Scarselli, N. B. Budanur, M. Timme, and B. Hof, Nature Communications 12, 2586 (2021).
  * D’Souza et al. (2019) R. M. D’Souza, J. Gómez-Gardeñes, J. Nagler, and A. Arenas, Advances in Physics 68, 123 (2019).
  * Cota et al. (2021) W. Cota, D. Soriano-Paños, A. Arenas, S. C. Ferreira, and J. Gómez-Gardeñes, New Journal of Physics 23, 073019 (2021).
  * Zhu et al. (2021) X. Zhu, Y. Liu, S. Wang, R. Wang, X. Chen, and W. Wang, Applied Mathematics and Computation 411, 126531 (2021).
  * Zhang et al. (2022) K. Zhang, Y. Han, M. Gou, and B. Wang, Physical Review E 105, 064308 (2022).


# Electric-field-driven conductance switching in encapsulated graphene
nanogaps

Eugenia Pyurbeeva<EMAIL_ADDRESS>School of Physics and Astronomy,
Queen Mary University of London, Mile End Road, London, E1 4NS, UK  Jacob L.
Swett  Department of Materials, University of Oxford, Oxford, OX1 3PH, UK
Qingyu Ye  School of Biological and Chemical Sciences, and Materials Research
Institute Queen Mary University of London ,Mile End Road, London E1 4NS, UK
Oscar W. Kennedy  London Centre for Nanotechnology and Department of Physics
and Astronomy, University College London, London, United Kingdom  Jan A. Mol
<EMAIL_ADDRESS>School of Physics and Astronomy, Queen Mary University of
London, Mile End Road, London, E1 4NS, UK

###### Abstract

Feedback-controlled electric breakdown of graphene in air or vacuum is a well-
established way of fabricating tunnel junctions, nanogaps, and quantum dots.
We show that the method is equally applicable to encapsulated graphene
constrictions fabricated using hydrogen silsesquioxane. The silica-like layer
left by hydrogen silsesquioxane resist after electron-beam exposure remains
intact after electric breakdown of the graphene. We explore the conductance
switching behavior that is common in graphene nanostructures fabricated via
feedback-controlled breakdown, and show that it can be attributed to atomic-
scale fluctuations of graphene below the encapsulating layer. Our findings
open up new ways of fabricating encapsulated room-temperature single-electron
nanodevices and shed light on the underlying physical mechanism of conductance
switching in these graphene nanodevices.

††preprint: AIP/123-QED

Electric breakdown of graphene Prins  _et al._ (2011); Lau  _et al._ (2014);
Sadeghi  _et al._ (2015) provides an inexpensive and facile way of creating
two-dimensional nanostructures, including tunnel junctions Prins  _et al._
(2011); Lau  _et al._ (2014); Sadeghi  _et al._ (2015); Limburg  _et al._
(2018) and quantum dots Harzheim  _et al._ (2020). While these structures have
mainly been studied in the context of electrodes for single-molecule devices
Prins  _et al._ (2011), they also hold potential as active electronic elements
in their own right as single electron transistors Barreiro  _et al._ (2012);
Puczkarski  _et al._ (2015); Harzheim  _et al._ (2020) and for phase-change
memory Pósa  _et al._ (2017). A key challenge in harnessing this potential is
the requirement to integrate electric breakdown with additional pre- and post-
processing fabrication steps Han  _et al._ (2014).

As graphene electronic devices approach practical applicability, encapsulation
can be beneficial, reducing the probability of damaging the graphene in
subsequent fabrication steps Li  _et al._ (2019). As 2D materials have a large
surface to volume ratio, surface contamination in the active region of the
device can result in large variations of conductance Ray (2016). While
environment-dependent conductance has succesfully been utilised to detect gas
molecules Ray (2016); Rani and Ray (2019) and DNA Mishra  _et al._ (2020), it
should be avoided in most applications to minimise device-to-device
variability.

Patterning graphene on a silica substrate with hydrogen silsesquioxane (HSQ)
enables direct integration of a silica-like encapsulation layer in the
electron beam lithography process, creating sandwiched graphene
nanostructures. We demonstrate that the electric breakdown process to
fabricate these nanostructures is equally efficient for sandwiched and open
graphene devices.

To benchmark the stability of encapsulated graphene devices, we compare the
observed conductance switching to that previously reported in open devices
Sarwat  _et al._ (2017); Pósa  _et al._ (2017). We attribute the electric-
field-induced switching to atomic-scale fluctuations within the silica-like
encapsulation. While this switching behavior is likely to limit the use of
encapsulated graphene nanostructures for some electronic applications, its
understanding can offer insights into the miniaturisation limits of graphene
electronics, while additionally, it opens a viable route towards the
development of nanoscale non-volatile memory elements Pósa  _et al._ (2017);
Fried  _et al._ (2020); Pósa  _et al._ (2021).

HSQ has been previously used for the fabrication of graphene devices Han  _et
al._ (2007); Chen  _et al._ (2007), however the developed resist is typically
removed after patterning with hydrofluoric acid for contact deposition or
further doping Chen  _et al._ (2007). Recent work has also shown quantum dot
formation through effective doping of graphene with varying exposure dose Wang
_et al._ (2019). It has been demonstrated that HSQ encapsulated devices have
an increased conductance and yield compared to those fabricated with
polymethyl methacrylate (PMMA) resist due to the graphene being protected
during post-exposure processing and the decreased contamination and doping of
graphene by resist residue Li  _et al._ (2019).

Figure 1:  a) 3D rendered AFM image of the fabricated device. A bowtie-shaped
graphene constriction overlaps two pre-fabricated gold electrodes. The
thickness of the residual HSQ layer is approximately 100nm. b) Conductance
dependencies on gate voltage of the non-electroburnt constrictions fabricated
with AR-N 7500 (graphene is uncovered) vs constriction fabricated with HSQ.
The latter shows a significantly lower shift of the Dirac point (conductance
minimum) from zero voltage.

We fabricate devices on heavily doped (0.005 Ohm-cm ) silicon substrates with
a 275 nm dry thermal oxide layer, on which electrodes have been deposited
using electron-beam lithography: first 5nm Cr + 15nm Au with electron beam
evaporation, followed by a photolithography step, using contact lithography
and 5nm Cr + 95nm Au deposited with electron beam evaporation. The wafer was
then completely covered with chemical vapour deposition grown graphene
(Graphenea). Using electron beam lithography, we patterned bowtie-shaped
graphene constrictions overlapping both electrodes with the width of in the
narrowest point of 300nm.

The AFM image of the resulting device is shown on figure 1a, where an
approximately 100 nm thick bowtie-shaped silica layer can be seen overlapping
with the two gold electrodes. The thickness of the residual silica-like layer
covering the graphene constriction is approximately 100 nm.

We measure the dependence of electrical conductance on the back-gate voltage,
applied to the doped silicon substrate, for a graphene constriction fabricated
with HSQ and compare it to the conductance of a constriction of the same
geometry fabricated with an organic negative resist (AR-N 7500) which is
removed after development. The results in Figure 1b show the device fabricated
with HSQ having a higher conductance and a smaller shift of the Dirac point
(corresponding to the conductance minimum) from zero gate voltage. While
absolute conductance is a poor measure of graphene contamination, as it is
affected by the geometry Qi  _et al._ (2014), the shift of the Dirac point is
a more reliable quantity. The Dirac point is typically shifted away from zero
gate voltage in graphene devices due to interaction with the substrate and
contamination and effective doping of graphene with the resist residue Wang
_et al._ (2019). In the device fabricated with HSQ this shift is several times
smaller than in the device fabricated with AR-N 7500 ($\sim 40$V of gate
voltage compared to $\sim 120$V) which indicates that the encapsulated
graphene is less contaminated. Both results agree with the previous studies Li
_et al._ (2019); Wang  _et al._ (2019).

Figure 2:  a) Current-voltage (IV) traces during the feedback-controlled
electroburning process. The first trace (shown in red) shows a dip in
conductance on the upwards ramp of the voltage, which is usual for electric
breakdown. The following ramps result in gradually decreasing conductance –
indicated by arrow. The inset shows sudden spikes in conductance immediately
before the gap formation that are attributed to due to the transition between
two-path and single-path conductance. b) IV trace of an electroburnt device
fitted to the Simmons model. The gap width is approximately 2nm. c) Current
stability diagram of a encapsulated device showing Coulomb blockade. d)
Differential conductance stability diagram of the same device.  Figure 3:
Compared height profiles of a device on which electric breakdown has not and
has been performed taken by AFM where the origin is the the centre of the
constriction. The profile cuts perpendicular to the constriction show that
after electroburning the device remains well-encapsulated. In the cuts
parallel to the constriction a slight dip can be seen in the device after
electric breakdown, coinciding with the centre of the constriction. Close-up
AFM images of the centre of the devices are also shown.

We perform feedback-controlled electric breakdown on the encapsulated graphene
constriction by repeatedly ramping up the voltage across the constriction up
to a point of a sudden drop in conductance, signifying the start of graphene
breakdown Lau  _et al._ (2014). After the conductance drop, the voltage is
quickly ramped down. Multiple voltage ramps lead to gradual narrowing of the
constriction, which can result in the formation of a nanogap with a
characteristic size of of several nanometres Limburg  _et al._ (2018) or a
quantum dot of similar characteristic size Barreiro  _et al._ (2012); Harzheim
_et al._ (2020); Fried  _et al._ (2020), which has been attributed to a small
graphene island remaining in the gap. However, due to the size of the system,
the exact structure of the dot is not known.

The electric breakdown process for the devices fabricated with HSQ is very
similar to the non-encapsulated devices. We performed feedback-controlled
breakdown on both at room temperature in ambient conditions and were able to
use identical parameters – positive and negative ramp rate (5–8 V/s and
2300–5000V/s ), maximum voltage (10V) and target resistance (500 M$\Omega$).
The breakdown success rate was comparable to that of non-encapsulated graphene
devices – 64% exhibited non-linear current-voltage characteristics indicative
of nanogaps compared to 71% in Lau  _et al._ (2014).

A common feature of the electric breakdown process is a significant
conductance enlargement in the end of the last cycle, which is attributed to a
quantum interference effect due to the transition between multi-path and
single path conductivity in the remaining carbon chains bridging the gap
Sadeghi  _et al._ (2015). This effect is also present in encapsulated devices
as shown in Figure 2a. In many cases, several peaks were present, which we
attribute to the repeated formation and destruction of the final carbon
filament.

We fit the current-voltage (IV) traces of several devices after electric
breakdown to the Simmons model Simmons (1963); Zimbovskaya (2013); Limburg
_et al._ (2018) (see figures 2b, 4b,c) and estimate the width of the gap to be
approximately 2 nm, which agrees with the results for devices fabricated with
AR-N 7500 El Abbassi  _et al._ (2017).

In addition to tunnel junctions, some devices show Coulomb diamonds indicative
of the formation of a quantum dot Gehring  _et al._ (2016) – see figures 2
c,d. The estimated lower bound for the quantum dot addition energy extracted
from 2 c,d is 200 meV, which is an order of magnitude greater than the
characteristic thermal energy $k_{\rm{B}}T$for room temperature (approximately
25 meV).

The large addition energy of quantum dots is one of the advantages of
feedback-controlled electric breakdown. In order for Coulomb blockade to be
observed, the energy of the electrostatic interaction of two electrons on the
quantum dot has to be much greater than the characteristic thermal energy
Barreiro  _et al._ (2012):

$\frac{e^{2}}{4\pi\varepsilon_{0}d}\gg k_{\rm{B}}T$ (1)

where $d$ is the diameter of the quantum dot. For Coulomb blockade to be
present at room temperature, this requires $d\ll 50$ nm and therefore $d$ in
the order of nanometres, which is close to the state-of-the-art in current
nanofabrication techniques Manfrinato  _et al._ (2013, 2017). As a result, the
size restriction makes electroburning an attractive, cheap alternative method
of fabricating quantum dots demonstrating room-temperature Coulomb blockade
Fried  _et al._ (2020).

To study the effect of electric breakdown on the encapsulation layer, we have
imaged devices using atomic force microscopy both before and after electric
breakdown and found that the 100 nm thick silica-like top layer is not removed
by the electric breakdown of the graphene underneath, despite the high
temperatures reached during breakdown – in excess of 1000K Pósa  _et al._
(2017); Evangeli  _et al._ (2021). Figure 3 shows height profiles on the
constrictions. After electric breakdown a slight depression in HSQ height is
present in the centre of the constriction, where the temperature is expected
to be the highest due to the highest current density, however with a height
decrease of 5 nm against the layer thickness of more than 100 nm, the device
is still fully encapsulated after the constriction breakdown (figure 3). The
depression is likely the result of melting or densification of the
encapsulating layer.

Figure 4:  a) Current-voltage traces of a device switching between two
conductance states. Two critical bias voltages can be seen: approximately
0.17V for a low to high switch and approximately 0.24V for the reverse
switching. b) Simmons model fits to the two conductance states in figure a.
The high state has a gap width several Ånarrower than the low state,
indicating geometric restructuring at the atomic scale. c) IV trace of an
electroburnt device before the voltage-driven restoration of ordered graphene
and Simmons model fit to it, indicating a gap width of 2nm. d) Conductance-
voltage traces for a graphene device before feedback-controlled electroburning
(top) and two electroburnt devices after a voltage-driven switch to a high
conductance state. The y-axis shows the typical conductances (in units of the
typical graphene device conductance) and the traces are shifted vertically for
ease of presentation.

We observe two types of conductance switching in the encapsulated devices,
both driven by bias voltage ramping. The first one is repeated switching
between two conductance states, both of which are in the tunneling regime. The
second one is consistent with that observed in Pósa  _et al._ (2017); Sarwat
_et al._ (2017) – conductance changes of several orders of magnitude from the
tunnelling regime to bulk conductance, which can only be reversed by repeating
the electric breakdown process. The large conductance changes have been
attributed to either a phase transition, forming a conductive Si crystal in
the substrate Yao  _et al._ (2010); He  _et al._ (2012); Pósa  _et al._ (2017)
or geometric restructuring of graphene, with the formation of carbon chains
and carbon nanotubes closing the gap Sarwat  _et al._ (2017). Both effects are
expected to occur due to the large field strength in the electroburnt gaps –
in the order of $10^{9}$V/m.

An example of the first type of switching is shown in Figure 4(a), with a
device switching between two similar conductive states (shown in green and
blue lines) more than 30 times during bias voltage sweeps from -0.3V to +0.3V.
The device switches from the low (green) to the high (blue) conductive state
when the bias voltage reaches approximately 0.17V, while the reverse switching
is only activated at bias voltages above 0.24V. Following the analysis in Pósa
_et al._ (2021), we have fitted both conductance states at voltage ranges
before switching is activated to the Simmons model (see figure 4b) and find
that the while both states have gap widths of about 2nm, the high state has a
gap width approximately 3Å lower, while the contact area is 20 times smaller
than in the low state. We note that this analysis relies on a continuum model
for tunneling, and is not fully applicable to the atomic scale. However, the
results are indicative for reproducible atomic scale reaarangements in the
nanogap.

In addition to atomic scale fluctuations, we observe conductance changes of
several orders of magnitude. The observed conductance changes are only
reversible by subsequent electric breakdown. We explore the nature of this
“breakdown reversal” by studying the dependence of the conductance of the
device after the “reversal” on the bias voltage. Prior to electroburning a
device conducts through bulk graphene and this function exhibits a dip in
conductance at zero bias Somphonsane  _et al._ (2017), due to weak
localisation Mccann  _et al._ (2006); Fal’ko  _et al._ (2007) and electron-
electron interactions Kozikov  _et al._ (2010), which is dramatic at low
temperatures, and is still observable in our sample at liquid nitrogen
temperature (Figure 4c).

The zero-bias conductance dip can be used as indication of conductance being
mediated by bulk graphene. Figure 2d shows the typical conductance-voltage
trace of a device before the breakdown of the constriction (top) and two types
of conductance-voltage traces of devices after breakdown reversal. The typical
conductance of a “reversed” device observed is $10^{-3}$ of the original
conductance of the graphene constriction and the conductance trace exhibits no
dip.

In rare cases, we have observed the device switching to a state with even
higher conductance: $10^{-2}$ of that of the original graphene constriction,
with a conductance-voltage trace demonstrating the same dip at zero bias
voltage as bulk graphene before electroburning. Thus, breakdown reversal
usually destroys bulk graphene structure. We have fitted the IV trace of a
device, that showed bulk graphene behaviour after the “breakdown reversal”,
just before the change in conductance (Figure 4c) and found the gap width to
be 2 nm – typical of these devices.

We believe that the first type of conductance switching we observed, repeated
switching between two states in the tunnelling regime is best described by
graphene restructuring. From the Simmons fits of both conductance state we
find that the conductance change is due to atomic-scale fluctuations, where
the increase of conductance is associated with both a decrease of the fitted
gap width and contact area. As the two states are highly reproducible (at the
atomic scale), repeated reforming of graphene in the gap, for instance,
formation of a protuberance (see Figure 5a,b) is more likely than the
formation and destruction of a silicon crystal identical on the atomic scale.

The large-scale conductance changes, however, can be explained by either
graphene restructuring, with the forming carbon chains closing the gap, or the
formation of a silicon crystal. Due to the size of the nanogap and the
encapsulation, this question cannot easily be solved by microscopy. It has
been demonstrated that the presence of the zero-bias conductance dip in bulk
exfoliated graphene depends on the deposition method Morozov  _et al._ (2006),
and its disappearance in some samples is attributed to corrugations in
graphene due to imperfect transfer. Therefore, the absence of the zero-bias
conductance dip of the devices after breakdown reversal can be attributed to
either carbon chains reforming in the gap in a disordered way, or conductance
being mediated by a silicon crystal. Similarly, the “reversed” devices that
demonstrated the zero-bias conductance dip can either have carbon bonds that
reformed in an ordered way, forming bulk-like crystaline graphene, (since the
conductance of a graphene nanoribbon is proportional to the width down to
nanometre-sized widths Qi  _et al._ (2014), the width of the graphene flake
closing the gap is a hundredth of the initial constriction width –
approximately 5 carbon chains that are closely packed together to form sheet
graphene), or a very conductive silicon crystal, so that the conductance dip
observed is from the quantum effects in the wider parts of the constriction.

Finally, we note that due to the encapsulation, a new mechanism of graphene
restructuring may be present in our devices, that does not exist in “open”
graphene nanogaps. As graphene breakdown occurs underneath a silica layer,
carbon atoms cannot escape the system (with the formation of CO2, for
instance), and remain in the silica. Therefore, conductance switching may
originate from the formation of new bond with the existing carbon sites in the
encapsulating layer (Figure 5b,d) without large-scale movement of the carbon
atoms, unlike the mechanism described in Sarwat  _et al._ (2017) (Figure
5a,c). Additionally, it is plausible that the restriction of dynamics to two
dimensions favours the reformation of ordered graphene (Figure 5c,d).

Figure 5:  a) The proposed mechanism of conductance switching in a non-
encapsulated device. The final row of carbon atoms “unzips” under an electric
field and can form a chain closing the gap. b) A proposed mechanism of
conductance switching in an encapsulated device. As carbon atoms remain in the
system after electric breakdown, carbon-carbon bonds can form with atoms in
the gap. c) Multiple carbon chains formed in a non-encapsulated gap. As the
dynamics is three-dimensional, a disordered structure is likely to form. d)
Multiple carbon chains formed in an encapsulated gap. Under two-dimensional
dynamics the formation of ordered graphene from adjacent carbon chains is more
likely than in the 3D dynamics.

To conclude, we investigated the feasibility of fabricating graphene
constrictions for further feedback-controlled electric breakdown using HSQ
resist and find that the presence of the passivation layer of silica-like
material remaining after exposure and development of the resist, does not
affect the breakdown process and both nanogaps and quantum dots can be formed.
In turn, the Joule heat dissipated during feedback-controlled breakdown does
not affect the encapsulating silica layer in a detrimental way.

As a result, we are able to fabricate nanometre-scale graphene nanogaps and
quantum dots with addition energies large enough to show room-temperature
Coulomb blockade, using a common negative resist. This approach offers
technological benefits as the devices being covered by a layer of SiO2 of
thickness which can be varied through initial spin coating. These include
protection against environmental conditions and the possibility of placing
additional electrodes over the silica layer, allowing the placement of gate
electrodes much closer to the quantum dot than a back or a side-gate. Side
gate electrodes can’t be placed closer than the original width of the
constriction Puczkarski  _et al._ (2015), approximately 100 nm, while the HSQ
layer can be down to 10nm thick. Close placement of gate electrodes reduces
the voltage needed to gate a single-electron transistor compared to using a
back gate. Additionally, in thermoelectric devices, reducing the distance
between the heater and the device (made possible by the top heater placement)
solves the problem of the large dispersed heater powersGehring  _et al._
(2019).

We have explored conductance switching in HSQ-encapsulated devices and found
that they demonstrate reproducible bias-voltage driven conductance changes in
the tunnelling regime due to atomic-scale fluctuations in the nanogap, as well
as previously described switches from tunnelling to bulk conductance. We
attribute the fluctuations in the encapsulated nanogap to the geometric
restructuring of graphene under electric field. We believe that studying the
physical mechanism behind the conductance switching in nanodevices can provide
insights on both the scaling limits of nanoelectronics and nanoscale memory
elements.

## References

  * Prins  _et al._ (2011) F. Prins, A. Barreiro, J. W. Ruitenberg, J. S. Seldenthuis, N. Aliaga-Alcalde, L. M. K. Vandersypen, and H. S. J. van der Zant, Nano Lett. 11, 4607 (2011), arXiv:1110.2335 .
  * Lau  _et al._ (2014) C. S. Lau, J. A. Mol, J. H. Warner, and G. A. D. Briggs, Phys. Chem. Chem. Phys. 16, 20398 (2014).
  * Sadeghi  _et al._ (2015) H. Sadeghi, J. A. Mol, C. S. Lau, G. A. D. Briggs, J. Warner, and C. J. Lambert, Proc. Natl. Acad. Sci. 112, 2658 (2015), arXiv:1601.00217 .
  * Limburg  _et al._ (2018) B. Limburg, J. O. Thomas, G. Holloway, H. Sadeghi, S. Sangtarash, I. C. Y. Hou, J. Cremers, A. Narita, K. Müllen, C. J. Lambert, G. A. D. Briggs, J. A. Mol, and H. L. Anderson, Adv. Funct. Mater. 28, 1 (2018).
  * Harzheim  _et al._ (2020) A. Harzheim, J. K. Sowa, J. L. Swett, G. A. D. Briggs, J. A. Mol, and P. Gehring, Phys. Rev. Res. 2, 013140 (2020), arXiv:1906.05401 .
  * Barreiro  _et al._ (2012) A. Barreiro, H. S. Van Der Zant, and L. M. Vandersypen, Nano Lett. 12, 6096 (2012).
  * Puczkarski  _et al._ (2015) P. Puczkarski, P. Gehring, C. S. Lau, J. Liu, A. Ardavan, J. H. Warner, G. A. D. Briggs, and J. A. Mol, Appl. Phys. Lett. 107, 133105 (2015).
  * Pósa  _et al._ (2017) L. Pósa, M. El Abbassi, P. Makk, B. Sánta, C. Nef, M. Csontos, M. Calame, and A. Halbritter, Nano Lett. 17, 6783 (2017).
  * Han  _et al._ (2014) S. J. Han, A. V. Garcia, S. Oida, K. A. Jenkins, and W. Haensch, Nat. Commun. 5, 1 (2014).
  * Li  _et al._ (2019) P. Li, D. Collomb, and S. Bending, Mater. Lett. 257, 126765 (2019).
  * Ray (2016) S. J. Ray, Sensors Actuators, B Chem. 222, 492 (2016).
  * Rani and Ray (2019) S. Rani and S. J. Ray, Carbon N. Y. 144, 235 (2019).
  * Mishra  _et al._ (2020) S. Mishra, S. Rani, and S. J. Ray, J. Appl. Phys. 128 (2020), 10.1063/5.0016104.
  * Sarwat  _et al._ (2017) S. G. Sarwat, P. Gehring, G. Rodriguez Hernandez, J. H. Warner, G. A. D. Briggs, J. A. Mol, and H. Bhaskaran, Nano Lett. 17, 3688 (2017).
  * Fried  _et al._ (2020) J. P. Fried, X. Bian, J. L. Swett, I. I. Kravchenko, G. A. D. Briggs, and J. A. Mol, Nanoscale 12, 871 (2020).
  * Pósa  _et al._ (2021) L. Pósa, Z. Balogh, D. Krisztián, P. Balázs, B. Sánta, R. Furrer, M. Csontos, and A. Halbritter, npj 2D Mater. Appl. 5, 1 (2021).
  * Han  _et al._ (2007) M. Y. Han, B. Özyilmaz, Y. Zhang, and P. Kim, Phys. Rev. Lett. 98, 1 (2007).
  * Chen  _et al._ (2007) Z. Chen, Y. M. Lin, M. J. Rooks, and P. Avouris, Phys. E Low-Dimensional Syst. Nanostructures 40, 228 (2007), arXiv:0701599 [cond-mat] .
  * Wang  _et al._ (2019) Z. Wang, Y. Yuan, X. Liu, J. Sun, M. Muruganathan, and H. Mizuta, ACS Nano 13, 7502 (2019).
  * Qi  _et al._ (2014) Z. J. Qi, J. A. Rodríguez-Manzo, A. R. Botello-Méndez, S. J. Hong, E. A. Stach, Y. W. Park, J. C. Charlier, M. Drndić, and A. T. Johnson, Nano Lett. 14, 4238 (2014).
  * Simmons (1963) J. G. Simmons, J. Appl. Phys. 34, 1793 (1963).
  * Zimbovskaya (2013) N. a. Zimbovskaya,  _Part Springer Tracts Mod. Phys. B. Ser._, Vol. 254 (2013) pp. 1–275.
  * El Abbassi  _et al._ (2017) M. El Abbassi, L. Pósa, P. Makk, C. Nef, K. Thodkar, A. Halbritter, and M. Calame, Nanoscale 9, 17312 (2017).
  * Gehring  _et al._ (2016) P. Gehring, H. Sadeghi, S. Sangtarash, C. S. Lau, J. Liu, A. Ardavan, J. H. Warner, C. J. Lambert, G. A. D. Briggs, and J. A. Mol, Nano Lett. 16, 4210 (2016).
  * Manfrinato  _et al._ (2013) V. R. Manfrinato, L. Zhang, D. Su, H. Duan, R. G. Hobbs, E. A. Stach, and K. K. Berggren, Nano Lett. 13, 1555 (2013).
  * Manfrinato  _et al._ (2017) V. R. Manfrinato, A. Stein, L. Zhang, C.-Y. Nam, K. G. Yager, E. A. Stach, and C. T. Black, Nano Lett. 17, 4562 (2017).
  * Evangeli  _et al._ (2021) C. Evangeli, E. McCann, J. L. Swett, S. Tewari, X. Bian, J. O. Thomas, G. A. D. Briggs, O. V. Kolosov, and J. A. Mol, Carbon N. Y. 178, 632 (2021).
  * Yao  _et al._ (2010) J. Yao, Z. Sun, L. Zhong, D. Natelson, and J. M. Tour, Nano Lett. 10, 4105 (2010).
  * He  _et al._ (2012) C. He, Z. Shi, L. Zhang, W. Yang, R. Yang, D. Shi, and G. Zhang, ACS Nano 6, 4214 (2012).
  * Somphonsane  _et al._ (2017) R. Somphonsane, H. Ramamoorthy, G. He, J. Nathawat, C. P. Kwan, N. Arabchigavkani, Y. H. Lee, J. Fransson, and J. P. Bird, Sci. Rep. 7, 1 (2017).
  * Mccann  _et al._ (2006) E. Mccann, K. Kechedzhi, V. I. Fala’Ko, H. Suzuura, T. Ando, and B. L. Altshuler, Phys. Rev. Lett. 97, 14 (2006), arXiv:0604015 [cond-mat] .
  * Fal’ko  _et al._ (2007) V. I. Fal’ko, K. Kechedzhi, E. McCann, B. Altshuler, H. Suzuura, and T. Ando, Solid State Commun. 143, 33 (2007).
  * Kozikov  _et al._ (2010) A. A. Kozikov, A. K. Savchenko, B. N. Narozhny, and A. V. Shytov, Phys. Rev. B - Condens. Matter Mater. Phys. 82, 1 (2010), arXiv:1004.0468 .
  * Morozov  _et al._ (2006) S. V. Morozov, K. S. Novoselov, M. I. Katsnelson, F. Schedin, L. A. Ponomarenko, D. Jiang, and A. K. Geim, Phys. Rev. Lett. 97, 7 (2006).
  * Gehring  _et al._ (2019) P. Gehring, M. van der Star, C. Evangeli, J. J. Le Roy, L. Bogani, O. V. Kolosov, and H. S. J. van der Zant, Appl. Phys. Lett. 115, 073103 (2019), arXiv:1907.02815 .