Theory of Shear Modulus in Glasses

TL;DR

This paper develops a linear response theory for shear deformations in glasses, linking shear modulus to stress and particle correlations, and confirms predictions through numerical simulations, revealing long-range correlations and sound amplification effects.

Contribution

It introduces a general linear response framework for glasses under shear, incorporating boundary effects and inhomogeneous responses, with numerical validation and application to sound propagation.

Findings

Shear modulus expressed via stress-force correlations.

Presence of long-range, long-lived force-displacement correlations.

Resonant amplification of transverse sound near mode frequency.

Abstract

We construct a linear response theory of applying shear deformations from boundary walls in the film geometry in Kubo's theoretical scheme. Our method is applicable to any solids and fluids. For glasses, we assume quasi-equilibrium around a fixed inherent state. Then, we obtain linear-response expressions for any variables including the stress and the particle displacements, even though the glass interior is elastically inhomogeneous. In particular, the shear modulus can be expressed in terms of the correlations between the interior stress and the forces from the walls. It can also be expressed in terms of the inter-particle correlations, as has been shown in the previous literature. Our stress relaxation function includes the effect of the boundary walls and can be used for inhomogeneous flow response. We show the presence of long-ranged, long-lived correlations among the fluctuations of the forces from the walls and the displacements of all the particles in the cell. We confirm these theoretical results numerically in a two-dimensional model glass. As an application, we describe propagation of transverse sounds after boundary wall motions using these time-correlation functions We also find resonant sound amplification when the frequency of an oscillatory shear approaches that of the first transverse sound mode.

Tables

Table 1. Table 1: Vibrational eigenmodes and projection coefficients.

$\lambda$		1	2	3	4	5	6	7	8
${\omega }_{\lambda }$		0.174	0.352	0.417	0.431	0.481	0.517	0.546	0.551
$Z_W^{\lambda }$		1.375	11.76	-0.226	-1.730	-7.920	-9.957	11.71	-13.64
$Z_A^{\lambda }$		37.83	-66.11	24.41	-13.08	231.3	-153.7	118.7	26.22

Equations

$${\mathcal{H}}^{\prime }=\mathcal{H}-{\gamma }_{\mathrm{ex}}(t)\mathcal{A}.$$

$$\psi (r)=\frac{1}{2}{s}_0{r}^2.$$

$$U=\frac{1}{2}\sum _{i,j}{\varphi }_{ij}+\sum _{k>N}{\psi }_k,$$

$${\dot{\mbox{\boldmath$J$}}}=-\nabla\cdot{\stackrel{{\scriptstyle\leftrightarrow}}{{\Pi}}}-\sum_{k>N}\delta({\mbox{\boldmath$r$}}-{\mbox{\boldmath$r$}}_{k})\nabla_{k}\psi_{k}.$$

$$\displaystyle\Pi_{\alpha\beta}^{\rm K}({\mbox{\boldmath$r$}})$$

$$\displaystyle\Pi_{\alpha\beta}^{\rm p}({\mbox{\boldmath$r$}})$$

$${\hat{\delta}}({\mbox{\boldmath$r$}},{\mbox{\boldmath$r$}}_{i},{\mbox{\boldmath$r$}}_{j})=\int_{0}^{1}d\lambda\delta({\mbox{\boldmath$r$}}-\lambda{\mbox{\boldmath$r$}}_{i}-(1-\lambda){\mbox{\boldmath$r$}}_{j}),$$

$$\int d{\mbox{\boldmath$r$}}\Pi_{z\alpha}({\mbox{\boldmath$r$}})-\int d{\mbox{\boldmath$r$}}z{\dot{J}}_{\alpha}({\mbox{\boldmath$r$}})=\sum_{k>N}z_{k}\nabla_{k}^{\alpha}\psi_{k},$$

$${\cal A}=\frac{H}{2}(F_{\rm bot}^{x}-F_{\rm top}^{x})=\frac{H}{2}\Big{(}\sum_{k\in{\rm top}}-\sum_{k\in{\rm bot}}\Big{)}\nabla_{k}^{x}\psi_{k}.$$

$$\mathcal{A}$$

$$\sum_{i}{\mbox{\boldmath$f$}}_{i}={\mbox{\boldmath$F$}}_{\rm top}+{\mbox{\boldmath$F$}}_{\rm bot}.$$

$$W^{\rm ih}=\int d{\mbox{\boldmath$r$}}\Pi_{zx}^{\rm ih}({\mbox{\boldmath$r$}})=\frac{1}{2}H(F_{\rm bot}^{x}-F_{\rm top}^{x})\quad(T\to 0),$$

$${\mbox{\boldmath$u$}}_{i}=(u_{i}^{x},u_{i}^{y},u_{i}^{z})={\mbox{\boldmath$r$}}_{i}-{\mbox{\boldmath$r$}}_{i}^{{\rm ih}},$$

$${\cal H}_{\rm ih}=\sum_{i}\frac{1}{2m_{i}}|{\mbox{\boldmath$p$}}_{i}|^{2}+\delta U,$$

$$\delta U=\frac{1}{4}\sum _{i,j,\alpha ,\beta }{s}_{ij}^{\alpha \beta }{u}_{ij}^{\alpha }{u}_{ij}^{\beta }+\frac{1}{2}\sum _{k>N,\alpha }{s}_0|u_k^{\alpha }{|}^2.$$

$$s_{ij}^{\alpha \beta }={\nabla }_i^{\alpha }{\nabla }_i^{\beta }{\varphi }_{ij}=(\frac{{\varphi }_{ij}^{\prime \prime }}{r_{ij}^2}-\frac{{\varphi }_{ij}^{\prime }}{r_{ij}^3})x_{ij}^{\alpha }{x}_{ij}^{\beta }+\frac{{\varphi }_{ij}^{\prime }}{r_{ij}}{\delta }_{\alpha \beta },$$

$$P_{\mathrm{ih}}\propto \exp (-{\mathcal{H}}_{\mathrm{ih}}/k_{B}T).$$

$$f_i^{\alpha }=-\frac{\partial }{\partial u_i^{\alpha }}\delta U=-\sum _{j,\beta }{h}_{ij}^{\alpha \beta }{u}_j^{\beta }.$$

$$h_{ij}^{\alpha \beta }=\sum _k{s}_{ik}^{\alpha \beta }{\delta }_{ij}-s_{ij}^{\alpha \beta }+s_0{\delta }_{\alpha \beta }{\delta }_{ij}{\theta }_{i-N},$$

$$(h^{-1}{)}_{ij}^{\alpha \beta }=⟨u_i^{\alpha }{u}_j^{\beta }{⟩}_{\mathrm{ih}}/k_{B}T,$$

$$-⟨\mathcal{B}{f}_i^{\alpha }{⟩}_{\mathrm{ih}}=\sum _{k,\beta }{h}_{ik}^{\alpha \beta }⟨u_k^{\beta }\mathcal{B}{⟩}_{\mathrm{ih}}=k_{B}T⟨\partial \mathcal{B}/\partial u_i^{\alpha }{⟩}_{\mathrm{ih}}.$$

$$\Pi^{\rm p}_{\alpha\beta}({\mbox{\boldmath$r$}})=\Pi_{\alpha\beta}^{\rm ih}({\mbox{\boldmath$r$}})+\delta\Pi^{\rm p}_{\alpha\beta}({\mbox{\boldmath$r$}})+\cdots.$$

$$\displaystyle\delta\Pi^{\rm p}_{\alpha\beta}({\mbox{\boldmath$r$}})=-\frac{1}{2}\sum_{i,j,\nu}w_{ij}^{{\alpha\beta}\nu}u_{ij}^{\nu}{\hat{\delta}}({\mbox{\boldmath$r$}},{\mbox{\boldmath$r$}}_{i},{\mbox{\boldmath$r$}}_{j})$$

$$\displaystyle+\sum_{i,j}\frac{{\phi_{ij}^{\prime}}{x_{ij}^{\alpha}}}{r_{ij}}({\mbox{\boldmath$u$}}_{i}\cdot\nabla)[(x_{\beta}-x_{j}^{\beta}){\hat{\delta}}({\mbox{\boldmath$r$}},{\mbox{\boldmath$r$}}_{i},{\mbox{\boldmath$r$}}_{j})].$$

$$w_{ij}^{\alpha \beta \nu }={\nabla }_i^{\nu }[{\varphi }_{ij}^{\prime }{x}_{ij}^{\alpha }{x}_{ij}^{\beta }/r_{ij}]=s_{ij}^{\alpha \nu }{x}_{ij}^{\beta }+{\nabla }_i^{\alpha }{\varphi }_{ij}{\delta }_{\beta \nu }.$$

$$\delta W_{\mathrm{p}}$$

$$w_i^{\alpha }=\sum _j{w}_{ij}^{zx\alpha }=\sum _j[s_{ij}^{z\alpha }{x}_{ij}+\frac{{\varphi }_{ij}^{\prime }}{r_{ij}}{z}_{ij}{\delta }_{x\alpha }].$$

$$\delta \mathcal{A}=\frac{1}{2}H(\delta F_{\mathrm{bot}}^x-\delta F_{\mathrm{top}}^x)$$

$$=\delta W_{\mathrm{p}}-\sum _i{z}_i{f}_i^x.$$

$$\delta F_{\mathrm{top}}^{\alpha }=-s_0\sum _{k\in \mathrm{top}}{u}_k^{\alpha },\delta F_{\mathrm{bot}}^{\alpha }=-s_0\sum _{k\in \mathrm{bot}}{u}_k^{\alpha }.$$

$${\cal H}^{\prime}_{\rm ih}=\sum_{i}\frac{1}{2m_{i}}|{\mbox{\boldmath$p$}}_{i}|^{2}+\delta U-\gamma_{\rm ex}(t)\delta{\cal A}$$

$$\delta \bar{\mathcal{B}}(t)={\chi }_{BA}(0){\gamma }_{\mathrm{ex}}(t)-{\int }_0^{t}ds{\chi }_{BA}(t-s){\dot{\gamma }}_{\mathrm{ex}}(s),$$

Full text

I Introduction

In glasses, the structural relaxation becomes exceedingly slow at low temperature $T$. The shear modulus $\mu$ is then well-defined for small deformations, though plastic events easily take place with increasing the applied strainLiu ; Lacks ; Maeda ; Yama1 ; Kawasaki ; Malo2 ; Barr3 . It is of great interest how $\mu$ in glasses depends on the disordered particle configuration. On the other hand, in crystals, the microscopic expressions for the elastic moduli can be derived under a homogeneously applied strain (or stress) in equilibriumHoover ; Ray ; Lutsko ; Pablo ; Hess . Such expressions are composed of a positive affine part and a negative nonaffine part, where the latter arises from the correlation of the stress fluctuations. If the moduli are homogeneous, they can be related to the variances of the thermal strain (or stress) fluctuations divided by $k_{B}T$Rahman ; Binder ; Gusev .

In glasses, $\mu$ is expressed in the same form as those in crystals in terms of the particle positions on timescales without plastic eventsMalo ; Malo1 ; Evans ; Wil ; Yoshi ; Ilg ; Barr1 ; Barr2 ; Saw ; Wit ; Wit1 ; Teren ; Za . To derive this expression, Maloney and LemaîtreMalo ; Malo1 examined the local minima of the potential energy under constraint of a fixed mean strain in the periodic boundary condition. Remarkably, the local values of $\mu$ exhibit mesoscopic inhomogeneityPabloP ; Miz1 . In fact, in glasses, the displacements (and suitably defined strains) in glasses are highly heterogeneous on mesoscopic scales under shearBarr1 ; Barr2 ; Malo2 ; Yama1 ; Kawasaki ; Liu ; Maeda ; Barr3 .

The linear response theory in statistical mechanics has a long historyHansen ; Onukibook ; Zwan . On the basis of Onsager’s theoryOnsager , GreenGreen derived time-evolution equations for gross variables, which was rigorously justified by ZwanzigZwan . Green then expressed the transport coefficients in fluids such as the viscosities and the thermal conductivity in terms of the time-correlation functions of the stress and the heat flux, respectively. These expressions also followed from the relaxation behaviors of the time-correlation functions of hydrodynamic variablesHansen ; Kada ; Zwan . KuboKubo studied linear response to mechanical forces, for which the Hamiltonian consists of the unperturbed one $\cal H$ and a small time-dependent perturbation as

[TABLE]

Here, $\gamma_{\rm ex}(t)$ is an applied force and ${\cal A}$ is its conjugate variable. Thus, there was a conceptual difference between the approaches to thermal and mechanical disturbances.

In this paper, we set up a Hamiltonian in the form of Eq.(1) for slight motions of the boundary walls in the film geometry, where the film thickness $H$ is much longer than the particle sizes. Here, $\gamma_{\rm ex}$ is the mean shear strain, and $\cal A$ is given by $H({F}^{x}_{\rm bot}-{F}^{x}_{\rm top})/2$, where ${F}^{x}_{\rm bot}$ and ${F}^{x}_{\rm top}$ are the tangential forces from the bottom and top walls to the particle system. For glasses, we can examine linear response for any variables assuming quasi-equilibrium around a fixed inherent state Malo ; Malo1 ; Evans ; Wil ; Yoshi ; Ilg . This is justified while jump motions do not occur among different inherent statesSastry ; Heuer ; Lacks ; Harro1 . For liquids, our theory yields Green’s expression for the shear viscosity with Kubo’s method in the low-frequency limit. It can further be used to analyze linear response in fluids near a moving wallslip .

In our theory, $\cal A$ in Eq.(1) can be expressed in terms of the particle positions near the walls. However, for nonvanishing $\mu$, the fluctuations of $\cal A$ are significantly correlated with those of all the particle displacements {\mbox{\boldmathu}}_{i} in the film due to the large factor $H$ in its definition. We shall even find a correlation between the fluctuations of ${F}^{x}_{\rm top}$ and ${F}^{x}_{\rm bot}$ proportional to $\mu/H$. On the other hand, in infinite glasses $(H\to\infty)$, the stress pair correlation decays algebraically in space (under the periodic boundary condition in simulations)Lema ; Fuchs ; Seme ; Harro1 ; Egami . We mention a similar effect in polar fluids, where the polarization pair correlation is dipolar in infinite systemsFelder but extends throughout the cell between metallic or polarizable wallsTakae .

We can also study propagation of sounds in glasses as a linear response to a small-amplitude wall motion. In our theory, its time-evolution can be described in terms of the time-correlation functions of the particle displacements {\mbox{\boldmathu}}_{i}(t) and ${\cal A}(0)$, where the quasi-equilibrium average is taken around a fixed inherent state. As in granular materialsJia ; Roux , we shall find rough wave fronts and random scattered waves. It is of general interest how thermal sound waves come into play in the time-correlation functions in films at low $T$.

The organization of this paper is as follows. In Sec. II, we will present the theoretical background of the linear response in glasses with respect to tangential motions of the boundary walls. In Sec.III, the linear response in supercooled and ordinary liquids will be briefly discussed. In Sec.IV, numerical results will be presented to confirm our theory in glasses. Additionally, a random elastic system will be treated in one dimension in Appendix A. Correlations among the displacements and the wall forces will be examined for homogeneous elastic moduli in Appendix B.

II Linear Response at a fixed inherent state

We consider a low-temperature glass composed of two species with particle numbers $N_{1}$ and $N_{2}$. The total particle number is $N=N_{1}+N_{2}$. We write the particle positions as {\mbox{\boldmathr}}_{i}=(x_{i},y_{i},z_{i}) and the momenta as {\mbox{\boldmathp}}_{i}=(p_{i}^{x},p_{i}^{y},p_{i}^{z}). We assume nearly rigid boundary walls at $z=\pm H/2$. These particles are confined in the cell region $-H/2<z<H/2$ along the $z$ axis, but the periodic boundary condition is imposed along the $x$ and $y$ axes with period $L$. The cell volume is $V=HL^{d-1}$. The lengths $H$ and $L$ are much longer than the particle diameters. Our results can be used both in two and three dimensions ($d=2$ and $3$), where the $y$ components are absent for $d=2$. In Appendix A, our theory will be presented in analytic forms in one dimension.

II.1 Applying shear from boundary walls

As illustrated in Fig,1, we induce shear deformations by motions of boundary wallsShiba ; Maeda , which are in the regions $-\ell_{\rm w}<z+H/2<0$ at the bottom and $0<z-H/2<\ell_{\rm w}$ at the top with $\ell_{\rm w}\ll H$. To each layer, $M$ particles are bound by spring potentials \psi(|{\mbox{\boldmathr}}_{k}-{\mbox{\boldmathR}}_{k}|), where {\mbox{\boldmathr}}_{k} are the positions of these particles and {\mbox{\boldmathR}}_{k} are the pinning centers fixed to the layers. We set $N<k\leq N+M$ at the top and $N+M<k\leq N+2M$ at the bottom. We assume the simple harmonic potential,

[TABLE]

Pinning becomes stronger with increasing the coefficient $s_{0}$. The bound particles also belong to either of the first or second species in the bulk and their density is equal to the bulk density, so $M=N\ell_{\rm w}/H$.

Particle pairs $i\in a$ and $j\in b$ (including the bound ones) interact via short-ranged potentials $\phi_{ab}(r_{ij})$, where r_{ij}=|{\mbox{\boldmathr}}_{i}-{\mbox{\boldmathr}}_{j}| and $a$ and $b$ denote the particle species ($1$ or $2$). As a result, the unbound particles do not penetrate into the boundary layers. For simplicity, we write $\phi_{ij}=\phi_{ab}(r_{ij})$ and \psi_{k}=\psi(|{\mbox{\boldmathr}}_{k}-{\mbox{\boldmathR}}_{k}|). At fixed {\mbox{\boldmathR}}_{k}, the total potential energy is given by

[TABLE]

where we sum over all the particles in the first term and the bound ones in the second term $(k>N$). We write the momentum density as {\mbox{\boldmathJ}}=\sum_{i}{\mbox{\boldmathp}}_{i}\delta({\mbox{\boldmathr}}-{\mbox{\boldmathr}}_{i}) using the $\delta$ function. Since the force {\mbox{\boldmathf}}_{i}=-\partial U/\partial{\mbox{\boldmathr}}_{i} on particle $i$ consists of the contributions from the particles and the walls, the time derivative {\dot{\mbox{\boldmathJ}}}=\partial{\mbox{\boldmathJ}}/\partial t is written as

[TABLE]

Here, $\stackrel{{\scriptstyle\leftrightarrow}}{{\Pi}}=\{\Pi_{\alpha\beta}\}$ is the microscopic stress tensorIrving at position $r$ and the second term is the force density from the walls. Hereafter, the over-dot denotes taking the time-derivative, $\alpha$ and $\beta$ represent the Cartesian coordinates, and we set $\nabla_{i}^{\alpha}=\partial/\partial x_{i}^{\alpha}$ and $\nabla_{i}=(\nabla_{i}^{x},\nabla_{i}^{y},\nabla_{i}^{z})$.

We divide the stress tensor into kinetic and potential parts as $\Pi_{\alpha\beta}=\Pi_{\alpha\beta}^{\rm K}+\Pi_{\alpha\beta}^{\rm p}$ withIrving ; Onukibook

[TABLE]

where we sum over all the particles and set $\phi^{\prime}_{ij}=d\phi_{ij}/dr_{ij}$ and $x_{ij}^{\alpha}=x_{i}^{\alpha}-x_{j}^{\alpha}$ (the $\alpha$ component of {\mbox{\boldmathr}}_{ij}={\mbox{\boldmathr}}_{i}-{\mbox{\boldmathr}}_{j}). In Eq.(6), we introduce the Irving-Kirkwood $\delta$ function by

[TABLE]

which satisfies \int d{\mbox{\boldmathr}}{\hat{\delta}}({\mbox{\boldmathr}},{\mbox{\boldmathr}}_{i},{\mbox{\boldmathr}}_{j})=1 and is nonvanishing only when $r$ is on the line segment connecting {\mbox{\boldmathr}}_{i} and {\mbox{\boldmathr}}_{j}. It follows the relation {\mbox{\boldmathr}}_{ij}\cdot\nabla\hat{\delta}({\mbox{\boldmathr}},{\mbox{\boldmathr}}_{i},{\mbox{\boldmathr}}_{j})=\delta({\mbox{\boldmathr}}-{\mbox{\boldmathr}}_{j})-\delta({\mbox{\boldmathr}}-{\mbox{\boldmathr}}_{i}), leading to \sum_{\beta}\nabla_{\beta}\cdot{{\Pi}}^{\rm p}_{\alpha\beta}=\sum_{i,j}\delta({\mbox{\boldmathr}}-{\mbox{\boldmathr}}_{i})\nabla_{i}^{\alpha}\phi_{ij} in Eq.(4). If we integrate $z{\dot{J}}_{\alpha}$ in a region containing all the particles, Eq.(4) gives a useful relation,

[TABLE]

where the right hand side arises from the wall potentials.

We next move the top wall by $\gamma_{\rm ex}H/2$ and the bottom wall by $-\gamma_{\rm ex}H/2$ along the $x$ axis. Here $\gamma_{\rm ex}$ is a small mean shear strain, which can depend on time. Then, in $U$ in Eq.(3), the positions {\mbox{\boldmathR}}_{k} are shifted by $\pm\gamma_{\rm ex}H/2$ along the $x$ axis and $\psi_{k}$ are changed by $\mp(\gamma_{\rm ex}H/2)(\partial\psi_{k}/\partial x_{k})$ to linear order in $\gamma_{\rm ex}$. Thus, $U$ is changed by $-\gamma_{\rm ex}{\cal A}$ with

[TABLE]

where $F_{\rm top}^{\alpha}=-\sum_{k\in{\rm top}}\nabla_{k}^{\alpha}\psi_{k}$ is the force on the particles from the top wall and $F_{\rm bot}^{\alpha}=-\sum_{k\in{\rm bot}}\nabla_{k}^{\alpha}\psi_{k}$ is that from the bottom wall. In Eq.(9), $\cal A$ consists of the contributions from the bound particles ($k>N$) and is amplified by the prefactor $H/2$. Hereafter, for any $a_{k}$, we write $(\sum_{k\in{\rm top}}-\sum_{k\in{\rm bot}})a_{k}=\sum_{k\in{\rm top}}a_{k}-\sum_{k\in{\rm bot}}a_{k}$. Previously, Shiba and one of the present authorsShiba applied a shear flow in the geometry in Fig.1 to binary particle systems in various states.

For $\alpha=x$, the right hand side of Eq.(8) is rewritten as ${\cal A}+\sum_{k>N}\delta z_{k}\nabla_{k}^{x}\psi_{k}$, where $\delta z_{k}=z_{k}\mp H/2$ at the top and the bottom layers, so it is nearly equal to $\cal A$ for $|\delta z_{k}|<\ell_{\rm w}\ll H$. Thus, $\cal A$ can also be expressed as a difference of two bulk integrals,

[TABLE]

Here, the kinetic parts cancel from the relation \int d{\mbox{\boldmathr}}z{\dot{J}}_{z}=d(\sum_{i}z_{i}p_{i}^{x})/dt=\sum_{i}[p_{i}^{z}p_{i}^{x}/m_{i}+z_{i}f_{i}^{x}]. In the second line, $f_{i}^{x}$ is the $x$ component of the force {\mbox{\boldmathf}}_{i}=-\partial U/\partial{\mbox{\boldmathr}}_{i} on particle $i$, so the sum $\sum_{i}z_{i}f_{i}^{x}$ is the off-diagonal virial. Since $U$ includes the wall potentials, the total force on the particles is given by

[TABLE]

Note that we can also move the top wall by $\gamma_{\rm ex}(H/2-a)$ and the bottom wall by $-\gamma_{\rm ex}(H/2+a)$, where $a$ is an arbitrary length. For example, the top wall is at rest for $a=H/2$. Then, ${\cal A}$ in Eq.(9) is changed to ${\cal A}^{\prime}=(H/2+a)F_{\rm bot}^{x}-(H/2-a)F_{\rm top}^{x}$ and $z_{i}$ in Eq.(10) is changed to $z_{i}-a$. Our expressions for $\mu$ will not depend on $a$ (see the sentences below Eqs.(44) and (50) and Eq.(51)).

II.2 Inherent state and quasi-equilibrium in glass

Glasses are nonergodic at low $T$, where the configuration changes are negligible during the observation in the limit of small applied strain. In the corresponding inherent state in the limit $T\to 0$, we write the particle positions as {\mbox{\boldmathr}}_{i}={\mbox{\boldmathr}}_{i}^{\rm ih}=(x^{\rm ih}_{i},y^{\rm ih}_{i},z^{\rm ih}_{i}), for which the mechanical equilibrium (${f}_{i}^{\alpha}=-\nabla_{i}^{\alpha}{U}={0}$) holds for all the particles. Thus, each inherent state is one of the local minima in the many-particle phase space. In addition, the inherent (residual) shear stress \Pi_{\alpha\beta}^{\rm ih}({\mbox{\boldmathr}})=\lim_{T\to 0}\Pi_{\alpha\beta}({\mbox{\boldmathr}}) is highly heterogeneous with long-range correlationsHarro ; Egami ; Lema ; Fuchs ; Seme . In our case, from Eqs.(9) and (10), the total shear stress is related to the force difference as

[TABLE]

while the total applied force {\mbox{\boldmathF}}_{\rm bot}+{\mbox{\boldmathF}}_{\rm top} vanishes as $T\to 0$ from Eq.(11). Fuereder and IlgIlg found a counterpart of Eq.(12) in the periodic boundary condition. See more discussions on \Pi_{zx}^{\rm ih}({\mbox{\boldmathr}}) in Sec.IIH.

We next consider the particle displacements,

[TABLE]

from the inherent positions. To leading order, the effective Hamiltoninan for {\mbox{\boldmathu}}_{i} is of the bilinear form,

[TABLE]

where {\mbox{\boldmathp}}_{i}=m_{i}{\dot{\mbox{\boldmathu}}}_{i} and the potential part $\delta U$ is given by

[TABLE]

The first term depends on $u_{ij}^{\alpha}=u_{i}^{\alpha}-u_{j}^{\alpha}$. We define

[TABLE]

where $\phi_{ij}^{\prime\prime}=d^{2}\phi_{ij}/dr_{ij}^{2}$ and the particle positions are at {\mbox{\boldmathr}}_{i}={\mbox{\boldmathr}}_{i}^{\rm ih}. For simplicity, $x_{i}^{\rm ih}$, $y_{i}^{\rm ih}$, and $z_{i}^{\rm ih}$ will be written as $x_{i}$, $y_{i}$, and $z_{i}$, when confusion will not occur.

The ${\cal H}_{\rm ih}$ describes the local vibrational motions in local potential minima and the collective acoustic modes on larger scales. We assume that the the observation time is much longer than the microscopic times but is much shorter than the structural relaxation time (see Sec.IIIA). Then, quasi-equilibrium should be attained at fixed {\mbox{\boldmathr}}_{i}^{\rm ih}, where ({\mbox{\boldmathu}}_{i},{\mbox{\boldmathp}}_{i}) obey the canonical distribution,

[TABLE]

Hereafter, the thermal average over this distribution will be written as \langle{\cdots}\rangle_{\rm ih}=\int\Pi_{i=1}^{N+2M}d{\mbox{\boldmathp}}_{i}d{\mbox{\boldmathu}}_{i}(\cdots)P_{\rm ih} for a fixed inherent state (isoinherent ensemble). For simplicity, we assume the Gaussian form of $P_{\rm ih}$ to obtain $\langle{u_{i}^{\alpha}u_{j}^{\beta}}\rangle_{\rm ih}\propto k_{B}T$. However, there is no difficulty to include the anharmonic potential terms, which can be important with increasng $T$ (see a remark below Eq.(42))Sastry .

To linear order, the force on particle $i$ is written as

[TABLE]

From Eq.(15) the Hessian matrix $\{h_{ij}^{\alpha\beta}\}$ is symmetric as

[TABLE]

where $\theta_{i-N}$ is 1 for $i>N$ and 0 for $i\leq N$. The last term is the contribution from the bound particles at the walls and is nonexistent in the periodic boundary condition.

The potential energy deviation in Eq.(15) assumes the symmetric bilinear form $\delta U=\sum_{i,j,\alpha,\beta}h_{ij}^{\alpha\beta}u_{i}^{\alpha}u_{j}^{\beta}/2$. Thus, the inverse matrix of $\{h_{ij}^{\alpha\beta}\}$ is related to the displacement variances as

[TABLE]

From $f_{i}^{\alpha}=k_{B}T\partial(\ln P_{\rm ih})/\partial u_{i}^{\alpha}$, any variable $\cal B$ satisfies

[TABLE]

For ${\cal B}=u_{i}^{\alpha}$ we find $\langle{u_{i}^{\alpha}f_{j}^{\beta}}\rangle_{\rm ih}=-k_{B}T\delta_{ij}\delta_{\alpha\beta}.$ See a similar relation in Eq.(69) in equilibrium.

Many authorsSch ; Liu ; Yama ; Gelin ; Barr1 ; Barr2 ; Malo ; Elliott have calculated the vibrational modes in glasses in the periodic boundary condition. We can examine them with boundary walls by including the last term in Eq.(19) (see Fig.2).

II.3 $\Pi_{\alpha\beta}^{\rm p}$ and $\cal A$

around an inherent state

Around an inherent state, the potential part of the stress tensor $\Pi_{\alpha\beta}$ in Eq.(6) is expanded as

[TABLE]

The first term is the inherent stress and the second term is the first-order deviation linear in {\mbox{\boldmathu}}_{i} given by

[TABLE]

In the first term the coefficients $w_{ij}^{{\alpha\beta}\nu}(=w_{ij}^{\beta\alpha\nu}=-w_{ji}^{{\alpha\beta}\nu}$) are defined in the inherent state by

[TABLE]

We derive the second term in Eq.(23) from the deviation of $\hat{\delta}$ in Eq.(7) using x_{ij}^{\beta}{\mbox{\boldmathu}}_{i}\cdot\nabla_{i}{\hat{\delta}}=-{\mbox{\boldmathu}}_{i}\cdot\nabla[(x_{\beta}-x_{j}^{\beta}){\hat{\delta}}], which vanishes upon space integration. In this paper, we consider the space integral of \delta\Pi_{zx}^{\rm p}({\mbox{\boldmathr}}) written as

[TABLE]

where we can replace $u_{ij}^{\nu}=u_{i}^{\nu}-u_{j}^{\nu}$ in the first line by $2u_{i}^{\nu}$ to obtain the second line and we define the coefficients,

[TABLE]

From Eqs.(9), (10), and (25) the deviation of $\cal A$ is written to linear order in the following two forms,

[TABLE]

Here, $\delta F_{\rm top}^{\alpha}$ is the deviation of the force from the top wall and $\delta F_{\rm bot}^{\alpha}$ is that from the bottom wall, so

[TABLE]

We can also write $\delta{\cal A}=s_{0}\sum_{k>N}z_{k}u_{k}^{x}$ for $H\gg\ell_{\rm w}$.

II.4 Linear response in quasi-equilibrium

For a fixed inherent state, we use Kubo’s methodKubo in the linear response theory for the perturbed Hamiltonian,

[TABLE]

around the quasi-equilibrium distribution $P_{\rm ih}$ in Eq.(17), where ${\delta U}$ and ${\delta{\cal A}}$ are given by Eqs.(15) and (27). The equations of motion are \dot{\mbox{\boldmathp}}_{i}={\mbox{\boldmathf}}_{i}+\gamma_{\rm ex}\partial(\delta{\cal A})/\partial{\mbox{\boldmathu}}_{i} including the perturbation. As a result, the phase-space distribution slightly deviates from $P_{\rm ih}$. For any variable ${\cal B}$, we consider its time-dependent average over this perturbed distribution. Its deviation due to $\gamma_{\rm ex}$ is given by

[TABLE]

where ${\dot{\gamma}}_{\rm ex}(t)=d\gamma_{\rm ex}(t)/dt$ and we assume ${\gamma}_{\rm ex}(t)=0$ for $t<0$. We define the response function,

[TABLE]

where the time-evolution is governed by the unperturbed Hamiltonian ${\cal H}_{\rm ih}$ in Eq.(14) (see Sec.IIG). The static susceptibility is given by the equal-time correlation in Eq.(30). Use of the second line of Eq.(27) gives

[TABLE]

where the first term arises from the affine change in $\cal B$ and the second term is due to the correlation between $\cal B$ and the deviation $\delta W_{\rm p}$ in the total shear stress.

In particular, for a stepwise shear strain $\gamma_{\rm ex}(t)=\gamma_{0}\theta(t)$ (being 0 for $t<0$ and $\gamma_{0}$ for $t>+0$), we find

[TABLE]

For ${\cal B}=\delta W_{\rm p}$, we obtain the stress relaxation function $G_{\rm ih}(t)=-\delta{\bar{W}}_{\rm p}(t)/\gamma_{0}$ for $t>0$. From Eq.(33) we find

[TABLE]

where $\mu$ is the average shear modulus in Eq.(38) below, so $G_{\rm ih}(0)=0$. Our $\chi_{BA}(t)$ and $G_{\rm ih}(t)$ exhibit oscillatory behavior and can be used only for $t$ before appreciable structural relaxation (see Fig.7). However, in the literature of rheologySaw ; Doi ; Wit1 , the stress relaxation function is the sum of $\mu$ and the stress time-correlation function in Eq.(61) (divided by $k_{B}T$).

II.5 Displacement and stress in steady shear strain

Let us make further calculations for a steady shear strain $\gamma_{\rm ex}$. For ${\cal B}=u_{i}^{\alpha}$, the first line of Eq.(27) gives the average dispacements. With the aid of Eq.(28) we find

[TABLE]

Here, all $u_{i}^{\alpha}$ are coupled with the force difference $\delta F_{\rm bot}^{x}-\delta F_{\rm top}^{x}$. See Sec.IIF for more discussions on $\delta F_{\rm bot}^{x}$ and $\delta F_{\rm top}^{x}$. Furthermore, Eq.(32) gives another expression,

[TABLE]

which consists of the affine and nonaffine partsMalo .

For {\cal B}=\delta\Pi_{zx}^{\rm p}({\mbox{\boldmathr}}), we define the local shear modulus by

[TABLE]

This {\hat{\mu}}({\mbox{\boldmathr}}) has the particle discreteness in an inherent state, so we need to integrate it in small squares or cubes in the cell to detect elastic heterogeneityMiz1 ; PabloP . The space average of {\hat{\mu}}({\mbox{\boldmathr}}) in the whole cell is simpler as

[TABLE]

Then, as $T\to 0$, Eqs.(35) and (38) yield

[TABLE]

which arises from the correlations between the bound and unbound particles.

From Eq.(36) we also find the well-known expression,

[TABLE]

where $\mu_{\infty}=\sum_{i}w_{i}^{x}z_{i}/V$ is the affine contribution,

[TABLE]

This $\mu_{\infty}$ was first introduced for fluids (see Eq.(57))Zwanzig . As $T\to 0$, $\langle{\mu_{\infty}}\rangle_{\rm ih}$ is equal to $\mu_{\infty}$ at the inherent positions {\mbox{\boldmathr}}_{i}={\mbox{\boldmathr}}_{i}^{\rm ih}, while the nonaffine part is proportional to the stress variance and is negative. As $T\to 0$, we find

[TABLE]

We can use Eq.(40) even if we include the anharmonic potential terms in the average $\langle{\cdots}\rangle_{\rm ih}$ (see a comment below Eq.(17)). For large $V$, the bulk contributions dominate over the surface ones in Eqs.(40)-(42). See Appendix A for counterparts of Eqs.(40)-(42) in one dimension.

From Eqs.(27) and (30) the average of $\delta{\cal A}$ is given by

[TABLE]

where we use $\langle{\delta{\cal A}\sum_{i}z_{i}f_{i}^{x}}\rangle_{\rm ih}/k_{B}T=-s_{0}\sum_{i>N}z_{i}^{2}$. Since $\langle{\delta{\cal A}^{2}}\rangle_{\rm ih}>0$, we require $s_{0}>2\mu V/H^{2}M$, which is well satisfied in our simulation in Sec.IV. From Eqs.(27), (28), and (43) we can also express $\mu$ as

[TABLE]

in terms of the variances of the bound particles. This relation is invariant with respect to the coordinate shift along the $z$ axis ($z_{i}\to z_{i}-a$ and $z_{j}\to z_{j}-a$) (which can be proved from Eqs.(47) and (48) below).

In Eq.(15) we replace $u_{k}^{x}$ by $u_{k}^{x}\mp H\gamma_{\rm ex}/2$ for $k>N$ to obtain the change in the potential energy deviation,

[TABLE]

up to of order $\gamma_{\rm ex}^{2}$, Minimization of $\delta U^{\prime}$ with respect to $u_{i}^{\alpha}$ gives $f_{i}^{\alpha}+\gamma_{\rm ex}\partial(\delta{\cal A})/\partial u_{i}^{\alpha}=0$. This leads to Eqs.(35) and (36) with the aid of the first and second lines of Eq.(27). Note that this derivation of Eq.(36) is equivalent to that of Maloney and LemaîtreMalo ; Malo1 . Let us then calculate the minimum value of $\delta U^{\prime}$ by seting $u_{i}^{\alpha}$ equal to ${\bar{u}}_{i}^{\alpha}$ in Eq.(35) and $\delta{\cal A}$ equal to $\delta{\bar{\cal A}}$ in Eq.(43). In terms of $\mu$ in Eq.(38), it is simply written as

[TABLE]

which is not obvious for inhomogeneous glassy systems.

We note the following. (i) For crystals, the shear moduli are expressed in the same form as in Eq.(40), where the thermal average is taken over the displacement fluctuations in a given crystal stateHess ; Ray ; Lutsko ; Hoover . For glasses, the average is taken over one inherent statesMalo ; Malo1 ; Evans ; Yoshi ; Ilg ; Saw . (ii) To examine elastic inhomogeneity in glasses, some authorsMiz1 ; PabloP divided the cell into small regions and integrated the average of the first term in Eq.(23) in each region (where $\hat{\delta}$ can be integrated in a simple form).

II.6 Fluctuations of forces from walls

We have introduced the forces from the walls to the particles in Eq.(9). We here examine the equal-time correlations among their thermal fluctuations and $u_{i}^{\alpha}$. From Eqs.(11) and (21) and $\sum_{j}h_{ij}^{\alpha\beta}=-s_{0}\theta_{i-N}\delta_{\alpha\beta}$, the total force deviation $\delta F_{\rm tot}^{\alpha}\equiv\delta F_{\rm top}^{\alpha}+\delta F_{\rm bot}^{\alpha}$ satisfies

[TABLE]

From Eq.(11) the left hand side of Eq.(47) is written as $\langle{u_{i}^{\alpha}\sum_{j}{\dot{p}}_{j}^{\beta}}\rangle_{\rm ih}=-\langle{{\dot{u}}_{i}^{\alpha}\sum_{j}p_{j}^{\beta}}\rangle_{\rm ih}$; then, Eq.(47) is obvious. From Eqs.(25) and (26) $\delta W_{\rm p}$ and $\delta F_{\rm tot}^{\alpha}$ are orthogonal as

[TABLE]

On the other hand, the $x$-component of the force difference $\delta F_{\rm top}^{x}-\delta F_{\rm bot}^{x}$ is proportional to $\delta{\cal A}$ as in Eq.(27), so its relations follow from Eqs.(35) and (43). Using Eq.(48) also, we find the cross correlation,

[TABLE]

See Appendices A and B for counterparts of Eq.(50) in simpler situations. Here, we argue that Eq.(50) is general for elastic films. Let us move the bottom layer by $-H\gamma_{\rm ex}$ with the top layer kept at rest (see the last paragraph of Sec.IIA). Then, $\delta{\cal A}$ is changed to $\delta{\cal A}^{\prime}=H\delta F_{\rm bot}^{x}$ in Eq.(29) and the average force from the top wall to the particles is given by $\langle{\delta F_{\rm top}^{x}\delta{\cal A}^{\prime}}\rangle_{\rm ih}/k_{B}T=\mu V/H$ in equilibrium, which coincides with Eq.(50). See Eq.(82) and Fig.8 for the time-correlation of the wall forces. In accord with this argument, Eqs.(38) and (49) give

[TABLE]

where $\delta{F}_{\rm bot}^{x}$ can be replaced by $-\delta{F}_{\rm top}^{x}$.

II.7 Linear dynamics around a fixed inherent state

Around a fixed inherent state, the dynamic equations follow from ${\cal H}_{\rm ih}$ in Eq.(14). They are rewritten as

[TABLE]

For $m_{1}\neq m_{2}$, we use the following scale changesElliott ; Yama ,

[TABLE]

where ${\hat{h}}_{ij}^{\alpha\beta}$ is the modified Hessian matrix. Its inverse is given by ${({\hat{h}}^{-1})}_{ij}^{\alpha\beta}=\sqrt{m_{i}m_{j}}{(h^{-1})}_{ij}^{\alpha\beta}$. We introduce the $d(N+2M)$ dimensional eigenvectors $e^{\lambda}_{i\alpha}$ satisfying

[TABLE]

where $0<\omega_{1}\leq\omega_{2}\leq\cdots$. The $e_{i\alpha}^{\lambda}$ are normalized as $\sum_{\lambda}e^{\lambda}_{i\alpha}e^{\lambda}_{j\beta}=\delta_{\alpha\beta}\delta_{ij}$ and $\sum_{i,\alpha}e^{\lambda}_{i\alpha}e^{\sigma}_{i\alpha}=\delta_{\lambda\sigma}$. Projection of ${\hat{u}}_{i}^{\alpha}$ on the eigenmodes yields the variables $s_{\lambda}$ as

[TABLE]

The potential energy deviation in Eq.(15) is written as

[TABLE]

from which we find $\langle{s_{\lambda}s_{\sigma}}\rangle_{\rm ih}=k_{B}T\delta_{\lambda\sigma}/\omega_{\lambda}^{2}$.

Now, from Eqs.(52)-(55), we obtain the dynamic equations $d^{2}s_{\lambda}(t)/dt^{2}=-\omega_{\lambda}^{2}s_{\lambda}(t)$. These are solved to give $s_{\lambda}(t)=s_{\lambda}(0)\cos(\omega_{\lambda}t)+{\dot{s}}_{\lambda}(0)\sin(\omega_{\lambda}t)/\omega_{\lambda}.$ with ${\dot{s}}_{\lambda}(0)=\sum_{i,\alpha}e^{\lambda}_{i\alpha}\sqrt{m_{i}}{\dot{u}}_{i}^{\alpha}(0)$ being linear in the velocities at $t=0$. Thus, averaging over $P_{\rm ih}$ yields

[TABLE]

For any variable $\delta{\cal B}(t)=\sum_{i,\alpha}b_{i}^{\alpha}u_{i}^{\alpha}(t)$ linear in $u_{i}^{\alpha}(t)$, we can express it as $\delta{\cal B}(t)=\sum_{\lambda}Z_{B}^{\lambda}s_{\lambda}(t)$, where

[TABLE]

Furthermore, if $\delta{\cal A}(t)=\sum_{\lambda}Z_{A}^{\lambda}s_{\lambda}(t)$, Eq.(57) gives

[TABLE]

At $t=0$, Eq.(59) is the equal-time correlation from $\sum_{\lambda}e^{\lambda}_{i\alpha}\omega_{\lambda}^{-2}e^{\lambda}_{j\beta}={({\hat{h}}^{-1})}_{ij}^{\alpha\beta}$. The expressions (57) and (59) exhibit oscillation even at long times without dissipationGelin ; Elliott . In particular, if $\delta{\cal B}=u_{i}^{\alpha}$ and $\delta{\cal A}$ is given by Eq.(27), Eqs.(35), (55), and (59) lead to another expression for the average induced displacements,

[TABLE]

II.8 Stress time-correlation in simulation

The stress time-correlation function has been calculated via molecular dynamics simulationOnukibook ; Harro ; Zwan ; Hansen ; Green . Using the integral of the shear stress W(t)=\int d{\mbox{\boldmathr}}\Pi_{zx}({\mbox{\boldmathr}},t) at time $t$ for large systems, we express it as

[TABLE]

As in usual simulations, $\langle{\cdots}\rangle_{\rm sim}$ includes the average over a long time interval of the initial time $t_{0}\in[0,t_{\rm sim}]$. It can also be over many simulation runs or over many inherent states in glasses. For a suitable ensemble, the integral of the inherent stress W^{\rm ih}=\int d{\mbox{\boldmathr}}\Pi^{\rm ih} can obey a distribution with $\langle{W^{\rm ih}}\rangle_{\rm sim}=0$ without applied strainHarro ; Ilg . Here, if $t$ and $t_{\rm sim}$ are sufficiently long at not very small $T$, $C(t)$ slowly decays to zero due to the configuration changes (see Sec,IIIA).

It is well-known that $C(t)$ decays from its initial value $C(0)$ to a well-defined plateau value $C_{\rm pl}$ after a microscopic time at low $T$. In our simulation in Sec.IV, this will be the case after averaging over $t_{0}$ even in a single run. Here, $C(0)$ and $C_{\rm pl}$ do not depend on the system size for large $V$. The $C_{\rm pl}$ is nearly independent of $T$ and is related to the inherent stress asIlg ; Harro ; Yoshi ; Evans ; Wil

[TABLE]

From Eq.(12) we also find

[TABLE]

On the other hand, the kinetic stress and the potential stress deviation $\delta W_{\rm p}(t)$ decay rapidly in Eq.(61). Thus,

[TABLE]

where $n(k_{B}T)^{2}$ is the kinetic contribution with $n=N/V$ and $\langle{{\delta W_{\rm p}}^{2}}\rangle_{\rm ih}$ is given in Eq.(42) for each inherent state.

From Eqs.(40) and (64) the average of $\mu$ is written as

[TABLE]

which holds at finite $T$. Since Eq.(11) indicates $F_{\rm bot}^{x}=-F_{\rm top}^{x}$ at $T=0$, Eqs.(62)-(65) yield the wall-force correlation including the inherent contribution,

[TABLE]

III Vanishing of shear modulus

III.1 Supercooled liquids

In supercooled liquids, the configuration changes occur appreciably on timescales longer than the bond breakage time $\tau_{b}$Kawasaki ; Yama1 ; Shiba . In each plastic event, some bonds are broken and some particles jump over distances longer than their diameters. As a result, the diffusion constants of the two componentsKawasaki become proportional to $\tau_{b}^{-1}$ with the coefficients independent of $T$. Therefore, the quasi-equilibrium distribution $P_{\rm ih}$ in Eq.(17) can be used on timescales shorter than $\tau_{b}$, where the mean-square displacements are smaller than the square of the particle diameters. At higher $T$, the second term in Eq.(65) increases such that $\langle{\mu}\rangle_{\rm sim}\to 0$ at a transition temperatureWil ; Evans ; Yoshi ; Saw ; Ilg ; Teren .

For shear rate ${\dot{\gamma}}_{\rm ex}$, we have a Newtonian regime for ${\dot{\gamma}}_{\rm ex}\tau_{b}<1$ and a shear-thinning regime for ${\dot{\gamma}}_{\rm ex}\tau_{b}>1$Kawasaki ; Yama1 ; Shiba , where the Newtonian viscosity is of order $\mu\tau_{b}$ and the nonlinear one is of order $\mu/{\dot{\gamma}}_{\rm ex}$. Note that the Green-Kubo formula for the former depends on the ensemble for deep supercoolingHarro ; Yip , where the average $\langle{\cdots}\rangle_{\rm sim}$ in Eq.(61) has to be still in some limited phase-space region.

III.2 Liquids

In liquids at higher $T$, thermal equilibration is rapidly achieved within experimental times. We can use the linear response theory around the equilibrium distribution $P_{\rm eq}\propto\exp(-{\cal H}/k_{B}T)$ with {\cal H}=\sum_{i}|{\mbox{\boldmathp}}_{i}|^{2}/2m_{i}+U. Here, the equilibrium average over $P_{\rm eq}$ is written as $\langle{\cdots}\rangle_{\rm eq}$.

If the boundaries are slightly moved, the perturbed Hamiltonian is ${\cal H}^{\prime}={\cal H}-\gamma_{\rm ex}(t){\cal A}$, where $\cal A$ is given in Eq.(9) with $\langle{{\cal A}}\rangle_{\rm eq}=0$. For any variable ${\cal B}$, its average response is thus expressed in the form of Eq.(30) with

[TABLE]

If {\cal B}(t)=J_{x}({\mbox{\boldmathr}},t)=\sum_{i}p_{i}^{x}(t)\delta({\mbox{\boldmathr}}-{\mbox{\boldmathr}}_{i}(t)) in Eq.(67), the response depends on $z$ and $t$. The resultant boundary flow profile will be investigated in future.

In liquids, the static shear modulus $\mu$ vanishes from $\langle{W{\cal A}}\rangle_{\rm eq}=0$, where W=\int d{\mbox{\boldmathr}}\Pi_{\rm xy} is the total shear stress. Using $\mu_{\infty}$ in Eq.(41) we rewrite it as

[TABLE]

where the left hand side is called the high-frequency shear modulusZwanzig . We can prove Eq.(68) using the relationZwanzig ,

[TABLE]

For a small oscillatory shear rate ${\dot{\gamma}}_{\rm ex}={\dot{\gamma}}_{0}\cos(\omega t)$, the average shear stress is the real part of $-{\dot{\gamma}}_{0}\eta^{*}(\omega)e^{{\rm i}\omega t}$. From Eq.(10) we find the complex shear viscosity,

[TABLE]

where {\cal G}_{x}(t)=\int d{\mbox{\boldmathr}}zJ_{x}({\mbox{\boldmathr}},t)=\sum_{j}z_{j}(t){p}_{j}^{x}(t). At $\omega=0$, the Green-Kubo formula surely holds for the viscosity $\eta=\eta^{*}(0)$ in the bulk. We can also use the argument below Eq.(50) to fluids. For liquids, Eq.(30) gives

[TABLE]

which should be compared with Eq.(50) for elastic films.

IV Numerical results in two dimensions

IV.1 Simulation method

We now present numerical results in two dimensions ($d=2$) in the $x$-$z$ plane. We have briefly described our model at the beginning of Sec.II. Our system is composed of two particle species with numbers $N_{1}=N_{2}=N/2=2000$ in the cell. The pairwise potentials are given by

[TABLE]

where $a$ and $b$ represent the particle species and we introduce $\epsilon$, $\sigma_{1}$, and $\sigma_{2}=1.4\sigma_{1}$. The potentials $\phi_{ab}(r)$ vanish for $r\geq r_{\rm c}=2.25(\sigma_{a}+\sigma_{b})$. The mass ratio is $m_{2}/m_{1}=1.96$. The cell lengths are $H=L=70.2\sigma_{1}$. The average particle density is $n=N/LH=0.81\sigma_{1}^{-2}$ and the average mass density is $\rho=1.2m_{1}\sigma_{1}^{-2}$. Hereafter, we will measure space, time, and temperature in units of $\sigma_{1}$, $t_{0}=\sigma_{1}(m_{1}/\epsilon)^{1/2}$, and $\epsilon/k_{B}$, respectively.

As in Fig.1, we attach two boundary layers with thickness $\ell_{\rm w}=H/16$. Each layer contains $M=250$ particles bound to pinning points {\mbox{\boldmathR}}_{j} by the potential (2). The total particle number is $N+2M=4500$. The {\mbox{\boldmathR}}_{j} were particle positions in a liquid stateShiba . The spring constant is chosen to be large at $s_{0}=200$. Then, the displacements of the bound particles from their inherent positions $u_{k}^{\alpha}$ ($k>N$) undergo thermal motions with amplitudes of order $(T/s_{0})^{1/2}=0.07T^{1/2}\ll 1$. Under applied strain, the motions of the bound particles are very small $(\propto s_{0}^{-1}$), which become even smaller with increasing the distance from the cell region. Our walls are nearly rigid and their surfaces are rough as in Fig.1, so slip motions do not occur for small wall motions.

We followed the following steps. (i) We started with a liquid at high $T$, lowered $T$ to 0.01 without crystallization, and waited for a time of $10^{3}$ to realize a glassy state, where we attached Nosé-Hoover thermostats in the cell and in the boundary layersShiba . (ii) We carried out some simulation runs at $T=0.01$ removing the thermostats. Results from these runs will be presented in Fig.3(c) and Fig.7(a). (iii) To seek an inherent state, we further cooled our system down to $T=10^{-5}$ keeping the thermostats without applied strain. Then, the particle positions became frozen. We use the resultant inherent state in Figs.2-9 without taking the ensemble average.

IV.2 Inherent state and eigenmodes

First, we examine the validity of Eqs.(62) and (63) in the limit $T\to 0$. Indeed, we obtain $(W^{\rm ih})^{2}/V=31.1$ and $(H^{2}/4V){(F_{\rm bot}^{x}-F_{\rm top}^{x})^{2}}=32.2$, while we find the plateau value $C_{\rm pl}=27.7$ from a simulation run at $T=0.01$. These three values are fairly close even for a single inherent state. For the inherent particle positions, we calculated the coefficients $w_{i}^{\alpha}$ in Eq.(26), the Hessian matrix $h_{ij}^{\alpha\beta}$ in Eq.(19), and its inverse $(h^{-1})_{ij}^{\alpha\beta}$ satisfying Eq.(20). The product of these two matrices is very close to the unit matrix $\delta_{\alpha\beta}\delta_{ij}$ with differences of order $10^{-16}$ for all the elements. The shear modulus $\mu$ is then given by $18.1$ from Eq.(39), $15.4$ from Eq.(40), $17.4$ from Eqs.(44) and (50), while it is $16.0$ from the stress-strain relation obtained from a simulation in Fig.3(c).

In Fig.2, we display the first eight eigenvectors $e_{i\alpha}^{\lambda}$, where the displacements of the bound particles are very small and are invisible here. The lowest eigenfrequency is $\omega_{1}=0.174$. This corresponds to the first transverse sound mode $e_{i\alpha}^{1}$, which is roughly proportional to $\delta_{\alpha x}\sqrt{m_{i}}\cos(\pi z_{i}/H)$ $(i\leq N)$. If we set $\omega_{1}=\pi c_{\perp}/H$, we obtain the transverse sound speed $c_{\perp}=3.9$. The same speed follows from $c_{\perp}=(\mu/\rho)^{1/2}$ for $\mu=16$ and $\rho=1.2$. The first longitudinal mode appears at $\lambda=5$ with $\omega_{5}=0.481$. We also obtain strongly localized modes for $\lambda=6$ and $7$, where clusters of large-amplitude oscillation are weakly connected to the bulkSch .

We project the variables $\delta W_{\rm p}$ in Eq.(25) and $\delta{\cal A}$ in Eq.(27) on the eigenmodes as $\delta W_{\rm p}=\sum_{\lambda}Z_{W}^{\lambda}s_{\lambda}$ and $\delta{\cal A}=\sum_{\lambda}Z_{A}^{\lambda}s_{\lambda}$ (see Sec.IIG). Here, Eq.(58) gives

[TABLE]

From Eq.(38), $\mu$ is expressed as

[TABLE]

which is calculated to be 17.7. In Table 1, we give $\omega_{\lambda}$, $Z_{W}^{\lambda}$, and $Z_{A}^{\lambda}$ ($\lambda\leq 8$) for the inherent state under consideration, where $Z_{A}^{\lambda}$ have large amplitudes as compared to $Z_{W}^{\lambda}$. In fact, we find $\langle{(\delta{\cal A})^{2}}\rangle_{\rm ih}/k_{B}T=\sum_{\lambda}(Z_{A}^{\lambda})^{2}/\omega_{\lambda}^{2}=2.5\times 10^{4}V$ from Eq.(43), while we calculate $\langle{(\delta{W}_{\rm p})^{2}}\rangle_{\rm ih}/k_{B}T=\sum_{\lambda}(Z_{W}^{\lambda})^{2}/\omega_{\lambda}^{2}=21.2V$. Notice that Eq.(26) also gives $Z_{W}^{\lambda}=-\sum_{i,j,\alpha}w_{ij}^{zx\alpha}E_{ij\alpha}^{\lambda}/2$, where $E_{ij\alpha}^{\Lambda}\equiv e_{i\alpha}^{\lambda}/{\sqrt{m_{i}}}-e_{j\alpha}^{\lambda}/{\sqrt{m_{j}}}$ are small for most adjacent $i$ and $j$. See Eq.(B8) for the eigenmode projection of $\delta{\cal A}$ in the continuum elasticity.

IV.3 Nonaffine displacements at static strain

In Fig.3, we display the normalized nonaffine displacements $({\bar{u}}_{i}^{x}/\gamma_{\rm ex}-z_{i},{\bar{u}}_{i}^{z}/\gamma_{\rm ex})$ for all the particles. They are highly heterogeneous on large scales for the unbound particles, as has been reported in the literaturePabloP ; Malo ; Malo2 ; Barr1 ; Barr2 ; Miz1 ; Liu ; Wit ; Wit1 . Those of the bound particles are small and invisible. To confirm the validity of our theory, we calculated these results (a) from the wall-particle correlations in Eq.(35), (b) from the particle-particle correlations in Eq.(36), and (c) from a single run of molecular dynamics simulation of applying a strain of $\gamma_{\rm ex}=10^{-3}$ at $T=0.01$. In (a) and (b), we use the inverse Hessian matrix and the results are the averages over $P_{\rm ih}$ in Eq.(17). In (c), use is made of the common inherent state, the time average is over a time interval of $t_{\rm sim}=5000$, and there is no irreversible motion. We can see good agreement of the results in (a), (b), and (c) from the three methods, which supports our theory. In particular, those in (a) and (c) are very close.

In Fig.4(a), we present the full normalized displacements $({\bar{u}}_{i}^{x}/\gamma_{\rm ex},{\bar{u}}_{i}^{z}/\gamma_{\rm ex})$ using Eq.(35), whose affine parts are conspicuous near the walls. Here, Eq.(60) gives

[TABLE]

which is equivalent to the second relation in Eq.(73). If we set ${{\bar{u}}_{i}^{\alpha}}/{\gamma_{\rm ex}}=z_{i}\delta_{\alpha x}$ in Eq.(75), we obtain the affine part of $Z_{A}^{\lambda}$. For $\lambda=1$, it is $-0.54$, so $Z_{A}^{1}$ mostly consisits of the nonaffine part since $Z_{A}^{1}\cong 38$ (see Eq.(B8)).

IV.4 Space-dependent dynamics

Next, we examine space-time-dependent effects. For $\delta{\cal B}=u_{i}^{\alpha}$ in Eq.(59), we define the response functions,

[TABLE]

for all the particles. At $t=0$, we have $\chi^{\alpha}_{i}(0)={\bar{u}}_{i}^{\alpha}/\gamma_{\rm ex}$ for a static strain in Fig.4(a). We also show $\chi_{i}^{\alpha}(t)$ at $t=200\pi/\omega_{1}$ in Fig.4(b), which retains no affine part but still keeps some space correlations. For a stepwise strain $\gamma_{\rm ex}(t)=\gamma_{0}\theta(t)$, which is zero for $t<0$ and is a constant $\gamma_{0}$ for $t>0$, Eq.(33) yields the subsequent evolution,

[TABLE]

In Fig.5, we show $\chi_{i}^{\alpha}(t)$ in the upper panels and $\chi_{i}^{\alpha}(0)-\chi_{i}^{\alpha}(t)$ in the lower panels for $\omega_{1}t/\pi=1/8,1/4,1/2$, and 1, In the initial stage, disturbances advance from the walls with the speed $c_{\perp}$ without noticeable changes in the center region. In (a)-(c), the initial affine correlations in Fig.4(a) disappear from the walls. In (a’)-(c’), shock-like transverse sounds propagate from the walls. Their fronts are irregular due to random scattering. In (d) and (d’), the sounds from the walls encounter at the center.

In Fig.6, we display the laterally averaged profiles,

[TABLE]

where $\theta(p)$ is the step function being 1 for $p>0$ and 0 for $p\leq 1$. We here remove the glassy irregularities to examine the acoustic behavior and the boundary relaxation along the $z$ axis. Setting $\Delta z=2.5$, we plot (a) ${\bar{\chi}}^{x}(z,t)$ and (b) ${\bar{\chi}}^{x}(z,0)-{\bar{\chi}}^{x}(z,t)$, where $\omega_{1}t/\pi=t/18=0,1/16,1/8,1/4,$ and $1/2$. In (a), the initial affine correlation soon disappears near the walls, whose timescale $t_{\rm w}$ is about $2$. In (b), the boundary values of ${\bar{\chi}}^{x}(z,0)-{\bar{\chi}}^{x}(z,t)$ at $z\cong\pm H/2$ change from 0 to the static values ($\sim\pm 40$) on the time $t_{\rm w}$. In the initial stage $t\ll t_{\rm a}=H/c_{\perp}\cong 18$, the expanding sound from each wall is of the form $g(t-\ell/c_{\perp})$, where $\ell$ is the distance from the wall and $g(t)$ is the boundary value being zero for $t\leq 0$. Since $t_{\rm w}\ll t_{\rm a}$, a shock wave is produced from each wall, whose front has a thickness of order $c_{\perp}t_{\rm w}\sim 10$. In (c) and (d), we show the profiles at $\omega_{1}t/\pi=1$ and 200, where complex oscillations still remain nonvanishing.

In agreement with our results, Jia et al.Jia observed a coherent ballistic pulse and speckle-like signals in sound propagation in a granular matter. Their observation was later reproduced in a simulation of granular matterRoux .

IV.5 Time-correlation functions

In Fig.7(a), we examine the stress time-correlation function at a fixed inherent state. We write it as

[TABLE]

In (a), this function decays to nearly zero on a rapid timescale of 0.2. Its oscillatory behavior is suppressed because of relatively small $Z_{W}^{\lambda}$ for not large $\lambda$ in Table 1, where the first term in Eq.(79) is $0.013\cos(\omega_{1}t)$. Notice that the large-scale sound modes do not contribute significantly to $\delta W_{\rm p}(t)$, because the difference $u_{ij}^{\alpha}=u_{i}^{\alpha}-u_{j}^{\alpha}$ appears for adjacent $i$ and $j$ in the first line of Eq.(25).

We also calculated the stress time-correlation function $C(t)$ in Eq.(61) at $T=0.01$ in a single run of molecular dynamic simulation with the same inherent state. Here, the time average was taken over the simulation time $t_{\rm sim}=5\times 10^{4}$. We used the stress integral $W(t)$ in Eq.(61) including the contributions from the bound particles, so there should be some boundary effect. Since Eq.(64) is predicted, Fig.7(a) gives

[TABLE]

where $C_{\rm pl}=27.7$ is the plateau value. In (a), its initial value $C^{\prime}(0)=33.1$ is larger than $C_{\rm ih}(0)=21.2$, but the two curves fairly agree for $t\gtrsim 0.1$.

In Fig.7(b), we plot the stress relaxation function $G_{\rm ih}(t)$ in Eq.(34), which can be expressed as

[TABLE]

In (b), this function exhibits complex oscillatory behavior with timescales shorter than $t_{\rm a}=c_{\perp}/H\sim 18$.

In Fig.8, we examine the time-correlation function of the forces from the walls. From Eq.(50), it is scaled as

[TABLE]

This function can be expressed in the form of Eq.(59) from Eqs.(28) and (58). Starting with ${\cal F}_{\rm w}(0)=\mu\cong 17$. ${\cal F}_{\rm w}(t)$ is nearly constant for $t<t_{\rm a}\sim 18$. Around $t\sim t_{\rm a}$, the sound from the bottom is reflected in the reverse direction at the top, resulting in a large drop in $\delta F_{\rm top}(t)$. For $t\gtrsim t_{a}$, it largely fluctuates due to scattered waves. See Appendix B for ${\cal F}_{\rm w}(t)$ for homogeneous $\mu$.

From the argument below Eq.(50) we recognize that ${\cal F}_{\rm w}(t)$ in Eq.(82) is measurable experimentally. That is, at $t=0$, we move the bottom layer by $-H\gamma_{0}$ in a stepwise manner keeping the top layer at rest. Then, from Eq.(33), the average force from the top wall to the particles is given by $\gamma_{0}[{\cal F}_{\rm w}(0)-{\cal F}_{\rm w}(t)]$ per unit area at time $t(>0)$. As in Fig.8, it should vanish before arrival of the sounds.

Somfai et al.Roux calculated a force signal on a wall after emission of a pulse strain from the opposite wall in a granular model. Their signals without damping resemble ours in Fig.8, but they decay to zero with increasing viscous dissipation. Wittmer et al.Wit1 introduced friction in the stress relaxation in a random elastic network.

IV.6 Resonance in periodic straining

We also suppose application of a sinusoidal strain $\gamma_{\rm ex}(t)=\gamma_{0}\theta(t)\sin(\omega t)$, which is zero for $t<0$. For $t>0$, substitution of this strain into Eq.(30) yields the average displacements ${\bar{u}}_{i}^{\alpha}(t)=\sum_{\lambda}e_{i\alpha}^{\lambda}{{\bar{s}_{\lambda}}(t)}/\sqrt{m_{i}}$. Here, ${{\bar{s}_{\lambda}}(t)}$ is the time-dependent average of the mode amplitude $s_{\lambda}(t)$ in Eq.(55). Then, ${{\bar{s}_{\lambda}}(t)}$ is a mixture of two oscillations,

[TABLE]

The term with the intrinsic frequency $\omega_{\lambda}$ remains nonvanishing, but it decays to zero in the presence of dissipation. Here, ${\bar{s}_{1}}(t)$ for the first mode grows as $(\omega_{1}^{2}-\omega^{2})^{-1}$ with increasing $\omega$ towards $\omega_{1}$. Thus, in Fig.9, ${\bar{s}_{1}}(t)$ is shown to increase with increasing $\omega$. However, as $\omega$ approaches $\omega_{1}$, the local strain increases near the walls and the linear theory becomes invalid,

The above displacement growth is a resonance effectLandau1 . Recently, we have performed a molecular dynamics simulation of a glass under periodic shearperiodic , where plastic events are proliferated at resonance frequencies. Wittmer et alWit1 studied resonance in a model network, where the growth is suppressed by the viscous friction.

V Summary and Remarks

In Sec.II, we have presented a linear response theory of applying a mean shear strain from boundary walls as in real experiments. Our theory is based on the relations (9) and (10) and is applicable to any solids and fluids in confined geometries. It can describe the linear dynamics in the bulk and near the boundary walls in terms of the appropriate time-correlation functions. It can be used even when the cell interior is inhomogeneous. As a first nontrivial application, this paper has mostly treated the linear response in glasses around a fixed inherent state.

In Sec.III, we have discussed the linear response in supercooled liquids with slow dynamics and in ordinary liquids with fast relaxations. In these states, the viscosity comes into play in the bulk. Our theory can further be used to study the boundary flow effects microscopically.

In Sec.IV, we have presented numerical results in a two-dimensional model glass on the basis of our theory. In Figs.3 and 4, the forces from the walls are correlated with the displacements of all the particles in the cell, resulting in heterogeneous responses. In Figs.5-9, time-dependent responses and time-correlation functions are strongly influenced by sound wave propagation and are very singular in glasses in the film geometry.

We make some remarks as follows.

(i) In real film systems, their temperature is regulated by heat transport between the walls and the interior particles. To realize this situation in molecular dynamics simulation, we can attach heat bathes to the boundary layers (not to the interior particles)Shiba . Sounds are then damped upon reflection at the walls, which leads to decays of the time-correlation functions.

(ii) Plastic events emit sounds leading to fast transport of the released potential energies throughout the system. Shiba and one of the present authorsShiba found that plastic events cause large oscillatory deviations in the force difference $F_{\rm top}^{x}(t)-F_{\rm bot}^{x}(t)$, where the heat bathes in the boundary layers damp such oscillations.

(iii) We can construct a microscopic theory of applying dilational strains by moving the walls along the $z$ axis, where propagation of the longitudinal sounds is crucial. In particular, the correlation between the normal components $\delta F_{\rm top}^{z}$ and $\delta F_{\rm bot}^{z}$ can be obtained if $\mu$ is replaced by $B+(2-2/d)\mu$ in Eq.(50), where $B$ is the bulk modulus.

(iv) In future, we should use a more realistic model of solid walls, where the forces from the particles in the walls to those in the cell are of great importanceslip . As discussed in Sec.IIIB, we can apply our theory to fluids to investigate the boundary flow profiles.

(v) There are a number of elastic systems with inhomogeneous elastic moduli on mesoscopic scalesHead ; Onukibook ; Hecke ; Luben ; Wit1 ; Bas , where nonaffine strains appear in applied stress. In gelsHead ; Onukibook , the crosslink structure is intrinsically random depending on the preparation condition. In multi-component metallic alloysOnukibook , the elasticity in the presence of precipitates is of great technological importance. which are harder or softer than the matrix.

Acknowledgements.

This work was partially supported by the Japan Society for the Promotion of Science Grants-in-Aid for Scientific Research (KAKENHI) (grants No. 16H04025, No. 16H04034, No. 16H06018, and 18H01188).

**Appendix A: One-dimensional elastic systems

**

Here, we examine random elastic systems in one dimension in our theoretical scheme. We apply a dilational strain $\epsilon_{\rm ex}$ from the walls at low $T$ to obtain an analytic expression for the dilational elastic modulus. We can mention a number of random network modelsWit1 ; Hecke ; Luben . Along the $x$ axis, the particles are at $x_{i}$ ($1\leq i\leq N)$ and the walls are at $R_{1}$ and $R_{N}$, where $R_{1}<x_{1}<\cdots<x_{N}<R_{N}$ and $N\gg 1$. We set $R_{N}=H/2$ and $R_{1}=-H/2$ with $H$ being the system length. The end particles 1 and N are bound to the walls by potentials $\psi(d_{0})$ and $\psi(d_{N})$, where $d_{0}=x_{1}-R_{1}$ and $d_{N}=R_{N}-x_{N}$. Particles $i$ and $i+1$ interact with potentials $\phi_{i}(d_{i})$, where $d_{i}=x_{i+1}-x_{i}$. The total potential energy is given by

[TABLE]

The potentials $\phi_{i}(d_{i})$ can be random depending on $i$. Hereafter, we set ${\sum_{i}}^{\prime}=\sum_{1\leq i<N}$. The potential part of the microscopic pressure is of the form,

[TABLE]

where $\phi_{i}^{\prime}(x)=d\phi_{i}(x)/dx$ and

[TABLE]

Here, $\theta(x)$ is the step function equal to 1 for $x>0$ and to 0 for $x\leq 0$. Thus, $\theta_{i}(x)$ is 1 in the interval $x_{i}<x<x_{i+1}$ and is zero outside it. From Eq.(7) we find ${\hat{\delta}}(x,x_{i},x_{i+1})=\theta_{i}(x)/d_{i}$, leading to Eq.(A2).

At fixed $H=R_{N}-R_{1}$, we assume that the mechanical equilibrium holds at $x_{i}=x_{i}^{{\rm ih}}$ as

[TABLE]

where $\psi^{\prime}(x)=d\psi(x)/dx$ and $d_{i}^{\rm ih}$, $d_{0}^{\rm ih}$, and $d_{N}^{\rm ih}(=d_{0}^{\rm ih})$ are the equilibrium values of $d_{i}$, $d_{0}$, and $d_{N}$, respectively. This state corresponds to the inherent state in glasses. From Eq.(A4) the potential part of the pressure is $\Pi^{\rm p}=p_{\rm ih}$ for $x_{1}^{\rm ih}<x<x_{N}^{\rm ih}$. We next consider small displacements $u_{i}=x_{i}-x_{i}^{\rm ih}$. At fixed $H$, the deviation of $U$ is bilinear in $u_{i}$ to leading order as

[TABLE]

where $\xi_{i}=d_{i}-d_{i}^{\rm ih}=u_{i+1}-u_{i}$, ${\sum_{i}}^{\prime}\xi_{i}=u_{N}-u_{1}$, and

[TABLE]

with $\phi_{i}^{\prime\prime}(x)=d^{2}\phi_{i}(x)/dx^{2}$ and $\psi^{\prime\prime}(x)=d^{2}\psi(x)/dx^{2}$. We assume $s_{i}>0$ and $s_{0}>0$. At low $T$, $u_{i}$ fluctuate thermally obeying the Gaussian distribution $P_{\rm ih}\propto\exp(-\delta U/k_{B}T)$. Hereafter, $\langle{\cdots}\rangle_{\rm ih}$ denotes this average.

To linear order, the deviation of $\Pi^{\rm p}(x)$ is written as

[TABLE]

where $\theta_{i}^{\rm ih}(x)=\theta(x-x_{i}^{\rm ih})-\theta(x,x_{i+1}^{\rm ih})$. The last term in Eq.(A7) vanishes in the range $x_{1}^{\rm ih}<x<x_{N}^{\rm ih}$. We then consider the correlation function for the thermal fluctuations of $\delta\Pi^{\rm p}$ defined by

[TABLE]

where $x$ and $x^{\prime}$ are in the range $[x_{1}^{\rm ih},x_{N}^{\rm ih}]$. Since $P_{\rm ih}$ is Gaussian, some calculations readily give

[TABLE]

The first term is nonvanishing only when $x$ and $x^{\prime}$ are in the same interval, so it is short-ranged. The second term is a constant ($\propto H^{-1})$, which arises from the global elastic coupling. As will be shown in Eqs.(A14) and (A16), $K$ has the meaning of the elastic constant given by

[TABLE]

where $\langle{s^{-1}}\rangle_{\rm ih}$ is the average of $s_{i}^{-1}$ over all the bonds.

To derive Eq.(A9) we can use the variance relations,

[TABLE]

where $c_{i}=k_{B}T/s_{i}$ and $e_{i}=K/Hs_{i}$. We obtain the displacements from $u_{i}=u_{0}+\sum_{1\leq k<i}\xi_{k}$. ThusLuben ,

[TABLE]

where $\ell_{ij}=\sum_{i\leq k<j}1/s_{k}$ for $i<j$. The two terms in Eq.(A12) grow with increasing $j-i$ but largely cancel for $\ell_{ij}>H/2K$. For $i=1$ and $j=N$, they nearly cancel as $\langle{(u_{N}-u_{1})^{2}}\rangle_{\rm ih}/k_{B}T=2(1-2e_{0})/s_{0}$. We also find the counterpart of Eq.(50) in the form,

[TABLE]

Next, we shift the top as $R_{N}\to R_{N}+\delta R_{N}$ keeping the bottom at rest with mean strain $\epsilon_{\rm ex}=\delta R_{N}/H$. The potential energy $U$ changes by $-p_{\rm ih}\delta R_{N}+\delta U^{\prime}$ with

[TABLE]

which corresponds to Eq.(46). Minimization of the first line of Eq.(A14) yields shifts of the particle positions given by $\xi_{i}=K{\epsilon_{\rm ex}}/{s_{i}}$, $u_{1}=\delta R_{N}-u_{N}=K{\epsilon_{\rm ex}}/{s_{0}}$, for which we obtain the second line.

We also use the linear response theoryKubo for a small strain. The perturbed Hamiltonian is {\cal H}_{\rm ih}^{\prime}=\sum_{i}|{\mbox{\boldmathp}}_{i}|^{2}/2m_{i}+\delta U-\epsilon_{\rm ex}\delta{\cal A}, where the first term is the kinetic energy, $\delta U$ is given by Eq.(A5), and

[TABLE]

where $f_{i}=-\partial(\delta U)/\partial u_{i}$. As in Eq.(38) we obtain

[TABLE]

From Eq.(A15), $K_{\rm a}$ is the affine part and $K_{\rm c}$ is the nonaffine part. From Eqs.(A7) and (A8) we find

[TABLE]

In $K_{\rm c}$, the first and second terms in Eq.(A9) yield $K_{\rm a}$ and $-K$, respectively, after double integration.

We notice that, if some bonds are very weak (with very small $s_{i}$), they can dominantly contribute to $\langle{s^{-1}}\rangle_{\rm ih}$ giving rise to a large reduction in $K$. In such cases, $K$ can be much smaller than $K_{\rm a}$.

**Appendix B: Thermal fluctuations in solid films in continuum elasticity

**

In the linear elasticity, we consider a solid film with homogeneous elastic moduli, where the displacement field {\mbox{\boldmathu}}({\mbox{\boldmathr}}) is well-defined. We assume that the thermal fluctuations of $u$ obey the distribution $\propto\exp(-F_{\rm el}/k_{B}T)$, where $F_{\rm el}$ is the elastic free energyGusev ; Binder ; Rahman . Here, the average over this distribution is written as $\langle{\cdots}\rangle$. We impose the rigid boundary condition at $z=\pm H/2$ and the periodic boundary condition along the $x$ and $y$ axes with period $L$. The film volume is $V=HS$ with $S=L^{d-1}$. .

First, we examine the equal-time correlations. For simplicity, we consider the lateral average of $u_{x}$ (the zero wavenumber component in the $x$-$y$ plane) given by

[TABLE]

where \int d{\mbox{\boldmathr}}_{\perp} is the integral in the $x$-$y$ plane. The normalized eigenfunctions are $e_{n}(z)=\sqrt{2/V}\sin[k_{n}(z+H/2)]$ $(n\geq 1$) with $k_{n}=\pi n/H$ for $U_{x}(\pm H/2)=0$. As in Eq.(55), we introduce the fluctuating variables $s_{n}$ by

[TABLE]

As in Eq.(57), the elastic free energy is expressed as

[TABLE]

where $\mu$ is the shear modulus and $U_{x}^{\prime}(z)=dU_{x}(z)/dz$. Then, we find $\langle{s_{n}s_{\ell}}\rangle=k_{B}T\delta_{n\ell}/\mu k_{n}^{2}$ and

[TABLE]

where we use the formula $\sum_{n\geq 1}\sin(np)\sin(nq)/n^{2}=q(\pi-p)/2$ for $0<q<p<\pi$. Differentiation of Eq.(B4) with respect to $z$ and $z^{\prime}$ gives the strain correlation,

[TABLE]

Here, the $\delta$-function appears, but it should be regarded as a function with a microscopic width in particle systems.

In Eq.(28) we have introduced the forces from the walls. In the continuum theory, we express them as

[TABLE]

To account for the particle discreteness, we assume the stress balance at $z=\pm H^{\prime}/2$, where $H^{\prime}=H-\ell_{\rm m}$ with $\ell_{\rm m}$ being a microscopic length. For $|z|<H^{\prime}/2$, we find $\langle{U_{x}(z)\delta F_{\rm top}^{x}}\rangle/k_{B}T=-z/H-1/2$ and $\langle{U_{x}(z)\delta F_{\rm bot}^{x}}\rangle/k_{B}T=z/H-1/2$ from Eq.(B4). These lead to the counterparts of Eqs.(35) and (47). If we define $\delta{\cal A}=H(\delta F_{\rm bot}^{x}-\delta F_{\rm top}^{x})/2$ as in Eq.(27), we find

[TABLE]

which consists of the affine displacement only. To be precise, $\langle{U_{x}(z)\delta{\cal A}}\rangle$ nearly vanishes in the narrow layers $H^{\prime}/2<|z|<H/2$. If we assume the linear combination $\delta{\cal A}=\sum_{n}Z_{A}^{(n)}s_{n}$ in terms of $s_{n}$ in Eq.(B2), the coefficients $Z_{A}^{(n)}$ are nonvanishing only for even positive $n$ as

[TABLE]

Next, we examine the time-correlations assuming the wave equation $\partial^{2}U_{x}/\partial t^{2}=c_{\perp}^{2}\partial^{2}U_{x}/\partial z^{2}$ without dissipation, where $c_{\perp}$ is the transverse sound speed. Fixing $z^{\prime}$ in the cell $|z^{\prime}|<H/2$, we have

[TABLE]

where $2g(z,z^{\prime})$ is equal to Eq.(B4) for $|z|<H/2$. We extend it outside the cell setting $g(z,z^{\prime})=-g(H-z,z^{\prime})=g(z+2H,z^{\prime})$. We can obtain Eq.(B9) if we replace $e_{n}(z)e_{n}(z^{\prime})$ by $e_{n}(z)e_{n}(z^{\prime})\cos(c_{\perp}k_{n}t)$ in the first line of Eq.(B4). Then, Eq.(B9) gives a periodic function of $t$ and the period is twice larger than the acoustic traversal time $t_{\rm a}=H/c_{\perp}$. From Eqs.(B6) and (B9), the function ${\cal F}_{\rm w}(t)$ defined in Eq.(82) is calculated as

[TABLE]

where the second term arises from impulses due to repeated reflections of transverse sounds without scattering. Thus, Eq.(B10) is consistent with Fig.8 for $t<t_{\rm a}$.