text:
# The Wavelength Dependence of Polarization in Active Galactic Nuclei

## 1 Introduction

The degree of polarization for thermal radiation in a magnetic field was obtained by J.C.Kemp in1970 (for review see Lang 1974). It is given by the formula

$$p=eB/2mc\omega $$
(1)
where e and m are the charge and mass of electron respectively, c is the speed of light, B is the strength of magnetic field, and w is the circular frequency of radiation.

According to equation (1), the wavelength dependence of a degree of polarization is determined by the magnetic field profile B(r), where r is the distance from a central energy source, and by the relationship between the frequency of radiation and the radius r. The last relationship is given by the condition for emission (Prigara 2003). To introduce the magnetic field profile we consider the gaseous disk model for radio source.

## 2 The gaseous disk model

It was shown recently (Prigara 2003) that thermal radio emission has a stimulated character. According to this conception thermal radio emission of non-uniform gas is produced by an ensemble of individual emitters. Each of these emitters is a molecular resonator the size of which has an order of magnitude of mean free path l of photons

$$l=\frac{1}{n\sigma }$$
(2)
where n is the number density of particles and $`\sigma `$ is the absorption cross-section.

The emission of each molecular resonator is coherent, with the wavelength

$$\lambda =l,$$
(3)
and thermal radio emission of gaseous layer is incoherent sum of radiation produced by individual emitters.

The condition (3) implies that the radiation with the wavelength $`\lambda `$ is produced by the gaseous layer with the definite number density of particles n .

The condition (3) is consistent with the experimental results by Looney and Brown on the excitation of plasma waves by electron beam (see Chen 1987, Alexeev 2003). The wavelength of standing wave with the Langmuir frequency of oscillations depends on the density as predicted by equation (2). The discreet spectrum of oscillations is produced by the non-uniformity of plasma and the reajustment of the wavelength to the length of resonator. From the results of experiment by Looney and Brown the absorption cross-section for plasma can be evaluated.

The product of the wavelength by density is weakly increasing with the increase of density. This may imply the weak dependence of the size of elementary resonator in terms of the wavelength upon the density or, equivalently, wavelength.

In the gaseous disk model, describing radio emitting gas nebulae (Prigara 2003), the number density of particles decreases reciprocally with respect to the distance r from the energy center

$$nr^1.$$
(4)
Together with the condition for emission (3) the last equation leads to the wavelength dependence of radio source size:

$$r_\lambda \lambda .$$
(5)
The relation (5) is indeed observed for sufficiently extended radio sources. For example, the size of radio core of galaxy M31 is 3.5 arcmin at the frequency 408 MHz and 1 arcmin at the frequency 1407 MHz (Sharov 1982).

## 3 The wavelength dependence of radio source size

In the case of compact radio sources instead of the relationship (4) the relationship

$$r_\lambda \lambda ^2$$
(6)
is observed (Lo et al. 1993, Lo 1982). This relationship may be explained by the effect of a gravitational field on the motion of gas which changes the equation (3) for the equation

$$nr^{1/2}.$$
(7)
The mass conservation in an outflow or inflow of gas gives nvr=const, where v is the velocity of flow. In the gravitational field of a central energy source the energy conservation gives

$$v=\left(v_0^2+c^2r_s/r\right)^{1/2}$$
(8)
where $`r_s`$ is the Schwarzschild radius. Therefore, at small values of the radius the equation (7) is valid, whereas at the larger radii we obtain the equation (4). The relationship between linear size and turnover frequency in gigahertz-peaked spectrum sources and steep-spectrum sources (Nagar et al. 2002) is a consequence of the wavelength dependence of radio source size. The turnover frequency is determined by the equation $`r_\nu =R`$, where R is the radius of a gaseous disk. The same equation determines a turnover frequency for planetary nebulae (Pottasch 1984; Siodmiak & Tylenda 2001).

To summarize, extended radio sources are characterized by the relation (5), and compact radio sources obey the relation (6).

## 4 The wavelength dependence of polarization

At first, consider extended radio sources for which the equation (5) is valid. Assuming the magnetic field profile $`Br^2`$, which is characteristic for the field of two-dimensional magnetic dipole, from equation (1) we obtain

$$pB\left(r_\lambda \right)\lambda \lambda ^1.$$
(9)
This relationship is qualitatively consistent with observational data. The degree of linear and circular polarization in active galactic nuclei typically decreases with the increase of wavelength (Bower et al. 1999; Bower et al. 2001). In the case of Cen A the law (9) is valid in the decimeter range.

For compact radio sources the law (6) is valid, and the magnetic field profile is also changed to $`Br^{3/2}`$ (Prigara 2003). Now the equation (1) gives

$$p\lambda ^2.$$
(10)
The degree of linear polarization in 3C 84 is p=0.03% at the frequency 4.8 GHz and p¿0.2% at the frequency 14.5 GHz (Bower et al. 1999), in agreement with the equation (10). The degree of linear polarization in Cen A in centimeter range also follows the equation (10).

It is known that the synchrotron theory predicts a change in the polarization position angle across the spectral peak of gigahertz-peaked source. No such a change was found for the six gigahertz-peaked sources (Mutoh et al. 2002).

## 5 The strength of magnetic field

According to equation (1) the field B=20 G is required to produce the degree of polarization p=0.2% at the frequency 14.5 GHz, as in the case of 3C 84 (Bower et al. 1999). If we assume, for compact radio sources, that $`r_\lambda =r_0`$ at the wavelength $`\lambda `$=1 mm, where $`r_0`$ is the radius of the photosphere of a central energy source (e.g., the Schwarzschild radius), then $`B_0=B_\lambda \left(\lambda /\lambda _0\right)^3=2\times 10^5G`$. Zakharov et al. (2003) have obtained an estimate $`B<10^{10}10^{11}G`$ for the Seyfert galaxy MCG-6-30-15.

The existence of regular magnetic fields in active galactic nuclei is supported by the geometry of radio jets following the lines of a magnetic field. In many cases jets have a significant curvature, sometimes up to 90 degrees or more (Kellermann et al. 1998). The stable handedness of circular polarization in Sgr A\* and M81\* requires stable global magnetic field components (Beckert & Falcke 2002, Bower, Falcke & Mellon 2002). This suggests the existence of regular magnetic fields in accretion flows too.

One of the possible explanations for line widths in active galactic nuclei is a magnetic broadening. A magnetic broadening of spectral lines is produced by the Zeeman splitting in a non-uniform magnetic field changing its strength in the emitting region. The magnetic broadening in AGN is supported by broader $`H_\alpha `$ lines in polarized emission than in total emission (Nagar et al. 2002; Hes, Barthel & Fosbury 1993). However, the magnetic broadening of broad emission lines in AGNs encounters some difficulties. First, a change in the continuum flux produces an earlier response in the red wing of the line than in the blue wing (e.g., Konigl 2003). Second, some AGNs show double-peaked broad emission lines (Strateva et al. 2003). Both these features are indicative of a convection in the disk corona. In particular, the double-peaked emission lines are observed in planetary nebulae (Pottasch 1984). Nevertheless, a magnetic broadening can be applied to the narrow emission lines in AGNs. It may also participate in the width of each component of a double-peaked emission line.

A magnetic broadening may be a solution to the well-known problem of line broadening in A, B and O stars (Bohm-Vitense 1980). The Ca II lines in B stars usually assumed to be interstellar alternatively may be interpreted as a result of a non-uniform shift of spectral lines (Ebbets 1980) together with the Zeeman splitting. The latter produces the several (6 to 10) components of Ca II lines.

The phase diagram for the matter in ultrastrong magnetic fields is virtually unknown. The phase transitions in sufficiently strong magnetic fields may be responsible both for hard gamma-ray flares (Kellermann et al. 1998) and cosmic rays of ultrahigh energies. The origin of the latter is not elucidated up to date.

The effect of the strong magnetic field on the matter is produced by the magnetic pressure. In magneto-hydrodynamics it is assumed that the magnetic pressure is proportional to the $`B^2`$. One-fluid magneto-hydrodynamics is, however, poorly justified (Kadomtsev 1988). Alternatively, the magnetic pressure in a hot plasma may be suggested to be proportional to the strength of a magnetic field as follows

$$p_B=eBv{}_{0}{}^{}/\sigma c,$$
(11)
where $`\sigma `$ is the collisional cross-section in the high-temperature limit (when $`e^2/kT<h^2/Zme^2`$), and $`v_0`$ is a constant.

In the case of flat-spectrum active galactic nuclei (Bower & Backer 1998; Nagar, Wilson & Falcke 2001; Ulvestad & Ho 2001) the density, temperature, and pressure profiles have a form $`nr^{1/2},Tr^1,P=nkTr^{3/2}`$. Since for compact radio sources $`Br^{3/2}`$, from equation (11) we obtain $`PP_B`$. The latter relationship confirms equation (11). Here the density, temperature and pressure profiles are the same as those in the convection-dominated accretion flow (CDAF) models (Nagar et al. 2001). The magnetic field profile is different since the expression for the magnetic pressure is changed. The ratio of magnetic to gas pressure will be fixed similar to the CDAF models.

## 6 Conclusions

The maser theory of thermal radio emission (Prigara 2003) can explain the wavelength dependence of polarization in active galactic nuclei, if we only assume the existence of sufficiently strong magnetic fields in AGN. The other indications of strong magnetic fields in AGN are a magnetic broadening of spectral lines and possibly the hard gamma-ray flares.

## References

Alexeev B.V. (2003) Usp. Fiz. Nauk 173, 145.

Beckert T., Falcke H. (2002) Astron. Astrophys. 388, 1106.

Bohm-Vitense E. (1980), in Stellar Turbulence, IAU Colloq. 51, eds. D.F.Gray, J.L.Linsky, Springer.

Bower G.C., Backer D.C. (1998) Astrophys. J., Lett. 507, L117.

Bower G.C., Falcke H., Mellon R.R. (2002) Astrophys. J., Lett. 578, L103.

Bower G.C., Wright M.C.H., Backer D.C., Falcke H. (1999) Astrophys. J. 527, 851.

Bower G.C., Wright M.C.H., Falcke H., Backer D.C. (2001) Astrophys. J., Lett. 555, L103.

Chen F.F. (1984), Introduction to Plasma Physics and Controlled Fusion, Vol.1: Plasma Physics, Plenum Press.

Ebbets D. (1980) in Stellar Turbulence, IAU Colloq. 51, eds. D.F.Gray, J.L.Linsky, Springer.

Hes R., Barthel P.D., Fosbury R.A.E. (1993) Nature 361, 326.

Kadomtsev B.B. (1988), Cooperative phenomena in plasmas, Nauka.

Kellermann K.I., Vermeulen R.C., Zensus J.A., Cohen M.H. (1998) Astron. J. 115, 1295.

Konigl A. (2003) in Active Galactic Nuclei: from Central Engine to Host Galaxy, ASP Conf. Series 290, eds. S.Collin, F.Combes, I.Shlosman, ASP.

Lang K.R. (1974), Astrophysical Formulae, Springer.

Lo K.Y. (1982) in AIP Proc. 83: The Galactic Center, eds. G.Riegler, R.Blandford, AIP.

Lo K.Y., Backer D.C., Kellermann K.I., Reid M., Zhao J.H., Goss M.H., Moran J.M. (1993) Nature 361, 38.

Mutoh M., Inoue M., Kameno S., Asada K., Fujisawa K., Uchida Y. (2002), astro-ph/0201144.

Nagar N.M., Wilson A.S., Falcke H. (2001) Astrophys. J., Lett. 559, L87.

Nagar N.M., Wilson A.S., Falcke H., Ulvestad J.S., Mundell C.G. (2002) in Issues in Unification of AGNs, ASP Conf. Series 258, eds. R.Maiolino, A.Marconi, N.Nagar, ASP.

Pottasch S.R. (1984), Planetary Nebulae, D.Reidel Reinhold Company.

Prigara F.V. (2003), Astron. Nachr., Vol. 324, No. S1, Supplement: Proceedings of the Galactic Center Workshop 2002 - The central 300 parsecs of the Milky Way, eds. A.Cotera, H.Falcke, T.R.Geballe, S.Markoff.

Sharov A.S. (1982), The Andromeda Nebula, Nauka.

Siodmiak N., Tylenda R. (2001) Astron. Astrophys. 373, 1032.

Strateva I., Strauss M.A., Hao L. (2003) in Carnegie Observatories Astrophysics Series, Vol.1: Coevolution of Black Holes and Galaxies, ed. L.C.Ho, Carnegie Observatories.

Tadhunter C., Wills K., Morganti R., Oosterloo T., Dickson R. (2001) Mon. Not. R. Astron. Soc. 327, 227.

Ulvestad J.S., Ho L.C. (2001) Astrophys. J., Lett. 562, L133.

Zakharov A.F., Kardashev N.S., Lukash V.N., Repin S.V. (2003) Mon. Not. R. Astron. Soc. 342, 1325.


format: text

text:
# Plastic Flow in Two-Dimensional Solids

## I Introduction

Plastic flow has long been studied in crystalline and amorphous solids and in glassy polymers. In crystals irreversible motions of dislocations give rise to plastic deformation and large strains produce high-density dislocations . The nonlinear flow properties are very sensitive to the amount of such defects and strongly dependent on the deformation history. Simulations of dynamics of dislocation lines have recently been performed but are still most difficult .

In amorphous solids at low temperature $`T`$ , salient features are asfollows. (i) Shear strains tend to be localized in narrow shear bands in plastic flow above a yield stress. The width of such shear bands is microscopic in the initial stage  but can grow to mesoscopic sizes , sometimes resulting in fracture. Shear bands were numerically realized at large shear strains in molecular dynamics (MD) simulations of two-dimensional (2D) two-component glasses  and in simulations of a 2D phenomenological stochastic model . (ii) As another aspect, in 3D MD simulations on model two-component glasses at low $`T`$, Takeuchi et al.  observed heterogeneities among mobile and immobile regions after application of shear strains. In 2D and 3D MD simulations on model supercooled binary mixtures above $`T_\mathrm{g}`$, similar or much more extended dynamic heterogeneities have been detected in quiescent states  and found to be sensitively suppressed by applied shear flow . (iii) Furthermore, at higher $`T(0.70.8T_\mathrm{g})`$ (where $`T_\mathrm{g}`$ is the glass transition temperature), shear deformation occurs quasi-homogeneously (still involving many particles in each configuration change), leading to highly viscous non-Newtonian behavior .

Glassy polymers also behave analogously , where shear bands appear above a yield stress at low $`T`$ , and highly viscous non-Newtonian flow and significant elongation of the chain shapes occur at elevated $`T`$ . In glassy polymers, the entropic stress arising from molecular orientations becomes significant at large strains and such systems behave like cross-linked rubbers .

Recently much attention has been paid to jamming rheology observed in sheared states of supercooled liquids, soft glassy materials such as dense microemulsions, or granular materials . In these systems, the thermal agitation effect is very small if the particle size is large, but universal constrained dynamics is realized under external forces. In supercooled liquids (at relatively high $`T`$)  and dense microemulsions (at effectively low $`T`$) , mesoscopic dynamic heterogeneity and strong shear-thinning behavior have been observed, but shear bands have not been identified. In granular materials (at effectively zero temperature), strain localization is most conspicuous .

In our recent work  we have constructed 2D nonlinear strain theory taking into account the underlying local periodic lattice structure, where the elastic energy density is periodic with respect to the shear and tetragonal strains. Then in our model plastic flow starts with appearance of slips. The simplest slips consist of two edge dislocations having opposite Burgers vectors with size $`a`$ ($`a`$ being the lattice constant) and they grow into mesoscopic shear bands as the applied strain is increased. Under uniaxial stretching , well developed shear bands were already numerically realized  and will be realized also in our simulations. These shear bands make angles of $`\pm \pi /4`$ with respect to the stretched direction in agreement with observations in various amorphous materials  including granular materials . These angles will be shown to minimize the elastic energy of incipient slips in this paper.

In crystalline solids, dislocation pairs in 2D or dislocation loops in 3D forming slip lines or surfaces should be nucleated at the inception of plastic deformation (in addition to preexisting dislocations). In amorphous solids, it has been controversial on whether dislocations themselves can be well-defined or not . It is not obvious how to characterize the local rearrangement processes as to their shapes and sizes. In MD simulations on binary mixtures, they have been detected as cluster-like objects with various visualization methods . If the size-ratio of the constituent two species is chosen such that crystallization is most suppressed, strong frustration occurs in the packing of large and small particles in jammed states. Then the local crystalline order can be well-defined only over short distances. However, in a 2D amorphous soap bubble raft, Argon and Kuo  observed that nucleation of a dislocation pair gave rise to a small-scale slip but such dislocations did not glide more than 2-3 bubble distances. It is worth noting that Deng et al.  found extended slip-like strain localization in 2D MD simulations with the size-ratio rather close to 1 (see comment (iv) in the last section). Notice that the displacement field around a slip is localized because the two constituent dislocations have opposite Burgers vectors. As a result, small-scale slips should be well-defined even in amorphous solids as long as the slip size does not much exceed the range of the local crystal structure.

The plastic flow phenomena are thus very complex, being influenced by many factors, but they are ubiquitous in various kinds of solid-like materials. The purpose of this paper is to present a well-defined Ginzburg-Landau model consistently taking account of nonlinear elasticity. A merit of this approach is that we can put emphasis on any aspects of the phenomena by controlling the parameters or changing the model itself. We will examine (i) the fundamental flow units, slips, in detail numerically and analytically and (ii) plastic flow numerically in simple shear and elongational (stretching) deformation. To make this paper simplest, as it is the first detailed exposition of our scheme, we will neglect (i) vacancies and interstitials, or a variable $`m`$ representing the local free-volume. A dynamic model including such an additional degree of freedom has already been presented in our previous work . We will also neglect (ii) the configurational frustration effect induced by the size difference between the two species. Introducing these two ingredients will constiture future development of our scheme.

This paper is organized as follows. In Section 2 we will present our dynamic model and explain our numerical method. In Section 3 we will discuss the simplest form of slips numerically obtained from our nonlinear elasticity theory. We will also derive some analytic expressions for the slip formation energy under general strain field starting with the Peach-Koehler theory  and compare them with numerical results. In our scheme stationary slip solutions exist for small externally applied strains, where the force balance is achieved in the presence of the Peierls potential energy for the dislocation position in crystals . Section 4 will present numerical results of the stress-strain relations and the patterns of the strains, the elastic energy density, and the displacement vector under applied strains. We will also give a method of dividing the elastic energy into an elastic part due to affine deformation and a defect part.

## II Model equations

### A Elastic energy

Recently Doi et al.  described plastic flow in a highly viscous 2D crystal phase of block copolymers in the framework of a 1D Frenkel-Kontrowa model  assuming 1D sliding motions. Hereafter we present a 2D Frenkel-Kontrowa model to describe plastic flow in 2D . In terms of the displacement vector $`𝒖=(u_x,u_y)`$ from a reference crystal state, we define the strain components as

$`e_1`$ $`=`$ $`_xu_x+_yu_y,`$ (2.1)
$`e_2`$ $`=`$ $`_xu_x_yu_y,`$ (2.2)
$`e_3`$ $`=`$ $`_xu_y+_yu_x,`$ (2.3)
where $`_x=/x`$ and $`_y=/y`$. We call $`e_1`$ the dilation strain, $`e_2`$ the tetragonal strain, and $`e_3`$ the shear strain. If we suppose a 2D triangular lattice with lattice constant $`a`$, the elastic energy should be invariant with respect to the rotations of the reference frame by $`\pm n\pi /3`$ ($`n=1,2,\mathrm{}`$). Due to this symmetry, the elasticity must be isotropic in the harmonic approximation , being characterized by the bulk and shear moduli, $`K_0`$ and $`\mu _0`$, but it depends on the orientational angle $`\theta `$ of one of the crystal axes with respect to the $`x`$ axis for large shear strains. Under rotation of the reference frame by $`\theta `$, the shear strains $`e_2`$ and $`e_3`$ are changed to $`e_2^{}`$ and $`e_3^{}`$, where

$`e_2^{}`$ $`=`$ $`e_2\mathrm{cos}2\theta +e_3\mathrm{sin}2\theta ,`$ (2.4)
$`e_3^{}`$ $`=`$ $`e_3\mathrm{cos}2\theta e_2\mathrm{sin}2\theta .`$ (2.5)
These relations are obtained from the orthogonal transformations, $`𝒓^{}=\stackrel{}{U}𝒓`$ and $`𝒖^{}=\stackrel{}{U}𝒖`$, where $`𝒓^{}=(x^{},y^{})`$ and $`𝒖^{}=(u_x^{},u_y^{})`$ represent the position and the displacement, respectively, in the new reference frame, and $`\stackrel{}{U}=\{U_{ij}\}`$ with $`U_{xx}=U_{yy}=\mathrm{cos}\theta `$ and $`U_{xy}=U_{yx}=\mathrm{sin}\theta `$. The $`e_2^{}`$ and $`e_3^{}`$ are the tetragonal and shear strains, respectively, in the new reference frame where the $`x^{}`$ axis is along one of the crystal axes.

The elastic energy is written as $`F_{\mathrm{el}}=𝑑𝒓f_{\mathrm{el}}`$ with the elastic energy density in the form,

$$f_{\mathrm{el}}=\frac{1}{2}K_0e_1^2+\mu _0\mathrm{\Phi }(e_3^{},e_2^{}),$$
(2.6)
which is independent of the rotation strain,

$$\omega =_xu_y_yu_x.$$
(2.7)
Note that $`e_1`$ and $`\omega `$ are invariant with respect to the rotation of the reference frame. The simplest form of $`\mathrm{\Phi }`$ is given by

$`\mathrm{\Phi }(e_3^{},e_2^{})`$ $`=`$ $`{\displaystyle \frac{1}{6\pi ^2}}[3\mathrm{cos}\pi (\sqrt{3}e_3^{}e_2^{})`$ (2.9)
$`\mathrm{cos}\pi (\sqrt{3}e_3^{}+e_2^{})\mathrm{cos}(2\pi e_2^{})].`$
This function is invariant with respect to the rotation $`\theta \theta +\pi /3`$, is a periodic function of $`e_3^{}`$ with period $`2/\sqrt{3}`$ for $`e_2^{}=0`$ (simple shear deformation), and becomes $`(e_2^2+e_3^2)/2`$ for small strains. In Fig.1 we display $`\mathrm{\Phi }(e_3,e_2)`$ supposing $`\theta =0`$, where one of the crystal axes is along the $`x`$ axis. We can see a hexagonal lattice structure in the $`e_3`$-$`e_2`$ plane. As a characteristic feature, Fig.2 shows that it is almost isotropic or is a function of $`e=(e_3^2+e_2^2)^{1/2}`$ only even for rather large $`e`$ ($`0.5`$) (in the unit cell at the origin). In fact, we have the Taylor expansion,

$`\mathrm{\Phi }(e_3,e_2)`$ $`=`$ $`{\displaystyle \frac{1}{2}}e^2{\displaystyle \frac{\pi ^2}{8}}e^4+{\displaystyle \frac{\pi ^4}{72}}e^6`$ (2.11)
$`+{\displaystyle \frac{\pi ^4}{720}}(e_2^2e_3^2)(e^416e_2^2e_3^2)+O(e^8)`$
If we set $`e_3=e\mathrm{cos}\chi `$ and $`e_2=e\mathrm{sin}\chi `$, the fourth term is rewritten as $`(\pi ^4/720)e^6\mathrm{cos}(6\chi )`$ and is known to keep the invariance $`\chi \chi +\pi /3`$. Anisotropy appears in the terms of order $`e^6`$, but the anisotropic fourth term is at most $`10\%`$ of the isotropic third term.

We examine elastic stability of homogeneously strained states. That is, we superimpose infinitesimal strains, $`\delta e_3`$ and $`\delta e_2`$, on $`e_3`$ and $`e_2`$ assumed to be homogeneous. The second order terms in $`\mathrm{\Phi }(e_3+\delta e_3,e_2+\delta e_2)`$ read

$$\delta ^{(2)}\mathrm{\Phi }=\frac{\mathrm{\Phi }_{33}}{2}(\delta e_3)^2+\mathrm{\Phi }_{23}\delta e_3\delta e_2+\frac{\mathrm{\Phi }_{22}}{2}(\delta e_2)^2,$$
(2.12)
where $`\mathrm{\Phi }_{\alpha \beta }=^2\mathrm{\Phi }/e_\alpha e_\beta `$ ($`\alpha ,\beta =2,3`$). In the stable regions the above second order contribution should be positive-definite or the two eigenvalues, $`\lambda _1`$ and $`\lambda _2`$, of the $`2\times 2`$ matrix $`\{\mathrm{\Phi }_{\alpha \beta }\}`$ should be both positive. In the appendix we will derive this linear stability criterion by solving the linearized version of our dynamic model to be presented in the next section. Fig.3 shows that the system is stable for $`e0.3`$, where $`\mathrm{\Phi }`$ depends almost only on $`e`$. This result readily follows if we neglect the anisotropy in $`\mathrm{\Phi }`$ by setting $`\mathrm{\Phi }(e_2,e_3)G(e^2)`$. Then some calculations yield

$$\mathrm{\Phi }_{\alpha \beta }2G^{}\delta _{\alpha \beta }+4G^{\prime \prime }e_\alpha e_\beta ,$$
(2.13)
where $`G^{}=dG(e^2)/de^2`$ and $`G^{\prime \prime }=dG^{}(e^2)/de^2`$. From (2.6) we find $`G^{}1/2`$ and $`G^{\prime \prime }\pi ^2/4`$. The determinant of $`\{\mathrm{\Phi }_{\alpha \beta }\}`$ becomes $`4G^{}(G^{}+2G^{\prime \prime }e^2)`$. Thus the stability condition becomes

$$e<(G^{}/|2G^{\prime \prime }|)^{1/2}0.3,$$
(2.14)
around the origin in the $`e_3`$-$`e_2`$ plane. The elastic instability with a negative eigenvalue ($`\lambda _1<0`$ or $`\lambda _2<0`$) causes rapid relaxation processes resulting in localized slips and stress release. In plastic flow the stability condition is mostly satisfied throughout the system (see Fig.11).

If we assume that small solid elements are rotated without shape changes with the local angular velocity $`(v_y/xv_x/y)/2`$, the orientation angle $`\theta `$ is related to the rotation strain as

$$\theta =\frac{1}{2}\omega +\theta _0,$$
(2.15)
where $`\theta _0`$ is the initial value independent of time $`t`$. This relation was assumed in our previous paper . However, it becomes not well defined when the three crystal axes are rotated differently. For example, let a crystal with the three axes at $`\theta =0,\pi /3,`$ and $`\pi /3`$ be affinely deformed by a simple shear deformation given by $`u_x=\gamma y`$ and $`u_y=0`$ ($`e_3=\gamma `$ and $`e_2=0`$). Then the the first axis along the $`x`$ direction is not rotated, while the other axes are rotated by

$$\delta \theta _\pm =\pm \left[\mathrm{tan}^1\left(\frac{\sqrt{3}}{1\pm \sqrt{3}\gamma }\right)\frac{\pi }{3}\right].$$
(2.16)
For $`\gamma 1`$, we have $`\delta \theta _\pm 3\gamma /4`$, so (2.10) holds for the average of the three rotation angles, $`(0+\delta \theta _++\delta \theta _{})/3=\gamma /2=\omega /2`$. Thus there remains ambiguity in (2.10) particularly at large strains $`e0.5`$ . In this paper, in view of the virtual isotropic behavior of our elastic energy for $`e0.5`$, we will simply set

$$\theta =0,$$
(2.17)
throughout the system. To check the appropriateness of this assumption, we have also performed simulations with various homogeneous $`\theta `$ held fixed. The resultant patterns and stress-strain relations are insensitive to $`\theta `$ and very similar to those obtained for $`\theta =0`$ to follow.

### B Dynamic equations

The elastic stress tensor should be defined for general strains. Let us change an arbitrary displacement $`𝒖`$ by an infinitesimal amount $`\delta 𝒖`$ as $`𝒖𝒖+\delta 𝒖`$. The incremental change of the elastic energy density is written as

$$\delta f_{\mathrm{el}}=\underset{ij}{}\sigma _{ij}\frac{}{x_j}\delta u_i.$$
(2.18)
This is just the definition of the elastic stress tensor $`\sigma _{ij}`$. Under (2.10)  it is symmetric and its components are expressed as

$`\sigma _{xx}`$ $`=`$ $`K_0e_1+\mu _0A_2,\sigma _{yy}=K_0e_1\mu _0A_2,`$ (2.19)
$`\sigma _{xy}`$ $`=`$ $`\sigma _{yx}=\mu _0A_3,`$ (2.20)
where $`A_\alpha =\mathrm{\Phi }(e_3,e_2)/e_\alpha `$ are written as

$`A_2`$ $`=`$ $`{\displaystyle \frac{1}{3\pi }}\left[\mathrm{cos}(\sqrt{3}\pi e_3)\mathrm{sin}(\pi e_2)+\mathrm{sin}(2\pi e_2)\right],`$ (2.21)
$`A_3`$ $`=`$ $`{\displaystyle \frac{1}{\sqrt{3}\pi }}\mathrm{sin}(\sqrt{3}\pi e_3)\mathrm{cos}(\pi e_2).`$ (2.22)
For small strains we have $`A_2e_2`$ and $`A_3e_3`$ to reproduce the isotropic, linear elastic theory.

We assume that the lattice velocity

$$𝒗=\frac{}{t}𝒖$$
(2.23)
obeys

$$\rho \frac{}{t}𝒗=\stackrel{}{\sigma }++\eta _0^2v+\underset{\mathrm{R}}{\overset{}{\sigma }},$$
(2.24)
where we introduce the shear viscosity $`\eta _0`$ but neglect the bulk viscosity . The $`\underset{\mathrm{R}}{\overset{}{\sigma }}=\{\sigma _{ij}^\mathrm{R}\}`$ is a symmetric random stress tensor and satisfies $`\sigma _{xx}^\mathrm{R}+\sigma _{yy}^\mathrm{R}=0`$, because the bulk viscosity is neglected, and is related to $`\eta _0`$ by

$$\sigma _{ij}^\mathrm{R}(𝒓,t)\sigma _{ij}^\mathrm{R}(𝒓^{},t^{})=2k_\mathrm{B}T\eta _0\delta (𝒓𝒓^{})\delta (tt^{}).$$
(2.25)
The mass density $`\rho `$ will be treated as a constant in (2.17). This is justified when the deviation $`\delta \rho =\rho \rho `$ is assumed to be much smaller than the average $`\rho `$ . In the usual linear elasticity theory  it is related to the dilation strain $`e_1`$ as $`\delta \rho \rho e_1`$, so we are assuming $`|e_1|1`$ (and coincidence of the lattice and mass velocities ). In our simulation in the plastic flow regime at $`\dot{\gamma }=10^3`$, for example, $`|e_1|`$ attains a maximum in a range of $`0.2`$-$`0.3`$ and the variance $`\sqrt{e_1^2}`$ increases up to about $`0.06`$.

Due to the presence of the random stress, $`𝒖`$ and $`𝒗`$ are random variables and (2.16) and (2.17) constitute nonlinear Langevin equations . Their equilibrium distribution attained in the unstrained condition is given by $`P\mathrm{exp}(F/k_\mathrm{B}T)`$, where the total free energy $`F`$ is the sum of the elastic energy $`F_{\mathrm{el}}`$ and the kinetic energy as

$$F=𝑑𝒓\left(f_{\mathrm{el}}+\frac{\rho }{2}𝒗^2\right).$$
(2.26)
Furthermore, if the random stress is omitted in (2.17) and the dynamic equations are treated as deterministic ones, the time derivative of $`F`$ is nonnegative-definite in the unstrained condition as

$$\frac{d}{dt}F=𝑑𝒓\underset{ij}{}\eta _0(_iv_j)^20,$$
(2.27)
where use is made of the relation,

$$\stackrel{}{\sigma }=\frac{\delta }{\delta 𝒖}F.$$
(2.28)
The stationary condition $`dF/dt=0`$ is attained for $`𝒗=\mathrm{𝟎}`$ and $`\stackrel{}{\sigma }=\mathrm{𝟎}`$. This is a condition to guarantee self-consistency of the dynamic equations which have stable equilibrium solutions.

In the appendix we will examine the linearized dynamic equations of our model around homogeneously strained states to obtain two sound modes in the stable region.

### C Dimensionless forms

We make our equations dimensionless measuring space in units of the lattice constant $`a`$ and time in units of

$$\tau _0=(\rho /\mu _0)^{1/2}a.$$
(2.29)
The stress components and the elastic energy density are measured in units of $`\mu _0`$ and the elastic energy in units of $`\mu _0a^2`$, while the strains remain unscaled. To avoid to introduce too many symbols, we will rewrite the scaled position vector $`a^1𝒓`$, time $`\tau _0^1t`$, displacement vector $`a^1𝒖`$, and velocity $`\tau _0a^1𝒗`$ simply as $`𝒓`$, $`t`$, $`𝒖`$, and $`𝒗`$ using the original notation. Then, in the dimensionless $`\stackrel{}{\sigma }`$, $`\mu _0`$ is replaced by 1 and $`K_0`$ by the ratio,

$$\lambda =K_0/\mu _0.$$
(2.30)
In terms of the scaled quantities the equilibrium distribution is written as

$$P\mathrm{exp}\left[\frac{1}{ϵ_{\mathrm{th}}}𝑑𝒓\left(\frac{\lambda }{2}e_1^2+\mathrm{\Phi }(e_3,e_2)+\frac{1}{2}𝒗^2\right)\right],$$
(2.31)
where

$$ϵ_{\mathrm{th}}=k_\mathrm{B}T/\mu _0a^2.$$
(2.32)
The parameter $`ϵ_{\mathrm{th}}`$ represents the degree of the thermal fluctuations (being proportional to $`T`$) and is an important parameter, for example, in describing the decay of metastable states by thermal agitations. If the dynamic equation (2.17) is made dimensionless, the dimensionless viscosity is given by

$$\eta _0^{}=\eta _0a^1(\rho \mu _0)^{1/2}.$$
(2.33)
In (2.18) the noise strength $`k_\mathrm{B}T\eta _0`$ is replaced by $`ϵ_{\mathrm{th}}\eta _0^{}`$.

### D Numerical method

We integrate (2.16) and (2.17) in the dimensionless units on a $`128\times 128`$ square lattice represented by $`(n,m)`$ ($`1n,m128`$) with $`x=n\mathrm{\Delta }x`$ and $`y=m\mathrm{\Delta }x`$. The modulus ratio $`\lambda `$ in (2.23) is set equal to 1. For simplicity, the mesh size $`\mathrm{\Delta }x`$ and the dimensionless viscosity in (2.26) are set equal to 1:

$$\mathrm{\Delta }x=1,\eta _0^{}=1.$$
(2.34)
The first relation means that the mesh size is just equal to the lattice constant, and the second one is rewritten as $`\eta _0=a(\rho \mu _0)^{1/2}`$.

We give our boundary conditions employed. At the bottom $`y=0`$, we set $`𝒖=\mathrm{𝟎}`$. At the top $`y=L_0=128`$, we set $`u_x=\gamma L_0`$ and $`u_y=0`$ in the presence of applied shear strain $`\gamma `$, and we set $`u_x=u_y=ϵL_0/2`$ in the presence of applied tetragonal strain $`ϵ`$. In the $`x`$ direction, we impose the periodic boundary condition, $`𝒖(x+L_0,y,t)=𝒖(x,y,t)`$. The unstrained condition is realized for $`\gamma =0`$ and $`ϵ=0`$.

We are interested in slips across which the atomic displacement is discontinuous by the lattice constant $`a`$. In this paper, by setting $`\mathrm{\Delta }x=a`$, we try to realize such singular objects numerically. For this purpose it is convenient to define the strains and tensors on the middle points $`(n+1/2,m+1/2)`$, while the vectors @ are defined at the lattice points $`(n,m)`$. For a vector component $`A`$ (say, $`A=u_x`$), $`_xA`$ and $`_yA`$ at $`(n+1/2,m+1/2)`$ are defined as $`[A(n+1,m+1)A(n,m)+A(n+1,m)A(n,m+1)]/2`$ and $`[A(n+1,m+1)A(n,m)A(n+1,m)+A(n,m+1)]/2`$, respectively, using $`A`$ at the four points $`(n+1/2\pm 1/2,m+1/2\pm 1/2)`$. In the same manner, we may construct the vector $`\stackrel{}{\sigma }`$ at $`(n,m)`$ using the stress components at the four points $`(n\pm 1/2,m\pm 1/2)`$. With this space discretization, slips consisting of a straight line segment become well-defined if their angle with respect to the $`x`$ or $`y`$ axis is $`0`$ or $`\pi /4`$. For other slip orientations, zigzag points appear along the slip line segment and an extra elastic energy becomes needed. On the other hand, if we define all the quantities on $`(n,m)`$, slip discontinuity takes place over a few lattice sizes, but macroscopic features such as the stress-strain relation remain almost unchanged. Furthermore, we will suppose simple shear or uniaxial deformation and, as will be shown in the next section, the preferred slip orientation angle is $`0`$ or $`\pi /4`$ with respect to the $`x`$ or $`y`$ axis. By this reason our simple numerical scheme seems to be allowable at least in this first attempt.

## III Slips

In our model system fundamental flow units in plastic deformation are slips composed of a pair of edge dislocations with opposite Burger vectors (dislocation dipoles). They are analogous to quantum vortex rings in superfluid helium . Their elastic structure far from the dislocation cores may well be described by the linear elasticity theory, but nonlinear elasticity theory is needed (i) to suppress the divergence of the stress at the cores and (ii) to stabilize the slips themselves when they adjust to the crystal structure. Slips are not in a stationary state in the linear elasticity theory in the absence of impurities etc. which can trap dislocations, as will be evident in (3.20). In our nonlinear theory, those along the $`x`$ axis ($`\theta =0)`$ can be in a stationary metastable state if their length is a multiple of the lattice constant. This is consistent with the Peierls-Nabarro theory , which takes into account the discreteness of the crystal structure and gives a periodic Peierls potential energy for the position of the dislocation center.

### A Slips in linear elasticity theory

To begin with, let us write the solution of an edge dislocation as $`𝒖=b𝒖^{\mathrm{Le}}=(bu_x^{\mathrm{Le}},bu_y^{\mathrm{Le}})`$, whose Burgers vector is assumed to be along the $`x`$ direction and is written as $`𝒃=(b,0)`$. The linear elasticity theory  gives

$`u_x^{\mathrm{Le}}`$ $`=`$ $`{\displaystyle \frac{1}{2\pi }}\left[\mathrm{tan}^1\left({\displaystyle \frac{y}{x}}\right)+{\displaystyle \frac{1}{2(1\nu )}}{\displaystyle \frac{xy}{x^2+y^2}}\right]`$ (3.1)
$`u_y^{\mathrm{Le}}`$ $`=`$ $`{\displaystyle \frac{1}{4\pi (1\nu )}}\left[(12\nu )\mathrm{ln}\sqrt{x^2+y^2}+{\displaystyle \frac{x^2}{x^2+y^2}}\right].`$ (3.2)
In the 3D theory the dislocation line is along the $`z`$ axis and $`\nu `$ is 3D Poisson’s ratio. In our 2D theory $`\nu `$ is related to $`\lambda `$ in (2.23) as

$$\nu =\frac{1}{2}\frac{1}{2\lambda }.$$
(3.3)
In our simulations $`\lambda =1`$ and hence $`\nu =0`$. In the linear theory a slip is a superposition of two edge dislocations with opposite Burgers vectors expressed as

$$𝒖^{\mathrm{Ls}\pm }=\pm \left[𝒖^{\mathrm{Le}}(x\frac{\mathrm{}}{2},y)𝒖^{\mathrm{Le}}(x+\frac{\mathrm{}}{2},y)\right],$$
(3.4)
where $`\mathrm{}`$ is the slip length and $`b=\pm 1`$ (both in units of $`a`$). The slip line segment is between the two points $`(\mathrm{}/2,0)`$ and $`(\mathrm{}/2,0)`$. The $`+`$ sign corresponds to a clockwise slip (type C), and the $``$ sign to a counterclockwise slip (type CC). The displacement vector around a slip is clockwise (counterclockwise) for type C (type CC) (as will be evident in Fig.5 below). Across the line segment $`|x|<\mathrm{}/2`$ and $`y=0`$, the displacement is discontinuous as

$$u_x^{\mathrm{Ls}\pm }(x,y+0)u_x^{\mathrm{Ls}\pm }(x,y0)=\pm 1.$$
(3.5)
The corresponding strains are written as

$`e_1^{\mathrm{Ls}\pm }`$ $`=`$ $`{\displaystyle \frac{\pm (2\nu 1)}{2\pi (1\nu )}}\left[{\displaystyle \frac{y}{x_+^2+y^2}}{\displaystyle \frac{y}{x_{}^2+y^2}}\right],`$ (3.6)
$`e_2^{\mathrm{Ls}\pm }`$ $`=`$ $`{\displaystyle \frac{1}{\pi (1\nu )}}\left[{\displaystyle \frac{x_+^2y}{(x_+^2+y^2)^2}}{\displaystyle \frac{x_{}^2y}{(x_{}^2+y^2)^2}}\right],`$ (3.7)
$`e_3^{\mathrm{Ls}\pm }`$ $`=`$ $`{\displaystyle \frac{\pm 1}{2\pi (1\nu )}}\left[{\displaystyle \frac{x_+(x_+^2y^2)}{(x_+^2+y^2)^2}}{\displaystyle \frac{x_{}(x_{}^2y^2)}{(x_{}^2+y^2)^2}}\right]`$ (3.9)
$`\pm \delta (y)\mathrm{\Theta }(\mathrm{}^2/4x^2),`$
where $`x_\pm =x\mathrm{}/2`$, and $`\mathrm{\Theta }(\zeta )`$ is the step function being equal to 1 for $`\zeta >0`$ and to $`0`$ for $`\zeta <0`$. There appears no $`\delta `$-function in the dilation and tetragonal strains. The strains diverge at the cores where the linear elasticity theory breaks down.

### B Slips in nonlinear elasticity theory

Next we numerically construct the corresponding slip solution in our nonlinear elasticity theory for integer $`\mathrm{}`$. That is, by starting with the linear solution $`𝒖^{\mathrm{Ls}\pm }`$ in (3.3) at $`t=0`$ and neglecting the random stress $`\stackrel{}{\sigma }_R`$, we integrate (2.16) and (2.17) to seek the steady solution $`𝒖^{\mathrm{s}\pm }`$ attained after transient relaxation. That is,

$$𝒖^{\mathrm{s}\pm }(x,y)=\underset{t\mathrm{}}{lim}𝒖(x,y,t),$$
(3.10)
where $`𝒖(x,y,0)=𝒖^{\mathrm{Ls}\pm }(x,y)`$. The core points of the linear solution, where the strains diverge, are placed at middle points $`(n+1/2,m+1/2)`$ of the mesh of integration at $`t=0`$. The limit $`𝒖^{\mathrm{s}\pm }`$ is a steady metastable solution satisfying the mechanical equilibrium condition,

$$\stackrel{}{\sigma }=\mathrm{𝟎}.$$
(3.11)
It nearly coincides with the linear solution $`𝒖^{\mathrm{Ls}\pm }`$ far from the dislocation cores ($`|x\pm \mathrm{}/2|^2+y^23`$ in our case) and keeps to satisfy the slip condition (3.4). The stain and stress components calculated from $`𝒖^{\mathrm{s}\pm }`$ are finite even in the core regions. In fact, in the presence of a slip along the $`x`$ direction, the maximum values attained by $`|e_1|`$, $`|e_2|`$, $`|e_3|`$, and $`|\sigma _{xy}|`$ in the core regions are about $`0.18`$, $`0.08`$, $`1.1`$, and $`0.1`$, respectively. In Fig.4 we show the $`x`$ component $`u_x=u_x^{\mathrm{s}+}`$ of a clockwise slip with $`\mathrm{}=20`$ in the unstrained condition, which is discontinuous by 1 across the slip segment.

The elastic energy of a slip in our nonlinear theory is then of central importance. We will neglect the Peierls potential energy for the time being. We generally assume that a slip line is oblique to the $`x`$ axis making an angle of $`\phi `$ and is under externally applied strains,

$$e_3=\gamma ,e_2=ϵ.$$
(3.12)
We assume $`e_1=0`$. The average displacement is given by $`u_x=\gamma y+ϵx/2`$ and $`u_y=ϵy/2`$. Here $`\mathrm{}`$ is the space average. The slip energy to create a single slip is defined by

$$F_{\mathrm{slip}}=FF_0=𝑑𝒓(f_{\mathrm{el}}f_{\mathrm{el}}^0),$$
(3.13)
where $`f_{\mathrm{el}}`$ is the elastic energy density in our nonlinear theory with one slip calculated numerically, and $`f_{\mathrm{el}}^0=\mathrm{\Phi }(\gamma ,ϵ)`$ is that in a homogeneously strained state. In the unstrained condition ($`\gamma =ϵ=0)`$, the expression $`F_{\mathrm{slip}}=\mathrm{ln}\mathrm{}/2\pi (1\nu )`$ (in units of $`\mu _0a^2`$) is well known .

Dislocation motions perpendicular to the slip line (climb motions) create a large number of defects ($`\mathrm{}`$) and the energy needed is very large, so they may be neglected . On the other hand, dislocation motions along the slip line (glide motions) involve displacements of a relatively small number of particles (of order 10 for a unit-length motion as can be seen in the inset of Fig.9) and hence play a major role in plastic flow. A force $`𝓕=(_x,_y)`$ acting on a dislocation under an applied stress $`\{\sigma _{ij}^{\mathrm{ex}}\}`$ may be calculated using the Peach-Koehler theory . For $`𝒃=(b,0)`$ along the $`x`$ axis, the components of the force are written as

$$_x=\sigma _{xy}^{\mathrm{ex}}b,_y=\sigma _{xx}^{\mathrm{ex}}b.$$
(3.14)
(i) For a slip along the $`x`$ axis, $`_y`$ is canceled by the force due to an additional elastic deformation in the surrounding medium. If the dislocation on the left is fixed and that on the right is moved by $`\delta \mathrm{}`$ along the $`x`$ axis, the change of $`F_{\mathrm{slip}}`$ is equal to $`_x\delta \mathrm{}`$. Here the Peierls force is neglected. Therefore,

$$\frac{}{\mathrm{}}F_{\mathrm{el}}=_x=\sigma _{xy}^{\mathrm{ex}}b.$$
(3.15)
The $`\sigma _{xy}^{\mathrm{ex}}`$ consists of the stress produced by the dislocation on the left and the externally applied stress. For small $`\gamma `$ we obtain

$$F_{\mathrm{slip}}=\frac{\mathrm{ln}\mathrm{}}{2\pi (1\nu )}\gamma \mathrm{},$$
(3.16)
where $``$ is for type C and $`+`$ is for type CC. (ii) For general angle $`\phi `$ of the slip with respect to the $`x`$ axis, we rotate the reference frame by $`\phi `$ to obtain the shear strain $`\gamma ^{}=\gamma \mathrm{cos}2\phi ϵ\mathrm{sin}2\phi `$ in the new reference frame from (2.2). Therefore,

$$F_{\mathrm{slip}}=\frac{\mathrm{ln}\mathrm{}}{2\pi (1\nu )}(\gamma \mathrm{cos}2\phi ϵ\mathrm{sin}2\phi )\mathrm{}.$$
(3.17)
If $`\phi `$ is varied in (3.13), $`F_{\mathrm{slip}}`$ is minimized for

$$\phi =\frac{1}{2}(n\pi \alpha )(n=0,1,2,\mathrm{}),$$
(3.18)
where $`n`$ is even for type C and odd for type CC, and $`\alpha `$ is determined by

$$\mathrm{cos}\alpha =\frac{\gamma }{\sqrt{\gamma ^2+ϵ^2}},\mathrm{sin}\alpha =\frac{ϵ}{\sqrt{\gamma ^2+ϵ^2}}.$$
(3.19)
For simple shear deformation with $`\gamma >0`$ and $`ϵ=0`$, the most favorable slip orientation with the lowest $`F_{\mathrm{slip}}`$ is $`\phi =0`$ for type C and $`\phi =\pi /2`$ for type CC. For uniaxial stretching with $`\gamma =0`$ and $`ϵ>0`$, it is given by $`\phi =\pi /4`$ for type C and $`\phi =\pi /4`$ for type CC, in agreement with the experiments .

As will be shown in the appendix, the slip directions determined by (3.14) are perpendicular to the wave vectors which minimize the angle-dependent sound velocity following from our linearized dynamic equations, or perpendicular to the softest directions for the sound modes. This coincidence is obtained with the aid of the isotropic behavior (2.8), so it is not a general result for anisotropic solids. We remark that previous theories of strain localization  are based on 1D analysis, where all the quantities vary only in one direction $`𝒏`$ normal to the plane of the band, and reduce to linear stability analysis for small amplitude perturbations for the wave vector in the direction of $`𝒏`$.

In Fig.5 we display the displacement vector $`𝒖`$ around type C and type CC slips calculated in the unstrained condition. The slips are oriented in the most favorable directions in shear deformation in (a) and in uniaxial stretching in (b). Away from the slips the directions of $`𝒖`$ continuously change to those of the macroscopic deformation supposed to be applied. In Fig.6 the slip energy $`F_{\mathrm{slip}}`$ of type C slips with $`\mathrm{}=10,20,`$ and 30 is shown as a function of the applied shear strain $`\gamma `$. For $`|\gamma |0.05`$ we confirm (3.12) with the $``$ sign. For larger $`\gamma `$ the linear relation $`F_{\mathrm{slip}}/\gamma `$ =const. does not hold. For $`\gamma <\gamma _{\mathrm{c1}}(0.1)`$ or for $`\gamma >\gamma _{\mathrm{c2}}(0.1)`$, a steady metastable solution becomes nonexistent and, as a result, the slip grows up to the system length or shrinks to vanish in the simulation. In Fig.6 the slip energy $`f_{\mathrm{el}}f_{\mathrm{el}}^0`$ with $`f_{\mathrm{el}}^0=\mathrm{\Phi }(\gamma ,0)`$ is displayed for the three values $`\gamma =0,0.065`$, and $`0.04`$. Interestingly, the elastic energy density in the middle region between the two dislocations at the ends is decreased for positive $`\gamma `$ and increased for negative $`\gamma `$, giving rise to the contribution $`\mathrm{}\gamma `$ in the slip energy. In the core region, $`f_{\mathrm{el}}`$ has two peaks and is rather insensitive to applied strains in our model.

Furthermore, in Fig.8 we numerically demonstrate coincidence of $`F_{\mathrm{slip}}/\gamma `$ and the space integral of $`\sigma _{xy}\gamma `$ for a single isolated slip in simple shear deformation. This relation may be obtained from (2.1). If $`\gamma `$ is increased by an infinitesimal amount $`\delta \gamma `$, the change of $`F_{\mathrm{slip}}`$ is written in terms of the incremental displacement $`\delta u_i`$ as

$$\delta F_{\mathrm{slip}}=𝑑𝒓\underset{ij}{}[\sigma _{ij}\sigma _{ij}^0]\frac{}{x_j}\delta u_i,$$
(3.20)
where $`\{\sigma _{ij}^0\}`$ is the stress in the homogeneous state and use is made of the relation for the space average $`\delta u_i/x_j=\delta _{ix}\delta _{jy}\delta \gamma `$. The deviation $`\delta u_i`$ consists of the applied displacement change $`\delta _{ix}y\delta \gamma `$ and the induced deviation $`\delta u_i^{}`$ localized near the slip. However, the contribution from $`\delta u_i^{}`$ vanishes in (3.16) from the mechanical equilibrium condition (3.7). Thus,

$$\frac{}{\gamma }F_{\mathrm{slip}}(\mathrm{},\gamma )=𝑑𝒓[\sigma _{xy}\sigma _{xy}^0].$$
(3.21)
In Fig.8 $`\gamma `$ is in the range $`|\gamma |0.1`$, so $`\sigma _{xy}^0\gamma `$ holds excellently for our elastic energy. Similarly, the counterpart of (3.17) in the case of uniaxial stretching is written as

$$\frac{}{ϵ}F_{\mathrm{slip}}(\mathrm{},ϵ)=\frac{1}{2}𝑑𝒓[\sigma _{xx}\sigma _{yy}\sigma _{xx}^0+\sigma _{yy}^0].$$
(3.22)
The relations (3.17) and (3.18) hold for any strain amplitudes as long as a steady slip solution can be obtained from (3.6), while (3.12) or (3.13) is valid only for very small strains ($`0.05`$ in our case).

### C Peierls potential energy

We continue to consider a slip along the $`x`$ axis under simple shear deformation, but the slip length $`\mathrm{}`$ here can be noninteger. For general $`\mathrm{}`$ we modify (3.12) as

$$F_{\mathrm{slip}}=\frac{\mathrm{ln}\mathrm{}}{2\pi (1\nu )}\gamma \mathrm{}+U_{\mathrm{PN}}(\mathrm{}),$$
(3.23)
where $`U_{\mathrm{PN}}(\mathrm{})`$ represents the Peierls potential energy being zero for integer $`\mathrm{}`$ and positive for noninteger $`\mathrm{}`$. The force acting on the slip along the glide direction is then given by

$$\frac{}{\mathrm{}}F_{\mathrm{slip}}=[2\pi (1\nu )]^1\frac{1}{\mathrm{}}\gamma +\frac{}{\mathrm{}}U_{\mathrm{PN}}(\mathrm{}).$$
(3.24)
This force vanishes for a stationary slip. In the realm of linear elasticity theory, where the crystal structure is smoothed out, the Peierls potential energy is nonexistent and we are led to the following conclusion: A type CC slip shrinks for any $`\mathrm{}(>1)`$, but a type C slip shrinks for $`\mathrm{}<\mathrm{}_c^\mathrm{L}`$ and expands for $`\mathrm{}>\mathrm{}_c^\mathrm{L}`$. Here we assume $`\gamma >0`$ and the critical length $`\mathrm{}_c^L`$ in the linear elasticity theory is determined by

$$\mathrm{}_\mathrm{c}^\mathrm{L}=1/[2\pi (1\nu )\gamma ].$$
(3.25)
If $`\mathrm{}`$ is fixed, we obtain a critical shear strain $`\gamma _\mathrm{c}^\mathrm{L}=1/[2\pi (1\nu )\mathrm{}]`$. In our nonlinear theory, however, this critical length (or strain) loses its physical relevance.

In our simulations (without the random stress) slips can be stationary suggesting the existence of the Peierls potential energy. To show this, we numerically create two type C slips along the $`x`$ axis obtained in the limit (3.6); one is in the range $`0x20`$ with length 20, and the other is in the range $`1x20`$ with length 21. Here the positions of the dislocation core on the left are different by $`1`$ but those on the right coincide. As a result, the corresponding displacements $`𝒖_{20}`$ and $`𝒖_{21}`$ are different only near the core region on the left, as can be seen in the inset of Fig.9. Note that the difference $`𝒖^{}=𝒖_{21}𝒖_{20}`$ is the displacement realized when the shorter slip grows into the longer one. Now, we calculate $`F`$ for the interpolated displacement,

$$𝒖_{\mathrm{}}=(1\alpha )𝒖_{20}+\alpha 𝒖_{21}=𝒖_{20}+\alpha 𝒖^{}.$$
(3.26)
As the slip energy at $`\mathrm{}=20+\alpha `$ ($`0<\alpha <1`$), $`F_{\mathrm{slip}}(\mathrm{})=FF_0`$ is determined as in (3.9). The $`F_{\mathrm{slip}}(\mathrm{})`$ here depends on the displacement path connecting $`𝐮_{20}`$ and $`𝐮_{21}`$.

In Fig.9 the resultant energy difference $`\mathrm{\Delta }F_{\mathrm{slip}}=F_{\mathrm{slip}}(\mathrm{})F_{\mathrm{slip}}(20)`$ is shown. From this calculation the Peierls potential cannot be determined uniquely but, if we define

$$U_{\mathrm{PN}}(\mathrm{})=F_{\mathrm{slip}}(\mathrm{})(1\alpha )F_{\mathrm{slip}}(20)\alpha F_{\mathrm{slip}}(21),$$
(3.27)
we can know its behavior in the range $`0<\alpha =\mathrm{}20<1`$ directly from Fig.9. We recognize that $`F_{\mathrm{slip}}(\mathrm{})`$ takes a maximum at $`\mathrm{}=\mathrm{}_{\mathrm{max}}`$ between $`20`$ and $`21`$ and takes local minima at $`\mathrm{}=20`$ and 21. It follows that the stable force-balance condition $`F_{\mathrm{slip}}/\mathrm{}=0`$ with $`^2F_{\mathrm{slip}}/\mathrm{}^2>0`$ holds at $`\mathrm{}=20`$ and 21. The maximum of $`U_{\mathrm{PN}}`$ is about 0.03 for $`\gamma =0`$. Fig.9 also indicates that $`\mathrm{}_{\mathrm{max}}21`$ as $`\gamma \gamma _{\mathrm{c1}}(0.09)`$ and $`\mathrm{}_{\mathrm{max}}20`$ as $`\gamma \gamma _{\mathrm{c2}}(0.12)`$. For integer $`\mathrm{}`$ and noninteger $`\mathrm{}^{}`$, the critical strains, $`\gamma _{\mathrm{c1}}(\mathrm{})`$ and $`\gamma _{\mathrm{c2}}(\mathrm{})`$, with respect to shrinkage and expansion satisfy the following: (i) $`F_{\mathrm{slip}}(\mathrm{}^{})/\mathrm{}^{}>0`$ for any $`\mathrm{}^{}`$ smaller than $`\mathrm{}`$ for $`\gamma <\gamma _{\mathrm{c1}}(\mathrm{})`$, (ii) $`F_{\mathrm{slip}}(\mathrm{}^{})/\mathrm{}^{}<0`$ for any $`\mathrm{}^{}`$ larger than $`\mathrm{}`$ for $`\gamma >\gamma _{\mathrm{c2}}(\mathrm{})`$, and (iii) $`F_{\mathrm{slip}}(\mathrm{})`$ is locally minimum at integer length $`\mathrm{}`$ for $`\gamma _{\mathrm{c1}}(\mathrm{})<\gamma <\gamma _{\mathrm{c2}}(\mathrm{})`$. These are consistent with the positions of the instability points in Fig.6.

## IV Plastic flow

In this section we induce deformation at a constant strain rate, $`\dot{\gamma }`$ or $`\dot{ϵ}`$, or cyclic shear deformation for $`t0`$ in the presence of the random stress tensor $`\stackrel{}{\sigma }_\mathrm{R}`$ with $`ϵ_{\mathrm{th}}=0.1`$ (except for the curve (b) in Fig.22 where $`ϵ_{\mathrm{th}}=0.25`$). At $`t=0`$, the values of $`𝒗`$ at the lattice points are Gaussian random numbers with variance $`0.01`$. The shear stress and the normal stress difference in the following figures are the space averages, $`\sigma _{xy}`$ and $`\sigma _{xx}\sigma _{yy}`$, respectively.

### A Shear deformation

In Fig.10 we show the stress-strain curves obtained by integration of (2.16) and (2.17) under a constant shear rate $`\dot{\gamma }`$. At small $`\gamma (t)=\dot{\gamma }t`$ these curves first follow the curve in the homogeneous case,

$$\sigma _{xy}A_3(\gamma )=\frac{1}{\sqrt{3}\pi }\mathrm{sin}(\sqrt{3}\pi \gamma ),$$
(4.1)
where $`A_3`$ was introduced in (2.15). For the curve of $`\dot{\gamma }=10^4`$ and that of $`\dot{\gamma }=10^3`$ (with the higher peak) we set $`𝒖=\mathrm{𝟎}`$ at $`t=0`$ (supposing a perfect crystal). They approach the elastic instability point $`\gamma =\sqrt{3}/60.289`$, where $`\sigma _{xy}/\gamma `$ tends to vanish and softening with respect to further shear deformation occurs. Then the shear stress drops sharply after the peak with catastrophic formation of slips. See the two snapshots of $`\delta e_3=e_3\gamma `$ in the figure. The orientations of the slips are the most favorable ones with $`\phi =0`$ and $`\pi /2`$ as determined in (3.14). For $`\dot{\gamma }=10^4`$ the slip formation is triggered at a smaller strain ($`0.25`$) and the spacing between the slips is a few times wider than in the case of $`\dot{\gamma }=10^3`$. In our simulation the slip spacing depends on the shear rate $`\dot{\gamma }`$ in the plastic flow regime with the other dimensionless parameters held fixed. For the other curve of $`\dot{\gamma }=10^3`$ (with the lower peak) we put four slips with length 20 at $`t=0`$, as will be illustrated in Fig.13 below in detail. This curve indicates that the overshoot in the stress-strain relation is weakened by the initially preexisting defects. In Fig.11 we display mechanically unstable points (dots) for $`\gamma =0.29`$, $`0.30`$, and $`0.40`$ at $`\dot{\gamma }=10^3`$, where at least one of the eigenvalues $`\lambda _1`$ and $`\lambda _2`$ of the matrix $`\{\mathrm{\Phi }_{\alpha \beta }\}`$ defined below (2.7) is negative. We recognize that the stability condition $`(\lambda _1>0`$ and $`\lambda _2>0`$) are satisfied at most points in plastic flow. Fig.12 displays a snapshot of the displacement vector $`𝒖`$ in a $`1/4`$ region of the total system at $`\gamma =0.4`$ with $`\dot{\gamma }=10^4`$. We can see a number of slips and bands (aggregates of slips here), where $`\phi =0`$ for type C and $`\phi =\pi /2`$ for type CC. The large horizontal shear band in the lower part is particularly conspicuous, where the band thickness and the discontinuity of $`u_x`$ across it are both increased up to about 3-4. We recognize that elementary slips (dipoles of edge dislocations) tend to be created around preexisting ones, yielding thicker shear bands. Similar thick shear bands have been observed in previous simulations .

In Fig.13 we follow time-development of $`\delta e_3`$ and $`\delta f_{\mathrm{el}}=f_{\mathrm{el}}f_{\mathrm{el}}`$ at $`\dot{\gamma }=10^3`$ in the presence of four slips at $`t=0`$, where the favorable slips (in black in the upper figure at $`\gamma =0`$) grow and the unfavorable ones (in white) shrink as $`\gamma `$ increases. In the lower snapshots the elastic energy density deviation $`\delta f_{\mathrm{el}}=f_{\mathrm{el}}f_{\mathrm{el}}`$ is shown, where the black dots represent the dislocation cores. At $`\gamma =0.176`$ we can see that the energy density between the two dislocations at the slip ends is decreased for the favorable slips and is increased for the unfavorable slips, in accord with (3.13) and Fig.7. At the plastic flow regime $`\gamma =0.387`$ the initial favorable slips grow into thick layers where the dislocation density is very high.

Next we apply a cyclic shear deformation, where $`\dot{\gamma }(t)=10^3`$ in the time regions $`nt_\mathrm{p}<t<(n+1/2)t_\mathrm{p}`$ and $`\dot{\gamma }(t)=10^3`$ in the time regions $`(n+1/2)t_\mathrm{p}<t<(n+1)t_\mathrm{p}`$. We choose $`t_\mathrm{p}=1000`$, so the maximum of $`\gamma (t)`$ is 0.5 and the minimum is 0. For the first two cycles $`(n=0`$ and 1), Fig.15 shows the stress-strain curve, while Fig.16 shows the average elastic energy density $`f_{\mathrm{el}}`$ with a snapshot of $`f_{\mathrm{el}}`$ at the point A where the average stress vanishes. As salient features, we notice (i) residual strain at vanishing stress, (ii) the shear stress becomes negative at the end of the first cycle, (iii) no overshoot in the stress and the elastic energy from the second cycle, and (iv) that the stress and the elastic energy take roughly constant values characteristic of well-developed plastic flow in the region $`0.25\gamma (t)<0.5`$ for increasing $`\gamma (t)`$ and in the region $`0<\gamma (t)0.25`$ for decreasing $`\gamma (t)`$. Thus we can see significant hysteresis behavior. In MD simulations of low $`T`$ glasses, similar stress-strain curves under step-wise strain rates have been obtained (but without overshoot behavior because of disordered initial states) .

We also notice that at the points A,B, and C, where the average stress vanishes as in Fig.14, the curves of $`f_{\mathrm{el}}`$ in Fig.15 are locally minimum. This suggests that the strain $`\gamma (t)`$ consists of an elastic strain $`\gamma _{\mathrm{el}}`$ and an slip strain

$$\gamma _\mathrm{s}=\gamma \gamma _{\mathrm{el}}.$$
(4.2)
Roughly speaking, the elastic strain outside the dislocation cores should give to the average stress, while the slip strain is caused by the jumps of $`u_x`$ across the slips. To be more quantitative, we define $`\gamma _{\mathrm{el}}`$ as

$$A_3(\gamma _{\mathrm{el}})=\sigma _{xy}.$$
(4.3)
The elastic energy density stored is then the sum of the elastic energy density and the defect energy density. We define the average defect energy density by

$$f_\mathrm{D}=f_{\mathrm{el}}\mathrm{\Phi }(\gamma _{\mathrm{el}},0),$$
(4.4)
where $`\mathrm{\Phi }`$ is defined by (2.5). For small $`\gamma _{\mathrm{el}}`$ we have $`\gamma _{\mathrm{el}}\sigma _{xy}`$ and $`\mathrm{\Phi }(\gamma _{\mathrm{el}},0)\gamma _{\mathrm{el}}^2/2`$. Hence, in the vicinity of the points where $`\sigma _{xy}=0`$, $`f_{\mathrm{el}}(t)`$ should take a minimum, provided that $`f_\mathrm{D}(t)`$ changes slowly there. To confirm these arguments, we plot $`f_\mathrm{D}(t)`$ in Fig.16. For $`\gamma (t)0.3`$ in the first cycle, however, it represents the elastic energy due to the inhomogeneous fluctuations of the local strains (mainly due to $`\delta e_3`$) with the peak height at 0.012 (not shown in the figure). After this initial period, $`f_\mathrm{D}(t)`$ is in a range of $`0.0040.005`$ reasonably representing the elastic energy due to the defects produced in plastic flow. The energy variance $`\sqrt{(\delta f_{\mathrm{el}})^2}`$ is also about $`0.005`$ in plastic flow obviously due to the discrete nature of dislocation cores. For $`\dot{\gamma }=10^4`$, on the other hand, $`f_\mathrm{D}(t)`$ is almost constant around $`0.002`$ in plastic flow.

### B Uniaxial stretching

In Fig.17 we show the normal stress vs strain at constant strain rate $`\dot{ϵ}`$ for $`t0`$. The characteristic features are very similar to those in Fig.10. At small $`ϵ(t)=\dot{ϵ}t`$ these curves first follow the curve in the homogeneous case,

$$\sigma _{xx}\sigma _{yy}2A_2(ϵ)=\frac{2}{3\pi }[\mathrm{sin}(\pi ϵ)+\mathrm{sin}(2\pi ϵ)],$$
(4.5)
where $`A_2`$ was introduced in (2.15). For the curve of $`\dot{ϵ}=10^4`$ and that of $`\dot{ϵ}=10^3`$ (with the higher peak) we set $`𝒖=\mathrm{𝟎}`$ at $`t=0`$. They approach the elastic instability point $`ϵ=\pi ^1\mathrm{cos}^1[(\sqrt{33}1)/8]0.298`$. After the catastrophic formation of slips, the slip orientation angles are $`\pm \pi /4`$ with respect to the stretched direction in agreement with the experiments  and the simulations . The discontinuity across the slip lines appears in the tetragonal strain $`e_2`$, which can be seen in the two snapshots of $`\delta e_2=e_2ϵ`$. In Fig.18 the displacement vector $`𝒖`$ in a $`1/4`$ region of the total system at $`ϵ=0.46`$ with $`\dot{ϵ}=10^3`$ is displayed, where the orientations of the slips are the most favorable ones with $`\phi =\pm \pi /4`$ as determined below (3.14). In Fig.19 four slips are placed at $`t=0`$, exactly as in Fig.13, and time-development of $`e_3`$, $`\delta e_2`$, and $`\delta f_{\mathrm{el}}`$ are followed. Here, however, these initial slip orientations are marginal in the sense that $`F_{\mathrm{slip}}`$ in (3.13) is independent of $`ϵ`$. We find that a dislocation pair is newly created near each end of the preexisting slips to glide in the favorable directions with $`\phi =\pm \pi /4`$. In the plastic flow regime at $`ϵ=0.4`$, they grow into thick layers containing high-density dislocations.

Next we apply cyclic uniaxial stretching, where $`\dot{ϵ}(t)=10^3`$ in the time regions $`nt_\mathrm{p}<t<(n+1/2)t_\mathrm{p}`$ and $`\dot{ϵ}(t)=10^3`$ in the time regions $`(n+1/2)t_\mathrm{p}<t<(n+1)t_\mathrm{p}`$. We choose $`t_\mathrm{p}=1000`$, so the maximum of $`ϵ(t)`$ is 0.5 and the minimum is 0. Fig.20 shows the stress-strain curve for the first four cycles, where it is nearly periodic from the second cycle. As in the shear deformation case we introduce the average elastic strain $`ϵ_{\mathrm{el}}`$ by

$$2A_2(ϵ_{\mathrm{el}})=\sigma _{xx}\sigma _{yy}$$
(4.6)
and the average defect energy density $`f_\mathrm{D}`$ by

$$f_\mathrm{D}=f_{\mathrm{el}}\mathrm{\Phi }(0,ϵ_{\mathrm{el}}).$$
(4.7)
where $`\mathrm{\Phi }`$ is defined by (2.5). The average slip strain is given by $`ϵ_\mathrm{s}=ϵϵ_{\mathrm{el}}`$. For small $`ϵ_{\mathrm{el}}`$ we have $`2ϵ_{\mathrm{el}}\sigma _{xx}\sigma _{yy}`$ and $`\mathrm{\Phi }(0,ϵ_{\mathrm{el}})ϵ_{\mathrm{el}}^2/2`$. Fig.21 shows the average elastic energy density $`f_{\mathrm{el}}`$ and the defect energy density $`f_\mathrm{D}`$. The latter is in a range of $`0.0080.009`$ in plastic flow. The inset of Fig.21 displays a snapshot of $`f_{\mathrm{el}}`$ at the point A in the first cycle, where the average normal stress difference vanishes. Comparing it with the snapshot of $`f_{\mathrm{el}}`$ in Fig.15 under shear deformation, we notice a considerable difference in the spatial anisotropy of the dislocation distribution between the two cases.

### C Strain-induced disordered states

In the plastic flow regime dislocations are proliferated and a structurally disordered state is realized. This effect may be called strain-induced disordering. We mention a simulation by Ikeda et al. , who applied a tensile strain to a 3D model to induce a change from a perfect crystal to an amorphous solid. In Fig.4 in our previous work  we switched off a shear flow (a) before the peak time of the stress, (b) just after the peak time, and (c) in well-developed plastic flow. Affine deformation in (a) was maintained, while no appreciable time evolution was detected after transients in (b) and (c). This means that the structurally disordered states are metastable.

In Fig.22, we show extremely slow relaxation of the average energy density $`f_{\mathrm{el}}(t)`$, where a shear flow with $`\dot{\gamma }=10^3`$ is stopped at time $`t_\mathrm{A}600`$ (at the point A in Fig.14) and the system is relaxed for a time period of $`t_\mathrm{w}=10^5`$. Here we set $`u_x=\gamma _\mathrm{A}L_0`$ and $`u_y=0`$ at the top $`y=L_0`$ with $`\gamma _\mathrm{A}0.40`$, together with $`𝒖=\mathrm{𝟎}`$ at the bottom. The dislocation distributions here closely resemble that in the inset of Fig.15 (are not identical because they are obtained from different runs). In the upper curve (a) the noise strength $`ϵ_{\mathrm{th}}`$ in (2.25) is $`0.1`$ as in the previous simulations in this section, where each maximum in the initial stage corresponds to an energy increase accompanied with a configuration change around a dislocation core. From (a) the typical energy increase is of order $`10^5N^20.1`$ for each event. However, there is no appreciable relaxation for $`t3\times 10^3`$. In (b) we increase the noise strength $`ϵ_{\mathrm{th}}`$ to $`0.25`$ to obtain a larger energy decrease in $`f_{\mathrm{el}}(t)`$ with a larger thermal noise superimposed. In these cases only a small number of configuration changes $`(10`$) occur around dislocation cores even for $`t_\mathrm{w}=10^5`$, so the effect of the structural relaxation is negligible on the macroscopic level (no aging effect). In fact, the stress-strain curves after switching-on of the shear flow at $`t=t_\mathrm{A}+t_\mathrm{w}`$ are almost independent of $`t_\mathrm{w}`$ if plotted as a function of $`\mathrm{\Delta }\gamma (t)=\dot{\gamma }(tt_\mathrm{A}t_\mathrm{w})`$. Interestingly, they exhibit a rounded peak at $`\mathrm{\Delta }\gamma =0.15`$ as can be seen in the inset. In addition, in the waiting time region ($`t_\mathrm{A}<t<t_\mathrm{A}+t_\mathrm{w})`$, the average $`\sigma _{xy}(t)`$ fluctuates in time around 0, but its noise amplitude is only of order $`0.002`$ for (a) and $`0.004`$ for (b) and the stress is nearly fixed.

## V Summary and concluding remarks

In summary, we have presented a nonlinear elasticity theory taking into account of periodicity of the elastic energy density with respect to the shear and tetragonal strains, $`e_3`$ and $`e_2`$. It has the symmetry of the 2D triangular lattice but is surprisingly isotropic in the $`e_3`$-$`e_2`$ plane. We summarize our main results together with some comments.

(i) We have numerically examined the slip structure as a function of its length $`\mathrm{}`$, its angle $`\phi `$, and applied strains $`\gamma `$ and $`ϵ`$. In external strain field, slips should appear in the orientations minimizing the slip energy $`F_{\mathrm{slip}}`$ and they should grow into shear bands observed in previous experiments under external load in various materials. The snapshots of the displacement vector in the plastic flow regime, Fig.12 for shear deformation and Fig.18 for uniaxial stretching, most unambiguously illustrate the physical processes taking place. We remark that the previous experiments have mostly been performed under uniaxial (biaxial or triaxial) deformation, so future experiments with shear deformation ($`u_x=\gamma y`$ and $`u_y=0`$) should be informative. On the other hand, in crystalline solids with strong crystal anisotropy, the slip planes are parallel to particular crystal planes . In our simulations we may obtain slips as steady solutions of our dynamic model due to the presence of the Peierls potential energy. However, they cannot be stationary for $`\gamma <\gamma _{\mathrm{c1}}0.1`$ or $`\gamma >\gamma _{\mathrm{c2}}0.1`$ in shear strain $`\gamma `$, where the potential minima disappear and the slips shrink or grow.

(ii) We have examined plastic flow by applying a constant strain rate at $`t=0`$. If there is no initial disorder, the stress-strain curve exhibits a pronounced overshoot with a peak stress of order $`10\%`$ of the shear modulus. However, in the presence of initial disorder, the overshoot is weakened or even erased, as revealed by simulations with initial four slips in Fig.10, under cyclic shear in Figs.14 and 20, and after staying at the zero-stress point A in the inset of Fig.22. In accordance with these findings, previous simulations performed with various initial states have demonstrated sensitive dependence of the overshoot behavior on the quenching conditions . In addition, a number of previous simulations have reported either of existence or nonexistence of the stress peak. In some real glassy systems including polymers, overshoot behavior has been widely observed . In an amorphous metal, the stress increased monotonically to a steady-state value for slow strain rate ($`\dot{ϵ}<10^3`$s<sup>-1</sup>), while a maximum appeared at $`ϵ0.06`$ for large strain rate ($`\dot{ϵ}>5\times 10^3`$s<sup>-1</sup>) .

(iii) As illustrated in Figs. 16 and 21 we have divided the total elastic energy into the affine part and the defect part following the definitions (4.4) and (4.7). In our model, once defects are created, the defect elastic energy is rather weakly dependent on the deformation history. On the basis of this division, we may easily understand the characteristic features of the stress-stain curves in the cyclic straining mentioned in the previous section. In future microscopic simulations on glassy materials, this kind of energy division should be informative at relatively small strains, where we are interested in the defect conribution to measurable quantities such as the specific heats. Here we mention microscopic calculations of the average potential energy per particle $`e`$ in supercooled states under shear , where $`e`$ was increased considerably by shear flow above the initial value in quiescent states. It is remarkable that, while $`e`$ in quiescent states sensitively depends on the quenching history (aging effect) , it becomes uniquely determined in shear flow for shear rates larger than the inverse structural relaxation time $`\tau _\alpha ^1`$. This is because sheared systems are effectively driven away from the glass transition . This effect would be consistent with our result that $`f_\mathrm{D}(t)`$ is kept nearly constant in plastic flow as in Figs.16 and 21.

(iv) In our simulations slips emerge as long straight lines, as shown in the snapshots of $`e_3`$ or $`e_2`$. In our model the crystal order is not broken over long distances, but if disorder is fully introduced, the glide motions of slips in particular directions should be much limited. In MD simulations of two component glasses, for example, such degree of disorder should be sensitive to the size-ratio of the two species. In fact, it was rather close to 1 in the simulation by Deng et al. , where nano-crystalline order was realized  because of rather weak frustration and long slips with atomic thickness emerged along the crystal axes. It is therefore informative to perform MD simulations with various size-ratios and examine the shapes of local configuration changes.

(v) As a special ingredient of our theory, our elastic energy is almost isotropic if the distance from the center of a unit cell in the $`e_3`$-$`e_2`$ plane is shorter than $`0.5`$ (see Figs.1-3). The snapshot at $`\gamma =0.4`$ in Fig.11 demonstrates that most spatial points are in the mechanically stable regions in plastic flow and hence are in the isotropic elasticity regions. Therefore, our results such as the orientations of well-developed shear bands should be applicable to those in amorphous materials (which are isotropic on large scales).

(vi) The dimensionless parameter $`\lambda `$ in (2.23) is related to Poisson’s ratio $`\nu `$ as in (3.2) and the dislocation energy depends on $`\nu `$. Our choice $`K_0/\mu _0=\lambda =1`$ or $`\nu =0`$ is rather unusual, since $`K_0`$ is usually considerably larger than $`\mu _0`$ in high density systems. The appropriateness of the other choices $`\eta _0^{}=1`$ in (2.26) and $`ϵ_{\mathrm{th}}=0.1`$ in Section 4 should also be examined in future.

(vii) As another aspect, we mention the effect of elastic interaction among dislocations. In our case slips are more easily created around preexisting ones, as already reported in Ref.. In the insets of Figs.15 and 21 we can see a tendency of aggregation of dislocation cores. In real 3D crystals a tangle of dislocations often appears as a characteristic feature of fatigued states  and has also been realized in 3D large-scale simulations . In the literature the inter-dislocation elastic interaction is believed to yield mesoscopic patterns in the dislocation distribution .

(viii) As already stressed in the introduction, in order to describe glass dynamics in a more satisfactory level, we should try to include a variable representing the local free-volume  and the configurational frustration effect induced by the size disparity between the two species. This generalization should be essential with increasing $`T`$ towards the glass transition. For example, we may predict a gradual diffusional increase of the free volume around dislocation cores , which induces configuration changes and presumably resulting in significant aging effects. (The critical strains $`\gamma _{\mathrm{c1}}`$ and $`\gamma _{\mathrm{c2}}`$ discussed below (3.23) and indicated in Fig.6 by arrows are sensitively decreased in magnitude by a small amount of the local free volume near the dislocation cores.) In the present study, as shown in Fig.22, we have found no appreciable aging effect.

## ACKNOWLEDGMENTS

This work is supported by Grants in Aid for Scientific Research from the Ministry of Education, Science and Culture. The author thanks R. Yamamoto for valuable discussions and A. Furukawa for his help in writing the figures. Thanks are also due to A.S. Argon for sending his papers and informative correspondence.

APPENDIX

Here we examine the linearized equations of (2.16) and (2.17) around a homogeneously strained state with $`e_1=0`$, $`e_3=\gamma `$, and $`e_2=ϵ`$. We neglect the random stress and assume that all the deviations depend on space and time as $`\mathrm{exp}(i𝒌𝒓+i\omega t)`$. Then the deviation of the displacement vector $`(\delta u_x,\delta u_y)`$ obeys

$`\mathrm{\Omega }\delta u_x`$ $`=`$ $`C_{xx}\delta u_x+C_{xy}\delta u_y,`$ (A.1)
$`\mathrm{\Omega }\delta u_y`$ $`=`$ $`C_{xy}\delta u_x+C_{yy}\delta u_y,`$ (A.2)
where

$$\mathrm{\Omega }=\rho \omega ^2/k^2i\eta _0\omega .$$
(A.3)
The frequency $`\omega `$ is expressed in terms of $`\mathrm{\Omega }`$ at small $`k`$ as

$$\omega =\pm (\mathrm{\Omega }/\rho )^{1/2}k+i(\eta _0/2\rho )k^2+O(k^3).$$
(A.4)
The first term is the oscillation frequency for $`\mathrm{\Omega }>0`$ and the second term viscous damping. The coefficients $`C_{\alpha \beta }`$ are expressed in terms of the coefficients $`\mathrm{\Phi }_{\alpha \beta }`$ in (2.7) as

$`C_{xx}=(K_0+\mathrm{\Phi }_{22})n_x^2+2\mathrm{\Phi }_{23}n_xn_y+\mathrm{\Phi }_{33}n_y^2,`$ (A.5)
$`C_{yy}=(K_0+\mathrm{\Phi }_{22})n_y^22\mathrm{\Phi }_{23}n_xn_y+\mathrm{\Phi }_{33}n_x^2,`$ (A.6)
$`C_{xy}=(K_0\mathrm{\Phi }_{22}+\mathrm{\Phi }_{33})n_xn_y+\mathrm{\Phi }_{23}(n_x^2n_y^2).`$ (A.7)
Here $`𝒏=k^1𝒌`$ is the unit vector representing the direction of the wave vector $`𝒌`$. From (A.1) it follows the relation,

$$(\mathrm{\Omega }C_{xx})(\mathrm{\Omega }C_{yy})=C_{xy}^2.$$
(A.8)
Let the angle of $`𝒏`$ be $`\phi +\pi /2`$ with respect to the $`x`$ axis; then,

$$n_x=\mathrm{sin}\phi ,n_y=\mathrm{cos}\phi .$$
(A.9)
After some calculations (A.5) is solved to give

$`\mathrm{\Omega }`$ $`=`$ $`{\displaystyle \frac{1}{2}}(K_0+\mathrm{\Phi }_{22}+\mathrm{\Phi }_{33})`$ (A.11)
$`\pm \sqrt{{\displaystyle \frac{1}{4}}K_0^2+A^2+B^2+K_0(A\mathrm{cos}4\phi +B\mathrm{sin}4\phi )},`$
where $`A=(\mathrm{\Phi }_{22}\mathrm{\Phi }_{33})/2`$ and $`B=\mathrm{\Phi }_{23}`$. As a function of $`\phi `$ the slowest mode is obtained if the combination $`A\mathrm{cos}4\phi +B\mathrm{sin}4\phi `$ takes the maximum $`(A^2+B^2)^{1/2}`$. The corresponding minimum of $`\mathrm{\Omega }`$ is given by

$$\mathrm{\Omega }_{\mathrm{min}}=\frac{1}{2}(\mathrm{\Phi }_{22}+\mathrm{\Phi }_{33})\sqrt{A^2+B^2}.$$
(A.12)
Notice that $`\mathrm{\Omega }_{\mathrm{min}}`$ is the smaller of the two eigenvalues of the matrix $`\{\mathrm{\Phi }_{\alpha \beta }\}`$. We may draw two conclusions. (i) As the instability point is approached, $`\mathrm{\Omega }_{\mathrm{min}}`$ tends to zero with softening of the sound speed $`(\mathrm{\Omega }_{\mathrm{min}}/\rho )^{1/2}`$ of the corresponding acoustic mode. (ii) If the approximate expression (2.8) is used, we have $`A=2G^{\prime \prime }(ϵ^2\gamma ^2)`$ and $`B=4G^{\prime \prime }ϵ\gamma `$ so that

$$A\mathrm{cos}4\phi +B\mathrm{sin}4\phi =C\mathrm{cos}(4\phi +2\alpha ),$$
(A.13)
where $`C=2|G^{\prime \prime }|(ϵ^2+\gamma ^2)`$ and $`\alpha `$ is defined by (3.15). Thus the minimum condition for $`\mathrm{\Omega }`$ yielding (A.8) is given by $`4\phi +2\alpha =2n\pi `$ and is equivalent to (3.14). That is, the slowest mode, which undergoes softening at the instability point, has a wave vector perpendicular to the favorable slip orientations given by (3.14). This is the case even far below the instability point. In addition, for the slowest mode, (A.1) gives

$$\delta u_y/\delta u_x=(\mathrm{\Omega }_{\mathrm{min}}C_{xx})/C_{xy}=\mathrm{cot}\phi ,$$
(A.14)
if $`\alpha `$ is eliminated using (3.14). Thus the deviation $`\delta 𝒖`$ is perpendicular to the wavevector $`𝒌`$ or the slowest mode is a transverse sound.

FIG.9. The slip energy difference $`\mathrm{\Delta }F_{\mathrm{slip}}(\mathrm{})=F_{\mathrm{slip}}(\mathrm{})F_{\mathrm{slip}}(20)`$ in the range $`20\mathrm{}21`$ obtained for the extrapolated displacement (3.22) for $`\gamma =0.04,0,0.08`$, and $`0.1`$. The position of the maximum is a decreasing function of $`\gamma `$. The maximum position approaches $`\mathrm{}=21`$ as $`\gamma \gamma _{\mathrm{c1}}0.09`$ and $`\mathrm{}=20`$ as $`\gamma \gamma _{\mathrm{c2}}0.12`$. The inset shows the displacement vector $`𝒖_{21}𝒖_{20}`$ needed for slip growth by unit length from $`\mathrm{}=20`$ to 21.

FIG.14. The stress-strain curve for cyclic shear deformation in the first two cycles at $`\dot{\gamma }=\pm 10^3`$ with period $`10^3`$. Once $`\gamma 0.3`$ in the first cycle, high-density dislocations are created. The shear stress vanishes at the three points A, B, and C ($`\times )`$.

FIG.20. The stress-strain curve for cyclic uniaxial stretching in the first four cycles at $`\dot{ϵ}=\pm 10^3`$ with period $`10^3`$. Once $`ϵ0.3`$ in the first cycle, high-density dislocations are created. In this case the stress-strain relation becomes nearly periodic from the second cycle.