Generic stress rectification in nonlinear elastic media

Stress propagation in nonlinear media is crucial in cell biology, where molecular motors exert anisotropic force dipoles on the fibrous cytoskeleton. While the force dipoles can be either contractile or expansile, a medium made of fibers which buckle under compression rectifies these stresses towards a biologically crucial contraction. A general understanding of this rectification phenomenon as a function of the medium's elasticity is however lacking. Here we use theoretical continuum elasticity to show that rectification is actually a very general effect in nonlinear materials subjected to anisotropic internal stresses. We analytically show that both bucklable and constitutively linear materials subjected to geometrical nonlinearities rectify small forces towards contraction, while granular-like materials rectify towards expansion. Using simulations, we moreover show that these results extend to larger forces. Beyond fiber networks, these results could shed light on the propagation of stresses in brittle or granular materials following a local plastic rearrangement.


I. CONSTITUTIVELY LINEAR MEDIA RECTIFY TOWARDS CONTRACTION
In this section, we consider only geometrical nonlinearities and set out to prove the inequality in Eq. (3). Firstly, the elastic energy density E is written with the Green-Lagrange strain tensor ε = (η + η T + η T η)/2, which depends quadratically on the displacement gradient η ij = ∂u i /∂x j . We further define the symmetric part of the displacement gradient U ij = (η ij + η ji )/2 which corresponds to the linearized strain. The stress tensor which naturally derives from E describes the surface force measured in the initial space with respect to the initial area: dF 0 /da = ∂E/∂ε, it is known as the second Piola-Kirchhoff stress. Then in order to find the Cauchy stress measured fully in the target space: σ = dF/dA, we need to transform the surface force and the area as dF = (1 + η) dF 0 and da = 1 + η T det (1 + η) dA, (S1) which ultimately gives the formula for the Cauchy stress in the main text [1]. Then, given the quadratic energy density E = κε 2 ii /2 + µ(ε 2 ij − 2 ii /d), the stress-strain relation displays a linear stress term σ L proportional to U, and a term which includes the geometrical nonlinearities σ G : Secondly, given the force balance equation f i = −∂σ ij /∂X j , the difference between the local and active coarse-grained stresses can be integrated by part to read where dv is the volume element in the initial space. Equation (S3) is known as the mean stress theorem [2,3]. Due to this relation, writing S = σ det(1 + η) and decomposing S = S L + S G as we did σ in Eq. (S2), the pressure difference reads P a − P l = − S ii /(V d). Due to our fixed boundary condition, the integral of the trace of the linear term S L = σ L vanishes and the trace of the nonlinear term is expressed in a closed form as Here, κ and µ are both positive for mechanical stability. Thus, for d ≥ 2, the geometrical term S G ii always gives a positive (contractile) contribution to the active stress. Indeed, using the eigenvalues λ i ∈ R of the symmetric matrix ε, we can rewrite ε 2 ij − ε 2 ii /d = i<j (λ i − λ j ) 2 /d, which is always non-negative. As a result S G ii is a sum of squares that is also non-negative, implying the inequality of Eq. (3): P a − P l 0.
Finally, we present an alternative derivation of this relation that highlights its frame indifference. The integrand of Eq. (3) can be rewritten using the deformation gradient Λ = 1 + η and the right Cauchy-Green deformation tensor we can write We see that the right-hand side Eq. (S6) is a sum of squares plus a term ∝ C ii = d + 4η ii + 4η 2 ij . Since η ii is integrated to zero due to the fixed boundary condition, the term in Eq. (S6) also gives a contractile contribution to P a − P l . Here, we derive typical values of κ 1 and µ 1 for granular media near the jamming transition. Let us consider a large equilibrium packing of bidisperse frictionless spherical grains in a 2D box with volume fraction φ. The grains interact through a harmonic potential V ∼ k δ 2 , where k is a spring constant, and δ the overlap divided by the sum of the two bead diameters, see Fig. S1. For volume fractions slightly above the jamming transition φ c ≈ 0.84, granular media display a strongly nonlinear elastic behavior. Indeed, multiple simulations and experiments [4,5] have shown that while the bulk modulus goes to a finite limit in φ + c and can thus be approximated by a constant, the shear modulus scales with φ − φ c and vanishes at the transition. This can be expressed as where K 0 , G 0 ≈ 0.2 and p ≈ 0.5. Based on this model, we impose an isotropic compression characterized by a displacement gradient tensor η ij = −η 0 δ ij /d on our granular material initially at φ c , where 0 < η 0 1. This results in a new volume fraction φ 0 ∼ φ c (1 + η 0 ). We then compute the elastic moduli κ and µ, and their nonlinear corrections κ 1 and µ 1 around this value of φ 0 . Let the bulk strain η ii = −η 0 + δη, where |δη| η 0 , corresponding to a volume fraction φ ∼ φ 0 − φ c δη. Then similarly to Eq. (6), the moduli are expressed as K ∼ κ(1 + κ 1 δη) and G ∼ µ(1 + µ 1 δη), where the parameters are derived from Eq. (S7): Therefore, while κ 1 vanishes, µ 1 diverges at the transition and scales as (φ 0 − φ c ) −1 . And close to the transition, around e.g. φ 0 − φ c = 0.001, 0.01 or 0.1, we find respectively µ 1 ≈ −400, −40 or −4 as in the main text.

III. COARSE-GRAINED STRESSES IN THE CIRCULAR GEOMETRY
In this section, we present the analytical calculations leading to the expressions of the coarse-grained stressesσ in Eq. (9) and of the rectification coefficient α of Eq. (10) of the main text. We first rewrite the coarse-grained stresses in the initial space where the calculations will be easier to handle in Sec. III A. Then in Sec. III B we present the Ansatz for the displacement field that allows us to solve the force balance condition. In Sec. III C, we show the expressions of the coarse-grained pressures and shear stresses with the stress and displacement fields. We finally display the detailed expressions of the coarse-grained stressesσ in the circular geometry as well as the expression of α in Sec. III D.
A. Rewriting the coarse-grained stresses in the initial space In order to calculate the coarse-grained stresses, we need to express them in the initial configuration, where the force density φ = f / det(1 + η) is related to the first Piola-Kirchhoff stress tensor τ = (1 + η) ∂E ∂ε . The force balance equation thus reads φ i = −∂ j τ ij and the coarse-grained stresses can be rewritten as where da is the outward-directed area element in the initial space. As before in Sec. I, we distinguish the linear term τ L (which is the same in the initial and target spaces τ L = σ L ), from the nonlinear term τ NL = τ G + τ C . This δ FIG. S1. Overlap between two interacting beads in a granular simulation.
3 last equality distinguishes the geometrical and constitutive nonlinearities. Given the non-harmonic energy density of Eq. (4), up to second order, the terms read (S10a) (S10b) (S10c)

B. Ansatz for the displacement field
In the initial space, the material is subjected to zero body force except at r in , where the stress is discontinuous. The fixed boundary at r out and the imposed displacement at r in additionally impose boundary conditions on the stress and displacement fields, resulting in the following system of equations: (S11) We solve this system perturbatively by expanding u, η, τ in the small scalar quantity The displacement gradient is hence written η = η L + η NL + O η 3 where the L and NL superscript refer to linear and quadratic (nonlinear) terms in η. This allows us to write the stress tensor as where τ L η L is of order 1, while the next two are of order 2. The linear displacement field u L is the solution of ∂ i τ L ji (η L ) = 0. This is solved by decomposing u L in the following Fourier modes due to the form of the imposed displacement: u L r /r in = e 0 ζ 0 (r) + e 2 ζ 2 (r) cos 2θ, u L θ /r in = e 2 ω 2 (r) sin 2θ.
(S14a) (S14b) Here ζ 0 (r), ζ 2 (r), ω 2 (r) are sums of r k , with k ∈ {−3, −1, 1, 3} and coefficients depending on the boundary conditions. Then, at the first nonlinear order, u NL is the solution of the linear equation ∂ i τ L ji (η NL ) = −∂ i τ NL ji (η L ) which is solved by expanding u NL as u NL r /r in = e 2 2 ξ 0 (r) + e 0 e 2 ξ 2 (r) cos 2θ + e 2 2 ξ 4 (r) cos 4θ, (S15b) where the ξ i (r), π i (r) are again sums of r k with k odd between −7 and +5. As a consequence, we obtain in Sec. III D the strain and stress fields up to second order in η.

C. Calculations of the coarse-grained stresses
We compute the coarse-grained active stressσ a and local stressσ l in the circular geometry by integrating the stresses in the material as in Eq. (S9). The coarse-grained stresses are expressed in Cartesian coordinates, while the stress field is more easily expressed in polar coordinates. In the following (and in this subsection only), we denote Cartesian indices x, y with Greek letters (µ, ν) and polar indices r, θ with Latin letters (i, j, k). The change-of-basis matrix between these two systems reads R = cos θ − sin θ sin θ cos θ . In the circular geometry of Fig. 1 where the active unit 4 produces a discontinuity in the stress at r in in the initial configuration, the coarse-grained stresses are expressed as follows: (S16b) where ∂ k τ ik denotes the stress divergence expressed in polar coordinates. Then, introducing β = (r out /r in ) 2 , the active pressure and shear stress rescaled so as to compensate for dilution read β(σ a xx +σ a yy ) = −2P a = β π τ rr (r out , θ), with similar expressions for the local pressure and shear stress: (S18a) (S18b)

D. Full expressions of the coarse-grained stresses
In Eqs. (9-10), we consider small displacements with two independent parameters e 0 , e 2 1. In the weakly nonlinear formalism, the simplest possible rectification requires that the term in e 0 be similar (and of opposite sign) to the term in e 2 2 . In this regime, given 1, the displacement parameters read e 0 = ẽ 0 and e 2 = √ ẽ 2 . As a result P x = P x and S x = √ S x for x ∈ {l, a}, where the tildes denote quantities of order one. Then to lowest order in , Eq. (9) can be rewritten asP (S19b) implying that the "∼" symbols of Eq. (10) denote equalities to lowest order in . The rectification behavior thus depends on the ratioP l /(αS 2 l ), i.e. onẽ 0 /ẽ 2 2 . We further display the complete expressions of the coefficients in Eq. (9) obtained after finding the strain field through force balance (see Sec. III B) and integrating the resulting stress via Eq. (S10), as in (S17) and (S18). In order to make sense of the cumbersome expressions of the A x , B x and C x , we introduce several quantities: In the end, we find As expected from the linear elasticity analysis of Sec. I, the coefficients in front of the linear e 0 and e 2 terms are identical for the local and boundary stresses, but discrepancies appear in the e 2 2 terms. This leads us to define α = (B a − B l )/C 2 l .

IV. BEHAVIOR OF α AND RECTIFICATION SATURATION RADIUS r *
To help better understand the lengthy expression of the rectification coefficient α = (B a − B l )/C 2 l in subsection III.D, we hereby discuss its dependence on the system size r out . As is apparent from Fig. 2(a), α increases with increasing r out for relatively small systems, then saturates as the size of the system goes to infinity. Indeed, away from the high-stress region close to the active unit, the stress decrease causes the nonlinearities to become negligible in front of the linear terms. Therefore, we examine the radius r * at which stress propagation switches from nonlinear to linear. To this end we define the system size parameter β = (r out /r in ) 2 . In the limit β → ∞, Eq. (S20) leads to , We introduce the value of the parameter β such that α is within 10% of its large-size limit through α(β10)−α∞ α∞ = 0.1. The square-root of β 10 , corresponding to the ratio of the radii, lies between 4 and 9, except for insignificant values  Fig. 2(a) showing the stabilization of the graphs as the system size increases for a Poisson's ratio ν = 1. The blue and yellow lines at constant |α|µ stay quite still between rout/rin = 8 and 10. The behavior is similar for ν = 0. of α, see Fig. S2. Fig. S3 also illustrates this behavior by showing the stabilization of the lines at constant α as the system size increase. Therefore, defining the rectification saturation radius r * such that β 10 = (r * /r in ) 2 , increasing r out past r * ∼ 10r in has little influence on the value of P a − P l , i.e. on the rectification effect. This indicates that the propagation is nonlinear only up to r * . In the study of stress propagation from multiple active units, one thus needs to compare this r * to the typical spacing between two active units.

V. THE RECTIFICATION DIAGRAM
Let us give further explanation to the shadings of the rectification diagram of Fig. 2(b), corresponding to different rectification regimes. In the circular geometry of Fig. 1, provided that α and P l have different signs, a change of sign of P a due to rectification can appear for all values of the local pressure P l as long as the local shear stress |S l | is large enough. Indeed, in Eq. (10) the sign switching of the active pressure P a (e.g. P a < 0 while P l > 0) requires This sets the boundary between the regions with light shading and the regions with intermediate shading. Then, the extreme case where all active stress components switch sign (e.g. P a ± S a < 0 while P l ± S l > 0) happens for i.e. for |P l | and |S l | both larger than 2/|α|. This extreme case corresponds to the dark regions of Fig. 2(b).

VI. FINITE-ELEMENT SIMULATIONS
This section provides further details on the finite element simulations used to produce Fig. 3 of the main text. Sec. VI A describes our simulation methods. In Sec. VI B, we discuss rectification in two additional fully nonlinear models.

A. Methods
We solve the set of equations (S11) via simulations with the finite element software Fenics [6] version 2019.2.0.dev0. We use a mesh with maximal size l = 0.01 for r in = 1 and r out = 2, and another one with l = 0.1 for r out = 10. They were both created with Gmsh version 4.4.1. In all figures, the error bars correspond to the differences between two meshes at l and l/10, which gives roughly 5% of S l or P l for all points. The meshes are created without enforcing rotational symmetry, which results in small non-zero values for the non-diagonal coefficients that should be zero in a continuum system (e.g.,σ l xy ), as shown in Eq. (8). However, we find that these values are smaller than 5% of the diagonal coefficients in all simulations. In the geometry of Fig. 1, if we apply too large a deformation at r in , we come into contact with the fixed boundary at r out , which poses some numerical issues. Therefore, we can only perform accurate simulations up to about η = |e 0 | + |e 2 | ∼ 0.6.

B. Additional data
In the main text, we studied two models with clear buckling and anti-buckling behaviors, which lead to a readily observable reversal of the active pressure sign, due to rectification. We also studied a standard neo-Hookean model of rubber in which this behavior is less pronounced. Here, we study another model which can mimic the shear-stiffening behavior of fiber networks. Consistent with analytical predictions, this system displays a smaller propensity for rectification, and we show that the predictions of Eq. (10) remain valid up to intermediate stress values.
In the fully nonlinear model with the elastic energy density of Eq. (11), we introduced the parameters a, b such that κ 1 = 1/2 − 3a and µ 1 = −3/2 − b. In Fig. 3, we showed the good agreement between the weakly nonlinear predictions of Eq. (10) and finite-element simulations for the bucklable and anti-bucklable models obtained by setting (a, b) = (−1/6, −5/2) and (3/2, 5/2) to obtain (κ 1 , µ 1 ) = (1, 1) and (−4, −4) respectively. As well as for the neo-Hookean model, obtained by setting a = b = 0 in Eq. (11), which has κ 1 = 1/2, µ 1 = −3/2. In Fig. S4(a) small stress regime for all considered values of the Poisson's ratio ν and the ratio of the boundary radius to the active unit radius r out /r in . We then investigate a variant of the neo-Hookean model which has the shear-stiffening behavior G ∝ σ 3/2 xy characteristic of fiber networks under large strains. Its elastic energy density reads [7,8] where J = det(1 + η) and I = Tr(1 + 2ε). Here, the shear strain threshold for the stiffening behavior corresponds to a fraction of 1/ √ c, the strain at which the shear stress diverges. This is such that the neo-Hookean model is recovered for c = 0. This model still has κ 1 = 1/2, µ 1 = −3/2, which corresponds to the same tendency to rectify towards contraction as in the neo-Hookean case. Indeed, c only affects higher order nonlinearities. As shown in Fig. S4(b,c) where c = 10, we recover Eq. (10) at small stress. But due to the shear stress divergence at finite shear strain in Eq. (S24), our predictions fail when S l /κ 1/ √ c. We finally compare the dependence of P a with P l at zero shear stress with the weakly nonlinear prediction P a ∼ P l . As displayed in Fig. S5, we recover the prediction for small stresses, but higher order nonlinearities induce significant deviations for P l /κ outside of [−0.3, 1]. In the end, we see that the predictions of Eq. (10) relating the active and local stress components fail when either the local pressure or the local shear stress become comparable to the bulk modulus κ. Our weakly nonlinear predictions additionally fail close to stress divergences, i.e. when the nonlinearities become too significant compared to the linear terms.