Moduli and modes in the Mikado model

Abstract

We determine how low frequency vibrational modes control the elastic shear modulus of Mikado networks, a minimal mechanical model for semi-flexible fiber networks. From prior work it is known that when the fiber bending modulus is sufficiently small, (i) the shear modulus of 2D Mikado networks scales as a power law in the fiber line density, G ∼ ρ^α+1, and (ii) the networks also possess an anomalous abundance of soft (low-frequency) vibrational modes with a characteristic frequency ω_κ ∼ ρ^β/2. While it has been suggested that α and β are identical, the preponderance of evidence indicates that α is larger than theoretical predictions for β. We resolve this inconsistency by measuring the vibrational density of states in Mikado networks for the first time. Supported by these results, we then demonstrate analytically that α = β + 1. In so doing, we uncover new insights into the coupling between soft modes and shear, as well as the origin of the crossover from bending- to stretching-dominated response.

---

Mikado networks¹ are a widely studied off-lattice model for the mechanics of semiflexible fiber networks.^2–9 Two-dimensional Mikado networks consist of slender rods with stretching modulus μ and bending modulus κ, which are randomly deposited in a box and crosslinked at each intersection point. At a phenomenological level, the model successfully reproduces important features of real fiber networks, including the affine to non-affine crossover in linear response (discussed below) and the onset of strain stiffening at a critical strain scale.^10–13 The model also serves as a bridge to other soft matter systems where floppiness and non-affinity play an important role, such as solid foams and confluent tissues,^14–16 soft sphere packings,^17–22 and random spring networks.^23–28

The goal of this article is to determine how low frequency (“soft”) vibrational modes control the shear modulus G in 2D Mikado networks. It is natural to expect the displacements of crosslinks during shear deformation to project strongly on soft modes, which implies a relationship to the modulus. This relationship has the potential to explain the strong non-affine fluctuations that are known to emerge in fiber networks, by connecting them to the floppy (zero frequency) vibrational modes of freely bending (κ = 0) networks, which are inherently non-affine. However, precisely how soft modes set the modulus remains unclear, as there is an apparent contradiction between existing results in the literature. While both G^1,3,4,7 and the typical soft mode frequency ω_κ^5,6 display anomalous power law scaling with the fiber line density ρ (defined below), the best estimates of their respective exponents do not appear to be compatible. In order to describe this discrepancy more precisely, we must first summarize what is known about moduli and modes separately.

Shear modulus. We begin with the shear modulus of bi-periodic networks of N fibers in an area A, as depicted in Fig. 1a. There are three independent, microscopic length scales that characterize such a network. The first is the fiber length [small script l]_f. The second is the mean length of a fiber segment between two crosslinkers. A Poissonian fiber deposition process results in a mean segment length [small script l]_s = π/(2ρ), where the fiber line density ρ ≡ N[small script l]_f/A is the total length of fiber per unit area. The third length [small script l]_b ≡ (κ/μ)^1/2 is mechanical; it represents the natural length over which a free fiber bends. When mechanical properties of networks are expressed in dimensionless form, they must be functions of two dimensionless ratios constructed from these three length scales. We choose [small script l]_b/[small script l]_f and ρ[small script l]_f.

Fig. 1 A Mikado network with ρ[small script l]_f = 40.0.

The most relevant features of G can be illustrated by varying the bending modulus κ while holding all other model parameters fixed (i.e., varying [small script l]_b/[small script l]_f at fixed ρ[small script l]_f), as shown using our own data in Fig. 2a (details below). If we look first at larger values of κ, then eventually the energetic penalty for bending will be so large that the fibers essentially do not bend. Numerics indicate that when deformation is stretching-dominated, the displacements of crosslinks follow the global deformation affinely.^1,3 The shear modulus is therefore given by its affine value G ≃ G_aff, which can be calculated exactly⁴ and scales as G_aff ∼ μρ/[small script l]_f in densely crosslinked networks. If instead we look at small κ, bending eventually becomes energetically favorable to stretching, the response becomes strongly non-affine, and G becomes much smaller than G_aff. Hence the shear modulus displays a crossover between stretching- and bending-dominated deformation, which correlate with affine and non-affine response, respectively. The bending-dominated regime is of particular interest, both because affine deformation is comparatively simple to model, and because non-affinity represents strong deviations from the mean response.

Fig. 2 (a) The shear modulus G of a Mikado network plotted versus ([small script l]_b/[small script l]_f)² for varying ρ[small script l]_f indicated in the legend (see text for definitions). The solid line denotes the affine shear modulus for ρ[small script l]_f = 48.9. (b) The same data as in (a), now collapsed by plotting G/G_aff as a function of ([small script l]_b/[small script l]_f)²(ρ[small script l]_f)^α with α = 7.2. The dashed line has a slope of 1, and the solid line denotes G = G_aff.

It is evident from Fig. 2a that the value ([small script l]_b/[small script l]_f)_* where the bending- to stretching-dominated crossover occurs depends itself on ρ[small script l]_f. If we assume power law scaling ([small script l]_b/[small script l]_f)_* ∼ (ρ[small script l]_f)^−α/2 for some undetermined exponent α, it is then natural to expect that the dimensionless shear modulus [capital G, script] ≡ G/G_aff can be written as a function of the variable z = ([small script l]_b/[small script l]_f)²(ρ[small script l]_f)^α. The scaling function has the form

where z₀ is a constant. The approach to [capital G, script] ≃ 1 for large z ensures that G ≃ G_aff is recovered, while the form [capital G, script] ∼ z for small z is necessary so that G is proportional to κ and independent of μ in the bending-dominated regime. The value of α has been estimated numerically by several groups.† Wilhelm and Frey found α ≈ 5.7.¹ Head et al.^3,4 estimated α ≈ 8 for networks with ρ[small script l]_f as large as 120 (comparable to Wilhelm and Frey), but suggested that α will approach a value of 7 when larger values of ρ[small script l]_f are accessed. Shahsavari and Picu⁷ simulated ρ[small script l]_f as large as 300 and found α ≈ 7. Head et al.⁴ also presented a mean field scaling argument that predicts α = 7; however, it does not invoke modes and so does not clarify their relationship to the modulus. As shown in Fig. 2b and described in greater detail below, we obtain good data collapse with our own numerical results for α ≈ 7.2.‡

For completeness, we note that Mikado networks undergo a rigidity percolation transition at a critical fiber line density ρ_c ≈ 5.7.^1,2 Like the prior studies discussed above, we estimate α using data well above ρ_c. G crosses over to a qualitatively different scaling regime as ρ → ρ_c from above. This regime is not the focus of the present study.

Soft modes. Let us now turn to soft vibrational modes in Mikado networks. If fibers are allowed to bend freely (κ = 0), then Mikado networks always possess so-called floppy modes.^2,5,8 These are non-rigid body motions that can be performed without doing work, and therefore correspond to vibrational modes with zero frequency. A finite bending stiffness stabilizes floppy motions, lifting them to finite frequency. As floppy motions are intrinsically non-affine, the number and vibrational frequencies of soft modes must be fundamentally related to non-affine fluctuations in mechanical response.

In two seminal papers, Heussinger and Frey (HF) developed what has come to be known as floppy mode theory (FMT).^5,6 The theory (i) explicitly constructs the floppy modes in freely bending Mikado networks, and (ii) uses these modes to build a set of trial modes, and thereby predict the characteristic soft mode frequency ω_κ at small but finite bending modulus using the variational method. Each trial mode is constructed by displacing one of the N fibers a distance δx along its axis. This motion is opposed by the rest of the network, and therefore carries an energetic cost W(δx) = (1/2)K_κδx². The average effective “spring constant” K_κ sets the characteristic frequency, ω_κ ∼ K_κ^1/2. By using the crosslink displacements from a floppy mode and the known distribution of segment sizes, HF found

where ω₀ = (μ/[small script l]_f)^1/2. Depending on the specific approximations that are made, FMT predicts β = 6 or β ≈ 5.75.

Rather than directly measuring the frequency ω_κ, HF tested their theory indirectly by measuring the shear modulus. To predict G, they equated the energy of the shear deformation, GAγ², with the energy NW(δx) stored by the soft modes under shear. However, it is not apparent what the typical axial displacement δx should be. In order to close this gap, they further assumed that the center of mass of each fiber displaces affinely. We shall refer to this as the affine fiber approximation (AFA). It implies that each fiber undergoes an axial displacement δx ∼ [small script l]_fγ relative to its environment. The shear modulus is then

Taken together, eqn (1)–(3) imply α = β. Based on this relation, Heussinger and Frey concluded that the numerical result α = 5.7 from Wilhelm and Frey¹ validates their theoretical predictions for β. However, as noted above, independent studies^3,4,7 conclude instead that α ≈ 7. In our opinion, the preponderance of evidence indicates that α is not equal to the FMT prediction for β.

Hence we have a puzzle. The discrepancy between α and the theoretically predicted value of β indicates that there is a gap in our understanding of non-affine response in the Mikado model – and, by extension, in fiber networks in general. Where is the missing physics? The persistence of this gap is due, at least in part, to the paucity of data for the vibrational density of states (and hence ω_κ) in Mikado networks and related models.²⁹ Since ω_κ has not been measured directly, one possibility is that β is larger than the HF theory predicts. We test this by presenting the first direct measurement of ω_κ in the Mikado model, as shown in Fig. 3a (details below). We obtain the best collapse of the data for β = 6.1, with reasonable collapse in the range 5.8 ≤ β ≤ 6.5. This result for β is in reasonable agreement with FMT, and also too small to account for the discrepancy between α and β.

Fig. 3 (a) The characteristic frequency ω_κ of soft vibrational modes in a Mikado network plotted versus ([small script l]_b/[small script l]_f)²(ρ[small script l]_f)^β with β = 6.1. See Fig. 2 for legend. The dashed line has a slope of 1. (b) The ratio of the dimensionless shear modulus shear modulus G/G_aff to the characteristic frequency ω_κ is not constant when plotted versus fiber line density, which implies α ≠ β.

The second possible explanation for the discrepancy, which we favor, is that the affine fiber approximation is incorrect, and hence that the relation α = β is unjustified. This view is supported by Fig. 3b, which plots the ratio (G/G_aff)/(ω_κ²/ω₀²). eqn (3) predicts this ratio to be independent of the fiber line density, but instead we find a roughly linear trend. In the following sections, we will present evidence that in fact

α = β + 1. (4)

This relation is in reasonable agreement with our numerical results in Fig. 2 and 3. More importantly, we will derive (4) in two different ways, neither of which invokes the AFA. This resolves the apparent discrepancy between moduli and modes in the Mikado model.

1 Shear response in the Mikado model

We now define the Mikado model in greater detail, and demonstrate how to calculate the shear modulus from the crosslink displacements in response to an imposed shear stress.

Mikado networks consist of N filaments of equal length [small script l]_f, which are deposited randomly in a biperiodic square system of linear size L and area A = L². Their positions and orientations are drawn from uniform probability distributions. If, during the deposition process, two filaments intersect, they form a crosslink around which the filaments can freely rotate. The crosslinks cannot be broken. In the initial condition, i.e. after deposition and prior to shear, the fibers are straight lines subdivided into co-linear segments by the crosslinks. We trim so-called “dangling ends,” i.e. segments at the end of a fiber that are attached via a single crosslink, because they neither stretch nor bend. We choose units of length and energy such that [small script l]_f = 1 and μ = 1 and fix A = 9[small script l]_f². Nevertheless, we continue to write [small script l]_f and κ in scaling relations for transparency.

The energy of a network configuation is given by

Here 〈ij〉 denotes the segment connecting crosslinks i and j, and 〈ijk〉 labels two consecutive segments 〈ij〉 and 〈jk〉 on the same fiber. The quantity δ[small script l]_〈ij〉 is the change in the length of segment 〈ij〉 from its value in the initial condition. The angle θ_〈ijk〉 measures the deviation from co-linearity of adjacent segments 〈ij〉 and 〈jk〉, and [small script l]_〈ijk〉 is the average of [small script l]_〈ij〉 and [small script l]_〈jk〉.

By construction, the network is initially in a stress-free state. When subjected to a shear stress σ, the box undergoes a shear strain γ and the crosslinks displace from their initial positions. We consider a simple shear in which the lattice vectors of the initially square unit cell become

where F̂ is the deformation gradient. The crosslink displacements can be expressed as a sum of affine and non-affine terms. Denoting the position of crosslink i as X⃑_i, its affine displacement is A⃑_i = F̂X⃑_i. Collecting the 2N_× Cartesian components of the crosslink positions in a vector |X〉 (and similarly for other quantities), the crosslink positions after shearing are |X〉 + |A〉 + |U〉, where |U〉 represents the non-affine part of the displacements. While each element of |U〉 is an independent degree of freedom, specifying the shear strain selects the value of |A〉 on each crosslink. Hence a shear displacement of the network is described by 2N_× + 1 degrees of freedom in total.

The response to shear can be determined by expanding the energy (5) in the components of |U〉 and γ. In the harmonic approximation, the result is

Here Ĥ is the “standard” Hessian, i.e. the 2N_× × 2N_× matrix of second derivatives of the energy with respect to each Cartesian component U_n of the crosslink displacements |U〉,Similarly, the vector |Ξ〉 is defined viaIt represents the force imbalance on each crosslink that results from affinely deforming the entire system (unit cell plus crosslinks). Finally, the affine modulus isWe refer to the entire matrix on the lefthand side of eqn (7) as the extended Hessian Ĥ^ext.

The shear modulus can be determined numerically by selecting a value of the shear stress σ and inverting eqn (7) numerically to find |U,γ〉. The modulus can then be read off from G = σ/γ. The results of this procedure, which were already discussed in the Introduction, are plotted in Fig. 2 after averaging over 100 samples per condition. To estimate the exponent α, we introduce a trial value α_trial and select the n data points (z_i,[capital G, script]_i) satisfying z(α_trial) < 0.1. According to eqn (1), the variable g_i = [capital G, script]_i/z_i should approach a constant value, independent of z_i. Deviations from this expectation are quantified by the cost function , which quantifies the relative amplitude of fluctuations about the mean ḡ. At α_trial = 7.2, the function obtains its minimum value E_min = 8.2 × 10⁻³, meaning data collapse to within less than 1% of the average. E remains within 50% of its minimum for 6.9 ≤ α_trial ≤ 7.6.

2 Density of states and soft modes

To make further progress, we now follow the approach developed in ref. 20,27,30, and consider the eigenmodes of Ĥ^ext. These serve as a convenient eigenbasis in which to express the quasistatic displacements of a network's nodes. Eigenmodes of the Hessian are also eigenmodes of the dynamical matrix for the case where each crosslink is assigned unit mass, which motivates our choice to describe them as (undamped) “vibrational modes.” Nevertheless, eigenvalues of the Hessian have units of stiffness (force/length), rather than (frequency)². The same Hessian can be used to model, e.g., overdamped relaxations in the presence of a solute, see ref. 20,27.

Each mode |U_n,Λ_n〉 satisfies the eigenvalue equation

Ĥ^ext|U_n,Λ_n〉 = ω_n²|U_n,Λ_n〉, (11)

where ω_n is the eigenfrequency of mode n. We then express the shear response |U,γ〉 as a superposition of the modes,To better understand the quantity Λ_n, it is useful to note that the shear strain isFor two modes with equal coefficients c, the mode with the larger Λ will contribute more strongly to the strain. As the sign of Λ can fluctuate, its square can be understood as a measure of how strongly the mode couples to shear.

We now solve for G = σ/γ by inserting eqn (12) in (7) and exploiting the orthonormality of modes to find

The restriction to positive frequencies excludes trivial rigid body modes. We will evaluate (14) for G in Section 4. But first, we focus on the statistics of modes.

It is useful to rewrite eqn (14) in the continuum limit. To do this we assume that Λ_n²N² is independent of N (which will be justified below) and define Λ(ω)/N² to be the average of Λ_n² in the interval [ω,ω + dω]). Similarly, we introduce the density of states D(ω), which is the probability density of ω in the same interval. We then have

The Mikado density of states D(ω) is plotted in Fig. 4 for one value of ρ and varying κ. It features two distinct bands. The lower band shifts with varying κ, while the upper band remains unchanged. We conclude that the lower band contains bending-dominated modes, and the upper band contains stretching-dominated modes.

Fig. 4 The Mikado density of states for ρ[small script l]_f = 20.0 and κ = 10⁻¹⁴…10⁻⁸ in decade steps. The frequencies ω_κ and ω_μ are indicated with arrows for the lowest value of κ. The dashed line has slope − 3.

The bending-dominated band is characterized by its peak frequency ω_κ, which we identify with the frequency scale studied by HF. We have estimated ω_κ by fitting log D to a parabola in log ω in the vicinity of the peak. These data are plotted in Fig. 3a. As noted in the Introduction, we find that ω_κ² scales with ([small script l]_b[small script l]_f)²(ρ[small script l]_f)^β. Estimating β in the same way as α above, we find the best collapse for β_trial = 6.1, with E_min = 0.01. The cost function remains within 50% of its minimum for 5.8 ≤ β_trial ≤ 6.5.

The stretching-dominated band is characterized by a low frequency shoulder at a frequency ω_μ. While ω_μ is independent of κ, it shifts to lower frequencies when the fiber line density is increased, as shown in Fig. 5a. We discuss ω_μ further in the following section. The large-ω tail of the stretching-dominated band is due to longitudinal vibrations ω_‖ ∼ (μ/[small script l])^1/2 of individual fiber segments with length [small script l]. These frequencies are spread out by the distribution of segment lengths P([small script l]), which is exponential due to the Poissonian deposition process. The tail is then D ∼ P[[small script l](ω_‖)]|d[small script l]/dω_‖| ∼ ω⁻³ (Fig. 4, dashed line).

Fig. 5 (a) Density of states for κ = 10⁻¹⁴ and varying fiber line density. (b) The same data, now with ω rescaled by the prediction of eqn (18) for ω_μ.

3 Two types of soft modes

We will demonstrate that α = β + 1 in two separate ways, neither of which invoke the AFA. In the following section, we will directly calculate G in terms of the peak frequency ω_κ in the bending-dominated regime. But first, we present a scaling analysis of the characteristic frequencies in the density of states. The main question will be determining when the bending- and stretching-dominated bands merge, i.e. when ω_κ becomes comparable to ω_μ. Therefore we must first determine how ω_μ scales with fiber line density.

The stretching-dominated band is independent of the bending modulus κ, and therefore persists in freely-bending networks with κ = 0. The vibrational modes of central force networks have been studied extensively,^{23–25,27,28,31–33} and we will now adapt these results to Mikado networks to determine the scaling of ω_μ.

The mean coordination Z of a Mikado network obeys an upper bound Z ≤ Z_c = 4, where Z_c is known as the central force isostatic point.§ It represents the lowest coordination where there are enough springs to constrain all possible non-rigid body motions of the network. It is straightforward to show that ΔZ ≡ Z_c − Z ∼ 1/(ρ[small script l]_f) in Mikado networks. Therefore densely crosslinked, freely bending Mikado networks are central force networks close to, but below, their rigidity transition.^2,8,29

Central force networks on either side of Z_c display an abundance of soft modes whose spatial character strongly resembles floppy motions,^24,33i.e. the relative motion of each pair of nodes i and j connected by a spring is predominantly transverse to the unit vector n̂_ij pointing from i to j. We shall refer to these as stretching-dominated soft modes, to distinguish them from the soft modes of FMT. The softest stretching-dominated soft modes have frequency ω_μ. When Z < Z_c, the density of states below ω_μ is gapped, i.e. there are no modes between ω_μ and the floppy modes at ω = 0.²⁴ The dependence of ω_μ on ΔZ can be calculated with variational arguments reminiscent of FMT,³³ or with effective medium theory.²⁴ The result is

ω_μ ∼ k_eff^1/2|ΔZ|, (16)

where k_eff is an effective spring constantHere we have used the fact that the mean segment length scales as 1/ρ. The typical stretching-dominated soft mode is thereforeThis prediction is verified in Fig. 5b, where the low frequency shoulder in the density of states is collapsed for varying ρ[small script l]_f by rescaling the frequency ω with the prediction of eqn (18).

Having determined the scaling of ω_μ, we can now motivate the crossover between bending- and stretching-dominated shear response. We expect network response to project strongly on the softest modes. Therefore the response to shear will be dominated by fiber bending when the density of states has a well-defined bending-dominated band. This condition is satisfied so long as ω_κ ≪ ω_μ. The crossover from bending- to stretching-dominated response will occur when ω_κ and ω_μ become comparable,

This defines the bending-to-stretching crossover in terms of the exponent β. In the Introduction, we showed that the crossover occurs when ([small script l]_b/[small script l]_f)²(ρ[small script l]_f)^α = z₀. As these two conditions must be equivalent, we conclude that α = β + 1. This is our first main result.

4 Directly relating the shear modulus to soft modes

We now derive the relation α = β + 1 by a different but compatible route. Namely, we calculate G in the bending-dominated regime directly from the sum rule (14).

We begin by rewriting the sum rule as a bound on G, by restricting the sum to the O(N) soft modes,

We expect G to follow this bound closely, because the soft modes are the lowest non-trivial eigenfrequencies. Assuming each soft mode has the same frequency ω_κ² with a corresponding coupling strength Λ_κ² gives

In order to evaluate G, we must determine Λ_κ². Acting on eqn (11) from the left with 〈0,1|, solving for Λ_n², and taking the thermodynamic limit gives

Assuming that correlations between the force imbalances |Ξ〉 and the displacements |U_n〉 of a soft mode are negligible, the squared inner product isThe modes are normalized, so their mean squared displacement |U⃑_n,i|² is O(1/N_×).¶

The quantity requires closer consideration. Recall that |Ξ〉 represents the force imbalance on each crosslink after an affine deformation of the network.|| Affine displacements do not bend fibers, so the forces they generate are only due to stretching. Let us consider a single fiber. Affine shear changes the length [small script l]_α of each segment by an amount δ[small script l]_〈ij〉 = A(θ)γ[small script l]_α, where A(θ) is a prefactor that depends only on the orientation of the fiber and is therefore the same for all segments on the same fiber. The force in segment α is f_〈ij〉 = μδ[small script l]_〈ij〉/[small script l]_〈ij〉 = μA(θ)γ. Note that the force is independent of the segment length [small script l]_〈ij〉. As a result, affine shear induces an identical force in each segment of a fiber. The segments are co-linear prior to shearing, and therefore the net force on a crosslink vanishes unless the crosslink is at one of the two ends of a fiber (or is attached to a dangling end, if they have not been trimmed). In other words, is zero on interior crosslinks and O(μ²) on terminal crosslinks. As there are O(N) terminal crosslinks, rather than O(N_×), the sum in eqn (23) can be rewritten

Replacing the subscript n in eqn (22) with κ to indicate a typical soft mode, we findThe fact that Λ_κ² scales inversely with N² guarantees that the shear modulus is intensive. Crucially, the coupling strength also depends on ρ[small script l]_f. Inserting (25) in (21), we findfrom which it follows that α = β + 1. This is our second main result.

5 Discussion

We have demonstrated how soft modes determine the shear modulus of densely crosslinked Mikado networks in the bending-dominated regime. In so doing, we have resolved an apparent discrepancy in the literature between the predictions of floppy mode theory and independent measurements of the shear modulus. We find no reason to question the predictions of FMT regarding the scaling of the characteristic frequency ω_κ in eqn (2). However, eqn (3) and the prediction that α = β, are not correct. As these relationships rely on the affine fiber approximation, the natural conclusion is that the AFA does not hold. Our results imply that the typical fiber displacement during shear is δx ∼ (ρ[small script l]_f)^(α−β)/2[small script l]_fγ = (ρ[small script l]_f)^1/2[small script l]_fγ, which is much larger than the value [small script l]_fγ predicted by the AFA. We have been unable to find an intuitive explanation for this scaling relation.

It is natural to ask why the early work of Wilhelm and Frey yielded an estimate α ≈ 5.7 that is significantly smaller than subsequent estimates. In our opinion, this is because they fitted data for ρ > 15 to the form G/G_aff ∼ Δρ^α, where Δρ = ρ − ρ_c (recall ρ_c ≈ 5.7). While the difference between using ρ and Δρ in the scaling relation becomes negligible sufficiently far above the percolation threshold, fitting to a power law in Δρ gives systematically lower estimates of α. We have shown that FMT predicts that G scales as a power of ρ in the bending-dominated regime.

Our work has generated two additional insights into Mikado network mechanics. First, bending-dominated response correlates with the existence of a well-separated band of bending-dominated soft modes. We have shown that Mikado networks also possess stretching-dominated soft modes. The bending-dominated regime ends, and the gap between bending- and stretching-dominated modes closes, when the characteristic frequency of stretching-dominated soft modes, ω_μ, becomes comparable to the bending-dominated soft mode frequency ω_κ. The frequency ω_μ is set by the network's proximity to the central-force isostatic state. Our results therefore indicate that the isostatic point plays a governing role in the bending- to stretching-dominated crossover in 2D Mikado networks. This is distinct from fiber networks in 3D, which can undergo a similar crossover while remaining well below the central force isostatic point.³⁴

A second insight is that bending-dominated soft modes couple to shear in a non-trivial way. This coupling is encoded in the system-size dependence of Λ_κ². Simple dimensional analysis reveals that Λ_κ² must be inversely proportional to the square of an extensive quantity. We have shown that this squared extensive quantity is an admixture of the number of fibers and the number of crosslinks, Λ_κ² ∼ 1/(NN_×). This is a consequence of the co-linearity of crosslinks along a fiber. This prediction is corroborated by numerical studies in which the crosslink positions in the initial condition are gradually perturbed, which qualitatively changes the scaling of G with fiber line density.⁶

Our results suggest several directions for future work. In addition to elastic moduli, the density of states can be used to calculate linear response properties such as the storage and loss moduli,²⁰ sound transmission,³⁵ and thermal expansion coefficients.²⁶ It can also predict the onset of nonlinear properties such as strain stiffening and shocks.^27,36,37 Finally, a detailed understanding of structure–property relationships can facilitate the design of network architectures with tunable properties, such as moduli and failure onset.^13,38–42

Conflicts of interest

There are no conflicts of interest to declare.

Acknowledgements

It is a pleasure to acknowledge helpful discussions with Claus Heussinger. This work was supported financially by the Nederlandse Organisatie voor Wetenschappelijk Onderzoek (Netherlands Organization for Scientific Research, NWO), as well as by NWO Exacte Wetenschappen (Physical Sciences) through the use of supercomputer facilities.

References

1. J. Wilhelm and E. Frey, Phys. Rev. Lett., 2003, 91, 108103.
2. M. Latva-Kokko, J. Mäkinen and J. Timonen, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2001, 63, 046113.
3. D. A. Head, A. J. Levine and F. MacKintosh, Phys. Rev. Lett., 2003, 91, 108102.
4. D. Head, A. Levine and F. MacKintosh, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2003, 68, 061907.
5. C. Heussinger and E. Frey, Phys. Rev. Lett., 2006, 97, 105501.
6. C. Heussinger, B. Schaefer and E. Frey, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2007, 76, 031906.
7. A. Shahsavari and R. Picu, Int. J. Solids Struct., 2013, 50, 3332–3338.
8. M. F. Vermeulen, A. Bose, C. Storm and W. G. Ellenbroek, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2017, 96, 053003.
9. D. Zhou, L. Zhang and X. Mao, Phys. Rev. Lett., 2018, 120, 068003.
10. M. L. Gardel, F. Nakamura, J. H. Hartwig, J. C. Crocker, T. P. Stossel and D. A. Weitz, Proc. Natl. Acad. Sci. U. S. A., 2006, 103, 1762–1767.
11. G. Koenderink, M. Atakhorrami, F. MacKintosh and C. F. Schmidt, Phys. Rev. Lett., 2006, 96, 138307.
12. R. Tharmann, M. M. Claessens and A. R. Bausch, Phys. Rev. Lett., 2007, 98, 088103.
13. P. H. Kouwer, M. Koepf, V. A. Le Sage, M. Jaspers, A. M. Van Buul, Z. H. Eksteen-Akeroyd, T. Woltinge, E. Schwartz, H. J. Kitto and R. Hoogenboom, et al. , Nature, 2013, 493, 651–655.
14. C. Heussinger and E. Frey, Phys. Rev. Lett., 2006, 96, 017802.
15. D. Bi, J. Lopez, J. M. Schwarz and M. L. Manning, Nat. Phys., 2015, 11, 1074–1079.
16. D. M. Sussman and M. Merkel, Soft Matter, 2018, 14, 3397–3403.
17. D. J. Durian, Phys. Rev. Lett., 1995, 75, 4780–4783.
18. C. S. O'Hern, L. E. Silbert, A. J. Liu and S. R. Nagel, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2003, 68, 011306.
19. C. Heussinger and J.-L. Barrat, Phys. Rev. Lett., 2009, 102, 218303.
20. B. P. Tighe, Phys. Rev. Lett., 2011, 107, 158303.
21. K. Baumgarten and B. P. Tighe, Soft Matter, 2017, 13, 8368–8378.
22. D. J. Koeze, L. Hong, A. Kumar and B. P. Tighe, Phys. Rev. Res., 2020, 2, 032047.
23. B. P. Tighe, Phys. Rev. Lett., 2012, 109, 168303.
24. G. Düring, E. Lerner and M. Wyart, Soft Matter, 2013, 9, 146–154.
25. M. Yucht, M. Sheinman and C. Broedersz, Soft Matter, 2013, 9, 7000–7006.
26. C. Buss, C. Heussinger and O. Hallatschek, Soft Matter, 2016, 12, 7682–7687.
27. K. Baumgarten and B. P. Tighe, Phys. Rev. Lett., 2018, 120, 148004.
28. M. Merkel, K. Baumgarten, B. P. Tighe and M. L. Manning, Proc. Natl. Acad. Sci. U. S. A., 2019, 116, 6560–6568.
29. E. Huisman and T. C. Lubensky, Phys. Rev. Lett., 2011, 106, 088301.
30. S. Dagois-Bohy, B. P. Tighe, J. Simon, S. Henkes and M. van Hecke, Phys. Rev. Lett., 2012, 109, 095703.
31. M. Wyart, L. E. Silbert, S. R. Nagel and T. A. Witten, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2005, 72, 051306.
32. W. G. Ellenbroek, Z. Zeravcic, W. van Saarloos and M. van Hecke, EPL, 2009, 87, 34004.
33. L. Yan, E. DeGiuli and M. Wyart, Europhys. Lett., 2016, 114, 26003.
34. C. Broedersz, M. Sheinman and F. MacKintosh, Phys. Rev. Lett., 2012, 108, 078102.
35. K. Saitoh and H. Mizuno, 2020, arXiv: 2008.09760.
36. S. Ulrich, N. Upadhyaya, B. van Opheusden and V. Vitelli, Proc. Natl. Acad. Sci. U. S. A., 2013, 110, 20929–20934.
37. R. Rens, C. Villarroel, G. Düring and E. Lerner, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2018, 98, 062411.
38. M. M. Driscoll, B. G.-g. Chen, T. H. Beuman, S. Ulrich, S. R. Nagel and V. Vitelli, Proc. Natl. Acad. Sci. U. S. A., 2016, 113, 10813–10817.
39. L. Zhang, D. Z. Rocklin, L. M. Sander and X. Mao, Phys. Rev. Mater., 2017, 1, 052602.
40. B. E. Vos, L. C. Liebrand, M. Vahabi, A. Biebricher, G. J. Wuite, E. J. Peterman, N. A. Kurniawan, F. C. MacKintosh and G. H. Koenderink, Soft Matter, 2017, 13, 8886–8893.
41. D. R. Reid, N. Pashine, J. M. Wozniak, H. M. Jaeger, A. J. Liu, S. R. Nagel and J. J. de Pablo, Proc. Natl. Acad. Sci. U. S. A., 2018, 201717442.
42. R. Rens and E. Lerner, Eur. Phys. J. E: Soft Matter Biol. Phys., 2019, 42, 1–11.

---

Footnotes

† The scaling function (1) follows the discussion by Shahsavari and Picu.⁷ Other authors rescale their data differently, but still determine exponents that can be compared to α directly. Namely, α is equal to μ − 1 from ref. 1, and 2 + 2/z from ref. 3,4.

‡ Our simulations access a narrower range in ρ[small script l]_f because, unlike prior work, we determine each network's eigenfrequencies, which is significantly more computationally expensive than calculating their shear moduli.

§ The value Z_c = 4 is the prediction of a mean field counting argument due to Maxwell. While not all random networks obey mean field counting, Mikado networks do.

¶ Normalization guarantees |U⃑_n,i|² ∼ O(1/N_×), provided the mode is extended. As the floppy modes in freely bending Mikado networks are localized (motion is restricted to a central fiber and the secondary fibers with which it is crosslinked), one might worry that the soft modes are localized, as well. However, this is not the case. When the bending modulus is small but finite, the secondary fibers trigger motion in ternary fibers, and so forth. HF demonstrated that the axial displacements of the non-central fibers are of the same order of magnitude as the central fiber.

|| It may seem odd that the force imbalance due to an affine shear plays a role in determining G in the bending regime, where deformations are strongly non-affine. This can be understood by thinking of shear as a two-step process. In the first step, the system is deformed affinely, which, due to disorder, generates force imbalances on nodes. In the second step, the nodes displace non-affinely to relieve these force imbalances.

---