Effects of spherical confinement and backbone stiffness on flexible polymer jamming

TL;DR

This study uses molecular simulations to explore how spherical confinement and backbone stiffness influence the jamming behavior of flexible polymers, revealing unique mechanical properties and boundary effects compared to repulsive spheres.

Contribution

It introduces a detailed simulation analysis of polymer jamming in confined spaces, highlighting the effects of backbone stiffness and curved boundaries on mechanical properties.

Findings

Jamming occurs isostatically similar to repulsive spheres.

Curved walls and backbone bonds create new vibrational modes.

Boundary effects reduce the bulk modulus in polymers.

Abstract

We use molecular simulations to study jamming of a crumpled bead-spring model polymer in a finite container and compare to jamming of repulsive spheres. After proper constraint counting, the onset of rigidity is seen to occur isostatically as in the case of repulsive spheres. Despite this commonality, the presence of the curved container wall and polymer backbone bonds introduce new mechanical properties. Notably, these include additional bands in the vibrational density of states that reflect the material structure as well as oscillations in local contact number and density near the wall but with lower amplitude for polymers. Polymers have fewer boundary contacts, and this low-density surface layer strongly reduces the global bulk modulus. We further show that bulk-modulus dependence on backbone stiffness can be described by a model of stiffnesses in series and discuss potential experimental and biological applications.

Figures

Tables

Table 1. Table 1: Curve-fitted parameters for Eq. ( 9 ).

$\log p$	$B_{\infty }$	$ϵ$
-1	0.562 $\pm$ 0.005	0.039 $\pm$ 0.004
-2	0.2676 $\pm$ 0.0006	0.110 $\pm$ 0.002
-3	0.1854 $\pm$ 0.0004	0.194 $\pm$ 0.002
-4	0.1663 $\pm$ 0.0004	0.247 $\pm$ 0.003
-5	0.1612 $\pm$ 0.0008	0.260 $\pm$ 0.006

Equations

$$V_0(r_{ij})$$

$$V_B(r_{kl})$$

$$V_W(r_i)$$

$$N_0-N_S=dN-N_C.$$

$$N_S=N_C-N_C^{\mathrm{iso}}+1.$$

$$z=\frac{1}{N}\sum _{i=1}^N{z}_i,$$

$$\frac{N_S}{N}=\frac{\Delta z}{2}\equiv \frac{z-z^{\mathrm{iso}}}{2},$$

$$\frac{2N_C}{N}-z=\frac{N_W}{N}\sim \frac{1}{L},$$

$$\underset{p\to 0^{+}}{\lim }{N}_C^{\mathrm{rigid}}=d{N}^{\mathrm{rigid}}-f(d)+1=N_C^{\mathrm{rigid},\mathrm{iso}},$$

$$k_{\mathrm{eff}}={\left(\frac{N_0^{\mathrm{eff}}}{k_0}+\frac{N_B^{\mathrm{eff}}}{k_B}\right)}^{-1},$$

$$B=B_{\infty }(1+ϵ/{\epsilon }_B{)}^{-1},$$

$$\mathcal{E}(\{\textbf{r}_{i}\})=\mathcal{E}_{0}+u_{i}^{\mu}\frac{\partial\mathcal{E}}{\partial r_{i}^{\mu}}\bigg{|}_{\{\textbf{r}_{i}^{0}\}}+\\ \frac{u_{i}^{\mu}u_{j}^{\nu}}{2}\frac{\partial^{2}\mathcal{E}}{\partial r_{i}^{\mu}\partial r_{j}^{\nu}}\bigg{|}_{\{\textbf{r}_{i}^{0}\}}+\cdots,$$

$$\Delta V_0(u_{\perp })$$

$$\Delta V_B(u_{\perp })$$

$$\Delta V_W(u_{\perp })$$

$$C_{i^{\prime }\mu }^i{U}_i^{\mu }=E_{i^{\prime }}.$$

$$N_0=\mathrm{nullity}(C_{i^{\prime }\mu }^i).$$

$$Q_i^{i^{\prime }\mu }{T}_{i^{\prime }}=-F_i^{\mu }.$$

$$N_S=\mathrm{nullity}(Q_i^{i^{\prime }\mu })=\mathrm{nullity}(C_i^{i^{\prime }\mu }).$$

$$N_0-N_S=dN-N_C.$$

$$D_{i\nu }^{j\mu }=\frac{1}{m}{C}_i^{i^{\prime }\mu }{K}_{i^{\prime }}^{j^{\prime }}{C}_{j^{\prime }\nu }^j=\frac{1}{m}{Q}_i^{i^{\prime }\mu }{K}_{i^{\prime }}^{j^{\prime }}{Q}_{j^{\prime }\nu }^j,$$

Full text

I Introduction

The phenomenon of jamming is observed in systems ranging from granular materials flowing down a chute To et al. (2001) to biomacromolecules Jose et al. (2012); Keller et al. (2014); Berndsen et al. (2014); Keller et al. (2016). Most theoretical understanding of jamming comes from ideal models of granular materials, e.g., simulations of repulsive disks and spheres in periodic boundary conditions (PBCs). In general, a packing of a biopolymer into a container involves, at least, two unavoidable additional features: backbone bonds that link the material into a polymer chain and container walls that influence the material near the boundary.

As a step toward closing the large gap between the most idealized models and experimental systems, we investigate the effects on jamming of unbreakable adhesive bonds and external spherical confinement (SC).

Jamming occurs as the constituents of a flowing material sufficiently constrain one another’s motion, leading to a configuration that resists applied stress. A central question is whether the material jams isostatically, that is, in precise balance of constraints and degrees of freedom, consistent with the boundary conditions. Simulations of frictionless repulsive disks and spheres in PBCs have shown that the onset of rigidity occurs as a jump in particle-particle coordination number from zero to twice the dimensionality: four for disks and six for spheres, which correspond to isostaticity O’Hern et al. (2002); O’Hern et al. (2003); Liu and Nagel (2010); Goodrich et al. (2012); Kasahara and Nakanishi (2004); Song et al. (2008); Silbert (2010); Vinutha and Sastry (2016). We show that jamming in our simulations with backbone bonds and a concave confining wall occurs isostatically in fundamentally the same way as for repulsive spheres.

Repulsive spheres typically jam at a packing fraction of about 64%, which corresponds to the density of the maximally random jammed (MRJ) state Berryman (1983); Torquato et al. (2000); Scott (1960); Bernal and Mason (1960); Scott and Kilgour (1969); Berryman (1983); Clarke and Wiley (1987); O’Hern et al. (2002); O’Hern et al. (2003). Unlike repulsive spheres, a bead-spring model polymer has “built-in” constraints provided by backbone bonds. Isostatic packings of freely jointed chains of tangent hard spheres can be obtained at $\phi^{\mathrm{MRJ}}$ using algorithms that eliminate the effects of connectivity and allow effective equilibration through chain-connectivity-altering Monte Carlo moves Karayiannis and Laso (2008a, b); Foteinopoulou et al. (2008); Karayiannis et al. (2009a, b, c, 2010); Hoy and Karayiannis (2013); Karayiannis et al. (2013); Foteinopoulou et al. (2015); Ramos et al. (2018). However, when connectivity is preserved, approximately tangent fully flexible bead-spring chains jam at about 2% below $\phi^{\mathrm{MRJ}}$ in PBCs with little system-size dependence and retain a significant fraction of unconstrained degrees of freedom Hoy (2017). Confinement of monomers also reduces the jamming density by inducing layering near the boundary Verman and Banerjee (1946); Brown and Hawksley (1946); Liu et al. (1974); Macpherson et al. (1987); Zou and Yu (1995); Seidler et al. (2000); Landry et al. (2003); Mittal et al. (2008); Desmond and Weeks (2009); Camenen et al. (2013); Ghosh et al. (2016). We present both the reduction in density due to SC alone using repulsive sphere packings and the further reduction due to the polymer backbone that links all particles together.

On the other hand, few studies of polymer packings in confined geometries address mechanical properties. Previous investigations have largely focused on chain conformation within the packing Dickman and Hall (1998); Bitsanis and Hadziioannou (1998); Zou et al. (2009); Lopatina et al. (2011); Carrillo and Sumpter (2014) and topological ordering of segments Nikoubashman et al. (2017). Long polymers with specified bond-bond angles typically coil during packaging in SC to minimize bending energy Kindt et al. (2001); Arsuaga et al. (2002); LaMarque et al. (2004); Ali et al. (2004); Spakowitz and Wang (2005); Marenduzzo et al. (2009, 2010) and thus exhibit boundary-induced layering Spakowitz and Wang (2005). In contrast, we use a crumpled flexible-chain model to avoid coiling Ali et al. (2004) and to focus on the role of backbone connectivity in distinguishing the polymer from the monomer systems.

In Sec. III, we explain the necessity of using direct constraint counting rather than coordination number to assess the onset of rigidity due to unique considerations of systems in external confinement.

In Sec. IV.1, we show that essentially the same understanding of states of self-stress (SSSs) and zero modes in repulsive sphere packings can be extended to the case of a polymer in SC. In Sec. IV.2, we provide the distribution of jamming densities in simulations of spherically confined polymers and compare to those of monomers in SC and in PBCs to isolate the effects of backbone bonds and the confining wall. We find boundary-induced order in local density and coordination (Sec. IV.3) and in the vibrational density of states with effects on band structure due to the confining wall and the backbone (Sec. IV.4). Finally, we show how the bulk modulus changes due to these structural differences between monomers and polymers as well as due to the polymer backbone stiffness (Sec. IV.5).

II Simulation details

To study jamming of flexible polymers, we use three-dimensional molecular dynamics simulations Plimpton (1995) of single chains, each composed of $256\leq N\leq 8192$ monodisperse frictionless spherical particles of diameter $\sigma$. Each particle represents a monomer along a polymer chain, and interactions are governed by the following potentials:

[TABLE]

Nonconsecutive monomers interact via the harmonic repulsive potential $V_{0}(r_{ij})$, where $r_{ij}$ is the distance between the centers of particles (also referred to as sites) $i$ and $j$, $\varepsilon_{0}$ is the characteristic energy, and $\theta(x)$ is the Heaviside step function. Consecutive monomers $k$ and $l$ are bound by the two-sided harmonic potential $V_{B}(r_{kl})$ so that backbone bonds have energy scale $\varepsilon_{B}$ and rest length $\sigma$. To induce jamming, the polymer is confined by a spherical wall centered at the origin according to the radial harmonic potential $V_{W}(r_{i})$, where $r_{i}$ is the radial coordinate of site $i$ and $R$ is the wall radius. The total potential energy $\mathcal{E}$ is the sum of all pairwise and wall potentials. We consider at each site an equal point mass $m$, which sets the mass scale, and energies will be reported in units of $\varepsilon_{0}$, distances will reported in units of $\sigma$, pressures will be reported in units of $\varepsilon_{0}/\sigma^{3}$, and frequencies will be reported in units of $\sqrt{\varepsilon_{0}/m\sigma^{2}}$.

Disordered configurations are generated by thermalizing the polymer chains at temperature $kT=0.003$ in a large confining sphere at packing fraction $\phi=N\left(\frac{\sigma}{2R}\right)^{3}=0.02$. Each thermal configuration is then quenched to $T=0$ using the FIRE algorithm Bitzek et al. (2006). We compress each quenched system in small increments of $0.001\leq\Delta\phi\leq 0.01$ (adjusted by system size) by decreasing $R$ and minimizing energy after each compression until a jammed configuration is obtained, indicated by a nonzero $\mathcal{E}$. We then expand or compress these configurations to within 1% of each target pressure $p\equiv-\partial\mathcal{E}/\partial V$ where the system volume $V=4\pi R^{3}/3$ is that bounded by the confining sphere. For each system of size $N$, at least 100 random configurations are prepared, and each of these is studied at a large range of target pressures $10^{-7}\leq p\leq 10^{-1}$, bond energies $0.1\leq\varepsilon_{B}\leq 10$, and wall energies $0.1\leq\varepsilon_{W}\leq 10$.

The same procedures are repeated for nonbonded monomers in SC (where $\varepsilon_{B}=0$) and in PBCs [where $\varepsilon_{W}=\varepsilon_{B}=0$ and $\phi=\frac{\pi N}{6}\left(\frac{\sigma}{2R}\right)^{3}$ in a cubic domain with side length $2R$].

III Isostaticity and coordination

We review the analysis of the mechanical constraints that resist deformations and cause jamming. This allows us to introduce the effects of confining walls and adhesive bonds. Here, we introduce the index theorem, and in Appendix B, we derive the theorem in detail and explain associated subtleties.

When interested in the linear response at low pressures, near jamming, we may consider the unstressed network of a given system by replacing all contacts (including backbone bonds and wall contacts) with unstretched harmonic springs in an analysis following Ref. Lubensky et al. (2015). The mapping to the spring system is exact in the limit of zero pressure, and each spring introduces one harmonic constraint.

Each contact $i^{\prime}\leq N_{C}$, where $N_{C}$ is the number of contacts, is replaced by a harmonic bond of rest length $r_{i^{\prime}}$ equal to $r_{ij}$, $r_{kl}$, or $R-r_{i}$ [referring to Eqs. (1)] depending on the interaction. A zero mode is a normal mode of the system that causes no springs to be extended or compressed and corresponds to a motion with zero stiffness. A SSS is a set of extensions and compressions assigned to the springs that results in zero net force at each site. The index theorem embodies the fact that each contact either reduces the number of zero modes or increases the number of SSSs Lubensky et al. (2015), which, for a $d$-dimensional system with $dN$ degrees of freedom, is

[TABLE]

Creating a rigid (i.e., having no floppy modes), $d$-dimensional packing of spheres requires the number of constraints to match or exceed the degrees of freedom to be constrained Maxwell (1864). Therefore, $N_{C}\geq dN-f(d)$, where $f(d)$ is the number of zero modes associated with rigid-body motions. PBCs permit $f(d)=d$ rigid translations whereas a frictionless $(d-1)$-spherical boundary permits $f(d)=\frac{1}{2}d(d-1)$ rigid rotations. By its strictest definition Lubensky et al. (2015), an isostatic system contains neither floppy modes nor SSSs [$N_{0}=f(d)$, $N_{S}=0$]; however, jammed packings necessarily have at least one SSS ($N_{S}\geq 1$) corresponding to a nonzero modulus O’Hern et al. (2003); Goodrich et al. (2012) so that the number of contacts of a jammed isostatic system is $N_{C}^{\mathrm{iso}}=dN-f(d)+1$. Each additional constraint added to such a system creates an additional SSS,

[TABLE]

Constraints in repulsive sphere packings are commonly characterized by the average coordination number,

[TABLE]

where $z_{i}$ is the number of contacts of particle $i$, but this is less appropriate in confinement. First, without external confinement, as in PBCs, all contacts are between two particles, so $z=\frac{2N_{C}}{N}$ is twice the contact density, and the relation between $N_{S}$ and $z$ is

[TABLE]

with $z^{\mathrm{iso}}=\frac{2N_{C}^{\mathrm{iso}}}{N}=2d-\frac{2f(d)}{N}$. However, in external confinement, each wall contact involves only one particle. Since the wall itself is not counted as a particle, wall contacts do not get double-counted, and the coordination number $z$ is lower than twice the contact density by an amount that decreases with system size,

[TABLE]

where $N_{W}$ is the number of wall contacts and $L\equiv N^{1/d}$ is the linear system size. Therefore, $z$ is twice the density of constraints only when each contact constrains two degrees of freedom.

Second, previous studies of monomers have removed rattlers in order to isolate the rigid subsystem so that $\Delta z$ is directly related to $N_{S}$ O’Hern et al. (2003); Goodrich et al. (2012). Due to unbreakable bonds, polymers instead contain particles called flippers, which are constrained only by backbone bonds and thus can freely move tangent to their neighbors Karayiannis and Laso (2008b). To accurately analyze the rigid subsystem of a confined polymer, an analogous computation of $\Delta z$ would require both a boundary correction and the removal of all flippers and the backbone bonds constraining them.

IV Results

IV.1 States of self-stress and zero modes

In Fig. 1, we compute $N_{S}$ [Eq. (15)]. Our results show that $N_{S}\geq 1$ as seen in the splitting of $N_{S}/N$ curves to $1/N$ in the low-$p$ limit. Jamming in our systems, even with adhesive bonds and confinement, therefore corresponds to the introduction of a single SSS. For the polymers, the SSSs may contain both extended and compressed backbone bonds; indeed, we find that $\approx 30\%$ of backbone bonds are extended near the jamming transition so the ratio of extended to compressed backbone bonds is $\approx 0.5$.

We see the power-scaling law $N_{S}/N\sim\Delta z\sim p^{1/2}$, the same as for spheres O’Hern et al. (2003); Goodrich et al. (2012) for both monomers and polymers ($\varepsilon_{B}=1$) in confinement. This may be contrasted with a perfect $d$-dimensional crystal in external confinement, which would contain $N_{S}\gtrsim L^{d-1}$ at $p\to 0^{+}$. The increasing number of SSSs involves an increasing number of sites ($N^{\mathrm{rigid}}$) and engaged contacts ($N_{C}^{\mathrm{rigid}}$) as the rigid subsystem grows.

To quantify the number of unconstrained motions, we compute $N_{0}$ [Eq. (13)]. We find that $N_{0}^{\mathrm{mono}}>N_{0}^{\mathrm{poly}}$ in the low-$p$ limit. For monomers, these are primarily rattlers, which have no constraints, so each contributes $d=3$ zero modes. For polymers, these are primarily flippers; the smaller number of zero modes reflects the extra constraints from the backbone bonds that constrain motion even on particles outside the rigid subsystem. Because flippers can occur at chain ends and may involve consecutive polymer sites, directly computing the precise number of flippers from $N_{0}$ requires distinguishing topologically distinct groups and is not necessary to see that about 1 to 2% of the degrees of freedom are unconstrained even at moderate pressures. The significant number of unconstrained motions is consistent with other realistic packing protocols Ness et al. (2017); Plaza-Rivera et al. (2017).

The fraction of rattlers (flippers) decreases with system size. In the high-$p$ limit, no rattlers (flippers) remain as all particles become sufficiently coordinated that the only remaining zero modes are those associated with $f(d)$ rigid rotations within the spherical container. As pressure increases, particles rearrange to allow the system to relax. Rearrangements only result in small-scale configurational changes, even though the chain spans the full system.

Next, we delete rattlers and flippers, isolating the $N^{\mathrm{rigid}}$ particles and $N_{C}^{\mathrm{rigid}}$ engaged contacts of the rigid subsystem. At all pressures, we find that the number of zero modes that remain is again $f(d)$, indicating that no other zero modes are present in the rigid subsystem. Therefore, from Eq. (2),

[TABLE]

and we find that the rigid subsystem jams isostatically.

IV.2 Packing fraction at jamming

For reference, we provide the fraction of systems that are jammed $f_{J}$ at packing fraction $\phi$ as well as the average packing fraction at jamming $\phi_{J}^{N}$ for $256\leq N\leq 8192$ (Fig. 2).

Monomers in PBCs jam near 64% as expected for MRJ states for all system sizes. Confinement shifts jamming distributions to lower densities and increases system-size dependence. $\phi_{J}^{N,\mathrm{mono}}<\phi^{\mathrm{MRJ}}$, in agreement with previous studies of confined monomers Zou and Yu (1995); Seidler et al. (2000); Desmond and Weeks (2009); Camenen et al. (2013). The deviation of $\phi_{J}$ from $\phi^{\mathrm{MRJ}}$ is almost 6% at $N=256$ and diminishes to less than 1% by $N=8192$.

Figures 2(b) and 2(c) show that the inclusion of unbreakable backbone bonds further reduces the jamming density to almost 10% below $\phi^{\mathrm{MRJ}}$ at $N=256$ and 4% below at $N=8192$.

The $\approx 4\%$ difference between $\phi_{J}^{N,\mathrm{mono}}$ and $\phi_{J}^{N,\mathrm{poly}}$ persists across system sizes, similar to the density shift seen in jamming of flexible thermal polymers in PBCs Hoy (2017). Backbone bond stiffness has no appreciable effect on $\phi_{J}^{N,\mathrm{poly}}$ of flexible polymers, so only $\varepsilon_{B}=1$ data are shown in Fig. 2.

In addition, we consider monomer packings generated from jammed polymer configurations by deleting the backbone bonds. Without the extended bonds, the packings are unstable, and jamming is reattained at densities similar to the monomer distributions in Fig. 2(a).

IV.3 Boundary-induced structure

We compute the average local number density $\bar{n}$ and average local coordination $\bar{z}$ by binning point masses at $\{\textbf{r}_{i}\}$ and their respective coordination values $\{z_{i}\}$ over distance from the boundary $R-r$ (Fig. 3). Here, each nonbonded contact, backbone bond, and wall contact involving particle $i$ is included as one contact in $z_{i}$. Density layering is significant near the boundary (and, as reflected in Fig. 2, reduces $\phi^{N}_{J}$). The global number density $n$ and coordination number $z$ are shown as dashed lines in Fig. 3. Oscillations occur in both $\bar{n}$ and $\bar{z}$, similar to previous density profiles of confined monomers determined in experiments Seidler et al. (2000); Ghosh et al. (2016) and simulations Liu et al. (1974); Camenen et al. (2013); Macpherson et al. (1987); Landry et al. (2003); Mittal et al. (2008); Desmond and Weeks (2009); Spakowitz and Wang (2005) as well as tangent hard-sphere chains Foteinopoulou et al. (2015). Both oscillatory periods are consistent with the height of a regular tetrahedron (3-simplex) $\sqrt{2/3}\sigma\approx 0.82\sigma$ and agree with the well-established polytetrahedral structure of jammed monomer Clarke and Jónsson (1993); Anikeenko and Medvedev (2007) and polymer Karayiannis et al. (2009b, c, 2010); Hoy and Karayiannis (2013); Karayiannis et al. (2013); Foteinopoulou et al. (2015); Hoy (2017); Ramos et al. (2018) states. We note that sharply peaked maxima (minima) in $\bar{n}$ ($\bar{z}$) are separated by broad rounded minima (maxima). This qualitative “inversion” of curves would suggest that sites of high-density layers are, perhaps unintuitively, less coordinated than sites in the low-density layers between them. This could be rationalized by considering that particles in high-density layers sit between two lower-density layers with which they have fewer contacts than particles in low-density layers that sit between two high-density layers. However, the structure is even more complex than this as $\bar{z}$ curves are also shifted to the right of their inverted $\bar{n}$ counterparts; qualitatively, this phase shift appears to be about one-quarter of the period.

Although the curves are similar for monomers and polymers, a first noticeable difference is the height of the initial narrow peak at $R-r=\sigma/2$, indicating $N_{W}^{\mathrm{mono}}>N_{W}^{\mathrm{poly}}$, which becomes important to the bulk modulus as considered in Sec. IV.5.1. Additionally, for polymers, the oscillation amplitude of $\bar{n}$ is noticeably less than that of monomers, indicating that polymers exhibit less-extreme layering. In contrast, $\bar{z}^{\mathrm{mono}}<\bar{z}^{\mathrm{poly}}$ at nearly all points because of backbone bonds retained by flippers, which lead to the higher global coordination $z$ of polymers than monomer systems with fully uncoordinated rattlers.

IV.4 Density of states

To investigate the vibrational density of states, we construct the dynamical matrix $D_{i\nu}^{j\mu}$ [Eq. (17)]. The set of eigenvectors $\{U_{i}^{\mu}\}$ of $D_{i\nu}^{j\mu}$ are the polarization vectors of the system’s normal modes, and the eigenvalues $\{\lambda\}=\{\omega^{2}\}$ are the squared frequencies of the normal modes Ashcroft and Mermin (1976). From $\{\omega\}$, we compute the density of states $D(\omega)$. Since there is little variation among system sizes, we present only $N=2048$ data.

IV.4.1 Boundary modes

We compute $D(\omega)$ in systems with wall potentials $0.1\leq\varepsilon_{W}\leq 10$ for monomers in SC (Fig. 4). We first note that peaks at $\omega=0$ represent zero modes due to rattlers and rigid rotations. The curves have the universal characteristic shapes seen previously in disordered systems in PBCs O’Hern et al. (2003), the so-called boson peak at small finite $\omega$. However, wall potentials induce $N_{W}$ boundary modes with typical frequencies of $\omega_{W}\equiv\sqrt{\varepsilon_{W}/m\sigma^{2}}$, resulting in additional pronounced peaks.

At large $\varepsilon_{W}$, the additional modes lead to a band gap in $D(\omega)$. In this case, modes with $\omega>3$ may be isolated, and we bin the total set of polarization magnitudes $\{|\textbf{u}_{i}|\}$ over $R-r$ to compute the average polarization $\langle|\textbf{u}|\rangle$ with respect to distance from the wall (Fig. 4, the inset). The boundary modes are almost entirely localized to the two layers of sites nearest the boundary, giving the two distinct peaks in $\langle|\mathbf{u}|\rangle$.

IV.4.2 Backbone modes

Next, we see the effect of backbone-bond stiffness on the density of states [Fig. 5(a)]. Backbone interactions lead to a broad band approximately centered at $\sqrt{2}\omega_{B}\equiv\sqrt{2\varepsilon_{B}/m\sigma^{2}}$ as identified in Refs. Jain and de Pablo (2004); Milkus et al. (2018). The broadness of this band may be contrasted with the narrower and more structured boundary-mode band in $D(\omega)$. Like the high-$\varepsilon_{W}$ boundary band in Fig. 4, the high-$\varepsilon_{B}$ backbone band’s separation from the bulk band suggests a degree of independence in mode structure, and the density of states of the full system can be broken down into contributions from all three sources.

In Fig. 5(b), we replot $D(\omega)$ when $\varepsilon_{B}=10$ for the system in Fig. 5(a) but set $\varepsilon_{0}=\varepsilon_{W}=0$ in our computation of $K_{i^{\prime}}^{j^{\prime}}$ (see Appendix B). Bulk and boundary bands vanish into the $\delta$-function peak of zero modes, but we observe almost no change in the backbone band, highlighting its independence from the bulk band. A universal feature of the polymer vibrational spectra, the broad backbone band is a feature of real globular proteins ben Avraham (1993); Na et al. (2016). For better resolution of its features, we compute the average curve from $>20$ systems. Several pronounced peaks appear in the backbone band, which are similar to the signatures of analytically derived modes in collections of short chains of length $N_{\mathrm{ch}}\leq 5$ in PBCs Milkus et al. (2018). The most pronounced peak is at $\omega=\sqrt{2}\omega_{B}$, which corresponds to the vibrational frequency of the 1-simplex (a single bond) as well as a normal mode of the general 3-simplex. There are also small peaks at $\omega=\sqrt{2\pm\frac{1}{2}}\omega_{B}$ and $\omega=\sqrt{2\pm\sqrt{\frac{1}{2}}}\omega_{B}$, which correspond to vibrational frequencies of regular 2- and 3-simplices, respectively.

IV.5 Bulk modulus

IV.5.1 Effect of backbone connectivity

Plotting the bulk modulus $B\equiv\phi\partial p/\partial\phi$ of monomers and polymers over a range of $10^{-7}\leq p\leq 10^{-1}$ (Fig. 6), we find a constant, nonzero limit $\lim_{p\to 0^{+}}=B_{0}$, consistent with the power-law scaling relation $B\sim p^{0}$ O’Hern et al. (2003); Liu and Nagel (2010). As pressure increases from zero, $B$ remains within 1% of $B_{0}$ until $p\sim 10^{-4}$ whereas over this range $N_{S}$ increases by orders of magnitude from $N_{S}=1$ in the system sizes considered here (Fig. 1). $B$ also varies with $N$, mostly due to variation of $\phi_{J}^{N}$ with system size (Fig. 2).

The bulk modulus is substantially ($\approx 40$%) higher for monomers than for polymers. Variation in the prefactor $\phi$ in the definition of $B$ accounts for only a small part of the difference; $\phi_{J}^{N,\mathrm{poly}}$ is only $\approx 4$% lower than $\phi_{J}^{N,\mathrm{mono}}$ [Fig. 2(c)]. Therefore, it must also be that $(\partial p/\partial\phi)^{N,\mathrm{mono}}>(\partial p/\partial\phi)^{N,\mathrm{poly}}$. Section IV.1 showed that the rigid subsystems are nearly equal in size between the two system types ($N^{\mathrm{rigid,mono}}\approx N^{\mathrm{rigid,poly}}$, $N_{C}^{\mathrm{rigid,mono}}\approx N_{C}^{\mathrm{rigid,poly}}$), so the difference in $B$ must be due to configurational differences.

Recall that monomer packings have stronger layering and far more wall contacts than polymers (Sec. IV.3). Only wall contacts couple the motion of the wall to the interior packing, and therefore we may expect $B$ to rise with the wall contact density $N_{W}/A\sim N_{W}(\phi/N)^{2/3}$. We consider the modulus scaled correspondingly $B_{\mathrm{scaled}}\equiv\frac{B}{N_{W}}\left(\frac{N}{\phi}\right)^{2/3}$ and plot the ratio $B_{\mathrm{scaled}}^{\mathrm{mono}}/B_{\mathrm{scaled}}^{\mathrm{poly}}$ in Fig. 6(c). We see that this ratio is approximately $1$ for all system sizes and pressures, demonstrating that the difference in $B$ is primarily due to $N_{W}$.

IV.5.2 Effect of backbone stiffness

We also investigate the effect on $B$ of the backbone stiffness by varying $\varepsilon_{B}$ into both low-stiffness and high-stiffness regimes at pressures $10^{-5}\leq p\leq 10^{-1}$, plotted in Fig. 7. At low pressures, the bulk modulus vanishes if $\varepsilon_{B}\to 0$ as the configuration without backbone bonds is undercoordinated for rigidity. The bulk modulus saturates to a constant as $\varepsilon_{B}\to\infty$; backbone bonds become essentially inextensible compared to other contacts, yet the material can still deform around an infinitely stiff backbone. (In the equivalent case of decreasing $\varepsilon_{0}$, recall that the units of $\varepsilon_{B}$ and $B$ are proportional to $\varepsilon_{0}$ so that $B$ decreases proportionally to $\varepsilon_{0}$.)

To motivate a simple curve-fitting relation, consider that the material is isostatic at jamming, so the existence of the bulk modulus is dependent on every contact, similar to the simple situation of springs all in series. Given that $B$ is measured by isotropically deforming the wall, we therefore consider a different system: a one-dimensional chain of $N_{0}^{\mathrm{eff}}$ springs of stiffness $k_{0}\equiv\varepsilon_{0}/\sigma^{2}$ (these represent both wall and nonbonded-particle interactions since we have set $\varepsilon_{0}=\varepsilon_{W}=1$) and $N_{B}^{\mathrm{eff}}$ springs of stiffness $k_{B}\equiv\varepsilon_{B}/\sigma^{2}$ (representing backbone interactions). The chain’s overall effective spring constant is

[TABLE]

which is proportional to the bulk modulus $B=\beta k_{\mathrm{eff}}/\sigma$, where $\beta$ is a dimensionless constant. Rearranging Eq. (8) in terms of $\varepsilon_{B}$ and $B_{\infty}=\lim_{\varepsilon_{B}\to\infty}B$ yields

[TABLE]

with $B_{\infty}=\beta k_{0}/\sigma N_{0}^{\mathrm{eff}}$ and $\epsilon=\sigma^{2}k_{0}N_{B}^{\mathrm{eff}}/N_{0}^{\mathrm{eff}}$. In natural units $\sigma=k_{0}=1$, the energy scale $\epsilon$ represents the ratio $N_{B}^{\mathrm{eff}}/N_{0}^{\mathrm{eff}}$.

We plot curve fits using Eq. (9) in Fig. 7, which agree well with data for $p\leq 10^{-2}$; curve-fitted values of $\epsilon$ and $B_{\infty}$ are given in Table 1. The upward deviation in our data at $p=10^{-1}$ for the lowest $\varepsilon_{B}$ is a result of extreme compression and overcoordination as second-nearest-neighbor interactions occur, which the fitting form is not meant to capture. Pressure effects diminish in the low-$p$ limit.

V Discussion and conclusions

We have analyzed jammed configurations of a flexible bead-spring polymer in SC. Despite the presence of adhesive backbone bonds and spherical confining walls, the conditions at jamming superficially carry over from the case of repulsive spheres in PBCs. After accounting for the rigid-body motions within the spherical container, wall contacts, and underconstrained particles (rattlers and flippers), we see that jamming occurs exactly at isostaticity and coincides with the emergence of a single SSS. Jamming occurs at somewhat reduced density compared to monomers, and, upon further compression, the number of SSSs scales as the square root of pressure as for monomers in SC.

The boundary causes layering in both local density and coordination, which are, unexpectedly, out of phase; qualitatively, curves for density and coordination are inverted in shape and phase-shifted $\approx\pi/2$. The boundary also introduces a narrow band of vibrational modes into the density of states with characteristic frequency scaling with the square root of the wall stiffness. At high wall stiffness, these modes are highly localized to the outermost two layers of sites. The independence of boundary modes from bulk modes extends to backbone modes; bands generated by high-stiffness backbone bonds are virtually unchanged after the removal of nonbonded and boundary potentials. Not only do these bands follow a universal pattern, but they also display peaks corresponding to regular low-dimensional simplices, indicating the possibility of inferring aspects of the internal structure from the vibrational spectrum.

The higher number of wall contacts in monomer packings raises the bulk modulus by $\approx 40\%$ compared to polymers. An explanation comes from a model of stiffnesses in series that scales with the wall contact number. A similar conceptual model motivates a fitting relation that describes the dependence of the bulk modulus on backbone stiffness and predicts its value in the limiting case of incompressible backbone bonds.

Although packing of a flexible-chain polymer is a highly idealized model of a biopolymer, several insights may apply immediately to experiment. The vibrational states convey information about the strength of confinement, the number of boundary constraints, and the backbone configuration, which could be exploited to study and potentially manipulate polymer structure. Our results may also apply to the cytoskeleton, the protein network that spans the cell from the nucleus to the cell membrane and accounts for cytoplasmic structure and rigidity. The number of contact points with the cell membrane may be strongly linked with cellular compressibility and membrane flexibility. This dependence could be measured experimentally, e.g., via atomic force microscopy. Further biological relevance could be found within the cell nucleus where our model may help elucidate the envelope’s influence on chromatin structure and mechanics. In addition, we hope this paper clarifies fundamental aspects of jamming with regard to internal constraints and the finite boundaries present in all real systems.

Future analysis may investigate the spatial structure of SSSs in SC, the origin of the apparent phase shift in local density and coordination, or the material elasticity at higher pressure and with higher-curvature walls where internal stresses and higher-order terms in the energy expansion are relevant (discussed in Appendix A). A fuller analysis may also consider nonbonded adhesive interactions, backbone-bending stiffness, dihedral stiffness, bond stresses, or finite temperature to yield more accurate models of real biopolymers.

Acknowledgements.

We are grateful for helpful discussions with A. Liu, C. Goodrich, and T. Lubensky. This research was funded, in part, by the U.S. DOE, Office of Basic Energy Sciences, Division of Materials Sciences and Engineering (Award No. DE-FG02-05ER46199). Computational support was provided by the LRSM HPC cluster at the University of Pennsylvania.

Appendix A Unstressed-network approximation

We first explain the approximation of low stresses, accurate at low pressures, that permits the index theorem analysis based on an unstressed spring network (Appendix B). In the following, we use Einstein notation, and sites are labeled by plain Roman indices, bonds are labeled by primed Roman indices, and Cartesian components are labeled by Greek indices.

We explicitly calculate the lowest-order energy terms to analyze stability of a static configuration $\{\textbf{r}_{i}^{0}\}$ after energy minimization, i.e., in force balance, with potential energy $\mathcal{E}_{0}=\mathcal{E}(\{\textbf{r}_{i}^{0}\})$. Let $\textbf{r}_{i}$ be the position of particle $i$ and $\textbf{u}_{i}=\textbf{r}_{i}-\textbf{r}_{i}^{0}$ be its displacement from its reference position. For small displacements, we can Taylor expand the energy,

[TABLE]

where terms proportional to $u_{i}^{\mu}$ are zero since we expand about a stable configuration.

All potentials in the simulation where nonzero have the form $\mathcal{V}(r)=\varepsilon(1-r/d)^{2}/2$, where $r=|\mathbf{r}|$ corresponds to displacements $|\mathbf{r}_{i}-\mathbf{r}_{j}|$, $|\mathbf{r}_{k}-\mathbf{r}_{l}|$, or $|\mathbf{r}_{i}|$ [referring to Eqs. (1)]. A displacement component parallel to the interaction direction $u_{\parallel}\equiv\textbf{u}\cdot\textbf{r}/r$ corresponds to stiffness $\kappa\equiv\partial^{2}\mathcal{V}/\partial u_{\parallel}^{2}=\varepsilon/d^{2}$. A component perpendicular $u_{\perp}\equiv|\textbf{u}-u_{\parallel}\textbf{r}/r|$ also has finite stiffness $\partial^{2}\mathcal{V}/\partial u_{\perp}^{2}=\kappa(1-d/r)$, which is positive for wall contacts and extended backbone bonds. Explicitly, the change in energy due to a small displacement perpendicular to $\mathbf{r}_{ij}$, $\mathbf{r}_{kl}$, or $\mathbf{r}_{i}$ is

[TABLE]

where $\theta(x)$ is omitted for brevity. The prefactor of the quadratic term is negative for overlapping monomers ($r_{ij}<\sigma$), so energy decreases in the perpendicular direction, and the particles tend to slip off one another. Small displacements perpendicular to extended backbone bonds ($r_{kl}>\sigma$) or tangential to the wall instead require an increase in energy, resulting in linear restoring forces. Prefactors vanish in the unstressed, i.e., zero-energy limit ($r=d$), so the energy costs of these motions appear only at $O(u_{\perp}^{4})$ and produce no linear response. Figure 8 illustrates the tangential curvature of $\mathcal{V}(r)$ for $r<d$ ($\partial^{2}\mathcal{V}/\partial u^{2}_{\perp}<0$), $r=d$ ($\partial^{2}\mathcal{V}/\partial u^{2}_{\perp}=0$), and $r>d$ ($\partial^{2}\mathcal{V}/\partial u^{2}_{\perp}>0$).

Although extended backbone bonds and wall contacts constrain tangential motion, we focus on the unstressed case (at jamming) where the mapping to unstretched springs is exact for harmonic analysis Lubensky et al. (2015). In that case, each interaction constrains motion only along the interaction direction. Therefore, at sufficiently low pressure, only relative motion in the direction normal to the contact contributes significantly to the linear response, but higher-order terms in the energy expansion can, in principle, affect jamming in some materials, e.g., they are seen to stabilize zero-frequency modes in packings of aspherical particles Donev et al. (2007); Mailman et al. (2009); VanderWerf et al. (2018). At rest length, nonzero contributions up to fourth order in the expansion come from the terms $\kappa u_{\parallel}u_{\perp}^{2}/6d$, $-\kappa u_{\parallel}^{2}u_{\perp}^{2}/12d^{2}$, and $\kappa u_{\perp}^{4}/8d^{2}$, so the curvature of confining walls and the adhesive regime of backbone bonds may contribute to higher-order stability at zero pressure. In principle, these terms may be able to stabilize zero modes in the packing; however, as stated in the main text, after rattlers and flippers have been deleted, no zero-frequency modes are present in our packings other than rigid-body motions, indicating that harmonic analysis accounts for all constraints in our sphere packings.

Appendix B Derivation of the index theorem

Site displacements form the $dN$-dimensional displacement vector $U_{i}^{\mu}$, where $\mu$ indexes the $d=3$ Cartesian components of each vector $\textbf{u}_{i}$. The linear operator $C_{i^{\prime}\mu}^{i}$, termed the compatibility matrix, maps $U_{i}^{\mu}$ to the $N_{C}$-dimensional bond elongation vector $E_{i^{\prime}}\equiv\frac{\partial r_{i^{\prime}}}{\partial r_{i}^{\mu}}U_{i}^{\mu}$,

[TABLE]

Since a zero mode is described by a set of displacements that causes no bond elongations, the null space of $C_{i^{\prime}\mu}^{i}$ is spanned by modes associated with both floppy modes and global rigid-body motions of which there are in total,

[TABLE]

Conversely, we may consider the resulting force on each site as the linear response to a tension vector $F_{i}^{\mu}\equiv-\frac{\partial r^{i^{\prime}}}{\partial r_{\mu}^{i}}T_{i^{\prime}}$. We then obtain the equilibrium matrix,

[TABLE]

Comparing with Eq. (12), we see that, in matrix form, $Q_{i}^{i^{\prime}\mu}$ is the transpose of $C_{i^{\prime}\mu}^{i}$.

In certain networks, the bonds may be placed under tension or compression while maintaining zero net force on each site, i.e., $Q_{i}^{i^{\prime}\mu}T_{i^{\prime}}^{S}=0$. Such a tensional state $T_{i^{\prime}}^{S}$ is referred to as a SSS and is contained in the null space of $Q_{i}^{i^{\prime}\mu}$. The number of SSSs in a system is thus given by

[TABLE]

From the rank-nullity theorem and given $\mathrm{rank}(C^{i}_{i^{\prime}\mu})=\mathrm{rank}(Q_{i}^{i^{\prime}\mu})$, we obtain the index theorem Lubensky et al. (2015),

[TABLE]

Finally, we note the connection to the dynamical matrix, defined as

[TABLE]

where $K_{i^{\prime}}^{j^{\prime}}\equiv\partial^{2}V(r_{i^{\prime}})/\partial r_{j^{\prime}}^{2}$ is the diagonal stiffness matrix.