# Self-supervised models of speech infer
universal articulatory kinematics

###### Abstract

Self-Supervised Learning (SSL) based models of speech have shown remarkable
performance on a range of downstream tasks. These state-of-the-art models have
remained blackboxes, but many recent studies have begun “probing” models like
HuBERT, to correlate their internal representations to different aspects of
speech. In this paper, we show “inference of articulatory kinematics” as
fundamental property of SSL models, i.e., the ability of these models to
transform acoustics into the causal articulatory dynamics underlying the
speech signal. We also show that this abstraction is largely overlapping
across the language of the data used to train the model, with preference to
the language with similar phonological system. Furthermore, we show that with
simple affine transformations, Acoustic-to-Articulatory inversion (AAI) is
transferrable across speakers, even across genders, languages, and dialects,
showing the generalizability of this property. Together, these results shed
new light on the internals of SSL models that are critical to their superior
performance, and open up new avenues into language-agnostic universal models
for speech engineering, that are interpretable and grounded in speech science.

Index Terms—  Self-Supervised Learning; Articulatory Kinematics; Cross-lingual
and Multilingual Speech Processing;

##  1 Introduction

Self-supervised learning (SSL) has revolutionized every field of machine
learning, by providing rich features of natural data without human-annotated
labels. Likewise, speech SSL models have been proven to be successful in
various speech downstream tasks [1]. To understand such utility, the internal
representation of speech SSL models has been scrutinized by probing analyses
for known speech and linguistic features, such as low-level acoustics,
phonetics, and lexical semantics [2, 3, 4, 5]. A comparative analysis by Cho
et al. (2023) [5] demonstrates that the state-of-the-art SSL models are highly
correlated with articulatory kinematics and the correlation score can indicate
the success of the SSL model in downstream tasks. This finding is extended to
developing a high-performance Acoustic-to-Articulatory inversion (AAI) model
[6].

Here, we test an intriguing hypothesis – speech SSL models infer the causal
articulatory processes that generate the speech acoustic signal. If this
hypothesis is true, it should be agnostic to language or dialect, as speech
acoustics are a result of the resonance associated with vocal tract shapes
that add spectral detail to the vocal source, i.e., vibration of the vocal
folds. The human species has a common canonical structure of vocal tract
anatomy and orofacial muscles, regardless of ethnic group, language, and
dialect.

We validate this hypothesis using two approaches: cross-lingual articulatory
probing and cross-speaker transferability analysis. First, we probe
articulatory representation in speech SSL models trained with different
languages, using a comprehensive set of public electromagnetic articulography
(EMA) datasets that cover up to 62 speakers from different language and
dialect groups. Then, we test how well each individual articulatory system can
be transferred to another system, by fitting an affine transformation.
Combining these two analyses elucidates the universality of the articulatory
representation in speech SSL models.

While EMA has been a major interface of studying articulatory phonology [7, 8,
9, 10], the data collection procedure is highly complicated with many physical
restrictions that have posed significant barriers in collecting large scale
dataset [11]. Moreover, a majority of the publically available datasets are
concentrated in a few common languages like English.111This is a significant
bottleneck given there are more than 7,000 languages around the world.
Therefore, the development of universal AAI models can reduce the burden of
collecting new EMA data for investigating diverse languages. This is
especially valuable for languages with low resources or near extinction where
only a short amount of audio is accessible. Our results suggest that the
articulatory representation in speech SSL models is transferable across
languages, making a critical step toward the universal AAI model.

Our major findings are:

  * •

Articulatory kinematics of different languages and dialects can be recovered
from speech SSL models with a simple linear projection, regardless of the
language for which SSL models are trained.

  * •

Individual articulatory systems can be aligned by affine transformations,
implying the existence of canonical basis of articulatory kinematics.

  * •

While the overall scores are surprisingly high, we find a preference for
language or dialect, which is evidence of language-specific articulatory
phonology.

##  2 Related Work

Speech SSL. Wav2vec-2.0 (W2V2) [12] and HuBERT (HB) [13] have been the most
successful SSL models of speech. Both models are trained by masked prediction
objective, where the Transformer encoder predicts randomly masked parts of
input features. As the raw audio waves are inadequate for target reference,
instead, the models predict self-driven units, either from vector quantization
(W2V2) or online clustering (H2B) of internal features.
Cross-lingual speech model. XLS-R [14] and MMS [15] are direct extensions of
W2V2, which aim to learn cross-lingual speech representation by training on
multiple languages. XLS-R is trained on 50K hours of speech from 53 languages,
and MMS is trained on 55K hours of speech from more than 1,000 languages. As
the models are trained across languages, they may learn articulatory
representation that is shared across languages.
Probing SSL representation. Several efforts have been made to explain which
features are encoded by SSL. Pasad et al. [2, 4] compare lexical semantics
with speech SSL models across layers, revealing that lexical semantics rise in
the later layer of SSL models. Shah et al. [3] train non-linear models to
probe audio, fluency, pronunciation, and text features, providing a detailed
description of information processing in SSL models. Furthermore, clustering
analysis suggests that the primary information encoded in SSL models is fine-
grained subphonemic units [16, 17]. Our work is a direct extension of Cho et
al. [5] that demonstrates evidence of high-fidelity articulatory kinematics in
SSL models.
AAI. Inversion models often suffer from individual differences in EMA data. A
major source of such variability is anatomical differences across speakers.
Furthermore, though key articulators are standardized, there is some level of
arbitrariness in the locations of the sensors. Together, such misalignment
across subjects has challenged developing a speaker independent AAI model. The
recent work by Wu et al. [6], mitigates the problem by utilizing relative
distances within articulators, which are more consistent across speakers than
EMA. This is well aligned with our finding that affine transformations can
align individual articulatory systems.

##  3 Methods

###  3.1 Speech SSL models trained for different languages

We use HB and W2V2 with large sizes (300M parameters) to extract speech
representations. We retrieved either of those models from Hugginface trained
with a specific language: English, Mandarin, Korean, French, or
Dutch.Furthermore, two cross-lingual W2V2 models are included: XLS-R [14] with
300M parameters trained on 53 languages, and MSS [15] with 1B parameters
trained on more than 1000 languages. For baseline reference, a set of acoustic
features (filter bank, Melspectrogram, and MFCC) are included as well.

###  3.2 Comprehensive set of multi-lingual, multi-dialect EMA

Our comprehensive analyses include five public EMA datasets that cover 62
speakers from three languages: English, Mandarin, and Italian (Table 1). For
English speakers, there are three regional dialect groups where speakers are
from the UK, USA, or China (Beijing or Shanghai). As the EMA data can be
highly noisy due to the known instability data collection procedure, we filter
out three Mandarin speakers from the original corpus. Note that the speakers
validated in the original dataset papers are included in the final dataset.
Each dataset has a frame-wise annotation of kinematic traces of six
articulators: low incisor (LI), upper lip (UL), lower lip(LL), tongue tip
(TT), tongue blade (TB), and tongue dorsum (TD)), on the midsagittal plane: X
(front-back) and Y (up-down). Each trace is normalized to be zero mean and
unit variance.222The normalization is done per channel and within the clip.
All of these data were collected while naturally speaking full sentences.

###  3.3 Probing SSL models with articulatory kinematics

We follow the previous articulatory probing approach [5] to measure the
correlation between SSL representation and articulatory kinematics. The
features from the 20 ms audio frame are extracted from each layer of each SSL
model and projected to the EMA space by fitting a linear inversion model. One
difference from the previous study is that a 6 Hz low-pass filter is applied
to remove high-frequency noise. The prediction performance is measured by the
correlation of prediction and ground truth reference and averaged over 5-fold
cross-validation.

###  3.4 Evaluating transferability between speakers

We hypothesize that the linear inversion models obtained from the probing
analysis are leveraging the same basis of the articulatory subspace that
resides in SSL representation. To test this hypothesis, we devise a
transferability metric, which measures the correlation between two
articulatory systems. Suppose we have two AAI models $f_{A}$ and $f_{B}$, for
different speakers, A and B. As there is variability in the vocal tract
anatomy across individuals, we trained another linear model, $g_{A\rightarrow
B}$, to align the articulatory systems, which is trained with L1
regularization weighted with $\alpha$ = 0.01 to impose sparsity in
coefficients. The transferability from speaker A to speaker B is then measured
as the correlation of the prediction by these two aligned inversion systems,
$\text{corr}(g_{A\rightarrow B}\circ f_{A},f_{B})$. For this analysis, linear
inversion models built upon the 17th layer of XLS-R are used and the data from
the target speaker (B) is used for training and testing. We only include
speakers with high inversion performance (corr $\geq 0.8$). We denote the A
and B speaker pair as the source and target pair.

Table 1: Summary of the public EMA datasets used in the analyses. Different language-dialect groups are denoted as EN.{UK, US, BJ, SH}: English speakers from the UK; the US; Beijing, or Shanghai, China, respectively. MAN: Mandarin speakers. IT: Italian speakers. The number of speakers used for the transferability analysis is denoted in ($\cdot$) if different. The gender distribution is balanced in MOCHA-TIMIT, HPRC, and EMA-MAE.  Corpus | Lang. | Spk # | m/spk
---|---|---|---
MNGU0 [18] | EN.UK | 1 | 75
MOCHA-TIMIT [19] | EN.UK | 7 | 27
HPRC [20] | EN.US | 8 | 59
EMA-MAE [21] |  | EN.US
---
EN.BJ
EN.SH
| 20 (18)
---
10 (9)
9 (5)
| 12
---
17
16
DKU-JNU-EMA [22] | MAN | 4 (2) | 20
MSPKA [23] | IT | 3 (2) | 47
Total | – | 62 (52) | –

##  4 Results

Fig. 1: (left) EMA prediction performance of probing SSL models from different
languages. The correlations are averaged across 12 EMA channels. Each
language-dialect group is denoted by a hatch pattern and each dot denotes an
individual speaker. Regardless of the language, the average performance
reaches over 0.8. (right) Performance comparison of English SSL versus
Mandarin SSL, each panel shows specific language-dialect groups. The diagonal
dashed lines denote identity lines. For English, native speakers (EN.UK/US)
prefer English models over Mandarin models, scattered slightly above the
diagonals, but speakers from China (EN.BJ/SH) show almost identical scores.

###  4.1 Articulatory probing shows high correlations agnostic to training
language of SSL models

Fig.1 shows the prediction performance of probing SSL models in §3.1 that are
trained in different languages. For specific SSL models, the performance of
the layer with the highest correlation is reported.333The layer-wise patterns
are omitted due to the limited space, but the patterns are consistent with
[5]. For some feature sets with multiple models, the models with higher scores
are selected.444English and Mandarin have both HB and W2V2. Despite the
substantial discrepancy between languages, the average correlations are
surprisingly high, reaching over 0.8 regardless of the training language
(Fig.1 (left)), which are far beyond the baseline acoustic features (corr =
0.66$\pm$0.040). This indicates that the SSL models naturally learn physical
relationships between raw speech signals and vocal tract articulations, which
are not bound to specific language, but rather universal properties of the
human vocal tract system. Amongst languages, the highest performance is
achievable by English SSL models (corr = 0.835$\pm$0.039) and cross-lingual
models (XLS-R: 0.830$\pm$0.039, MMS: 0.832$\pm$0.040). The performance
difference between MMS and the other two best models is not significant. Given
the scale of MMS in both training data and model size, this finding suggests
that scaling up has a minimal effect, and learning articulatory representation
is readily saturated in the SSL regime.

Probing models for individual speakers show variable performances due to the
noisy nature of the data collection, which is mostly induced by the unstable
attachment of sensors. Still, 44% of speakers show high performance above 0.85
correlations, and the best subject achieves 0.89 correlation.555One of the
Italian speakers using the 10th layer of HB trained on English (among 24
layers). This is remarkable given that we are applying a very simple linear
regression on the frozen SSL models that are never exposed to EMA data while
training.

###  4.2 Preference for language in articulatory probing

There is some noteworthy variability across languages. English models and
Mandarin models show higher correlations than models of other languages. As
most of the speakers are from either the US, UK, or China, this may indicate
some language specificity. However, this difference can not be fully
attributed to the language difference because the probing performance is also
correlated with in general representational quality of SSL model [5].

To further test language specificity in SSL models, we compare English SSL
models and Mandarin SSL models for each dialect group in English speakers
(Fig.1 (right)). The native English speakers (EN.UK/US) prefer English models
over Mandarin models (Fig.1 (right) top two panels), to a small extent yet
statistically significant (mean diff = 0.019; $p$ $<$ $1e-5$). However, the
English speakers from China (EN.BJ/SH) show almost identical probing
performances of these models (Fig.1 (right) bottom left panel), showing non
significant difference (mean diff = -0.002; $p$ = 0.29). As the language is
controlled to be English, the observed preference is attributed to differences
in articulatory phonology rather than other linguistic components (e.g.,
semantics or syntax). This indicates that there are some residual articulatory
habits originating from speaking Mandarin as a mother tongue language. For the
other language groups (MAN and IT), show relatively minor differences and the
number of speakers is too few to evaluate the significance.

###  4.3 Individual articulatory systems are mutually transferable with affine
transformations

Fig. 2: (far-left) Correlation matrix denoting transferabilities between
language-dialect groups. The scores are averaged over possible pairs between
groups. (mid-left) Distributions of correlations across dialects (blue) and
within dialects (orange). The English speakers from China (EN.SH+EN.BJ) and
the native English speakers from the US (EN.US) groups in the EMA-MAE dataset
are used. (mid-right) Distributions of correlations across gender (blue) and
within gender (orange). Four male and four female speakers from the HPRC
dataset are used. (far-right) The average absolute coefficients in the affine
transformations between the articulatory systems.

For each pair of speakers, we measure the transferability scores (§3.4)
between the articulatory systems obtained from the probing analyses. The far-
left panel of Fig.2 shows the average correlations between source and target
group pairs. The scores within groups are strikingly high, showing near or
above 0.85 average correlation (diagonal brown boxes). This suggests that the
articulatory systems of individual speakers are affine transformations of each
other. Such high correlations are also observable in the transformations
across language-dialect groups, especially between the Italian (IT) and the
native English speakers EN.UK/US which shows even higher scores than within-
group scores. The IT has in general high transferability to others, indicating
that the AAI models built upon one of these particular speakers can serve as
the closest universal AAI model. These high-fidelity affine transformations
between different speakers demonstrate a significant level of isometry in the
human articulatory system. This is strong evidence that speech SSL models
infer the canonical basis of the articulatory kinematics.

The coefficients in the affine transformations are mostly assigned to the same
articulators (red boxes in Fig.2 (far-right)). This suggests that the
transformations are geometric alignments that correct the variability induced
by the difference in vocal tract anatomy or the sensor locations. While most
of the six articulators show definite self-assignments, the tongue blade (TB)
is affected by other tongue parts, which is natural since TB bridges the
tongue tip and dorsum. Fig.3 shows transferability scores for each
articulator. Some channels, LIX, LIX, ULY, and LLX, are relatively difficult
to align, which reflects their innate subtlety in shaping phonetic
information.

###  4.4 The variance in transferability reflects language-specific
articulatory phonology

When we compare the within and across dialects transferability, the across
dialect transformations show significantly lower scores than those of the
within dialect cases (Fig.2 (mid-left)). However, when groups are divided by
gender, there is no significance between the within and across gender
transformations (Fig.2 (mid-right)). As gender is the most significant factor
of anatomical differences, this result suggests that the transferability
between articulatory systems is not affected by the vocal tract anatomy.
Therefore, the observed patterns across languages and dialects are more likely
due to differences in articulatory phonology. This is consistent with the
language preference in cross-lingual probing §4.2. To avoid variance induced
by the data collection site, we confined analyses on dialect groups and gender
groups to the EMA-MAE corpus and HPRC corpus, respectively.

Fig. 3: Transferability correlations of each articulator averaged over all
source-target pairs.

This language-specific pattern is also reflected in the transferability matrix
(Fig.2 (far-left) black box). While the transferability correlations from the
English groups (EN.*) to the Mandarin group (MAN) are generally lower than in
other cases, the Chinese English groups (EN.BJ/SH) are more transferable than
the native English speakers (EN.UK/US). Such a pattern is not observed in the
opposite cases (MAN $\rightarrow$ EN.*).

##  5 Discussion

We demonstrate that the speech SSL models are able to recover articulatory
kinematics with simple linear mapping, achieving high performance inversion
regardless of speakers, languages, and dialects. This is further verified by
finding affine transformations from one articulatory system to another. Our
findings provide strong evidence that there is a canonical basis of
articulatory phonology which is naturally emerging in self-supervised learning
of speech.

Our findings evoke another interesting hypothesis that articulatory kinematics
are not only the physical interface of speech but also the continuous
embedding representations of phonetics. As the SSL models are trained by the
masked prediction objective without any external labels, the resulting
representations are perceptual descriptions of how sounds are shaped in
natural speech data, which also shows high correspondence to the human
auditory perception [24]. Here, we demonstrate that the natural selection of
the perceptual embedding space of speech is a continuous, dynamic articulatory
space. Particularly, the group-level analyses in §4.4 reveal the language
specificity in articulatory systems but little evidence of gender specificity
that entails large anatomical differences. This observation further supports
the hypothesis of equivalence between articulatory kinematics and continuous
phonetic embeddings.

The articulatory information captured by EMA is not complete as limited to the
midsagittal plane, and upper bounded by the inevitable noise ceiling of the
sensors. Also, the source information such as vocal fold vibration is missing.
However, our findings indicate that a complete set of articulatory kinematics
may lie in SSL representations which are missing in EMA. This is an
interesting future work toward the unsupervised discovery of articulatory
kinematics without relying on EMA.

##  6 Acknowledgements

This research is supported by the following grants to PI Anumanchipalli — NSF
award 2106928, BAIR Commons-Meta AI Research, the Rose Hills Innovator
Program, and UC Noyce Initiative, at UC Berkeley. Special thanks to Shang-Wen
(Daniel) Li for discussions and for providing advice.

## References

  * [1] Shu-wen Yang, Po-Han Chi, Yung-Sung Chuang, Cheng-I Jeff Lai, Kushal Lakhotia, Yist Y Lin, Andy T Liu, Jiatong Shi, Xuankai Chang, Guan-Ting Lin, et al.,  “Superb: Speech processing universal performance benchmark,”  arXiv preprint arXiv:2105.01051, 2021.

  * [2] Ankita Pasad, Ju-Chieh Chou, and Karen Livescu,  “Layer-wise analysis of a self-supervised speech representation model,”  in 2021 IEEE Automatic Speech Recognition and Understanding Workshop (ASRU), 2021, pp. 914–921.
  * [3] Jui Shah, Yaman Kumar Singla, Changyou Chen, and Rajiv Ratn Shah,  “What all do audio transformer models hear? probing acoustic representations for language delivery and its structure,”  arXiv preprint arXiv:2101.00387, 2021.
  * [4] Ankita Pasad, Bowen Shi, and Karen Livescu,  “Comparative layer-wise analysis of self-supervised speech models,”  in ICASSP 2023-2023 IEEE International Conference on Acoustics, Speech and Signal Processing (ICASSP). IEEE, 2023, pp. 1–5.
  * [5] Cheol Jun Cho, Peter Wu, Abdelrahman Mohamed, and Gopala K Anumanchipalli,  “Evidence of vocal tract articulation in self-supervised learning of speech,”  in ICASSP 2023-2023 IEEE International Conference on Acoustics, Speech and Signal Processing (ICASSP). IEEE, 2023, pp. 1–5.
  * [6] Perter Wu, Li-wei Chen, Cheol Jun Cho, Shinji Watanabe, Louis Goldstein, Alan W Black, and Gopala K Anumanchipalli,  “Speaker-independent acoustic-to-articulatory speech inversion,”  in ICASSP 2023-2023 IEEE International Conference on Acoustics, Speech and Signal Processing (ICASSP). IEEE, 2023, pp. 1–5.
  * [7] Catherine P Browman and Louis Goldstein,  “Articulatory phonology: An overview,”  Phonetica, vol. 49, no. 3-4, pp. 155–180, 1992.
  * [8] Vikram Ramanarayanan, Louis Goldstein, and Shrikanth S Narayanan,  “Spatio-temporal articulatory movement primitives during speech production: Extraction, interpretation, and validation,”  The Journal of the Acoustical Society of America, vol. 134, no. 2, pp. 1378–1394, 2013.
  * [9] Jiachen Lian, Alan W Black, Louis Goldstein, and Gopala Krishna Anumanchipalli,  “Deep neural convolutive matrix factorization for articulatory representation decomposition,”  arXiv preprint arXiv:2204.00465, 2022.
  * [10] Peter Wu, Shinji Watanabe, Louis Goldstein, Alan W Black, and Gopala Krishna Anumanchipalli,  “Deep Speech Synthesis from Articulatory Representations,”  in Proc. Interspeech 2022, 2022, pp. 779–783.
  * [11] Teja Rebernik, Jidde Jacobi, Roel Jonkers, Aude Noiray, and Martijn Wieling,  “A review of data collection practices using electromagnetic articulography,”  Laboratory Phonology, vol. 12, no. 1, pp. 6, 2021.
  * [12] Alexei Baevski, Yuhao Zhou, Abdelrahman Mohamed, and Michael Auli,  “wav2vec 2.0: A framework for self-supervised learning of speech representations,”  Advances in neural information processing systems, vol. 33, pp. 12449–12460, 2020.
  * [13] Wei-Ning Hsu, Benjamin Bolte, Yao-Hung Hubert Tsai, Kushal Lakhotia, Ruslan Salakhutdinov, and Abdelrahman Mohamed,  “Hubert: Self-supervised speech representation learning by masked prediction of hidden units,”  IEEE/ACM Transactions on Audio, Speech, and Language Processing, vol. 29, pp. 3451–3460, 2021.
  * [14] Arun Babu, Changhan Wang, Andros Tjandra, Kushal Lakhotia, Qiantong Xu, Naman Goyal, Kritika Singh, Patrick von Platen, Yatharth Saraf, Juan Pino, et al.,  “Xls-r: Self-supervised cross-lingual speech representation learning at scale,”  arXiv preprint arXiv:2111.09296, 2021.
  * [15] Vineel Pratap, Andros Tjandra, Bowen Shi, Paden Tomasello, Arun Babu, Sayani Kundu, Ali Elkahky, Zhaoheng Ni, Apoorv Vyas, Maryam Fazel-Zarandi, et al.,  “Scaling speech technology to 1,000+ languages,”  arXiv preprint arXiv:2305.13516, 2023.
  * [16] Amitay Sicherman and Yossi Adi,  “Analysing discrete self supervised speech representation for spoken language modeling,”  in ICASSP 2023-2023 IEEE International Conference on Acoustics, Speech and Signal Processing (ICASSP). IEEE, 2023, pp. 1–5.
  * [17] Badr M Abdullah, Mohammed Maqsood Shaik, Bernd Möbius, and Dietrich Klakow,  “An information-theoretic analysis of self-supervised discrete representations of speech,”  arXiv preprint arXiv:2306.02405, 2023.
  * [18] Korin Richmond, Phil Hoole, and Simon King,  “Announcing the electromagnetic articulography (day 1) subset of the mngu0 articulatory corpus,”  in Twelfth Annual Conference of the International Speech Communication Association, 2011.
  * [19] Alan Wrench,  “Mocha: multichannel articulatory database,” http://www.cstr.ed.ac.uk/research/project/
artic/mocha.html, 1999.

  * [20] Mark Tiede, Carol Y Espy-Wilson, Dolly Goldenberg, Vikramjit Mitra, Hosung Nam, and Ganesh Sivaraman,  “Quantifying kinematic aspects of reduction in a contrasting rate production task,”  The Journal of the Acoustical Society of America, vol. 141, no. 5, pp. 3580–3580, 2017.
  * [21] An Ji, Jeffrey J Berry, and Michael T Johnson,  “The electromagnetic articulography mandarin accented english (ema-mae) corpus of acoustic and 3d articulatory kinematic data,”  in 2014 IEEE International Conference on Acoustics, Speech and Signal Processing (ICASSP). IEEE, 2014, pp. 7719–7723.
  * [22] Zexin Cai, Xiaoyi Qin, Danwei Cai, Ming Li, Xinzhong Liu, and Haibin Zhong,  “The dku-jnu-ema electromagnetic articulography database on mandarin and chinese dialects with tandem feature based acoustic-to-articulatory inversion,”  in 2018 11th International Symposium on Chinese Spoken Language Processing (ISCSLP). IEEE, 2018, pp. 235–239.
  * [23] Claudia Canevari, Leonardo Badino, and Luciano Fadiga,  “A new italian dataset of parallel acoustic and articulatory data,”  in Sixteenth Annual Conference of the International Speech Communication Association, 2015.
  * [24] Yuanning Li, Gopala K Anumanchipalli, Abdelrahman Mohamed, Junfeng Lu, Jinsong Wu, and Edward F Chang,  “Dissecting neural computations of the human auditory pathway using deep neural networks for speech,”  bioRxiv, pp. 2022–03, 2022.


Propagating degrees of freedom
in f(R) gravity

Yun Soo Myunga

Institute of Basic Science and Department of Computer Simulation,

Inje University, Gimhae 621-749, Korea

Abstract

We have computed the number of polarization modes of gravitational waves
propagating in the Minkowski background in $f(R)$ gravity. This is three of
two from transverse-traceless tensor modes and one from a massive trace mode,
which confirms the results found in the literature. There is no massless
breathing mode and the massive trace mode corresponds to the Ricci scalar. A
newly defined metric tensor in $f(R)$ gravity satisfies the transverse-
traceless (TT) condition as well as the TT wave equation.

PACS numbers: 04.30.-w, 04.50.Kd

<EMAIL_ADDRESS>

##  1 Introduction

The $f(R)$ gravity theory is considered as a representative theory of modified
gravities. The $f(R)$ gravity [1, 2, 3, 4] has much attentions as a strong
candidate for explaining the current accelerating universe [5, 6]. When
choosing the Hu-Sawicki model [7], the theory could give rise to the late time
cosmic acceleration without violating the gravity tests in the solar system
and without affecting high redshift physics. Very recently, the observational
constraint on this model were reported from weal lensing peak abundances [8].
Particularly, $f(R)=R+R^{2}/(6M^{2})$ gravity [9, 10, 11] has shown a strong
evidence for inflation to support recent Planck data [12]. An important
feature of this model indicates that the inflationary dynamics were driven by
the purely gravitational interaction $R^{2}$ and the scale of inflation is
linked to the mass parameter $M^{2}$. This theory could thus provide a unified
picture of both inflation in the early universe and the accelerated expansion
at later times. In addition, black hole [13, 14, 15] and traversal wormhole
solutions [16, 17] have been found within $f(R)$ gravity in recent years. The
recent detection of gravitational waves by the LIGO Collaboration [18] is
surely a milestone in gravitational waves research and opens new perspectives
in the study of Einstein gravity (general relativity) and astrophysics. Hence,
it is meaningful to explore gravitational waves in the modified theory of
gravity, especially in $f(R)$ gravity. The observation of the polarization
modes of gravitational waves will be a crucial tool to obtain valuable
information about the black holes and the physics of the early universe.

It is well-known that the Einstein gravity with two polarization degrees of
freedom (DOF) is distinguished from the metric $f(R)$ gravity with three DOF
[19]. Importantly, it is worth noting that the Einstein equation derived from
$f(R)$ gravity contains fourth-order derivative terms. A simple way to avoid a
difficulty dealing with the fourth-order equation is to transform the $f(R)$
gravity into a scalar-tensor theory which is surely a second-order theory.

Very recently, it was reported that a polynomial $f(R)$ model could provide
two additional scalars of a massive longitudinal mode (perturbed Ricci scalar:
$R^{(1)}$) and a massless transverse mode (breathing mode: $\hat{h}_{b}$), in
addition to the two TT tensor modes ($\hat{h}^{+},\hat{h}^{\times}$) [20]. A
breathing mode seems to be overlooked in the literature because of the
assumption that the application of the Lorentz gauge implies the TT wave
equation. Also, it was insisted that four DOF found in [20] is consistent with
the result obtained from the Newman-Penrose (NP) formalism. However, the
presence of a breathing mode contradicts to the well-known fact in the
literature that the $f(R)$ gravity involves three DOF of a massive
longitudinal mode and two spin-2 modes. Hereafter, we wish to call this as the
issue of DOF in $f(R)$ theories.

In this work, we wish to point out that the $f(R)$ gravity still involves
three DOF by investigating the fourth-order equation composed of a second-
order tensor and a fourth-order scalar.

It seems that there is no breathing mode because the perturbed Ricci scalar
$R^{(1)}$ is related closely to the trace ‘$h$’ of perturbed metric tensor.
Hence, the allocation of the Ricci scalar as a newly scalar represents the the
trace of metric tensor. Also, we wish to remind the reader that the Ricci
scalar equation is not an independent equation and is not separated from the
perturbed Einstein equation because it comes out just from taking the trace of
the latter equation. This implies that the Ricci scalar is an emergent mode
from $h$. Furthermore, it is instructive to note that in the TT gauge, there
is a close connection between the metric perturbation and the linearized
Riemann tensor, implying that $\delta R_{itjt}=-\ddot{h}^{\rm TT}_{ij}/2$
[21]. This gauge is very convenient because it fixes all local gauge freedoms.
But, it might be unclear that there exists a close relation between the metric
perturbation and the NP formalism unless one chooses the TT gauge. One could
not naively choose the TT gauge in the perturbed $f(R)$ gravity because of
$h\not=0$, whereas the Lorentz gauge is easily implemented to eliminate the
gauge DOF. However, one might choose the TT gauge to obtain a massless spin-2
in the perturbed $f(R)$ gravity when one introduces a newly metric
perturbation $\tilde{h}_{\mu\nu}$.

The organization of our work is as follows. In the section 2, we briefly
describe the $f(R)$ gravity and its scalar-tensor theory and derive two sets
of perturbed equations around the Minkowski background in the section 3. Sec.4
is focused on obtaining the number of propagating DOF when one chooses the
Lorentz gauge. Finally, we will discuss our result which shows that there is
no breathing mode in the section 5.

##  2 $f(R)$ gravity and its scalar-tensor theory

Instead of a polynomial model of $f(R)=R+\alpha R^{2}+\beta R^{3}+\cdots$ [22,
20], we start with a specific $f(R)$ gravity (Starobinsky model [9])

$\displaystyle S_{\rm f}=\frac{M^{2}_{\rm P}}{2}\int
d^{4}x\sqrt{-g}f(R),\leavevmode\nobreak\ \leavevmode\nobreak\
f(R)=R+\frac{R^{2}}{6M^{2}},$ (1)

where the $R^{2}$-term was originally motivated by one-loop correction to
Einstein gravity. Here the mass parameter $M^{2}$ is chosen to be a positive
value, which is consistent with the stability condition of
$f^{\prime\prime}(0)>0$ [2]. This model of $f(R)$ gravity is enough to find
the propagating DOF around the Minkowski background. The Einstein equation
takes the form

$\displaystyle
R_{\mu\nu}f^{\prime}(R)-\frac{1}{2}g_{\mu\nu}f(R)+\Big{(}g_{\mu\nu}\nabla^{2}-\nabla_{\mu}\nabla_{\nu}\Big{)}f^{\prime}(R)=0,$
(2)

where the prime (′) denotes the differentiation with respect to its argument.

On the other hand, one might represent (1) as a scalar-tensor theory by
introducing an auxiliary field $\psi$ [10]

$S_{\rm A}=\int d^{4}x\sqrt{-g^{J}}\Big{[}\frac{M^{2}_{\rm
P}}{2}R+\frac{M_{\rm P}}{M}R\psi-3\psi^{2}\Big{]},$ (3)

where the superscript $J$ means the Jordan frame. Varying $S_{\rm A}$ with
respect to $\psi$ provides

$\psi=\frac{M_{\rm P}}{6M}R$ (4)

which means that the Ricci scalar is treated as an independent scalar degree
of freedom. Plugging (4) into $S_{\rm A}$ again leads to the original $f(R)$
gravity (1) exactly.

Making use of the conformal transformation and redefining the scalar field
($\psi\to\phi$)

$g^{J}_{\mu\nu}\to\frac{1}{1+\frac{2\psi}{MM_{\rm P}}}g^{E}_{\mu\nu}\to
e^{-\sqrt{\frac{2}{3}}\frac{\phi}{M_{\rm P}}}g^{E}_{\mu\nu},$ (5)

one arrives at the Starobinsky model in the Einstein frame [11]

$S_{\rm S}=\int d^{4}x\sqrt{-g^{E}}\Big{[}\frac{M^{2}_{\rm
P}}{2}R-\frac{1}{2}\partial_{\mu}\phi\partial^{\mu}\phi-V_{\rm
S}(\phi)\Big{]}$ (6)

with the Starobinsky potential

$V_{\rm S}(\phi)=\frac{3M^{2}_{\rm
P}M^{2}}{4}\Big{[}1-e^{-\sqrt{\frac{2}{3}}\frac{\phi}{M_{\rm P}}}\Big{]}^{2}.$
(7)

At this stage, we note that the conformal transformation (5) is a purely
classified transformation of coordinates and results in one frame are
classically equivalent to the ones obtained other frame. Hence, it is
plausible that the number of DOF in the scalar-tensor theory (6) is three
because of two TT tensor modes and one scalar mode. From (6), the Einstein and
scalar equations are derived as

$\displaystyle G_{\mu\nu}=\frac{1}{M^{2}_{\rm
P}}T^{\phi}_{\mu\nu},\leavevmode\nobreak\ \leavevmode\nobreak\
T^{\phi}_{\mu\nu}=\partial_{\mu}\phi\partial_{\nu}\phi-\frac{1}{2}g_{\mu\nu}\Big{[}(\partial\phi)^{2}+V_{\rm
S}\Big{]},$ (8) $\displaystyle\nabla^{2}\phi-V^{\prime}_{\rm
S}=0,\leavevmode\nobreak\ \leavevmode\nobreak\ V^{\prime}_{\rm
S}=\sqrt{\frac{3}{2}}M_{\rm P}M^{2}e^{-\sqrt{\frac{2}{3}}\frac{\phi}{M_{\rm
P}}}\Big{[}1-e^{-\sqrt{\frac{2}{3}}\frac{\phi}{M_{\rm P}}}\Big{]}.$ (9)

The above describes a process of [$R^{2}\to
R\psi-3\psi^{2}\to-(\partial\phi)^{2}-V$] briefly.

##  3 Two sets of perturbed equations

We introduce the metric perturbation around the Minkowski background to find
out the propagating DOF

$\displaystyle g_{\mu\nu}=\eta_{\mu\nu}+h_{\mu\nu}.$ (10)

The Taylor expansions around $R=0$ are employed to define the linearized Ricci
scalar $\delta R(h)$ as [23]

$\displaystyle f(R)$ $\displaystyle=$ $\displaystyle f(0)+f^{\prime}(0)\delta
R(h)+\cdots,$ (11) $\displaystyle f^{\prime}(R)$ $\displaystyle=$
$\displaystyle f^{\prime}(0)+f^{\prime\prime}(0)\delta R(h)+\cdots$ (12)

with $f(0)=0$, $f^{\prime}(0)=1$, and $f^{\prime\prime}(0)=1/3M^{2}$. We note
that $\delta R(h)$ will be used here instead of $R^{(1)}$ in [20]. The
perturbed (linearized) equation to (2) is given by the fourth-order coupled
equation

$\displaystyle\delta
R_{\mu\nu}(h)+\frac{1}{3M^{2}}\Big{[}\eta_{\mu\nu}\Big{(}-\frac{3M^{2}}{2}+\square\Big{)}-\partial_{\mu}\partial_{\nu}\Big{]}\delta
R(h)=0,\leavevmode\nobreak\ \leavevmode\nobreak\ \square=\partial^{2},$ (13)

where the linearized Ricci tensor and scalar are given by

$\displaystyle\delta
R_{\mu\nu}(h)=\frac{1}{2}\Big{[}\partial^{\rho}\partial_{\mu}h_{\nu\rho}+\partial^{\rho}\partial_{\nu}h_{\mu\rho}-\square
h_{\mu\nu}-\partial_{\mu}\partial_{\nu}h\Big{]},$ (14) $\displaystyle\delta
R(h)=\partial^{\rho}\partial^{\sigma}h_{\rho\sigma}-\square h.$ (15)

When using (14) and (15), the linearized equation (13) becomes a second
(fourth)-order differential equation with respect to $h_{\mu\nu}(h)$.
Obviously, it is not a tractable equation. Furthermore, its trace equation
leads to the linearized Ricci scalar equation

$(\square-M^{2})\delta R=0.$ (16)

Introducing the linearized Einstein tensor $\delta G_{\mu\nu}=\delta
R_{\mu\nu}-\eta_{\mu\nu}\delta R/2$, Eq.(13) takes a compact form

$\displaystyle\delta
G_{\mu\nu}(h)+\frac{1}{3M^{2}}\Big{[}\eta_{\mu\nu}\square-\partial_{\mu}\partial_{\nu}\Big{]}\delta
R(h)=0.$ (17)

We note that the Bianchi identity is satisfied when acting $\partial^{\mu}$ on
(17).

On the other hand, two linearized equations from (8) and (9) together with
$\phi=0+\varphi$ take the simple forms with $\delta T^{\phi}_{\mu\nu}=0$ and
$\delta R=0$

$\displaystyle\delta R_{\mu\nu}(h)=0,$ (18)
$\displaystyle(\square-M^{2})\varphi=0.$ (19)

We note that Eqs.(16) and (19) are the same when replacing $\delta R$ by
$\varphi$, but the fourth-order coupled equation (13) is quite different from
the linearized Einstein equation (18). This indicates that (13) can be reduced
to two decoupled second-order equations (18) and (19) if one employs the
conformal transformation and redefinition of scalar appropriately after
choosing (3), leading to a canonical scalar action with the Starobinsky
potential in the Einstein frame. In the scalar-tensor theory approach, one
assigns the perturbed Ricci scalar to an independent scalar $\varphi$.
Instead, it does not have the trace of metric perturbation $h$.

##  4 Propagating DOF with the Lorentz gauge

In order to take into account the propagating DOF in $f(R)$ gravity, it is
convenient to separate the metric tensor $h_{\mu\nu}$ into the traceless part
$h^{T}_{\mu\nu}$ and the trace part $h$ as

$h_{\mu\nu}=h^{\rm T}_{\mu\nu}+\frac{h}{4}\eta_{\mu\nu}$ (20)

with $h^{{\rm T}\mu}_{\leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \mu}\equiv h^{\rm T}=0$. This splitting is meaningful
because the issue of DOF in $f(R)$ theories is related to the presence of $h$.

First of all, let us choose the Lorentz (harmonic) gauge to eliminate the
gauge DOF

$\partial_{\mu}h^{\mu\nu}=\frac{1}{2}\partial^{\nu}h\to\partial_{\mu}h^{{\rm
T}\mu\nu}=\frac{1}{4}\partial^{\nu}h.$ (21)

Here, we note that the transverse condition of $\partial_{\mu}h^{{\rm
T}\mu\nu}=0$ cannot be achieved in $f(R)$ gravity because of $h\not=0$. Then,
the linearized Ricci tensor and scalar are given by

$\delta R(h)_{\mu\nu}\equiv\delta R^{\rm
T}_{\mu\nu}+\frac{\eta_{\mu\nu}}{4}\delta R=-\frac{1}{2}\Big{[}\square h^{\rm
T}_{\mu\nu}+\frac{\eta_{\mu\nu}}{4}\square h\Big{]},\leavevmode\nobreak\
\leavevmode\nobreak\ \delta R(h)=-\frac{1}{2}\square h,$ (22)

where the last equation indicates that the linearized Ricci scalar exists iff
$h\not=0$ under the Lorentz gauge. This implies that $\delta R$ cannot be
defined without $h$. That is, if $h=0$, $\delta R=0$.

Now, the fourth-order equation (13) leads to

$\square
h_{\mu\nu}+\frac{1}{3M^{2}}\Big{[}\eta_{\mu\nu}\Big{(}-\frac{3M^{2}}{2}+\square\Big{)}-\partial_{\mu}\partial_{\nu}\Big{]}\square
h=0.$ (23)

The other form of (23) takes the form

$\square h^{\rm
T}_{\mu\nu}+\frac{1}{3M^{2}}\Big{[}\eta_{\mu\nu}\Big{(}-\frac{3M^{2}}{4}+\square\Big{)}-\partial_{\mu}\partial_{\nu}\Big{]}\square
h=0.$ (24)

If $\square h\not=0$, its trace equation takes the form

$(\square-M^{2})\square h=0$ (25)

which is actually the same equation as in (16). Here, Eq.(25) implies

$(\square-M^{2})h=0$ (26)

because $h$ could represent a massive (scalar) graviton mode in $f(R)$
gravity. The other case of $\square h=0$ is not allowed since if $\square
h=0$, one could not derive the trace equation (25) itself. Importantly, this
issue should be carefully treated because the massless mode satisfying
$\square h=0$ may correspond to the breathing mode, which is the main subject
of this work. In general, it seems that the solution of Eq.(25) is given by
the sum of the massive mode and massless mode which are independent with each
other. However, the massless mode which is a solution to $\square h=0$ does
not exist in $f(R)$ gravity. We stress that $h$ plays the role of a
propagating massive mode instead of $\delta R$. Here is the additional reason
to understand why the massless mode (breathing mode) cannot survive in $f(R)$
gravity. If one requires $\square h=0[\delta R=0$, via (22)], Eq.(24) reduces
to $\square h^{\rm T}_{\mu\nu}=0$, which is just the tensor equation in
Einstein gravity when choosing the Lorentz gauge. It is worth noting that the
last fourth-order term of (24) indicates a feature of the perturbed $f(R)$
gravity clearly. If one chooses $\square h=0$, this term disappears.
Therefore, we clarify that the massless scalar mode does not exist.

Acting $\partial^{\mu}$ on (24) leads to (21), which implies that the Lorentz-
gauge condition is satisfied in the perturbed equation level. Plugging (25)
into (23), we have

$\square
h_{\mu\nu}-\frac{1}{3M^{2}}\square\Big{[}\frac{M^{2}}{2}\eta_{\mu\nu}+\partial_{\mu}\partial_{\nu}\Big{]}h=0.$
(27)

Substituting (25) into (24) implies a fourth-order equation

$\square h^{\rm
T}_{\mu\nu}+\frac{1}{3M^{2}}\square\Big{[}\frac{M^{2}}{4}\eta_{\mu\nu}-\partial_{\mu}\partial_{\nu}\Big{]}h=0$
(28)

which shows clearly that the traceless metric perturbation $h^{\rm
T}_{\mu\nu}$ is closely coupled to the trace of metric perturbation $h$. It is
clear that the trace mode $h$ cannot be decoupled from the traceless tensor
mode $h^{\rm T}_{\mu\nu}$. This is the origin of difficulty met when taking
into account DOF arisen from the $f(R)$ gravity. Interestingly, Eq.(28) can be
transformed to the Ricci tensor-Ricci scalar equation

$\delta R^{\rm
T}_{\mu\nu}+\frac{1}{3M^{2}}\Big{[}\frac{M^{2}}{4}\eta_{\mu\nu}-\partial_{\mu}\partial_{\nu}\Big{]}\delta
R=0,$ (29)

which indicates that the traceless Ricci tensor is coupled to the Ricci
scalar.

At this stage, we observe that Eq.(28) may become a massless propagating
tensor equation for $\tilde{h}_{\mu\nu}$ as

$\tilde{h}_{\mu\nu}\equiv h^{\rm
T}_{\mu\nu}+\frac{1}{3}\Big{(}\frac{\eta_{\mu\nu}}{4}-\frac{\partial_{\mu}\partial_{\nu}}{M^{2}}\Big{)}h,\leavevmode\nobreak\
\leavevmode\nobreak\ \square\tilde{h}_{\mu\nu}=0,$ (30)

which suggests a way of defining a massless spin-2 in $f(R)$ gravity. That is,
plugging $h^{\rm T}_{\mu\nu}$ in (30) into (28) leads to
$\square\tilde{h}_{\mu\nu}=0$. We find from (20) that $h^{\rm T}_{\mu\nu}$
differs from $h_{\mu\nu}$ by $h$. Especially, the splitting of $h$ in
$\tilde{h}_{\mu\nu}$ is nontrivial, which reflects a feature of $f(R)$
gravity. This is compared to that of $h_{\mu\nu}$ in Einstein gravity. For
this purpose, we may express $\tilde{h}_{\mu\nu}$ as

$\displaystyle\tilde{h}_{\mu\nu}$ $\displaystyle=$ $\displaystyle
h_{\mu\nu}-\frac{1}{3}\Big{(}\frac{\eta_{\mu\nu}}{2}+\frac{\partial_{\mu}\partial_{\nu}}{M^{2}}\Big{)}h$
(31) $\displaystyle=$
$\displaystyle\bar{h}_{\mu\nu}-\frac{1}{3}\Big{(}\eta_{\mu\nu}-\frac{\partial_{\mu}\partial_{\nu}}{M^{2}}\Big{)}\bar{h},$
(32)

where $\bar{h}_{\mu\nu}=h_{\mu\nu}-\eta_{\mu\nu}h/2$ is the trace-reversed
metric perturbation $(\bar{h}=-h)$ [21]. Here, the Lorentz gauge is given by
$\partial^{\mu}\bar{h}_{\mu\nu}=0$. We may rewrite (23) in term of
$\bar{h}_{\mu\nu}$

$\square\bar{h}_{\mu\nu}-\frac{1}{3M^{2}}\Big{[}\eta_{\mu\nu}\square-\partial_{\mu}\partial_{\nu}\Big{]}\square\bar{h}=0,$
(33)

which is surely the same equation found in Ref.[19]. Substituting
$\bar{h}_{\mu\nu}$ defined in (32) into that in (33) arrives at
$\square\tilde{h}_{\mu\nu}=0$ which is the same equation as in (30).

Using the trace equation (26) and the Lorentz gauge (21), we may impose the TT
condition for $\tilde{h}_{\mu\nu}$

$\partial_{\mu}\tilde{h}^{\mu\nu}=0,\leavevmode\nobreak\ \leavevmode\nobreak\
\tilde{h}^{\mu}\leavevmode\nobreak\ _{\mu}=0,$ (34)

which indicates that $\tilde{h}_{\mu\nu}$ is a newly tensor mode defined in
$f(R)$ gravity. This may imply that $f(R)$ gravity accommodates three DOF of
two from $\tilde{h}_{\mu\nu}$ and one from $h$. We note that if $h=0$,
$\tilde{h}_{\mu\nu}$ reduces to $h^{\rm T}_{\mu\nu}$ and to $h^{\rm
TT}_{\mu\nu}$ finally, leading to Einstein gravity.

Considering a gravitational wave that propagates in the $z$ direction, Eq.(30)
together with (34) exhibits the fact that gravitational waves have two
polarization components. Explicitly, Eq.(34) implies

$\tilde{h}_{tt}=\tilde{h}_{ti}=0,\leavevmode\nobreak\ \leavevmode\nobreak\
\tilde{h}^{i}\leavevmode\nobreak\ _{i}=0,\leavevmode\nobreak\
\leavevmode\nobreak\ \partial^{i}\tilde{h}_{ij}=0,$ (35)

where the last two expressions correspond to the TT gauge.
$\tilde{h}_{ij}^{\rm TT}=\tilde{h}_{ij}^{\rm TT}(t-z)$ is a valid solution to
the TT wave equation $\square\tilde{h}_{ij}^{\rm TT}=0$. The TT gauge
condition of $\partial_{z}\tilde{h}_{zj}^{\rm TT}=0$ implies
$\tilde{h}_{zj}^{\rm TT}(t-z)$=constant and however, this component should be
zero to satisfy a condition of the asymptotic flatness: $\tilde{h}_{\mu\nu}\to
0$ as $z\to\infty$. The remaining non-zero components of $\tilde{h}_{ij}^{\rm
TT}$ are given by $\tilde{h}_{xx}^{\rm TT},\leavevmode\nobreak\
\tilde{h}_{xy}^{\rm TT},\leavevmode\nobreak\ \tilde{h}_{yx}^{\rm TT}$ and
$\tilde{h}_{yy}^{\rm TT}$. Requiring the symmetry and traceless condition
leads to the two independent components

$\tilde{h}_{xx}^{\rm TT}=-\tilde{h}_{yy}^{\rm
TT}\equiv\tilde{h}^{+}(t-z),\leavevmode\nobreak\ \leavevmode\nobreak\
\tilde{h}_{xy}^{\rm TT}=\tilde{h}_{yx}^{\rm TT}\equiv\tilde{h}^{\times}(t-z).$
(36)

Now, considering (26), one find the trace solution [24]

$h=h^{0}e^{ik_{\mu}x^{\mu}}\to h(t-v_{G}z)$ (37)

where
$v_{G}=k/\omega=\sqrt{\omega^{2}-M^{2}}/\omega<1(\omega^{2}=M^{2}+k^{2})$ is
the group velocity of a massive scalar graviton.

Finally, we obtain $h_{\mu\nu}$ as the solution to (27) with (31)

$h_{\mu\nu}(t,z)=\tilde{h}^{+}(t-z)e^{(+)}_{\mu\nu}+\tilde{h}^{\times}(t-z)e^{(\times)}_{\mu\nu}+\frac{1}{3}\Big{(}\frac{\eta_{\mu\nu}}{2}-\frac{k_{\mu}k_{\nu}}{M^{2}}\Big{)}h(t-v_{G}z).$
(38)

The other solution $\bar{h}_{\mu\nu}$ as the solution to (33) with (32) takes
the form

$\bar{h}_{\mu\nu}(t,z)=\tilde{h}^{+}(t-z)e^{(+)}_{\mu\nu}+\tilde{h}^{\times}(t-z)e^{(\times)}_{\mu\nu}+\frac{1}{3}\Big{(}\eta_{\mu\nu}+\frac{k_{\mu}k_{\nu}}{M^{2}}\Big{)}h(t-v_{G}z),$
(39)

which is the same solution found in Ref.[19]. This encodes that three DOF of
($\tilde{h}^{+},\tilde{h}^{\times},h$) are found from the $f(R)$ gravity.

In order to compare (38) and (39) with the Ricci scalar-like solution [24, 25,
20], we write down its solution by replacing $\delta R=-\square h/2$ with
$h_{f}$

$h_{\mu\nu}^{\rm
R}(t,z)=A^{+}(t-z)e^{(+)}_{\mu\nu}+A^{\times}(t-z)e^{(\times)}_{\mu\nu}+h_{f}(t-v_{G}z)\eta_{\mu\nu}.$
(40)

Even though ($A^{+},A^{\times},h_{f}$) are similar to
($\tilde{h}^{+},\tilde{h}^{\times},h$), the last term of solution
$h_{\mu\nu}^{\rm R}$ (40) differs from that of $h_{\mu\nu}$ (38).

On the other side of the scalar-tensor theory, Eq.(18) together with $\delta
R=0$ reduces to

$\delta R^{\rm T}_{\mu\nu}=0\to\square h^{\rm T}_{\mu\nu}=0.$ (41)

Also, as was shown in Eq.(19), the scalar mode $\varphi$ is decoupled
completely from the $h^{\rm T}_{\mu\nu}$. Requiring the transverse condition
of $\partial_{\mu}h^{{\rm T}\mu\nu}=0$ (Lorentz condition with $h=0$) leads to
two DOF of $(h^{+},h^{\times})$, which describe the general relativity. Hence,
it is obvious that the scalar-tensor theory has three DOF (one scalar DOF+ two
tensor DOF).

##  5 Discussions

First of all, we note that three DOF of ($\tilde{h}^{+},\tilde{h}^{\times},h$)
(38) are found from analyzing the perturbed $f(R)$ gravity. Here, we did not
introduce the Ricci scalar mode ($\delta R=R^{(1)}$) separately because it is
closely related to the trace of metric tensor $h$. We have solved the fourth-
order coupled equation (27) together with the trace equation (26) directly.

We have found that there is no breathing mode in $f(R)$ gravity. The four DOF
including breathing mode have been obtained in [20] by assuming that the
traceless condition of $h=0$ cannot be imposed on the perturbed $f(R)$
gravity, after counting the Ricci scalar mode. The authors in [20] have
discovered the breathing mode ($\hat{h}_{b}$) from the condition of $h\not=0$
when the background spacetime is not Minkowski. The approach used in [20] was
based on the observation that $R^{(1)}$ is considered as a different mode from
$h$ initially [25]. This might lead to overcounting of DOF. However, noting an
expression of $\delta R=-\frac{\square h}{2}$ (22) in the Lorentz gauge
implies that $\delta R$ is closely connected to $h$.

One might attempt to argue from (25) that the massless mode satisfying
$\square h=0$ may correspond to the breathing mode. In general, it seems that
the solution to the fourth-order equation (25) is given by the sum of the
massive mode and massless mode which are independent with each other. However,
the massless mode which is a solution to $\square h=0(h\not=0)$ does not exist
in $f(R)$ gravity since $\square h=0$ [via (22)] means $\delta R=0$ in the
Lorentz gauge.

Consequently, we have clarified the the issue of DOF in $f(R)$ theories. The
number of polarization modes of gravitational waves in $f(R)$ gravity is still
‘three’ in Minkowski spacetime, which is consistent with the results in the
literature (especially for [19]). Also, we would like to mention that the DOF
counting of $f(R)$ theories should be independent of propagating spacetime.

## Acknowledgement

This work was supported by the 2016 Inje University research grant.

## References

  * [1] S. Nojiri and S. D. Odintsov, eConf C0602061 (2006) 06 [Int. J. Geom. Meth. Mod. Phys. 4, 115 (2007)] [arXiv:hep-th/0601213].
  * [2] T. P. Sotiriou and V. Faraoni, Rev. Mod. Phys. 82, 451 (2010) [arXiv:0805.1726 [gr-qc]].
  * [3] A. De Felice and S. Tsujikawa, Living Rev. Rel. 13, 3 (2010) [arXiv:1002.4928 [gr-qc]].
  * [4] S. ’i. Nojiri and S. D. Odintsov, Phys. Rept. 505, 59 (2011) [arXiv:1011.0544 [gr-qc]].
  * [5] S. Perlmutter et al. [Supernova Cosmology Project Collaboration], Astrophys. J. 517 (1999) 565 [arXiv:astro-ph/9812133].
  * [6] A. G. Riess et al. [Supernova Search Team Collaboration], Astron. J. 116 (1998) 1009 [arXiv:astro-ph/9805201].
  * [7] W. Hu and I. Sawicki, Phys. Rev. D 76, 064004 (2007) doi:10.1103/PhysRevD.76.064004 [arXiv:0705.1158 [astro-ph]].
  * [8] X. Liu et al., Phys. Rev. Lett. 117, no. 5, 051101 (2016) doi:10.1103/PhysRevLett.117.051101 [arXiv:1607.00184 [astro-ph.CO]].
  * [9] A. A. Starobinsky, Phys. Lett. B 91, 99 (1980). doi:10.1016/0370-2693(80)90670-X
  * [10] B. Whitt, Phys. Lett. B 145, 176 (1984). doi:10.1016/0370-2693(84)90332-0
  * [11] S. Ferrara, A. Kehagias and A. Riotto, Fortsch. Phys. 62, 573 (2014) doi:10.1002/prop.201400018 [arXiv:1403.5531 [hep-th]].
  * [12] P. A. R. Ade et al. [Planck Collaboration], arXiv:1502.02114 [astro-ph.CO].
  * [13] T. Multamaki and I. Vilja, Phys. Rev. D 74, 064022 (2006) doi:10.1103/PhysRevD.74.064022 [astro-ph/0606373].
  * [14] A. de la Cruz-Dombriz, A. Dobado and A. L. Maroto, Phys. Rev. D 80, 124011 (2009) Erratum: [Phys. Rev. D 83, 029903 (2011)] doi:10.1103/PhysRevD.83.029903, 10.1103/PhysRevD.80.124011 [arXiv:0907.3872 [gr-qc]].
  * [15] T. Moon, Y. S. Myung and E. J. Son, Gen. Rel. Grav. 43, 3079 (2011) doi:10.1007/s10714-011-1225-3 [arXiv:1101.1153 [gr-qc]].
  * [16] F. S. N. Lobo and M. A. Oliveira, Phys. Rev. D 80, 104012 (2009) doi:10.1103/PhysRevD.80.104012 [arXiv:0909.5539 [gr-qc]].
  * [17] K. A. Bronnikov, M. V. Skvortsova and A. A. Starobinsky, Grav. Cosmol. 16, 216 (2010) doi:10.1134/S0202289310030047 [arXiv:1005.3262 [gr-qc]].
  * [18] B. P. Abbott et al. [LIGO Scientific and Virgo Collaborations], Phys. Rev. Lett. 116, no. 6, 061102 (2016) doi:10.1103/PhysRevLett.116.061102 [arXiv:1602.03837 [gr-qc]].
  * [19] E. Barausse and T. P. Sotiriou, Phys. Rev. Lett. 101, 099001 (2008) doi:10.1103/PhysRevLett.101.099001 [arXiv:0803.3433 [gr-qc]].
  * [20] H. R. Kausar, L. Philippoz and P. Jetzer, Phys. Rev. D 93, no. 12, 124071 (2016) doi:10.1103/PhysRevD.93.124071 [arXiv:1606.07000 [gr-qc]].
  * [21] E. E. Flanagan and S. A. Hughes, New J. Phys. 7, 204 (2005) doi:10.1088/1367-2630/7/1/204 [gr-qc/0501041].
  * [22] Y. S. Myung, Phys. Rev. D 84, 024048 (2011) doi:10.1103/PhysRevD.84.024048 [arXiv:1104.3180 [gr-qc]].
  * [23] Y. S. Myung, T. Moon and E. J. Son, Phys. Rev. D 83, 124009 (2011) doi:10.1103/PhysRevD.83.124009 [arXiv:1103.0343 [gr-qc]].
  * [24] S. Capozziello, C. Corda and M. F. De Laurentis, Phys. Lett. B 669, 255 (2008) doi:10.1016/j.physletb.2008.10.001 [arXiv:0812.2272 [astro-ph]].
  * [25] C. P. L. Berry and J. R. Gair, Phys. Rev. D 83, 104022 (2011) Erratum: [Phys. Rev. D 85, 089906 (2012)] doi:10.1103/PhysRevD.85.089906, 10.1103/PhysRevD.83.104022 [arXiv:1104.0819 [gr-qc]].


∎

e1e-mail<EMAIL_ADDRESS>e2e-mail<EMAIL_ADDRESS>e3e-mail:
<EMAIL_ADDRESS>
11institutetext: Centro de Física Teórica e Computacional, Faculdade de
Ciências, Universidade de Lisboa, 1749-016 Lisboa, Portugal 22institutetext:
Centro de Física das Universidades do Minho e do Porto, Campus de Gualtar,
4710-057 Braga, Portugal 33institutetext: Departamento de Física, Faculdade de
Ciências, Universidade de Lisboa, 1749-016 Lisboa, Portugal

Three-dimensional shell-like structures can be obtained spontaneously at the
microscale from the self-folding of 2D templates of rigid panels. At least for
simple structures, the motion of each panel is consistent with a Brownian
process and folding occurs through a sequence of binding events, where pairs
of panels meet at a specific closing angle. Here, we propose a lattice model
to describe the dynamics of self-folding. As an example, we study the folding
of a pyramid of $N$ lateral faces. We combine analytical and numerical Monte
Carlo simulations to find how the folding time depends on the number of faces,
closing angle, and initial configuration. Implications for the study of more
complex structures are discussed.

# Lattice model for self-folding at the microscale

T. S. A. N. Simõese1,addr1,addr2  H. P. M. Meloe2,addr1  N. A. M.
Araújoe3,addr1,addr3

(Received: date / Accepted: date)

††journal:

##  1 Introduction

Folding is considered a promising strategy to design shape-changing materials
with fine tuned mechanical, electrical, and optical properties
shyu2015kirigami ; blees2015graphene ; liu2018nano ; zhang2015mechanically .
The idea is to act over a single degree of freedom to promote a reversible
change in shape, as in the case of the Miura-Ori fold miura1985method ;
mahadevan2005self ; schenk2013geometry . The folding trajectories are usually
deterministic and the final folded configuration unique. However, with
technology putting pressure to go to smaller and smaller scales as, for
example, in encapsulation and soft robotics fernandes2012self ;
shim2012buckling ; filippousi2013polyhedral ; felton2014method , new
challenges are posed. At the microscale, thermal fluctuations are no longer
negligible and so the folding process is stochastic. Experimental and
numerical results show that the final folded configuration is no longer
unique, but rather dependent on the experimental conditions
pandey2011algorithmic ; araujo2018finding ; dodd2018universal ;
azam2009compactness .

Thermal fluctuations have also advantages. They can drive the folding process
and, if the templates are properly designed, there is no need to control any
degree of freedom. Recently, the spontaneous folding of a pyramid from a
planar template consisting of rigid micron-size panels connected through
flexible hinges has been studied using molecular dynamics simulations Melo2020
. The numerical results suggest that, when thermal fluctuations dominate, the
motion of individual panels is well described by a Brownian process and
folding evolves through a sequence of binding events between pairs of faces,
which correspond to first-passage processes.

The first binding event occurs when the first pair of lateral faces meet at
the closing angle $\phi$. The characteristic time is the first binding time
$T_{F}$ and can be estimated by mapping the problem into a two-dimensional
first-passage process Melo2020 . The remaining $N-2$ faces move freely until
they also have the same closing angle. Therefore, the subsequent binding
events can all be mapped to a one-dimensional first-passage processes Melo2020
. When the last face binds, the folding process is complete. The time between
the first and last binding is called the last binding time $T_{L}$. The total
folding time $T$ is then $T=T_{F}+T_{L}$.

Due to the stochastic nature of the folding process, the timescale of binding
is typically much larger than the one of thermal jiggling of the individual
panels. Thus, to observe the complete folding, one needs very long molecular
dynamics simulations. In order to access the folding timescales efficiently,
here, we propose a lattice model. We also consider the folding of a pyramid
with $N$ lateral faces. We show that the dynamics of folding can be mapped
into a set of random walks on a lattice, each one representing a lateral face,
and the binding events occur when two random walks visit a predefined set of
lattice sites simultaneously. We show that our model is consistent with
existing results and reproduces the relation between folding time and number
of faces $N$, reported previously Melo2020 . The proposed lattice model allows
us to explore the folding dynamics for different structures and initial
conditions. The model is generic and can be generalized to simulate the
folding of more complex structures, by including new interaction zones and/or
changing the interaction rules.

The paper is organized as follows. In sect. 2, we introduce the lattice model
and we define the first and last binding times. In sect. 3 we present
numerical and analytical results for the total folding time. We draw some
conclusions in sect. 4.

Figure 1:  Schematic representation of the folding of a pyramid of three
lateral faces and closing angle $\phi$, starting from a flat template, and the
corresponding dynamics on the lattice model. The horizontal axis represents
the angle $\theta$ of each face at a certain time $t$, in the domain
$\theta\in\left[0,\pi\right]$ with reflective boundaries, and the vertical
axis is the time $t$. The motion of each face is mapped into the movement of a
random walker in a 1D. The first binding time $T_{F}$ is the time at which the
first pair of walkers meet at the closing angle $\phi$ and bind together. The
folding time $T$ is when all walkers meet at the closing angle, thus
completing the folding process. The last binding time $T_{L}$ is defined to be
the time interval between first binding and the folding time, so
$T_{L}=T-T_{F}$. The excluded volume interaction between walkers and with
boundaries is described by elastic collisions as depicted. The corresponding
configuration for four times ($t=0$, $0<t<T_{F}$, $t=T_{F}$ and $t=T$) is also
shown.

##  2 Model

We map the spontaneous folding of a regular pyramid with $N$ lateral faces
into a set of random walks on a one-dimensional lattice. Each face $i$ is a
random walk and the angle $\theta_{i}$ relative to the base is given by the
position on a 1D lattice of $l$ sites. To study the folding of a regular
pyramid without misfolding, we consider that the domain of the lattice ranges
from $\theta\in\left[0,\pi\right]$, such that the site $j$ corresponds to an
angle in the range $\left[j\frac{\pi}{l},(j+1)\frac{\pi}{l}\right]$, see model
scheme in Fig. 1.

We consider a rejection-free kinetic Monte Carlo scheme andersen2019practical
; stamatakis2014kinetic ; sickafus2007radiation , where, at each time step, a
randomly selected walker tries to hop to one of its neighbors. Time is
incremented by $\Delta t=\frac{\Delta\theta^{2}}{2D}\frac{1}{N_{F}}$, where
$\Delta\theta=\frac{\pi}{l}$ is the lattice spacing, $D$ is the angular
diffusion coefficient of the lateral faces, and $N_{F}$ is the number of
walkers that are still free. For simplicity, we consider an unbiased folding,
and thus the hopping probability is the same in both directions.

The binding between any two lateral faces occurs when two random walkers first
meet at a specific site corresponding to the closing angle $\phi$.
Additionally, when two faces are in the region beyond the closing angle
($\theta>\phi$) there is a geometrical constraint that they can not overpass
each other. To include this interaction, we assume that in the region
$\theta\in\left]\phi,\pi\right]$, random walkers interact through excluded
volume: If a walker tries to hop into an already occupied site, it goes back
to the original site.

##  3 Results

For each number of random walkers $N$, we performed $10^{4}$ independent
samples, starting with all walkers at the origin of the lattice, which
correspond to the open template shown in Fig. 1. Figure 2 shows the average
folding time $\langle T\rangle$, where we recover the non-monotonic dependence
on the number of lateral faces $N$, reported previously Melo2020 . This
behavior stems from a different functional dependence on $N$ of the mean first
binding time $\langle T_{F}\rangle$ and the mean last binding time $\langle
T_{L}\rangle$, as discussed below.

###  3.1 First binding time

We first estimate the dependence of $\left\langle T_{F}\right\rangle$ on $N$.
For a pyramid with $N$ lateral faces, the time that the first two faces bind
can be estimated in the following way. Let us define $g(t)$ as the probability
density function that a pair of Brownian particles meet at an angle $\phi$ at
time $t$ for the first time. Assuming there is no correlation between the
movement of the $N_{p}=\binom{N}{2}=\frac{N(N-1)}{2}$ pairs, for each sample,
$T_{F}=min(t_{1},t_{2},t_{3}...t_{N_{p}})$, where $t_{i}$ is a random number
drawn from the distribution $g(t)$, for each pair $i$. From the theory of
order statistics arnold1992first  the expected time for the first binding is,

$\left\langle
T_{F}(N_{p})\right\rangle=N_{p}\int_{0}^{\infty}{tg(t)G(t)^{N_{p}-1}dt},$ (1)

where $G(t)=\int_{t}^{\infty}g(t^{\prime})dt^{\prime}$ is the survival
function. Using partial integration we obtain,

$\left\langle T_{F}(N_{p})\right\rangle=\int_{0}^{\infty}{e^{N_{p}\ln
G(t)}dt}.$ (2)

Equation 2 shows that, since $\ln G(t)\leq 1$, the first binding time $T_{F}$
is a monotonic decreasing function of $N$, as also observed with the lattice
model (see Fig. 2). The value of $\left\langle T_{F}\right\rangle$ for
$N\to\infty$ can be estimated for different geometries and initial conditions
yuste1996order ; weiss1983order ; holcman2014narrow ; singer2006narrow ;
yuste1997escape . For all faces starting at $\theta_{i}(0)=0$ (planar initial
condition), it was shown that $\left\langle T_{F}\right\rangle\sim
1/\ln(N_{p})$ weiss1983order ; basnayake2019asymptotic , which is in agreement
with the simulation results shown in Fig. 3(a).

Figure 2: Dependence of the folding time on the number of lateral faces $N$.
$\langle T_{F}\rangle$, $\langle T_{L}\rangle$ and $\langle T\rangle$ are the
average first binding, last binding and folding time, respectively. Time is in
units of Brownian time, i.e. $\frac{\pi^{2}}{2D}$, where $D$ is the diffusion
coefficient of each face. We started from a planar condition, where
$\theta_{i}=0$ for all faces $i$, and a fixed closing angle to
$\phi=\frac{2}{3}\pi$. Results are averages over $10^{4}$ independent samples
on a lattice of $l=181$ sites and the error bars are given by the standard
error.

For the random initial condition, all walkers start at a lattice site selected
uniformly at random. Since there is a chance of two faces to start close to
$\phi$, we can assume that $G(t)$ decays rapidly with $t$, and for large
values of $N$ the most significant contribution for the integral in Eq. 2 is
from $t$ near zero. Therefore, using the approximation $\ln
G(t)\approx\ln(G_{0})+\frac{G^{\prime}_{0}}{G_{0}}t$, where $G_{0}=G(t=0)=1$
and $G^{\prime}_{0}=\frac{dG}{dt}(t=0)$, the first binding time can be
calculated as

$\left\langle T_{F}(N)\right\rangle\approx-\frac{2}{N(N-1)G^{\prime}_{0}}.$
(3)

Since $G^{\prime}_{0}<0$, Eq. 3 implies that the first binding time decays
with $1/N(N-1)$ in perfect agreement with the simulation results shown in Fig.
3(a).

Figure 3: Dependence on the initial conditions. We consider two initial
conditions: all faces start at the origin (planar) or uniformly distributed
(random). (a) First binding time as a function of $N$ for the two initial
conditions. The solid and dashed lines are given by
$T_{F}=\omega_{P_{0}}+\omega_{P}/\ln(N(N-1)/2)$ where
$\omega_{P_{0}}=-0.0632\pm 0.0003$ and $\omega_{P}=1.199\pm 0.002$ for the
planar, and $T_{F}=2\omega_{R}/N(N-1)$ with $\omega_{R}=2.830\pm 0.007$ for
the random, all obtained by least mean square fitting. (b) Last binding time
as a function of $N$. Asymptotically for both initial conditions time scales
logarithmically with $N$, as predicted by Eq. 5. The solid and dashed lines
given by $T_{L}=\tau_{P_{0}}+\tau_{P}\ln(N-2)$ where $\tau_{P_{0}}=0.196\pm
0.002$ and $\tau_{P}=0.3662\pm 0.0006$ for the planar, and
$T_{L}=\tau_{R_{0}}+\tau_{R}\ln(N-2)$ where $\tau_{R_{0}}=0.171\pm 0.002$ and
$\tau_{R}=0.3217\pm 0.0006$ for the random, all obtained by least square
fitting. (c) Total folding time. All times are in units of Brownian time and a
fixed closing angle of $\phi=\frac{2}{3}\pi$. Results are averages over
$10^{4}$ independent samples on a lattice of $l=181$ sites and the error bars
are given by the standard error.

###  3.2 Last binding time

After the first closing, the two faces that form a bond stay fixed at the
angle $\phi$, while the remaining $N-2$ faces move until they reach the same
angle for the first time. To determine $\left\langle T_{L}\right\rangle$, we
calculate the average time for the binding of the last face. From a first-
passage time probability $f(t)$, we have
$T_{L}=max(t_{1},t_{2},t_{3}...t_{N-2})$, where the average time is given by

$\left\langle
T_{L}(N)\right\rangle=(N-2)\int_{0}^{\infty}{tf(t)\left[1-F(t)\right]^{(N-2)-1}dt},$
(4)

where $F(t)=\int_{t}^{\infty}f(t^{\prime})dt^{\prime}$. We assume an uniform
distribution of $\theta_{i}<\phi$ for all the $N-2$ free faces at the moment
of the first binding. For large values of $t$, the distribution of first-
passage times is well described by $f(t)\approx e^{-t/\tau_{L}}$, with
$\tau_{L}=4\phi^{2}/D\pi^{2}$ redner2001guide , and the average time for the
last binding is

$\left\langle T_{L}(N)\right\rangle=\tau_{L}\sum_{i=1}^{N-2}\frac{1}{i}.$ (5)

For large $N$, $\left\langle
T_{L}(N)\right\rangle\approx\tau_{L}\ln(N-2)+\gamma\tau_{L}$, where $\gamma$
is the Euler-Mascheroni constant. In Fig. 3(b) we show that the last binding
time grows with the logarithm of $N$ as predicted by our analytical
calculations for both initial conditions.

For low values of $N$, the first binding time is large and, therefore, after
the first binding the $N-2$ free faces are uniformly distributed, $T_{L}$ is
the same for the planar and random initial condition (see Fig. 3(b)). However,
for large values of $N$, we have a fast first binding and the remaining $N-2$
free faces will still be close to $\theta_{i}=0$ for the planar initial
condition. This is consistent with the fact that $T_{L}$ is larger for the
planar than the random initial condition (see Fig. 3(b)).

Since the average folding time is the sum between $\langle T_{F}\rangle$ and
$\langle T_{L}\rangle$, the different functional dependencies for these two
quantities explains the non-monotonic behavior of $\langle T\rangle$ shown in
Fig. 3(c) for both initial conditions, in line with previous results obtained
using molecular dynamic simulations Melo2020 .

###  3.3 Closing angle

Since the overall dynamics of all faces is diffusive, we would expect that the
time scale that controls the binding process is the Brownian time
$\phi^{2}/2D$. Figure 4 shows simulation results for the average first (a) and
last (b) binding times as a function of $\phi$, for different values of the
number of faces $N$. The solid and dashed lines are given by fitting the
simulation data to a power law $\langle T_{F}\rangle=A_{F}\phi^{m_{F}}$ and
$\langle T_{L}\rangle=A_{L}\phi^{m_{L}}$, using least mean square method, with
fitting parameters $(A_{F},m_{F})$ and $(A_{L},m_{L})$. However, the results
show that the binding time tends to be proportional to $\phi^{2}$ only for
large values of $N$. In particular, for $N=2$ and $N=3$, we see that there is
a slight deviation from the quadratic scaling. This deviation is due to the
fact that the first binding is a two-dimensional first-passage process where
the binding site is located inside the domain, and the time scale is
controlled by the size of the boundary. However, for large values of $N$ the
path for the first binding will be close to a straight line
lawley2020probabilistic , recovering the timescale $\phi^{2}/2D$ as shown in
Fig. 4(c).

Figure 4: Time dependence with the closing angle $\phi$. (a) First binding
time $\langle T_{F}\rangle$ as a function of $\phi$ for $N=\\{2,3,100\\}$. (b)
Last binding time $\langle T_{L}\rangle$ as a function of $\phi$ for
$N=\\{3,100\\}$. In (a) and (b) the data was fitted by a power law $\langle
T_{L/F}\rangle\sim\phi^{m}$, represented by the dashed lines. From the theory
of first-passage processes, the time scale should be proportional to the
square of the domain size, however since the adsorbing site is not on the
boundary of the lattice we can have deviations from that scaling,
corresponding to $m\neq 2.0$. We show in (c) that for both the first and last
binding times, $m$ converges to $2.0$ as we increase $N$, since for large $N$
the first face to reach $\phi$ closely follows a straight path
lawley2020probabilistic . All times are in units of Brownian time and
$\phi\in\left[\frac{19\pi}{36},\frac{35\pi}{36}\right]$. Results are averages
over $10^{4}$ independent samples with planar initial conditions on a lattice
of $l=181$ sites and the error bars are given by the standard error.

##  4 Discussion and conclusions

We proposed a lattice model to simulate the spontaneous folding of a pyramid
at the microscale, driven by thermal fluctuations. We map the angular motion
of each face into a random walk on a lattice. The pairwise binding of faces
corresponds to having two walkers at a specific lattice site. We recover a
recent result obtained from Molecular Dynamics simulations, namely, that the
average folding time is a non-monotonic function of the number of lateral
faces $N$. The folding dynamics involves two types of first-passage processes:
The first binding, where the first two faces bind together; and the last
binding, where the last face binds. We show that the first binding corresponds
to a two-dimensional first-passage process, where two random walkers have to
bind for the first time at a specific lattice site. After that, the remaining
free walkers bind through a sequence of one-dimensional first-passage
processes. The total folding time is the sum of the first and last binding
times.

Describing folding as a sequence of binding events, we demonstrated that the
characteristic time of the first binding has to decrease with $N$, while for
the last binding it increases. It is the balance between these two processes
that leads to a non-monotonic dependence on the number of faces. We show that
this non-monotonic dependence is robust to changes on the initial conditions
or value of the closing angle. Although the initial condition could affect
both first and last binding times, it only affects significantly the first
binding time and only for large $N$. For a small number of lateral faces the
binding times slightly deviates from a quadratic dependence on $\phi$, however
for large $N$ the quadratic dependence is recovered, as expected for a
diffusive process.

The lattice model allows to evaluate not only the folding time, but also the
yield of more complex structures by extending the domain and including new
closing angles. For example, an additional closing angle at $-\phi$
corresponds to the possibility of folding in two different sides of the base
of the pyramid, with possible misfolding for a pyramid with $N>3$, where the
final folded configuration combines faces at $\phi$ and $-\phi$. Analytical
techniques developed for first-passage process with multiple adsorbing
boundaries can be used to calculate the yield as a function of $N$.

###### Acknowledgements.

We acknowledge financial support from the Portuguese Foundation for Science
and Technology (FCT) under Contracts No. PTDC/FIS-MAC/28146/2017
(LISBOA–01–0145–FEDER–028146), UIDB/00618/2020, and UIDP/00618/2020.

## References

  * (1) T.C. Shyu, P.F. Damasceno, P.M. Dodd, A. Lamoureux, L. Xu, M. Shlian, M. Shtein, S.C. Glotzer, N.A. Kotov, Nat. Mater. 14, 785 (2015).
  * (2) M.K. Blees, A.W. Barnard, P.A. Rose, S.P. Roberts, K.L. McGill, P.Y. Huang, A.R. Ruyack, J.W. Kevek, B. Kobrin, D.A. Muller, et al., Nature 524, 204 (2015).
  * (3) Z. Liu, H. Du, J. Li, L. Lu, Z.Y. Li, N.X. Fang, Sci. Adv. 4, eaat4436 (2018).
  * (4) Y. Zhang, Z. Yan, K. Nan, D. Xiao, Y. Liu, H. Luan, H. Fu, X. Wang, Q. Yang, J. Wang, et al., Proc. Natl. Acad. Sci. U.S.A. 112, 11757 (2015).
  * (5) K. Miura, Inst. Space Astronaut. Sci. Rep. 618, 1 (1985).
  * (6) L. Mahadevan, S. Rica, Science 307, 1740 (2005).
  * (7) M. Schenk, S.D. Guest, Proc. Natl. Acad. Sci. U.S.A. 110, 3276 (2013).
  * (8) R. Fernandes, D.H. Gracias, Adv. Drug Deliv. Rev. 64, 1579 (2012).
  * (9) J. Shim, C. Perdigou, E.R. Chen, K. Bertoldi, P.M. Reis, Proc. Natl. Acad. Sci. U.S.A. 109, 5978 (2012).
  * (10) M. Filippousi, T. Altantzis, G. Stefanou, M. Betsiou, D.N. Bikiaris, M. Angelakeris, E. Pavlidou, D. Zamboulis, G. van Tendeloo, RSC Adv. 3, 24367 (2013).
  * (11) S. Felton, M. Tolley, E. Demaine, D. Rus, R. Wood, Science 345, 644 (2014).
  * (12) S. Pandey, M. Ewing, A. Kunas, N. Nguyen, D.H. Gracias, G. Menon, Proc. Natl. Acad. Sci. U.S.A. 108, 19885 (2011).
  * (13) N.A.M. Araújo, R.A. da Costa, S.N. Dorogovtsev, J.F.F. Mendes, Phys. Rev. Lett. 120, 188001 (2018).
  * (14) P.M. Dodd, P.F. Damasceno, S.C. Glotzer, Proc. Natl. Acad. Sci. U.S.A. 115, E6690 (2018).
  * (15) A. Azam, T.G. Leong, A.M. Zarafshar, D.H. Gracias, PLOS ONE 4, e4451 (2009).
  * (16) H.P.M. Melo, C.S. Dias, N.A.M. Araújo, Commun. Phys. 3, 154 (2020).
  * (17) M. Andersen, C. Panosetti, K. Reuter, Front. Chem. 7, 202 (2019).
  * (18) M. Stamatakis, J. Condens. Matter Phys. 27, 013001 (2014).
  * (19) K.E. Sickafus, E.A. Kotomin, B.P. Uberuaga, _Radiation effects in solids_ , vol. 235 (Springer Science & Business Media, 2007).
  * (20) B.C. Arnold, N. Balakrishnan, H.N. Nagaraja, _A first course in order statistics_ , vol. 54 (Siam, 1992).
  * (21) S.B. Yuste, K. Lindenberg, J. Stat. Phys. 85, 501 (1996).
  * (22) G.H. Weiss, K.E. Shuler, K. Lindenberg, J. Stat. Phys. 31, 255 (1983).
  * (23) D. Holcman, Z. Schuss, SIAM Rev. 56, 213 (2014).
  * (24) A. Singer, Z. Schuss, D. Holcman, J. Stat. Phys. 122, 491 (2006).
  * (25) S.B. Yuste, Phys. Rev. Lett. 79, 3565 (1997).
  * (26) K. Basnayake, Z. Schuss, D. Holcman, J. Nonlinear Sci. 29, 461 (2019).
  * (27) S. Redner, _A guide to first-passage processes_ (Cambridge University Press, 2001).
  * (28) S.D. Lawley, J.B. Madrid, J. Nonlinear Sci. 30, 1207 (2020).