format: text

text:
# 1 Introduction

## 1 Introduction

When a system is in contact with a heat bath, in the most common situation it relaxes toward the equilibrium state rather rapidly, typically in a manner characterized by exponential decay, except at a phase transition point. For instance, a hot particle in cold water generally loses its energy very quickly. However, does such rapid relaxation of a high-energy state hold generally? Consider, for example, a system composed of several elements that is in contact with a heat bath only through its outer part. A protein molecule in water represents such an example, since the hydrophobic part of a protein molecule will be segregated from an aqueous solution . In such a situation, the relaxation process of the internal part naturally depends on the dynamics of the whole system, including the temporal and spatial relationship between the internal and external parts. Indeed, it has been reported that proteins excited into high-energy states, exhibit very slow relaxation, compared with that seen in standard dissipation processes in systems of similar size . Such long-term energy storage may be necessary for proteins to work functionary. For instance, consider a thermal ratchet model  for motor proteins that work in a heat bath. In the model, the existence of both high and low temperature domains in proteins are assumed, to generate a directive motion. Then, for such a mechanism to work, it should be answered how the proteins maintain a high energy (i.e., keep an effectively high temperature state) at a localized domain within them.

In connection with this long-term energy storage in proteins, two of the present authors  have studied a model of a coupled pendulum in partial contact with a heat bath. It was found that when a pendulum that is not in contact with a heat bath is highly excited, the energy dissipates only slowly, approximately with logarithmic decay. This slow relaxation seems to be consistent with the Boltzmann-Jeans conjecture  for very slow relaxation in Hamiltonian systems without contact to a heat bath. In such systems, the relaxation time increases exponentially of some positive power of the injected energy .

When all elements in a Hamiltonian system are in full contact with a heat bath, the mechanism allowing for slow relaxation may be destroyed completely. On the other hand, when a system contacts a heat bath only through an external part, the internal part can evolve according to its ‘own’ dynamics. For this reason, the internal Hamiltonian dynamics may have a strong influence on the relaxation process, and this can lead to a slow relaxation process and long-term energy storage. It is thus important to study how a system in partial contact with a heat bath insulates its internal part from the effect of the heat bath.

To study such behavior of Hamiltonian systems in partial contact with a heat bath, we adopt a simple model, introduced in §2, consisting of two coupled elements: The first element is in direct contacted with a heat bath to allow for the dissipation of the energy, while the second element interacts only with the first element, and not directly with the heat bath. We study relaxation from a far-from-equilibrium state by preparing the second element in a highly excited state with a large kinetic energy. The first element, on the other hand, remains always in near-equilibrium with the heat bath. It should be noted that, without the heat bath, the presently studied Hamiltonian system is always integrable and does not relax, even when prepared in a highly excited state. Therefore the relaxation of the second element to thermal equilibrium is due to the indirect effect of the heat bath. In §3, the relaxation equation for the underdamped limit is obtained, while the case of the overdamped limit is studied in §4. The validity of these equations is confirmed by numerical simulations discussed in §5. Relaxation in the intermediate damping case is examined in §6. The results of numerical experiments for this case are also reported there. From these analytical and computational studies, two kinds of power law are obtained for the underdamped and overdamped cases, respectively. It is found that these power laws are independent of the temperature.

## 2 Model System

To study relaxation of a system with internal degrees of freedom that is in partial contact with a heat bath, we consider a Hamiltonian system with $`2`$ degrees of freedom of the general form

$$H=\frac{p_1^2}{2}+\frac{p_2^2}{2}+U(\theta _1\theta _2),$$
(1)
where $`U(\theta )`$ is a periodic function of period $`2\pi `$ satisfying $`|U|<\mathrm{}`$. Considering the system to be in partial contact with a heat bath, and adopting the standard representation of the heat bath, we take the motion of the elements to be given by the following Langevin equations:

$$\{\begin{array}{ccc}\hfill \dot{p}_1& =& \frac{U}{\theta _1}\gamma p_1+\xi (t),\hfill \\ & & \\ \hfill \dot{p}_2& =& \frac{U}{\theta _2}.\hfill \end{array}$$
(2)
Here the random force $`\xi (t)`$ due to the heat bath consists of a Gaussian white noise satisfying $`\xi (t)=0`$ and $`\xi (t_1)\xi (t_2)=2\gamma T\delta (t_2t_1)`$, with $``$ representing the temporal average. (We use units in which $`k_B=1`$). The parameter $`\gamma `$ is the dissipation coefficient characterizing the dissipation to the heat bath, and $`T`$ is the temperature of the heat bath. With such contact to a heat bath, it is generally expected that a system of this kind will relax to equilibrium. The change of the total energy $`E`$ due to dissipation to the heat bath can be written

$$\dot{E}=\gamma p_1^2+p_1\xi (t),$$
(3)
and we see that $`E`$ remains constant if there is no contact ($`\gamma =\xi =0`$) with the heat bath.

Using the set of variables $`\mathrm{\Theta }\theta _1+\theta _2,\varphi \theta _1\theta _2`$ and $`Pp_1+p_2=\dot{\mathrm{\Theta }}`$, we rewrite the Hamiltonian as

$$H=\frac{P^2}{4}+\frac{\dot{\varphi }^2}{4}+U(\varphi ).$$
(4)
When there is no heat bath (i.e., for $`\gamma =\xi =0`$), the equation of motion is completely separable in these variables: $`\dot{P}=0`$, $`\ddot{\varphi }=2U/\varphi `$. Thus we see that in this case this system is completely integrable, and the total momentum $`P`$ is a conserved quantity.

In the following sections, we study this system analytically to elucidate the relaxational behavior in certain limits. The results of numerical simulations of some models are given in §5 and §6.

## 3 Underdamped Limit

In this section we consider (2) in the underdamped limit, i.e. with a small value of $`\gamma `$. With the presence of this weak damping term, the total energy $`E`$ and total momentum $`P`$, which are conserved quantities for $`\gamma =0`$, become slow variables.

### 3.1 Derivation of relaxation equation

Let us consider the situation in which the system is initially far from equilibrium, with very large $`|p_2||p_1|`$. The rotation of $`\varphi `$, accordingly, is much faster than its change due to dissipation. Then, the terms $`U/\theta _1`$ and $`U/\theta _2`$ in (2) change so fast that their effects on the dynamics of $`p_1`$ and $`p_2`$ are averaged out over a time scale relevant for these slower variables. Here we have to distinguish the orders of characteristic time scales. We have the following relations : (the time scale of $`\xi (t)`$) $``$ (the time scale of rotation $`2\pi |\dot{\varphi }|^1`$) $``$ (the time scale of relaxation $`\gamma ^1`$). Now, we introduce a slower time scale $`\tau `$ satisfying $`2\pi |\dot{\varphi }|^1\tau \gamma ^1`$, and consider averages over this time scale, $`\stackrel{~}{X}(t)\frac{1}{\tau }_t^{t+\tau }Xdt`$. Using variables averaged over this time scale, the mode of fast rotation is eliminated.

The equations of motion for the slow time scale are obtained by integrating (2) from $`t`$ to $`t+\tau `$. This yields

$$\{\begin{array}{ccc}\hfill \frac{\mathrm{d}\stackrel{~}{p}_1}{\mathrm{d}t}& =& \gamma \stackrel{~}{p}_1\stackrel{~}{U}^{}+\stackrel{~}{\xi },\hfill \\ & & \\ \hfill \frac{\mathrm{d}\stackrel{~}{p}_2}{\mathrm{d}t}& =& \stackrel{~}{U}^{},\hfill \end{array}$$
(5)
where $`U^{}=\mathrm{d}U/\mathrm{d}\varphi =U/\theta _1=U/\theta _2`$, $`\stackrel{~}{U}^{}\frac{1}{\tau }_t^{t+\tau }dU(\varphi (t^{}))/d\varphi dt^{}`$ and $`\stackrel{~}{\xi }(t)\stackrel{~}{\xi }(t^{})=2\gamma T\delta (tt^{})`$, with $`\stackrel{~}{\xi }(t)=0`$. The $`\delta `$-function here is defined with respect to the time scale $`\tau `$.

In the limit $`|\dot{\varphi }|\mathrm{}`$, the averaged variables can simply be replaced by their uniform averages over $`\varphi `$, defined generally by $`\overline{X}=_0^{2\pi }X(\varphi )d\varphi /(2\pi )`$. Then we obtain

$$\{\begin{array}{ccc}\hfill \frac{\mathrm{d}\overline{p_1}}{\mathrm{d}t}& =& \gamma \overline{p_1}+\overline{\xi },\hfill \\ & & \\ \hfill \frac{\mathrm{d}\overline{p_2}}{\mathrm{d}t}& =& 0,\hfill \end{array}$$
(6)
where $`\overline{U^{}}=0`$ for periodic potentials. Then, as the lowest-order estimation in this limit, we have

$$\stackrel{~}{p}_1\stackrel{~}{p}_1(0)\mathrm{e}^{\gamma t},\stackrel{~}{p}_2=\stackrel{~}{p}_2(0),$$
(7)
starting in a state far from equilibrium with large values of $`|p_1|`$ and $`|p_2|`$. Beginning in such a state, after a time scale of $`O(\gamma ^1)`$, the first element approaches thermal equilibrium ($`p_1^2T`$), while, within this lowest-order approximation, the second element continues to satisfy $`p_2=p_2(0)`$ and does not evolve toward equilibrium.

To observe the relaxation process of the second element, we need the next order estimation. Since $`\dot{\varphi }`$ is finite and can depend on $`\varphi `$, in general the value of a time-averaged variable $`\stackrel{~}{X}`$ deviates from that of the corresponding phase-averaged variable $`\overline{X}`$. We now derive the correction to the equations of motion that this deviation introduces to first order in $`|\dot{\varphi }|^1`$, and thereby obtain the relaxational equations for the second element and the total energy $`E`$.

After the exponential decay derived by (7) occurring on a time scale $`O(\gamma ^1)`$, $`p_1`$ is approximately equal to its equilibrium value while $`p_2`$ is still large. Then the relations

$$|p_1|\sqrt{T}p_2P$$
(8)
hold, where $`P`$ is the total momentum of the system and is large enough for the fast rotational motion of the second element to be maintained. From (8), we straightforwardly obtain $`\dot{\varphi }P`$. The total energy $`E`$ of the system can then be expressed by the variables $`P`$, $`\varphi `$ and $`p_1`$ as

$$E=\frac{P^2}{2}+U(\varphi )p_1P+p_1^2,$$
(9)
where $`U(\varphi )`$, $`p_1P`$ and $`p_1^2`$ are much smaller than $`P^2/2`$.

Since $`|\dot{\varphi }|`$ is sufficiently large and $`\gamma `$ is sufficiently small, $`E`$ and $`P`$ are almost conserved over the time scale $`\tau `$. To consider the slow evolution of $`E`$ and $`P`$, we eliminate the fast variable $`\varphi `$, using the time averaged variables $`\stackrel{~}{X}`$. Within the slow time scale, the value of $`E`$ can then be expressed as

$$E=\frac{P^2}{2}+\stackrel{~}{U}\stackrel{~}{p}_1P+\stackrel{~}{p}_{1}^{}{}_{}{}^{2},$$
(10)
where $`\frac{1}{\tau }_t^{t+\tau }p_1^2dt^{}`$ is approximated by $`\stackrel{~}{p}_{1}^{}{}_{}{}^{2}`$ on the time scale $`\tau \gamma ^1`$. From the relation between (4) and (10), the value of $`\dot{\varphi }`$ is estimated as

$`\dot{\varphi }`$ $`=`$ $`\sqrt{4(EU(\varphi ))P^2}=P\sqrt{1{\displaystyle \frac{4\stackrel{~}{p}_1}{P}}{\displaystyle \frac{4(\delta U\stackrel{~}{p}_1^2)}{P^2}}}`$ (11)
$`=`$ $`P+2\stackrel{~}{p}_1+{\displaystyle \frac{2\delta U}{P}}+{\displaystyle \frac{4\stackrel{~}{p}_1\delta U}{P^2}}+O({\displaystyle \frac{\delta U^2}{P^3}},{\displaystyle \frac{T\delta U}{P^3}}),`$
where $`\delta UU\stackrel{~}{U}`$. The last expansion is valid if the conditions $`|\delta U|P^2`$ and $`TP^2`$ are satisfied. The momenta of the two elements are then given by

$$p_1=\frac{P+\dot{\varphi }}{2}=\stackrel{~}{p}_1+\frac{\delta U}{P}+\frac{2\stackrel{~}{p}_1\delta U}{P^2}+O(\frac{\delta U^2}{P^3},\frac{T\delta U}{P^3}),$$
(12)
This result reveals that both elements oscillate rapidly around their averages $`\stackrel{~}{p}_1`$ and $`\stackrel{~}{p}_2`$ ($`=P\stackrel{~}{p}_1`$). We see that the amplitude of oscillation becomes smaller as the second element rotates faster (i.e., for larger $`p_2`$ and $`P`$).

In order to obtain the relaxation equation for $`E`$, we substitute (12) into (3). By considering only the case with a sufficiently low temperature, and neglecting terms of $`O(TP^2)`$, we obtain

$$\dot{E}=\frac{\gamma \delta U^2}{P^2}+\frac{2\stackrel{~}{p}_1\delta U}{P}+\frac{\delta U\xi }{P}+\stackrel{~}{p}_1(\gamma \stackrel{~}{p}_1+\xi )+O(\frac{\delta U^2}{P^3},\frac{T}{P^2}).$$
(13)
This equation includes fast variables, while the total energy $`E`$ is a slow variable. It is thus appropriate to average out the fast change in this equation. Doing so, we obtain

$$\frac{\mathrm{d}E}{\mathrm{d}t}=\frac{\gamma }{P^2}\{\frac{1}{\tau }_t^{t+\tau }\delta U^2dt^{}\}+\frac{1}{P}\{\frac{1}{\tau }_t^{t+\tau }\delta U\xi (t^{})dt^{}\}+\frac{\mathrm{d}}{\mathrm{d}t}(\frac{\stackrel{~}{p}_1^2}{2}+\stackrel{~}{U})+O(\frac{\delta U^2}{P^3},\frac{T}{P^2}),$$
(14)
by using (5). To first order in $`|\dot{\varphi }|^1`$, we can replace the time average of $`\delta U^2`$ by the cycle average $`\overline{\delta U^2}`$:

$$\frac{1}{\tau }_t^{t+\tau }\delta U(\varphi (t^{}))^2dt^{}\frac{1}{2\pi }_0^{2\pi }\delta U(\varphi )^2d\varphi \overline{\delta U^2}.$$

Applying this approximation, we obtain

$$\frac{\mathrm{d}E}{\mathrm{d}t}=\frac{\gamma \overline{\delta U^2}}{P^2}+\frac{1}{\tau }_t^{t+\tau }\delta U\xi (t^{})dt^{}+\frac{\mathrm{d}}{\mathrm{d}t}(\frac{\stackrel{~}{p}_1^2}{2}+\stackrel{~}{U})+O(\frac{\delta U^2}{P^3},\frac{T}{P^2}).$$
(15)
Since $`E=\stackrel{~}{U}+\stackrel{~}{p}_1^2/2+\stackrel{~}{p}_2^2/2`$ and $`|\stackrel{~}{p}_2||\stackrel{~}{p}_1|`$ ($`\sqrt{T}`$), (15) can be rewritten in terms of $`\stackrel{~}{p_2}`$:

$$\frac{\mathrm{d}}{\mathrm{d}t}(\frac{\stackrel{~}{p}_2^4}{4})=\gamma \overline{\delta U^2}+\stackrel{~}{p_2}\{\frac{1}{\tau }_t^{t+\tau }\delta U\xi (t^{})dt^{}\}+O(\frac{\delta U^2}{\stackrel{~}{p}_2},T).$$
(16)
This equation describes the relaxation to thermal equilibrium for the second element. To consider only the long-time ($`O(\tau )`$) behavior at the current order of approximation, it can be written as

$$\frac{\mathrm{d}}{\mathrm{d}t}(\frac{p_2^4}{4})=\gamma \overline{\delta U^2}+p_2\delta U\xi (t).$$
(17)
### 3.2 Relaxation at T=0

When $`T=0`$, $`\xi (t)`$ vanishes, and $`\stackrel{~}{p}_1`$ is negligible (or, more precisely speaking, it is $`O(P^3)`$). Then we obtain the following relations from (11) and (12)

$$\dot{\varphi }=P+\frac{2\delta U}{P},p_1=\frac{\delta U}{P},p_2=P\frac{\delta U}{P}.$$

The relaxation equation for the total energy $`E`$, (15), is then simplified as

$$\frac{\mathrm{d}E}{\mathrm{d}t}=\frac{\gamma \overline{\delta U^2}}{2E},$$
(18)
for $`\overline{U}E`$ (that is $`EP^2/2`$). The relaxation equation for the second element is now given by

$$\frac{\mathrm{d}}{\mathrm{d}t}(\frac{p_2^4}{4})=\gamma \overline{\delta U^2}.$$
(19)
These relaxation equations reveal that a coupled system in partial contact with a heat bath exhibits non-exponential, slow relaxation after excitation to a high-energy state. Indeed, from (18) and (19), the form of relaxation is obtained as

$`E^2`$ $`=`$ $`E_0^2\gamma \overline{\delta U^2}t,`$ (20)
$`p_2(t)^4`$ $`=`$ $`p_2(0)^44\gamma \overline{\delta U^2}t,`$ (21)
where $`E_0`$ and $`p_2(0)`$ are the initial total energy and momentum of the second element at $`t=0`$, which is defined as a time at which the first element reaches thermal equilibrium. Thus, the characteristic time $`t_r`$ for the relaxation of $`p_2`$ is given by

$$t_r=\frac{p_2(0)^4}{4\gamma \overline{\delta U^2}}=\frac{E_0^2}{\gamma \overline{\delta U^2}}.$$
(22)
With an increase in the injected energy (or momentum), the relaxation time increases in proportion to $`p_2(0)^4`$. This increase is slower than the $`\mathrm{exp}(p_2(0))`$ form obtained from the Boltzmann-Jeans conjecture , but even here, there is a rapid increase of $`t_r`$ with the injected energy. Comparing the characteristic time (22) to that of the standard dissipation processes exhibited by systems in full contact with a heat bath (where $`t_r\gamma ^1`$), it is confirmed that partial contact indeed results in long-time relaxation behavior.

### 3.3 Relaxation at finite temperature

Using the expression (17), we derive a relaxation equation for finite temperature. The solution of (17) can be written

$$\frac{p_2(t)^4}{4}=\frac{p_2(0)^4}{4}\gamma \overline{\delta U^2}t+_0^tp_2\delta U\xi (t^{})dt^{}$$
(23)
for $`t2\pi |p_2|^1`$. We next take an ensemble average (indicated by $``$) of (23), in order to see the averaged behavior for the long-time profile of the relaxation. Although the ensemble average of the third term, $`p_2\delta U\xi `$, cannot necessarily be decomposed into the product of the average of each quantity, such a decomposition is valid at the current order of approximation, as discussed in Appendix 1. Hence we get

$$\frac{p_2(t)^4}{4}\frac{p_2(0)^4}{4}\gamma \overline{\delta U^2}t+_0^tp_2(t^{})\delta U(\varphi (t^{}))\xi (t^{})dt^{}=\frac{p_2(0)^4}{4}\gamma \overline{\delta U^2}t.$$
(24)
Since $`Ep_2^2/2`$ for $`|p_2||p_1|`$, this result leads to

$$E^2=E_0^2\gamma \overline{\delta U^2}t.$$
(25)
Thus the relaxations of the mean values of $`E^2`$ and $`p_2^4`$ exhibit the same forms as these in the case with $`T=0`$. The validity of the approximations used to obtain these forms have been confirmed numerically, as discussed in §5, where (25) to be valid over a sufficiently long interval.

Furthermore, we can approximate the evolution of the variance of $`E^2`$ (or $`p_2^4`$) around the mean value,

$$(E^2E^2)^2_0^t_0^tp_2(t^{})p_2(t^{\prime \prime })\delta U(t^{})\delta U(t^{\prime \prime })\xi (t^{})\xi (t^{\prime \prime })dt^{}dt^{\prime \prime }=2\gamma T_0^tp_2(t^{})^2\delta U(t^{})^2dt^{},$$
(26)
where we have again used the decomposition discussed in Appendix 1. If the time $`t`$ is not too large, this variance can be approximated by

$$(E^2E^2)^22\gamma Tp_2(0)^2\overline{\delta U^2}t4\gamma TE_0\overline{\delta U^2}t.$$
(27)
This form for the increase of the variance has also be shown to agree with numerical results, as discussed in §5.

From (25), we obtain the mean relaxation time $`t_r`$ as

$$t_r\frac{E_0^2}{\gamma \overline{\delta U^2}}.$$
(28)
We thus find that the value of the temperature has little effect on the nature of the relaxation from a high-energy state, as long as the temperature is sufficiently small.

## 4 Overdamped Limit

When the dissipation is sufficiently large, the inertial force can be neglected. In this overdamped limit, the equations of motion (2) can be written

$$\{\begin{array}{ccc}\hfill p_1& =& \gamma ^1\frac{U}{\theta _1}+\gamma ^1\xi (t)+O(\gamma ^2),\hfill \\ & & \\ \hfill \dot{p}_2& =& \frac{U}{\theta _2}.\hfill \end{array}$$
(29)
Then, for sufficiently large $`\gamma `$, the equation for the energy dissipation is obtained straightforwardly as

$$\dot{E}=\gamma p_1^2+p_1\xi (t)=\gamma ^1U^2+\gamma ^1U^{}\xi (t)+O(\gamma ^2),$$
(30)
where $`U^{}\mathrm{d}U/\mathrm{d}\varphi =U/\theta _1=U/\theta _2`$. Solving this equation, we have

$$E=E_0\gamma ^1_0^tU_{}^{}{}_{}{}^{2}dt^{}+\gamma ^1_0^tU^{}\xi (t^{})dt^{},$$
(31)
with the initial energy $`E_0`$ at $`t=0`$.

We now consider the case $`\gamma P/(2\pi )`$. In other words, that in which the time scale of dissipation from the first element is much faster than the rapid rotation of the second element. Here, we are interested in relaxation from a high-energy state (with very large $`|p_2|`$) taking place over a much longer time scale than the rotation time ($`2\pi P^1`$). In this case, the integration over time can be approximated by the cycle average: $`\frac{1}{t}_0^tU_{}^{}{}_{}{}^{2}dt^{}\frac{1}{2\pi }_0^{2\pi }U_{}^{}{}_{}{}^{2}d\varphi \overline{U^2}.`$ By using this approximation, (31) can be replaced by

$$E=E_0\gamma ^1\overline{U_{}^{}{}_{}{}^{2}}t+\gamma ^1_0^tU^{}\xi (t^{})dt^{}.$$
(32)
An intuitive explanation for this approximation is as follows: In the limit $`P\mathrm{}`$, the change of $`\varphi `$ is expected to be uniform in time. Then relaxation on the slower time scale can be approximated by using the uniform average of $`U^{}(\varphi )`$ over $`\varphi `$, to lowest order in $`\gamma ^1`$ and $`P^1`$.

### 4.1 Relaxation at T=0

For $`T=0`$, there is no random force, and the relaxation is obtained directly as

$$E=E_0\gamma ^1\overline{U_{}^{}{}_{}{}^{2}}t.$$
(33)
Since $`Ep_2^2/2`$, this leads to relaxation of the second element given by

$$p_2(t)^2=p_2(0)^22\gamma ^1\overline{U_{}^{}{}_{}{}^{2}}t,$$
(34)
in which $`p_2(0)`$ is the initial momentum.

From (34), the characteristic time $`t_r`$ for the relaxation is obtained as

$$t_r=\frac{\gamma p_2(0)^2}{2\overline{U^2}}=\frac{\gamma E_0}{\overline{U^2}}.$$
(35)
Thus in the present case, the dependence of $`t_r`$ on $`p_2(0)`$ is different from that in the underdamped case, where $`t_rp_2(0)^4`$. We remark that even in this overdamped case, the relaxation is much slower than in the case of full contact with the heat bath ($`t_r\gamma ^1`$). The crossover between the two limiting cases is discussed on the basis of numerical computations in §6.

### 4.2 Relaxation at finite temperature

Taking an ensemble average of (32), we consider the relaxation of the mean value of $`E`$:

$$EE_0\gamma ^1\overline{U_{}^{}{}_{}{}^{2}}t+\gamma ^1_0^tU^{}\xi (t^{})dt^{}=E_0\gamma ^1\overline{U_{}^{}{}_{}{}^{2}}t.$$
(36)
Here $`U^{}(\varphi )\xi (t^{})`$ has been decomposed as $`U^{}(\varphi )\xi (t^{})`$, as is justified at the present order of approximation. This relation can be intuitively explained as follows. The quantity $`U^{}(\varphi )`$ changes due to the Brownian motion of $`\varphi `$ in a manner that is independent of $`E`$, because the dependence of $`U^{}(\varphi )`$ on $`E`$ appears only at a higher order of $`E^1`$. Equation (36) leads to an equation for the relaxation of the second element:

$$p_2(t)^2p_2(0)^22\gamma ^1\overline{U_{}^{}{}_{}{}^{2}}t.$$
(37)
It is noted that this form is independent of the temperature of the heat bath.

Fluctuations around the average relaxation path, which depend on the temperature, can be straightforwardly estimated from (32) as follows:

$$(EE)^2\gamma ^2_0^t_0^tdt^{}dt^{\prime \prime }U^{}(t^{})U^{}(t^{\prime \prime })\xi (t^{})\xi (t^{\prime \prime })2\gamma ^1T\overline{U_{}^{}{}_{}{}^{2}}t.$$
(38)
Thus, the distribution of $`E`$ becomes broader with time, and this variance is proportional to the temperature $`T`$ of the heat bath. The results represented by (36) and (38) agree well with the numerical results given in §5.

The ensemble average for the relaxation time $`t_r`$ is obtained from eq.(36)

$$t_r\frac{\gamma E_0}{\overline{U^2}}.$$
(39)
This is again independent of the temperature. Thus the discussion given in §4.1 is also valid at finite temperature.

## 5 Numerical Experiments

In this section, we check the validity of our analytical approximations, using results obtained by numerically integrating the Langevin equation (2) directly. In order to demonstrate the generality of our results, we adopted three forms of the potential $`U`$ with $`0U1`$ ($`|\delta U|1`$):

$`U_1(\theta )`$ $`=`$ $`{\displaystyle \frac{1}{2}}(1\mathrm{cos}\theta ),`$
$`U_2(\theta )`$ $`=`$ $`{\displaystyle \frac{1}{4}}\{2(\mathrm{cos}\theta +\mathrm{cos}g\theta )\},`$
$`U_3(\theta )`$ $`=`$ $`{\displaystyle \frac{1}{\mathrm{cosh}(\alpha \mathrm{cos}\theta )}}{\displaystyle \frac{1}{\mathrm{cosh}(\alpha )}}.`$
Here, $`g`$ is the golden mean $`(1+\sqrt{5})/2`$, and the parameter $`\alpha `$ represents the steepness of the potential $`U_3(\theta )`$. $`U_1(\theta )`$ and $`U_3(\theta )`$ are periodic functions with period $`2\pi `$, while $`U_2(\theta )`$ is quasi-periodic with two incommensurate frequencies. $`U_3(\theta )`$ possesses a steeper form as $`\alpha `$ increases, and it approaches a train of $`\delta `$-functions with period $`2\pi `$ as $`\alpha \mathrm{}`$. We set $`\alpha =10`$ to realize a sufficiently steep shape of the potential. The values of $`\overline{\delta U^2}`$ and $`\overline{U^2}`$ are given in Table I for each potential form.

The relaxation at a finite temperature $`T`$ was studied by taking an ensemble average. We integrated (2) for almost $`1000`$ samples with random seeds for each set of parameters, and thereby obtained a distribution $`\rho _t(x)`$ of the variable $`x(t)`$, and thus an ensemble average $`x(t)`$ over this distribution. For every sample in a given ensemble, we used the same initial conditions (given below), and thus the differences among paths for individual samples resulted only from the effect of the random force.

In this section, we study the scaling form for the relaxation of the total energy $`E`$ (not $`p_2`$) by considering variations of $`\gamma `$, $`U(\theta )`$ and $`T`$. For this purpose, we chose an initial conditions at a given value of $`T`$ with initial energy $`E_0`$ as

$$\theta _1(0)=\theta _2(0)=p_1(0)=0,p_2(0)=\sqrt{2E_0T}.$$
(40)
Here, $`E_0`$ is a constant value, independent of $`T`$. By using these initial conditions, the total energy $`E`$ was evaluated with the initial energy $`E_0`$ for a number of temperatures, after a time to allow for the relaxation of the first element. After this relaxation, it is obtained that $`\stackrel{~}{p}_1^2/2T/2`$ and $`\stackrel{~}{p}_2^2/2E_0T/2\overline{U}`$ (i.e., $`E=\stackrel{~}{p}_1^2/2+\stackrel{~}{p}_2^2/2+\overline{U}E_0`$).

### 5.1 Underdamped case

We first give the numerical results for $`T=0`$. Figure 1 displays the energy relaxation for the three kinds of potentials described in Table I. For the underdamped case, the values $`\gamma _u`$ in Table I were used for each potential. The relaxation profile obtained from the numerical experiments is fit well by $`E^2=E_0^2\gamma \overline{\delta U^2}t`$, for $`\delta U^2`$ corresponding to each potential form. Thus the numerical results confirm the analytical expression (20).

As a next step, we examine the relaxation process at finite temperature. In Fig.2(a), ensemble averages $`E^2`$ are displayed as functions of time $`t`$ during the relaxation process. This figure includes numerical results from six different set-ups, using the two potential forms $`U_1`$ and $`U_3`$ each with three different temperature $`T=0.01,0.1`$ and $`0.5`$. For a small time $`t`$, the means $`E^2`$ exhibit the same scaling relation as in the $`T=0`$ case: $`E^2E_0^2\gamma \overline{\delta U^2}t`$. This is consistent with (25). As the time increases, the mean begins to deviate slightly from the line defined by (25). At higher temperature, there is a tendency to faster relaxation.

In these computations, the variance of $`E^2`$ (i.e. $`(E^2E^2)^2`$) was found to increase linearly with time, as is seen in Fig.2(b). This figure includes the variances for the above mentioned six cases, as well as those for two more cases with a different value of $`E_0`$ at $`T=0.01`$ for $`U_1`$ and $`U_3`$. The last two sets were included to examine the effect of changing the energy $`E_0`$. From these results, the variance is shown to satisfy the scaling relation $`(E^2E^2)^2/E_0T4\gamma \overline{\delta U^2}t`$, as expected from the analytical result in §5 (27). Note the slight deviation upward at large time.

### 5.2 Overdamped case

In this subsection, we report the numerical results for the relaxational behavior in the overdamped case. Here the parameter values $`\gamma _o`$ in Table I were used. When $`T=0`$, we found that $`E(t)`$ beginning from a high-energy state relaxes as in Fig.3. By considering a variety of potential forms (see Fig.3), the result $`EE_0\gamma ^1\overline{U^2}t`$ is obtained, confirming the analytical result, (33).

For finite temperature $`T`$, the ensemble average was studied numerically, and the same relaxation form, $`EE_0\gamma ^1\overline{U^2}t`$, is found to hold within the range of temperatures considered, as displayed in Fig.4(a). This figure includes results for six parameter sets: $`U_1`$ and $`U_3`$ each with $`T=0.01,0.1,0.5`$.

The variance of an ensemble exhibits a linear increase with time, as seen in Fig.4(b). From the numerical results, it obeys the scaling form $`(EE)^2/T2\gamma ^1\overline{U^2}t`$, which agrees with the analytical relation (38). We see that there is slightly larger deviation from the expected scaling relation here compared to the small deviation exhibited by the mean value in Fig.4(a).

## 6 Crossover Regime

In this section we study the dissipation of energy in intermediate values of damping between the overdamped and underdamped limits. In the preceding sections, the relaxation of the ensemble average of $`E`$ (or $`E^2`$) was shown to be independent of the temperature. For this reason, we confine our study of the crossover regime to the case $`T=0`$. The result we obtain here is expected to hold at a finite temperature also.

Recall the relaxational behavior at $`T=0`$ obtained for the underdamped and overdamped limit:

$$\frac{\mathrm{d}E}{\mathrm{d}t}=\{\begin{array}{cc}\frac{\gamma \overline{\delta U^2}}{2}E^1,\hfill & P\gamma ,\hfill \\ & \\ \gamma ^1\overline{U^2},\hfill & P\gamma .\hfill \end{array}$$
(41)
Here $`EP^2/2`$ and $`\overline{U}E`$. In Appendix 2, we obtain more general approximations of $`p_1`$ and $`\dot{E}`$, for any values of $`\gamma `$ and $`P`$. Here we give an intuitive description of the change in the relaxational behavior with the decrease of $`P`$ for a fixed value of $`\gamma `$.

The first element is driven by the second element and damped by the frictional effect of the heat bath. When $`\gamma `$ is sufficiently small, and $`P`$ is sufficiently larger than $`O(\gamma )`$, this damping due to friction is negligible, and the dynamics of the first element are approximately determined solely by its interaction with the second element. This leads to $`p_1\delta U/P`$. From this we obtain the upper form in (41). As the momentum $`P`$ decreases, the interaction between the first and the second element becomes stronger, and this is accompanied by the growth of $`p_1`$. Due to this input of energy to the first element, it rapidly dissipates energy to the heat bath. In such a situation, the value of $`p_1`$ cannot be estimated as $`\delta U/P`$, due to the large friction. Instead, it is estimated as $`\gamma ^1U^{}`$ for sufficiently small values of $`P`$. The relaxation equation in this case thus assumes the lower form in (41).

From the above picture, we expect a smooth interpolation between the two relaxation forms in (41). The effect of the heat bath on the first element is comparable to that of the second element when the energy $`E`$ ($`P^2/2`$) satisfies the condition $`\gamma \overline{\delta U^2}E^1/2\gamma ^1\overline{U^2}`$. Taking this into consideration, we define the crossover energy for the relaxation behavior,

$$E_c\frac{\gamma ^2\overline{\delta U^2}}{2\overline{U^2}}.$$
(42)
This represents the energy at which the two relaxational forms are equal. By introducing $`t_c`$ defined by

$$t_c\frac{\gamma ^3\overline{\delta U^2}}{2\overline{U^2}^2},$$
(43)
the normalized energy $`\epsilon E/E_c`$, and the normalized time $`\tau t/t_c`$, (41) can be written as

$$\frac{\mathrm{d}\epsilon }{\mathrm{d}\tau }=\{\begin{array}{cc}2\epsilon ^1,\hfill & \epsilon 1,\hfill \\ & \\ 1,\hfill & \epsilon 1.\hfill \end{array}$$
(44)
These equations include neither $`\gamma `$ nor $`U`$ explicitly. Hence, the value of $`\mathrm{d}\epsilon /\mathrm{d}\tau `$ is expected to depend only on the value of $`\epsilon `$, regardless of the value of $`\gamma `$ and the shape of $`U`$.

The relaxation form $`\epsilon `$ obtained numerically is displayed in Fig.5 for the case $`T=0`$. To obtain this figure, we computed (2) at $`T=0`$ for various values of $`\gamma `$ and plotted the relation between $`\epsilon `$ and $`\mathrm{d}\epsilon /\mathrm{d}\tau `$ after the first element reached thermal equilibrium. We see that the value of $`\mathrm{d}\epsilon /\mathrm{d}\tau `$ becomes proportional to $`\epsilon ^1`$ as $`\epsilon `$ increases beyond $`1`$, while $`\mathrm{d}\epsilon /\mathrm{d}\tau `$ approaches $`1`$ and becomes independent of the value of $`\epsilon `$ for smaller $`\epsilon `$, as expressed by (44). This figure shows that the relaxation form changes smoothly between the two limiting regimes. We also remark that the nature of potential form $`U(\theta )`$ seems to affect the sharpness of the crossover near $`E=E_c`$.

## 7 Summary and Discussion

In this paper, we have studied the relaxational behavior from a far-from-equilibrium (highly excited) state for a system consisting of internal and external elements in which only the external element is in contact with a heat bath. In spite of the simplicity of the model, our system exhibits a power-law-type slow relaxation when the potential is bounded, i.e., for $`|U|<\mathrm{}`$. This slow relaxation is observed when the inner part of the system is excited to a high-energy state, with energy much greater than $`|\delta U|`$. The slow relaxation of the internal part results from its indirect contact with the heat bath. In spite of the system’s contact with the heat bath, guaranteeing its approach to equilibrium, the observed slow relaxation is in a strong contrast with the rapid exponential decay exhibited by systems in full contact with a heat bath. We therefore find that a Hamiltonian system in partial contact with a heat bath can maintain its own internal dynamics for a much longer time than one in full contact with a heat bath. If a system has an internal part further separated from the heat bath, the Hamiltonian dynamics of this part can be preserved much longer.

Suppose the momentum of the second element $`p_2`$ is constant in time. Then the present model becomes $`\ddot{\varphi }+\gamma \dot{\varphi }+2U^{}(\varphi )=\gamma p_2+\xi (t)`$. The Brownian motion of a particle in a periodic potential $`U(\varphi )`$ in thermal equilibrium has been analyzed using the Fokker-Planck equation for this model  in the contexts of Josephson tunneling junction , super-ionic conductors , and two-mode lasers . In contrast, we are interested in the autonomous relaxation of $`p_2`$ from a highly excited state, through the internal dynamics obeying (2). This autonomous change of $`p_2`$, due to the internal dynamics of the composite system, brings about slow relaxation characterized by power laws, even though the system is neither subject to special conditions, such as being placed at a phase transition point, nor exists in a glassy state involving many degrees of freedom.

While our interest in Hamiltonian systems in partial contact with a heat bath originates in our desire to gain an understanding of the long-term maintenance of high-energy states, such systems have also been studied in relationship with thermal conduction , for example, in the situation that two ends of a chain are in contact with two heat baths at different temperatures. According to several simulations performed on systems of this type, the Fourier law of thermal conduction is not always obeyed. It appears that certain conditions concerning the internal dynamics, as well as a suitable relation of the system with the heat baths, are required in addition to ergodicity (or mixing), in order for the Fourier law to be obeyed. These studies also suggest that parts of systems in indirect contact with heat baths can maintain over extended periods of time the characteristic dynamics of isolated Hamiltonian systems.

In a class of Hamiltonian systems in the absence of a heat bath, very slow relaxation to equilibrium was suggested in the form of the ‘Boltzmann-Jeans conjecture’ , and such behavior was reported in Refs. . Long-term storage of a large amount of energy within a degree of freedom through this mechanism has also been studied in a Hamiltonian system with chaotic behavior , in which localized energy is transferred to other degrees of freedom only after a very long time. In that case, the energy dissipates logarithmically in time. We have reported in this paper that situations exist in which localized energy in a system dissipates slowly, according to power laws, and that long-term energy storage is possible even if the system is integrable and in contact with a heat bath. From this point of view, our result can be regarded as an extension of the Boltzmann-Jeans conjecture to include dissipation, while the relaxational behavior itself is modified from the logarithmic form of the Boltzmann-Jeans conjecture. As mentioned in the Introduction, long-term energy storage has been reported in some protein molecules , and the relevance of this storage to molecular motors (for muscle contraction) has been discussed. To apply the present ideas regarding slow relaxation to proteins, we need to examine relaxation from a high-energy states for a system with more degrees of freedom. This problem will be studied in a future publication . The main question to be addressed is whether the relaxational behavior exhibits a power-law-form or a logarithmic form when Hamiltonian systems more complex than that considered in the present study are in partial contact with a heat bath.

This work is supported by Grants-in-Aid for Scientific Research from the Ministry of Education, Science and Culture of Japan (11CE2006, 11837004 and 11740218), and the REIMEI Research Resources of Japan Atomic Energy Research Institute.

## 8 Appendix 1

In this appendix, an expression for the ensemble average $`p_2^4`$ is derived from (17). Since this equation is a Langevin equation with multiplicative noise, the drift term of the corresponding Fokker-Planck equation is not necessarily $`\gamma \overline{\delta U^2}`$. We carry out the Kramers-Moyal expansion  for the evolution of the variable $`zp_2^4`$ to obtain the Fokker-Planck equation for the probability $`P(z,t)`$ in the form $`P(z,t)/t=D^{(1)}P/z+(1/2)/zD^{(2)}/zP`$. Computing straightforwardly from (17), we obtain the drift term

$$D^{(1)}=\gamma \overline{\delta U^2}+\frac{\delta U^2\gamma T}{\sqrt{2}p_2^2}+o(\frac{T}{p_2^2})$$

and the diffusion term

$$D^{(2)}=p_2^2\delta U^2\gamma T+o(\frac{T}{p_2^2}),$$

where the second term of $`D^{(1)}`$ comes from the multiplicative noise.

Since we confine our study to the situation with sufficiently large $`p_2`$ and sufficiently small $`T`$, this second term of $`D^{(1)}`$ can be ignored. Thus we have $`D^{(1)}=\gamma \overline{\delta U^2}`$ to the order of approximation of the current analysis. Hence, the ensemble average $`p_2\delta U\xi `$ can be decomposed as $`p_2\delta U\xi `$.

## 9 Appendix 2

In this appendix, we give another derivation of the relaxation equation and a more precise expression of the crossover condition at $`T=0`$. Solving the equation of motion for the first element ( $`\dot{p}_1=\mathrm{d}U(\varphi (t))/\mathrm{d}\varphi \gamma p_1`$), we obtain

$$p_1(t)=_0^{\mathrm{}}d\tau \frac{\mathrm{d}U(\varphi (t\tau ))}{\mathrm{d}\varphi }\mathrm{e}^{\gamma \tau }.$$
(45)
We write the periodic potential $`U(\varphi )`$ as a Fourier series:

$$\{\begin{array}{ccc}U\hfill & =\hfill & \underset{j}{}a_j\mathrm{cos}(\omega _j(\varphi \varphi _j)),\hfill \\ U^{}\hfill & =\hfill & \underset{j}{}a_j\omega _j\mathrm{sin}(\omega _j(\varphi \varphi _j)),\hfill \end{array}$$
(46)
where we assume $`\overline{U}=0`$ without loss of generality. Then, assuming that $`\varphi (t)=\dot{\varphi }t`$, (45) can be written,

$$p_1=\underset{j}{}a_j\omega _j_0^{\mathrm{}}d\tau \mathrm{sin}(\omega _j\dot{\varphi }(t\tau )+\varphi _j)\mathrm{e}^{\gamma \tau }=\underset{j}{}\frac{\dot{\varphi }\omega _j^2U_j+\gamma U_j^{}}{(\dot{\varphi }\omega _j)^2+\gamma ^2}.$$
(47)
Here $`U_j`$ and $`U_j^{}`$ are the $`j`$th components of the expressions in (46). This solution implies that

$$p_1=\dot{\varphi }^1U,\dot{E}=\gamma \dot{\varphi }^2\overline{U^2}(\dot{\varphi }\omega _j\gamma \text{for all}j),$$

and

$$p_1=\gamma ^1U,\dot{E}=\gamma ^1\overline{U^2}(\dot{\varphi }\omega _j\gamma \text{for all}j).$$

These equations are equivalent to those obtained in §6, but in the presents derivation, the crossover condition is more explicitly represented as $`\dot{\varphi }\omega _j\gamma `$. Since $`\overline{U_j^2}=\omega _j^2\overline{U_j^2}`$, this crossover condition can be expressed as

$$\dot{\varphi }^2\gamma ^2\frac{\overline{U_j^2}}{\overline{U_j^2}},$$

which is consistent with the condition discussed in §6.


text:
# Disordering Transitions in Vortex Matter: Peak Effect and Phase Diagram

\[

## Abstract

Using numerical simulations of magnetically interacting vortices in disordered layered superconductors we obtain the static vortex phase diagram as a function of magnetic field and temperature. For increasing field or temperature, we find a transition from ordered straight vortices to disordered decoupled vortices. This transition is associated with a peak effect in the critical current. For samples with increasing disorder strength the field at which the decoupling occurs decreases. Long range, nonlinear interactions in the c-axis are required to observe the effect.

\]

Vortex matter in superconductors exhibits a remarkably rich variety of distinct phases due to the numerous competing interactions . These phases can strongly affect the response of the system, such as the critical current $`J_c`$ and magnetization. In highly anisotropic superconductors, such as BSCCO, a pronounced fish-tail or peak effect, in which $`J_c`$ shows a sharp increase, is observed as a function of increasing field . This peak can be interpreted as occurring when the increasing field weakens the interlayer coupling of vortex pancakes due to geometric constraints, and a transition occurs from weakly pinned 3D line vortices to decoupled 2D pancake vortices which can more easily adjust their positions to maximize the pinning . The peak has also been proposed to arise from plasticity, proliferation of in-plane defects, dynamical effects, or matching effects .

There is mounting experimental evidence that the peak effect is associated with a sharp transition in the vortex lattice from an ordered state to a disordered state. In BSCCO, neutron scattering  and muon lifetime  experiments provide evidence that a transition from an ordered 3D vortex arrangement to a disordered or decoupled arrangement is associated with the peak effect. Additional evidence from plasma-resonance  and magneto-optical studies  point to the first-order nature of this transition. In YBCO a rapid increase in $`J_c`$ as a function of magnetic field is observed and is thought to indicate a transition from an ordered to a disordered state. History and memory effects found near this peak indicate that this transition is also first order , suggesting that the physics of the second peak is similar in BSCCO and YBCO. In recent muon-spin measurements in YBCO, the order-disorder transition has been interpreted to be a 3D-2D transition of the vortex lattice . As a function of temperature, in YBCO a peak effect can be observed very near or at the melting line . Transformer measurements in this regime provide evidence of vortex cutting in the liquid state, suggesting that the breakdown of the 3D nature of vortex lines is relevant in this case as well.

A key question is what type of sharp transition is responsible for the observed peak in $`J_c`$, and whether the mechanism of the peak effect is the same going through the melting line as through the second peak. In this paper we argue that it is a decoupling transition along the vortex line, in the $`c`$-direction, which is responsible. We present results from a simulation of magnetically interacting vortices in which we demonstrate a robust peak in $`J_c`$ as a function of magnetic field and produce an H-T phase diagram in good agreement with experiment.

We consider a 3D layered superconducting model material containing an equal number of pancake vortices in each layer, interacting magnetically. This simulation differs from those used previously by several other groups in that it treats the long-range interactions along the length of the vortex line exactly, and neglects the Josephson coupling, approximating highly anisotropic materials. Previous models have treated the inter-plane interactions as simple (linear) elastic connections .

The overdamped equation of motion for vortex $`i`$ is $`𝐟_i=_{j=1}^{N_v}_i𝐔(\rho _{ij},z_{ij})+𝐟_i^{vp}+𝐟_d+𝐟^T=𝐯_i`$. The total number of pancakes per layer is $`N_v`$, and $`\rho _{ij}`$ and $`z_{ij}`$ are the distance between vortex $`i`$ and vortex $`j`$ in cylindrical coordinates. We impose periodic boundary conditions in $`x`$ and $`y`$ directions . The interaction energy between pancakes is

$`𝐔(\rho _{ij},0)=2dϵ_0\left((1{\displaystyle \frac{d}{2\lambda }})\mathrm{ln}{\displaystyle \frac{R}{\rho }}+{\displaystyle \frac{d}{2\lambda }}E_1\right),`$ (1)
$`𝐔(\rho _{ij},z)=s_m{\displaystyle \frac{d^2ϵ_0}{\lambda }}\left(\mathrm{exp}(z/\lambda )\mathrm{ln}{\displaystyle \frac{R}{\rho }}+E_2\right),`$ (2)
where $`R`$ is the maximum in-plane distance, $`E_1=_\rho ^{\mathrm{}}𝑑\rho ^{}\mathrm{exp}(\rho ^{}/\lambda )/\rho ^{}`$, $`E_2=_\rho ^{\mathrm{}}𝑑\rho ^{}\mathrm{exp}(\sqrt{z^2+\rho ^2}/\lambda )/\rho ^{}`$, $`ϵ_0=\mathrm{\Phi }_0^2/(4\pi \lambda )^2`$, $`d=0.005\lambda `$ is the interlayer spacing, and $`\lambda `$ is the London penetration depth. We model the pinning as $`N_p`$ short range attractive parabolic traps that are randomly distributed in each layer. The pinning interaction is $`𝐟_i^{vp}=_{k=1}^{N_p}(f_p/\xi _p)(𝐫_i𝐫_k^{(p)})\mathrm{\Theta }((\xi _p|𝐫_i𝐫_k^{(p)}|)/\lambda )`$, where the pin radius is $`\xi _p`$, the pinning force is $`f_p`$, and $`f_0^{}=ϵ_0/\lambda `$. Thermal fluctuations are represented as Gaussian noise of width $`f_T`$, with $`<f_i^T(t)>=0`$ and $`<f_i^T(t)f_j^T(t^{})>=2\eta k_BT\delta _{ij}\delta (tt^{})`$. The parameter $`s_m`$ is used to vary the coupling strength between planes, but except where noted it is set to $`s_m=1.0`$. To vary the applied magnetic field $`H`$, we fix the number of vortices in the system and change the system size (vortex density $`n_v`$). The pin density remains fixed, but in all cases $`N_p>N_v`$. We consider systems of $`L=8`$ to 32 layers containing from 1 to 80 vortex pancakes per layer. The system size ranges from $`0.6\lambda \times 0.6\lambda `$ to $`141\lambda \times 141\lambda `$.

The pancakes in our model can behave in two possible ways. They may align along the c-axis into well-defined vortex lines (“3D”). Alternatively, the pancakes may break apart in the c-direction and move independently in each plane (“2D”). We can cross between the two types of behavior in the presence of disorder by varying the inter-plane coupling strength, $`s_m`$. The vortices are 3D at high coupling, while 2D at low coupling. We find that due to the long range interactions and the nonlinearity, the transition between these two states is sharp and has many first-order characteristics, including strong hysteresis as well as superheating and supercooling effects. We have studied this coupling strength crossover in detail in previous papers , and have shown that it is associated with a large change in the critical current $`f_c`$, as indicated in Fig. 1(a). To quantify the alignment of the vortices between layers, we also plot the correlation function in the $`z`$-direction, $`C_z=1\mathrm{\Theta }(a_0/2|𝐫_{i,L}𝐫_{j,L+1}|)|𝐫_{i,L}𝐫_{j,L+1}|2/a_0`$, where $`a_0`$ is the vortex lattice constant. When $`C_z=1.0`$, the pancakes are aligned into 3D lines, whereas a low value of $`C_z`$ indicates that the pancakes are decoupled.

It is expected that changing the magnetic field in the system will cause a decoupling transition for the same reason that changing the coupling strength for interlayer interactions does. This is because the vortices come closer together in plane as the field increases while the spacing between planes remains constant. Thus, at higher fields the relative coupling between planes becomes weaker. Here, we show that the same decoupling transition previously observed as a function of coupling strength also occurs as a function of magnetic field.

In Fig. 1(b-d) we plot $`f_c`$ as a function of vortex density (magnetic field) $`n_v`$, at a fixed temperature of $`f_T=0.0005`$ for several different samples containing different numbers of vortices. A clear peak or fish-tail effect is present. In Fig. 1(c), a sample containing $`N_v=4`$ is shown. Here, for low fields ($`n_v<0.1`$), $`f_c`$ is high and the vortices are uncorrelated in the $`z`$ direction as well as disordered in plane. At these low fields the vortex-pin interactions dominate. $`f_c`$ drops rapidly with field in this regime due to the

increasing in-plane vortex lattice stiffness, which causes the pinning to be less effective . For intermediate fields, $`0.1<n_v<0.2`$, $`f_c`$ is low and the vortices form an ordered 3D structure as indicated by a near unity correlation, $`C_z1`$. The 3D vortex lines in this regime are poorly pinned by the point pinning. At $`n_v=0.1`$ there is a sharp decoupling transition as is apparent in the abrupt drop of $`C_z`$. Simultaneously, $`f_c`$ increases rapidly to $`f_c=0.8f_p`$. At higher fields $`f_c`$ gradually decreases due to the increasing in-plane interactions, which continue to stiffen the vortex lattice and thereby weaken the effectiveness of the pinning. We have observed similar peaks in the critical current for samples containing $`N_v=12`$, $`30`$, and $`80`$ vortices. Data for the $`N_v=80`$ sample is shown in Fig. 1(d). The exact value of $`n_v`$ at which the peak falls is affected by the pinning strength and can be adjusted within a wide range, as indicated in Fig. 1.

To show that the peak effect observed here arises from the onset of plasticity of the vortices between layers, and not from in-plane plasticity, we conduct simulations of a rigid vortex lattice, represented by a single vortex with periodic boundary conditions. In the $`N_v=1`$ simulation in-plane defects are not permitted. Thus, the vortex line can disorder along the c-axis only. The fact that we can observe the second peak in a sample containing only a rigid lattice, as in Fig. 1(b), indicates that the observed peak is not being generated by the proliferation of defects in the plane. This, along with our previous results on the existence of a transition as a function of $`c`$-axis coupling strength, provides evidence that it is a change in $`c`$-axis correlation that causes striking changes in $`f_c`$.

We can now probe the effect of temperature on the

observed transition. As illustrated in Fig. 2(a-b), the field at which the 3D-2D transition occurs remains almost constant with $`T`$, in agreement with experiment . $`f_c`$ in the decoupled state drops as the temperature is raised, but the jump in $`f_c`$ persists. In Fig. 2(c) we show the peak effect in the average vortex velocity in the direction of drive, $`V_x=<v_x>`$, as a function of temperature for fixed values of $`n_v`$. Starting at $`f_T=0`$, we apply a constant driving force of $`f_d=0.1f_0^{}`$ such that the vortices are moving and ordered. As $`f_T`$ is increased there is a transition to a 2D decoupled vortex arrangement as seen in the drop in $`C_z`$. This coincides with a sharp drop in the vortex velocities $`V_x`$ since the decoupled vortices experience stronger effective pinning. As $`f_T`$ is further increased, the effectiveness of the pinning gradually decreases and $`V_x`$ gradually increases back to the free flow level. The $`f_c`$ measurements at different temperatures and the resistivity measurements as a function of temperature exhibit the same features as experimental data taken near the peak effect regime.

Next, we construct a phase diagram for our model in the H-T plane. Since we wish to stress that all qualitative features of the phase diagram result from the c-axis interactions, we consider a system containing a rigid vortex lattice. We first locate the thermal melting line for $`N_v=1`$ in the absence of pinning. This is shown as open circles on the phase diagram. We find a reentrance in the melting, shown in more detail in the top inset of Fig. 3, which agrees with theoretical predictions for the melting . In the presence of pinning, we measure $`f_c`$ as a function of field (as in Fig. 1) at several different temperatures. We plot the 3D-2D transition field as solid squares. As was shown in Fig. 2, the location of the transition is insensitive to temperature at low

temperature. We find, however, that the temperature weakens the pinning enough that for $`f_T>1`$ a transition in $`f_c`$ is no longer observed. Instead $`f_c`$ is small everywhere. The transition is still observable in $`C_z`$, and data taken from $`C_z`$ is shown above $`f_T=1`$, indicating that the transition begins to drop in field at higher temperatures until it merges with the clean melting line.

The field at which the second peak occurs is lowered in samples with stronger pinning, as indicated in Fig. 3 where the second peak line for a sample with $`N_v=1`$ and stronger pinning is plotted. This is consistent with experiments in which the peak shifts to lower fields in samples that have had their pinning strength increased artificially.

We have also constructed a phase diagram for a sample with $`N_v=4`$ and much stronger pinning, shown in the lower inset of Fig. 3. In the disordered system the 3D-2D transition is relatively flat in $`f_T`$ for $`f_T<0.3`$ and then begins to turn down as the clean melting line is approached. We also observe that the clean system shows a reentrant behavior for low fields due to the reduced vortex interactions. In the disordered sample this reentrant line is pushed up in field. We see reentrance of the second peak line in good agreement with recent measurements . In still larger systems containing $`N_v=80`$ vortices, we have taken two slices across the phase diagram as a function of field and temperature, shown in Fig. 1(d) and

the inset of Fig. 2(c), respectively, and find decoupling transitions consistent with those shown above.

To show that it is the long range interactions and the nonlinearity in the c-axis that cause the sharpness of the transition, we replace the long range interplane interactions by a model in which the pancakes are coupled by elastic nearest-neighbor springs. This model resembles one considered earlier . The interlayer interaction is $`𝐟_{i,j}=s_m(𝐫_{i,L+1}2𝐫_{i,L}+𝐫_{i,L1})`$, where $`s_m`$ can be used to vary the coupling strength.

As shown in Fig. 4, for varying numbers of vortices, we find a smooth change in depinning force for the same parameters and over the entire range of $`n_v`$ considered in the long range interacting case. For $`N_v=1`$, the interlayer interaction of the spring is not affected by varying $`n_v`$ so we vary $`s_m`$ instead. For $`N_v=4`$, 12, and 80, we set $`s_m=1`$ and vary $`n_v`$ directly. In each case we find a smooth decrease of $`f_c`$ with field. We do not observe any indications that a model containing springs and lacking long range interaction and nonlinearity in the z-direction can show a statistically significant peak in the critical current, as was recently suggested .

In conclusion, using numerical simulations of magnetically interacting vortex pancakes in three dimensional layered superconductors, we have obtained static vortex phase diagrams as a function of magnetic field and temperature. We have demonstrated a unique relationship between c-axis vortex correlations and the critical current $`J_c`$. Specifically, we have identified c-axis correlation transitions to be responsible for the peaks in $`J_c`$ as well as the fish-tail effect in magnetization measurements. Our simulations of rigid in-plane vortex lattices, which also exhibit the peak and fish-tail phenomena, support the claim that changes in c-axis (as opposed to in-plane) correlations are responsible for the observed anomalies. Finally, we have addressed recent suggestions that the peak effect in $`J_c`$ can be observed in a model system where inter-plane interactions are modeled by linear springs between pancakes. We only observe the critical current peak when long-range interactions are present, again suggesting that the peak and fish-tail effect are directly related to the nonlinear and nonlocal nature of the inter-plane vortex interactions that lead to decoupling transitions.

We acknowledge useful discussions with D. Domínguez, A. Kolton, A. Koshelev, W.K. Kwok, V. Vinokur, and E. Zeldov. This work was supported by CLC and CULAR (LANL/UC), by NSF-DMR-9985978, and by the Director, Office of Adv. Scientific Comp. Res., Div. of Math., Information and Comp. Sciences, U.S. DoE contract DE-AC03-76SF00098.


text:
# Baryonic and Non-Baryonic Dark Matter

### Abstract

Cosmological nucleosynthesis calculations imply that there should be both non-baryonic and baryonic dark matter. Recent data suggest that some of the non-baryonic dark matter must be “hot” (i.e. massive neutrinos) and there may also be evidence for “cold” dark matter (i.e. WIMPs). If the baryonic dark matter resides in galactic halos, it is likely to be in the form of compact objects (i.e. MACHOs) and these would probably be the remnants of a first generation of pregalactic or protogalactic Population III stars. Many candidates have been proposed - brown dwarfs, red dwarfs, white dwarfs or black holes - and at various times each of these has been in vogue. We review the many types of observations which can be used to constrain or exclude both baryonic and non-baryonic dark matter candidates.

### 1. Introduction

Evidence for dark matter has been claimed in many different contexts. There may be local dark matter in the Galactic disc, dark matter in the halos of our own and other galaxies, dark matter associated with clusters of galaxies and finally - if one believes that the total cosmological density has the critical value - smoothly distributed background dark matter. Since dark matter probably takes as many different forms as visible matter, it would be simplistic to assume that all these dark matter problems have a single solution. The local dark matter is almost certainly unrelated to the other ones and, while the halo and cluster dark matter would be naturally connected in many evolutionary scenarios, there is a growing tendency to regard the unclustered background dark matter as different from the clustered component. This is because the latest supernovae measurements indicate that the cosmological expansion is accelerating, which means that the total density must be dominated by some form of “exotic” energy with negative pressure (possibly a cosmological constant).

The combination of the supernovae and microwave background observations suggest that the density parameters associated with the exotic (unclustered) component and the ordinary (positive pressure) matter component are $`\mathrm{\Omega }_X=0.8\pm 0.2`$ and $`\mathrm{\Omega }_M=0.4\pm 0.1`$ respectively (Turner 1999). This is compatible with the total density parameter being $`1`$, as expected in the inflationary scenario, but does not definitely require it. This paper will focus exclusively on the matter component of the cosmological density, with particular emphasis on the question of how much of this is baryonic and non-baryonic. We will also be interested in the relative contributions of hot and cold non-baryonic matter. As emphasized by Turner (1999), each of these components seem to be required and it is remarkable that all their densities are within one or two orders of magnitude of each other. Why this should be remains a mystery but we will not explore this further here.

### 2. Evidence for Baryonic and Non-Baryonic Dark Matter

The main argument for both baryonic and non-baryonic dark matter comes from Big Bang nucleosynthesis calculations. This is because the success of the standard picture in explaining the primordial light element abundances only applies if the baryon density parameter lies in the range (Copi et al. 1995)

$$0.007h^2<\mathrm{\Omega }_B<0.022h^2$$
(1)
where $`h`$ is the Hubble parameter in units of $`100`$ km s-1 Mpc-1. A few years ago the anomalously high deuterium abundance measured in some intergalactic Lyman-$`\alpha `$ clouds suggested that the nucleosynthesis value for $`\mathrm{\Omega }_B`$ could be lower than indicated by eqn (1) (Carswell et al. 1994, Songaila et al. 1994, Rugers & Hogan 1996, Webb at al. 1997). However, the evidence for this has always been strongly disputed (Tytler et al. 1996) and recent studies of the quasars Q1937 and Q1009 suggest that the deuterium abundance is $`3.3\times 10^5`$ (Burles & Tytler 1998). This corresponds to the range

$$0.018h^2<\mathrm{\Omega }_B<0.020h^2,$$
(2)
which is much tighter than indicated by eqn (1) (Burles et al. 1999).

Both eqns (1) and (2) suggest that the upper limit on $`\mathrm{\Omega }_B`$ is well below 1, which implies that no baryonic candidate could provide the matter density required by large-scale structure observations. This conclusion also applies if one invokes inhomogeneous nucleosynthesis since one requires $`\mathrm{\Omega }_B<0.09h^2`$ even in this case (Mathews et al. 1993). On the other hand, the upper limit on $`\mathrm{\Omega }_B`$ allowed by eqns (1) and (2) almost certainly exceeds the density of visible baryons $`\mathrm{\Omega }_V`$. A careful inventory by Persic & Salucci (1992) shows that the density in galaxies and cluster gas is

$$\mathrm{\Omega }_V(2.2+0.6h^{1.5})\times 10^30.003$$
(3)
where the last approximation applies for reasonable values of $`h`$. A recent review by Hogan (1999) has updated this estimate of $`\mathrm{\Omega }_V`$: he replaces the factors in brackets in eqn (3) by 2.6 for spheroid stars, 0.86 for disc stars, 0.33 for neutral atomic gas, 0.30 for molecular gas and 2.6 for the ionized gas in clusters (assuming $`H_o=70`$ km s-1 Mpc-1). Although this is still well below the nucelosynthesis bound, this may not account for all components. Some of the baryons may be in an intergalactic medium; for example, Cen & Ostriker (1999) find that - in the context of the “Cold Dark Matter” scenario - half the mass of baryons should be in warm ($`10^510^7`$K) intergalactic gas at the present epoch. Another possibility is that the missing baryons could be in gas in groups of clusters. This has been emphasized by Fukugita et al. (1998), who argue that plasma in groups (which would be too cool to have appreciable X-ray emission) could provide all of the cosmological nucleosynthesis shortfall. Indeed the review by Hogan (1999) suggests that the ionized gas in groups could have a density parameter as high as $`0.014`$. However, it must be stressed that this estimate is based on an extrapolation of observations of rich clusters and there is no direct evidence for this.

Which of the dark matter problems identified in Section 1 would be baryonic or non-baryonic? We have already seen that the dark matter in clusters would need to be mainly non-baryonic. On the other hand, the dark matter in galactic discs could certainly be baryonic: even if all discs have the 60% dark component envisaged for the Galaxy by Bahcall et al. (1992), this only corresponds to $`\mathrm{\Omega }_d0.001`$, well below the nucleosynthesis bound. Indeed the disc dark matter would need to be baryonic in order to dissipate enough to form a disc. However, the Bahcall et al. claim is strongly disputed by Kuijken & Gilmore (1989) and Flynn & Fuchs (1995), In fact, recent Hipparcos observations suggest that the dark disc fraction is below 10% (Crézé et al. 1998).

The more interesting question is whether the halo dark matter could be baryonic. If the Milky Way is typical, the density associated with halos would be

$$\mathrm{\Omega }_h0.03h^1(R_h/100\text{kpc}),$$
(4)
where $`R_h`$ is the (uncertain) halo radius, so the upper limit in eqn (1) implies that all the dark matter in halos could be baryonic providing $`R_h<70h^1`$kpc. This is marginally possible for our galaxy (Fich & Tremaine 1991), in which case the dark baryons could be contained in the remnants of a first generation of “Population III” stars (Carr 1994). This corresponds to the “Massive Compact Halo Object” or “MACHO” scenario and has attracted considerable interest as a result of the LMC microlensing observations. Even if halos are larger than $`70h^1`$kpc, and studies of the kinematics of other spirals by Zaritsky et al. (1993) suggest that they could be as large as $`200h^1`$kpc, a considerable fraction of their mass could still be in stellar remnants.

On theoretial grounds one would expect the halo dark matter to be a mixture of WIMPs and MACHOs. For since the cluster dark matter must be predominantly cold, one would expect at least some of it to clump into galactic halos. The relative densities of the two components would depend sensitively on the formation epoch of the Population III stars. If they formed pregalactically, one would expect the halo ratio to reflect the cosmological background ratio (viz. $`\mathrm{\Omega }_B/\mathrm{\Omega }_M0.1`$). However, if they formed protogalactically, the ratio could be larger since the baryons could have dissipated and become more centrally concentrated than the WIMPs. In order to distinguish between the pregalactic and protogalactic scenarios, it is important to gain independence evidence about the formation epoch of the putative MACHOs. At moderate redshifts one can obtain a lower limit to the baryon density by studying Lyman-$`\alpha `$ clouds. The simulations of Weinberg et al. (1997) suggest that the density parameter of the clouds must be at least $`0.017h^2`$ at a redshift of 2 and this is already close to the upper limit given by eqn (1). By today some of these baryons might have been transformed into a hot intergalactic medium or stars but this suggests that there was little room for any dark baryons before $`z=2`$. On the other hand, we will see in Section 4 that background light constraints require that any massive Population III stars must have formed much earlier than this.

### 3. Non-Baryonic Candidates

Not only is the presence of non-baryonic dark matter required by cosmological nucleosynthesis arguments, but it also expected in the Big Bang scenario. This is because any particles which exist in nature must have existed profusely in the early Universe once the temperature exceeded their rest mass $`m_X`$. Some of these particles should therefore have survived as “relics” once the temperature fell below the value $`T_F`$ at which their interactions froze-out. If the particles were relativistic at freeze-out ($`T_F>m_X`$), their current number density should be comparable to that of the microwave background photons, so one expects a relic density parameter $`\mathrm{\Omega }_X(m_X/100eV)`$. On the other hand, if they were non-relativistic at freeze-out ($`T_F<m_X`$), the current density declines with the particle mass ($`\mathrm{\Omega }_Xm_X^{1.9}`$) due to Boltzmann suppression. This gives two mass windows in which $`\mathrm{\Omega }_X`$ can be interesting: $`M_X10100`$ eV and $`M_X1030`$ GeV. Since the current velocity dispersion of the particles should be $`<v^2>^{1/2}T/m_X`$, which decreases with mass, particles in the smaller and larger ranges are termed “hot” and “cold” respectively. Particles can clump into bound objects of radius $`R`$ and velocity dispersion $`\sigma `$ only for

$$m_X>30(R/10kpc)^{1/2}(\sigma /100km/s)^{1/4}eV$$
(5)
and this gives a lower limit of 4 eV for clusters, 20 eV for galactic halos and 300 eV for dwarf galaxies.

Large-scale structure observations also require that the Universe be dominated by non-baryonic matter. This is because, in a purely baryonic universe, the density fluctuations required at decoupling to generate this structure would induce larger microwave background anisotropies than observed. However, this argument does not specify the nature of the non-baryonic dark matter precisely. The original “hot dark matter” (HDM) scenario, proposed in the early 1980s, assumed that it was in the form of neutrinos but that would be incompatible with the presence of dark matter in dwarf galaxies. This was superseded in the mid-1980s by the “cold dark matter” (CDM) scenario and this remained popular for almost a decade when it was found that it did not fit the latest large-scale structure data or give the observed gas fraction in clusters. In the mid-1990s models invoking a mixture of hot and cold dark matter (CHDM) became topical but these still had problems explaining the observed abundances of clusters, bright galaxies and Lyman-$`\alpha `$ clouds. Currently we have seen that a Universe containing both a cosmological constant and cold dark matter ($`\mathrm{\Lambda }CDM`$) is in vogue and some people have considered a hot component ($`\mathrm{\Lambda }CHDM`$) as well.

Eventually the combination of large-scale structure data and microwave background anisotropy measurements could determine the relative contributions of all the different types of dark matter very precisely (Gawiser & Silk 1998). However, as mentioned in the Introduction, it already seems clear that each of the components is needed and that their densities are all comparable. Of course, these arguments provide only indirect evidence for non-baryonic dark matter. More definitive evidence can come only from attempts to detect the hot or cold particles directly. There have been several recent interesting developments in this context, as we now discuss.

Massive Neutrinos. There are now three experiments which indicate a non-zero neutrino mass and hence a non-zero cosmological neutrino density \[see Caldwell (1999) for a recent review\]. (1) Studies of the ratio of muon and electron neutrinos produced by the decay of pions in the atmosphere suggest that there are oscillations between muon and tau neutrinos, corresponding to a mass difference $`\mathrm{\Delta }m_{\tau \mu }^22\times 10^3eV^2`$. (2) The deficit of electron-neutrinos coming from the Sun implies that these neutrinos must be oscillating into another neutrino species (labelled “X”) with a mass difference, if attributed to the MSW effect, given by $`\mathrm{\Delta }m_{eX}^210^5eV^2`$. (3) LSND claims to have detected oscillations between muon and electron neutrinos in neutrino beams coming from accelerator experiments, corresponding to a mass difference $`10eV^2>\mathrm{\Delta }m_{\mu e}^2>0.2eV^2`$. This is a controversial result but, if true, it implies $`m(\nu _\mu )>0.5eV`$ and hence $`\mathrm{\Omega }_\nu >0.005`$, which is already comparable to the visible density given by eqn (3). Since the first experiment implies $`\mathrm{\Delta }m_{\tau e}^2>>10^5eV^2`$, the particle X required to explain the solar neutrino deficit cannot be the tau nuetrino and so would most naturally be identified with the muon neutrino. In this case,

$$m(\nu _\tau )\mathrm{\Delta }m_{\tau \mu }0.05eV,m(\nu _\mu )\mathrm{\Delta }m_{\mu e}0.005eV,m(\nu _e)<<0.005eV.$$
(6)
This is consistent with the “see-saw” mechanism and suggests that the total neutrino density is very low. Another possibility, and the only one compatible with the LSND result, is that $`m(\nu _\tau )m(\nu _\mu )`$, so that one has mass degeneracy and $`m(\nu _\tau )>>\mathrm{\Delta }m_{\tau \mu }`$. In this case $`\mathrm{\Omega }_\nu `$ could be large and is limited only by the large-scale structure observations (i.e. $`\mathrm{\Omega }_\nu <0.2`$). However, the solar neutrino problem could not then be atributed to electron neutrinos oscillating into muon neutrinos since $`\mathrm{\Delta }m_{\mu e}`$ would be too large. One would therefore need to identify X with a new “sterile” neutrino, whose mass is very close to that of the electron neutrino, so that $`\mathrm{\Delta }m_{\mu e}>>\mathrm{\Delta }m_{Xe}`$. It has been claimed that the sterile neutrino may be useful for astrophysical purposes (eg. in explaining the abundance of r-process elements) but this remains controversial.

WIMPs. Supersymmetry theory requires the existence of a large number extra “sparticles” with a mass in the range $`0.11`$ TeV. \[See Roszkowski (1999) and Ellis (2000) and for recent reviews.\] The Lightest Supersymmetric Particle (LSP) should be stable and observations require that it has neither electric nor strong interactions. Three candidates have been suggested: the (spin 0) sneutrino, the (spin 3/2) gravitino, and the (spin 1/2) neutralino. However, the LEP results imply that the sneutrino mass must exceed $`1`$ TeV, which gives a negligible density, while the gravitino would have a mass of order $`1`$ keV, corresponding to “warm” dark matter. The latter would be inconsistent with large-scale structure observations, although another scenario has recently been proposed in which 100 GeV gravitinos are produced in the context of leptogenesis. The neutralino is therefore the favoured WIMP candidate and this would need to be some mixture of the photino, Higgsino and zino. LEP implies that the mass must be in the range $`30600`$ GeV, corresponding to $`\mathrm{\Omega }_Xh^20.1`$. One would expect there to be a continual flux of halo neutralinos passing through the Earth and efforts to detect these has been underway for more than a decade. Recently the DAMA group claims to have made a positive detection, associated with the annual modulation of the flux as the Earth goes around the Sun. They find a mass in the range $`4576`$ GeV but it should be stressed that this depends sensitively on the halo model.

Axions. The axion is a particle associated with the breaking of Peccei-Quinn symmetry at a temperature $`f_a10^{12}`$GeV and it aquires a mass $`m_a`$ due to instanton effects at a temperature $`T\mathrm{\Lambda }_{QCD}`$. In the original thermal production scenario, the associated parameters are

$$\mathrm{\Omega }_a\theta ^2(m_a/10^6eV)^{1.2},m_a\mathrm{\Lambda }_{QCD}^2/f_a10^5(f_a/10^{12})^1eV$$
(7)
where $`\theta `$ is the “misalignment” angle (Preskill et al. 1983). Several astrophysical constraints require $`10^6eV<m_a<10^2eV`$ and microwave cavity searches focus on the mass range $`10^610^5eV`$. Later it was realized that a network of axion strings could form at the Peccei-Quinn epoch (Davis & Shellard 1989) and axion domain walls could form at the QCD epoch (Shellard 1990). The decay and interaction of these strings and walls could produce an axion density in the range $`\mathrm{\Omega }_a(1001000)(m_a/10^6eV)^{1.2}`$, which is considerably larger than indicated by eqn (7). Various refinements of these scenarios arise in the inflationary scenarios, depending on whether inflation occurs at a higher or lower temperature than $`f_a`$ (Shellard & Battye 1998). The combination of the Supersymmetry and Peccei-Quinn symmetry also implies the existence of an axino. In the simplest scenario this has a mass of order a keV and therefore acts like warm dark matter. However, in a more recent scenario the photino decays into an axino and a photon and this predicts a much larger axino density (Covi et al. 1999).

Superheavy Relics. The existence of WIMPs is an inevitable prediction of the “bottom-up” approach of supersymmetry theory. However, one can also adopt various “top-down” approaches (eg. heterotic string theory and minimal supergravity) and these predict the existence of a plethora of particles with much larger mass. For example, Kaluza-Klein theory predicts many helicity states with mass near the Planck mass (Ellis et al. 1990) and the non-thermal decay of inflatons would produce a large density of GUT-mass particles (Kolb 2000). Such particles are sometimes described as “wimpzillas” and their decay might provide a possible explanation of the Ultra-High Energy (UHE) cosmic rays with $`E>10^{19}`$GeV. The origin of such cosmic rays is currently unknown but they must originate within 30-50 Mpc (because of the GZK effect) and they do not seem to be associated with individual galaxies or gamma-ray bursts. Another possible explanation would be the decay of topological defects. In any case, observations of UHE cosmic rays already place interesting constraints on such scenarios (Ziaeepour 1999).

Indirect Detections. There are various ways in which dark matter particles might be detected indirectly. For example, the annihilations of WIMPs in the Galactic halo would produce a flux of positron and antiproton cosmic rays (Bergstrom 1999). Balloon experiments have detected an excess of such particles but it is difficult to identify them specifically with WIMPs. They might also be secondary particles produced by the interaction of primary cosmic rays with interstellar protons and there is much uncertainty associated with the propagation of cosmic rays through the Galaxy and their modulation by the solar wind. WIMP-annihilations would also induce a a Galactic $`\gamma `$-ray background, with a characteristic enhancement of a $`\gamma `$-line near the Galactic centre in the CDM scenario. Such a background may have already been detected with the EGRET satellite (Chary & Wright 1998).

### 4. Baryonic Candidates

Over the last few decades, a number of observational claims have been made from time to time which seemed to indicate the existence of MACHOs. Many of these observations turned out to be spurious but they served the useful purpose of stimulating theorists to study cosmological scenarios involving baryonic dark matter and to try to predict their various signatures. Currently one cannot state definitely that MACHOs exist but one can place important constraints on the possible candidates. In this section we will discuss these candidates in turn, focussing particularly on brown dwarfs, red dwarfs, white dwarfs and black holes.

Background Radiation from VMOs. Stars in some range above $`8M_{}`$ would leave neutron star or black hole remnants but neither of these would be plausible candidates for either the disc or halo dark matter because their precursors would have unacceptable nucleosynthetic yields. Stars larger than $`200M_{}`$ are termed “Very Massive Objects” or VMOs and might collapse to black holes without these nucleosynthetic consequences (Carr et al. 1984). However, during their main-sequence phase, such VMOs would be expected to generate of a lot of background light. By today this should have been shifted into the infrared or submillimetre band, as a result of either redshift effects or dust reprocessing, so one would expect a sizeable infrared/submillimetre cosmic background (Bond et al. 1991). Over the last few decades there have been several reported detections of such a background but these have usually turned out to be false alarms. COBE does now seem to have detected a genuine infrared background (Fixsen et al. 1998) but this can probably be attributed to ordinary Population I and II stars. In any case, the current constraints on such a background strongly limit the density of any VMOs unless they form at a very high redshift. For this reason massive Population III stars would need to be pregalactic rather than protogalactic.

Dynamical Effects of SMOs. Stars larger than $`10^5M_{}`$ \- termed “Supermassive Objects” or SMOs - would collapse directly to black holes without any nucleosynthetic or background light production. However, supermassive black holes would still have noticeable dynamical effects and these have been investigated by many authors. The constraints on black holes in our own disc - due to the disruption of open clusters - and in our own halo - due to the heating of disc stars, the disruption of globular clusters and dynamical friction effects - have been discussed in detail by Carr & Sakellariadou (1998). Although it has been claimed that there is positive evidence for some of these effects, such as disc heating (Lacey & Ostriker 1985), the interpretation of this evidence is not clear-cut. It is therefore more natural to regard these dynamical effects as merely imposing an upper limit on the density of black holes or indeed any other type of compact object. It should be stressed that many of these limits would also apply if the black holes were replaced by “dark clusters” of smaller objects, a scenario which has been explored by many authors (Carr & Lacey 1987, Ashman 1990, Kerins & Carr 1994, De Paolis et al. 1995, Moore & Silk 1995).

Halo Light from Red Dwarfs. A particularly interesting recent development was the discovery of red light around NGC 5907 by Sackett et al. (1994), apparently emanating from low mass stars with a density profile like that of the halo. This detection was confirmed in V and I by Lequeux et al. (1996) and in J and K by James & Casali (1996). However, the suggestion that the stars might be of primordial origin (with low metallicity) was contradicted by the results of Rudy et al. (1997), who found that the color was indicative of low mass stars with solar metallicity. Furthermore, it must be stressed that the red light light has only been observed within a few kpc and no NIR emission is detected at 10-30 kpc (Yost et al. 1999). Both these points go against the suggestion that the red light is associated with MACHOs. Recently it has been suggested that the red light seen in NGC 5907 is more likely to derive from a disrupted dwarf galaxy, the stars of which would naturally follow the dark matter profile (Lequeux et al. 1998), or to be a ring left over from a disrupted dwarf spheroidal galaxy (Zheng et al. 1999). In any case, NGC 5907 does not seem to be typical since ISO observations of four other edge-on bulgeless spiral galaxies give no evidence for red halos (Gilmore & Unavane 1998).

Microlensing of Macrolensed Quasars by LMOs. If a galaxy is suitably positioned to image-double a quasar, then there is also a high probability that an individual halo object will traverse the line of sight of one of the images and this will give intensity fluctuations in one but not both images (Gott 1981). The effect would be observable for objects bigger than $`10^4M_{}`$ but the timescale of the fluctuations $`40(M/M_{})^{1/2}`$y would make them detectable only for $`M<0.1M_{}`$, so this method works only for “Low Mass Objects” or LMOs. There is already evidence of this effect for the quasar 2237+0305 (Irwin et al. 1989), the observed timescale for the variation in the luminosity of one of the images possibly indicating a mass below $`0.1M_{}`$ (Webster et al. 1991). However, because the optical depth is high, the mass estimate is very uncertain and a more recent analysis suggests that it could be in the range $`0.110M_{}`$ (Lewis et al. 1998), in which case the lens could be an ordinary star. The absence of this effect in the quasar 0957+561 has also been used to exclude MACHOs with mass in the range $`10^710^3M_{}`$ from making up all of the halo of the intervening galaxy, although the precise limit has some dependence on the quasar size (Schmidt & Wambsganss 1998).

Microlensing of Quasars by Jupiters. More dramatic but rather controversial evidence for the microlensing of quasars comes from Hawkins (1993, 1996, 1999), who has been monitoring 300 quasars in the redshift range $`13`$ for nearly 20 years using a wide-field Schmidt camera. He finds quasi-sinusoidal variations with an amplitude of 0.5 magnitudes on a timescale 5 y and attributes this to lenses with mass $`10^3M_{}`$. The crucial point is that the timescale decreases with increasing redshift, which is the opposite to what one would expect for intrinsic variations, although this has been disputed (Alexander 1995, Baganoff & Malkan 1995). The timescale also increases with the luminosity of the quasar and he explains this by noting that the variability timescale should scale with the size of the accretion disc (which should itself correlate with luminosity). A rather worrying feature of Hawkins’ claim is that he requires the density of the lenses to be close to critical (in order that the sources are transited continuously), so he has to invoke primordial black holes which form at the quark-hadron phase transition (Crawford & Schramm 1982). However, this requires fine-tuning since the fraction of the Universe going into black holes at this transition must only be about $`10^9`$. Walker (1999) has proposed that Hawkins’ lenses might also be jupiter-mass gas clouds.

Microlensing of Stars in the LMC. Attempts to detect microlensing by objects in our own halo by looking for intensity variations in stars in the Magellanic Clouds and the Galactic Bulge have now been underway for a decade (Paczynski 1996). This method is sensitive to lens masses in the range $`10^710^2M_{}`$ but the probability of an individual star being lensed is only $`\tau 10^6`$, so one has to look at many stars for a long time (Paczynski 1986). The duration and likely event rate are

$$P0.2(M/M_{})^{1/2}y,\mathrm{\Gamma }N\tau P^1(M/M_{})^{1/2}y^1$$
(8)
where $`N10^6`$ is the number of stars. The MACHO group currently has 13-17 LMC events and the durations span the range $`34230`$ days (Alcock et al. 2000). For a standard halo model, the data suggest an average lens mass of around $`0.5M_{}`$ and a halo fraction of 0.2, with the 95% confidence ranges being $`0.150.9M_{}`$ and $`0.080.5`$. The mass is comparable with the earlier estimates but the fraction is somewhat smaller (Alcock et al. 1997). This might appear to indicate that the MACHOs are white dwarfs but one could also invoke a less conventional candidate. For example, primordial black holes forming at the quark-hadron phase transition might have the required mass (Jedamzik 1997) and the microlensing implications of this scenario are discussed by Green (2000). Perhaps the most important result of the LMC searches is that they eliminate many candidates. Indeed the combination of the MACHO and EROS results already excludes objects in the mass range $`5\times 10^70.002M_{}`$ from having more than 0.2 of the halo density (Alfonso et al. 1997).

White Dwarfs. A few years ago white dwarfs (WDs) were regarded as rather implausible dark matter candidates. One required a very contrived IMF, lying between $`2M_{}`$ and $`8M_{}`$, in order to avoid excessive production of light or metals (Ryu et al. 1990); the fraction of WD precursors in binaries would be expected to produce too many type 1A supernovae (Smecker & Wyse 1991); and the halo fraction was constrained to be less than 10% in order to avoid the luminous precursors contradicting the upper limits from galaxy counts (Charlot & Silk 1995). The observed WD luminosity function also placed a severe lower limit on the age of any WDs in our own halo (Tamanaha et al. 1990). More recent constraints - from CNO production (Gibson & Mould 1997), helium and deuterium production (Field et al. 2000) and extragalactic background light limits (Madau & Pozzetti 1999) - have strengthened these limits (Graff et al. 1999). At least some of these limits must be reconsidered in view of recent claims by Hansen (1998) that metal-poor old WDs with hydrogen envelopes could be much bluer and brighter than previously supposed, essentially because the light emerges from deeper in the atmosphere. This suggestion has been supported by HST observations of Ibata et al. (1999), who claim to have detected five candidates of this kind. The objects are blue and isolated and show high proper motion. They infer that they are $`0.5M_{}`$ hydrogen-atmosphere WDs with an age of around $`12`$ Gyr. Three such objects have now been identified spectroscopically (Hodgkin et al. 2000, Ibata et al. 2000), so this possibility must be taken very seriously. However, this does not circumvent the nucleosynthetic arguments against WDs.

Brown Dwarfs. Objects in the range $`0.0010.08M_{}`$ would never burn hydrogen and are termed “brown dwarfs” (BDs). They represent a balance between gravity and degeneracy pressure. It has been argued that objects below the hydrogen-burning limit may form efficiently in pregalactic or protogalactic cooling flows (Ashman & Carr 1990, Thomas & Fabian 1990) but the direct evidence for such objects remains weak. While some BDs have been found as companions to ordinary stars, these can only have a tiny cosmological density and it is much harder to find isolated field BDs. The best argument therefore comes from extrapolating the initial mass function (IMF) of hydrogen-burning stars to lower masses than can be observed directly. The IMF for Population I stars ($`dN/dmm^\alpha `$ with $`\alpha <1.8`$) suggests that only 1% of the disc could be in BDs (Kroupa et al. 1995). However, one might wonder whether $`\alpha `$ could be larger, increasing the BD fraction, for zero-metallicity stars. Although there are theoretical reasons for entertaining this possibility, earlier observational claims that low metallicity objects have a steeper IMF than usual are now discredited. Indeed observations of Galactic and LMC globular clusters (Elson et al. 1999) and dwarf spheroidal field stars (Feltzing et al. 1999) suggest that the IMF is universal with $`\alpha <1.5`$ at low masses (Gilmore 1999). This implies that the BD fraction is much less than 1% by mass. However, it should be stressed that nobody has yet measured the IMF in the sites which are most likely to be associated with Population III stars. We have seen that the LMC microlensing results would now seem to exclude a large fraction of BDs in our own halo. Although Honma & Kan-ya (1998) have presented 100% BD models, these would require falling rotation curves and most theorists would regard these as rather implausible.

### 5. Conclusions

1) Cosmological nucleosynthesis calculations suggest that there is both baryonic and non-baryonic dark matter. From a theoretical perspective, the latter is better motivated since one cannot be sure that the dark baryons are inside galaxies. The more conservative conclusion would be that they are contained in an intergalactic medium or in gas within groups or clusters of galaxies.

2) Over the years there have been several observational claims which seem to indicate the existence of MACHOs but these have usually turned out to be false alarms. Currently the only positive evidence comes from microlensing observations and the LMC results suggest that white dwarfs may be the best MACHO candidate. However, it must be stressed that the mass estimate upon which this inference is based is sensitive to assumptions about the halo model and the large number of arguments which have been voiced against white dwarfs in the past cannot be brushed aside too cavalierly.

3) There are still many viable non-baryonic candidates. As with the baryonic candidates, there have been a number of flase alarms in the past. However, the most recent experimental results require that there be at least some hot dark matter in the form of neutrinos and there are also clams to have detected cold dark matter in the form of WIMPs.

4) We have seen that there are many important constraints on both baryonic and non-baryonic dark matter candidates from a wide variety of other astrophysical effects. Until there is a definite detection, therefore, the best strategy is to proceed by eliminating candidates, on the Sherlock Holmes principle that whatever candidate remains, however implausible, must be correct.

1. Alcock, C. et al. 1997, ApJ 486, 697

2. Alcock, C. et al. 2000, astro-ph/0001272

3. Alexander, T. 1995, MNRAS 274, 909

4. Alfonso, A. et al. 1997, ApJ 99, L12

5. Ashman, K.A. 1990, MNRAS 247, 662

6. Ashman, K.A., Carr, B.J. 1991, MNRAS 249, 13

7. Baganoff, F.K., Malkan, M.A. 1995, ApJ 444, L13

8. Bahcall, J.N., Flynn, C., Gould, A. 1992, ApJ 389, 234

9. Bergstrom, L. 1999, in Cosmo-98, ed. D.Caldwell, in press (astro-ph/9902172)

10. Bond, J.R., Carr, B.J., Hogan, C.J. 1991, ApJ 367, 420

11. Burles, S., Tytler, D. 1998, ApJ 499, 699; 507, 732

12. Burles, S. et al. 1999, Phys. Rev. Lett. 82, 4176

13. Caldwell, D. O. 1999, in The Identification of Dark Matter, ed. N.J.C. Spooner & V. Kudryavtsev, p 527 (World Scientific)

14. Carr, B.J. 1994, ARA&A 32, 531

15. Carr, B.J., Lacey, C.G. 1987, Ap.J 316, 23

16. Carr, B.J., Sakellariadou, M. 1998, ApJ 516, 195

17. Carr, B.J., Bond, J.R., Arnett, W.D. 1984, ApJ 277, 445

18. Carswell R.F. et al. 1994, MNRAS 268, L1

19. Cen, R., Ostriker, J.P. 1999, ApJ 514, 1

20. Charlot, S., Silk, J. 1995, ApJ 445, 124

21. Chary, R., Wright, E.L., 1998, in 3rd Stromlo Symposium, ed. B.K.Gibson et al. (ASP Conference Series) (astro-ph/9811324)

22. Copi, C.J., Schramm, D.N., Turner, M.S. 1995, Phys.Rev.Lett. 75, 3981

23. Covi, L., Kim, J.E., Roszkowski, L. 1999, Phys. Rev. Lett. 82, 4180.

24. Crézé, M. et al. 1998, A & A 329 920

25. Crawford, M., Schramm, D.N. 1982, Nature 298, 538

26. Davis, R.L., Shellard, E.P.S. 1989, Nuc. Phys. B. 324, 167

27. De Paolis, F. et al. 1995, A & A 295, 567; 299, 647

28. Ellis, J. 2000, Phys. Scripta T85, 221

29. Ellis, J., Lopez, J., Nanopoulos, D.V. 1990, Phys.Lett. B. 247, 257

30. Elson, R.A. et al. 1999, in Stellar Populations in the Magellanic Clouds, in press

31. Feltzing, S., Gilmore, G., Wyse, R.F.G. 1999, ApJ 516, L17

32. Fich, M., Tremaine, S. 1991, ARA&A 29, 409

33. Field, B.D., Freese, K., Graff, D.S. 2000, ApJ 534, 265

34. Fixsen, D.J. et al. 1998, ApJ 508, 123

35. Flynn, C., Fuchs, B. 1995, MNRAS 270, 471

36. Fukugita, M., Hogan, C.J., Peebles, P.J.E. 1998, ApJ 503, 518

37. Gawiser, E., Silk, J. 1998, Science 280, 1405

38. Gibson, B.K., Mould, J.R. 1997, ApJ 482, 98

39. Gilmore, G. 1999, in The Identification of Dark Matter, ed. N.J.C. Spooner & V. Kudryavtsev, p 121 (World Scientific)

40. Gilmore, G., Unavane, M. 1998, MNRAS 301, 813

41. Gott, J.R. 1981, ApJ 243, 140

42. Graff, D.S., Freese, K., Walker, T.P., Pinsonneault, M.H. 1999, ApJ 523, L77

43. Green, A. 2000, ApJ 537, 708

44. Hansen, B.M.S. 1998, Nature, 394 860

45. Hawkins, M.R.S. 1993, Nature 366, 242

46. Hawkins, M.R.S. 1996, MNRAS 278, 787

47. Hawkins, M.R.S. 1999, in The Identification of Dark Matter, ed. N.J.C. Spooner & V. Kudryavtsev, p 206 (World Scientific)

48. Hodgkin, S.T. et al. 2000, Nature 403, 57

49. Hogan, C.J. 1999, in Inner Space/Outer Space II; astro-ph/9912107

50. Honma, M., Kan-ya, Y. 1998, ApJ 503, L139

51. Ibata, R.A., Richer, H.B., Gilliland, R.L., Scott, D. 1999, ApJ 524, L95

52. Ibata, R.A. et al. 2000, astro-ph/0002138

53. Irwin, M.J. et al. 1989, AJ 98, 1989

54. James, P., Casali, M.M. 1996, Spectrum 9, 14

55. Jedamzik, K. 1997, Phys.Rev.D. 55, R5871

56. Kerins, E., Carr, B.J. 1994, MNRAS 266, 775

57. Kolb, E.W. 2000, Phys. Scripta T85, 231

58. Kuijken, P., Gilmore, G. 1989, MNRAS 239, 571; 605, 651

59. Kroupa, P., Tout, C., Gilmore, G. 1993, MNRAS 262, 545

60. Lacey, C.G., Ostriker, J.P. 1985, ApJ 299, 633

61. Lequeux, J. et al. 1996, A & A 312, L1

62. Lequeux, J. et al. 1998, A & A 334, L9

63. Lewis, G.F. et al. 1998, MNRAS 295, 573

64. Madau, P., Pozzetti, L. 2000, MNRAS 312, L9

65. Mathews, G.J., Schramm, D.N., Meyer, B.S. 1993, ApJ 404, 476

66. Paczynski, B. 1986, ApJ 304, 1; 308, L43

67. Paczynski, B. 1996, ARA&A 34, 419

68. Persic, M., Salucci, P. 1992, MNRAS 258, 14P

69. Preskill, J., Wise, M.B., Wilczek, F. 1983, Phys. Lett. B 120, 127

70. Roszkowski, L. 1999, in Cosmo-98, ed. D.Caldwell, in press (hep-ph/9903467)

71. Rudy, R.P. et al. 1997, Nature 387, 159

72. Rugers, H., Hogan, C.J. 1996, ApJ 459, L1

73. Ryu, D., Olive, K.A., Silk, J. 1990, ApJ 353, 81

74. Sackett, P.D. et al. 1994, Nature 370, 441

75. Schmidt , R., Wambsganss, J. 1998, A & A 335, 379

76. Shellard, E.P.S. 1990, in Formation and Evolution of Cosmic Strings, ed. Gibbons, G.W., Hawking, S.W. & Vachaspati, V. (Cambridge University Press)

77. Shellard, E.P.S., Battye, R.A. 1998, Phys. Rep. 307, 227.

78. Smecker, T.A., Wyse, R.F.G. 1991, ApJ 372, 448

79. Songaila, A. et al. 1994, Nature 368, 599

80. Tamanaha, C.M., Silk, J., Wood, M.A., Winget, D.E. 1990, ApJ 358, 164

81. Thomas, P., Fabian, A.C. 1990, MNRAS 246, 156

82. Turner, M.S. 1999, in Physics in Collision, ed. M.Campbell & T.M.Wells (World Scientific)

83. Tytler, D., Fan, X.M., Burles, S. 1996, Nature 381, 207

84. Walker, M.A. 1999, MNRAS 306, 504; 308, 551

85. Webb, J.K. et al. 1997, Nature 388, 250

86. Webster, R.L. et al. 1991, AJ 102, 1939

87. Weinberg, D.H. et al. 1997, ApJ 490, 564

88. Yost, S.A. et al. 1999, ApJ in press; astro-ph/9908364

89. Zaritsky, D. et al. 1993, ApJ 405, 464

90. Zheng, Z. et al. 1999, AJ 117, 2757

91. Ziaeepour, H. 1999, in The Identification of Dark Matter, ed. N.J.C. Spooner & V. Kudryavtsev, p 106 (World Scientific)

## Entry Form for the Proceedings

### 6. Title of the Paper

Baryonic and Non-Baryonic Dark Matter.

$``$

### 7. Author(s)

I am the author of this paper.

* + Full Name: Bernard John Carr
    + First Name: Bernard
    + Middle Name: John
    + Surname: Carr
    + Initialized Name: B. J. Carr
    + Affiliation: Queen Mary & Westfield College, University of London
    + E-Mail: B.J.Carr@qmw.ac.uk
    + Ship the Proceedings to: Astronomy Unit, Queen Mary & Westfield College, Mile End Road, London E1 4NS, England