[datatype=bibtex, overwrite=true] [fieldsource=booktitle,
match=.*Interspeech.*, replace=Proc. Interspeech] [fieldsource=journal,
match=.*INTERSPEECH.*, replace=Proc. Interspeech] [fieldsource=booktitle,
match=.*icassp_inpress.*, replace=ICASSP (in press)] [fieldsource=booktitle,
match=.*International.*Conference.*on.*Acoustics,.*Speech.*and.*Signal.*Processing.*,
replace=Proc. ICASSP] [fieldsource=booktitle,
match=.*International.*Conference.*on.*Learning.*Representations.*,
replace=ICLR] [fieldsource=booktitle,
match=.*International.*Conference.*on.*Machine.*Learning.*, replace=ICML]
[fieldsource=booktitle,
match=.*Automatic.*Speech.*Recognition.*and.*Understanding.*, replace=Proc.
ASRU] [fieldsource=booktitle, match=.*Spoken.*Language.*Technology.*,
replace=Proc. SLT] [fieldsource=booktitle,
match=.*Speech.*Synthesis.*Workshop.*, replace=Proc. SSW]
[fieldsource=booktitle, match=.*workshop.*on.*speech.*synthesis.*,
replace=Proc. SSW] [fieldsource=booktitle,
match=.*Advances.*in.*neural.*information.*processing.*, replace=Proc. NIPS]
[fieldsource=booktitle,
match=.*Advances.*in.*Neural.*Information.*Processing.*, replace=Proc. NIPS]
[fieldsource=booktitle, match=.*Workshop.*on.* Applications.* of.*
Signal.*Processing.*to.*Audio.*and.*Acoustics.*, replace=Proc. WASPAA]
[fieldsource=booktitle,
match=.*International.*Conference.*on.*Language.*Resources.*and.*Evaluation.*,
replace=Proc. LREC] [fieldsource=journal,
match=.*Spontaneous.*Speech.*Processing.*and.*Recognition, replace=Proc. SSPR]
[fieldsource=publisher, match=.+, replace=] [fieldsource=month, match=.+,
replace=] [fieldsource=location, match=.+, replace=] [fieldsource=address,
match=.+, replace=] [fieldsource=organization, match=.+, replace=]
[fieldsource=doi, match=.+, replace=] [fieldsource=url, match=.+, replace=]
[fieldsource=editor, match=.+, replace=] The University of Tokyo, Japan.

# Do learned speech symbols follow Zipf’s law?

###### Abstract

In this study, we investigate whether speech symbols, learned through deep
learning, follow Zipf’s law, akin to natural language symbols. Zipf’s law is
an empirical law that delineates the frequency distribution of words, forming
fundamentals for statistical analysis in natural language processing. Natural
language symbols, which are invented by humans to symbolize speech content,
are recognized to comply with this law. On the other hand, recent
breakthroughs in spoken language processing have given rise to the development
of learned speech symbols; these are data-driven symbolizations of speech
content. Our objective is to ascertain whether these data-driven speech
symbols follow Zipf’s law, as the same as natural language symbols. Through
our investigation, we aim to forge new ways for the statistical analysis of
spoken language processing.

Index Terms—  speech analysis, Zipf’s law, generative spoken language model,
speech representation

##  1 Introduction

Zipf’s law, a well-known empirical principle, delineates the frequency of
occurrence of elements within a dataset [1]. Specifically, when the occurrence
frequency of an element ranks as the $k$-th highest within a dataset, it
equates to $1/k$ of the frequency of the most occurring element. This law is
observed to be applicable across various data domains, with natural language
symbols (e.g., words), follow this pattern [2]. To illustrate, the third most
frequent word in an English document, “and”, appears approximately one-third
as often as the most frequent word,
“the”111https://www.cs.cmu.edu/~cburch/words/top.html. In the practice of
analyzing the occurrence frequency of natural language symbols within a text
corpus, one can denote the frequency rank as $r$ and its respective occurrence
frequency as $f_{r}$. Consequently, the following relationship is given by
Zipf’s law:

$\displaystyle f_{r}=ar^{-\eta},$ (1)

where, $a$ and $\eta$ are model parameters. The text corpus follows Zipf’s law
when $\eta\approx 1$, and follows a power law otherwise. In essence, Zipf’s
law is a specific type of power law. When following a power law, the log-log
plot of rank against frequency appears linear.

By examining the adherence to or deviation from Zipf’s law, one can analyze
the distinct characteristics of a text corpus. This analysis finds notable
applications in natural language processing, as illustrated by the following
examples [3]:

  * Infants’ language acquisition: The vocabulary that 2- to 4-year-olds acquire tends to gravitate towards high-frequency words, a phenomenon indicated by a rank-frequency distribution that is convex, rather than linear [4].

  * Writing system variations across languages: Writing systems significantly differ from one language to another, ranging from sound-based phonograms to semantic-based logograms. Tanaka observed that as notation shifts from phonogram to logogram, the rank-frequency distribution transitions to follow a power law; phonographic languages display a convex, while logographic languages demonstrate a linear on a plot [3].

  * Quantification of communication effort: Zipf’s law is also recognized as the principle of least effort [5]. This principle stipulates that the utilization of frequently used words minimizes speaking and listening efforts during human communication. Viewing communication through this perspective allows for the potential measurement of the naturalness exhibited in communication, including emergent communication [6] and machine-generated communication [7].

Meanwhile, recent advancements in deep learning, such as self-supervised
learning, have facilitated the discovery of discrete symbol representations
[8, 9, 10]. These representations are learned in a data-driven manner from
speech. In contrast to traditional signal-processing-based representations,
such as mel-spectrograms, the symbol sequences learned through these methods
captures a wealth of phonetic and semantic contents. While natural language
symbols were crafted by humans to encode the contents of speech, the learned
speech symbols can be seen as their data-driven counterparts, that means,
symbols crafted in a data-driven manner.

With this context in mind, we are led to a research question:

“Does Zipf’s law, which holds true for natural language symbols, also apply to
the learned speech symbols?”

Verifying this hypothesis could potentially unveil the capability to extend
the statistical analysis techniques, traditionally employed in natural
language analysis, to the spoken language analysis. Moreover, since this
methodology bypasses the necessity for transcriptions, it might forge a path
towards textless analysis applicable to a diverse range of sounds,
encompassing non-verbal vocalizations and non-speech sounds (further details
are elaborated in Section 5).

In this study, we conduct experiments to address the posed question, utilizing
the generative spoken language model (GSLM) [8], a variant of speech symbol
representation methodology. Initially, we perform a foundational experiment to
ascertain whether the learned speech symbols follow Zipf’s law. Leveraging a
speech corpus comprising paired text and speech, we analyze symbol rank-
frequency distributions by correlating them with the accompanying text.
Subsequently, we explore whether our textless analysis identifies deviation in
non-textual contents of speech. Specifically, we aim to capture variations in
language fluency between native and non-native speech utterances. Our efforts
endeavor to pave a way for comprehensive speech analysis.

##  2 Generative spoken language model (GSLM)

Fig. 1: Generative spoken language model.

Figure 1 illustrates the GSLM [8], an analysis-synthesis system operating
through discrete speech symbols. This system is structured into three modules:
speech2unit, unit language model, and unit2speech. However, this study focuses
solely on the utilization of the speech2unit. This specific module integrates
a pre-trained self-supervised learning (SSL) model, featuring technologies
such as contrastive predictive coding [11], wav2vec2.0 [12], HuBERT [13], in
conjunction with a $k$-means clustering model. The SSL model remains fixed,
while the $k$-means clustering model is trained using a speech corpus. Within
this setup, the speech2unit module transforms a speech waveform into a
sequence of discrete symbols.

Symbols are generated at intervals of $20$ ms, this interval generally shorter
than the phoneme duration. Consequently, the speech2unit frequently predicts
identical symbols. To prevent this redundancy, we have opted to consolidate
consecutive identical symbol sequences into a single representation. For
instance, if the speech2unit outputs a sequence like $[3,3,3,50,200,200]$
(where the numbers indicate symbol indices), we simplify this to $[3,50,200]$.
It’s important to note that the speech2unit, particularly the $k$-means
clustering model, is highly sensitive to language of input speech. Hence, it’s
imperative that the model and the encoded speech align in terms of language.

##  3 Methodology

We delineate two methodologies employing GSLM speech symbols to explore Zipf’s
law.

Examining the applicability of Zipf’s law to speech symbols. This methodology
is conducted utilizing text-speech pairs. In the case of text, we employ
natural language symbols, represented either as words or as character
$n$-grams. Words are extracted through morphological analysis to discern the
original form of the words, while character $n$-grams involve sequences of $n$
consecutive characters. For speech, we obtain speech symbols as outlined in
Section 2, calculate speech symbol $n$-grams, that is, sequences of $n$
consecutive speech symbols. We then calculate the ratio of the lengths of
natural language symbol sequences to speech symbol sequences to determine the
value of $n$. To describe this intuitively, $n$ signifies the average number
of speech symbols corresponding to a single natural language symbol. We then
assess whether Zipf’s law holds true for speech symbol $n$-grams within text
adhering to the law.

Identifying non-textual deviations from Zipf’s law. As explained in Section 1,
identifying deviations from this law can be a method to pinpoint non-standard
word usage. In this methodology, we explore to detect non-standard speech
patterns based on rank-frequency distribution. By comparing distributions of
standard and non-standard speech, we aim to identify deviations in non-
standard speech from the standard speech.

##  4 Experimental evaluation

###  4.1 Experimental condition

We utilized HuBERT [13], which was trained on LibriSpeech [14], as the SSL
model in the GSLM speech2unit. The language-specific $k$-means clustering
models were trained using JSUT/JVS [15], J-KAC [16], and J-MAC [17] for
Japanese, and LibriSpeech for English. The number of classes was set to $200$
for both languages. The models, implemented using the fairseq toolkit [18],
are publicly available222https://huggingface.co/nonmetal/gslm-japanese
(Japanese)[19]333https://github.com/facebookresearch/fairseq/tree/main/examples/textless_nlp/gslm
(English). The speech sampling frequency was set at $16$ kHz, and speech
symbols were extracted every $20$ ms. Given that high-frequency and low-
frequency items often deviate from Zipf’s law[3], we estimated model
parameters $a$ and $\eta$ based only on the top $0.1$ % to $10$ % of the
frequencies. These model parameters were determined using the least square
method. For comparative analysis, we fixed $\eta$ at $1.0$ while estimating
only $a$, indicating a strict adherence to Zipf’s law in the rank-frequency
distribution. To reduce the data size of figures, we thinned the data to be
plotted.

To verify Zipf’s law in Section 4.2, we utilized approximately $7,600$
Japanese utterances from JSUT [15] and $13,000$ English utterances from
LJSpeech [14]. For word tokenization, MeCab444https://taku910.github.io/mecab/
and NLTK555https://www.nltk.org/ served as morphological analyzers. In the
character $n$-gram analysis, we used a set of characters encompassing
Chinese/Japanese characters and marks for Japanese, and lowercase alphabets,
symbols (e.g., ”:”, ”?”), and whitespace for English. The average number of
characters per word was $1.6$ (ja) and $5.1$ (en), while the average number of
speech symbols per character stood at $5.7$ (ja) and $1.9$ (en). Moreover, the
average number of speech symbols per word was $8.9$ (ja) and $9.0$ (en). The
value of $n$ in both character $n$-gram and speech symbol $n$-gram was set to
the ceiling of these values, for instance, $\lceil 1.6\rceil=2$ for Japanese
characters.

For identifying non-textual deviation in Section 4.3, we employed native and
Japanese-accented English utterances from UME-
ERJ666https://research.nii.ac.jp/src/en/UME-ERJ.html. The corpus contained
about $20$ native and $200$ non-native speakers, each reading approximately
300–500 English sentences. These sentences varied between speakers but
maintained a balanced phoneme distribution. Non-native speakers were assigned
a language fluency score on a five-point scale. Based on these scores, we
categorized non-native speakers into three groups: low-
($\mathrm{score}<3.0$), mid- ($3.0\leq\mathrm{score}<3.5$), and high-level
($\mathrm{score}\geq 3.5$). These groups consisted of $74$, $66$, and $46$
speakers, respectively. We aggregated speech symbols separately for each non-
native speaker group and the native speakers. To equalize data size across
groups, we randomly selected $10,000$ utterances per group. For symbol
encoding, we used the English GSLM speech2unit.

###  4.2 Verifying Zipf’s law of speech symbols

We verify the law via three steps: word, character $n$-gram, and speech symbol
$n$-gram.

####  4.2.1 Word

First, we verify that the words in the corpora we used adhere to Zipf’s law.
Figure 2 illustrates the rank-frequency distributions in both Japanese and
English. These distributions appear to be linear and the value of $\eta$ is
close to $1.0$. Therefore, we can say that the words follow Zipf’s law, and
that the corpora used have the universal statistics in terms of the word
distribution. A minor observation, as highlighted in previous studies [3], is
that the high-frequency ($\mathrm{rank}<20$) and low-frequency
($10^{3}<\mathrm{rank}$) items deviate from the regression line.

Fig. 2: Word rank-frequency distributions. left: Japanese, right: English.

####  4.2.2 Character $n$-gram

Next, we explore another natural language symbol: the character $n$-gram.
Figure 3 shows these distributions. When analyzed for the same value of $n$,
differences between languages become evident. In Japanese, the distribution
attains linear at $n=2$, whereas in English, it displays convex. The English
distribution gradually transitions to a linear shape as $n$ increases, and at
$n=6$, it remains convex but closely resembles a linear trend. At $n=2$ for
Japanese and $n=6$ for English that correspond to a word, the distributions of
character $n$-grams more closely follow a power law, rather than Zipf’s law.

Fig. 3: Character $n$-gram rank-frequency distributions. left: Japanese,
right: English. $n=2$ in Japanese and $n=6$ in English correspond a word.

####  4.2.3 Speech symbol $n$-gram

Finally, we explore whether speech symbols follow Zipf’s law by comparing the
results with those of words and character $n$-grams. Figure 4 shows these
distributions.

Difference between languages. When analyzed for the same value of $n$, we find
that the distributions are almost identical. These results suggest that the
distributions of consecutive speech symbols, which indicate the frequency of
speech segment use, are language-independent, at least between Japanese and
English. The distributions for English are slightly shifted upwards compared
to those for Japanese. This discrepancy is due to the difference in data size;
the English corpus is twice as large as the Japanese corpus.

Comparison in $n$ corresponding to a character. $n=6$ for Japanese and $n=2$
for English corresponds to a character. We observe that distributions where
$n$ corresponds to a character maintain shapes similar to those of character
1-gram. Specifically, the distribution for $n=6$ in Japanese is nearing
linearity (but not $\eta=1.0$), while the distribution for $n=2$ in English is
convex (Although not illustrated in Figure 3, the distributions of character
1-grams resemble those of the 2-gram distributions: they are linear for
Japanese and convex for English.). This difference can be traced back to
differences in writing systems, as discussed in Section 1. Namely, Japanese
and English are more akin to logographic and phonographic languages,
respectively. The closer a symbol is to representing semantic information
rather than phonetic one, the more linear the distribution becomes. Our
findings suggest that this trend holds true for speech symbols and that the
speech symbol $n$-gram might statistically reflect the sound or meaning of a
character. The ability to attain character statistics without relying on
characters also facilitates analysis based on Zipf’s law (and power law)
without involving characters.

Comparison in $n$ corresponding to a word. $n=9$ signifies for both Japanese
and English. As $n$ escalates, the distributions gravitate towards linear (not
$\eta=1.0$). We theorize that this phenomenon occurs because, as noted
previously, $n$-grams tend to represent semantic information as $n$ increases.
$\eta$ not being $1.0$, but the distribution at $9$-grams is linear, mirroring
the word distribution observed in Figure 2. This suggests that our methodology
holds promise for statistically analyzing features associated with words
without utilizing words themselves.

Fig. 4: Speech symbol $n$-gram rank-frequency distributions. left: Japanese,
right: English. $n=6$ in Japanese and $n=2$ in English correspond to a
character. $n=9$ corresponds a word.

###  4.3 Identifying non-textual deviation

Fig. 5: Speech symbol $n$-gram frequency of native and non-native speech.
“low,” “mid,” and “high” indicate English fluency level of non-native
speakers. For clear illustration, we choose line plots rather than scatter
plots.

Figure 5 displays the speech symbol $n$-gram ($n={1,3,5,7}$) frequencies in
both native and non-native speakers. We delve into the distinctions between
native and non-native speakers based on the data presented in the figure.

Linearization by increasing $n$. The distribution tends to linearize as $n$
increases, a trend observed in both native and non-native speakers. Since the
corpus utilized in this experiment employs the same reading text for both
groups, language fluency does not influence the linguistic content.
Consequently, the linearization is anticipated to be driven by the textual
information (particularly the semantic information hypothesized in Section
4.2), rather than language fluency.

Deviation of non-native speakers. When we examine the deviations of non-native
speakers compared to native speakers, we notice deviations correlated with
language fluency. Specifically, non-native speakers tend to use high-frequency
symbols more frequently. These findings indicate that Zipf’s law (and the
power law) can be employed to discern differences between standard and non-
standard speech with regards to language fluency. Further exploration is
required to ascertain the specific nature of these deviations. Intriguingly,
speakers with high-level proficiency diverge more from native speakers. This
counters the intuitive expectation that individuals with low-level fluency
would deviate more significantly.

##  5 Conclusion

In this paper, we investigated whether speech symbols follow Zipf’s law.
Through our experiments, we determined that: 1) speech symbol $n$-grams
corresponding to a word follow a power law rather than Zipf’s law, and 2) non-
textual deviations in non-standard speech can be identified through the power
law of the rank-frequency distributions.

Our research paves the way for textless analysis methods applicable to various
audio data. The following are possible directions for future research.

  * •

Language development. As outlined in Section 1, the Zipf’s (or power) law can
potentially aid in analyzing infants’ vocabulary acquisition. Despite the
existing challenges in robust automatic speech recognition of infants’ voices
[20, 20], our approach can encode voices into symbols without transcriptions.
This enables statistical analyses utilizing voices exclusively.

  * •

Animal voices and non-speech audio. Zipf’s law has been demonstrated to apply
to symbolic audio, including animal calls [21] and music scores [22]. While
current research often relies on human-invented symbols and annotations, our
method could potentially extend to audio sources beyond human speech without
using the human-invented symbols and annotations. General-purpose audio
representation models [23] could be employed for this purpose.

  * •

Emergent speech communication. The analysis of emergent languages constitutes
a significant area of research, with objectives including 1) fostering
language-based machine-machine communication (i.e., communication between
artificial intelligence) [7, 24], and 2) evaluating the extent to which these
communications mirror statistics of human-to-human communications [6, 25].
Given that Zipf’s law (and related principles, such as the Zipf’s law of
abbreviation [26]) embodies the principle of least effort in human
communication, it may offer a way to explore machine-machine communication
[27]. Our method holds the potential to forge a way in examining speech-based
machine-machine communication, potentially benefitting developments of
exploring human spoken language emergence and human-machine communication.

## References

  * [1] George Kingsley Zipf “Human Behaviour and the Principle of Least Effort”, 1949
  * [2] Steven T. Piantadosi “Zipf’s word frequency law in natural language: A critical review and future directions” In _Psychonomic Bulletin & Review_ 21.5, 2009, pp. 1112–1130
  * [3] Kumiko Tanaka-Ishii “Statistical Universals of Language: Mathematical Chance vs. Human Choice [Japanese Edition]”, 2021
  * [4] Lieven Elena, Salomo Dorothé and Tomasello Michael “Two-year-old children’s production of multiword utterances: A usage-based analysis” In _Cognitive Linguistics_ 20.3, 2009
  * [5] Guido M. Linders and Max M. Louwerse “Zipf’s law revisited: Spoken dialog, linguistic units, parameters, and the principle of least effort” In _Psychon Bull Rev_ , 2023
  * [6] Satwik Kottur et al. “Natural Language Does Not Emerge ‘Naturally’ in Multi-Agent Dialog” In _Proceedings of the 2017 Conference on Empirical Methods in Natural Language Processing_ , 2017, pp. 2962–2967
  * [7] Serhii Havrylov and Ivan Titov “Emergence of Language with Multi-Agent Games: Learning to Communicate with Sequences of Symbols” In _Proceedings of the 31st International Conference on Neural Information Processing Systems_ , NIPS’17, 2017, pp. 2146–2156
  * [8] Kushal Lakhotia et al. “On Generative Spoken Language Modeling from Raw Audio” In _Transactions of the ACL_ 9, 2021, pp. 1336–1354
  * [9] Joonyong Park et al. “How Generative Spoken Language Modeling Encodes Noisy Speech: Investigation from Phonetics to Syntactics” In _INTERSPEECH_ , 2023
  * [10] Zalán Borsos et al. “AudioLM: A Language Modeling Approach to Audio Generation”, arXiv 2209.03143, 2022
  * [11] Aäron Oord, Yazhe Li and Oriol Vinyals “Representation Learning with Contrastive Predictive Coding”, arXiv preprint, arXiv:1807.03748, 2018
  * [12] Alexei Baevski et al. “wav2vec 2.0: A Framework for Self-Supervised Learning of Speech Representations” In _Proc. NIPS_ 33, 2020, pp. 12449–12460
  * [13] Wei-Ning Hsu et al. “HuBERT: Self-Supervised Speech Representation Learning by Masked Prediction of Hidden Units” In _IEEE/ACM Transactions on Audio, Speech, and Language Processing_ 29, 2021, pp. 3451–3460
  * [14] Vassil Panayotov et al. “Librispeech: An ASR corpus based on public domain audio books” In _Proc. ICASSP_ , 2015, pp. 5206–5210
  * [15] Shinnosuke Takamichi et al. “JSUT and JVS: Free Japanese voice corpora for accelerating speech synthesis research” In _Acoustical Science and Technology_ 41.5, 2020, pp. 761–768
  * [16] Wataru Nakata et al. “Audiobook Speech Synthesis Conditioned by Cross-Sentence Context-Aware Word Embeddings” In _Proc. SSW_ , 2021, pp. 211–215
  * [17] Shinnosuke Takamichi et al. “J-MAC: Japanese multi-speaker audiobook corpus for speech synthesis” In _Proc. Interspeech_ , 2022, pp. 2358–2362
  * [18] Myle Ott et al. “fairseq: A Fast, Extensible Toolkit for Sequence Modeling” In _Proceedings of the Conference of the North American Chapter of the Association for Computational Linguistics_ , 2019, pp. 48–53
  * [19] Joonyong Park et al. “Analysis of Degraded Noisy Voice and Application to Other Languages Using Generative Spoken Language Model” in Japanese In _Annual meeting of acoustic society of Japan_ , 2023
  * [20] Christopher Gebauer et al. “Exploiting Diversity of Automatic Transcripts from Distinct Speech Recognition Techniques for Children’s Speech” In _Proc. INTERSPEECH_ , 2023, pp. 4578–4582
  * [21] Arik Kershenbaum et al. “Shannon entropy as a robust estimator of Zipf’s Law in animal vocal communication repertoires” In _Methods in Ecology and Evolution_ 12.3, 2021, pp. 553–564 eprint: https://besjournals.onlinelibrary.wiley.com/doi/pdf/10.1111/2041-210X.13536
  * [22] Juan I. Perotti and Orlando V. Billoni “On the emergence of Zipf ’s law in music” In _Physica A: Statistical Mechanics and its Applications_ 549, 2020, pp. 124309
  * [23] Daisuke Niizumi et al. “BYOL for Audio: Self-Supervised Learning for General-Purpose Audio Representation” In _International Joint Conference on Neural Networks (IJCNN)_ , 2021
  * [24] Angeliki Lazaridou, Alexander Peysakhovich and Marco Baroni “Multi-Agent Cooperation and the Emergence of (Natural) Language” In _ICLR_ , 2017
  * [25] Laura Harding Graesser, Kyunghyun Cho and Douwe Kiela “Emergent Linguistic Phenomena in Multi-Agent Communication Games” In _Proceedings of the 2019 Conference on Empirical Methods in Natural Language Processing and the 9th International Joint Conference on Natural Language Processing (EMNLP-IJCNLP)_ , 2019, pp. 3700–3710
  * [26] George Kingsley Zipf “The psycho-biology of language: An introduction to dynamic philology”, 2013
  * [27] Ryo Ueda and Koki Washio “On the Relationship between Zipf’s Law of Abbreviation and Interfering Noise in Emergent Languages” In _Proceedings of the 59th Annual Meeting of the Association for Computational Linguistics and the 11th International Joint Conference on Natural Language Processing: Student Research Workshop_ , 2021, pp. 60–70


# Stochastic resetting in a networked multiparticle system with correlated
transitions

Oriol Artime  CHuB Lab, Fondazione Bruno Kessler, Via Sommarive 18, 38123 Povo
(TN), Italy<EMAIL_ADDRESS>
###### Abstract

The state of many physical, biological and socio-technical systems evolves by
combining smooth local transitions and abrupt resetting events to a set of
reference values. The inclusion of the resetting mechanism not only provides
the possibility of modeling a wide variety of realistic systems but also leads
to interesting novel phenomenology not present in reset-free cases. However,
most models where stochastic resetting is studied address the case of a finite
number of uncorrelated variables, commonly a single one, such as the position
of non-interacting random walkers. Here we overcome this limitation by framing
the process of network growth with node deletion as a stochastic resetting
problem where an arbitrarily large number of degrees of freedom are coupled
and influence each other, both in the resetting and non-resetting (growth)
events. We find the exact, full-time solution of the model, and several out-
of-equilibrium properties are characterized as function of the growth and
resetting rates, such as the emergence of a time-dependent percolation-like
phase transition, and first-passage statistics. Coupled multiparticle systems
subjected to resetting are a necessary generalization in the theory of
stochastic resetting, and the model presented herein serves as an
illustrative, natural and solvable example of such a generalization.

  *

##  1 Introduction

The mechanism of randomly resetting a dynamical variable to a particular set
of values is widespread in nature, such as in quantum mechanics, stochastic
thermodynamics, chemical reactions and movement ecology [1]. Stochastic
resetting can be convenient for a number of reasons, and rooted in diverse
origins, e.g., given by evolutionary advantage or designed to fulfill specific
purposes. Such a mechanism can lead to the optimization of some task with
respect to the reset-free scenario, e.g., in search strategies [2], enzymatic
catalysis [3], rare event sampling [4], or Internet congestion reduction
strategies [5], among many others [1]. From a fundamental viewpoint, in recent
years the focus has been put on unveiling the properties of the relaxation to
the stationary state, of the first-passage times and of the survival
statistics in systems following diverse dynamics that undergo stochastic
resetting [6, 7, 8, 9, 10, 11]. Of particular note are the emergence of a
nonequilibrium steady state and the existence of an optimal value of the
resetting rate that minimizes the mean first-passage time (MFPT) to a target,
properties that usually hold true for any embedding dimension [12, 13, 14],
for several underlying dynamics followed by the particle [2, 15, 16, 17] and
for most inter-resetting time statistics [18, 19, 20, 21, 22]. See [1] for a
recent review.

Beyond particles diffusing in regular and heterogeneous topologies, resetting
can be actually worth exploring in any dynamical system with a stochastically
evolving state [23, 24, 9, 10]. An interesting example due to its theoretical
and practical relevance is complex network growth undergoing stochastic
resetting. Indeed, networks are abstract representations of complex systems,
where nodes represent individual units and edges encode the interactions among
these units [25], e.g., neurons and synapses in the brain, people and
acquaintance ties in a social network, or web pages and hyperlinks in the
world wide web. Apart from characterizing the topological structure of
empirical networks [26, 27, 28, 29, 30, 31, 32, 33, 34], and study the
behavior of dynamical models on top of them [35, 36, 37, 38, 39, 40], it is
crucial to understand how simple growing mechanisms yield the emergence of
topological features present in real interconnected systems. Diverse
mechanisms are able to predict or explain a plethora of these features: link
rewiring leads to small-worldness [41], node addition with preferential
attachment linking leads to scale-free networks [42, 43, 44, 45], or
heterogeneous node fitness leads to winner-takes-all (or Bose-Einstein
condensation) phenomena [46].

Users in a social network can deactivate their accounts, websites in the WWW
can be deleted, and neurons in the brain can be injured and stop functioning.
In fact, many empirical networks might shrink, and even crash. Hence, a
realistic ingredient to take into consideration in modeling network evolution
is node removal [47, 48, 49, 50, 51, 52, 53, 54, 55, 56]. It turns out that
this effect can be framed as a general stochastic resetting problem that
provides a conceptual advancement with respect to more traditional
multiparticle stochastic resetting models, namely, the fact that the
stochastic events on a particle influence the state of an arbitrary and
variable number of other particles. Recently, there have been some
contributions that take into account the particle interactions in the reset
probability, but in an indirect way, see e.g., [57, 58, 59, 60, 61]. In this
article we show that such a direct interaction can be mathematically handled
and we derive some consequences from it: the emergence of a macroscopic
network structure via a time-dependent phase transition, and the emergence of
an inflection point in the mean first-passage time.

The remainder of this paper is structured as follows. First, we explain in
detail the model of network growth with node removal. We then write a master
equation accounting for its probabilistic description, which is exactly
solved. Afterward, the time-dependent percolation phase transition is examined
and its critical point is derived. Finally, the analysis of the first-passage
distribution and the mean first-passage time are considered. We close the
article drawing some conclusions.

##  2 Model definition

We consider a set of $N$ interconnected nodes, each one characterized by its
degree $k$, i.e., the number of undirected connections to its neighbors. Two
processes compete in the formation of the network. On the one hand, a link
between two nodes picked uniformly at random is added at rate $\alpha N/2$. In
the absence of resetting, the average degree would grow at rate $\alpha$, and
if $N\to\infty$ this corresponds to a kinetic formulation of the Erdős-Rényi
model with parameter (ratio between the number of links $M(t)$ at time $t$ and
the number of possible links between $N$ pair of nodes) $p=\alpha t/N$, where
$t$ is time [62]. On the other hand, with rate $rN$, a node selected uniformly
at random is removed, along with all of its edges, and a new node with no
connections is added into the network. Notice the coupling induced by the
edges in both the growth and removal process: the increase of a degree unit
occurs simultaneously in two nodes, and a deletion of a node with degree $k$
implies a degree unit loss in other $k$ nodes. See Fig. 1 for a sketch of this
process. With this model we prioritize the analytical understanding of the
coupled multiparticle stochastic resetting process over a faithful description
of network evolution. The incorporation of some realistic effects is possible
—e.g., the higher chances of well-connected nodes to acquire new links, the
correlations between the degree and the removal probability, a variable system
size $N$, or the fact that resetting event could make a node to lose not
necessarily all its connections—, but they make the analytical treatment much
more limited.

##  3 Time-dependent degree distribution

We can describe the temporal behavior of the model by studying the
probabilities $p_{k}(t)$, the so-called degree distribution, that give the
fraction of nodes with degree $k$ at time $t$ [49]. These probabilities follow
the master equations

$\mathop{}\\!\frac{\mathrm{d}p_{k}}{\mathrm{d}t}=\alpha p_{k-1}-\alpha
p_{k}-rp_{k}+r(k+1)p_{k+1}-rkp_{k}+\delta_{k,0}r.$ (1)

The first two terms are associated to the uniformly random addition of a link,
while the rest correspond to the resetting process. In specific, the third
term encodes the direct resetting of a node of degree $k$, while the fourth
and the fifth terms correspond to the loss of one link due to the removal of a
neighbor. The last term stands for the incoming flux of nodes that are reset
to degree $0$. Negative degrees are not allowed, so the boundary condition is
$p_{-1}(t)=0$. As initial condition, we choose $p_{k}(0)=\delta_{k,0}$, that
is, initially there are no links in the network.

The exact full-time solution to the master equations (1) can be obtained by
means of the generating function. Introducing the time-dependent degree
generating function $g(z,t)=\sum_{k=0}^{\infty}z^{k}p_{k}(t)$, if we multiply
Eqs. (1) by $z^{k}$ and sum over all degrees, we obtain

$\mathop{}\\!\frac{\partial g}{\partial t}=\left[\alpha
z-(\alpha+r)\right]g+r(1-z)\mathop{}\\!\frac{\partial g}{\partial z}+r,$ (2)

with the conditions $g(1,t)=1$, which derives from the normalization of the
degree distribution, and $g(z,0)=1$, which derives from the initial condition.
Employing the method of characteristics, the solution reads

$g(z,t)=\frac{1}{1-z}\left[e^{\alpha
z/r}\mathcal{G}\left((1-z)e^{-rt}\right)+\frac{r}{\alpha}\right],$ (3)

where $\mathcal{G}(x)=(x-r/\alpha)\exp\left[\alpha/r\,(x-1)\right]$. To obtain
the time-dependent degree distribution, we need to rewrite Eq. (3) as a power
series of the auxiliary variable $z$. To do so, we use the well-known
expansions $(1-x)^{-1}=\sum_{k=0}^{\infty}x^{k}$ and
$e^{x}=\sum_{k=0}^{\infty}x^{k}/k!$, along with the relation

$\sum_{k=0}^{\infty}z^{k}\sum_{m=0}^{\infty}\frac{(xz)^{m}}{m!}=e^{x}\sum_{k=0}^{\infty}z^{k}Q(k+1,x),$
(4)

where $Q(a,b)\equiv\gamma(a,b)/\Gamma(a)$ is the regularized Gamma function,
with $\Gamma(a)$ the complete Gamma function and

$\gamma(a,b)=\int_{b}^{\infty}\mathrm{d}x\,x^{a-1}e^{-x}$ (5)

the upper incomplete Gamma function. After some simple algebra, we arrive at
our desired solution

$p_{k}(t)=\frac{r}{\alpha}\left[1-Q(k+1,\,c(t))\right]+e^{-c(t)-rt}\frac{c(t)^{k}}{k!},$
(6)

where we have introduced the function $c(t)=\frac{\alpha}{r}(1-e^{-rt})$ to
ease the notation.

Figure 1: In $(a)$, stylized time evolution of the number of links of node
$i$, $k_{i}(t)$, and of the total number of links in the system, $M(t)$. The
vertical drops in $M$ correspond to the resetting event of a node losing all
of its links. In the first drop, one of $i$’s neighbor is reset, while in the
second drop $i$ itself is reset. The last reset event does not affect on $i$
because it does not occur in its immediate vicinity. In $(b)$, resetting event
on a node and the impact on its neighbors, which lose one degree unit.

Although the full-time solution (6) is a complicated expression, we can
extract some useful information under certain limits. Let us consider first
the long-time limit, where we have the stationary degree distribution

$p_{k}^{\text{st}}=\frac{r}{\alpha}\left[1-Q\left(k+1,\frac{\alpha}{r}\right)\right].$
(7)

By virtue of the relation

$1-Q\left(k+1,\frac{\alpha}{r}\right)=e^{-\alpha/r}\sum_{m=k+1}^{\infty}\frac{\left(\alpha/r\right)^{m}}{m!},$
(8)

which can be proved by integration by parts, Eq. (7) becomes

$p_{k}^{\text{st}}=\frac{r}{\alpha}e^{-\alpha/r}\sum_{m=k+1}^{\infty}\frac{\left(\alpha/r\right)^{m}}{m!}\approx\left(\frac{\alpha}{r}\right)^{k}\frac{e^{-\alpha/r}}{(k+1)!},$
(9)

where in the last step we keep only the smallest value of the sum, which is at
the same time the largest contribution. Using the Stirling’s formula for the
factorial, we obtain the asymptotic behavior of the stationary degree
distribution

$p_{k}^{\text{st}}\approx
e^{-\alpha/r}\left(\frac{e\,\alpha}{r}\right)^{k}(k+1)^{-3/2-k},$ (10)

which decays much faster than an exponential [49]. In Fig. 2$(a)$ we verify
that both the analytical time-dependent degree distribution and its stationary
values coincide with the values coming from the simulations of the model.

Another interesting limit to explore is the one of small reset rate $r$. On
the one hand, if we set $r=0$ and solve the master equations (1), the time-
dependent degree distribution is $p_{k}(t)=(\alpha t)^{k}e^{-\alpha t}/k!$.
That is, in the reset-free scenario there is no stationary state, and the
network grows indefinitely. On the other hand, if we take the limit of small
resetting rates in the time-dependent solution (6), we get

$\frac{(\alpha t)^{k}e^{-\alpha t}}{k!}+r\left[\frac{t(\alpha
t-k-2)}{2}\frac{(\alpha t)^{k}e^{-\alpha
t}}{k!}+\frac{1-Q(k+1,t\alpha)}{\alpha}\right]+\mathcal{O}\left(r^{2}\right).$
(11)

We see that only when $r=0$ Eq. (11) reduces to the reset-free degree
distribution. As far as the resetting rate becomes finite, this no longer
holds true and the system is able to tend to a well-defined stationary state.
Thus, the introduction of a resetting rate, however small it may be, induces
the emergence of a nonequilibrium steady-state.

Figure 2: In $(a)$, temporal evolution of the degree distribution. We show the
curves for several degrees, indicated in the caption, together with their
stationary value, computed from Eq. (7). Markers come from simulations and
solid lines come from Eq. (6). In $(b)$, temporal evolution of the mean degree
and its variance. Markers come from simulations and lines from the theory. In
all simulations we use $N=1000$, $\alpha=1$, $r=0.5$ and averages are computed
over $200$ independent realizations.

Additional information of the process of network growth under stochastic
resetting can be gained with the temporal evolution of the moments of the
degree distributions $\langle k^{n}\rangle(t)\equiv\sum_{k}k^{n}p_{k}(t)$.
Although they could be computed directly from the full-time expression (6), a
significantly easier approach is to solve the differential equation governing
them. We can solve exactly up to arbitrary order $n$ because the differential
equations for the moments are closed. For example, after some trivial algebra,
for the mean degree we obtain

$\mathop{}\\!\frac{\mathrm{d}\langle k\rangle}{\mathrm{d}t}=\alpha-2r\langle
k\rangle,$ (12)

with initial condition $\langle k\rangle(0)=0$. The solution reads

$\langle k\rangle(t)=\frac{\alpha}{2r}\left(1-e^{-2rt}\right),$ (13)

and it is displayed in Fig. 2$(b)$, along with the standard deviation, which
can be trivially calculated.

##  4 Percolation transition

The time-dependent degree distribution and their moments allow us to study
network-wide connectivity properties during the growth process. As an
illustration, we herein consider the emergence of the giant component, which
is the order parameter of the percolation phase transition and is of crucial
importance, because connectivity is usually assumed the first proxy for
network functionality. Let us denote $u(t)$ the probability that, at time $t$,
a node is not in the giant component via one of its links. If that node has
$k$ connections, the probability to belong to the giant component is then
$1-u(t)^{k}$. Averaging over the degree distribution, we obtain the size of
the giant component

$S(t)=\sum_{k=0}^{\infty}p_{k}(t)\left(1-u(t)^{k}\right)=1-g(u,t),$ (14)

where $g(z,t)$ is the time-dependent degree generating function used to solve
Eq. (1). To solve Eq. (14), we need to know what is the value of $u(t)$. This
can be obtained in a recursive manner, by noting that the condition for a node
to not belong to the giant component following one of its links is that the
node at the end of the link we are following does not belong to the giant
component via any of its other $k-1$ neighbors, i.e., $u(t)^{k-1}$. To account
for the network heterogeneity, this latter quantity needs to be averaged in a
similar fashion to the computation of the giant component $S(t)$, although
without the use of the degree distribution. Indeed, let $q_{k}(t)$ denote the
probability that a node at the end of a link has $k$ connections at time $t$.
It then follows that $q_{k}(t)$ is proportional to $kp_{k}(t)$. The
proportional factor must ensure the normalization of this distribution, that
is, $q_{k}(t)=kp_{k}(t)/\langle k\rangle(t)$. Hence, we reach to our self-
consistent equation for $u(t)$,

$u(t)=\frac{1}{\langle
k\rangle(t)}\sum_{k=0}^{\infty}k\,p_{k}(t)\,u(t)^{k-1}=\frac{\partial_{u}g(u,t)}{\partial_{u}g(1,t)}.$
(15)

Note that $u(t)=1$ is always a solution to Eq. (15) and, if plugged in Eq.
(14), we obtain $S(t)=0$. This absence of the giant component is of course
consistent with the physical meaning of $u(t)$, and it could have been guessed
a priori. When a second, non-trivial solution $u(t)\neq 1$ exists, we find a
finite $S(t)$. The appearance of this second solution marks the existence of a
percolation-like phase transition, and occurs when the derivatives of both
sides of Eq. (15), evaluated at $u(t)=1$, are equal. In our case, the
condition for criticality is given by the combination of values $(t,\alpha,r)$
that satisfy

$1=\frac{2\alpha}{3r}\frac{1-2e^{-rt}+e^{rt}}{1+e^{rt}}.$ (16)

Figure 3: In $(a)$, size of the giant component as a function of time and the
reset rate. The growth rate is $\alpha=1$. Separating the percolating and non-
percolating phases, there is the critical line $t_{c}(\alpha,r)$, from Eq.
(16). In $(b)$, temporal evolution of the size of the giant component for
different values of $\alpha$ and $r$. Solid lines come from the analytical
approximations Eqs. (14) and (15), markers come from simulations. In all
curves, we have considered $N=20000$ and averages over $10$ realizations.

We graphically investigate the time-dependent percolation framework developed
above, in Fig. 3$(a)$. We show the size giant component $S$ as a function of
time and of the resetting parameter, for fixed $\alpha=1$. We see that as
resetting becomes more and more probable, the giant component needs more time
to emerge and its stationary value tends to be smaller. There is a value of
resetting rate, $2\alpha/3$, above which the giant component never emerges, no
matter how long we wait. That is, the minimum average degree (at stationarity)
for a network to display a percolating structure is $3/4$. This is somehow
counterintuitive at first sight, because we would expect, on average, at least
one link per node in order to find a connected subgraph traversing a large
fraction of the network. Indeed, for the Erdős-Rényi model the giant component
is found only if $\langle k\rangle>1$. This conundrum is solved by looking at
higher moments of the degree distribution: if $\alpha$ is large enough, or $r$
is small enough —they need to satisfy $\alpha/r>3/2$— the fluctuations from
the mean value become relevant and create nodes with large enough degrees, so
the giant component can emerge even if, on average, nodes have less than one
neighbor. It is not difficult to verify that the Moloy-Reed criterion [63] for
the existence of the giant component precisely leads to $\langle
k\rangle^{\text{st}}>3/4$.

To further evince the non-trivial impact of the growth and resetting rate on
the percolation transition, in Fig. 3$(b)$ we show the time evolution of the
size of the giant component for several pairs $(\alpha,r)$. We see that if
their ratio is the same, the curves tend to the same stationary value, as
expected from the long-time limit of the degree distribution and the mean
degree. However, the absolute values of $(\alpha,r)$ do have a strong impact
on the critical point and on the time scale to reach the stationary value, as
it can be appreciated in the ordering of the curves. We see, moreover, that
theory and simulations match very well in all the cases.

##  5 First-passage statistics

A final point that we address is the first-passage properties of this coupled
multiparticle system. Let $q_{k}(t)$ be the probability that at time $t$ a
randomly selected node has degree $k$ without having arrived at degree
$k^{*}>k$, and let $a_{k^{*}}(t)$ be the probability that at time $t$ a
randomly selected node has arrived at degree $k^{*}$. Mathematically, this is
equivalent to setting an absorbing boundary at $k^{*}$, so there is no outflow
probability from that state. The master equations for these quantities are

$\displaystyle\mathop{}\\!\frac{\mathrm{d}q_{k}}{\mathrm{d}t}=$
$\displaystyle\alpha q_{k-1}-\alpha q_{k}-rq_{k}+r(k+1)q_{k+1}-rkq_{k}+$
$\displaystyle+\delta_{k,0}r(1-a_{k^{*}})\quad\text{for $0\leq k\leq
k^{*}-2$}.$ (17)
$\displaystyle\mathop{}\\!\frac{\mathrm{d}q_{k}}{\mathrm{d}t}=$
$\displaystyle\alpha q_{k-1}-\alpha q_{k}-rq_{k}-rkq_{k}\quad\text{for
$k=k^{*}-1$}.$ (18)
$\displaystyle\mathop{}\\!\frac{\mathrm{d}a_{k}}{\mathrm{d}t}=$
$\displaystyle\alpha q_{k-1}\quad\text{for $k=k^{*}$}.$ (19)

The initial conditions are $q_{k}(0)=\delta_{k,0}$. Clearly, the first-passage
time probability to the target degree $k^{*}$ is given by $f_{k^{*}}(t)=\alpha
q_{k^{*}-1}$, with the mean first-passage time being then $\langle
t_{k^{*}}\rangle=\int\mathop{}\\!\mathrm{d}t\,tf_{k^{*}}$. We show in Fig.
4$(a)$ that the first-passage distribution obtained from the numerical
solution of Eqs. (17)–(19) match well with the results from simulations.

It is a non-trivial task to obtain the solution of the above set of equations.
In order to proceed and find some analytical results regarding the first-
passage statistics, we apply some approximations whose validity will be later
checked. In particular, we transform the set of discrete master equations
(17)–(19) into a single Fokker-Planck equation, approximating the degree as a
continuous variable (see Appendix for the derivation). If we define $q(k,t)$
as the continuous version of $q_{k}(t)$, that is, $q(k,t)$ is a probability
density function that gives the probability of finding a randomly chosen node
with degree $k\in\left[k,k+\mathop{}\\!\mathrm{d}k\right]$ at time
$t\in\left[t,t+\mathop{}\\!\mathrm{d}t\right]$ that has not yet arrived at
$k^{*}$, we obtain

$\mathop{}\\!\frac{\partial q}{\partial
t}=-\mathop{}\\!\frac{\partial}{\partial
k}\left[v(k)q(k,t)\right]+\mathop{}\\!\frac{\partial^{2}}{\partial
k^{2}}\left[D(k)q(k,t)\right]-rq(k,t)+r\delta(k),$ (20)

with the drift coefficient $v(k)=(\alpha-rk)$ and the diffusion coefficient
$D(k)=(rk+\alpha)/2$. The boundary conditions are $q(k^{*},t)=0$ and
$\partial_{k}q(k,t)|_{k=0}=0$, and the initial condition is
$q(k,0)=\delta(k)$. The survival probability is defined such that

$S(k^{*},t)=\int_{0}^{k^{*}}\,\mathop{}\\!\mathrm{d}k\,q(k,t),$ (21)

that relates to the mean first-passage time as

$\langle
T_{k^{*}}\rangle=-\int_{0}^{\infty}\mathop{}\\!\mathrm{d}t\,t\mathop{}\\!\frac{\partial
S(k^{*},t)}{\partial t}=\tilde{S}(k^{*},s=0),$ (22)

where
$\tilde{S}(k^{*},s)=\int_{0}^{\infty}\mathop{}\\!\mathrm{d}te^{-st}S(k^{*},t)$
is the Laplace transform of the survival probability.

Figure 4: In $(a)$, first-passage time distribution to target degrees
$k^{*}=4,5,\ldots,11$. Solid lines correspond to $f_{k^{*}}$ computed from the
numerical integration of Eqs. (17)–(19), while markers come from simulations.
The growth and reset parameters $\alpha=2$ and $r=0.5$. In $(b)$, mean first-
passage time computed from the analytical expression obtained in the Appendix
(solid lines) and from simulations (markers). From left to right, curves
corresponding to $(\alpha,r)=(2,2.5)$, $(2,0.5)$ and $(2,0.2)$. In all
simulations, we use networks with $N=100$ nodes.

In the Appendix, we obtain a solution for the $\tilde{S}$ in terms of special
functions. However, a direct inspection of Equation (20) already provides
information that would hardly be obtained from the discrete master equations
or grasped by analyzing directly the analytical solution. Given the growth and
reset parameters, we see that the drift term vanishes at
$k_{\text{ip}}=\alpha/r$. Those nodes with degrees $k<k_{\text{ip}}$
experience a positive drift toward $k_{\text{ip}}$, while, on the contrary,
those with degrees $k>k_{\text{ip}}$ suffer a negative drift toward
$k_{\text{ip}}$. Consequently, we anticipate different behaviors depending on
whether the target degree $k^{*}$ is larger or smaller than $k_{\text{ip}}$.
In particular, in the regime $k^{*}<k_{\text{ip}}$, the creation of links
dominates over the resetting of nodes, and the arrival to the target degree
will be relatively fast. In the regime $k^{*}>k_{\text{ip}}$, however, the
resetting dominates and the system does not create connections quick enough in
order for a node to easily reach the target degree. Thus, we expect a
considerable increase in the mean first-passage time. In fact, because the
drift is degree-dependent, we expect that the larger the difference
$k^{*}-k_{\text{ip}}>0$, the more difficult will be to arrive to the target
degree and the longer it will take. Likewise, the larger the difference
$k_{\text{ip}}-k^{*}>0$, the faster the target will be attained, hence
observing a decrease in the MFPT. Note, however, that we should not wait for a
non-trivial combination of $(\alpha,r)$ that minimizes the mean first-passage
time, because our system is now bounded and the non-resetting transitions are
unidirectional. Hence, it is natural to suppose that $\langle
T_{k^{*}}\rangle$ shows a concave-to-convex inflection point at
$k^{*}=k_{\text{ip}}$ if $\alpha>r$. When $\alpha<r$ the inflection point will
disappear and the mean first-passage time will be purely convex over all
traget degrees $k^{*}$. In Fig. 4$(b)$, we verify that the analytical
approximation for the MFPT reproduces well the results from simulations.
Indeed, we find an overall good agreement, especially in the concave region.
When we enter in the convex regime, where resetting dominates, we observe that
some slight differences between the prediction and simulations arise, that can
be associated with the failure of the mean-field approach in capturing the
full complexity of the growth and reset network process.

##  6 Conclusions

In summary, we have approached a central problem in network science by
reformulating it as a stochastic resetting process. In doing so, we have
provided insights that are relevant to both areas of research. On the one
hand, network growth under node removal advances our knowledge of stochastic
resetting, providing an analytically amenable model characterized by many-body
variable interactions. In it, the particles’ state evolution is non-trivially
coupled in the degree space: growth events can be seen as a two-particle
coupling, with a simultaneous increase of one degree unit; reset events can be
seen as a $(k+1)$-particle coupling, where one node in state $k$ resets and
$k$ other nodes lose one degree unit. On the other hand, we have studied
several out-of-equilibrium properties of network growth with node removal, as
it is customary in stochastic resetting problems.

In particular, we have obtained an exact expression for the time-dependent
degree distribution, which has allowed us to elucidate network-wide properties
such as the existence of a connected component occupying a macroscopic portion
of the network. In addition, we have studied first-passage statistics, finding
that our system does not display a minimum in the mean first-passage time from
the origin to a target state $k^{*}>0$, typical of many processes with
resetting [1]. Instead, we find a monotonously increasing function that,
however, can present an inflection point depending on the growth and resetting
rates.

There are several directions for further research. One deals with the
inclusion of more realistic dynamical rules both in the network formation and
in the resetting process. This would partially change the mathematical form of
the master equations and the challenge would be related to solving more
complicated equations that could lead to new physical insights. On these
lines, it would be interesting to test the performance of different growth and
reset mechanisms in reproducing the evolution of empirical networks. Another
research avenue deals with a more complete characterization of the very same
model presented in this article. For example, when computing first-passage
statistics we could be interested in knowing the first-passage time of the
$i$th fastest node arriving to a target. The cluster dynamics would be very
interesting to be explored as well, as network growth with node removal can be
seen as a highly non-trivial aggregation-fragmentation process for which
writing a Smoluchowski-like equation and grasping analytical insights remain
as a considerable technical challenge. Finally, we acknowledge that correlated
transitions like the ones presented here are still an under-researched
characteristic in stochastic resetting models. Therefore, devising new
mathematically tractable models that incorporate this feature and that lend
themselves to be eventually compared to empirical data is a necessary step
toward a more complete, solid and useful theory of stochastic resetting.

## Acknowledgments

The author thanks M. De Domenico and D. Kiziridis for their comments on the
manuscript.

## Appendix: Continuous approach

To gain insight into the first-passage problem, we proceed by applying some
approximations, namely, we are treating the discrete degrees $k$ as a
continuous variable. We will see that even for small target degrees $k^{*}$,
hence a small effective interval, the approximation works well.

Our starting point are the master equations (17)–(19). Let $q(k,t)$ be
continuous version of $q_{k}(t)$, that is, $q(k,t)$ is a probability density
function that gives us the probability of finding a node with degree
$k\in\left[k,k+\mathop{}\\!\mathrm{d}k\right]$ at time
$t\in\left[t,t+\mathop{}\\!\mathrm{d}t\right]$. We notice that we can write

$\mathop{}\\!\frac{\partial q}{\partial t}=\Omega^{+}(k-\delta k)q(k-\delta
k,t)+\Omega^{-}(k+\delta k)q(k+\delta
k,t)-\left[\Omega^{+}(k)+\Omega^{-}(k)+\rho(k)\right]q(k,t),$ (A.1)

where $\Omega^{+}(k)=\alpha$ is the growing rate, while $\Omega^{-}(k)=rk$ and
$\rho(k)=r$ correspond to the resetting rates, the former due to a neighbor
deletion and the latter due to the deletion of the node itself. For
convenience, we leave $\delta k$ undetermined during the calculations and set
$\delta k=1$ at the end. Expanding both the rates and the pdf up to order
$(\delta k)^{2}$, we arrive at

$\mathop{}\\!\frac{\partial q}{\partial
t}=-\mathop{}\\!\frac{\partial}{\partial
k}\left[v(k)q(k,t)\right]+\mathop{}\\!\frac{\partial^{2}}{\partial
k^{2}}\left[D(k)q(k,t)\right]-\rho(k)q(k,t),$ (A.2)

where the drift and diffusion terms are, respectively, $v(k)=\delta
k(\Omega^{+}(k)-\Omega^{-}(k))$ and $D(k)=\delta
k^{2}(\Omega^{+}(k)+\Omega^{-}(k))/2$. We have not yet incorporated the
information given in master equations for the boundary degrees $k=0$ and
$k=k^{*}$. The former can be manually introduced as a delta-like source of
probability at the resetting point, $k_{\text{reset}}$,while the former takes
the form of an absorbing boundary condition. Putting all the pieces together,
we are dealt with

$\mathop{}\\!\frac{\partial q}{\partial t}=-\delta
k\mathop{}\\!\frac{\partial}{\partial k}\left[(\alpha-
rk)q(k,t)\right]+\frac{\delta
k^{2}}{2}\mathop{}\\!\frac{\partial^{2}}{\partial
k^{2}}\left[(rk+\alpha)q(k,t)\right]-rq(k,t)+r\delta(k-k_{\text{reset}}),$
(A.3)

with $q(k^{*},t)=0$, $\partial_{k}q(k,t)|_{k=0}=0$ and
$q(k,0)=\delta(k-k_{0})$.

In this continuous-degree approach, the survival probability reads
$S(k^{*},k_{0},t)=\int_{0}^{k^{*}}\,\mathop{}\\!\mathrm{d}k\,q(k,t)$, where
$k_{0}$ is the initial degree. Treating $k_{0}$ as a variable, the survival
probability satisfies the backward equation

$\mathop{}\\!\frac{\partial S}{\partial t}=(\alpha-
rk)\mathop{}\\!\frac{\partial S}{\partial k_{0}}+\frac{\delta
k^{2}}{2}(rk+\alpha)\mathop{}\\!\frac{\partial^{2}S}{\partial
k_{0}^{2}}-rS(k^{*},k_{0},t)+rS(k^{*},k_{\text{reset}},t).$ (A.4)

with $S(k^{*},k^{*},t)=0$ and $\partial_{k_{0}}S(k^{*},k_{0},t)|_{k_{0}=0}=0$.
Introducing the Laplace transform
$\tilde{S}(k^{*},k_{0},s)=\int_{0}^{\infty}\mathop{}\\!\mathrm{d}te^{-st}S(k^{*},k_{0},t)$,
we obtain

$(\alpha-rk)\mathop{}\\!\frac{\partial\tilde{S}}{\partial
k_{0}}+\frac{rk+\alpha}{2}\mathop{}\\!\frac{\partial^{2}\tilde{S}}{\partial
k_{0}^{2}}-(r+s)\tilde{S}(k^{*},k_{0},s)=-1-r\tilde{S}(k^{*},k_{\text{reset}},t).$
(A.5)

The general solution to this equation is

$\tilde{S}(k^{*},k_{0},s)=\frac{r\tilde{S}(k^{*},k_{\text{reset}},t)+1}{s+r}+c_{1}\mathcal{U}(1,k_{0})+c_{2}\mathcal{L}(0,k_{0}),$
(A.6)

where to ease the notation we have introduced the functions
$\mathcal{U}(x,y)=U(x+r/\alpha,x-1+4\alpha/r,2y+2\alpha/r)$, being
$U(\cdot,\cdot,\cdot)$ the confluent hypergeometric function of the second
kind, and $\mathcal{L}(x,y)=L_{-(x+1)-s/r}^{x-1+4\alpha/r}(2y+\alpha/r)$,
being $L_{a}^{b}(\cdot)$ the generalized Laguerre polynomial. The derivative
is

$\mathop{}\\!\frac{\partial\tilde{S}}{\partial
k_{0}}=-\frac{2}{r}(r+s)\,c_{1}\,\mathcal{U}(2,k_{0})-2\,c_{2}\,\mathcal{L}(1,k_{0}).$
(A.7)

Applying the boundary conditions, after some easy algebra we obtain

$\displaystyle c_{1}$
$\displaystyle=\frac{r\tilde{S}(k^{*},k_{\text{reset}},t)+1}{s+r}\frac{B_{2}}{A_{2}B_{1}-A_{1}B_{2}}$
(A.8) $\displaystyle c_{2}$
$\displaystyle=\frac{r\tilde{S}(k^{*},k_{\text{reset}},t)+1}{s+r}\frac{A_{2}}{A_{1}B_{2}-A_{2}B_{1}},$
(A.9)

with

$\displaystyle A_{1}$ $\displaystyle=\mathcal{U}(1,k^{*}),$ (A.10)
$\displaystyle B_{1}$ $\displaystyle=\mathcal{L}(0,k^{*}),$ (A.11)
$\displaystyle A_{2}$ $\displaystyle=-\frac{2}{r}(r+s)\mathcal{U}(2,0),$
(A.12) $\displaystyle B_{2}$ $\displaystyle=-2\mathcal{L}(1,0).$ (A.13)

Note that all these quantities are constant with respect to the “spatial”
variable $k_{0}$ but of course do depend explicitly on the growth and reset
rates, $\alpha$ and $r$, as well as the target degree $k^{*}$ and the time in
Laplace domain $s$. Now setting $k_{0}=k_{\text{reset}}=0$, we obtain the
Laplace transform of the probability that a randomly chosen node starts with
no connections and arrives for first time to reach degree $k^{*}$. Taking the
limit $s\to 0$, we have an analytical expression for the Laplace transform of
the survival probability, from which we can obtain the mean first-passage time
of a randomly chosen node of an initially empty network to reach degree
$k^{*}$.

Note that if we relax the condition $k_{\text{reset}}=0$, the value and
relative ordering of $k_{\text{reset}}$, $k^{*}$, and $k_{\text{ip}}$ (the
degree value at which the drift vanishes, see main text) will impact in the
behavior of the first-passage distribution and its moments. Despite being an
interesting theoretical exercise, these cases do not bring new
phenomenological insights with respect to the case $k_{\text{reset}}=0$. For
this reason, a careful report of all the combinations falls outside the scope
of the present article, and we stick to the case of nodes losing all their
links in the resetting events.

## References

## References

  * [1] Evans M R, Majumdar S N and Schehr G 2020 Journal of Physics A: Mathematical and Theoretical 53 193001
  * [2] Kusmierz L, Majumdar S N, Sabhapandit S and Schehr G 2014 Physical Review Letters 113 220602
  * [3] Rotbart T, Reuveni S and Urbakh M 2015 Physical Review E 92 060101
  * [4] Villen-Altamirano M, Villen-Altamirano J et al. 1991 Queueing, performance and Control in ATM 71–76
  * [5] Maurer S M and Huberman B A 2001 Journal of Economic Dynamics and Control 25 641–654
  * [6] Evans M R and Majumdar S N 2011 Physical Review Letters 106 160601
  * [7] Evans M R and Majumdar S N 2011 Journal of Physics A: Mathematical and Theoretical 44 435001
  * [8] Mukherjee B, Sengupta K and Majumdar S N 2018 Physical Review B 98 104309
  * [9] Basu U, Kundu A and Pal A 2019 Physical Review E 100 032136
  * [10] Grange P 2021 Journal of Physics A: Mathematical and Theoretical
  * [11] Sandev T, Domazetoski V, Kocarev L, Metzler R and Chechkin A 2022 Journal of Physics A: Mathematical and Theoretical 55 074003
  * [12] Evans M R and Majumdar S N 2014 Journal of Physics A: Mathematical and Theoretical 47 285001
  * [13] Riascos A P, Boyer D, Herringer P and Mateos J L 2020 Physical Review E 101 062147
  * [14] Bressloff P C 2021 Journal of Statistical Mechanics: Theory and Experiment 2021 063206
  * [15] Evans M R and Majumdar S N 2018 Journal of Physics A: Mathematical and Theoretical 51 475003
  * [16] Masoliver J 2019 Physical Review E 99 012121
  * [17] Gupta D 2019 Journal of Statistical Mechanics: Theory and Experiment 2019 033212
  * [18] Pal A and Reuveni S 2017 Physical Review Letters 118 030603
  * [19] Chechkin A and Sokolov I M 2018 Physical Review Letters 121 050601
  * [20] Nagar A and Gupta S 2016 Physical Review E 93 060102
  * [21] Kuśmierz Ł and Toyoizumi T 2019 Physical Review E 100 032110
  * [22] Masoliver J and Montero M 2019 Physical Review E 100 042103
  * [23] Gupta S, Majumdar S N and Schehr G 2014 Physical Review Letters 112 220601
  * [24] Durang X, Henkel M and Park H 2014 Journal of Physics A: Mathematical and Theoretical 47 045002
  * [25] Newman M E 2003 SIAM Review 45 167–256
  * [26] Broder A, Kumar R, Maghoul F, Raghavan P, Rajagopalan S, Stata R, Tomkins A and Wiener J 2000 Computer Networks 33 309–320
  * [27] Colizza V, Flammini A, Serrano M A and Vespignani A 2006 Nature Physics 2 110–115
  * [28] Peixoto T P 2014 Physical Review X 4 011047
  * [29] Noldus R and Van Mieghem P 2015 Journal of Complex Networks 3 507–542
  * [30] Broido A D and Clauset A 2019 Nature Communications 10 1–10
  * [31] Voitalov I, van der Hoorn P, van der Hofstad R and Krioukov D 2019 Physical Review Research 1 033034
  * [32] De Domenico M, Solé-Ribalta A, Cozzo E, Kivelä M, Moreno Y, Porter M A, Gómez S and Arenas A 2013 Physical Review X 3 041022
  * [33] Artime O, Benigni B, Bertagnolli G, d’Andrea V, Gallotti R, Ghavasieh A, Raimondo S and De Domenico M 2022 Multilayer Network Science: From Cells to Societies Elements in Structure and Dynamics of Complex Networks (Cambridge University Press)
  * [34] Battiston F, Cencetti G, Iacopini I, Latora V, Lucas M, Patania A, Young J G and Petri G 2020 Physics Reports 874 1–92
  * [35] Arenas A, Díaz-Guilera A, Kurths J, Moreno Y and Zhou C 2008 Physics Reports 469 93–153
  * [36] Pastor-Satorras R, Castellano C, Van Mieghem P and Vespignani A 2015 Reviews of Modern Physics 87 925
  * [37] Porter M A and Gleeson J P 2016 Frontiers in Applied Dynamical Systems 4
  * [38] Artime O, Ramasco J J and San Miguel M 2017 Scientific Reports 7 1–10
  * [39] Hens C, Harush U, Haber S, Cohen R and Barzel B 2019 Nature Physics 15 403–412
  * [40] Artime O and De Domenico M 2020 New Journal of Physics 22 093035
  * [41] Watts D J and Strogatz S H 1998 Nature 393 440–442
  * [42] Barabási A L and Albert R 1999 Science 286 509–512
  * [43] Dorogovtsev S N, Mendes J F F and Samukhin A N 2000 Physical Review Letters 85 4633
  * [44] Dorogovtsev S N and Mendes J F F 2001 Physical Review E 63 025101
  * [45] Krapivsky P L, Redner S and Leyvraz F 2000 Physical Review Letters 85 4629
  * [46] Bianconi G and Barabási A L 2001 Physical Review Letters 86 5632
  * [47] Dorogovtsev S N and Mendes J F F 2000 EPL (Europhysics Letters) 52 33
  * [48] Sarshar N and Roychowdhury V 2004 Physical Review E 69 026101
  * [49] Moore C, Ghoshal G and Newman M E 2006 Physical Review E 74 036121
  * [50] Srinivasan U T, Dunne J A, Harte J and Martinez N D 2007 Ecology 88 671–682
  * [51] Saavedra S, Reed-Tsochas F and Uzzi B 2008 Proceedings of the National Academy of Sciences 105 16466–16471
  * [52] Chung F and Lu L 2004 Internet Mathematics 1 409–461
  * [53] Bauke H, Moore C, Rouquier J B and Sherrington D 2011 The European Physical Journal B 83 519–524
  * [54] Kalay Z and Ben-Naim E 2014 Journal of Physics A: Mathematical and Theoretical 48 045001
  * [55] Crane E, Freeman N and Tóth B 2015 Electronic Journal of Probability 20 1–33
  * [56] Zhang X, He Z and Rayman-Bacchus L 2016 Journal of Statistical Physics 162 842–854
  * [57] Pelizzola A, Pretti M and Zamparo M 2021 Europhysics Letters 133 60003
  * [58] Miron A and Reuveni S 2021 Physical Review Research 3 L012023
  * [59] Miron A 2022 arXiv preprint arXiv:2207.07114
  * [60] Bertin E 2022 Journal of Physics A: Mathematical and Theoretical 55 384007
  * [61] Krapivsky P L, Vilk O and Meerson B 2022 Physical Review E 106(3) 034125
  * [62] Krapivsky P L, Redner S and Ben-Naim E 2010 A kinetic view of Statistical Physics (Cambridge University Press)
  * [63] Molloy M and Reed B 1995 Random Structures & Algorithms 6 161–180