text:
# 1 INTRODUCTION

## 1 INTRODUCTION

The protein folding problem is one of the most challenging problems in computational biophysics. The difficulty comes from the fact that the number of possible conformations for each protein is astronomically large. Simulations by conventional methods such as Monte Carlo (MC) or molecular dynamics (MD) algorithms in canonical ensemble will necessarily get trapped in one of many local-minimum states in the energy function. In order to overcome this multiple-minima problem, many methods have been proposed (for a review, see, e.g., Ref. ).

One way to alleviate the difficulty is to perform a simulation in a generalized ensemble where each state is weighted by a non-Boltzmann probability weight factor so that a random walk in potential energy space may be realized. The random walk allows the simulation to escape from any energy barrier and to sample much wider configurational space than by conventional methods. Monitoring the energy in a single simulation run, one can obtain not only the global-minimum-energy state but also canonical ensemble averages as functions of temperature by the single-histogram  and multiple-histogram  reweighting techniques.

One of the most well-known generalized-ensemble methods is perhaps multicanonical algorithm (MUCA)  (for a recent review, see Ref. ). MUCA was first introduced to the molecular simulation field in Ref. . Since then MUCA has been extensively used in many applications in protein and related systems (for a review, see, e.g., Ref. ).

While a simulation in multicanonical ensemble performs a free 1D random walk in potential energy space, that in simulated tempering (ST)  performs a free random walk in temperature space (for a review, see, e.g., Ref. ). This random walk, in turn, induces a random walk in potential energy space and allows the simulation to escape from states of energy local minima. ST has also been introduced to the protein folding problem .

The generalized-ensemble method is powerful, but in the above two methods the probability weight factors are not a priori known and have to be determined by iterations of short trial simulations. This process can be non-trivial and very tedius for complex systems with many local-minimum-energy states. Therefore, there have been attempts to accelerate the convergence of the iterative process for MUCA  (see also Ref. ).

In the replica-exchange method (REM) , the difficulty of weight factor determination is greatly alleviated. (REM is also referred to as multiple Markov chain method  and parallel tempering . For recent reviews with detailed references about the method, see, e.g., Refs. .) In this method, a number of non-interacting copies (or replicas) of the original system at different temperatures are simulated independently and simultaneously by the conventional MC or MD method. Every few steps, pairs of replicas are exchanged with a specified transition probability. REM has also been introduced to the protein folding problem . We further developed a multidimensional REM which is particularly useful in free energy calculations .

However, REM also has a computational difficulty: As the number of degrees of freedom of the system increases, the required number of replicas also greatly increases, whereas only a single replica is simulated in MUCA or ST. This demands a lot of computer power for complex systems. Our solution to this problem is: Use REM for the weight factor determinations of MUCA or ST, which is much simpler than previous iterative methods of weight determinations, and then perform a long MUCA or ST production run. The methods are referred to as the replica-exchange multicanonical algorithm (REMUCA)  and the replica-exchange simulated tempering (REST) . We have introduced a further extension of REMUCA, which we refer to as multicanonical replica-exchange method (MUCAREM) . In MUCAREM, the multicanonical weight factor is first determined as in REMUCA, and then a replica-exchange multicanonical production simulation is performed with a small number of replicas (for a review of all these new methods, see Ref. ).

In this article, we describe the six generalized-ensemble algorithms mentioned above. Namely, we first describe the three familiar methods: MUCA, ST, and REM. We then present the three new algorithms: REMUCA, REST, and MUCAREM. The effectiveness of these methods is tested with a 2-dimensional Potts model and protein systems.

## 2 METHODS

In the regular canonical ensemble with a given inverse temperature $`\beta 1/k_BT`$ ($`k_B`$ is the Boltzmann constant), the probability distribution of potential energy $`E`$ is given by

$$P_B(E;T)n(E)W_B(E;T)n(E)e^{\beta E},$$
(1)
where $`n(E)`$ is the density of states. Since the density of states $`n(E)`$ is a rapidly increasing function of $`E`$ and the Boltzmann factor $`W_B(E;T)`$ decreases exponentially with $`E`$, the probability distribution $`P_B(E;T)`$ has a bell-like shape in general. However, it is very difficult to obtain canonical distributions at low temperatures with conventional simulation methods. This is because the thermal fluctuations at low temperatures are small and the simulation will certainly get trapped in states of energy local minima.

Multicanonical algorithm (MUCA)  is one of the most well-known generalized-ensemble algorithms. In the “multicanonical ensemble” the probability distribution of potential energy is defined as follows:

$$P_{mu}(E)n(E)W_{mu}(E)\mathrm{constant}.$$
(2)
Because the multicanonical weight factor $`W_{mu}(E)`$ is (proportional to the inverse of the density of states and) not a priori known, one has to determine it for each system by iterations of trial simulations. See, for instance, Ref.  for details of the method to determine the MUCA weight factor $`W_{mu}(E)`$.

After the optimal MUCA weight factor is obtained, one performs a long MUCA simulation once. By monitoring the potential energy throughout the simulation, one can find the global-minimum-energy state. Moreover, by using the obtained histogram $`N_{\mathrm{mu}}(E)`$ of the potential energy distribution $`P_{\mathrm{mu}}(E)`$, the expectation value of a physical quantity $`A`$ at any temperature $`T=1/k_\mathrm{B}\beta `$ can be calculated from

$$<A>_T=\frac{{\displaystyle \underset{E}{}}A(E)n(E)e^{\beta E}}{{\displaystyle \underset{E}{}}n(E)e^{\beta E}},$$
(3)
where the best estimate of the density of states is given by the single-histogram reweighting techniques (see Eq. (2)) :

$$n(E)=\frac{N_{\mathrm{mu}}(E)}{W_{\mathrm{mu}}(E)}.$$
(4)
In simulated tempering (ST)  temperature itself becomes a dynamical variable, and both the configuration and the temperature are updated during the simulation with a weight:

$$W_{\mathrm{ST}}(E;T_m)=e^{\beta _mE+a_m},$$
(5)
where we discretize the temperature in $`M`$ different values, $`T_m`$ ($`m=1,\mathrm{},M`$). Without loss of generality we can order the temperature so that $`T_1<T_2<\mathrm{}<T_M`$. The lowest temperature $`T_1`$ should be sufficiently low so that the simulation can explore the global-minimum-energy region, and the highest temperature $`T_M`$ should be sufficiently high so that no trapping in a local-minimum-energy state occurs. The parameters $`a_m`$ are chosen so that the probability distribution of temperature is flat:

$$P_{\mathrm{ST}}(T_m)=𝑑En(E)W_{\mathrm{ST}}(E;T_m)=𝑑En(E)e^{\beta _mE+a_m}=\mathrm{constant}.$$
(6)
Hence, in simulated tempering the temperature is sampled uniformly. A free random walk in temperature space is realized, which in turn induces a random walk in potential energy space and allows the simulation to escape from states of energy local minima.

The parameters $`a_m`$ are not known a priori and have to be determined by iterations of short simulations. See, for instance, Ref.  for details of the method to determine the ST weight factor $`W_{\mathrm{ST}}(E;T_m)`$.

Note that from Eq. (6) we have

$$e^{a_m}𝑑En(E)e^{\beta _mE}.$$
(7)
The parameters $`a_m`$ are therefore “dimensionless” Helmholtz free energy at temperature $`T_m`$ (i.e., the inverse temperature $`\beta _m`$ multiplied by the Helmholtz free energy).

A simulation of ST is realized by alternately performing the following two steps . Step 1: A canonical MC or MD simulation at the fixed temperature $`T_m`$ is carried out for a certain MC or MD steps. Step 2: The temperature $`T_m`$ is updated to the neighboring values $`T_{m\pm 1}`$ with the configuration fixed. The transition probability of this temperature-updating process is given by the Metropolis criterion (see Eq. (5)):

$$w(T_mT_{m\pm 1})=\mathrm{min}(1,e^\mathrm{\Delta }),$$
(8)
where

$$\mathrm{\Delta }=\left(\beta _{m\pm 1}\beta _m\right)E\left(a_{m\pm 1}a_m\right).$$
(9)
After the optimal ST weight factor is determined, one performs a long ST simulation once. From the results of this production run, one can obtain the canonical ensemble average of a physical quantity $`A`$ as a function of temperature from Eq. (3), where the density of states is given by the multiple-histogram reweighting techniques  as follows. Let $`N_m(E)`$ and $`n_m`$ be respectively the potential-energy histogram and the total number of samples obtained at temperature $`T_m=1/k_\mathrm{B}\beta _m`$ ($`m=1,\mathrm{},M`$). The best estimate of the density of states is then given by

$$n(E)=\frac{{\displaystyle \underset{m=1}{\overset{M}{}}}g_m^1N_m(E)}{{\displaystyle \underset{m=1}{\overset{M}{}}}g_m^1n_me^{f_m\beta _mE}},$$
(10)
where

$$e^{f_m}=\underset{E}{}n(E)e^{\beta _mE}.$$
(11)
Here, $`g_m=1+2\tau _m`$, and $`\tau _m`$ is the integrated autocorrelation time at temperature $`T_m`$. Note that Eqs. (10) and (11) are solved self-consistently by iteration  to obtain the dimensionless Helmholtz free energy $`f_m`$ and the density of states $`n(E)`$.

The system for replica-exchange method (REM)  consists of $`M`$ non-interacting copies, or replicas, of the original system in canonical ensemble at $`M`$ different temperatures $`T_m`$ ($`m=1,\mathrm{},M`$). We arrange the replicas so that there is always one replica at each temperature. Then there is a one-to-one correspondence between replicas and temperatures. Let $`X=\{\mathrm{},x_m^{[i]},\mathrm{}\}`$ stand for a state in this generalized ensemble. Here, the superscript $`i`$ and the subscript $`m`$ in $`x_m^{[i]}`$ label the replica and the temperature, respectively. The state $`X`$ is specified by the $`M`$ sets of coordinates $`q^{[i]}`$. A simulation of REM is then realized by alternately performing the following two steps  (for details of the molecular dynamics version of REM, see Ref. ). Step 1: Each replica in the canonical ensemble at a fixed temperature is simulated simultaneously and independently for a certain number of MC or MD steps. Step 2: A pair of replicas, say $`i`$ and $`j`$, which are at neighboring temperatures, say $`T_m`$ and $`T_{m+1}`$, respectively, are exchanged: $`X=\{\mathrm{},x_m^{[i]},\mathrm{},x_{m+1}^{[j]},\mathrm{}\}X^{}=\{\mathrm{},x_m^{[j]},\mathrm{},x_{m+1}^{[i]},\mathrm{}\}`$. The transition probability of this replica exchange is given by the Metropolis criterion:

$$w(XX^{})=\mathrm{min}(1,e^\mathrm{\Delta }),$$
(12)
where

$$\mathrm{\Delta }\left(\beta _{m+1}\beta _m\right)\left(E\left(q^{[i]}\right)E\left(q^{[j]}\right)\right).$$
(13)
The replica-exchange multicanonical algorithm (REMUCA)  and replica-exchange simulated tempering (REST)  overcome both the difficulties of MUCA and ST (the weight factor determination is non-trivial) and REM (a lot of replicas, or computation time, is required).

In REMUCA  we first perform a short REM simulation (with $`M`$ replicas) to determine the MUCA weight factor and then perform with this weight factor a regular MUCA simulation with high statistics. The first step is accomplished by the multiple-histogram reweighting techniques . Let $`N_m(E)`$ and $`n_m`$ be respectively the potential-energy histogram and the total number of samples obtained at temperature $`T_m=1/k_\mathrm{B}\beta _m`$ of the REM run. The density of states $`n(E)`$ is then given by solving Eqs. (10) and (11) self-consistently by iteration . Once the estimate of the density of states is obtained, the multicanonical weight factor can be directly determined from Eq. (2).

In REST , just as in REMUCA, we first perform a short REM simulation (with $`M`$ replicas) to determine the ST weight factor and then perform with this weight factor a regular ST simulation with high statistics. The first step is accomplished by the multiple-histogram reweighting techniques , which give the dimensionless Helmholtz free energy $`f_m`$ (see Eqs. (10) and (11)). Once the estimate of the dimensionless Helmholtz free energy $`f_m`$ are obtained, the simulated tempering weight factor can be directly determined by using Eq. (5) where we set $`a_m=f_m`$ (compare Eq. (7) with Eq. (11)).

The formulations of REMUCA and REST are simple and straightforward, but the numerical improvement is great, because the weight factor determination for MUCA and ST becomes very difficult by the usual iterative processes for complex systems.

While multicanonical simulations are usually based on local updates, a replica-exchange process can be considered to be a global update, and global updates enhance the sampling further. Here, we present a further modification of REMUCA and refer to the new method as multicanonical replica-exchange method (MUCAREM) . In MUCAREM the final production run is not a regular multicanonical simulation but a replica-exchange simulation with a few replicas in the multicanonical ensemble. Because multicanonical simulations cover much wider energy ranges than regular canonical simulations, the number of required replicas for the production run of MUCAREM is much less than that for the regular REM, and we can keep the merits of REMUCA (and improve the sampling further).

## 3 RESULTS

We now present the results of our simulations based on the algorithms described in the previous section.

The first example is a spin system. We studied the 2-dimensional 10-state Potts model . The lattice size was $`34\times 34`$. This system exhibits a first-order phase transition . In Figure 1 we show the probability distributions of energy at three tempeartures (above the critical temperature $`T_C`$, at $`T_C`$, and below $`T_C`$). At the critical temperature we observe two peaks in the distribution, indicating that the system indeed undergoes a first-order phase transition.

In Figure 2 we show how the iterative procedure  for the MUCA weight factor determination converges. We see that a flat distribution in the entire energy range was obtained after 960,000 MC sweeps. Note that the convergence slows down drastically near the global-minimum-energy region (step from 300,000 MC sweeps to 960,000 MC sweeps).

In Figure 3 we show the results of our new method for the MUCA weight factor determination. We first made a REM simulation of 10,000 MC sweeps (for each replica) with 32 replicas (Figure 3(a)). Using the obtained energy distributions, we determined the (preliminary) MUCA weight factor by the REMUCA procedure as described in the previous section. Because the trials of replica exchange are not accepted near the critical temperature for first-order phase transitions, the probability distributions in Figure 3(a) for the energy range from $`1.5`$ to $`1.0`$ fails to have sufficient overlap, which is required for successful application of REM. This means that the MUCA weight factor, or density of states, in this energy range thus determined is of “poor quality.” With this MUCA weight factor, however, we made iterations of three MUCAREM simulations of 10,000 MC sweeps (for each replica) with 8 replicas (Figures 3(b), 3(c), 3(d)). In Figure 3(b) we see that the distributions are not completely flat, reflecting the poor quality in the phase-transition region. This problem is rapidly rectified as iterations continue, and the distributions are completely flat in Figure 3(d), which gives an optimal MUCA weight factor in the entire energy range. The details including the comparisons with the new method in Ref.  will be published elsewhere .

The second example is a protein system. We first illustrate how effectively generalized-ensemble simulations can sample the configurational space compared to the conventional simulations in the canonical ensemble. It is known by experiments that the system of a 17-residue peptide fragment from ribonuclease T1 tends to form $`\alpha `$-helical conformations. We have performed both a canonical MC simulation of this peptide at a low temperature ($`T=200`$ K) and a multicanonical MC simulation . In Figure 4 we show the time series of potential energy from these simulations.

We see that the canonical simulation thermalizes very slowly. On the other hand, the MUCA simulation indeed performs a random walk in potential energy space covering a very wide energy range. Four conformations chosen during this period (from 120,000 MC sweeps to 300,000 MC sweeps) are shown in Figure 5 for the MUCA simulation. The MUCA simulation indeed samples a wide conformational space.

The last example is a penta peptide, Met-enkephalin, whose amino-acid sequence is: Tyr-Gly-Gly-Phe-Met. In Figure 6, we show the average potential energy of Met-enkephalin in gas phase as a function of temperature that was calculated by the single- and multiple-histogram reweighting techniques from the four generalized-ensemble algorithms, MUCA, REMUCA, MUCAREM, and REST . The results are in good agreement.

## 4 CONCLUSIONS

In this article we have described the formulations of the three well-known generalized-ensemble algorithms, namely, multicanonical algorithm (MUCA), simulated tempering (ST), and replica-exchange method (REM). We then introduced three new generalized-ensemble algorithms that combine the merits of the above three methods, which we refer to as replica-exchange multicanonical algorithm (REMUCA), replica-exchange simulated tempering (REST), and multicanonical replica-exchange method (MUCAREM).

With these new methods available, we believe that we now have working simulation algorithms for spin systems and protein systems.

## 5 Acknowledgements

Our simulations were performed on the HITACHI and other computers at the Research Center for Computational Science, Okazaki National Research Institutes. This work is supported, in part, by a grant from the Research for the Future Program of the Japan Society for the Promotion of Science (JSPS-RFTF98P01101).


format: text

text:
# Electric Field Screening by Pairs in the Presence of Returning Positrons

## 1. Introduction

High energy pulsed emission from rotation powered pulsars are accounted for by presuming an electric potential drop along a magnetic flux in the magnetosphere, as typically seen in ‘polar cap models’ or ‘outer gap models’. The available potential drop produced by the neutron star is enough for the emission. In particular, for young pulsars acceleration potential is only a small portion of the total electromotive force. Rather important is how we understand localization of the electric field. In some mechanism, the actual charge density deviates from the Goldreich-Julian (GJ) charge density to produce field-aligned electric field. It is sometimes naively assumed that if pairs are created with a density higher than the GJ density, then produced field-aligned electric field is screened out and shut down the accelerator localized in a small region. In this work, we treat such a screening process by pairs and find out conditions to obtain a finite potential drop.

Let us consider a steady model for the field-aligned accelerator which has a finite potential drop with the electric field screened at the both ends. For definite sign of charge, let us assume electrons are accelerated outwards, i.e., the electric field points toward the star. Gamma-rays emitted by the electrons convert into pairs beyond a certain surface called the pair production front (PPF); beyond PPF pairs are assumed to be created continuously in space. A schematic structure is shown in Fig. 1, and the electric potential may be given by, in non-dimensional form,

$$\frac{d^2\varphi }{dx^2}=\frac{j}{v_1/c}+\frac{B_z}{B}+\overline{\rho }_{}+\overline{\rho }_+$$
(1)
where the first negative term is due to the primary electron beam, the second term is the normalized GJ density, and the last two terms appear due to pairs only beyond PPF.

In the ‘standard’ polar cap model, the negative space charge (cathode) is produced by non-relativistic electrons near the surface (when $`v_1c`$). If the length of the acceleration region along field lines is larger than the transverse scale, negative charge appearing on the side wall also contributes. Such an electric field can be screened by two possible reasons; one is when the GJ term (the second term) dominates the primary electronic charge (the first term), and the other is when pairs are created.

In the former case, it is obvious that the primary current should be less than a value, typically of order of GJ value, depending how slowly $`B_z`$ decreases as compared with $`B`$. As is well-known, field lines curving toward the rotation axis, $`B_z/B`$ term increases and is favorable to the case (Scharlemann, Arons, & Fawley; 1978, Arons, & Scharlemann; 1979). In order to shutdown the potential drop, the current density should be adjusted to be a critical value determined by the field geometry. However, the current density is determined globally, especially, through interaction with the wind zone, so the actual current density is likely to be different from the critical current density. The screening is nothing to do with the pair creation in this case. If pair creation takes place, the strict condition on the current density is slightly relaxed; as will be shown below some of pair positron returns back to the star and left pair electrons behind, which provide negative space charge in the ‘anode region’. Therefore if the actual primary current is less than the critical value, i.e., if the region needs negative charge to adjust the GJ (positive) term, then pair plasma helps to adjust the space charge to close the accelerator. If the situation is opposite, namely, if current density is much higher than the critical value, there is no way to close the accelerating electric field, even in the away curved field lines.

I would like to argue that such super critical current density seems realistic. The critical charge density is more or less GJ. The GJ current is obtained by some sort of dimensional analysis: the current derived form the magnetic moment $`\mu `$, angular velocity $`\mathrm{\Omega }`$ and other physical constants is $`I=\mu ^2\mathrm{\Omega }^4/c`$. This estimate for the circulating current in the magnetosphere will be correct. However, the GJ current density $`B\mathrm{\Omega }/2\pi c`$ is just $`I`$ divided by the polar cap area. It is not likely that current density is distributed more or less uniformly over the polar caps, and much likely that it is localized in annular rings or clumps, hinted by radio profiles and also by terrestrial aurora. An expected current density can be by factor of 10 or 100 larger than the GJ current density. If one assume such a super GJ current density, it is certain that strong accelerating electric field appears as shown by Shibata (1997) regardless of the field geometry.

The second case that the screening is due to pairs may be sometimes assumed. Pair polarization is an effect that screens out the electric field. The pair positrons are decelerated to non-relativistic speed soon after their birth, while the pair electrons are accelerated, and as a result of continuity, a positive space charge appears to reduce the electric field. However, some positrons reflected backward to the star by un-screened electric field. In this paper, we calculate the positive space charge produced by pair polarization and obtain the electric field strength which can be screened.

## 2. Model

We assume one dimensional primary stream along a field line, and pairs, which is created beyond the PPF more or less uniformly. Pair density is parameterized by ’multiplication factor’ $`M(x)`$, by which the primary electronic flux is multiplied to give pair flux at the point $`x`$. The Poisson equation for the electric potential and equations of motion and of continuity for electrons and positrons are solved self-consistently for a steady state solution.

## 3. Results

We obtain the electric field strength which can be screened, given a pair creation rate. It is found that returning of positrons makes the screening difficult seriously, because the pair electrons left behind the returning positrons produce negative space charge in the screening region where the positive space charge is required. As a result, the thickness of the screening is restricted to be as small as the braking distance $`\mathrm{\Delta }smc^2/e|E_{}|`$ for which positrons become non-relativistic, where $`E_{}`$ is the electric field strength just before the PPF. We confirmed the previous result of Shibata et al. (1998) that the electric field which can be screened by pair polarization is fairly small:

$$\frac{E_{}^2}{8\pi }<mc^2\left(\frac{\mathrm{\Omega }B}{2\pi ce}\right)\zeta ^{}j\mathrm{\Delta }M_{\mathrm{screen}}$$
(2)
where $`\mathrm{\Delta }M_{\mathrm{screen}}`$ is the pair multiplication factor within $`\mathrm{\Delta }s`$, $`\zeta ^{}`$ is a numerical factor of order of unity.

If the primary current density is of order of the Goldreich-Julian (GJ) value, the required pair multiplication factor per one primary electron is enormously large and cannot be realized in the conventional pair creation models. A previous belief that pair creation with a pair density higher than the GJ density immediately screens out the electric field is unjustified.

Some mechanism to salvage this difficulty should be found. As has been mentioned in §1, the toward curvature provides an effective positive charge, so that the required multiplication factor is considerably reduced. In this case, over-screening by the positive GJ term is canceled by the pair electrons left behind the returning positrons (this is the opposite of pair polarization). Harding & Muslimov (2001) is the case, but the current is restricted to be sub-GJ. In a more likely case with super-GJ current density, toward curvature will not be helpful. In any case, physics in focused field-aligned current, e.g., effects of inhomogeneity across the magnetic field lines, two stream instability in a high current density with unscreened electric field, has yet been properly studied.

Another possible way out is to include frictional forces between various components of charged currents. Since friction on positrons pulls them outwards along with electrons, returning fraction of positrons will be reduced, which will make screening easier. Although there are some ideas for physical processes of friction (two stream instability, production of positroniums and others), it is not clear whether the frictional force can be strong enough to lift positrons and screen the electric field. Study on the frictional interaction under unscreened electric field is strongly demanded.

More detailed analysis and results will be published in near future (Shibata, Miyazaki & Takahara 2001).

### Acknowledgments

We thank D. Melrose and Q. Luo for useful comments and discussions. This work is supported in part by Grant-in-Aid for Scientific Research from the Ministry of Education, Science, and Culture of Japan Nos. 12640229.

## References

Arons, J., Scharlemann, E.T., 1979, ApJ, 231, 854

Harding, A.K., Muslimov, A.G., 2001, ApJ, 556, 987

Miyazaki, J.& Takahara, F., 1997, MNRAS, 290, 49

Scharlemann, E.T., Arons, J., Fawley, W.M., 1978, ApJ, 222, 297

Shibata, S., 1997, MNRAS, 287, 262

Shibata, S., Miyazaki, J. & Takahara, F., 1998, MNRAS, 295, L53

Shibata, S., Miyazaki, J. & Takahara, F., 2001, submitted to MNRAS


text:
# Glueball Properties at Finite Temperature

## 1 INTRODUCTION

Finite temperature QCD, including the quark gluon plasma(QGP) physics, is one of the most interesting subjects. Even in the hadronic phase, a lot of effective models suggest the changes in the hadronic properties. However, lattice QCD studies of thermal hadron properties have been inadequate until quite recent. This is due to the difficulty in measuring hadronic two-point correlators at finite temperature. At high temperature, due to the shrink in the temporal extention and the consequent decrease in the number of data, the mass measurement from the temporal correlations is difficult. Hence, finite temperature mass shifts have been studied through the spatial correlations, which, however, is afflicted with the mixture of large Matsubara frequencies . Recently, the use of anisotropic lattice is established , which has the finer temporal lattice spacing $`a_t`$ than the spatial one $`a_s`$. In this way, available data increases, and the accurate mass measurement directly from the temporal correlations becomes possible . Here, we report the thermal properties of the $`0^{++}`$ glueball based on the quenched SU(3) anisotropic lattice QCD.

## 2 SMEARING METHOD

We consider the glueball correlator as $`G(t)=O(t)O(0)`$, $`O(t)={\displaystyle \underset{\stackrel{}{x}}{}}\left(\text{}\stackrel{~}{O}(\stackrel{}{x},t)\stackrel{~}{O}\right)`$ where $`\stackrel{~}{O}(\stackrel{}{x},t)=\text{Re}\text{Tr}(P_{12}(\stackrel{}{x},t)+P_{23}(\stackrel{}{x},t)+P_{31}(\stackrel{}{x},t))`$ denotes an interpolating field for the $`0^{++}`$ glueball. $`P_{ij}(\stackrel{}{x},t)\mathrm{SU}(3)`$ denotes the plaquette operator in the $`i`$-$`j`$-plain. We apply the spectral representation to $`G(t)`$ to have $`G(t)/G(0)=C_ne^{E_nt},`$ where $`C_n=|n|O|0|^2/|k|O|0|^2`$ with $`E_n`$ being the energy of the $`n`$-th excited state $`|n`$. Here, $`|0`$ denotes the vacuum and $`|1`$ denotes the ground state glueball. $`C_n`$ is non-negative with $`C_n=1`$. It is known that $`G(t)`$ receives large contributions from the excited states and the ground state contribution is small. As a consequence, the extracted mass always behaves much heaver than the ground-state mass. The problem originates from the fact that the “size” of the ordinary plaquette operator is $`O(a_s)`$, which is much smaller than the physical size $`R`$ of the glueball. Hence, the problem can be resolved by generating physically extended glueball operators, which is achieved by the smearing method . This extended operator is refered to as the smeared operator, which is obtained by replacing the ordinary link variables $`U_i(s)`$ in the plaquette operator $`P_{ij}(s)`$ by the fat link variables $`\overline{U}_i(s)\mathrm{SU}(3)`$. $`\overline{U}_i(s)`$ is defined so as to maximize $`\mathrm{ReTr}\left(\overline{U}_i^{}(s)V_i(s)\right)`$, where $`V_i(s)\alpha U_i(s)+{\displaystyle \underset{ji,\pm }{}}U_{\pm j}(s)U_i(s\pm \widehat{j})U_{\pm j}^{}(s+\widehat{i})`$ with $`U_\mu (s)U_\mu ^{}(s\widehat{\mu })`$. $`\alpha `$ is a real parameter. The summation involves only the spatial directions to avoid the temporal nonlocality. Note that $`\overline{U}_i(s)`$ holds the same gauge transformation properties as $`U_i(s)`$. We refer to the fat link defined in this way as the first fat link $`U_i^{(1)}(s)`$. The $`n`$-th fat link $`U_i^{(n)}(s)`$ is defined iteratively by $`U_i^{(n)}(s)\overline{U}_i^{(n1)}(s)`$ starting from $`U_i^{(1)}(s)\overline{U}_i(s)`$. The plaquette operator constructed with $`U_i^{(n)}(s)`$ is refered to as the $`n`$-th smeared operator. We next consider the size of the $`n`$-th smeared operator. By using the linearization and the continuum approximation, we obtain the diffution equation as

$$\frac{}{n}\varphi _i(n;\stackrel{}{x})=D\mathrm{}\varphi _i(n;\stackrel{}{x}),D\frac{a_s^2}{\alpha +4}.$$
(1)
$`\varphi _i(n;\stackrel{}{x})`$ describes the distribution of the gluon field $`A_i(\stackrel{}{x},t)`$ in the $`n`$-th smeared plaquette. In the $`n`$-th smeared plaquette located at the origin $`\stackrel{}{x}=\stackrel{}{0}`$, the gauge field $`A_i(\stackrel{}{x},t)`$ is distributed in the Gaussian form as $`e^{\stackrel{}{x}^2/(4Dn)}/(4\pi Dn)^{3/2}`$. Hence the size of the operator is estimated as

$$R\sqrt{\stackrel{}{x}^2}=\sqrt{6Dn}=a_s\sqrt{\frac{6n}{\alpha +4}}.$$
(2)
The original aim of the smearing method is the accurate mass measurement by maximizing the ground-state overlap. However, it can be also used to give a rough estimate of the physical glueball size. In fact, once the maximum overlap is achieved with some $`n`$ and $`\alpha `$, the glueball size can be estimated with Eq.(2).

## 3 SU(3) LATTICE QCD RESULT

We use the SU(3) anisotropic lattice plaquette action as $`S_G={\displaystyle \frac{\beta }{N_c}}{\displaystyle \frac{1}{\xi _0}}{\displaystyle \underset{s,i<j3}{}}\mathrm{ReTr}(1P_{ij}(s))+{\displaystyle \frac{\beta }{N_c}}\xi _0{\displaystyle \underset{s,i3}{}}\mathrm{ReTr}(1P_{i4}(s))`$ with the plaquette operator $`P_{\mu \nu }(s)`$ in the $`\mu `$-$`\nu `$-plane. The lattice parameter and the bare anisotropy parameter are fixed as $`\beta =6.25`$, $`\xi _0=3.2552`$, respectively. These parameters reproduce the spatial lattice spacing as $`a_s^1=2.272(16)`$ GeV, and the temporal one as $`a_t^1=9.088(64)`$ GeV. The renormalized anisotropy is $`\gamma a_s/a_t=4`$ . Here, the scale unit is determined by reproducing the string tension as $`\sqrt{\sigma }=427`$ MeV from the on-axis data of the static inter-quark potential. The pseudo-heat-bath algorithm is used to update the gauge field configurations on the lattice of the sizes $`20^3\times N_t`$ with various $`N_t`$ listed in Table 1. For each temperature, we pick up gauge field configurations every 100 sweeps for measurements, after skipping more than 20,000 sweeps of the thermalization. The numbers of gauge configurations used in our calculations are summarized in Table 1. From the analysis of the Polyakov loop, we find $`T_c260`$ MeV. To enhance the ground-state contribution, we adopt the smearing method with $`\alpha =2.1`$.

In Fig.1, the $`0^{++}`$ glueball correlator $`G(t)/G(0)`$ for $`N_{\mathrm{smear}}=40`$ at $`T=245`$ MeV is shown. The statistical errors are estimated with the jackknife analysis. The solid line denotes the best single hyperbolic cosine fitting as $`G(t)/G(0)=C(e^{m_\mathrm{G}ta_t}+e^{m_\mathrm{G}(N_tt)a_t})`$, performed in the interval indicated by vertical dotted lines. In the most suitable smearing number $`N_{\mathrm{smear}}^{\mathrm{best}}`$, the ground-state overlap $`C`$ is maximized and the mass $`m_G`$ is minimized, which indicates the achievement of the ground-state enhancement. In practical calculations, the maximum overlap and the mass minimization are achieved at almost the same $`N_{\mathrm{smear}}`$, and both of these two conditions would work as an indication of the maximal ground-state enhancement. Here, we take the maximum ground-state overlap condition.

In Table 1, we summarize the SU(3) lattice QCD result of the lowest $`0^{++}`$ glueball mass $`m_\mathrm{G}(T)`$ at various temperatures. We observe that about 20 % mass reduction of the lowest $`0^{++}`$ glueball near $`T_c`$ as $`m_\mathrm{G}(T)=1.25\pm 0.1`$ GeV for $`0.8T_cTT_c`$ in comparison with $`m_\mathrm{G}(T0)1.51.7`$ GeV . To estimate the glueball size, we search $`N_{\mathrm{smear}}^{\mathrm{best}}`$ which achieves the maximum overlap as $`C=C^{\mathrm{max}}`$. Considering about 0.5 % errors of $`C`$, we define the “smearing window” $`N_{\mathrm{smear}}^{\mathrm{window}}`$ where $`C0.995C^{\mathrm{max}}`$ is satisfied. From Eq.(2) with $`N_{\mathrm{smear}}^{\mathrm{window}}`$, we estimate the glueball size as $`R0.50.6`$ fm both at low and high temperatures. In Table 1, we summarize the lowest $`0^{++}`$ glueball mass $`m_\mathrm{G}(T)`$, the ground-state overlap $`C^{\mathrm{max}}`$, uncorrelated $`\chi ^2/N_{\mathrm{DF}}`$, $`N_{\mathrm{smear}}^{\mathrm{best}}`$, $`N_{\mathrm{smear}}^{\mathrm{window}}`$, the number of gauge configurations $`N_{\mathrm{conf}}`$ and the glueball size $`R`$.

## 4 SUMMARY

We have studied the $`0^{++}`$ glueball properties at finite temperature using SU(3) anisotropic quenched lattice QCD with more than 1,000 gauge configurations at each temperature. From the temporal correlation analysis with the smearing method, we have observed about 20 % mass reduction of the lowest $`0^{++}`$ glueball as $`m_G(T)=1.25\pm 0.1`$ GeV for $`0.8T_c<T<T_c`$, while no significant change has been observed for the glueball size as $`R0.50.6`$ fm.


text:
# Bose-Einstein Correlations in WW Events at LEP2

## 1 Introduction

Bose-Einstein interference is observed in hadronic Z-decay as an enhanced production of identical bosons at small four-momentum difference . There is no reason why such an interference should not be present within hadronic W-decay (intra-W BE interference), as well. Indeed, the four LEP collaborations measure intra-W Bose-Einstein correlations (BEC), as a consequence of intra-W BE interference, see  and references therein. Furthermore, since in fully hadronic WW events ($`\mathrm{W}^+\mathrm{W}^{}\mathrm{q}_1\overline{\mathrm{q}}_2\mathrm{q}_3\overline{\mathrm{q}}_4`$) the W-decay products overlap in space-time at LEP2 energies, it may be natural to expect interference also between identical bosons originating from different W’s (inter-W BE interference). Together with colour reconnection , inter-W BE interference not only forms a potential bias in the determination of the W-mass, but also may provide a laboratory to measure the space-time development of this overlap. In this contribution, the latest preliminary results of the four LEP experiments on inter-W BEC are presented.

## 2 Analysis

The largely model independent method used to study inter-W BEC is based on that proposed in.  If the two W’s would decay independently, the two-particle density in fully hadronic WW events, $`\rho _2^{\mathrm{WW}}`$, would be given by

$`\rho _2^{\mathrm{WW}}(p_1,p_2)`$ $`=`$ $`\rho _2^{\mathrm{W}^+}(p_1,p_2)+\rho _2^\mathrm{W}^{}(p_1,p_2)+\rho _1^{\mathrm{W}^+}(p_1)\rho _1^\mathrm{W}^{}(p_2)+\rho _1^\mathrm{W}^{}(p_1)\rho _1^{\mathrm{W}^+}(p_2)`$ (1)
$`=`$ $`2\rho _2^\mathrm{W}(p_1,p_2)+2\rho _1^\mathrm{W}(p_1)\rho _1^\mathrm{W}(p_2),`$
where the superscripts, e.g., WW, indicate which W’s decay hadronically. Experimentally, $`\rho _2^{\mathrm{WW}}`$ is measured in the fully hadronic WW events and $`\rho _2^\mathrm{W}`$ in the semi-hadronic events, i.e., events where one W decays hadronically and the other W decays into a charged lepton and its corresponding neutrino. To measure the product of the single-particle densities, a two-particle density, $`\rho (p_1,p_2)_{\mathrm{mix}}^{\mathrm{W}^+\mathrm{W}^{}}`$, is constructed by pairing particles originating from two different semi-hadronic WW events. By construction these pairs of particles are uncorrelated.

Since BEC are large at small values of the four-momentum difference, $`Q`$, between two identical particles, and since the statistics are too poor (about 10k WW events are selected per experiment at LEP2, at $`160<\sqrt{s}<209\mathrm{GeV}`$) to do a multi-dimensional analysis, the two-particle densities are measured in this one-dimensional distance measure: $`\rho _2(Q)\frac{1}{N_{\mathrm{ev}}}\frac{\mathrm{d}n_{\mathrm{pairs}}}{\mathrm{d}Q}`$, where $`N_{\mathrm{ev}}`$ is the number of events that is used and $`n_{\mathrm{pairs}}`$ is the number of particle pairs.

Deviations from the hypothesis that the two W’s decay independently can be found in deviations from Eq. (1). Therefore, the following test statistics are used:

$`\mathrm{\Delta }\rho \rho _2^{\mathrm{WW}}2\rho _2^\mathrm{W}2\rho _{\mathrm{mix}}^{\mathrm{W}^+\mathrm{W}^{}}\mathrm{and}D{\displaystyle \frac{\rho _2^{\mathrm{WW}}}{2\rho _2^\mathrm{W}+2\rho _{\mathrm{mix}}^{\mathrm{W}^+\mathrm{W}^{}}}}.`$ (2)
To diminish artificial distortions due to event mixing, non-BEC and various detector effects, also the double ratio $`D^{}D/D_{\mathrm{MC},\mathrm{noBE}}`$ is used, where $`D_{\mathrm{MC},\mathrm{noBE}}`$ is a MC sample without BEC, or at least without inter-W BEC. Evidently, $`\mathrm{\Delta }\rho 0`$ and $`D=D^{}1`$ in the case of inter-W BEC.

## 3 Results

Fig. 2 shows the distributions of $`D`$ and $`D^{}`$ for like-sign track pairs, for the L3 data collected in the period 1998-2000. It is clear that the data show flat distributions in both cases (the same conclusion is obtained for the $`\mathrm{\Delta }\rho (Q)`$ distribution, which is not shown), indicating no evidence for inter-W BEC and $`DD^{}`$ (i.e., only a small model dependence). Also shown in the figure are the predictions of a MC with the BE32 algorithm  and tuned to the Z-data, for the scenario where correlations between all particles are allowed (BEA) and for the scenario where only intra-W BEC are allowed (BES). It is clear that the BEA scenario is disfavoured, while BES describes the data. In particular, $`D^{}(Q)`$ is fitted by $`(1+\epsilon Q)(1+\mathrm{\Lambda }\mathrm{exp}(k^2Q^2))`$, where $`\epsilon `$, $`\mathrm{\Lambda }`$ and $`k`$ are the fit parameters. The parameter $`\mathrm{\Lambda }`$ measures the strength of inter-W BEC. The result of the fit is $`\mathrm{\Lambda }=0.008\pm 0.018(\mathrm{stat})\pm 0.016(\mathrm{syst})`$. This can be interpreted as evidence against BEC between identical particles originating from different W’s.

Analogously, Fig. 2 shows $`D(Q)`$ and $`\mathrm{\Delta }\rho (Q)`$ for the DELPHI data. The same conclusions are obtained as for the L3 data: there is no evidence for inter-W BEC and a MC with inter-W BEC is clearly disfavoured. In particular, $`D(Q)`$ is fitted by $`\gamma (1+\epsilon Q)(1+\mathrm{\Lambda }\mathrm{exp}(kQ))`$. Fixing $`k`$ at 1.01 fm (computed from a MC with inter-W BEC, using BE32), one obtained $`\mathrm{\Lambda }=0.037\pm 0.055(\mathrm{stat})\pm 0.055(\mathrm{syst})`$, consistent with no inter-W BEC.

The ALEPH collaboration also computed the $`\mathrm{\Delta }\rho (Q)`$ and $`D(Q)`$ test statistics. Qualitatively, the same observations are made as by L3 and DELPHI: the data show flat distributions, indicating no evidence for inter-W BEC. A MC with BEA is disfavoured and the BES scenario describes the data. No numbers are computed, yet. The ‘standard’ analysis of ALEPH differs from the method explained in Sect. 2. It starts by defining a correlation function $`R_2\rho _2/\rho _0`$, where $`\rho _0`$ is the two-particle density that would occur in a world without BEC. This correlation function is computed for the selected fully hadronic WW events and for two MC scenarios, BEA and BES, using the BE3 algorithm  (the MC is tuned on a light-quark Z-decay sample). One observes that the BES scenario describes the data, whereas BEA is disfavoured. Although this result agrees with the results from the method explained in Sect. 2, it should be realized that there is a strong model dependence, due to the MC that is needed to compute $`\rho _0`$, and that $`R_2`$ is also sensitive to intra-W BEC. These are two disadvantages with respect to the method of Sect. 2.

Unfortunately, the OPAL collaboration has no new results since 1999. A new analysis, using the method of Sect. 2 is in progress, but is not finalized, yet. The results of the ‘old’ analysis are based on a completely different method and due to the lack of statistics the results are inconclusive, but not in contradiction with the others.

## 4 Summary and Interpretation of the Results

A common method is used by L3, DELPHI and ALEPH (and will be used by OPAL), to search for inter-W BEC. No evidence is found for the existence of these correlations and MC models allowing inter-W BEC are (greatly) disfavoured. The picture is now more consistent than a year ago:  There is a generally accepted common method and the results from all experiments are more or less consistent. This provides a basis to combine results of the four LEP experiments into one single number in the near future.

Provided that there is no colour reconnection, the Lund string model  can only give the BE effect between bosons coming from the same string. Thus, the data, that do not show any indication for inter-W BEC, are no problem for this model.

The positive outcome of the BE analyses, is that the effect on the W mass measurement due to inter-W BEC, is small. At the time that W production started at LEP2, the worst case scenario predicted an effect on the measurement of the W mass up to 100 $`\mathrm{MeV}`$ Before the results shown in this contribution, the L3 collaboration estimated the error to be 60 $`\mathrm{MeV}`$ Now, it has been reduced to 20 $`\mathrm{MeV}`$ The error is expected to go down even further, when the LEP results are combined.


text:
# Mode-Coupling Theory of Colloids with Short-range Attractions

## 1 Introduction

The experimental study of the glass transition in colloidal systems  has been a very important test case to assess the validity of theories concerning the formation of an amorphous solid from super-cooled liquids. In particular the mode-coupling theory (MCT) has been useful when applied to colloids, modeled as spherical particles interacting through a hard sphere purely repulsive potential . In this case the MCT predicts the existence of a critical volume fraction $`\varphi `$ where the system undergoes an ergodic-nonergodic transition, which was observed experimentally using quasi-elastic dynamic light scattering . From the physical point of view the MCT describes in a fairly accurate way the so-called cage effect, i.e. the fact that at high densities molecular motions of particles are constrained by the presence of the surrounding ones which form a cage around each particle. Prior to reaching the MCT transition, we may think of the particle vibrating within their cages at short time-scales, and escaping from their cages at somewhat longer time-scales. The two time scales show up in a well-defined way, when the liquid system gets closer to the critical threshold, as distinct relaxation regimes separated by a plateau region which becomes more and more extended as criticality is approached. Close and above the plateau, there is a power law and the plateau ends with another power law prior to enter into the so-called full $`\alpha `$ regime. The power-law region is called $`\beta `$ correlator in the framework of $`MCT`$. The $`\alpha `$ decay is quite well phenomenologically described by a stretched exponential.

We have briefly sketched the main features of the MCT for colloids considered as hard spheres and proceed to consider the novel consequences of adding an attractive contribution to the inter-particle potential. In this paper we mainly consider the effects of using a hard core followed by a square well potential in order to mimic the interactions between colloidal particles. In real systems this can be obtained, for example, by covering the surface of the particles with a polymer coating or with depletion interaction . A description of such systems using the MCT approach has revealed the existence of a set of new and interesting phenomena, that we briefly summarise in what follows . In the temperature $`T`$ and volume fraction $`\varphi `$ plane MCT predicts two lines of transition from liquid to glass.One of the lines extends to high temperatures, it can be traced essentially to the repulsive part of the potential and tends asymptotically for high temperatures to the value $`\varphi `$ corresponding to the critical value for hard spheres. We will call it the repulsive glass transition line (RGL). On lowering the temperature the transition line moves toward higher values of $`\varphi `$ and gives rise to a re-entrant behaviour, i.e. one can pass from a liquid to a glassy state either by lowering or by raising $`T`$. In other words the liquid phase tends to exist in a region which extends more deeply into the glassy phase compared to the hard spheres case. The origin of this unusual effect is that if the interaction is short ranged enough, there can, in this temperature and density regime, be partial cancellation of the repulsive and attractive interactions leading to a phenomenon for glasses, not unlike that of a theta point for polymers.

The other glass transition line, originating in a relatively well-defined energy scale of the well-depth, is almost parallel to the $`\varphi `$ axis and extends on one side of the binodal to the other, until it eventually crosses the repulsive transition line. We call it the attractive glass transition line (AGL). For sufficiently narrow well-widths and on the high $`\varphi `$ end the AGL enters the glass region and terminates in an endpoint. Thus the two sides of the line separate two different types of glassy system that we will characterise through their mechanical properties. The endpoint of the AGL is a higher order glass transition point (called $`A_3`$), which will be shown to have special dynamical properties, namely the relaxation slows down dramatically on approaching it. Upon increasing the width of the potential well, the length of the AGL shrinks and its extension into the glass region reduces to zero, and as a consequence the endpoint turns into a higher transition point (called $`A_4`$). We will characterise this rather complex behaviour of the system, where repulsion and attraction compete, using two typical dynamical quantities, the shear viscosity, which characterises the liquid phase, and the elastic shear modulus for the amorphous glass.

Aspects of the MCT have been tested in many different cases, both experimentally and with computer molecular dynamics . The most extensive and accurate experimental check has been performed in the liquid-glass transition of colloidal systems, treated as hard sphere systems, and studied with dynamic light scattering . The agreement with MCT is quantitatively satisfactory . In a similar fashion computer simulations have been used to study the glass transition in simple model systems, as diverse as Lennard-Jones binary systems , the SPC/E model for water , orthoterphenil  and silica . In all these cases evidence shows that these systems undergo a kinetic glass transition, and the molecular dynamics is well accounted for by the idealised MCT of super-cooled liquids. On the other hand there is extensive evidence, mainly due to experimental results in colloidal systems, that cannot be interpreted in terms only of hard-core potentials.

Dense systems of colloidal particles characterised by a hard core and strong short-ranged attractions have been realized experimentally by adding polymers to either a suspension of colloidal hard spheres , in solutions of sterically stabilised particles when decreasing the solvent quality , and in copolymer micellar systems when changing the temperature . Such systems were also studied by Monte Carlo simulations . The facts that cannot be simply explained using hard-core potentials are the following.

(i) An amorphous material can be formed by increasing the attraction strength even though the volume fraction is kept well below the value of the hard sphere glass transition .

(ii) In mixtures of colloids and polymers, melting of the glass states is obtained by increasing the strength of a short range attraction by the addition of small polymers .

(iii) Using solvents of decreasing quality  in polymer coated colloidal particles the long time limit of the density time correlation function at small wave vectors are much larger than in hard-sphere systems.

(iv) Viscoelastic measurements for intermediate frequencies found strongly concentration dependent elastic moduli .

(v) A system in which an anomalous dynamical behaviour has been reported is a polymer solution (called L64) where structural arrest accompanied by a logarithmic-like decay of the density correlator is observed .

The paper is organised as follows. In Section 2 we report the details of the calculation of the structure factor for a square-well system (SWS), the only input needed for the mode-coupling dynamics, which is described in Section 3. The resulting phase diagram is described in Section 4, the time dependent density correlation functions in Section 5 and finally the various phases are characterised by their mechanical properties in Section 6. In Section 7 we report our conclusions.

## 2 Calculation of the structure factor

In the framework of the MCT the only quantity that is needed in order to calculate the dynamical properties of a super-cooled liquid approaching the glass transition is the structure factor. It is well known that, given an inter-particle potential satisfying the Ornstein-Zernike equation for the space-dependent pair correlation function $`g(r)`$, one needs a closure approximation in order to solve the equation. Many approximations of this type have been proposed and solved in the case of simple model potential . In the case of the square-well potential the Wiener-Hopf method as formulated by Baxter  can be used. The interaction potential $`V(r)`$ for particles a distance $`r`$ apart, is obtained with a hard-core repulsion for $`r<d`$, and the negative attractive value $`u_0`$ within the range $`d<r<sD`$.

The structure factor is specified by three control parameters, i.e. the packing fraction $`\varphi `$ related to the hard cores, the temperature $`T`$, and the relative width $`ϵ=1d/D`$ of the attraction shell. In the Baxter approach the structure factor $`S_q`$ is expressed in terms of the Fourier transform $`\stackrel{~}{Q}(q)`$ of the factor function $`Q(r)`$, a continuous real function defined for $`r0`$,

$`S_q^1`$ $`=`$ $`\stackrel{~}{Q}(q)\stackrel{~}{Q}^{}(q),`$ (1)
$`\stackrel{~}{Q}(q)`$ $`=`$ $`12\pi \rho {\displaystyle _0^{\mathrm{}}}𝑑r\mathrm{exp}(iqr)Q(r).`$ (2)
The funcion $`Q(r)`$ is related to the direct correlation function $`c(r)`$  and both functionss vainsh beyond the distance $`D`$. For $`0rD`$, the following equation for $`Q(r)`$ holds

$$rc(r)=Q^{}(r)+2\pi \rho _r^R𝑑sQ^{}(s)Q(sr).$$
(3)
In addition for $`r>0`$

$$rh(r)=Q^{}(r)+2\pi \rho _0^R𝑑s(rs)h(|rs|)Q(s).$$
(4)
For the SWS, $`g(r)=0`$ is fulfilled for $`0<r<d`$, and therefore, using $`h(r)=g(r)1`$, equation (4) splits into three sub-equations. The result for the middle part, $`\mathrm{\Delta }rd`$, where $`\mathrm{\Delta }=Dd`$, is very simple since the formula known from the theory for the hard-sphere system is reproduced

$$Q^{}(r)=ar+b,$$
(5)
where the coefficients $`a`$ and $`b`$ are introduced by

$$a=12\pi \rho _0^{d+\mathrm{\Delta }}𝑑sQ(s),b=2\pi \rho _0^{d+\mathrm{\Delta }}𝑑ssQ(s).$$
(6)
Defining $`G(r)=rg(r)`$ one finds, for small distances, $`0r\mathrm{\Delta }`$, where $`\mathrm{\Delta }=Dd`$,

$$Q^{}(r)=ar+b2\pi \rho _{r+d}^{d+\mathrm{\Delta }}𝑑sG(sr)Q(s),$$
(7)
and for the attraction shell, $`drD`$, one obtains

$$Q^{}(r)=ar+bG(r)+2\pi \rho _0^{rd}𝑑sG(rs)Q(s).$$
(8)
A closure approximation for $`c(r)`$ must be introduced into equation (3) in order to complete the system of equations (3) and (8). Using the Percus-Yevick (PY) approximation , according to which

$$c(r)=g(r)\left[1exp\left(\frac{V(r)}{k_BT}\right)\right]$$
(9)
inside the hard core in equation (3) and equation (8) leads to the approximation valid for $`drd+\mathrm{\Delta }`$

$`e^{u_0/k_BT}G(r)=`$ $`ar+b2\pi \rho {\displaystyle _r^{d+\mathrm{\Delta }}}𝑑sQ^{}(s)Q(sr)`$ (10)
$`+2\pi \rho {\displaystyle _0^{rd}}𝑑sG(rs)Q(s).`$
Equations (8) and (10) for $`Q(r)`$ and $`G(r)`$ are solved numerically choosing on each of the three intervals of the variable a grid of equally spaced points $`r_n`$, where $`n=1,2,\mathrm{},1000`$. The integral in equation (2) is then performed to obtain $`\stackrel{~}{Q}(q)`$ and hence $`S_q`$. The typical values for the factor function $`Q(r)`$ are shown in figure 1, while figure 2 shows the corresponding structure factors $`S_q`$ for various temperatures, given in term of the well depth $`u_0`$ as $`k_bT/u_0`$.

In a first attempt to solve the problem associated with the attractive part of the inter-particle potential the Baxter potential has been used , i.e. the limit of a square-well in which the range of the potential vanishes while its depth increases in such a way that their product remains constant . The mean spherical approximation has also been used in conjunction with an attractive Yukawa potential in an attempt to explain colloidal gelation . We may note that more recently other methods of generating the structure factors have been explored . The main results reported here are found not to depend on these details, and we therefore here only comment on the square well potential.

## 3 The mode-coupling equations

The MCT of the ideal liquid-glass transition is capable of describing the ergodic to non-ergodic transition in the normalised density correlation function

$$\varphi _q(t)=\frac{\rho _\stackrel{}{q}^{}(t)\rho _\stackrel{}{q}}{|\rho _\stackrel{}{q}|^2}.$$
(11)
In the long time limit $`\varphi _q(t)`$ shows a discontinuous variation from a vanishing value in a super-cooled liquid to a finite value different from zero $`f_q`$ in a glassy state, on varying the external control parameters. In the latter state the system experiences a structural arrest. In real systems the transition point is never reached since new dynamical mechanisms set in which bypass the ergodicity breakdown. The MCT equations are

$$\ddot{\varphi _q}(t)+\mathrm{\Omega }_q^2\varphi _q(t)+\nu _q\dot{\varphi _q}(t)+\mathrm{\Omega }_q^2_0^tm_q(tt^{})\dot{\varphi _q}(t^{})𝑑t^{}=0$$
(12)
where the mode-coupling memory functional $`m_q`$ is given by

$$m_q(f)=\frac{1}{2}\frac{d^3k}{(2\pi )^3}V_{\stackrel{}{q},\stackrel{}{k}}f_kf_{|\stackrel{}{q}\stackrel{}{k}|}.$$
(13)
The mode-coupling vertices are determined by the density $`\rho `$, the structure factor $`S_q`$ and the direct correlation function $`c_q=(11/S_q)/\rho `$

$$V_{\stackrel{}{q},\stackrel{}{k}}S_qS_kS_{|\stackrel{}{q}\stackrel{}{k}|}\frac{\rho }{q^4}\left[\stackrel{}{q}\stackrel{}{k}c_k+\stackrel{}{q}(\stackrel{}{q}\stackrel{}{k})c_{|\stackrel{}{q}\stackrel{}{k}|}\right]^2.$$
(14)
The two quantities $`\mathrm{\Omega }_q`$ and $`\nu _q`$ are respectively the characteristic frequency of the phonon-type motions of the fluid, and a term that describes instantaneous damping arising from the fast contribute to the memory function. They are defined as

$`\mathrm{\Omega }_q`$ $`=`$ $`{\displaystyle \frac{q^2k_BT}{mS(q)}}`$
$`\nu _q`$ $`=`$ $`\nu _1q^2`$
and $`\nu _q=1`$ in our calculations. In the long time limit $`t\mathrm{}`$, the density correlators $`\varphi _q(t)`$ tend to a value

$$f_q=\frac{\rho _q^{}(0)\rho _q(\mathrm{})}{|\rho _q(0)|^2}$$
(15)
the non-ergodicity factor, or Debye-Waller factor. The MCT equations (12) in the static limit give rise to the bifurcation relation

$$\frac{f_q}{1f_q}=\frac{1}{2}\frac{d^3k}{(2\pi )^3}V(𝐪,𝐤)f_kf_{|𝐪𝐤|}$$
(16)
It is clear that $`f_q=0`$ is a solution of equations (16) and it corresponds to an ergodic state of the system in which the correlations decay for long time. The correlators tend instead to a finite value different from zero if the system is kinetically arrested. This loss of ergodicity for $`\varphi `$ is interpreted as the transition to a kinetic glassy state. Therefore for some critical values of the thermodynamic parameters, density and temperature in our case, bifurcations of the solutions of the asymptotic equations appear that produce non-zero solutions. The bifurcations can be multiple, up to the number of control parameters of the system. Thus, when a bifurcation gives rise to more than two solutions of equations (16), there will exist multiple solutions with finite non-ergodicity factors. In these cases one speaks of $`A_k`$-type bifurcations, with $`k=2,3,4,\mathrm{}`$ In this case, MCT predicts that only the state corresponding to the largest value of $`f_q`$ is a stable solution of the equations .

## 4 Phase diagram

The phase diagram resulting from the solution of the MCT equations in the large time limit is shown in figure 3. We report the liquid-glass transition lines for some values of the width parameter $`ϵ`$. At this point it is worth reviewing the points made earlier about the experimental observations in the introduction. Indeed, the calculated phase diagram exhibits the phenomena of points (i) and (ii), and in later sections we shall show that all the other points are satisfactorily described. We point out the salient features of the phase diagram below.

(i) For high enough temperatures the repulsive glass lines, RGL, tend to converge to a value approaching the critical volume fraction for hard-sphere systems $`\varphi _c0.516`$, corresponding to the fact that for large $`T`$ the short-range part of the inter-particle potential weakly affects the dynamics. In this case the microscopic effect described by the MCT is the excluded volume cage effect, the impossibility of the molecules to move freely due to the presence of the neighbouring particles.

(ii) On lowering the temperature a liquid stabilising effect due to the attractive forces sets in, so that the liquid-glass lines tend to extend to larger values of the volume fraction. This behaviour gives also rise to the possibility of the glass melting when lowering the temperature.

(iii) The reentrant behaviour characterises the interplay of the repulsive and the attractive inter-particle forces; the system can become a glass both when lowering or raising the temperature.

(iv) For even lower values of the temperature another phase line, the AGL, appears and it is almost parallel to the $`\varphi `$ axis. In this case the particle motions are difficult since, due to the attractive interactions, the particle tend to form bonds.

(v) The two sections of the phase lines, RGL and AGL, intersect at a finite angle for small values of the parameter $`ϵ`$. The AGL extends for higher volume fractions beyond the crossing point and finally ends in a point.

(vi) The points on the glass lines we have discussed so far all correspond to bifurcations of type $`A_2`$, while the end point of the AGL relates to a higher order bifurcation $`A_3`$.

(vii) The extension of the AGL that develops in the glass region leads to the coexistence of two different types of glasses. In order to represent the differences between the two glasses we present in figure 4 the shape of the non-ergodicity factor $`f_q`$ crossing the glass-glass transition. The results are now consistent with the experimental observation (iii) in the introduction. The width of the $`f_q`$ is related to the localisation length of the particles in the glass. If the width is larger the particles are more localised. It is evident that the colloidal particles are more localised in the attractive glass then in the repulsive one. It is also possible to show that the localisation of the length in the repulsive glass remains more or less unchanged decreasing the temperature, whereas for the attractive one the particles become more localised. This seems to confirm the idea that the repulsive glass is dominated by the so-called cage effect: when the system gets to high packing fraction the particles start to be blocked by their neighbours, at certain critical packing fraction $`\varphi _c`$ each particle is trapped in a cage formed by the surrounding particles. Indeed this is a purely geometrical effect and it should not depend on temperature.

(viii) On increasing the relative square well width parameter $`ϵ`$ the glass-glass line tends to shrink and finally coincides with the point of intersection of the RGL and AGL when $`ϵ0.04`$. This end-point is a higher bifurcation point of type $`A_4`$. We note that there is no difference in the structure factor across this curve. Only the late stage dynamics, reflected in the values of the non-ergodicity factors, may be used to differentiate between these glasses.

## 5 The intermediate scattering function

We have studied the dynamical behaviour of the system in different regions of the phase diagram, in particular in the proximity of the higher order bifurcation points $`A_3`$ and $`A_4`$. We will mainly comment on the case where the well width parameter is $`ϵ=0.03`$. The wave vectors are given in units of the particles diameter. In figure 5 we present different intermediate scattering functions for $`q=21.75`$, obtained by solving MCT equations, approaching the glass transition at $`\varphi =0.504449`$ and $`T_c0.8000`$. The relevant cut across the attractive glass curve is shown in figure 3 using crosses, and a transition point. For temperatures well above the glass transition, i.e. $`T=0.85`$, $`\varphi _q(t)`$ presents the typical liquid behaviour: after a first short-time decaying due to microscopic dynamics, the correlation function starts to relax toward ergodicity and no sign of a critical slowing down is evident. When the system becomes close enough to the transition temperature ($`T=0.80010`$ in figure 5) the system starts to freeze and the typical two step relaxation scenario starts to be evident. We note that we are here crossing that transition curve which we believe corresponds to the process of gelation at high density. We see, therefore, that we may expect the characteristic behaviour of a glass transition in the slow dynamics near gelation at high density.

Once the system is below the transition temperature (the case $`T=0.80000`$ in figure 5), $`\varphi _q(t)`$ corresponding to the gel does not relax to zero anymore, the system is non-ergodic. In this situation there is only part of the $`\beta `$-relaxation toward a non-zero constant value, i.e. the non-ergodicity parameter $`f_q0`$. Perhaps it is useful now to discuss some of the MCT predictions for the asymptotic behaviour of the intermediate scattering function (see for example  and references therein). In the time region where $`t_0t\tau _\alpha `$ the so-called factorisation theorem holds

$$\varphi _q(t)=f^c(q)+h_qG(t)$$
(17)
$`t_0`$ being a typical macroscopic time, $`\tau _\alpha `$ the $`\alpha `$ relaxation time scale and $`f_q^c`$ the critical non ergodicity parameter. Equation 17 implies a universal behaviour in which wave-vectors and time factorize. $`G(t)`$ is a scaling function which describes the whole relaxation pattern. Near the glass transition the behaviour of the function $`G(t)`$ may be calculated analytically in the proximity of the $`\alpha `$ and $`\beta `$ relaxation. It can be shown that the asymptotic behaviour can be expressed in terms of a rescaled time $`\tau =t/t_\sigma `$ where $`t_\sigma `$ is a time scale that depends crucially on the distance to the transition, i.e. on the separation parameter $`\sigma `$ . Close to the glass transition, for $`T>T_c`$, $`G(\tau )`$ is given by

$`G(\tau 1)1/\tau ^a`$ (18)
$`G(\tau 1)\tau ^b`$ (19)
where the two exponents are related to the exponent parameter $`\lambda `$, which is obtained from a stability matrix , via the relation $`\lambda =\mathrm{\Gamma }^2(1a)/\mathrm{\Gamma }(12a)=\mathrm{\Gamma }^2(1+b)/\mathrm{\Gamma }(1+2b)`$, which implies $`0<a<1/2`$ and $`0<b<1`$. The scaling behaviour of the correlation function changes in the proximity of higher bifurcation points. In particular it is possible to show that near an $`A_3`$ point $`G(t)1/\mathrm{ln}(t/t_1)`$ whereas near an $`A_4`$ it goes like $`G(t)1/\mathrm{ln}^2(t/t_1)`$, where $`t_1`$ is a microscopic time scale . This is quite a remarkable behaviour and represents a rather strong prediction of the theory. There have been reports of such behaviour in some system , and it will be interesting to see if this phenomenon is discovered in a variety of different colloidal systems.

## 6 Mechanical properties

Two particularly interesting quantities to investigate in the vicinity of the glass transition are the shear viscosity and the elastic shear modulus. Those two quantities can be measured in real systems and such measurements can provide a good insight in the glass transition phenomena and its relations with MCT. The complex shear viscosity $`\eta ^{}(\varphi ,T,\omega )`$ can be calculated in terms of the normalised time correlation function of the density fluctuations $`\varphi _k(t)`$

$$\eta ^{}(\varphi ,T,\omega )=\frac{k_BT}{60\pi ^2}_0^{\mathrm{}}𝑑te^{i\omega t}_0^{\mathrm{}}𝑑kk^4\left[\frac{dlnS_k}{dk}\varphi _k(t)\right]^2$$
(20)
and is related to the complex shear modulus $`G^{}(\varphi ,T,\omega )`$ by

$$G^{}(\varphi ,T,\omega )=i\omega \eta ^{}(\varphi ,T,\omega ).$$
(21)
It is possible to take the $`\omega 0`$ limit in equation 20

$$\eta ^{}(\varphi ,T,\omega =0)=\frac{k_BT}{60\pi ^2}_0^{\mathrm{}}𝑑kk^4\left[\frac{dlnS_k}{dk}f_k\right]^2$$
(22)
and this gives the static shear viscosity. We note that this formula, in that it depends on the square of the correlation function, is different in structure from that used by Weitz and co-workers, where the stress-relaxation function is linearly related to the density correlation function . This difference would lead, in principle, to substantially different predictions for the loss modulus, and in particular around its minimum. It seems likely that the idea the stress modulus is linear in the correlation function is more suited to visco-elastic properties dominated by large domains, and their changing surface area under shear , whilst the dependence on the square of the correlation function is more suited to dense fluids where the properties are dominated by loss of mobility due to caging effects.

All the quantities in the equations (20), (21) and (22) can be evaluated for the present system, and therefore it is possible to use the numerical solution of equation (12) in order to solve them. This has been accomplished performing the integrals with standard numerical integration over $`2000`$ $`q`$-vector varying between $`0`$ and $`500`$. In figure 6 the static shear modulus has been represented for three different packing fractions ($`\varphi =0.539672`$, $`\varphi =0.544052`$ and $`\varphi =0.580000`$) as a function of the temperature. These three cases correspond to crossing the glass-glass line, crossing the $`A_3`$ point and a density where the repulsive and attractive glasses have become indistinguishable. In the first case it is clearly possible to distinguish between the two glasses. For low temperatures there is strong dependence of the elastic viscosity on the temperature, in particular the system becomes more and more rigid on decreasing the temperature. When the system crosses the glass-glass transition there is a discontinuity in the elastic response which clearly indicates that the structure is changed. Increasing further the temperature the elastic behaviour does not change so much anymore and for large temperature the system behaves like a hard spheres suspension. In this case the glass is originated by the cage effect and consequently the particle are forced to move inside a fixed volume that does not change with temperature. It is then evident how the differences between the two glasses are manifested by mechanical properties. We have noted that the difference in the shear modulus as the $`A_3`$ point is approached is described a by a power law, and it is also possible to show that in this regime the differences in mechanical properties of the two systems are due to the different shape of the non-ergodicity parameter $`f_q`$, and hence the long time residual motions in the gel, and not to the contribution of the equilibrium structure factor .

We have performed this calculation also for the dynamical quantities defined by (20) and (21). It is evident that, if the system is in the glassy phase we have for the long time behaviour of the density correlators $`\varphi _q(t\mathrm{})=f_q0`$, and consequently the integral in (20) is unbounded. Thus the solution posses a zero frequency pole. In order to obtain convergent solutions for a glass we have replaced $`\varphi _q(t)`$ with $`\widehat{\varphi }_q(t)=\varphi _q(t)f_q`$ which decays to zero at infinite time. In figure 7a and figure 7b the real and the imaginary part of the complex viscosity, $`\eta ^{}(\omega )`$ and $`\eta ^{\prime \prime }(\omega )`$, are presented for a fixed volume fraction, $`\varphi =0.504449`$, for different temperatures, approaching and crossing the glass transition. In figures 8a and 8b $`G^{}(\omega )`$ and $`G^{\prime \prime }(\omega )`$ are shown for the same points as figure 7. Since figures 7 and 8 contain the same information we discuss only figure 8. $`G^{}(\omega )`$ represents the elastic response of the material: the higher its value the stiffer the material at that frequency scale. The imaginary part of the shear modulus describes the viscous behaviour of the system and so the dissipation. Both quantities can be easily measured experimentally, for example as the response of a sample to small oscillatory shears which weakly perturb the system . Approaching the glass transition we note that the $`G^{}(\omega )`$ begins to develop a shoulder at low frequency. This effect is due to the slowing down of the dynamics of the system, in other words the formation of the plateau in the $`\varphi _q(t)`$ is responsible for the formation of a region where $`G^{}(\omega )`$ varies only slightly with frequency. The range of such a region tends to increase approaching $`T_g`$. Indeed a similar behaviour has been observed experimentally in measurements of linear viscoelasticity in a colloidal suspension with hard sphere interaction . If the system is in a glassy state the value of the real part of the static shear modulus $`G^{}(\omega =0)`$ is finite, indicating that the system is solid and consequently it presents an elastic behaviour. In the liquid phase, however, the system does not show any elastic behaviour so the $`G^{}(\omega =0)`$ tends to zero. It is then clear that at the glass transition the system abruptly changes its behaviour, presenting a singularity in the static shear modulus. In figure 8b the behaviour of the imaginary part of the shear modulus is represented. At high frequency the curves for different temperature are the same, showing the microscopic dynamics which is the same for all the temperatures. Approaching the glass transition a second maximum starts to emerge. Such a maximum represents the $`\alpha `$-relaxation and it moves toward low frequency on decreasing the distance to the glass transition temperature $`T_g`$. The minimum of $`G^{\prime \prime }(\omega )`$ corresponds to the plateau region in time, and the power law behavior in frequency on both sides has been already observed in hard spheres systems .

## 7 Conclusions

In this paper we have outlined some of the experimental observations that have been made of particle gels in colloidal systems where we know there to be a strong short-ranged potential. We have shown that all of these features may be interpreted within the paradigm that such particle gels represent a new type of attractive glass. However, having made this correspondence, we then find a number of predictions of new phenomena that have yet to be confirmed experimentally, or are just being so. These include the re-entrant behaviour for the phase diagram, the logarithmic dynamic behaviour near this re-entrant regime, and in the general the characteristic alpha relaxational behaviour near gelation. We have calculated the mechanical manifestations of these phenomena also. Thus we have shown how the glass-glass transition may be studied using the shear modulus, and that the particle gel must be expected to be much stiffer than the colloidal glass, but that the differences between them vanishes near the $`A_3`$ point. We have also studied, to our knowledge for the first time, the frequency dependent modulus and loss modulus for these systems, using the paradigm that the system is undergoing a glass transition. MCT type treatments seem to be suited to study colloidal systems of this type because, at least of the repulsive system, the glass transition seems to lack some of the fluctuations of molecular liquids. We also can calculate phase diagrams, and mechanical and scattering properties from a single theory. We expect that, in future, dense particle gels will be better fitted to the theory present here, than previous phenomenological treatments.

There remains much to be done in many directions. It now begins to seems likely that the glass paradigm is suitable for interpretation of particle gelation. However we must await further detailed experimental work to see just how extensive the agreement with MCT type treatments will be. We may also speculate that this type of idea will have much broader significance that particle gels. There is every reason to hope that polymer gelation and protein gelation may also find a description in similar terms. Certainly this is a rich area for future exploration.

P.T. thanks Sow-Hsin Chen, to whom this paper is dedicated, for the invaluable example and the scientific collaboration lasting for more than three decades, and still active. K.D. and F.S. also thanks Sow-Hsin Chen for many stimulating conversations over a period of years, and spanning a number of complex fluid fields. We all appreciate his continued example in commitment to the scientific endeavour after so many years. The authors also acknowledge the collaboration with the group of Prof. Wolfgang Göetze, together which part of the work reported here was peformed. The work in Rome is supported by PRIN-2000-MURST and PRA-HOP-INFM, and the work both in Rome and in Dublin is supported by COST P1.

## References


text:
# On topological defects

Abstract The language and methods of algebraic topology, particularly homotopy theory, have been extensively used in the study of the identification, the classification and the evolution of defects. Topological methods provide the means for the detection of gross features, such as holes in a manifold, and therefore one would clearly anticipate that the defects, if identifiable by such means, to be linked with that kind of configurations. In physical media, however, defects usually manifest themselves in ways which have no distinct and direct association with patterns of topological value along the lines of our previous reference. In this paper we will postulate that the defects do not themselves correspond to topological features. Instead, they are forced to exist in order to redress the violation of physical rules happening as a result of the topological nature of our system. We will, thus, assume that the topology of our system changes with the introduction of defects. Therefore, the relevent to our study topological methods will be applied to a well defined physical system only in the absence of defects.
An ordered medium corresponds to some real smooth manifold $`M`$ where a function $`f`$ is defined which assigns to every point of the medium an order parameter. Furthermore, over our medium, $`f`$ should be a continuous function. For example, if $`f:^3`$ is a continuous function throughout $`^3`$ then our medium is $`^3`$ and the order parameter is a complex number. In our discussion, we will refer to that function $`f`$ by the name order parameter, even though this title might be somehow confusing111one, for example, might consider the space of functions as that relevant to an ordered medium and thus start all sorts of investigations on a different sphere.

Excluding the trivial case where the order parameter is constant throughout the medium (which is called thus uniform), we will focus our attention to non uniform media where the function, through connected space, varies continuously apart perhaps (depending on the specific configuration) at isolated regions .

Suppose $`f`$ to be an order parameter. Take $`f:\stackrel{~}{M}`$ to be everywhere continuous apart from some region $`M_d\stackrel{~}{M}`$. Obviously, according to our choice of $`f`$ we can force the existence of a certain $`M_d`$ in $`\stackrel{~}{M}`$. We could consider the $`M_d`$ region as a hole in $`\stackrel{~}{M}`$ so that $`f`$ is kept always continuous on a certain area of its domain. The latter we will denote as $`M`$: $`\stackrel{~}{M}M_d`$. In this way, the order parameter will be assigning topological properties to our medium $`M`$. However, this is an ambiguous technique in case $`f`$ is defined throughout $`\stackrel{~}{M}`$ since we will be just ignoring certain values of the order parameter and thus be lead to a non equivalent situation. As we will see, discontinuities in $`f`$ can be appropriately smoothed out so that any measures of the kind outlined already will be unnecessary.

Primarily, our contemplation to consider order parameters with discontinuities might seem quite far fetched to someone having in mind the situation in physical systems where the order parameter is always continuous. Its relevance will become evident if we introduced the notion of a defect. We will postulate that a defect corresponds to an appropriate change in the order parameter, applied only to those regions where originally it was discontinuous, having as a main aim to restore the continuity of that function throughout the area where it was originally defined. Thus, the discontinuities will indicate where the introduction of a defect is needed. In this way our medium will become the same as the domain of the order parameter diminishing the possibility of the previously mentioned ambiguity.

As we will see, for a certain order parameter behavior, the defect induces an invertible change to that function. Thus, if we were to lift the defect we would surely recover the underlying discontinuity. However, there are situations for which the discontinuity happens artificially, and thus, unless we remember the exact formula for the order parameter, any change intended to make that function continuous can not be reversed. In physical systems where the order parameter is everywhere continuous, defects are, possibly, already there and therefore the continuity has been appropriately restored. If, though, we were to remove those defects and recover the order parameter configuration, we would surely see those discontinuities for which this procedure is invertible.

It is, thus, inevitable to ask, how one identifies defects in continuous configurations. Moreover, how one knows how the order parameter should behave if one were to take them away. Here, we, thus, have to introduce the notion of a preferred order parameter. Since our medium is a physical system, it should be described by certain equations the solution to which is our order parameter. For various systems there exist many solutions to the relevant equations (and here one should consider all possible solutions including those having the order parameter constant). Each such order parameter corresponds to a particular energy distribution over our medium. We will assume that the solution preferred will be the least energetic one. There could be situations where several order parameters will be equally good in terms of their associated energy distribution. We will postulate that our physical system does not distinguish between them.

Suppose that for a system in equilibrium (where we can safely identify preferred order parameters) one finds that the corresponding order parameter does not, over all, have a preferred configuration. If, in this particular case one finds that there are some regions in our medium where the order parameter deviates from the surrounding preferred configuration then one should think of these areas as possible candidates for defects. If, further, one lifts those configurations and seeks for a good replacement among those functions having some physically explained priority then, in case such an investigation turns out to be in vain, one can safely link these areas with a defect. Here we have assumed that the way one discards certain - otherwise preferred - order parameters for a given region has to do with whether they violate the continuity requirement for the total order parameter defined over all the manifold. Thus, order parameters which occur in contrast to the natural tendency of the system, should underline the existence of defects.

To make the point clearer consider the situation where there are two different order parameters which give the same energy distribution over all the medium. Suppose that $`f_1=c_1`$ and $`f_2=c_2`$ with $`c_1`$, $`c_2`$ constants. Our medium is one dimensional. Take that the over all our manifold order parameter is

$$f(x)=\{\begin{array}{cc}f_1\hfill & \text{if }x>0\hfill \\ c\hfill & \text{for }x=0\hfill \\ f_2\hfill & \text{if }x<0\hfill \end{array}$$

Obviously, at $`x=0`$ there will be a discontinuity of $`f`$ regardless of the actual value of $`c`$. Thus, that discontinuity will be hidden if we introduced a defect which would take us continuously from $`c_1`$ to $`c_2`$ or vice versa. The order parameter configuration which we could choose should satisfy the relevant equations for our system. This example, for those familiar with our subject, monitors the situation of the one dimensional kink. In that system, $`c_1`$ and $`c_2`$ correspond to the two minima of the relevant potential. We took that, in the absence of a defect, those two values for the order parameter would create a discontinuity at the $`x=0`$ point of our one dimensional manifold. Thus, introducing a defect there, would smooth out our function in accordance to the relevant equations, and form a continuous transition between those two minima.

Alternatively, if our relevant order parameter was

$$f(x)=\{\begin{array}{cc}f_1\hfill & \text{for }x>0\hfill \\ c\hfill & \text{at }x=0\hfill \\ f_1\hfill & \text{for }x<0\hfill \end{array}$$

for some constant $`cf_1`$ then that discontinuity at $`x=0`$, if removed by forcing, for example, $`f(0)=f_1`$, would not be recoverable since there would be no trace left to indicate what type of measure we took, and indeed if we took any measure, in order to make the order parameter continuous. In that particular situation we will assume that the domain of the order parameter and $`M`$ are already the same even though there is a formula dependent discontinuity at $`x=0`$. In such cases the $`M_d`$ could be considered to be the empty set. Surely, there are other ways of making the order parameter be continuous throughout $``$. We will consider that the process of smoothing-out should happen in a way that the least changes to $`f`$ would take place.

Hereafter we will consider that if the order parameter, around any point $`p`$ of the manifold, has such a configuration so that the $`limf`$ at $`p`$ will be well defined, then

$$f(p)=\underset{p^{}p}{lim}f(p^{})$$

regardless of the original value of the order parameter there. In the converse case where the $`limf`$ at $`p`$ does not return a unique value, then in order to restore the continuity - our primary aim - we need to smooth-out our function in an appropriate way around the point $`p`$. Our assumption here is that the order parameter will be exhibiting discontinuities only in the absence of defects. Thus, $`f`$ is continuous because of the appearance of those configurations.

The idea of invertible comes as a direct consequence of what has been already stated. If for a given $`f`$ and a given $`\stackrel{~}{M}`$ we find that on some $`M_d`$ there are discontinuities of that function being exhibited then, our first step is to call $`M:\stackrel{~}{M}M_d`$ our relevant domain. Then, we will try and see if we can apply any preferred order parameter configuration on $`M_d`$ so that $`M\stackrel{~}{M}`$ in the sense that the new order parameter will be a continuous function throughout a medium ($`M`$) equal its original domain of definition ($`\stackrel{~}{M}`$). If we are unable to find any such preferred order parameter that will smooth-out our $`f`$ completely then the introduction of a new order parameter configuration will be an invertible process. Otherwise, the additional alterations could not be recovered.

Moreover, if we are given an order parameter configuration, take, for instance, the former case described above (the kink), then we will be unable to find if there are any defects there unless we are provided with the information as to which are the order parameters that the system would prefer to have.

Assume the planar spin case, also quoted by Mermin. We remind the reader that the system of planar spins corresponds to a flat $`2D`$ medium at each point of which a $`2D`$ vector is being assigned via the appropriate order parameter. Mermin claims that if we removed the information about the order parameter from inside a certain disk on that plane, we would still be able to recover a discontinuity at the center of that region in case the relevant vector winds by $`2\pi n`$ with ($`n1`$) around it. That is obviously true if one is provided with the additional information that our order parameter should be of a non zero magnitude throughout that disk. If we were allowed to put to zero the magnitude of the order parameter (thus, the preferred order parameter could take the zero value) it is not necessary that we would be able to recover any discontinuity at the center. Thus, the identification of areas where a defect is needed is the result of the application of topological methods to a physical system which, we assumed, tends to acquire only preferred order parameters.

It might seem, though, oxymoron the fact that the system does not take care of the continuity of the order parameter itself, even though any discontinuous behavior of that function cannot be physically allowed. It is inevitable, thus, to require an answer to the question of how the system manages to create a situation against its natural tendency. One, thus, needs to clarify whether the order parameter takes a particular form because there is some intrinsic value to the possibly created discontinuities or whether this is just a trick to explain how the system chooses a state associated with a certain continuous function $`f`$.

First of all, in order to create a discontinuity of the kind we already outlined we need to have a system which has more than one preferred order parameters. Furthermore, to simplify matters, we will assume that those functions are bounded from above and below and that their range of values correspond to entirely different sets. Suppose now that our medium $`\stackrel{~}{M}`$ can be divided into independent sections denoted as $`S_k`$’s, for which we have

$$\underset{k=1}{\overset{m}{}}S_k\stackrel{}{\mathrm{\_}}\stackrel{~}{M}$$

and

$$B_{kk^{}}=S_kS_k^{}$$

a part of the common boundary of $`S_k`$ and $`S_k^{}`$ or, otherwise, the empty set $`\mathrm{}`$. We will take that all the $`S_k`$’s around a certain $`S_{k_0}`$ belong to the latter’s environment and do not affect the physics of that set. Thus, the $`S_{k_0}`$ can be thought of as an independent physical system subject to no, by the environment imposed, boundary conditions. If that happens for each $`S_k`$ of $`\stackrel{~}{M}`$ then the order parameter will be allowed to evolve autonomically in each different section of the manifold.

Even if our manifold $`\stackrel{~}{M}`$ is separated into physically independent regions, we will assume that the parameterization in each of those areas does not happen independently. Thus, there is a global chart $`\psi `$ assigning coordinates to each point of $`\stackrel{~}{M}`$. Further, the identification of preferred order parameters will result from minimizing the relevant potential. If the latter is given in terms of manifold coordinates, then the order parameter configuration will be parameterization dependent something which we wouldn’t like to have. The reason is that we want to allow the system to choose any order parameter at each $`S_k`$ for any given $`\psi `$. Thus, the potential will have to be given in terms of $`f`$, the order parameter. It is essential to stress that by minimizing that potential we will be finding order parameter values which can be applied in different ways to our manifold. That is to say that the way the order parameters we will be finding, should vary with manifold coordinates is not determined and therefore the realization of the potential over some $`S_k`$ can have many different forms for a given $`\psi `$.

Suppose that

$$\underset{k=1}{\overset{m}{}}S_k\text{=}\stackrel{~}{M}$$

and

$$B_{kk^{}}=S_kS_k^{}$$

a part of the common boundary of $`S_k`$ and $`S_k^{}`$. The boundary of every $`S_k`$ is being thought of as a wall that prevents information about the order parameter on its one side to cross to the other side. Assume that the way to create such a situation is by externally forcing the system to become fragmented.

Taking that the manifold is separated into independent regions each of which corresponds to a physical system in its own right seems accommodating. Since in each of those regions the order parameter is a continuous function, each physical system will be related to a function with no discontinuities within it. Thus, the whole $`\stackrel{~}{M}`$ cannot be considered as one physical system. The latter will be implemented when the boundaries between the various $`S_k`$’s fall and defects appear which smooth-out any discontinuities at the connecting borderlines.

This approach lets us assume that those discontinuities are happening not because of some intrinsic property of our physical system but instead because of external intervention. They are, thus, forced to appear because of the from the outside imposed fragmentation of our medium. This can be related to what Polturak , Carmi & Koren did in their recent experiment forcing a predetermined domain structure. In this way the defects will be created after that imposed structure is left to naturally evolve and smooth-out, in an appropriate way, any initial discontinuities.

Another way we can create discontinuities is by assuming that the medium consists of disconnected regions being at some distance $`D`$ away from each other. Thus,

$$\underset{k=1}{\overset{m}{}}S_k\stackrel{~}{M}$$

and

$$B_{kk^{}}=S_kS_k^{}=\mathrm{}$$

At each $`S_k`$, being thought of as an independent physical system, the order parameter takes a preferred configuration. One can postulate that those regions either grow or come closer so that they become correlated. Thus, the one feels the presence of the other, say after a distance $`d`$. Since then the order parameter in each of those regions cannot acquire a form independently. It could be possible that the form of the order parameter, beyond that point, will not anymore correspond to a preferred one for either of those regions. The defect, thus, will be appearing due to the interaction between those areas.

This situation could be related to what has been performed already in superfluid $`{}_{}{}^{3}He`$. Using the fact that $`{}_{}{}^{3}He`$ is an excellent neutron absorber, two experiments, one in Helsinki and another in Grenoble subject it to neutron bombardment to heat a small region of the superfluid. Small bubbles of normal fluid appeared which, by rapidly cooling down, evolved independently from the surrounding superfluid. The order parameter of the surrounding superfluid $`{}_{}{}^{3}He`$ could not follow the changing temperature front fast enough. Consequently, internal regions of the hot volume transit into the superfluid phase independently in accordance with the Kibble/Zurek cosmological scenario.

Here, thus, we have postulated that the main ingredient for the production of defects is to create independent regions within the manifold where the order parameter acquires some preferred form. Moreover, we assumed that in each one of those regions the phase transition takes place individually. Hence, the phase transition by itself will not create those domains in our manifold.

E.D.M.K. would like to thank Dr. R. J. Rivers for fruitful discussions and support.