text:
# 1 Introduction

## 1 Introduction

Liquids undercooled below the freezing temperature without crystallization shows a sharp increase in their transport properties like viscosity and acquire solid like properties e.g. has a finite shear modulus. The dynamics of the dense liquids are studied in terms of generalized transport coefficients, expressed in the form of memory functions which account for the correlated motion of the fluid particles at high density. Frequency and wave vector dependence of these quantities, relevant for behavior over different time and length scales are studied through models of the liquid obtained from the statistical mechanics of a many particle system. Much progress has been made in this field through a self-consistent approximation of the memory function in terms of the slowly decaying modes that exist as a consequence of microscopic conservation laws . This has been termed in the literature as the mode coupling theory (MCT) . A key result of the simple mode coupling model  for the memory function, is a two step relaxation process involving the power law decay of correlation over intermediate time scales crossing over to the stretched exponential behavior in the long time, termed as the $`\alpha `$-relaxation. With increase of density, the liquid become very viscous and relaxation slows down strongly. The basic driving mechanism for the enhancement of the viscosity comes from a feed back to the transport coefficient from the coupling of the slowly decaying density correlation functions. In the simple model, long time correlation of density fluctuations freeze completely beyond a critical density and the fluid undergoes a dynamic transition to a nonergodic phase. An extended version  of the mode coupling theory with proper coupling of the density and current modes demonstrated that there is a final exponential relaxation of the density correlation function over asymptotic time scales. The final decay process which shifts to longer time scales with density can be interpreted as a process similar to the diffusion of free volumes in the amorphous system whereby ergodicity is maintained.

In the theoretical formulation of the dynamics of the supercooled liquids often extra slow modes  in addition to the usual ones that follow from simple hydrodynamics, has been considered. Extension of the simple MCT has been done  along these lines to include the extra slow mode of defect density that develop in the amorphous solid like state. In Ref.  and  the displacement field $`\stackrel{}{u}`$ was used to describe the extra hydrodynamic mode due to the solid like nature of the system. Ref.  demonstrated that due to the nonlinear couplings of the density fluctuations and the displacement field $`\stackrel{}{u}`$ in a solid like state, the renormalized expression for transport coefficients include terms linear in density correlation function. Such terms in the memory function appears over time scales in which the defect correlation is taken to be a constant. Due to this phenomenological inclusion of the extra slow mode of the defect density, a line of transition is obtained instead of a single critical density. In this formulation of the extended MCT, temperature or density dependent exponents of power law relaxation  are obtained. Over the longest time scale the defect correlation as well as the density correlation function do however decay to zero and thus restores the ergodic behavior. In a recent work  we have studied the coupled dynamics over longest time scale, through a schematic model dropping all wave vector dependencies. The asymptotic behavior for both these correlation functions have been considered there to study the relaxation over various time regimes. However, in order to study the transport behavior with respect to change in the thermodynamic parameters like density, the wave vector dependence of the model has to be taken into account. In the present work we consider the wave vector dependent model for the dense supercooled liquid and the coupling of the density fluctuations to the extra slow mode of defect density are included. The dynamics of the density correlation function is considered up to the time scales where the final ergodicity restoring mechanism present. It was demonstrated  from a non perturbative approach that the self consistent feed back mechanism allows a non-zero value of the cutoff function. However, explicit evaluation of the strength of the function is extremely difficult beyond the one loop order. On the other hand in certain simple systems like colloids very slow relaxation of the structure function has indeed been observed. This indicate that the cutoff function is effective over very long time scales. In the present calculation we model this through a self-consistent picture where the defect density and the density fluctuations over longest time scales behave in a similar way. We use here this criteria to close the dynamical equation for the density correlation function - instead of actually evaluating cutoff function self consistently. The cutoff function is estimated to give decay of the structure factor that agrees with experimental results on hard sphere like systems like Colloids. We then compute the mode coupling integral to obtain the renormalized transport coefficient as function of density. The coupling of the slowly decaying density fluctuations to the defect density in the super cooled state is a key aspect of this model and allows a physical picture for the ergodicity restoring mechanism.

The plan of the paper is as follows. In the section II we give a brief description of the model studied and the approximations involved while in section III we describe our results for the transport coefficients. We end the paper with a small discussion.

## 2 Transport Coefficients at High Density

### 2.1 Longitudinal Viscosity

In this section we describe briefly the model used in computing the longitudinal viscosity or the sound attenuation constant. The key quantity of interest is the dynamic structure factor or the density auto correlation function (DACF). The Fourier-Laplace transform of the DACF $`\psi _q(t)`$ normalized with respect to its equal time value , defined as,

$$\psi _q(z)=i_0^{\mathrm{}}e^{izt}\psi _q(t)𝑑tIm(z)>0$$
(1)
In the basic formulation of the MCT, this quantity is expressed in terms of the generalized memory function $`\mathrm{\Gamma }_q(z)`$,

$$\psi _q(z)=\frac{z+i\mathrm{\Gamma }_q(z)}{z^2\mathrm{\Omega }_q^2+i\mathrm{\Gamma }_q(z)[z+i\gamma _q(z)]}.$$
(2)
$`\mathrm{\Omega }_q`$ corresponds to a characteristic microscopic frequency for the liquid state dynamics. The memory function which is the generalized longitudinal viscosity, consists of the bare and mode coupling part, $`\mathrm{\Gamma }_q(z)=\mathrm{\Gamma }_q^0+\mathrm{\Gamma }_q^{mc}(z)`$. $`\mathrm{\Gamma }_q^0`$ is related to bare or short time dynamics with uncorrelated collisions. The mode coupling contribution $`\mathrm{\Gamma }_q^{mc}(z)`$ signify the correlated motion in the dense liquid and can be expressed in terms of the coupling of density fluctuations dominant at supercooled densities. Several models has been used to compute the Generalized memory function for the supercooled state. The simplest model . express the memory function in terms of a bilinear product of the density correlation functions i.e. $`\mathrm{\Gamma }^{mc}(t)`$ $`=`$ $`c_2\psi ^2(t)`$ . The model that contains a linear term of the density correlation function in addition to the quadratic one, i.e. $`\mathrm{\Gamma }^{mc}(t)`$ $`=`$ $`c_1\psi (t)+c_2\psi ^2(t)`$, has been referred to in the literature  $`\varphi _{12}`$ model. In the wave vector independent form this model gave rise to the stretched exponential relaxation. In later works it was demonstrated that by considering the solid like nature of the amorphous state at high density such a linear term can be obtained in the renormalization of the transport coefficients. The set of hydrodynamic variables were now extended to include the slowly relaxing defect density. The resulting non linearities in the equations of Generalized Hydrodynamics gives rise to the above memory function in terms the DACF. In order to consider the very dense system approaching the glass transition we will use here this model for the dynamical behavior. Following Ref. , the mode coupling part of the generalized transport coefficient is given by,

$$\mathrm{\Gamma }_q^{mc}(t)=\frac{d\stackrel{}{k}}{n_0(2\pi )^3}S(q)[\stackrel{~}{V}^{(1)}(q,k)\psi _k(t)\varphi _{k_1}(t)+\stackrel{~}{V}^{(2)}(q,k)\psi _k(t)\psi _{k_1}(t)]$$
(3)
with $`\stackrel{}{k_1}=\stackrel{}{q}\stackrel{}{k}`$. Here $`\varphi _{k_1}(t)`$ is the normalized defect density auto correlation function and gives the contribution to the transport coefficient coming from the coupling of the defect density with the density fluctuations. The mode coupling vertices with full wave vector dependence are given by

$$\stackrel{~}{V}^{(1)}(q,k)=2y\stackrel{~}{U}(q,k)S(k)S(k_1)+\kappa n_0c(k)S(k)$$
(4)
and

$$\stackrel{~}{V}^{(2)}(q,k)=\frac{1}{2}\stackrel{~}{U}^2(q,k)y\stackrel{~}{U}(q,k)S(k)S(k_1)$$
(5)
The vertex function $`\stackrel{~}{U}`$ in (5) is given by,

$$\stackrel{~}{U}(q,k)=\frac{n_0}{q}[\widehat{q}.kc(k)+\widehat{q}.(k_1)c(k_1)]$$
(6)
where $`n_0`$ is the particle number density. It is assumed that the defects are moving in a metastable double well potential and this well depth is characterized by $`y`$. The defect density is assumed here to be a variable similar to mass density and is weakly interacting with the latter. $`\kappa `$ is a dimensionless parameter characterizing the coupling of the defect density with the particle density in the Free energy functional corresponding to the stationary state around which fluctuations are considered. For $`\kappa `$ and $`y`$ both equal to zero, (3) reduces to the simple MCT result, . If the defect correlation is assumed to be very long lived compared to the time scales over which initial power law relaxation persists, $`\varphi _{k_1}(t)`$ in equation (3) can be taken as a constant and the so called $`\varphi _{12}`$ model results on ignoring all wave vector dependencies. The density dependence of the power law exponents and the resulting dynamic instability in the model has been studied in Ref. . Here we will consider the time scale over which both density correlation and the defect correlation decays. The quantity $`\gamma _q(z)`$ in the R.H.S of eqn. (2) has crucial implications for the asymptotic dynamics . If $`\gamma `$ is ignored then a sharp transition in the supercooled liquid to an ideal glassy phase results beyond a critical density. The Laplace transform of the density correlation function develops a $`1/z`$ pole. This is a widely studied  model for the dynamics of supercooled liquids. However, with the presence of $`\gamma _q`$ at high density when $`\mathrm{\Gamma }`$ gets large, the pole shifts to $`1/(z+\gamma _q)`$. It has been demonstrated  that in the small $`q`$ and $`\omega `$ limit, the quantity $`\gamma _q`$ $``$ $`q^2\gamma `$, where $`\gamma `$ remains finite, instead of self consistently going to zero. This implies a diffusive decay of the density correlation restoring ergodicity in the longtime limit. The quantity $`\gamma `$ that provides a mechanism that cuts off the sharp transition of the fluid to an ideal glassy phase is $`O(k_BT)`$  to leading order in the perturbation theory. It is a consequence of the coupling of the density and current correlation in the compressible fluid. Formal expression for $`\gamma `$ was obtained in self consistent form in ref.  in terms of current correlation function. However beyond the one loop order explicit evaluation of the strength of the cutoff function becomes extremely complicated. In the present work, we estimate this by interpreting it as similar to slow diffusion of free volumes in the dense liquid. We use the simple model of defect correlation as $`\varphi _q(t)=e^{\delta (q)t}`$ and replace the $`\gamma (q)`$ in equation (2) by $`\delta (q)`$ to have a closed equation for the density correlation function. The density correlation function is obtained by solving equation (2) in the time space . The quantity $`\delta `$ is treated as a parameter and is adjusted to obtain agreement in the asymptotic behavior of the dynamic structure factor as seen in experiments. For this purpose the data from light scattering on colloids by Pusey et al. is  is considered. The solution of the MCT equations provide the result for the relaxation of the density correlation function with time scaled w.r.t the microscopic time scale for the Hard Sphere system - namely the Enskog time. The present analogy of the Hard Sphere system with the colloids is used only with regard to the nature of the cooperativity. The unit of time $`\tau `$ in terms of which the experimental results are presented, then needs to be related to the characteristics microscopic time scale of the Hard Sphere model i.e. $`\sqrt{\beta m}\sigma `$. $`m`$ is the mass of the liquid particle and $`\beta `$ is the Boltzmann factor. The latter is related to the Enskog time  for the system. The MCT models for the dynamics we consider here does not describe the short time dynamics or the bare transport coefficient for Colloids - the relevant comparison is over the long time scales or the nature of the correlated motion of individual units. These correlated motions are quantified in terms of the mode coupling terms. Thus the comparison of the theoretical data with the results of the experiment is done upto an undetermined scaling factor of the short time scales in the two systems. The decay of density correlation as a function of time is evaluated over a suitable wave number grid. The generalized transport coefficients expressed in terms of the coupling of density correlations from the mode coupling contribution is then evaluated. The longitudinal viscosity is computed relative to its bare or short time contribution. The dependence of the Generalized Viscosity can be probed in the present model from various approach - namely the time dependence as well as the wave vector dependence of the long time limit. The detailed results from the numerical solutions are presented in the next section.

### 2.2 Long time Diffusion

The density fluctuation in a dense fluid at finite wave number close to the peak of the static structure factor follows a diffusive behavior. Once again similar to the viscosity the diffusion constant has short time value $`D_s`$ and the long time $`D_L`$ behaviors. If all kinds of correlated motions are completely ignored then the diffusion coefficient is related to the bare transport coefficient. At high densities the cooperative nature of the dynamics is essential in computing the diffusion coefficient. This is obtained by including in the memory function the mode coupling contribution. The short time dynamics depends on the microscopic nature of the system. We use this diffusion constant to define the unit of time namely $`\tau =\sigma ^2/D_s`$ The mode coupling part is computed through the self-consistent expression in terms of the density correlation function. For the DACF, the model described in the first part of this section is used. Once again instead of focussing on the microscopic model we compute the ratio of the longtime to short time diffusion coefficient. This ratio of diffusion coefficients for the density fluctuations at finite wave number is expressed in terms of the memory function through the formula

$$\frac{D_L(q)}{D_s}=\frac{1}{s(q)\stackrel{~}{\mathrm{\Gamma }}_q}$$
(7)
where $`\stackrel{~}{\mathrm{\Gamma }}_q=_0^{\mathrm{}}𝑑\tau \mathrm{\Gamma }_q^{mc}(\tau )`$, is the long time limit of the memory function given by equation (3) and s(q) is the static structure factor. The colloidal time ($`\tau `$) is related to the Enskog time $`t_E`$ through a constant $`\mathrm{\Delta }`$, i.e. $`\tau =\mathrm{\Delta }t_E`$. The numerical results for the transport coefficients obtained from these formulations w.r.t to density and wave number variations are described in the next section.

## 3 Results for the transport coefficients

The dynamic density-density correlation function $`\psi _q(t)`$ is obtained as a solution of the nonlinear equation obtained by inverse transforming equation (2). The generalized transport coefficient $`\mathrm{\Gamma }_q(t)`$ has two parts namely the bare and the mode coupling part. The bare quantities in the Hard Sphere system for which the dynamical equations are solved are obtained in terms of standard models of kinetic theory and the unit of time is expressed in terms of the Enskog time . The integro-differential equation obtained from equation (2) is solved for $`\psi _q(t)`$ self-consistently over a wave vector grid of $`N=200`$ and an upper cutoff of $`\mathrm{\Lambda }\sigma =30`$ using numerical integration. Here we have used the Percus-Yevick  solution with Verlet-Wiess  correction for the hard sphere static structure factor. The constant $`\mathrm{\Delta }`$ relating the two time scales, colloidal and the microscopic time, is estimated from the dynamic structure factor obtained from colloidal experiments. The colloidal time ($`\tau `$) is estimated as $`2.9\times 10^5`$ times Enskog time ($`t_E`$). As discussed in the previous section the cutoff function $`\gamma `$ in equation (2) is estimated through a parameter $`\delta `$ which is inversely related to the relaxation time for the defects. The parameter $`\delta `$ is adjusted to get an agreement with the light scattering data, reported in Fig. 5(a) of Ref.  for the dynamic structure factor, in the asymptotic limit. The light scattering data supplied in Ref.  was thus used to obtain the corresponding value of the parameter $`\delta `$ at that density. In Fig. 1 the $`\delta `$ values thus obtained in units of $`\sigma ^2/\tau _H`$ are plotted as a function of packing fraction for the dense fluid. Here $`\tau _H=10t_E`$ is the hard sphere unit of time with respect to which the time is scaled. For intermediate densities, the $`\delta `$ value follows a power law behavior approaching divergence at $`\eta 0.536`$ as shown by dotted line. With increase of density the viscosity smoothes off to a slower enhancement. The parameter $`\delta `$ as estimated from the asymptotic behavior of the dynamic stucture factor of the colloids, is then used to compute the density correlation function from the model equations self consistently. This is also used to compute the memory function over very long time scales and the transport coefficients in the supercooled state is thereby computed in terms of the Mode Coupling contributions accounting for the correlated motions at higher densities. The long wavelength limit of the memory function $`\mathrm{\Gamma }_q(t)`$, reduces to the form,

$$\underset{q0}{lim}\mathrm{\Gamma }_q(t)=\frac{1}{\beta m}\frac{d\stackrel{}{k}}{(2\pi )^2}\{\stackrel{~}{V}^{(1)}(k)\psi _k(t)\varphi _{k_1}(t)+\stackrel{~}{V}^{(2)}(k)\psi _k(t)\psi _{k_1}(t)\}$$
(8)
with the vertex functions given by,

$$\stackrel{~}{V}^{(1)}(q0,k)\stackrel{~}{V}^{(1)}(k)=2c(k)\{2yS(k)+\kappa \}S(k)+\frac{2}{3}c_1(k)S^2(k)$$
(9)
and in a similar way,

$$\stackrel{~}{V}^{(2)}(k)=\left[c^2(k)+\frac{1}{5}c_{1}^{}{}_{}{}^{2}(k)+\frac{2}{3}c(k)c_1(k)2yc(k)\frac{2y}{3}c_1(k)\right]n_0S^2(k).$$
(10)
where $`c_1(k)=kc^{}(k)`$. In Fig. 2 the generalized memory function $`\mathrm{\Gamma }(t)`$ in the long wavelength limit is plotted as a function of time for the packing fraction $`\eta =0.542`$. Here the time $`(\tau _H)`$ is expressed in the units of the Enskog time, $`t_E`$. The corresponding value of $`\delta `$ is $`1.0\times 10^7`$ in units of $`(\sigma ^2/\tau _H)`$ and the metastability parameters in the MCT equation are $`\kappa =0.01,y=0.12`$ respectively. The metastability parameter $`y`$ determines the potential well used for the defect density. In figure 3 the form of the well used in the calculations for figure 2 is shown. We have expressed the result for $`\mathrm{\Gamma }(t)`$ in a dimensionless quantity in terms of the unit $`(\beta m)^1`$. $`\mathrm{\Gamma }(t)`$ also follows the two step relaxation process. At intermediate time scales long compared to the microscopic time scales, it follows a power law decay finally crossing over to the stretched exponential decay at longer times. The variation of the stretching exponent $`\beta `$ as a function of packing fraction is shown in Fig. 4. With increase in density $`\beta `$ decreases implying increase or stretching of relaxation of the DACF at higher densities. This critical slowing down of the system is estimated through the zero frequency limit of the longitudinal viscosity. In Fig. 5 longitudinal viscosity scaled with respect to the corresponding bare transport coefficient, $`\mathrm{\Gamma }(z=0)`$ is plotted as a function of packing fraction. Over the low density region the viscosity shows a power law behavior approaching a divergence around $`\eta =0.538`$. The corresponding exponent is $`1.1`$. With increase of density an abrupt increase in the viscosity by approximately $`4`$ orders of magnitude is observed. Thus at the onset of the glassy regime, a relatively strong enhancement in viscosity is obtained from the self-consistent mode coupling models. We also present the wave number dependence of the viscosity to illustrate the nature of transport over different length scales. Indeed an analysis of the self consistent expression for the density correlation function indicates that the relevant quantity for long time decay is inverse of $`S(q)\mathrm{\Gamma }_q`$. This is plotted in Fig. 6 scaled with respect to the bare time $`v_o\sigma `$ , where $`v_o`$ is the thermal velocity of the particles. Here it is shown as a function of the wave number $`q\sigma `$ for $`\eta =0.542`$. The values for the metastability parameters $`\kappa `$ and $`y`$ are same as stated for figure 2. The result shows that the time scale indicating decay rates for different wave numbers is slowest near the peak of the structure factor and for small wave number it approaches a constant value. It should be noted that this is the final decay process in the case when one completely ignores the cutoff mechanism. Thus with increase of density this will go into a delta function type peak. However the existence of the cutoff mechanism ensures the proper ergodic behavior through the diffusive process over longest time scales  At very small wavelength i.e. large wave numbers the transport coefficients are same as calculated from the short time properties of the liquid, i.e., $`\mathrm{\Gamma }_q1`$. This indicate that over small length scales the effect of correlated motion is insignificant even at very high density and the mode coupling has negligible implications on the transport properties. With increasing wavelength the viscosity increases and is maximum at $`\lambda \mathrm{}`$ or $`q\sigma 0`$. The signature of the structure of the liquid are seen in terms of very small maximas and minimas (see Fig. 6) over the whole wave vector range. Fig 7 shows the ratio of long time to short time diffusion coefficient at the peak of the structure factor as a function of density of the supercooled system. It shows a power law decay of diffusion at intermediate densities which shows a divergence around $`\eta =0.538`$ with an exponent close to unity. As the density increases further the decrease in the diffusion coefficient slows down due to the presence of cutoff effect in the present model. The diffusion coeffients are calculated in colloidal units with $`\mathrm{\Delta }=2.9\times 10^5`$. To study the dynamics at the microscopic time scales we have calculated the ratio of the long time to short time diffusion coefficient as a function of wave vector. Fig 8, presents the behavior of the Diffusion Coefficient over a wave vector range around the diffraction peak. The results shown correspond to the packing fraction $`\eta `$=.542. We observe that the diffusion coeffient is minimum at the peak of the structure factor where maximum number of nearest neighbours are present. We observe qualitatively similar behavior as observed by Pusey et al. . The values of the metastability parameter used are same as mentioned above for figure 2. It should be noted that it is only over the finite wave vector range one can describe the density fluctuations following a diffusive mechanism even in the models where there is no cutoff mechanism of hydrodynamic diffusion is considered. Of course in this case it will only be vaild upto densities below the ideal transition point beyond which transport coefficient at all wave numbers would diverge in the ideal nonergodic state. However, with cutoff mechanism taken into account in the present work, as a motion of defect density - one is not restricted to densities below the ideal glass instability only. Fig. 8 considers the effect on the diffusive process at finite length scale within the self consistent formulation of the problem where the cutoff mechanism is also being taken into account.

## 4 Discussion

The standard mode coupling theory of glass transition predicts a sharp dynamic transition with diverging transport coefficients. However in realistic system such a freezing does not occur through substantial slowing down of the ralaxation as seen in the light scattering data of colloids.

In the present work we considered a phenomenological model for the dynamics of the supercooled system where structure factor for the hard sphere system was used as an input for thermodynamic properties. The density correlation function is the key quantity in the mode coupling model. In this formulation the self-consistent form of the cutoff function used is obtained in terms of the coupling to current correlation functions to lowest order in perturbation theory. The present model includes the phenomenological description of slowly decaying defects in the amorphous structures. We interpret the cutoff mechanism of final decay of density correlation function as being related to the diffusion of free volume or defect in the amorphous system. The strength of the cutoff function is estimated by seeking an agreement with the relaxation of dynamic structure factor seen in experimental studies on colloidal systems. The correlation of the density fluctuations in the presence of very slowly decaying defect correlations contribute to the transport coefficients calculated using the feed back mechanism. The present theory of course does not compute the bare transport coefficient which is related to the short time scale dynamical behavior. The mode coupling contributios are computed relative to the short time or the bare properties. The longitudinal viscosity in the long time limit is computed and shows an enhancement by four orders of magnitude over small density range with an approach to power law divergence around a critical range of density $`\rho =1.028`$ which finally smooths off to a slower decay. The wave vector dependence of the transport coefficient is also considered. The decay rate over finite wave vector shows a minimum at the peak of the structure factor.

The relaxation of density fluctuations over the intermediate length scales follows a diffusive behavior . Due to effects of correlated motions in the fluid particles at high density the diffusion coefficient is changed over long times. We compute this by taking into account the Mode Coupling effects due to slow decay of density correlations. Thus the ratio of the long time to short time diffusion constant is calculated showing that the diffusion process slows down. However the diffusion process is fater than the final diffusive decay of defect density fluctuations given by $`\delta `$ in the $`q0`$ limit. The latter is related to the ergodicity restoring mechanism in the fluid and correspond to the largest time scale in the dynamics.

The present work demonstrates that the very slow decay of DACF seen in colloids can be understood in terms of the extended MCT and the corresponding values of the transport coefficients like longitudinal viscosity or long time diffusion constant is calculated in terms of the short time properties of the system.

## Acknowledgements

SPD acknowledges support from the NSF grant INT9615212. SS acknowledges the financial support from the UGC, India.

## Figure Captions

### Figure 1

Plot of $`\delta `$ in units $`\sigma ^2/\tau _H`$ Vs packing fraction.

### Figure 2

The normalized generalized viscosity $`\mathrm{\Gamma }(t)`$ in units of $`(\beta m)^1`$ as a function of time at $`q\sigma =0.0`$, $`\eta =0.542`$, $`y=0.15`$, $`\kappa =0.01`$ and $`\delta =1.0\times 10^7`$.

### Figure 3

Variation of depth of the potential with $`n/n^{}`$ for $`\eta _c`$ 0.542. In the figure $`h^{}(n)`$ represent the dimensionless quantity $`h(n)\beta ϵn\sigma ^3`$. The metastability parameters used are $`\kappa =0.01`$ and $`y=0.12`$.

### Figure 4

Variation of $`\alpha `$-relaxation stretching exponent $`\beta `$ with packing fraction.

### Figure 5

Plot of the longitudinal viscosity $`\mathrm{\Gamma }`$ scaled w.r.t to the bare transport coefficient as a function of packing fraction $`\eta `$.

### Figure 6

variation of $`\mathrm{\Gamma }_q^1`$ in units of $`v_0\sigma `$ with wave number $`q\sigma `$ at $`\eta =0.542`$.

### Figure 7

Variation the longtime to short time diffusion constant $`(D_L/D_s)`$ as a function of packing fraction.

### Figure 8

Wave vector dependence of the ratio of longtime to short time diffusion constant,
$`D_L(q^{})/D_s\times 10^5`$ near the peak of the structure factor.


format: text

text:
# References

On the Spin-Orbital Structure of the “Zero” Magnetization
Spin-Precessing Modes of the Superfluid $`{}_{}{}^{3}HeB`$

G. Kharadze and N. Suramlishvili

Andronikashvili Institute of Physics, Georgian Academy of Sciences, 6 Tamarashvili St., 380077, Tbilisi, Georgia

The spin-orbital configurations of the coherently precessing spin modes characterized by a small value of the magnetization $`(MM_0=\chi H_0)`$ is considered. Various regimes are analysed depending on the transverse rf-field strength.

PACS: 67.57.Lm

In the superfluid phases of liquid $`{}_{}{}^{3}He`$ the unusual spin-precessing modes can be exited at the magnitude $`M=|\stackrel{}{M}|`$ of the magnetization essentially different from its equilibrium value $`M_0=\chi H_0`$ (here $`\chi `$ is the magnetic susceptibility and $`H_0`$ denotes the strength of an applied static magnetic field). Such a possibility is realized, in particular, for $`{}_{}{}^{3}HeB`$ at $`M=M_0/2`$ (or $`M=2M_0`$), as has been pointed out in Ref. and demonstrated experimentally \[2-4\] (for the details see, also, Ref.). These homogeneously precessing spin-modes are stabilized at the local minima of time-averaged dipole-dipole potential $`U_D`$ and are characterized by the specific (non-Leggett) orbital configurations of the order parameter.

For the superfluid $`B`$-phase

$$U_D=\frac{2}{15}\chi _B\left(\frac{\mathrm{\Omega }_B}{g}\right)^2\left(Tr\widehat{R}\frac{1}{2}\right)^2,$$
(1)
where $`\mathrm{\Omega }_B`$ is the longitudinal NMR frequency, $`g`$-the gyromagnetic ratio for $`{}_{}{}^{3}He`$ nuclei and the orthogonal matrix $`\widehat{R}`$ describes 3D relative rotations of the spin and orbital degrees of freedom. By introducing the triples of Euler angles $`(\alpha _S,\beta _S,\gamma _S)`$ and $`(\alpha _L,\beta _L,\gamma _L)`$, which parametrize rotations in the spin and orbital spaces, respectively, it can be shown that $`U_D=U_D(s_Z,l_Z,\alpha ,\gamma )`$, where $`s_Z=\mathrm{cos}\beta _\mathrm{S}(\mathrm{l}_\mathrm{Z}=\mathrm{cos}\beta _\mathrm{L})`$ is the projection of the Cooper pair spin quantization axis (orbital momentum quantization axis) along the direction of the applied static magnetic field $`\stackrel{}{H}_0,\alpha =\alpha _S\alpha _L`$ and $`\gamma =\gamma _S\gamma _L`$. Measuring the energy density in units of $`\chi H_0^2`$ and noticing that $`U_D`$ is a periodic function of $`\alpha `$ and $`\gamma `$, we have:

$$\frac{U_D}{\chi _BH_0^2}=\epsilon f(s_Z,l_Z,\alpha ,\gamma )=\epsilon \underset{kl}{}f_{kl}(s_Z,l_Z)e^{i(k\alpha +l\gamma )}$$
(2)
with $`\epsilon (\mathrm{\Omega }_B/\omega _0)^2`$ (here $`\omega _0=gH_0`$).

In the case of a strong magnetic field ($`\omega _0\mathrm{\Omega }_B`$) the angular variables $`\alpha `$ and $`\gamma `$ perform rapid rotations within a long time scale $`1/\mathrm{\Omega }_D`$ and the result of the time-averaging procedure of Eq.(2) essentially depends on the possible presence of a slow combination $`\varphi _{kl}=k\alpha +l\gamma `$. Since in the strong-field case ($`\epsilon 1`$) $`\dot{\alpha }\omega _0`$ and $`\dot{\gamma }gM/\chi `$, the phase $`\varphi _{kl}`$ turns out to be a slow variable at the resonance condition $`M/M_0=k/l`$. Addressing the explicit expressions of the Fourier coefficients $`f_{kl}`$ for $`{}_{}{}^{3}HeB`$ (see Ref.) it is concluded that, along with the conventional resonance at $`M=M_0`$ with $`k=l=\pm 1,\pm 2`$, the two other (unconventional) resonances are realized at $`M=M_0/2(M=2M_0)`$ with $`2k=l(k=2l)`$. This happens in the case of $`l_Z1`$ (non-Leggett orbital configuration) where $`f_{12}=f_{21}0`$. In particular, at $`M=M_0/2`$ the time-averaged dipole-dipole potential

$$\overline{f}=1+2s_Z^2l_Z^2+(1s_Z^2)(1l_Z^2)+\frac{2}{3}\sqrt{1s_Z^2}\sqrt{1l_Z^2}(1+s_Z)(1+l_Z)\mathrm{cos}\varphi _{12}$$
(3)
and the half-magnetization (HM) spin-precessing mode is trapped at the local minima of Eq.(3) .

In Refs., along with the above-mentioned HM spin-precessing mode, another unusual spin-precessing state with $`MM_0`$ has been observed experimentally. This “zero” magnetization mode is not within the category of the resonances with a slow phase. Instead, the “zero” magnetization spin-precessing mode can be stabilized at the balance of the dissipative energy losses and the transverse rf-field energy pumping, as has been discussed in Ref.. The experimentally realized “zero” magnetization spin-precessing mode seems to be characterized by the Cooper pairs orbital configuration close to $`l_Z=0`$, although in Ref. the case of the Leggett configuration ($`l_Z=1`$) was adopted. The computer simulation was used to resolve this controversy (see, e.g., ). Here we apply an analytical approach to the same question.

In order to construct the stationary spin-precessing state with $`MM_0`$ we address the equations describing the evolution of the spin density $`\stackrel{}{S}=\stackrel{}{M}/g`$. In the strong magnetic field case the spin dynamics is governed by the two pairs of canonically conjugate variables ($`S_Z,\alpha `$) and ($`S,\gamma `$) subject to a set of the equations

$$\dot{\alpha }=1+\epsilon \frac{f_D}{S_Z}+\frac{}{S_Z}\sqrt{S^2S_Z^2}h_{}\mathrm{cos}\theta \epsilon \kappa \frac{\mathrm{S}^2}{\sqrt{\mathrm{S}^2\mathrm{S}_\mathrm{Z}^2}}\left(\frac{\mathrm{f}_\mathrm{D}}{\alpha }\mathrm{s}_\mathrm{Z}\frac{\mathrm{f}_\mathrm{D}}{\gamma }\right),$$
(4)
$$\dot{\gamma }=S+\epsilon \frac{f_D}{S}\frac{S}{\sqrt{S^2S_Z^2}}h_{}\mathrm{cos}\theta \epsilon \kappa \frac{\mathrm{S}^2}{\sqrt{\mathrm{S}^2\mathrm{S}_\mathrm{Z}^2}}\left(\frac{\mathrm{f}_\mathrm{D}}{\gamma }\mathrm{s}_\mathrm{Z}\frac{\mathrm{f}_\mathrm{D}}{\alpha }\right),$$
(5)
$$\dot{S_Z}=\epsilon \frac{f_D}{\alpha }\sqrt{S^2S_Z^2}h_{}\mathrm{sin}\theta +\epsilon ^2\kappa (\mathrm{S}^2\mathrm{S}_\mathrm{Z}^2)\left(\frac{\mathrm{f}_\mathrm{D}}{\mathrm{S}_\mathrm{Z}}\right)^2,$$
(6)
$$\dot{S}=\epsilon \frac{f_D}{\gamma }.$$
(7)
This equations are put into a dimensionless form with the time measured is the units of $`1/\omega _0`$ and the spin density - in the units of $`S_0=M_0/g`$. The presence of the transverse rf field $`\stackrel{}{H}_{}(t)=H_0(\widehat{x}\mathrm{cos}\varphi (\mathrm{t})+\widehat{\mathrm{y}}\mathrm{sin}\varphi (\mathrm{t}))`$ is taken into account $`(h_{}=H_0/H_0,\theta =\alpha \varphi )`$. The dissipation in the spin dynamics is characterized by a phenomenological parameter $`\kappa `$ (as in Ref.).

Now we focus on the case of the spin precession in the regime with $`S1`$. In this situation $`\dot{\gamma }1`$ but if $`\sqrt{\epsilon }S1`$ the angular variable will change faster then $`S_Z`$ and $`S`$. Then, in order to describe the evolution of the spin system in $`\epsilon `$-approximation (the Van der Pol picture), we have to use the time-averaged dipole-dipole potential $`\overline{f}_D`$ with respect to both angular variables $`\alpha `$ and $`\gamma `$ independently. In this non-resonance case

$$\overline{f}=f_{00}=1+2s_Z^2l_Z^2+(1s_Z^2)(1l_Z^2),$$
(8)
and the disipationless dynamics of $`\alpha `$ and $`\gamma `$ is governed by the equations

$$\dot{\alpha }=1+\epsilon \frac{f_{00}}{S_Z}+\frac{S_Z}{\sqrt{S^2S_Z^2}}h_{}\mathrm{cos}\theta ,$$
(9)
$$\dot{\gamma }=S+\epsilon \frac{f_{00}}{S}\frac{S}{\sqrt{S^2S_Z^2}}h_{}\mathrm{cos}\theta ,$$
(10)
with $`S_Z`$ and $`S`$ being constants.

The stationary solutions of Eqs.(9) and (10) are found according to the conditions

$$\frac{f}{S_Z}=\frac{f}{S}=0,$$
(11)
where the free energy

$$f=\epsilon f+\frac{1}{2}(S\omega _\gamma )^2+(\omega _\alpha 1)S_Z\sqrt{S^2S_Z^2}h_{}\mathrm{cos}\theta $$
(12)
with $`\omega _\alpha =\dot{\alpha }|_{st}`$ and $`\omega _\gamma =\dot{\gamma }|_{st}`$.

In what follows we consider the case where the first (dipole-dipole) term in Eq.(12) dominates over the rest part of the free energy. Than it is concluded that the stationarity conditions (11) are reduced to the equation

$$\frac{f_{00}}{\beta _s}=2\sqrt{1s_Z^2}s_Z(13l_Z^2)=0.$$
(13)
In the situation where $`l_Z`$ is free to adjust to the minimal value of $`f_{00}`$, Eq.(13) is to be supplemented by an analogous condition

$$\frac{f_{00}}{\beta _L}=2\sqrt{1l_Z^2}l_Z(13s_Z^2)=0.$$
(14)
From Eqs.(13) and (14) it is readily concluded that the minimal value of $`f_{00}(=1)`$ is realized at the following spin-orbital configurations:

$`a)s_Z=\pm 1,l_Z=0,`$ (15)
$`b)s_Z=0,l_Z=\pm 1.`$ (16)
The case with $`s_Z=l_Z=1`$ corresponds to the maximal value of $`f_{00}`$.

The experimentally observed “zero” magnetization spin-precessing mode is developed near the spin-orbital configuration $`s_Z=1`$ and $`l_Z=0`$  and we have to verify whether it corresponds to the dynamical regime with $`\sqrt{\epsilon }<S1`$. For this purpose the set of equations for $`S`$ and $`S_Z`$ with dissipative terms is to be addressed. Using the results of Ref. and puting $`s_Z1\frac{1}{2}\beta _S^2`$, for the case $`l_Z=0`$ the following set of equations is obtained:

$`\dot{S}=4\kappa \epsilon ^2\left(1{\displaystyle \frac{\beta _S^2}{2S}}\right),`$ (17)
$`\dot{\beta _S}={\displaystyle \frac{\kappa \epsilon ^2}{S\beta _S}}\left({\displaystyle \frac{77}{9}}+{\displaystyle \frac{\beta _S}{S}}\right)+h_{}\mathrm{sin}\theta ,`$ (18)
with the stationary solution

$$S=a\left(\frac{\epsilon ^2}{h_{}}\right)^{2/3}.$$
(19)
where $`a=\frac{1}{2}\left(\frac{190\kappa }{9sin\theta }\right)^{2/3}`$. It can be verified that this solution is locally stable. Since in Eq.(19) the coefficient $`a`$ is of the order of unity, the initially imposed condition on the value of $`S(\sqrt{\epsilon }S1)`$ will be fulfilled for $`h_{}`$ being within the limits

$$\epsilon ^2h_{}\epsilon ^{5/4}.$$
(20)
For the case with $`\epsilon 10^3`$, wich corresponds to the situation realized experimentally, the value of $`h_{}=H_{}/H_010^5`$ fits Eq.(20). This consideration shows that the “zero” magnetization spin-precessing mode is realized at the spin-orbital configuration $`s_Z1,l_Z0`$, in accordance with experimental observations.

Now we turn to the quastion of stability of the “zero” magnetization spin-precessing mode near the orbital state $`l_Z=1`$, which realizes the maximum value of the time-averaged dipole-dipole potenyial (8). Again addresing the set of Eqs. for $`S`$ and $`s_Z`$ (see Ref.) this time for the case of a spin-orbital configuration with $`s_Z1`$ and $`l_Z1`$, we obtain that for $`\sqrt{\epsilon }S1`$

$$\dot{S}=\frac{160}{9}\kappa \epsilon ^2,$$
(21)
$$\dot{\beta _S}=\frac{\kappa \epsilon ^2}{S\beta _S}\left(\frac{77}{9}+\frac{\beta _S}{S}\right)+h_{}\mathrm{sin}\theta .$$
(22)
It can be easily seen that this set of Eqs. has no stationary solutions.

It is interesting to analyze another possible regime of “zero” magnetization spin-precessing state with $`S\sqrt{\epsilon }`$. In this case $`\dot{\gamma }\sqrt{\epsilon }`$ and the rate of the time evolution of angular variable $`\gamma `$ is still faster then the temporal variations of $`S`$ and $`S_Z`$. On the other hand, $`\dot{s_Z}=(\dot{S_Z}s_Z\dot{S})/S\sqrt{\epsilon }`$ and the only fast variable upon which depends the dipole-dipole potential is $`\alpha `$. In the considered situation

$$f_D=\underset{k}{}f_k(s_Z,l_Z,\gamma )e^{ik\alpha },$$
(23)
where the time-averaged value

$$f_0=1+2s_Z^2l_Z^2+4s_Zl_Z\sqrt{1s_Z^2}\sqrt{1l_Z^2}\mathrm{cos}\gamma +(1\mathrm{s}_\mathrm{Z}^2)(1\mathrm{l}_\mathrm{Z}^2)(1+\mathrm{cos2}\gamma ).$$
(24)
After having minimized $`f_0`$ with respect to the slow variable $`\gamma `$ it can be shown that $`f_0(s_Z,l_Z)`$ is highly degenerate with respect to the spin and orbital variables: the minimum of $`f_0(s_Z,l_Z)`$ is realized within the circle $`s_Z^2+l_Z^21`$. Considering again the orbital state $`l_Z=0`$ with the stationary value $`\gamma _{st}=\pi /2`$ it can be shown that the time evolution of $`s_Z`$ is governed by an equation

$$\dot{s}_Z=\frac{4\kappa \epsilon ^2}{S(1s_Z^2)}(s_Z^4s_Z^3+s_Z^2+4)h_{}\sqrt{1s_Z^2}\mathrm{sin}\theta .$$
(25)
From Eq.(25) it follow that near $`s_Z=\pm 1`$ the stationary value of $`\beta _S`$ is given by

$$\beta _S2\left(\frac{2\kappa \epsilon ^2}{3S}\frac{1}{h_{}\mathrm{sin}\theta }\right)^{1/3}.$$
(26)
It is easy to show that the solution (26) for the case $`s_Z1`$ is locally stable (in contrast to the case with $`s_Z1`$). This spin-orbital configuration $`(s_Z1,l_Z0)`$ is one of the stable spin-precessing states at $`S\sqrt{\epsilon }`$. The phase diagram of the “zero” magnetization in the degeneracy domain $`s_Z^2+l_Z^21`$ will be presented in a separate publication.

We are indepted to E.Sonin for his valuable comments.


text:
# The Hamiltonian of the V15 Spin System from first-principles Density-Functional Calculations

## Abstract

We report first-principles all-electron density-functional based studies of the electronic structure, magnetic ordering and anisotropy for the V15 molecular magnet. From these calculations, we determine a Heisenberg Hamiltonian with five antiferromagnetic and one ferromagnetic exchange couplings. We perform direct diagonalization to determine the temperature dependence of the susceptibility. This Hamiltonian reproduces the experimentally observed spin $`S`$=1/2 ground state and low-lying $`S`$=3/2 excited state. A small anisotropy term is necessary to account for the temperature independent part of the magnetization curve.

With the continued interest in the fabrication and optimization of miniaturized magnetic devices future design considerations will require an understanding of nanoscale magnetic systems. In order to transition such materials into simple devices it is necessary to be able to explain how interactions such as spin-spatial coupling and spin-spin exchange effects may couple collectively to create a seemingly single-spin system. Further, it is necessary to determine the temperature range at which such systems will indeed behave collectively. In general the properties of a nanoscale system of coupled spins depend directly on the strength of the exchange-parameters and on the size and sign of the anisotropy energy due to spin-orbit coupling. While these parameters are generally determined by the transition metal atoms, the ligands and other nonmagnetic atoms are responsible for stabilizing the array of spins. Requisite to a complete computational understanding of such spin systems is the ability to account for the strong ligand-metal interactions and to determine whether the behavior of a given spin system is mainly mediated by the anisotropy, by spin-spin coupling or by a combination of the two.

Recently, the Mn12-Acetate and Fe8 molecules have attracted considerable interest because they behave as high-single-spin systems (total spin $`S`$=10) at temperature ranges on the order of 20-60K. Due to their large magneto-molecular anisotropy energy these systems retain their moment orientation at reasonably high temperatures and exhibit the phenomena of resonant tunneling of magnetization at well defined magnetic fields.

The $`\mathrm{K}_6[\mathrm{V}_{15}\mathrm{As}_6\mathrm{O}_{42}(\mathrm{H}_2\mathrm{O})]8\mathrm{H}_2\mathrm{O}`$ molecular crystal, first synthesized by Müller and Döring, represents a transition-metal spin system in the same size regime as the Mn12 and Fe8 molecular crystals. In contrast to Mn12 and Fe8 molecules, the V15 molecule is thought to behave as a weakly anisotropic magnet composed of 15 spin $`s`$=1/2 particles which couple together to form a molecule with a total spin $`S`$=1/2 ground state. Besides the fundamental interest in understanding quantum effects in these nanomagnets they might be also relevant for implementations of quantum computers. Calculations on such correlated systems present a challenge to mean-field frameworks such as density-functional theory because it is often not possible to construct a single collinear reference state which preserves the inherent symmetry of the system and has the correct spin quantum numbers.

This work utilizes an efficient coupled multilevel analysis which relies on fitting density-functional energies to mean-field Heisenberg or Ising energies in order to determine the exchange parameters. The approximate exchange parameters gleaned from the first $`N`$ Ising configurations were used to find the next lowest energy Ising configuration and subsequently to improve the parameterization of the exchange parameters. “Self Consistency” of this approach is determined when the predicted Ising levels are unchanged by the addition of data from new Ising configurations. The coupling of the density-functional method to a classical Ising representation allowed us to determine the exchange parameters by considering only several spin configurations. Further, the resulting ground-state spin configuration within density-functional theory exhibits the correct spin projection of $`1/2`$. With the exchange parameters determined, we diagonalize the complete many-body Heisenberg Hamiltonian to calculate the susceptibility and spin correlation functions for comparison with experiment. The many-body basis is complete, so all states are allowed including non-collinear spin arrangements and quantum disordered phases .

Starting from X-ray data we generated several unit cells and isolated a single $`\mathrm{K}_6[\mathrm{V}_{15}\mathrm{As}_6\mathrm{O}_{42}(\mathrm{H}_2\mathrm{O})]`$ unit. In order to optimize the geometry within the quasi-$`D_3`$ symmetry of the V15 molecule, we initially replaced the statistically oriented water molecule at the center of the molecule by a neon atom and used an $`S`$=3/2 spin configuration which does not break the crystallographic symmetry. The Ising configuration of this molecule consists of three aligned spin-$`\frac{1}{2}`$ V atoms in the triangle and equivalent upper and lower hexagons composed of a ring of antiferromagnetically (AF) coupled spin-$`\frac{1}{2}`$ V atoms.

The geometry of the molecule was then optimized within the all-electron density-functional methodology using the generalized-gradient approximation (GGA). The calculations were performed with the Naval Research Laboratory Molecular Orbital Library (NRLMOL). Calculations on 49 geometrical configurations were performed during the conjugate-gradient relaxation of the molecule. Subsequent calculations show that the geometrical, electronic and magnetic properties of this system are unaffected by the presence or type of inert moiety enclosed within the void. Using this geometry, we performed eleven additional calculations on different spin configurations (See Table I) to determine the six exchange parameters ($`J^{}`$s) of the Heisenberg Hamiltonian

$$H=J_{ij}𝐒_i𝐒_j,$$
(1)
as well as the spin configuration of the density-functional ground state. The $`J^{}`$s used in the above Hamiltonian are defined according to Fig. 1. As shown in Table I, we have included high-spin configurations (XI, XII, XIII), which generally have some symmetry as well as lower-spin non-symmetric configurations. The energy for the high-spin $`S=15/2`$ ferromagnetic (FM) state (XIII) of 873 meV is predominantly caused by a large AF exchange coupling ($`J`$) between the most closely bonded hexagonal V atoms. However, the 113 meV splitting between the $`S`$=9/2 and $`S`$=15/2 states (XII and XIII) shows that there is a reasonably strong AF coupling, approximately 18 meV on average, between the triangular and hexagonal atoms. All of the data displayed in Table I has been used to determine the exchange parameters from a least square fit to the mean-field solution of the Heisenberg Hamiltonian (1). The fit is very good with errors ranging from 0.1-1.55 meV. The fit leads to exchange parameters of $`J`$ = 290.3 meV, $`J^{}`$= -22.7 meV, $`J^{\prime \prime }`$=15.9 meV, $`J1`$ = 13.8 meV, $`J2`$ = 23.4 meV and $`J3`$ = 0.55 meV, where positive numbers correspond to AF and negative ones to FM interactions. The ferromagnetic interaction $`J^{}`$ is a surprising result  and deserves further discussion since it is qualitatively different from earlier assumptions based on entirely AF interactions. Ferromagnetic coupling is possible without polarizing the oxygens through a 4’th order process similar to super-exchange. In super-exchange, the intermediate state has the lowest $`d`$-orbital on the V doubly occupied with up and down electrons. However, electrons can also hop to higher energy $`d`$-orbitals on the V’s. In this case both parallel and antiparallel spins are allowed without violating the Pauli exclusion principle, but Hund’s rule coupling prefers parallel alignment. The super-exchange process (same $`d`$-orbital) completely excludes the process with same-spin electrons while the FM process (different $`d`$-orbitals) merely favors FM alignment. Thus a FM coupling is obtained if the V-O hopping matrix elements into the higher $`d`$-orbital are significantly larger than the matrix elements for the hopping of O electrons into the lowest energy $`d`$-orbital. The occurrence of such interactions are possible in a low-symmetry system such as the one studied here.

Even with this FM interaction, our spin Hamiltonian yields an $`S`$=1/2 ground state composed largely of Ising configurations similar to the one depicted in Fig. 1. This Ising configuration was predicted from the $`J`$’s from earlier fits to DFT energies and corresponds to the ground state DFT configuration (I).

We have fully diagonalized the Heisenberg Hamiltonian (1). Using all symmetries, the largest irreducible many-body subspace has dimension 2145. We find a spin-1/2 Kramer doublet as the ground state with a low-lying spin-3/2 quadruplet as shown in Fig. 1. The rest of the spectrum is well separated from these eight states. The large value of $`J`$ binds the spins in the hexagons into singlets. The low-energy physics arises from the inner triangle spins interacting with each other both directly and with an effective coupling through the hexagons, yielding the doublet-quadruplet spectrum . There are two important energy scales in the spectrum: $`\mathrm{\Delta }`$, the gap between the doublet and quadruplet, and $`J`$, the energy at which the singlets in the hexagon break and the molecule starts to behave as more than three spins.

The low-energy effective interaction between inner triangle spins proceeds via $`J1`$ and $`J2`$ (which frustrate each other) to a hexagon singlet, via $`J^{}`$ and $`J^{\prime \prime }`$ to a neighboring singlet, and finally via $`J1`$ and $`J2`$. A larger contribution to $`\mathrm{\Delta }`$ comes from a direct interaction through $`J3`$ mediated by hopping through O and As levels. Thus simple perturbation theory  yields

$$\mathrm{\Delta }=\frac{3}{4}\frac{(J2J1)^2(J^{\prime \prime }J^{})}{J^2}+\frac{3}{2}J3.$$
(2)
Comparing our calculated susceptibility with experiment , we find the low-temperature behavior indicates our doublet-quadruplet gap $`\mathrm{\Delta }10`$K to be significantly larger than the experimental value of $`\mathrm{\Delta }3.7`$K, while the high-temperature behavior shows our calculated value of $`J`$ is too large. Both of these discrepancies can be explained almost entirely by $`J^{}`$s that are too large within the density-functional-based treatment, as known for other vanadium systems .

Agreement with experiment for the low temperature $`\mathrm{\Delta }`$ can be achieved by dividing all $`J^{}`$s by a constant factor of 2.9. Fig. 2 shows our calculated result with rescaled $`J^{}`$s compared to experimental data from Chiorescu. A uniform scaling of our calculated exchange parameters is not able to obtain the right low and high temperature behavior at the same time. The high temperature behavior could be improved by futher reducing $`J`$ with corresponding adjustments to the other couplings to keep $`\mathrm{\Delta }`$ constant.

A set of AF interactions  also fits the experimental results. In fact any set of parameters with the correct values of $`J`$ and $`\mathrm{\Delta }`$ given by Eq. (2) will fit the effective moment well. To support future efforts aimed at experimentally confirming our relative parameters, we have calculated the spin-spin correlation function $`C_S=4S_i^qS_j^q`$, with $`q`$ representing the arbitrary quantization axis. Since $`J`$ is always largest, the spins tend to form a singlet along this bond, yielding $`C_S(J)1`$ (AF correlation). The three inner triangle spins form a non-collinear arrangement with $`C_S0.33`$ between all pairs of spins. The direct $`J3`$ coupling reduces the correlations between triangle and hexagon spins. We find $`C_S(J1)+0.02`$ and $`C_S(J2)0.02`$ (FM and AF correlations, respectively), while the coupling constants of Ref.  yield $`C_S(J1)+0.12`$ and $`C_S(J2)0.14`$. This prediction should be measurable with neutron scattering.

Chiorescu et al. observe that rotation of the spin projection is achieved without encountering a barrier and Dobrovitski et al. posit that the V15 molecule is indeed a low anisotropy system. As shown below the existence of either easy-plane or easy-axis anisotropy will shift the $`M_s`$=3/2-1/2 Zener-Landau tunneling transitions that have been observed by Chiorescu et al..

Recently, Pederson and Khanna have developed a new method for accounting for second-order anisotropy energies. This method relies on a simple albeit exact method for spin-orbit coupling and a second-order perturbative treatment of the spin-orbit operator to determine the dependence of the total energy on spin projection. Initial applications to Mn12 lead to a density-functional-based second-order anisotropy energy of 55.7K which is in essential agreement with the experimentally deduced values of 55.6K. We have generalized this methodology to systems with arbitrary symmetry and have calculated the anisotropy energy for several different spin configurations of the V15 molecule.

We have calculated the anisotropy energy for the lowest Ising configurations with one, three, nine and fifteen unpaired electrons, as given in Table I. In all cases we find that the V15 possesses weak easy-plane anisotropy. This result ensures that anisotropy effects will not change the total spin of the V15 ground state. Examination of the expression for the second order anisotropy energy in Ref. , shows that such energies do not necessarily scale as the square of the total moment. Indeed, as shown in Table I, we find that the anisotropic effects are in fact only weakly dependent on the total spin and that the energy of the Ms=S state increases by approximately $`\delta `$=0.8K to 1.8K. Chiorescu et al. show that the broadening of the Zener-Landau tunneling fields decreases with temperature for the $`|1/2,1/2`$ to $`|1/2,1/2`$ transitions but are independent of temperature for the $`|3/2,3/2`$ to the $`|3/2,1/2`$ transitions. This behavior is exactly what is expected from a sample containing weak easy-plane spin anisotropy. At sufficiently low temperatures, only the $`S=3/2`$ and $`S=1/2`$ states are relevant, and the field-dependent crossing of these states depends on whether the magnetic field is parallel or perpendicular to the easy-plane. The broadening is proportional to the difference of the magnetic anisotropy energy for the different spin configurations involved. Due to the small anisotropy in V15 the effect will be small. Although it is not possible to translate the DFT obtained anisotropy energies directly to the quantum mechanical many-spin ground state discussed here, we obtain from these energies tunnel field broadenings between 0.1 T to 0.48 T which envelope the experimentally observed field broadening of about 0.2 T. In powdered samples, the small easy-plane anisotropy would lead to a broadening in the tunneling field and in single crystals the effect would change the tunneling fields as a function of field orientations.

To summarize, we have performed accurate all-electron density-functional calculations on the V15 cluster as a function of geometry and spin configuration. By dynamically coupling the mean-field density-functional approach to exact-diagonalization of a many-spin Heisenberg representation, we have efficiently determined the lowest density-functional configurations and the entire Heisenberg spin excitation spectrum. Our calculations suggest that the small experimentally observed orientational dependence of the tunneling field for the $`M_s`$=3/2 to $`M_s`$=1/2 is a signature of configuration dependent magnetic anisotropy in this molecule. The method used here is general and allows one to characterize both systems which are potentially useful for magnetic storage (Mn12, Fe8) and systems which show quantum coherence, such as the one studied here.

This work was supported in part by ONR grant N0001498WX20709 and N0001400AF00002. Computations were performed at the DoD Major Shared Resource Centers. We thank I. Chiorescu for enlightening discussions and for supplying the data used in Fig. 2.


text:
# 1 Introduction

## 1 Introduction

The talk is dedicated to the experimental selection of the events with total disintegration of nuclei or with central collisions of ones. In paper  it was experimentally shown that the processes with total disintegration of nuclei correspond to the central collisions. Therefore it can generalize the rezults of the experiments on the total disintegration of nuclei and the central collisions. First of all it is necessary to remind that at present there are many papers in which the processes of nuclear fragmentation  and the processes of total disintegration of nuclei  are considered as a critical phenomenon. Therefore it is possible to suppose that if there are cases corresponding to the critical phenomena among the events with the total disintegration of nuclei then the points of the regime changes in the behaviour of some characteristics of secondary particles ($`a_i`$) depending on the centrality degree of collisions - $`Q`$ could be observed. Some experimental data obtained in nuclear-nuclear interactions at high energy demonstrate the existence of the points of the regime changes in the behaviour of the $`a_i=f(Q)`$ distributions. For example in fig. 1 the averege multiplicity of relativistic chaged particles depending on $`Q`$ is shown for the $`{}_{}{}^{28}Si_{14}+Em`$ reactions at the energies 3.7 and 14.6 GeV per nucleon. To determine the $`Q`$ a number of charged projectile fragments ($`Z_{f_i}`$) were used. The figures were obtained from paper . The points of the regime change are observed in these dependences. These points were used by the authors to select the events with central collisions of nuclei.

In fig. 2 are shown the averege values of pseudorapidity - $`\eta (\eta =lntg(\theta /2))`$ for s-particles(the particles with $`\beta >0.7`$) depending on the number of $`g`$-particles( the particles with $`\beta 0.7`$ ) for $`pEm`$-reactions at the momenta of $`p_0=`$4.5; 24.0;50.0;67.0 and 200.0 GeV/c. This figure was obtained from paper . The dashed line in the figure corresponds to the cascade-evaporation model calculations. The points of the regime changes in these distributions are also seen. The cascade-evaporation model calculations cannot describe numerically these distributions.

In fig. 3 the $`Q`$-dependences of the absolute values of one-particle correlation function $`R(pt,Y)`$ are shown for the protons in the $`{}_{}{}^{12}CC`$-interactions at the momentum of 4.2 A GeV/c. The data was obtained from paper . Here $`Q`$ is a number of protons in an event. The points of the regime changes in these figures are also seen. Thus the results clearly demonstrate that there exist the points of the regime change of the behavior of $`a_i=f(Q)`$ distribution which could fix the region of the central collisions of nuclei. We believe that these points could be used to select the events with the central collisions of nuclei. So to find the points of the regime changes it is possible to investigate the $`a_i=f(Q)`$ distributions’ behaviour. The $`Q`$ could be determinated as:

\- a number of fragments($`n_f`$);

\- a number of protons in an event($`n_p`$);

\- an energy flow of particles at the emission angulars $`\theta 0(E_{ZDC})`$ or $`\theta 90(E_t)`$ ones.

For further confirmation of the results on the existence the regime change points and for clarification of the reasons of these points appearance we studied the influence of nuclear fragmentation process on the behaviour of the events number depending on $`Q`$. Supposing these points beeing connected with the appearence of a critical phenomenon the fragment number change could also have a critical character with the $`Q`$ increase. Because the intermediate nuclear formations (for example nuclear cluster) could be a source of the nuclear fragmentation.

## 2 Experiment

We used 20407 events of $`{}_{}{}^{12}CC`$\- reactions at the momentum of 4.2 A GeV/c obtained from the 2-m propane bubble chamber of LHE, JINR(the methodical details are described in papers ).To determine the values of the $`Q`$ two variants were considered. In the first variant the values of the $`Q`$ were determined as

$$Q=n_++N_pn_{},$$
(1)
here $`n_+,n_{}`$ and $`N_p`$ are the numbers of identified $`\pi ^+,\pi ^{}`$-mesons and protons respectively. In that determination the $`Q`$ is a number of all the protons in an event without taking into account a remainder of nuclei. In the second variant the Q values were determined as

$$Q=N_+n_{},$$
(2)
here $`N_+`$ are charges of all the positively charged particales in an event including the nuclear fragments. In that determination the $`Q`$ is a summary charge of an event.

## 3 Results

The distributions of the events number depending on the $`Q`$ are shown in fig.4a,b. In fig.4a the empty starlets correspond to the first variant and the full starlets - to the second variant. From fig. 4a it is seen that the fragments number to determine the $`Q`$ being included the form of distributions sharply changes and has a two-steps structure.

In fig.4b are shown the $`Q`$-dependences of the events number for the calculation data obtained from the quark-gluon string model  (QGSM) without the nuclear fragments. The empty starlets correspond to the cases in which the stripping protons were not taken into account and the full starlets correspond to the cases in which the stripping protons were included. It is seen that the form of the distribution strongly differs from the experimetal one in fig.4a. There is no two-steps structure in this figure. Therefore we can assert that this difference is connected with the existence of fragments in experimental events.
Thus, the results demonstrate that the influence of nuclear fragmentation process in the behaviour of the events number depending on the $`Q`$ has a critical character. To explain this result we suppose that it could be connect with the existence of nuclear clusters.It is possible to think that with the increasing of centrality degree the probability of claster formation grows but further increasing the Q ( in the region of high Q) leads to the big clasters decay on nuclear fragments and then on free nucleons(see fig.5). It could be a reason for observation of two-step structure in the distributions. The first step connected with the formation of a cluster and second one with its decay.

## 4 Conclusion

1. The experimental results obtained in the nuclear-nuclear collisions at high energy clearly demonstrate the points of the regime change of the behaviour of some characteristics of econdary particles depending on the centrality degree of collisions. It could mean that there are cases corresponding to the critical phenomena among the events with the central collisions of nuclei.

2. For $`{}_{}{}^{12}CC`$-interaction the behaviour of the number of events, depending on $`Q`$ also depends on the number of fragments and has a two–steps form which is not reproduced by the calculated data in the framework of the QGS model without nuclear fragments. This result could be explained by the existence of nuclear clusters.

Captions of figures

1. Q-dependences of multiplicity of relativistic charged particles in Si+Em reactions at the energies of a) 3.7 GeV per nucleon; b) 14.6 GeV per nucleon.

2. $`n_g`$-dependences of the averege values of $`\eta =lntg(\theta /2))`$ for $`n_s`$-particles.

3. Q-dependences of the absolute values of one particle correlation function for the protons in $`{}_{}{}^{12}CC`$-interaction.

4. Q-dependences of the event number of $`{}_{}{}^{12}CC`$-interactions: a) experiment; b) quark-gluon string model.

5. Fig.5

Figure 4a,b.


text:
# Energy Dependence of the Contribution of Pion Exchange to Large-Rapidity-Gap Events in Deep Inelastic Scattering

We study the energy dependence of the contribution of pion exchange to large-rapidity-gap events in deep inelastic scattering. The results show that this contribution can be quite significant at low energy and that the LRG events observed by E665 collaboration in $`\mu Xe`$ and $`\mu D`$ interactions at 490 $`GeV`$ can be reasonably well described in terms of meson exchange. We also show that the distribution of the maximum rapidity for all hadrons is quite different from that for charged hadrons only and that the former exhibits also shoulder-like structure for events at 490 $`GeV`$ similar to that at HERA.

Large-rapidity-gap (LRG) events in deep inelastic scattering have been studied extensively at HERA both experimentally  and theoretically \[References,References-References\]. It has been observed that there exist a class of events in which there are no particles or energy depositions in the forward part of the detector, i.e., there is a large rapidity gap in the final state. Such events take about $`10\%`$ of the total deep inelastic events in the kinematic range $`10^4<x<10^2`$ and $`5<Q^2<120`$ $`GeV^2`$. (Here $`x`$ is the Bjorken-$`x`$, $`Q^2=q^2`$ and $`q`$ is the four-momentum transfer carried by the virtual photon.) The experimental data obtained by ZEUS  and H1  collaborations show that the distribution of the maximum pseudo-rapidity $`\eta _{max}`$ for the produced hadrons has a clear shoulder-like structure which signifies the occurrence of the LRG events.

The existence of the LRG shows that, in such events, the exchanged object between the virtual photon and the incoming proton must be colorless (c.f. Fig. 1a given later in this paper). It has been shown that different features of LRG events can be well described  in terms of pomeron ($`IP`$ ) exchange in Regge phenomenology. At the same time, it was also well known that one-meson exchange such as one-pion exchange process contributes significantly  to deep inelastic scattering $`epeX`$. One-pion exchange itself takes about $`10\%`$ of the whole $`epeX`$ events in the abovementioned kinematic region . Therefore, it was expected  that meson exchange might also contribute to the LRG events. This can be studied using Monte Carlo event generators. Such study showed that the contribution is negligible  at HERA energy. It has completely no contribution to the characteristic shoulder-like structure in $`\eta _{max}`$ distribution.

It is interesting to note that LRG events have also been observed  by E665 collaboration at FNAL in $`\mu Xe`$ and $`\mu D`$ fixed target experiments at 490 $`GeV`$. In these experiments, LRG events were defined as the events with a rapidity gap $`\mathrm{\Delta }y^{}>2`$, where $`\mathrm{\Delta }y^{}`$ is the difference between the rapidity of the target nucleon before the scattering and the lowest rapidity of the charged hadrons in the event in $`\gamma ^{}`$-nucleon c.m. frame. The distribution for the probability of $`\mathrm{\Delta }y^{}`$ greater than a given value i.e., the probability distribution $`P(\mathrm{\Delta }y^{})`$, rather than the probability density, has been given. It is interesting to note that the $`P(\mathrm{\Delta }y^{})`$ obtained by E665 does not have the same shoulder-like structure as that for $`\eta _{max}`$ distribution at HERA. Having in mind that $`\mathrm{\Delta }y^{}`$ and $`\eta _{max}`$ are essentially the same except for the small difference in reference system, we are naturally led to the following questions: Why is there such a large difference between the $`\mathrm{\Delta }y^{}`$ distribution obtained by E665 collaboration and the $`\eta _{max}`$ distribution obtained by H1 and ZEUS collaborations? How is the energy dependence of the rapidity distribution in events where $`IP`$ or $`\pi `$ is exchanged? Can meson exchange give a significant contribution to the LRG events observed by E665 collaboration?

These are the questions that we would like to study in this note. For explicity, we take one-pion exchange as an example. Similar effects should also exist for other mesons. We now start with a qualitative analysis. In fact, from the following simple qualitative analysis, we already expect that the contribution of $`IP`$ or $`\pi `$ exchange to LRG events at E665 energy may be quite different from that at HERA energy. We use the notations as those shown in Fig.1a and recall that the invariant mass of the total hadronic system and that of the hadronic subsystem $`X`$ are given by,

$$W^2(P+q)^2=Q^2(1/x1)+M^2,$$
(1)
$$M_X^2(q+q_\kappa )^2=Q^2(\xi /x1)+t,$$
(2)
respectively. Here $`tq_\kappa ^2`$ is the square of the four-momentum transfer between the proton and the virtual photon; $`\xi (qq_\kappa )/(qP)`$ can be interpreted as the fraction of the momentum carried by the exchanged object $`\kappa `$ (which represents $`IP`$ or $`\pi `$) from the incident proton in the infinite momentum frame. Independent of the reference frame one uses, the width of the rapidity distribution of the hadronic subsystem $`X`$ is proportional to $`\mathrm{ln}M_X^2`$, i.e., $`\mathrm{\Delta }y_X\mathrm{ln}M_X^2`$; and the width of the rapidity distribution of the total hadronic system is proportional to $`\mathrm{ln}W^2`$. This is illustrated in Fig.1b. We see that, for a given $`W`$, the smaller $`M_X`$ is, the larger the rapidity gap between the hadronic subsystem $`X`$ and the scattered nucleon is. Hence, $`M_XW`$ is a necessary condition for LRG to appear. From Eq. (2), we see that $`M_X`$ increases with increasing $`\xi `$ and decreases with increasing $`|t|`$ at fixed $`Q^2`$ and $`x`$. This implies that $`\xi `$ can not be large to ensure $`M_XW`$. In the experiments at HERA, the kinematic range is $`5<Q^2<120`$ $`GeV^2`$ and $`0.0001<x<0.01`$, and the typical value for $`M_X`$ is several ($`17`$, say) $`GeV`$. The corresponding $`\xi `$ value is of the order of 0.01 or less. This means that only events where $`\xi `$ is very small contribute to the LRG events at HERA. Hence, the contribution of $`\pi `$ or $`IP`$ exchange to LRG events at HERA is determined by its contribution to deep inelastic scattering at small $`\xi `$.

The $`\xi `$-dependences of the contributions of $`IP`$ and $`\pi `$ exchange to deep inelastic scattering $`epeX`$ at given $`x`$ and $`Q^2`$ are given by their contributions to the “diffractive structure function” $`F_2^{D(4)}(x,Q^2,\xi ,t)`$. To make a qualitative comparison of the two contributions, we use the approximately valid factorization theorem, i.e. take $`F_{2(\kappa )}^{D(4)}(x,Q^2,\xi ,t)=`$ $`f_{\kappa /p}(\xi ,t)F_2^\kappa (x/\xi ,Q^2)`$. Here, $`f_{\kappa /p}`$ is the corresponding flux factor; $`F_2^\kappa `$ is the structure function of $`\kappa `$ ($`IP`$ or $`\pi `$). In the case of $`IP`$ exchange, $`f_{IP/p}(\xi ,t)`$ can be parameterized as

$$f_{IP/p}(\xi ,t)=\frac{9\delta ^2}{4\pi ^2}[F(t)]^2\xi ^{12\alpha _{IP}(t)},$$
(3)
where $`\delta ^2=3.24`$ $`GeV^2`$, $`\alpha _{IP}(t)=1+ϵ+\alpha ^{}t`$, $`ϵ0.085`$, $`\alpha ^{}=0.25`$ and the elastic form factor $`F(t)`$ is given by:

$$F(t)=\frac{4m_p^22.8t}{4m_p^2t}\frac{1}{(1t/0.7)^2},$$
(4)
where $`m_p`$ is the mass of proton. The pion flux factor can be derived from the pion cloud model. For the case $`p\pi ^{}N`$, we have,

$$f_{\pi /p}(\xi ,t)=3.257\xi \frac{t}{(m_\pi ^2t)^2}\mathrm{exp}\left(\frac{m_\pi ^2t}{1.21\xi }\right),$$
(5)
where both $`m_\pi ^2`$ and $`t`$ are taken in unit of GeV2. Take,

$$F_2^{IP}(x,Q^2)=3x(1x)/2,$$
(6)
and the SMRS-P2-parameterization of $`F_2^\pi (x,Q^2)`$, we obtained their contributions to $`F_2^{D(4)}(x,Q^2,\xi ,t)`$ at $`t=|t|_{min}`$ $`m_p^2\xi ^2/(1\xi )`$ in Fig.2. From the figure, we explicitly see that $`F_{2(IP)}^{D(4)}(x,Q^2,\xi ,|t|_{min})F_{2(\pi )}^{D(4)}(x,Q^2,\xi ,|t|_{min})`$ in the small $`\xi `$ (say, $`\xi <0.05`$) region. Since the LRG events at HERA come mainly from this small $`\xi `$ region, we reach the conclusion that the contribution of $`IP`$ exchange to LRG events at HERA is dominant and that the contribution of $`\pi `$ exchange can be neglected.

The situation is, however, quite different in the E665 experiment. Here, events in the kinematic range $`1<Q^2<100`$ $`GeV^2`$ and $`0.002<x<0.3`$ were selected. Compared with those at HERA, $`Q^2`$ is much smaller and $`x`$ is larger. Thus, to obtain the same $`M_X`$ as that at HERA, $`\xi `$ should be much larger. In fact, for typical $`M_X`$ around $`17`$ $`GeV`$, $`\xi `$ is of the order of $`10^110^2`$ and can even be significantly larger than 0.1. From Fig.2, we see that $`F_{2(IP)}^{D(4)}(x,Q^2,\xi ,|t|_{min})`$ decreases with increasing $`\xi `$, but $`F_{2(\pi )}^{D(4)}(x,Q^2,\xi ,|t|_{min})`$ increases very rapidly with increasing $`\xi `$. As a result, their difference becomes very small in the region of $`\xi 0.05`$, and the latter can even be larger than the former for large $`\xi `$ ($`>0.1`$, say). Thus the contribution of pion exchange to the LRG events obtained in the E665 experiments should be quite significant compared with that from Pomeron.

Now, we explicitly calculate the rapidity distribution of final hadrons in events where $`IP`$ or $`\pi `$ exchange takes place. Presently this can only be carried out using Monte Carlo events generator. There exist a number of Monte Carlo event generators, such as POMPYT and RAPGAP, which simulate the processes shown in Fig.1a. Here, it is envisaged that the incoming proton ‘emits’ a $`IP`$ or a $`\pi ^{}`$, and the $`IP`$ or $`\pi ^{}`$ then collides with the virtual photon emitted by incident lepton to produce the sub-hadronic-system $`X`$ shown in Fig.1a. Various options for the effective $`IP`$ flux and the parton densities in $`IP`$ are available in the programs. Using such Monte Carlo program we can easily calculate the contribution of $`IP`$ or $`\pi `$ exchange to the LRG events at different energies or in different reference frames. Both POMPYT and RAPGAP are slave systems, which must be called by our own steering program. Thus we first simulate the events at HERA to check our steering program then apply it to E665 energy.

We simulate the events where the abovementioned $`IP`$ or $`\pi `$ exchange takes place using POMPYT or RAPGAP and the usual DIS events using LEPTO and obtain the $`\eta _{max}`$ distribution for each class of events respectively. We calculated them using different options for $`IP`$ structure functions. We found out that both the results from POMPYT and those from RAPGAP with different options for $`IP`$ structure function are essentially the same. Adding the different contributions together with the corresponding weights which measure the relative contribution of each type of process to the inclusive process $`epeX`$, we can obtain the $`\eta _{max}`$ distribution for final hadrons in deep inelastic scattering at HERA energy. The contribution of $`\pi `$ exchange to deep inelastic scattering $`epeX`$ in the HERA kinematic region was estimated in \[References\]. The results depend on the parameterization of pion structure function, but they are all of the order of $`10\%`$ of all DIS events at HERA energy. In Fig.3a, we show the results that we obtained by adding 12% from $`IP`$ or $`\pi `$ exchange with 88% usual DIS events from LEPTO (the lower solid and dotted lines). Here, in obtaining these results, all hadrons that can be observed by the H1 detector, i.e., those with $`3.8<\eta <3.65`$ and energy higher than 400 $`MeV`$, are taken into account. Since our purpose is to study the energy dependence of the contribution of $`IP`$ or $`\pi `$ exchange to LRG events, we do not simulate the detector effects at HERA. The results show that the usual DIS events can give a reasonable account of the shape of the $`\eta _{max}`$ distribution for values above 1.5 and that the usual DIS and $`IP`$ exchange together can well describe the $`\eta _{max}`$ distribution at HERA. The results also show that the $`\pi `$ exchange has a very small contribution to the LRG events at HERA.

Subsequently, we apply the method to the $`\mu Xe`$ and $`\mu D`$ fixed target scattering at 490 $`GeV`$ beam energy, and calculate the rapidity distribution of final hadrons in events where $`IP`$ or $`\pi `$ exchange takes place. Since our purpose is to study the energy dependence of the contribution of $`IP`$ or $`\pi `$ exchange to the LRG events, we do not take the nuclear effects into account. This means that we simply regard $`\mu D`$ or $`\mu Xe`$ scattering in the E665 kinematic region as $`\mu N`$ scattering in the same kinematic region. We select the events in the same kinematic range as that chosen by E665 collaboration  and obtain the probability distributions of $`\mathrm{\Delta }y^{}`$ in events where $`IP`$ (upper solid line) or $`\pi `$ (upper dotted line) exchange takes place compared with the E665 $`\mu D`$ or $`\mu Xe`$ data shown in Fig.4a. From the results we see clearly that, as we expected in the abovementioned qualitative analysis, $`\pi `$ exchange can have a very significant contribution  to the LRG events observed by E665 collaboration. We found out also, to reproduce the E665 data, we need a rather large contribution of $`IP`$ and/or $`\pi `$ exchange. We estimated the contribution of pion-exchange to $`\mu p\mu X`$ using the pion flux factor given in Eq.(5) and different parameterizations of pion structure functions. We found out that the results are slightly different if different parameterizations are used. But they are all of the order of $`1020\%`$ in the E665 kinematic regions. In Fig.4a, we show the results obtained by adding 20% $`IP`$ exchange with 80% from LEPTO (lower solid line) and those obtained by adding 20% $`\pi `$ exchange with 80% from LEPTO (lower dotted line). We see in particular that pion exchange contributes significantly to the LRG events but cannot account for all of them. However, what we discussed till now is only the contribution from the case $`p\pi ^{}N`$, which is an explicit and calculable example of different meson exchange processes. Similar contributions should be expected from other mesons which cannot be calculated presently because of the lack of the corresponding flux factors and the structure functions. To show what we may expect if all different meson exchange processes are taken into account, we simply add more contributions from $`\pi `$ exchange to the whole events. Hence, in Fig.4a, we show also the results obtained by adding 40% from $`\pi `$ exchange with 60% from LEPTO (dash-dotted line). We see that the results agree reasonably well with the data.

From Fig.4a, we also see that the shape of the $`\mathrm{\Delta }y^{}`$ distribution does not have a clear shoulder-like structure at the E665 energy as that in $`\eta _{max}`$ distribution at HERA (see Fig.3a). We are therefore led to the question about the reason of the disappearance of such characteristic structure for LRG events. We note that, besides the energy is lower, only the charged hadrons are taken into account by E665 whereas all the hadrons are taken into account at HERA. Thus we take also all the hadrons with energy higher than 400 $`MeV`$ into account and re-calculate $`\mathrm{\Delta }y^{}`$ distribution at E665 energy. The obtained results are shown in Fig.4b. It is interesting to see that, compared with that for charged hadrons only, the $`\mathrm{\Delta }y^{}`$ distribution obtained from LEPTO is much narrower whereas those obtained from RAPGAP for $`\pi `$ or $`IP`$ exchange remain essentially the same. Their difference becomes much larger. Adding them together with the corresponding weights mentioned above, we obtain the total $`\mathrm{\Delta }y^{}`$ distribution, which has now a significant shoulder-like structure in the case of $`IP`$ exchange. But the shoulder structure is not clear if only $`\pi `$ exchange is involved. This is similar to that at HERA. These results clearly show that it is much more efficient to distinguish events where $`IP`$ or $`\pi `$ is exchanged from the usual DIS events by studying the $`\mathrm{\Delta }y^{}`$ distribution for all the hadrons than that for charged hadrons only. To check whether this is also true at other energies, we re-calculate the $`\eta _{max}`$ distribution for the charged hadrons only at HERA energy. The obtained result is shown in Fig.3b. The results show similar effect, i.e., compared with those for all the hadrons (see Fig.3a), the $`\eta _{max}`$ distribution for the charged hadrons from usual DIS or $`\pi `$ exchange events at HERA is much wider whereas that for $`IP`$ exchange remains essentially the same. The contamination from the fluctuation in usual DIS to LRG events would be much higher if one would study charged hadrons only.

We note also that one can use different variables to describe LRG events, such as the $`\mathrm{\Delta }y^{}`$ used by E665, $`\eta _{max}`$ used at HERA or $`\eta _{max}^{}`$ and $`y_{max}^{}`$ in the hadronic c.m. system. We also studied the question of which one is more efficient by calculating the corresponding distributions in the scattering of the 490 $`GeV`$ lepton beam off the fixed target proton. The obtained results show no significant difference between these variables, all of them can give good description to the occurrence of LRG events.

In summary, using the Monte Carlo event generators POMPYT, RAPGAP and LEPTO, we showed that $`\pi `$ exchange has a significant contribution to the LRG events in the $`\mu Xe`$ and $`\mu D`$ fixed target scattering at 490 $`GeV`$ beam energy. Taking all the contributions from different meson exchange processes into account, we should obtain a reasonably well description of the corresponding E665 data at that energy. The distribution of the maximum rapidity for all hadrons is quite different from that for charged hadrons only and the former exhibits also shoulder-like structure for events at 490 $`GeV`$ beam energy similar to that at HERA.

We are in debated to C. Boros who took part in the early stage of this research. We thank Li Shi-yuan, Xie Qu-bing, and Xu Qing-hua for discussions. This work was supported in part by the National Science Foundation (NSFC) and the Education Ministry of China.

Figure captions

Fig.1: (a) Schematic of the semi-inclusive process $`epeNX`$ with the exchange of a colorless object $`\kappa ^{}`$ which can be a $`IP`$ or a $`\pi ^{}`$; and (b) Diagram illustrating the range of the rapidities of the hadrons in the subsystem $`X`$ and that for all the hadrons including the outgoing nucleon $`N`$. Here we see in particular that the necessary condition for LRG to appear is $`\mathrm{ln}W^2\mathrm{ln}M_X^2`$.

Fig.2: $`\xi `$-dependence of $`IP`$ or $`\pi `$ exchange contribution to the “diffractive structure function” $`F_{2(\kappa )}^{D(4)}(x,Q^2,\xi ,t)`$ at $`t=|t|_{min}`$.

Fig.3: Distribution of $`\eta _{max}`$ for hadrons in DIS events at HERA. The data were taken from \[References\]. The total number of the Monte Carlo events is normalized to the data. In obtaining the results in (a), all hadrons which can be observed by H1 detector (i.e. those with $`3.8<\eta <3.65`$ and energy higher than 400 MeV) are taken into account. In (b), only the charged hadrons are taken into account.

Fig.4: Probability $`P(\mathrm{\Delta }y^{})`$ of events with rapidity difference between the most backward hadron and the incoming nucleon in the $`\gamma ^{}p`$ c.m. frame to be greater than $`\mathrm{\Delta }y^{}`$ in the scattering of 490 $`GeV`$ lepton off fixed proton target. In obtaining the results in (a), only the charged hadrons with energy higher than 400 $`MeV`$ are taken into account. The data are obtained by E665 collaboration \[References\] in $`\mu Xe`$ and $`\mu D`$ scattering at 490 GeV. In (b), all the hadrons with energy higher than 400 $`MeV`$ are taken into account.


text:
# Information transport by sine-Gordon solitons in microtubules

\[

## Abstract

We study the problem of information propagation in brain microtubules. After considering the propagation of electromagnetic waves in a fluid of permanent electric dipoles, the problem reduces to the sine-Gordon wave equation in one space and one time dimensions. The problem of propagation of information is thus set.

PACS number(s): 87.17.Aa 87.17.-d 89.70.+c 89.75.Fb 89.75.-k

\]

Information is a central question in the understanding of the mechanisms regulating the brain. Questions as bounds on information, capacity of communication channels are of extreme relevance and a theory of information and communication is of utter value . Thus, how much information can be stored by a cube is an important question. Foreseeable technology making use of atomic manipulation would suggest an upper bound of around $`10^{20}`$ bits. But as technology takes advantage of unforeseen paradigms, this number could grow up. Could the bound grow without limit? With black hole thermodynamics some definite answers are forthcoming . In Quantum Field Theory the question developed unexpectedly in the framework of black hole thermodynamics and Quantum Gravity . Some knowledge was already known from Shannon’s Theory of Information as well .

Shannon imagined a system capable of storing information by virtue of it possessing many distinguishable states. A state $`a`$ is not known a priori but its probability $`p_a`$ is known. The measure of uncertainty corresponds to an entropy $`S`$ which formally coincides with the corresponding statistical interpretation of entropy:

$$S=K\underset{a}{}p_alnp_a$$
(1)
On the other hand, it has been conjectured  that the brain microtubules might be an active component of the brain functioning. It is thus natural to consider electromagnetic waves moving in that cavity as transporting and carrying information.

With these matters in mind, we consider here the effective electromagnetic wave obtained when the second quantized electromagnetic field interacts with the permanent dipole moment of the vicinal water in brain microtubules. The second quantized electromagnetic field shall be given by the usual development in frequency components . The electric dipoles can be seen as two-component spinors described as effective spin-fields. Therefore, we have, for the effective Hamiltonian,

$$H=H_{QED}\mu \underset{j=1}{\overset{N}{}}[\stackrel{}{E}_{tr}(\stackrel{}{x}_j,t)\stackrel{}{s}]+ϵ\underset{j=1}{\overset{N}{}}s_z,$$
(2)
where $`\stackrel{}{E}_{tr}`$ is the electric field transversal to the wave propagation direction and $`\stackrel{}{s}`$ is a spin field describing the electric dipole moment degree of freedom. The last term represents an effective interaction of the $`z`$-component of the electric dipole with an average electric field, that is, it represents a two-energy eigenstates system. The value of $`ϵ50`$cm$`{}_{}{}^{1}=6.3\times 10^3`$eV$`=1.0\times 10^{14}`$erg has been claimed in .

Such a problem has been considered by . We derive some results which have not been explicitely obtained in . We suppose that the electromagnetic field has a fast dependence on $`zct`$ and a slow dependence on $`z`$ and $`t`$, allowing us to write the expansion

$$\stackrel{}{E}(\stackrel{}{x},t)=\stackrel{}{E}_{tr}^n(z,t)e^{ik_n(zct)}$$
(3)
Using the equations of motion derived from the Hamiltonian (2), that is, Maxwell equations with sources, we arrive at

$$\frac{E^\pm }{z}+\frac{1}{c}\frac{E^\pm }{t}=\pm i\frac{2\pi ϵ\mu }{\mathrm{}V}s^{}$$
(4)
This is a quantum equation of motion. However we do not have any practical means to either measure the variable $`s`$, or take care of its detailed dynamics, therefore we take its quantum average. Such an average is easily obtained due to the simple description of $`s`$ in terms of Pauli matrices, leading to a result written in terms of the exponential of the field $`\theta `$, defined by

$$\theta ^\pm (z,t)=\frac{\mu }{\mathrm{}}_0^tE^\pm (z,t)_{qu}𝑑u,$$
(5)
where we take the quantum average $`_{qu}`$. Following uch a procedure in equation (4) leads to the semiclassical equation of motion

$$\frac{^2\theta ^\pm }{t\sigma }=\frac{4\pi ϵN\mu ^2}{\mathrm{}^2V}\mathrm{sin}\theta ^\pm $$
(6)
where $`N/V`$ is the number of dipoles (molecules) per unit volume, and $`\sigma =t+\frac{z}{c}`$. Above, the indices $`\pm `$ correspond to the usual combinations of the transversal direction, and we supposed also that the longitudinal direction does not propagate. This is a variant of the well known sine-Gordon equation. The one-soliton solution is given by the expression

$$E=\frac{\mathrm{}}{\mu }AsechA(t\frac{z}{\nu _0})$$
(7)
where the angular frequency characteristic of the model is

$$A=\sqrt{\frac{2\pi ϵ\mu ^2N\nu _0}{\mathrm{}^2V(c\nu _0)}}$$
(8)
and $`\nu _0`$ is the velocity of the soliton.

In order to understand the propagation of information in such a device, we follow  and consider small perturbations around the soliton, which amounts to solving the equation

$$\omega ^2\eta +i\omega \frac{\eta }{t}A_0^2\mathrm{cos}\theta _0\eta =0$$
(9)
where $`A_0=\sqrt{2\pi ϵ\frac{N}{V}}\frac{\mu }{\mathrm{}}3.1\times 10^{14}`$s-1. Plugging in back the solution $`\theta _0=`$ $`4\mathrm{arctan}\mathrm{exp}[Az/\nu _0]`$, that is, $`\mathrm{cos}\theta _0=12\mathrm{s}\mathrm{e}\mathrm{c}\mathrm{h}^2(Az/\nu _0)`$, we obtain the equation

$$i\omega \frac{\eta }{t}=2A_0^2\mathrm{sech}^2(Az/\nu _0)\eta $$
(10)
where we chose the boundary conditions such that $`\omega =A_0`$. The only solution is

$$\eta =exp[2i\sqrt{\frac{\nu _0}{c}}\mathrm{tanh}\frac{Az}{\nu _0}]$$
(11)
Let us discuss the physics behind the problem and the consequences for the constants appearing in the solution. First, the constant $`ϵ`$ is a free parameter and represents the energy of a dipole in the vicinal water. It is of the order of magnitude of difference of two molecular energy levels. The study of vibration in water indicates the value $`ϵ50`$cm$`{}_{}{}^{1}=6.3\times 10^3`$eV$`=1.0\times 10^{14}`$erg . The constant $`\mu `$ represents the permanent electric dipole moment of the water, which is the electron charge times $`0.2\times 10^8`$cm, that is, in CGS units, $`\mu 6.8\times 10^{18}`$. Finally, the number of molecules per unit volume is easily obtained for the water, it is of the order of $`0.3\times 10^{23}`$cm-3.

In order to fix the velocity of the wave, we integrate the solitonic electric field imposing that it is the unit synaptic potential coming from the quantum of transmitter producing a postsynaptic potential, typically $`0.5`$ to $`1.0`$mV, as discussed in , who has proposed it to be a quantum unit of potential in such a context. We have

$$VE𝑑z=\frac{\pi }{2}\frac{\mathrm{}\nu _0}{\mu },$$
(12)
thus obtaining, for the velocity parameter, the value $`\nu _01.4\times 10^4`$cm/s, or $`\frac{\nu _0}{c}0.5\times 10^6`$. With this result for $`\nu _0`$ we obtain for the constant $`A`$ the result $`A2.2\times 10^{11}`$s-1. Estimating the time to send information as that necessary to pass the bulk of the soliton (7), we get a rough estimate for the frequency of the waves as $`\nu \frac{A}{6}3.7\times 10^{10}`$s-1.

On the other hand, taking the average electric field in the brain as corresponding to the quantum unit of electric potential value as discussed above, namely $`1`$mV, divided by the lenght of the typical microtubule, that is, $`10`$nm we are led, for the average electric field, to the value

$$E_{ave}10^1V2\times 10^2\times 10^9m3statV/cm$$
(13)
Now, the constants $`A`$ and $`E_{ave}`$ are related by

$$A=E_{max}\frac{\mu }{\mathrm{}}=\sqrt{2}E_{ave}\frac{\mu }{\mathrm{}}4\times 10^{10}s^1$$
(14)
This is compatible with the previous value for $`A`$, giving us some confidence on the result.

The corresponding wavelenght is $`\lambda =\frac{c}{2\pi A}3`$mm. This corresponds to the order of magnitude of the pineal gland. Whether this is just a coincidence or whether it has a deeper meaning is a question that deserves further study. Moreover, it corresponds to a typical frequency already obtained for phonon transitions in the brain, allowing for new theoretical models of the interaction of electromagnetism with the biological cells !

Furthermore, there are bounds on information storage. Theoretically, in a problem of completely different character, one arrives at the maximum entropy a cache can hold, with the result

$$I_{max}<\frac{2\pi R}{\mathrm{}c\mathrm{ln}2}$$
(15)
where $`R`$ is the overall radius of the object under study and $``$ its energy.

In the theory of solitons, the number of possible information-holding configurations based on the soliton equals the number of quanta that might populate the first excited level. To this number we must add unity to account for the soliton configuration itself. So, the possible configurations within an energy budget $``$ above the soliton energy is $`N()=1+[[/\mathrm{}\omega ]]`$ ($`[[x]]`$ stands for the integral part of $`x`$), then

$$I_{max}=ln(1+[[/\omega _1]])log_2ebits$$
(16)
consistent with the bound (15). In our case, for a microtubule, we have $`R10^4`$cm. Taking the energy $``$ as corresponding to a quantum of energy $`\mathrm{}A_0`$, from the source of the electromagnetic field, we find $`=\mathrm{}A_00.2`$eV. In such a case,

$$I_{max}1$$
(17)
while the bound (15) corresponds to $`I_{max}<6`$.

It would not be too original to call such an information a quantum information unit sent via the microtubuiles, in view of similar considerations, in a diferent context by Gabor . In the present case, the soliton, formed by the interaction of Quantum Electrodynamics with the electric dipole moment of the background water in the one-dimensional device offered by microtubules is the natural way chosen by nature to send information bits.

Another interesting point is the fact that the frequency parameters which showed up naturally in the course of the computations have natural interpretations in terms of brain structures. The frequency $`A_03.1\times 10^{14}`$s-1 is compatible with the size of the microtubules. But the real frequency $`\nu \frac{A_0}{6}\sqrt{\frac{\nu _0}{c}}3.7\times 10^{10}`$s-1 is compatible with the transition period observed for the socalled conformational changes connected with tubulin dimer protein (namely $`10^9`$ to $`10^{11}`$s-1) .

This is a further example of the application of Quantum Field Theory to general aspects of matter interactions in complex systems. It is clear that this is not a way of achieving comprehension of the complexity aspects of such a sophisticated system, but it certainly provides valuable tools for working in this field as well.

Acknowledgements: this work has been partially supported by CNPq (Conselho Nacional de Desenvolvimento Científico e Tecnológico) and FAPESP (Fundação de Amparo à Pesquisa do Estado de São Paulo). B.Maroufi would like to thank the Instituto de Física, Universidade de São Paulo, for the hospitality. E.A. thanks Drs. I. Prates de Oliveira and S. F. de Oliveira for discussions concerning aspects of the brain functions.