\section{Supplemental Material for \\ ``Athermal Fracture of Elastic Networks:\\ How Rigidity Challenges the Unavoidable Size-Induced Brittleness''}

\renewcommand{\thefigure}{S\arabic{figure}}
\noindent

\section{Simulation methods}
We generate networks by diluting regular lattices consisting of $L \times L$ nodes with spacing $l_0$ and different topologies. A fraction $p$ of all possible nearest-neighbor bonds is present. The associated average network connectivity is therefore $\langle z \rangle = p \, z_{\text{max}}$, with $z_{\text{max}}$ depending on network topology (see also Fig.~\ref{SfigB}). All bonds are harmonic springs with stiffness $\mu$ and rest length $l_0$ and break irreversibly when exceeding the threshold $\lambda$ (same for all springs). By setting $l_0=1$ and $\mu=1$, we can express all quantities in reduced units, and $\lambda$ can be equivalently considered as deformation or force threshold. The energy of the system reads $\mathcal{H}=\frac{\mu}{2l_0}\sum_{\langle ij \rangle} {(l_{ij}-l_0)}^2$, where $\langle ij \rangle$ indicates bonded pairs of nodes and $l_{ij}$ the distance between them. We deform the networks uniaxially under athermal and quasistatic loading conditions, using the FIRE algorithm for the energy minimization~\cite{bitzek2006}. The tolerance $F_{\text{RMS}}$ is typically set at $10^{-5}$ for $L \leq 256$, but it is varied for larger sizes as discussed in the text. In our case, $F_{\text{RMS}}$ is the maximum numerical tolerance for both the system root-mean-square force~\cite{bitzek2006} and the maximum force on each single node. The networks are deformed uniaxially in the $y$-direction by employing a very small strain increment $\Delta \epsilon=0.001$, after which the energy minimization procedure is performed to obtain the equilibrium (off-lattice) node positions. Based on the mechanically equilibrated configuration, we identify the bonds that exceed $\lambda$ (if any) and remove the weakest bond defined as the one that exceeds the threshold the most, i.e. the one with the largest $(l_{ij}-l_0)/\lambda$. In most of the simulations, we additionally allow the strain increment to be reduced by up to two order of magnitude ($\Delta \epsilon=0.00001$) if the most stressed bond is very close to $\lambda$, to ensure that we are in the quasistatic condition and only one bond at a time can be broken. Typically, we use periodic boundary conditions in all directions. However, following Ref.~\cite{zhang2017}, we also perform simulations in which the nodes on the top and bottom boundaries are free to move only in the $x$-direction and move affinely in the $y$-direction, without noticing any appreciable qualitative difference. Furthermore, as detailed below, we also perform simulations by employing different boundary conditions in the $x$-direction. We confirmed that the conclusions of our study are robust with respect to all the variations on the simulation protocol.


\section{Observables and data analysis}
We define the mechanical response of the network by considering the stress along the direction of the deformation. In particular, we use the $yy$-component of the virial stress tensor $\sigma \equiv \sigma_{yy}^{\text{vir}} = \frac{1}{d V} \sum_{\langle ij \rangle} f_{ij,y} \, r_{ij,y}$ , where the sum runs over all the bonded pairs of nodes $\langle ij \rangle$, $\mathbf{f}_{ij}$ is the force acting on node $i$ due to $j$, $\mathbf{r}_{ij}$ is the vector connecting the two nodes, $d$ is the dimensionality of the network and $V$ is its (instantaneous) volume. Using a different definition for the macroscopic stress (e.g. $\sigma_{xx}^{\text{vir}}+\sigma_{yy}^{\text{vir}}$) does not change the conclusions of the paper. A network is defined rigid if $\sigma$ is greater than a small value ($10^{-4}$) and $\epsilon_{\text{R}}$ is the first strain for which this condition is satisfied. $\sigma_{\text{c}}$ is the maximum in the stress response and $\epsilon_{\text{p}}$ the associated strain. The work of extension $W_{\text{f}}$ is calculated as the integral of the stress-strain curve until fracture (network broken in two parts). $\Delta \sigma_{\text{max}}$ is the maximum stress drop and it is calculated following a procedure highlighted in Ref.~\cite{bonfanti2018}. It consists of (i) calculating the derivative of the stress-strain curve; (ii) identifying the strain interval for which the derivative is greater than a certain threshold (here we use zero, but we checked that our results are robust using different thresholds); (iii) for each strain interval, calculate the stress drop as the difference in the stress associated to the extremal values of the interval; (iv) the maximum stress drop is $\Delta \sigma_{\text{max}}$. In most of the cases, this procedure gives exactly the same result as just considering the difference in stress between two consecutive strains; whereas in few cases the values are slightly different but the overall trends are consistent with both definitions.

All quantities are averaged over many configurations to ensure proper statistical sampling of the disorder. Errors are calculated as standard deviation divided by the number of sampled configurations (standard error on the mean) and are shown when they are larger than the symbols displayed in the graphs. For example, for the size scaling analysis for networks close to the isostatic point ($p=0.65$ and $p=0.70$) and for small $\lambda$ we employ more than 2000 configurations for small system sizes ($L \leq 32$); 1000 for $L=50$; 500 for $L=64$; 200 for $L=128,256$; 50 for $L=512$; and $\sim 10$ configurations for larger system sizes ($L=640,768,1024$) that are statistically sufficient due to the self-averaging nature of the fracture process in such large systems. In general, for all networks we average quantities over at least 200 configurations when $L=128$ and a larger number (typically 500) for $L=64$ or smaller sizes. Results for different network geometries are also averaged over approximately 200 configurations for 2D networks and over 20 to 100 configurations for 3D networks, depending on system size.

\section{Additional results}

\subsection{Fixed system size}
In Fig.~\ref{SfigA}, some additional results for triangular networks with fixed system size $L=128$ and different $p$ and $\lambda$ are shown. In particular, we note in Fig.~\ref{SfigA}(a) that for fixed $\lambda$ the work of extension $W_f$, i.e. the energy needed to break the material (calculated as the integral of the stress-strain curve), shows a maximum around the isostatic point ($p_{iso}\simeq0.66$) for small $\lambda$, analogously to the trend for $n_f$ (shown in Fig.~1(c) of the main text). This has been observed also in Ref.~\cite{zhang2017} and implies that networks close to the isostatic point can be counterintuitively considered tougher than networks with more bonds (when $\lambda \leq 0.10$). However, as shown in Fig.~\ref{SfigA}(b), when considering the maximum stress $\sigma_{\text{c}}$ as a function of connectivity we observe that there is no maximum and the more intuitive monotonic dependence on $p$ --- the more the bonds, the stronger the material --- is observed. This discrepancy between the dependence on connectivity of two quantities describing resistance to fracture is an effect of the size-dependent brittleness demonstrated in this study: at moderate system size, the networks are ductile and significant energy is dissipated in bond-breaking events occurring after the peak stress. To further quantify the tunability of the fracture process by varying $p$ and $\lambda$ at fixed system size, we show the peak strain $\epsilon_p$ in panel (c) and the maximum stress drop in panel (d).

\subsection{Different topologies}
In Fig.~\ref{SfigB}, we show $n_{\text{f}}$ as a function of the network average connectivity $\langle z \rangle=pz_{\text{max}}$ for different network topology and dimensionality, and small $\lambda$. We observe that the position of the maximum in $n_{\text{f}}$ depends not only on dimensionality but also on topology and orientation. For example, despite the fact that the diluted square (SQ) and the 2D diamond (D2) networks have the same connectivity properties, SQ is more brittle (fewer bonds are broken) because it features more bonds already aligned with the direction of deformation. Analogously, because of geometric effects (a large number of system-spanning bonds at small connectivity) in the honeycomb (HC) and in the 3D diamond (D3) networks, a larger number of bonds is broken when approaching the limit of undiluted lattices. Only for the diluted triangular (TR) and face-centred-cubic (FCC) networks, that are sufficiently isotropic, the maximum of $n_{\text{f}}$ roughly corresponds to the isostatic point $z_{\text{iso}}$. Our results enrich the data reported in Refs.~\cite{zhang2017,driscoll2016,berthier2019}. We can summarize that damage is maximized close to rigidity points, either of isostatic or geometric nature, when $\lambda$ is small.

\subsection{Size scaling}
In Fig.~\ref{SfigC}, additional results from the size scaling are shown. To extract the damage fractal dimension $d_{\text{f}}$ and generate the plot of Fig.~2(b) of the main text, we simulate at least six different system sizes (typically seven, up to $L=256$) and we fit $N_{\text{f}} \sim N_{\text{in}}^{d_{\text{f}}/2}$ where $N_{\text{in}}$ is the number of intact initial springs. This is analogous to $N_{\text{f}} \sim L^{d_{\text{f}}}$ but sometimes better fits are obtained with the previous expression since the values on the x-axis span more orders of magnitude. As shown in the main text, we notice that the fractal dimension $d_{\text{f}}$ is (much) greater than 1 for a large range of $p$ and $\lambda$. This is a distinct feature of our elastic network model compared to the scalar random fuse models~\cite{alava2006} or models assuming linear elasticity~\cite{nukala2005}, where the number of broken bonds at peak load scales \emph{at most} linearly with the system size. In our case $d_f$ refers to the exponent of the total broken bonds at failure and not to the exponent of broken bonds at peak load. However, since the fracture becomes more brittle (i.e. bonds break only before and at peak load) the exponent for the broken bonds at peak load is actually larger than $d_f$. We also show how the maximum stress $\sigma_{\text{c}}$ and peak strain $\epsilon_p$ decrease with increasing system size. Inspired by previous works~\cite{duxbury1987,manzato2012,roy2017b}, we fitted the decay of $\sigma_{\text{c}}$ with three functional forms: (i) logarithmic decay $\propto 1/\ln(L)$; (ii) Duxbury-Leath-Bale (DLB)~\cite{duxbury1987,manzato2012} $\propto 1/\sqrt{\ln(L^2)}$; (iii) power-law decay $\propto 1/L^{\beta}$. We found that the power-law decay fits our data the best, with an exponent $\beta$ that clearly depends on $p$ and $\lambda$. For example, we obtain $\beta=0.73$ for $p=0.65$, $\lambda=0.03$ (Fig.~2(d) of main text) and $\beta=0.89$ for $p=0.70$, $\lambda=0.03$ (Fig.~\ref{SfigC}(d)). Differently from previous models where it is assumed that the material strength (i.e. $\sigma_c$) is zero in the thermodynamic limit~\cite{alava2006}, in our case we have an offset value corresponding to a finite network strength $\sigma_c^{(\infty)} \neq 0$ that depends on $p$ and $\lambda$. Our data suggest that the spring network model studied here is qualitatively different from scalar models that have been so far considered as the golden standard for statistical physics studies of fracture.

\subsection{Universality of fracture abruptness trend}
In Figs.~\ref{SfigD} and~\ref{SfigE}, we show that the non-monotonic size-dependence of fracture abruptness $\Delta \sigma_{\text{max}}$ is universal with respect to network topology and robust with respect to the simulation procedure. In Fig.~\ref{SfigD}(a), we compare the trend for the same set of initial configurations simulated with a protocol where the small strain increment is kept fixed at $\Delta \epsilon=0.001$ and one where it can be adjusted to much smaller values (up to $\Delta \epsilon=0.00001$) if the most stressed bond is very close to $\lambda$. A minimum is observed in both cases for similar values of $L$. A slightly smaller $\Delta \sigma_{\text{max}}$ is observed for large $L$ when $\Delta \epsilon$ can be reduced by two orders of magnitude. Next, we check if the choice of boundary conditions in the $x$-direction, transverse to the deformation direction, has some qualitative effect on the fracture abruptness. We tested three different boundary conditions: (i) periodic (mimicking an infinitely large system); (ii) free, i.e. without bonds crossing the $x$-direction (the system is finite in that direction); (iii) with periodic bonds and with adjustable size in the $x$-direction such that the area is kept constant during the uniaxial deformation. As shown in Fig.~\ref{SfigD}(b), the qualitative trends are robust and independent of the fundamentally different boundary conditions employed (the first one was chosen for the rest of our study since it seems the most computationally efficient). Furthermore, we confirm that the non-monotonic trend persists when changing network topology. As already highlighted in the main text, the minimum and the maximum of $\Delta \sigma_{\text{max}}$ as a function of $L$ depend on network properties. This can be appreciated in Fig.~\ref{SfigD}(c) where we report selected cases for various 2D and 3D networks. Finally, preliminary results confirm the non-monotonic trend also for shear deformation, supporting the conclusion that it is a universal size-dependent behavior of diluted networks. As shown in Fig.~\ref{SfigE}, we also succeeded in manually collapsing the curves of $\Delta \sigma_{\text{max}}$ obtained for networks close to (just below and slightly above) the isostatic point and various $\lambda$ into a single master curve. The rescaling $L \rightarrow \alpha L$ and $\Delta \sigma_{\text{max}} \rightarrow \beta \Delta \sigma_{\text{max}}$ was performed by imposing that the different curves approximately share the same x-value for the minimum. We obtained values for $\alpha$ that are inversely proportional to $\lambda$, suggesting that the ``critical'' system size $L^* \sim 1/\alpha$ at which brittleness kicks in decreases by making the springs stronger. On the other hand, $\beta$ increases more than linearly with $\lambda$, consistently with the idea that due to inhomogeneous stress distribution, the increase of the springs strength has non-linear consequences on the overall network.

\subsection{Effect of tolerance in energy minimization}
Finally, in Fig.~\ref{SfigF} we report the size scaling of the ductility interval $\Delta \epsilon_{\text{duct}}$, defined as the strain difference between $\epsilon_p$ and the strain corresponding to the maximum stress drop, for two different values of the tolerance $F_{\text{RMS}}$ and various networks. Brittle fracture (vanishing $\Delta \epsilon_{\text{duct}}$) is observed in the athermal limit, whereas a non-monotonic size-dependence of the ductility interval is observed for larger $F_{\text{RMS}}$. The system size at which ductility is increasing again (the minimum in $\Delta \epsilon_{\text{duct}}$ as a function of $L$) does strongly depend on $\lambda$ and therefore seems to be coupled with the network rigidity. If a variation of the tolerance in the energy minimization protocol $F_{RMS}$ is interpreted as a variation of the thermal effects in the system (that would keep the system not exactly at the energy minimum), these results might suggest that by increasing temperature it is possible to observe a more ductile fracture occurring via many large cracks. Future research should further investigate the intriguing possibility that the minimum temperature at which such ductile behavior is first observed, and the material rigidity are strongly coupled. However, to precisely assess the role of temperature, more traditional simulation methods for thermal systems should be employed.

\clearpage

\begin{figure*}[h!t]
	\begin{center}
		\includegraphics[width=0.9\textwidth]{./SfigA.pdf}
	\end{center}
	\caption{Additional results for triangular networks with $L=128$ as a function of $p$ and various $\lambda$ (see legend in panel (b)). (a) Work of extension $W_{\text{f}}$, defined as the integral of the stress-strain curve. (b) Maximum stress $\sigma_{\text{c}}$. (c) Strain $\epsilon_{\text{p}}$ corresponding to the stress peak as a function of $p$ for various $\lambda$. Inset: strain $\epsilon_{\text{R}}$ at which the network is rigid upon uniaxial extension as a function of $p$. (d) Maximum stress drop $\Delta \sigma_{\text{max}}$ versus $p$ for different $\lambda$.  }
	\label{SfigA}
\end{figure*}

\begin{figure*}[h!t]
	\begin{center}
		\includegraphics[width=0.75\textwidth]{./SfigB.pdf}
	\end{center}
	\caption{Fraction of broken bonds as a function of averaged network connectivity $\langle z \rangle$ for (a) 2D and (b) 3D networks with the different architectures shown on the right, same $\lambda$ and similar system size. Maxima are observed close to rigidity points that depend on topology, geometry and dimensionality.}
	\label{SfigB}
\end{figure*}


\begin{figure*}[h!t]
	\begin{center}
		\includegraphics[width=0.9\textwidth]{./SfigC.pdf}
	\end{center}
	\caption{(a) Number of broken bonds at failure $N_{\text{f}}$ as a function of number of initial intact springs $N_{\text{in}}$ for various combination of $p$ and $\lambda$. Dashed lines are best power-law fits from which the damage fractal dimension $d_{\text{f}}$ is extracted. (b) Size scaling of broken bonds for triangular network close from above to the isostatic point with small $\lambda$. Black circles indicate single-run values, red squares averages and dashed line best power-law fit. (c) Size scaling of strain at peak stress $\epsilon_{\text{p}}$ for networks in the proximity of the isostatic point with small $\lambda=0.03$. Symbols represent averages and lines are best fits. Data are well fitted by a power-law decay with an offset $\epsilon_{\text{p}} = \epsilon_{\text{p}}^{(\infty)} + A/L^{\beta}$, with $\epsilon_{\text{p}}^{(\infty)}$ the strain at peak in the thermodynamic limit. For $p=0.65$, we obtain $\epsilon_{\text{p}}^{(\infty)}=0.050$ and $\beta=0.47$. For $p=0.70$, we obtain $\epsilon_{\text{p}}^{(\infty)}=0.016$ and $\beta=0.24$. (d) Size scaling of the network strength $\sigma_{\text{c}}$ (log-log plot). Symbols correspond to averages, error bars show the standard deviations and lines are the best fits. Three functional forms were used: logarithmic decay (dotted black line), DLB decay (dashed red line), and power-law decay $\sigma_{\text{c}} = \sigma_{\text{c}}^{(\infty)} + A/L^{\beta}$ (green solid line). The latter fits the data the best, we obtain $\sigma_{\text{c}}^{(\infty)}=0.0053$ and $\beta=0.89$. }
	\label{SfigC}
\end{figure*}

\begin{figure*}[h!t]
	\begin{center}
		\includegraphics[width=0.8\textwidth]{./SfigD.pdf}
	\end{center}
	\caption{(a) Size scaling of maximum stress drop for two different simulation protocols. Black circles represent data from simulations where the strain step can be adjusted (reduced) to smaller steps when there are bonds very close to their thresholds (see section on simulation methods); whereas blue squares are data from simulations where the strain step is kept fixed at a small value. (b) Size scaling of maximum stress drop for four networks with $p$ and $\lambda$ as indicated in the legend, and three different simulation protocols for which periodic boundary conditions in the $x$ direction (normal to the deformation direction) are present (solid lines) or not (dashed lines), or for which the $x$ direction is adjusted to keep the system area constant (dotted lines). (c) Size scaling of fracture abruptness for diluted 2D honeycomb (HC), 2D square (SQ) and 3D FCC (inset) networks with averaged connectivity $\langle z \rangle$ and $\lambda$ as indicated.}
	\label{SfigD}
\end{figure*}

\begin{figure*}[h!t]
	\begin{center}
		\includegraphics[width=0.9\textwidth]{./SfigE.pdf}
	\end{center}
	\caption{Rescaled maximum stress drop $\Delta \sigma_{\text{max}}$ as a function of rescaled system size $L$ for (a) $p=0.65$ (see also Fig.~3(d) of main text) and (b) $p=0.70$, color-coded according to $\lambda$. Symbols represent averages and error bars are standard deviations. (c) Both data sets combined (circles for $p=0.65$, squares for $p=0.70$). The values used for the rescaling parameters $\alpha$ and $\beta$ are shown in the insets as a function of $\lambda$. $\beta$ grows (more than linearly) with increasing $\lambda$. $\alpha$ is related to the transition system size $L^*$ discussed in the main text and decreases by two order of magnitudes in the simulated range of $\lambda$. }
	\label{SfigE}
\end{figure*}

\begin{figure*}[h!t]
	\begin{center}
		\includegraphics[width=0.9\textwidth]{./SfigF.pdf}
	\end{center}
	\caption{Size scaling of the ductility interval $\Delta \epsilon_{\text{duct}} = \epsilon_{\Delta\sigma_{\text{max}}} - \epsilon_{\text{p}}$, where $\epsilon_{\Delta\sigma_{\text{max}}}$ is the strain at which the maximum stress drop is observed (see also Fig.~4 of the main text), for simulations with (a) $F_{\text{RMS}}=10^{-5}$ and (b) $F_{\text{RMS}}=10^{-6}$ for different $p$, $\lambda$ as indicated in the legend. For small $F_{\text{RMS}}$ (panel (b)), it is evident that fracture becomes more brittle upon increasing system size. For larger $F_{\text{RMS}}$ (panel (a)), that are nevertheless three orders of magnitude smaller than $\lambda$, a non-monotonic trend is observed and networks with large $L$ reacquire ductility. At a first glance, the minima in $\Delta \epsilon_{\text{duct}}$ seem to roughly correspond with the minima in $\Delta \sigma_{\text{max}}$ (see Fig.~3(a) of main text), but further studies are needed to confirm this correspondence. }
	\label{SfigF}
\end{figure*}

\end{document}