text:
# Monopolar Optical Orientation of Electronic Spins in Semiconductors

## 1 Introduction

Absorption of circularly polarized light in semiconductors may result in spin polarization of photoexcited carriers. This phenomenon of optical orientation is well known for interband transitions in semiconductors . At interband excitation with circularly polarized light transitions from the valence to the conduction band are allowed only if the angular momentum is changed by $`\pm 1`$. These selection rules lead to the spin orientation of carriers with the sign and degree of polarization depending on the light helicity. While optical orientation at interband excitation has been widely studied, it is not obvious that the absorption of infrared radiation due to intraband optical transitions can result in a spin polarization.

In this paper we show that free carrier absorption of circularly polarized radiation due to both, indirect optical transitions for bulk semiconductors or quantum well structures, and direct transitions between size-quantized subbands also leads to spin orientation of free electrons. This optical orientation has not been considered previously and may be referred to as ‘monopolar spin orientation’ because only one type of carriers is excited, electrons in $`n`$-type materials, holes in $`p`$-type materials. We present theoretical and experimental results on the monopolar optical orientation of intraband absorption by free carriers in $`n`$-type zinc-blende-structure bulk semiconductors and low dimensional structures. We show that the monopolar optical orientation of electrons can be obtained if an admixture of the $`\mathrm{\Gamma }_7`$ and $`\mathrm{\Gamma }_8`$ valence band states to the conduction band wave functions is taken into account.

## 2 Monopolar optical spin orientation in bulk semiconductors

Absorption of infrared light by free electrons in $`n`$-type semiconductors (Drude-like absorption) occurs by indirect intraband transitions where momentum conservation law is satisfied due to acoustic phonons, optical phonons, static defects etc. The intraband absorption of the light involving both electron-photon interaction and scattering is described by second-order processes with virtual intermediate states. The compound matrix elements for such kind of transitions with the initial and final states in the same band has the standard form

$$M_{cs^{}𝒌^{}cs𝒌}=\underset{\nu }{}\left(\frac{V_{cs^{}𝒌^{},\nu 𝒌}R_{\nu 𝒌,cs𝒌}}{E_{\nu 𝒌}E_{c𝒌}\mathrm{}\omega }+\frac{R_{cs^{}𝒌^{},\nu 𝒌^{}}V_{\nu 𝒌^{},cs𝒌}}{E_{\nu 𝒌^{}}E_{c𝒌}\pm \mathrm{}\mathrm{\Omega }_{𝒌𝒌^{}}}\right).$$
(1)
Here $`E_{c𝒌}`$, $`E_{c𝒌^{}}`$ and $`E_\nu `$ are the electron energies in the initial $`|c,s,𝒌`$, final $`|c,s^{},𝒌^{}`$ and intermediate $`|\nu `$ states, respectively, $`s`$ is the spin index, $`𝒌`$ is the electron wavevector, $`R`$ is the matrix element of electron interaction with the electromagnetic wave, $`V`$ is the matrix element of electron-phonon or electron-defect interaction, and $`\mathrm{}\mathrm{\Omega }_{𝒌𝒌^{}}`$ is the energy of the involved phonon. At static defect assisted scattering $`\mathrm{}\mathrm{\Omega }=0`$. The sign $`\pm `$ in Eq. (1) correspond to emission and absorption of phonons. A dominant contribution to the light absorption due to indirect transitions in the conduction band is caused by processes with intermediate states in the same band (see Fig. 1). The corresponding matrix element has the form

$$M_{cs^{}𝒌^{}cs𝒌}\frac{V_{cs^{}𝒌^{},cs𝒌}}{\mathrm{}\omega }\frac{m_0}{\mathrm{}}𝒆\left(\frac{E_{c𝒌^{}}}{𝒌^{}}\frac{E_{c𝒌}}{𝒌}\right),$$
(2)
where $`m_0`$ is the free electron mass, and $`𝒆`$ is the unit vector of the electric field polarization. The absolute value of the matrix element in Eq. (2) is independent of the degree of circular polarization of radiation, $`P_{circ}`$. Hence, the intraband transitions with intermediate states in the same band do not contribute to the optical orientation.

Spin orientation caused by intraband absorption of circularly polarized light can be obtained considering processes with intermediate states in the valence band and taking into account its spin-orbit splitting. Fig. 2 demonstrates schematically the spin orientation at intraband absorption of right handed circularly polarized light ($`\sigma ^+`$). Because of the selection rules for interband optical matrix elements, the electron transitions with spin reversal $`1/2+1/2`$ are possible via intermediate states in the light-hole and spin-orbit split subbands, while the opposite processes, $`+1/21/2`$ are forbidden. We assume the photon energy $`\mathrm{}\omega `$ to exceed the typical energy of equilibrium electrons, $`\mathrm{}\omega k_BT`$ for a nondegenerate electron gas at the temperature $`T`$, or $`\mathrm{}\omega \epsilon _F`$ for the case of degenerate statistics with the Fermi energy $`\epsilon _F`$. Considering acoustic-phonon-assisted light absorption as an example, the following equation for the spin generation rate is obtained

$$\dot{𝑺}=\frac{1}{6}\frac{\mathrm{\Xi }_{cv}^2}{\mathrm{\Xi }_c^2}\frac{\mathrm{\Delta }_{so}^2}{E_g(E_g+\mathrm{\Delta }_{so})(3E_g+2\mathrm{\Delta }_{so})}K\widehat{𝒆}IP_{circ}.$$
(3)
Here $`\mathrm{\Xi }_{cv}`$ and $`\mathrm{\Xi }_c`$ denote the interband and intraband deformation potential constants, respectively, $`E_g`$ stands for the energies of the band gap and $`\mathrm{\Delta }_{so}`$ for the valence band spin-orbit splitting, $`I`$ is the light intensity, and $`\widehat{𝒆}`$ is the unit vector in the light propagation direction. The factor $`K`$ in Eq. (3) is the coefficient of the light absorption under phonon-assisted intraband optical transitions. It is determined by the dominant processes with intermediate states in the same band, and has the form

$$K=\frac{4\alpha }{3n_\omega }\left(\frac{\mathrm{\Xi }_c}{\mathrm{}\omega }\right)^2\frac{k_BT}{\rho v_s^2}\left(\frac{2m^{}\omega }{\mathrm{}}\right)^{1/2}N_e,$$
(4)
where $`\alpha `$ and $`n_\omega `$ are the fine structure constant and the refraction index of the medium, respectively, $`\rho `$ is the mass density of the crystal, $`v_s`$ is the sound velocity, $`m^{}`$ is the effective mass of the conduction electron, and $`N_e`$ is the carrier concentration.

It should be emphasized that the spin generation depends strongly on the energy of spin-orbit splitting of the valence band, $`\mathrm{\Delta }_{so}`$. It is due to the fact that the contributions to the matrix element giving rise to the optical orientation from the valence band $`\mathrm{\Gamma }_8`$ and spin-orbit split band $`\mathrm{\Gamma }_7`$ have the opposite signs. An estimation for GaAs material shows that the typical value of spin generated per one absorbed photon is of order $`10^6`$ at a photon energy of $`\mathrm{}\omega =10`$ meV.

## 3 Monopolar optical spin orientation in QWs

### 3.1 Experimental technique

Experimentally, the monopolar spin orientation has been investigated on MBE (001)-grown $`n`$-GaAs/AlGaAs QW of width $`d_W=7`$ nm and $`n`$-In0.2Ga0.8As QWs of 7.6 nm width. Samples with free carrier densities of about $`210^{11}`$ cm-2 were studied in the temperature range from liquid helium to room temperature. A pair of ohmic contacts was centered on opposite sample edges along the direction $`x[1\overline{1}0]`$ (see in Fig. 3).

A high power pulsed mid-infrared (MIR) TEA-CO2 laser and a far-infrared (FIR) NH3 laser were used as radiation sources delivering 100 ns pulses with radiation power $`P`$ up to 100 kW. Several lines of the CO2 laser between 9.2 $`\mu `$m and 10.6 $`\mu `$m and of the NH3-laser  between $`\lambda `$ = 76 $`\mu `$m and 280 $`\mu `$m were chosen for excitation in the MIR and FIR range, respectively. Excitation of samples by FIR radiation with photon energy less than the separation of size-quantized subbands leads to absorption caused by indirect optical transitions in the lowest subband (Drude absorption). The MIR radiation induces direct optical transitions between the first and the second subband of QWs.

The laser light polarization was modified from linear to circular using a Fresnel rhombus and quartz $`\lambda /4`$ plates for MIR and FIR radiation, respectively. The helicity of the incident light was varied according to $`P_{circ}=\mathrm{sin}2\phi `$ where $`\phi `$ is the angle between the initial plane of linear polarization and the optical axis of the polarizer. Spin polarization has been investigated making use of the circular photogalvanic effect (CPGE)  and the spin-galvanic effect (SGE) . The experimental procedure is sketched in Fig 3. For investigation of the spin-galvanic effect an in-plane magnetic field $`B`$ up to 1 T has been applied. The current $`j`$ generated by polarized light in the unbiased structures was measured via the voltage drop across a 50 $`\mathrm{\Omega }`$ load resistor in a closed circuit configuration. The voltage was recorded with a storage oscilloscope. The measured current pulses of 100 ns duration reproduce the temporal structure of the laser pulses.

### 3.2 Intra-subband transitions in QWs

In quantum well structures absorption of infrared radiation may be achieved by indirect intra-subband optical transitions and, for photon energies being in resonance with the energy distance between size quantized subbands, by direct transitions between these subbands. To obtain absorption caused by indirect transitions we used FIR radiation in the range from 76 $`\mu `$m to 280 $`\mu `$m (corresponding photon energies are from 16 meV to 4.4 meV). The experiments were carried out on GaAs and InAs QWs with the energy separation $`\mathrm{\Delta }E=E_2E_1`$ between $`e1`$ and $`e2`$ size-quantized subbands equal to 120 meV and 112 meV, respectively. Therefore the photon energies used are much smaller than $`\mathrm{\Delta }E`$ and absorption is caused by indirect intra-subband optical transitions. With illumination of (001)-grown GaAs and InAs QWs at oblique incidence of FIR radiation a current signal proportional to the helicity $`P_{circ}`$ has been observed (see Fig. 4a) indicating the circular photogalvanic effect . At normal incidence of radiation, where the CPGE vanishes, the spin-galvanic current  is also observed applying an in-plane magnetic field (see Fig. 4b). Both effects are due to spin orientation. Therefore the observation of the CPGE and the spin-galvanic effect gives clear evidence that the absorption of far infrared circularly polarized radiation results in spin orientation. We note, that monopolar orientation has also been observed for $`p`$-type QW structures, but this is out of the scope of the present paper.

Monopolar optical orientation caused by indirect intra-subband transitions in QWs can be obtained in the similar way as that for bulk semiconductors (see previous section). For this particular mechanism of the monopolar optical orientation, the spin generation rate has been derived for two scattering processes: acoustic phonons and elastic scattering on static defects.

Considering acoustic-phonon-assisted processes with virtual intermediate states in the heavy-hole $`hh`$, light-hole $`lh`$, and spin-orbit split $`so`$ subband, we obtain for the spin generation rate

$$\dot{𝑺}_{ph}=\frac{1}{6}\frac{\mathrm{\Xi }_{cv}^2}{\mathrm{\Xi }_c^2}\frac{\mathrm{\Delta }_{so}^2}{E_g(E_g+\mathrm{\Delta }_{so})(3E_g+2\mathrm{\Delta }_{so})}\left[\widehat{𝒆}_{}+\frac{d_W}{3}\sqrt{\frac{2m^{}\omega }{\mathrm{}}}\widehat{𝒆}_z\right]\eta _{ph}IP_{circ},$$
(5)
where $`d_W`$ is the width of a QW, $`\eta _{ph}`$ is the fraction of the energy flux absorbed in a QW under normal incidence. It is given by

$$\eta _{ph}=\frac{3\pi \alpha }{n_\omega }\left(\frac{\mathrm{\Xi }_c}{\mathrm{}\omega }\right)^2\frac{k_BT}{\rho av_s^2}N_e,$$
(6)
where $`N_e`$ is the concentration of the two-dimensional electron gas.

In the case of elastic scattering on static impurities, the expression for spin generation rate has the form

$$\dot{𝑺}_{imp}=\frac{2}{3}\frac{\mathrm{\Delta }_{so}^2}{E_g(E_g+\mathrm{\Delta }_{so})(3E_g+2\mathrm{\Delta }_{so})}\left[\frac{V_{}^2}{V_0^2}\widehat{𝒆}_{}+\frac{V_z^2}{V_0^2}\widehat{𝒆}_z\right]\eta _{imp}IP_{circ}.$$
(7)
Here $`V_0`$ is the intraband matrix element of scattering, the factors $`V_z`$ and $`V_{}`$ describe impurity-induced mixing of the conduction band Bloch function $`S`$ with the valence band Bloch functions $`Z`$ and $`X,Y`$ respectively. The fraction of the energy flux absorbed in the QW at impurity-assisted intrasubband optical transitions under normal incidence is given by

$$\eta _{imp}=\frac{2\pi \alpha }{n_\omega }\left(\frac{V_0}{\mathrm{}\omega }\right)^2N_eN_d,$$
(8)
where $`N_d`$ is the concentration of scatterers.

### 3.3 Direct inter-subband transitions in QWs

Absorption of radiation in the range from 9 $`\mu `$m to 11 $`\mu `$m in our GaAs and InAs samples is dominated by resonant direct inter-subband optical transitions between the first and the second subband. Fig. 5b shows the resonance behaviour of absorption measured in GaAs QWs by making use of transmission Fourier spectroscopy in a multiple-reflexion waveguide geometry (see inset Fig. 5b). Applying MIR radiation of the CO2 laser, which causes direct transitions in GaAs and InAs QWs, the circular photogalvanic current at oblique incidence (Fig. 5a) and the spin-galvanic current at normal incidence of radiation (Fig. 5b) have also been observed. The wavelength dependence of the spin-galvanic effect obtained between 9.2 $`\mu `$m and 10.6 $`\mu `$m repeats the spectral behaviour of direct inter-subband absorption. This unambiguously demonstrates that in this case the spin orientation of $`n`$-type QWs is obtained by inter-subband transitions.

We would like to emphasize that spin sensitive inter-subband transitions in $`n`$-type QWs have been observed at normal incidence when there is no component of the electric field of the radiation normal to the plane of the QWs. Generally it is believed that inter-subband transitions in $`n`$-type QWs can only be excited by infrared light polarized in the growth direction $`z`$ of the QWs . Furthermore such transitions are spin insensitive and, hence, do not lead to optical orientation. Since the argument, leading to these selection rules, is based on the effective mass approximation in a single band model, the selection rules are not rigorous.

In order to explain the observed spin orientation as well as the absorption of light polarized in the plane of the QW we show that a $`𝐤𝐩`$ admixture of valence band states to the conduction band wave functions has to be taken into account. Calculations yield that inter-subband absorption of circularly polarized light propagating along $`z`$ induces only spin-flip transitions resulting in 100% optical orientation of photoexcited carriers. In this geometry the spin generation rate for resonant inter-subband transitions has the form

$$\dot{S}_z=\frac{128\alpha }{9n_\omega }\frac{\mathrm{\Delta }_{so}^2(2E_g+\mathrm{\Delta }_{so})^2E_1}{E_g^2(E_g+\mathrm{\Delta }_{so})^2(3E_g+2\mathrm{\Delta }_{so})^2}\frac{\mathrm{}^2N_e}{m^{}}\delta (\mathrm{}\omega E_1+E_2)IP_{circ}$$
(9)
where the function $`\delta `$ describes the resonant behaviour of the inter-subband transitions.

## 4 Summary

In conclusion, our results demonstrate that in $`n`$-type bulk semiconductors as well as QW structures monopolar spin orientation can be achieved applying circularly polarized radiation with photon energies less than the fundamental energy gap. Spin orientation has been observed by indirect intra-subband absorption as well as by inter-subband transitions. It is shown that monopolar spin orientation in $`n`$-type materials becomes possible if an admixture of valence band states to the conduction band wave function and the spin-orbit splitting of the valence band are taken into account. We emphasize that the spin generation rate under monopolar optical orientation depends strongly on the energy of spin-orbit splitting of the valence band, $`\mathrm{\Delta }_{so}`$. It is due to the fact that the valence band $`\mathrm{\Gamma }_8`$ and the spin-orbit split band $`\mathrm{\Gamma }_7`$ contribute to the matrix element of spin-flip transitions with opposite signs.

## 5 Acknowledgments

We thank D. Weiss and W. Wegscheider for many helpful discussions. Financial support from the DFG, the RFFI, the Russian Ministry of Science and the NATO linkage program is gratefully acknowledged.

## References


format: text

text:
# Boundary Lubrication: Squeeze-out Dynamics of a Compressible 2D Liquid

## I Introduction

Sliding friction is one of the oldest problems in physics, and has undoubtedly a huge practical importanceP1a ; P1b ; P1c . In recent years, the ability to produce durable low-friction surfaces and lubricants has become an important factor in the miniaturization of moving components in technologically advanced devices. For such applications, the interest is focused on the stability under pressure of thin lubricant films, since the complete squeeze-out of the lubricant from an interface may give rise to cold-welded junctions, resulting in high friction and catastrophically large wear.

It has been shown both experimentally and theoretically that when simple fluids (quasi-spherical molecules and linear hydrocarbons) are confined between atomically flat surfaces at microscopic separations, the behavior of the lubricant is mainly determined by its interaction with the solids that induce layering in the perpendicular direction Israel3a ; Israel3b ; Gao97a ; Gao97b ; Demi96a ; Demi96b ; Kleina ; Kleinb . The thinning of the lubrication film under applied pressure occurs step-wise, by expulsion of individual layers. These layering transitions appear to be thermally activated PT ; Persson2000 ; zilberman2001x1 ; zilberman2001x2 . Under strong confinement conditions, some lubricant fluids become solid-like Israel3a ; Israel3b ; Gao97a ; Gao97b ; Demi96a ; Demi96b ; Kleina ; Kleinb . Other fluids, notably water KleinW ; GranickW , remain liquid-like up to the last layer that can be removed upon squeezing. This is related to the expansion of water upon freezing.Jagla ; PerJCP , and should also hold for other liquids which expand upon freezing.

The phenomenology of layering transitions in 2D-solid-like boundary lubrication has been studied in Ref. PerJCP ; Perssona ; Perssonb . It has been shown in a series of computer simulations that for solid-like layers, the layering transitions are sometimes initiated by a disordering transition, after which the lubricant behaves in a liquid-like manner for the rest of the squeeze-out process. Since the typical lateral extension in surface force apparatus (SFA) experiments is of order 10-100 $`\mu `$m (much greater than atomic dimensions), it is reasonable to expect that during the layering transition the squeeze-out can be described in the framework of 2D continuum fluid mechanics.

For the first time such squeeze-out layering transitions were quite recently observed in a chain alcohol, C11H23OH (Ref. Mugele ; Mugelesubm ), by imaging the lateral variation of the gap between the two anvil surfaces as a function of time. These experiments addressed the $`n=10`$ transition. More recently, in a refined experimental setup, Mugele et al. MugelToBe  were able to image several layering transitions ($`nn1`$, $`n=5,4,3`$) of the silicon oil OMCTS (spherical molecule, diameter $`10`$ Å) in great detail.

The basic theory of 2D squeeze-out dynamics was outlined in Ref. PT . Initially the system is trapped in a meta-stable state at the initial film thickness. Squeeze-out starts by a thermally activated nucleation process in which a density fluctuation forms a small hole, of critical radius $`R_c10`$ Å. Once formed, a 2D pressure difference $`\mathrm{\Delta }p`$ develops between the boundary line separating the squeezed out region from the rest of the system, and the outer (roughly circular) boundary line of the contact area, thus driving out the rest of the 2D fluid. The origin of $`\mathrm{\Delta }p`$ is the elastic relaxation of the confining solids at the boundary line as is illustrated in Fig. 1.

All earlier analytical studies of squeeze out have assumed that the lubricant behaves as an incompressible 2D liquid. While that assumption is quite good for many practical situations, recent computer simulationsPerJCP  have shown that, at least at high squeezing pressures, strong density fluctuations may occur in the lubrication film. For example, in Ref. PerJCP  it was found that during the layering transition $`n=21`$, while islands of (temporarily) trapped bilayer ($`n=2`$) were removed by being squeezed into the monolayer, the density of the monolayer film was much higher in the region close to the trapped $`n=2`$ islands than further away. This resulted in a 2D pressure gradient in the film which induced a flow of the lubricant molecules away from the trapped islands. This kind of situations clearly calls for a consideration of the finite compressibility of the film.

In this paper we study the dynamics of the expulsion of the last liquid monolayer of molecules confined between two surfaces by solving the two-dimensional (2D) Navier-Stokes equation for a compressible liquid monolayer. We find that the squeeze-out is characterized by the parameter $`g_0P_0/\rho c^2`$, where $`P_0`$ is the average perpendicular (squeezing) pressure, $`\rho `$ the liquid (3D) density and $`c`$ the longitudinal sound velocity in the monolayer film. When $`g_01`$ the 2D liquid can be considered as incompressible, in which case the results of the earlier treatment are reproduced. We present numerical results for the squeeze-out time, and for the time-dependence of the radius of the squeezed-out region, for several values of the parameter $`g_0`$. The main changes due to compressibility occur right at the onset of the squeeze-out process, and just before its completion.

## II Theory

We assume the lubricant film to be in a 2D-liquid-like state, and the squeeze-out to be described by the 2D Navier-Stokes equations. For a compressible 2D liquid these equations take the form

$$\frac{n}{t}+(n𝐯)=0$$
$`(1)`$
$$\frac{𝐯}{t}+𝐯𝐯=\frac{1}{mn}p$$

$$+\nu _1^2𝐯+\nu _2𝐯\overline{\eta }𝐯$$
$`(2)`$
where $`mn`$ and $`𝐯`$ are the local 2D mass density and the velocity of the fluid, $`p`$ is the 2D-pressure, $`\nu _1`$ and $`\nu _2`$ are viscosities. The last term in (2), i. e. $`\overline{\eta }𝐯`$, describes the drag force acting on the fluid as it slides relatively to the solid walls.slidingfriction  The magnitude of the friction also depends on the nature of the solid walls, e. g. on the amplitude of the atomic corrugation and on the structure of the solid walls (amorphous vs. crystalline, commensurate vs. incommensurate). As long as the lubricant can be considered as a 2D fluid, as in the recent squeeze out experiments by Mugele et al., all the details of the solid walls are properly taken into account by the friction coefficient $`\overline{\eta }`$. In principle $`\overline{\eta }`$ depends on the normal pressure too, but in the measurements of Mugele et al. this dependence is slight enough that it appears to be negligiblezilbermannew .

Simple dimensional arguments (see Ref. zilberman2001x2 ) show that one can usually neglect the nonlinear and the viscosity terms in (2), and that one can also assume the velocity field to change so slowly that the time derivative term can be neglected too. Thus,

$$p+mn\overline{\eta }𝐯=0.$$
$`(3)`$
In what follows we will assume that the squeeze-out nucleated in the center ($`r=0`$) of the contact area, and spread circularly towards the periphery. The 2D-pressure $`p`$ at the outer boundary $`r=R`$ of the contact area takes the constant value $`p_0`$ (the spreading pressure) while it takes a higher value $`p_1`$ at the inner boundary towards the $`n=0`$ area (Fig. 1). In fact, if $`P(r)`$ is the perpendicular pressure acting in the contact area, then $`p_1(r)=p_0+P(r)a`$, where $`a`$ is the thickness of the monolayer (see Ref. PT ). We may in most practical applications assume a Hertz contact pressure

$$P(r)=\frac{3}{2}P_0\left(1\frac{r^2}{R^2}\right)^{1/2}$$
$`(4)`$
Finally, in order to have a complete set of equations we must specify the relation between the 2D pressure $`p`$ and the 2D density $`n`$. We will assume that

$$p=p_0+mc^2(nn_0)$$
$`(5)`$
where the compressibility is $`B=1/mc^2`$, $`c`$ being the longitudinal 2D sound velocity. Here $`n_0`$ is the 2D lubricant density at the periphery $`r=R`$ of the contact area.

Writing the velocity field as $`𝐯=\widehat{r}v(t)`$, equations (1) and (3) takes the form

$$\frac{n}{t}+\left(\frac{}{r}+\frac{1}{r}\right)(nv)=0$$
$`(6)`$
$$\frac{p}{r}+mn\overline{\eta }v=0$$

Using (5) the last equation takes the form

$$\frac{n}{r}+\frac{\overline{\eta }}{c^2}(nv)=0$$
$`(7)`$
Combining (6) and (7) gives

$$\frac{n}{t}\frac{c^2}{\overline{\eta }}\left(\frac{}{r}+\frac{1}{r}\right)\frac{n}{r}=0$$
$`(8)`$
The density $`n`$ satisfies the boundary conditions

$$n(R,t)=n_0$$
$`(9)`$
$$n(r_1(t),t)=n_0+\frac{a}{mc^2}P(r_1(t))$$
$`(10)`$
where $`r=r_1(t)`$ is the equation for the inner boundary line. Finally, we note that the velocity $`v(r_1(t),t)`$ of the 2D-liquid at the inner boundary $`r=r_1(t)`$ must equal the radial velocity $`\dot{r}_1(t)`$ of the boundary line. Thus, if we put $`r=r_1(t)`$ and $`v(r_1(t),t)=\dot{r}_1(t)`$ in (7) we get

$$\frac{n}{r}(r_1(t),t)=\frac{\overline{\eta }}{c^2}n(r_1(t),t)\dot{r}_1$$
$`(11)`$
Let us at this point introduce dimensionless variables. If we measure the radius $`r`$ in units of $`R`$, time $`t`$ in units of $`\tau =\overline{\eta }R^2/c^2`$ and density $`n`$ in units of $`n_0`$ we get

$$\frac{n}{t}\left(\frac{}{r}+\frac{1}{r}\right)\frac{n}{r}=0$$
$`(12)`$
and the boundary conditions becomes

$$n(1,t)=1$$
$`(13)`$
$$n(r_1(t),t)=1+g(t)$$
$`(14)`$
and

$$\frac{n}{r}(r_1(t),t)=[1+g(t)]\dot{r}_1$$
$`(15)`$
where

$$g(t)=\frac{a}{mc^2n_0}P(r_1(t))$$
$`(16)`$
If we assume that $`P(r)`$ is of the Hertz form, then

$$g(t)=g_0\frac{3}{2}\left(1[r_1(t)]^2\right)^{1/2}$$
$`(17)`$
where

$$g_0=\frac{aP_0}{mc^2n_0}.$$
$`(18)`$
Thus, the theory depends only on a single parameter $`g_0`$. Note that $`n_0/a\rho `$, where $`\rho `$ is the 3D density of the liquid, and with the typical values $`\rho 1000\mathrm{kg}/\mathrm{m}^3`$ and $`c700\mathrm{m}/\mathrm{s}`$ we get $`\rho c^2500\mathrm{MPa}`$. In the experimental studies by Mugele and Salmeron (as well as in most other Surface Force Apparatus studies) the average squeezing pressure $`P_0500\mathrm{MPa}`$ which implies that the liquid can be considered as incompressible and the theory developed elsewhere can be usedPT ; zilberman2001x1 ; zilberman2001x2 . However, in many practical situations the pressure $`P_0`$ might be similar to the yield stress of the solids which for metals is typically of order $`1000\mathrm{MPa}`$. In these cases it is necessary to include the finite compressibility of the lubricant in order to accurately describe the squeeze-out dynamics. In Sec. III we shall present numerical results based on the above equations for a range of values of the parameter $`g_0`$.

Let us first calculate the squeeze-out time in the limit when the compressibility $`B=\mathrm{}`$. In this limit, all the adsorbates from the region $`r<r_1(t)`$ will be piled-up right at the boundary line $`r=r_1`$. If we consider a small angular section $`\mathrm{\Delta }\varphi `$ then the driving force acting on the boundary line is $`F=r_1\mathrm{\Delta }\varphi [p_1p_0]=r_1\mathrm{\Delta }\varphi P(r_1)a`$. This must balance the frictional drag force which is $`Nm\overline{\eta }\dot{r}_1`$ where the number of adsorbates $`N=\pi r_1^2(\mathrm{\Delta }\varphi /2\pi )n_0`$. Thus we get

$$\dot{r}_1r_1=\frac{2P(r)a}{n_0m\overline{\eta }}$$
$`(19)`$
Assume first that $`P=P_0`$ is constant. Thus, (19) gives

$$\dot{r}_1r_1=\frac{R^2}{2T}$$
$`(20)`$
where

$$T=\frac{mn_0\overline{\eta }R^2}{4P_0a}$$
$`(21)`$
is the squeeze out time for an incompressible 2D fluid, see Ref. PT . Integrating (20) gives

$$\frac{r_1(t)}{R}=\left(\frac{t}{T}\right)^{1/2}$$
$`(22)`$
so that the squeeze-out time, $`T^{}`$, in the limit $`B=\mathrm{}`$ is the same as for an incompressible 2D-fluid, $`T^{}=T`$. This suggests that the squeeze-out time is independent of the compressibility of the 2D-liquid, which our numerical simulations presented below indeed show to be the case.

Next, let us assume that the perpendicular pressure is of the Hertz form, appropriate for a fluid between curved surfaces. In that case (19) gives

$$\dot{r}_1r_1=\frac{R^2}{2T}\frac{3}{2}\left[1\left(\frac{r_1}{R}\right)^2\right]^{1/2}$$

or, with $`r_1^2/R^2=\xi `$,

$$_0^{r_1^2/R^2}\frac{d\xi }{(1\xi )^{1/2}}=\frac{3t}{2T}$$

Performing the integral gives

$$\left[1\left(\frac{r_1}{R}\right)^2\right]^{1/2}=1\frac{3t}{4T}$$

Thus, the squeeze-out time for the Hertzian contact pressure is $`T_\mathrm{H}^{}=4T/3`$, and the time dependence of $`r_1(t)`$ is given by

$$\frac{r_1}{R}=\left(\frac{3t}{2T}\right)^{1/2}\left(1\frac{3t}{8T}\right)^{1/2}$$

or

$$\frac{r_1}{R}=\left(\frac{2t}{T_\mathrm{H}^{}}\right)^{1/2}\left(1\frac{t}{2T_\mathrm{H}^{}}\right)^{1/2}$$

The squeeze-out time for an incompressible fluid with a Hertzian contact pressure is $`T_\mathrm{H}=(4T/3)(2\mathrm{ln4})0.8183T`$ which is a factor $`2\mathrm{ln4}0.6137`$ smaller than for the $`B=\mathrm{}`$ limiting case. Thus, for a Hertzian contact pressure the squeeze-out time does depend on the compressibility of the 2D lubrication film.

## III Numerical results

We consider first the case of a spatially constant squeezing pressure, $`P(r)=P_0`$, appropriate for a fluid between flat surfaces. In this case numerical calculations show that the squeeze-out time is independent of the compressibility, i.e., independent of the parameter $`g_0`$. However, the time dependence of the squeeze-out radius $`r_1(t)`$ does depend on $`g_0`$. In Fig. 2(a) we show this time dependence for four cases, $`g_0=0`$, $`0.2`$, $`0.5`$ and $`2`$. In Fig. 2(b) we show the adsorbate density profile at the time point when the squeeze-out radius $`r_1=0.3R`$, for the same four $`g_0`$-values as in (a).

Let us now assume the Hertzian squeezing pressure profile, $`P(r)=P_\mathrm{H}(r)`$. In this case the squeeze-out time depends on the compressibility, increasing from $`0.8183T`$ to $`1.3333T`$ as the compressibility increases from $`B=0`$ to $`\mathrm{}`$. In Fig. 3 we show the squeeze-out time as a function of $`g_0`$. In Fig. 4 we show the time dependence of $`r_1(t)`$ for four cases, $`g_0=0`$, $`0.2`$, $`0.5`$ and $`2.0`$. In Fig. 4(b) we show the adsorbate density profile at the time point when the squeeze-out radius $`r_1=0.3R`$, for the same four $`g_0`$-values as in (a). The results of Fig. 2 and Fig. 4 show that a compressibility parameter $`g_0`$ as small as 0.2 should produce a measurable difference in the squeeze-out evolution at all times. The main effect of compressibility appears at the beginning and at the end of the squeeze-out process. Initially, compressibility favors piling up of fluid at the squeeze-out boundary, which can as a result expand more rapidly, compared with the case of an incompressible fluid. On the other hand, when the hole approaches the boundary of the contact region, the squeezing-out speed of the compressed fluid is smaller, due to its increased density and, consequently, friction. In the case of uniform squeezing pressure these two effects compensate exactly, leading to a total squeeze-out time independent of the compressibility. For a Hertzian pressure distribution, the squeeze-out time instead increases with increasing compressibility: the initial speed up is overcompensated by the enhanced friction at the periphery of the contact area.

## IV Summary

The continuum mechanics theory of squeeze-out has been solved numerically for a 2D compressible liquid. We considered both a constant normal stress, and a Hertzian normal stress, and assumed a centro-symmetric squeeze-out. We found that the squeeze-out is completely characterized by the compressibility parameter $`g_0P_0/\rho c^2`$, where $`P_0`$ is the average perpendicular squeezing pressure, $`\rho `$ the (3D) liquid density and $`c`$ the longitudinal sound velocity in the monolayer film. When $`g_01`$ the 2D liquid can be considered incompressible, and the earlier results are reproduced. We presented numerical results for the squeeze-out time, and for the time dependence of the squeeze-out radius, for a grid of values of the parameter $`g_0`$. For a constant squeezing pressure, the squeeze-out time was found to be independent of the compressibility parameter $`g_0`$, while for a Hertzian contact pressure it increased slightly with increasing $`g_0`$. It is hoped that these theoretical results will soon be submitted to experimental check.

###### Acknowledgements.

B.P. thanks BMBF for a grant related to the German-Israeli Project Cooperation “Novel Tribological Strategies from the Nano-to Meso-Scales”, and ICTP/SISSA in Trieste where the main part of this work was carried out. Work at SISSA was sponsored by INFM PRA NANORUB, by MIUR through COFIN 2001, and by EU, contract ERBFMRXCT970155 (FULPROP).


text:
# Linear temperature dependence of conductivity in the “insulating” regime of dilute two-dimensional holes in GaAs

## Abstract

The conductivity of extremely high mobility dilute two-dimensional holes in GaAs changes linearly with temperature in the insulating side of the metal-insulator transition. Hopping conduction, characterized by an exponentially decreasing conductivity with decreasing temperature, is not observed when the conductivity is smaller than $`e^2/h`$. We suggest that strong interactions in a regime close to the Wigner crystallization must be playing a role in the unusual transport.

In recent years, there has been great interest in the transport of low density two-dimensional (2D) charge carrier systems. In contradiction to the the scaling theory of localizationloc  which predicted that all states in 2D in the absence of electron-electron interactions are localized, experimental observations of “metallic” behavior showing increasing conductivity ($`\sigma `$) with decreasing temperature ($`T`$) ($`d\sigma /dT<0`$) and an apparent metal-insulator transition (MIT) as the carrier density is lowered have been reported on many low density, low disorder 2D systems.mit  Although there have been many experimental and theoretical studies trying to understand this metallic behavior, there is still no consensus on its origin to date.

On the other hand, the “insulating” behavior ($`d\sigma /dT>0`$) for lower densities in these low disorder systems has not been studied as extensively as the metallic behavior. In highly disordered systems, the transport in the insulating regime is described by variable range hoppingmott ; es  (VRH) among the localized states, and the $`T`$-dependence of $`\sigma `$ is expressed as $`\sigma (T)=\sigma _0exp[(T_0/T)^x]`$, where $`\sigma _0`$ is a prefactor and $`T_0`$ is a characteristic temperature. The exponent $`x`$ depends on the density of states (DOS) at the Fermi energy. For a constant DOS in the absence of interactions, $`x=1/3`$ corresponding to the Mott VRH,mott  and for a DOS which has a Coulomb gap at the Fermi energy due to interactions, $`x=1/2`$ corresponding to the Efros-Shklovski VRH.es  In clean systems with low disorder, however, the origin of the “insulating” behavior ($`d\sigma /dT>0`$) in low density regime, where the interactions between carriers are strong, could be significantly different from that described by simple localization. In the ideal case, it was predicted that the ground state of the system in the low density limit is a Wigner crystal,wigner  which can be pinned by even a tiny amount of impurities in the system.

We have recently studied the MIT of extremely high mobility dilute 2D holes in GaAs and presented a detailed analysis for the $`T`$-dependence of $`\sigma `$ on the “metallic” side of the transition.noh  The metallic behavior persisted up to $`r_s=57`$, much higher than 37 predicted for Wigner crystallization ($`r_s`$ is the interaction energy to Fermi energy ratio, given by $`r_s=(p\pi )^{1/2}m^{}e^2/4\pi \mathrm{}^2ϵ`$, where $`p`$ is the hole density, $`m^{}`$ effective mass, and $`ϵ`$ the dielectric constant). With the persistence of metallic behavior up to such a large $`r_s`$ value, the question arises regarding the nature of the “insulating” behavior observed in this 2D system for $`r_s`$ reaching 80. To address this question, we present the $`T`$-dependence of $`\sigma `$ of the same sample on the “insulating” side of the MIT in this report. We show that the conductivity exhibits a linear dependence on temperature and does not follow the VRH behavior, and discuss its implications with regard to strong interaction effects close to or in the Wigner crystallization regime.

The sample used in this study is a heterojunction insulated-gate field-effect transistor (HIGFET) made on a (100) surface of GaAs.noh ; kane  A metallic gate, separated by an insulator (AlGaAs) from the semiconducting GaAs, is used to induce the 2D holes at the interface between the GaAs and AlGaAs. The mobility ($`\mu `$) of our sample measured at $`T=65`$ mK reaches $`1.8\times 10^6`$ cm2/Vs for $`p=3.2\times 10^{10}`$ cm-2, which is the highest achieved for 2D holes in this low density regime. This high mobility allows us to measure the temperature dependence of conductivity down to very low densities reaching $`p=1.5\times 10^9`$ cm-2, with $`r_s`$ near 80, and the MIT is observed at a critical density of $`p_c=3\times 10^9`$ cm-2, i.e. $`r_s=57`$. The experiment was done in a dilution refrigerator with a base temperature around 65 mK. The resistivity of the 2D holes were measured by a low-frequency lock-in technique with frequency as low as 0.1 Hz to maintain the out-of-phase signal less than 5 %. Excitation current of as small as 50 pA was used to ensure that the effect of Joule heating is negligible. In addition, the drive current was varied at base temperature to make sure that the linear response was measured.

Figure 1 shows the $`T`$-dependence of $`\sigma `$ for five 2D hole densities varying from $`5\times 10^9`$ cm-2 ($`r_s=44`$) to $`1.5\times 10^9`$ cm-2 ($`r_s=80`$). Data for higher densities deep inside the metallic region have been presented in Ref noh, . The top two traces are on the metallic side with $`\sigma >e^2/h`$ in the entire $`T`$ range. For $`p=5\times 10^9`$ cm-2, $`\sigma `$ decreases with increasing $`T`$ from 0.065 K to 0.15 K. This decrease of $`\sigma `$ with increasing $`T`$ ($`d\sigma /dT<0`$) is typical of the metallic behavior, the physical origin of which has been the focus of much recent research on 2D MIT. For $`T>0.15`$ K, $`\sigma `$ increases with increasing $`T`$ and the increase becomes linear for $`T>0.25`$ K. This linear dependence of $`\sigma `$ on $`T`$ continues for $`T`$ exceeding $`T_F`$, the Fermi temperature of the 2D hole gas (marked by the vertical bar in Fig. 1), which is 0.36 K for this density. Das Sarma and Hwangdas\_sarma  explained this change in $`T`$-dependence of $`\sigma `$ from $`d\sigma /dT<0`$ to $`d\sigma /dT>0`$ as a quantum-to-classical crossover. In the classical regime, when $`TT_F`$, scattering of the 2D carriers by charged impurities would give rise to linear increase of $`\sigma `$ with increasing $`T`$. We note that the linear dependence starts below $`T_F`$ in our data, and the crossover happens at a temperature about one third of $`T_F`$.

For $`p=4\times 10^9`$ cm-2 ($`T_F=0.29`$ K), similar $`T`$-dependence of $`\sigma `$ is observed, although the metallic behavior is much weaker, and the quantum-to-classical crossover occurs at lower $`T`$ ($`0.1`$ K) as expected. The linear dependence sets in at $`T=0.22`$ K. However, for $`p=3\times 10^9`$ cm-2 ($`r_s=57`$), the third trace from the top, $`d\sigma /dT>0`$ is observed in all our $`T`$ range. The metallic behavior ($`d\sigma /dT<0`$) is not seen in the data down to the lowest $`T`$, and for $`T<0.13`$ K $`\sigma `$ becomes smaller than $`e^2/h`$, suggestive of an insulator, although strictly linear dependence is not seen for $`T<0.1`$ K.

The bottom two traces show insulating behavior in that $`\sigma `$ is observed to decrease with decreasing $`T`$, albeit only linearly. For $`p=2.3\times 10^9`$ cm-2 ($`r_s=65`$), $`\sigma <e^2/h`$ for $`T<0.22`$ K and the linear dependence extends down to the base temperature. This linear dependence appears to be the same as that seen in the high hole density metallic side ($`5\times 10^9`$ cm-2 and $`4\times 10^9`$ cm-2), and the same linear dependence is seen for $`\sigma <e^2/h`$ as well as for $`\sigma >e^2/h`$. Further reducing the density down to $`1.5\times 10^9`$ cm-2 ($`r_s=80`$) does not change the $`T`$-dependence very much, but simply decreases the magnitude of $`\sigma `$, making the $`T`$ range where $`\sigma <e^2/h`$ wider.

In Fig. 2, we show the $`T`$-dependence of $`\sigma `$ for $`p=1.5\times 10^9`$ cm-2 in more detail. Clearly, $`\sigma `$ depends linearly on $`T`$ (as shown by the solid line) in the entire $`T`$ range, and $`\sigma <e^2/h`$ in almost the same range. When we tried to fit the data with the VRH model with exponent $`x`$ equal to either 1/3 or 1/2, none of them produced a good fit, as shown in the insets. This suggests that the transition to an “insulating” behavior as we lower the 2D hole density in our sample, determined from the sign of $`d\sigma /dT`$, is not a disorder induced transition. While one can argue that further decreasing the density will eventually lead to an exponentially localized behavior, strong interactions (with $`r_s=80`$) which must be present for this density make it unlikely that the transport can be described by single particle localization.

For a 2D electronic system with such large $`r_s`$ value, formation of a Wigner crystal should be considered. When the Wigner crystal is pinned by disorder, the transport is expected through thermal activation of the crystal over the pinning potential. In this case, $`\sigma `$ should decrease exponentially with decreasing $`T`$. Our data, on the contrary, shows a linear $`T`$-dependence in the $`\sigma <e^2/h`$ insulating regime. Here, we should note that reducing the carrier density not only increases $`r_s`$, but also decreases $`T_F`$, thus decreases the quantum-to-classical crossover temperature. In our data, for hole densities $`p=5\times 10^9`$ cm-2 and $`p=4\times 10^9`$ cm-2, the linear dependence of $`\sigma `$ on $`T`$ predicted for the $`TT_F`$ classical regime actually extends down to a $`T`$ even below $`T_F`$, suggesting that the crossover from classical to the degenerate regime seems to occur much below $`T_F`$. In any case, for a classical 2D electronic system, the relevant interaction parameter is $`\mathrm{\Gamma }`$, which is given by $`\mathrm{\Gamma }=e^2\sqrt{\pi n}/4\pi ϵk_BT`$, with $`n`$ being the carrier density. This parameter is equivalent to $`r_s`$ in the degenerate regime, with the kinetic energy replaced by $`k_BT`$. Studies of electrons on liquid Heliumgrimes ; mehrotra ; morf ; gann  have shown that transition to the Wigner crystal is around $`\mathrm{\Gamma }=127`$. In our sample, for the lowest hole density, $`p=1.5\times 10^9`$ cm-2, $`T_F=0.11`$ K. $`\mathrm{\Gamma }`$ is between 80 and 27 for $`T`$ from 0.11 K to 0.32 K, and equals 146 at $`T=0.06`$ K. Therefore, Wigner crystallization is not likely to occur in most of our $`T`$ range until $`T`$ reaches close to the base temperature. However, a more recent electron-on-helium experimentdjerfi  has reported residual signature of the Wigner crystal up to a much higher temperature with $`\mathrm{\Gamma }`$ down to 46. Thus, we should expect the 2D holes in our sample to be strongly correlated and the Wigner crystallization physics be important. The $`T`$-dependence arising from scattering of the 2D hole gas by charged impurities, as considered by Das Sarma and Hwang, even though remarkably similar to that observed in our experiment, cannot be the adequate explanation for the insulating regime.

We should also mention that localization corrections to the conductivity in the weakly interacting classical regime do not explain our data. Experimental study on weak localization in the classical regime has been performed for electrons on solid hydrogen surface.adams  In this case, the correction to the conductivity is proportional to $`T_F/T`$, and our data are not consistent with this expectation. The interaction correction is proportional to $`(T_F/T)^2`$. This term is negligible compared with that from the weak localization, and does not explain our data either.

Finally, the linear dependence of $`\sigma `$ on $`T`$ observed in the insulating regime of our data must result from strong interaction effects. In this respect, we think that the recent theory by Spivakspivak  may be relevant to our data. In his theory, the metallic state is a state in which a small concentration of Wigner crystal droplets are embedded in the Fermi liquid, and the insulating state corresponds to one in which a small fraction of Fermi liquid is present in a mostly crystalline phase. In the metallic side, the resistivity increases linearly in $`T`$ in the low temperature degenerate regime as the size of the Wigner crystal droplets grows. In the insulating side, he argues that the crystal forms a supersolid, which is not pinned by disorder, and $`\sigma `$ remains finite as $`T`$ goes to zero. At high temperatures, Wigner crystal droplets melt and the system behaves classically. In this case, the resistivity is inversely proportional to $`T`$ and $`\sigma `$ is linear in $`T`$. We note that this predicted $`T`$-dependence is identical to that by Das Sarma and Hwang, but it is based on strong correlation effects close to the Wigner crystallization. The prediction that $`\sigma `$ remains finite at $`T`$=0 in the insulating side is also significant, since in our data the linear dependence appears to continue below $`T_F`$ with finite offset.

In summary, we have studied the temperature dependence of the conductivity in the insulating side of the MIT for extremely high mobility dilute 2D hole system with $`r_s`$ reaching 80. A linear dependence of conductivity on temperature has been observed in the insulating side. We suggest that strong interaction effects in both degenerate and classical regime close to the Wigner crystallization should be considered to provide an understanding for this unusual transport.

This work has been supported by the NSF and the DOE at Princeton University.


text:
# Dark Energy and the Fate of the Universe

## I Introduction

Recent observations indicate that the Universe is accelerating, and it is spatially flat ($`\mathrm{\Omega }_{\mathrm{tot}}=\mathrm{\Omega }_M+\mathrm{\Omega }_D=1`$). Approximately $`0.3`$ of the total energy density of the universe $`\rho _010^{120}M_p^410^{29}`$ g/cm3 consists of ordinary matter ($`\mathrm{\Omega }_M0.3`$), and $`0.7`$ of the energy density corresponds to dark energy ($`\mathrm{\Omega }_D0.7`$) supernova ; Bond . There are three basic scenarios describing the evolution of the universe filled by dark energy.

i) From a purely phenomenological point of view, the simplest possibility is that the dark energy $`\rho _D`$ is represented by a positive vacuum energy (cosmological constant) $`\mathrm{\Lambda }0.7\rho _0`$. If this is the case, the universe will reach de Sitter (dS) regime and expand exponentially for an indefinitely long time, $`ae^{Ht}`$. This possibility seems much more natural than the other two possibility to be discussed below.

ii) It may also happen that the dark energy is the energy of a slowly changing scalar field $`\varphi `$ with equation of state $`p_D=w\rho _D`$, $`w1`$. In most of the models of dark energy it is assumed that the cosmological constant is equal to zero, and the potential energy $`V(\varphi )`$ of the scalar field driving the present stage of acceleration, slowly decreases and eventually vanishes as the field rolls to $`\varphi =\mathrm{}`$, see e.g. Dolgov:gh . In this case, after a transient dS-like stage, the speed of expansion of the universe decreases, and the universe reaches Minkowski regime.

iii) It is also possible that $`V(\varphi )`$ has a minimum at $`V(\varphi )<0`$, or that it does not have any minimum at all and the field $`\varphi `$ is free to fall to $`V(\varphi )=\mathrm{}`$. In this case the universe eventually collapses, even if it is flat Linde1986 . The simplest way to understand this unusual effect is to analyse the Friedmann equation $`\left(\frac{\dot{a}}{a}\right)^2=\rho /3`$ (in units $`M_p=1`$). The positive energy density of a normal matter, as well as the positive kinetic energy density of the scalar field, tend to decrease in an expanding universe. At some moment, the total energy density $`\rho `$, including the negative contribution $`V(\varphi )<0`$, vanishes. Once it happens, the universe, in accordance with the equation $`\left(\frac{\dot{a}}{a}\right)^2=\rho /3`$, stops expanding and enters the stage of irreversible collapse.

The last possibility for a while did not attract much attention. There was no specific reason to expect that the present regime of acceleration is going to end, and there was even less reason to believe that the universe is going to collapse any time soon.

Unfortunately, despite many attempts, we were unable to obtain a good theoretical description of the models i) and ii) in the context of M-theory. We will discuss this issue in Sect. 2. Meanwhile, in Kallosh:2002gf ; Kallosh:2002gg  it was found that one can describe the present state of acceleration of the universe in a broad class of models based on N=8 extended supergravity (the theory closely related to M-theory). However, the universe described by these models typically collapses within the time $`t_{\mathrm{collapse}}`$ comparable to its present age $`t_014`$ billion years.

In the beginning, this seemed to be a model-specific result that should be taken seriously only if one can construct fully realistic models of elementary particles based on extended supergravity. However, this result is valid not only in for the models based on extended supergravity but for many phenomenological models of dark energy based on the minimal N=1 supergravity Kallosh:2002gf .

In this paper we will argue that this result is even more general. It can be obtained in almost every model of dark energy, either based on supergravity or not, if one takes into account the possibility that the effective potential of the field $`\varphi `$ may have a minimum at $`V(\varphi )<0`$ or may be unbounded from below.

We will first review the possibilities to describe the accelerating universe starting with M/string theory and supergravity. Then we will describe our general argument.

## II Dark energy in M-theory and extended supergravity

### II.1 M-theory with compactification

The standard approach to the description of our world in M/string-theory is based on the assumption that our space-time is 11 or 10 dimensional, but 7 or 6 of these dimensions are compactified. The finitness of the volume of the compactified space is required so that the original D=11 or D=10 theory can be related to the resulting D=4 theory. One can describe the scale of compactification for example, in string theory by introducing scalar moduli field $`\varphi `$, which appears as a coefficient $`e^{\sqrt{6}\varphi }`$ in front of the potential energy $`V`$ in 4D. (We use the units $`M_p=1`$, where $`M_p^2=\frac{1}{8\pi G}`$.) Here $`V`$ is the value of the potential at which other scalar fields of the higher-dimensional theory are stabilized. However, the mechanism of stabilization of the compactified space is still unknown. In the absence of such mechanism, the term $`e^{\sqrt{6}\varphi }V`$ leads to the runaway behaviour $`\varphi \mathrm{}`$, which implies decompactification. During this process, the energy density $`e^{\sqrt{6}\varphi }V`$ falls down very quickly. A similar result is valid if we consider D=10 string theory as a result of compactification of D=11 supergravity. The bottom line is: during the cosmological evolution in D=4 the runaway moduli represent decompactification of all internal dimensions of M-theory.

The scale factor of the universe $`a(t)`$ in the theories with exponential potentials of the type of $`e^{\lambda \varphi }V`$ grows as $`t^{2/\lambda ^2}`$. In order to describe acceleration of the universe in such theories one would need to have $`\lambda <\sqrt{2}`$. In the theory with the potential $`e^{\sqrt{6}\varphi }V`$ the universe can only decelerate.

In application to string cosmology these observations imply that until we learn how to stabilize the compactified space, we cannot describe the accelerating universe approaching dS regime, as well as the accelerating universe with the energy density slowly approaching zero.

### II.2 M-theory with non-compactification

One way to avoid the problems discussed above is to consider the models with a non-compact internal space Kallosh:2001gr ; Kallosh:2002wj . One may start with D=11 or D=10 supergravity with internal space with an infinite volume and relate it to N=8 supergravity in D=4 Hull:1988jw . This is called ‘non-compactification.’ In this approach the connection between the original D=10 or D=11 theory and our D=4 world is more complicated than in the usual case of dimensional reduction. However, one can study these theories directly in D=4. These theories are interesting because they have maximal amount of supersymmetry, related to D=11 and/or D=10 supergravities with non-compact internal spaces.

The number of such models successfully describing an accelerated universe is very limited, due to the maximal amount of supersymmetries. Some of these theories have dS solutions and can describe dark energy Kallosh:2002gf ; Kallosh:2002wj . These dS solutions correspond to the extrema of the effective potentials $`V(\varphi )`$ for some scalar fields $`\varphi `$. An interesting and very unusual feature of these scalars in all known theories with $`N2`$ is that their mass squared is quantized in units of the Hubble constant $`H_0`$ corresponding to dS solutions: $`\frac{m^2}{H_0^2}=n`$, where $`n`$ are some integers of the order 1. This property was first observed in Kallosh:2001gr  for a large class of extended supergravities with unstable dS vacua, and confirmed and discussed in detail more recently in Kallosh:2002wj  with respect to a new class of $`N=2`$ gauged supergravities with stable dS vacua Fre:2002pd .

The meaning of this result can be explained in the following way. The simplest potential for a scalar field $`\varphi `$ in this theory has the form

$$V(\varphi )=\mathrm{\Lambda }\left(2\mathrm{cosh}\sqrt{2}\varphi \right).$$
(1)
Usually the potential near its extremum can be represented as $`V(\varphi )=\mathrm{\Lambda }+m^2\varphi ^2/2`$, where $`\mathrm{\Lambda }`$ and $`m^2=V^{\prime \prime }(0)`$ are two free independent parameters. However, in extended supergravities with $`\mathrm{\Lambda }>0`$ one always has $`m^2=V^{\prime \prime }(0)=nV(0)/3=n\mathrm{\Lambda }/3`$, where $`n`$ are integers (we are using units $`M_p=1`$) Kallosh:2001gr ; Kallosh:2002wj . Taking into account that in dS space the Hubble constant is given by $`H_0^2=\mathrm{\Lambda }/3`$, one has, for $`|\varphi |1`$, $`V(\varphi )=\mathrm{\Lambda }(1+n\varphi ^2/6)=3H_0^2(1+n\varphi ^2/6)`$. In particular, in all known versions of $`N=8`$ supergravity dS vacuum corresponds to an unstable maximum, $`m^2=6H^2`$ Kallosh:2001gr ; Kallosh:2002wj , i.e. at $`|\varphi |1`$ one has

$$V(\varphi )=\mathrm{\Lambda }(1\varphi ^2)=3H^2(1\varphi ^2).$$
(2)
One can easily verify that the simplest potential (1) satisfies this rule. The main property of this potential is that $`m^2=V^{\prime \prime }(0)=2V(0)=2\mathrm{\Lambda }=6H_0^2`$. One can show that a homogeneous field $`\varphi 1`$ with $`m^2=6H_0^2`$ in the universe with the Hubble constant $`H_0`$ grows as follows: $`\varphi (t)=\varphi _0\mathrm{exp}cH_0t`$, where $`c=(\sqrt{33}3)/21.4`$. Consequently, in the universe with the energy density dominated by $`V(\varphi )`$ it takes time $`t0.7H_0^1\mathrm{ln}\varphi _0^1`$ until the scalar field rolls down from $`\varphi _0`$ to the region $`\varphi 1`$, where $`V(\varphi )`$ becomes negative. Once it happens, the universe rapidly collapses.

Note, that the present age of the universe $`t_0`$ is approximately equal to $`H_0^1`$, and the total time of the development of the instability leading to the global collapse of the universe is given by $`t0.7H_0^1\mathrm{ln}\varphi _0^1`$, which is also of the same order as $`H_0`$, unless $`\varphi _0`$ is exponentially small. This explains the main result of Ref. Kallosh:2002gf ; Kallosh:2002gg : the universe described by this class of theories is going to collapse within the time $`t_{\mathrm{collapse}}`$ comparable to its present age $`t_014`$ billion years, see Fig. 1.

### II.3 Two-field models of dark energy in N=8 supergravity

There are also 3 different cosmological models of dark energy based on N=8 supergravity, with the following potentials of two fields, $`\varphi `$ and $`\sigma `$. We will list them here, leaving a detailed description of their cosmological implications for the future publication Shmakova

$$V(\varphi ,\sigma )=\mathrm{\Lambda }e^{\sqrt{\frac{2}{3}}\sigma }\left(3\mathrm{cosh}(2\sqrt{2}\varphi )\right).$$
(3)
This is an improved version of Townsend’s model of M-theory quintessence Townsend:2001ea . Townsend used this potential at $`\varphi =0`$ and found the exponential potential $`e^{\sqrt{\frac{2}{3}}\sigma }`$ which may describe the current acceleration with $`a(t)t^3`$. The complete form of the potential shows that the full potential of this theory is unbounded from below and unstable with respect to the generation of the field $`s`$. The speed of the development of the instability is determined by the curvature of the potential in the $`\varphi `$-direction, $`V^{\prime \prime }=4V(0)`$, i.e. $`m_\varphi ^212H_0^2`$. This allows for the existence of a stage of accelerated expansion, but eventually the universe collapses, just like in the theory (1) discussed above.

The second model is

$$V(\varphi )=\mathrm{\Lambda }e^{\sqrt{\frac{2}{7}}\sigma }\left(24e^{\sqrt{\frac{1}{21}}\varphi }8e^{2\sqrt{\frac{3}{7}}\varphi }3e^{8\sqrt{\frac{1}{21}}\varphi }\right).$$
(4)
The $`\varphi `$ part of the potential has a dS maximum at $`\varphi =\sqrt{\frac{3}{7}}\mathrm{log}2`$. It may describe the current acceleration during the slow growth of the field $`t`$: $`a(t)t^7`$. The instability with respect to the field $`\varphi `$ eventually develops and the universe collapses.

The last model of this type has a saddle point dS solution at $`\varphi =0,\sigma =0`$

$$V(\varphi )=\mathrm{\Lambda }\left(2\mathrm{cosh}\sigma \mathrm{cosh}\sqrt{2}\varphi \right).$$
(5)
At the point $`\sigma =0`$ this model reduces to the model (1) which describes dark energy and eventual collapse of the universe.

A lot of work should be done to incorporate usual matter fields and construct realistic cosmology in N=8 supergravity. However, it is quite encouraging that there is a class of models based on N=8 supergravity which can describe the present stage of acceleration of the universe. All of these models share the same property: at some moment expansion of the universe stops and the universe collapses within time $`t_{\mathrm{collapse}}`$ comparable to its present age $`t_014`$ billion years.

## III A general class of dark energy models

The main reason for the coincidence of the two different time scales, the present age of the universe and the time until the Big Crunch, is the relation $`|m^2|H^2`$. But this relation appears not only in the extended supergravity. It is often valid for the moduli fields in $`N=1`$ supergravity Dine:1983ys . In particular, it is valid in the simplest Polónyi-type toy model for dark energy in $`N=1`$ supergravity Kallosh:2002gf . If one does not fine-tune the value of the cosmological constant in this model to be equal to zero, one has two equally compelling options: $`V(\varphi )>0`$ and $`V(\varphi )<0`$ in the minimum of the effective potential. In the first case, the universe enters dS regime of eternal expansion. In the second case, the universe collapses. The time until the global collapse depends on the parameters of the model and the initial conditions, but typically it is of the same order as $`t_014`$ billion years, just as in the extended supergravity Kallosh:2002gf .

Another interesting model is the axion quintessence Frieman:1995pm ; Choi:1999xn . In the M-theory motivated version of this model proposed in Choi:1999xn  one has $`V(\varphi )\mathrm{\Lambda }(\mathrm{cos}(\varphi /f)+C)`$, where the value of constant $`C`$ depends on the details of the model. For $`C=0`$, $`f=O(M_p)`$ one finds $`m^2=V^{\prime \prime }(0)=O(H_0^2)`$. According to Kallosh:2002gf ; Jan , this version of the axion quintessence model can successfully describe the present stage of acceleration of the universe, but, just like the $`N=8`$ models, it leads to a global collapse of the universe in the future within the typical time $`t_{\mathrm{collapse}}t_014`$ billion years.

In fact, the crucial relation $`|m^2|H^2`$ is valid in most of the dark energy models. It is the standard inflationary slow-roll condition, which should be valid at the present stage of the late-time inflation/acceleration of the universe. Sometimes this slow-roll condition can be violated Kallosh:2001gr ; Linde:2001ae , but typically one cannot have a prolonged stage of acceleration of the universe for $`|m^2|H^2`$. As a result, one can give a simple argument suggesting that the coincidence of the two different time scales, the present age of the universe and the time remaining until the Big Crunch, is a generic property of many models of dark energy.

Let us make the simplest assumption that the expansion of the universe at $`t>t_0`$ can be approximately described by the simple power-law equation $`a(t)a(t_0)\left(\frac{t+c}{t_0+c}\right)^r,`$ where $`c`$ is some constant. The Hubble constant in the universe with $`a(t)(t+c)^r`$ is given by $`\frac{r}{t+c}`$, which means that the total energy density is equal to $`3H^2=\frac{3r^2}{(t+c)^2}`$. The parameter $`c`$ is determined by the continuity condition for $`H`$: $`\frac{r}{t_0+c}=H_0t_0^1.`$ The last part of this equation is the observed relation between the present value of $`H`$ and the age of the universe in the simplest $`\mathrm{\Lambda }`$CDM model. This gives $`c(r1)t_0`$.

The energy density will become $`N^2`$ times smaller than now at the time $`tt_0(N1)rt_0`$. In particular, the energy density will become 2 times smaller at the time $`tt_00.4rt_0`$ from now, and it will become 9 times smaller at $`tt_0=2rt_0`$. How large is this time interval?

Acceleration of the universe implies that $`r>1`$. Taking for definiteness $`r=2`$, one finds that the density of the universe will decrease 2 times at the moment $`tt_00.8t_010`$ billion years from now, and it will drop 9 times $`50`$ billion years from now.

Now let us take any model of dark energy with $`r=2`$ and add to the scalar potential $`V(\varphi )`$ a tiny negative cosmological constant $`\mathrm{\Lambda }=\rho _0/N^2`$ with $`N^21`$. The evolution of the universe up to the present moment will not change significantly. One can make this model even better by multiplying $`V(\varphi )`$ by some factor $`B>1`$ to ensure that the present value of $`\mathrm{\Omega }_D`$ remains equal $`0.7`$ even after we add a negative constant $`\mathrm{\Lambda }=\rho _0/N^2`$ to $`V(\varphi )`$. This model will remain a viable model of dark energy. However, after some time $`tt_02(N1)t_0`$ the value of $`V(\varphi )`$ will drop down more than N times, the total energy density of the universe will vanish, the universe will stop expanding, and soon after that it will collapse. If $`N`$ is not too large, the universe collapses at the time $`t_{\mathrm{collapse}}=O(t_0)`$.

One could think that the argument given above works only for the marginal situation, when $`a(t)t^r`$ with $`r=O(1)`$. However, the final result is very general. Our procedure of modification of the potential (subtraction of a constant and a subsequent compensation of the decrease of $`\mathrm{\Omega }_D`$ via the multiplication of $`V(\varphi )`$ by a constant $`B>1`$) works until the moment when the theory no longer represents dark energy because the potential becomes too steep. This happens for any theory of dark energy, even if the original potential was extremely flat (which corresponds to $`r1`$). Once it happens, the field rolls down within the time $`O(H^1)`$. The rolling field rapidly approaches the region of negative $`V(\varphi )`$, which typically leads to the collapse of the universe within the time $`O(H^1)`$.

This suggests that many models of dark energy considered in the literature have viable counterparts that can be obtained from the original models by adding a negative cosmological constant $`\mathrm{\Lambda }\rho _0`$. Unless the absolute value of this extra term is many orders of magnitude smaller than $`\rho _0`$, these models will describe the universe collapsing within the time comparable to the present age of the universe.

## IV Example: Dark energy with an exponential potential

As an example illustrating the general argument given in the previous section, let us consider dark energy described by the scalar field $`\varphi `$ with an exponential potential

$$V(\varphi )=\mathrm{\Lambda }e^{\lambda \varphi }.$$
(6)
We already mentioned that this theory describes an accelerating universe for $`\lambda <\sqrt{2}`$. However, in the realistic situation this condition can be slightly relaxed. One should take into account that the post-inflationary universe rapidly expanded, being dominated by hot matter and then by cold dark matter. This expansion does not allow the field $`\varphi `$ to move until the energy density of matter becomes sufficiently small. Therefore in the beginning the kinetic energy of the scalar field is very small, so it has the same equation of state as the cosmological constant. Then the field starts moving slowly, loosing its energy at a much slower paste than CDM. Eventually, the universe enters the stage with $`\mathrm{\Omega }_D=0.7`$ (the present time), and then $`\mathrm{\Omega }_D`$ continues to grow. As a result, one can have $`\mathrm{\Omega }_D=0.7`$ and $`w<0.6`$ at present for $`\lambda 1.7`$ LopesFranca:2002ek ; Kallosh:2002gf .

In our investigation we will assume, without loosing the generality, that the initial value of the field was $`\varphi =0`$ (one can always rescale the field and the potential). Then one should find such value of the parameter $`\mathrm{\Lambda }`$ that the universe enters the stage with $`\mathrm{\Omega }_D=0.7`$ at the same time when its Hubble constant acquires its present value $`H_010^{60}M_p`$. This requires fine-tuning, but this is the same fine-tuning that hampers all models of dark energy LopesFranca:2002ek ; Kallosh:2002gf .

Now we will consider a class of potentials of a more general type, containing a constant negative contribution, as suggested in the previous section:

$$V(\varphi )=\mathrm{\Lambda }_C\left(e^{\lambda \varphi }C\right).$$
(7)
Here $`C`$ is some positive constant. The constant $`\mathrm{\Lambda }_c`$ should be found anew for each new value of $`C`$, just as we did for the case $`C=0`$ described above. This can be done, and the cosmological evolution of this model can be easily studied using the methods described in Kallosh:2002gf . Here we will only present the results of our investigation of a model with $`\lambda =1`$.

Figure 3 shows the evolution of the scale factor of the universe for 4 different values of parameter $`C`$: $`C=0,\mathrm{\hspace{0.17em}0.1},\mathrm{\hspace{0.17em}0.2},\mathrm{\hspace{0.17em}0.3}`$ and $`0.4`$. For all of these cases one can find such parameters $`\mathrm{\Lambda }_C`$ that the value of the Hubble constant at $`t=t_0`$ coincides with its present value (the curves $`a(t)`$ have the same derivative at $`t=0`$ in Figure 3).

All of the models with $`C=0,\mathrm{\hspace{0.17em}0.1},\mathrm{\hspace{0.17em}0.2},\mathrm{\hspace{0.17em}0.3}`$ with the potentials shown in Figure 2 can represent dark energy with $`\mathrm{\Omega }_D=0.7`$ at the present moment (the point $`z=0`$ in Figure 4). However, the largest value of $`\mathrm{\Omega }_D`$ in the case $`C=0.4`$ is $`0.58`$, and the equation of state $`w`$ in this case blows up near $`z=0`$, see Figure 5.

As we see, the universe with $`C=0`$ (thick red line in Figures 2-6) always continues its accelerated expansion. However, in all other cases the universe collapses within the time ranging from $`1.3t_018`$ billion years (for $`C=0.4`$) to $`7.5t_0105`$ billion years (for $`C=0.1`$), in agreement with the argument given in the previous section. The universe with $`C=0.4`$ never accelerates, so this model is ruled out by the existing observations.

Thus, for every model $`V(\varphi )=\mathrm{\Lambda }e^\varphi `$ successfully describing dark energy there exist many other models $`V(\varphi )=\mathrm{\Lambda }_C\left(e^\varphi C\right)`$ with $`0<C0.3`$ which provide an equally good description of the present stage of acceleration, but lead to the global collapse of the universe within the next $`10^{10}10^{11}`$ years.

Similar results can be obtained for many other models of dark energy, e.g. for the theories $`V(\varphi )=\mathrm{\Lambda }_C\left(e^{\lambda \varphi }C\right)`$ with all possible values of $`\lambda `$, or for the models with the inverse power law potentials.

## V Discussion

Recent discovery of acceleration of the universe is one of the major challenges for the modern theory of fundamental interactions. After many attempts to explain why the cosmological constant must be zero, the theorists switched to the new paradigm and started trying to explain why it should be positive and why, consequently, the universe should expand forever. The first attempts to do so in the context of M/string theory and extended supergravity revealed many problems described, e.g., in Gibbons:85 ; Hellerman:2001yi . Then we learned that one can describe acceleration of the universe in extended supergravity, but in all models based on extended supergravity, with the exception of the N=2 model of Fre:2002pd , the regime of acceleration is unstable Kallosh:2001gr . Typically it ends by a global collapse of the universe within the time comparable with the present age of the universe, $`t_{\mathrm{collapse}}1030`$ billion years Kallosh:2002gf ; Kallosh:2002gg .

In this paper we have shown that the possibility of a global collapse is not specific to supergravity but is, in fact, quite generic. For every model of dark energy describing eternally expanding universe one can construct many closely related models which describe the present stage of acceleration of the universe followed by its global collapse.

This does not mean that we are making a doomsday prediction. None of the existing theoretical models of dark energy look particularly natural and attractive. It may happen that eventually we will find good theoretical models describing an eternally accelerating universe, or conclude, on the basis of anthropic considerations, that dark energy should change in time extremely slowly, so that the collapse will occur exponentially far away in the future Vil2002 . We hope to return to this question in the future publications. However, in the absence of a compelling theory of dark energy one may also consider a more humble approach and try to compare predictions of various models of dark energy with observations, see e.g. Copeland .

If the universe is going to collapse, then in the beginning of this process the speed of expansion of the universe should gradually slow down, see Fig. 6. This is accompanied by the rapid growth of the parameter $`w`$ at small $`z`$, as shown in Fig. 5. We find it quite significant that some of the models predicting global collapse can be already ruled out by the existing observational data. For example, all models $`\mathrm{\Lambda }_C\left(e^\varphi C\right)`$ predicting the global collapse within the next 18 billion years, do not describe the present stage of acceleration, and therefore contradict the recent cosmological observations. Thus, even though the observational data cannot rule out the general possibility of the global collapse in the distant future, they can help us to put strong constraints on the time of the possible Big Crunch.

It is a pleasure to thank R. Bond, L. Kofman, J. Kratochvil, E. Linder, S. Prokushkin and M. Shmakova, for useful discussions. This work was supported by NSF grant PHY-9870115. The work by A.L. was also supported by the Templeton Foundation grant No. 938-COS273.


text:
# Subharmonics and Aperiodicity in Hysteresis Loops

## Abstract

We show that it is possible to have hysteretic behavior for magnets that does not form simple closed loops in steady state, but must cycle multiple times before returning to its initial state. We show this by studying the zero-temperature dynamics of the 3d Edwards Anderson spin glass. The specific multiple varies from system to system and is often quite large and increases with system size. The last result suggests that the magnetization could be aperiodic in the large system limit for some realizations of randomness. It should be possible to observe this phenomena in low-temperature experiments.

Hysteresis in magnetic systemsmasters  is the prototype for understanding history dependent behavior in all complex systems.colloids ; cdw ; alloys  Although the usual example of this is for ferromagnets, it can exist in other magnetic systems such as anti-ferromagnetsafm  or spin glasses.Zimanyi ; Helmut

For magnetic systems, when one cycles adiabatically between a large positive and negative field $`H`$, the magnetization $`M`$ forms a closed loop in the $`HM`$ plane because of saturation. There is an implicit expectation that even if the field were cycled repeatedly between two moderate values, the magnetization would settle down to a steady state closed loop, after some initial transient response. In this letter, we show that this expectation is often wrong. Many cycles of the external field are necessary to bring the magnetization back to its initial value, and that as the system size grows, so does the number of cycles needed.

We do this for Ising spin glass models: both short range and long range models violate the usual assumption of simple closed loops. Even in steady state, the magnetization returns to its initial value only after multiple cycles of the field. Another way of looking at this phenomenon, probably more relevant to real experiments, is that with a steady state oscillatory field the magnetization shows a subharmonic component. The number of cycles it takes for the magnetization to return varies from system to system, depending on the realization of the randomness, and can be very large. Our estimates suggest that this effect should be observable in nano-particles of spin glass material at very low temperatures.

This result might appear to be rather surprising as there are many cases when the strong result of “return point memory”Sethna ; RPM  can be shown: if a system starts at a field $`H_1`$ and ends at a field $`H_2`$, its final state is independent of the time dependence of H as long as H is never outside the interval $`[H_1,H_2]`$. In the case above this would imply at most one transient cycle before the steady state closed loop was reached. An example of such a system is a random field Ising ferromagnet where an elegant proof of return point memory was given.Sethna  However the proof does not apply to cases where antiferromagnetic interactions are also present, such as a spin glass. We show below that in this case there is a much richer class of memory effects.

Earlier work has reported subharmonic response in the magnetization of ferromagnetic materials subjected to oscillating magnetic fields. engineers  It is a specific example of the general phenomenon of ferroresonance. However, this occurs when the magnetic material is part of an electrical circuit (e.g. as the core of an inductor) that is driven at finite frequency. Furthermore the intrinsic response of the ferromagnetic material is expected to show return point memory. From a dynamical systems viewpoint, this is simply the observation that a driven nonlinear oscillator can respond at a frequency different from the driving frequency, possibly depending on the initial conditions of the system. There has also been theoretical work on finite frequency hysteresis loops in ferromagnets.Nattermann ; Dhar1 ; Sides  In contrast, we only consider $`H`$ fields that vary sufficiently slowly to be adiabatic and our non-return point effect is intrinsic to the material.

We consider the standard Edwards-Anderson spin glass Hamiltoniansg1 ; sg2  with Ising spins and nearest neighbor interactions in three dimensions,

$$=\underset{i,j}{}J_{i,j}S_iS_jH\underset{i}{}S_i.$$
(1)
The coupling strengths $`J_{i,j}`$ are chosen to be uniform random numbers between $`1`$ and $`1`$. The Ising spin variables $`S_i`$ take values $`\pm 1`$. Because we expect to apply this model to granular materials, we employ free boundary conditions. Other boundary conditions have been verified to yield similar results. We use single spin-flip dynamics at zero temperature, and vary the magnetic field $`H`$ adiabatically. Thus as the magnetic field is changed, a spin only flips if it is energetically favorable for it to do so. Once a spin flips, it can render other spins unstable. If more than one spin is unstable, it is important to have an algorithm that chooses the order that spins flip in a deterministic manner, to avoid introducing randomness into the system through the dynamics. Accordingly we flip the spin that lowers the energy the most. This process is repeated until all the spins are stable, after which the magnetic field can be further varied.foot:dyn

We start the system at large $`H`$ with all the spins positive and lower $`H`$ to $`H_{min}`$ adiabatically. The field is then raised to $`H_{max}`$. Thereafter, the field is cycled repeatedly between these two extremal values. After an initial transient period, the system reaches steady state. In this steady state we determine the number of cycles of $`H`$ that are required for the system to repeat its configuration. This varies randomly, depending on the coupling strengths $`J_{i,j}`$.

Clearly if $`H_{min}`$ and $`H_max`$ are too large in magnitude, the magnetization will saturate at the extremal points of the $`H`$ cycles and the behavior seen is trivial. On the other hand, if the range of $`H`$ is too small, hysteretic effects are minimal and behavior of the kind above will be hard to observe. By trying various values of the extremal fields, we find that it is best to choose them to be about half of the saturation field. Consequently we work with $`H_{max}=H_{min}=1.4`$.

Figure 1 shows an example of a system with linear dimension $`L=8,`$ i.e. with $`8^3`$ spins, where the steady state behavior is a two-cycle, i.e. two cycles of the magnetic field correspond to one cycle in the magnetization. Even though the gap between the two halves of the magnetization loop is much smaller than the width of the loop itself, with an oscillatory magnetic field, it would be straightforward to detect the subharmonic response. Although more complicated cycles are often found, they are not shown in the figure for the sake of clarity.

More generally, $`m`$ cycles of the magnetic field can be required for the (steady state) magnetization to return to its initial value. The distribution of $`m`$ for systems of different sizes is shown in Table 1. One observes that the larger the system size, the more likely it is to find long cycles. In fact for our largest size, with $`16^3`$ spins, for about 1% of the systems we were unable to detect any periodicity even when checking for $`m`$ up to 128. The fact that for the remaining systems, where periodicity was seen, the biggest value of $`m`$ detected was 42, suggests the existence of truly aperiodic behaviorperiodicity\_footnote .

Although Table 1 shows an increase in typical cycle size as the system is made larger, to some extent this is not caused by the dynamics becoming more complex, but because of the way that the cycle length is defined. For instance if the system were divisible into two independent halves, with a 6 and 7 cycle respectively, the full system would need 42 cycles to return to its original configuration. As the system is made larger, the number of such subsystems could proliferate, driving up the apparent cycle length. To eliminate this kind of behavior, it is useful to examine the power spectrum of these systems. One measures the magnetization at $`H_{min}`$ and $`H_{max}`$, the extrema of the cycles of $`H`$. If the magnetization goes through an $`m`$-cycle, these measurements will have a periodicity of $`2m,`$ while corresponding measurements of $`H`$ would oscillate between $`H_{min}`$ and $`H_{max}`$, exhibiting a periodicity of 2. The discrete Fourier transform of the magnetization measurements is taken and the power spectrum is calculated. In the case above, the spectrum will show peaks only at multiples of $`f=1/6`$ and $`f=1/7`$, where $`f=1`$ is the frequency of the $`H`$ field, whereas if the system were truly indivisible, peaks at multiples of $`f=1/42`$ would be observed. In addition, the power spectrum is crucial to detecting these effects experimentally, as we have discussed in connection with figure 1. Although in an experimental situation, one would measure the power spectrum for the full time dependent magnetization $`M(t)`$ instead of considering only its values at the extrema of $`H`$, the discrete Fourier transform computed here should be a good indicator of the experimental results.

Figure 2 shows the power spectrum for $`L=16`$. This is an average of the spectra for all the different $`m`$-cycles, weighted according to their probabilities (given in Table 1. Sharp peaks are seen at discrete frequencies, with the dominant peak at $`f=1/2`$. Figure 3 is an enlarged plot of the power spectrum at low frequencies for $`L=4`$ and $`L=16.`$ Although small, the low frequency power is clearly greater for $`L=16.`$ In addition, the peaks go down to lower frequencies, demonstrating that the dynamics get more complex as $`L`$ is increased in addition to the trivial growth in cycle length discussed in the previous paragraph.

In Figure 2 we have not plotted the two trivial peaks, at $`f=0`$ and $`1`$. These are approximately $`10^5`$ times the strength of the peak at $`f=1/2`$. As one would expect from Table. 1, the corresponding power spectrum for $`L=4`$ has much fewer peaks. However, the strength of the peaks is larger: the peak at $`f=1/2`$ is approximately $`10^4`$ times the strength of the peaks at $`f=0`$ and $`1.`$ More generally, as $`L`$ is varied, the strength of the trivial peaks will be proportional to the number of spins, i.e. $`L^3`$. Although the limited range of $`L`$ that we have been able to simulate does not allow us to determine the dependence on $`L`$ of the other peaks, for large $`L`$ we might expect any system to behave as a collection of reasonably independent regions, and the strength of the low frequency peaks to be proportional to $`L^{3/2}`$. However, this argument is hazardous, since even a small coupling between different regions of the system might affect subtle details of the dynamics. On the basis of our numerics, we make the more conservative claim, that the strength of the low frequency peaks increases with $`L`$, while their relative strength (compared to the peaks at $`f=0`$ and 1) decreases with $`L.`$ In a phase-locked experimental setup, the absolute rather than the relative strength of these peaks is important.

We now turn to experimental considerations. Nanoparticles of spin glass material, such as CuMn, with of order $`16^3`$ separate spins would have a linear dimension of a few nanometers, which could be readily manufactured. In fact it should also be possible to reduce the grain size somewhat, although as discussed in the previous paragraph the subharmonic content is reduced. In order to eliminate possible problems from interaction between the grains, it would advisable to coat or embed the grains in a neutral material to isolate them sufficiently from each other. With $`N_g`$ independent grains, the strength of the signal would be proportional to $`\sqrt{N_g}.`$

As discussed earlier, the field should be varied over a range of order half the saturation field, which would be $`1T`$. The time dependence of the magnetic field should be sufficiently slow so as to evolve the spins adiabatically, and to prevent any significant eddy current heating. For the small grains that we consider, a frequency of the order of $`10Hz`$ should be adequate.

The critical temperature $`T_c`$ of $`CuMn`$ is of order $`1K`$; in order to see the effects that we have discussed, which are very small compared to the overall hysteresis loop, it is necessary to essentially eliminate thermal noise over the time scales we are interested in. With activated dynamics, it should be possible to achieve this if the temperature is of order $`10mK`$.

It should be possible to satisfy all these conditions in an experiment without too much difficulty. We note that recent work has found super spin glass behavior in nanoparticle multilayer composites ssg  and in frozen ferrofluids nordblad , where the effective individual spins are large and $`T_c`$ is correspondingly high. While this would render the phenomena we wish to observe accessible at higher temperatures, one should be careful because the saturation field increases with $`T_c`$.

From a theoretical perspective, it is of interest to ask whether the results we have found here are generic for systems with hysteresis. As mentioned earlier in this paper, there are strong constraints on many ferromagnetic models that force them to have return point memory, precluding the possibility of multi-cycle hysteresis loops. Even with the spin glass model that we have considered in this paper, we find numerically that in one dimensional systems, when $`H`$ is cycled repeatedly, the hysteresis loop always settles down to a one-cycle. However, steady state multi-cycles are seen in two dimensions, at rates comparable to those for three dimensional systems.

From a dynamical systems approach, it would be natural to expect a system with many interacting degrees of freedom to exhibit complex dynamical behavior, possibly even chaos. In order to verify whether this is indeed the source of the results reported here, we have carried out similar numerical simulations on the Sherrington Kirkpatrick model SK  for spin glasses, as an example of a coupled non-linear system with several degrees of freedom (the spins). We find that one needs a minimum of $`5`$ spins to find multi-cycle steady state magnetization. In this case, one finds only three cycles and no others aside from the trivial closed loop case. The fraction of systems showing three cycles for $`5`$ spins is low, approximately $`3\times 10^4`$. However this fraction rises steeply with the number of spins. One has to be careful in applying general results from dynamical systems, since the dynamics we consider here are of a special kind, being driven by energy minimization. For instance, with a time-independent magnetic field and any initial configuration, the system must evolve towards a fixed point. Further work in this direction may explain the essential ingredients necessary to see the behavior reported in this paper.

Our results for the SK model suggest that systems which are described by a Landau theory with five or more components to the order parameter might, in the correct parameter regime, show multi-cycle hysteresis behavior. This would have the advantage that the effect would show up in intensive quantitites. As discussed earlier in this paper, for the spin-glass system we have considered here the non-return point effects for the average magnetization per spin should vanish in the large size limit. Zimanyi

In conclusion, we have shown in this paper that hysteresis loops have much richer behavior than commonly assumed. For Ising spin glass models with zero temperature dynamics, a slowly varying cyclical magnetic field produces a subharmonic response in the magnetization. Our estimates indicate that this effect should be observable in low temperature experiments. Broader implications for other systems have been discussed.

We thank David Belanger, Peter Young and Steve Ford for useful discussions.