text:
# A Micromechanical Model of Hardening, Rate sensitivity and Thermal Softening in BCC Single Crystals

## 1 Introduction

The present paper is concerned with the development of a micromechanical model of the hardening, rate-sensitivity and thermal softening of bcc crystals. We place primary emphasis on the derivation of closed-form analytical expressions describing the macroscopic behavior of the crystals amenable to implementation as constitutive relations within a standard finite-element code. In developing the model, we follow the well-established paradigm of micromechanical modeling, consisting of: the identification of the dominant or rate-limiting ‘unit’ processes operating at the microscale; the identification of the macroscopic forces driving the unit processes; the analysis of the response of the unit processes to the macroscopic driving forces; and the determination of the average or macroscopic effect of the combined operation of all the micromechanical unit processes (see, e. g.,  for recent overviews and discussions pertaining to micromechanics and multiscale modeling of materials).

In formulating the present model we specifically consider the following unit processes: double-kink formation and thermally activated motion of kinks; the close-range interactions between primary and forest dislocations, and the subsequent formation of jogs; the percolation motion of dislocations through a random array of forest dislocations introducing short-range obstacles of different strengths; dislocation multiplication due to breeding by double cross-slip; and dislocation pair annihilation. We believe that this forms an ‘irreducible’ set of unit processes, in that each of these processes accounts for–and is needed for matching–salient and clearly recognizable features of the experimental record. For instance, consideration of dislocation multiplication and annihilation leads to a predicted stage I-II transition strain which decreases with temperature and increases with strain rate, as observed experimentally.

We bring an assortment of analysis tools to bear on each of the unit processes under consideration. As already mentioned, the choice of tools is to a large extent conditioned by our desire to derive closed-form analytical expressions for all constitutive relations. The motion of dislocations through an otherwise defect-free lattice is assumed to be thermally activated and controlled by the nucleation of kink pairs. A analysis of this process based on transition-state theory yields the effective Peierls stress of the lattice as a function of temperature and strain rate. The short-range interactions between a primary and secondary dislocation are assumed to result in the acquisition by both dislocations of a jog. The process by which the primary dislocation unpins from a forest obstacle is assumed to be thermally activated and the activation energy barrier is identified with the jog-formation energy. These assumptions result in temperature and rate-dependent obstacle strengths which are in good qualitative and quantitative agreement with experiment.

The percolation motion of dislocations through random arrays of point obstacles is studied in detail using statistical methods. The analysis presented here generalizes the analysis of Cuitiño and Ortiz  so as to account for multi-species distributions of finite-strength obstacles. Finally, we model dislocation multiplication as the result of two competing effects: the proliferation of dynamic sources by double cross-slip; and dislocation pair annihilation. We model cross-slip as a thermally activated process and develop an elastic model, similar to that of Huang et al. , which estimates the probability that the trajectories of two parallel screw dislocations collide, resulting in the annihilation of the pair. The model predicts that multiplication and annihilation balance out at sufficiently large strains, leading to saturation.

We validate the micromechanical model by recourse to detailed comparisons with the uniaxial tension tests on Ta single crystals of Mitchell and Spitzig . The model is found to capture salient features of the behavior of Ta crystals such as: the dependence of the initial yield point on temperature and strain rate; the presence of a marked stage I of easy glide, specially at low temperature and high strain rates; the sharp onset of stage II hardening and its tendency to shift towards lower strains, and eventually disappear, as the temperature increases or the strain rate decreases; the parabolic stage II hardening at low strain rates or high temperatures; the stage II softening at high strain rates or low temperatures; the trend towards saturation at high strains; the temperature and strain-rate dependence of the saturation stress; and the orientation dependence of the hardening rate.

## 2 General Framework

Our objective is to formulate a model of the hardening of bcc crystals which is well-suited to finite-deformation large-scale finite-element calculations and, therefore, we couch the model within a nonlinear kinematics framework. To this end, we adopt a conventional multiplicative elastic-plastic kinematics of the form

$$𝑭=𝑭^e𝑭^p$$
(1)
where $`𝑭`$ is the deformation gradient, $`𝑭^e`$ is the elastic lattice distortion and rotation, and $`𝑭^p`$ is the plastic deformation, which represents the net effect of crystallographic slip and leaves the lattice undistorted and unrotated. In addition, the plastic deformation is volume preserving. In materials such as metals, the elastic response is ostensibly independent of the internal processes and the free energy density decomposes additively as

$$A=W^e(𝑭^e,T)+W^p(T,𝑭^p,𝑸)$$
(2)
where $`T`$ is the absolute temperature and $`𝑸`$ is some suitable set of internal variables. The function $`W^e`$ determines the elastic response of the metal, e. g., upon unloading, whereas the function $`W^p`$ describes the hardening of the crystal. Physically, $`W^p`$ measures the energy stored in the crystal due to the plastic working of the material.

The examples of validation reported subsequently probe the hardening behavior of the material due to crystallographic slip over a range of moderate temperatures and strain rates. We therefore eschew a number of issues which become important at high pressures and elevated temperatures, but which play a limited role otherwise, including the volumetric equation of state (EoS), the pressure dependence of yield and elastic moduli, and others. First-principles calculations of the EoS and elastic moduli of bcc metals up to high pressures and temperatures may be found in . These results may be used as a basis for extending the present theory to conditions such as arise in shocked materials. These extensions notwithstanding, a simple form of the elastic energy density appropriate for present purposes is

$$W^e(𝑭^e,T)=\frac{1}{2}c_{ijkl}(T)(ϵ_{ij}^e\alpha _{ij}T)(ϵ_{kl}^e\alpha _{kl}T)$$
(3)
where $`𝒄`$ is the tensor of elastic moduli,

$$\mathit{ϵ}^e=\mathrm{log}\sqrt{𝑪^e}=\frac{1}{2}\mathrm{log}(𝑭_{}^{e}{}_{}{}^{T}𝑭^e)$$
(4)
is the logarithmic elastic strain, and $`𝜶`$ is the thermal expansion tensor. The elastic properties of a cubic crystal are fully described by the three Voigt constants $`c_{11}`$, $`c_{12}`$ and $`c_{44}`$ and a scalar thermal expansion coefficient $`\alpha `$. For low temperatures, the temperature dependence of the moduli may be taken to be linear in $`T`$ to a first approximation, leading to the relation:

$$c_{ijkl}(T)c_{ijkl}^0Tc_{ijkl}^1$$
(5)
This linear dependence of the elastic moduli is observed experimentally up to moderate temperatures , and can be justified within the framework of statistical mechanics .

Plastic deformations in single crystals are crystallographic in nature. Following Rice , we adopt a flow rule of the form

$$\dot{𝑭}^p=\left(\underset{\alpha =1}{\overset{N}{}}\dot{\gamma }^\alpha 𝒔^\alpha 𝒎^\alpha \right)𝑭^p$$
(6)
where $`\gamma ^\alpha `$ is the slip strain, and $`𝒔^\alpha `$, and $`𝒎^\alpha `$ are orthogonal unit vectors defining the slip direction and slip-plane normal corresponding to slip system $`\alpha `$. The collection $`𝜸`$ of slip strains may be regarded as a subset of the internal variable set $`𝑸`$. A zero value of a slip rate $`\dot{\gamma }^\alpha `$ signifies that the corresponding slip system $`\alpha `$ is inactive. The flow rule (6) allows for multiple slip, i. e., for simultaneous activity on more than one system over a region of the crystal. The vectors $`\{𝒔^\alpha ,𝒎^\alpha \}`$ remain constant throughout the deformation and are determined by crystallography. For bcc crystals, we consider the 24 slip systems listed in Table 1 in the Schmid and Boas nomenclature.

Slip on the $`\{112\}`$ systems is known to be asymmetric at low temperatures: slip is easier when the applied stress is such that dislocations move in the twinning direction (e. g., ). However, over the temperature range of interest here the experimental evidence  suggests that the extent of the twinning/anti-twinning asymmetry is relatively small, e. g., of the order of 20 MPa for Mo at 150 K , and, for simplicity, we will neglect it to a first approximation.

In the examples presented in Section 7, the constitutive equations are integrated in time using the variational update of Ortiz and Stainier . The variational formulation of the rate problem proposed by Ortiz and Stainier  additionally furnishes a convenient avenue for the superposition of the various unit mechanisms analyzed subsequently. Our working assumption is that the dissipation rates arising from these mechanisms are *additive*. This assumption in turn implies that the critical resolved shear stress $`\tau _c^\alpha `$ for the operation of system $`\alpha `$ may be computed as the sum of unit-process contributions, i. e.,

$$\tau _c^\alpha =\underset{\mathrm{processes}}{}\tau _c^\alpha (\mathrm{process})$$
(7)
where the sum extends to all unit processes. In the present model, these contributions are the Peierls stress of the lattice and the forest-obstacle resistance. Other dissipation mechanisms, such as phonon drag, may in principle be superposed likewise.

## 3 Dislocation mobility

In this section we begin by considering the thermally activated motion of dislocations within an *obstacle-free* slip plane. Under these conditions, the motion of the dislocations is driven by an applied resolved shear stress $`\tau `$ and is hindered by the lattice resistance, which is weak enough that it may be overcome by thermal activation. The lattice resistance is presumed to be well-described by a Peierls energy function, which assigns an energy per unit length to dislocation segments as a function of their position on the slip plane.

In bcc crystals, the core of screw dislocation segments relaxes into low-energy non-planar configurations . This introduces deep valleys into the Peierls energy function aligned with the Burgers vector directions and possessing the periodicity of the lattice. At low temperatures, the dislocations tend to adopt low-energy configurations and, consequently, the dislocation population predominantly consists of long screw segments. In order to move a screw segment normal to itself, the dislocation core must first be constricted, which requires a substantial supply of energy. Thus, the energy barrier for the motion of screw segments, and the attendant Peierls stress, may be expected to be large, and the energy barrier for the motion of edge segments to be comparatively smaller. For instance, Duesbery and Xu  have calculated the Peierls stress for a rigid screw dislocation in Mo to be 0.022$`\mu `$, where $`\mu `$ is the $`111`$ shear modulus, whereas the corresponding Peierls stress for a rigid edge dislocation is 0.006$`\mu `$, or about one fourth of the screw value. This suggests that the rate-limiting mechanism for dislocation motion is the thermally activated motion of kinks along screw segments ().

Consider now a screw dislocation segment, possibly pinned at both ends by an obstacle pair, subjected to a resolved shear stress $`\tau >0`$. At sufficiently high temperatures a double-kink may be nucleated with the assistance of thermal activation (e. g., , and the subsequent motion of the kinks causes the screw segment to effectively move forward, Fig. 1. Neglecting entropic effects and invoking Peach-Koehler’s formula, the forward and backward activation enthalpies $`\mathrm{\Delta }G^+`$ and $`\mathrm{\Delta }G^{}`$, respectively, for the formation of a double kink may be approximated as

$$\mathrm{\Delta }G^\pm E^{\mathrm{kink}}b\tau L^{\mathrm{kink}}l_P$$
(8)
where $`E^{\mathrm{kink}}`$ is the energy of formation of a kink-pair, $`L^{\mathrm{kink}}`$ is the length of an incipient double kink, and $`l_P`$ is the distance between two consecutive Peierls valleys. For bcc crystals, $`l_P=\sqrt{2/3}a`$ if the slip plane is $`\{110\}`$, $`l_P=\sqrt{2}a`$, if the slip plane is $`\{112\}`$, and $`l_P=\sqrt{8/3}a`$ if the slip plane is $`\{123\}`$, where $`a`$ is the cubic lattice size .

The formation energy $`E^{\mathrm{kink}}`$ cannot be reliably estimated from elasticity since is is composed mostly of core region. It can, however, be accurately computed by recourse to atomistic models. For instance, for Mo at zero stress Xu and Moriarty  have found formation energies $`E^{\mathrm{kink}}`$ of the order of 1 eV for kinks separated by a distance greater than $`L^{\mathrm{kink}}=15b`$. The core structure, gamma surfaces, Peierls stress, and kink-pair formation energies associated with the motion of $`a/2111`$ screw dislocations in Ta and Mo  have also been calculated by Moriarty et al.

We may expect a proliferation of kinks at a stress of the order of

$$\tau _0=\frac{E^{\mathrm{kink}}}{bL^{\mathrm{kink}}l_P}$$
(9)
for which $`\mathrm{\Delta }G=0`$. For Mo, Xu and Moriarty  have computed $`\tau _0`$ to be of the order of a few GPa. The activation free enthalpy for double kink nucleation may thus be rewritten as

$$\mathrm{\Delta }G^\pm (\tau )=E^{\mathrm{kink}}\left(1\frac{\tau }{\tau _0}\right),0\tau <\tau _0$$
(10)
A slightly more elaborate empirical formula that is widely used to fit activation energies was proposed by Kocks et al.  and the corresponding parameters have been determined for Ta by Tang al. , but this enhancement will not be pursued here.

Transition state theory predicts that a dislocation segment effects $`\nu _D\text{e}^{\beta G^+}`$ and $`\nu _D\text{e}^{\beta G^{}}`$ jumps per unit time in the positive and negative directions, respectively. Here, $`\beta =1/k_BT`$, $`k_B`$ is Boltzmann’s constant, $`T`$ is the absolute temperature, and $`\nu _D`$ is the attempt frequency which may be identified with the Debye frequency to a first approximation. Since the length of the individual jumps is $`l_P`$, the mean velocity of the dislocations then follows as

$$v=2l_P\nu _D\mathrm{e}^{\beta E^{\mathrm{kink}}}\mathrm{sinh}\left(\beta E^{\mathrm{kink}}\frac{\tau }{\tau _0}\right)$$
(11)
An application of Orowan’s formula then gives the slip strain rate as

$$\dot{\gamma }=\dot{\gamma }_0^{\mathrm{kink}}\mathrm{e}^{\beta E^{\mathrm{kink}}}\mathrm{sinh}\left(\beta E^{\mathrm{kink}}\frac{\tau }{\tau _0}\right)$$
(12)
where

$$\dot{\gamma }_0^{\mathrm{kink}}=2b\rho l_P\nu _D$$
(13)
is a reference strain rate, and $`\rho `$ is the dislocation density. In writing (12), we have taken into account the possibility of thermally activated jumps in the reverse direction, with forward and backward activation enthalpies $`\mathrm{\Delta }G^+`$ and $`\mathrm{\Delta }G^{}`$, respectively. For slip in the positive direction, $`\dot{\gamma }0`$, (12) may be inverted to give:

$$\frac{\tau _P}{\tau _0}=\frac{1}{\beta E^{\mathrm{kink}}}\mathrm{asinh}\left(\frac{\dot{\gamma }}{\dot{\gamma }_0^{\mathrm{kink}}}\mathrm{e}^{\beta E^{\mathrm{kink}}}\right)$$
(14)
where $`\tau _P`$ may be regarded as a temperature and rate-dependent *effective* Peierls stress.

Figure 2 illustrates the dependence of the effective Peierls stress on temperature and rate of deformation. The Peierls stress decreases ostensibly linearly up to a critical temperature $`T_c`$, beyond which it rapidly tends to zero. These trends are in agreement with the experimental observations of Wasserbäch  and Lachenmann and Schultz . The critical temperature $`T_c`$ increases with the strain rate. In particular, in this model the effect of increasing (decreasing) the strain rate has an analogous effect to decreasing (increasing) the temperature, and vice-versa, as noted by Tang et al. . In the regime of very high strain-rates ($`\dot{\gamma }>10^5s^1`$), effects such as electron and phonon drag become important and control the velocity of dislocations .

## 4 Forest hardening

In the forest-dislocation theory of hardening, the motion of dislocations, which are the agents of plastic deformation in crystals, is impeded by secondary –or ‘forest’– dislocations crossing the slip plane. As the moving and forest dislocations intersect, they form jogs or junctions of varying strengths  which, provided the junction is sufficiently short, may be idealized as point obstacles. Moving dislocations are pinned down by the forest dislocations and require a certain elevation of the applied resolved shear stress in order to bow out and bypass the pinning obstacles. The net effect of this mechanism is macroscopic hardening. Detailed numerical simulations of a dislocation line propagating through forest dislocations have been carried out by Foreman and Makin , and by Kocks , and more recently by Tang et al. . Analytical treatments of the model have been given by Kocks , Ortiz and Popov , and Cuitiño and Ortiz  for the case of infinitely strong obstacles. A phase-field model of the forest hardening mechanism has been proposed by Ortiz  and by Cuitiño et al. . Here we extend the statistical analysis of Cuitiño and Ortiz  to the case of several species of obstacles of finite strength.

Because of the random nature of the interactions, the motion of dislocations through a distribution of obstacles is best described in statistical terms . We begin by treating the case of infinitely strong obstacles. In this case, pairs of obstacles pin down dislocation segments, which require a certain threshold resolved shear stress $`s`$ in order to overcome the obstacle pair. Since the distribution of point obstacles within the slip plane is random, it follows that $`s`$ is itself a random variable. We shall let $`\stackrel{~}{f}_0^\alpha (s,t)`$ denote the probability density function of two-point barrier strengths on slip system $`\alpha `$ at time $`t`$. The time dependence of $`\stackrel{~}{f}_0^\alpha (s,t)`$ is a consequence of the variation in forest dislocation density.

In order to determine the precise form of $`\stackrel{~}{f}_0^\alpha (s,t)`$, we begin by noting that the Peierls energy landscape of bcc crystals strongly favors either screw or edge segments . In addition kinks, or points of change of direction of the dislocation line, carry a non-negligible amount of energy. It therefore follows that the lowest-energy configuration of unstressed dislocation segments spanning an obstacle pair is a step of the form shown in Fig. 3. Under these conditions, the bow-out mechanism by which a dislocation segment bypasses an obstacle pair may be expected to result in the configuration shown in Fig. 3 (bold line). If the edge-segment length is $`l_e`$, a displacement $`da_e`$ of the dislocation requires a supply of energy equal to $`2U^{\mathrm{screw}}da_e+b\tau _P^{\mathrm{edge}}l_eda_e`$ in order to overcome the Peierls resistance $`\tau _P^{\mathrm{edge}}`$ and to extend the screw segments. The corresponding energy release is $`b\tau l_eda_e`$. Similar contributions result from a displacement $`da_s`$ of the screw-segment of length $`l_s`$. Therefore, for the bow-out of the dislocation to be energetically possible we must have

$$b\tau l_eda_e+b\tau l_sda_s2U^{\mathrm{screw}}da_e+b\tau _P^{\mathrm{edge}}l_eda_e+2U^{\mathrm{edge}}da_s+b\tau _P^{\mathrm{screw}}l_sda_s$$
(15)
As already noted, in bcc crystals the core energy $`U^{\mathrm{screw}}`$ per unit length of the screw segments is smaller than the core energy $`U^{\mathrm{edge}}`$ per unit length of the edge segments (e. g., ). Conversely, the Peierls resistance $`\tau _P^{\mathrm{screw}}`$ to the glide of screw segments is larger than the Peierls resistance $`\tau _P^{\mathrm{edge}}`$ for edge segments (e. g., ). Retaining dominant terms only, Eq. (15) simplifies to

$$\tau s=\tau _P^{\mathrm{screw}}+\frac{2U^{\mathrm{edge}}}{bl_s}$$
(16)
We may further identify the Peierls stress $`\tau _P^{\mathrm{screw}}`$ with (14), whence it follows that $`s`$ comprises a thermally-activated and rate-sensitive term $`\tau _P`$ and an athermal and rate-insensitive term $`2U^{\mathrm{edge}}/bl_s`$.

We shall assume that the point obstacles are randomly distributed over the slip plane with a mean density $`n^\alpha `$ of obstacles per unit area. We shall also assume that the obstacle pairs spanned by dislocation segments are nearest-neighbors in the obstacle ensemble. Thus, if $`r`$ is the distance between the obstacles in an obstacle pair, then the circle of radius $`r`$ centered at either obstacle contains no other obstacles. Under these conditions, the probability density of $`r`$ is given by

$$\stackrel{~}{f}_0^\alpha (r,t)=2\pi n^\alpha r\mathrm{exp}(\pi n^\alpha r^2)$$
(17)
where the dependence of $`\stackrel{~}{f}_0`$ on time stems from the time-dependence of the obstacle density $`n^\alpha `$.

In order to deduce the probability density of $`l_s`$, and by extension of $`s`$, we note that, given a point obstacle located at the origin, the probability of finding another obstacle in the element of area $`rdrd\theta `$ is:

$$d\stackrel{~}{P}_0^\alpha (r,t)=n^\alpha r\mathrm{exp}(\pi n^\alpha r^2)drd\theta $$
(18)
Changing variables to Cartesian coordinates, $`x=r\mathrm{cos}\theta `$, $`y=r\mathrm{sin}\theta `$, (18) may be recast in the form:

$$d\stackrel{~}{P}_0^\alpha =n^\alpha \mathrm{exp}[\pi n^\alpha (x^2+y^2)]dxdy$$
(19)
The frequency of obstacle pairs with screw-segment length $`l_s`$ is, therefore:

$$\stackrel{~}{f}_0^\alpha (l_s,t)=2_{\mathrm{}}^{\mathrm{}}n^\alpha \mathrm{exp}[\pi n^\alpha (x^2+l_s^2)]𝑑x=\sqrt{n^\alpha }\mathrm{exp}(\pi n^\alpha l_s^2)$$
(20)
The corresponding probability of finding an obstacle pair of strength $`s`$ follows from (20) and (16) as:

$$\stackrel{~}{f}_0^\alpha (s,t)=\frac{4U^{\mathrm{edge}}\sqrt{n^\alpha }}{b(s\tau _P^\alpha )^2}\mathrm{exp}\left[\frac{\pi n^\alpha }{b^2}\left(\frac{2U^{\mathrm{edge}}}{s\tau _P^\alpha }\right)^2\right],s\tau _P^\alpha $$
(21)
and the associated distribution function is:

$$\stackrel{~}{P}_0^\alpha (s,t)=_{\tau _P}^s\stackrel{~}{f}_0^\alpha (\xi ,t)𝑑\xi =\left[1\mathrm{erf}\left(\frac{2\sqrt{\pi n^\alpha }U^{\mathrm{edge}}}{b(s\tau _P^\alpha )}\right)\right]$$
(22)
We note that the Peierls stress $`\tau _P^\alpha `$ depends on the slip system $`\alpha `$ through its dependence on the slip-strain rate $`\dot{\gamma }^\alpha `$.

It is interesting to note that the probability density $`\stackrel{~}{f}_0^\alpha (s,t)`$ of obstacle-pair strengths just derived for bcc crystals differs markedly from those which are obtained for fcc crystals by a similar argument , namely:

$$\stackrel{~}{f}_0^\alpha (s,t)=\frac{2\pi n^\alpha U^2}{b^2s^3}\mathrm{exp}\left(\frac{n^\alpha \pi }{b^2}\frac{U^2}{s^2}\right)$$
(23)
where $`U`$ is the dislocation core energy per unit length of dislocation. This difference owes to the different bow-out configurations for the two crystal classes and the comparatively larger values of the Peierls stress in bcc crystals. Thus, the Peierls stress of fcc crystals is generally quite small and, as in the derivation of (23) is often neglected entirely to a first approximation. For bcc crystals, the effective Peierls stress $`\tau _P^\alpha `$ decreases with increasing temperature and decreasing rate of deformation, Eq. (14), and hence the behavior of bcc crystals may be expected to be closer to that of fcc crystals under those conditions, as noted by Tang et al. . Conversely, the hardening behavior of bcc crystals may be expected to differ sharply from that of fcc crystals at low temperatures and high rates of deformation.

The function $`\stackrel{~}{f}_0^\alpha (s,t)`$ just derived provides a complete description of the distribution of the obstacle-pair strengths when the point obstacles are of infinite strength and, consequently, impenetrable to the dislocations. Next we extend the preceding analysis to the case of finite obstacle strengths. Let $`s^{\alpha \beta }`$ be the strength of the jogs or junctions formed by dislocations of systems $`\alpha `$ and $`\beta `$. Now consider an obstacle pair in which the weakest point obstacle corresponds to a forest dislocation of type $`\beta `$ and, therefore, has strength $`s^{\alpha \beta }`$. The probability that the strength of the obstacle pair be $`s`$ is, therefore,

$$\stackrel{~}{f}^\alpha (s|s^{\alpha \beta },t)=\frac{\stackrel{~}{f}_0^\alpha (s,t)}{\stackrel{~}{P}_0^\alpha (s^{\alpha \beta },t)}\left[1H(ss^{\alpha \beta })\right]$$
(24)
where $`H(s)`$ is the Heaviside function. The probability that an obstacle on system $`\alpha `$ be of type $`\beta \alpha `$ is $`n^{\alpha \beta }/n^\alpha `$, where $`n^{\alpha \beta }`$ is the number of obstacles of type $`\beta `$ per unit area of the slip plane $`\alpha `$, and

$$n^\alpha (t)=\underset{\beta \alpha }{}n^{\alpha \beta }(t)$$
(25)
is the total obstacle density on slip system $`\alpha `$.

Next we note that the probability that the weakest of the two obstacles forming a obstacle pair be of type $`\beta `$ is:

$$p^{\alpha \beta }=\frac{n^{\alpha \beta }}{n^\alpha }\left(\frac{n^{\alpha \beta }}{n^\alpha }+2\underset{s^{\alpha \beta ^{}}>s^{\alpha \beta }}{}\frac{n^{\alpha \beta ^{}}}{n^\alpha }\right),\beta \alpha ,\beta ^{}\alpha $$
(26)
It is readily verified that

$$\underset{\beta \alpha }{}p^{\alpha \beta }=\left[\underset{\beta \alpha }{}\left(\frac{n^{\alpha \beta }}{n^\alpha }\right)\right]^2=1$$
(27)
as required. Finally, the probability that an obstacle pair have strength $`s`$ follows as

$$\stackrel{~}{f}^\alpha (s,t)=\underset{\beta \alpha }{}p^{\alpha \beta }(t)\stackrel{~}{f}^\alpha (s|s^{\alpha \beta },t)$$
(28)
This probability distribution function jointly accounts for the strength of the obstacle pairs due to line tension and to the obstacle strength. Evidently, $`\stackrel{~}{f}^\alpha (s,t)`$ is supported in the interval $`[0,s_{\mathrm{max}}]`$, where $`s_{\mathrm{max}}`$ is the maximum obstacle strength. Making use of (24), the probability density function (28) can be rewritten in the form

$$\stackrel{~}{f}^\alpha (s,t)=\underset{\beta \alpha }{}p^{\alpha \beta }(t)\stackrel{~}{f}_0^\alpha (s,t)\frac{1H(ss^{\alpha \beta })}{\stackrel{~}{P}_0^\alpha (s^{\alpha \beta },t)}$$
(29)
and the associated distribution function becomes

$$\stackrel{~}{P}^\alpha (s,t)=\underset{\beta \alpha }{}p^{\alpha \beta }(t)\left\{\stackrel{~}{P}_0^\alpha (s,t)\frac{1H(ss^{\alpha \beta })}{\stackrel{~}{P}_0^\alpha (s^{\alpha \beta },t)}+H(ss^{\alpha \beta })\right\}$$
(30)
We assume that $`n^{\alpha \beta }`$, the number of obstacles of type $`\beta `$ per unit area of the slip plane $`\alpha `$, scales with the dislocation density $`\rho ^\beta `$ according to the relation:

$$n^{\alpha \beta }=a^{\alpha \beta }\rho ^\beta $$
(31)
where the coupling constants $`a^{\alpha \beta }`$ are regarded here as purely geometrical parameters. A simple geometrical argument based on counting intersections of randomly distributed lines with a slip plane gives :

$$n^{\alpha \beta }=\frac{2}{\pi }\sqrt{1(𝒎^\alpha 𝒎^\beta )^2}\rho ^\beta $$
(32)
where $`𝒎^\alpha `$ and $`𝒎^\beta `$ are the unit normals to slip planes of types $`\alpha `$ and $`\beta `$ and $`\rho ^\beta `$ is the dislocation line density per unit volume in slip system $`\beta `$. A comparison between (31) and (32) suggests writing:

$$a^{\alpha \beta }=a_0\frac{2}{\pi }\sqrt{1(𝒎^\alpha 𝒎^\beta )^2}$$
(33)
where $`a_0<1`$ is an ‘efficiency’ factor which accounts for the tendency of dislocations to tangle and form loops, which in turn tends to lower the number of slip-plane crossings.

Let $`\tau ^\alpha (t)`$ now be the resolved shear stress acting on the slip system $`\alpha `$ at time $`t`$. Assume for now that $`\tau ^\alpha (t)`$ increases monotonically from zero at $`t=0`$. Evidently, for dislocations to be stable at time $`t`$ they must face barriers of strengths $`s`$ in excess of $`\tau ^\alpha (t)`$. As $`\tau ^\alpha (t)`$ is increased to $`\tau ^\alpha (t)+\dot{\tau }^\alpha (t)dt`$, the dislocation segments held at barriers of strengths in the range $`\tau ^\alpha (t)s\tau ^\alpha (t)+\dot{\tau }^\alpha (t)dt`$ are dislodged and move forward until they reach barriers of strength $`s\tau ^\alpha (t)+\dot{\tau }^\alpha (t)dt`$. This motion of dislocations results in an net increase in the plastic deformation. The dislocation jumps between obstacles are assumed to be instantaneous. This idealization is justified when the duration of the flights is much smaller than the characteristic time of variation of the loads.

As noted by Ortiz and Popov , the information needed to describe the dislocation motion is fully contained in the probability density function $`f^\alpha (s,t)`$, which represents the fraction of dislocation length facing obstacle pairs of strength $`s`$ at time $`t`$. The function $`f^\alpha (s,t)`$ evolves in time due to the process of redistribution of dislocation line described above. Initially, though, the dislocations may be assumed to be randomly distributed over their slip plane, and $`f^\alpha (s,0)=\stackrel{~}{f}^\alpha (s,0)`$. At later times, $`f^\alpha (s,t)`$ must vanish identically for $`0s<\tau ^\alpha (t)`$ in the rate independent limit.

A kinetic equation for the evolution of $`f^\alpha (s,t)`$ was derived by Ortiz and Popov  using standard tools of non-equilibrium statistical mechanics. Cuitiño and Ortiz  were able to obtain an analytical solution to this kinetic equation for the case of monotonic loading and an arbitrary time variation of the density of point obstacles. Remarkably, the solution takes the simple closed form

$$f^\alpha (s,t)=\frac{\stackrel{~}{f}^\alpha (s,t)}{1\stackrel{~}{P}^\alpha (\tau ^\alpha (t),t)}H\left(s\tau ^\alpha (t)\right)$$
(34)
Solution (34) implies that, under the conditions of the analysis, the probability density $`f^\alpha (s,t)`$ remains proportional to $`\stackrel{~}{f}^\alpha (s,t)`$ over the current admissible range $`[\tau ^\alpha (t),s_{\text{max}}]`$.

Let $`\rho ^\alpha (t)`$ denote the current dislocation length per unit volume for the slip system $`\alpha `$. The dislocation density released during an increment of the resolved shear stress from $`\tau ^\alpha (t)`$ to $`\tau ^\alpha (t)+\dot{\tau }^\alpha (t)dt`$ gives rise to an incremental plastic strain

$$d\gamma ^\alpha (t)=b\rho ^\alpha (t)f(\tau ^\alpha (t),t))\dot{\tau }^\alpha (t)dt\overline{N}(t)\overline{l}(t)$$
(35)
where $`\overline{l}(t)`$ is the average distance between obstacles, and $`\overline{N}(t)`$ is the average number of jumps the dislocation segments make before attaining stable positions.

To a good approximation, $`\overline{l}(t)`$ can be identified with the average distance between point obstacles. Assuming that the obstacles are randomly distributed over the slip plane with density $`n^\alpha (t)`$, a straightforward derivation from (17) gives

$$\overline{l}(t)=l(t)=\frac{1}{2\sqrt{n^\alpha (t)}}$$
(36)
Next we compute the average number of jumps $`\overline{N}(t)`$. Evidently, the probability that an unstable segment becomes arrested after the first jump is equal to the probability that the first barrier encountered is of a strength $`s\tau ^\alpha (t)`$. This probability is $`1\stackrel{~}{P}^\alpha (\tau ^\alpha (t),t)`$. The probability that the segment goes beyond the first barrier is $`\stackrel{~}{P}^\alpha (\tau ^\alpha (t),t)`$. Likewise, the probability that a segment gets arrested at the second barrier encountered is $`\stackrel{~}{P}^\alpha (\tau ^\alpha (t),t)[1\stackrel{~}{P}^\alpha (\tau ^\alpha (t),t)]`$, and the probability that it goes beyond is $`\stackrel{~}{P}^\alpha {}_{}{}^{2}(\tau ^\alpha (t),t)`$, and so on. Hence, the average number of jumps between barriers taken by an unstable segment is

$$\begin{array}{cc}\hfill \overline{N}(t)& =[1\stackrel{~}{P}^\alpha (\tau ^\alpha (t),t)]+2\stackrel{~}{P}^\alpha (\tau ^\alpha (t),t)[1\stackrel{~}{P}^\alpha (\tau ^\alpha (t),t)]\hfill \\ & +3\stackrel{~}{P}^\alpha {}_{}{}^{2}(\tau ^\alpha (t),t)[1\stackrel{~}{P}^\alpha (\tau ^\alpha (t),t)]+\mathrm{}=\frac{1}{1\stackrel{~}{P}^\alpha (\tau ^\alpha (t),t)}\hfill \end{array}$$
(37)
Interestingly, if $`\tau ^\alpha =0`$, any moving dislocation segment is sure to be arrested at the first obstacle it encounters, and $`\overline{N}=1`$, in agreement with (37). Likewise, if $`\tau ^\alpha (t)>s_{\mathrm{max}}`$, then the segment never reaches a stable barrier and $`\overline{N}\mathrm{}`$, as predicted by (37).

Substituting (37) and (34) with $`s=\tau ^\alpha (t)`$ into (35) we finally obtain

$$\dot{\gamma }^\alpha (t)=\gamma _\text{c}^\alpha (t)\frac{\stackrel{~}{f}^\alpha (\tau ^\alpha (t),t)}{[1\stackrel{~}{P}^\alpha (\tau ^\alpha (t),t)]^2}\dot{\tau }^\alpha (t)$$
(38)
where the characteristic plastic strain $`\gamma _\text{c}^\alpha (t)`$ is defined as

$$\gamma _\text{c}^\alpha (t)=\frac{b\rho ^\alpha (t)}{2\sqrt{n^\alpha (t)}}$$
(39)
Equation (38) defines a relation of the form

$$\dot{\gamma }^\alpha (t)=\frac{\dot{\tau }^\alpha (t)}{h^\alpha (t)}$$
(40)
where

$$h^\alpha (t)=\frac{1}{\gamma _\text{c}^\alpha (t)}\frac{[1\stackrel{~}{P}^\alpha (\tau ^\alpha (t),t)]^2}{\stackrel{~}{f}^\alpha (\tau ^\alpha (t),t)}$$
(41)
is the hardening modulus of the slip system $`\alpha `$. By way of example, in the particular case of obstacles of uniform strength, $`s^{\alpha \beta }=s_{\text{max}}`$, the hardening modulus takes the form:

$$\begin{array}{c}h^\alpha (t)=h_c^\alpha (t)\left(\frac{\tau ^\alpha (t)\tau _P^\alpha }{\tau _c^\alpha }\right)^2\frac{1}{\sqrt{\pi }}\mathrm{exp}\left[\left(\frac{\tau _c^\alpha }{\tau ^\alpha (t)\tau _P^\alpha }\right)^2\right]\hfill \\ \hfill \frac{[\stackrel{~}{P}_0^\alpha (s_{\mathrm{max}},t)\stackrel{~}{P}_0^\alpha (\tau ^\alpha ,t)]^2}{\stackrel{~}{P}_0^\alpha (s_{\mathrm{max}},t)}(1H(\tau ^\alpha s_{\mathrm{max}}))\end{array}$$
(42)
where

$$\tau _\text{c}^\alpha (t)=2\sqrt{\pi n^\alpha }\frac{U^{\mathrm{edge}}}{b}\text{and}h_\text{c}^\alpha (t)=\frac{\tau _\text{c}^\alpha (t)}{\gamma _\text{c}^\alpha (t)}$$
(43)
are a characteristic shear stress and plastic modulus, respectively. Equation (42) predicts an initial infinite hardening modulus at $`\tau =\tau _P`$. The hardening modulus subsequently decreases monotonically to zero as $`\tau `$ approaches $`s_{\text{max}}`$.

As shown in , the above relations can be extended simply to account for elastic unloading, provided that $`\tau ^\alpha (t)`$ does not change sign at any time during the loading history. This extension consists in defining the current flow stress $`g^\alpha (t)`$ of slip system $`\alpha `$ as the maximum previously attained value of $`\tau ^\alpha (t)`$. The hardening relations are then rewritten in the form

$$\dot{\gamma }^\alpha (t)=\frac{\dot{g}^\alpha (t)}{h^\alpha (t)}$$
(44)
and

$$h^\alpha (t)=\frac{1}{\gamma _\text{c}^\alpha (t)}\frac{[1\stackrel{~}{P}^\alpha (g^\alpha (t),t)]^2}{\stackrel{~}{f}^\alpha (g^\alpha (t),t)}$$
(45)
It should be noted that, since $`\dot{g}^\alpha (t)=0`$ whenever either $`\tau ^\alpha (t)<g^\alpha (t)`$ or $`\dot{\tau }^\alpha (t)0`$, relation (44) implicitly accounts for elastic unloading. To account for the initial distribution of dislocations and obstacles, we take $`g^\alpha (t)=g_0`$.

## 5 Dislocation intersections

In this section we proceed to estimate the obstacle strengths which enter in the preceding analysis of forest hardening. The interaction between primary and secondary dislocations may result in a variety of reaction products, including jogs and junctions . Experimental estimates of junction strengths have been given by Franciosi and Zaoui  for the twelve slip systems belonging to the family of $`\{111\}`$ planes and $`[110]`$ directions in fcc crystals, and by Franciosi  for the twenty-four systems of types $`\{211\}`$ $`[111]`$ and $`\{110\}`$ $`[111]`$ in bcc crystals. The strength of some of these interactions has recently been computed using atomistic and continuum models . Tang et al. have numerically estimated the average strength of dislocation junctions for Nb and Ta crystals .

For purposes of the present theory, we specifically concern ourselves with short-range interactions between dislocations which can be idealized as point defects. For simplicity, we consider the case in which each intersecting dislocation acquires a jog. The energy of a pair of crossing dislocations is schematically shown in Fig. 4 as a function of some convenient reaction coordinate, such as the distance between the dislocations. The interaction may be repulsive, resulting in an energy barrier, or attractive, resulting in a binding energy, Fig. 4. In the spirit of an equilibrium theory, here we consider only the final reaction product, corresponding to a pair of jogged dislocations at infinite distance from each other, and neglect the intermediate states along the reaction path. In addition, we deduce the strength of the obstacles directly from the energy supply required to attain the final state, i. e. the jog-formation energy. Despite the sweeping nature of these assumptions, the predicted saturation strengths in multiple slip are in good agreement with experiment (cf Section 7), which lends some empirical support to the theory.

We estimate the jog formation energy as follows. Based on energy and mobility considerations already discussed, we may expect the preponderance of forest dislocations to be of screw character, and the mobile dislocation segments to be predominantly of edge character. We therefore restrict our analysis to intersections between screw and edge segments. The geometry of the crossing process is schematically shown in Fig. 5. Each dislocation acquires a jog equal to the Burgers vector of the remaining dislocation. The energy expended in the formation of the jogs may be estimated as

$$E_{\alpha \beta }^{\mathrm{jog}}|𝒃^\beta |\left[U^{\alpha \beta }U^{\alpha \alpha }\mathrm{cos}\theta ^{\alpha \beta }\right]+|𝒃^\alpha |\left[U^{\beta \alpha }U^{\beta \beta }\mathrm{cos}\theta ^{\beta \alpha }\right]$$
(46)
Here, $`\alpha `$ designates the slip system of the moving edge segment, $`\beta `$ the slip system of the forest screw dislocation, $`U^{\alpha \alpha }=U^{\mathrm{edge}}`$ is the energy per unit length of an edge segment in the slip system $`\alpha `$, $`U^{\beta \beta }=U^{\mathrm{screw}}`$ is the energy per unit length of a screw segment in slip system $`\beta `$, $`U^{\alpha \beta }`$ is the energy per unit length of a segment in slip system $`\alpha `$ aligned with $`𝒃^\beta `$, $`U^{\beta \alpha }`$ is the energy per unit length of a segment in slip system $`\beta `$ aligned with $`𝒃^\alpha `$, $`\theta ^{\alpha \beta }`$ is the angle between the edge direction in system $`\alpha `$ and $`𝒃^\beta `$, and $`\theta ^{\beta \alpha }`$ is the angle between the screw direction in system $`\beta `$ and $`𝒃^\alpha `$.

We additionally assume that the screw orientation defines a cusp in the dependence of the dislocation line energy on segment orientation, and that small deviations of a segment from a pure screw character raise its energy to a level comparable to $`U^{\mathrm{edge}}`$. This gives the energy estimate

$$U^{\alpha \beta }=U^{\beta \alpha }=\{\begin{array}{cc}U^{\mathrm{screw}}\hfill & \text{if }𝒃^\alpha =𝒃^\beta \hfill \\ U^{\mathrm{edge}}\hfill & \text{otherwise}\hfill \end{array}$$
(47)
Inserting these energies into (46) gives

$$E_{\alpha \beta }^{jogs}\{\begin{array}{cc}bU^{\mathrm{screw}}\left[1r\mathrm{cos}\theta ^{\alpha \beta }\right]\hfill & \text{if }𝒃^\alpha =𝒃^\beta \hfill \\ bU^{\mathrm{screw}}\left[2r\mathrm{cos}(\theta ^{\beta \alpha })r\mathrm{cos}\theta ^{\alpha \beta }\right]\hfill & \text{otherwise}\hfill \end{array}$$
(48)
where $`r=U^{\mathrm{edge}}/U^{\mathrm{screw}}`$ is the ratio of screw to edge dislocation line energies. For Ta, recent atomistic calculations  give a value of $`r=1.77`$. The resulting jog formation energies for the complete collection of pairs of $`\{211\}`$ and $`\{110\}`$ dislocations are tabulated in Table 2.

The net effect of jog formation on hardening may be ascertained as follows. Consider the case in which a primary system $`\alpha `$ contains forest obstacles of a single species, corresponding to secondary dislocations in the slip system $`\beta `$. Let $`n^{\alpha \beta }`$ be the density of forest obstacles per unit area of the primary plane $`\alpha `$. As a primary dislocation sweeps through an area $`A`$, the energy expended in forming jogs with forest dislocations of the $`\beta `$-type is $`E_{\alpha \beta }^{\mathrm{jog}}n^{\alpha \beta }A`$. On the other hand, the potential energy released as a result of the motion of the dislocation follows from the Peach-Koehler formula as $`b\tau ^\alpha A`$, where $`\tau ^\alpha `$ is the resolved shear stress acting on the primary system $`\alpha `$. Hence, the forest obstacles of type $`\beta `$ become ‘transparent’ to the motion of primary dislocations when $`b\tau ^\alpha AE_{\alpha \beta }^{\mathrm{jog}}n^{\alpha \beta }A`$, or

$$\tau ^\alpha \frac{1}{b}E_{\alpha \beta }^{\mathrm{jog}}n^{\alpha \beta }$$
(49)
Since the jog energies scale with the elastic moduli, they may be expected to reduce to zero at the melting temperature, which to some extent accounts for the observed thermal softening. In addition, the jog energies are small enough for thermal activation to be operative at the level of individual obstacles. A derivation entirely analogous to that leading to Eq. (14) yields, in this case,

$$\frac{s^{\alpha \beta }}{s_0^{\alpha \beta }}=\frac{1}{\beta E_{\alpha \beta }^{\mathrm{jog}}}\mathrm{asinh}\left(\frac{\dot{\gamma }^\alpha }{\dot{\gamma }_0^\alpha }\mathrm{e}^{\beta E_{\alpha \beta }^{\mathrm{jog}}}\right)$$
(50)
where

$$s_0^{\alpha \beta }=\frac{E_{\alpha \beta }^{\mathrm{jog}}}{b\overline{l}^\alpha L^{\mathrm{junct}}}$$
(51)
and

$$\dot{\gamma }_0^\alpha =2\rho ^\alpha b\overline{l}^\alpha \nu _D$$
(52)
The lengths $`\overline{l}^\alpha `$ and $`L^{\mathrm{junct}}`$ describe the geometry of the junction, as illustrated in Fig. 6, and are of the order of a few Burgers vectors.

It follows from these developments that, at low temperatures, the energies collected in Table 2 provide a measure of the corresponding obstacle strengths. We may also recall that the saturation flow stress of a slip system is directly related to the strengths of the obstacles. Consequently, information about obstacle strengths may be inferred from saturation flow stresses in crystals deformed in double slip. Measurements of this type have been performed by Franciosi for $`\alpha `$-Fe . The general trends exhibited by Table 2 are consistent with Franciosi’s data, as well as with the recent analytical model of Lee et al. . In addition, the saturation flow stresses resulting from the model developed above are in good quantitative agreement with those measured experimentally, cf Section 7, which lends empirical support to jog-formation as a plausible mechanism underlying short-range obstacle strength and saturation in bcc crystals.

## 6 Dislocation evolution

The density of forest obstacles depends directly on the dislocation densities in all slip systems of the crystal. Therefore, in order to close the model we require a equation of evolution for the dislocation densities. Processes resulting in changes in dislocation density include production by fixed sources, such as Frank-Read sources, breeding by double cross slip and pair annihilation (see  for a review; see also ). Although the operation of fixed Frank-Read sources is quickly eclipsed by production due to cross slip at finite temperatures, it is an important mechanisms at low temperatures. The double cross-slip, fixed Frank-Read sources and pair annihilation mechanisms are next considered in turn.

### 6.1 Breeding by double cross-slip

The importance of breeding by double cross-slip as a dislocation generation mechanism in crystals was emphasized by Johnston and Gilman . In this mechanism, the screw section of a moving dislocation migrates to a parallel plane by double cross-slip, thus creating a pair of sessile segments which pin down the dislocation and act in a manner similar to a Frank-Read source.

The rate at which dislocation length is produced by this mechanism can be estimated as follows. Let $`N_{\mathrm{CS}}^\alpha `$ denote the number of dynamic sources per unit volume induced by double cross-slip on slip system $`\alpha `$. Additionally, let $`L^{\mathrm{sat}}`$ denote the dislocation length emitted by each source prior to saturation. The rate of dislocation production per unit volume is, therefore:

$$\dot{\rho }_{\mathrm{CS}}^\alpha =L^{\mathrm{sat}}\dot{N}_{\mathrm{CS}}^\alpha $$
(53)
Assuming that cross-slip is thermally activated (see, e. g.,  for recent calculations of pathways and energy barriers in copper), the rate $`\dot{N}_{\mathrm{CS}}^\alpha `$ at which dynamic sources are generated may in turn be computed as follows. Let $`\overline{L}`$ be the mean-free path between cross-slip events, $`L^{\mathrm{cross}}`$ the length of the screw segment effecting the double cross glide, and $`E^{\mathrm{cross}}`$ the energy barrier for cross slip. The frequency of cross glide attempts per unit volume is, therefore, $`\rho ^\alpha v^\alpha /(\overline{L}L^{\mathrm{cross}})`$, where $`v^\alpha `$ is the mean dislocation velocity. According to transition-state theory, of these attempts a fraction $`\mathrm{e}^{\beta E^{\mathrm{cross}}}`$ is successful. This gives:

$$\dot{N}_{\mathrm{CS}}^\alpha =\frac{\rho ^\alpha v^\alpha }{\overline{L}L^{\mathrm{cross}}}\mathrm{e}^{\beta E^{\mathrm{cross}}}$$
(54)
Indeed, theoretical  and experimental  investigations suggest that the breeding rate due to cross-glide is proportional to the mean dislocation speed. Using Orowan’s formula, (54) may be rewritten as

$$\dot{N}_{\mathrm{CS}}^\alpha =\frac{\dot{\gamma }^\alpha }{b\overline{L}L^{\mathrm{cross}}}\mathrm{e}^{\beta E^{\mathrm{cross}}}$$
(55)
Inserting this relation into (53) gives

$$\dot{\rho }_{\mathrm{CS}}^\alpha =\frac{L^{\mathrm{sat}}\dot{\gamma }^\alpha }{b\overline{L}L^{\mathrm{cross}}}\mathrm{e}^{\beta E^{\mathrm{cross}}}$$
(56)
Following Essmann and Rapp , we shall additionally assume that the mean free path is inversely proportional to the dislocation density, which gives:

$$b\dot{\rho }_{\mathrm{CS}}^\alpha =\lambda _{\mathrm{CS}}\sqrt{\rho ^\alpha }\dot{\gamma }^\alpha $$
(57)
where we have written

$$\lambda _{\mathrm{CS}}=\frac{L^{\mathrm{sat}}}{L^{\mathrm{cross}}}\mathrm{e}^{\beta E^{\mathrm{cross}}}$$
(58)
An improvement on this model is to account for the resolved shear stress acting on the cross-slip plane in the computation of the activation energy, resulting in a so-called non-Schmid effect, but this possibility will not be pursued here in the interest of simplicity.

The double cross-slip mechanism is geometrically similar to the double kink formation alluded to earlier. By virtue of this similarity, we may expect that $`L^{\mathrm{cross}}L^{\mathrm{kink}}`$ and $`E^{\mathrm{cross}}E^{\mathrm{kink}}`$ to a first approximation.

### 6.2 Multiplication by fixed Frank-Read sources

An infinite amount of dislocation multiplication can be sustained by the fixed Frank-Read sources as long as newly formed loops can expand and move away from the pinning points. The precise details of this mechanism, which was independently proposed by Frank and Read in the early 50, are well documented and can be consulted elsewhere. In here, we limit our analysis to obtained an estimate for the dislocation multiplication based on this mechanism. Considering that the obstacles in the slip plane $`\alpha `$ serve as pinning points to operate the Frank-Read sources, the increase of the dislocation population is given by

$$\dot{\rho }_{\mathrm{FR}}^\alpha =L_{\mathrm{FR}}\dot{N}_{\mathrm{FR}}^\alpha $$
(59)
where $`N_{\mathrm{FR}}^\alpha `$ is the number of active sources per unit volume and $`L_{\mathrm{FR}}^\alpha `$ is the length of the emitted segment by the source and $`n^\alpha `$ is the obstacle density given by Eq. (25). The effective rate of activation of Frank-Read sources in a given system $`\alpha `$ can be estimated by the number of intersections of the gliding loops with obstacles in that system weighted by the efficiency of the source. In order to account for the possibility of the newly formed loops to glide away from the obstacles, we introduce the factor $`\sqrt{\rho ^\alpha /n^\alpha }`$ which is the ratio between mean distance between the obstacles and the mean dislocation spacing. Then,

$$\dot{N}_{\mathrm{FR}}^\alpha =\sqrt{\frac{\rho ^\alpha }{n^\alpha }}n^\alpha v^\alpha \rho ^\alpha .$$
(60)
Introducing Orowan’s formula and assuming that the length of dislocation segments is proportional to the mean distance between obstacles, the multiplication rate due to fixed Frank-Read sources can be written as

$$b\dot{\rho }_{\mathrm{FR}}^\alpha =\lambda _{\mathrm{FR}}\sqrt{\rho ^\alpha }\dot{\gamma }^\alpha $$
(61)
Unlike breeding by cross-glide, dislocation multiplication by fixed Frank-Read sources is not thermally activated process, and thus, it remains operative even at low temperatures. At finite temperatures however both mechanisms contribute to the total dislocation multiplication rate.

It should pointed out that both Eqs. (57) and (61) predict that $`\rho ^\alpha `$ grows as $`(\gamma ^\alpha )^2`$. This rate of growth is indeed observed in many crystals during the intermediate stages of dislocation multiplication.

### 6.3 Attrition by pair annihilation

Dislocation densities are often observed to attain a saturation density at sufficiently large strains . This saturation stage arises as a result of the competition between dislocation multiplication mechanisms such as double cross-slip and pair annihilation. For instance, spontaneous annihilation is observed in metals such as copper for screw dipole heights less than 1 nm . Pair annihilation is mainly the result of the cross-slip of screw segments of opposite sign .

Huang et al. , have studied the dynamic stability of short-range linear-elastic interactions between two dislocations of parallel line vectors which glide on two parallel slip planes in bcc crystals. Here we develop a similar but somewhat simpler linear-elastic model of dislocation pair annihilation. Thus, our goal is to estimate the frequency with which two parallel screw segments moving on parallel planes will acquire converging trajectories leading to their mutual annihilation. We note that, for simplicity, our analysis is restricted to annihilation between pairs of dislocations belonging to the same slip system.

Consider a screw segment moving under the action of an applied shear stress $`\tau `$. The segment follows a path which brings it in close proximity to a second immobile screw segment, Fig. 7. Let $`(x,y)`$ be coordinates centered at the fixed dislocation such that $`x`$ points in the direction of motion of the incoming dislocation and $`y`$ points in the direction of the cross-slip plane, Fig. 7. The interaction force per unit length exerted on the moving segment along the cross slip plane is :

$$f_y=K\frac{y}{r^2}$$
(62)
where $`K`$ is the pre-logarithmic factor for a screw segment, and

$$r=\sqrt{x^2+y^2+2xy\mathrm{cos}\alpha }$$
(63)
where $`\alpha `$ is the angle subtended by the $`x`$ and $`y`$ directions, Fig. 7. For bcc crystals, $`\alpha =2\pi /3`$. For isotropic crystals, $`K=\mu b^2/4\pi `$. In view of Eq. (11), the equation of motion for the incoming screw dislocation is:

$`\dot{y}`$ $`=`$ $`2l_P\nu _D\mathrm{e}^{\beta E^{\mathrm{cross}}}\mathrm{sinh}\left[\beta E^{\mathrm{cross}}{\displaystyle \frac{(Ky/br^2)}{\tau _0}}\right]`$ (64)
for $`y>\kappa _0>0`$, where the characteristic distance

$$\kappa _0=\frac{K}{b\tau _0}$$
(65)
is determined by the condition $`f_y=b\tau _0`$. For $`y<\kappa _0`$ dislocation motion is no longer a thermally activated process and therefore Eq. 64, which is based on the transition-state theory, does no longer apply. In this regime, where the interaction force exceeds the Peierls barrier, dislocation mobility increases significantly.

Depending on the value of $`y`$ as $`x\mathrm{}`$, the moving segment bypasses or is captured by the second segment. The precise calculation of the trajectories of the moving dislocation requires numerical computation. However, the essence of the annihilation mechanism can be captured by the following simple argument. Consider a dislocation initially at rest at position $`(x=0,y_0)`$. Subsequently, the motion of the dislocation for $`y\kappa _0`$ is governed by the equation:

$$\dot{y}=2l_P\nu _D\mathrm{e}^{\beta E^{\mathrm{cross}}}\mathrm{sinh}\left(\beta E^{\mathrm{cross}}\frac{(K/by)}{\tau _0}\right)$$
(66)
This equation is separable, and the time required for the annihilation process follows simply as:

$$t=_{\kappa _0}^{y_0}\left\{2l_P\nu _D\mathrm{e}^{\beta E^{\mathrm{cross}}}\mathrm{sinh}\left(\beta E^{\mathrm{cross}}\frac{(K/by)}{\tau _0}\right)\right\}^1𝑑y$$
(67)
In order to facilitate this calculation, the time required to travel from $`y=\kappa _0`$ to $`y=0`$ is neglected in the previous equation due to the much higher dislocation mobility in this regime and the small value of $`\kappa _0`$, of the order of few $`b`$ for BCC crystals. Additionally, we may simply set $`\mathrm{sinh}(x)x`$ to a first approximation, with the result:

$$t=\frac{\mathrm{e}^{\beta E^{\mathrm{cross}}}}{2l_P\nu _D\beta E^{\mathrm{cross}}}\frac{y_0^2\kappa _0^2}{2\kappa _0}$$
(68)
On the other hand, the time which the incoming dislocation spends at distances of the order of $`y_0`$ to the receiving dislocation is

$$t\frac{b\rho ^\alpha }{\dot{\gamma }^\alpha }y_0$$
(69)
For annihilation to be possible the annihilation time (68) must be less than the time (69) which the incoming dislocation spends in the immediate vicinity of the target dislocation. This yields the condition:

$$y_0\kappa =\kappa _0\left(A+\sqrt{A^2+1}\right)$$
(70)
where

$$A=\mathrm{e}^{\beta E^{\mathrm{cross}}}\beta E^{\mathrm{cross}}\dot{\gamma }_0^{\mathrm{cross}}/\dot{\gamma }^\alpha $$
(71)
is a factor depending on the strain rate and temperature,

$$\dot{\gamma }_0^{\mathrm{cross}}=2b\rho l_P\nu _D$$
(72)
is a reference slip-strain rate, and $`\kappa `$ may now be regarded as an effective pair annihilation distance. The cut-off value $`\kappa _c`$, corresponding to the maximal pair annihilation distance, is the effective screening distance which can be set equal to the mean distance between dislocations.

A simple expression which interpolates between the extreme values of $`\kappa `$ is

$$\frac{1}{\kappa }=\frac{1}{\kappa _c}+\frac{1}{\kappa _0\left(A+\sqrt{A^2+1}\right)}$$
(73)
It follows that the critical pair-annihilation distance $`\kappa `$ decreases with increasing strain rate and decreasing temperature. Thus, at high strain rates the dislocation velocities are high and the probability of being captured by another dislocation diminishes accordingly. Additionally, an increase in temperature increases the dislocation mobility and speeds up the annihilation process, which results in an attendant increase in annihilation rates. As will be demonstrated in Section 7, these trends are important in order to capture the temperature and strain-rate dependence of the stage I-II transition in Ta.

The rate of dislocation attrition due to pair annihilation may finally be estimated as follows. Over a differential of time $`dt`$, the number of annihilation events per unit volume is

$$dN^\alpha =\frac{1}{2L^{\mathrm{cross}}}(\rho ^\alpha )^2\kappa v^\alpha dt$$
(74)
The attendant loss of dislocation length per unit volume is

$$d\rho ^\alpha =2L^{\mathrm{cross}}dN^\alpha $$
(75)
Combining these relations and using Orowan’s formula, the rate of pair annihilation may finally be expressed in the form

$$b\dot{\rho }^\alpha =\kappa \rho ^\alpha \dot{\gamma }^\alpha $$
(76)
### 6.4 Dislocation multiplication rate

Combining (57), (61)) and (76), the total rate of change of the dislocation density follows as

$$b\dot{\rho }^\alpha =\lambda \sqrt{\rho ^\alpha }\dot{\gamma }^\alpha \kappa \rho ^\alpha \dot{\gamma }^\alpha $$
(77)
where $`\lambda =\lambda _{\mathrm{CS}}+\lambda _{\mathrm{FR}}`$. Evidently, $`\dot{\rho }^\alpha =0`$ upon the attainment of the saturation density:

$$\rho ^{\mathrm{sat}}=\left(\frac{\lambda }{\kappa }\right)^2$$
(78)
Dividing through by $`\rho ^{\mathrm{sat}}`$, (78) may be recast in the form:

$$\frac{\dot{\rho }^\alpha }{\rho ^{\mathrm{sat}}}=\left(\sqrt{\frac{\rho ^\alpha }{\rho ^{\mathrm{sat}}}}\frac{\rho ^\alpha }{\rho ^{\mathrm{sat}}}\right)\frac{\dot{\gamma }^\alpha }{\gamma ^{\mathrm{sat}}}$$
(79)
where

$$\gamma ^{\mathrm{sat}}=\frac{b}{\kappa }$$
(80)
is a saturation slip strain. Equation (79) can be integrated, yielding:

$$\frac{\rho ^\alpha }{\rho ^{\mathrm{sat}}}=\left[1+\left(\sqrt{\frac{\rho _0^\alpha }{\rho ^{\mathrm{sat}}}}1\right)\mathrm{exp}\left(\frac{1}{2}\frac{\gamma ^\alpha }{\gamma ^{\mathrm{sat}}}\right)\right]^2$$
(81)
The rate equation (79) expresses a competition between the dislocation multiplication and annihilation mechanisms. For slip strains $`\gamma ^\alpha \gamma ^{\mathrm{sat}}`$, the multiplication term dominates and, as noted previously, the dislocation density $`\dot{\rho }^\alpha `$ grows as $`(\dot{\gamma }^\alpha )^2`$. By contrast, when $`\gamma ^\alpha \gamma ^{\mathrm{sat}}`$ the rates of multiplication and annihilation balance out and saturation sets in. After saturation is attained, the dislocation density remains essentially unchanged. It should be carefully noted, that, in view of (73), the saturation slip strain $`\gamma _{\mathrm{sat}}`$ is a function of temperature and strain rate. In particular, $`\gamma _{\mathrm{sat}}`$ decreases with increasing temperature and decreasing strain rate. Since the stage I-II transition strain scales roughly with $`\gamma _{\mathrm{sat}}`$, we expect these trends to be exhibited by the transition strain itself, in accordance with experimental observation.

## 7 Comparison with experiment

We proceed to validate the theory against the uniaxial tests on Ta single crystals of Mitchell and Spitzig . In these tests, 99.97%-pure Ta specimens were loaded in tension along the $`[213]`$ crystallographic axis, at various combinations of temperature and strain rate. In particular we considered temperatures ranging from $`296K`$ to $`573K`$, and strain rates ranging from $`10^1s^1`$ to $`10^5s^1`$. The numerical procedure employed for the integration of the constitutive equations has been described elsewhere . The constitutive update is fully implicit, with the active systems determined iteratively so as to minimize an incremental work function. All stress-strain curves are reported in terms of nominal stress and engineering strain.

The material property set used in calculations is collected in Table 3. The elastic moduli $`C_{11}`$, $`C_{12}`$ and $`C_{44}`$ were obtained by fitting to the tables of Simmons and Wang . The ratio between edge and screw dislocation-line energies ($`U^{\text{edge}}/U^{\text{screw}}`$) is taken from the atomistic calculations of Wang et al. . The remaining parameters have been obtained by fitting to the experimental data of Mitchell and Spitzig . Note that we do not give any value for $`L^{\mathrm{sat}}`$, since it has no measurable influence in the range of temperatures considered here. It is expected, however, to play a role at higher temperatures.

Figs. 8 and 9 show the predicted and measured stress-strain curves for a $`[213]`$ Ta crystal over a range of temperatures and strain rates. It is evident from these figures that the model captures salient features of the behavior of Ta crystals such as: the dependence of the initial yield point on temperature and strain rate; the presence of a marked stage I of easy glide, specially at low temperature and high strain rates; the sharp onset of stage II hardening and its tendency to shift towards lower strains, and eventually disappear, as the temperature increases or the strain rate decreases; the parabolic stage II hardening at low strain rates or high temperatures; the stage II softening at high strain rates or low temperatures; the trend towards saturation at high strains; and the temperature and strain-rate dependence of the saturation stress.

The theory reveals useful insights into the mechanisms underlying these behaviors. For instance, since during state I the crystal deforms in single slip and the secondary dislocation densities are low, the Peierls resistance dominates and the temperature and strain-rate dependency of yield owes mainly to the thermally activated formation of kinks and crossing of forest dislocations. It is interesting to note that during this stage the effect of increasing (decreasing) temperature is similar to the effect of decreasing (increasing) strain rate, as noted by Tang et al. . The onset of stage II is due to the activation of secondary systems. The rate at which these secondary systems harden during stage I depends on the rate of dislocation multiplication in the primary system. This rate is in turn sensitive to the saturation strain $`\gamma ^{\mathrm{sat}}`$, which increases with strain rate and decreases with temperature. As a result, the length of the stage I of hardening is predicted to increase with strain rate and decrease with temperature, as observed experimentally. Finally, the saturation stress is mainly governed by the forest hardening mechanism and, in particular, by the strength of the forest obstacles. This process is less thermally activated than the Peierls stress, since the corresponding energy barriers are comparatively higher. Consequently, the stress-strain curves tend to converge in this regime, in keeping with observation.

The apparent softening observed in simulation results at the lowest temperature ($`296K`$) and the highest strain rate ($`10^1s^1`$) is actually an effect of the boundary conditions, which allow for some level of rotation of the specimen. Since in those cases, the material hardening is relatively low (stage I only), this geometrical softening dominates the apparent macroscopic behavior. In the other cases, the activation of several systems limits the extend of the rotations reducing the effects of the macroscopic hardening. In order to simulate more precisely the experimental boundary conditions, a model of the entire specimen allowing for a non-homogeneous deformation field should be considered.

The effect of temperature and strain-rate on hardening is also illustrated by the evolution of slip activity in the primary and secondary slip systems, $`D1`$ and $`D4^{\prime \prime }`$, respectively. Figs. 10 and 11 show the evolution of slip strains and dislocation densities as a function of temperature. It is evident from these figures that dislocation density saturation in the primary system occurs earlier as the temperature is increased, resulting in the activation of the secondary system and in the onset of stage II at lower strains. A similar effect is observed when the strain rate is decreased, Figs. 12 and 13.

The effect of loading direction on the hardening rate is illustrated in Fig. 14. The figure shows the stress-strain curves obtained by loading the crystal in the ,  and  directions at a strain rate of $`10^3`$ s<sup>-1</sup> and a temperature of 373 K. Each of these loading directions results in the activation of a different set of slip systems. As may be observed in the figure, the higher initial yield stress in the  direction relative to the baseline  direction, and the initial negative hardening rate, are fairly well captured by the model. The experimental curve exhibits a subsequent upturn, most likely due to the activation of additional secondary systems, which is not captured by the model. Loading in the  direction results in the activation of a large number of systems from the outset and the rate of hardening is correspondingly high. The model appears to over-predict the rate of hardening.

However, as was noted in , the stress-strain curve for crystals loaded in high-symmetry orientations are extremely sensitive to small misalignments in the loading axis, which accounts for the large experimental scatter characteristic of those orientations. This extreme sensitivity is due to the fact that small deviations from a high-symmetry loading axis, of the order of a degree or less, may result in the activation of a different set of slip systems, which may in turn have a large effect on the hardening rates. This pathological behavior is illustrated in Fig. 15, which shows the stress-strain curve obtained by randomly offsetting less than one degree the nominal loading axis . The resulting reduction in the hardening rate is quite remarkable. In addition, the experimental curve falls between the bounds of the stress-strain curves for the nominal and perturbed directions. This comparison exemplifies the need to take experimental scatter into account when assessing the fidelity of models, specially where crystals loaded along directions of high symmetry are concerned.

## 8 Summary and conclusions

We have developed a micromechanical model of the hardening, rate-sensitivity and thermal softening of bcc crystals. The model is predicated upon the consideration of an ‘irreducible’ set of unit processes, consisting of: double-kink formation and thermally activated motion of kinks; the close-range interactions between primary and forest dislocation, leading to the formation of jogs; the percolation motion of dislocations through a random array of forest dislocations introducing short-range obstacles of different strengths; dislocation multiplication due to breeding by double cross-slip; and dislocation pair-annihilation. Each of these processes accounts for–and is needed for matching–salient and clearly recognizable features of the experimental record. In particular, on the basis of detailed comparisons with the experimental data of Mitchell and Spitzig , the model is found to capture: the dependence of the initial yield point on temperature and strain rate; the presence of a marked stage I of easy glide, specially at low temperature and high strain rates; the sharp onset of stage II hardening and its tendency to shift towards lower strains as the temperature increases or the strain rate decreases; the initial parabolic hardening followed by saturation within the stage II of hardening; the temperature and strain-rate dependence of the saturation stress; and the orientation dependence of the hardening rates.

The choice of analysis tools which we have brought to bear on the unit processes of interest, e. g., transition-state theory, stochastic modeling, and simple linear-elastic models of defects and their interactions, is to a large extent conditioned by our desire to derive closed-form analytical expressions for all constitutive relations. As noted throughout the paper, many of the mechanisms under consideration are amenable to a more complete analysis by recourse to atomistic or continuum methods. However, at this stage of development, direct simulation methods, be it atomistic or continuum based, tend to produce unmanageable quantities of numerical data and rarely result in analytical descriptions of effective behavior. The daunting task of post-processing these data sets and uncovering patterns and laws within them which can be given analytical expression is as yet a largely unfulfilled goal of multiscale modeling.

This larger picture notwithstanding, one concrete and workable link between micromechanical models and first-principles calculations concerns the calculation of material constants. A partial list relevant to the present model includes: energy barriers and attempt frequencies for double-kink formation, kink migration, dislocation unpinning, cross-slip, and pair annihilation; dislocation-line and jog energies; and junction strengths. Other properties which have yielded to direct calculation include the volumetric equation of state (EoS), the pressure dependence of yield, and the pressure and temperature dependence of elastic moduli. References to recent work concerned with the calculation of these material properties have been given throughout the paper. As noted earlier, these results provide a suitable basis for future extensions of the present model to higher temperatures, pressures and strain-rates.

## Acknowledgments

The support of the DOE through Caltech’s ASCI Center for the Simulation of the Dynamic Response of Materials is gratefully acknowledged. LS also wishes to acknowledge support from the Belgian National Fund for Scientific Research (FNRS).

## List of Figures

Schematic of the double-kink mechanism.

Temperature dependence of the effective Peierls stress for various strain rates. Note that the typical order of magnitude of $`\dot{\gamma }_0^{\mathrm{kink}}=10^6s^1`$.

Bow-out mechanism for a dislocation segment bypassing an obstacle pair.

Schematic of energy variation as a function of a reaction coordinate during dislocation intersection and crossing.

Schematic of jog formation during dislocation intersection.

Schematic of a dislocation line overcoming a junction.

Interaction between two parallel screw dislocations of opposite sign moving on parallel slip planes.

Temperature dependence of stress-strain curves for $`[213]`$ Ta single crystal ($`\dot{ϵ}=10^3`$ s<sup>-1</sup>).(a) Experimental data of Mitchell & Spitzig . (b) Predictions of the model.

Strain-rate dependence of stress-strain curves for $`[213]`$ Ta single crystal ($`T`$=373 K). (a) Experimental data of Mitchell & Spitzig . (b) Predictions of the model.

Slip activity as a function of temperature ($`\dot{ϵ}=10^3`$ s<sup>-1</sup>).

Evolution of dislocation densities as a function of temperature ($`\dot{ϵ}=10^3`$ s<sup>-1</sup>).

Slip activity as a function of strain-rate ($`T`$=373 K).

Evolution of dislocation densities as a function of strain-rate ($`T`$=373 K).

Orientation dependence of stress-strain curves: comparison with experimental results (Mitchell & Spitzig ).

Effect of a small perturbation in the tensile orientation on the stress-strain curve.


format: text

text:
# Worm algorithms for classical statistical models

## Abstract

We show that high-temperature expansions may serve as a basis for the novel approach to efficient Monte Carlo simulations. “Worm” algorithms utilize the idea of updating closed path configurations (produced by high-temperature expansions) through the motion of end points of a disconnected path. An amazing result is that local, Metropolis-type schemes may have dynamical critical exponents close to zero (i.e., their efficiency is comparable to the best cluster methods). We demonstrate this by calculating finite size scaling of the autocorrelation time for various universality classes.

Metropolis scheme  is usually the most universal and easy to program approach to Monte Carlo simulations. However, its advantages are virtually canceled out near phase transition points. It is believed that any scheme based on local  Metropolis-type updates connecting system configurations into Markovian chain is inefficient at the transition point because its autocorrelation time, $`\tau `$, scales as $`L^z`$, where $`L`$ is the system linear dimension and $`z`$ is the dynamical critical exponent which is close to $`2`$ in most systems .

An enormous acceleration of simulations at the critical point has been achieved with the invention of cluster algorithms by Swendsen and Wang . However, the original method and its developments (both classical and quantum)  are essentially non-local schemes, and we are not aware of any exception from this rule.

In this Letter we propose a method which essentially eliminates the critical slowing down problem and yet remains local. The corner stone of our approach is the possibility to introduce the configuration space of closed paths. Closed-path (CP) configurations may be then sampled very efficiently using Worm algorithm (WA) introduced in Ref.  for quantum statistical models in which closed trajectories naturally arise from imaginary-time evolution of world lines. In classical models the CP representation derives from high-temperature expansions for a broad class of lattice models (see, e.g. Ref. ). In 2D, another family of WA may be introduced by considering domain-wall boundaries as paths.

We note, that our approach is based on principles which differ radically from cluster methods and, most probably, has another range of applicability. For one thing, the CP representation is most suitable for the study of superfluid models by having direct Monte Carlo estimators for the superfluid stiffness (through the histogram of winding numbers ) which are not available in the standard site representation.

In what follows we first recall how high-temperature expansions work by employing Ising model as an example (still, trying to keep notations as general as possible). We then explain how WA is used to update the path configuration space. Next, we discuss specific implementations of WA for $`|\psi |^4`$-, XY-, and $`q=3`$ Potts models, and comment on the special property of 2D models which allows an alternative CP parameterization of the configuration space. The efficiency of the new method is studied by looking at autocorrelation properties for six different universality classes. It is found that for 2D and 3D Ising models, 2D and 3D XY-models, and Gaussian model the $`\tau (L)`$ scaling is consistent with the law $`\tau (L)=\tau _0+c\mathrm{ln}(L)`$, i.e., its critical exponent is close to zero. For the two-dimensional $`q=3`$ Potts model our data are consistent with the power law with $`z0.55`$ which means that the Li-Sokal  bound $`z>\alpha /\nu `$ ($`\alpha `$ and $`\nu `$ are the specific heat and correlation radius critical exponents) derived for the Swendsen-Wang algorithm is seemingly applicable to our method as well.

Since high-temperature expansions for various models can be found in standard texts (see, e.g., Ref. ) we briefly remind the procedure for the Ising model

$$\frac{H}{T}=\beta \underset{b=<ij>}{}s_is_j,$$
(1)
where $`\beta =J/T`$ is the dimensionless nearest-neighbor coupling parameter between spin variables $`s_i=\pm 1`$ and index $`b=<ij>`$ refers to the simple cubic/square lattice bonds (we will also use another notation: $`b=(i,\nu )`$, in which $`\nu `$ enumerates bonds containing site $`i`$). Since $`H`$ is additive, the corresponding Gibbs exponent factories in terms of exponents for each bond. Expanding each exponent in Taylor series allows one to perform summation over site variables and to arrive at an expansion in powers of $`\beta `$. The partition function, for example, takes the form

$$Z=\underset{\{s_i\}}{}\underset{b=<ij>}{}\left(\underset{N_b=0}{\overset{\mathrm{}}{}}\frac{\beta ^{N_b}(s_is_j)^{N_b}}{N_b!}\right)=\underset{CP}{}W_{CP}.$$
(2)
Note that summation includes all possible CPs, both connected and disconnected, with self-intersections and overlaps. Graphically, each elementary factor in the order-by-order expansion over the bond Hamiltonian can be represented by a line ascribed to the corresponding bond , and the closed path constraint is imposed by symmetry: the sum $`_{s_i}s_i^{k_i}`$ is non-zero only if the number of bond lines connected to the site $`i`$ is even. For future purposes it is convenient to denote the number of bond lines $`N_b=0,1,2,\mathrm{}`$ and call it the “bond state”. Then the “site state” is $`k_i=_\nu N_{i,\nu }`$, and the configuration weight, $`W_{CP}`$, is given by

$$W_{CP}=\left(\underset{b}{}\frac{\beta ^{N_b}}{N_b!}\right)\left(\underset{i}{}Q(k_i)\right),Q(k)=\underset{s_i}{}s_i^k.$$
(3)
\[$`Q(k)2`$ for the Ising model.\]

The spin-spin correlation function $`G(i_1i_2)=s_{i_1}s_{i_2}`$ (we understand site indices as vectors) by definition equals to $`g(i_1i_2)/Z`$ where $`g(i_1i_2)=_{\{s_i\}}s_{i_1}s_{i_2}e^{H/T}`$. Configurations for $`g`$ differ from those contributing to the partition function in only one respect: there are two special sites $`i_1`$ and $`i_2`$ with odd number of bond lines connected to them, see Fig. 1. The corresponding site states are given by $`k_i=1+_\nu N_{i,\nu }`$. These points are the only places where the path may be disconnected. We will abbreviate the configuration space for $`g(i)`$ as CPg, and note, that weights $`W_{CP_g}`$ are given by the same expression (3).

The above consideration is actually a perfect setup for the Monte Carlo simulation since configurations for $`Z`$ and for $`g(0)`$ have identical bond elements (for the Ising model they even have equal weights, i.e. $`g(0)=Z`$). Thus, if an ergodic Metropolis process is sampling g-configurations according to their weights, then the Monte Carlo estimators for $`g(i)`$ and $`Z`$ are just

$$g(i|CP_g)=\delta _{i_2i_1,i},Z(CP_g)=\delta _{i_2i_1,0}.$$
(4)
The spin susceptibility is given by $`\chi =_ig(i)/Z`$, and energy may be computed in two ways: (i) as the nearest-neighbor sum $`E=(L^d/2)\beta _{|i|=1}g(i)/Z`$, and (ii) using direct estimator, $`E=/Z`$, which is non-zero only for the CP configurations contributing to Z:

$$(CP)=\beta \underset{b}{}N_b.$$
(5)
Our results for $`\tau `$ were obtained using this estimator.

The WA emerges as an idea to update $`g`$-configurations through the space motion of the end points $`i_1`$ and $`i_2`$ only . The algorithm itself is extremely simple and consists of two elementary updates which are proposed in the context of a given configuration: if $`i_1=i_2`$, then with probability $`p_0`$ we suggest to “move” both end points to site $`i_0`$ selected at random among all $`L^d`$ lattice cites, and with probability $`p_1=1p_0`$ to “shift” the end point $`i_1`$ to one of the neighboring sites by selecting the direction of the shift at random; if $`i_1i_2`$ we always suggest to “shift” $`i_1`$.

For the Ising model the move-update is automatically self-balanced and its acceptance ratio is unity, otherwise it would involve the ratio of $`Q`$-functions

$$P_{\mathrm{mv}}(ij)=\frac{Q(k_j+2)}{Q(k_j)}\frac{Q(k_i2)}{Q(k_i)}.$$
(6)
In contrast, shift-updates also increase/decrease the bond state by “drawing/erasing” bond lines, and we select which way to proceed with probability $`1/2`$. The shifts $`ij`$ with $`N_{b=<ij>}N_b+1`$ are balanced by shifts $`ji`$ with $`N_bN_b1`$. By standard rules  we have acceptance ratios as

$`P_{\mathrm{sh}}(ij,N_bN_b+1)`$ $`=`$ $`r{\displaystyle \frac{\beta }{(N_b+1)}}{\displaystyle \frac{Q(k_j+2)}{Q(k_j)}},`$ (7)
$`P_{\mathrm{sh}}(ij,N_bN_b1)`$ $`=`$ $`r{\displaystyle \frac{N_b}{\beta }}{\displaystyle \frac{Q(k_i2)}{Q(k_i)}},`$ (8)
where

$`1/(2p_1)`$ $`\text{ if }i_1=i_2=i,`$ (9)
$`r=`$ $`2p_1`$ $`\text{ if }i_1=i,i_2=j,`$ (10)
$`1`$  otherwise (11)
The Ising model has a special property which can be used in practical applications to enhance the efficiency by truncating the configuration space. Namely, using identity $`e^{\beta s_is_j}\mathrm{cosh}\beta [1+\mathrm{tanh}\beta s_is_j]`$ one can restrict the bond summation to include terms $`N_b=0,1`$ only. The corresponding modifications of the scheme are straightforward.

We now turn to the $`|\psi |^4`$-model described by the Hamiltonian ($`\psi _i`$ is a complex variable)

$$\frac{H}{T}=\underset{<ij>}{}(\psi _i\psi _j^{}+c.c.)+\underset{i}{}[\mu |\psi _i|^2U|\psi _i|^4],$$
(12)
the two limiting cases of which are the Gaussian model, $`U=0`$, and the XY-model, $`U=\mathrm{}`$. The procedure of factorizing the partition function expansion in terms of the bond Hamiltonian is exactly as before. The only new ingredient which is not present in the Ising model, is that for each bond we have two different terms: $`\psi _i\psi _j^{}`$ and $`\psi _i^{}\psi _j`$, and therefore the expansion has to be performed for each of them separately. Graphically, this can be captured by drawing lines with arrows, thus specifying each term by the arrow direction. Correspondingly, the bond state is defined by two numbers $`(N_b^{(1)},N_b^{(2)})`$ which tell how many lines and in which direction go along this bond. Since integrals $`𝑑\psi e^{\mu |\psi |^2U|\psi |^4}\psi ^m(\psi ^{})^m^{}=\delta _{m,m^{}}Q(k=2m)`$ are non-zero only for $`m=m^{}`$, we conclude that the configuration space for $`Z`$ is defined on CP and site states are given by $`k_i=_\nu \left(N_{i,\nu }^{(1)}+N_{i,\nu }^{(2)}\right)=`$even. \[The $`Q(k)`$ integrals are easily tabulated prior to the simulation.\]

In close analogy with the previous case, the configuration space for the two-point correlation function $`g(i_1i_2)=_id\psi _i\psi _{i_1}\psi _{i_2}^{}e^{H/T}`$, is given by CPg containing two special sites with $`k_i=1+_\nu \left(N_{i,\nu }^{(1)}+N_{i,\nu }^{(2)}\right)`$. WA updates consist of “moves” which change the location of end points when $`i_1=i_2`$ and “shifts” of $`i_1`$ which increase/decrease numbers $`N_{i_1,\nu }^{(1)}`$ and $`N_{i_1,\nu }^{(2)}`$. Acceptance ratios are given by Eqs. (6-8) where $`N_b`$ is understood as one of the numbers describing the bond state. An estimator for the winding number (which has exactly the same meaning as in worldline quantum Monte Carlo ) of the CP configuration is given by the difference between the bond numbers

$$M_\alpha =L^1\underset{b_\alpha }{}\left(N_{b_\alpha }^{(1)}N_{b_\alpha }^{(2)}\right),$$
(13)
where $`\alpha =x,y,z,\mathrm{}`$, and notation $`b_\alpha `$ is used to specify bonds connecting sites in direction $`\alpha `$.

The $`q=3`$ Potts model (which we considered to compare with the exhaustive cluster method study by Li and Sokal ) is described by

$$\frac{H}{T}=\beta \underset{<ij>}{}(\delta _{s_i,s_j}1),s_i=0,1,2.$$
(14)
Introducing phase variable $`\phi _i=(2\pi /3)s_i`$ we can rewrite Eq. (14) (up to a constant term) as

$$\frac{H}{T}=\beta ^{}\underset{<ij>}{}(e^{i(\phi _i\phi _j)}+c.c.),\beta ^{}=\beta /3.$$
(15)
With this form one may immediately proceed with the WA along the lines described above for the $`|\psi |^4`$-model. However, more efficient scheme results from the identity $`e^{2\beta ^{}\mathrm{cos}(\phi _i\phi _j)}=A[1+\gamma \mathrm{cos}(\phi _i\phi _j)]`$, where $`A=(e^{2\beta ^{}}+2e^\beta ^{})/3`$ and $`\gamma =2(e^{2\beta ^{}}e^\beta ^{})/3A`$, which allows to restrict summation over bond states to just three values $`N_b=(0,0),(1,0),(0,1)`$.

We would like to comment that high-temperature expansions are not the only procedure to arrive at the CP configuration space and WA. To underline this point we note that in lattice models with discrete site variables one can unambiguously specify the state of the system (up to symmetry transformations) by drawing domain boundaries, which, in 2D, may be considered as a CP configuration. For the Ising model the configuration weight is simply $`W_{CP}=_be^{\beta (2N_b1)}`$ where the bond state takes values $`N_b=0,1`$.

To implement WA for efficient sampling of domain boundaries we formally enlarge the configuration space to include CPg configurations with two end points and proceed exactly as discussed above. The only difference is that now “open” boundaries with $`i_1i_2`$ have no physical meaning and serve just for updating purposes. Our results for the 2D Ising model in Fig. 2 were obtained using domain-wall representation.

In Fig. 2 we present data for the energy autocorrelation time. In all cases, except $`q=3`$ Potts model in 2D, they are consistent with the linear in $`\mathrm{ln}L`$ behavior. Of course, we may not exclude that logarithmic scaling is actually an intermediate length-scale behavior and at some larger $`L`$ it crosses over to the power law, $`\tau (L)L^z`$, with small $`z`$. To get an estimate (upper bound) for the possible dynamic exponent we mention in the figure caption the slope $`z(L)=d\mathrm{ln}(\tau )/d\mathrm{ln}(L)`$ at the largest simulated $`L`$. For the $`q=3`$ Potts model our data scale as $`\tau L^{0.55}`$ for $`L>64`$, although the dynamic exponent is showing a systematic decrease with $`L`$. This has to be compared with the value $`z=0.515`$ for the Swendsen-Wong algorithm . It seems that the Li-Sokal lower bound for the dynamic exponent $`z>\alpha /\nu `$ also applies to our method. For all models discussed here $`\alpha /\nu `$ is very small with exception of the $`q=3`$ Potts model were $`\alpha /\nu =0.4`$.

To summarize, we found that WA Monte Carlo schemes working with closed path representations of classical statistical models are as efficient as the best cluster methods in reducing the problem of critical slowing down near transition points. The new approach is general enough to find applications other than discussed in this Letter, and may be used to study quantities like superfluid density for which other methods may not have direct estimators.

We thank J. Machta for valuable discussions. This work was supported by the National Science Foundation under Grant DMR-0071767.