text:
# Frequency-independent effective rest mass of photons in the 2TDLM model

## I Introduction

The effective mass of electrons and holes in a semiconductor possesses the physical significance: it contains details and information on the interaction between electrons and lattices, and provides us with a vivid and convenient description of motion of electrons and electron-lattice couplings in a solid state. Moreover, this effective mass can be measured experimentally. So, the effective mass of electrons in a semiconductor is regarded as an essential concept, which is helpful for describing the electron wave propagating inside the semiconductor (and the solid state). Likewise, it is reasonably believed that the effective rest mass of photons in media may also be physically meaningful. In the present paper, two kinds of effective rest mass (i.e., the frequency-dependent and frequency-independent rest mass) of photons will be suggested. The definition of frequency-dependent effective rest mass of photons in media will be merely briefly discussed and its deficiencies will also be demonstrated from the physical point of view. To eliminate the drawback of the frequency-dependent effective rest mass, we will suggest a new definition of effective rest mass, i.e., the frequency-independent photon rest mass in electromagnetic media and consider it in more detail. It is shown that the frequency-independent effective rest mass of photons can be extracted from the optical refractive index squared $`n^2\left(\omega \right)`$. As an illustrative example, we obtain the expression for the frequency-independent effective rest mass in the two time derivative Lorentz material (2TDLM) model.

Recently, metamaterials such as left-handed media, has captured attention of a large number of researchers in various fields, e.g., materials science, condensed matter physics, optics and electromagnetism as well. One of the metamaterials models proposed by Ziolkowski and Auzanneau is called the two time derivative Lorentz material (2TDLM) model, which is a generalization of the standard Lorentz material model and encompasses the permittivity and permeability material responses experimentally obtained, which may be seen, for example, in the reference of Smith et al..

This paper is organized as follows: first we consider the frequency-dependent effective rest mass of photons in media. The formulation of the frequency-independent effective rest mass is presented in Sec. III. In Sec. IV, the frequency-independent effective rest mass in the 2TDLM model and left-handed media is read off from $`n^2\left(\omega \right)`$. In Sec. V, we conclude with some remarks on the frequency-independent effective rest mass of photons in media.

## II Frequency-dependent effective rest mass of photons in media

We have studied the wave propagation problem (i.e., polarization of plane wave, geometric phases, time evolution of photon wavefunction, helicity reversal) of photons field in media (optical fiber) at the quantum level, where we have not considered the dispersion and absorption properties of media. In the present paper, these electromagnetic and optical properties of media contributing to the wave propagation of light field will be taken into account. It is readily verified that the interaction of light wave with media can give rise to an effective rest mass of photons, which may be illustrated as follows:

if the concept of effective rest mass is applicable to the photons in media, then in terms of the Einstein-de Broglie relation, a photon possessing $`m_{\mathrm{eff}}`$ agrees with the dispersion relation $`\left(\mathrm{}k\right)^2c^2=\left(\mathrm{}\omega \right)^2m_{\mathrm{eff}}^2c^4`$ with $`\mathrm{}`$ and $`c`$ being the Planck’s constant and the speed of light in a vacuum, respectively. Thus one can arrive at $`k^2c^2=\omega ^2\left(1\frac{c^4}{\mathrm{}^2}\frac{m_{\mathrm{eff}}^2}{\omega ^2}\right)`$ or

$$\frac{\omega ^2}{k^2}=\frac{c^2}{1\frac{c^4}{\mathrm{}^2}\frac{m_{\mathrm{eff}}^2}{\omega ^2}}.$$
(1)
Since the relation between the frequency $`\omega `$ and wave vector $`k`$ of light wave propagating inside a dispersive medium may be $`\frac{\omega ^2}{k^2}=\frac{c^2}{n^2}`$ with $`n`$ being the optical refractive index. So, with the help of (1) and dispersion relation $`\frac{\omega ^2}{k^2}=\frac{c^2}{n^2}`$, one can conclude that the effective rest mass $`m_{\mathrm{eff}}`$, refractive index $`n`$ and wave frequency $`\omega `$ together satisfy the following relation

$$\frac{m_{\mathrm{eff}}^2c^4}{\mathrm{}^2}=(1n^2)\omega ^2.$$
(2)
Thus we obtained the expression for the effective rest mass squared $`m_{\mathrm{eff}}^2`$. For Lorentz dispersive media, the refractive index squared is of the form

$$n^2=ϵ=1+\frac{Ne^2}{ϵ_0m_\mathrm{e}}\frac{1}{\omega _0^2\omega ^2i\gamma \omega },$$
(3)
and the effective rest mass squared is therefore written

$$m_{\mathrm{eff}}^2=\frac{\mathrm{}^2}{c^2}\left(\frac{\mu _0Ne^2}{m_\mathrm{e}}\frac{\omega ^2}{\omega ^2+i\gamma \omega \omega _0^2}\right),$$
(4)
where use is made of the relation $`c^2=\frac{1}{ϵ_0\mu _0}`$. Here $`\gamma `$ and $`\omega _0`$ respectively denote the damping frequency and resonance frequency of the electric dipole oscillators; $`N`$ and $`m_\mathrm{e}`$ represent the electron number per unit volume and electron mass, respectively. Note that this effective rest mass of photons in electromagnetic media originates from the interaction term $`\mu _0𝐉𝐀`$ in the Lagrangian density, where $`𝐉`$ and $`𝐀`$ respectively denote the electric current density and the magnetic vector potentials. Since the electric current acted upon by external electromagnetic wave behaves like the self-induced charge current, $`𝐉`$ is proportional to $`𝐀`$, for instance, in Lorentz materials, the electric current density in the presence of electromagnetic wave is given

$$𝐉=Ne𝐯=\frac{Ne^2\omega ^2}{m_\mathrm{e}}\frac{1}{\omega _0^2\omega ^2i\gamma \omega }𝐀.$$
(5)
Thus the interaction term $`\mu _0𝐉𝐀`$ is changed into a mass term $`\frac{m_{\mathrm{eff}}^2c^2}{\mathrm{}^2}𝐀^2`$ of electromagnetic fields in the Lagrangian density ( see in the Appendix for the detailed discussions on this problem ). Apparently, if $`\gamma =0`$ and $`\omega _0^2=0`$, then the expression (4) for $`m_{\mathrm{eff}}^2`$ is reduced to the rest mass squared that is independent of $`\omega `$. We are certain that in this case ($`\gamma =\omega _0=0`$), $`m_{\mathrm{eff}}`$ in (4) possesses a form of plasma-type effective rest mass $`m_{\mathrm{eff}}=\frac{\mathrm{}\omega _\mathrm{p}}{c^2}`$ with $`\omega _\mathrm{p}^2=\frac{\mu _0c^2Ne^2}{m_\mathrm{e}}`$. It is readily verified that $`L=\frac{\mathrm{}}{m_{\mathrm{eff}}c}`$ is just the Londons’ penetration length (i.e., the magnetic field can penetrate only the scale $`L`$ within the surface of the superconductor). Thus it is concluded without fear that the above formulation for effective rest mass of photons in media is self-consistent.

However, it should be noted that in Lorentz dispersive materials, $`m_{\mathrm{eff}}^2`$ in (4) depends strongly on $`\omega `$, and sometimes possesses an imaginary part (due to the damping parameter and resonance frequency) and the physical meanings of $`m_{\mathrm{eff}}^2`$ is therefore not very apparently seen. In what follows we will define a physically meaningful effective rest mass of photons, which is independent of the frequency $`\omega `$, and will explain a method by which the frequency-independent effective rest mass of photons in the 2TDLM media is extracted from $`n^2\left(\omega \right)`$. Note that the above non-physical disadvantage of frequency-dependent mass in (4) will be avoided in the frequency-independent effective rest mass.

## III Extracting frequency-independent effective mass from refractive index squared

In this framework of definition of frequency-independent effective rest mass, the photon should be regarded as acted upon by a hypothetical force field, namely, it behaves like a massive de Broglie particle moving in a force field. It is verified in the following that, for the case of de Broglie particle in a force field, one can also consider its “optical refractive index” $`n`$. According to the Einstein- de Broglie relation, the dispersion relation of the de Broglie particle with rest mass $`m_0`$ in the potential field $`V\left(𝐱\right)`$ agrees with

$$\left(\omega \varphi \right)^2=k^2c^2+\frac{m_0^2c^4}{\mathrm{}^2}$$
(6)
with $`\varphi =\frac{V}{\mathrm{}}`$. It follows that

$$\omega ^2\left[1\frac{m_0^2c^4}{\mathrm{}^2\omega ^2}2\frac{\varphi }{\omega }+\left(\frac{\varphi }{\omega }\right)^2\right]=k^2c^2,$$
(7)
which yields

$$\frac{\omega ^2}{k^2}=\frac{c^2}{1\frac{m_0^2c^4}{\mathrm{}^2\omega ^2}2\frac{\varphi }{\omega }+\left(\frac{\varphi }{\omega }\right)^2}.$$
(8)
Compared with the dispersion relation $`\frac{\omega ^2}{k^2}=\frac{c^2}{n^2}`$, one can arrive at

$$n^2=1\frac{m_0^2c^4}{\mathrm{}^2\omega ^2}2\frac{\varphi }{\omega }+\left(\frac{\varphi }{\omega }\right)^2,$$
(9)
which is the square of “optical refractive index” of the de Broglie particle in the presence of a potential field $`\varphi `$.

If a photon is permitted to possess a frequency-independent effective rest mass $`m_{\mathrm{eff}}`$, then the expression (9) for the case of photons in media is rewritten

$$n^2\left(\omega \right)=1\frac{m_{\mathrm{eff}}^2c^4}{\mathrm{}^2\omega ^2}2\frac{\varphi }{\omega }+\left(\frac{\varphi }{\omega }\right)^2.$$
(10)
In what follows, an approach to frequency-independent effective mass extracted from the optical refractive index squared $`n^2\left(\omega \right)`$ is presented. It is assumed that $`n^2\left(\omega \right)`$ can be rewritten as the following series expansions

$$n^2\left(\omega \right)=\underset{k}{}\frac{a_k}{\omega ^k},$$
(11)
then one can arrive at $`n^2\left(\omega \right)\omega ^2=_k\frac{a_k}{\omega ^{k2}}`$, namely,

$$n^2\left(\omega \right)\omega ^2=\mathrm{}+\frac{a_3}{\omega }+a_2+a_1\omega +a_0\omega ^2+a_{+1}\omega ^3+\mathrm{}.$$
(12)
Compared (12) with (10), it is apparent that the frequency-independent effective rest mass squared $`m_{\mathrm{eff}}^2`$ of the photon can be read off from $`n^2\left(\omega \right)`$, namely, the constant term on the right-handed side on (12) is related only to $`m_{\mathrm{eff}}^2`$, i.e.,

$$a_2=\frac{m_{\mathrm{eff}}^2c^4}{\mathrm{}^2}.$$
(13)
Generally speaking, the optical refractive index $`n^2\left(\omega \right)`$ is often of complicated form, particularly for the metamaterials. So, we cannot extract frequency-independent $`m_{\mathrm{eff}}^2`$ immediately from $`n^2\left(\omega \right)`$. In the following, this problem is resolved by calculating the limit value of $`n^2\left(\omega \right)\omega ^2`$, where $`\omega `$ tends to $`\mathrm{}`$ or zero. We first consider the large frequency behavior of $`n^2\left(\omega \right)\omega ^2`$. It follows from (12) that $`n^2\left(\omega \right)\omega ^2`$ at high frequencies (i.e., $`\omega \mathrm{}`$) is of the form

$$n^2\left(\omega \right)\omega ^2a_2+a_1\omega +a_0\omega ^2+a_{+1}\omega ^3+a_{+2}\omega ^4+\mathrm{}.$$
(14)
This, therefore, means that $`n^2\left(\omega \right)\omega ^2`$ at high frequencies ($`\omega \mathrm{}`$) can be split into the following two terms

$$n^2\left(\omega \right)\omega ^2\left(\mathrm{constant}\mathrm{term}\right)+\left(\mathrm{divergent}\mathrm{term}\right).$$
(15)
So, we can obtain $`m_{\mathrm{eff}}^2`$ from the following formula

$$\underset{\omega \mathrm{}}{lim}\left[n^2\left(\omega \right)\omega ^2\left(\mathrm{divergent}\mathrm{term}\right)\right]=\frac{m_{\mathrm{eff}}^2c^4}{\mathrm{}^2}.$$
(16)
Note that the constant term in $`n^2\left(\omega \right)\omega ^2`$ can also be extracted by considering the low-frequency (zero-frequency) behavior of $`n^2\left(\omega \right)`$. It is certain that in this case the obtained constant term is the same as (16). Apparently $`n^2\left(\omega \right)\omega ^2`$ at low-frequencies (i.e., $`\omega 0`$) behaves like

$$n^2\left(\omega \right)\omega ^2\mathrm{}+\frac{a_5}{\omega ^3}+\frac{a_4}{\omega ^2}+\frac{a_3}{\omega }+a_2,$$
(17)
where the sum term ($`\mathrm{}+\frac{a_5}{\omega ^3}+\frac{a_4}{\omega ^2}+\frac{a_3}{\omega }`$) is divergent if $`\omega `$ approaches zero. Thus, subtraction of the divergent term from $`n^2\left(\omega \right)\omega ^2`$ yields

$$\underset{\omega 0}{lim}\left[n^2\left(\omega \right)\omega ^2\left(\mathrm{divergent}\mathrm{term}\right)\right]=\frac{m_{\mathrm{eff}}^2c^4}{\mathrm{}^2}.$$
(18)
In the next section, we will calculate the frequency-independent effective rest mass of photons in 2TDLM model by making use of the method presented above.

## IV Frequency-independent effective rest mass of photons in 2TDLM model

According to the definition of the 2TDLM model, the frequency-domain electric and magnetic susceptibilities are given

$$\chi ^e\left(\omega \right)=\frac{\left(\omega _\mathrm{p}^e\right)^2\chi _\alpha ^e+i\omega \omega _\mathrm{p}^e\chi _\beta ^e\omega ^2\chi _\gamma ^e}{\omega ^2+i\omega \mathrm{\Gamma }^e+\left(\omega _0^e\right)^2},\chi ^m\left(\omega \right)=\frac{\left(\omega _\mathrm{p}^m\right)^2\chi _\alpha ^m+i\omega \omega _\mathrm{p}^m\chi _\beta ^m\omega ^2\chi _\gamma ^m}{\omega ^2+i\omega \mathrm{\Gamma }^m+\left(\omega _0^m\right)^2},$$
(19)
where $`\chi _\alpha ^{e,m}`$, $`\chi _\beta ^{e,m}`$ and $`\chi _\gamma ^{e,m}`$ represent, respectively, the coupling of the electric (magnetic) field and its first and second time derivatives to the local electric (magnetic) dipole motions. $`\omega _\mathrm{p}`$, $`\mathrm{\Gamma }^e`$, $`\mathrm{\Gamma }^m`$ and $`\omega _0`$ can be viewed as the plasma frequency, damping frequency and resonance frequency of the electric (magnetic) dipole oscillators, respectively. So, the refractive index squared in the 2TDLM model reads

$$n^2\left(\omega \right)=1+\chi ^e\left(\omega \right)+\chi ^m\left(\omega \right)+\chi ^e\left(\omega \right)\chi ^m\left(\omega \right).$$
(20)
In order to obtain $`m_{\mathrm{eff}}^2`$ from $`n^2\left(\omega \right)`$, according to the formulation in the previous section, the high-frequency ($`\omega \mathrm{}`$) behavior of $`\chi ^e\left(\omega \right)\omega ^2`$ should first be taken into consideration and the result is

$$\chi ^e\left(\omega \right)\omega ^2i\omega \omega _\mathrm{p}^e\chi _\beta ^e+\omega ^2\chi _\gamma ^e.$$
(21)
Note that $`i\omega \omega _\mathrm{p}^e\chi _\beta ^e+\omega ^2\chi _\gamma ^e`$ is the divergent term of $`\chi ^e\left(\omega \right)\omega ^2`$ at large frequencies. Subtracting the divergent term from $`\chi ^e\left(\omega \right)\omega ^2`$, we consider again the high-frequency behavior ($`\omega \mathrm{}`$) and obtain another divergent term $`i\omega ^3\mathrm{\Gamma }^e\chi _\gamma ^e`$, i.e.,

$$\chi ^e\left(\omega \right)\omega ^2\left(i\omega \omega _\mathrm{p}^e\chi _\beta ^e+\omega ^2\chi _\gamma ^e\right)=\frac{\omega ^2\left(\omega _\mathrm{p}^e\right)^2\chi _\alpha ^e\omega ^2\omega _\mathrm{p}^e\mathrm{\Gamma }^e\chi _\beta ^e+i\omega \omega _\mathrm{p}^e\left(\omega _0^e\right)^2\chi _\beta ^ei\omega ^3\mathrm{\Gamma }^e\chi _\gamma ^e\omega ^2\left(\omega _0^e\right)^2\chi _\gamma ^e}{\omega ^2+i\omega \mathrm{\Gamma }^e+\left(\omega _0^e\right)^2}i\omega ^3\mathrm{\Gamma }^e\chi _\gamma ^e.$$
(22)
Thus, after subtracting all the divergent terms from $`\chi ^e\left(\omega \right)\omega ^2`$, we obtain

$$\underset{\omega \mathrm{}}{lim}\left[\chi ^e\left(\omega \right)\omega ^2\left(i\omega \omega _\mathrm{p}^e\chi _\beta ^e+\omega ^2\chi _\gamma ^e\right)i\omega ^3\mathrm{\Gamma }^e\chi _\gamma ^e\right]=\left[\left(\omega _\mathrm{p}^e\right)^2\chi _\alpha ^e\omega _\mathrm{p}^e\mathrm{\Gamma }^e\chi _\beta ^e\left(\omega _0^e\right)^2\chi _\gamma ^e+\left(\mathrm{\Gamma }^e\right)^2\chi _\gamma ^e\right].$$
(23)
In the same fashion, we obtain

$$\underset{\omega \mathrm{}}{lim}\left[\chi ^m\left(\omega \right)\omega ^2\left(i\omega \omega _\mathrm{p}^m\chi _\beta ^m+\omega ^2\chi _\gamma ^m\right)i\omega ^3\mathrm{\Gamma }^m\chi _\gamma ^m\right]=\left[\left(\omega _\mathrm{p}^m\right)^2\chi _\alpha ^m\omega _\mathrm{p}^m\mathrm{\Gamma }^m\chi _\beta ^m\left(\omega _0^m\right)^2\chi _\gamma ^m+\left(\mathrm{\Gamma }^m\right)^2\chi _\gamma ^m\right]$$
(24)
for the magnetic susceptibility $`\chi ^m\left(\omega \right)`$. In what follows we continue to consider the term $`\chi ^e\left(\omega \right)\chi ^m\left(\omega \right)`$ on the right-handed side of (20). In the similar manner, the infinite-frequency limit of $`\chi ^e\left(\omega \right)\chi ^m\left(\omega \right)\omega ^2`$ is given as follows

$`\underset{\omega \mathrm{}}{lim}\chi ^e\left(\omega \right)\chi ^m\left(\omega \right)\omega ^2`$ $`=`$ $`\chi _\gamma ^m\left[\left(\omega _\mathrm{p}^e\right)^2\chi _\alpha ^e\omega _\mathrm{p}^e\mathrm{\Gamma }^e\chi _\beta ^e\left(\omega _0^e\right)^2\chi _\gamma ^e+\left(\mathrm{\Gamma }^e\right)^2\chi _\gamma ^e\right]`$ (25)
$``$ $`\chi _\gamma ^e\left[\left(\omega _\mathrm{p}^m\right)^2\chi _\alpha ^m\omega _\mathrm{p}^m\mathrm{\Gamma }^m\chi _\beta ^m\left(\omega _0^m\right)^2\chi _\gamma ^m+\left(\mathrm{\Gamma }^m\right)^2\chi _\gamma ^m\right].`$ (26)
Hence, insertion of (23)-(26) into (16) yields

$`{\displaystyle \frac{m_{\mathrm{eff}}^2c^4}{\mathrm{}^2}}`$ $`=`$ $`\left(1+\chi _\gamma ^m\right)\left[\left(\omega _\mathrm{p}^e\right)^2\chi _\alpha ^e\omega _\mathrm{p}^e\mathrm{\Gamma }^e\chi _\beta ^e\left(\omega _0^e\right)^2\chi _\gamma ^e+\left(\mathrm{\Gamma }^e\right)^2\chi _\gamma ^e\right]`$ (27)
$`+`$ $`\left(1+\chi _\gamma ^e\right)\left[\left(\omega _\mathrm{p}^m\right)^2\chi _\alpha ^m\omega _\mathrm{p}^m\mathrm{\Gamma }^m\chi _\beta ^m\left(\omega _0^m\right)^2\chi _\gamma ^m+\left(\mathrm{\Gamma }^m\right)^2\chi _\gamma ^m\right].`$ (28)
Thus the $`\omega `$-independent effective rest mass squared $`m_{\mathrm{eff}}^2`$ is extracted from $`n^2\left(\omega \right)`$. Apparently, $`m_{\mathrm{eff}}^2`$ in Eq. (28) does not depend on $`\omega `$. The remainder on the right-handed side of Eq. (10) is $`n^2\left(\omega \right)1+\frac{m_{\mathrm{eff}}^2c^4}{\mathrm{}^2\omega ^2}`$, which is equal to $`\left(\frac{\varphi }{\omega }\right)^22\frac{\varphi }{\omega }`$. For simplicity, we set

$$\mathrm{\Delta }\left(\omega \right)=n^2\left(\omega \right)1+\frac{m_{\mathrm{eff}}^2c^4}{\mathrm{}^2\omega ^2},$$
(29)
and then from the equation

$$\left(\frac{\varphi }{\omega }\right)^22\frac{\varphi }{\omega }=\mathrm{\Delta }\left(\omega \right),$$
(30)
one can readily obtain the expression for the hypothetical potential field $`\varphi \left(\omega \right)`$ is $`\varphi \left(\omega \right)=\omega \left[1\sqrt{1+\mathrm{\Delta }\left(\omega \right)}\right]`$, or

$$V\left(\omega \right)=\mathrm{}\omega \left[1\sqrt{1+\mathrm{\Delta }\left(\omega \right)}\right],$$
(31)
where $`\varphi =\frac{V}{\mathrm{}}`$ has been previously defined in Sec. III.

Many metamaterials which have complicated expressions for electric and magnetic susceptibilities such as (19) are generally rarely seen in nature and often arises from the artificial manufactures. More recently, one of the artificial composite metamaterials, the left-handed medium which has a frequency band (GHz) where the electric permittivity ($`ϵ`$) and the magnetic permeability ($`\mu `$) are simultaneously negative, has focused attention of many authors both experimentally and theoretically. In the left-handed medium, most phenomena as the Doppler effect, Vavilov-Cherenkov radiation and even Snell’s law are inverted. In 1964Note that, in the literature, many authors mention the year when Veselago suggested the left-handed media by mistake. They claim that Veselago proposed the concept of left-handed media in 1968. On the contrary, the true fact is as follows: Veselago’s excellent paper was first published in July, 1964 \[Usp. Fiz. Nauk 92 517-526\]. In 1968, this original paper was translated into English by W. H. Furry and published again in the journal of Sov. Phys. Usp.., Veselago first considered many peculiar optical and electromagnetic properties, phenomena and effects in this medium and referred to such materials as left-handed media, since in this case the propagation vector $`𝐤`$, electric field $`𝐄`$ and magnetic field $`𝐇`$ of light wave propagating inside it form a left-handed system. It follows from Maxwellian curl equations that such media having negative simultaneously negative $`ϵ`$ and $`\mu `$ exhibit a negative index of refraction, i.e., $`n=\sqrt{ϵ\mu }`$. Since negative refractive index occurs only rarely, this medium attracts attention of many physicists in various fields. In experiments, the negative $`ϵ`$ and $`\mu `$ can be respectively realized by using a network (array) of thin (long) metal wires and a periodic arrangement of split ring resonators. A combination of the two structures yields a left-handed medium. In general, the dielectric parameter $`ϵ\left(\omega \right)`$ and magnetic permeability $`\mu \left(\omega \right)`$ in the isotropic homogeneous left-handed medium are of the form

$$ϵ\left(\omega \right)=1\frac{\omega _\mathrm{p}^2}{\omega \left(\omega +i\gamma \right)},\mu \left(\omega \right)=1\frac{F\omega ^2}{\omega ^2\omega _0^2+i\omega \mathrm{\Gamma }}$$
(32)
with the plasma frequency $`\omega _\mathrm{p}`$ and the magnetic resonance frequency $`\omega _0`$ being in the GHz region. $`\gamma `$ and $`\mathrm{\Gamma }`$ stand for the electric and magnetic damping parameters, respectively. Compared (32) with the electric and magnetic susceptibilities (19), one can arrive at

$`\chi _\alpha ^e=1,\omega _\mathrm{p}^e=\omega _\mathrm{p},\chi _\beta ^e=\chi _\gamma ^e=0,\omega _0^e=0,\mathrm{\Gamma }^e=\gamma ,`$ (33)
$`\chi _\gamma ^m=F,\chi _\alpha ^m=\chi _\beta ^m=0,\omega _0^m=\omega _0,\mathrm{\Gamma }^m=\mathrm{\Gamma }.`$ (34)
So, it is readily verified with the help of (28) that the $`\omega `$-independent effective rest mass of photons inside left-handed media is written as follows

$$\frac{m_{\mathrm{eff}}^2c^4}{\mathrm{}^2}=\left(1F\right)\omega _\mathrm{p}^2+F\left(\omega _0^2\mathrm{\Gamma }^2\right).$$
(35)
It follows that Eq. (35) is a restriction on the electromagnetic parameters $`F,\omega _\mathrm{p},\omega _0,\mathrm{\Gamma }`$ and $`F`$ in the electric permittivity and magnetic permeability (32), since $`m_{\mathrm{eff}}^20`$. Insertion of the experimentally chosen values of the electromagnetic parameters in the literature into Eq. (35) shows that this restriction condition is satisfied. For instance, in Ruppin’s work, $`F,\omega _\mathrm{p},\omega _0,\gamma ,\mathrm{\Gamma }`$ and $`F`$ were chosen $`F=0.56,\omega _\mathrm{p}=10.0\mathrm{GHz},\omega _0=4.0\mathrm{GHz},\gamma =0.03\omega _\mathrm{p},\mathrm{\Gamma }=0.03\omega _0`$, which indicates that the inequality $`m_{\mathrm{eff}}^20`$ is satisfied.

## V Concluding remarks

This paper deals mainly with a problem of how to extract the frequency-independent effective rest mass of photons in media from the optical refractive index squared $`n^2\left(\omega \right)`$. It is of physical interest that the expressions (28) and (35) for $`m_{\mathrm{eff}}^2`$ are the mathematical restrictions on the electromagnetic parameters in the electric permittivity and magnetic permeability of artificial composite metamaterials. Apparently, it is certain that the two essential requirements (28) and (35) are related close to the Kramers-Kronig dispersion relation.

To close this paper, we consider briefly the potential physical significance of $`m_{\mathrm{eff}}^2`$ in electromagnetic media. Since it is independent of frequency of light wave, this effective rest mass of photons possesses the explicit physical meanings. It follows from (28) that all the parameters in $`ϵ\left(\omega \right)`$ and $`\mu \left(\omega \right)`$ contribute to $`m_{\mathrm{eff}}^2`$ and that $`m_{\mathrm{eff}}^2`$ therefore consists of various interactions between light fields and electromagnetic materials. If $`m_{\mathrm{eff}}^2`$ can be obtained by experimental measurements, then we are able to investigate the parameters in $`ϵ\left(\omega \right)`$ and $`\mu \left(\omega \right)`$ as well as the electromagnetic interaction inside media in more detail. It should be noted that, for the plasma-type media such as the superconductor and the array of long metal wires (which possesses a negative permittivity in the GHz range), the above two effective rest mass (i.e., frequency-dependent and frequency-independent rest mass) are equivalent, namely, the frequency-dependent effective rest mass is reduced to the frequency-independent effective rest mass. In the former case of effective rest mass (frequency-dependent), the photon in media is regarded as a free massive de Broglie particle. But, its effective rest mass lacks explicit physical meanings since it is dependent strongly on the frequency $`\omega `$. Whereas in the latter case of effective rest mass (frequency-independent), the photon in media may be considered a massive de Broglie particle acted upon by a hypothetical potential field. Hence, the frequency-independent effective rest mass has somewhat explicit physical meanings and may therefore have potential applications to the investigation of light propagation in electromagnetic media whose permittivity and permeability are of the complicated form. We hold that both theoretical and experimental analysis of this frequency-independent effective rest mass of photons in metamaterials deserve further considerations.
Acknowledgements I thank Sai-Ling He for useful discussions on the wave propagation in left-handed media. This work is supported in part by the National Natural Science Foundation of China under the project No. $`90101024`$.
Appendix

In this Appendix, we introduce the effective rest mass of electromagnetic wave from the fundamental principle ( based on the Lagrangian density of electrodynamics ). Note that in Londons’ electromagnetics for superconductivity, Ginzberg-Landau superconductivity theory and Higgs mechanism, there are also the viewpoints of effective rest mass. Here we consider the effective rest mass of electromagnetic wave in superconducting media ( or electron plasma ) and dispersive/absorptive materials. It is well known that the Lagrangian density of electromagnetics is

$$\mathrm{}=\frac{1}{4}F_{\mu \nu }F^{\mu \nu }\frac{1}{ϵ_0c^2}\rho \phi +\mu _0𝐉𝐀.$$
(36)
The equation of motion of an electron in superconducting media ( or electron plasma ) acted upon by the external electromagnetic wave is $`m_\mathrm{e}\dot{𝐯}=e𝐄`$ with $`𝐄=\frac{}{t}𝐀`$, where dot on the velocity $`𝐯`$ denotes the derivative of $`𝐯`$ with respect to time $`t`$ and $`𝐀`$ stands for the magnetic vector potentials of the applied electromagnetic wave. One can therefore arrive at

$$\frac{\mathrm{d}}{\mathrm{d}t}\left(m_\mathrm{e}𝐯+e𝐀\right)=0,$$
(36)
namely, the canonical momentum $`m_\mathrm{e}𝐯+e𝐀`$ of the electron is conserved. Set $`m_\mathrm{e}𝐯+e𝐀=0`$ and then $`m_\mathrm{e}𝐯=e𝐀`$. So, the electric current density is of the form

$$𝐉=Ne𝐯=\frac{Ne^2}{m_\mathrm{e}}𝐀$$
(37)
with $`N`$ being the electron number in per unit volume. It follows from (36) that the interaction term $`\mu _0𝐉𝐀`$ in Lagrangian density is rewritten as follows

$$\mu _0𝐉𝐀=\frac{\mu _0Ne^2}{m_\mathrm{e}}𝐀^2=\frac{Ne^2}{ϵ_0m_\mathrm{e}c^2}𝐀^2,$$
(38)
where use is made of the relation $`c^2=\frac{1}{ϵ_0\mu _0}`$. It is apparently seen that since electron current acts as the self-induced charge current, the interaction term $`\mu _0𝐉𝐀`$ is therefore transformed into the mass term $`\frac{m_{\mathrm{eff}}^2c^2}{\mathrm{}^2}𝐀^2`$ of electromagnetic fields. Hence, we have

$$\frac{m_{\mathrm{eff}}^2c^2}{\mathrm{}^2}=\frac{Ne^2}{ϵ_0m_\mathrm{e}c^2}$$
(39)
and the effective mass squared

$$m_{\mathrm{eff}}^2=\frac{\mathrm{}^2}{c^2}\left(\frac{Ne^2}{ϵ_0m_\mathrm{e}c^2}\right).$$
(40)
It is readily verified that $`L=\frac{\mathrm{}}{m_{\mathrm{eff}}c}`$ is just the Londons’ penetration length ( i. e., the magnetic field can penetrate only the scale $`L`$ within the surface of the superconductor ).

In what follows we take into account the case of dispersive and absorptive materials. Here we only consider the motion of electrons, and the result presented in the following holds also for ions in media. The equation of motion of electrons acted upon by the time harmonic electromagnetic wave in these electromagnetic materials is written

$$\ddot{𝐱}+\gamma \dot{𝐱}+\omega _0^2𝐱=\frac{e}{m_\mathrm{e}}𝐄_0\mathrm{exp}\left[i\omega t\right].$$
(41)
Here $`𝐱`$, $`\gamma `$ and $`\omega _0`$ denote the electron spatial displacement vector, damping frequency and resonance frequency of the electric dipole oscillators, respectively. One can obtains

$$𝐱=\frac{e}{m_\mathrm{e}}\frac{1}{\omega _0^2\omega ^2i\gamma \omega }𝐄_0\mathrm{exp}\left[i\omega t\right]$$
(42)
from Eq. (41). Further calculation yields

$$\ddot{𝐱}=\frac{e\omega ^2}{m_\mathrm{e}}\frac{1}{\omega _0^2\omega ^2i\gamma \omega }𝐄_0\mathrm{exp}\left[i\omega t\right]=\frac{e\omega ^2}{m_\mathrm{e}}\frac{1}{\omega _0^2\omega ^2i\gamma \omega }\left(\frac{}{t}𝐀\right),$$
(43)
where use is made of $`𝐄_0\mathrm{exp}\left[i\omega t\right]=\frac{}{t}𝐀`$. Both the electron velocity ( $`𝐯=\dot{𝐱}`$ ) and the electric current density J are therefore obtained, i. e.,

$$𝐯=\frac{e\omega ^2}{m_\mathrm{e}}\frac{1}{\omega _0^2\omega ^2i\gamma \omega }𝐀$$
(44)
and

$$𝐉=Ne𝐯=\frac{Ne^2\omega ^2}{m_\mathrm{e}}\frac{1}{\omega _0^2\omega ^2i\gamma \omega }𝐀.$$
(45)
Hence, similar to the above procedure (38)-(40) ( i. e., the effective mass term arises from the interaction term ), here the effective mass can be obtained from

$$\frac{m_{\mathrm{eff}}^2c^2}{\mathrm{}^2}=\frac{\mu _0Ne^2\omega ^2}{m_\mathrm{e}}\frac{1}{\omega _0^2\omega ^2i\gamma \omega }$$
(46)
or

$$m_{\mathrm{eff}}^2=\frac{\mathrm{}^2}{c^2}\left(\frac{\mu _0Ne^2}{m_\mathrm{e}}\frac{\omega ^2}{\omega ^2+i\gamma \omega \omega _0^2}\right).$$
(47)
Apparently, if $`\gamma =0`$ and $`\omega _0^2=0`$, then the expression (47) for $`m_{\mathrm{eff}}^2`$ is reduced to (40), which is independent of $`\omega `$. It is certain that in this case ( $`\gamma =\omega _0=0`$ ), $`m_{\mathrm{eff}}`$ in (47) possesses the form of plasma-type effective mass ( frequency-independent rest mass ).


format: text

text:
# Long-range Static Directional Stress Transfer in a Cracked, Nonlinear Elastic Crust

## 1 Introduction

There are many evidences that faults and earthquakes interact, as suggested by calculations of stress redistribution , elastodynamic propagation of ruptures using laboratory-based friction law , simplified models of multiple faults , as well as general constraints of kinematic and geometric compatibility of the deformations . Maybe the simplest mechanism for earthquake interaction involves stress re-distribution, both static  and dynamical  associated with a given earthquake modeled as a set of dislocations or cracks. In this simple mechanical view, earthquakes cast stress shadows in lobes of stress unloading  and increase the probability of rupture in zones of stress increase , according to the laws of linear elasticity. These elastic stress transfer models are useful for their conceptual simplicity and are increasingly used. Notwithstanding their extended use, the calculations of stress transfer have large uncertainties stemming from (i) the usually poorly known geometry of the rupture surfaces, (ii) the unconstrained homogeneity and amplitude of the stress drop and/or of the slip distribution on the fault plane, (iii) the use of simplified models of the crust (3D semi-infinite, or thin elastic plate, or plate coupled to a semi-infinite visco-elastic asthenosphere, etc.), and (iv) the unknown direction and amplitude of the absolute stress field that pre-existed before the event, including its possible spatial inhomogeneity.

Such elastic stress transfer models seem unable to account for a growing phenomenology of long-range earthquake interactions. For instance, many large earthquakes have been preceded by an increase in the number of intermediate sized events over very broad areas . The relation between these intermediate sized events and the subsequent main event has only recently been recognized on a large scale because the precursory events occur over such a large area that they do not fit prior definitions of foreshocks . In particular, the $`11`$ earthquakes in California with magnitudes greater than $`6.8`$ in the last century are associated with an increase of precursory intermediate magnitude earthquakes measured in a running time window of five years . What is strange about the result is that the precursory pattern occured with distances of the order of $`300`$ to $`500km`$ from the futur epicenter, i.e. at distances up to ten times larger that the size of the futur earthquake rupture. Furthermore, the increased intermediate magnitude activity switched off rapidly after a big earthquake in about half of the cases. This implies that stress changes due to an earthquake of rupture dimension as small as $`35km`$ can influence the stress distribution to distances more than ten times its size. These observations of earthquake-earthquake interactions over long times and large spatial separations have been strengthened by several other works on different catalogs using a variety of techniques . These results defy usual mechanical models of linear elasticity and one proposed explanation is that seismic cycles represent the approach to and retreat from a critical state of a fault network . Within the critical earthquake concept, the anomalous long-range interactions between earthquakes reflect the increasing stress-stress correlation length upon the approach of the critical earthquake. Another explanation involves dynamical stress triggering  (see however ). Additional seismic, geophysical, and hydrogeological observations  cannot be accounted for by using models derived from the elastic stress transfer mechanism. In particular, standard poro-elastic models underestimate grossly the observed amplitudes of hydrogeological anomalous rises and drops in wells at large distances from earthquakes.

## 2 Mechanical Model of the Earth’s Crust

Here, we investigate the hypothesis, and its implications for the above observations, that the crust is a nonlinear elastic medium characterized by an asymmetric response to compressive versus extensive perturbations around the lithostatic stress. We call this a “threshold nonlinearity.” This nonlinearity stems from a mechanically-justified argumentation based on the fact that the Earth’s crust at seismogenic depth is crisscrossed by joints, cracks or faults at many different scales filled with drained fluid in contact with delocalized reservoirs at pressures close to the lithostatic pressure. It has been argued that rock permeability and thus microcracking adjusts itself, so that fluid pressure is always close to rock pressure irrespective of the extend of hydration/dehydration . One possible mechanism for this involves a time-dependent process that relates fluid pressure, flow pathways and fluid volumes .

### 2.1 Presence and role of fluids

Indeed, a lot of data collectively support the existence of significant fluid circulation to crustal depths of at least $`1015km`$. Much attention has been devoted to the role of overpressurized fluid . It is more and more recognized that fluids play an essential role in virtually all crustal processes. Ref. reviews the historical development of the conciousness among researchers of the ubiquitous presence and importance of fluids within the crust. Numerous examples exist that demonstrate water as an active agent of the mechanical, chemical  and thermal processes that control many geologic processes that operate within the crust . The bulk of available information on the behavior of fluids comes from observations of exposed rocks that once resided at deeper crustal levels. In any case, present day surface exposed metamorphic rocks indicate that, at all crustal levels, fluids have been present in significant volume. Because the porosity of metamorphic rock is probably less than $`1\%`$, the high volume of calculated fluid necessary to produce observed chemical changes suggests that fluid must have been replenished thousands of times. It has also been proposed that gold-quartz vein fields in metamorphic terranes provide evidence for the involvement of large volumes of fluids during faulting and may be the product of seismic processes . Water is also released from transformed minerals. For instance, Montmorillonite changes to illite with a release of free water from the clay structure at approximately the same depth as the first occurrence of the anomalous pore pressure . This is the most commonly discussed example of hydration and dehydration of minerals changing the fluid mass and the pore pressure. Fluids have been directly sampled at about $`11km`$ by the Soviets at the Kola Peninsula drillhole.

Many observations suggest that there are massive crustal fluid displacements correlated with seismic events. Among them, one can cite the fault-valve mechanism  or the migration and diffusion of aftershocks . Several mechanisms have been proposed for the mechanical effect of fluids to decrease compressive lithostatic stresses: mantle-derived source of fluids can maintain overpressure within a leaky fault ; laboratory sliding experiments on granite show that the sliding resistance of shear planes can be significantly decreased by pore sealing and compaction which prevent the communication of fluids between the porous deforming shear zone and the surrounding material . Physico-chemical processes such as mineral dehydration during metamorphism may provide large fluid abundance over large areas . Such fluid presence or migration implies that cracks may open and close even at seismogenic depth, justifying the relevance of asymetric nonlinear elasticity for such depths as we elaborate below.

### 2.2 Physical mechanism for the threshold nonlinearity

From a mechanical point of view, the threshold nonlinearity we invoke is different from the ubiquitous nonlinearity of rocks. In the later, nonlinearity becomes important only under large deformations. Such nonlinear elasticity of rocks is well-documented from its nonlinear wave signatures . In contrast, the nonlinearity we invoke is revealed for almost arbitrary small perturbations by the difference between compressive versus extensive perturbations around a mean lithostatic stress field. Many crustal rocks have a Young’s modulus depending on confining pressure, in particular with a Young’s modulus in tension smaller than the Young’s modulus in compression in a ratio from $`1/2`$ to $`1/10`$ .

The highly damaged upper crust, as described above, would behave as a standard continuous elastic medium if almost all cracks are kept closed and are thus completely transparent to the applied tectonic stress. In a totally dry crust, this would occur at depths larger than about $`1km`$. Near the surface, cracks can open under sufficiently tensile stresses, so that the rock (if we neglect crack growth phenomena) will display low elastic modulii (what we will call later a soft state). Under compressive (or weakly tensile) stresses, those same cracks will close so that the rock will display elastic modulii close to those of uncracked material (the hard state). This simple compressive/tensile asymmetry changes strongly the behaviour of the rock: the rheology is still elastic (in the sense that it is reversible) but it is nonlinear (a change of sign of applied stress does not only changes the sign of strain, it also changes its modulus). In addition, if we now take account of the possible slow growth of cracks under tensile stress, the nature of the nonlinearity will not change, but its amplitude will.

Let us now consider the influence of fluids whose presence are pervasive in the crust, as summarized in section 2.1. Let us assume that the crust is saturated with fluids and that the network of cracks is sufficiently dense so that it behaves as a drained medium. For the sake of simplicity, we will also assume that the crust is characterized by a homogeneous spatial distribution of cracks, so that the permeability is uniform and isotropic. We will also suppose that there exists a horizontal interface at depth $`z_{seal}`$ where permeability is $`0`$. The last ingredient is that the part of the crust below $`z_{\mathrm{seal}}`$ is connected to a reservoir of fluid. Many indirect observations suggest indeed the presence of large sources of fluids . Then, above depth $`z_{\mathrm{seal}}`$, water will be at hydrostatic pressure $`P_h`$ (lower than lithostatic pressure $`P_l`$) so that cracks will be able to close if the tensile applied tectonic stress is smaller than $`P_lP_h`$ in modulus (cracks remain closed if tectonic stress is compressive). Below $`z_{\mathrm{seal}}`$, the scenario is different. Water trapped in cracks is now at lithostatic pressure. If there is no applied tectonic stress, fluid pressure inside cracks compensates exactly the lithostatic pressure and the medium is exactly at the hard/soft boundary. The net stress acting on any crack’s lips is $`0`$, so that the crack doesn’t grow. If the applied tectonic stress is compressive, cracks will close and fluids are expelled towards the reservoir: the crust is in the hard phase. If the applied tectonic stress is tensile, cracks will open and the crust is in the soft phase. If a crack becomes unstable, pressure drops within it, so that it tends to close to re-establish initial fluid pressure. If a crack grows very slowly, pressure within the crack will stay constant with time. We will from now on neglect the possible slow crack growths and assume that fluid pressure is constant in each crack and that the crack network geometry does not show any evolution with time or applied stress. Thus, below $`z_{\mathrm{seal}}`$, the applied stress threshold to get from the hard to the soft state is $`0`$. This explains our terminology of a “threshold nonlinearity.”

The remaining of the paper will investigate stress redistributions associated with earthquakes occurring below $`z_{\mathrm{seal}}`$. The existence of $`z_{\mathrm{seal}}`$ is difficult to prove for the real crust. It would correspond, for instance, to a depth where the crack network geometry changes abruptly, for instance at the boundary between the sedimentary basin and cristalline rocks. We can also propose an alternate model in which the vertical permeability below $`z_{\mathrm{seal}}`$ is much lower than the horizontal one, and would obtain the same kind of soft/hard transition. Note that much more complicated scenarii involving the chemistry of both fluids and rock matrix could be taken into account , but we choosed to neglect them in this purely mechanical paper.

We then propose a simple central-force spring model in two dimensions to study the behaviour of such an asymmetric nonlinear medium when subjected to an infinitesimal internal source of stress or strain. We show that the spatial structure of the stress transfer associated with such a perturbation evolves with the strength of the nonlinearity (defined as the ratio of the springs’ stiffness in tensile and in compressive states). Those changes are quantified in terms of the symmetry of the resulting stress field and of the decay rate of the amplitude of the stress perturbation with distance from the source.

## 3 Numerical Model

Solving theoretical problems of nonlinear elasticity proves to be very tough, even for the simple asymmetry of our threshold nonlinearity. We thus choosed to solve a couple of simple problems related to seismology using numerical modelling. Stress, strain and material rigidity being $`2^{nd}`$ or $`4th`$ rank tensors, and as there is an obvious and complex feedback between strain and rigidity in our model, we choosed to use a simple spring model to illustrate the concept and major consequences of the threshold nonlinear elasticity.

A plate of size $`L`$ by $`L`$ is discretized onto a regular grid of mesh size $`a`$. In the following, we choose $`L=2000km`$ and $`a=10km`$. Each elementary cell is defined by $`4`$ nodes, each node being shared between $`4`$ different neighbouring cells (except on the boundary of the plate). Figure 1 shows the mechanical structure defined for each cell: each node is connected to its two nearest neighbours by springs of stiffness $`K_1`$ (those springs indeed define the $`4`$ edges of the cell). As those are central force springs, the shear modulus of such a cell is indeed $`0`$. To get shear elasticity , $`2`$ springs are added along the cell’s diagonals, such that each node is also connected to its next-nearest neighbours. The stiffness of those diagonal springs is $`K_2`$. Once the plate is discretized with such cells, it can be shown that the plate behaves as an isotropic elastic medium if and only if we have $`K_2=K_1/2`$ . The two independent elastic modulii of the plate can then be shown to be $`\lambda =\mu =K_1/2`$ for the two Lamé coefficients, thus yielding $`E=(4/3)K_1=(8/3)\mu `$ for the Young modulus and $`\nu =1/3`$ for the Poisson coefficient (note that we are dealing with a pure $`2D`$ model). Of course, many other geometries of the elementary cell are possible, but we had to choose square cells which help to handle more easily with boundary conditions used in the problems we want to study.

The relationships we just defined are true only if the springs are symmetric, i.e., their stiffnesses is the same under tensile or compressive states. The next and last step to model our nonlinear threshold elastic rheology is to impose that the stiffness of each spring can vary with its length. Thus, if a spring is shortened, its stiffness will be, say, $`K`$. If a spring is lengthened, its stiffness is lower and taken equal to $`\alpha K`$, with $`\alpha 1`$. Each spring represents an effective volume filled with a uniform and isotropic distribution of cracks with sizes smaller than the representative mesh size. Real damage in the crust is of course much more complex, with anisotropic, space and scale dependence. These complications are neglected in our first exploration. We can obtain an order of magnitude estimate of the density of fluid-filled faults associated with a given asymmetric coefficient $`\alpha `$, using the effective medium calculations in . To simplify, let us assume $`\lambda =\mu `$. Then, $`\alpha =\lambda _d/\lambda _1=1/(1+5d/2)`$, where $`\lambda _d`$ is the Lamé coefficient of the damaged material and $`d=N(\mathrm{}/L)^3`$ is the density of faults (assumed identical) of radius $`\mathrm{}`$ in a cube of volume $`L^3`$. $`N`$ is the number of faults in that volume. For instance, we need about $`11`$ faults of size $`L/3`$ to get $`\alpha =0.5`$. Such estimate must however be taken with caution since the effective medium calculation $`\alpha =1/(1+5d/2)`$ is valid only for small crack densities, while any piece of rock and the real crust are crisscrossed by many faults at many length scales, most of them being healed at varying degrees. We think that values of $`\alpha `$ significantly smaller than $`1`$ should thus not be excluded. It is also probably that $`\alpha `$ is not uniform within the crust and can be expected to reflect the past history of deformations and ruptures.

The ratio $`\alpha `$ of the extensive over compressive elastic coefficient can also be seen as equivalent to $`1D_m`$ in damage mechanics, where $`D_m`$ is the scalar damage variable. If the spring isn’t damaged, then $`D=0`$, so that $`\alpha =1`$ and the stiffness is the same under tensile and compressive states. If the spring is totally damaged (near failure), $`D_m`$ is close to $`1`$, so that $`\alpha 0`$ and the stiffness of the spring in tensile state vanishes, while it is still $`K`$ if the spring is compressed. Under arbitrary loading conditions, some springs in the plate will be in tensile state, while others will be in compressive state, so that it is difficult to analytically compute the stiffness tensor of the whole plate when $`\alpha <1`$. In addition, the stiffness tensor may feature more than $`2`$ independent modulii. This justifies the use of an iterative numerical method as described in the next section.

## 4 Numerical Method

The method we use belongs to the so-called iterative ‘type-writer’ methods. The first step consists in defining boundary conditions, i.e. to fix displacements and/or applied forces on set of nodes. As we are dealing with statics, such forces and displacements are held constants throughout the numerical process. We then consider, say, the node on the top left corner of the mesh and, according to the forces/displacements applied to this node and his nearest and next-nearest neighbours, we can compute the net force acting on that node. As we are dealing with a statics problem, the net force acting on that node at equilibrium should vanish. If the force vector we computed for that node is not $`\stackrel{}{0}`$, then we move it in the direction of the force vector to decrease the net force. We store the new position of this node and get to its right neighbour and follow the same scheme. We thus sweep the mesh line by line down to the node on the bottom right, and iterate the same operations again from the node on top left. As iterations accumulate, the net force acting on each node decreases, and we stop the process once all nodes are subjected to a force whose modulus is under a given threshold. This threshold is chosen such that the modulus of the incremental displacement necessary to decrease the net force modulus on each node is of the order of the accuracy of the computer, namely about $`10^{12}`$. To avoid any spurious result due to the direction of type-writing, iterations alternatively begin on each of the four corners of the plate on either horizontal or vertical directions, right or left, up or down. Once the ‘type-writer’ is stopped, displacements of all nodes relative to their initial positions are stored. Those positions also allow to compute, for a given node, the forces exerted on it by each of its neighbours. Those last quantities allow to compute the full stress tensor at that node , which is also stored. Computations of forces transmitted by a spring from one node to another take account of the spring stiffness which, as described in the previous section, depends on the state (stretched or shortened) of the spring through the value of $`\alpha `$ which is kept constant throughout the network.

Starting from given boundary conditions, we solve a statics problem, which means that we do not take account of neither wave propagation nor fluid migration. We just compute the final equilibrium solution. This point will be discussed at the end of the paper when considering the application of this model to real Earth data.

## 5 Earthquake Modeling

In this paper, we want to examine the main differences between the stress field pattern generated by an earthquake in our nonlinear threshold elastic medium and in a standard linear elastic medium. The first step is thus to define mesh parameters representative of the real crust, the second one is to define what is an earthquake source in such a model.

### 5.1 Parametering the Earth’s crust

Our numerical model is strictly a $`2D`$ one, as dealing with the third dimension would lead to serious memory and computing time problems, and we would have to handle boundary problems such as the free surface and coupling with lower viscous layers. We can however choose mesh parameters such that, using relationships linking plane elasticity to $`3D`$ elasticity, we can model realistic crust properties (however neglecting boundary problems). The size of the plate as well of the cells has been previously given. We fixed the (virtual) thickness of the plate to be $`10km`$, so that it corresponds roughly to the thickness of the seismogenic zone within the crust. We then choosed $`K_c=510^{14}Nm^1`$ (where $`K_c`$ is the value of $`K_1`$ when springs are compressed), so that the $`3D`$ elastic modulii become $`E=6.2510^{10}Pa`$, $`\lambda =\mu =2.510^{10}Pa`$, and $`\nu =0.25`$, which are close to usual modulii measured in rock mechanics experiments. The value of the assymmetry parameter $`\alpha `$ will be varied through several numerical experiments from $`1`$ (standard isotropic elasticity) down to $`0.01`$ (strong asymmetry of the elastic response under extension versus compression).

### 5.2 Earthquake Source modelling

Earthquake source theory has until now been theoretically studied within the framework of linear elasticity. This allows one to use very powerful tools such as the representation theorem and Green functions. An earthquake can then be viewed equivalently as a displacement discontinuity across a fault plane, or a distribution of double-couples and dipoles of forces along the same plane in a continuous medium .

In the case of a fault of finite dimensions, if we assume that the stress drop is uniform along the fault, then we deal with a crack problem. If we assume that the displacement discontinuity across the fault is uniform, then we deal with a dislocation problem. At distances from the fault much larger than its size, and in the case of linear elasticity, both models yield the same spatial patterns of stress and displacement fields, which are linked by Hooke’s law. This will be illustrated below in our diagrams obtained for the symmetric elasticity case $`\alpha =0`$.

In the nonlinear case, it is easy to show that the representation theorem fails to apply, and then so does the Green function concept. This stems from the fact that the principle of linear superposition fails in the presence of nonlinearity. It follows that even the most simple earthquake source problem has to be defined either as a crack or a dislocation problem, and both problems should give different stress and displacement patterns at long wavelengths. We will, in this preliminary work, study only pointwise sources, which will allow simple comparisons with elementary solutions obtained from linear elasticity.

We will thus study two cases:

* an initially stress-free medium within which a single cell (located at its center) is subjected to the following stress field tensor: $`\sigma _{xx}=\sigma _{yy}=0`$ while $`\sigma _{xy}0`$ (and is hereafter called pointwise shear stress load or crack model) - this model is reminiscent of the standard dynamical model of an event in standard linear elasticity ,
* an initially stress-free medium within which a single cell is subjected to a pure shear strain field: $`ϵ_{xx}=ϵ_{yy}=0`$ while $`ϵ_{xy}0`$ (hereafter called pointwise shear strain load or dislocation problem) - this model rather views the event as a shear displacement discontinuity.

In both cases, the corresponding infinitesimal planar defect (the source of the earthquake) suffers from undeterminacy and is oriented either along the $`\stackrel{}{O}_x`$ direction (plane $`P_x`$) or along the $`\stackrel{}{O}_y`$ direction (plane $`P_y`$). In the first case, the slip discontinuity is dextral, and it is sinistral in the second case.

### 5.3 Quantitative source parameters

The small scale of our mechanical model is that of a cell, and this is thus the smallest scale we have to deal with in order to model an earthquake source. Figure 2 shows the source cell and $`4`$ different vectors originating from each of its $`4`$ corners.

In the case of the pure shear stress load model, each vector represents a force applied to the corresponding node. All forces have the same modulus so that the stress tensor within the central cell indeed corresponds to the one we defined above. The same set of forces is applied whatever the value of $`\alpha `$.

In the case of the pure shear strain load model, each vector represents a displacement applied to the corresponding node. All displacements have the same modulus so that the strain tensor within the central cell indeed corresponds to the one we defined above. The same set of displacements is applied whatever the value of $`\alpha `$.

In order to ensure that we can compare results obtained from both types of boundary conditions, we have to fullfil a very simple condition: the stress and displacement fields must be identical in the classic linear case, i.e. when $`\alpha =1`$. In the pure shear stress load model, we imposed the modulus of each applied force equal to $`F=7.110^{14}N`$, so that is corresponds to an event of scalar moment $`M_0=7.110^{18}Nm`$, i.e., of magnitude $`6.5`$. We computed the displacement field for $`\alpha =1`$, and observed that the magnitude of the displacement at each node of the source cell was $`0.526640588m`$. To be perfectly consistent, we thus impose this displacement amplitude at each source cell node in the case of the pure shear strain load model.

## 6 Displacement field at the source

In the pointwise dislocation case, displacements are held constant at the source whatever the value of $`\alpha `$. In the pointwise crack model, only forces are kept constant as $`\alpha `$ is changed, and displacements are expected to vary as the asymmetry increases (i.e. $`\alpha `$ decreases). We have already pointed out that both boundary conditions assume that the mechanical defect is either parallel to $`P_x`$ or to $`P_y`$. Indeed, we will show that the displacement (hereafter named $`u_n`$) normal to each of these conjugate defects is of the same type (opening) while the shear displacement (hereafter $`u_s`$) along each of them is of different type: dextral along $`P_x`$ and sinistral along $`P_y`$. That said, we will focus only on the modulii of those displacements.

How do we obtain $`u_n`$ and $`u_s`$? Indeed, according to the orientations of the conjugate plane defects, we have to compute displacements at points $`A,B,C`$ or $`D`$, whereas our model provides solutions at nodes $`1,2,3`$ and $`4`$ (see Figure 2). Displacements at points $`A`$ to $`D`$ are thus computed through bilinear interpolations within the cell. We then define $`u_n=u_x(B)=u_x(D)=u_y(C)=u_y(A)`$ and $`u_s=u_x(A)=u_y(B)=u_x(C)=u_y(D)`$ as a result of the symmetries of the system.

Figure 3 shows the variations of both $`u_n`$ and $`u_s`$ with $`\alpha `$ (in fact, it shows the values of displacements discontinuities across the crack, i.e. $`2u_n`$ and $`2u_s`$). For $`\alpha =1`$, we find that $`u_n`$ is very close to $`0`$, so that the displacement along the pointwise crack lips is pure shear. As $`\alpha `$ decreases from $`1`$ to close to $`0`$, the shear displacements increase by a factor of about $`4.5`$. This is perfectly understandable, as some springs are sollicited in tension, leading to a decrease in their stiffness. This decrease, in the presence of constant forces, implies that displacements increase.

More surprising is the behaviour of normal displacements, which increase drastically as $`\alpha `$ decreases, tending to be about $`2.7m`$ when $`\alpha `$ tends to $`0`$ (i.e. about half the shear displacements). For $`\alpha =0.1`$ we have $`u_n/u_s`$ close to $`1/3`$. Moreover, $`u_n`$ and $`u_s`$ are both positive, which signifies that, under the shear stress load assumption, the defect opens when the medium is asymmetric. In seismological words, this means that the static moment tensor of the source has a non-vanishing trace and can thus be decomposed into an isotropic part and a deviatoric one. Despite the observation that most earthquake sources are thought to be well modelled by the deviatoric part alone, a few catalogues report isotropic components. Several mechanisms have been invoked to explain a non-vanishing isotropic component of the seismic moment. The standard explanation for non-double-couple components relies on the fault zone irregularity . Some earthquakes with non-double-couple mechanisms have been claimed not to be explained solely by such a composite rupture . It is then important to note that the seismic moment tensor reported in catalogues is a variable that quantifies static displacements at the source. A static dilational strain at the source can occur even when the dynamic representation of the source is a pure double couple of forces (yielding a stress tensor with zero trace).

Figure 3 also reports the shear strain that can be measured at the source cell as well as the relative dilation of its surface $`dS/S`$. Both quantities of course corroborate the previous results on $`u_s`$ and $`u_n`$.

## 7 The stress field

We will now focus on the structure of the stress field generated by our pointwise earthquakes. The stress field will be studied at large wavelengths, i.e. at distances from the source cell of more than a few cell sizes. This gives the rate of stress decay with distance from the source and thus the range of interactions between events. We will in the following implicitly assume that the plate is affected by many other faults that are locked and oriented in the same direction (say, along direction $`P_x`$ with potential dextral displacement along the fault). The source cell is one of such fault producing an event. We then study the effect of this event on all other faults in the plate.

### 7.1 Spatial patterns

Figures 4 to 6 show the variation of the shear stress component $`\sigma _{xy}`$ within the plate near the source cell. A positive variation signifies that this stress component increases. On each figure, the panels on the left represent patterns obtained with the pure shear stress load hypothesis, while the panels on the right represent patterns obtained with the pure shear strain load hypothesis. Each row corresponds to a different value of $`\alpha `$, i.e., to a different degree of elastic asymmetry between extensive and compressive deformation.

In the case $`\alpha =1`$, both boundary conditions yield exactly the same spatial pattern, as expected. This is in agreement with the fact that, in linear elasticity, both boundary conditions are equivalent. We obtain $`8`$ lobes of identical shapes within which stress amplitude alternates from positive to negative (red and blue lobes, repectively).

As $`\alpha `$ decreases, the symmetry of the patterns decreases: some lobes are rotating, some are widening while others are narrowing. We will study stress variation in those lobes in a subsequent section. The most striking observation is however that the spatial structure of the patterns is the same for both types of boundary conditions for a given value of $`\alpha `$, notwithstanding a significantly smaller stress amplitude in the pure shear strain load boundary condition.

Figures 7 to 9 show the variation of the stress component $`\sigma _{xx}`$ within the plate. The same kind of comments apply here as for $`\sigma _{xy}`$. This is also the case for component $`\sigma _{yy}`$ which is shown in Figures 10 to 12.

We have thus introduced a mechanical asymmetry at the microscopical level (i.e. at the spring scale), for all springs’orientations corresponding to an isotropic asymmetry, which translates into a loss of symmetry at the macroscopical scale. We shall quantify this loss of symmetry more precisely in the next section.

### 7.2 Decay of the stress amplitude away from the source

If we consider the center of the source cell as the origin of our frame, every point within the plate can be located using polar coordinates $`(r,\theta )`$, where $`r`$ is the distance to the origin (the source cell center), and $`\theta `$ is an azimuth measured clockwise from the $`\stackrel{}{O}_x`$ axis. We can then, for a fixed $`\theta `$, look at the decay rate with $`r`$ of the modulus of any of the stress tensor components. To avoid problems due to the finite size of the source and of the whole plate, we quantify the decay of the modulus of any stress component by a power-law of the type $`r^\gamma `$, within a distance interval bracketed within a few cell to a few tens of cell sizes. We then repeat the same computation for different values of $`\theta `$. We then change the value of the asymmetry factor $`\alpha `$ and obtain the corresponding dependences $`\gamma (\theta )`$ for each value of $`\alpha `$.

Figures 13 to 15 show the variation of $`\gamma `$ with $`\theta `$, for different values of $`\alpha `$, quantifying the decay of the $`\sigma _{xy}`$ component. Each frame features two curves, one corresponding to the pure shear stress condition at the source (in red), the other one to the pure shear strain condition (in blue).

When $`\alpha =1`$, $`\gamma `$ is very close to $`2`$, which is the theoretical value predicted by planar elasticity. There are some small fluctuations around that value, the largest ones being obtained for values of the azimuth corresponding to a change of sign of the stress, i.e. where the stress itself almost vanishes. In those peculiar directions, the determination of the exponent is very unstable.

When $`\alpha `$ decreases, $`\gamma `$ values can reach values very different from $`2`$. Some of those values correspond to azimuths where stress vanishes (and are thus spurious), but others reflect genuine consequences of the nonlinearity of the medium. One can see that exponents can thus reach values larger than $`2`$ (reflecting a very rapid decay and thus short interaction range), but that they can also get down to values around or lower than $`1`$, leading to very large interaction ranges. When $`\alpha =0.01`$, one should not consider negative values of $`\gamma `$ too seriously as such negative values would imply that the stress increases with distance. The increase is occurring only over a finite distance range and gives way to a decrease at larger distance. The measured exponent is thus only valid at very short distances and is not an asymptotic value. We are unable for those cases to quantify accurately the value of the asymptotic $`\gamma `$ due to finite size effects. For that $`\alpha `$ values, one should thus consider that the asymptotic value of $`\gamma `$ varies en general between $`0`$ and $`1`$.

The found values of $`\gamma `$ as a function of $`\theta `$ also reveal that exponents do not vary significantly with the conditions imposed at the source, which constitutes another surprise. However, stresses obtained in the pure shear strain condition are lower than in the pure shear stress condition, as the prefactor of the power-law decay is found smaller than for the pure shear stress case.

Figures 16 to 18, as well as Figures 19 to 21 show the same results for components $`\sigma _{xx}`$ and $`\sigma _{yy}`$. For both components, and for $`\alpha =1`$, we recover the theoretical value $`\gamma =2`$ for any $`\theta `$. As $`\alpha `$ decreases, the exponents can take very different values, including some which imply very long range decay. We observe again that the exponents do not vary with the type of loading at the source. We also checked that the explaination for negative values of exponents was the same as for $`\sigma _{xy}`$. Thus, for low $`\alpha `$ values, $`\gamma `$ decreases to values between $`0`$ and $`1`$.

## 8 Discussion

This idea of mechanical asymmetry, and/or of the feedback between local damage and stress decay from perturbative sources is not new but, to the best of our knowledge, it is the first time that it is implemented in a real $`2D`$ plane elastic problem applied to Earth mechanics and seismotectonics. For example,  gives the analytical solution for the stress field and for the dependence $`\gamma (\alpha )`$ in a nonlinear asymmetric elastic medium in the case of antiplane mode III loading, which is thus the scalar equivalent to the problem studied here. In the antiplane case, there is only one stress component and a single exponent $`\gamma (\alpha )`$. In this antiplane case, it can be shown analytically that the exponent $`\gamma (\alpha )`$ is indeed decreasing from the value $`2`$ for $`\alpha =1`$ to smaller values as $`\alpha `$ decreases. But there is not dependence on azimuth for this scalar case.

The existence of an asymmtry in the crust elasticity has been proposed on the basis of observations of the Manyi ($`Mw=7.6`$) earthquake . Using SAR interferometry data, Peltzer et al.  interpreted the mismatch between the displacement across each side of the left-lateral strike-slip fault as due to a mechanical asymmetry between dilational and compressional quadrants. Using a first-order perturbative calculation, they estimated a coefficient $`1/4\alpha 1/2`$ to explain the observed displacement asymmetry. However, they did not consider the possibility the asymmetry could modify the long range decay rate of the stress field. This in turn can modify the future seismic history in the neighborhood of this event. Their computation showed that the asymmetric effect was probably confined in the very shallow part of the crust which, if true, implies that the stress transfer at seismogenic depths after this event obeys standard linear elastic solutions.

The fact that this model should be relevant for shallow crustal mechanics is rather obvious (as shown from the previous field example as well as from the short discussion at the beginning of this paper). The important question is to check if this model also holds (at least in limited spatial domains) at depth. In that case, stress transfer case-studies should take account of this asymmetry effect, which can greatly enhance the distance at which a given event can trigger another one. Testing this hypothesis is not simple, as the rheology we assumed is nonlinear, which means that the effect of successive events can not be simply added. Stress field evolutions with time may then have much sharper transitions in space and time than predicted by models involving linear elasticity, a behaviour reminiscent of the mechanics of granular media . The consequence is that, to compare with the standard stress transfer mechanism , we need to know in details the state of stress within the crust prior to an event to map predicted stress transfer lobes onto aftershocks location catalogs.

Another possibility for testing our hypothesis would be to study the statistics of seismic moment tensors of events, as we saw that, if we assume that the source can be describe dynamically as a pure double couple of forces, this tensor should display a non-vanishing isotropic part. However, the seismic moment tensor is computed from seismic wave observations (i.e. of dynamical nature), and not from static displacements in situ at depth. Morover, the time at which the solutions we computed really hold depends on the diffusion properties of fluids in rock. This is why such a way of testing would indeed imply to compute the whole time-dependent dynamical asymmetric poro-elastic solution to really propose quantitative results allowing a reliable comparison.

Other data that could be used to test the model are the SAR data, in the spirit of the work of . Interpretation of SAR data done after a large event could prove the pertinence of our model, provided that stress field evolution after the event is not modified by other events or by the far-field loading of the plate.

Proving or refuting this model is of prime importance for the understanding of the spatio-temporal patterns of earthquakes, including those preceeding a large and potentially destructive event. We already discussed the fact that it could improve the potential of prediction methods based on concepts such as stress transfer. But it has even deeper implications on the hope of predicting large events from the behavior of the statistics of populations of shocks preceeding that event. In another numerical work, Ref.  studied the progressive damage of a fault plane before its macroscopic rupture. Their model employs cellular automaton techniques to simulate tectonic loading, rupture events and strain redistribution. Note that in that case, strain is equivalent to stress. The elastodynamic Green function for stress/strain redistribution is taken to vary as $`1/r^p`$, where $`p`$ is a parameter which is varied. The systems displays two different regimes depending on the $`p`$ value. For $`p2`$, large events are preceeded by a clear power-law acceleration of energy release of the system, together with the growth of strain energy correlations. This is of course reminiscent of the critical earthquake hypothesis . For $`p`$ larger than $`2`$, the trend of energy release before a large event is linear. This means that the lower $`p`$ is, the more predictable is the large event, using time series of precursory energy release. Their model does not map exactly to ours,, but we could expect that if $`\alpha `$ is under a certain threshold (still to be determined), then $`\gamma `$ would be low enough for the critical earthquake scenario to apply, making large events predictable from time series. In the other hand, if in some areas $`\alpha `$ is above that threshold, then the local tectonic domain would belong to the other regime, and large events would be unpredictable. If the predictability of large events relies, as suggested by , on the exponent of the Green function of stress transfer, then we have pointed out a very simple physical mechanism allowing to tune that exponent. Large scale fluctuations of fluid pressure from lithostatic to infra-lithostatic could then explain why in some cases large events are preceeded by strain energy release acceleration, while the opposite holds in other cases.


text:
# A Molecular Hydrodynamic Theory of Supercooled Liquids and Colloidal Suspensions under Shear

## Abstract

We extend the conventional mode-coupling theory of supercooled liquids to systems under stationary shear flow. Starting from generalized fluctuating hydrodynamics, a nonlinear equation for the intermediate scattering function is constructed. We evaluate the solution numerically for a model of a two dimensional colloidal suspension and find that the structural relaxation time decreases as $`\dot{\gamma }^\nu `$ with an exponent $`\nu 1`$, where $`\dot{\gamma }`$ is the shear rate. The results are in qualitative agreement with recent molecular dynamics simulations. We discuss the physical implications of the results.

Recently, there has been an explosion of interest in understanding out of equilibrium properties in supercooled liquids. In general, the nonequilibrium behavior of a glassy system is characterized by a violation of the fluctuation-dissipation theorem and the absence of time translation invariance. In the particular case where the system is subjected to a homogeneous, steady shear flow, time translation invariance is recovered. This simpler nonequilibrium situation is of interest for two reasons. First, understanding the rheological properties of complex fluids such as colloidal suspensions and polymers at a microscopic level is important for the design and control of new materials. On a more fundamental level, it has recently been suggested that for supercooled liquids there are fundamental connections between the standard thermodynamic control variables of temperature and density in the equilibrium case and steady state shear out of equilibriumliu1998b . A major goal of this work is to develop a theory that provides an explicit microscopic connection between the temperature, density and shear rate in a supercooled liquid.

Dense colloidal suspensions are known to exhibit weak shear thinning behavior. Such behavior is predicted for simple liquids as well, but the effect is too small to observe at temperatures well above the glass transition. For supercooled liquids, however, the situation is different. Recent numerical simulations have revealed anomalous rheological behavior in supercooled liquids. Yamamoto and Onukiyamamoto1998c  and Berthier and Barratberthier2002  have simulated supercooled liquids under strong stationary shear flow and observed a characteristic shear dependence of the structural relaxation time and the shear viscosity, $`\tau _\alpha ,\eta \dot{\gamma }^\nu `$, where $`\tau _\alpha `$ is the structural relaxation time, $`\eta `$ is the shear viscosity, $`\dot{\gamma }`$ is the shear rate and the exponent $`\nu `$ is empirically found to range between $`2/3`$ and $`1`$. An abstract schematic approach based on the exactly solvable p-spin spin glass has been proposed and studied by Berthier, Barrat and Kurchanberthier2000b . This model predicts $`\nu =2/3`$ in agreement with the lower bound found in the simulations of Berthier and Barratberthier2002 . Since this model is schematic, it cannot be used to understand in detail the microscopic relationship between fluid structure and dynamics as a function of thermodynamic control variables and external driving.

In this letter, we shall generalize the mode-coupling theory developed to describe the fluctuations in an equilibrium state to that of a system under a stationary shear flow. Our starting point is generalized fluctuating hydrodynamics. Using several approximations, we obtain a closed nonlinear equation for the sheared generalization of the intermediate scattering function. In this letter, we shall consider both normal liquids as well as the overdamped Brownian behavior of a colloidal suspension in the absence of hydrodynamic interactions. Numerical results will only be presented for the Brownian case, but the more general results derived here could be used to make quantitative contact with recent molecular dynamics simulations.

Consider the shear flow given by

$$𝐯_0(𝐫)=𝚪𝐫=(\dot{\gamma }y,0,0),$$
(1)
where $`(𝚪)_{\alpha \beta }=\dot{\gamma }\delta _{\alpha x}\delta _{\beta y}`$ is the velocity gradient matrix. The hydrodynamic fluctuations for density $`\rho (𝐫,t)`$ and the velocity field $`𝐯(𝐫,t)`$ obey the following set of Langevin equationskirkpatrick1986c ; indrani1995 .

$$\begin{array}{cc}& \frac{\rho }{t}=(\rho 𝐯),\hfill \\ & \frac{(\rho 𝐯)}{t}+(\rho \mathrm{𝐯𝐯})=\rho \frac{\delta }{\delta \rho }\zeta _0(𝐯𝐯_0)+𝐟_R,\hfill \end{array}$$
(2)
where $`\zeta _0`$ is the collective friction coefficient for colloidal particles and $`𝐟_R(𝐫,t)`$ is the random force. The $`\zeta _0`$ term is specific for the colloidal case. In the case of atomic liquids, the friction term should be replaced by a stress term which is proportional to the gradient of velocity field multiplied by the position dependent shear viscosities. Both cases, however, lead to the same dynamical behavior at long time scales. We neglect the weak $`𝐫`$ dependence of the friction coefficient that could arise from hydrodynamic interactions between colloidal particles in the Brownian case. The first term in the right hand side of the equation for the momentum is the pressure term and $``$ is the total free energy in a stationary state. Here we assume that the free energy is well approximated by that of the equilibrium form and is given by a well-known expression;

$$\begin{array}{cc}\hfill \beta & \text{d}𝐫\rho (𝐫)\left\{\mathrm{ln}\rho (𝐫)1\right\}\hfill \\ & \frac{1}{2}\text{d}𝐫_1\text{d}𝐫_2\delta \rho (𝐫_1)c(r_{12})\delta \rho (𝐫_2),\hfill \end{array}$$
(3)
where $`\beta =1/k_\text{B}T`$ and $`c(r)`$ is the direct correlation function. Under shear, it is expected that $`c(r)`$ will be distorted and should be replaced by a nonequilibrium, steady state form $`c_{\text{noneq}}(𝐫)`$, which is an anisotropic function of $`𝐫`$. The effect of shear on the static correlation functions has been studiedonuki1997 ; ronis1986  and is found that the distortion is characterized by the Péclet number Pe$`=\dot{\gamma }\sigma ^2/D_0`$, where $`\sigma `$ is the diameter of the particle and $`D_0=k_\text{B}T/\zeta _0`$ is the diffusion constant. Thus, if the Péclet number is small, the above assumption is valid. If necessary, the full anisotropic steady state structure may be used. By linearizing eq.(2) around the stationary state as $`\rho =\rho _0+\delta \rho `$ and $`𝐯=𝐯_0+\delta 𝐯`$, where $`\rho _0`$ is the average density, we obtain the following equations,

$$\begin{array}{cc}& \left(\frac{}{t}𝐤𝚪\frac{}{𝐤}\right)\delta \rho _𝐤(t)=ikJ_𝐤(t),\hfill \\ & \left(\frac{}{t}𝐤𝚪\frac{}{𝐤}+\widehat{𝐤}𝚪\widehat{𝐤}\right)J_𝐤(t)\hfill \\ & =\frac{ik}{m\beta S(k)}\delta \rho _𝐤(t)\hfill \\ & \frac{1}{m\beta }_𝐪i\widehat{𝐤}𝐪c(q)\delta \rho _{𝐤𝐪}(t)\delta \rho _𝐪(t)\frac{\zeta _0}{m}J_𝐤(t)+𝐟_{R𝐤}(t),\hfill \end{array}$$
(4)
where $`c(q)`$ is the Fourier transform of $`c(r)`$, $`\widehat{𝐤}𝐤/|𝐤|`$, $`J_𝐤(t)=\rho _0\widehat{𝐤}\delta 𝐯_𝐤(t)`$ is the longitudinal momentum fluctuation, and $`_𝐪\text{d}𝐪/(2\pi )^d`$ for a $`d`$-dimensional system. Note that our approximate equation does not contain coupling to transverse momentum fluctuations even in the presence of shear.

In order to construct equations for the appropriate correlations from the above expressions, an approximate symmetry is necessary. In the presence of shear, translational invariance is violated. In other words, correlations of arbitrary fluctuations, $`f(𝐫,t)`$ and $`g(𝐫,t)`$, do not satisfy $`f(𝐫,t)g(𝐫^{},0)f(𝐫𝐫^{},t)g(\mathrm{𝟎},0)`$. Instead, we shall assume that the following symmetry is validonuki1997 ,

$$f(𝐫,t)g(𝐫^{},0)=f(𝐫𝐫^{}(t),t)g(\mathrm{𝟎},0),$$
(5)
where we defined the time-dependent position vector by $`𝐫(t)\mathrm{exp}[𝚪t]𝐫=𝐫+\dot{\gamma }ty\widehat{𝐞}_x`$, where $`\widehat{𝐞}_x`$ is an unit vector oriented along the $`x`$-axis. In wave vector space, this is expressed as

$$\begin{array}{cc}\hfill f_𝐤(t)g_𝐤^{}^{}(0)& =f_𝐤(t)g_{𝐤(t)}^{}(0)\times \delta _{𝐤(t),𝐤^{}}\hfill \\ & =f_{𝐤^{}(t)}(t)g_𝐤^{}^{}(0)\times \delta _{𝐤,𝐤^{}(t)},\hfill \end{array}$$
(6)
where $`𝐤(t)=\mathrm{exp}[{}_{}{}^{t}𝚪t]𝐤=𝐤+\dot{\gamma }tk_x\widehat{𝐞}_y`$, where $`{}_{}{}^{t}𝚪`$ denotes the transpose of $`𝚪`$ and $`\delta _{𝐤,𝐤^{}}(2\pi )^dV^1\delta (𝐤𝐤^{})`$ for a system of volume $`V`$.

Eq.(5) states that the fluctuations satisfy translational invariance in a reference frame flowing with the shear contours. This approximation holds for long wavelengths where the direct interactions between particles are not important. On the other hand, for correlations between particles separated by molecular length scales, this is not generally true. The validity of this approximation for molecular length scales should be systematically examined in the future.

Using this approximation, it is straightforward to construct the mode-coupling equations for the appropriate correlation functions. We shall derive the equation for the intermediate scattering function defined by

$$F(𝐤,t)\frac{1}{N}\delta \rho _{𝐤(t)}(t)\delta \rho _𝐤^{}(0),$$
(7)
where $`N`$ is the total number of the particles in the system. Note that the wave vector in $`\delta \rho _𝐤(t)`$ is now replaced by a time-dependent one $`𝐤(t)`$.

Eq.(4) has a nonlinear term

$$R_𝐤(t)=\frac{1}{m\beta }_𝐪i(\widehat{𝐤}𝐪)c(𝐪)\delta \rho _{𝐤𝐪}(t)\delta \rho _𝐪(t).$$
(8)
This term can be renormalized with the definition of a generalized friction coefficient following the standard procedure of derivation of the mode-coupling equationskirkpatrick1986c ; zwanzig2001 . To lowest order in the fluctuations, we have

$$\begin{array}{cc}& \left(\frac{}{t}𝐤𝚪\frac{}{𝐤}+\widehat{𝐤}𝚪\widehat{𝐤}\right)J_𝐤(t)=\frac{ik}{m\beta S(k)}\delta \rho _𝐤(t)\hfill \\ & \frac{1}{m}_{\mathrm{}}^t\text{d}t^{}\text{d}𝐤^{}\zeta (𝐤,𝐤^{},tt^{})J_𝐤^{}(t^{})+𝐟_{R𝐤}^{}(t),\hfill \end{array}$$
(9)
where $`\zeta (𝐤,𝐤^{},t)`$ is the generalized friction coefficient and $`𝐟_{R𝐤}^{}(t)`$ is a corresponding random force. $`\zeta (𝐤,𝐤^{},t)`$ is given by the sum of the bare friction coefficient and the mode-coupling term as

$$\zeta (𝐤,𝐤^{},t)=\zeta _0\times 2\delta (t)+\delta \zeta (𝐤,𝐤^{},t)$$
(10)
with the mode-coupling contribution given by

$$\delta \zeta (𝐤,𝐤^{},t)=\frac{m^2\beta }{N}R_𝐤(t)R_𝐤^{}^{}(0).$$
(11)
Substituting eq.(8) into the above expression, we obtain

$$\begin{array}{cc}& \delta \zeta (𝐤,𝐤^{}t)=\frac{1}{\beta N}_𝐪_𝐪^{}\widehat{𝐤}𝐪c(𝐪)\widehat{𝐤}^{}𝐪^{}c(𝐪^{})\hfill \\ & \times \delta \rho _{𝐤𝐪}(t)\delta \rho _𝐪(t)\delta \rho _{𝐤^{}𝐪^{}}^{}(0)\delta \rho _𝐪^{}^{}(0).\hfill \end{array}$$
(12)
This involves a four point correlation function. Using the Gaussian approximation, this can be decomposed into a product of two-point correlation functions as

$$\begin{array}{cc}& \delta \rho _{𝐤𝐪}(t)\delta \rho _𝐪(t)\delta \rho _{𝐤^{}𝐪^{}}^{}(0)\delta \rho _𝐪^{}^{}(0)\hfill \\ & N^2F(𝐤(t)𝐪(t),t)F(𝐪(t),t)\hfill \\ & \times \left\{\delta _{𝐤(t),𝐤^{}}\delta _{𝐪(t),𝐪^{}}+\delta _{𝐤(t),𝐤^{}}\delta _{𝐪(t),𝐤^{}𝐪^{}}\right\},\hfill \end{array}$$
(13)
where use has been made of the translational invariance, eq.(6).

Substituting eq.(13) back to eq.(12), we obtain

$$\delta \zeta (𝐤,𝐤^{},t)=\delta \zeta (𝐤,t)\delta _{𝐤(t),𝐤^{}}$$
(14)
with

$$\begin{array}{cc}& \delta \zeta (𝐤,t)=\frac{\rho _0}{\beta }_𝐪\left\{\widehat{𝐤}𝐪c(𝐪)+\widehat{𝐤}\left(𝐤𝐪\right)c(𝐤𝐪)\right\}\hfill \\ & \times \widehat{𝐤}(t)𝐪(t)c(𝐪(t))F(𝐤(t)𝐪(t),t)F(𝐪(t),t)\hfill \\ & =\frac{\rho _0}{2\beta }_𝐪𝒱(𝐤,𝐪)𝒱(𝐤(t),𝐪(t))F(𝐤(t)𝐪(t),t)F(𝐪(t),t),\hfill \end{array}$$
(15)
where $`𝒱(𝐤,𝐪)`$ is the vertex function given by

$$𝒱(𝐤,𝐪)=\widehat{𝐤}\left\{𝐪c(𝐪)+\left(𝐤𝐪\right)c(𝐤𝐪)\right\}.$$
(16)
From these results and eq.(9), the equation for the correlation function,

$$C(𝐤,t)\frac{1}{N}J_{𝐤(t)}(t)n_𝐤^{}(0)$$

is given by

$$\begin{array}{cc}& \frac{\text{d}C(𝐤,t)}{\text{d}t}\hfill \\ & =\widehat{𝐤}(t)𝚪\widehat{𝐤}(t)C(𝐤,t)\frac{ik(t)}{m\beta S(k(t))}F(𝐤,t)\hfill \\ & \frac{1}{m}_0^t\text{d}t^{}\delta \zeta (𝐤(t),tt^{})C(𝐤,t^{}).\hfill \end{array}$$
(17)
Note that in the above equation, the differential operator $`𝐤𝚪/𝐤`$ disappears because

$$\frac{\text{d}C(𝐤,t)}{\text{d}t}=\frac{C(𝐤,t)}{t}𝐤(t)𝚪\frac{}{𝐤(t)}C(𝐤,t).$$
(18)
Likewise, the continuity equation (the first term in eq.(4)) can be written as

$$\begin{array}{cc}\hfill \frac{\text{d}F(𝐤,t)}{\text{d}t}& =ik(t)C(𝐤,t).\hfill \end{array}$$
(19)
This equation together with eq.(17) comprises the closed set of the mode-coupling equations for $`F(𝐤,t)`$ and $`C(𝐤,t)`$ under shear.

For colloidal suspensions the relaxation time of the momentum fluctuations is of the order of $`\tau _m=m/\zeta _0`$ and is much shorter than the relaxation time for density fluctuations which is of the order of or longer than $`\tau _d=\sigma ^2/D_0`$. For the time scale of interest, $`\widehat{𝐤}(t)𝚪\widehat{𝐤}(t)`$ as well as the inertial term can be neglected in the equation for the momentum fluctuations since $`\dot{\gamma }\tau _m1`$ in realistic situations. Thus, the equation for the momentum fluctuations may be written as

$$\begin{array}{cc}\hfill 0=& \frac{ik(t)}{\beta S(k(t))}F(𝐤,t)\zeta _0C(𝐤,t)\hfill \\ & _0^t\text{d}t^{}\delta \zeta (𝐤(t),tt^{})C(𝐤,t^{}).\hfill \end{array}$$
(20)
Substituting this back into eq.(19), we arrive at

$$\begin{array}{cc}\hfill \frac{\text{d}F(𝐤,t)}{\text{d}t}=& \frac{D_0k(t)^2}{S(k(t))}F(𝐤,t)\hfill \\ & _0^t\text{d}t^{}M(𝐤(t),tt^{})\frac{\text{d}F(𝐤,t^{})}{\text{d}t^{}},\hfill \end{array}$$
(21)
where

$$\begin{array}{cc}\hfill M(𝐤,t)=& \frac{\rho _0D_0}{2}\frac{k}{k(t)}_𝐪𝒱(𝐤,𝐪)𝒱(𝐤(t),𝐪(t))\hfill \\ & \times F(𝐤(t)𝐪(t),t)F(𝐪(t),t).\hfill \end{array}$$
(22)
Eqs.(21) and (22) are the major result of this letter. In the absence of the shear, they reduce to the conventional mode-coupling equationsgotze1992 .

In order to study the shear thinning effect in the supercooled state, we shall solve eqs.(21) and (22) numerically. Solving this equation is more difficult than solving the corresponding equation in the equilibrium case because the wave vectors are distorted by shear and the system is not isotropic. For simplicity, we shall consider a hypothetical two-dimensional colloidal suspension which is simple to handle numerically but still undergoes an ergodic-nonergodic transition below a certain density. The shear flow occurs in the $`x`$ direction. We have chosen the following form of the static structure factor $`S(k)`$:

$$S(k)=S_{\text{PY}}(k+k_0,\alpha \rho _0)f(kk_c),$$
(23)
where $`S_{\text{PY}}(k,\rho )`$ is the static structure factor for a hard-sphere system at the density $`\rho `$ obtained from the Percus-Yevick closure, $`k_0`$ and $`\alpha `$ are parameters which were chosen in such a way that $`S(k)`$ is short-ranged and has broader peak. $`f(kk_c)`$ is a cut-off function which makes $`S(k)`$ approach unity smoothly for wave vectors larger than the cut-off $`k_c`$. The choice of $`S(k)`$ mimics the shape of $`S(k)`$ of real systems although it does not satisfy sum-rule restrictions. In our calculation, we chose $`k_0=4.0`$, $`\alpha =32`$ and $`k_c=4.0`$. For this system, the ergodic-nonergodic transition occurs around a “density” $`\rho _c\sigma ^21.2\times 10^2`$ in the absence of shear. In Figure 1, we show the behavior of $`F(𝐤,t)`$ for $`\rho \sigma ^2=1.15\times 10^2`$, slightly below $`\rho _c`$. The wave vector is $`𝐤=(0,2\sigma )`$. Since $`k_x=0`$, the expression for $`F(𝐤,t)`$ is equivalent to that in the quiescent state, $`F_{\text{eq}}(𝐤,t)N^1\delta \rho _𝐤(t)\delta \rho _𝐤^{}(0)`$. Thus, there is no direct effect from the convection term but, due to the nonlinear coupling through the mode-coupling term $`M(𝐤,t)`$, a strong shear dependence of the relaxation time can be seen. The dependence of the structural relaxation time $`\tau _\alpha (\dot{\gamma })`$ on the shear rate $`\dot{\gamma }`$ is estimated from the value where $`F(𝐤,\tau _\alpha )=\text{e}^1`$. For the particular case of $`\rho \sigma ^2=1.15\times 10^2`$ we find the power law $`\tau _\alpha \dot{\gamma }^\nu `$ with $`\nu 0.8`$ for Pe$`10^3`$ and $`\tau _\alpha `$ saturates to the equilibrium value at Pe$`10^3`$. This is similar to the behavior reported in recent molecular dynamics simulations. For values of $`\rho \sigma ^2>1.2\times 10^2`$ (not shown), we find that the exponent $`\nu `$ saturates at a higher value, in agreement with the simulations of Berthier and Barratberthier2002 .

The physical picture that emerges from the molecular hydrodynamic theory developed here is simple. The shear flow perturbs and randomizes the coupling between different modes. Physically, this perturbation dissipates the cage that transiently immobilizes particles. Mathematically, this is reflected through the time dependence of the vertex, which vanishes as $`t\mathrm{}`$. This simple picture illustrates the essence of the mode-coupling approach to the shear thinning effect in simple supercooled systems. Note that even for fluctuations orthogonal to the direction of flow, thinning occurs due to the coupling of fluctuations in all directions. In this sense, the picture of cage breakup in a supercooled liquid due to external flow is quite different from that of dynamic critical phenomena under shear, in which the faster relaxation occurs solely because the fluctuations are stretched out by the shear flow and pushed to larger wave vectors where faster relaxation occurs.

In this letter, we have derived an approximate mode-coupling theory for a supercooled liquid under steady shear flow. The most important assumption is the use of approximate translational invariance, eq.(6). This allows one to derive a nonlinear integro-differential equation for $`F(𝐤,t)`$ similar to the one for the equilibrium state. The numerical analysis for a hypothetical two-dimensional colloidal suspension has been carried out and a typical behavior of $`F(𝐤,t)`$ was shown to be consistent with recent simulations. The relaxation time is found to have the strong shear dependence. More systematic and thorough analysis of the numerical solutions of the mode-coupling equations are left for future work111 As this manuscript was being written for publication, we became aware of a parallel effort by Fuchs and Cates \[cond-mat/0204628;cond-mat/0207530\]. While the omission of several details in their calculation make a detailed comparison between the two approaches impossible at this time, there are many similarities in the final expressions, and the physical picture is also quite similar. Differences do exist in the expressions, which most likely lead to different estimates of the exponent $`\nu `$. It will be most useful in the future to make a detailed comparison of the two theories since the differences in approach (fluctuating hydrodynamics verses the projection operator formalism) and the final results may lead to a better understanding of the limitations of several assumptions used to here and by Fuchs and Cates to describe how an external flow disrupts jamming..

The authors acknowledge support from NSF grant #0134969. The authors would like to express their gratitude to Dr. Ryoichi Yamamoto for suggesting this problem during his stay at Harvard, and for useful discussions.

FIGURE CAPTIONS