text:
# Phonon deficit effect and solid state refrigerators based on superconducting tunnel junctions.

## I Introduction

Solid state micro–refrigerators ($`SSMR`$) have many advantages over conventional cooling devices: they have a long life-time, they are compact and simple to apply. These $`SSMRs`$ utilize motion of electrons rather than atoms and molecules. The flow of electric current through a semiconductor–semiconductor contact is accompanied by a transfer of heat (the Peltier effect) from one semiconductor to the other. Micro-refrigerators based on the Peltier effect are working at temperatures above $`150K`$. Solid-state micro-refrigerators operating below $`150K`$ are a scientific and technical challenge. A possible approach to this problem is to consider electron kinetics in superconductor–insulator–normal metal ($`SIN`$) or superconductor–insulator–superconductor ($`SIS`$) tunnel junctions caspek ; Martinis . As we will demonstrate in this article, another even deeper understanding of the same effects is possible from the concept of the phonon deficit effect ($`PDE`$). The phonon deficit effect was initially predicted GulianZharkov2 ; GulianZharkov1  to occur when superconductivity is enhanced in high frequency electromagnetic fields ($`UHF`$eliash . Later the $`PDE`$ was recognized to exist in more general cases GulBook ; GulianZharkov3 . It is related to the violation of detailed balance of kinetic processes in superconductors in the presence of an external field. If a superconducting film is immersed into a heatbath and the electron system goes over to an excited nonequilibrium state due to an external supply of energy, then under appropriate conditions the superconducting film absorbs phonons from the heatbath rather than emitting phonons.

The situation is essentially similar in the case of thin tunnel junctions of superconductors. Phonon emission of a symmetric $`SIS`$ tunnel junction has been considered in detail in Ref. GulianZharkov3  and an asymmetric junction has been considered in Refs. GulBook ; GulianSergoyan ; GulianMelkonyan . In the sub–threshold tunneling regime ($`V<2\mathrm{\Delta }/e`$, where $`\mathrm{\Delta }`$ denotes the gap) the resulting phonon flux is shown in Fig.. One observes the phonon emission in the frequency windows $`0<\omega _q<2\mathrm{\Delta }/\mathrm{}`$ and $`\omega _q>3\mathrm{\Delta }/\mathrm{}`$. In the same time there is the phonon absorption in the frequency window $`2\mathrm{\Delta }/\mathrm{}<\omega _q<3\mathrm{\Delta }/\mathrm{}`$. The analogy between this case and the superconductor in the $`UHF`$ field \[see Ref. GulBook , Figure 6.2\] strikingly demonstrates the universality of the mechanism of the phonon deficit effect.

Experiments on refrigeration are usually performed on the $`SIN`$ tunnel junctions. So called “cooling by extraction” is recognized as the operational principle of these type of refrigerators. At temperatures much below the critical temperature $`T_c`$ of the superconducting electrode ($`S`$) there are few thermally excited quasiparticles above the superconducting gap $`\mathrm{\Delta }`$, so that tunneling from $`S`$ to the normal metal electrode ($`N`$) can be neglected. Meanwhile, in the normal metal electrode ($`N`$) there are enough electron-like and hole-like thermal excitations near the Fermi energy $`E_F`$ which can tunnel to $`S`$ and effectively modify the electron distribution in $`N`$. The rate of electronic heat exchange with the lattice is proportional to $`T_e^5T_{ph}^5`$, where $`T_e`$ is electronic temperature and $`T_{ph}`$ is the temperature of phonons. The cooling process is a balance between phonon heat flow, quasiparticle extraction and Joule heating. The value of the applied potential $`V`$ on the tunnel junction determines if cooling or heating occurs depending on which process is dominant in the heat transfer balance.

Here we have studied the behavior of $`SIN`$ tunnel junction using the phonon deficit effect. The outline of this paper is as follows. In Section II we describe the phonon heatbath model and decoupling of electrons from phonons. In Section III we solve kinetic equations for the $`SIN`$ tunnel junction. We calculate the net phonon flux to the $`SIN`$ tunnel junction when electron and phonon systems are maintened at different temperatures. We show that the rule for heat exchange $`PT_e^5T_{ph}^5`$ follows from the $`PDE`$ in $`SIN`$ tunnel junctions. In analyzing the $`PDE`$ we propose the phonon filters. We present the results of our numerical simulations, which can be utilised to enhance the refrigeration effect. The constraction of phonon filters in accordance with the calculated phonon absorbtion spectra is discussed in Section IV. Conclusions are given in Section V.

## II Phonon heatbath model and decoupling of electrons from phonons

Experimentally junction electrodes are thin films placed on a substrate. These thin films should be considered as a system of two interacting subsystems: electrons and phonons. Nonequilibrium dynamics of phonons and electrons in superconductors is described by a coupled system of kinetic equations GulBook . In the simplest case, when the lattice phonons play the role of a thermostat, the kinetic equations for electron (hole) excitations and phonons in a metal can be written as in the following form GulBook

$`{\displaystyle \frac{n_{\pm ϵ}}{t}}=Q^f\left(n_{\pm ϵ}\right)+I^{elph}(n_{\pm ϵ},N_{\omega _q}^i),`$ (1)
$`{\displaystyle \frac{N_{\omega _q}^i}{t}}=L\left(N_{\omega _q}^i\right)+I^{phel}(N_{\omega _q}^i,n_{\pm ϵ}).`$ (2)
Here $`n_{\pm ϵ}`$ denotes the distribution functions of nonequilibrium electron-like ($`n_ϵ`$), respectively hole-like ($`n_ϵ`$) excitations in the film. $`N_{\omega _g}^i`$ represents the distribution function of internal phonons, $`I^{elph}(n_{\pm ϵ},N_{\omega _q}^i)`$ stands for the electron–phonon collision integral, and $`I^{phel}(N_{\omega _q}^i,n_{\pm ϵ})`$ denotes the phonon–electron collision integral. $`Q^f\left(n_{\pm ϵ}\right)`$ is an external source of nonequilibrium quasiparticles and $`f`$ represents external fields to be specified later. $`L(N_{\omega _q}^i)`$ is an operator describing the interaction of phonons with an external heatbath.

An important problem in a dynamical theory of metals is the study of nonequilibrium stationary states, i.e. when $`n_{\pm ϵ}/t=0`$, $`N_{\omega _q}^i/t=0`$. For thin films in an external field the dynamical problem can be simplified by considering the phonon heatbath model GulBook ; Eliashberg  and also the ”geometric-acoustical” approximation of phonon propagation. In this approximation the distribution function of internal phonons is given by $`N_{\omega _q}^i=N_{\omega _q}^0`$, where $`N_{\omega _q}^0`$ is the distribution function of phonons at thermodynamic equilibrium. The distribution of electrons and holes corresponds to a nonequilibrium state. In a steady state, using the phonon heatbath model, the following system of equations holds

$`Q^f\left(n_{\pm ϵ}\right)`$ $`=I^{elph}(n_{\pm ϵ},N_{\omega _q}^0),`$ (3)
$`L\left(N_{\omega _q}^0\right)`$ $`=I^{phel}(N_{\omega _q}^0,n_{\pm ϵ}).`$ (4)
These equations are uncoupled. Eq.(3) can be used to obtain stationary steady state solutions for the distribution functions of nonequilibrium excitations (electron and hole branches. Eq.(4) becomes an identity, which describes the phonon exchange (phonon flux) between the thin film and its environment GulBook .

## III $`SIN`$ tunnel junction - kinetic equations in the steady state

We consider now a massive superconductor (injector $`S^{}`$), having a superconducting gap $`\mathrm{\Delta }^{}`$, joined with a thermostat held at temperature $`T`$. A thin film having a thickness $`d\xi _0`$ (superconductor or normal metal) is attached to the massive superconductor by means of a thin oxide layer. It has the same temperature (see Fig.). In the dissipative steady state, nonequilibrium electron ($`n_ϵ`$) and hole ($`n_ϵ`$) excitations in the thin film ($`S`$ or $`N`$) are described by the following kinetic equation (see Eqs. (1) and (3) )

$$u_ϵ\frac{n_{\pm ϵ}}{t}=Q(n_{\pm ϵ})+I^{elph}(n_{\pm ϵ},N_\omega ^0)=0.$$
(5)
Here $`u_ϵ=|ϵ|\theta (ϵ^2\mathrm{\Delta }^2)/\sqrt{ϵ^2\mathrm{\Delta }^2}`$ denotes the BCS density of states, $`ϵ`$ is the quasiparticle energy and $`Q(n_{\pm ϵ})Q^f\left(n_{\pm ϵ}\right)`$ describes the injector. In a particular geometry the external field $`f`$ represents the distribution function of quasiparticles in the electrode denoted by $`S^{}`$. The explicit form of the electron–phonon collision integral $`I^{elph}(n_ϵ,N_{\omega _q}^0)`$ and $`Q(n_{\pm ϵ})`$ are given in Refs. GulBook ; GulianZharkov1 . The dissipative steady state solutions of Equation (5) give the electron and the hole distribution functions in the thin film electrode. The corresponding phonon fluxes emitted to (or absorbed from) its environment in a frequency interval $`d\omega _q`$ centered at frequency $`\omega _q`$ can be obtained from (4).

The number of phonons absorbed per time unit in the volume $`\mathrm{\Omega }`$ at frequency $`\omega _q`$ in the spectral interval $`\mathrm{d}\omega _q`$ is given by

$$\text{d}\dot{N}_{\omega _q}=L\left(N_{\omega _q}^0\right)\rho _0(\omega _q)\text{d}\omega _q=\rho _0(\omega _q)I^{phel}\left(N_{\omega _q}^0\right)\text{d}\omega _q.$$
(6)
Here $`\rho _0(\omega _q)=\mathrm{\Omega }\omega _q^2/(2\pi ^2u^3)`$, $`u`$ is the velocity of sound and $`\dot{X}X/t`$. The absorbed power per volume $`\mathrm{\Omega }`$ at frequency $`\omega _q`$ in the spectral interval $`\mathrm{d}\omega _q`$ is given by

$$\frac{\mathrm{d}P}{\mathrm{d}\omega _q}=\mathrm{}\omega _q\rho _0(\omega _q)I^{phel}(N_{\omega _q}^0).$$
(7)
Introducing the operator $`I(N_{\omega _q}^0)`$ via

$$I(N_{\omega _q}^0)=\frac{16\pi ϵ_F}{\lambda \omega _D}\omega _q^3I^{phel}\left(N_{\omega _q}^0\right),$$

allows to express the absorbed power (7) in dimensionless form

$`{\displaystyle \frac{\mathrm{d}P}{\mathrm{d}x}}={\displaystyle \frac{\mathrm{\Omega }\lambda \omega _D}{16\pi ϵ_Fu^3}}{\displaystyle \frac{(k_BT_c^{})^5}{\mathrm{}^3}}I(N_x^0)=P_0I(N_x^0)`$
$`P_0={\displaystyle \frac{\mathrm{\Omega }\lambda \omega _D}{16\pi ϵ_Fu^3}}{\displaystyle \frac{(k_BT_c^{})^5}{\mathrm{}^3}},`$ (8)
where we abbreviate $`x=\mathrm{}\omega _q/(k_BT_c^{})`$ and $`T_c^{}`$ denotes the critical temperature of the injector.

In the following we use the units $`k_B=\mathrm{}=e=1`$ if not stated otherwise. Also, recall that the primed variables correspond to the injector material.

### III.1 Numerical results and analysis

To simplify further the analysis, we need to analyze properties of equation (5). The electron–phonon collision integral $`I^{elph}`$ is proportional to the parameter $`\gamma T^3/\omega _D^2`$ (damping) which characterizes the strength of the relaxation processes. The parameter $`\nu `$ describes the intensity of the particle and hole injection $`Q(n_{\pm ϵ})`$. The typical values of $`\nu `$ are small compared to $`\gamma `$. When $`\nu /\gamma `$ tends to zero, the particle and hole distribution functions go over to thermal equilibrium distribution functions and both, $`I^{elph}(n_{\pm ϵ},N_{\omega _q}^0)`$ and $`I^{phel}(N_{\omega _q}^0,n_{\pm ϵ})`$ vanish. At finite values of $`\nu /\gamma `$ the particle and hole distribution functions differ from the thermal equilibrium functions and thus phonon fluxes arise. As observed in Refs. GulianSergoyan ; GulianZharkov1 ; GulianMelkonyan  both negative and positive phonon fluxes occur between the thin film and the thermal reservoir. The net phonon flux integrated over all frequencies can be either negative or positive.

Figure  shows the power $`P(V)/P(0)`$ absorbed by the normal metal thin film electrode. It is shown as a function of the potential applied across the $`S^{}IN`$ tunnel junction for different values of the parameter $`\nu /\gamma `$. All curves have a minimum (maximal absorption) at $`V<\mathrm{\Delta }^{}(T)`$. When the potential difference $`V`$ across the tunnel junction approaches the superconducting gap $`\mathrm{\Delta }^{}(T)`$, the phonon absorption decreases. The absorbed power increases with the increasing values of $`\nu /\gamma `$. The situation changes dramatically in the high temperature regime $`T/T_c^{}>0.5`$. When the value of the potential applied across the tunnel junction approaches the superconducting gap $`\mathrm{\Delta }^{}`$, then the $`S^{}IN`$ tunnel junction starts to emit phonons to the environment. This is shown in Fig.\[2b\].

In Fig.  we present the dynamics of phonon flux in the vicinity of the minima of Figs \[2a,b\]. As expected from Fig.\[2a\], the phonon absorption spectrum is spread out over the whole frequency range of acoustic phonons at the temperature $`T/T_c^{}=0.5`$ and $`\mathrm{\Delta }^{}V>0.18`$. At temperature $`T/T_c^{}=0.8`$ there is both phonon emission and phonon absorption at high and low frequency regions, respectively.

The phonon emission spectrum of a $`SIN`$ tunnel junction is different from the $`S^{}IS`$ case (see Fig.). The $`SIS`$ tunnel junction has two emission and one absorption frequency regions. In $`SIN`$ tunnel junctions the low frequency emission region is absent. Also, as temperature decreases the phonon emission decreases faster than the absorption.

### III.2 $`PDE`$ in the $`SIN`$ tunnel junctions:”Equilibrium” limit

An analytic expression for the power absorbed by the $`SIN`$ tunnel junction when electron and phonon subsystems are maintained at temperature $`T_e`$ and $`T_{ph}`$, respectively, can be obtained by integrating expression (III). The heatflow from the phonon subsystem to the electron subsystem maintained at different temperatures is given then by (for details see Appendix )

$$P=\alpha \left[\left(\frac{T_e}{T_c^{}}\right)^5\left(\frac{T_{ph}}{T_c^{}}\right)^5\right],$$
(9)
where $`\alpha =4P_0\mathrm{\Gamma }(5)\zeta (5)`$, $`\zeta (x)`$ denotes Riemann’s Zeta function and $`\mathrm{\Gamma }(n)`$ is the Gamma function. A comparison of Equation (9) with the Equation (19) of Ref. PBAllen  shows that $`PT_e^5T_{ph}^5`$ is compatible with the $`PDE`$ in $`SIN`$ tunnel junctions \[see Refs. FCWellstood ; Pekola \]. This result shows that the usual approach caspek ; caspek1 ; Martinis ; Pekola  takes into account only the net phonon flux. In this approach details like absorption at low frequency region and emission at high frequency region (see Fig. 3 ) are inaccessible. However, the details of phonon absorption spectra are important and may open new possibilities for cooling by $`SIN`$ tunnel junctions as discussed in the next Section.

## IV Microrefrigerator based on the phonon deficit effect

Why is it important to know the spectral dependence of the phonon fluxes? How can the phonon-deficit-effect influence the microrefrigerator design? The answers to these questions can be found in spectral-selective phonon filters.

To clarify this point, we should calculate the energy, which the nonequilibrium phonons yield to the heatbath. This quantity is equal to the integral of the spectral function plotted in Fig.. This integral having a net negative value is likely to be higher in the case of an asymmetric $`S^{}IS`$ junction GulianSergoyan ; GulianMelkonyan  and also in asymmetric $`S^{}IN`$ junctions. Fig. displays the nonequilibrium phonon flux from the $`S^{}IS`$ tunnel junction, $`\mathrm{\Delta }_S<<\mathrm{\Delta }_S^{}`$.

As follows from Figs.,  and  the net flux integrated over energy is negative GulianSergoyan ; GulianSergoyan2 ; GulianMelkonyan . To enhance the $`PDE`$ effect and achieve better cooling, one should separate nonequilibrium negative and positive fluxes across the interface between the superconductor and the attached heat reservoir. Let us consider placing a phonon spectral filter between superconductor $`S`$ and reservoir 2 (Fig.).

Suppose we can construct a filter AGulian ; GulianMelkonyan1  which is transparent to phonons in a window of frequencies where the flux in Fig. is negative. Such a filter would permit the flow of phonons ($`\mathrm{}\omega _𝐪2\mathrm{\Delta }`$) from the reservoir 2 to the superconductor $`S`$. At the same time the filter will absorb phonons below and above that frequency window, including those coming from the heated barrier. Another filter with complementary (i.e. the band pass) transmission properties may be placed between reservoir 1 and the electrode $`S^{}`$. It will prevent phonons coming from $`S^{}`$ (i.e., the reservoir 1) from being absorbed in $`S`$. This will make the absorption from the reservoir 2 more efficient. Then the system shown in Fig. will work as a refrigerator, by coolingthe reservoir 2. This cooling effect cannot reduce the temperature much below the critical temperature of superconductor $`S^{}`$ (Fig.). However, there are different classes of superconductors with critical temperatures covering a wide range of temperatures from about $`150K`$ to very low temperatures. Then, cooling cascades can be organized in such a way to generate substantial cooling. We suggest that such designs of a $`PDE`$-based microrefrigerator could be carefully considered for the purpose of practical use. In the next section we will discuss some designs.

### IV.1 Design of phonon filters

From the analysis in the preceding subsection it follows that the design of appropriate phonon filters is of primary importance. Using the analogy between reflection of a plane electromagnetic wave at the interface between two optical media with different refractive indexes $`n_1`$ and $`n_2`$ and the acoustic wave reflecting at the boundary of two elastic media with different acoustic impendance $`Z_1`$ and $`Z_2`$, Narayanamurti et al. Narayanamurti ; Koblinger  have developed appropriate phonon filters for phonon spectrometery.

Here we consider the simplest case of a phonon filter GulianMelkonyan ; GulianMelkonyan1  for the case when acoustic phonons propagate through a superlattice perpendicular to interfaces. In Fig. the superlattice system is shown.

The system consists of a sequence of different layers: material 1 with thickness $`d_1`$ and impedance $`Z_1`$, material 2 with thickness $`d_2`$ and impedance $`Z_2`$, etc. If one repeats this configuration $`p`$ times then the superlattice will have a period $`D=d_1+d_2+\mathrm{}+d_n`$ and $`L=Np`$ layers. The interfaces are perpendicular to the wave vector of the incident wave. The superlattice is placed between two media with the acoustic impedance $`Z_{medium}`$ and $`Z_{substrate}`$. Methods of calculation for the transmission and reflection coefficients in a linear stratified medium are described in Refs. Pierce ; Richards ; Shelkunoff ; ShunLien . We have presented results of our calculation in Fig..

One should note that all filters have an oscillatory behavior for the transmission coefficient in the band pass region. This is typical for band pass filters. In Tab. we list impedance and thickness of the superlattices used in our calculations. Fig. presents curves for the transmission spectra of the filter between reservoir 2 and superconductor $`S`$ with different acoustic impedance, layer thickness and layer number. Comparing with Fig. we observe large phonon fluxes at frequencies about $`\omega _{bfq}=\mathrm{\Delta }/\mathrm{}`$ and $`\omega _{bfq}=3\mathrm{\Delta }/\mathrm{}`$. The filter shown in Fig.\[6d\] may work well for those frequencies. The filter between reservoir 1 and superconductor $`S^{}`$ is presented in Ref GulianMelkonyan1 . The phonon filters for the $`SIN`$ tunnel junctions can be designed by a similar analysis.

### IV.2 Phonon filters and refrigeration

The proposed $`PDE`$ refrigerator is suited to refrigerate an object attached to it. The effect of phonon filters can be estimated by the value of $`P/P_0`$ from an experiment carried out by Pekola et al Pekola  where a membrane attached to the $`SIN`$ tunnel junction was cooled. Our estimate for the net phonon absorption $`P/P_0`$ by Equation (9) varies from $`0`$ to $`10^2`$ depending on temperatures $`T_e`$ and $`T_{ph}`$. Table 2 presents our calculated values for the nonequilibrium tunnel junctions (see also Figure 3). We find qualitative agreement between experiment and calculation for the value of absorbed power.

The phonon filters will increase the absorbed power by cutting off the emission window (see Figure 3). Table 3 presents the absorbed power after placing a low pass filter (see Figures 3 and 7) between the $`SIN`$ tunnel junction and an object. Phonon absorption is increased at $`T/T_c=0.5`$. At $`T/T_c=0.8`$ we see that $`SIN`$ tunnel junction with phonon filter will absorb phonons (the minus sign in Table 3 means absorption).

Fig.  shows the temperature dependence of the total phonon absorption $`P/P_0`$ for three different casses. The solid line is the phonon absorption by a $`SIN`$ tunnel junction without phonon filter. The dotted and dashed curves present the filtered phonon absorption for two different low pass filters with transmittance $`1`$ at $`\omega <\mathrm{\hspace{0.17em}2.95}T_c`$ and $`\omega <\mathrm{\hspace{0.17em}3.25}T_c`$. Thus the phonon absorption is incresed in a wide temperature region.

We recall that the refrigeration schema presented above is suited to refrigerate an object attached to it. But as follows from the phonon absorption spectra an analogous design with high pass filters can prevent the phonon absorption by $`SIN`$ tunnel junction and thus reduce the heating of electrons by phonons (see Figures ).

## V Conclusions

The $`SSMR`$ performance depends on various external and internal parameters. The $`PDE`$-based approach presented here and the usual apprach PBAllen ; FCWellstood ; Pekola  give the similar expression (Eq. 9) for the net cooling in case of $`SIN`$ tunnel junction. Thus one can conclude that the cooling of a membrane in the experiment by Pekola et al. Pekola  is a direct implication of the phonon deficit effect. The $`PDE`$approach offers explicit spectral description for participating phonon fluxes. It provides receipes to overcome unwanted heating channels and shows new details which have been omitted in a simplified ”cooling by extraction” description.

We have analyzed the “phonon deficit effect” in $`SIN`$ tunnel junctions. In contrast to $`SIS`$ tunnel junctions there is a temperature region where the $`SIN`$ tunnel junctions essentially absorb low frequency phonons. As temperature is increased, the $`SIN`$ tunnel junctions also emit phonons to the environment. Depending on the value of the ratio $`\nu /\gamma `$ there is net phonon emission or absorption. The temperature dependence of the net phonon absorption is consistent with experimental observations Pekola ; caspek ; Martinis : there is heating at high and low temperatures and refrigeration at intermediate temperatures.

To enhance the cooling power of the $`SIS`$($`N`$) and tunnel junctions and overcome the phonon heating at high temperatures, phonon filters can be effective. We performed calculations for the microrefrigerator phonon filters. Those filters have superlattice design and can be obtained using contemporary nanotechnology methods. They can enhance the cooling power of $`SIS`$ and $`SIN`$ tunnel junctions.

###### Acknowledgements.

H.K. has been supported by NSERC Canada. G.M. is grateful for support by the “Programme québécois de bourses d’excellance” of the governement of Québec.

*

## Appendix A Phonon flux between electrons and Phonons

In this appendix, we give the derivation of the expression (9).

The general expression for the phonon flux is given by (compare with the formula (7), $`\mathrm{}=1`$,$`k_B=1`$)

$$P=_0^{\omega _D}\omega _q\rho _0(\omega _q)I^{phel}(N_{\omega _q}^0)d\omega _q.$$
(10)
Here $`I^{phel}(N_{\omega _q}^0)`$ is the phonon-electron coillision integral and it is given by

$`I^{phel}\left(N_{\omega _q}\right)`$ $`=`$ $`{\displaystyle \frac{\pi \lambda \omega _D}{8ϵ_F}}{\displaystyle _0^{\mathrm{}}}\mathrm{d}ϵ{\displaystyle _0^{\mathrm{}}}\mathrm{d}ϵ^{}\{𝒮_1\delta (ϵ+ϵ^{}\omega _q)`$ (11)
$`+2𝒮_2\delta (ϵϵ^{}\omega _q)\},`$
where $`𝒮_1`$ and $`𝒮_2`$ are elementary collison processes. When the distribution function of electrons and holes is Fermi function then $`n_ϵ=n_ϵ`$. In this case $`𝒮_1`$ and $`𝒮_2`$ can be written as

$$𝒮_1(ϵ,ϵ^{},\omega _q)=4\left[\left(N_{\omega _q}+1\right)n_ϵn_ϵ^{}N_{\omega _q}\left(1n_ϵ\right)\left(1n_ϵ^{}\right)\right],$$
(12)
$$𝒮_2(ϵ,ϵ^{},\omega _q)=4\left[\left(N_{\omega _q}+1\right)n_ϵ\left(1n_ϵ^{}\right)N_{\omega _q}\left(1n_ϵ\right)n_ϵ^{}\right].$$
(13)
After substitutions of the expressions (12), (13) and (11) in the expression (10), the integration by delta function gives

$`P`$ $`=`$ $`{\displaystyle \frac{4\pi \lambda \omega _D}{8ϵ_F}}{\displaystyle _0^{\omega _D}}\mathrm{d}\omega _q\omega _q\rho _0(\omega _q)[{\displaystyle _0^{\omega _q}}𝒮_1(ϵ,\omega _qϵ,\omega _q)\mathrm{d}ϵ`$ (14)
$`+{\displaystyle _{\omega _q}^{\omega _D+\omega _q}}2𝒮_2(ϵ,ϵ\omega _q,\omega _q)\mathrm{d}ϵ].`$
Finally evaluating these integrals by noting that $`T\omega _D`$ and using the integral representation of Rieman zeta function $`\zeta (n+1)\mathrm{\Gamma }(n+1)=_0^{\mathrm{}}x^n/(exp(x)1)𝑑x`$ one obtains

$$P=4\frac{\mathrm{\Omega }}{16\pi }\frac{\lambda \omega _D}{ϵ_Fu^3}\mathrm{\Gamma }(5)\zeta (5)\left(T_e^5T_{ph}^5\right).$$
(15)


format: text

text:
# Dephasing from interactions and spin disorder

## Abstract

We calculate the dephasing rate of the electrons in the presence of interactions and elastic spin disorder scattering. In the frame of a self-consistent diagrammatic treatment, we obtain saturation of the dephasing rate in the limit of zero temperature for spin-orbit disorder in 2 dimensions. This result is in agreement with relevant experiments.

The dephasing rate $`\tau _\varphi ^1`$ provides a measure of the loss of coherence of the carriers. The saturation of the dephasing rate at low temperature seen in some samples  has attracted a vigorous interest.

In the presence of spin-less disorder, the Cooperon (particle-particle diffusion correlator) is given by

$$C_0(q,\omega )=\frac{1}{2\pi N_F\tau ^2}\frac{1}{Dq^2i\omega +\tau _\varphi ^1}.$$
(1)
$`D`$ is the diffusion coefficient, $`N_F`$ is the density of states at the Fermi level and $`\tau ^1`$ the total impurity scattering rate. We will work in the diffusive regime $`ϵ_F\tau >1`$ ($`\mathrm{}=1`$), $`ϵ_F`$ being the Fermi energy.

Previous studies  have focused on the calculation of $`\tau _\varphi ^1`$ in the absence of spin-scattering disorder. Here we calculate $`\tau _\varphi ^1`$ in the presence of spin-scattering disorder, whereupon the Cooperon becomes spin-dependent. The relevant terms $`C_i`$ are shown in fig. 1. We start by giving the explicit form of these $`C_{1,2}^o`$ without a dephasing rate. For the case $`\tau _{so}^1>0,\tau _S^1=0`$ \- with $`\tau _{so}^1`$ the spin-orbit impurity scattering rate and $`\tau _S^1`$ the magnetic impurity scattering rate - they are given by

$$C_{1,2}^o=b_{1,2}\left[\frac{1}{Dq^2i\omega }\frac{1}{Dq^2i\omega +3/(4\tau _{so})}\right],$$
(2)
with $`b_1=(3\tau _{so}/\tau 2)/(2u),b_2=1/u,u=4\pi N_F\tau ^2`$. We should emphasize that the impurity scattering considered is elastic.

To calculate the dephasing rate, we solve the appropriate Bethe-Salpeter equations for all three Cooperons $`C_i(q,\omega )`$, $`i=0,1,2`$. These are shown schematically in fig. 2:

$`C_0=C_0^o+C_0^oWC_0,`$ (3)
$`C_1=C_1^o+C_1^o[FC_1+RC_2]+C_2^o[RC_1+FC_2],`$ (4)
$`C_2=C_2^o+C_2^o[FC_1+RC_2]+C_1^o[RC_1+FC_2].`$ (5)
In fig. 3 we show explicitly the components of the self-energy terms $`\mathrm{\Sigma }=W,F,R`$. These quantities are given by

$`W=\mathrm{\Sigma }_0+d_0\mathrm{\Sigma }_1,F=\mathrm{\Sigma }_1+d_1\mathrm{\Sigma }_0,R=d_2\mathrm{\Sigma }_2.`$ (6)
Here, $`d_0=d_1=\{1/(ϵ_F\tau )+4\pi \}/(2\pi ^2ϵ_F\tau _{so})`$ and $`d_2=2\{1/(ϵ_F\tau _{so})+2\pi \tau /\tau _{so}\}/\pi `$. The terms containing the factors $`d_i`$ \- with a spin impurity line either looping around the Cooperon or crossing it - provide the coupling between the spin-independent and the spin-dependent Cooperons, and they turn out to be crucial for the saturation of the dephasing rate obtained below. In the spirit of ref.  we obtain

$$\mathrm{\Sigma }_i\mathrm{\Sigma }_i(q=0,\omega =0)\underset{q}{}_{\mathrm{}}^{\mathrm{}}𝑑x\frac{C_i(q,x+i0)\text{Im}V(q,x+i0)}{\mathrm{sinh}(x/T)},$$
(7)
where

$$V(q,\omega )=\frac{v_q}{1+v_q\mathrm{\Pi }(q,\omega )},\mathrm{\Pi }(q,\omega )=\frac{N_F}{2}\left[\frac{Dq^2}{Dq^2i\omega }+\frac{Dq^2+\tau _{sp}^1}{Dq^2i\omega +\tau _{sp}^1}\right].$$
(8)
Here $`v_q=2\pi e^2/q`$ is the bare Coulomb interaction and $`\tau _{sp}^1=(4/3)(\tau _S^1+\tau _{so}^1)`$ the total spin scattering rate. We make the approximation  $`_{\mathrm{}}^{\mathrm{}}𝑑xF(x)/\mathrm{sinh}(x/T)T_T^T𝑑xF(x)/x`$. The bosonic modes with energy greater than $`T`$ manifestly do not contribute to the self-energy and the dephasing process.

Since the coefficients $`d_i1`$ for $`ϵ_F\tau _{so}1`$, we can decouple a $`2\times 2`$ system of equations involving only $`C_{1,2}`$, to facilitate the solution of eqs. (3-5). Subsequently, $`C_{1,2}`$ are given by

$$C_{1,2}(q,\omega )=S_{1,2}\left[\frac{1}{Dq^2i\omega +r_{}}\frac{1}{Dq^2i\omega +r_+}\right].$$
(9)
Here $`r_+=(b_1b_2)(\mathrm{\Sigma }_2\mathrm{\Sigma }_1),r_{}=(b_1+b_2)(\mathrm{\Sigma }_2+\mathrm{\Sigma }_1)`$, $`S_i=a_o\{b_i(Dq^2i\omega )+M_i\}/X,X=2(b_1\mathrm{\Sigma }_2+b_2\mathrm{\Sigma }_1)`$ and $`a_o=8N_F\tau ^4/D`$.

The self-energies are given by

$$\mathrm{\Sigma }_i=\frac{Ta_o}{X}\left\{(M_ib_ir_{})\mathrm{ln}\left(\frac{z_or_{}}{Tr_{}}\right)(M_ib_ir_+)\mathrm{ln}\left(\frac{z_or_+}{Tr_+}\right)\right\},$$
(10)
with $`i=1,2`$, $`z_o=\tau ^1`$ and $`M_1=(b_2^2b_1^2)\mathrm{\Sigma }_1,M_2=(b_2^2b_1^2)\mathrm{\Sigma }_2`$.

Taking

$$b_1=b_2(\tau _{so}=\frac{4}{3}\tau ),$$
(11)
we obtain a simplified version of these equations:

$$\mathrm{\Sigma }_1=Ta_ob_1\mathrm{ln}\left(\frac{z_o+4b_1\mathrm{\Sigma }_1}{T+4b_1\mathrm{\Sigma }_1}\right),$$
(12)
and $`\mathrm{\Sigma }_2=\mathrm{\Sigma }_1`$. In the limit $`T0`$ we obtain the solution

$$\mathrm{\Sigma }_1=z_o/(4b_1)+O(e^{1/T}).$$
(13)
This finite solution for $`\mathrm{\Sigma }_{1,2}`$ for $`T0`$ is a generic fact for any finite value of $`\tau _{so}^1`$.

In 1-D we obtain similarly

$$\mathrm{\Sigma }_i=\frac{Ta_{1D}}{X}\left\{\frac{M_ib_ir_{}}{\sqrt{r_{}}}\left[f(z_o,r_{})f(T,r_{})\right]\frac{M_ib_ir_+}{\sqrt{r_+}}\left[f(z_o,r_+)f(T,r_+)\right]\right\},$$
(14)
with $`f(x,y)=\mathrm{ln}\{(\sqrt{x}\sqrt{y})/(\sqrt{x}+\sqrt{y})\}`$ and $`a_{1D}=2\pi \sqrt{D}a_o`$. In 3-D we have

$`\mathrm{\Sigma }_i={\displaystyle \frac{Ta_{3D}}{X}}\left\{b_iX\sqrt{z_o}+(M_ib_ir_{})\sqrt{r_{}}\left[f(z_o,r_{})f(T,r_{})\right](M_ib_ir_+)\sqrt{r_+}\left[f(z_o,r_+)f(T,r_+)\right]\right\},`$ (15)
with $`a_{3D}=2a_o/\pi `$. Both the 1-D and 3-D cases differ from 2-D in that no finite solution exists for $`\mathrm{\Sigma }_{1,2}`$ in the limit $`T0`$. This is due to the presence of the factors $`\sqrt{r_\pm }\sqrt{\mathrm{\Sigma }_i}`$ in the right-hand side of the respective eqs. for $`\mathrm{\Sigma }_i`$.

It turns out that

$$\tau _\varphi ^1b_0d_0\mathrm{\Sigma }_1=\frac{2[4\pi +1/(ϵ_F\tau )]}{\pi ^2\tau (ϵ_F\tau _{so})(3\tau _{so}/\tau 2)},$$
(16)
with $`b_0=1/(2\pi N_F\tau ^2)`$. Hence $`\tau _\varphi ^1`$ saturates in the limit $`T0`$ in 2-D. This is to be contrasted with the absence of spin disorder, where it has been shown that $`\tau _\varphi ^1T\mathrm{ln}T0`$ , and this fact yields eq. (16). Moreover, we would like to emphasize that other processes in $`\mathrm{\Sigma }`$, which are first order in the interaction $`V(q,\omega )`$, e.g. with $`V`$ crossing diagonally the Cooperon, do not modify qualitatively the result in eq. (16).

Now, the total correction to the conductivity can be written as

$$\delta \sigma _{tot}=\delta \sigma _o+\delta \sigma _{sp},$$
(17)
where the first term is the spin independent one - involving $`C_0`$ \- and the second term the spin dependent one - involving $`C_{1,2}`$. The latter is expected to saturate always in the low $`T`$ limit. The form usually fit to experiments  for the former is

$$\delta \sigma _o=\sigma _0\left[\mathrm{ln}(y)\mathrm{\Psi }(y+\frac{1}{2})\right],$$
(18)
where $`\mathrm{\Psi }`$ is the digamma function, $`y=1/4eHD\tau _\varphi `$, $`H`$ the magnetic field and $`\sigma _0=e^2/(2\pi ^2)`$. Obviously, this expression saturates for a saturating $`\tau _\varphi ^1`$, and vice-versa for the experimental determination of $`\tau _\varphi ^1`$.

A number of experiments show saturation of $`\tau _\varphi ^1`$ in the zero temperature limit. The samples in which saturation is observed are made of elements with a high atomic number, which induces a substantial spin-orbit scattering. Moreover, several of the samples, in which the dephasing saturation is observed, are truly 2-dimensional, such as the wires in refs. , the quantum dots in refs.  etc. We believe that the observed saturation can be understood in the frame of our results above.

I have enjoyed useful discussions/correspondence with Peter Kopietz and Pritiraj Mohanty.

e-mail : kast@iesl.forth.gr


text:
# Ordering of the lamellar phase under a shear flow

## Abstract

The dynamics of a system quenched into a state with lamellar order and subject to an uniform shear flow is solved in the large-$`N`$ limit. The description is based on the Brazovskii free-energy and the evolution follows a convection-diffusion equation. Lamellae order preferentially with the normal along the vorticity direction. Typical lengths grow as $`\gamma t^{5/4}`$ (with logarithmic corrections) in the flow direction and logarithmically in the shear direction. Dynamical scaling holds in the two-dimensional case while it is violated in $`D=3`$.

Many systems in nature exhibit lamellar order. One example is a diblock copolymer melt where chains of type A and B covalently bonded end-to-end in pairs segregate at low temperatures with A-B junctions forming a stack of lamellae . In ternary mixtures a lamellar phase is stable with ordered sheets of surfactant separating (say) oil and water domains . Lamellar order is also observed in Raleigh-Bénard cells, where convective rolls form above the convective threshold . Further examples include smectic liquid crystals , dipolar fluids with long-range interactions  and chemically reactive binary mixtures . A theoretical model for the general description of the lamellar-disordered phase transition was proposed by Brazovskii  who showed the first-order character of the transition induced by fluctuations.

Lamellar phases under an applied shear flow show a very rich behavior which is relevant for many applications . A variety of transitions in morphology and orientation occur as shear rate and temperature are changed. Stable configurations of lamellae lying along the flow with the normals differently oriented have been observed  and analyzed evaluating the effects of the flow on the fluctuation spectra .

Non-stationary properties are far less considered . In this Letter we study the effects of a simple planar shear flow on the ordering of the lamellar phase in a system quenched from an initially high temperature disordered state. We present the first analytical results on the kinetics of this system by solving the Brazovskii model in the limit of an infinite number of components of the order parameter. This is one of the few methods allowing an explicit solution for phase ordering systems . A similar approach for fluids under shear flow has been used in .

The behavior of quenched binary mixtures without imposed flows is characterized by dynamical scaling: The structure factor obeys $`C(\stackrel{}{k},t)=R(t)^Df[kR(t)]`$ where $`R(t)t^\alpha `$ is the typical domain size and $`D`$ is the space dimensionality . In the case of a fluid with lamellar order, if the order parameter is not conserved as in the Swift-Hohenberg model for Raleigh-Bénard convection , regimes exhibiting dynamical scaling have been found . In models with conserved order parameter, when diffusion is the only segregating physical mechanism, simulations show the entanglement of the fluid into frozen intertwined structures . In this case the effects of hydrodynamical modes are crucial for reaching order on large scales .

When shear is applied, our results show that lamellae grow preferentially with the perpendicular orientation , namely along the plane formed by the flow and the shear (velocity gradient) directions. Their typical size, obtained from the second momentum of the structure factor, grows as $`\gamma t^{5/4}\sqrt{\mathrm{ln}t}`$ in the flow direction and as $`\sqrt{\mathrm{ln}t}`$ in the shear direction. Surprisingly we find that dynamical scaling is obeyed in two dimensions but not in $`D=3`$. Our results concern the cases of conserved and not conserved order parameter, and apply to most of the systems mentioned above.

We consider the Brazovskii free energy

$$\{\varphi \}=𝑑\stackrel{}{r}\left\{\frac{r}{2}\varphi ^2+\frac{u}{4}\varphi ^4\frac{b}{2}|\varphi |^2+\frac{1}{2}(^2\varphi )^2\right\}$$
(1)
where $`\varphi (\stackrel{}{r},t)`$ is the order parameter field and $`u>0`$. With a negative value of $`r`$ and $`b<0`$ the system orders in one of the two minima of the local potential. However, when $`b>0`$ interfaces are favoured and a modulated state with wavevector $`k_M=\sqrt{b/2}`$ is stable .

When a flow is imposed, the kinetics can be described by the convection-diffusion equation

$$\frac{\varphi }{t}+\stackrel{}{}(\varphi \stackrel{}{v})=\mathrm{\Gamma }_p\frac{\delta }{\delta \varphi },p=0,2$$
(2)
where

$$\mathrm{\Gamma }_0=\mathrm{\Gamma },\mathrm{\Gamma }_2=\mathrm{\Gamma }^2$$
(3)
and $`\mathrm{\Gamma }`$ is a mobility coefficient. $`p=0`$ describes systems with non-conserved order parameter (NCOP) and corresponds to the Swift-Hohenberg equation ; the case $`p=2`$ is for conserved order parameter (COP) and applies for example to copolymer melts. For uniform shear flow in the $`x`$-direction $`v_x=\gamma y`$, $`\gamma `$ being the shear rate. Eq.(2) neglects thermal fluctuations and possible effects due to differences in viscosities between the two components . The complete description should take into account the coupling of Eq. (2) with the Navier-Stokes equation. However, the study of Eq. (2) is a prerequisite for any more general theory and usually applies to intermediate temporal regimes, as recognized for phase separation of simple binary mixtures .

Eq. (2) can be studied analytically generalizing the field $`\varphi `$ to a vector order parameter with $`N`$ components and taking the large-$`N`$ limit . One obtains the following equation in Fourier space for the structure factor $`C(\stackrel{}{k},t)=\varphi (\stackrel{}{k},t)\varphi (\stackrel{}{k},t)`$

$`{\displaystyle \frac{C(\stackrel{}{k},t)}{t}}`$ $``$ $`\gamma k_x{\displaystyle \frac{C(\stackrel{}{k},t)}{k_y}}=`$ (4)
$``$ $`2\mathrm{\Gamma }k^p\left[r+uS(t)bk^2+k^4\right]C(\stackrel{}{k},t)`$ (5)
where we have dropped the component indices due to internal symmetry. $`S(t)`$ has to be computed self-consistently through

$$S(t)=_{|k|<\mathrm{\Lambda }}\frac{d\stackrel{}{k}}{(2\pi )^D}C(\stackrel{}{k},t),$$
(6)
$`\mathrm{\Lambda }`$ being an ultraviolet cut-off. For a symmetric mixture quenched from a high temperature homogeneous phase an appropriate initial condition is $`C(\stackrel{}{k},0)=\mathrm{\Delta }`$, where $`\mathrm{\Delta }`$ is a constant. Eq. (5) can be integrated yielding

$$C(\stackrel{}{k},t)=\mathrm{\Delta }e^{2\mathrm{\Gamma }_0^t𝑑\tau \stackrel{}{𝒦}^p(\tau )\left[𝒦^4(\tau )b𝒦^2(\tau )+r+uS(t\tau )\right]}$$
(7)
where $`\stackrel{}{𝒦}(\tau )\stackrel{}{k}+\gamma \tau k_x\stackrel{}{e}_y`$. Defining

$$Q(t)_0^t𝑑\tau [r+uS(t\tau )]$$
(8)
the analysis can be carried out asymptotically through the ansatz

$$Q(t)=k_M^4t\frac{1}{2\mathrm{\Gamma }k_M^p}(2\mathrm{ln}t\mathrm{ln}\mathrm{ln}t+\mathrm{ln}v_p),$$
(9)
where $`v_p`$ is a constant to be determined. Eq.(9) will be justified a posteriori by proving the solution of the self-consistency problem. Specifically, the large time value $`S_{\mathrm{}}`$ calculated through (6), inserted into (8), must satisfy (9) asymptotically, namely

$$S_{\mathrm{}}=(k_M^4r)/u$$
(10)
We now briefly illustrate this for COP in $`D=3`$  setting $`\mathrm{\Gamma }=1=u=b=r`$.

The symmetry of the problem suggests the use of the cylindrical set of variables $`(k_x,k_{},\theta )`$, where $`k_y=k_{}\mathrm{cos}\theta `$ and $`k_z=k_{}\mathrm{sin}\theta `$. Naive power counting applied to Eq. (7) indicates that wave vectors in the flow direction scale asymptotically as $`k_xt^{3/2}`$ while in the perpendicular directions $`k_{}k_Mt^{1/2}`$. Then, introducing the new variables $`X=k_xL_{}(t)`$, $`Q=(k_{}k_M)L_{}(t)`$, with $`L_{}(t)=2\sqrt{2/3}k_M^2\gamma t^{3/2}`$ and $`L_{}(t)=2\sqrt{2t}k_M^2`$, and taking into account Eq. (9), Eq. (7) reads as

$$C(\stackrel{}{k},t)=v_2\mathrm{\Delta }\frac{t^2}{\mathrm{ln}t}e^{[f(X,Q,\theta )+g(X,Q,\theta ,t)]}$$
(11)
with

$$f=Q^2+\sqrt{3}QX\mathrm{cos}\theta +X^2\mathrm{cos}^2\theta $$
(12)
and $`g(X,Q,\theta ,t)=_{j=1}^{16}t^{j/2}g_j`$, where $`g_j`$ are polynomials in $`X,Q,\theta ,\mathrm{log}t`$. For long times $`g`$ can be shown to provide a relevant contribution $`g_2\frac{9}{160}\frac{X^4}{k_M^6}\frac{3}{4k_M^6}X^2\mathrm{ln}t`$ only in the region $`Q0,\mathrm{cos}\theta 0`$. Then the integration over $`Q`$ and $`\theta `$ can be performed yielding

$$S(t)\frac{\sqrt{3}\mathrm{\Delta }}{32k_M^3\pi ^{\frac{3}{2}}}\frac{v_2}{\gamma \mathrm{ln}t}_{\mathrm{}}^{\mathrm{}}𝑑Xe^{\frac{X^2}{8}\frac{g_2}{t}}I_0(X^2/8)$$
(13)
where $`I_0(z)`$ is a Bessel function. The integral of Eq. (13), evaluated by the asymptotic expansion of $`I_0(z)`$ at large $`z`$, behaves as $`\mathrm{ln}t/\sqrt{\pi }`$. Hence at large times $`S(t)`$ approaches the constant $`S_{\mathrm{}}=\frac{\sqrt{3}v_2\mathrm{\Delta }}{32\pi ^2k_M^3\gamma }`$. A comparison with Eq.(10) fixes $`v_2=\frac{32\pi ^2k_M^3\gamma (k_M^4+1)}{\sqrt{3}\mathrm{\Delta }}`$ and verifies the ansazt.

Due to the presence of $`g`$, the function $`C(\stackrel{}{k},t)`$ of Eq.(11) cannot be cast in a scaling form. However, different scaling behaviors are obeyed in the regions where the functions $`f`$ or $`g`$ can be respectively neglected. In particular, for $`Q0,\mathrm{cos}\theta 0`$, where $`g`$ dominates, $`C(\stackrel{}{k},t)`$ scales in the $`x`$-direction with respect to the length $`l_{}(t)=(2/5)^{1/4}\sqrt{k_M}\gamma t^{5/4}`$ different from $`L_{}(t)`$. In $`D=2`$, on the other hand, where the contribution of the function $`g`$ is always negligible, one finds the scaling form $`C(\stackrel{}{k},t)=v_2\mathrm{\Delta }t^2e^{[f(X,Q,0)]}`$.

Ordering properties are usually inferred from the momenta of the structure factor. We define $`R_x=[𝑑\stackrel{}{k}C(\stackrel{}{k},t)/𝑑\stackrel{}{k}k_x^2C(\stackrel{}{k},t)]^{1/2}`$ and similarly for the shear and vorticity directions. In the absence of dynamical scaling definitions of $`R_x`$ based on different momenta of $`C(\stackrel{}{k},t)`$ may lead to different results. However in this case it can be proven that changing the order of the momentum does not change the growth exponent but only the logarithmic correction. ¿From Eq. (11) we find $`R_x\gamma t^{5/4}\sqrt{\mathrm{ln}t}`$, $`R_y\sqrt{\mathrm{ln}t}`$ and $`R_zk_M^1`$. The behavior of $`R_x`$ shows the relevance of the scaling with respect to $`l_{}(t)`$ for the ordering of the system. The growth of $`R_y`$ indicates that lamellae order preferentially in the plane $`xy`$. In $`D=2`$ we find $`R_x\gamma t^{3/2}`$ and $`R_yk_M^1`$. The same results are obtained with non-conserved order parameter. The behavior of $`R_x,R_y,R_z`$ resulting from the numerical integration of Eq. (5) with an adaptive grid algorithm is shown in Fig. 1. After the initial isotropic evolution, the shear-induced anisotropy becomes evident for values of the strain $`\gamma t`$ larger than one and agreement with the analytical behavior is observed.

It is also useful to illustrate the behavior of the structure factor (11) shown in Fig. 2. The maxima of $`C(\stackrel{}{k},t)`$ are located at $`k_x^2=5\frac{1}{k_M^4\gamma ^2}\frac{\mathrm{ln}t}{t^3},k_y=0,k_z^2=k_M^2`$ for COP and at $`k_x=0,k_y^2+k_z^2=k_M^2`$ for NCOP. The shape of $`C(\stackrel{}{k},t)`$ on the planes $`k_x=0`$ and $`k_z=0`$ is qualitatively similar for COP and NCOP. At $`k_x=0`$, in particular, the wavevectors are not affected by the flow and the structure factor has a circular structure with radius $`k_M`$. At $`k_y=0`$ or $`k_z=k_M`$ the behavior of $`C(\stackrel{}{k},t)`$ depends on the conservation law. The maxima of $`C(\stackrel{}{k},t)`$ developed with NCOP are splitted with COP into a pair of symmetric narrowing peaks. This pattern is typical of the case with conserved dynamics and has been observed in other segregating systems under shear flow .

Finally we turn to the study of the rheological properties. In phase separation of binary mixtures the flow acts against the surface tension inducing stretching of domains followed by break-up processes and burst of small bubbles . An excess viscosity $`\mathrm{\Delta }\eta `$ is measured : it reaches a maximum generally for $`\gamma t1`$ and later decays. In our case we measure the evolution of the shear stress $`\sigma _{xy}`$ deducible in a general way from (1) in terms of the structure factor; it is given by $`\sigma _{xy}(t)=\frac{d\stackrel{}{k}}{(2\pi )^D}k_xk_y(2k^2b)C(\stackrel{}{k},t)`$ . An excess viscosity can be introduced as $`\mathrm{\Delta }\eta (t)=\sigma _{xy}(t)/\gamma `$ . Using Eq. (11) we find $`\mathrm{\Delta }\eta \gamma ^1t^2`$. The numerical results of Fig. 3 show that, after reaching a maximum, the decay of $`\mathrm{\Delta }\eta `$ agrees with the analytical behavior. The other rheological indicators can be shown to behave similarly.

In conclusion we have studied the ordering kinetics of a lamellar phase in shear flow by solving the dynamics of the convective-diffusion equation in the large-N limit. Regarding the debated question of the stable orientation, our results are in agreement with the expectation of the stability of the perpendicular phase at high shear rates . In $`D=3`$ a logarithmic growth law is found along the shear direction while $`R_x\gamma t^{5/4}\sqrt{\mathrm{ln}t}`$. Interestingly, the same exponent $`5/4`$ is found for simple binary fluids . Our results show that the scaling depends on dimensionality and is violated in $`D=3`$. This is at variance with the case without shear where the same approach gives a scaling form independently of $`D`$ . Violation of scaling is quite uncommon . In this context this phenomenon may be related to the existence of two lengths $`L_{}`$, $`l_{}`$ growing with different exponents.

A natural question is the relevance of the behavior of the large-N model to the physical case with $`N=1`$. As without flow, symmetry arguments suggest an analogy between the coarsening behavior of lamellar phases and that of vectorial models and one would argue that the large-N results give reliable indications also for the physical case ; in particular the large-N exponent is the one expected asymptotically . Simulations of the scalar case with shear could elucidate this point; however strong finite size effects make the evaluation of the growth exponent difficult (see  and refs. therein). We also mention that it will be important to study how the presence of hydrodynamics affects the picture provided in this Letter.

We thank Mario Pellicoro and Marco Zannetti for useful discussions. F.C. and G.G. acknowledge M.Cirillo and R. Del Sole for hospitality in Rome University and support by the PRA-HOP 1999 INFM and MURST (PRIN 1999).


text:
# Limits on Supersymmetric Dark Matter From EGRET Observations of the Galactic Center Region

## I Introduction

Observations by a variety of experiments have revealed that a great deal of the mass of our universe is dark and cold . Despite this growing body of evidence, we are still ignorant of the nature of dark matter.

One of the most promising dark matter candidates is the lightest neutralino in supersymmetric models . In most models, the lightest supersymmetric particle (LSP) is stable by the virtue of R-parity . Often, this particle is a neutralino, $`\chi ^0`$, the partner of the photon, Z-boson and neutral higgs bosons. This candidate is attractive due to the fact that it is electrically neutral, not colored and naturally has the approprite annihilation cross section and mass to provide a cosmologically interesting relic density.

Many methods have been proposed to search for evidence of supersymmetric dark matter. These include experiments which hope to measure the recoil of dark matter particles elastically scattering off of a detector (direct searches) , experiments which hope to observe the products of dark matter annihilation (indirect searches) and, of course, collider experiments . Indirect searches include searches for neutrinos from the Sun, Earth or Galactic center , positrons  or anti-protons  from the Galactic halo and gamma rays from the Galactic center and halo .

Methods of indirect detection which involve the Galactic center depend strongly on the distribution of dark matter in the Galaxy. At this time, there is a great deal of debate and speculation over the merits of various Galactic dark matter halo models. Numerical simulations favor models with strong cusps in the central region, such as the Navarro, Frenk, White (NFW) and Moore, et. al. models . These models predict increasing dark matter density as one approaches the Galactic center, $`\rho 1/r^\gamma `$, where $`\gamma `$ is 1.0 for the NFW case and 1.5 for the Moore case.

There have been arguments made, based on observations, in the favor of flat density core models. These distributions, although possible, are probably not capable of producing observable signals from dark matter annihilation, and are not discussed in this letter for this reason.

Models with strong density spikes at the center of the halo have recently receive some attention . In these models, cuspy halos generate spikes as a result of adiabatic acretion of matter into the central Galactic black hole.

Finally, if halo distributions are clumpy, rather than smooth, it would be possible that less dark matter would be present in the central region, and the dark matter signal diminished.

In this letter, we will show results for smooth, cuspy, halo distributions of the NFW and the Moore profiles. These distributions (as well as spikey models) are especially interesting to gamma ray experiments, as they provide signals from dark matter annihilation which appear as point sources. The angular distribution of events is proportional to the dark matter density squared integrated over the line of sight. A strongly cusped distribution produces the vast majority of the annihilation signal in an angular region much smaller than the point spread function of EGRET, and hence are indistinguishable from a point source.

## II EGRET Point Source Location Analysis

EGRET, the Energetic Gamma Ray Experiment Telescope, launched on the Compton Gamma Ray Observatory in 1991, is sensitive to gamma rays in the range of approximately 30 MeV to 30 GeV. During its operation, EGRET’s observations included an exposure of approximately $`2\times 10^9\mathrm{cm}^2\mathrm{sec}`$ in the direction of the Galactic center. This paper presents an analysis of only the gamma rays $`>`$ 1 GeV because the continuum spectrum of gamma rays from neutralino annihilation peaks at higher energies and because the point spread function of EGRET improves with energy.

Previous searches for point sources in EGRET data, such as the 3EG catalog , used a single mean point spread function for each observed gamma-ray above 1 GeV and spatially binned the data in square bins of sides 0.5 degrees. However, the EGRET point spread function significantly improves with increasing gamma-ray energy with 68% of 1 GeV gamma-rays reconstructed within 1.3 degrees of the true direction as compared to 68% of the 10 GeV gamma rays within 0.4 degrees. Our analysis uses the point spread function as determined by the preflight calibration  for 6 energy bins above 1 GeV, and does not degrade the reconstructed gamma-ray direction to the nearest 0.5 degree bin. We use a spatially unbinned maximum likelihood analysis to determine the best localization of a point source. The diffuse Galactic background is from the same model used for the production of the 3EG catalog. Before addressing the Galactic center region, we tested our method on well known sources such as the Vela pulsar and the Crab pulsar.

We found the position with maximum likelihood for the EGRET source near the Crab pulsar to be l=184.52,b=-5.79. The known location of the Crab is l=184.56, b=-5.78 which is within the 95% confidence region as determined by of our analysis. The 3EG catalog lists the location of this source as l=184.53, b=-5.84, with the known location well outside of the 95% confidence contour.

The results for the Vela pulsar are similar. We found the maximum likelihood at l=263.53, b=-2.82 with a known location of l=263.55, b=-2.79 again within the 95% contour. The 3rd EGRET catalog lists Vela at l=263.52, b=-2.86 and the known location is well outside of the 99% confidence contour.

When our technique is applied to the Galactic center region, we find a point source located at l=0.19, b=-0.08. The Galactic center is excluded as the source beyond the 99.9% confidence level, see figure 1. The 3EG 95% confidence region includes the Galactic center as shown by the circle in figure 1. We find that if this source, modeled with a differential power law spectrum with slope determined by the maximum likelihood technique to be -2.2, is included in the background of the region, the 95% confidence upper limit on the number of gamma rays from a point source at the Galactic center is 10 to 100 (depending on the spectrum of the source). By contrast, the source at l=0.19, b=-0.08 is a bright EGRET source of 370 gamma-rays. The identification of this source is unknown , but this new localization agrees well with a postulated source of inverse Compton gamma rays from the electrons which create the Galactic center radio arc .

## III Calculating The Gamma Ray Flux from the Galactic Center

We calculated, for a variety of supersymmetric models, the number of events EGRET would have been expected to have observed, as a function of the halo model. We only consider those models which do not violate accelerator limits , including b to s$`\gamma `$  and invisible Z decay width measurements. Furthermore, we require that the relic density of the LSP be $`0.05<\mathrm{\Omega }_\chi h^2<0.2`$. We calculate the neutralino relic density using the full cross section, including all resonances and thresholds, and solving the Boltzmann equation numerically . Coannihilations with Charginos and Neutralinos are included. We then calculate the LSP annihilation cross section, mass and resulting gamma ray spectrum, for a given halo model. The results are shown in figures 2 and 3.

The general supersymmetric parameter space, even for the minimal supersymmetric standard model, consists of more than 100 free parameter and, therefore, must be simplifed to do any practical calculations. We considered a 7-dimensional parameter space consisting of the gaugino mass parameter, $`M_2`$, the non physical mass $`\mu `$, the ratio of higgs vacuum expectation values, $`\mathrm{tan}\beta `$, a universal SUSY mass scale, $`M_{\mathrm{SUSY}}`$, the pseudoscalar higgs mass $`m_A`$, and the couplings $`A_t`$ and $`A_b`$.

We parametrized the continuum gamma-ray spectrum as a function of the LSP mass and calculated the 95% exclusion confidence levels which could be placed by the EGRET data. The parametrization depends on the neutralino annihilation branching fractions, but varies little in the majority of models. This exclusion contour is shown in figures 2 and 3 as a solid line. We did not consider the $`\gamma \gamma `$ or $`\gamma Z`$ line emission as these are generally above the energy range sensitive to EGRET.

We also considered the ability of the future experiment, GLAST, to probe the galactic center for dark matter annihilations. With larger area and better angular resolution, GLAST, will be capable of testing many more models than EGRET. Furthermore, GLAST, with sensitivity to energies as high as $``$ 1 TeV, can test models with heavier LSP neutralinos somewhat more easily than EGRET. In figures 2 and 3, the expected sensitivity of GLAST, after 3 years of observation, is shown as a dashed line.

## IV Conclusions

Our analysis of the EGRET data in the Galactic center region indicates an off-center point source, excluded beyond 99.9% as the Galactic center. Considering this source as background, we found no evidence of a point source at the Galactic center and determined the 95% confidence upper limits on the flux of gamma rays as a function of the WIMP mass.

We compared these limits to the flux predicted for a variety of supersymmetric models and galactic halo models. We find that for very cuspy (or spikey) halo models, such as the Moore et. al. profile, the majority of viable supersymmetric models are excluded by our limit. We show that the GLAST experiment will have the sensitivity to further constrain the Galactic halo profile of neutralino dark matter.

Acknowledgments: This work was supported in part by a DOE grant No. DE-FG02-95ER40896 and NASA grant No. NAG5-9712 and in part by the Wisconsin Alumni Research Foundation.


text:
# 1 Inequalities in the Bonabeau model, for density 0.4. The straight line symbolizes the absence of hierarchies (all 𝑞=1/2) for the lower density 0.3. Only at the higher density do the probabilities 𝑞 to win a fight concentrate on new very low and very high values.

Sociophysics Simulations

Dietrich Stauffer, Institute for Theoretical Physics

Cologne University, D-50923 Köln, Euroland

Abstract: Check if we need one

Introduction:

Physicists since a long time have tried to apply their skills to fields outside physics, with more or less success. Econophysics is a big fashion at present, though partly based on the wrong belief that economists did not make empirical observations of financial markets or made no Monte Carlo simulations where rational judgement is replaced by random decisions. Also Frederick Soddy ventured into economics after his 1921 Nobel prize in chemistry. Quantum mechanics co-inventor Erwin Schrödinger more than half a century ago wrote a book asking “What is life ?”, long before today’s interests in biophysics and bioinformatics. Sociophysics has been around for at least three decades, with or without that name . The present review summarizes some of the more recent simulations in sociophysics and is clearly biased by the personal preferences and experience of its author. In June 2002, I visited the first conference devoted only to Sociophysics, organized by F. Schweitzer and K.G. Treutsch, www.ais.fhg.de/ frank/. We ignore here car traffic simulations , scale-free networks , social percolation  and active Brownian particles  since they were reviewed recently. Instead we look at the models of Bonabeau et al , Sznajd  and similar consensus models . We concentrate on simple models which take about one page of Fortran program, available from stauffer@thp.uni-koeln.de. We thus update similar summaries published before .

Hierarchies:

How come someone is born into nobility, and others are members of the proletariat. Some scientists are given tenure, and others have to leave academia The elites of all countries and all times has always some explanations, like the Grace of God having them put into the upper levels of society. Statistical Physicists, of course, assume these hierarchies to arise from randomness. (The illusion that everything is random is a professional disease ”morbus Boltzmann” among these physicists, just as silicosis = black lung affects mine corkers.)

If nobility is connected with ownership of the land, then it cannot develop easily in a nomadic society, while sedentary societies may have ground property. The peasants then can become slaves of the nobility owning the piece of territory on which the peasants work. Thus sedentary societies with agricultural fields may develop stronger hierarchies than nomadic tribes with just a few light goods to be carried around. The Hollywood movie “Dances with Wolves” is an example how the industrialized world may imagine beautiful nomadic paradise to have been.

The problem therefore is to develop a model giving rise to strong hierarchies at high population densities and weak hierarchies at low population densities. Such a model, with a sharp first-order phase transition at some critical density, was proposed by Bonabeau et al  and followed-up by others . People diffuse on a square lattice to nearest neighbours (Margittai Neumann neighbourhood, not Moore neighbourhood). When one person tries to move to a place already occupied by another person, the two have a fight, the winner takes the contested place, and the loser moves to (or stays at) the other place. Initially, the probability to win or lose is 50 percent. But after some time, memories of past fights and their outcomes accumulate and are stored in the history $`h(i)`$ for each person $`i`$, where $`h`$ is the number of victories minus the number of defeats of this person. After each iteration, $`h`$ is diminished by, say, 10 percent of its current value to take into account that memories fade away.

The probability $`q`$ for $`i`$ to win against $`k`$ is assumed as

$$q=1/(1+\mathrm{exp}(\sigma [h(k)h(i)])$$
$`(1)`$
like a Fermi function: More victories in the past mean a higher chance to win now. Here $`\sigma `$, the inverse width of the transition at Fermi level, is the rms fluctuation in the winning probabilities:

$$\sigma =(<q^2><q>^2)^{1/2}$$
$`(2)`$
and thus is zero if everybody is equal ($`q=1/2`$), and is $`1/\sqrt{6}`$ if the values for $`q`$ are distributed homogeneously between zero and unity. Thus $`\sigma `$ measures the inequality in this society. And this inequality re-enters into the probabilities to win and lose, thus enhancing existing inequalities. (For the first ten Monte Carlo steps per site, $`\sigma =1`$ to allow a build-up of hierarchies.)

The simulations show either a $`\sigma `$ going after some initial positive values rapidly to zero with time (density below 32 percent), or a $`\sigma `$ staying at or above 0.25 for density above 32 percent. This first-order phase transition actually has a complicated history , and only the present status is summarized here. Fig.1 shows the distribution of $`q`$ values near the end of a simulation with and without hierarchy.

Even in a strongly hierarchical society, revolutions can happen, and then a new hierarchy builds up. Thus averaging for one individual over very long time would give an average winning probability $`q`$ near 1/2, for all people. The proper way to distinguish between hierarchical and non-hierarchical societies is thus the above snapshot method, where the inequality $`\sigma `$ is determined from all $`q`$ values at one given moment, via eq.(2).

F. Schweitzer, private communication, has criticized the symmetry built into this model, where the probability to belong to the leaders (high $`q`$) is the same as the one to belong to the followers (low $`q`$). A monarchy needs one king with many subjects. This is only one of many questions which are still open for further modelling.

Consensus:

A democracy needs both a stable government and a viable opposition. But for a selection comittee trying to fill a professor position at the university it is much nicer if finally a consensus is established about who are the best candidates. Thus consensus may be a good or a bad thing, depending on the application. This section deals with some models of opinion formation where, depending on details, a complete consensus is or is not found .

Imagine that the possible opinions of a large set of people $`i=1,2,\mathrm{},N`$ are described by a real number $`S`$ between zero and unity; it politics this correponds to the traditional left-right classification. Now these people talk to each other and try to convince each other of their opinions. Normally people with vastly different opinions will not agree on any compromise but when the opinions are simular then a compromise is possible. For example, if two people $`i`$ and $`k`$ have similar opinions $`S_i`$ and $`S_k`$ such that their difference $`|S_iS_k|`$ is smaller than some small limit $`ϵ`$ (”bounded confidence”), then a reasonable compromise is that both accept the average $`(S_i+S_k)/2`$ as their new opinions. This rule can be generalized to several people agreeing within $`ϵ`$ with opinion $`S_i`$, or to a vector of several opinion variables per person, instead of merely one $`Si`$ .

With simultaneous updating of 625 initially random opinions,  found dozens of final opinions for $`ϵ=0.01`$, two opinions (”polarization”) for $`ϵ=0.15`$, and one opinion (”consensus”) for $`ϵ=0.25`$. Thus, the more tolerant people are, the higher is the chance for consensus, certainly a plausible result showing that the model is reasonable. Similar results were found by the French group , giving about $`1/(2ϵ)`$ different final opinions in a similar model. In both cases, everybody could interact with everybody, like at a long conference reception where nobody can sit and thus everybody walks around a lot.

In the opposite extreme, people site on a lattice and talk only with their nearest neighbours. And in between is the case of slow diffusion when after every chat with a neighbour everybody tries to move to an empty neighbour place, similar to . Computationally it is simpler to model opinions as discrete integers $`S=1,2,3,\mathrm{}q`$ such that only neighbouring opinions $`S\pm 1`$ can influence opinion $`S`$. Thus geometric space is a square lattice, and opinion space is a one-dimensional chain. The role of the previous parameter $`ϵ`$ is now played by $`1/q`$.

This discrete dynamics was used in particular for the Sznajd model  where, however, opinion $`S_i`$ is not influenced by its neighbours (as in Ising models ), but instead it influence them. Thus in a simple version, one site $`i`$ forces its opinion $`S_i`$ onto all those nearest neighbours on the lattice which have the opinion $`S_i\pm 1`$. Information thus flows from inside out, instead of the usual information flow from outside inwards. The people are now more like missionaries trying to convince others, not negotiators looking for a compromise or opportunists accepting the majority opinion of their neighbours (voter models).

Most of the research on Sznajd models uses, however, the original principle of ”united we stand, divided we fall”. Then two neighbours with the same opinion convince their six neighbours on the square lattice. If the two people in the middle have divided opinions, they do not convince anybody. Similarly, the children obey perhaps their parents of both parents agree with each other; if mother says something different from father, the children are more free to do what they want. This model, which has also been simulated in one and three dimensions as well as on the triangular lattice, always leads to a consensus when $`q=2`$, i.e. when all opinions are similar to each other in the sense of the above $`S\pm 1`$ rule. (If we allow $`q>2`$ and relax the bounded confidence restriction $`S\pm 1`$ by letting a pair convince all neighbours irrepective of the difference in opinion, then again always a consensus is found.) Since the Sznajd model was reviewed in  we now summarize mainly the more recent results.

The difference between simultaneous and random sequential updating of opinions  is quite crucial for the Sznajd model . If simultaneously two pairs of neighbours tell me that I should vote the way they want, and these two pairs have different opinions, then I am frustrated and do not change my opinion. This frustration effect makes a consensus very difficult, one needs a very large initial majority for one opinion to find a consensus for this opinion ($`q=2`$, square lattice). Thus formal committee meetings have less chance of success than informal encounters spread over a longer time interval.

What if the Sznajd agents diffuse slowly on a half-filled square lattice under bounded confidence, i.e. $`S`$ convinces only $`S\pm 1`$, with $`q>2`$ opinions. For $`q=5`$, usually most people at the end will have adopted the centrist opinion 3, some the extremist opinions 1 and 5, but opinions 2 and 4 have died out completely. More interesting are $`q=4`$ opinions: One of the more centrist opinions, like 2, wins over nearly everybody at the end, after having eaten up all neighbouring opinions 1 and 3. A small opposition of opinion 4 usually remains left. Thus if you want to stay with the winner, observe the evolution of votes: The one who is leading half-way through the race will get most votes at the end, the one who is on second place at half-time will get nothing, while the third- ranked opinion has a small set of followers at the end. The fourth-ranked opinion dies out soon. Fig.2 shows two typical examples on a $`101\times 101`$ square lattice, with 4 and 5 opinions. (For clarity the opinions in the $`q=4`$ case are plotted not at 1, 2, 3, 4 but at 1.5, 2.5, 3.5 and 4.5; the 39 votes at 4.5 are barely visible on this scale.)

Thus far opinions were exchanged only among neighbours. There are, however, also mass media which try to influence us through advertising. This can be simulated  by assuming for $`q=2`$ that at every iteration after the exchange with neighbours every person is flipped into opinion 1, if opinion 1 is the one which advertises. Then for large lattices already a small amount of advertising is sufficient to convince the whole square lattice, if initially the two opinions were shared equally often.

Physicists like to have a Hamiltonian (usually meaning just energy) for their models; for Sznajd this was achieved only in one dimension . Physicists also are accustomed that long range forces proportional to 1/distancez facilitate phase transitions. The Sznajd model with interactions only between nearest neighbours on the square lattice has such a phase transition for $`q=2`$: That opinion which initially has a tiny majority, at the end gets all votes. The introduction of long-range forces  does not change this; in fact, $`z=1`$ seems to make the transition less pronounced than $`z=4`$.

Summary:

In the above examples, the whole human being is reduced to a simple number, which represents her opinion (consensus models) or his history (hierarchies). This is, of course, a great simplification, and cognitive scientists may dislike it. Actually, more complicated models of human behaviour have been simulated extensively as neural networks, and one could apply these neural network techniques onto the above models: Does a person recognize before a fight the enemy through associative memory ? How much external influence is needed to move from one intrinsic fixed point of the neural network to another and thus to change opinion (consensus models). The above models ignore these details just as Kepler’s laws how the Earth circles the Sun ignore the whole structure of the Earth. Clearly, the Earth is not point-like, Kepler knew that, and geographers endanger their employment if they treat the Earth as a point mass. But for the purpose of describing celestian motions, the model of a point-like Earth is good and lead to the development of theoretical mechanics by Newton and others later.

Humans, in contrast to the Earth, should have intelligence and thus the above analogy with the Earth may be inappropriate. We make our own decision whether we smoke, drink beer, or go on a diet; and all these decisions may influence our health and age of death. Clearly, we do not make these decisions randomly. Nevertheless, averaged over a large population, experts have constructed life tables, which entered into health and life insurance business and seem to work reasonably even though they assume humans as dying randomly. The first life tables were constructed centuries ago by Halley, for whom a famous ”comet” is named. In a similar sense, the above computer simulations may give us information on averages over many people, but not the fate of one specific person. For example, the Sznajd model was used to simulate the distribution of votes among candidates in Brazilian elections in general, but could not predict how many votes one named candidate in one specific election got.

I hope this small selection of examples encourages the readers to enter this field and to invent their own simulation models. Thanks are due to A. Maksymowicz and K. Kułalkowski of AGH in Kraków, Poland for hospitality when this review was drafted there.