text:
# Shear yielding of amorphous glassy solids: Effect of temperature and strain rate

## I Introduction

Deformation processes and plasticity in amorphous materials such as metallic or polymeric glasses have recently received a lot of attention Larson1999 ; Rottler2001 ; STZ ; Lemaitre2001 ; Hasan1995 . These materials are used in many load-bearing applications. However, an understanding of their yield and flow properties is hampered by the absence of long-range order and easily identifiable mechanisms that mediate the deformation such as dislocation motion in crystals. A similar situation is encountered in “soft” glassy materials such as foams, pastes, and colloidal suspensions, which are also characterised by a liquid-like structure and long relaxation times. In fact, it has been suggested recently that these very different materials can be viewed as particular realizations of a jammed state Liu2001 , which implies that their mechanical behavior could be described in a common framework.

Much insight into the mechanical behavior of structural glasses has been gained from molecular simulations of simple glass-forming liquids or polymers, where particles interact through a Lennard-Jones potential. For example, Falk and Langer STZ  studied two-dimensional shear deformation of a mixture of such particles and found localization of plastic events in so-called shear transformation zones. Barrat and Berthier Barrat2000  studied the steady state flow of a similar model and analyzed its relation to the fluctuation-dissipation theorem in out-of-equilibrium situations.

In a recent paper Rottler2001 , we have studied the onset of shear yielding of amorphous polymer glasses under multiaxial loading conditions. A pressure-modified von Mises criterion Ward1983  accurately describes the maximum shear yield stress as a function of the applied stress for different temperatures. However, in these simulations the bulk polymer was deformed at a single constant strain rate. The present work investigates the effect of strain rate on the stress at the onset of shear yielding as well as on the steady state flow stress of glassy materials in simple shear. In order to relate the onset of yield to steady shear, we not only discuss the macroscopic material response, but also perform an analysis of the localized dynamics of the stress distribution in the deforming solid. Such a program is particulary suited to test predictions of models of plasticity Eyring1936 ; STZ ; Lemaitre2001 ; SGR ; Berthier2000  and may lead to a deeper understanding of the nature of plastic deformation.

In the following section, we briefly summarize traditional and more recent models of viscoplasticity and rheology of glassy materials. Sections III and IV discuss the molecular models and the simulation results, respectively. Section V critically reviews how the theoretical ideas of Section II describe the simulation results.

## II Models of plasticity and rheology

### II.1 Rate and temperature dependence: Eyring model

The simplest model that makes a prediction for the rate and temperature dependence of shear yielding is the rate-state Eyring model Larson1999 ; Eyring1936  of stress-biased thermal activation. Structural rearrangement is associated with a single energy barrier $`E`$ that is lowered or increased linearly by an applied stress $`\sigma `$. This defines transition rates of the form

$$R_\pm =\nu _0\mathrm{exp}\left[\frac{E}{k_BT}\right]\mathrm{exp}\left[\pm \frac{\sigma V^{}}{k_BT}\right],$$
(1)
where $`\nu _0`$ is an attempt frequency and $`V^{}`$ is a constant called the “activation volume.” In glasses, the transition rates are negligible at zero stress. Thus at finite stress one needs to consider only the rate $`R_+`$ of transitions in the direction aided by stress. The plastic strain rate $`\dot{ϵ}_{pl}`$ will be proportional to $`R_+`$, $`\dot{ϵ}_{pl}=cR_+`$. Solving for the stress $`\sigma `$, one obtains

$$\sigma =\frac{E}{V^{}}+\frac{k_BT}{V^{}}\mathrm{ln}\left[\frac{\dot{ϵ}_{pl}}{c\nu _0}\right].$$
(2)
Eq. (2) contains only a single relaxation time scale and predicts an apparent yield stress that varies logarithmically with the strain rate and a prefactor that depends linearly on temperature. Despite its simplicity, experimental results Ward1983  are often fitted to Eq. (2), and the value of $`V^{}`$ is associated with a typical volume required for a molecular shear rearrangement.

### II.2 Modern approaches

Modern phenomenological approaches pay tribute to the complexity of glassy systems through several extensions. First, it has been realized that assuming a single energy barrier for rearrangements is an oversimplified description of glassy materials Hasan1995 . One can therefore introduce a distribution of barriers and add additional time scales. Second, any theory that attempts to predict a full stress-strain curve must contain some information about the internal state of the system as a function of time or strain. Extensions therefore consider dynamical internal state variables. In the following, we describe particular realizations of these ideas.

#### II.2.1 Shear transformation zone theory

Falk and Langer used molecular dynamics simulations very similar to the present study to identify local plastic rearrangements under shear. They formulated a theory of viscoplasticity STZ  based on the concept of “shear transformation zones” (STZ), bistable (mesoscopic) regions that transform under shear between $`\pm `$-states. One then considers the dynamics of an ensemble of STZ with number density $`n_\pm `$ on a mean-field level, which determines the plastic strain rate

$$\dot{ϵ}_{pl}=A_0(R_+n_+R_{}n_{}),$$
(3)
where $`A_0`$ is a constant. A key difference from the Eyring model resides in the form of the transition rates $`R_\pm `$, which are assumed to be free-volume (entropically) activated rather than thermally activated, i.e.

$$R_\pm =R_0\mathrm{exp}\left[\frac{v_0\mathrm{exp}[\sigma /\overline{\mu }]}{v_f}\right],$$
(4)
where $`v_0`$ is a characteristic free volume required for a STZ flip and $`\overline{\mu }`$ a characteristic stress scale required for a molecular rearrangement. The role of temperature is played by a “free volume” $`v_f`$ per particle. The authors motivate this with the observation that in a solid at very low temperature, energy barriers should be very large compared to thermal energies and thus, as in granular systems, thermal activation over these barriers should be negligible. The population densities themselves evolve according to the rate equation

$$\dot{n}_\pm =R_{}n_{}R_\pm n_\pm +\sigma \dot{ϵ}_{pl}(A_cA_an_\pm ),$$
(5)
where the last term introduces creation and annihilation processes of STZ’s proportional to the work of plastic deformation $`\sigma \dot{ϵ}_{pl}`$. The STZ equations can be solved analytically in special steady state situations and otherwise solved numerically. They were shown STZ  to have both a jammed solution for which $`\dot{ϵ}_{pl}=0`$ and a flowing solution once $`\sigma `$ exceeds a true yield stress $`\sigma _y`$. The shear rate rises linearly as sigma increases above $`\sigma _y`$ as in a Bingham fluid. Numerical stress-strain curves, hysteresis experiments and creep tests were shown to be accurately reproduced by the model after the adjustment of several fit parameters.

Lemaître Lemaitre2001  recently extended STZ theory with concepts from the physics of granular media. He treated the free volume parameter $`v_f`$ as a dynamical state variable and proposed the following time evolution:

$$\dot{v}_f=R_1\mathrm{exp}\left[\frac{v_1}{v_f}\right]+A_v\sigma \dot{ϵ}_{pl}.$$
(6)
This expression is motivated by slow density relaxations in granular materials that decrease the free volume. The activation factor $`\mathrm{exp}[v_1/v_f]`$ describes the probability for volume fluctuations larger than a characteristic volume $`v_1`$. Note that the “activation barrier” for compaction $`v_1`$ differs a priori from the barrier for shear transformation $`v_0`$. The second term refers to creation of free volume (shear induced dilatancy) again due to plastic deformation. A ”linearized” version of this theory produces a power law relation between shear rate and shear stress, $`\sigma \dot{ϵ_{pl}}^{\kappa 1/\kappa +1}`$, where $`\kappa =v_1/v_0`$. The full theory can yield more complicated functional forms with a true yield stress Lemaitre2001 .

#### II.2.2 Soft Glassy Rheology model

The soft glassy rheology (SGR) model of Sollich et al. SGR  is an extension of a trap model for glasses originally proposed by J.-P. Bouchaud Monthus1996  with stress acting as an external drive. It was designed to describe the flow behavior of foams, dense emulsions, pastes and slurries. However, it is very similar to the previous models and should also be relevant to the materials of interest here. Small volume elements are assumed to yield with a rate $`\mathrm{\Gamma }_0\mathrm{exp}[(Ekl^2/2)/x]`$, where $`l`$ is the local strain and $`k`$ an elastic constant. The role of temperature is replaced by a “noise” temperature $`x`$, which is assumed to describe the effect of structural rearrangements in a mean-field spirit. Note that the model describes stress-assisted yielding as in the other two models, but stress enters quadratically and not linearly. Structural disorder is modeled with an exponential distribution of yield energies $`E`$, $`\rho (E)=\mathrm{exp}(E/x_g)/x_g`$, where $`x_g`$ is usually set to one. This introduces a distribution of relaxation timescales. After yielding of a volume element, a new yield energy is drawn from $`\rho (E)`$ and the local strain $`l`$ rises again from 0 according to a macroscopic shear rate $`\dot{\gamma }`$. The time evolution of the probability $`P(l,E;t)`$ describing the ensemble of volume elements can be obtained from a master equation.

Like the STZ theory, the predictions arising from the SGR equations are richer than the simple Eyring model. The exponential distribution of traps induces a dynamical glass transition, and the system exhibits aging for $`x<1`$. Analysis has mainly focused on the steady state situation under constant shear rate $`\dot{\gamma }`$, which is the generic experiment used to determine the mechanical properties of soft glassy materials. Salient predictions are: a Newtonian fluid flow $`\sigma \dot{\gamma }`$ for $`x>2`$ and a power law fluid $`\sigma \dot{\gamma }^{x1}`$ for $`1<x<2`$ SGR . In the glassy phase, a scaling of the form $`\sigma \sigma _y\dot{\gamma }^{1x}`$ is predicted. The nature of the noise temperature $`x`$ (not the true thermodynamic temperature) and the pre-exponential factor $`\mathrm{\Gamma }_0`$ (thermal) remain largely unspecified.

#### II.2.3 Microscopic Approaches

The descriptions of nonlinear rheology and plasticity described so far are appealing because of their simplicity, but much of the physics is put in “by hand”. In the most recent literature, efforts are made to derive the response of driven glassy systems from microscopic considerations. Berthier et al. Berthier2000  consider a generic driven glassy system by calculating the correlation and response functions from the microscopic Langevin equations in a mean-field approximation. Based on this approach, they suggest a “two-time-scale scenario,” in which the slow time scales associated with structural relaxations are accelerated by the drive, while the fast degrees of freedom (phonons) remain at the thermodynamic temperature. This concept can be extended to introduce an “effective temperature” different from the thermodynamic temperature. A relevant conclusion for the present work is that in their calculations for $`T>T_c`$ tc-comm , the slow relaxation time $`t_\alpha `$ decreases with the drive, $`t_\alpha \dot{\gamma }^{2/3}`$. Since the relaxation time determines the viscosity, this leads to power law shear thinning, $`\sigma \dot{\gamma }^n`$ with exponent $`n=1/3`$. Below $`T_c`$, their numerical results indicate that $`n`$ is only very weakly temperature dependent.

## III Molecular Simulations

We perform three-dimensional molecular dynamics (MD) simulations for two model glasses. A bead-spring model is used for polymers. Beads of mass $`m`$ interact via a conventional 6-12 Lennard-Jones (LJ) potential. All results will be expressed in terms of the characteristic length $`a`$, energy $`u_0`$ and time $`\tau _{\mathrm{LJ}}=\sqrt{ma^2/u_0}`$ of this potential. Unless otherwise noted, the potential is truncated at $`r=r_c=1.5a`$ for computational convenience. We construct linear polymer chains by connecting adjacent beads with the finite extensible nonlinear elastic (FENE) bond potential Kremer1990 . This polymer model has been used extensively to study polymer melt dynamics Puetz2000  and was also used in our previous study of yield conditions Rottler2001 . The number of beads per chain used below is usually 256, but the chain length and entanglement effects are not important for the small strains considered in Section IV.1.

In addition to the polymer, a binary mixture composed of 80% A-particles and 20% B-particles without covalent bonds is also studied. LJ interaction parameters were set to the values employed in previous studies that aimed at verifying predictions of mode-coupling theory for supercooled liquids Kob1995 , and studied aging Kob1997  or dynamical heterogeneities Kob1997b  during the glass transition. A very similar system with slightly different parameters was used by Falk STZ  to study deformation and plasticity in amorphous metals in two dimensions.

Both the polymer and binary mixture models enter an amorphous glassy state without crystallization upon cooling. For the polymer model, the glass transition temperature $`T_g0.35\pm 0.05u_0/k_B`$, while $`T_g`$ is smaller for the binary system. These values are for $`r_c=1.5a`$, and slightly higher values of $`T_g`$ are obtained with larger $`r_c`$ Bennemann1998 . Here, we focus on a temperature range between $`T=0.01u_0/k_B`$ and $`T=0.3u_0/k_B`$. Unless noted, the temperature is controlled with a Langevin thermostat (damping rate $`1\tau _{\mathrm{LJ}}^1`$). The equations of motion are solved using the velocity Verlet algorithm with a timestep of $`dt=0.0075\tau _{\mathrm{LJ}}`$.

The simulation cell contained 32768 LJ beads. Previous studies of the yield behavior Rottler2001  have shown that an increase of the system size beyond this size leads to slightly lower values of the shear yield stress, but the generic behavior is unchanged. In order to minimize statistical fluctuations while varying the strain rate, we use the same initial state for all runs at a given temperature. The glassy states were prepared by a quench from a fluid temperature of $`T=1.3u_0/k_B`$ to a glassy temperature of $`T=0.3u_0/k_B`$ at constant volume over a time interval of order 1000$`\tau _{\mathrm{LJ}}`$. The density was chosen so that the hydrostatic pressure was zero at this temperature. Lower temperatures were then reached by cooling at constant pressure.

In studies of initial yield behavior, periodic boundary conditions are applied in all directions to eliminate edge effects. The original cell is cubic and at zero hydrostatic pressure. Tensile or compressive strains are imposed on one or more axes by rescaling the periods $`L_i`$. We use true strain rates $`\dot{ϵ}_i=L_i^1dL_i/dt`$ between $`\dot{ϵ}_i=10^6\tau _{\mathrm{LJ}}^1`$ and $`\dot{ϵ}_i=10^3\tau _{\mathrm{LJ}}^1`$. Note that since $`\tau _{\mathrm{LJ}}3\mathrm{ps}`$, the strain rates employed are much higher than typical experimental strain rates. As discussed below, sound propagation is too slow for stress to equilibrate across the system at the highest strain rates. However, most simulations are slow enough that loading proceeds nearly quasistatically.

Steady state flow cannot be investigated by deforming the simulation box in the above manner, because some box periods would soon decrease to a single molecular diameter. One approach is to use Lees-Edwards boundary conditions Allen1987 . Here, we choose a different route and replace the periodic boundary conditions in the 3-direction with two rigid walls composed of 2 layers of an fcc(111) crystal for simulations of simple shear. Otherwise, the simulation cell has the same dimensions and size, and the wall beads are strongly coupled to the sheared glass so that slip at the interface is prohibited. By moving one wall parallel to the solid at constant velocity $`v_x`$, a steady state shear profile is imposed, and we can measure the shear stress $`\tau `$ as a function of the average strain $`\gamma =tv_x/h`$, where $`h=32a`$ is the wall separation and $`t`$ is the elapsed time. In these simulations, the Langevin thermostat is only coupled to the irrelevant (y) direction perpedicular to the flow (x) and velocity gradient (z) directions Thomson1990 .

## IV Results

### IV.1 Onset of shear

In ref. Rottler2001  we analyzed the onset of shear deformation in polymer glasses. Results for a wide range of multiaxial loading conditions, temperatures and potential parameters were consistent with the pressure-modified von Mises criterion. In this criterion, the driving force for shear is the deviatoric stress $`\tau _{\mathrm{dev}}=\left((\sigma _1\sigma _2)^2+(\sigma _2\sigma _3)^2+(\sigma _3\sigma _1)^2\right)^{1/2}/3`$, where the $`\sigma _i`$ are the three eigenvalues of the stress tensor. Rather than having a sudden onset at a finite strain, irrecoverable deformation was observed at arbitrarily small strains. The most robust definition of the yield stress was $`\tau _{\mathrm{dev}}^y`$, the peak value of the deviatoric stress as a function of strain. This quantity increases linearly with pressure in agreement with the pressure-modified von Mises criterion.

The results in ref. Rottler2001  were obtained for a single value of the strain rate. Here we focus on the effect of varying strain rate on the shear yield stress. The initially cubic simulation cell is expanded in one direction at a constant strain rate $`\dot{ϵ}_z=L_z^1dL_z/dt`$, and volume $`V`$ is conserved by maintaining $`L_x=L_y=\sqrt{V/L_z}`$. Fig. 1 shows the deviatoric stress averaged over the entire simulation cell as a function of strain at three different strain rates. As can be seen, the behavior of the bead spring model (a) and the binary LJ glass (b) is qualitatively similar. These curves also closely resemble experimental stress-strain curves for e. g. polymeric glasses Hasan1995 .

For all systems the initial response is nearly elastic, i.e. the stress rises almost linearly with strain. In the quasistatic limit, the initial slope gives the elastic modulus. Results for strain rates of $`3\times 10^4\tau _{\mathrm{LJ}}^1`$ and below collapse onto this quasistatic behavior. As the strain rate rises to $`10^3\tau _{\mathrm{LJ}}^1`$ and beyond, the initial slope grows. The reason is that stress is no longer able to equilibrate across the system. It is well known that the elastic modulus of a heterogeneous system is over-estimated by applying a uniform strain (the Voight limit Voight1928 ). Our algorithm imposes a uniform strain at each step by rescaling the cell dimensions and particle coordinates and the resulting stress will be too high when the strain rate becomes too large. A rough estimate for the characteristic time for stress equilibration through the system is $`L/c10\tau _{\mathrm{LJ}}`$, where $`c`$ is the speed of sound. When the strain is stopped suddenly we find an exponential stress relaxation with a characteristic time $`15\tau _{\mathrm{LJ}}`$. At a strain rate of $`\dot{ϵ_i}=10^3\tau _{\mathrm{LJ}}^1`$, the time to reach a strain of 2% is $`20\tau _{\mathrm{LJ}}`$ and becomes comparable to the above estimates. We conclude from Fig. 1 that the shear rate is slow enough to produce a quasistatic deformation for $`\dot{ϵ_i}<10^3\tau _{\mathrm{LJ}}^1`$, but not at the highest shear rates.

As the strain increases, even the curves for lower strain rates split apart and saturate at different maximum heights $`\tau _{\mathrm{dev}}^y`$. The maxima are broad and centered at strains of about 6% in the binary LJ glass and 8% in the polymer. The fine structure on the curves corresponds to individual plastic yield events that are discussed further below. This structure decreases with increasing system size and temperature as the fraction of the system involved with typical events decreases. Results for larger systems were consistent with those shown here, but could not be extended over as wide a range of shear rates and other parameters.

Fig. 2 summarizes values for $`\tau _{\mathrm{dev}}^y`$ obtained from the maximum of curves such as those in Fig. 1 as a function of strain rate. Four different temperatures were studied, with the highest value, $`T=0.3u_0/k_B`$, close to $`T_g`$ and $`T=0.01u_0/k_B`$ far away. This covers a much wider range than usually explored in experiments. For a given temperature and most strain rates, both models give nearly straight lines in a semilogarithmic plot. The salient observation that can be made in this figure is that all curves are nearly parallel, and temperature merely changes the offset value on the stress axis.

As suggested by the Eyring-model, an obvious way to describe the rate dependence is through a logarithmic fit of the form $`\tau _{\mathrm{dev}}^y=\tau _0+s\mathrm{ln}\dot{ϵ}`$ (solid lines). The observation of parallel curves then implies a nearly constant prefactor $`s`$ in front of the logarithm, and the variation with temperature is described by $`\tau _0`$. Often the rate dependence of the shear yield stress of complex materials is also described with a power-law added to a constant, i. e. $`\tau _{\mathrm{dev}}^y=\tau _0+r\dot{ϵ}^n`$ (dashed lines). As can be seen in Fig. 2, such fits provide an equally good description of the data. Due to the small variations of $`\tau _{\mathrm{dev}}^y`$, a determination of the exponent $`n`$ via best fits is very unreliable. Since the curves for different $`T`$ are nearly parallel, there is no reason to expect $`n`$ to vary. The curves in Fig. 2 are for $`n=0.2`$, which gave the smallest variation of the prefactor $`r`$ ($`<10\%`$) with temperature, but other choices between 0.1 and 0.3 are also possible. Since $`\tau _0>0`$ for $`T0.2u_0/k_B`$, the data is not consistent with a pure power law.

Neither functional form provides a good fit to results at the highest strain rate $`10^3\tau _{\mathrm{LJ}}^1`$. We have already shown that the elastic behavior changes at this shear rate because stress can not equilibrate throughout the system. The plastic response will also be affected by the dynamics of stress distribution, and new behavior is expected to set in at this rate. Therefore, these points were not included in the fits shown in Fig. 2. Note that experiments are always at much lower strain rates, while most simulations of shear have been done at $`\dot{ϵ_i}=10^5`$ to $`10^1\tau _{\mathrm{LJ}}^1`$. A clear distinction between a logarithm and a power-law dependence would thus only be possible by measuring the yield stress at lower rates, but unfortunately the simulation time becomes prohibitively large.

The fact that the curves of Fig. 2 are parallel implies that replotting the shear yield stress as a function of temperature also leads to parallel behavior. Fig. 3 shows that the trend of $`\tau _{\mathrm{dev}}^y`$ with $`T`$ is in fact linear, and the slope is nearly rate independent. The rate merely changes the offset value.

In order to investigate whether the rate dependence is in any way influenced by the methodology, we have expanded our studies and considered other model parameters, different thermostat methods and loading states. In particular, we obtained data analogous to that shown in Fig. 2 for a longer range of the Lennard-Jones potential $`r_c=2.2`$a, Nosé-Hoover and Langevin thermostats with different rates, and uniaxial strain as opposed to the volume-conserving shear in Fig. 1. We also considered a system that was quenched into the glassy state 10 times faster than the other systems to investigate cooling rate effects. All these different situations show essentially the same robust scenario: the rate dependence is nearly independent of temperature with a constant offset that decreases linearly as $`T`$ increases. The slope $`s`$ is generally unaffected within the noise, but is roughly twice as large for uniaxial strain, and slightly larger for rapidly quenched initial states. Again, power-law fits could also be used to describe the data, and a fixed exponent suffices to describe all temperatures.

### IV.2 Steady shear

Most rheological measurements of soft glassy materials focus on the steady state stress and not on the initial transient maximum of the stress-strain curve. Thus it is interesting to compare the rate dependence for steady state shear to the above results. The steady state can not be studied with the entangled bead-spring model used above <sup>2</sup><sup>2</sup>2The bead-spring model could be used when the chain length is very short., because the shear stresses at accessible strain rates would soon break covalent bonds. However previous studies with short chains Baljon  and the binary system Yamamoto1998 ; Lacks2001 ; Barrat2000  show similar trends with decreasing $`T`$. At high $`T`$ ($`0.7u_0/k_B`$ for the binary system), Newtonian behavior is observed up to the highest practical shear rates. At lower temperatures, where the liquid is in a supercooled state, there is a crossover from Newtonian behavior at low rates to power law shear thinning $`\sigma \dot{\gamma }^n`$ at high rates. The crossover rate goes to zero as $`T`$ decreases to $`T_g`$, and at lower temperatures the shear stress is nearly independent of strain rate Baljon .

We focus here on the binary (80/20) LJ system at $`T0.3u_0/k_B`$, varying $`T`$ from just above $`T_g`$ to well below. As described in Section III, shear was imposed by confining the system between walls separated by $`h`$ along the z-direction and translating one at a fixed speed $`v_x`$. Fig. 4 shows typical behavior of the stress as a function of the average shear strain $`\gamma `$. There is an initial maximum at $`\gamma 0.1`$ that corresponds to the transient behavior studied in the previous section. The stress then decreases rapidly to a steady state plateau value. This strain softening becomes more pronounced with decreasing temperature.

The average value from the steady state region is plotted vs. shear rate in Fig. 5. At $`T=0.3u_0/k_B`$ (lowest curve), the solid exhibits glassy behavior at high rates, but is Newtonian at lower rates, which indicates that the temperature is still above $`T_g`$. For all temperatures below $`T_g`$ the curves are nearly parallel, as for the results for the transient maximum shear stress in Fig. 2. A logarithmic rate dependence may be fitted at small rates, but a power-law added to a constant clearly provides a much better fit over the whole range of rates. The fits for the temperatures $`T0.2u_0/k_B`$ shown in Fig. 5 use $`n=0.3`$, for which the prefactor $`r`$ to the power-law deviates from unity by less than 10%. Our results for the glassy state are very similar to observations reported by Varnik et al. Varnik2002 , who used the same model with a larger cutoff value $`r_c=2.5a`$ for the LJ potential and only considered $`T=0.2u_0/k_B`$.

The rate dependence of the steady shear stress is very similar to that of the yield stress (Fig. 2). In both cases, there is a rapid rise above the logarithmic fits at rates of $`10^3\tau _{\mathrm{LJ}}`$ and above. Although in Fig. 5 the system has had time to reach a steady state, one may wonder whether there is still a change in behavior near $`\dot{\gamma }=10^3\tau _{\mathrm{LJ}}^1`$. At higher rates the stress may not relax between the local yield events discussed below, leading to a more rapid rise in mean stress. The results were not extended to higher shear rates because the temperature in unthermostatted directions begins to rise slightly above the set temperature even at $`\dot{\gamma }=10^2\tau _{\mathrm{LJ}}^1`$.

A plot of the shear stress vs. rate as in Fig. 5 assumes that the local strain rate is equal to the average implied by $`v_x/h`$. Fig. 6 shows the velocity profile for two different temperatures and several different shear rates. The first few layers always move with the walls due to the strong coupling. At $`T=0.2u_0/k_B`$, the profile in the center region is homogeneous and here the assumption that $`v_s/h`$ equals the shear rate is satisfied. However, at $`T=0.01u_0/k_B`$ shear is localized in the lower 60% of the simulation cell, and the upper part of the glass moves at the constant wall speed. This is a clear indication of shear banding. The local shear rate is now larger than the average value, but since the change is only by about a factor of two, the data points in Fig. 4 are not dramatically affected.

In their simulations at $`T=0.2u_0/k_B`$, Varnik et al. Varnik2002  also found a transition from homogeneous flow to shear banding as the shear rate dropped below $`10^3\tau _{\mathrm{LJ}}^1`$. In our simulations, shear is homogeneous at this temperature, but the shorter cutoff in our model implies a slightly lower $`T_g`$. The increasing amount of shear localization with decreasing temperature is consistent with a stress peak that becomes more pronounced as T is lowered (see Fig. 4). The region with negative slope on the stress-strain curve signals a mechanical instability and shear localization. This localization is inhibited in the simulations that were used to study the transient stress maximum, since they used periodic boundary conditions in all directions. Varnik et al. Varnik2002  noted that simulations with Lees-Edwards boundary conditions Berthier2002  also suppressed shear banding, but that the shear stresses obtained from simulations with walls and periodic boundary conditions are similar.

### IV.3 Analysis of the dynamics of the local stress distribution

One of the great advantages of molecular simulations is the ability to reach beyond a measurement of the macroscopic response function and in addition obtain information about the local dynamics and rearrangements in nonequilibrium situations. Such information is essential to construct a clear picture of the underlying microscopic processes.

Here, we follow this route by decomposing the total simulation cell into small volume elements of size $`78a^3`$ and measuring the local stress tensor in those regions. Since the density is close to $`1a^3`$, these regions contain typically 7-8 particles. A particle number of that size should be sufficient to constitute a locally transforming region as envisioned in several of the theoretical models described in Section II. Some of the models suggested that local rearrangements replace thermal noise as the source for activation over barriers, which motivates a study of the stress changes experienced by local regions.

Our results for general stress states Rottler2001  have shown that the deviatoric shear stress $`\tau _{\mathrm{dev}}`$ is the relevant stress tensor invariant that describes shear deformation. We therefore calculate the local change $`\mathrm{\Delta }\sigma _{ij}`$ in the stress tensor in every volume element during a small time interval and then study the size distribution of changes in the scalar variable $`\mathrm{\Delta }\tau _{\mathrm{dev}}`$ calculated from $`\mathrm{\Delta }\sigma _{ij}`$ (note that $`\mathrm{\Delta }\tau _{\mathrm{dev}}`$ is always positive, since it refers to the deviatioric stress of $`\mathrm{\Delta }\sigma _{ij}`$) . In general, one expects $`\mathrm{\Delta }\tau _{\mathrm{dev}}`$ to fluctuate even in the undriven case, and one should expect to find a stationary distribution of stress jumps over not too long timescales. Since the glass is out of equilibrium, it is by construction not stationary and will exhibit aging phenomena, etc. Such long timescales, however, are not explored here.

Fig. 7 shows an example of such stationary distributions for four different temperatures. The jumps correspond to an average change in local deviatoric stress in small regions over a time difference of $`7.5\tau _{\mathrm{LJ}}`$. The narrowest distribution is found at the lowest temperature $`T=0.01u_0/k_B`$ and the distribution widens as the temperature increases to $`T=0.3u_0/k_B`$. The cause for the stress jumps is obviously thermal motion of the particles, and only small excursions about their positions can occur in the glass (cage effect). The part of the distribution at small values of $`\mathrm{\Delta }\tau _{\mathrm{dev}}`$ can be fitted to a Gaussian, but deviations become visible at large values.

We are now able to determine how the distribution of stress jumps in the undriven system changes when the system is macroscopically strained. Fig. 8 shows several distributions at the highest ($`T=0.3u_0/k_B`$) and lowest ($`T=0.01u_0/k_B`$) temperatures considered. $`\mathrm{\Delta }\tau _{\mathrm{dev}}`$ was calculated for the same time difference $`7.5\tau _{\mathrm{LJ}}`$. For each temperature, four distributions are shown, which correspond to four different strains between zero and the maximum shear yield strain. Also shown for comparison are the distributions in steady shear.

Several observations can be made in this figure. Note first that at $`T=0.3u_0/k_B`$, the unstrained distribution is not changed dramatically under external driving. Only the tail at large jumps gets modified. This temperature is above $`T_g`$ and the shear rate is low enough that the deviation from Newtonian behavior is small. This is in sharp contrast to the situation at $`T=0.01u_0/k_B`$, where strain produces a dramatic increase in the number of large jumps. The tail of the distribution can be fit to an exponential form $`\mathrm{exp}(\mathrm{\Delta }\tau _{\mathrm{dev}}/\tau _c)`$ that extends to larger stresses as the strain increases to the peak in the stress-strain curve. The steady state stress results are close to those at a transient strain of 4.5%. The steady state probability curves decay with the characteristic stresses $`\tau _c=0.18u_0/a^3`$ $`(T=0.01u_0/k_B)`$ and $`\tau _c=0.12u_0/a^3`$ $`(T=0.3u_0/k_B)`$.

## V Interpretation and comparison to models of plasticity

The above results for the strain rate dependence of shear yielding offer an opportunity to test the theoretical models described in Section II. The Eyring model Eq. (2) predicts a logarithmic dependence of $`\tau _{\mathrm{dev}}^y`$, and the prefactor $`s`$ is given by $`k_BT/V^{}`$. Since we found typical values of $`s`$ between 0.02 and 0.03 for all temperatures, our result can only be reconciled with the Eyring model if one allows for huge variations of the “activation volume” $`V^{}`$ between $`0.3a^3`$ and $`10a^3`$. $`V^{}`$ is a phenomenological fit parameter, but is typically interpreted as a characteristic volume for a local shear event. It should then be at least of order the volume per particle, i.e. of order $`a^3`$ or larger. The implied linear change in $`V^{}`$ with $`T`$ would also be inconsistent with the observed linear temperature dependence of $`\tau _{\mathrm{dev}}^y`$. From Eq. (2) the stress would vary as $`E/V^{}E/T`$, which is inconsistent with the data. Besides, the usefulness of the Eyring expression is reduced if $`V^{}`$ is allowed to be temperature dependent.

Modern theories offer a different interpretation of the rate dependence. Lemaître’s extension Lemaitre2001  of the STZ theory STZ  to include free volume relaxation produces complex rate dependence that could describe our results. However, the issue of temperature dependence has not been addressed in this approach and we have not examined the issue of free volume relaxation here.

A pure power-law shear thinning behavior was predicted in the calculations of Berthier et al. Berthier2000 . Starting from very different considerations, these authors found that the exponent $`n`$ should have a generic value of 1/3 at a transition temperature $`T_c`$. Below $`T_c`$, the exponent should in principle be temperature dependent, but their numerical results indicated only a very weak dependence. Both results are consistent with our findings near $`T_g`$, but not at lower temperatures where we find a constant plus power-law behavior.

Below the glass transition temperature, the SGR model SGR  predicts such constant plus power-law behavior of the flow curve with an exponent of the form $`1x`$, where $`x`$ is an effective noise temperature in units where the glass transition temperature $`x_g=1`$. In steady shear, we found good agreement with this functional form and an exponent $`n=0.3\pm 0.1`$. The authors of the SGR model go to great length in providing an interpretation for $`x`$. They argue that in a sheared state, the energy for surmounting a barrier for rearrangement is not only provided by the thermal energy, but also by the energy released from rearrangements elsewhere in the material. This energy diffuses through the material and provides an effective thermal bath in a mean field sense. The energy released in such a rearrangement must therefore be of order of the typical yield energies, which implies that $`x`$ should be of order unity. When the yield energies are much larger than typical thermal energies, $`x`$ will be independent of the true thermodynamic temperature $`T`$. Our finding of an exponent independent of $`T`$ could thus be rationalized in this framework.

Our analysis of the distribution of local stress jumps could lend additional support to the concepts behind the SGR model. Changes in the local deviatoric stress as shown in Fig. 8 can activate yield events. We found that the distribution of stress jumps could be fitted to an exponential distribution with a characteristic decay stress that varied less than 30% as $`T`$ changed from $`0.3u_0/k_B`$ to $`0.01u_0/k_B`$. At small values of $`\mathrm{\Delta }\tau _{dev}`$, the distribution retains the equilibrium form. This result suggests that the drive generates internal dynamics on a common scale despite very different thermodynamic temperatures and might provide a more microscopic justification of the background “effective noise temperature” proposed in the SGR model.

The shear energy released by a yielding local region should indeed trigger additional yield events at other locations in the material. However, the ensuing dynamics may be more complicated than suggested by mean-field approaches. Many models of material breakdown include load redistribution mechanisms. In fiber bundle models Kloster1997 , for instance, one finds avalanche behavior that precedes total failure. The avalanche size distribution follows a power law. When shearing a material, however, a yielded region will typically reemerge with a different yield energy as assumed in the SGR model. In ref. Zapperi1997 , the simpler situation of an elastic network with random breaking thresholds of the links was subjected to an external drive. Damaged elements were replaced with new ones with a different breaking threshold. Power-law behavior in the size and duration of failure events was found. Recent work on molecular Gagnon2001  as well as simplified models Baret2002  shows that load redistribution in the yielding glass mediated by elastic interactions between rearranging regions can lead to avalanche behavior.

## VI Conclusions

The (transient) maximum deviatoric shear stresses for two model glasses, binary LJ mixtures and polymers, were shown to exhibit strikingly similar trends with rate and temperature. The rate dependence was remarkably constant from $`T_g`$ to temperatures 30 times lower. The entire flow curve shifted linearly to lower stress as $`T`$ increased. At low rates, the rate dependence of the peak stress could be described by a logarithm or constant plus power-law, where the prefactor and/or exponent did not vary with temperature. A more rapid rise in stress was observed when strain rate was increased to $`10^3\tau _{\mathrm{LJ}}^1`$. Deviations in the elastic response also set in at this strain rate, indicating that stress could not relax throughout the system. This is not surprising given that the time to strain to 1% is comparable to that of sound propagation back and forth across the system.

The stress for steady shear flow was calculated for the LJ mixture. Curves for different temperatures were also nearly parallel, shifting rigidly to lower stresses with increasing temperature. Below strain rates of $`10^3\tau _{\mathrm{LJ}}`$, the dependence of the flow stress on rate could be described by a logarithm. The entire flow curve could be fit by a constant plus power law with a temperature independent exponent $`n=0.3\pm 0.1`$. However, as for the transient stress, the stress begins to rise more rapidly at $`10^3\tau _{\mathrm{LJ}}^1`$. It is interesting to ask whether the emergence of a power-law vs. logarithmic behavior is related to the overlap of shear rate and stress equilibration timescales. Although the system is in steady state, regions of a few atoms undergo a substantial yield event after of order 1% strain increments. At strain rates of $`10^3\tau _{\mathrm{LJ}}^1`$ and above, the time between these yield events will be too short for stress relaxation in our system. This point should be kept in mind when comparing molecular simulations to experiments at lower strain rates.

The Eyring model was shown to be incompatible even with logarithmic fits to the shear stress at low rates. It predicts that the prefactor of the logarithm should scale as $`kT/V^{}`$, while the observed prefactor is nearly independent of temperature. Fixing this discrepancy by requiring $`V^{}`$ to scale linearly with $`T`$ leads to unphysically small values of $`V^{}`$ and is inconsistent with the linear drop in yield stress with temperature at a given rate. The Eyring model has been very helpful in analyzing experimental data, but typically over a narrow range of temperatures. It would be interesting to extend these measurements to very low temperatures and to higher shear rates. Note that typical room temperature experimental values for $`V^{}`$ in polymers correspond to 3 or 4 repeat units Ward1983 , which is of the same order as the values we find at $`T=0.1`$ or $`0.2u_0/k_B`$.

The insensitivity of the rate dependence to temperature changes indicates that thermal activation is not dominant at the rates studied. Analysis of the individual particle trajectories offers an explanation for this. Even in the limit of zero temperature we find that the trajectories are exponentially sensitive to changes in rate and other parameters. The system does not travel along a single path through the energy landscape at different rates, but gets deflected between many possible paths by small perturbations.

Both the STZ and SGR models describe complex dynamics in systems where temperature is unimportant. Instead, activation is due to another internal state variable, either $`v_f`$ or $`x`$, that couples to the external drive via a feedback mechanism. The original STZ theory gives a simple linear rise in stress with rate above the yield stress STZ . However recent generalizations Lemaitre2001  can produce more complex rate dependence like that found here. The SGR model predicts a constant plus power-law behavior that is also consistent with our simulations and can account for a temperature independent exponent $`n`$. It remains to be seen whether the generating mechanism of structural disorder (SGR) or free volume dynamics (STZ) provides a more useful description of glassy dynamics.

An external drive introduces a new timescale into the problem that couples to the structural rearrangements. It will therefore modify the associated relaxation time scales and alter the unperturbed glassy dynamics, e.g to stop aging and induce rejuvenation Liu2001 . At zero temperature, this timescale should be the only one present in the system, while other competing timescales will arise at finite temperature. Associated with this timescale is the concept of an effective temperature Barrat2000 ; Ono2001  that also appears in the SGR model. While this picture is clearly a useful starting point, the analysis of the jump distribution has shown that the internal dynamics of the flowing glassy solid is much more intricate and strongly calls for extensions of analytical models beyond mean-field concepts. Such a treatment should describe more accurately the load redistribution mechanisms and the propagation of released shear energy on a microscopic level. A preliminary analysis has shown that local plastic events occur in avalanches over a wide range of length and time scales.

Finally, we note that the above analytic models of glassy rheology and viscoplasticity are all scalar models that cannot describe self-organization on larger length scales such as shear bands. Tensorial versions of the STZ theory are being discussed to address banding and necking Langer2001 . These are promising approaches that should continue to benefit from insight gained from molecular simulations, in particular the load redistribution mechanisms of shear yielding zones.

## VII Acknowledgements

We thank J. S. Langer and M. L. Falk for useful discussions. 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:
# A grid of synthetic ionizing spectra for very hot compact stars from NLTE model atmospheres

## 1 Introduction

Post-AGB stars display the hottest stage of stellar evolution just before they enter the white dwarf cooling sequence. Effective temperatures up to 180 kK have been determined (e.g. Werner et al. 1997) by means of spectral analyses based on NLTE (non local thermodynamic equilibrium) model atmosphere techniques.

A sufficient condition for the use of LTE models in spectral analyses is, that the collisional rates between atomic states are much larger than the corresponding radiative rates (cf. Kudritzki & Hummer 1990). This is e.g. fulfilled in the case of “cool” ($`T_{\mathrm{eff}}`$$`\stackrel{<}{}\mathrm{\hspace{0.17em}40}\mathrm{kK}`$) DA white dwarfs (Napiwotzki 1997). However, there are always NLTE effects in any star which are easily detectable if high-resolution high-S/N spectra, especially of the high energy wavelength range (UV, EUV, X-ray), are evaluated (e.g. Rauch 1997 and references therein) and these effects are particularly important for hot stars like e.g. central stars of planetary nebulae (CSPN).

Observations in the UV – X-ray wavelength range have shown that this region is generally dominated by metal opacities. Since the flux maxima of hot stars lie in the EUV/X-ray wavelength (Fig. 1), and the flux in this region is crucial to the ionization of ambient matter, e.g. a surrounding planetary nebulae (PN), the employment of realistic fluxes calculated from fully line-blanketed NLTE model atmospheres as ionizing spectra in photoionization models is very important. However, blackbody fluxes are still often in use.

In the early 80’s of the last century, the calculation of NLTE models was rudimental and hampered by the existing insufficient numerical methods and computational capacity (cf. Rauch 1997). Following Pottasch et al. (1978), who has shown that, based on intermediate band photometry within 1500 - 3300 Å, some CSPN have their radiation well represented by blackbodies (and others not so well), Kaler (1983) addressed that the major problem in PN analyses is the model selection. He invoked then that a blackbody has at least the virtue of simplicity.

The implementation of the “accelerated lambda iteration” (ALI) technique in the NLTE model-atmosphere codes (Werner 1986) provided an efficient method to calculate realistic synthetic spectra of hot stars. However, at that time it had been impossible to calculate extensive grids of fully line-blanketed NLTE models due to the lack of reliable atomic data (cf. Rauch 2003). Thus, Kwok (1994) repeated Kaler’s (1983) statement, that due to the variety of CSPN parameters, the blackbody assumption is probably not worse than any other specific set of models. In order to demonstrate the strong deviations at energies higher than 13.6 eV and 54.4 eV (I and He II ground state thresholds, respectively), we show in Fig. 1 a comparison of blackbody fluxes to NLTE model atmosphere fluxes. The weaknesses of a blackbody are obvious. Even in a H+He model with $`T_{\mathrm{eff}}=110\mathrm{k}\mathrm{K}`$, $`\mathrm{log}g=7`$, and a solar He/H ratio, the model flux at energies higher and lower than the He II ground state cannot be matched simultaneously with one blackbody flux. Instead, blackbody fluxes with 90 kK and 130 kK are needed, respectively. At higher energies ($`\lambda <160\mathrm{\AA }`$), the situation gets even worse.

NLTE model atmosphere fluxes have been successfully used, e.g. by Rauch et al. (1994, 1996) for the construction of consistent models of PN and CSPN in the cases of the PNe K 1-27 and LoTr 4. This approach — consistent models of PN and CSPN — appears to be a promising way to improve PN analyses.

Recently, Armsdorfer et al. (2002) selected the spatially resolved, round PN NGC 2438 as a test example to check for effects of CSPN models on PN shell modeling. For their study they used the grid parameters $`T_{\mathrm{eff}}=110\mathrm{k}\mathrm{K}`$ and $`\mathrm{log}g=7`$ which are closest to the result of Rauch et al. (1999, $`T_{\mathrm{eff}}=114\pm 10\mathrm{k}\mathrm{K}`$, $`\mathrm{log}g=6.62\pm 0.22`$) from spectral analysis of the CSPN of NGC 2438. They found that the use of a blackbody as ionizing source leads e.g. to an underestimation of the \[O III$`\lambda \mathrm{\hspace{0.17em}5007}`$Å line and an overestimation of the He II $`\lambda \mathrm{\hspace{0.17em}4686}`$Å line strength and, hence, an overestimation of the CSPN temperature, compared to the use of a model from the H-Ca grid (Sect. 1).

Rauch (1997) has investigated in detail on the impact of light metals (Li – Ca) on NLTE model atmospheres of compact hot stars (see this paper for a more detailed description of advances in NLTE calculations and the consideration of metal opacities). He found that CNO opacities reduce the flux strongly at energies higher than the oxygen absorption edges (Fig. 1) while the opacities of light metal reduce the flux even more, but at higher energies. Rauch concluded that the metal opacities have a great influence on the structure of model atmospheres. For the calculation of realistic emergent fluxes in the EUV and X-ray range, the consideration of metals is indispensable. Consequently, a first grid of models which are composed out of all elements from hydrogen to calcium had been calculated with solar (\[X\]111\[ \]: log (abundance / solar abundance) = \[Y\] = \[Z\] = 0) and halo (\[X\] = \[Y\] = 0, \[Z\] = $`1`$) abundance ratios. This allows to interpolate or extrapolate slightly in the metal abundances \[Z\]. The emergent fluxes are available in tabular form at http://astro.uni-tuebingen.de/$``$rauch.

The influence of the iron-group elements was demonstrated e.g. by Dreizler & Werner (1993), Lanz & Hubeny (1995), and Deetjen et al. (1999). Thus, parallel to the calculation of the H-Ca grid, test calculations have been performed in order to check the impact of the iron-group opacities (Rauch 2000, Rauch & Deetjen 2001, Rauch 2002). The consideration has an influence on the atmospheric structure and thus, also on the strengths of the absorption edges of light metals. Since significant changes in the flux level in the UV – X-ray region occur (Fig. 2), the iron-group opacities should not be neglected. Thus, a second grid with all elements from hydrogen up to the iron group has been calculated. This is described in the following.

## 2 NLTE model atmospheres

Our NLTE models (Werner 1986, 1989; Werner & Dreizler 1999, Werner et al. 2003) are plane-parallel and in hydrostatic and radiative equilibrium. Hence, our models neglect wind effect onto the spectrum, which can be significant close to the Eddington limit (Gabler et al. 1989).

All elements from hydrogen to the iron-group elements can be considered simultaneously. More than 300 levels (limited by numerical accuracy in the used code-version) can be treated in NLTE with more than 1 000 line transitions calculated in detail. Millions of lines of the iron-group elements tabulated in Kurucz (1996) and data from the Opacity Project (Seaton et al. 1994) are accounted for using a statistical approach and an opacity sampling method (Deetjen 1999, http://astro.uni-tuebingen.de/$``$deetjen).

For the construction of the model atoms used in our calculations, we follow Rauch (1997). H, He, C, N, O, and Ne are represented by small model atoms, the other species are modeled to fit into the capacity limit of the code. As a consequence, many of the ions of F, Na - K, (Table 1) are too small to include radiative transitions, not even the resonance lines. We performed some test calculations with larger model ions for selected atoms in order to include at least all resonance lines of these ions. We found that these lines are very weak, so they do not have any significant influence on the atmospheric structure. Since our aim is not a precise spectral analysis of a given star but the calculation of emergent fluxes as ionizing spectra, this is a good approximation.

In contrast to Rauch (1997), Ca is not treated as a single species but combined with the iron-group elements into a generic model atom in order to consider all Ca lines given by Kurucz (1996). Within the generic model atom, Ca - Ni are represented with solar abundance ratios. Depending of $`T_{\mathrm{eff}}`$, the model atoms have individually been adjusted due to the dominant ionization stages. The statistics of representative model atoms are summarized in Table 1. For a more detailed description of the used atomic data and the consideration of the iron-group elements, see Rauch & Deetjen (2003).

We calculated a grid of models within $`T_{\mathrm{eff}}=50\mathrm{1\hspace{0.17em}000}\mathrm{kK}`$ and $`\mathrm{log}g=59`$ (Fig. 3), with solar and halo abundances (Fig. 4). We use $`T_{\mathrm{eff}}^{\mathrm{mod}}=\left(4\pi /\sigma _0^{\mathrm{}}_\nu 𝑑\nu \right)^{1/4}`$ for the convergence criterion ($`_\nu `$ is the Eddington flux). If this is closer than 2% to the target $`T_{\mathrm{eff}}`$, the model calculation is stopped in order to stay within a reasonable cpu time on the CRAY computers. We mention that it may take about one day to calculate a single model.

Flux tables have subsequently been calculated within 5-2000Å binned to 0.1Å intervals. $`T_{\mathrm{eff}}^{\mathrm{mod}}`$ as well as the model’s $`T_{\mathrm{eff}}`$ and $`\mathrm{log}g`$ are given in the tables’ headers (Table 2). These tables are directly suited to be used with the photoionization code CLOUDY (Ferland 2000, http://www.pa.uky.edu/$``$gary/cloudy/, versions C96 Beta 4 and later). The complete flux tables, based on the frequency grids used in the model calculations, are also available in case that a normalization to a certain (e.g. V) magnitude is required.

In order to judge the necessity of a detailed consideration of iron-group opacities in calculations of ionizing fluxes which shall be used in photoionization models of a PN, we measured the changes of the strengths of the I and He II absorption edges in the H-Ca and H-Ni fluxes. We calculated

$$s=\frac{_{\lambda _0\mathrm{\Delta }\lambda }^{\lambda _0}F_\lambda ^{\mathrm{H}\mathrm{Ni}}d\lambda }{_{\lambda _0}^{\lambda _0+\mathrm{\Delta }\lambda }F_\lambda ^{\mathrm{H}\mathrm{Ni}}d\lambda }/\frac{_{\lambda _0\mathrm{\Delta }\lambda }^{\lambda _0}F_\lambda ^{\mathrm{H}\mathrm{Ca}}d\lambda }{_{\lambda _0}^{\lambda _0+\mathrm{\Delta }\lambda }F_\lambda ^{\mathrm{H}\mathrm{Ca}}d\lambda }$$
(1)
with $`\lambda _0=911.5\mathrm{\AA }`$ and $`\lambda _0=227.8\mathrm{\AA }`$ for I and He II, respectively, and $`\mathrm{\Delta }\lambda =50\mathrm{\AA }`$. While the changes of the I absorption edge are generally less than 10%, the strengths of the He II absorption edge may change by a factor of about $`23`$ at lower temperatures. Depending on $`\mathrm{log}g`$, there is a maximum of $`s`$, e.g. $`s=2.3`$ at $`T_{\mathrm{eff}}=70\mathrm{k}\mathrm{K}`$ and $`\mathrm{log}g=5`$. $`s`$ decreases towards lower and higher $`T_{\mathrm{eff}}`$ and is about $`1\text{ }\pm \text{ }0.2`$ at $`T_{\mathrm{eff}}`$$`>120\mathrm{kK}`$.

A similar grid for hydrogen-deficient stars (Fig. 5) has been calculated. For this grid we selected the “typical” abundance pattern (He:C:N:O=33:50:2:15 by mass) of PG 1159 stars (Werner 2001). This allows to estimate abundance effects onto the emergent flux.

## 3 Conclusions

Presently available NLTE model atmosphere codes have arrived at a very high level of sophistication (Mihalas 2003, Werner et al. 2003). The accelerated lambda iteration (ALI) method represents a powerful tool to calculate realistic metal-line blanketed atmospheres in a reasonable time, even on nowadays available PCs — these got such a tremendous number-crunching power in the last two years that they can now easily compete with CRAY and other (super)computers. The statistical approach which has been used to construct a generic model atom for all elements from Ca to Ni (Sect. 2) has been improved and can now be employed to construct model atoms for other light metals in order to consider all of their lines, provided by Kurucz (1996) and the Opacity Project (Seaton et al. 1994), in detail. Thus, there is no need to reflect on blackbodies or simple model atmospheres as ionizing sources for photoionization models.

To summarize, the use of metal-line blanketed NLTE model atmosphere fluxes, which consider at least the CNO opacities in detail, as ionizing spectra for photoionization models is highly recommended. The here presented fully metal-line blanketed NLTE model atmosphere fluxes which consider all elements from hydrogen up to the iron group are well suited for this purpose.

However, the user of these model atmosphere fluxes should be aware that they should only be used in a PN analysis if the CSPN abundances are in agreement with the model abundances. The use of fluxes from models with solar abundance ratios will spoil any PN analysis if the CSPN is hydrogen deficient because the flux shortwards of the He II ground state edge can be different by about one dex (Fig. 5). A preceding spectral analysis of the CSPN which yields $`T_{\mathrm{eff}}`$, $`\mathrm{log}g`$, and abundance ratios and an individually calculated small “grid” of NLTE model atmosphere fluxes may improve a subsequent PN analysis — and allows to construct consistent models of PN and CSPN (Sect. 1). This is also a prerequisite for reliable hydrodynamical models of the PN.

###### Acknowledgements.

I would like to thank Klaus Werner for many helpful discussions and careful reading of the manuscript, Jochen L. Deetjen for his continuous support to calculate the iron-group opacities, the Rechenzentrum der Universität Kiel, where the model atmosphere computations were performed on CRAY computers from 1997 to 2002, for generously awarding us the necessary computational time, and my colleagues from the Institut für Theoretische Physik und Astrophysik, Abteilung Astrophysik, at Kiel University for their patience and cooperation when the hugest batch classes on the CRAYs had been blocked for a long time. This research was supported by the DLR under grants 50 OR 9705 5 and 50 OR 0201.


text:
# Instabilities in the transient response of muscle

## I Introduction

The transient response of muscle to a sudden adjustment of its extension, or to an abrupt change in load, has been one of the most important sources of information about the mechanism of contraction for over three decades. Ever since the pioneering work of Huxley and Simmons, (1971), experimental data on transients has informed theoretical models of the interaction between myosin and actin (Huxley and Simmons,, 1971; Hill,, 1974; Eisenberg et al.,, 1980; Chen and Brenner,, 1993; Huxley,, 2000), providing a more detailed picture than could be obtained from the force-velocity relation (Hill,, 1939; Huxley,, 1957) alone. The reason is that the actomyosin interaction involves several processes which occur on different time scales, and these individual components can be resolved during the transient response (Huxley and Simmons,, 1971; Ford et al.,, 1977, 1985, 1986; Brenner,, 1991; Brenner et al.,, 1995). The quickest process is the elastic deformation of the myosin crossbridges that link the thick and thin filaments. Rapid transitions between two or more different bound states of the myosin molecule are thought to be the next fastest events, while detachment and reattachment of myosin heads occur on a slower time scale.

In an experiment to determine the isometric transient response, a muscle fiber is held at both ends to prevent it from contracting. The muscle is then suddenly shortened (or stretched) by a fixed amount, and the tension $`T`$ that it generates is measured. Immediately after the imposed change of length, the tension shifts from the isometric value $`T_0`$ to a new value, which is termed $`T_1`$. But shortly afterwards (typically within $`2\mathrm{ms}`$), the tension adjusts to a new value, termed $`T_2`$. Subsequently, it gradually reverts to the original isometric value $`T_0`$, and the entire transient response is usually completed in a fraction of a second. It is generally accepted that the initial response $`T_1`$ corresponds to the mechanical deformation of cross-bridges and provides a direct measure of their elasticity (Huxley and Simmons,, 1971). The interpretation of $`T_2`$ is rather more controversial. It is often attributed to force generation by the working stroke of bound myosin molecules (Huxley and Simmons,, 1971; Hill,, 1974; Eisenberg et al.,, 1980; Huxley and Tideswell,, 1996; Brenner et al.,, 1995; Duke,, 1999, 2000) and this interpretation has recently gained support from X-ray interference techniques applied to shortening fibers (Irving et al.,, 2000; Piazzesi et al., 2002a, ). But alternative models suggest that the force regeneration might be due, in part, to the rapid binding of new myosin heads to the thin filament (Brenner,, 1991; Howard,, 2001), or that it might involve the activation of the second myosin head (Huxley and Tideswell,, 1997).

In this article, we wish to address a fundamental problem connected with the interpretation of force transients. The present theories are all based on the consideration of a single pair of filaments, i.e. one filament containing myosin molecules, interacting with one actin filament. The dynamics of this filament pair is generalized to that of a whole muscle fiber by assuming that all filament pairs in a fiber behave in exactly the same way. This assumption is certainly justified as long as there are no static or dynamic instabilities in the system. However, the possibility of such instabilities has been known for a long time (Hill,, 1974). Moreover, a stochastic model of the actomyosin cycle, based on the swinging lever-arm hypothesis, has shown that instabilities do arise when values of parameters such as the lever-arm displacement and the crossbridge elasticity are chosen to provide effective energy transduction (Duke,, 1999, 2000). Such instabilities would give rise to a region of negative slope in the $`T_2`$ curve of a single filament pair. Several reasons have been advanced for the absence of any negative slope in the experimentally determined $`T_2`$ curve. Huxley and Simmons, (1971) argued that the power stroke is sub-divided into several small steps, and fixed the step size so that the $`T_2`$ curve had zero slope for limitingly small changes of length. In the model proposed by Eisenberg et al., (1980), the flatness of the $`T_2`$ curve was explained by a broad distribution of cross-bridge strain after attachment, combined with a specific strain-dependence of the transition rates to ensure the proper occupancies of the two bound states. A further explanation involved the compliance of the filaments and the distribution of binding sites on the thin filament in addition to a sub-divided power stroke (Huxley and Tideswell,, 1996). Duke, (1999) has suggested that the flat $`T_2`$ curve of a muscle fiber can arise despite an instability in the dynamics of a single pair of filaments, owing to the symmetry of a sarcomere. We investigate this possibility further in this article.

## II Methods

Simulation of the stochastic evolution of the system was performed using the Gillespie kinetic Monte Carlo algorithm, which works as follows. In each simulation step the rates of all possible transitions are calculated. The time until the next event is chosen as a random number with an exponential distribution and the expectation value given by the inverse of the sum of all rates. The event itself is chosen randomly with a statistical weight proportional to its rate. In the situations with stiff (non-compliant) filaments and continuous binding sites the transition rates can be factorized into factors that depend only on the backbone position (which are the same for all motors in a group) and factors that only depend on the binding position of a motor (which does not change with time unless that motor undergoes a transition). This allowed us to use a very efficient ($`n\mathrm{log}n`$) algorithm based on binary trees.

We assumed complete mechanical relaxation of the system in each step, i.e. the strain of all elastic elements is equilibrated before the next transition takes place. The structure involving sarcomeres, filaments and myosin heads was described as a circuit of elements with given resting lengths and compliances. The strain of every cross-bridge was calculated, given the constraint of fixed total length of the system (isometric conditions), or of fixed force acting on the ends of the system (isotonic conditions).

The $`T_2`$ transients were always determined $`2\mathrm{ms}`$ after the stretch/release. All other parameter values are summarized in Table 1.

## III Motors acting between a single pair of filaments

### III.1 Isometric transient in a swinging lever-arm model

We model the chemical cycle of a myosin head as described in Refs. (Duke,, 1999, 2000) and shown in Fig. 1. A molecule can exist in either a detached state (D) or in one of the attached states, A1 (with ADP$``$Pi) or A2 (with ADP). The state A1 is also referred to as the weak-binding state (Brenner,, 1991). The chemical transition between states A1 and A2 is concomittant with a conformational change of the molecule, in which the lever arm moves through displacement $`d`$. From state A2 the head undergoes another conformational change with a lever-arm shift $`\delta `$, associated with ADP release. In the ATP-rich physiological environment, this step is quickly followed by detachment of the head and we therefore skip the transitional state (A3) in our model. The attachment and detachment rates, which are summarized in Table 1, determine the shape of the force-velocity relation, but have little influence on the transient response to length steps, which mainly depend on the power-stroke displacement $`d`$ and the elastic constant $`K`$ of the myosin cross-bridge. We assume that the transitions A1$``$A2 take place on a faster time scale than the detachment and reattachment of heads. Therefore, once a head has bound with strain $`\xi `$ in the state A1, its state can be described as a statistical ensemble of the states A1 and A2 with probabilities given by the Boltzmann factors

$$P_2(\xi )=\left(1+\mathrm{exp}\frac{\mathrm{\Delta }G_{\mathrm{stroke}}+\frac{1}{2}Kd^2+Kd\xi }{k_BT}\right)^1,$$
(1)
where $`\mathrm{\Delta }G_{\mathrm{stroke}}`$ is the free energy change associated with phosphate release.

Now consider a filament containing $`N`$ myosin motors interacting with a actin filament which is held to prevent it from sliding. If $`N_A`$ heads are bound to the actin with a distribution of strains $`\mathrm{\Phi }(\xi )`$ (normalized to 1), then the isometric force that they generate is

$$T_0=N_A\mathrm{\Phi }(\xi )K(\xi +P_2(\xi )d)𝑑\xi .$$
(2)
If a stretch $`\mathrm{\Delta }x`$ is suddenly applied to the pair of filaments, the force will first change to

$$T_1(\mathrm{\Delta }x)=T_0+N_AK\mathrm{\Delta }x,$$
(3)
as the cross-bridges are deformed. In the next instance, transitions between states A1 and A2 occur and the probability distribution re-equilibrates to $`P_2(\xi +\mathrm{\Delta }x)`$. As a result the force adjusts to

$$T_2(\mathrm{\Delta }x)=N_A\mathrm{\Phi }(\xi )K(\xi +\mathrm{\Delta }x+P_2(\xi +\mathrm{\Delta }x)d)𝑑\xi .$$
(4)
As a first approximation, we can neglect the distribution of strains and assume $`\xi =0`$ for all myosin heads, which allows us to calculate the $`T_2`$ curves analytically. Their shape depends on two dimensionless parameters,

$`ϵ_0`$ $`={\displaystyle \frac{Kd^2}{2k_BT}},`$ (5)
$`ϵ_1`$ $`={\displaystyle \frac{Kd^2}{2\mathrm{\Delta }G_{\mathrm{stroke}}}},`$ (6)
which provide a measure of the energy stored in the elastic element when the head performs a power stroke (A1$``$A2). The parmeter $`ϵ_0`$ measures this energy relative to the thermal energy, while $`ϵ_1`$ compares it to the magnitude of the chemical free-energy change $`\mathrm{\Delta }G_{\mathrm{stroke}}`$ that accompanies the power-stroke transition. Note that if $`ϵ_1<1`$, the power stroke can occur immediately after the myosin head binds to the thin filament; but if $`ϵ_1>1`$, the conformational change is energetically inhibited initially, and it is only likely to occur once the thin filament has been moved forward by the action of other motors. On grounds of efficiency of energy transduction, we expect the value of $`ϵ_1`$ to be as high as possible, while ensuring that there is sufficient chemical energy to drive the power stroke, i.e. $`ϵ_11`$ (Duke,, 1999). The power-stroke displacement $`d`$ has been measured directly in a number of single-molecule experiments, and the values obtained range between $`5`$ and $`10\mathrm{nm}`$ (reviewed by Tyska and Warshaw,, 2002). Data on the crossbridge elasticity $`K`$ are less reliable, because there are inevitably other sources of compliance in the system (Veigel et al.,, 1998). But indirect evidence is provided by the energetic efficiency of muscle, which peaks at about $`50\%`$ (Kushmerick and Davies,, 1969; Barclay,, 1998; Piazzesi et al., 2002b, ). Thus a lower estimate for the energy stored in the elastic element following the power stroke is one half of the free-energy change accompanying the hydrolysis of an ATP molecule $`\mathrm{\Delta }G_{\mathrm{ATP}}20\mathrm{k}_\mathrm{B}\mathrm{T}`$, which leads to $`ϵ_010`$. Taking into account that not all energy stored in the spring can be converted to mechanical work one obtains a better estimate $`ϵ_020`$ (Duke,, 1999) (which is consistent with $`ϵ_11`$ if most of the energy of hydolysis is used to power the stroke $`\mathrm{\Delta }G_{\mathrm{stroke}}\mathrm{\Delta }G_{\mathrm{ATP}}`$). This corresponds to an elastic constant $`K=2.5\mathrm{pN}/\mathrm{nm}`$ if a power-stroke distance of $`d=8\mathrm{nm}`$ is assumed.

The form of the $`T_2`$ curve depends critically on the values of these dimensionless parameters.

If $`ϵ_0>2`$, an interval with a negative slope (a hysteresis) exists (Hill,, 1974). This means that under isotonic (constant load) conditions the system can be bistable. The location of the interval of negative slope depends on the value of $`ϵ_1`$. If the hysteresis spans the origin, then the state with $`\mathrm{\Delta }x=0`$ will be unstable under isotonic conditions. Putting $`\frac{dT_2}{d(\mathrm{\Delta }x)}|_{\mathrm{\Delta }x=0}<0`$ we see that this occurs if

$$\begin{array}{c}\left(1\frac{1}{ϵ_0}\mathrm{ln}\left(ϵ_01\sqrt{ϵ_0\left(ϵ_02\right)}\right)\right)^1<ϵ_1\hfill \\ \hfill <\left(1\frac{1}{ϵ_0}\mathrm{ln}\left(ϵ_01+\sqrt{ϵ_0\left(ϵ_02\right)}\right)\right)^1\end{array}$$
(7)
With a value of $`ϵ_11`$ the $`T_2`$ curve always has a negative slope around the stationary point. For example, with $`ϵ_0=20`$, Eq. (7) yields $`0.85<ϵ_1<1.22`$. However, it should be bourne in mind that this calculation does not take into account the distribution of strains on cross-bridges within the ensemble, and so the actual range might deviate slightly from this estimate.

An example of the $`T_2`$ transient for a group of myosin heads between two firmly clamped filaments is shown in Fig. 2.

### III.2 Group of stiffly coupled motors under near-isometric conditions

For a single pair of filaments, true isometric conditions imply that both filaments are held fixed. However, experimental conditions are usually near-isometric, in that the average velocity is zero, but the filaments still have the freedom to move. This situation occurs, for example, if both filaments are held at their ends by a flexible spring (in which case the motors cause the filaments to slide until the stall force is reached), or if the filaments suffer a constant load which is precisely adjusted to prevent net sliding. In this case the conditions are actually isotonic (or nearly so) and the instability discussed in the previous section implies that there is no steady state of the system with zero velocity when Eq. (7) is satisfied. Instead, as shown in Fig. 3, the pair of filaments oscillates in near-isometric conditions. The mechanism that generates the oscillations is the following. An ensemble of bound motors pulling against a constant load has two stable configurations: one with a majority of motors in state A2, the other with the majority in state A1. Because the detachment of motors is faster from state A2 than from state A1, the total number of bound motors decreases when the system is in the first configuration. As a result, the load per motor increases until the remaining motors cannot support the load any more, whereupon the system flips into the second configuration. Subsequently, the number of bound motors starts to grow again, and when the load per motor falls below a critical level the system flips back to the first configuration. Repetition of the cycle gives rise to an oscillation whose asymmetry reflects the differences in attachment and detachment rates.

Could the isometric transient be modified if it is measured in near-isometric conditions, instead of exactly isometric conditions? Fig. 2 shows the average $`T_2`$ response in the situation where the filaments are held at the ends by an external spring. In the non-oscillating case ($`ϵ_1=1.33`$, dashed line) the difference between this $`T_2`$ curve and the one in the strictly isometric state is only quantitative. The most important result is that the range of negative slope in $`T_2`$ remains. On the other hand, in the oscillating case ($`ϵ_1=0.8`$, dot-dashed line) the negative slope is flattened out because of the broader distribution of strains on cross-bridges.

### III.3 Effect of filament elasticity

The compliance of the thick and thin filaments is also expected to affect the $`T_2`$ curves. We model it by introducing a linear elasticity in the backbone connecting the heads (which is, to first order, equivalent to an elasticity in the track the motors are running on) (Vilfan et al.,, 1998). With $`\gamma `$ being the linear modulus of the filaments, which has the measured value $`\gamma =44000\mathrm{pN}`$ (Kojima et al.,, 1994), the spring constant of a filament segment of length $`L`$ is $`K_f=L^1\gamma `$. In (Vilfan et al.,, 1998) a linear two-state model was used to study the effect of filament elasticity on force-velocity relations, showing that they reduce the isometric force if the compliance of a filament segment between two bound motors becomes comparable to the motor spring constant $`K`$. Thin filaments are just stiff enough to prevent significant losses due to this effect. However, the filament compliance can have a significantly bigger effect on the transient response. The result of a simulation is shown in Fig. 4. Although the filament compliance almost halves the total stiffness of the system, it leaves the major part of the interval of negative slope unaffected.

### III.4 Effect of discrete binding sites

So far we have assumed that a myosin head can bind anywhere on the actin filament with equal probability. However, in reality the binding sites on actin are $`c=5.5\mathrm{nm}`$ apart. In addition, the actin filament has a helical structure with a half-pitch of about $`38\mathrm{nm}`$ (7 monomers). We therefore propose that the binding rate to a particular site is proportional to the Boltzmann factor determined by the spring distortion energy, composed of a longitudinal and a lateral component

$$k_{\mathrm{bind}}(\xi ,\varphi )=k_{\mathrm{bind}}\mathrm{exp}\left[\frac{K\xi ^2+K_A\varphi ^2}{2k_BT}\right]$$
(8)
with $`\xi =x(i+7j)c`$ and $`\varphi =\pi \frac{i}{7}`$, where $`j`$ denotes the repeat on the helix and $`i`$ the consecutive number of the actin monomer. The total binding rate at a given position $`x`$ then reads

$$k_{\mathrm{bind}}(x)=k_{\mathrm{bind}}\underset{i,j}{}\mathrm{exp}\left[\frac{K(x(i+7j)c)^2+K_A(\pi i/7)^2}{2k_BT}\right].$$
(9)
These curves have been measured experimentally using S1 myosin heads by Steffen et al., (2001) and the data were fitted with $`K_A=15\mathrm{pN}/\mathrm{nm}`$, the value which we use here. The position-dependence of the binding rate for these parameters is shown in Fig. 5. A similar distribution has also measured for myosin V (Veigel et al.,, 2002).

The dot-dashed line in Fig. 4 shows the $`T_2`$ response when taking into account discrete binding sites and the elasticity of actin filaments. As noted by Huxley and Tideswell, (1996), both the filament compliance and the discrete binding sites contribute to the flattening of the $`T_2`$ curve. However, our simulation indicates that their combined effect is not sufficient to cancel the negative slope with the parameters used here. We therefore conclude that the absence of negative slope in the measured $`T_2`$ curves cannot be explained in terms of a single-filament model, but requires taking into account the action of serially connected sarcomeres.

## IV Motors in muscle sarcomeres

### IV.1 Redistribution of filaments

The situation changes essentially if we take into account that the experiments are performed on whole muscle fibers, where several hundreds of sarcomeres are connected in series, as shown in Fig. 6d. The tension in each half-sarcomere must be the same, but the extension of each half-sarcomere can differ, provided that their sum equals the total stretch of the muscle fiber. To see what happens in this situation, it is instructive first to look at a minimum unit consisting of two filament pairs, joined back-to-back as shown in Fig. 6b. This unit represents two thin filaments and a thick filament, the basic building block of a sarcomere. A sudden stretch of the unit will first displace both halves equally, according to the elasticity of cross-bridges. But in the next phase, when the distribution of states A1 and A2 equilibrates, the position of the thick filament can become unstable. This means that a slight fluctuation in one direction will cause the motors pulling the thick filament in that direction to increase their force and those on the other side to decrease it, and the thick filament will jump sideways to one of the stable points. A schematic example is shown in Fig. 7.

The result of a simulation using different numbers of filament pairs, connected in series (the situation illustrated in Fig. 6c), is shown in Fig. 8. The interval of negative slope is increasingly flattened out as the number of pairs increases. Apart from the cancellation of the hysteresis, a shift in the curves is also observable, corresponding to the transition between the strictly isometric conditions and those where a filament is held under isotonic conditions before the stretch, cf. Fig. 2. The effect of redistribution can best be seen if one plots the distribution of stretches of individual half-sarcomeres against the total stretch per half-sarcomere (Fig. 9). For large stretch amplitudes the distribution is concentrated on the diagonal, in agreement with the expectation that all filaments experience the same stretch. For two filament pairs in series, the distribution splits in the intermediate range. The two peaks represent two stretches for which the single-filament $`T_2`$ curves produce the same force. Their position is determined by the condition that the average stretch must equal the value on the $`x`$-axis. With four filament pairs in series, a similar situation occurs, but with three different regimes. First three out of four pairs get displaced in one direction while the fourth is displaced in the other direction. Then this ratio changes to two against two and finally to one against three, before the regime with all filaments being subject to the same displacement settles in. With higher numbers of filaments in series (Fig. 9c shows the situation for 32) this pattern becomes increasingly continuous. In the regime where the negative slope stretches over the isometric point, the individual filaments perform small oscillations (Fig. 3) while the whole muscle is held at a constant length. These oscillations blur the distribution further (Fig. 9d).

$`T_2`$ curves for four different values of $`K`$ are shown in Fig. 10.

### IV.2 Isometric state

It is instructive to take a closer look at the parameters of the globally isometric state (by globally we mean that the total length of a fiber composed of many sarcomeres is constant, rather than the positions of individual filaments) in the two cases: one where the individual filaments oscillate and one where their positions just fluctuate. All data were determined in a situation where a large number (32 in the simulation) of filament pairs were connected in series and the total length held constant (situation (c) in Fig. 6). One expects the oscillating filaments to consume more ATP than the non-oscillating filaments. The results are summarized in Table 2. An exact comparison with experimental data is difficult, because the fraction of attached heads in the isometric state is not known reliably, and also because the experiments were carried out at different temperatures and with different types of myosin. For example, the value measured by Barsotti and Ferenczi, (1988), $`0.4`$ ATP molecules per second per head (with $`25\%`$ of the heads attached this gives $`1.6\mathrm{s}^1`$ per attached head), is compatible with the first scenario. However, other measurements give higher values (Ebus and Stienen,, 1996), more consistent with the second case. We therefore conclude that from the experimental values of the ATPase rate in the isometric state neither scenario can be excluded.

At high values of $`ϵ_1`$, the model shows an interesting feature. In the isometric state the number of cross-bridges in the state A2 can be very low, although the generated force per attached cross-bridge is as high as $`5\mathrm{pN}`$, in agreement with single-filament measurements (Kawai et al.,, 2000). This allows the muscle to support a load under isometric conditions with little ATP consumption (Duke,, 2000). It might at first sound paradoxical that most of the isometric force is generated by the pre-power-stroke state A1. This is possible because the myosin heads bind with a stochastic distribution of strains. Because those with a negative strain are more likely to undergo the power-stroke and then detach, the remaining ensemble produces a positive force. In this aspect, the power-stroke actually serves to eliminate negatively strained cross-bridges rather than direct generation of positive force. This notion is in agreement with recent experiments, which have shown that high phosphate concentration does reduce the isometric force significantly (Cooke and Pate,, 1985), but it does not have a visible effect on the conformation of the myosin heads or even their catalytic domains (Baker et al.,, 1999). But let us stress again that this holds only for the isometric state. The force in a contracting muscle originates mainly from the state A2.

### IV.3 Isotonic response

Another important class of experiments which provides information on the actomyosin interation involves the isotonic transient. Here the applied force is initially set at the value of the stalling load $`T_0`$, so that the fiber is prevented from contracting. The force is then suddenly changed and held constant at a different value, while the length of the fiber is recorded. Some early experiments showed that following a small step change of load, damped oscillations were imposed on the steady contraction or extension of the fiber (Podolsky,, 1960; Granzier et al.,, 1990). Such oscillations are particularly clear in recent experiments on single muscle fibers (Edman and Curtin,, 2001).

One possible cause of this oscillatory response has previously been suggested on the basis of the stochastic model of the actomyosin interaction used in this article. When $`ϵ_1>1`$, the chemical cycles of myosin motors on the same filament can become synchronized at loads close to the stalling force (Duke,, 1999, 2000). A pair of filaments then slides in a step-wise fashion under isotonic conditions. But during steady shortening, the motors on different filaments within the same muscle fiber operate out of phase, so that there is no macroscopic manifestation of the steps. However, an abrupt change in the load can cause the synchronization of a large fraction of the bound motors, whereupon the steps do become observable (Duke,, 1999). Because the correlation of the motors soon decays, the macroscopic steps fade and a damped oscillation is seen.

A much stronger oscillatory response is seen in the regime where individual pairs of filaments perform oscillations in the near-isometric state, as described in Sect. III.2. The synchronized oscillations can then be very pronounced following a small decrease in the load, as shown in Fig. 11a. On the other hand, no damped oscillations are observed after a larger drop in the load, e.g. to $`T_0/3`$ (Fig. 11b), because the individual filament pairs immediately move out of the hysteretic regime. These properties are in agreement with recent experiments on the isotonic response of single muscle fibers carried out by Edman and Curtin, (2001).

In steady isotonic conditions, another kind of instability can arise due to a hysteresis in the force-velocity relationship. The possibility of such an instability was first discussed in the context of a two-state ratchet model (Jülicher and Prost,, 1995) and subsequently in a kinetic cross-bridge model with a strain-dependent detachment rate (Vilfan et al.,, 1999). With the parameters we use this hysteresis, though existent, covers a rather small range of velocities and only results in a small inflexion in the force-velocity curve (Duke,, 1999, 2000).

## V Discussion

In the past, the isometric transient model of muscle has been modelled by considering the dynamics of a single pair of filaments. A significant problem with this approach is that, with a general choice of parameters in a cross-bridge model, the $`T2`$ curve is not flat (Hill,, 1974) — it typically has either a positive or a negative slope for limitingly small step displacements. Thus, in order to reproduce the experimental $`T2`$ curve of muscle, some of the parameters have needed to be finely adjusted (Huxley and Simmons,, 1971; Eisenberg et al.,, 1980), and limitations have been imposed on the size $`d`$ of the power-stroke and the rigidity $`K`$ of the cross-bridges (Huxley and Tideswell,, 1996). In this article we have investigated the isometric transient response using a stochastic model of the actomyosin cycle. When values of the parameters $`d`$ and $`K`$ are chosen to explain other characteristics of muscle, like the force-velocity relation and the efficiency of energy transduction, the $`T_2`$ curve of a single pair of filaments displays a region of negative slope, which cannot be eliminated by factors including filament compliance and the discrete nature of the binding sites on the actin filament. However, we argue that the symmetric structure of sarcomeres must be taken into consideration when computing the isometric transient response of muscle. Following a step change of length of a muscle fiber, a redistribution of sarcomere lengths occurs within the fiber. Some half-sarcomeres contract, while others extend. This redistribution always eliminates the negative slope, leading to a flat $`T_2`$ curve. It is a generic feature of unstable elements connected in series, and does not require any special values of model parameters.

Our model shows that there are two different regimes of the microscopic dynamics in near-isometric conditions. For a range of values of the parameter $`ϵ_1`$ close to unity, the isometric point falls in the interval where the slope of the $`T_2`$ curve is negative. In this case individual filament pairs oscillate with small amplitude. In the other regime, where the $`T_2`$ curve has positive slope at the isometric point, the individual filaments are stationary, apart from stochastic fluctuations. The macroscopic manifestations of these two regimes differ in few respects. Because the oscillations of different filament pairs have different phases, oscillatory motion is not normally observable on the scale of a whole muscle fiber in steady conditions. However, a sudden change of load can synchronize the oscillations and thereby make them visible. The existence of damped oscillations in the isotonic transient response of single muscle fibers (Edman and Curtin,, 2001) therefore argues in favor of the oscillating regime. We note, however, that damped oscillations can also be a manifestation of step-wise shortening (Duke,, 1999), which can exist in both regimes. Because the efficient transduction of energy demands $`ϵ_11`$, which is close to value of this parameter at the boundary of the two regimes, it is possible that both regimes exist depending on conditions such as the myosin isoform, phosphate concentration, pH, ionic strength and temperature. Indeed, measurements by Edman and Curtin, (2001) show a dependence of the oscillation decay on the solution pH and on muscle fatigue. Further experiments, in which conditions are systematically varied, could shed more light on the mechanism of oscillation.

As a final remark, we emphasize that according to our model, the $`T_2`$ curve of an individual pair of filaments differs from that of a muscle fiber. Recently, assays have been developed to measure the force-velocity relation of a single filament within a half-sarcomere (Kawai et al.,, 2000). In order to test our predictions, it would worthwhile to develop such high-precision techniques to measure the transient response of a single filament.

## Acknowledgment

A.V. would like to acknowledge support from the European Union through a Marie Curie Fellowship (No. HPMFCT-2000-00522) and from the Slovenian Office of Science (Grant No. Z1-4509-0106-02). T.D. acknowledges support from the Royal Society.