# Siloxane molecules: Nonlinear elastic behavior and fracture characteristics

Tianchi Li  [  [  Eric R. Dufresne  [  Martin Kröger  Polymer Physics,
Department of Materials, ETH Zurich, CH–8093 Zurich, Switzerland  Stefanie
Heyden  [<EMAIL_ADDRESS>
###### Abstract

Fracture phenomena in soft materials span multiple length- and timescales.
This poses a major challenge in computational modeling and predictive
materials design. To pass quantitatively from molecular- to continuum scales,
a precise representation of the material response at the molecular level is
vital. Here, we derive the nonlinear elastic response and fracture
characteristics of individual siloxane molecules using molecular dynamics (MD)
studies. For short chains, we find deviations from classical scalings for both
the effective stiffness and mean chain rupture times. A simple model of a non-
uniform chain of Kuhn segments captures the observed effect and agrees well
with MD data. We find that the dominating fracture mechanism depends on the
applied force scale in a non-monotonic fashion. This analysis suggests that
common polydimethylsiloxane (PDMS) networks fail at crosslinking points. Our
results can be readily lumped into coarse-grained models. Although focusing on
PDMS as a model system, our study presents a general procedure to pass beyond
the window of accessible rupture times in MD studies employing mean first
passage time theory, which can be exploited for arbitrary molecular systems.

###### keywords:

siloxane molecules, nonlinear elastic response, fracture characteristics

A]Soft and Living Materials, Department of Materials, ETH Zurich, CH–8093
Zurich, Switzerland A]Soft and Living Materials, Department of Materials, ETH
Zurich, CH–8093 Zurich, Switzerland A]Soft and Living Materials, Department of
Materials, ETH Zurich, CH–8093 Zurich, Switzerland Magnetism and Interface
Physics, Department of Materials, ETH Zurich, CH–8093 Zurich, Switzerland
A]Soft and Living Materials, Department of Materials, ETH Zurich, CH–8093
Zurich, Switzerland

##  1 Introduction

Most things in life start small. This basic concept also applies to the
failure of soft materials, emerging from the rupture of interatomic bonds.
Predicting the fracture journey that follows becomes a question of failure
mechanisms and lengthscales 1, 2, 3. In view of failure mechanisms, Lake-
Thomas theory has formed our understanding of how much energy it takes to
break an elastic chain 4. When a crack propagates within a stretched elastic
material, each repeat unit within chains crossing the fracture plane stores
energy. The resultant fracture energy should thus reflect the elastic energy
stored within the entire chain instead of pure single bond scission. Recent
works on tough hydrogels hint at a more complicated picture, in which network
characteristics such as entanglements have a crucial effect on fracture 5, 6.

Multiple lengthscales form the basis of the classical fracture mechanics
picture. In the ideally brittle limit, dissipation and material failure occur
on the scale of the atomistic separation length. In soft tough materials, the
characteristic lengthscale in the continuum limit is the so-called _elasto-
adhesive_ length. This lengthscale is typically microscopic and represents the
region of nonlinear elastic deformation around a macroscopic crack tip 3. It
can be coupled to molecular failure processes at small scales 4, 7, as well as
energy dissipation at the mesoscale 8, 9, 10 and macroscopic effects such as
crack blunting 11, 12. In a recent work, scale-free cavity growth at constant
driving pressure was accessed in the mesoscopic region 13. In this picture, no
well defined crack tip exists and corresponding process zones for the
calculation of fracture energies becomes obsolete.

Fracture in soft solids thus displays manifold characteristics that are
deviating from classical theories. To get further insight into what governs
these deviations, multiple length- and timescales need to be bridged, which
poses a major computational challenge. Here, we address this challenge by
providing a detailed description of the nonlinear elastic response of
molecular building blocks up to fracture. These building blocks can then
provide starting grounds for higher level coarse grained models.

Previous studies on the force-extension relation and fracture of individual
molecules encompass both experimental- and computational investigations.
Experimental studies include atomic force microscopy (AFM) 14, 15, 16, 17,
optical tweezers 18, 19, 20, as well as magnetic tweezers 21, 22, 23. Due to
its large accessible force range up to $\mathcal{O}($nN$)$ and high
resolution, AFM has been widely adopted 24. Investigations using AFM comprise
a wide spectrum, ranging from proteins 25, 26, DNA 27, 28, polysaccharides 29,
poly(ethylene glycol) 30 and poly(methacrylic acid) 31 to polydimethylsiloxane
(PDMS) 32.

Using computational methods, ab initio molecular dynamics (AIMD) 33, 34
simulations allow for an on-the-fly computation of electronic structures based
on quantum mechanics. While bond fracture can be modeled in this setting, high
computational costs limit AIMD studies to $\mathcal{O}($nm$)$ and
$\mathcal{O}($ps$)$ 35, 36, 37. At higher length- and timescales, steered
molecular dynamics (MD) simulations have emerged as the primary method in
studying the force-extension behavior of molecules 38, 39, 40, 41. Classical
MD methods are amenable of treating system sizes of several hundreds of
nanometers and time scales on the order of nanoseconds. However, atomic
interactions are typically modeled _via_ empirical interatomic potentials,
which require a predefined atomic connectivity remaining unchanged throughout
simulations, such that fracture of interatomic bonds cannot be described. As
an alternative, bond-order based force fields were developed to bridge the gap
between ab-initio and empirical force fields. Here, we derive the quasi-static
force-extension and rupture properties of single molecules up to
$\mathcal{O}(100\,\text{nm})$ and $\mathcal{O}($ns$)$ by enriching all-atom
steered molecular dynamics simulations with a bond-order based force field
(_ReaxFF_) 42, 43, 44, 45, 46, with the help of the LAMMPS software package.47
Unlike classical atomistic bond potentials, ReaxFF allows for different atomic
bonding states, such that fracture events can be captured. Simulations thus
reduce the gap between length- and time scales accessible using ab initio
computational methods and experimental approaches. We focus on PDMS as a model
system as used in previous studies 13, for which both linear PDMS and
crosslinked PDMS are investigated.

##  2 Molecular Dynamics Studies

###  2.1 Nonlinear elastic response

Prior to failure, the static molecular response is governed by entropic
elasticity at extensions well below the unstretched contour length, and
enthalpic elasticity at higher extensions. The exact shape of this nonlinear
elastic force-extension relation depends on the specific molecular structure
under investigation. To derive the nonlinear elastic response of siloxane
molecules, PDMS-$n$ molecules of varying polymerization degree $n$ are
created. Figure 1a illustrates the chemical structure of PDMS-$n$ as an
example. Each PDMS-$n$ molecule is embedded in a simulation box, which is set
up both with and without solvent molecules. When solvent molecules are
present, periodic boundary conditions are applied. We use hexamethyldisiloxane
(HMDSO) molecules as a solvent, as interactions between HMDSO and PDMS do not
alter the rupture behavior of PDMS (compared to interactions with itself) 37.
In comparison, trace amounts of water were found to lower the maximally
attained rupture stretch 37.

All simulations are performed with the parameter set specifically trained and
optimized for PDMS 43. Without solvent molecules, the number of degrees of
freedom is $3n_{a}$, with $n_{a}$ being the number of atoms in PDMS-$n$.
$n_{a}$ scales linearly with polymerization degree $n$. For simulations in
which solvent molecules are present, the overall system size increases by
$3n_{s}$ degrees of freedom based on $n_{s}$ solvent molecules. There is no
upper constraint on $n_{s}$. Its lower bound is set by the requirement of
generating sufficiently large RVE’s for subsequent steered molecular dynamics
runs, preventing self-interactions. For simulations of PDMS-27 up to fracture,
$n_{s}/n_{a}\sim 60$.

amounts to $60\times 3n_{a}$ degrees of freedom. Here, the multiplicative
factor of $60$ stems from the presence of solvent molecules. In all
simulations, a timestep of $0.1$ fs is applied. Systems are relaxed in an NPT
ensemble at ambient conditions ($T=300$ K, $p=1$ atm). Following relaxation, a
constant repulsive force $F$ between the 2 terminal Si atoms is applied in an
NVT ensemble. In this ensemble, volume $V$ is held constant, such that the
equilibrated end-to-end distance is purely based on the applied force
(rescalings of the simulation box are prohibited). Results are compared to a
displacement controlled setting, in which both terminal Si atoms are held
constant at fixed end-to-end distance $R$ and the exerted force is recorded.
Figure 1c shows the nonlinear elastic reponse of PDMS-27 up to fracture. The
choice of boundary condition does not influence the force-extension relation
in both entropic ($R\ll L$) and enthalpic ($R>L$) regimes, where the
unstretched contour length $L$ marks the crossover point. For comparison with
classical polymer models, the inset of Figure 1c highlights the divergence of
the freely jointed chain model (FJC) 21 for an end-to-end chain distance $R$
approaching the unstretched contour length $L=6.4\,$nm. In contrast, the
elastic freely jointed chain model (EFJC) 48 captures the nonlinear elastic
force-extension relation also within the enthalpic regime $R>L$. The change of
slope at large forces is encoded in the $F$-dependent bond potential of mean
force $U(b;F)$, as investigated in more detail in Section 2.2.

Figure 1: (a) Chemical structure of PDMS-$n$ illustrating $(n-1)$ repeating
units. (b) Non-uniform chain of springs model consisting of two types of Kuhn
springs. (c) Force-extension relation $F(R)$ of a single PDMS-27 molecule
measured at $T=300$ K and $P=1$ atm in the presence of HMDSO solvent
molecules. Red/black data points are obtained in a displacement/force
controlled setting, respectively. Inset showing the small extension regime
($R<L$). The FJC model (pink solid upper line) and elastic FJC model (blue
solid lower line) both capture the entropic regime $R\ll L$, while the latter
shows better agreement in the enthalpic region $R>L$. (d) Evolution of
$F_{0}=L\,\partial F/\partial R$, obtained at $T=1$ K within the elastic
regime as shown in the inset for $n=2$ up to $n=376$. MD data is captured well
by a model of non-uniform chain of springs, consisting of two types of Kuhn
springs as sketched in (b).

The force-extension curve at large deformation is linear. As expected from the
EFJC model, the force-extension relation of a single polymer chain is given as

$R(F)=L\left(1+\frac{F}{F_{0}}\right){\cal
L}(\xi),\qquad\text{where}\;\;\xi=\frac{FL_{k}}{k_{\text{B}}T},\quad
F_{0}=KL.$ (1)

Equation (1) represents a classical FJC model with an added elastic extension
$(1+F/F_{0})$. Within the entropic regime, elasticity is modeled via the
Langevin function ${\cal L}(\xi)=\coth(\xi)-\xi^{-1}$. With the added elastic
extension, the EFJC model introduces the effective Hookean spring constant $K$
as an additional elastic parameter within the enthalpic regime.

For computational efficiency, solvent molecules are removed for the
determination of $F_{0}$, and $T=1$ K is chosen to reduce thermal noise. All
other simulations are performed at $T=300$ K. Note that in this study, we
focus on the enthalpic regime, in which temperature effects on the mechanical
response become negligible. Insensitivity of $F_{0}$ towards both temperature
and solvent molecules in the enthalpic regime is tested for short oligomers
(see Figures S3 and S4 in the Appendix). We find that both temperature and
solvent molecules do not influence $F_{0}$ in the ethalpic regime. For small
forces $F\ll F_{0}$ in the entropic limit, we have $R(F)\simeq F/\tilde{K}$,
with $\tilde{K}=3k_{B}T/(LL_{k})$ the elastic Hookean spring constant within
the entropic regime. For large forces $F\gg F_{0}$ in the enthalpic limit,
$R(F)\simeq F/K$, as the Langevin function ${\cal L}(\xi)$ approaches unity
for $\xi\gg 1$. $L_{k}$ is obtained from fitting the EFJC model to the
measured force-extension curve at $T=300$ K, see Figure 1c. A comparison to
other hypothetical Kuhn segments (which consistently overpredict the
unstretched contour length $L$) is shown in Figure S2 in the Appendix).
Further support giving an independent estimate of $L_{k}$ from analyzing the
Si-Si vector correlation function is given in supplementary Section S2. Here,
$L_{k}=5.5\pm 0.7$ Å is the Kuhn length of the polymer chain corresponding to
the mean end-to-end distance of a Si-O-Si-O-Si triplet, which in the following
is abbreviated as Si-Si-Si.

For a low polymerization degree of $n=2$, the resultant slope $F_{0}$ when
plotting force $F$ versus stretch $R/L$ attains its maximum value
$F_{0}=31.35\,$nN (see inset in Figure 1d). $F_{0}$ decreases with increasing
$n$ in a nonlinear fashion as illustrated in Figure 1d. Characteristic forces
are calculated for a large range of polymerization degrees
$n\in[2,3,4,5,7,10,15,20,50,100,376]$. This differs from a classical model of
$n$ identical springs in series, for which $K\propto 1/n$ and $L\propto n$. To
determine the source of this deviation, we track bond length distributions at
fixed repulsive force between terminal Si atoms. We find that within the
enthalpic regime, internal triplet distances are shorter than terminal ones.
This difference in triplet distance distributions can be related to
restrictions in bond angles and dihedrals at terminal atoms. Endowing internal
Si-Si-Si Kuhn segments with spring stiffness $k_{i}$ and equilibrium length
$l_{i}$, whereas terminal Kuhn segments possess spring stiffness $k_{t}$ and
equilibrium length $l_{t}$, gives the overall Hookean spring constant $K$ and
contour length $L$ as

$K=\left(\frac{2}{k_{t}}+\frac{n/2-2}{k_{i}}\right)^{-1},\qquad
L=2l_{t}+\left(\frac{n}{2}-2\right)l_{i}.$ (2)

Here, $n/2$ denotes the total number of Kuhn segments, where $n$ is the
polymerization degree. $k_{t}=47.3\,\pm 0.1$ nN/nm and $l_{t}=0.52\pm
0.01\,$nm are directly computed from the bond length distribution of PDMS-4,
which only consists of two terminal Kuhn segments ($k_{t}=2K$ and
$l_{t}=L/2$). Using a fit to simulation results for PDMS-5 to calculate the
remaining free parameters, we find that $k_{i}=40.82\pm 0.01\,$nN/nm and
$l_{i}=0.48\pm 0.01\,$nm of the internal Kuhn segment. As shown by the solid
black line in Figure 1d, this model of a non-uniform chain of springs agrees
well with MD data and captures the effect of terminal springs at small $n$, as
well as convergence of $F_{0}$ to a plateau at large $n$.

To summarize, the elastic response of PDMS oligomers (both within entropic and
enthalpic regimes) is characterized by the Si triplet length $L_{k}$ based on
two reasons: First, the entropic part of the force-extension curve suggests
$L_{k}$ to be identical to the extension of a Si-Si-Si triplet. Second,
$n$-dependencies of $K$ and $L$ are consistently captured only if the number
($n/2$) of Si-Si-Si triplets is used in Equation (2). In sharp contrast, the
fracture behavior to be discussed next will be dominated by the $F$-dependent
characteristics of single covalent atomic bonds.

###  2.2 Fracture characteristics

At the molecular scale, fracture is stochastic. An intuitive question to ask
is ’Where and when does a network tend to break?’. Here, we try to
quantitatively answer this question for PDMS in terms of mean rupture times
and preferred fracture modes.

To distinguish between rupture of the PDMS-$n$ backbone (_chain scission_) and
rupture at crosslinking sites (_crosslink failure_), we take into account two
different structures: PDMS-$n$ as used in the previous Section, as well as two
PDMS-4 molecules linked via a crosslinking site -$CH_{2}$-$CH_{2}$-. Chain
scission thus stems from the rupture of Si-O bonds, while crosslink failure
results from rupturing Si-C bonds. The accessible window of mean rupture times
in MD studies lies in the range of $10^{-2}-10^{0}$ ns. The upper limit is set
by computational feasibility, while the lower limit depends on the molecular
vibration frequency of the polymer chain below which inertial effects dominate
the response.

To extend beyond this rupture time window and determine mean bond rupture
times $\tau(F)$ on longer timescales, we use a statistical extrapolation
scheme. Our approach renders a close analogy to the calculation of mean first
passage times for chemical processes with a single reaction coordinate 49, for
the thermal or enforced breakage of discrete one-dimensional chains 50, Morse-
chains 51 and biomolecules 52. Its derivation is provided in the Supplementary
Information. We proceed in the following way: Stationary equilibrium Si-O and
Si-C bond length probability densities $p(b)$ of PDMS-4, PDMS-376 and linked
PDMS-4 are measured at different levels of constant force (_cf._ Figure 2).
Using $p(b)$, we calculate mean chain rupture times, for which we need to pass
from rupture times of single bonds to those of chains with $2n$ bonds.

Figure 2: (a and b) Stationary nonequilibrium MD probability densities of Si-O
and Si-C bond lengths on linked PDMS-4 at various force levels (simulations
performed in the presence of solvent molecules). Grey curves denote
intermediate force levels. (c and d) Radial bond potential of mean force
calculated from bond length distributions highlighted in (a,b). Black curves
represent polynomial fits of order 4 (Si-O) and 2 (Si-C). All simulations are
performed at $T=300$ K.

At fixed force, directly measured single bond probability densities $p(b)$
serve to define an effective potential, with

$p(b)=\frac{\exp(-U/k_{\text{B}}T)}{\int_{0}^{b_{r}}\exp(-U/k_{\text{B}}T)db}.$
(3)

Here, $b_{r}$ is the rupture bond length. For later calculations of mean chain
rupture times, we need a functional form of the effective bond potential $U$.
In the following, using the notation $U(b;F)$, we emphasize that the potential
is a function of bond length $b$, while its parameters (in this case $b_{1}$)
depend on $F$. Numerically solving for $U$ from Equation (3) at different
levels of applied force, we see that $U$ needs to satisfy the following
properties: It should have the generic form of a double-well potential (fourth
order polynomial) with minima corresponding to two different equilibrium bond
lengths (cf. Figure S7). Here, we denote $b_{1}$ as the equilibrium bond
length corresponding to the first minimum, and choose $U(b_{1})=0$ for
convenience. Most importantly, we observe a nearly $F$-independent parabolic
shape of $U^{\prime\prime}(b)$ about its second minimum at $b\simeq b_{2}$,
i.e., $U^{\prime\prime}(b)=c_{2}+k_{2}(b-b_{2})^{2}$ (cf. Figure S7). Since
$U^{\prime\prime}(b)$ exhibits a parabolic and $F$-independent shape and
location, the corresponding parameters $b_{2}$, $c_{2}$, $k_{2}$ can be
treated as $F$-independent constants. This finding forms the basis of
rendering our statistical extrapolation scheme feasible, since $b_{1}(F)$
remains as the only force-dependent fitting parameter. In addition, we only
observe a weak dependence of $b_{1}$ on $F$, which forms the basis for an
extrapolation to higher force regimes. These observations lead to a functional
form of $U(b;F)$ as

$\displaystyle U(b;F)$ $\displaystyle=$
$\displaystyle\int_{b_{1}}^{b}\int_{b_{1}}^{x}\left[c_{2}+k_{2}(y-b_{2})^{2}\right]\,dy\,dx$
(4) $\displaystyle=$
$\displaystyle\frac{(b-b_{1})^{2}\\{6c_{2}+k_{2}[3b_{1}^{2}+6b_{2}^{2}-8b_{1}b_{2}+2(b_{1}-2b_{2})b+b^{2}]\\}}{12}.$

Bond length distributions (Figures 2a and 2b) are fitted to Equations (3), (4)
with fitting parameters given in Table SI. The resulting radial bond
potentials of mean force are illustrated in Figure 2c and 2d, with
corresponding polynomials of order 4 (Si-O bond) and order 2 (Si-C bond, for
which $b_{2}=k_{2}=0$). Note that fitting deviations in Si-C potentials are
based on restricting $b_{1}(F)$ to be the only force-dependent parameter. The
largest possible instantaneous Si-O an Si-C bond length value in stable chain
configurations is given in Table SI. As shown in Figure 2c, an increasing non-
linearity develops with increasing tension for Si-O bonds. This is rooted in
bond angle potentials losing their dominance within the energy landscape due
to the externally enforced alignment. At this point, the remaining
interactions (dihedral, Si-C, Si-H, O-H) come into play.

With an expression for $U(b;F)$ at hand, we proceed with the calculation of
mean chain rupture times, passing from rupture times of single bonds to those
of chains with $2n$ bonds. Furthermore, it needs to be verified that a theory
neglecting inertia effects captures the attendant fracture characteristics.

Neglecting inertia effects, the mean rupture time of a single bond is
calculated from the Fokker-Planck equation as 53, 54, 55

$\tau(F)=\frac{\zeta}{k_{\text{B}}T}\int_{b_{1}(F)}^{b_{r}}\int_{0}^{z}\frac{\Psi(y)}{\Psi(z)}\,dy\,dz,\qquad\text{with}\quad\Psi(b)=\exp\left[-\frac{U(b;F)}{k_{\text{B}}T}\right].$
(5)

This is a purely theoretical limit, as atomistic simulations (PDMS chains
consist of multiple Si-O and Si-C bonds) measure $\tau_{n}$. $\zeta$ is an a
priori unknown friction coefficient, which will be determined later by
matching theoretical rupture times $\tau_{n}(F)$ with those obtained from MD
simulations. Utilizing the Fokker-Planck approach, Equation (S-2) (or
equivalently Brownian Dynamics simulations via Equation (S-1)) can be used to
explore the rupture time distribution $p(t_{r};F)$ of a single bond, which is
nearly mono-exponential (apart from a small dip at $t_{r}\rightarrow 0$).

In order to pass to the rupture time distribution $p_{n}(t_{r};F)$ of a chain
with polymerization degree $n$ (which thus contains $2n$ bonds), we assume
independent bonds. The probability of a chain (_i.e._ at least one of its
assumed identical bonds) rupturing during time interval $t$ after onset of $F$
at time $t=0$ is

$P_{n}(t;F)=1-\left[1-\int_{0}^{t}p(t_{r};F)dt_{r}\right]^{2n}=1-e^{-2nt/\tau(F)}.$
(6)

The term in parentheses in Equation (6) denotes the probability of an
individual bond staying intact until time $t$. The probability distribution
for rupture times $t_{r}$ of $n$-chains (PDMS-$n$) is thus
$p_{n}(t_{r};F)=(d/dt_{r})P_{n}(t_{r};F)=2ne^{-2nt_{r}/\tau(F)}/\tau(F)$, from
which the mean chain rupture time $\tau_{n}(F)$ follows as
$\tau_{n}(F)=\tau(F)/2n$. We compare these theoretical expressions to those
measured in MD simulations. Figure 3a shows measurements on PDMS-4. MD
measurements show a mono-exponential shape of $p_{4}(t_{r};F)$, which is in
agreement with the Fokker-Planck prediction for a single bond and the
assumption of independent bonds in chains of higher polymerization degree.

Figure 3b illustrates $\tau_{n}(F)$ at three different constant stretching
forces $F$, for which each data point is the average of $10000$ independent
samples. For large $n$, $\tau_{n}(F)$ approaches the expected $\propto 1/n$
limit. Deviations from this scaling for short chains are reminiscent of the
non-uniform chain of springs effect highlighted in Figure 1b. Equivalent to
the functional form given in Equation (2), we have

$\frac{\tau_{n}(F)}{\textrm{ns}}=\frac{1}{a_{0}+(n-1)/a},\qquad\frac{\tau(F)}{\textrm{ns}}=\lim_{n\rightarrow\infty}\frac{2n\tau_{n}(F)}{\textrm{ns}}=2a$
(7)

Fitting parameters at $F=6.25$ nN are obtained as $a_{0}=-24\pm 1$ and
$a=0.108\pm 0.004$, while at $F=6.11$ nN, $a_{0}=-9.1\pm 0.4$ and $a=0.28\pm
0.01$. Single bond mean rupture times depicted in Figure 3c are calculated as
$\tau^{\rm Si-O}(F)=\lim_{n\rightarrow\infty}2n\tau_{n}(F)=2a$ ns. For Si-C
bonds (which are present twice in linked PDMS-4), force levels of
$F=5.21\,$nN, $F=5.28\,$nN and $F=5.35\,$nN are investigated. This force range
in MD already spans two decades in single-bond mean rupture time $\tau(F)$.
The resulting single bond mean rupture time is computed as $\tau^{\rm
Si-C}(F)=2\tau^{\rm linked-PDMS-4}(F)$.

By matching the measured $\tau_{n}$ with the theoretically predicted one, we
can furthermore determine $\zeta$ (which is the shape-preserving, force-
independent vertical shift required to match measurement and theory). Solid
lines in Figure 3c illustrate the $\tau(F)$ resulting from the Fokker-Planck
equation (5). This solution allows to extend beyond the rupture time window
accessible in MD studies (shaded region in Figure 3c). We find that the
lifetime of a single representative Si-O bond is much longer than that of a
Si-C bond, with an increasing gap for larger $F$ (based on the significant
difference in potentials at high forces).

Figure 3: (a) Measured rupture time probability density $p_{4}(t_{r};F)$ of
PDMS-4 at constant stretching force $F=6.46$ nN. Green solid line: Mono-
exponential used in the calculation of mean rupture time $\tau_{n}(F)$. (b)
Evolution of mean chain rupture time $\tau_{n}(F)$ with increasing
polymerization degree $n$ at three different forces. Each data point is the
average of $10000$ simulations. Solid lines are obtained from fitting the
functional form given in Equation (7). For large $n$, $\tau_{n}(F)\propto 1/n$
as expected. Deviations for short chains are related to the effect of non-
uniform chain of springs highlighted in Figure 1d. (c) Single-bond mean
rupture time $\tau(F)$ for Si-O (obtained from (7)) and Si-C (average of
individual MD simulations). Solid lines: Solution to the mean first passage
time problem (5) based on $U$ (cf. Figure 2b-c). The shaded region highlights
the rupture time window accessible in MD studies. All simulations are
performed at $T=300$ K.

To determine preferred fracture mechanisms, we note that with increasing
polymerization degree $n$, $\tau_{n}^{\rm Si-O}(F)$ decreases, while
$\tau_{n}^{\rm Si-C}(F)$ is constant (based on the constant number of Si-C
bonds when focusing on the single chain level, see Figure 4b). Rupture of Si-O
bonds (chain scission) and Si-C bonds (crosslink failure) thus becomes
comparable at a crossover polymerization degree $n_{c}$. Figure 4b displays a
comparison of mean chain rupture times $\tau_{n}(F)$ at different levels of
applied force. Tracking the crossover polymerization degree $n_{c}$ allows to
identify two different failure regimes: For $n<n_{c}$, crosslink failure is
anticipated, while chain scission is the preferred failure mode for $n>n_{c}$.
Figure 4c highlights the effective rupture time (taking into account both Si-O
and Si-C bonds) as contour lines as a function of $n$ and $F$. Again, the
crossover polymerization degree $n_{c}$ differentiates a region dominated by
crosslink failure (green) from a regime dominated by chain scission (blue).

We observe a strengthening effect in Si-O bonds with increasing force, which
is reminiscent of phenomena observed in systems involving catch bonds, e.g.,
membrane-to-surface adhesion 56, myosin and actin 57, or signaling receptors
and their ligands 58, 59. This strengthening emerges as a ’re-entrant’ effect
of crosslink failure for polymerization degrees $n>10^{3}$, which can be
related to the higher order structure of Si-O bond potentials (see the
nonlinearity developing at higher forces, Figures 2, S5 and S9). Si-O bonds
are stable at low forces, at which the failure of crosslinking junctions is
the dominating fracture mechanism. At intermediate forces, chain scission
dominates. At high forces, at which the increasing stiffness of the second
minimum in Si-O potentials comes into play, fracture characteristics are
dominated by crosslink failure again. With increasing $n$, the force regime
dominated by chain scission grows, which is in keeping with Equation (6).

Typical siloxane materials used in the laboratory setting are highlighted in
Figure 4c in terms of polymerization degree $n$. Single molecules in _Sylgard
184_ and _DMS-V31_ are entirely dominated by crosslink failure. In contrast,
individual molecules in _Sylgard 186_ feature a much higher polymerization
degree. As such, their fracture behavior strongly depends on the applied
force, with crosslink failure in the low force regime transitioning to chain
scission at $F>2.5\,$nN.

chain scissioncrosslink failure Figure 4: Crosslink failure versus chain
scission. (a) We model constituent molecules of monodisperse strands, testing
for chain scission along the PDMS backbone (blue) versus failure at
crosslinking junctions (green). (b) Evolution of crossover polymerization
degree $n_{c}$ (dashed black line) with increasing force $F$. Intersections of
$\tau_{n}(F)$ for Si-O along the backbone and $\tau_{2}(F)$ for Si-C at two
crosslinking junctions are illustrated at three different force levels. (c)
Effective rupture time of a network strand (contour lines) as a function of
$n$ and $F$. Green: Region dominated by failure of crosslinking junctions.
Blue: Chain scission dominated region. (A) _Sylgard 184_ ($n=78$), (B)
_DMS-V31_ ($n=376$), (C) _Sylgard 186_ ($n=64848$) 60. All predictions apply
to $T=300$ K.

In the above, we model rupture time distributions using an inertia-free
Fokker-Planck approach. To justify this approach, it remains to investigate
the absence of solvent molecules on rupture time distributions
$p_{n}(t_{r},F)$. While the dynamics of bond lengths exhibits inertia effects,
we find the rupture time distribution to be unaffected by the presence/absence
of explicit solvent molecules, as shown in Figure S8. Without solvent
molecules, inertia effects are maximal and friction is absent. The presence of
solvent molecules provides additional noise and stochastic collisions, such
that inertia effects are diminished. The rupture time distribution is thus
also unaffected by the degree of suppression of inertia effects. Furthermore,
we observe a nearly mono-exponential rupture time probability distribution, in
which the majority of bonds does not break during the first oscillation but at
later times. This is in stark contrast to the dominance of inertial effects,
for which bonds that are still intact after the first oscillation would never
fail (the rupture criterion would not be fulfilled in subsequent oscillations
if it was not already fulfilled in the first oscillation).

##  3 Discussion and Outlook

This work characterizes the nonlinear elastic- and fracture behavior of PDMS.
The nonlinear elastic repsonse of siloxane oligomers is captured well using
the EFJC model, both within entropic- and enthalpic regimes. At low
polymerization degrees $n$, we find a deviation from a classical scaling of
effective stiffness $K\propto 1/n$, which can be captured with a simple model
of non-uniform chains of springs, with each spring constituting a Kuhn
segment.

Passing to the inelastic behavior of PDMS, we focus on rupture times of both
Si-O bonds (present in the backbone of siloxane oligomers), as well as Si-C
bonds (present at crosslinking sites). When calculating mean chain rupture
times of siloxane oligomers, we again observe a deviation from classical
scalings ($\tau_{n}(F)\propto 1/n$) at low polymerization degrees. Similar to
the nonlinear elastic part, a model including the ’end effect’ on mean chain
rupture times agrees well with simulations.

We find that the lifetime of single Si-O bonds is much longer than that of a
Si-C bond. We define a crossover polymerization degree $n_{c}$ at which the
dominating rupture mechanism of the crosslinked network passes from crosslink
failure (rupture of Si-C bonds) to chain scission (rupture of Si-O bonds on
PDMS-$n$). To pass to long timescales, we use a statistical extrapolation
scheme in order to understand bond fracture within the network. Surprisingly,
we find that the dominating fracture mechanism depends on the applied force
scale in a non-monotonous fashion. The non-monotonic dependence of $n_{c}$ on
$F$ is rooted in nonlinearities in the Si-O bond potential and the
corresponding stiffness dependence on force at the two minima.

For single molecules in typical siloxane materials used in the laboratory
setting, such as _Sylgard 184_ and _DMS-V31_ , our analysis suggests breakage
exclusively at crosslinks. For individual molecules in _Sylgard 186_ featuring
higher polymerization degrees however, the attendant failure mode depends on
the magnitude of applied force. While crosslink failure dominates at low
forces, chain scission is expected for $F>2.5\,$nN. It furthermore bears
mentioning that our analysis focuses on the single chain level. Possible
inhomogeneities in force distribution based on network topology (and the
resultant changes in fracture characteristics) furnish an important point for
future studies. These could aid in elucidating network properties passing from
the single chain to the continuum level, focusing on the influence of
polymerization degree and crosslink density.

Our results provide building blocks, which can be readily used in coarse-
grained higher scale models. As an example, network models of nonlinear
springs could be easily tailored to siloxane systems by implementing a
material model corresponding to the nonlinear elastic response derived in this
work. In this setting, crosslinking molecules could be lumped into nodes
representing crosslinking sites. Tracking attendant forces upon deformation
would then allow for the implementation of fracture criteria corresponding to
those derived in this work.

Finally, our work provides a generic procedure to pass beyond the window of
accessible rupture times in MD studies, which can be applied for arbitrary
molecular systems.

SH and TL gratefully acknowledge funding via the SNF Ambizione grant
PZ00P2186041. We are furthermore thankful for many insightful discussions with
our Soft $\&$ Living Materials peers Dr. Robert Style and Dr. Nicolas Bain.

The following file is available free of charge.

  * •

SI.pdf: (i) Characteristic force at finite temperature, (ii) Kuhn length,
(iii) parameters of the double well potential, (iv) mean bond rupture time,
(v) influence of boundary conditions, (vi) bond length distributions and
potentials, (vii) rupture times with and without HMDSO solvent, (viii) mean
chain rupture times for different polymerization degrees.

## References

  * Zhao 2014 Zhao, X. Multi-scale multi-mechanism design of tough hydrogels: building dissipation into stretchy networks. _Soft Matter_ 2014, _10_ , 672–687
  * Bai et al. 2019 Bai, R.; Yang, J.; Suo, Z. Fatigue of hydrogels. _Europ. J. Mech. A_ 2019, _74_ , 337–370
  * Long et al. 2021 Long, R.; Hui, C.-Y.; Gong, J.; Bouchbinder, E. The Fracture of Highly Deformable Soft Materials: A Tale of Two Length Scales. _Annu. Rev. Condens. Matter Phys._ 2021, _12_ , 71–94
  * Lake and Thomas 1967 Lake, G.; Thomas, A. The strength of highly elastic materials. _Proc. R. Soc. A_ 1967, _300_ , 108–119
  * Kim et al. 2021 Kim, J.; Zhang, G.; Shi, M.; Suo, Z. Fracture, fatigue, and friction of polymers in which entanglements greatly outnumber cross-links. _Science_ 2021, _374_ , 212–216
  * Nian et al. 2022 Nian, G.; Kim, J.; Bao, X.; Suo, Z. Making Highly Elastic and Tough Hydrogels from Doughs. _Adv. Mater._ 2022, _34_ , e2206577
  * Baumberger et al. 2006 Baumberger, T.; Caroli, C.; Martina, D. Solvent control of crack dynamics in a reversible hydrogel. _Nat. Mater._ 2006, _5_ , 552–5
  * Brown 2007 Brown, H. A Model of the Fracture of Double Network Gels. _Macromolecules_ 2007, _40_ , 3815–3818
  * Tanaka 2007 Tanaka, Y. A local damage model for anomalous high toughness of double-network gels. _Europhys. Lett._ 2007, _78_ , 56005
  * Zhang et al. 2015 Zhang, T.; Lin, S.; Yuk, H.; Zhao, X. Predicting fracture energies and crack-tip fields of soft tough materials. _Extreme Mech. Lett._ 2015, _4_ , 1–8
  * Hui et al. 2003 Hui, C.-Y.; Jagota, A.; Bennison, S.; Londono, J. Crack blunting and the strength of soft elastic solids. _Proc. R. Soc. A_ 2003, _459_ , 1489–1516
  * Seitz et al. 2009 Seitz, M. E.; Martina, D.; Baumberger, T.; Krishnan, V. R.; Hui, C.-Y.; Shull, K. R. Fracture and large strain behavior of self-assembled triblock copolymer gels. _Soft Matter_ 2009, _5_ , 447–456
  * Kim et al. 2020 Kim, J. Y.; Liu, Z.; Weon, B. M.; Cohen, T.; Hui, C.-Y.; Dufresne, E. R.; Style, R. W. Extreme cavity expansion in soft solids: Damage without fracture. _Sci. Adv._ 2020, _6_ , eaaz0418
  * Binnig et al. 1986 Binnig, G.; Quate, C. F.; Gerber, C. Atomic Force Microscope. _Phys. Rev. Lett._ 1986, _56_ , 930–933
  * Xu et al. 2002 Xu, Q.; Zhang, W.; Zhang, X. Oxygen bridge inhibits conformational transition of 1,4-linked alpha-D-galactose detected by single-molecule atomic force microscopy. _Macromolecules_ 2002, _35_ , 871–876
  * Gunari et al. 2006 Gunari, N.; Schmidt, M.; Janshoff, A. Persistence length of cylindrical brush molecules measured by atomic force microscopy. _Macromolecules_ 2006, _39_ , 2219–2224
  * Yang et al. 2018 Yang, P.; Song, Y.; Feng, W.; Zhang, W. Unfolding of a Single Polymer Chain from the Single Crystal by Air-Phase Single-Molecule Force Spectroscopy: Toward Better Force Precision and More Accurate Description of Molecular Behaviors. _Macromolecules_ 2018, _51_ , 7052–7060
  * Ashkin et al. 1986 Ashkin, A.; Dziedzic, J. M.; Bjorkholm, J. E.; Chu, S. Observation of a single-beam gradient force optical trap for dielectric particles. _Opt. Lett._ 1986, _11_ , 288–290
  * Wang et al. 1997 Wang, M.; Yin, H.; Landick, R.; Gelles, J.; Block, S. Stretching DNA with optical tweezers. _Biophys. J._ 1997, _72_ , 1335–1346
  * Rocha et al. 2018 Rocha, M. S.; Storm, I. M.; Bazoni, R. F.; Ramos, E. B.; Hernandez-Garcia, A.; Stuart, M. A. C.; Leermakers, F.; de Vries, R. Force and Scale Dependence of the Elasticity of Self-Assembled DNA Bottle Brushes. _Macromolecules_ 2018, _51_ , 204–212
  * Smith et al. 1992 Smith, S. B.; Finzi, L.; Bustamante, C. Direct Mechanical Measurements of the Elasticity of Single DNA Molecules by Using Magnetic Beads. _Science_ 1992, _258_ , 1122–1126
  * Strick et al. 1996 Strick, T.; Allemand, J.; Bensimon, D.; Bensimon, A.; Croquette, V. The Elasticity of a Single Supercoiled DNA Molecule. _Science_ 1996, _271_ , 1835–7
  * del Rio et al. 2009 del Rio, A.; Perez-Jimenez, R.; Liu, R.; Roca-Cusachs, P.; Fernandez, J. M.; Sheetz, M. P. Stretching Single Talin Rod Molecules Activates Vinculin Binding. _Science_ 2009, _323_ , 638–641
  * Butt et al. 2005 Butt, H.-J.; Cappella, B.; Kappl, M. Force measurements with the atomic force microscope: Technique, interpretation and applications. _Surf. Sci. Rep._ 2005, _59_ , 1–152
  * Rief et al. 1997 Rief, M.; Gautel, M.; Fern, J. M.; Gaub, H. E. reversible unfolding of individual titin immunoglobulin dimains by AFM. _Science_ 1997, 1109–1112
  * Fisher et al. 1999 Fisher, T.; Oberhauser, A.; Carrion-Vazquez, M.; Marszalek, P.; Fernandez, J. The study of protein mechanics with the atomic force microscope. _Trends Biochem. Sci._ 1999, _25_ , 379–84
  * Rief et al. 1999 Rief, M.; Clausen-Schaumann, H.; Gaub, H. E. Sequence-dependent mechanics of single DNA molecules. _Nat. Struct. Biol._ 1999, _6_ , 346–349
  * Bustamante et al. 2000 Bustamante, C.; Smith, S. B.; Liphardt, J.; Smith, D. Single-molecule studies of DNA mechanics. _Curr. Opin. Struct. Biol._ 2000, _10_ , 279–285
  * Rief et al. 1997 Rief, M.; Oesterhelt, F.; Heymann, B.; Gaub, H. E. Single molecule force spectroscopy on polysaccharides by atomic force microscopy. _Science_ 1997, _275_ , 1295–1297
  * Oesterhelt et al. 1999 Oesterhelt, F.; Rief, M.; Gaub, H. E. Single molecule force spectroscopy by AFM indicates helical structure of poly(ethylene-glycol) in water. _New J. Phys._ 1999, _1_ , 6
  * Ortiz and Hadziioannou 1999 Ortiz, C.; Hadziioannou, G. Entropic elasticity of single polymer chains of poly(methacrylic acid) measured by atomic force microscopy. _Macromolecules_ 1999, _32_ , 780–787
  * Schwaderer et al. 2008 Schwaderer, P.; Funk, E.; Achenbach, F.; Weis, J.; Bräuchle, C.; Michaelis, J. Single-Molecule Measurement of the Strength of a Siloxane Bond †. _Langmuir_ 2008, _24_ , 1343–9
  * Kresse and Hafner 1993 Kresse, G.; Hafner, J. Ab initio molecular dynamics for liquid metals. _Phys. Rev. B_ 1993, _47_ , 558–561
  * Sun 1995 Sun, H. Ab initio calculations and force field development for computer simulation of polysilanes. _Macromolecules_ 1995, _28_ , 701–712
  * Lupton et al. 2005 Lupton, E.; Nonnenberg, C.; Frank, I.; Achenbach, F.; Weis, J.; Bräuchle, C. Stretching siloxanes: An ab initio molecular dynamics study. _Chem. Phys. Lett._ 2005, _414_ , 132–137
  * Lupton et al. 2009 Lupton, E.; Achenbach, F.; Weis, J.; Bräuchle, C.; Frank, I. Origins of Material Failure in Siloxane Elastomers from First Principles. _ChemPhysChem_ 2009, _10_ , 119–23
  * Lupton et al. 2006 Lupton, E.; Achenbach, F.; Weis, J.; Bräuchle, C.; Frank, I. Modified Chemistry of Siloxanes under Tensile Stress: Interaction with Environment. _J. Phys. Chem. B_ 2006, _110_ , 14557–63
  * Lu et al. 1998 Lu, H.; Isralewitz, B.; Krammer, A.; Vogel, V.; Schulten, K. Unfolding of Titin Immunoglobulin Domains by Steered Molecular Dynamics Simulation. _Biophys. J._ 1998, _75_ , 662–671
  * Lu and Schulten 1999 Lu, H.; Schulten, K. Steered molecular dynamics simulations of force-induced protein domain unfolding. _Proteins: Struct. Function Bioinf._ 1999, _35_ , 453–463
  * Isralewitz et al. 2001 Isralewitz, B.; Baudry, J.; Gullingsrud, J.; Kosztin, D.; Schulten, K. Steered molecular dynamics investigations of protein function. _J. Molec. Graph. Model._ 2001, _19_ , 13–25
  * Han et al. 2022 Han, Z.; Hilburg, S. L.; Alexander-Katz, A. Forced Unfolding of Protein-Inspired Single-Chain Random Heteropolymers. _Macromolecules_ 2022, _55_ , 1295–1309
  * Van Duin et al. 2001 Van Duin, A.; Dasgupta, S.; Lorant, F.; Goddard, W. ReaxFF: A reactive force field for hydrocarbons. _J. Phys. Chem. A_ 2001, _105_ , 9396–9409
  * Chenoweth et al. 2005 Chenoweth, K.; Cheung, S.; van Duin, A.; Goddard, W.; Kober, E. Simulations on the Thermal Decomposition of a Poly(dimethylsiloxane) Polymer Using the ReaxFF Reactive Force Field. _J. Amer. Chem. Soc._ 2005, _127_ , 7192–202
  * Newsome et al. 2012 Newsome, D. A.; Sengupta, D.; Foroutan, H.; Russo, M. F.; Duin, A. C. T. Oxidation of Silicon Carbide by O2 and H2O: A ReaxFF Reactive Molecular Dynamics Study, Part I. _J. Phys. Chem. C_ 2012, _116_ , 16111–16121
  * Soria et al. 2017 Soria, F. A.; Zhang, W.; van Duin, A. C. T.; Patrito, E. M. Thermal Stability of Organic Monolayers Grafted to Si(111): Insights from ReaxFF Reactive Molecular Dynamics Simulations. _ACS Appl. Mater. Interf._ 2017, _9_ , 30969–30981
  * Soria et al. 2018-10-18 Soria, F. A.; Zhang, W.; Paredes-Olivera, P. A.; van Duin, A. C. T.; Patrito, E. M. Si/C/H ReaxFF Reactive Potential for Silicon Surfaces Grafted with Organic Molecules. _J. Phys. Chem._ 2018-10-18, _122_
  * Plimpton 1995 Plimpton, S. Fast Parallel Algorithms for Short-Range Molecular Dynamics. _J. Comput. Phys._ 1995, _117_ , 1–19
  * Smith et al. 1996 Smith, S. B.; Cui, Y.; Bustamante, C. Overstretching B-DNA: The Elastic Response of Individual Double-Stranded and Single-Stranded DNA Molecules. _Science_ 1996, _271_ , 795–799
  * Preston et al. 2021 Preston, R. J.; Gelin, M. F.; Kosov, D. S. First-passage time theory of activated rate chemical processes in electronic molecular junctions. _J. Chem. Phys._ 2021, _154_ , 114108
  * Razbin et al. 2019 Razbin, M.; Benetatos, P.; Moosavi-Movahedi, A. A. A first-passage approach to the thermal breakage of a discrete one-dimensional chain. _Soft Matter_ 2019, _15_ , 2469–2478
  * Puthur and Sebastian 2002 Puthur, R.; Sebastian, K. Theory of polymer breaking under tension. _Phys. Rev. B_ 2002, _66_ , 024304
  * Berezhkovskii and Makarov 2019 Berezhkovskii, A. M.; Makarov, D. E. On distributions of barrier crossing times as observed in single-molecule studies of biomolecules. _Biophys. Rep._ 2019, _1_ , 100029
  * Risken 1996 Risken, H. _The Fokker-Planck Equation_ ; Springer: Berlin, 1996
  * Kampen 2007 Kampen, N. V. _Stochastic Processes in Physics and Chemistry, 3rd Ed._ ; North Holland: Amsterdam, The Netherlands, 2007
  * Gardiner 1985 Gardiner, C. W. _Handbook of Stochastic Methods, 2nd Ed._ ; Springer, Berlin: Berlin, 1985
  * Dembo et al. 1988 Dembo, M.; Torney, D.; Saxman, K.; Hammer, D. The reaction-limited kinetics of membrane-to-surface adhesion and detachment. _Proc. R. Soc. B_ 1988, _234_ , 55–83
  * Guo and Guilford 2006 Guo, B.; Guilford, W. Mechanics of actomyosin bonds in different nucleotide states are tuned to muscle contraction. _Proc. Natl. Acadm. Sci. USA_ 2006, _103_ , 9844–9
  * Liu et al. 2014 Liu, B.; Chen, W.; Evavold, B.; Zhu, C. Accumulation of dynamic catch bonds between TCR and agonist peptide-MHC triggers T cell signaling. _Cell_ 2014, _157_ , 357–68
  * Das et al. 2016 Das, D.; Mallis, R.; Duke-Cohan, J.; Hussey, R.; Tetteh, P.; Hilton, M.; Wagner, G.; Lang, M.; Reinherz, E. Pre-T Cell Receptors (Pre-TCRs) Leverage Vbeta Complementarity Determining Regions (CDRs) and Hydrophobic Patch in Mechanosensing Thymic Self-ligands. _J. Biolog. Chem._ 2016, _291_ , 25292–305
  * Flowers and Switzer 1978 Flowers, G. L.; Switzer, S. T. Background material properties of selected silicone potting compounds and raw materials for their substitutes. 1978, United States Government document
  * Evmenenko et al. 2006 Evmenenko, G.; Mo, H.; Kewalramani, S.; Dutta, P. Conformational rearrangements in interfacial region of polydimethylsiloxane melt films. _Polymer_ 2006, _47_ , 878–882

## Supplementary Information

Siloxane molecules: Nonlinear elastic behavior and fracture characteristics

Tianchi Li1, Eric R. Dufresne1, Martin Kröger2,3, Stefanie Heyden1∗)
1) Soft and Living Materials, Department of Materials, ETH Zurich, CH–8093
Zurich, Switzerland

2) Polymer Physics, Department of Materials, ETH Zurich, CH–8093 Zurich,
Switzerland

3) Magnetism and Interface Physics, Department of Materials, ETH Zurich,
CH–8093 Zurich, Switzerland
∗) Corresponding author<EMAIL_ADDRESS>(S.H.)

###  S1 Characteristic force at finite temperature

Figure S1: Characteristic force (equivalent to Fig. 1(d) in the main
manuscript) at $T=300$ K. The end effect (decreasing $F_{0}$ with increasing
$n$) is also observed at finite temperature. Added noise due to entropic
contributions impede the characterization of $F_{0}(n)$ within the enthalpic
region, for which we still recover $F_{0}\propto(n+1)/n$. In combination with
a study of bond length distributions along the chain (for all forces in the
enthalpic part, the two terminal Kuhn springs (Si-Si-Si) are longer than all
internal springs), this leads us to the non-uniform chain of springs model.

###  S2 Kuhn length

Figure S2: Measured equilibrium contour length (blue circles) in comparison to
different models of Kuhn segments (Si-O, Si-Si, and Si-Si-Si). Using force-
extension simulations at $T=1$ K, the measured equilibrium contour length
(specified in units of nm) follows from the intercept of the enthalpic part on
the horizontal axis at $F=0$. Following the EFJC model, the equilibrium
contour length is the sum of all Kuhn segment lengths. We find that taking
Si-O or Si-Si as a Kuhn segment overestimates the measured equilibrium contour
length, while taking Si-Si-Si as the elementary Kuhn unit is in good agreement
with measured values.

We obtained an independent estimate of the Kuhn length of the PDMS chains by
analyzing the Si-Si vector correlation function $C(i)=\langle{\bf
u}_{j}\cdot{\bf u}_{j+i}\rangle$, where ${\bf u}_{j}$ denotes the unit vector
parallel to the vector connecting the $j$th and $(j+1)$th Si atom along the
PDMS backbone, and the average is taken over all $j\in\\{1,..,n\\}$ within an
ensemble of equilibrium PDMS-$n$ chains. For the FRC model, $\ln
C(i)=-\ell/L_{p}$, where $L_{p}$ is the persistence length, and $\ell\approx
2.932\pm 0.003$ Å the measured average distance between adjacent Si atoms. We
obtain $L_{p}=2.67\pm 0.07$ Å ($n=10$), $L_{p}=2.86\pm 0.09$ Å ($n=20$),
$L_{p}=2.69\pm 0.06$ Å ($n=50$), $L_{p}=2.77\pm 0.05$ Å ($n=100$). The Kuhn
length $L_{k}$ is twice as large as the persistence length, confirming
$L_{k}\approx 5.5$ Å. This is in agreement with prior literature estimates 61.

###  S3 Parameters of the double well potential

| $b_{r}$ | $k_{2}$ | $c_{2}$ | $b_{2}$ | $b_{1}$
---|---|---|---|---|---
System | [nm] |  [kg.nm-2.s-2] |  [kg.s-2] | [nm] |  [nm] $\qquad\tilde{F}\equiv F/$ nN
PDMS 4 (8 Si-O bonds) | 0.22 | $579866$ | $-55.9156$ | $0.186$ | $0.15471+0.00085126\,\tilde{F}+0.0002648\,\tilde{F}^{2}$
PDMS 376 (752 Si-O bonds) | 0.22 | $621285$ | $-41.419$ | $0.185$ | $0.15471+0.00085126\,\tilde{F}+0.0002648\,\tilde{F}^{2}$
linked PDMS 4 (2 Si-C bonds) | 0.255 | $0$ | $186$ | $0$ | $0.20042+0.0029917\,\tilde{F}+0.00031558\,\tilde{F}^{2}$
linked PDMS 4 (1 C-C bond) |  (∗)0.181 | $0$ | $800$ | $0$ | $0.151-0.0005\,\tilde{F}$
Table SI: Parameters of the effective bond potential $U(b;F)$ given by Eq.
(4). All results have been produced using the stated values, whose error is
less than 3% for $k_{2}$ and $c_{2}$, and less than 1% for $b_{1}$ and
$b_{2}$. For PDMS 4 (4 Si-O bonds) we obtained $\tau_{4}(F=90$
kcal.mol-1.Å${}^{-1}=6.2531$ nN$)=0.221\pm 0.002$ ns from atomistic simulation
at $T=300$ K. This value is reproduced via Brownian dynamics using $\zeta=105$
ng/m2. For linked PDMS 4 (2 Si-C bonds) $b_{r}=0.255$ nm and we obtained
$\tau_{2}(F=90$ kcal.mol-1.Å${}^{-1})=0.70\pm 0.04$ ps and $\tau_{1}(F=80$
kcal.mol-1.Å${}^{-1})=0.94\pm 0.05$ ps and $\tau_{1}(F=75$
kcal.mol-1.Å${}^{-1})=51.3\pm 0.8$ ps from atomistic simulation at $T=300$ K.
This value is reproduced via Brownian dynamics using $\zeta=0.002$ ng/m2. (∗)
literature estimate.

###  S4 Mean bond rupture time

Consider a Brownian bond whose time-dependent length $b(t)$ resides within the
interval $[b_{0},b_{r}]$. The bond is assumed to change its length resulting
from three types of forces: Deterministic forces due to a one-dimensional
potential $U(b;F)$ of mean force, which we determine from atomistic simulation
at constant force $F$, a frictional force (friction coefficient $\zeta$)
resulting from the surrounding medium, and a stochastic force whose strength
is governed by the fluctuation-dissipation theorem. The Langevin equation for
the bond length $b$ thus reads 53

$\frac{d}{dt}b=-\frac{1}{\zeta}\frac{dU(b;F)}{db}+\sqrt{\frac{2k_{\text{B}}T}{\zeta}}\eta(t),$
(S-1)

where $\eta(t)$ represents uncorrelated white noise, $\langle\eta(t)\rangle=0$
and $\langle\eta(t)\eta(t^{\prime})\rangle=\delta(t-t^{\prime})$. We further
consider an adsorbing boundary at $b=b_{r}$ (the rupture bond length) and a
reflecting boundary at $b=b_{0}$. Let the conditional probability
$p_{2}(b^{\prime},t|b,0)$ distribution capture the probability that a bond,
whose length is $b$ at time $0$, assumes length $b^{\prime}$ at a later time
$t^{\prime}\geq 0$, with $b,b^{\prime}\in[b_{0},b_{r}]$. Inline with our
assumptions, $p_{2}(b^{\prime},t|b,0)$ solves an adjungated Fokker-Planck
equation corresponding to the Langevin Eq. (S-1) 55, 54

$\displaystyle-\frac{\partial}{\partial t}p_{2}(b^{\prime},t|b,0)$
$\displaystyle=$ $\displaystyle-\left[-\frac{1}{\zeta}\frac{\partial
U(b;F)}{\partial b}\frac{\partial}{\partial
b}+\frac{k_{\text{B}}T}{\zeta}\frac{\partial^{2}}{\partial
b^{2}}\right]p_{2}(b^{\prime},t|b,0)$ (S-2)

subject to initial condition $p_{2}(b^{\prime},0|b,0)=\delta(b^{\prime}-b)$
and the abovementioned constraints. Then

$G(b_{r},t|b)=\int_{b_{0}}^{b_{r}}p_{2}(b^{\prime},t|b,0)\,db^{\prime}$ (S-3)

is the probability that $b(t)$ resides within the interval $[b_{0},b_{r}]$ at
time $t$. The $G$ is thus not normalized except in the limit
$b_{r}\rightarrow\infty$, i.e., one has $G(\infty,t|b)=1$. Further
$G(b_{r},0|b)=1$ since $p_{2}(b^{\prime},0|b,0)=\delta(b^{\prime}-b)$, and
$G(b_{r},\infty|b)=0$, since the bond length exceeds $b_{r}$ with a nonzero
probability. One can write down an equation for $G$ based on the equation
(S-2) for $p_{2}$ 53. Due to the boundary conditions for $p_{2}$, the boundary
conditions for $G$ read $G(b_{r},0|b)=1$ for $b\in[b_{0},b_{r}]$ and
$G(b_{r},0|b)=0$ otherwise, and

$\left.G(b_{r},t|b)\right|_{b=b_{r}}=0,\qquad\left.\frac{\partial}{\partial
b}G(b_{r},t|b)\right|_{b=b_{0}}=0.$ (S-4)

Because we are interested in the mean rupture time, we introduce the fraction
$f(b_{r},t|b)$ of bonds that reach $b_{r}$ (and thus leave the interval
$[b_{0},b_{r}]$) within the time interval $[t,t+dt]$. One has

$-dG(b_{r},t|b)=-\partial_{t}G(b_{r},t|b)dt\equiv f(b_{r},t|b)dt,\qquad
f(b_{r},t|b)=-\partial_{t}G(b_{r},t|b).$ (S-5)

The quantity

$T_{1}(b_{r},b)=\int_{0}^{\infty}tf(b_{r},t|b)dt=-\int_{0}^{\infty}t\partial_{t}G(b_{r},t|b)dt=\int_{0}^{\infty}G(b_{r},t|b)dt$
(S-6)

is the mean bond rupture time. Higher moments can also be calculated with the
cumulative distribution function $G(b_{r},t|b)$ at hand. The equation for $G$
can now be used to write down coupled equations for the moments

$T_{j}(b_{r},b)=\int_{0}^{\infty}t^{n}f(b_{r},t|b)=j\int_{0}^{\infty}t^{j-1}G(b_{r},t|b)\qquad(j\geq
1,T_{0}=1).$ (S-7)

The equation for the $j$th moment reads

$\left[\frac{k_{\text{B}}T}{\zeta}\frac{\partial^{2}}{\partial
b^{2}}-\frac{1}{\zeta}\frac{\partial U(b;F)}{\partial
b}\frac{\partial}{\partial b}\right]T_{j}(b_{r},b)=-jT_{j-1}(b_{r},b)\qquad
j=1,2,..$ (S-8)

and the boundary conditions for $T_{j}(b_{r},b)$ are

$T_{j}(b_{r},b_{r})=0,\qquad\left.\frac{\partial}{\partial
b}T_{j}(b_{r},b)\right|_{b=b_{0}}=0.$ (S-9)

For $j=1$ the above Eq. (S-8) reduces to

$\left[\frac{k_{\text{B}}T}{\zeta}\frac{\partial^{2}}{\partial
b^{2}}-\frac{1}{\zeta}\frac{\partial U(b;F)}{\partial
b}\frac{\partial}{\partial b}\right]T_{1}(b_{r},b)=-1.$ (S-10)

This ordinary differential boundary problem is solved by

$T_{1}(b_{r},b)=\int_{b}^{b_{r}}dz\,\frac{1}{\Psi(z)}\int_{b_{0}}^{z}\frac{\Psi(y)}{D}\,dy,$
(S-11)

with

$\Psi(z)=\exp\left[\int^{z}\frac{-\frac{1}{\zeta}\frac{\partial
U(b;F)}{\partial
b}}{k_{\text{B}}T/\zeta}db\right]=\exp\left[-\frac{U(z;F)}{k_{\text{B}}T}\right].$
(S-12)

If we average the mean rupture time over all possible initial lengths of the
bond,

$\overline{T_{1}}(b_{r})=\frac{\int_{b_{0}}^{b_{r}}T_{1}(b_{r},b)p_{0}(b)\,db}{\int_{b_{0}}^{b_{r}}p_{0}(b)\,db},$
(S-13)

where $p_{0}(b)$ is the density distribution of the initial value. Within the
manuscript we denote the mean rupture time by $\tau(F)$ with

$\tau(F)=T_{1}(b_{r},b=b_{1}),$ (S-14)

to highlight its dependency on $F$, because $b_{r}$ is a bond type-specific
constant, because we are not considering higher moments than the first moment,
and because we choose the bond to reside at $t=0$ in its energetic minimum,
located at $b=b_{1}(F)$. Recall that our potential has the features
$U(b_{1})=U^{\prime}(b_{1})=0$ and $U^{\prime\prime}(b_{1})>0$. Moreover, we
use $b_{0}=0$ as the reflecting boundary, noting that the precise choice does
not matter as $U$ tends to diverge at $b\rightarrow 0$ due to excluded volume
interactions. We checked that the $\tau(F)$ calculated semi-analytically via
Eq. (S-11) with Eq. (S-12) (numerical integration of the double-integral) is
exactly identical with the mean rupture time obtained via Brownian dynamics of
the Langevin equation (S-1) with reflecting boundary at $b_{0}=0$, adsorbing
boundary at $b_{r}$, and initial condition $b(0)=b_{1}(F)$, if results are
extrapolated to infinitely small time step.

###  S5 Influence of boundary conditions

Explanation of fix-move, fix-force, fix-smd and fix-spring. fix-smd is the
primary method for chain stretching used within the manuscript. Only Figure
1(c) is obtained using fix-move and fix-force, while all other figures are
obtained with fix-smd.

fix-move: move the positions of 2 terminal Si atoms of a PDMS chain with
constant velocity $V$ (i.e. stretch a chain with constant velocity $V$). When
$V$ is set to 0, we can fix the extension $R$ of a PDMS chain and measure the
force $F$ to obtain force-extension relation.

fix-force: apply a constant force $F$ on a terminal Si atom while fixing the
position of the other terminal Si atom (i.e. stretch a polymer chain with
constant force $F$). After equilibrium is reached, we register $R$ to measure
force extension.

fix-smd: apply a constant repulsive force $F$ between 2 terminal Si atoms.
After equilibrium is reached, we register $R$ to measure force extension.

fix-spring: use 2 springs to stretch the 2 terminal Si atoms of a PDMS chain.
We can measure the $R$ and $F$ (forces in 2 springs) to obtain force-
extension. However due to the oscillation of springs, the measurement error is
extremely large so no results presented in this manuscript are measured with
this method.

Figure S3: Investigation of the influence of boundary conditions, temperature,
and the presence/absence of solvent molecules on the force-extension relation
of PDMS-27. We find that the choice of boundary condition does not influence
the obtained force-extension relation (fix-spring is omitted, as it is not
suitable in force-extension measurements due to large measurement errors
linked to spring oscillations). fix-move measurements at $T=$ 1K (red
triangles) show a slightly smoothed-out transition between entropic- and
enthalpic regimes, which stems from reduced entropic contributions. Figure S4:
Investigation of the influence of boundary conditions on the observed ’end
effect’. We focus on three different molecules: PDMS-2, PDMS-20, as well as a
PDMS-2 segment (Si-Si-Si) located at the center of PDMS-20. Different boundary
conditions do not influence the force-extension relation within the enthalpic
part. Furthermore, stretching an isolated PDMS-2 molecule is equivalent of
stretching PDMS-2 at the center of PDMS-20. The difference between data
plottes as $\times$ and $\square$ is the chain length between the 2 stretching
points.

###  S6 Bond length distributions and potentials

(a) (b)
(c) (d) (e)

Figure S5: Bond length distribution on (a) PDMS-4, (b) PDMS-376 and (c-e)
linked PDMS-4 obtained from atomistic MD. Simulations are performed at $T=300$
K.

(a) (b)
(c) (d) (e)

Figure S6: Effective bond potentials calculated from the bond length
distributions on (a) PDMS-4, (b) PDMS-376 and (c-e) linked PDMS-4, shown in
Fig. S1, along with the Si-O and Si-C fit functions (solid black lines) stated
in the manuscript. Figure S7: Potential of mean force and its derivatives,
$U(b)$, $U^{\prime}(b)$ and $U^{\prime\prime}(b)$ for Si-O at different levels
of relatively strong applied force $F>3$ nN (MD simulations performed in the
presence of solvent molecules). Since $U^{\prime\prime}$ exhibits parabolic
and $F$-independent shape and location, the corresponding parameters $b_{2}$,
$c_{2}$, $k_{2}$ in the 4th order polynomial are treated as $F$-independent
constants.

###  S7 Rupture times

Here, we provide evidence that the measured rupture times are basically
unaffected by the presence of HMDSO solvent molecules. In the absence of
solvent, the friction coefficient $\zeta$ is implicitly captured by the
employed thermostat. This finding allows us to simulate the exponential tail
of the rupture time distribution in the absence of solvent (Fig. 3a). This
renders computatiosn feasible (as it is two orders of magnitude cheaper than
the full atomistic simulation of the solvated PDMS chain). Shown in Fig. S8 is
the cumulative fraction of ruptured PDMS-6 chains versus time both in the
presence and absence of solvent molecules.

Figure S8: Cumulative rupture time distribution for PDMS-6 at $F=6.39$ nN.
Red: Simulation in the presence of HMDSO solvent molecules. Black: Simulations
without solvent. Results obtained by averaging over 100 (with) and 10000
(without solvent) independent start configurations. Simulations are performed
at $T=300$ K.

###  S8 Mean chain rupture times for different polymerization degrees

Figure S9: Equivalent of Figure 3c for different polymerization degrees $n$.
At low $n$ (blue dashed line), crosslinking junctions are weaker than Si-O
bonds along the backbone (corresponding to the dashed line at position (A) in
Figure 4c. A double crossover between Si-C and Si-O bonds is observed at high
$n$ (dashed red line). This corresponds the re-entrant effect observed in
Figure 4c.