text:
# Growth, microstructure, and failure of crazes in glassy polymers

## I Introduction

The failure of glassy polymers such as polystyrene (PS) or polymethylmethacrylate (PMMA) under external stresses occurs either through shear deformation or through crazing Haward1997 ; Ward1983 . While shear yielding occurs essentially at constant volume, crazing has a strong dilational component, and the volume of the material increases to several times its original value before catastrophic fracture occurs. Crazing is a failure mechanism unique to entangled polymeric materials and usually precedes a crack tip (see Fig. 1). The fundamental and technological importance of crazes is that they are in part responsible for the large fracture energy $`G_c`$ of polymer glasses Wool1995 ; Pocius1997 ; Brown1991 ; Rottler2002a  that makes them useful load-bearing materials. They control the crack tip advance and require a large amount of energy dissipation up to the point of catastrophic failure. Crazes can reach several $`\mu m`$ in width and consist of an intriguing network of fibrils and voids that spans the entire deformed region.

Despite the frequent appearance of crazes, there is still comparatively little theoretical understanding about the conditions and mechanisms of craze nucleation, growth and ultimate breakdown Haward1997 ; Ward1983 ; Wool1995 ; Pocius1997 ; Krupenkin1999 ; Kramer1990 . In this paper, we present an extensive set of nonequilibrium molecular dynamics (MD) simulations that address these various phenomena. In this approach, polymers are modeled on a coarse-grained scale that takes into account van-der-Waals (vdW) and covalent interactions without specific reference to chemical detail. The effect of chain length, temperature $`T`$, widening velocity $`v`$ and vdW interaction strength on the craze structure can be studied over a wide range of parameters. The molecular simulations allow insight into microscopic details not accessible to experiments and offer an opportunity to test and develop theoretical models of crazing.

A fundamental limitation on molecular level treatments is of course the finite system size. The largest volumes accessible at present are $`100nm^3`$, while the craze spans many $`\mu m`$. We are thus limited to a study of craze widening in a small representative region and cannot include e.g. the entire crack tip. Craze tip advance processes Kramer1990  are beyond the scope of the present work.

Several aspects of craze physics have already been adressed with simulations in previous papers. Baljon and Robbins Baljon1996 ; Baljon2001  demonstrated the importance of chain length for the onset of craze growth. Rottler et al. studied the elastic properties and fracture stresses of fully evolved crazes and used them in combination with linear fracture mechanics to calculate the macroscopic fracture energy of glassy polymers that fail by crazing Rottler2002a . Rottler and Robbins also investigated how polymer entanglements affect the craze structure on a microscopic level and argued that they “jam” the expansion of the glass under tension Rottler2002b .

This paper extends the previous work and is organized as follows. In Section II, we briefly summarize the key experimental observations and review existing theoretical models of crazing. Section III gives the technical details of the molecular models used in this study. We then analyze results for craze nucleation (Section IV), growth (Section V), microstructure (Section VI), and failure (Section VII), and compare our findings to previous models and experiments. Final conclusions are offered in Section VIII.

## II Craze phenomenology and theory

### II.1 Experiments

Crazes have been studied experimentally for more than 30 years Argon1977 ; ArgonHan1977 ; Lauterwasser1979 ; Kramer1983 ; Kramer1990 . The techniques most commonly used to analyse the craze structure are transmission electron microscopy (TEM) Washiyama1992 ; Michler1986 , low angle electron diffraction (LAED) and small angle x-ray scattering (SAXS) Paredes1979 ; Brown1981 ; Brown1987 . Comprehesive reviews of theoretical and experimental results have been presented by Kramer and Berger Kramer1983 ; Kramer1990  and Creton et al. Creton2001 .

The density in the undeformed polymer $`\rho _i`$ and craze $`\rho _f`$ is obtained from TEM measurements. The increase in volume during craze formation, or extension ratio $`\lambda \rho _i/\rho _f`$, is found to have a characteristic value for a given polymer that is independent of molecular weight. Typical values of $`\lambda `$ for different polymers range from two to seven.

Real space images of the craze show that the polymers are bundled into fibrils that merge and split to form an intricate network. The fibrils are highly aligned with the applied tensile stress and vary in diameter and length. However, this complex structure is normally idealized as a set of uniform vertical cylinders connected by short cross-tie fibrils Kramer1990 . The characteristic fibril diameter $`D`$ and separation $`D_0`$ (see Fig. 2) are then determined from a Porod analysis of scattering experiments (see Section V). Measured values range between $`330`$ nm for $`D`$ and $`2050`$ nm for $`D_0`$ Kramer1983 ; Kramer1990 ; Salomons1999a ; Salomons1999b . For example, for polystyrene one obtains $`D6`$nm and $`D_020`$nm Haward1997 .

Nucleation of crazes ArgonHan1977  occurs preferentially near defects in the polymer. These produce large local tensile stresses that lead to the formation of microvoids that evolve into a craze. Once nucleated, the craze grows in length and width (see Fig. 1). It is well established Kramer1990  that the craze widens by drawing material from the dense polymer into new fibrils. This deformation is confined to a narrow active zone at the interface between the dense polymer and craze. The width of the active zone $`h`$ (see also Fig. 2) is usually between $`D`$ and $`D_0`$. Craze widening is a steady-state process, in which a constant “drawing stress” $`S`$ ranging between $`20`$ and $`100`$MPa is applied. Typical experimental values are 35 MPa (polystyrene) Kramer1983  and 70 MPa (polymethymethacrylate) Brown1990 . The value of $`S`$ is of the same order as the shear yield stress of the polymer, and is found to increase with the entanglement density.

### II.2 Theory

A theory of crazing has to explain the molecular origin of the craze structure and the interdependencies of the various quantities measured in experiments. Despite a wealth of experimental data on crazing, there is currently no theoretical description that addresses all aspects of craze physics. The following models have been proposed to explain the extension ratio $`\lambda `$ and the relationship between fibril spacing $`D_0`$ and drawing stress $`S`$.

#### II.2.1 The extension ratio $`\lambda `$

The extension ratio $`\lambda `$ has been successfully explained by a simple scaling argument, that relates $`\lambda `$ to the microscopic entanglement network in the polymer glass. Entanglements arise in dense polymeric systems from the topological constraints that the chains impose upon each other. The mobility of the chains is greatly restricted, because they cannot pass through each other. The starting point for the present argument is the assumption that the glass inherits these entanglements from the melt, where an entanglement molecular weight is given by the plateau modulus under shear, $`G_N^{(0)}`$:

$$M_e=\rho 4RT/5G_N^{(0)}.$$
(1)
This result can be derived from the microscopic tube model Doi1986 , which relates the rheological response of the polymer melt to the deformation of a tube to which the polymer chain is confined. With repeat units of weight $`M_0`$, one can define a typical number of steps $`N_e=M_e/M_0`$ (entanglement length) between entanglements along the polymer backbone.

These entanglements are assumed to act like permanent chemical crosslinks during crazing, which implies that the expansion ends when segments of length $`N_e`$ are fully stretched. The initial separation of entanglement points is $`d_i=(l_pl_0N_e)^{1/2}`$, according to standard random walk (RW) scaling, where $`l_0`$ is an elementary step length and $`l_p`$ the persistence length. The length of this segment rises from $`d_i`$ to a maximum final length $`d_f=\lambda _{\mathrm{max}}d_i=N_el_0`$, and thus

$$\lambda _{\mathrm{max}}=(N_el_0/l_p)^{1/2}.$$
(2)
Experimentally, Eq. (2) is well confirmed, but Section V.4 shows that the picture motivating this expression is oversimplified.

#### II.2.2 The drawing stress $`S`$

The value of the drawing stress $`S`$ has traditionally been related to the craze microstructure ($`D`$, $`D_0`$) via capillary models Kramer1983 ; Kramer1990 . In these models, the polymer in the active zone is treated as a viscous fluid with a surface tension $`\mathrm{\Gamma }`$ and a viscosity $`\eta `$. Figure 2 shows an idealized picture of the craze geometry, where craze formation is modelled as the propagation of void fingers with a characteristic spacing $`D_0`$ into the strain-softened fluid. The applied stress $`S`$ required to advance the interface has a dissipative contribution arising from a suitable flow law (e.g. power-law fluid) and an energy penalty contribution due to the surface tension. The tension is $`S`$ in the polymer glass and the Laplace pressure $`2\mathrm{\Gamma }/(D_0/2)`$ at the ceiling of the finger, where $`D_0/2`$ is the characteristic radius of curvature (see Fig. 2). By estimating the width of the active zone as $`hD_0/2`$, Kramer calculated a stress gradient between glassy polymer and the finger void ceiling,

$$\sigma \frac{\mathrm{\Delta }\sigma }{h}\frac{S4\mathrm{\Gamma }/D_0}{D_0/2}.$$
(3)
Since $`\sigma `$ is proportional to the interface velocity, he then predicted that the system will select a value of

$$D_08\mathrm{\Gamma }/S,$$
(4)
which maximizes the stress gradient between finger ceiling and bulk polymer and thus will lead to the fastest propagation velocity of the fingers.

More recently, Krupenkin and Fredrickson Krupenkin1999  have formulated a theory of craze widening that is similar in spirit to Kramer’s arguments and also equates the craze widening stress with a viscous and a surface tension contribution. However, these authors suggest a different interpretation of $`\mathrm{\Gamma }`$. They introduce an effective surface tension that begins to rise above the vdW value when the finger radius rises above the rms spacing between entanglement lengths $`d_i`$. This ansatz is motivated by the idea that expanding the random walk between entanglements generates an additional energy penalty. An upper bound to $`\mathrm{\Gamma }`$ is provided by the energy required for chain breaking, which sets in once the finger radius exceeds the maximum elongation between entanglement points, $`l_0N_e`$. By minimizing the finger propagation stress, they conclude that the fibril spacing will always be

$$D_0d_i,$$
(5)
independent of surface tension. In their model, the fibril spacing is determined exclusively by the entanglement network.

## III Simulations and Molecular Models

We study craze formation by performing molecular dynamics simulations of a standard coarse-grained polymer model Puetz2000 , where each linear polymer contains $`N`$ spherical beads of mass $`m`$. Models of this kind have a long tradition in polymer research and have verified theories of polymer dynamics Doi1986  in the melt. They have recently been employed by other researchers to study failure in network polymer adhesives Stevens2001  and end-grafted polymer chains between surfaces Sides2001 .

In this bead-spring model, van der Waals interactions between beads separated by a distance $`r`$ are modeled with a truncated Lennard-Jones potential:

$$V_{\mathrm{LJ}}(r)=4u_0\left[(a/r)^{12}(a/r)^6(a/r_c)^{12}+(a/r_c)^6\right]$$
(6)
for $`rr_c`$, where $`u_02040`$ meV and $`a0.81.5`$ nm are characteristic energy and length scales Kremer1990 . A simple analytic potential Sides2001

$$V_{\mathrm{br}}(r)=k_1(rR_0)^3(rR_1)$$
(7)
is used for covalent bonds between adjacent beads along the chain. The form of this potential was chosen to allow for covalent bond breaking, which is not possible with other standard bond potentials such as the popular FENE potential Puetz2000 . Bonds are permanently broken when $`r`$ exceeds $`R_0=1.5a`$. The constant $`R_1=0.7575a`$ was chosen to set the equilibrium bond length $`l_0=0.96a`$, which is the “canonical” value for the bead-spring model with the FENE potential Puetz2000 . This allows us to use results from previous studies, most importantly the entanglement length. The constant $`k_1`$ determines the ratio of the forces at which covalent and van der Waals bonds break. We find that this ratio is the only important parameter in the covalent potential and set it to 100 based on data for real polymers Sides2001 ; Stevens2001 , which implies $`k_1=2351u_0/a^4`$. Tests with other analytical forms of the bond potential showed no appreciable impact on our results as long as the bonds break before the chains can pass through each other.

In order to vary the entanglement length, we include a bond-bending potential Sides2001 ; Faller2000

$$V_B=b\underset{i=2}{\overset{N1}{}}\left(1\frac{(\stackrel{}{r}_{i1}\stackrel{}{r}_i)(\stackrel{}{r}_i\stackrel{}{r}_{i+1})}{|(\stackrel{}{r}_{i1}\stackrel{}{r}_i)||(\stackrel{}{r}_i\stackrel{}{r}_{i+1})|}\right)$$
(8)
that stiffens the chain locally and increases the radius of gyration. Here, $`\stackrel{}{r}_i`$ denotes the position of the $`i`$th bead along the chain, and $`b`$ characterizes the stiffness. We will consider two cases here referred to as flexible $`(b=0)`$ and semiflexible $`(b=1.5u_0)`$ polymers. The corresponding entanglement lengths are $`N_e^{\mathrm{fl}}70`$ and $`N_e^{\mathrm{sfl}}30`$ beads, respectively Kremer1990 ; Puetz2000 ; Faller2000 .

We consider three temperatures $`T=0.01u_0/k_B`$, $`T=0.1u_0/k_B`$ and $`T=0.3u_0/k_B`$, where the last temperature is close to the glass transition temperature. The amount of adhesive interaction between beads is varied by changing the range $`r_c`$ of the LJ potential from $`1.5a`$ to $`2.2a`$.

The equations of motion are solved using the velocity Verlet algorithm with a timestep of $`dt=0.0075\tau _{\mathrm{LJ}}`$, where $`\tau _{\mathrm{LJ}}=\sqrt{ma^2/u_0}`$ is the characteristic time given by the LJ energy and length scales. Periodic boundary conditions are employed in all directions to eliminate edge effects. The temperature is controlled with a Nosé-Hoover thermostat (thermostat rate $`1\tau _{\mathrm{LJ}}^1`$), and the thermostat is only employed perpendicular to the direction of craze growth. Simulations with a Langevin thermostat showed no appreciable difference between the two methods.

In all simulations of crazing, an initial isotropic state in a cubic simulation cell of edge length $`L`$ is created using standard techniques Kremer1990 . Polymer chains are constructed as ideal RWs with a suitably chosen persistence length $`l_p`$. $`l_p`$ is fixed by matching the radius of gyration of the chains to the equilibrium value in the melt, and the values are $`l_p^{\mathrm{fl}}=1.65a`$ and $`l_p^{\mathrm{sfl}}=2.7a`$ for flexible and semiflexible chains, respectively. Subsequently, the interaction potentials are imposed and the system is cooled at constant volume from a melt temperature $`T_m=1.3u_0/k_B`$, to the desired run temperature.

All runs begin at zero hydrostatic pressure. Strains $`ϵ_{ii}`$ are then imposed by rescaling the simulation box periods $`L_i`$ and all particle coordinates proportionately Allen1987 . This allows arbitrary stress states to be studied in Section IV.

## IV Criteria for cavitation and craze nucleation

The loading conditions on the polymer glass determine whether it will fail initially by shear yielding or the formation of voids and cavities. In general, strong triaxial tensile stresses will favor cavitation. Cavitation and crazing are closely related, because crazes usually require the initial formation of microvoids Kramer1990 . We therefore first address the initial failure of the polymer glass through either shear yielding or cavitation, and later discuss the formation of crazes.

The loading conditions that lead to shear yielding in many experimental polymers Kody1997 ; Raghava1973  are most accurately represented by the pressure-modified von Mises yield criterion. It is formulated in terms of simple stress invariants, the hydrostatic pressure $`p=(\sigma _1+\sigma _2+\sigma _3)/3`$ and the deviatoric or octahedral shear stress $`\tau _{\mathrm{oct}}=\left((\sigma _1\sigma _2)^2+(\sigma _2\sigma _3)^2+(\sigma _3\sigma _1)^2\right)^{1/2}/3`$, where the $`\sigma _i`$ denote the three principal stress components. The pressure-modified von Mises criterion states that yield will occur at an octahedral yield stress $`\tau _{oct}^y`$ given by

$$\tau _{oct}^y=\tau _0+\alpha p,$$
(9)
where $`\tau _0`$ is the yield stress at zero hydrostatic pressure and $`\alpha `$ is a dimensionless constant. Its physical motivation is that the elastic free energy stored in shear deformation is proportional to $`\tau _{\mathrm{oct}}^2`$ and failure should occur when this energy exceeds a threshold that rises slowly with $`p`$.

In ref. Rottler2001, , we examined a much larger range of stress states than in previous experimental studies and showed that the pressure-modified von Mises criterion provides a good description of shear yield in our bead-spring model. Data for two extremal temperatures are replotted in Fig. 3 along with solid lines showing fits to Eq. (9). Shear yield was observed to the right of the dot-dashed line, and these data points follow Eq. (9) quite accurately. To the left of the line cavitation was observed. The deviation from the von Mises fits is very sharp, and $`\tau _{\mathrm{oct}}`$ quickly drops to zero. The values of $`\tau _{\mathrm{oct}}^c`$ where cavitation occurs are well described by a straight lines

$$\tau _{oct}^c=\tau _0^c+\alpha ^cp$$
(10)
with new constants $`\tau _0^c`$ and $`\alpha ^c`$. This new “cavitation criterion” can be motivated in analogy to the von Mises criterion by assuming that the elastic free energy $`F_V`$ associated with volume changes must reach a critical value for cavitation to occur. $`F_V`$ is proportional to $`p^2`$, which gives a criterion of the form $`p=p_0`$. One can then assume that shear components in the stress tensor aid cavitation in a linear fashion, i.e. $`p=p_0+\tau _{\mathrm{oct}}/\alpha _c`$, which can be rearranged to give Eq. (10) with $`\tau _0^c=\alpha _cp_0`$.

No clear experimental consensus exists about the stress state required for crazing, partly because of the importance of surface defects in nucleating crazes. However, several criteria for craze nucleation were proposed almost 30 years ago. They all try to take into account the critical role of tensile stress components. Sternstein et al. Courtney1990  suggested a craze yield criterion of the form

$$\tau _{\mathrm{max}}\frac{1}{2}|\sigma _i\sigma _j|_{\mathrm{max}}=A+B/p,$$
(11)
where $`A`$ and $`B`$ are constants that depend on temperature. With respect to our criterion Eq. (10), $`p`$ has been replaced by $`1/p`$ and $`\tau _{\mathrm{oct}}`$ by the largest difference between any two stress components. Bowden and Oxborough Haward1997  formulated a similar criterion, where $`\tau _{\mathrm{max}}`$ is replaced by $`\sigma _1\nu \sigma _2\nu \sigma _3`$ and $`\nu `$ is Poisson’s ratio for the polymer glass. This expression is another possibility to describe the shear components of the stress state, and it reduces to $`\tau _{\mathrm{max}}`$ when $`\nu =1/2`$ and $`\sigma _2=\sigma _3`$. The Sternstein and Bowden and Oxborough expressions could in principle also be fitted to the rather narrow range of pressure in Fig. 3 where cavitation occurs. However, we are unaware of a convincing physical motivation for the $`1/p`$ term, which leads to obvious analytical problems at small $`p`$. In addition, the experimental results that motivated Eq. (11) are sensitive to surface defects ArgonHan1977 .

The above considerations pertain to the initial mode of failure of the polymer glass at strains typically less than 10%. However, crazing is a large strain deformation with strains of several hundred percent. Although we find voiding to be a neccessary precursor to crazing, it is not guaranteed that a loading state that leads to cavitational failure according to Eq. (10) will ultimately produce stable crazes. Likewise, we have observed that an initial failure through shear deformation can still lead to later void formation and crazing. One should thus strictly call Eq. (10) a cavitation failure criterion and not a craze yielding criterion.

## V Growth of Crazes

In order to induce crazing, we enforce cavitation by expanding the periodic simulation box in the z-direction at constant velocity while maintaining the simulation box periods in the perpendicular $`xy`$ plane. This leads to an initial stress state where all three principal stresses are tensile. The initial voids formed during cavitation expand upon further straining, but their growth rapidly becomes arrested Baljon2001 . Instead of forming new voids, additional material is drawn out of the uncavitated polymer, and stable craze growth occurs. In our simulations, growth continues until all material in the simulation box is converted into the craze.

### V.1 Images of crazes

A good impression of the crazing process can be obtained by inspecting the snapshots of the simulation cell shown in Figs. 4 \- 6. Each slice has a lateral width of $`64a`$, and three different strains are shown. In all images, the chain length $`N=512`$. Previous studies Baljon2001  had shown that $`N`$ has to be twice the entanglement length or greater in order to form stable crazes. For shorter chains, the material cavitates, but then rapidly fails due to chain pullout. In the following, we only consider chains with $`N2N_e`$.

Note first that in all cases, there is a sharp interface between dense polymer and crazed material. This narrow ’active zone’ is one of the key features of craze phenomenology found in experiment. In the craze, the polymer chains have merged into fibrils that are strongly aligned. However, the structure is quite complicated, as there are many lateral connections between fibers.

One can also observe that the fine structure of the crazes in the three sequences varies greatly. Fig. 4 with flexible chains at the low temperature of $`T=0.1u_0/k_B`$ and the weak adhesive interaction (cutoff distance $`r_c=1.5a`$) shows many thin fibrils, whereas the fibrils in Fig. 6 at the higher temperature of $`T=0.3u_0/k_B`$ and the stronger adhesive interaction $`r_c=2.2a`$ are much thicker in diameter. These trends are not surprising, because increased chain mobility at higher temperatures and stronger adhesive interactions should drive the system to larger fibril diameters, which minimize the surface area.

### V.2 The drawing process and stress-strain curves

A second characteristic feature of craze growth is that deformation occurs at a constant plateau or drawing stress $`S`$. This plateau can be easily identified in the stress-strain curves shown in Fig. 7. The curves can be separated into three different regimes. In regime I, the stress rises to a peak of $`2.6u_0/a^3`$ and then drops when the polymer yields by cavitation. Following cavitation, the stress rapidly relaxes and remains at the plateau value $`S`$ in regime II, the growth regime. Regime II is much shorter in the semiflexible case Fig. 7(b) than in the flexible case Fig. 7(a) (note different lateral scales). Regime II ends when the strain $`L_z/L`$ reaches the extension ratio $`\lambda `$. At this point, all the material in the simulation cell has been converted into the craze, and additional deformation strains the entire craze uniformly. As a consequence, the stress rises again in regime III. This regime finally ends in catastrophic failure either through chain disentanglement or chain scission (see Section VII).

Note first that neither the peak stress at cavitation nor the value of $`S`$ depends on the chain length $`N`$. The curves for different $`N`$ in Fig. 7 only split apart after completion of craze growth when $`L_z/L`$ reaches $`\lambda `$ and the entire craze is strained. Baljon and Robbins Baljon2001  showed that the peak stress remained constant for much shorter chains, but that regime II only appeared when $`N`$ was $`2N_e`$ or longer. Another important fact to note is that $`S`$ is independent of system size. For example, values of $`S`$ in systems ranging between 32768 and 1048576 beads are the same within a few percent. The biggest change with increasing system size is that temporal fluctuations in $`S`$ decrease.

In Fig. 8(a), we analyze trends of $`S`$ with $`T`$ and $`r_c`$. The drawing stress decreases linearly with increasing temperature and increases with increasing adhesive interactions (i.e. increasing $`r_c`$). Fig. 8 (b) shows that $`S`$ varies logarithmically with the widening velocity $`v`$ over two orders of magnitude, which is indicative of a thermally activated process. For the subsequent figures, we choose $`v=0.06a/\tau _{\mathrm{LJ}}`$, which is at the upper end of the logarithmic regime Baljon2001 . Similar behavior is also found for the shear yield stress of glassy polymers Rottler2001 ; Rottler2003 .

### V.3 Crazing under plane stress conditions

The results of Section IV show that cavitation only occurs when all three principal stresses are tensile. Many experimental crazes grow in a thin film geometry under plane stress conditions. However, in these experiments the craze is often prenucleated or nucleates near a defect ArgonHan1977 . This situation can also be mimicked in our simulations. To this end, the periodic boundary conditions in the $`x`$ direction were replaced with free boundaries, so that the solid is free to relax in that direction. Initial failure is now nucleated by placing 1000 purely repulsive LJ beads in the center plane of the simulation cell located at $`z=L_z/2`$ <sup>1</sup><sup>1</sup>1This is only done in Fig. 9. All other simulations use fully 3D periodic boundary conditions without repulsive beads.. This weakens the solid locally and constrains the location of initial failure, while not affecting subsequent craze growth.

Fig. 9 shows three snapshots of a craze in this geometry. As in experiments, necking is observed at the craze-bulk interface. Although $`\sigma _{\mathrm{xx}}`$ vanishes in the rest of the film, the neck produces strong tensile stresses in all three directions in the active zone. The craze grows in the same fashion as in the simulations with 3D periodic boundary conditions. Since the latter yield better statistics for the craze structure, we have focused on this methodology for our analysis.

### V.4 The extension ratio

The extension ratio $`\lambda `$ can be calculated from the average densities of crazed and uncrazed material. Fig. 10 shows how the density drops from the initial value $`\rho _i`$ to $`\rho _f`$ in the craze. As can be seen, $`\rho _f`$ is higher for the semiflexible chains, which have a smaller value of $`N_e30`$. Remarkably, we find that $`\lambda `$ is a function of $`N_e`$ only and decreases with decreasing $`N_e`$. For instance, while increasing $`T`$ and $`r_c`$ produces dramatic coarsening of the fibril structure in Fig. 6 relative to Fig. 4, $`\lambda `$ is unchanged. We obtain values of $`\lambda _{\mathrm{fl}}=6.0\pm 0.6`$ and $`\lambda _{\mathrm{sfl}}=3.5\pm 0.3`$ independent of $`N`$, $`T`$, and adhesive interaction strength.

In order to understand the dependence of the macroscopic quantity $`\lambda `$ on $`N_e`$, we analyze the structural changes in the polymer glass during deformation on a microscopic level (see also ref. Rottler2002b ). Figure 11(a) shows the average final position of beads in the completely evolved craze as a function of their initial positions along the direction of the expansion (z-axis). The average was taken over all beads with initial heights in a bin of width $`1`$a. Although the strain rate is strongly localized during the craze process, the ultimate displacement profile is linear, $`z_f=\lambda z_i`$.

To measure deviations from a purely affine (uniform) deformation, we evaluated the rms variation $`\delta z`$ in $`z_f`$ for beads in each bin. This quantity is indicated by error bars in Fig. 11(a). Note that the variation in each bin is very reproducible. We find that $`\delta z`$ is nearly independent of $`T`$ and $`r_c`$ and has values on the order of 19a and 9a for flexible and semiflexible chains, respectively.

Since no strain is applied in the perpendicular x and y directions, one would assume that there is on average no displacement in these directions. That this is indeed the case is shown in Fig. 12, which repeats the analysis of Fig. 11 for the x-direction. Average final bead positions are identical to initial positions, but there are lateral variations $`\delta x`$ that are indicated by error bars. These lateral displacements allow chains to gather in fibrils at the initial density to minimize surface area. Unlike the vertical displacements $`\delta z`$, these lateral displacements depend strongly on $`T`$ and $`r_c`$. For example, $`\delta x2.5`$a for the fine structure shown in Fig. 4, where many thin fibrils can be seen, while $`\delta x5.6`$a for the much coarser structure of Fig. 6. In general, $`\delta x`$ correlates with the spacing between fibrils as discussed in Section VI and is less than $`d_i`$. Krupenkin and Fredrickson Krupenkin1999  suggested that $`d_i`$ provides an upper bound for the lateral chain deformations.

We now examine changes in the conformation of individual chains. In the initial state, the polymer chains exhibit an ideal random walk (RW) structure inherited from the melt. The average end-to-end vector $`𝐑^\mathrm{𝟐}`$ thus scales with the number of covalent bonds connecting two beads $`\mathrm{\Delta }N`$ as $`𝐑^2=l_pl_0\mathrm{\Delta }N`$. The component along each direction is $`1/3`$ of that value since the initial state is isotropic. Fig. 12(b) shows this initial scaling behavior for $`\mathrm{\Delta }x^2`$ (dashed line) and that $`\mathrm{\Delta }x^2`$ is not affected by crazing (solid line).

After an affine deformation by $`\lambda `$ along $`z`$, one would have an anisotropic RW with no change in $`\mathrm{\Delta }x`$ or $`\mathrm{\Delta }y`$, but $`\mathrm{\Delta }z^2=\lambda ^2l_pl_0\mathrm{\Delta }N/3`$. Fig. 11(b) shows the actual behavior (solid lines) of $`\mathrm{\Delta }z^2`$ in the craze. At large scales, it exhibits the expected scaling for an affine deformation (dashed lines). However, the separation between beads is fixed by the length of the covalent bonds, so the deformation of individual polymers along $`z`$ cannot be purely affine. At small scales, the linear scaling behavior of $`\mathrm{\Delta }z^2`$ crosses over into a quadratic behavior, which indicates that the polymer has been pulled taut on this scale. The typical number of beads in such a straight segment $`\stackrel{~}{N}_{\mathrm{st}}`$ can be calculated by letting $`(\stackrel{~}{N}_{\mathrm{st}}l_0)^2=\mathrm{\Delta }z^2=\lambda ^2l_pl_0N_{\mathrm{st}}/3`$ at the crossover point, which yields $`\stackrel{~}{N}_{\mathrm{st}}=\lambda ^2l_p/3l_0`$. Inserting the observed values of $`\lambda `$, $`l_p`$, and $`l_0`$, we arrive at $`\stackrel{~}{N}_{\mathrm{st}}^{\mathrm{fl}}=21\pm 4`$ and $`\stackrel{~}{N}_{\mathrm{st}}^{\mathrm{sfl}}=12\pm 2`$, respectively. These values are comparable to the values of $`\delta z`$ found in Fig. 11(a). On this scale, the deformation is non-affine.

The length of taut sections can also be determined by direct analysis of the chain geometry. To this end, we calculate the angle between every covalent bond and the z-axis and label a bond as pointing up (down) if the angle is within $`45^{}`$ of the z (-z) axis. We then count the number $`N_{\mathrm{st}}`$ of consecutive up (down) steps. The probability $`P(N_{\mathrm{st}})`$ of finding a straight segment containing $`N_{st}`$ steps is shown in Fig. 13(a). For both flexible and semiflexible chains, the distribution develops an exponential tail. Like $`\lambda `$, this tail is independent of $`N`$, $`T`$, and $`r_c`$. The characteristic length scales that arise from these tails are $`\stackrel{~}{N}_{\mathrm{st}}^{\mathrm{fl}}21`$ and $`\stackrel{~}{N}_{\mathrm{st}}^{\mathrm{sfl}}13`$, in good agreement with the prediction from the RW argument. Fig. 13(b) shows that very similar length scales arise from an equivalent analysis of the decay of the correlation function for the z-component of successive bonds.

In section II.2.1, we introduced the standard scaling argument Eq. (2) that relates extension ratio and entanglement length, which has been verified experimentally with great success. In our cases, it predicts $`\lambda _{\mathrm{max}}^{\mathrm{fl}}=6.5`$ and $`\lambda _{\mathrm{max}}^{\mathrm{sfl}}=3.5`$, which agree with the observed values of $`\lambda `$. However, the argument was motivated by the idea that segments between entanglements become fully stretched and thus appears to be at odds with the finding of an average straight segment length of only $`N_e/3`$ rather than $`N_e`$. This discrepancy is resolved by realizing that since the deformation is uniaxial, only the projection of $`d_i`$ onto the z-axis, $`d_i\mathrm{cos}(\mathrm{\Theta }_z)`$, is expanded, where $`\mathrm{\Theta }_z`$ is the angle between $`d_i`$ and the z-axis. The average projection is thus only $`1/\sqrt{3}`$ of the total length. Indeed, it was noted already in earlier work Kramer1983  that, due to this geometric factor, $`\lambda `$ should be $`\sqrt{3}\lambda _{\mathrm{max}}`$ for fully stretched chains. However, this result is little cited since $`\lambda \lambda _{\mathrm{max}}`$ in many systems and, until our work, there was no reason to expect the length of straight segments to be $`N_e/3`$.

The emergence of the length scale $`N_e/3`$ is a consequence of the random nature of the entanglement mesh. Clearly, all strands would be expanded simultaneously by the same factor in a regular mesh as reported in a simulation study by Stevens Stevens2001 . By contrast, in the polymer glass only the segments that are initially aligned with the stretching direction become fully stretched. These fully stretched segments are able to prevent further extension, because the entanglements act like chemical crosslinks. Barsky and Robbins have confirmed the equivalence between entanglements and crosslinks by adding permanent crosslinks randomly to the system Barsky2000 . The length between constraints then decreases from $`N_e`$, and $`\lambda _{\mathrm{max}}`$ decreases accordingly. They found $`\lambda \lambda _{\mathrm{max}}`$ in all cases and that the average stretched length $`N_{\mathrm{st}}`$ remains at $`1/3`$ of the distance between constraints.

The success of the scaling argument Eq. (2) and the constancy of the extension ratio imply that there is no appreciable loss of entanglements in our simulations during craze growth. Chains do not disentangle once $`N>2N_e`$, and chain scission (see also section VII) is not observed during growth for any choice of parameters in our model.

### V.5 Energy dissipation and stress transfer during crazing

The work done in transforming a volume $`dV`$ of polymer into a craze is $`\delta W=S(\lambda 1)dV`$. This work can either increase the potential energy $`dU`$ or be dissipated as heat $`\delta Q`$. The division between energy and heat is hard to determine experimentally, but simulations with short chains found that both contributions were substantial Baljon1996 . We have measured $`\delta W`$ and the energy change directly in our simulations and calculated $`\delta Q`$ using the first law of thermodynamics: $`\delta Q=dU\delta W`$. $`dU`$ can be calculated directly from the bead positions and interaction potentials. Table 1 shows the fraction $`\delta Q/\delta W`$ of dissipated total work for a number of large systems. In all cases, a large percentage, $`80\%`$, of the total work is dissipated, and only $`20\%`$ is stored as potential energy.

Since the craze drawing stress varies logarithmically with velocity (see Fig. 8(b)), these percentages could change with velocity. However, we find that $`dU`$ also decreases with decreasing velocity, and there is no measurable change in the percentage of work converted to heat over at least two orders of magnitude in velocity.

Stress in the craze can also be partitioned into two components that originate either from van der Waals (LJ) interactions Eq. (6) or from covalent interactions Eq. (7). The two contributions are very different in the uncrazed and crazed material. In the undeformed polymer, the tensile stress is mainly carried by the vdW bonds. This stress is transferred in large part to the covalent bonds as they pass through the active zone. Evidence for this is provided in Fig. 14, which shows the covalent and LJ contributions to the total stress as a function of height in the widening direction. Panel (a) displays the density profile in order to identify dense polymer regions (high density) and craze regions (low density). In Panel (b) one observes that in the dense region all the tension is carried by the van der Waals bonds and the covalent bonds are under slight compression. In the craze, between $`6095\%`$ of the total stress (see Table 1) is carried by the covalent bonds, and the van der Waals bonds only contribute a small fraction.

### V.6 Problems with the capillary models

The results presented so far reveal serious difficulties with the surface tension models discussed in the Introduction. The first evidence of this comes from the observation that $`S`$ is independent of system size. In our smallest simulations of lateral width $`32a`$ the simulation box only contains a few fibrils at $`T=0.1u_0/k_B,r_c=2.2a`$. If $`S`$ were controlled by Eq. (4), one should expect that the simulation box would need to contain a statistically significant number of fibrils of spacing $`D_0`$ for $`S`$ to reach its steady state value. However, the value of $`S`$ does not fluctuate as the lateral dimensions are increased to $`64a`$ or $`128a`$.

The second and more severe problem concerns the distribution of stress in the craze. The surface tension model assumes that all the stress is carried at the interfaces of fibrils or in viscous stress in the active zone. However, Table 1 and Fig. 14 show that almost all the stress in the fibrils is carried by the covalent bonds, while the surface tension is entirely associated with broken vdW bonds and small entropic contributions.

In the following, we make an alternative proposition to relate craze microstructure and drawing stress. This proposition is based on the observation that the values of $`S`$ and $`\lambda `$ obtained from our simulations obey the equality

$$S_{\mathrm{fl}}\lambda _{\mathrm{fl}}=S_{\mathrm{sfl}}\lambda _{\mathrm{sfl}}=S_0(T,r_c,v).$$
(12)
This can be verified for each $`T`$ and $`r_c`$ using $`\lambda _{\mathrm{fl}}=6.0`$, $`\lambda _{\mathrm{sfl}}=3.6`$ and values of S from Fig. 8. Since the fraction of area occupied by the fibrils is $`1/\lambda `$, $`S_0`$ is the local stress within the fibrils.

It is perhaps surprising that the value of local stress needed to draw fibrils is independent of the entanglement length. In order to further test Eq. (12), it would be desirable to consider additional values of $`N_e`$ and thus $`\lambda `$. Unfortunately, reliable values for $`N_e`$ exist for only a few values of $`b`$. It is, however, not necessary to know $`N_e`$ for the present purpose, since both $`S`$ and $`\lambda `$ can be measured directly from the craze simulation. Moreover the entanglement length should only depend upon the chain statistics Fetters1994 , and glassy states with arbitrary statistics can easily be created. We confirmed that simulations with the same persistence length in the undeformed glass gave the same values of $`\lambda `$ and $`S`$ independent of whether the bond-bending potential (Eq. (8)) was included. Fig. 15(a) compares stress strain curves for $`b=0u_0`$ at four values of the initial persistence length. Increasing $`l_p`$ lowers $`\lambda `$ and increases $`S`$. However, Fig. 15(b) shows that all curves can be collapsed if $`\sigma _{zz}`$ is scaled by $`\lambda `$ and the extension by $`\lambda 1`$. This confirms that $`S_0`$ is the stress that controls craze growth and only depends on the van der Waals interactions. An experimental version of this test would be difficult, since it is hard to change $`N_e`$ for real polymers without changing the chemistry as well.

It is interesting to compare values of $`S_0`$ to the stresses required for shear yielding and cavitation, which are also independent of chain statistics. Table 2 compares these three stresses for two ranges of the LJ potential. Here $`b=0u_0`$, but the bond-bending potential has little effect on the values. The three stresses are clearly correlated, decreasing with increasing temperature and decreasing $`r_c`$. The local fibril stress is always about twice the cavitation stress and ranges from 7 to 12 times $`\tau _0`$. The implication is that local stress for drawing material into fibrils $`S_0`$ is related to the bulk yield stresses, but it is difficult to determine the relative role of shear and cavity growth. It is interesting to note that the experimental values of $`S_0/\tau _0`$ for PS and PMMA are about 5 Kramer1983 ; Quinson1997 , and it would be useful to have values of this ratio for other polymers. A comparison to $`p_{\mathrm{cav}}`$ would also be interesting, but its value is sensitive to system size, strain rate and inhomogeneities and it is difficult to measure experimentally.

Given that previous results for $`N_e`$ Kremer1990 ; Puetz2000 ; Faller2000  are consistent with values inferred from the extension ratio (Eq. 2), our results for $`\lambda `$ as a function of $`l_p`$ allow a rapid estimation of $`N_e`$. Table 3 presents results for a wide range of $`l_p`$ and shows that the product $`N_el_p`$ is constant within our errorbars. Fetters et al. Fetters1994  have presented a model for the relation between chain statistics and $`N_e`$ that predicts

$$N_e\frac{N^3}{R^2^3}$$
(13)
where $`R^2=l_pl_0N`$ denotes the average end-to-end vector of the polymer chain and the proportionality constant only depends upon density. This implies $`N_el_p^3`$, while our data is clearly consistent with a simple inverse relation $`N_el_p^1`$. Equation (13) describes many experimental polymers, but it is difficult to change $`l_p`$ without changing all the other parameters in the equation. The flexible model ($`b=0u_0`$) is known to be quantitatively inconsistent with Equation (13) Puetz2000 , which has been one motivation for studies of more rigid models. It would be interesting to have additional values of $`l_p`$ from melt simulations to test whether the inverse relation between $`N_e`$ and $`l_p`$ found here holds more generally, and, if so, to understand its origin.

### V.7 Width of the active zone

At the interface between dense polymer and craze, polymer chains are locally mobilized and brought into the new fibril structure. The region in which this motion takes place is called the ’active zone’. In Fig. 2 the height of the active zone $`h`$ was defined as the distance between undeformed polymer layer and the void ceiling, and this layer was assumed to behave like a strain-softened fluid. The main drop in density should occur over the height of the void ceiling of order $`D_0/2`$ Kramer1990 . In this section, we compare this simple picture to our simulations.

Fig. 10 shows typical results for the density profile near the craze boundary. For both flexible and semiflexible polymers, the density drops over a region of width $`20a`$. The average strain rate $`\dot{ϵ}`$ must be localized in the same region, since $`\dot{ϵ}=\mathrm{ln}\rho /t`$. If the active zone advances at velocity $`v`$ then $`\dot{ϵ}=\pm v\mathrm{ln}\rho /z`$ (the sign depends on whether the top or bottom interface is growing). In Fig. 16(a) we present $`\dot{ϵ}`$ as a function of position along $`z`$. Averages were taken over layers of height 1a. The curves shown correspond to the 4 different values of $`l_p`$ used in the creation of the polymer chains. As discussed above, this varies the entanglement length $`N_e`$. Curves were shifted by $`z_0`$, which corresponds to the center of the peak. At this point, the density is roughly halfway between $`\rho _i`$ and $`\rho _f`$. Curves for very different $`N_e`$ essentially overlap with a width of $`10a`$ at half maximum and $`20a`$ at the base. The latter corresponds to the range of rapid density change in Fig. 10.

The width of the region over which beads are mobilized was determined from the relative diffusion as a function of height. Fig. 16(b) shows the standard deviation of displacements in the lengthening direction $`\sigma (\mathrm{\Delta }z)`$. The curves peak at the same location as the curves in Fig. 16(a), but are more asymmetric. In the direction of the dense polymer glass $`\sigma (\mathrm{\Delta }z)`$ and $`\dot{epsilon}`$ fall to zero over a comparable range. By contrast, $`\sigma (\mathrm{\Delta }z)`$ shows a long exponential tail into the craze with a characteristic decay length $`l_{\mathrm{AZ}}`$ that varies with entanglement length. The fit values of $`l_{\mathrm{AZ}}`$ indicate that there is a definite trend to larger values as $`N_e`$ increases, and $`l_{\mathrm{AZ}}`$ tends to be somewhat smaller than $`N_e`$. This result is not surprising, because $`N_e`$ is the longest length scale over which particle mobilization should occur. Standard deviations of the lateral displacements $`\mathrm{\Delta }x`$ and $`\mathrm{\Delta }y`$ are smaller, but show essentially the same decay lengths.

The above analysis indicates that while the mobility of the beads is constrained by entanglements, the regions of localization of strain rate and the density drop are related to the craze microstructure (see Section VI) . Typical values for $`D`$ and $`D_0`$ are given in Table 4. From this, $`D_0/210a`$, which compares well to the width of the strain localization peak at half maximum. Experimentally, the width of the active zone has been measured by a gold decoration technique Kramer1990 . It was concluded that it lies between $`D`$ and $`D_0`$, which agrees with our results.

## VI Microstructure of Crazes

Another fascinating aspect of crazes is their complicated microstructure. Figs. 4 \- 6 give an impression of the range of length scales appearing in the voided fibril network. Clearly, the picture of cylindrical fibers and void fingers (Fig. 2) is an oversimplification. It is nevertheless helpful to build more realistic models starting from this simple scenario.

### VI.1 Structure factor

Experimentally, the standard measurement of the craze microstructure is done via scattering experiments. The scattering intensity in these measurements is proportional to the structure factor

$$S(𝐤)=\frac{1}{N}\underset{j,k}{}e^{i𝐤(𝐫_j𝐫_k)}$$
(14)
times the form factor for the monomers. Since the craze structure has azimuthal symmetry, one decomposes the wavevector $`𝐤`$ into components parallel and perpendicular to the fibrils. Contour plots of $`S(k_{},k_{})`$ for two crazes are shown in Fig 17. The microstructure was varied by changing the cutoff length $`r_c`$ and chain flexibility. Both patterns are asymmetric with the intensity decaying much faster in the direction parallel to the fibrils than perpendicular to them.

Most experimental setups integrate over $`k_{}`$ using slit collimation and measure the integrated structure factor $`S(k_{})=_{\mathrm{}}𝑑k_{}S(𝐤)`$. In Fig. 18(a), we plot $`S(k_{})`$ as a function of the magnitude of $`k_{}`$ for the same systems shown in Fig. 17. At large wavevectors, the curves rise to a peak at $`2\pi /a`$ (not shown), which corresponds to the typical separation of two beads. This length scale is so short that it is usually not resolved in typical experimental scattering plots shown in e. g. Brown1987 . The characteristic feature of $`S(k_{})`$ is found at smaller $`k_{}`$ in form of a power-law regime with exponent $`3`$. The extent of this scaling regime is bound at large wavevectors by the small scale cutoff provided by the interparticle spacing and at small wavevectors by the distance to the next fibril. The power-law regime is more pronounced for the crazes shown in Fig. 18(b), which were created at higher temperatures where $`D_0`$ is bigger.

### VI.2 Interpretation of scattering data

The traditional interpretation of craze scattering data begins with idealizing a craze fibril as a straight cylinder of diameter $`D`$ and length $`l`$ along $`z`$. The scattering intensity is then proportional to the squared magnitude of the form factor for such a cylinder Brown1987 ,

$$F(k_{})=\frac{\mathrm{\Delta }\rho _{el}\pi D^2l^{1/2}}{2}\frac{J_1(\pi Dk_{})}{\pi Dk_{}},$$
(15)
where $`\mathrm{\Delta }\rho _{el}`$ denotes the electron density and $`J_1`$ is the first-order Bessel function. Due to the asymptotic behavior of $`J_1(x)=(2/\pi x)^{1/2}\mathrm{cos}[x3\pi /4]+𝒪(x^1)`$ for large arguments, the scattering intensity of a single cylinder will exhibit an oscillating power-law behavior, $`|F(k_{})|^2k_{}^3`$. This is also called Porod scattering.

In general the fibrils do not have a single diameter, but rather a diameter distribution $`P(D)`$. One can introduce an average scattering intensity

$$I_0(k_{})=F(k_{})^2=_{D_{\mathrm{min}}}^{D_{\mathrm{max}}}𝑑DP(D)F(k_{})^2.$$
(16)
by averaging the form factor over the diameter distribution and neglecting correlations between fibrils. The main effect is to smooth out the oscillations so that a straight power-law tail results. The average diameter $`D=DP(D)𝑑D`$ can be obtained via a Porod analysis, in which one determines the prefactor $`\alpha `$ to the power-law tail, $`S(k_{})=\alpha k_{}^3`$. This can be related to $`D`$ Brown1981 ; Paredes1979  through,

$$D=\frac{Q}{\pi ^3(11/\lambda )\alpha },$$
(17)
where $`Q=𝑑k_{}2\pi k_{}S(k_{})`$ is a scattering invariant. Values for $`D`$ obtained from this formula are collected in Table 4.

The craze fibrils do not all have the same distance from each other, but have in general varying distances that can be described by a radial distribution function $`g(r)`$. This will lead to interference effects in the scattering intensity that can be described by an interference function $`j(k_{})=I(k_{})/I_0(k_{})1`$. The interference function is related to $`g(r)`$ by Brown1987

$$g(r)=1+\frac{\pi D^2\lambda }{4}_0^{\mathrm{}}2\pi k_{}j(k_{})J_0(2\pi rk_{})𝑑k_{}.$$
(18)
As a result, the power-law tail will be modified at small wavevectors. In particular, the first peak in $`g(r)`$, corresponding to a typical fibril separation $`D_0`$, should translate into a maximum in $`I(k_{})`$.

The Porod scattering law $`|F(k_{})|^2k_{}^3`$ is well confirmed at higher temperatures in Fig. 18(b), while the power-law regime is shorter at lower temperatures. Visual inspection of the craze images suggests that the straight cylinder approximation is not so well justified in this case. Fibrils branch more often and intersect the z-axis at varying angles. At the higher temperature, the chains are more mobile and can align more easily, but they are still not ideal cylinders. We note furthermore that the curves shown in the log-log plot of Fig. 18 do not exhibit a clear maximum (a maximum would be more easily identifiable in a linear plot normally used for experimental results). This would suggest that the ordering of fibrils is mostly random without a clear characteristic separation. However, our statistics are limited by the system size at these large length scales.

### VI.3 Distributions of fibril diameter and spacing from real space analysis

In a previous analysis of scattering data Brown1981  for polystyrene and polycarbonate crazes, Eqs. (15)-(18) were used to extract the diameter distribution $`P(D)`$ and $`g(r)`$ by means of a detailed fitting procedure. The craze images shown above suggest that $`P(D)`$ is rather broad, and wide distributions for $`P(D)`$ were found for both materials, with a significant increase in breadth for polycarbonate Brown1987 .

Here, we access these distributions by direct geometrical analysis of the bead positions. To this end, we bin the particle positions onto a square grid with gridsize $`1a`$ normal to the widening direction and take lateral cross-sections of height $`1.5a`$. As illustrated in Fig. 19, a fibril now appears as a 2-dimensional connected cluster, whose area $`A`$ is taken to be the sum of the areas of the occupied squares. We define $`D=\sqrt{4A/\pi }`$.

Fig. 20 shows the resulting distributions of $`D`$. As expected, $`P(D)`$ is very broad. The distributions for flexible and semiflexible chains are very similar at small $`D`$ for a fixed temperature, but differ for larger $`D`$. The tail of the distributions could be fitted to an exponential function, but our statistics are too small for a conclusive statement.

Mean values of the diameter are given in Table 4 together with the standard deviations of the distributions. The large values of the latter suggest that $`D`$ has to be used with care when describing the craze microstructure. Previously, Baljon and Robbins reported similar values of $`D=7a`$, $`\sigma (D)=11a`$ for flexible chains at $`T=0.3u_0/k_B,r_c=1.5a`$ Baljon2001 .

The value of $`D`$ obtained from the Porod analysis is always larger than the value from the cluster (real-space) analysis. Both values rise with increasing adhesive interaction and increasing temperature as expected. Taking $`a0.8`$nm, diameters in our model crazes would correspond to a range of $`D=1.46.3`$nm, which is at the small end of the experimental range. The reason is that an artificially small value of the surface tension and a high widening velocity $`v`$ are used here. Both lower $`D`$, which allows us to use smaller system sizes.

An estimate for the mean fibril spacing $`D_0`$ can be obtained by equating the area per fibril, $`\pi D_0^2/4`$, to the inverse areal density $`1/n`$, i.e. $`D_0=2\sqrt{1/n\pi }`$. The areal density was obtained by counting the number of separate fibrils per cross-section. Values for $`D_0`$ are also given in Table 4 and translate into a range between 10 nm and 25 nm. The higher numbers are comparable to experiment and are obtained with $`r_c=2.2a`$, which produces more realistic surface energies.

In order to obtain the radial distribution function of the fibrils, we continued the analysis described above and calculated the center of mass for each 2D cluster. The positions given by this procedure were used to calculate $`g(r)`$ in Fig. 21. In general, these functions have very little structure. There is a size exclusion minimum at the origin, and the curves have a weak first maximum around $`10a`$. As the fibrils become thicker, the location of the maximum shifts to larger values. Qualitatively similar curves were obtained from experiment Brown1987 , which confirms the basically random nature of fibril positions. The height of the maximum is too small to be reflected in the scattering intensity.

## VII Stress Distribution and Craze Breakdown

In regime III of the stress-strain curve of Fig. 7 the entire volume of the simulation cell has been converted to a craze. Elongation past the extension ratio causes uniform straining of the craze and eventually leads to craze failure. Studies of this regime are directly relevant to crack propagation in glassy polymers (Fig. 1). The stress in the craze region rises from $`S`$ at the active zone to a maximum value $`S_{\mathrm{max}}`$ at the crack tip. The elastic properties of the craze determine the rate at which the stress rises with distance, and $`S_{\mathrm{max}}`$ determines how large the craze region can become before the crack propagates Brown1991 . These properties were recently obtained from MD simulations and combined with continuum theory to predict the macroscopic fracture energy Rottler2002a . Here we focus on the microscopic stress distribution and its relation to $`S_{\mathrm{max}}`$.

### VII.1 Disentanglement versus chain scission

The craze can fail by two different mechanisms that depend on the chain length $`N`$: short chains can disentangle, while very long chains fail through chain scission Rottler2002a . Both limiting behaviors and the crossover between them can be addressed through our simulations. As can be seen in Fig. 7, short chains of length $`N=128`$ form crazes that grow at the constant plateau stress $`S`$, but the stress drops monotonically to zero upon straining past $`\lambda `$. For longer chains, the stress $`\sigma _{zz}`$ rises to a maximum value $`S_{\mathrm{max}}`$ that exceeds $`S`$.

Values for $`S_{\mathrm{max}}`$ were systematically obtained as a function of normalized chain length $`N/N_e`$ from curves such as those shown in Fig. 7. Fig. 22 summarizes the breaking stresses for the craze fibrils normalized by the breaking stress in the limit of very long chains $`S_{\mathrm{}}`$. $`S_{\mathrm{max}}`$ is zero for $`N<2N_e`$, since stable crazes do not form for such short chains. $`S_{\mathrm{max}}/S_{\mathrm{}}`$ first rises roughly linearly with $`N/N_e`$, and then saturates at unity for chain lengths longer than about $`10N_e`$. The saturation coincides with the observation of significant amounts of chain scission. Interestingly, the data seems to collapse onto a single curve (solid line).

Note that the maxima of the stress-strain curves in Fig. 7 are reached at strains of $`6`$ and $`10`$ for flexible and semiflexible chains, respectively. These values are close to $`\sqrt{3}\lambda _{\mathrm{max}}`$, which implies that at the breaking point the chains are pulled fully taut between entanglement points. This was confirmed by direct analysis of the craze microstructure.

#### VII.1.1 Chain end relaxation

In order to understand the crossover regime and the competition between the two failure mechanisms, it is useful to study the distribution of tension along a given chain. Fig. 23 shows the tension as a function of distance from the chain end at several stages of craze breaking. Since the chain ends are identical, symmetry was used to improve statistics. In the unstrained craze (lowest curves), both flexible and semiflexible systems exhibit a constant stress in the center of the chain, but a relaxation toward the free ends. The characteristic length scales for this relaxation were extracted by fitting an exponential decay to the transition region. The values of the decay lengths $`N_{\mathrm{end}}^{\mathrm{fl}}=21`$ and $`N_{\mathrm{end}}^{\mathrm{sfl}}=13`$ are comparable to the characteristic length of $`N_e/3`$ for stretched segments, but are not universal. Stronger adhesive interactions were found to increase $`N_{\mathrm{end}}^{(\mathrm{s})\mathrm{fl}}`$.

Upon straining the craze, the tension in the center of the chains and the values of $`N_{\mathrm{end}}^{(\mathrm{s})\mathrm{fl}}`$ rise. At the breaking point (last curves), the end relaxation extends over a length scale comparable to the entanglement length.

These results help to formulate a simple argument for the universal curve plotted in Fig. 22. The average distance of an entanglement point from the chain end is $`N/4`$. We assume that the probability of disentanglement decreases exponentially with distance from the chain end, as suggested by the tension relaxation curves. The characteristic length scale $`N_{\mathrm{end}}`$ at the breaking point in these curves was on the order of $`N_e`$. Because of the above, we expect this length scale to be the characteristic decay length for the probability of disentanglement, and thus postulate a disentanglement probability of the form $`\mathrm{exp}[(N2Ne)/4N_e]`$. Here, $`N`$ was reduced by $`2N_e`$, since for this chain length the disentanglement probability is one and the chain is free on either side. The maximum stress can now be written as the limiting value of $`S_{\mathrm{}}`$ times the probability for non-disentanglement, which gives

$$S_{\mathrm{max}}/S_{\mathrm{}}=1\mathrm{exp}[(N2Ne)/4N_e].$$
(19)
Fig. 22 shows that this curve agrees well with the data.

### VII.2 Global tension distribution

The parameter governing chain scission and thus the value of $`S_{\mathrm{}}`$ is the distribution of tension in the polymer craze. In a previous paper Rottler2002b , we reported that this distribution is characterized by an exponential tail at large tensile forces, in analogy to jammed systems such as granular materials Liu2001 . This distribution is shown in Fig. 24 for flexible and semiflexible chains of length $`N=512`$ and several strain states. The tensile (positive) part of the distribution is well fitted by $`1/f\mathrm{exp}[f/f]`$, where only the positive tensions are included in the average $`f`$. Note that $`f`$ is the same for flexible and semiflexible chains at the plateau. The distribution with the steepest slope (smallest $`f`$) corresponds to the fully developed craze. Additional curves with higher $`f`$ correspond to stressed crazes at the same strain with respect to the unstrained craze. Note that the semiflexible and flexible crazes have the same values of $`f`$ at each strain. This is related to the fact that the stress $`S_0`$ in the fibrils is independent of $`N_e`$ (see Section V.6). The last curve shows the tension distribution at the breaking point where $`\sigma _{\mathrm{zz}}`$ is largest (see also Fig. 7). The effect of straining the craze only changes $`f`$, and all curves could be collapsed after normalizing by $`f`$.

We note that in our simulations, no scission or disentanglement occurs during craze growth proper. The fraction of bonds that break at a given average tension is $`P_{br}(f)=_{f_c}^{\mathrm{}}\mathrm{exp}[f/f]𝑑f=\mathrm{exp}[f_c/f]`$. In our simulations $`f_c=100f_{\mathrm{LJ}}`$ as described in Section III. The onset of scission can be estimated using a simple scaling argument: $`N_{bonds}P_{br}(f)1`$. From this we estimate an average value of the tension at breaking

$$f=f_c/\mathrm{ln}N_{bonds}.$$
(20)
For $`N_{bonds}262144`$, this implies that scission will occur when $`f8.0f_{\mathrm{LJ}}`$, which was confirmed by direct inspection of the chains at the corresponding strains. Such high tensions only occur when the craze is strained past the extension ratio. The largest value of $`f`$ observed with the present model during craze growth was $`4.5f_{\mathrm{LJ}}`$ and occurred at very low temperature $`T=0.01u_0/k_B`$ and $`r_c=2.2a`$.

The degree of chain scission in experimental crazes is still a matter of debate, but it appears likely that at least some chains do experience scission. The absence of scission in the present study is most likely due to the low monomeric friction coefficient of the bead-spring model. As the above argument showed, a relatively modest increase in the average tension will quickly lead to appreciable scission. More realistic polymer models should be able to capture this effect. An increase in system size will likewise raise the number of broken bonds. For a typical value of $`f=3f_{\mathrm{LJ}}`$, one bond would break for every $`10^{15}`$ bonds. Note that the exponential tension distribution leads to a logarithmic size dependence Eq. (20) and allows for sequential bond breaking. The fibrils are thus much weaker than implied by the common simple assumption that all bonds carry the same tension and break when $`f=f_c`$.

## VIII Summary and Conclusions

This paper presented molecular dynamics simulations of craze nucleation, widening and breakdown. Initial failure of the LJ polymer glass occurred through shear in biaxial loading. Only when all three principal stresses were tensile did cavitation and craze formation occur. However, once past the nucleation phase, plane stress conditions are sufficient for continuing craze growth. Cavitational failure could be fitted to a cavitation criterion of the form $`\tau _{dev}^c=\tau _0^c+\alpha ^cp`$.

Craze widening proceeds in the simulations by a clearly identifiable fibril drawing process. This interpretation is also well-supported by experiments. The resulting craze microstructure is compellingly similar to TEM images of experimental crazes, and the length scales quantified by $`D`$ and $`D_0`$ are within experimental limits. The simulations clearly establish the connection between extension ratio and entanglement length. In the glass, disentanglement is prevented and the entanglements act like chemical crosslinks. A microscopic analysis of the length of stretched chain segments has shown that unlike the case of a regular mesh, only a few segments are fully expanded to the entanglement length, and the average extension is only $`N_e/3`$. The factor 1/3 arises from averaging over all angles that a given segment can form with the stretching direction.

Another salient finding of this study is the exponential distribution of tension in the craze. The presence of large stress fluctuations makes chain scission much more likely than e. g. a Gaussian distribution or uniform loading. Since force distributions of this kind are also often seen in conventional “jammed” systems such as foams, colloids and granular media, we have suggested Rottler2002b  that a craze can be viewed as a system that jams under tension.

The highly nonequilibrium nature of the force distribution, and the strong concentration of stress in the covalent backbone bonds, formed the basis for our critique of the conventional capillary model of craze widening. The polymer glass is not a viscous fluid in the active zone, and the hydrodynamic description does not apply. The picture suggested by our simulations is that crazing is a form of localized shear deformation, but with a much greater mobilization of material than in standard shear yielding. The very similar rate and temperature dependence of $`S`$ is another indication of the close relation between the processes. Based on trends observed in the simulations, we have suggested that the local stress in the fibrils $`S_0=S\lambda `$ is independent of the entanglement length. $`S_0`$ varies with temperature and strength of adhesive interaction in a manner very similar to the yield stresses for shear and cavitation. Establishing a precise connection between these characteristic stresses should be a most interesting direction for future work.

A detailed analysis of the microstructure of crazes was also presented. The calculated structure factor is similar to measured scattering intensities. As in these experiments, a Porod analysis was used to extract a measure of the mean fibril diameter $`D`$ from the structure factor. While the extension ratio depends only on $`N_e`$, the mean fibril diameter depends on many factors. The value of $`D`$ increases with increasing $`T`$ and with increasing strength of the van der Waals interactions. Chain stiffness has less effect, although $`D`$ is larger for flexible chains than semiflexible chains at high temperatures.

The distribution of fibril diameters was determined from the real space structure of the crazes. The average fibril diameter from this method was always smaller than that determined from the structure factor. The distribution was also very wide with a variance that exceeded the mean and a tail extending to many times the mean. The radial distribution function for the fibrils shows almost no correlation, merely an exclusion minimum near the origin. Fibrils merge and split with each other directly, rather than being joined by smaller cross-tie fibrils.

The simulations described here capture the generic features of experiments on many different polymers and provide previously inaccessible information about the dynamics and microstructure. However, they are unable to address quantitative behavior of specific polymers. Future studies with chemically realistic potentials will be of great value, but require orders of magnitude more computational effort.

## IX Acknowledgements

We are indebted to E. J. Kramer and H. R. Brown for very insightful discussions of this work. Financial support from the Semiconductor Research Corporation (SRC) and NSF grant No. DMR0083286 is gratefully acknowledged. The simulations were performed with LAMMPS 2001 Lammps2001 , a molecular dynamics package developed by Sandia National Laboratories.


format: text

text:
# Chaos in Turbulence Driven by the Magnetorotational Instability

## 1 Introduction

Astrophysical accretion disks are able to evolve because angular momentum is extracted from fluid elements and transported outward. This is effected by the presence of nonvanishing radial-azimuthal components of the Maxwell and Reynolds stress tensors, produced by magnetohydrodynamic (MHD) turbulence driven by the magnetorotational instability (MRI) (Balbus & Hawley, 1998). As there is no analytic theory of MHD turbulence at hand (nor is there one in sight), large-scale numerical simulation has been the main avenue of progress toward understanding its properties.

Many numerical studies make use of the local “shearing box” approximation (Hawley, Gammie & Balbus, 1995). The shearing box is an invaluable tool for studying local flow dynamics in detail, and for resolving turbulent flow with the largest possible dynamical range. MHD turbulence, like all true turbulence, should be chaotic, as quantified formally by a measured positive Lyapunov exponent. What is less clear, but of great astrophysical significance, is whether long term flow averages are even well-defined. To put the question most starkly, imagine macroscopically identical disks, with fluid perturbations that vary by a tiny amount. The fine scale description of their internal turbulence will surely differ, but will quantities such as the stress tensor components converge to the same values in the long term? This question goes directly to the heart of the standard $`\alpha `$ disk formalism (Shakura & Sunyaev, 1973) which assumes that disk transport may be described by a spatially constant or slowly varying $`\alpha `$ parameter. (The quantity $`\alpha `$ is defined as the ratio of radial-azimuthal component of the stress tensor $`T_{R\varphi }`$ to the gas, or gas plus radiation, pressure.) How well supported is this assumption?

This paper begins to study these complex issues by examining the chaotic nature of the MHD turbulence. The presentation is organized as follows. In §2 we briefly review the shearing box. In §3 we present the results of experiments designed to reveal how the measured properties of turbulent flow are related to both computational and physical input parameters. In §4, we carry out a series of experiments that demonstrate qualitatively that the MHD turbulence is chaotic, and in §5 this is quantified by computing the Lyapunov exponents for a set of simulations. Our conclusions are summarized in §6.

## 2 Shearing Box System

The shearing box system, described in Hawley, Gammie & Balbus (1995), is designed to represent a very local section of an accretion disk, viewed in corotating coordinates with angular frequency $`\mathrm{\Omega }`$. Starting with the full set of dynamical equations in cylindrical coordinates, the equations are locally expanded about a fiducial cylindrical radius $`R`$, with $`(dx,dy,dz)`$ corresponding to cylindrical coordinates $`(dR,Rd\varphi ,dz)`$. The computational domain is then a Cartesian box, but with the rotational inertial forces (Coriolis and centrifugal) retained. The centrifugal force nearly balances gravity, leaving a remnant tidal force linear in $`x`$. All other forces are directly retained. For this system, the ideal MHD equations of motion become,

$`{\displaystyle \frac{𝐯}{t}}+\left(𝐯\right)𝐯={\displaystyle \frac{1}{\rho }}\left(P+{\displaystyle \frac{B^2}{8\pi }}\right)+{\displaystyle \frac{\left(𝐁\right)𝐁}{4\pi \rho }}`$
$`2𝛀\times 𝐯+2q\mathrm{\Omega }^2x\widehat{𝐱}\mathrm{\Omega }^2z\widehat{𝐳}.`$ (1)
The mass conservation equation

$$\frac{\rho }{t}+\left(\rho 𝐯\right)=0,$$
(2)
the adiabatic internal energy equation,

$$\frac{\rho ϵ}{t}+\left(\rho ϵ𝐯\right)+P𝐯=0,$$
(3)
and the induction equation

$$\frac{𝐁}{t}=\times \left(𝐯\times 𝐁\right),$$
(4)
retain their usual forms. The internal energy per particle $`ϵ`$ is defined by

$$P=\rho ϵ\left(\gamma 1\right).$$
(5)
For an isothermal gas, the energy equation is replaced with

$$P\rho ^1=\mathrm{constant}.$$
(6)
The variables in these equations have their usual meanings. We shall also introduce the constant $`q`$, which is a parameter describing the local radial dependence of the angular frequency, i.e., $`q=d\mathrm{ln}\mathrm{\Omega }/d\mathrm{ln}R`$. For a Keplerian angular momentum distribution, $`q=3/2`$. In this study we have ignored vertical gravity, and have thus dropped the $`\mathrm{\Omega }^2z\widehat{𝐳}`$ term from equation (2).

These equations are solved using the time-explicit, operator-split finite differencing algorithm of the ZEUS code for hydrodynamics (Stone & Norman, 1992a) and MHD (Stone & Norman, 1992b; Hawley & Stone, 1995). The shearing box computational domain extends in $`x`$, $`y`$, and $`z`$ respectively from $`L_x/2`$ to $`+L_x/2`$, 0 to $`L_y`$, 0 to $`L_z`$. The boundary conditions are periodic in both the $`y`$ and $`z`$ directions, and “shearing periodic” in $`x`$, meaning we equate

$$f(x,y,z)=f(x+L_x,yq\mathrm{\Omega }L_xt,z).$$
(7)
The azimuthal velocity has an additional correction to offset the angular velocity difference between the inner and outer radial boundaries:

$$v_y(x,y,z)=v_y(x+L_x,yq\mathrm{\Omega }L_xt,z)+q\mathrm{\Omega }L_x.$$
(8)
The initial state is one of uniform density and pressure, and the velocity profile is $`𝐯=q\mathrm{\Omega }x\widehat{𝐲}`$. In this paper, our initial magnetic field configurations are either a uniform toroidal magnetic field, or a vertical magnetic field varying sinusoidally in the radial $`x`$ direction. The magnetic field strength is described by the standard plasma $`\beta `$ parameter

$$\beta \frac{8\pi P}{B^2},$$
(9)
the ratio of the gas to magnetic pressures. The MRI is triggered by seeding the initial equilibrium state with small pressure and velocity fluctuations.

## 3 Saturation Amplitude in the Shearing Box: a Survey

One of the goals of shearing box simulations has been to try to understand what determines the sustained saturation amplitude of the MRI-driven turbulence. In accretion disk applications, the focus is often on the $`T_{r\varphi }`$ component of the combined Reynolds and Maxwell stress tensor,

$$T_{r\varphi }=\left(B_xB_y/4\pi +\rho v_x\delta v_y\right).$$
(10)
Here $`\delta v_y`$ is the azimuthal velocity fluctuation, obtained by subtracting the equilibrium shear flow. This can be directly measured in shearing box simulations. The volume-averaged stress is usually normalized to the averaged gas pressure, a quantity equivalent to the (Shakura & Sunyaev, 1973) $`\alpha `$ parameter.

Even with a system as simple as the shearing box, there are many possible factors, both computational as well as physical, that could in principle govern the saturation amplitude. These include the initial field strength and geometry, box size and aspect ratio, value of $`\mathrm{\Omega }`$, background gas pressure, and numerical resolution. Many of these are physically significant only in relation to others. For example, $`\mathrm{\Omega }`$ is an arbitrary time scale in the system. The shearing box system is independent of $`\mathrm{\Omega }`$ so long as $`q`$ remains constant and the physical box size changes proportionately. The formal scaling invariance is that if

$$\mathrm{\Omega }c\mathrm{\Omega },$$
(11)
where $`c`$ is an arbitrary constant, then the equations remain unchanged if

$$(x,t)(x/c,t/c).$$
(12)
We will make use of this fundamental property in §4 when we explore chaotic behavior.

For the cases of an initial uniform field oriented along either the $`y`$ or $`z`$ axes, Hawley, Gammie & Balbus (1995) found that the saturation amplitude of the turbulence depended upon both the box size and initial field strength. When the field had a vanishing mean value over the computational box volume, the final outcome was much less sensitive to the initial field energy. Here, we present a series of zero mean field simulations designed to test how the box size affects the overall level of the turbulence. The baseline calculation has $`q=1.5`$ corresponding to a Keplerian background, and a box size of $`L_x=1`$, $`L_y=2\pi `$, and $`L_z=1`$, set so that $`L_z=c_s^2/\mathrm{\Omega }`$ where $`c_s^2=P/\rho `$. The actual values used are $`P=10^6`$ and $`\mathrm{\Omega }=0.001`$. The gas is adiabatic with $`\gamma =5/3`$. The initial field is vertical and varies sinusoidally in the $`x`$ direction, $`B_z\mathrm{sin}(x/L_x)`$. The field strength set so that the volume average magnetic energy corresponds to $`\beta =800`$. The maximum vertical field has $`\beta =400`$, which gives a fastest growing wavelength for the MRI of $`\lambda _{\mathrm{max}}=8\pi /15^{1/2}(v_A/\mathrm{\Omega })=0.46`$. The grid resolution is $`32\times 64\times 32`$ zones in $`(x,y,z)`$. This setup corresponds to the so-called “fiducial box” used in Hawley, Gammie & Balbus (1995, 1996) and Brandenburg et al. (1995).

In addition to the fiducial run, we computed four additional simulations with different box dimensions. In three simulations we doubled each of the box dimensions in turn, and in the final simulation of this set we quadrupled the vertical size. To maintain a constant resolution number of grid zones per unit length, the number of grid zones used was increased appropriately. Figure 1 shows the magnetic energy per unit volume evolving over 30 orbits for each run in this series of models. While one must be cautious extrapolating from this limited time baseline, the results generally suggest that increasing the vertical dimension is the most important effect over the long term. The temporal variations in magnetic energy within a given run are, however, quite large, often well in excess of the running average.

To isolate the effect of the vertical size, we carried out a second set of simulations in which $`L_z`$ is varied, from $`L_z=1/4`$ to $`L_z=4`$, and all other factors are held fixed, including the grid zone size $`\mathrm{\Delta }z`$. The initial magnetic field is as described above. Here we used an isothermal equation of state, $`\gamma =1.0`$ which keeps the box from heating due to dissipation of the turbulence. Figure 2 shows the time evolution of the magnetic energy and even over a limited time baseline clearly reveals a general trend toward increasing magnetic energy with increasing $`z`$ box size, again with significant variation. The most striking conclusion is that there is a lower bound to the vertical size required for amplification to occur: the field energy declines dramatically in the $`L_z=1/4`$ model, and even the initial linear growth rate is reduced. The fastest growing vertical wavelengths for initial field strengths with $`\beta <1350`$ no longer fit inside the box, and evidently losses win over the reduced power input at the top of the turbulent cascade.

For box sizes $`L_z=1/2`$ and larger the turbulence is sustained over the evolution period of 30 orbits. However, other trends are visible. First, the amplitude of the initial spike of linear growth saturation decreases with increasing $`L_z`$. Second, the amplitude of the magnetic energy variations in the turbulent state is reduced with increasing box size.

These history plots demonstrate that it is not a simple matter to quantify the saturation amplitude of the turbulence in a shearing box, or to determine unambiguously the influence of physical parameters in setting that amplitude. It is not obvious, in general, over what length of time one should average to obtain a characteristic amplitude; nor is it even obvious that such a baseline must always exist. Large fluctuations, even in volume-averaged data, are the rule. This is the overarching and defining characteristic of the MHD turbulence: it is chaotic.

## 4 Onset of Chaos

It is clear that changing the physical parameters of a shearing box run can result in different levels of turbulence. Here we show that almost exactly the same physical parameters, but with tiny variations in numerical value, leads to macroscopically different levels of turbulence. The sensitivity of the system to infinitesimal perturbations is a key feature of chaos.

The first set of simulations consists of three shearing boxes with different $`\mathrm{\Omega }`$ values, but which have been rescaled so as to be physically equivalent. All three have $`q=1.5`$, $`\beta =800`$, $`\gamma =1.0`$, $`B_z\mathrm{sin}(x/L_x)`$, and $`32\times 64\times 32`$ grid zones. One run uses the standard box size of $`1\times 2\pi \times 1`$ with $`\mathrm{\Omega }=0.001`$. The other two have $`\mathrm{\Omega }=0.0031`$ and $`\mathrm{\Omega }=0.000507`$ with the physical box sizes scaled accordingly.111Care must be taken so that the rescaling does not correspond to a simple shift by an integer number of bits in the machine registers. An identical sinusoidal velocity perturbation was applied to the initial state in each box. In principle, since the length to time ratio is preserved between these different shearing box systems, the resulting evolution should be identical. In the computational experiment, however, these scaled systems evolved very differently.

Figure 3 shows the volume averaged magnetic energy density histories for the three runs. The initial linear growth and saturation phase is indistinguishable. Significant deviations appear after orbit 8, and at orbit 20 the energies are widely different. Since these simulations correspond to identical physical conditions, one might expect that, at the very least, time-averaged values would coincide. Figure 4 illustrates how the running time averages of the volume-averaged stress parameter $`\alpha `$ for each run evolves with time. The running time average is defined

$$\alpha (t)=\frac{1}{(tt_o)}_{t_o}^t\alpha 𝑑t,$$
(13)
and the expectation is that this should converge to a steady value characteristic of the turbulence. The question is, over what time interval should this happen? As figure 4 shows, the fluctuations in the stress are so large that the running time averages $`\alpha (t)`$ do not coincide. Since a very large excursion in magnetic energy is visible in figure 3 at 20 orbits, we computed the running average from 45 to 120 orbits. But over this interval a unique characteristic $`\alpha `$ does not emerge even for a single simulation.

The next experiment further demonstrates the chaotic nature of the turbulence by taking the data halfway through an extended run and randomly perturbing the flow velocities with a gaussian rms value of 0.001%. The subsequent evolution is then compared with the extended run. The parameters for this test were $`\mathrm{\Omega }`$ = 0.0005, $`q=1.5`$, $`\beta =400`$, $`\gamma =1.0`$, $`B_z\mathrm{sin}(x/L_x)`$, grid size of $`32\times 64\times 32`$, and an appropriately scaled box equivalent to the standard box. The results are shown in figure 5. The arrow indicates the time when the perturbations were added, and the perturbed simulation was then run on to 60 orbits in time. The plot of magnetic energy reveals a visible divergence by orbit 30, and a substantial difference by orbit 35. This is the defining feature of chaos: substantially different macroscopic behavior created by infinitesimal perturbations.

To compare the behavior of models with different initial field topologies, we perform an experiment with an initial uniform toroidal magnetic field for two boxes, one with $`\mathrm{\Omega }=0.001`$ the other with $`\mathrm{\Omega }=0.0031`$. Other physical parameters were the same as in the vertical field experiments. Figure 6 displays the volume-averaged magnetic energy density histories for these runs. Again, the initial linear growth phases of the two runs were identical; diverging behavior is a nonlinear phenomenon. The amplitude of the variations is not as great as in the vertical field case. Nevertheless, as figure 7 shows, the running time averages of $`\alpha `$ required 100 orbits to converge.

As a contrast to the chaotic behavior of MHD turbulence, we also ran several purely hydrodynamic shearing box systems. Hydrodynamic shearing boxes are linearly stable when $`q<2`$ (the Rayleigh criterion), and nonlinearly stable when $`q<2ϵ`$, where $`ϵ`$ is a number on order $`0.01`$ (Hawley, Balbus, & Winters, 1999). All of the $`q=1.5`$ hydrodynamical simulations gave similar results. The volume-averaged kinetic energies die out in an identical fashion for equivalent problems that use different $`\mathrm{\Omega }`$ values. Small fluctuations exist in the form of waves, and are typical of decaying hydrodynamic turbulence. These waves were reproduced almost identically from simulation to simulation—the measured differences were consistent with round-off error. Additional experiments with different $`q`$ parameters (up to $`q=1.95`$) were carried out, and all of these stable hydrodynamical experiments remained identical in their macroscopic properties. The numerical differences were small, not increasing, and consistent with machine precision. These results are consistent with the conclusion that the observed chaos is indeed a property of the MHD turbulence, not of the shearing box itself or the numerics.

### 4.1 Lyapunov Exponents

The experiments provide compelling, but qualitative evidence that MRI-driven MHD turbulence is chaotic. A more quantitative parameterization of chaotic behavior is afforded by the use of Lyapunov exponents. A system with $`N`$ degrees of freedom has $`N`$ Lyapunov exponents. The largest positive exponent represents the average rate of divergence of two initially close evolution paths in phase space. The evolution paths are traced out by monitoring the state vector of the system, which we define as

$$𝐯=(v_x,v_y,v_z,\frac{B_x}{\sqrt{4\pi \rho }},\frac{B_y}{\sqrt{4\pi \rho }},\frac{B_z}{\sqrt{4\pi \rho }}).$$
(14)
The density itself is generally of secondary importance, and therefore is tracked only indirectly in the Alfvén velocity. Close evolution paths imply that there are small differences between the state vectors for each of the two trajectories. If an initially small difference diverges exponentially between two state vectors, the system is formally chaotic. This corresponds to a positive Lyapunov exponent.

Kurths & Brandenburg (1991) demonstrated chaos in MHD solar convection simulations in three-dimensions and computed estimates of the Lyapunov exponents. We follow their procedure to perform a similar analysis for MRI-driven turbulence. We begin with an evolution in which MHD turbulence has fully developed. Then, a second evolution path is created by randomly perturbing the velocity components of the first system’s state vector. The two state vector trajectories are then evolved for about 4 orbits, and, at regular time intervals the fractional state vector difference between the two paths was calculated:

$$\delta =\frac{|𝐕_𝐩𝐕|}{|𝐕|}$$
(15)
where $`𝐕_𝐩`$ and $`𝐕`$ are the perturbed and original state vectors. Data from the last two orbits in the state vector difference history are then fit to an exponential from which we estimate the value of the largest Lyapunov exponent.

In the first experiment we employ a standard shearing box, namely $`1\times 2\pi \times 1`$, grid size of $`32\times 64\times 32`$, with an isothermal equation of state, $`\gamma =1.0`$, sinusoidal vertical field, $`B_z\mathrm{sin}(x/L_x)`$, of average magnetic energy $`\beta =800`$, and a background Keplerian shear, $`q=1.5`$. This was evolved until MHD turbulence had fully developed. Perturbations were then applied to a parallel computation and the state vector difference was calculated as a function of time. The resulting Lyapunov exponent is $`0.458\omega _{\mathrm{max}}`$, where $`\omega _{\mathrm{max}}=0.75\mathrm{\Omega }`$ is the maximum unstable growth rate for the linear phase of magnetorotational instability (Balbus & Hawley, 1998). All of our experiments yield a Lyapunov exponent comparable to $`\omega _{\mathrm{max}}`$.

That the Lyapunov exponent would be on order of the maximum MRI growth rate is not surprising. It is precisely the linearly unstable, exponentially growing MRI that is feeding the turbulence, and driving exponential divergence of the state vector. To examine this more carefully, a series of experiments was performed on shearing boxes with a variety of background shear $`q`$ parameters and for both sinusoidal vertical and toroidal initial field configurations. Because the maximum linear growth rate of the MRI is $`q/2`$ (Balbus & Hawley, 1998), the ensemble of simulations spans an interesting range of MRI growth rates. We would expect the Lyapunov exponent to be proportional to $`q`$ as well. Figure 8 displays the volume averaged, percent state vector difference histories for these runs. The curves are labeled by their $`q`$ values, with solid lines for the vertical field runs and dashed lines for the toroidal field runs. The corresponding first Lyapunov exponent values normalized by $`\omega _{\mathrm{max}}`$ are $`0.521`$, $`0.458`$, $`0.422`$ and $`0.583`$ for the vertical field cases, in order of descending $`q`$ value, and $`0.484`$ and $`0.644`$ for the toroidal field $`q=1.5`$ and $`q=1.3`$ runs. Clearly, the first Lyapunov exponents are positive, and of the same magnitude, normalized by $`\omega _{\mathrm{max}}`$. It should also be noted that first Lyapunov exponents were calculated at many points in time in the simulations, always with similar results.

In summary, the Lyapunov exponents in MRI-driven MHD turbulence simulations are all positive, and, when normalized by $`\omega _{\mathrm{max}}`$, they all lie within in a range of 0.4 to 0.6. There is a slight trend of larger Lyapunov exponent for larger $`q`$ value (non-normalized). Previous experiments also found stronger overall levels of turbulence with larger $`q`$ values (Hawley, Balbus, & Winters, 1999); in a sense, a larger Lyapunov exponent is “more turbulent.” Finally, the chaos parameters of the turbulence are independent of the initial magnetic field configuration. Initial vertical magnetic fields have similar Lyapunov exponent values as initial toroidal magnetic fields.

## 5 Conclusions

Shearing box simulations of the MRI constitute an excellent numerical laboratory in which to study chaos in turbulent flow. Their compactness and simple boundary conditions make them a very convenient system to study, but they also require care to interpret. For example, the variance of the the flow fluctuations, which may be of direct astrophysical interest because of its connection with radiative emission, is a function of the box size adopted (cf. figure 2).

In this paper, we have demonstrated the extreme sensitivity of MHD turbulence to infinitesimal deviations in the flow. This was done by several different methods: showing that invariant scaling laws fail when implemented numerically, for both vertical and toroidal initial fields, and externally imposing tiny perturbations on an established turbulent flow and following the growing deviations in the subsequent evolution of the original and the perturbed system. Estimates of the largest Lyapunov exponent in a variety of turbulent flows with different field geometries yielded values near the characteristic growth rate of the linear MRI. This is an indication that the linear physics of the instability plays an active role in defining the highly nonlinear turbulent dynamics of these flows. One way this could come about would be if the energy input into a Kolmogorov-like cascade was essentially the linear MRI.

The most important practical consequence of this behavior is that the nongaussian statistical properties of chaotic flows severely limit the extent to which $`\alpha `$ modeling can be used uncritically. Though Maxwell and viscous stress have some formal properties in common (Balbus & Papaloizou, 1999), the averaging procedure necessary for a semi-local treatment of the turbulence is a very delicate matter. A time base of hundreds of orbits is clearly necessary to establish a meaningful estimate of the characteristic stress. In astrophysical systems, especially those in transience, it may not be possible to ascribe an instantaneous $`\alpha `$ value to the stress, and there may be no recourse other than detailed numerical modeling.

## Acknowledgements

We acknowledge support under NSF grant AST-0070979, and NASA grant NAG5-9266. Some of the simulations described here were carried out on computational platforms at the NSF-supported National Center for Supercomputing Applications and the San Diego Supercomputer Center.