Persistently Auxetic Materials: Engineering the Poisson Ratio of 2D

Jul 18, 2016 - size membranes with unclamped boundary conditions have positive Poisson's ratio due to spontaneous non-zero mean curvature, which can ...
0 downloads 0 Views 2MB Size
Download PDF
Persistently Auxetic Materials: Engineering the Poisson Ratio of 2D Self-Avoiding Membranes under Conditions of Non-Zero Anisotropic Strain Zachary W. Ulissi,†,‡ Ananth Govind Rajan,† and Michael S. Strano*,† †

Department of Chemical Engineering, Massachusetts Institute of Technology, Cambridge, Massachusetts 02139, United States ABSTRACT: Entropic surfaces represented by fluctuating two-dimensional (2D) membranes are predicted to have desirable mechanical properties when unstressed, including a negative Poisson’s ratio (“auxetic” behavior). Herein, we present calculations of the strain-dependent Poisson ratio of self-avoiding 2D membranes demonstrating desirable auxetic properties over a range of mechanical strain. Finitesize membranes with unclamped boundary conditions have positive Poisson’s ratio due to spontaneous non-zero mean curvature, which can be suppressed with an explicit bending rigidity in agreement with prior findings. Applying longitudinal strain along a singular axis to this system suppresses this mean curvature and the entropic out-of-plane fluctuations, resulting in a molecular-scale mechanism for realizing a negative Poisson’s ratio above a critical strain, with values significantly more negative than the previously observed zero-strain limit for infinite sheets. We find that auxetic behavior persists over surprisingly high strains of more than 20% for the smallest surfaces, with desirable finite-size scaling producing surfaces with negative Poisson’s ratio over a wide range of strains. These results promise the design of surfaces and composite materials with tunable Poisson’s ratio by prestressing platelet inclusions or controlling the surface rigidity of a matrix of 2D materials. KEYWORDS: auxetic materials, 2D materials, 2D membranes, negative Poisson ratio

M

lattices possessing small bending rigidity with free boundary conditions.12−14 Experimental realizations of these predictions for microscopic surfaces, which could be used as fillers in composites,15,16 have shown limited success. In principle, graphene should exhibit auxetic behavior, as it can be modeled as a 2D self-avoiding fixed-connectivity membrane (e.g., some previous studies have used finite-element-based models to uncover slightly auxetic behavior of graphene17 and hexagonal boron nitride,18 but only for strains up to 0.01%). However, the significant bending rigidity of graphene suppresses most out-ofplane fluctuations leading to no measurable transverse expansion in detailed ab initio simulations.19,20 The difficulties in manipulating the Poisson ratio of graphene illustrates an important problem that has been overlooked in the more general literature of 2D entropic surfaces: How do we design materials with a negative Poisson’s ratio over a range of relevant strains, and what properties affect how persistent these

aterials with a negative Poisson’s ratio, known as auxetics,1 expand in transverse directions under an applied uniaxial tensile strain. This property is highly desirable in several applications, such as impact mitigation, where material failure is caused by a thinning of the materials in the impact region, and for sealants, where it is desirable for materials to expand to fill regions as anisotropic pressure is applied. Recently, the auxetic nature of graphene oxide2 was also reported to enable its use in water desalination membranes.3 A negative Poisson ratio can be engineered at the macroscopic scale with carefully designed lattice structure,4−6 using auxetic foams,7 or using carbon nanotube sheets,8,9 but there are limited ways to create molecular structures with a tunable Poisson’s ratio. At the microscopic scale, it is well established that fluctuating membranes should, in the zero-strain limit, possess a negative Poisson ratio.10 Applying anisotropic in-plane tensile strain to a self-avoiding two-dimensional (2D) membrane suppresses out-of-plane roughness, which allows the material to expand in the transverse direction. This effect has been confirmed with simulations in the zero-strain limit for simple interconnected beads with periodic boundary conditions11 and for triangular © 2016 American Chemical Society

Received: April 13, 2016 Accepted: July 7, 2016 Published: July 18, 2016 7542

DOI: 10.1021/acsnano.6b02512 ACS Nano 2016, 10, 7542−7549

Article

www.acsnano.org

Article

ACS Nano properties are? The increasing complexity of nanoscale devices suggests that there may be many new materials which could be engineered to exhibit auxetic properties. For example, possible routes to surfaces with desirable functionality include 2D DNA origami surfaces,21 interlocking molecular structures,22 and 2D materials more flexible than graphene. In fact, a recent study proposed the introduction of atomic vacancies into the graphene lattice to promote out-of-plane structure that could lead to auxetic behavior.23 In this work, we present simulations of the strain-dependent behavior of 2D self-avoiding fixed-connectivity surfaces and compare the results with a simple model studied in the zero strain limit in the previous literature. We show that membranes with free boundary conditions and without an explicit bending rigidity actually possess a positive Poisson ratio at zero-strain due to a non-zero mean curvature of the surface. The possibility for suppression of this curvature leads to auxetic behavior (possession of a negative Poisson ratio), and this persists over a large range of applied strains for small surfaces. Larger surfaces or those with more significant bending rigidity show a smaller window for a negative Poisson ratio. At very large strains, the membrane connectivity becomes limiting, and the membrane begins to contract, placing an upper limit on the window for auxetic behavior. Previous studies have considered in detail the behavior of 2D sheets in the limit of zero applied strain or with uniform isotropic strain. These studies often described the out-of-plane fluctuations in 2D membranes using a mean-field theoretical approach. Note that, we will show later that these theoretical approaches fail to quantitatively capture the simulated mechanical behavior of 2D entropic surfaces, and instead we propose a semiempirical descriptor for the out-of-plane fluctuations. The mean-field theoretical approach is formulated in terms of the Hamiltonian for a strained continuum elastic surface with small out-of-plane fluctuations in the harmonic approximation, coupled to an external stress source:25 F[u , z , ταβ ] =

1 2

This definition is valuable for considering how the local strain contributes to the macroscopic strain for the entire membrane and is related to the angle between the local surface tangent and the mean plane of the membrane. The above system has been studied in detail, especially in the 24 case of no external strain, i.e., uext Key measurables that αβ = 0. have been described for 2D self-avoiding membranes with fixed-connectivity based on this Hamiltonian include the rootmean-square (RMS) out-of-plane fluctuations ⟨z2⟩1/2, the radius of gyration Rg, and the induced bending rigidity from selfavoidance κ. The most interesting property of these systems is the prediction of a universal Poisson ratio σ, defined as the ratio of transverse to longitudinal strain. In this study, we consider two macroscopic definitions of the Poisson ratio:

∫ d2r[κ(∇2 z)2 + 2μuαβ 2 + λuαα 2 + ταβuαβ]

Here, u is the internal strain tensor, z is the membrane height, ext ταβ is the in-plane stress tensor (ταβ = λδαβuext αβ + 2μuαβ ) (α and ext β indicate the coordinate axes), u is the external strain tensor, κ is the bending rigidity, λ is the first Lamé coefficient, μ is the shear modulus, and r is a position vector indicating the location on the 2D sheet. Note that the Einstein summation convention is followed by which each repeated subscript is assumed to have a summation over it. This theoretical model is considerably more simplified than the system considered later using molecular simulations. The internal strain tensor of a thin membrane, to first order, is: 1 (∂αuβ + ∂βuα + ∂αz ∂βz) 2

(2)

where, ∂α indicates the partial derivative with respect to coordinate α and uα denotes the displacement in the α direction due to the strain. The internal strain is important for describing the state of the system. We can also define a spatial 2D strain for the membrane projected onto the mean surface: P uαβ =

1 (∂αuβ + ∂βuα) 2

∂uyy uxx ∂uxx uyy

σM = −

δW /W δL /L

(4)

The first definition σL is based on the strain definition presented in eq 2 and represents the local Poisson ratio of the membrane surface including out-of-plane fluctuations. The macroscopic Poisson ratio σM is defined as the negative of the ratio of the center-of-mass distance between the membrane edges in the longitudinal/strained direction L and the transverse/free direction W. Hence, σM is the Poisson ratio we would observe for the entire sheet if we only probed it through interactions at the membrane edges and is the most important for consideration of practical experiments. It is important to note that all of these definitions are highly straindependent, as their properties depend on out-of-plane fluctuations that can be suppressed or enhanced through external strain. Interestingly, the local Poisson ratio σL has been predicted to be approximately −1/3 in the limit of no external strain. This follows analytical theory using renormalization group26 and numerical simulations of finite-sized membranes.11,13 The Poisson ratio in this limit of zero strain appears to be nearly universal and independent of the precise connectivity geometry or tethering, under the presence of periodic boundary conditions or explicit bending rigidity. In this study, we report strain-dependent calculations of the Poisson ratio in 2D self-avoiding fixed-connectivity membranes. These membranes exhibit a transition from previously reported universal limits to a positive Poisson ratio at a size-dependent critical strain. Further, we show that this critical strain decreases with increasing membrane size, suggesting that maintaining a small membrane size is actually beneficial for maintaining a negative Poisson ratio over a large regime. In addition, we show that other experimentally-relevant measurables, such as the RMS out-of-plane fluctuations, undergo nonlinear transformations at similar strains. We show that this transition is necessary, as the behavior transitions from being dominated by entropic fluctuations of the membrane to enthalpic interactions in the strained lattice.

(1)

uαβ =

σL = −

RESULTS AND DISCUSSION Simulation Methods. We performed molecular dynamics (MD) simulations of self-avoiding fixed-connectivity membranes in the canonical (NVT) ensemble with varying external anisotropic strain. For easy comparison to previous studies, we adopted interaction potentials from a previous study in the

(3) 7543

DOI: 10.1021/acsnano.6b02512 ACS Nano 2016, 10, 7542−7549

Article

ACS Nano literature.11 N × N lattices of beads (N = 28, 40, 100, 150) were connected in a triangular mesh (Figure 1). All beads interacted UFENE

⎧ ⎡ ⎛ r ⎞2 ⎤ 1 2 ⎪ ⎢ ⎪− kR 0 ln 1 − ⎜ ⎟ ⎥ r ≤ R 0 ⎢⎣ ⎝ R 0 ⎠ ⎥⎦ =⎨ 2 ⎪ ⎪ r > R0 ∞ ⎩

(6)

As in previous studies,11 parameters were fixed to be k = 6ε/ σ , R0 = 1.5σ, and T = ε/kBT. We chose these force field scalings to relate to previous quantitative simulations in the literature. Changing the force field outside of these simple scalings would change the results quantitatively, but would require orders of magnitude more simulations to map out the effect of every force field parameter on our final results. However, we do not believe that the choice of these scalings would affect our conclusions qualitatively. Periodic boundary conditions used previously11 were not amenable to the computation of anisotropic strain, so the edges of the lattice were kept unclamped to match the experimental conditions for such systems. Forces were applied to the edge beads in the longitudinal direction. The same force was applied to each bead along each edge, with opposite edges having opposite forces. The direction of the applied forces was updated periodically based on the orientation of the membrane. Self-avoidance in fluctuating membranes leads to an induced bending rigidity which stabilizes the flat phase of these membranes, preventing crumpling. However, many physical membranes of 2D materials like graphene have explicit bending rigidity as well. To investigate the effect of the surface rigidity on out-of-plane fluctuations and the Poisson ratio, three quadratic improper dihedral terms were added for each bead (Figure 1), of the form: 2

Figure 1. Example geometry for a 6 × 6 mesh. Forces were added to edge beads in the longer (longitudinal) direction. Example improper dihedrals, used to increase the bending rigidity of the sheet, are also indicated for a single central atom, as shown with the red/green/blue connections.

via a truncated self-avoiding Lennard-Jones (LJ) potential of the form: ⎧ ⎡⎛ ⎞12 ⎛ ⎞6 1⎤ σ σ ⎪ ⎪ 4ε⎢⎜ ⎟ − ⎜ ⎟ + ⎥ r ≤ 21/6σ ⎝r⎠ 4⎦ ULJ = ⎨ ⎣⎝ r ⎠ ⎪ ⎪ ⎩ 0 r > 21/6σ

(5)

Bonded beads were connected with a finite extensible nonlinear elastic (FENE)27 bond of the form:

V (χ ) = kimp(χ − χ0 )2

(7)

Figure 2. Surface behavior under longitudinal extension for 40 × 40 membranes with varying bending rigidity. (A) Transverse strain as a function of applied longitudinal strain, showing an initial shrinking of the membrane in the transverse direction followed by a maximum. Dashed lines are approximations of the local gradient for each set of calculations. (B) Cartoon of the effect of applying longitudinal strain for a sheet without explicit rigidity, showing that longitudinal extension causes the surface to flatten in the longitudinal direction, while becoming more curved in the transverse direction. (C) Mean curvature as a function of the applied longitudinal strain for surfaces of various bending rigidity. Solid lines indicate the mean surface curvature in the longitudinal direction, while dashed lines indicate mean curvature in the transverse direction. Open symbols are averages for ensemble simulations (constant applied force), and filled symbols are averages taken by subsampling ensemble results. Lines are smoothing splines included only as visual guides. 7544

DOI: 10.1021/acsnano.6b02512 ACS Nano 2016, 10, 7542−7549

Article

ACS Nano

Figure 3. Strain-dependent Poisson ratio for several surface sizes and relation to the out-of-plane fluctuations. (A) Macroscopic Poisson ratio as a function of applied strain for various surface sizes, showing an initially positive value, a local minima at ∼10%, followed by a return to positive values. (B) Transverse extension of the same surfaces. Circles represent ensemble averages for each simulation condition (i.e., constant end forces), with error bars generated using the bootstrap method. Small points represent averages of collections from each ensemble, showing variation within each condition. Solid lines are smoothing splines included for visual guides. Dashed lines are approximations of the local gradient for each set of calculations. (C) Mean curvature as a function of the applied longitudinal strain for surfaces of various sizes (kimp = 0). Solid lines indicate the mean surface curvature in the longitudinal direction, while dashed lines indicate mean curvature in the transverse direction. Open symbols are averages for ensemble simulations (constant applied force), and filled symbols are averages taken by subsampling ensemble results. Lines are smoothing splines included only as visual guides.

where kimp represents the strength of the added rigidity, χ is the dihedral angle, and χ0 = π is the equilibrium angle ensuring that the equilibrium surface is flat. kimp was varied from [0ε,1000ε] for the N = 40 sheet. Details about the precise simulation method and equilibration protocols are contained in the Methods section. Effect of Bending Rigidity on Auxetic Behavior. In our simulations, all systems initially contracted in response to longitudinal strain, in contrast to previous studies, due to the non-zero mean curvature of the membrane in absence of explicit bending rigidity or periodic boundary conditions. Accordingly, anisotropic straining of the membrane leads to an initial transverse contraction of the surface, as shown in Figure 2A for an example case of a 40 × 40 system, for the various imposed bending rigidities. At approximately 5% strain for the 40 × 40 surface without any applied bending rigidity (kimp = 0), the transverse strain reaches a minimum. After that point, further straining the membrane leads to an increase in the transverse strain. In this region, the surface expands in the transverse direction as entropic out-of-plane fluctuations are suppressed, thereby the Poisson ratio (i.e., negative of the ratio of transverse to longitudinal strain) is negative, as explained in detail later. Further straining causes the membrane to lose its auxetic property (at about 20% strain in Figure 2A), leading to a transition from a negative Poisson ratio to a positive one. Adding an explicit bending rigidity to the surface through the application of an improper dihedral with force constant kimp leads to a reduction in the maximum in the transverse extension with respect to the applied longitudinal strain. Increasing kimp above 100ε suppresses the initial contraction by reducing the mean curvature of the unstrained membrane, resulting in a positive Poisson ratio at zero strain. Although finite-size self-avoiding 2D membranes are asymptotically flat (they do not crumple entropically into

spheres), they do contain non-zero mean curvature over the scale of the membrane, as illustrated by the cartoon of a 40 × 40 membrane without explicit bending rigidity, in Figure 2B. This non-zero mean curvature of the membrane arises due to the random thermal (entropic) fluctuations of the atoms and depends on the bending rigidity of the membrane. To quantify the membrane curvature, the principal curvatures of the membrane in the longitudinal and transverse directions, calculated using the surfature module in MATLAB, were spatially averaged over the entire surface. From Figure 2B, at uxx = 0, there is a non-zero mean curvature in both longitudinal and transverse directions. At small extensions, the curvature in the longitudinal direction is quickly suppressed, but the transverse curvature increases slightly. For large extensions, curvature in both directions is suppressed, and the surface becomes more flat. The extent of this curvature depends on the surface parameters (bead size, bonding, etc.), the size of the sheet, and whether there is an explicit bending rigidity present. Previous simulation efforts did not observe this effect because periodic boundary conditions were imposed11 or significant bending rigidity was maintained that suppressed this curvature.12−14 Figure 2C depicts the mean curvature in the transverse (solid lines) and longitudinal (dashed lines) directions for membranes with varying explicit bending rigidity after averaging out small-scale spatial variation using a Gaussian spatial filter. For all membranes, applying longitudinal strain to the membrane immediately leads to a reduction in curvature in the longitudinal direction, as shown by the solid lines in Figure 2C. However, for the transverse curvature (dashed lines in Figure 2C), there exist two different trends: (i) for membranes with high bending rigidity, the transverse curvature also decreases with initial strain, while (ii) for membranes with low bending rigidity, the transverse curvature actually increases slightly with initial strain. Further longitudinal straining for low 7545

DOI: 10.1021/acsnano.6b02512 ACS Nano 2016, 10, 7542−7549

Article

ACS Nano bending rigidity membranes leads to a decrease in the transverse curvature, but at a rate slower than in the longitudinal direction. Therefore, adding bending rigidity by adjusting kimp, the force constant for the improper dihedrals as described above, reduces the out-of-plane fluctuations and mean curvature of the membrane. Doing so reduces both the minimum Poisson ratio and the longitudinal strain at which the system transitions back to a positive Poisson ratio. These results show that the explicit bending rigidity of a 2D material plays a crucial role in determining the auxetic behavior of a material, thereby providing avenues for introducing a negative Poisson’s ratio by reducing the explicit bending rigidity. Effect of Membrane Size on Auxetic Behavior. The strain-dependent Poisson ratio for the various membrane sizes considered, shown in Figure 3A, was calculated from the spline fits of the strain−strain relations in Figure 3B. The Poisson ratio is positive at zero-strain for membranes of all sizes with no explicit bending rigidity, with larger membranes having larger zero-strain Poisson ratio. This is significantly different from previous estimates of the zero-strain Poisson ratio for the same model system due to the difference in boundary conditions and resulting non-zero mean curvature in the membrane, as discussed above. Most of the sheets transition to a negative Poisson’s ratio at approximately 5% applied longitudinal strain and reach a minimum Poisson ratio significantly more negative than those reported in the zero-strain limit, with larger membranes having lower minimum Poisson ratios. This continues until the membrane is essentially flat in the longitudinal direction, and further strain forces neighboring beads closer together, leading to a transverse contraction and a return to a positive Poisson’s ratio. Note that, because the 2D membranes considered here are anisotropic (even within the xy plane), they are not subject to the classical thermodynamic limit of −1 < ν < 0.5, which is true only for isotropic materials. Indeed, in Figure 3A, one can see that the Poisson ratio exceeds 0.5 at extremely high strains. Figure 3C depicts the mean curvature in the transverse (solid lines) and longitudinal (dashed lines) directions for membranes with varying sizes (kimp = 0) after averaging out small-scale spatial variation using a Gaussian spatial filter. For all membrane sizes, applying longitudinal strain to the membrane immediately leads to a reduction in curvature in the longitudinal direction, as shown by the solid lines in Figure 3C. On the other hand, the transverse curvature actually increases slightly with initial strain and then decreases at higher strain. The above results show that the strain-dependent Poisson’s ratio can be engineered by adjusting the membrane size (Figure 3) or by adjusting the bending rigidity of the membrane (Figure 2). In both situations, the change in Poisson ratio is a direct effect of adjusting the out-of-plane fluctuations at low strains that allow the system to expand when flattened. The return to a positive Poisson ratio occurs in all membranes when the mean out-of-plane fluctuations is reduced below a critical value of approximately ⟨z2⟩1/2 ≈ 0.4, as shown in Figure 4, after which there are essentially no more out-of-plane entropic fluctuations to remove. This continues until most out-of-plane fluctuations have been reduced (⟨z2⟩1/2 < 0.4) corresponding to approximately 10−20% applied longitudinal strain, depending on system size. As applied strain increases beyond this critical level, interbead stretching takes over, leading to a contraction. This effect can be seen more clearly in Figure 4; as the out-ofplane fluctuations are reduced from around ⟨z2⟩1/2 = 2.1, the

Figure 4. Transverse strain as a function of the out-of-plane fluctuations, demonstrating that the minimum in the Poisson ratio is related to the out-of-plane fluctuations. When all of these fluctuations have been removed, i.e. ⟨z2⟩ → 0, the transverse direction begins to contract.

transverse strain increases, until the critical point at which ⟨z2⟩1/2 = 0.4. At very large applied longitudinal strains (⟨z2⟩1/2 < 0.4 in Figure 4), the deformation of the triangular lattice leads to a reduction in the transverse membrane extent. This is a simple geometric effect arising due to the surface connectivity and depends on the precise connectivity and interbead potential. For an infinitely stiff network in the unstrained direction (i.e., nondeformable interbead distances in the transverse direction), it is evident that the transverse membrane extent will reduce to accommodate for the increased longitudinal extent. Comparison of Results to Previous Theory. For comparison to previous studies, we also report the size scaling of the radius of gyration Rg ∝ Lν, ν = 1.02 ± 0.0004, and the size-scaling roughness exponent for the out-of-plane fluctuations, ⟨z2⟩ ∝ L2ζ, ζ = 0.79 ± 0.07. The roughness scaling is slightly larger than previous studies11,12 due to the difference in boundary conditions. The initial decrease in RMS height (Figure 4) can be understood using continuum elastic theory in the harmonic approximation, represented by the Hamiltonian in eq 1.25 The RMS height fluctuation follows as: ⟨z 2⟩ =

∫ dqG (22ππ|q)| , 0

G0(q) =

2

kBT 2

(

2

ext ext q κq + λuαα + 2μuαβ

qαqβ |q|2

)

(8)

where G0 is the Fourier component of the height−height correlation function and uext αβ is the external strain applied to the membrane. For purely longitudinal external strain (uyy,uxy = 0), G0 simplifies to G0(q) =

kBT 2

2

q (κq +

ext uxx (λ

ext + μ) + μuxx cos(2θ ))

(9)

where θ is the angle between the Fourier mode q and the x-axis. Integrating over the magnitude of the Fourier component q in eq 8, starting at a minimum Fourier component inversely proportional to the system size qmin = 2π/L yields 7546

DOI: 10.1021/acsnano.6b02512 ACS Nano 2016, 10, 7542−7549

Article

ACS Nano

Figure 5. (A) Straining the membrane initially decreases the out-of-plane fluctuations, until approximately 10% strain. Theory for each membrane size is a fit to eq 10, demonstrating the theory is a poor approximation at large strains. (B) Relation of the maximum transverse extension to an empirically-derived descriptor, (⟨z2⟩1/2/R0.88 g )uxx=0, which describes the amount of out-of-plane fluctuations in the unstrained state that can be suppressed by applying longitudinal strain. The empirical exponent of 0.88 for Rg collapses the results obtained by changing the membrane size or adding explicit bending rigidity to the system. ⟨z 2⟩ =

=

∫0 ∫0

2π

dθ

∫q

∞

dq

min

2π

q kBT ext ext (2π )2 q2(κq2 + uxx (λ + μ) + μuxx cos(2θ))

Table 1. Theoretical Bending Rigidity As a Function of the Improper Dihedral, For a Given Membrane Size of N = 40

⎡ ⎛ u ext (λ + μ) u extμ cos(2θ) ⎞ ⎤ ⎟⎥ + xx 2 ⎢ kBT log⎜1 + xx 2 qminκ qminκ ⎠⎥ ⎝ dθ ⎢ 2 ext ⎥ ⎢ 8π uxx (λ + μ + μ cos(2θ)) ⎥ ⎢ ⎦ ⎣

(10)

Equation 10 allows for a theoretical estimate of the mean out-of-plane fluctuations ⟨z2⟩, given the values of the effective bending rigidity and the Lamé coefficients. For this purpose, eq 10 was used to derive two important parameter groupings, evaluated in the limit of zero strain:

improper dihedral, kimp (N = 40)

bending rigidity, κ

0ε 10−1ε 100ε 101ε 102ε

0.5352ε 0.6339ε 1.3218ε 4.3469ε 9.5203ε

Table 2. Theoretical Bending Rigidity As a Function of the Membrane Size, For No Applied Improper Dihedral

d log⟨z 2⟩ L2(λ + μ) 2 lim z , lim ⟨ ⟩ = = − ext ext uxx →0 16π 3κ uxxext → 0 duxx 8π 2κ kBTL2

(11)

Using these two zero-strain limits allowed for the estimation of the effective bending rigidity and the combination λ + μ, under the assumption that λ = μ. With these values, eq 10 was integrated numerically over a range of imposed external strains for various membrane sizes as shown in Figure 5A, using dashed lines. It can be seen that the theory appears to deviate from the simulation results at strains of just a few percent. This is because the mean-square roughness decreases faster than the theory suggests, due to a reduction in the non-zero mean curvature of the unstressed membrane with free boundary conditions, which is not accounted for in the theory. At sufficiently extensive strains, the effective bending rigidity decreases, leading to deviation from the theoretical limit of constant rigidity and elastic coefficients. The induced bending rigidity κ from the self-avoiding character necessarily weakens at increasing strain to become anisotropic, causing the membrane to collapse in the transverse direction. Table 1 lists the calculated theoretical bending rigidities as a function of the membrane size, N (for kimp = 0), and Table 2 lists the theoretical bending rigidities as a function of the improper dihedral, kimp (for N = 40). Data in Tables 1 and 2 indicate that the induced bending rigidity increases with increasing membrane size and with increasing improper dihedral, which is consistent with the smaller windows for negative Poisson ratio with increasing membrane size and with increasing improper dihedral, as discussed above.

membrane size (kimp = 0)

bending rigidity, κ

28 × 28 40 × 40 100 × 100 150 × 150

0.4961ε 0.5352ε 0.8048ε 0.9908ε

The transition from negative to positive Poisson ratio at large strains can be correlated to the out-of-plane fluctuations ⟨z2⟩1/2 and the radius of gyration Rg of the membrane at zero-strain. A simple metric [⟨z2⟩1/2/R0.88 g ]uxx=0 predicts the maximal transverse extension under the addition of bending rigidity or varying system size, as shown in Figure 5B. The scaling exponents for ⟨z2⟩ and Rg given above suggest that this transition point should scale as ⟨z2⟩1/2/R0.88 ∼ L−0.11, suggesting that even membranes g sizes several orders of magnitude larger will have negative Poisson ratios over substantial strain ranges. This could be verified with further work, but simulating larger membranes takes significantly more computational resources. Comparison to Real 2D Materials. The bending rigidities of real, covalently-bonded 2D materials, such as hexagonal boron nitride (hBN), graphene, and molybdenum disulfide (MoS2) are 0.98, 1.18, and 9.61 eV, respectively.31,32 Note that these materials possess a honeycomb lattice (which is a combination of two triangular lattices), as opposed to the single triangular lattice considered here, and therefore the bending rigidities cannot be directly compared. Our choice of the triangular lattice was motivated by its use in previous simulations studies. Nevertheless, typical values of the improper dihedral potential for the above-mentioned 2D materials will likely lie beyond 10ε. Therefore, our simulations put into 7547

DOI: 10.1021/acsnano.6b02512 ACS Nano 2016, 10, 7542−7549

Article

ACS Nano

The direction of the forces was updated every 1000 steps based on the vector defined by the center of each edge, after projection onto a 2D surface by using the two largest components of the SVD of the bead positions. The system was integrated using a time step of 10−3σ m/ε . Note that, in our work, the lattice spacing, the bead mass, and the energy scale were all set to unity and represent the scales in the system. Simulations were performed for systems with size N = 28, 40, 100, 150. Forces were applied to reach longitudinal extension from approximately 0−60%. For each simulation, the first 107 steps were discarded for equilibration. Simulation statistics were collected every 105 steps thereafter, which was the minimum to provide uncorrelated samples as determined by the autocorrelation function. For larger systems, multiple replicas with random initial velocities were used to generate extra samples. For most samples, >1010 steps were achieved, resulting in ∼105 independent samples. In all cases, error estimates on variables were generated using the bootstrap method with 103 samples.

context what is actually realizable in experiments. Indeed, while our simulations reveal that previous theoretical work was unable to uncover the finite strain auxetic behavior of 2D entropic surfaces, the bending rigidity dependence of the current results suggests that, going forward, if elusive auxetic behavior is desired, innovative approaches like engineering defects or designing 2D-polymer-like structures will be required.

CONCLUSIONS In summary, we have performed large-scale, GPU-accelerated MD simulations to uncover the strain-dependent Poisson ratio of entropic 2D surfaces, with varying sizes and bending rigidities. Unlike previous theoretical efforts, we have considered systems with unclamped boundary conditions, under conditions of non-zero anisotropic strain. These simulations show that exotic auxetic (negative Poisson ratio) properties of entropic 2D surfaces are more interesting than the universal zero-strain limits proposed in previous literature and that the properties persist over a surprisingly large range of applied strains. The most-negative Poisson ratios observed in this study are considerably more negative than the zero-strainlimit simulations suggested. Adjusting the rigidity and the size of the membrane both affect the strain range over which auxetic behavior persists, and this suggests that there should be noticeable strain-dependent effects for practical surfaces. Sufficiently stiff membranes can completely suppress these regions of auxetic behavior, explaining why attempts to use graphene membranes as auxetic materials have failed. Our simulations also uncover that larger membranes possess a smaller window for auxetic behavior. Finally, we show that previous simplified mean-field theoretical approaches, formulated in the harmonic approximation, are not able to capture the realistic mechanical behavior of 2D membranes. Instead, we provide a metric for the persistence range based on a correlation consisting of the out-of-plane fluctuations and the radius of gyration. Going forward, these results suggest that another path toward making materials with a specific Poisson ratio is to prestrain the materials to the desired point. This could be achieved by embedding the membranes in a second matrix that can be adjusted.

AUTHOR INFORMATION Corresponding Author

*E-mail: [email protected]. Present Address ‡

Department of Chemical Engineering, Stanford University, Stanford, CA 94305, United States Notes

The authors declare no competing financial interest.

ACKNOWLEDGMENTS This work was supported by a grant from the Army Research Office, supported via award no. 64655-CH-ISN to the Institute for Soldier Nanotechnologies. Z.W.U. acknowledges Department of Energy CSGF support (DOE grant DE-FG0297ER25308). This work used the Extreme Science and Engineering Discovery Environment (XSEDE), which is supported by National Science Foundation grant no. ACI1053575. REFERENCES (1) Evans, K. E.; Alderson, A. Auxetic Materials: Functional Materials and Structures from Lateral Thinking! Adv. Mater. 2000, 12, 617−628. (2) Talyzin, A. V.; Solozhenko, V. L.; Kurakevych, O. O.; Szabó, T.; Dékány, I.; Kurnosov, A.; Dmitriev, V. Colossal Pressure-Induced Lattice Expansion of Graphite Oxide in the Presence of Water. Angew. Chem., Int. Ed. 2008, 47, 8268−8271. (3) Boukhvalov, D. W.; Katsnelson, M. I.; Son, Y.-W. Origin of Anomalous Water Permeation through Graphene Oxide Membrane. Nano Lett. 2013, 13, 3930−3935. (4) Yang, W.; Li, Z.-M.; Shi, W.; Xie, B.-H.; Yang, M.-B. Review on Auxetic Materials. J. Mater. Sci. 2004, 39, 3269−3279. (5) Caddock, B. D.; Evans, K. E. Microporous Materials with Negative Poisson’s Ratios. I. Microstructure and Mechanical Properties. J. Phys. D: Appl. Phys. 1989, 22, 1877−1882. (6) Alderson, A.; Evans, K. E. Molecular Origin of Auxetic Behavior in Tetrahedral Framework Silicates. Phys. Rev. Lett. 2002, 89, 225503. (7) Lakes, R. Foam Structures with a Negative Poisson’s Ratio. Science 1987, 235, 1038−1040. (8) Hall, L. J.; Coluci, V. R.; Galvao, D. S.; Kozlov, M. E.; Zhang, M.; Dantas, S. O.; Baughman, R. H. Sign Change of Poisson’s Ratio for Carbon Nanotube Sheets. Science 2008, 320, 504−507. (9) Coluci, V. R.; Hall, L. J.; Kozlov, M. E.; Zhang, M.; Dantas, S. O.; Galvão, D. S.; Baughman, R. H. Modeling the Auxetic Transition for Carbon Nanotube Sheets. Phys. Rev. B: Condens. Matter Mater. Phys. 2008, 78, 115408. (10) Nelson, D. R.; Piran, T.; Weinberg, S. Statistical Mechanics of Membranes and Surfaces; World Scientific Publishing Co.: Singapore, 2004.

METHODS MD simulations were performed on graphical processing units (GPUs) in the HOOMD-BLUE 1.0 package,28−30 enabling several orders of magnitude more simulation capability than previous studies. This capability was necessary to analyze surface behavior under a variety of longitudinal strains. The use of GPUs placed an upper bound on the system size that could be considered, since the volume of the simulation cell is directly proportional to the number of regional neighbor lists needed for the simulation. To maximize possible system size, the membrane was placed in a simulation cell with a large aspect ratio (e.g., a width/length proportional to the area of the membrane, and fixed cell height of 100 bead radii corresponding to relatively small out-of-plane fluctuations). To prevent the membrane from rotating out of plane and interacting with neighbors via the periodic boundary conditions, every 103 time steps the membrane and velocities were rotated such that the smallest basis vector from the singular value decomposition (SVD) was aligned to the z-axis, keeping the membrane extent in the x- and y-directions with out-of-plane fluctuations in the z-direction. Temperature control was implemented with the default HOOMD NVT integrator with a time scale of 103 simulation steps. The momentum of the simulation was zeroed every 103 steps. 7548

DOI: 10.1021/acsnano.6b02512 ACS Nano 2016, 10, 7542−7549

Article

ACS Nano (11) Zhang, Z.; Davis, H. T.; Kroll, D. M. Molecular Dynamics Simulations of Tethered Membranes with Periodic Boundary Conditions. Phys. Rev. E: Stat. Phys., Plasmas, Fluids, Relat. Interdiscip. Top. 1996, 53, 1422−1429. (12) Bowick, M. J.; Catterall, S. M.; Falcioni, M.; Thorleifsson, G.; Anagnostopoulos, K. N. The Flat Phase of Crystalline Membranes. J. Phys. I 1996, 6, 1321−1345. (13) Bowick, M.; Cacciuto, A.; Thorleifsson, G.; Travesset, A. Universal Negative Poisson Ratio of Self-Avoiding Fixed-Connectivity Membranes. Phys. Rev. Lett. 2001, 87, 148103. (14) Bowick, M. J.; Cacciuto, A.; Thorleifsson, G.; Travesset, A. Universality Classes of Self-Avoiding Fixed-Connectivity Membranes. Eur. Phys. J. E: Soft Matter Biol. Phys. 2001, 5, 149−160. (15) Wei, G.; Edwards, S. F. Poisson Ratio in Composites of Auxetics. Phys. Rev. E: Stat. Phys., Plasmas, Fluids, Relat. Interdiscip. Top. 1998, 58, 6173−6181. (16) Suter, J. L.; Coveney, P. V. Materials Properties of Clay Nanocomposites: Onset of Negative Poisson Ratio in Large-Scale Molecular Dynamics Simulation. Soft Matter 2009, 5, 3896. (17) Scarpa, F.; Adhikari, S.; Srikantha Phani, A. Effective Elastic Mechanical Properties of Single Layer Graphene Sheets. Nanotechnology 2009, 20, 065709. (18) Boldrin, L.; Scarpa, F.; Chowdhury, R.; Adhikari, S. Effective Mechanical Properties of Hexagonal Boron Nitride Nanosheets. Nanotechnology 2011, 22, 505702. (19) Liu, F.; Ming, P.; Li, J. Ab Initio Calculation of Ideal Strength and Phonon Instability of Graphene under Tension. Phys. Rev. B: Condens. Matter Mater. Phys. 2007, 76, 064120. (20) Gui, G.; Li, J.; Zhong, J. Band Structure Engineering of Graphene by Strain: First-Principles Calculations. Phys. Rev. B: Condens. Matter Mater. Phys. 2008, 78, 075435. (21) Han, D.; Pal, S.; Nangreave, J.; Deng, Z.; Liu, Y.; Yan, H. DNA Origami with Complex Curvatures in Three-Dimensional Space. Science 2011, 332, 342−346. (22) Ravirala, N.; Alderson, A.; Alderson, K. L. Interlocking Hexagons Model for Auxetic Behaviour. J. Mater. Sci. 2007, 42, 7433−7445. (23) Grima, J. N.; Winczewski, S.; Mizzi, L.; Grech, M. C.; Cauchi, R.; Gatt, R.; Attard, D.; Wojciechowski, K. W.; Rybicki, J. Tailoring Graphene to Achieve Negative Poisson’s Ratio Properties. Adv. Mater. 2015, 27, 1455−1459. (24) Falcioni, M.; Bowick, M. J.; Guitter, E.; Thorleifsson, G. The Poisson Ratio of Crystalline Surfaces. Europhys. Lett. 1997, 38, 67−72. (25) Roldán, R.; Fasolino, A.; Zakharchenko, K. V.; Katsnelson, M. I. Suppression of Anharmonicities in Crystalline Membranes by External Strain. Phys. Rev. B: Condens. Matter Mater. Phys. 2011, 83, 174104. (26) Le Doussal, P.; Radzihovsky, L. Self-Consistent Theory of Polymerized Membranes. Phys. Rev. Lett. 1992, 69, 1209−1212. (27) Kremer, K.; Grest, G. S. Dynamics of Entangled Linear Polymer Melts: A Molecular-Dynamics Simulation. J. Chem. Phys. 1990, 92, 5057. (28) Anderson, J. A.; Lorenz, C. D.; Travesset, A. General Purpose Molecular Dynamics Simulations Fully Implemented on Graphics Processing Units. J. Comput. Phys. 2008, 227, 5342−5359. (29) Glaser, J.; Nguyen, T. D.; Anderson, J. A.; Lui, P.; Spiga, F.; Millan, J. A.; Morse, D. C.; Glotzer, S. C. Strong Scaling of GeneralPurpose Molecular Dynamics Simulations on GPUs. Comput. Phys. Commun. 2015, 192, 97−107. (30) HOOMD-Blue, http://codeblue.umich.edu/hoomd-blue. (31) Singh, S. K.; Neek-Amal, M.; Costamagna, S.; Peeters, F. M. Thermomechanical Properties of a Single Hexagonal Boron Nitride Sheet. Phys. Rev. B: Condens. Matter Mater. Phys. 2013, 87, 184106. (32) Jiang, J.-W.; Qi, Z.; Park, H. S.; Rabczuk, T. Elastic Bending Modulus of Single-Layer Molybdenum Disulfide (MoS2): Finite Thickness Effect. Nanotechnology 2013, 24, 435705.

7549

DOI: 10.1021/acsnano.6b02512 ACS Nano 2016, 10, 7542−7549