text:
# Confined water in the low hydration regime.

## I  Introduction

Thin films of water confined in restricted geometries play an important role in many different phenomena in geology, chemistry, biology. Important processes taking place in minerals, plants, biological membranes are determined by the presence of water confined in nanopores inside the system or at contact with the surface. debenedetti ; robinson ; israe  In biology for instance the enzymatic activity of globular proteins or the function of biological membranes are due to the interaction with water debenedetti , and the decrease of the hydration level may inhibit their biological functionality. biol  It is still difficult to predict accurately how the properties of water are modified when it is at contact with the various substrates or confined in different geometries. Nevertheless recent theoretical and computational work on water confined in hydrophobic rigid environments indicate strong changes in its thermodynamical behavior. truskett ; meyer

A strong reduction of the tetrahedral order with consequent distortion of the hydrogen bond network has been evidenced by neutron scattering experiments with isotopic substitution performed on water confined in Vycor mar1 ; mar2 , by experiments performed on freezing of water confined in an hydrophilic environment of mica raviv  and in water confined in vermiculite clay. bergman  These findings are in agreement with recent computer simulation results. brovchenko ; jmoliq ; gr-jcp  Particularly relevant are the observations in experiments and computer simulation of the formation of distinct layer of water with different structural and dynamical properties: an interfacial layer close to the substrate, where water molecules show a slow dynamics, and a layer far from the surface with almost bulk like properties. melni ; kremer ; barut ; mckenna ; noi-epl ; noi-prl ; noi-jcp  The layering effect is expected to be relevant in many phenomena associated with confined water. There are experimental evidences of a slow relaxation of water at contact with proteins doster ; bizzarri1 ; bizzarri2 ; bizzarri3 ; bizzarri4 , while signatures of a similar behavior has been found in water confined in Vycor for low hydrations. venturini

As long as a detailed description of the microscopic properties of water at interface with complex materials like biological macromolecules is difficult to be obtained it could be appropriate to use model systems with well defined geometry and interaction between water and substrate. We present here the results obtained by a molecular dynamics (MD) study of a model of water confined in Vycor at low hydration levels. By computer simulation we can freely vary temperature and hydration level and make an accurate study of the changes induced in the properties of water by the confinement. This work is part of a detailed study of the structural and dynamical properties of confined water. jmoliq ; gr-jcp ; noi-epl ; noi-prl ; noi-jcp

The low hydration regime is somehow more important to be explored at least for two main reasons: (i) the effect of the substrate on the confined water is expected to be enhanced; (ii) the experimental observations are the result of an averaging over different water layers and an accurate analysis of the low hydration regime is the only way to gain some insight on the water-substrate interaction.

An accurate layer analysis enables us to examine the modifications of the main structural properties of confined water and to make connection with recent experimental and computer simulation results which concern water at contact with proteins.

In Sec. II we briefly describe our simulation method, in Sec. III, IV and V we present the main results of our calculation. Sec. VI is devoted to conclusions.

## II  Computer simulation of confined water

The computer simulation is performed in the microcanonical ensemble on water molecules described by the simple point charge/extended (SPC/E) site model potential inserted in a cylindrical cavity of radius $`20`$ Å. The cavity has been carved in a cubic cell of silica glass obtained by computer simulation, as described in previous work. jmoliq ; gr-jcp ; noi-jcp  The system is modeled to represent water confined in Vycor, the internal surface of the cylindrical cell is corrugated and composed by silicon atoms (Si), bridging oxygens (BO) which are bonded to two silicons and non-bridging oxygens (NBO) connected only to one silicon. The NBO are saturated with acidic hydrogens (AH). The interaction of water with the substrate is modeled by assigning different charges to the sites representing the atoms of Vycor. A Lennard-Jones interaction is assumed between the oxygens of water and the BO and NBO of the confining system. jmoliq  Periodic boundary conditions are applied along the axis of the cylinder assumed as the $`z`$ direction. The motion is confined in the $`xy`$ plane, where the distance from the axis is defined in terms of the radius $`R=\sqrt{x^2+y^2}`$. During the simulation the substrate is kept rigid.

The density of the full hydration in the experiments mar1  $`0.0297`$ Å<sup>-3</sup> corresponds in our system to roughly $`N_W=2600`$ molecules. We performed our simulation at three different hydration levels $`N_W=500`$ which corresponds to $`19`$ % of hydration, $`N_W=1000`$ ($`38`$ %) and $`N_W=1500`$ ($`58`$ %). In Fig. 1 the snapshots of the configurations at the three different hydration levels show the hydrophilic character of the substrate.

## III  Density profiles

In the density profiles reported in Fig. 2 it is evident for all the hydrations the formation of a well defined first layer of molecules close to the surface, while a second layer grows up with the hydration level. Few molecules penetrate and are trapped inside the substrate ($`R>20`$ Å). We notice that at the highest hydration ($`N_W=1500`$) the density in the double layer structure increases above the value of density of water at ambient conditions. The decrease of temperature has very little effect on the density profiles, apart from the sharpening of the second layer at the hydration $`N_W=1500`$.

It is clear from Fig. 1 that the arrangement of the molecules close to the Vycor surface is not uniform and the formation of clusters of different sizes is observed, this is even more evident in Fig. 3 where we report an histogram of the distribution of the molecules in the $`xy`$ plane along a circumference at distance $`R=17.5`$ Å where the density profile reaches the maximum. The cluster sizes fluctuate up to hundred molecules with an average value of $`40`$ molecules, this is in agreement with neutron diffraction experiments agamalian  where clusters of the same size are observed at similar value of the hydration.

## IV Layer analysis of the hydrogen bond

Figure 4 shows the average number of hydrogen bond(HB) per molecule along the pore radius for the different hydration levels. We assume that two neighbors water molecules are hydrogen bonded if the $`HO\mathrm{}O`$ angle is less than $`30^o`$jedlov  Two water molecules are considered neighbors if their O-O separation is less than $`3.35`$ Å. The number of water-water HB $`N_{ww}^{HB}(R)`$ decreases approaching the surface, while inside the double layer structure it grows up the number of HB of the water with the Vycor surface $`N_{wv}^{HB}(R)`$. The water molecules form preferentially hydrogen bonds with the bridging oxygens. The arrangement of the clusters close to the Vycor surface is mainly determined by this type of bonds. The number of HB between water molecules and the acidic H atoms is negligible on the scale of the figure. The $`N_{wv}^{HB}(R)`$ curves do not show substantial differences as function of hydration level of the pore.

In Fig. 5 the radial profile of the water-Vycor bonds is compared with the density profile of the bridging oxygens at the surface; by considering that the bond length is approximately $`2`$ Å there is a clear correspondence between the two profiles. It is evident that the bridging oxygens substitute as acceptors the water oxygens approaching the surface since the total number of HB $`N_{tot}^{HB}(R)=N_{ww}^{HB}(R)+N_{wv}^{HB}(R)`$ , also reported in Fig. 4, remains practically constant. The maximum number of water-water HB grows with increasing hydration level since more water molecules become available for the formation of bonds. The maximum is located around $`R=16`$ Å where the density profile shows a minimum (see Fig. 2), in agreement with the fact that the hydrogen bonds increase with decreasing density at room temperature and upon supercooling. geiger

At lower temperatures the number of water-water and water-Vycor hydrogen bonds increases but the profiles are similar to those shown in Fig. 4 at $`T=300`$ K.

In Fig. 6 we report the radial profiles of the average number of water-water hydrogen bonds compared with the profiles of the average number of neighbors $`N_{nn}`$. The ratio between the two quantities $`N_{ww}^{HB}(R)/N_{nn}(R)`$ estimates the number of molecules forming hydrogen bonds between the neighbors. It can be considered as a measure of the rigidity of the hydrogen bond network and it is reported in Fig. 7 for selected values of the hydration and temperature. The ratio is smaller close to the surface ($`R>16`$ Å) where the hydrogen bonds are made mainly with the Vycor surface. Below $`16`$ Å each molecule is able to make hydrogen bonds with almost all its neighbors.

The network becomes more rigid on lowering the temperature (see Fig. 7) since the number of hydrogen bonds increases while the number of neighbors is almost insensible to the variation of temperature, like in the density profiles in Fig. 2. This might be possibly related to the recently proposed fragile to strong transition for strongly confined water. bergman

We can conclude that the water molecules rearrange themselves at decreasing temperature in order to maximize the number of hydrogen bonds. We observe moreover that the rigidity increases at decreasing hydration (see inset of Fig. 7) with some important consequences for the residence times, as discussed later.

We investigated how the short range order present in bulk liquid water is modified upon confinement by looking at the distribution of the hydrogen bond angles. Fig. 8 shows the distribution of $`cos(\theta )`$, where $`\theta `$ is the angle between the intramolecular bond O-H and the intermolecular O$`\mathrm{}`$O vector, for $`N_W=1000`$ and different temperatures. The water-water linear hydrogen bonds are still favored like in bulk water, there is however a more pronounced tail toward large angle with respect to bulk water. By lowering the temperature the distribution peak becomes sharper indicating an increase in the rigidity of the bonds. The distortion in the angular distribution with respect to the bulk is due to the water molecules close to the pore surface, as clearly appears in the inset of Fig. 8, where for $`N_W=500`$ the separated contributions from the layer $`0<R<16`$ Å and $`16<R<20`$ Å are shown. For the molecules close to the surface the peak of the angular distribution is depressed while there is a large increase of the tail at high angle, the angular distribution of the molecules closer to the interior of the pore is similar to the one of the bulk water.

The distortion of the hydrogen bonds is confirmed by experimental studies on water confined in Vycor mar1 ; mar2  and at contact with other surfaces. raviv

We studied by layer analysis the local arrangement of the nearest neighbor molecules by looking to the distribution of the angle $`\gamma `$ between the two vectors joining the oxygen atom of a water molecule with the oxygen atoms of the two closest water neighbors. This distribution in the computer simulation of bulk water at ambient condition has a well defined peak around $`\gamma 109^o`$ ($`\mathrm{cos}(\gamma )0.326`$), this peak is a signature of the tetrahedral order present in liquid water at short range distance. There is also a secondary peak at $`\gamma 54^0`$ corresponding to interstitial neighbor molecules. These findings are in agreement with experimental results on water at ambient conditions.jedlov  In Fig. 9 we report a layer analysis of the distribution function of $`\mathrm{cos}(\gamma )`$ of water in Vycor for different hydrations. The molecules in the layer $`0<R<16`$ Å for $`N_W=1500`$ and $`N_W=1000`$ show a distribution very similar to what found in the simulation of bulk water, the distributions become sharper at lower temperature. For $`N_W=500`$ the peak shifts toward $`90^o`$, the value at which the maximum is located for $`T=220`$ K. This behavior for the lowest hydration is in agreement with experimental results indicating that in confined water freezing in a cubic crystalline phase could be favored. rrgrenoble

As a matter of fact a good short range order is present for all the hydrations in the layer $`0<R<16`$ Å where the water molecules are able to form hydrogen bonds with almost all their neighbors as discussed above (see Fig. 6).

In the layer closest to the Vycor surface $`16<R<18`$ Å the tetrahedral arrangement is completely lost due to the formation of hydrogen bonds between water molecules and the substrate. The temperature has very little effect on these distributions. In the layer $`18<R<20`$ Å the specific shape of the distribution is not very significant and indicates the strong distortion of the hydrogen bonds formed by the molecules found in the inner surface of the pore.

## V Residence time

The dynamics of the molecules of confined water is strongly determined by the interaction with the substrate. In Fig. 10 is reported the trajectory of the center of mass of a molecule which resides at the beginning in the first layer at roughly $`2`$ Å from the surface for the case $`N_W=500`$. At ambient temperature the molecule oscillates mainly between the two layers and persists in a position for fractions of nanoseconds, at lower temperature the molecule oscillates for all the time explored around the same position and does not go out from its initial layer.

Fig. 11 shows how the average residence time of the water molecules changes along the pore radius for different hydrations and temperatures. The residence time for each region is calculated as

$$t_{res}=\underset{i=1}{\overset{N_W}{}}\underset{k=1}{\overset{q_i}{}}\mathrm{\Delta }t_{ik}/\underset{i=1}{\overset{N_W}{}}q_i$$
(1)
where $`\mathrm{\Delta }t_{ik}`$ is the fraction of time that the molecule $`i`$ spent during the $`k`$-th visit in the region. For all hydrations oscillations are present which are particularly marked for $`N_W=500`$. The residence time reaches the maximum values at the surface. The structure of the residence time profiles around the two main peaks is reminiscent of the double layer structure of the density profiles (see Fig. 2). Apart for the layer attached to the Vycor surface the molecules have longer residence time in the regions where the density profiles have local maximal values, while the residence time minimum is found roughly at the minimum of the density distribution between the two layers. The water molecules stay for a shorter time in the range between the two layers where the maximum of the water-water HB profile is located. So in the region of minimum density the HB formed are not very stable.

The maxima of the residence time are enhanced by decreasing the hydration. As discussed in Sec. IV this effect can be explained in terms of an increase of the rigidity of the hydrogen bond network at decreasing hydration.

Close to the substrate where an elevated number of water-Vycor hydrogen bonds $`N_{wv}^{HB}`$ are present a comparison of the average residence time in the zone $`16<R<20`$ Å with $`N_{wv}^{HB}`$ shows that the two quantities are proportional in the range of temperature explored

$$t_{res}(T)=𝒞_NN_{wv}^{HB}(T)$$
(2)
where the constant $`𝒞_N`$ increases at decreasing hydration level going from $`10`$ for $`N_W=1500`$ to $`20`$ for $`N_W=500`$.

The distributions of the residence times in the different layers show peculiar features, as illustrated in Fig. 12 for $`N_W=1500`$ at room temperature. The distributions $`𝒫(t_{res})`$ for the layers more close to the surface are well described by a power law

$$𝒫(t)=at^\mu $$
(3)
where for the molecules in layer $`16<R<20`$ Å $`\mu =1.51`$, while by considering the layer $`14<R<20`$ Å $`\mu =1.36`$. This behavior is a signature of a non brownian motion of the particles and it is observed at all the hydration levels explored in this work. Non brownian diffusion has been found recently in neutron scattering experiment and MD simulation on water at contact with a globular protein bizzarri3  where the distribution of the residence times of the water molecules in the layers close to the protein surface is fitted to Eq. (3) with exponents $`\mu `$ in fairly good agreement with the exponents found in the present work.

For water molecules far from the protein surface it is found bizzarri3  that the brownian regime is established with an exponential decays of the residence time distribution. In the inset of Fig. 12 the residence time distribution is reported for the internal layer $`0<R<16`$ Å of our system, also in our case the decay at long time can be fitted with an exponential

$$𝒫(t)=be^{Bt}$$
(4)
## VI Summary and Conclusions

The structural properties of water confined in a cavity representing the average properties of the pores of Vycor glass have been studied by computer simulation. Water molecules are attracted by the hydrophilic surface and form a double layer structure of $`4÷5`$ Å. In this region it is found that the local properties of water are strongly dependent on the hydration level and the distance from the substrate.

The hydrogen bond network is deformed with respect to bulk water, in agreement with experimental results. mar1 ; mar2 ; raviv  At a distance greater than $`4`$ Å from the surface for the highest hydrations investigated ($`N_W=1500`$ and $`N_W=1000`$) the water molecules are arranged in a local tetrahedral order similar to bulk water, while approaching the substrate the water-Vycor hydrogen bonds substitute the water-water hydrogen bonds and the tetrahedral order is completely lost. At the lowest studied hydration ($`N_W=500`$) the local order far from the surface is compatible with nucleation of a cubic crystalline phase at low temperature. Experimental evidences for nucleation of cubic ice has been indeed found for water under confinement. rrgrenoble

At a given hydration the average number of water-water HB increases with decreasing temperature while the number of nearest neighbors remains constant, this implies that the hydrogen bond network becomes more rigid upon supercooling. The effect is accentuated by decreasing the hydration. This findings are consistent with the recently proposed fragile to strong transition for strongly confined water. bergman

The analysis of the residence time of the water molecules along the pore radius evidenced connections between this quantity and the density profiles, the water molecules persist for longer time in the regions of the maxima of the double layer structure and very close to the Vycor surface.

The residence time distributions show an anomalous non brownian behavior in the layers close to the substrate with a power law decay for long time. On the contrary the distributions decay exponentially in the layers far from the surface. Experimental and computer simulation studies on water at contact with a globular protein bizzarri3  show a similar behavior. This agreement indicates that there are useful analogies between the behavior of water confined in Vycor and the behavior of water at contact with proteins. This is confirmed by a more refined analysis of the dynamical properties of confined water in the low hydration regime. lowhdyn

## VII Acknowledgments

We thank M. A. Ricci for critical reading of the manuscript.


format: text

text:
# Chern-Simons Theory for Magnetization Plateaus of Frustrated 𝐽₁-𝐽₂ Heisenberg model

## Abstract

The magnetization curve of the two-dimensional spin-1/2 $`J_1`$-$`J_2`$ Heisenberg model is investigated by using the Chern-Simons theory under a uniform mean-field approximation. We find that the magnetization curve is monotonically increasing for $`J_2/J_1<0.267949`$, where the system under zero external field is in the antiferromagnetic Néel phase. For larger ratios of $`J_2/J_1`$, various plateaus will appear in the magnetization curve. In particular, in the disordered phase, our result supports the existence of the $`M/M_{\mathrm{sat}}=1/2`$ plateau and predicts a new plateau at $`M/M_{\mathrm{sat}}=1/3`$. By identifying the onset ratio $`J_2/J_1`$ for the appearance of the 1/2-plateau with the boundary between the Néel and the spin-disordered phases in zero field, we can determine this phase boundary accurately by this mean-field calculation. Verification of these interesting results would indicate a strong connection between the frustrated antiferromagnetic system and the quantum Hall system.

Due to quantum and frustration effects, rich physics can appear in the frustrated quantum spin systems at zero external field.frust  Exciting behavior was also observed recently in several cases with an external magnetic field. For instance, the recently discovered two-dimensional (2D) $`S=1/2`$ spin-gap material SrCu2(BO3)2, which can be described by the Shastry-Sutherland model,SS81  exhibits several plateaus at $`M/M_{\mathrm{sat}}=1/3`$, $`1/4`$, and $`1/8`$ in its magnetization curve, where $`M`$ ($`M_{\mathrm{sat}}=1/2`$) is the (saturating) magnetization per site.Kageyama  The origin of these plateaus and the nature of the corresponding spin states are under intense debate.Miyahara ; Momoi ; Muller-Hartmann ; Fukumoto  Recently, by mapping onto spinless fermions carrying one quantum of statistical flux and under a mean-field approximation, Misguich et al. show that the original spin model can be related to a generalized Hofstadter problem, where the spin excitation gaps that produce the observed magnetization plateaus arise from some of the Landau level gaps in the integer quantum Hall effect for the fermions on a lattice.Misguich  For realistic values of the exchange constants, they obtain an excellent quantitative fit to the observed magnetization curve, which demonstrates the success of their approach.

Another prototype of a realistic frustrated two-dimensional system, which has been recently realized experimentally in $`\mathrm{Li}_2\mathrm{VOSiO}_4`$ and $`\mathrm{Li}_2\mathrm{VOGeO}_4`$ compounds,carretta  is the so-called $`J_1`$-$`J_2`$ Heisenberg model with the Hamiltonian

$$H=\underset{i,j}{}J_{ij}\stackrel{}{S}_i\stackrel{}{S}_jB\underset{i}{}S_i^z,$$
(1)
where the exchange couplings $`J_{ij}`$ are equal to $`J_1`$ when $`i,j`$ are nearest neighbors on the square lattice; $`J_{ij}`$ are equal to $`J_2`$ when $`i,j`$ are connected by a diagonal bond. The external magnetic field $`B`$ is applied along the $`z`$-axis. Both couplings are antiferromagnetic, i.e. $`J_{1,2}>0`$, and the spins $`S_i=1/2`$. This model has been the object of intense investigation through years.Chandra ; Read ; series1 ; ed ; Schulz ; series2 ; Capriotti ; lozovik93 ; Gluzman ; Mutter97 ; ZHP ; honecker00 ; Mutter01  At $`B=0`$, for small $`J_2/J_1`$ where the frustration is weak, the system exhibits Néel ordering described by a wave vector $`(\pi ,\pi )`$. When $`J_2/J_1`$ is large enough, the ground state is dominated by the next-nearest-neighbor interactions and has a collinear order described by $`(\pi ,0)`$ or $`(0,\pi )`$. There is also a general consensus on the disappearance of the magnetic ordering at $`0.38<J_2/J_1<0.6`$, while the identification of the ground state is still a subject of much controversy.Chandra ; Read ; series1 ; ed ; Schulz ; series2 ; Capriotti

Just like the case of the Shastry-Sutherland model, the plateaus in the magnetization curve are predicted for the $`J_1`$-$`J_2`$ model,lozovik93 ; Gluzman ; Mutter97 ; ZHP ; honecker00 ; Mutter01  while the situation is even more controversial. Strong evidence for a plateau at $`M/M_{\mathrm{sat}}=1/2`$ in the region $`0.5J_2/J_10.65`$ has been recently reported by Honecker et al.ZHP ; honecker00  However, Fledderjohann and MütterMutter01  do not find a plateau at $`M/M_{\mathrm{sat}}=1/2`$ in the region of the quantum disorder phase; but instead they find some indications for a plateau-like structure at $`M/M_{\mathrm{sat}}=2/3`$. Most of the previous studies on the plateau are obtained by the numerical calculations on a small clusters (the typical number of lattice sites in these works is about $`6\times 6=36`$). As discussed in Ref. Mutter01, , the plateau structures in the magnetization curve can depend sensitively on the system size. Therefore, these results may be plagued with the finite-size effects. For example, some of the predicted plateaus may be an artifact of the special lattice geometry, and the boundary conditions used may frustrate the order which tends to develop. Thus the precise determination of the positions and widths of the plateaus in the frustrated Heisenberg model is indeed a very delicate problem, and a better theoretical understanding of the magnetic order and of the mechanisms which create the plateaus is needed.

To avoid the possible finite-size effects, we apply the Chern-Simons (CS) theory for the magnetization plateausMisguich  to the 2D spin-1/2 $`J_1`$-$`J_2`$ Heisenberg model. Because of its success for the Shastry-Sutherland model, this approach should give us reasonable results in the present case. We find that the magnetization curve is monotonically increasing for $`J_2/J_1<0.267949`$, where the system at zero field is in the antiferromagnetic Néel phase. Besides, various plateaus will appear in the the magnetization curve for larger ratios of $`J_2/J_1`$. In particular, in the disordered phase, our result supports the existence of the $`M/M_{\mathrm{sat}}=1/2`$ plateau and predicts a new plateau at $`M/M_{\mathrm{sat}}=1/3`$. Furthermore, we note that, by identifying the onset value of $`J_2/J_1`$ for the appearance of the 1/2-plateau with the critical value of the phase transition between the Néel and the spin-disordered phases in zero field, this phase boundary, which was determined earlier by heavy numerical means,series2 ; Capriotti  can be reproduced accurately by this mean-field calculation.

The Hamiltonian in Eq. (1) can be rewritten as

$`H=H_{xy}+H_zB{\displaystyle \underset{i}{}}S_i^z,`$
$`H_{xy}={\displaystyle \frac{1}{2}}{\displaystyle \underset{i,j}{}}J_{ij}\left(S_i^+S_j^{}+S_j^+S_i^{}\right),`$ (2)
$`H_z={\displaystyle \underset{i,j}{}}J_{ij}S_i^zS_j^z,`$
where $`H_{xy}`$ and $`H_z`$ are the spin-flip and the Ising parts of the Hamiltonian at zero field. According to Ref. Misguich, , $`H_{xy}`$ and $`H_z`$ are treated in different ways. First, the Ising part is approximated by a simple uniform mean-field decoupling,

$$H_z4(J_1+J_2)M\underset{i}{}S_i^z2(J_1+J_2)M^2N,$$
(3)
where $`N`$ is the number of lattice sites. Thus $`H_z`$ gives a contribution to the total energy with a simple dependence on the magnetization. Second, for $`H_{xy}`$, one can exactly map the spin operators to spinless fermions attached with a flux tube carrying one flux quantum of statistical CS magnetic field,CS ; ywg93 ; lrf94  where $`S_i^z+1/2`$ corresponds to the occupation number $`n_i`$ of site $`i`$. Under a mean-field treatment such that the flux tubes are smeared out into a uniform background magnetic field, the flux per square plaquette $`\varphi `$ is then tied to the density of fermions and thus to the magnetization $`M`$ of the spin system :

$$\frac{\varphi }{2\pi }=n=M+\frac{1}{2}.$$
(4)
Because of this flux, each energy band splits to subbands with a complicated structure. Therefore, the present spin system can be identified with a Hofstadter problemhofstadter76  for fermions moving on a square lattice with nearest-neighbor and next-nearest-neighbor hoppings.HK ; LCH  For a given $`M`$ \[or $`\varphi =2\pi (M+1/2)`$\], the mean-field ground state is obtained by filling the lowest energy subbands with fermions until their density satisfies $`n=\varphi /2\pi `$. The one-body problem from $`H_{xy}`$ can be straightforwardly analyzed for rational values of $`\varphi /2\pi `$. For $`\varphi /2\pi =p/q`$ ($`p`$ and $`q`$ are mutually prime integers), there are $`q`$ subbands and the ground state is the Slater determinant with the lowest $`p`$ subbands being completely filled. The energy of the filled subbands leads to another contribution to the total energy. Hence the total energy per site $`E(M)`$ as a function of the magnetization becomes

$$E(M)=\frac{1}{N}\underset{\alpha =1,\mathrm{},p}{}\underset{\stackrel{}{k}}{}ϵ_\stackrel{}{k}^{(\alpha )}+2(J_1+J_2)M^2.$$
(5)
Here $`ϵ_\stackrel{}{k}^{(\alpha )}`$ is the eigenenergy of the $`\alpha `$-th subband with the wavevector $`\stackrel{}{k}`$ being restricted to the magnetic Brillouin zone. The magnetization can be obtained as a function of $`B`$ by minimizing $`E(M)BM`$. It is clear from Eq. (5) that, without the contribution from $`H_{xy}`$, the magnetization $`M`$ will be linearly proportional to $`B`$, and there is no magnetization plateau. Therefore, the appearance of magnetization plateaus is related to certain features of the Hofstadter spectrum.

The Hofstadter diagrams for $`J_2/J_1`$=0.2 and 0.3 are shown in Fig. 1, where the lower bold line marks the Fermi level (highest occupied state) and the upper one marks the lowest unoccupied level. A jump of the Fermi energy as a function of $`M`$ in Fig. 1(b) leads to discontinuity of the slope of the function $`E(M)`$. These jumps for various $`M`$ give rise to plateaus in the magnetization curve. They are closely related to the occurrence of band-crossing when the value of $`J_2/J_1`$ is varied. For example, in Fig. 1(a) before the upper two subbands at $`\varphi /2\pi =2/3`$ touch at $`J_2/J_1=0.267949`$, the Fermi energy is continuous and there is no magnetization plateau at $`M/M_{\mathrm{sat}}=1/3`$ (see Fig. 2). However, in Fig. 1(b), the “pockets” enclosed by the bold lines are separated at $`\varphi /2\pi =2/3`$. It can be seen that the Fermi level to the right of the contact point is below (above) the band gap before (after) band-crossing. As discussed earlier, the Fermi level marks the position of the $`p`$-th subband \[see Eq. (5)\]. Therefore, apparently a subband associated with some flux slightly larger than 2/3 is shifted above the energy gap after band-crossing. This shift of a fine subband due to the crossing of the broader subbands at $`\varphi /2\pi =2/3`$ was studied earlier in the context of the Hofstadter spectrum.LCH  It is closely related to the jump of the integer-valued Hall conductances of the broader subbands.LCH ; chang  The emergence of the magnetization plateaus for the spin system thus has an interesting connection with the change of the integer-valued Hall conductances induced by band-crossing.

To justify the present approach, it is important to check whether the plateau states are robust against fluctuations around the mean-field solutions. It was showed that the Gaussian fluctuations of the CS gauge field are massless only when the TKNN integertknn82  describing the quantized Hall coefficient of the fermions on the frustrated lattice becomes unity.ywg93  In that case, the Gaussian fluctuations induce instabilities for the mean-field ground states. We have computed the TKNN integers numerically (for example, see Fig. 1) and found the Gaussian fluctuations to be massive. Therefore, the plateaus are not destroyed by quantum fluctuations.

Magnetization curves for various $`J_2/J_1`$ ratios are shown in Fig. 2. We note that the saturation field $`B_{\mathrm{sat}}`$ can be computed exactly by identifying the energy $`E_f`$ of the fully polarized state with the (exact) minimum energy $`E_{1s}^{\mathrm{min}}`$ of the states with one spin flipped, $`E_f=E_{1s}^{\mathrm{min}}`$. Thus $`B_{\mathrm{sat}}/J_1=4`$ for $`J_2/J_1<1/2`$; $`B_{\mathrm{sat}}/J_1=2+4J_2/J_1`$ for $`J_2/J_1>1/2`$. Our findings agree with these exact results near the saturation.

In the Néel phase, it is expected that the spins cant gradually from the antiparallel configuration toward the parallel configuration until the magnetization saturates at the saturation field $`B_{\mathrm{sat}}`$. The magnetization curve obtained from the present approach is consistent with this expectation: it is featureless all the way to full saturation when $`J_2/J_1`$ is small (see the curves for $`J_2/J_1=0`$ and 0.2 in Fig. 2). Upon increasing $`J_2/J_1`$, plateaus emerge and the magnetization curves become more complex. In particular, because of the band-crossing at $`\varphi /2\pi =2/3`$ when $`J_2/J_1=0.267949`$ (see Table I of Ref. LCH, ), a plateau at $`M/M_{\mathrm{sat}}=1/3`$ is found \[see Eq. (4)\]. This is an unexpected result, especially for $`0.267949<J_2/J_10.38`$ where the ground state in the absence of an external field is in the Néel phase. In the previous finite-size studies,lozovik93 ; Gluzman ; Mutter97 ; ZHP ; honecker00 ; Mutter01  there is no indication for the appearance of this plateau. However, the system sizes and the boundary conditions they used forbid the appearance of the $`M/M_{\mathrm{sat}}=1/3`$ plateau, therefore the possibility of this plateau cannot be ruled out. Furthermore, when $`J_2/J_1=0.382683`$, another band-crossing in the Hofstadter spectrum occurs at $`\varphi /2\pi =3/4`$ (see Table I of Ref. LCH, ), and a plateau at $`M/M_{\mathrm{sat}}=1/2`$ ensued. It is interesting to find that the critical value for the appearance of the $`M/M_{\mathrm{sat}}=1/2`$ plateau agrees quite well with the critical point of the quantum phase transition at zero field between the Néel and the quantum disorder phases.series2 ; Capriotti  This reinforces our confidence on the present CS mean-field approach. As mentioned before, while the appearance of a plateau in the quantum disorder phase had been predicted, the value of the plateau is still under debate.ZHP ; honecker00 ; Mutter01  The controversy may come from the subtle finite-size effects in their investigations. Since the present CS theory is free from the finite-size effects, we give a strong support for the existence of the $`M/M_{\mathrm{sat}}=1/2`$ plateau.

More complex structures in the magnetization curves appear when $`J_2/J_1`$ is further increased. For example, when $`J_2/J_1=1/2`$, a series of plateaus at $`M/M_{\mathrm{sat}}=n/(n+2)`$ is found, which corresponds to the band-crossing at some magic numbers $`\varphi /2\pi =(n+1)/(n+2)`$. Irregular plateau structures are found for even higher values of $`J_2/J_1`$ (see the $`J_2/J_1=0.7`$ curve in Fig. 2). This behavior is quite similar to the case of the triangular lattice, where many plateaus are predicted under the uniform CS mean-field approximation.Misguich  In the case of the triangular lattice, it is shown that, when the nonuniform mean-field solutions are used, only the main plateau ($`M/M_{\mathrm{sat}}=1/3`$ in that case) survives and other mini-plateaus disappear. We believe that the same situation will happen in the present study of the $`J_1`$-$`J_2`$ model. That is, when going beyond the uniform mean-field approximation, many of the mini-plateaus in the magnetization process may disappear and only the main plateaus with simple fractions of $`M/M_{\mathrm{sat}}`$ survive. Moreover, as shown in Fig. 2, the spin gap (the plateau at $`M/M_{\mathrm{sat}}=0`$) opens at $`J_2/J_1=0`$, instead of opening at the correct critical coupling $`J_2/J_10.38`$. This feature comes from the fact that the Néel state of the square-lattice antiferromagnet is not correctly described by the present uniform mean-field approximation.

In conclusion, the magnetization curve of the 2D spin-1/2 $`J_1`$-$`J_2`$ Heisenberg model is studied by using the Chern-Simons theory under a uniform mean-field approximation. In the disordered phase, our result supports the existence of the $`M/M_{\mathrm{sat}}=1/2`$ plateau and predicts a new plateau at $`M/M_{\mathrm{sat}}=1/3`$. Moreover, various plateaus appear in the magnetization curve both in the disordered and the collinear phases, which could be an artifact of our mean field approximation. More work is needed to be conclusive at this point. We note that it is experimentally accessible to confirm our results. As mentioned before, the 2D spin-1/2 $`J_1`$-$`J_2`$ Heisenberg model with $`J_2/J_11`$ has been recently realized experimentally in $`\mathrm{Li}_2\mathrm{VOSiO}_4`$ ($`J_1+J_28.2`$ K) and $`\mathrm{Li}_2\mathrm{VOGeO}_4`$ ($`J_1+J_28.2`$ K) compounds.carretta  If the $`g`$-factor is taken to be 2, the corresponding saturation fields will be approximately 18 T for $`\mathrm{Li}_2\mathrm{VOSiO}_4`$ and 13 T for $`\mathrm{Li}_2\mathrm{VOGeO}_4`$. In both cases, these values of the magnetic fields can be reached experimentally. Thus the full magnetization curve could be mapped out to test our results. Verification of the magnetization plateaus would indicate a strong connection between the frustrated antiferromagnetic system and the quantum Hall system.

M.C.C. is supported by the National Science Council of Taiwan under Contract No. NSC 90-2112-M-003-022. M.F.Y. acknowledges financial support by the National Science Council of Taiwan under Contract No. NSC 90-2112-M-029-004.


text:
# Flow Measurements via Two-particle Azimuthal Correlations in Au + Au Collisions at √𝑠_{𝑁𝑁}=130 GeV.

## Abstract

Two particle azimuthal correlation functions are presented for charged hadrons produced in Au + Au collisions at RHIC ($`\sqrt{s_{_{NN}}}=130`$ GeV). The measurements permit determination of elliptic flow without event-by-event estimation of the reaction plane. The extracted elliptic flow values ($`v_2`$) show significant sensitivity to both the collision centrality and the transverse momenta of emitted hadrons, suggesting rapid thermalization and relatively strong velocity fields. When scaled by the eccentricity of the collision zone $`\epsilon `$, the scaled elliptic flow shows little or no dependence on centrality for charged hadrons with relatively low $`p_T`$. A breakdown of this $`\epsilon `$ scaling is observed for charged hadrons with $`p_T>`$ 1.0 GeV/c for the most central collisions.

The primary goal of current relativistic heavy ion research is the creation and study of nuclear matter at high energy densities. Open questions include the detailed properties of such excited matter, as well as the existence of a transition to the quark-gluon plasma (QGP) phase. Such a phase of deconfined quarks and gluons has been predicted to survive for $`310`$ fm/c in Au + Au collisions at the Relativistic Heavy Ion Collider (RHIC), and several experimental probes have been proposed for its possible detection and study. Elliptic flow constitutes an important observable because it is thought to be driven by pressure built up early in the collision, and therefore can reflect conditions existing in the first few fm/c. Elliptic flow leads to an anisotropy in the azimuthal distribution of emitted particles. A Fourier decomposition of this distribution

$$\frac{dN}{d(\varphi \mathrm{\Phi }_R)}(1+\underset{n=1}{\overset{\mathrm{}}{}}2v_n\mathrm{cos}(n(\varphi \mathrm{\Phi }_R))),$$
(1)
provides a characterization of the elliptic flow via the second Fourier coefficient $`v_2`$. Here, $`\varphi `$ is the azimuthal angle of an emitted particle and $`\mathrm{\Phi }_R`$ is the azimuth of the reaction plane, defined by the beam direction and the impact parameter vector. The apparent reaction plane is determined from the azimuthal correlations between emitted particles, and its dispersion correction from the azimuthal correlations between two “subevents”. An alternative technique for elliptic flow analysis, is the Fourier decomposition of the pair-wise distribution in the azimuthal angle difference ($`\mathrm{\Delta }\varphi =\varphi _1\varphi _2`$) between pairs of emitted particles:

$$\frac{dN}{d\mathrm{\Delta }\varphi }(1+\underset{n=1}{\overset{\mathrm{}}{}}2v_n^2\mathrm{cos}(n\mathrm{\Delta }\varphi )).$$
(2)
In this case the magnitude of the elliptic flow is characterized by the square of the second Fourier coefficient in Eq. (1), i.e. $`v_2^2`$. These methods of analysis can be taken as equivalent since (i) the correlation between every particle and the reaction plane induces correlations among the particles, and (ii) correlating two subevents amounts to summing two-particle correlations. The results in this letter have several advantages over elliptic flow measurements already performed at the same beam energy by the STAR and PHOBOS collaborations. First, two particle correlations circumvent the need for full azimuthal detector acceptance. Second, it allows the determination of elliptic flow without event-by-event estimation of the reaction plane and the associated corrections for its dispersion. Third, this method can serve to minimize many important systematic uncertainties (detector acceptance, efficiency, etc) which may influence the accuracy of elliptic flow measurements.

Elliptic flow is predicted and found to be negative for beam energies $`4`$ AGeV and positive for higher beam energies in Au + Au collisions. Recent theoretical investigations have made predictions for the centrality dependence of the scaled elliptic flow $`A_2v_2/\epsilon `$  where $`\epsilon `$ is the eccentricity or initial spatial anisotropy of the “participant” nucleons in the collision zone. The eccentricity $`\epsilon `$ shows an essentially linear variation with impact parameter $`b`$, for $`0.2b_{max}b0.8b_{max}`$ in Au + Au collisions ($`b_{max}14`$ fm). For central collisions ( b$``$ 5 - 6 fm), it is predicted that higher energy densities are produced and rapid kinetic equilibration in the QGP phase leads to a characteristic rise in $`A_2`$. In addition, elliptic flow for high $`p_T`$ particles has been proposed as an observable sensitive to the energy loss of scattered partons in a QGP phase .

The colliding Au beams ($`\sqrt{s_{_{NN}}}=130`$ GeV) used in these measurements have been provided by the Relativistic Heavy Ion Collider at Brookhaven National Laboratory (BNL). Charged reaction products were detected in the east and west central arms of the PHENIX detector. Each of these arms subtends 90o in azimuth $`\varphi `$, and $`\pm 0.35`$ units of pseudo-rapidity $`\eta `$. The axially symmetric magnetic field of PHENIX (0.5 T) allowed for the tracking of particles with $`p_T200`$ MeV/c in the fiducial volume of both arms. The Drift Chamber (DC) and a layer of Pad Chambers (PC1) located at radii of 2 m and 2.5 m respectively, in each arm, served as the primary tracking detector for these measurements. A second layer of Pad Chambers (PC3), positioned at 5 m in the east arm, was employed to confirm the trajectory of charged particles which traversed both the DC and PC1. The Zero Degree Calorimeters (ZDC), were used in conjunction with the Beam-Beam Counters (BBC) to provide the position of the vertex along the beam direction as well as a trigger for a wide range of centrality selections.

The present data analysis uses two-particle azimuthal correlation functions to measure the distribution of the azimuthal angle difference ($`\mathrm{\Delta }\varphi =\varphi _1\varphi _2`$) between pairs of charged hadrons. Following an approach commonly exploited in interferometry studies, a two-particle azimuthal correlation function can be defined as follows

$$C(\mathrm{\Delta }\varphi )=\frac{N_{cor}(\mathrm{\Delta }\varphi )}{N_{uncor}(\mathrm{\Delta }\varphi )},$$
(3)
where $`N_{cor}(\mathrm{\Delta }\varphi )`$ is the observed $`\mathrm{\Delta }\varphi `$ distribution for charged particle pairs selected from the same event, and $`N_{uncor}(\mathrm{\Delta }\varphi )`$ is the $`\mathrm{\Delta }\varphi `$ distribution for particle pairs selected from mixed events. Events were selected with a collision vertex position, $`20<`$ z $`<20`$ cm, along the beam axis. Mixed events were obtained by randomly selecting each member of a particle pair from different events having similar multiplicity and vertex position. In order to suppress an over-efficiency in finding two tracks at close angles, hadron pairs within 1 cm of each other in the DC were removed from both the $`N_{cor}(\mathrm{\Delta }\varphi )`$ and $`N_{uncor}(\mathrm{\Delta }\varphi )`$ distributions. Event centralities were obtained via a series of cuts in the space of BBC versus ZDC analog response; they reflect percentile selections of the total interaction cross section of 6.8 barns. Estimates for the impact-parameter and the eccentricity, $`\epsilon (<y^2><x^2>)/(<y^2>+<x^2>)`$ were also made for each of these selections following the model detailed in Ref.. Here, $`<\mathrm{}>`$ represents the spatial average (weighted by the density) of participant nucleons over the transverse plane of the collision zone. Systematic uncertainties associated with the determination of $`\epsilon `$ are estimated to be $``$ 7%.

Correlation functions were obtained via two separate methods. In the first, charged hadron pairs were formed by selecting each particle from a common $`p_T`$ range (fixed-$`p_T`$ method). In the second hadron pairs were formed by selecting one member from a fixed $`p_T`$ range and the other from outside this range (assorted-$`p_T`$ method). Within statistical errors, both methods of analysis yield similar results for the $`p_T`$ range presented below, as would be expected for a system dominated by collective motion.

An important prerequisite for reliable flow extraction from PHENIX data is to establish whether or not the $`180^o`$ azimuthal coverage of the detector results in significant distortions to the correlation function. To this end, detailed simulations of the detector response, acceptance and efficiency, have been performed for simulated data incorporating specific amounts of flow (parameterized by $`v_2`$). The results from these simulations indicated no significant distortion to the correlation functions due to the PHENIX acceptance. On the other hand, small distortions to the correlation function (for $`\mathrm{\Delta }\varphi 25^o`$) as well as an incomplete recovery of $`v_2`$ could be attributed to background contributions.

These background contributions principally affect the extraction of $`v_2`$ from the correlation function in two ways. The distortion to the correlation function at small relative angle introduces a small systematic distortion when fit with a Fourier function \[c.f. Eq. (2)\]. A good representation of the data was obtained with the fit function $`C(\mathrm{\Delta }\varphi )=\lambda exp(0.5(\mathrm{\Delta }\varphi /\sigma )^2)+a_1(1+2v_2^2cos(2\mathrm{\Delta }\varphi ))`$, where the Gaussian term is used to characterize the background distortion at small angles. In addition, there is an isotropic background of false tracks which are predominantly misidentified as high $`p_T`$ particles. These contributions can be efficiently suppressed in the east central arm of PHENIX, by requiring a relatively stringent association between tracks found in the DC and their associated hits in PC3. Using such a procedure, the fraction of background tracks has been evaluated as a function of $`p_T`$ and used to correct $`v_2`$. Corrections range from $``$ 10% at low $`p_T`$ to $``$ 25% at 2 GeV/c with a systematic uncertainty of 5%. The current analysis is restricted to the range $`0.3<p_t2.5`$ GeV/c to maintain this relatively small systematic uncertainty.

Figures 1a - d show representative $`\mathrm{\Delta }\varphi `$ correlation functions obtained for charged hadrons detected in the pseudo-rapidity range $`0.35<\eta <0.35`$. Correlation functions for relatively central events (centrality = 20 - 25%, b $``$ 7.0 fm) are shown for hadrons with $`0.3<p_T<2.5`$ GeV/c and $`0.5<p_T<2.5`$ GeV/c in Figs. 1a and c respectively. The same $`p_T`$ selections have been made for the correlation functions shown in Figs. 1b and d but for more peripheral collisions ( centrality = 40 - 45%, b $``$ 9.6 fm) as indicated. Figs. 1a - d show a clear anisotropic pattern which is essentially symmetric about $`\mathrm{\Delta }\varphi =90^o`$. There is also a visible increase of this anisotropy with increasing impact parameter and $`p_T`$. These trends are all consistent with those expected for in-plane elliptic flow.

The magnitude of elliptic flow and the mechanism for its development, can be shown to be related to (a) the geometry of the collision zone, (b) the initial baryon and energy density developed in this zone, and (c) the detailed nature of the equation of state for the created nuclear matter. Since differential flow measurements can serve to provide important insights for disentangling these separate aspects, we show the results of such measurements in Figs. 2 and 3. Fig. 2 shows $`v_2`$ as a function of centrality for several $`p_T`$ selections; $`0.40<p_T<0.60`$ (diamonds), $`0.60<p_T<1.00`$ (squares), and $`1.0<p_T<2.5`$ (circles) GeV/c respectively. Fig. 3 compares the differential flow $`v_2(p_T)`$, for several centralities as indicated.

Figs. 2 and 3 both show relatively large differential flow values which increase with increasing impact parameter and the $`p_T`$ of emitted hadrons. The separate effects of spatial asymmetry and the response of the collision zone to the generated pressure are also evident in Fig. 3. That is, $`v_2`$ not only increases with increasing impact parameter for a fixed $`p_T`$, but also increases with increasing $`p_T`$ for a fixed centrality selection. Trivially, the magnitude of flow should go to zero for very small and very large impact parameters. Similarly, its magnitude can be expected to be zero for $`p_T=0`$. It is interesting that Fig. 3 indicates an essentially linear rise of $`v_2`$ with $`p_T`$ for each of the centrality selections presented. Such a trend can not be accounted for via simple geometric considerations alone. However, it is compatible with model calculations which assume a strong transverse velocity field. This suggests the presence of strong dynamically driven transverse flow at RHIC. The magnitude and trends for $`v_2`$ shown in Figs. 2 and 3 are consistent with other elliptic flow measurements at RHIC.

Figure 4 aims to disentangle the geometric and dynamical ($`p_T`$) contributions to the elliptic flow over a broad range of centralities or energy densities. To do this, we plot $`A_2`$ as a function of centrality to obtain the dynamical contributions . This evaluation is performed for two $`p_T`$ selections ($`0.40<p_T<0.60`$ and $`1.0<p_T<2.5`$ GeV/c) which give rise to relatively low and high $`p_T`$ values respectively. The underlying idea is that this ratio should remove the geometric dependence of $`v_2`$, while the $`p_T`$ selections may provide greater sensitivity to different time scales and energy densities associated with the expanding system.

Figure 4 shows an increase in the magnitude of $`A_2`$ with increasing $`p_T`$. This increase can be attributed to the dynamical response of the created system, resulting from the generated pressure gradients. For hadrons of $`0.4<p_T<0.6`$ the observed centrality dependence of $`A_2`$ is compatible with $`\epsilon `$ scaling. However, a breakdown of this scaling is observed for hadrons with $`1.0<p_T<2.5`$. Such a trend may point to a change in the particle production mechanism or the possibility that pressures larger than those predicted by current hydrodynamic calculations  are being produced in the most central collisions at RHIC. It is also interesting to note that the species composition of the charged particle spectra changes dramatically between the two $`p_T`$ ranges used in Fig. 4.

To summarize, we have measured two-particle azimuthal correlation functions for charged hadrons produced in Au + Au collisions at RHIC ($`\sqrt{s_{_{NN}}}=130`$ GeV). The integral, differential and scaled elliptic flow values extracted from these measurements indicate strong sensitivity to the collision centrality and the transverse momenta of emitted hadrons. The centrality dependence of $`v_2`$ suggests that the high-energy-density nuclear matter created at RHIC, efficiently translates the initial spatial asymmetry into a similar asymmetry in momentum space. The $`p_T`$ dependence is consistent with the development of strong transverse velocity fields in the collision zone. The centrality dependence of $`A_2`$ for hadrons in the range $`0.4<p_T<0.6`$ is compatible with $`\epsilon `$ scaling. However, a breakdown of this scaling is observed for hadrons with $`1.0<p_T<2.5`$. Such a trend could result from a number of effects, the most intriguing of which is a possible change in the equation of state. Additional experimental signatures and model calculations will undoubtedly be necessary to test the detailed implications of these results. Nevertheless, the results presented here clearly show that two particle azimuthal correlation measurements provide an important probe for the high-energy-density nuclear matter created at RHIC.

We thank the staff of the Collider-Accelerator and Physics Departments at BNL for their vital contributions. We acknowledge support from the Department of Energy and NSF (U.S.A.), MEXT and JSPS (Japan), RAS, RMAE, and RMS (Russia), BMBF, DAAD, and AvH (Germany), VR and KAW (Sweden), MIST and NSERC (Canada), CNPq and FAPESP (Brazil), IN2P3/CNRS (France), DAE and DST (India), KRF and CHEP (Korea), the U.S. CRDF for the FSU, and the US-Israel BSF.


text:
# First-Order Phase Transitions in Frustrated Spin Systems

## 1 Introduction

Recently antiferromagnetic quantum spin systems with strong frustration have been studied intensively. A typical example is $`\mathrm{SrCu}_2(\mathrm{BO}_3)_2`$ where the spin-1/2 magnetic Cu ions are located on the orthogonal-dimer lattice (the so-called Shastry-Sutherland model). Various interesting phenomena were observed for this compound such as dispersionless triplet excitations, correlated hopping of two magnons, plateaus in the magnetization curve. It is suggested that novel properties in this compound are closely related to frustration due to the competing antiferromagnetic couplings. Besides such spin systems, metallic systems with frustration have also attracted much attention recently. For instance, it is suggested that heavy-fermion behavior observed in $`\mathrm{LiV}_2\mathrm{O}_4`$  may be caused by geometrical frustration inherent in the pyrochlore-lattice structure.

One of the remarkable properties common in such frustrated systems is that strong frustration enhances the competition among several eigenstates with distinct character as a candidate for the ground state, giving rise to first-order quantum phase transitions when we vary the exchange couplings, the magnetic field, the chemical potential, etc. In order to address the role of geometrical frustration in quantum phase transitions, in this paper, we investigate a simple one-dimensional (1D) version of the frustrated spin model with orthogonal-dimer structure. By means of the series expansion technique, the exact diagonalization (ED) and the density matrix renormalization group (DMRG), we demonstrate that strong frustration indeed causes a wide variety of first-order quantum phase transitions.

## 2 Quantum spin system with orthogonal-dimer structure

In this section, we study first-order transitions in the 1D orthogonal-dimer spin chain.

The Hamiltonian we shall deal with is the standard $`S=1/2`$ antiferromagnetic Heisenberg model in a magnetic field $`H`$, $`=J_{ij}𝐒_i𝐒_jHS_i^z`$, on the orthogonal-dimer chain shown in Fig. 1. We will use the normalized parameters $`j=J^{}/J`$ and $`h=H/J`$ for convenience.

We first summarize the ground state properties of the orthogonal-dimer spin chain. The remarkable point characteristic of the orthogonal-dimer system is that a direct product of local dimer-singlet states indicated by the solid line in Fig. 1 is always an exact eigenstate of the Hamiltonian. Therefore, the dimer-singlet state should be the exact ground state up to a certain critical value of $`j`$.

On the other hand, in the large $`j`$ region, another singlet state is stabilized, which is adiabatically connected to the isolated plaquette-singlet denoted by the broken line in Fig. 1. As seen from Fig. 2 (a), a first-order quantum phase transition occurs between the dimer and the plaquette phases at the critical point, $`j_c=0.81900`$.

Let us now turn to the spin excitation spectrum. In the dimer phase $`(j<j_c)`$, the orthogonal-dimer structure forbids free hopping of a triplet excitation created over the local dimer state. On the other hand, in the plaquette phase $`(j>j_c)`$ there are two kinds of low-energy magnetic excitations. One is a dispersive triplet excitation, which is formed by simply breaking the plaquette singlet state. The corresponding excitation energy is shown as the solid line in Fig. 2 (b). The other is an unusual four-fold degenerate excitation shown as the dashed line, which is formed by breaking the plaquette singlet state first, and then making a local dimer state accompanied by two free spins. Since the latter excitation becomes the lowest-excitation in the region $`j_c<j<0.87`$, as shown in Fig. 2, it plays an essential role for the phase transitions in the vicinity of the critical point, as will be shown below. In this way, the unusual four-fold degenerate excitation found here possesses intermediate properties between those typical for the dimer phase and the plaquette phase.

We have seen so far that several different kinds of low-energy excitations coexist in our frustrated model. It is thus interesting to see how these excitations affect the magnetization process. To this end, we focus here on the magnetization curve for the ratio of the exchange couplings $`j=0.94`$, where the system belongs to the plaquette phase. As shown in Fig. 3, the lowest excitation is a plaquette-triplet excitation (solid line), and another four-fold dispersionless excitation (broken line) lies slightly above it.

The coexistence of two-kinds of distinct excitations influences the magnetization curve considerably. As $`h`$ increases beyond the critical field $`(h_c=\mathrm{\Delta })`$, the magnetization should develop with $`M(hh_c)^{1/2}`$, since it is dominated by the triplet excitation whose dispersion relation is quadratic near the bottom $`(q=\pi )`$. With slight increase in $`h`$, however, the magnetization suddenly jumps to the $`1/4`$-plateau. This first-order transition results from the coexistence of two distinct levels in low-energy excitations. Therefore, it is instructive to remark that the 1/4-plateau is not generated by the crystallization of the lowest triplet excitation, but of the four-fold degenerate excitation. Beyond the 1/4-plateau, we again encounter another first-order phase transition accompanied by a jump in $`M`$.

As seen in this section, nontrivial behavior shows up reflecting the competing interactions in the present model. We wish to observe below what happens for such 1D spin systems, when we introduce the charge degree of freedom, which is important to discuss a metallic system with strong frustration.

## 3 Correlated electron system with strong frustration

In this section, we extend our discussions to a hole-doped system. We here employ another orthogonal-dimer model  shown in Fig. 4. A nice feature in this model is that it can describe not only the orthogonal-dimer spin system but also the linked-tetrahedra system, which plays an important role for understanding low-energy properties of the pyrochlore system. In the undoped case, a first-order quantum phase transition between the dimer phase and the plaquette phase occurs at the critical point $`j_c=0.71`$.

We shall deal with a hole-doped system by introducing the t-J model with electron hopping, $`t`$ and $`t^{}`$, and clarify how frustration affects the nature of the metal-insulator transition. Let us first consider the hole-doping effects on the dimer phase. For simplicity, we set $`t/J=t^{}/J^{}=4.0`$.

In Fig. 5 (a), the electron density $`n`$ calculated by DMRG is shown as a function of the chemical potential $`\mu `$ for $`j=0.4`$. In this figure, $`n=1`$ corresponds to half-filling at which the system is in the insulating phase with singlet-dimer state. It is seen that hole-doping smoothly drives the system to a metallic phase. At quarter filling ($`n=1/2`$), we encounter another insulating phase, which is stabilized by the formation of a CDW state.

In order to see the nature of spin correlations in the metallic phase, we compute the short-range spin correlation functions in the rung and chain directions, which are defined as (see Fig. 4),

$`\overline{S}_1`$ $`=`$ $`{\displaystyle \frac{1}{2}}(<S_1^zS_2^z>+<S_3^zS_4^z>),`$
$`\overline{S}_2`$ $`=`$ $`{\displaystyle \frac{1}{4}}(<S_1^zS_3^z>+<S_1^zS_4^z>+<S_2^zS_3^z>+<S_2^zS_4^z>).`$
In Fig. 5 (b), the spin correlation functions calculated by DMRG are shown as a function of the chemical potential $`\mu `$. Note that we have $`\overline{S}_1=1/4`$ and $`\overline{S}_2=0`$ at half filling, which accord with those expected for the isolated dimers. It is seen that the spin correlation along the rung, $`\overline{S}_1`$, rapidly decreases down to $`1/4`$ as the system approaches half filling from the metallic side. On the other hand, $`\overline{S}_2`$ becomes almost zero along the chain direction, reflecting the fact that the system is an assembly of decoupled dimers at half filling. Therefore, near half filling, the metallic state is considered as a resonating state composed of such dimer pairs.

In contrast to the dimer phase, strong frustration dramatically changes the nature of the metal-insulator transition in the plaquette phase.

In Fig. 6 (a), we show the electron density $`n`$ as a function of $`\mu `$ for $`j=1.0`$, where the system is in the plaquette-singlet phase at half filling. It should be noticed that as $`\mu `$ decreases from half filling ($`n=1`$), a first-order phase transition suddenly occurs to the metallic phase at $`\mu 2.2`$, which is accompanied by the discontinuity in the electron density. To clarify the nature of the metal-insulator transition in more detail, we again discuss the short-range spin correlation function \[see Fig. 6 (b)\]. In the insulating phase ($`n=1`$) for $`\mu >2.2`$, the ground state is in the plaquette-singlet phase, which is characterized by $`\overline{S}_1>0`$ and $`\overline{S}_2<0`$ in accordance with the numerical results. On the other hand, once the first-order transition occurs, both of $`\overline{S}_1`$ and $`\overline{S}_2`$ have negative values in the metallic phase. Note that $`\overline{S}_1`$ decreases with negative values while $`\overline{S}_2`$ gets very small as the system approaches half-filling from the metallic side. This implies that the dimer correlation along the rung is largely enhanced in the metallic phase close to half filling. This dimer-dominant state competes with the plaquette-singlet state realized at half filling, triggering the first-order phase transition.

## 4 Summary

We have investigated the first-order quantum phase transitions in antiferromagnetic spin chains with orthogonal-dimer structure. It has been clarified that there appears rich structure in the excitation spectrum, which is caused by strong frustration. Such nontrivial excitations determine the nature of the quantum phase transitions, when we change the interactions, the magnetic field, etc. We have also studied the hole-doping effect on another orthogonal-dimer chain with linked-tetrahedra structure. It has been shown that the first-order metal-insulator transition is triggered by strong geometrical frustration. We have checked in a preliminary calculation that such a first-order metal-insulator transition also appears in the model of §2, implying that it may be common in this class of the orthogonal-dimer models.

## Acknowledgements

This work was partly supported by a Grant-in-Aid from the Ministry of Education, Science, Sports and Culture of Japan. A part of computations was done at the Supercomputer Center at the Institute for Solid State Physics, University of Tokyo and Yukawa Institute Computer Facility. A. Kawaguchi is supported by the Japan Society for the Promotion of Science.


text:
# Limiting temperatures and the Equation of State of Nuclear Matter

## I Acknowledgements

We appreciate helpful discussions with C. A. Gagliardi, J. Hardy, S. Shlomo and D. H. Youngblood. This work was supported by the U S Department of energy under Grant DE-FG03-93ER40773 and by the Robert A. Welch Foundation.


text:
# First results from dark matter search experiment with LiF bolometer at Kamioka Underground Laboratory

## 1 Introduction

Weakly Interacting Massive Particles (WIMPs) are thought to be the most promising dark matter candidates and many experiments have been performed for direct detection of WIMPs. The DAMA experiment reported having detected WIMPs by an annual modulation signal with their 10kg NaI detectors. The reported WIMP parameter region is $`M_{\mathrm{WIMP}}=(52_8^{+10})\mathrm{GeV}`$ and $`\sigma _{\mathrm{WIMP}\mathrm{p}}^{\mathrm{SI}}=(7.2_{0.9}^{+0.4})10^6\mathrm{pb}`$. On the other hand, the CDMS Collaboration and EDELWEISS Collaboration reported having excluded most of the parameter space claimed by the DAMA experiment within the standard supersymmetry(scalar couplings) and the halo models. These experiments aim to detect WIMPs via the spin-independent interactions. In 2001, the DAMA group published a paper in which they interpreted the annual modulation signal in a spin-independent and spin-dependent mixed framework. According to Ref. , $`\sigma _{\mathrm{WIMP}\mathrm{p}}^{\mathrm{SD}}`$ could be as high as 0.6pb for $`3080\mathrm{G}\mathrm{e}\mathrm{V}`$ if we assume certain sets of parameters. These results suggest that WIMPs could be detected via spin-dependent interaction and the dark matter search experiments should be performed both via spin-independent and spin-dependent interactions.

We, the Tokyo group, had developed a lithium fluoride (LiF) bolometer to detect WIMPs via spin-dependent interaction in the surface laboratory then performed pilot runs at the shallow depth site. We have improved our LiF bolometer since we installed the detector system in 2000 at Kamioka Observatory (2700m.w.e) and recently we made a dark matter search measurement from 2001 through 2002.

In this paper, we report our first results from the dark matter search experiment performed at Kamioka Observatory with new limits in the $`a_\mathrm{p}a_\mathrm{n}`$ (WIMPs-nucleon couplings) plane.

## 2 Experimental set-up and measurements

The laboratory for this dark matter search experiment was newly caved near the Super-Kamiokande detector at Mozumi Mine of Kamioka Mining and Smelting Co. located at Kamioka-cho, Gifu, Japan. The laboratory is at the depths of 2700 m.w.e and the muon flux is reduced about five orders of magnitude compared to the surface laboratory. As the radon concentration in the mine air is high and has a seasonal change (40 (winter)$``$1200 (summer) Bq/$`\mathrm{m}^3`$), the fresh air from the outside of the mine is blown into the laboratory with a rate of 1 $`\mathrm{m}^3/\mathrm{min}`$. The radon concentration in the laboratory is less than 50 Bq/$`\mathrm{m}^3`$.

The detector used in this dark matter search experiment is a LiF bolometer array consisting of eight LiF bolometers with a total mass of 168g. Each bolometer consists of a 2$`\times `$2$`\times `$2 cm3 LiF crystal and a neutron transmutation doped (NTD) Ge thermistor (1.5$`\times `$1$`\times `$1 mm3 ) glued to the LiF with GE varnish. We refer to each bolometer as D1, D2,…, and D8. Four of the LiF crystals are purchased from BICRON Co. Ltd. (expressed as BIC in Table 1) while the others are from OHYO KOKEN KOGYO Co. Ltd. (expressed as OKEN in Table 1). NTD3 thermistors, which have been used for most of the previous measurements, are used for six out of the eight bolometers, while NTD2 and NTD4 are used for one bolometer each. Characteristics of the NTD thermistors are shown in Ref. . Radioactive contamination in the crystals were previously checked with a low-background Ge spectrometer and was found to be less than 0.2 ppb for U, 2 ppb for Th, and 20 ppm for K. Surfaces of the LiF crystals were etched with perchloric acid in order to eliminate remaining grains of the polishing powder. Features of the eight bolometers are summarized in Table 1.

Schematic drawing of the bolometer array is shown in Fig. 1. The array consists of two stages, each stage holding four bolometers. The bolometers are fixed with spring pins and Delrin balls. Thermal contacts between the LiF crystals and the oxygen free high conductivity (OFHC) copper holder are achieved by annealed OFHC copper ribbons of 10 mm width and 0.1 mm thickness. The etching process and the assembly were performed in a clean bench of class 100. The bolometer array is encapsulated in a low-radioactivity lead (older than 200 years) of 2cm thickness. The OFHC copper holder and the inner shield were etched with nitric acid in order to eliminate the surface contamination.

The bolometer array and the inner shields are suspended with three Kevlar cords from the bottom of the dilution refrigerator which is cooled down to below 10mK. Thermal connection to the dilution refrigerator is achieved by four annealed OFHC copper ribbons of 10 mm width, 120 mm length, and 0.1 mm thickness. Each Kevlar cord is 8 cm long and eliminates the microphonic noise due to the helium liquefier. The dilution refrigerator is constructed with the low-radioactivity materials assayed with a low-background Ge spectrometer in advance.

The refrigerator is shielded with 10 cm of OFHC copper, 15cm of lead, $`1\mathrm{g}/\mathrm{cm}^2`$ of boric acid sheet, and 20 cm of polyethylene as shown in Fig. 1. Materials for the passive shields were previously radio-assayed with a low-background Ge spectrometer. An air tight plastic bag is set between the copper shield and the cryostat. The plastic bag is filled with nitrogen gas evaporated from liquid nitrogen. The radon concentration was less than 0.2 Bq/$`\mathrm{m}^3`$ throughout the measurement period, which makes negligible background to this dark matter search experiment.

Each thermistor is DC biased through a 300 M$`\mathrm{\Omega }`$ load resistor. The bias current of each bolometer is listed in Table 2. The voltage change across the thermistor is fed to the voltage sensitive amplifier placed at the 4K cold stage. A junction field effect transistor (J-FET), NEC 2SK163, is used for the cold stage amplifier with a gain of 10. The signal from the cold stage amplifier is in turn amplified by the voltage amplifier placed at the room temperature with a gain of 4.86 $`\times 10^2`$. The signal from the voltage amplifier undergoes through a low pass filter with a cut-off frequency of 11 Hz and is digitized with a sampling rate of 1 kHz. All the data are recorded continuously without any online trigger for a complete offline analysis. The data size for an hour measurement is 74M bytes.

The measurement was performed from November 22, 2001 through January 12, 2002. Eight of the eight bolometers worked, however, D1 and D5 showed poor signal to noise ratios because of the thermistors (NTD2 and NTD4). We analyzed the data obtained with the other six bolometers.

We applied an offline level trigger to the data averaged over $`\pm `$10 msec. Trigger levels are determined so that the trigger rates do not exceed 1 Hz. The trigger levels and the trigger rates are shown in Table 2. We determined the energy of each event by fitting the scaled template signal to the raw data. The electric noise events and the events that are thought to have occurred within the thermistors are eliminated by the wave form analysis. Then the events above the threshold energy that is determined as the energy with efficiency=0.8 are selected as real bolometer signals. The trigger and event selection efficiencies are estimated by adding the scaled template signals to the raw data as described in Ref. .Finally, the real bolometer events in coincidence with two or more detectors within $`\pm `$ 10 msec are discarded since WIMPs hardly interact in two or more detectors simultaneously. Energy resolutions of the bolometers are estimated by the same method as the efficiency estimation. The threshold energies used for the event selection and estimated energy resolution near the threshold are listed in Table 2.

The calibrations are carried out by irradiating $`{}_{}{}^{60}\mathrm{Co}`$ and $`{}_{}{}^{137}\mathrm{Cs}`$ gamma-ray sources suspended between the heavy shield and the cryogenics. Compton edges are used and good agreements with the simulations are seen. Good linearities are seen as shown in Ref. . Linearities down to 60keV are confirmed previously. Recorded pulse heights to a given signal for six detectors are shown as ’gain’ in Table 2. Gains are within the range of factor 7 which seems to be reasonable, considering that the assembly conditions are not exactly the same for the six bolometers. Stability of the detector response and linearities are checked every three or four days and the gains are found to be stable within $`\pm `$ 5$`\%`$.

## 3 Results

The spectra obtained with the six bolometers are shown in Fig. 2. One of the spectra (D3) obtained in the pilot run at Nokogiriyama underground cell is shown for comparison. Live times of the six bolometers are listed in Table 2. The total exposure is 4.1 kg$``$days. Background levels of D6, D7, and D8 have decreased by nearly one order of magnitude compared to the pilot run. D4 suffers from electronics noises, therefore the background level of D4 is higher than those of other detectors.

From the obtained spectra, we derived the spin-dependent WIMP-proton elastic scattering cross section ($`\sigma _{\mathrm{WIMP}\mathrm{p}}^{\mathrm{SD}}`$) limits. The limits are shown in Fig. 3.We calculated the limits in the same manner as used in Ref.  and Ref. . Here we used the astrophysical and nuclear parameters listed in Table 3. We used the earth velocity which corresponds to the velocity in the measurement period. We used the $`\lambda ^2J(J+1)`$ values calculated by assuming the odd group shell model. For the WIMPs heavier than 50 GeV, the limits have been improved by factor 5 from the pilot run. Since we determined the energy threshold (efficiency=0.8) higher than the pilot run, we did not improve the cross section limits for the lighter WIMPs compared to the pilot run. The best limit is 23pb for $`M_{\mathrm{WIMP}}=`$40GeV.

We also derived the spin-dependent WIMP-proton and WIMP-neutron elastic scattering cross section limits from $`{}_{}{}^{7}\mathrm{Li}`$ and $`{}_{}{}^{19}\mathrm{F}`$ independently in the same manner as described in Ref. , or using

$$\sigma _{\mathrm{WIMP}\mathrm{p}(\mathrm{N})}^{\mathrm{SDlim}}=\sigma _{\mathrm{WIMP}\mathrm{N}}^{\mathrm{SDlim}}\frac{\mu _{\mathrm{WIMP}\mathrm{p}}^2}{\mu _{\mathrm{WIMP}\mathrm{N}}^2}\left(\frac{C_{\mathrm{p}(\mathrm{N})}^{\mathrm{SD}}}{C_\mathrm{p}^{\mathrm{SD}}}\right)^1$$
(1)
and

$$\sigma _{\mathrm{WIMP}\mathrm{n}(\mathrm{N})}^{\mathrm{SDlim}}=\sigma _{\mathrm{WIMP}\mathrm{N}}^{\mathrm{SDlim}}\frac{\mu _{\mathrm{WIMP}\mathrm{n}}^2}{\mu _{\mathrm{WIMP}\mathrm{N}}^2}\left(\frac{C_{\mathrm{n}(\mathrm{N})}^{\mathrm{SD}}}{C_\mathrm{n}^{\mathrm{SD}}}\right)^1,$$
(2)
where $`\sigma _{\mathrm{WIMP}\mathrm{p}(\mathrm{N})}^{\mathrm{SDlim}}`$ and $`\sigma _{\mathrm{WIMP}\mathrm{n}(\mathrm{N})}^{\mathrm{SDlim}}`$ are the spin-dependent WIMP-proton and WIMP-neutron cross section limits obtained on the assumption that all events are due to WIMP-proton and WIMP-neutron elastic scatterings in the nucleus N, respectively, $`\sigma _{\mathrm{WIMP}\mathrm{N}}^{\mathrm{SDlim}}`$ is the WIMP-nucleus cross section limit which is directly set by the experiments, $`\mu _{\mathrm{WIMP}\mathrm{p}}`$ and $`\mu _{\mathrm{WIMP}\mathrm{n}}`$ are the WIMP-proton and WIMP-neutron reduced masses, respectively, $`\mu _{\mathrm{WIMP}\mathrm{N}}`$ is the WIMP-nucleus reduced mass, $`C_{\mathrm{p}(\mathrm{N})}^{\mathrm{SD}}`$ and $`C_{\mathrm{n}(\mathrm{N})}^{\mathrm{SD}}`$ are the proton and neutron contributions to the total enhancement factor of nucleus N, respectively, and $`C_\mathrm{p}^{\mathrm{SD}}`$ and $`C_\mathrm{n}^{\mathrm{SD}}`$ are the enhancement factors of proton and neutron themselves, respectively. Ratios of the enhancement factors calculated by assuming the odd group shell model are shown in Table 4.

The enhancement factor and the cross section are related by the following equation:

$`\sigma _{\mathrm{WIMP}\mathrm{N}}^{\mathrm{SD}}`$ $`=`$ $`4G_\mathrm{F}\mu _{\mathrm{WIMP}\mathrm{N}}^2C_\mathrm{N}^{\mathrm{SD}}`$ (3)
$`=`$ $`4G_\mathrm{F}^2\mu _{\mathrm{WIMP}\mathrm{N}}^2{\displaystyle \frac{8}{\pi }}(a_\mathrm{p}S_{\mathrm{p}(\mathrm{N})}+a_\mathrm{n}S_{\mathrm{n}(\mathrm{N})})^2{\displaystyle \frac{J+1}{J}},`$
where $`\sigma _{\mathrm{WIMP}\mathrm{N}}^{\mathrm{SD}}`$ is the WIMP-nucleus cross section, $`G_\mathrm{F}`$ is the Fermi coupling constant, $`C_\mathrm{N}^{\mathrm{SD}}`$ is total enhancement factor of nucleus N, $`a_\mathrm{p}`$ and $`a_\mathrm{n}`$ are the WIMP-proton and WIMP-neutron couplings, respectively, $`S_{\mathrm{p}(\mathrm{N})}`$ and $`S_{\mathrm{n}(\mathrm{N})}`$ are the expectation values of the proton and neutron spins within the nucleus N, respectively, and $`J`$ is the total nuclear spin. $`S_{\mathrm{p}(\mathrm{N})}`$, $`S_{\mathrm{n}(\mathrm{N})}`$, and $`J`$ can be calculated with shell models and are shown in Table 4.

From the obtained WIMP-proton and WIMP-neutron elastic scattering cross section limits, we derived the limits in the $`a_\mathrm{p}a_\mathrm{n}`$ plane. We again followed Ref. . The allowed region in the $`a_\mathrm{p}a_\mathrm{n}`$ plane is defined as follows:

$$\left(\frac{a_\mathrm{p}}{\sqrt{\sigma _{\mathrm{WIMP}\mathrm{p}(\mathrm{N})}^{\mathrm{SDlim}}}}\pm \frac{a_\mathrm{n}}{\sqrt{\sigma _{\mathrm{WIMP}\mathrm{n}(\mathrm{N})}^{\mathrm{SDlim}}}}\right)^2<\frac{\pi }{24G_\mathrm{F}^2\mu _{\mathrm{WIMP}\mathrm{p}}^2},$$
(4)
where the relative sign inside the square is determined by the sign of $`S_{\mathrm{p}(\mathrm{N})}/S_{\mathrm{n}(\mathrm{N})}`$.

Obtained limits in the $`a_\mathrm{p}a_\mathrm{n}`$ plane for WIMPs with mass of 50 GeV are shown in Fig.4 (a). The regions between the two parallel lines are the allowed regions from 7Li and 19F alone, respectively. It is seen that a single nuclide is not sufficient to set a finite limit in the $`a_\mathrm{p}a_\mathrm{n}`$ plane. Since we have 7Li and 19F as the targets, we can take the intersection of the two regions as the combined allowed region from the LiF detector. This is the merit of the detector consisting of more than two nuclides with sufficient spins. We show the allowed regions in the $`a_\mathrm{p}a_\mathrm{n}`$ plane obtained by this experiment for light WIMPs (Fig.4 (b)) and heavy WIMPs (Fig.4 (c)). The limits set by this experiment for $`M_{\mathrm{WIMP}}`$=30GeV are $`|a_\mathrm{p}|<32`$ and $`|a_\mathrm{p}|<133`$.

We compared our results with those of other experiments. As the DAMA group had not shown the results from 23Na and 127I independently, we were not able to compare our results with the spin-dependent annual modulation allowed regions which is shown in Fig. 3. Instead, we compared our results with those of the UKDMC experiments. Since UKDMC shows the limits from 23Na and 127I independently , we can calculate the limits in the $`a_\mathrm{p}a_\mathrm{n}`$ plane from 23Na and 127I and combine them without any assumptions. It should be noted that limits in the $`a_\mathrm{p}a_\mathrm{n}`$ plane from the NaI detector are largely affected by the sensitivity ratio of 23Na and 127I; the threshold of the detector and the quenching factor of 127I largely change the limits. Furthermore, accurate calculation of the annual modulation allowed region requires all of the obtained experimental data with time information, which we do not have. Therefore, we did not take the risk of deriving the limits or allowed region in the $`a_\mathrm{p}a_\mathrm{n}`$ plane from the experiments of which independent results from 23Na and 127I are not shown.

In Fig. 5, the allowed regions for the (a): sufficiently light WIMPs for which 127I loses the sensitivity and NaI effectively becomes a detector with one nuclide, (b): WIMPs with mass near the center value of DAMA’s annual modulation signal, and (c): WIMPs heavy enough for which 23Na and 127I both have good sensitivities, are shown. In the three plots, the limits from Na and I of UKDMC are shown in dotted and dash-dotted lines, respectively, NaI combined limits of UKDMC are shown in thin solid lines and LiF combined limits from this work are shown in thick solid lines. UKDMC published the recalculated data in which only NaI combined limits are shown, and we derived the limits in the $`a_\mathrm{p}a_\mathrm{n}`$ plane on the assumption that the sensitivity ratios of 23Na and 127I are the same as those shown in Ref. . These new limits from Ref.  are shown in dashed lines just for reference.

It is seen that our results can set limits complementary to the results obtained with NaI detectors. For the light WIMPs (lighter than 20 GeV for most of the detectors) to which 127I loses sensitivity, this is the only experiment that can set a finite allowed region in the $`a_\mathrm{p}a_\mathrm{n}`$ plane.

For the WIMPs with a mass of 50 GeV, which is near the center value of DAMA’s annual modulation signal, we exclude more than two thirds of the parameter space allowed by the UKDMC experiments. As the limits from the UKDMC experiments are close to the DAMA’s spin-dependent allowed region, as shown in Fig. 3, we should also be able to compare our limits in the $`a_\mathrm{p}a_\mathrm{n}`$ plane. We, however, have to wait until DAMA publishes relevant data to calculate the allowed region in the $`a_\mathrm{p}a_\mathrm{n}`$ plane.

For the WIMPs with a mass of 100 GeV, we only exclude a part of the UKDMC allowed region. However, it must be stressed that though our $`\sigma _{\mathrm{WIMP}\mathrm{p}}^{\mathrm{SD}}`$ limits are about two orders of magnitude worse than the limits of UKDMC as shown in Fig. 3 , we are able to tighten their limits in the $`a_\mathrm{p}a_\mathrm{n}`$ plane. This is because the spin factors of 19F are large and the $`S_{\mathrm{p}(\mathrm{N})}/S_{\mathrm{n}(\mathrm{N})}`$ sign of 19F is different from those of 7Li, 23Na, and 127I.

## 4 Discussions and prospects

Though we are able to tighten the limits from the UKDMC experiment in the $`a_\mathrm{p}a_\mathrm{n}`$ plane, the background is still high for the direct detection of WIMPs. We studied the background sources and found that potassium contamination on the surface of the LiF crystals, tritium produced by the neutron capture of $`{}_{}{}^{6}\mathrm{Li}`$ during the development in the surface laboratory, and the uranium and thorium contamination in the OFHC copper holder and in the inner lead shield dominate the background in the energy region between 10 keV and 90 keV. We are planning to etch the LiF crystals more deeply in order to eliminate the contamination on the surface of the LiF crystals. We also plan to use NaF crystals as a bolometer in order to eliminate the effects of tritium. It is also expected that NaF crystals give more stringent limits in the $`a_\mathrm{p}a_\mathrm{n}`$ plane since $`{}_{}{}^{19}\mathrm{F}`$ and $`{}_{}{}^{23}\mathrm{Na}`$ set limits orthogonal to each other in the $`a_\mathrm{p}a_\mathrm{n}`$ plane as shown in Fig. 4 and Fig. 5. We are planning to install active shields between the LiF crystals and the OFHC copper holders to reduce the background from the OFHC copper holder and the inner shield. With these improvements, we expect to improve the sensitivity of our detector by more than one order of magnitude.

## 5 Conclusions

We have performed the first underground dark matter experiment at Kamioka Observatory. With a total exposure of 4.1 kg$``$days, we have obtained $`\sigma _{\mathrm{WIMP}\mathrm{p}}^{\mathrm{SD}}`$ limits factor 5 improved from the pilot run at Nokogiriyama underground cell. We calculated the limits in the $`a_\mathrm{p}a_\mathrm{n}`$ plane and showed that our results tightened the limits from the UKDMC experiment. For the WIMPs with a mass of 50 GeV, we have excluded more than two thirds of the parameter space in the $`a_\mathrm{p}a_\mathrm{n}`$ plane allowed by the UKDMC experiment. As DAMA’s spin-dependent annual modulation allowed region in Fig. 3 is close to UKDMC’s limits, we should also be able to compare our limits in the $`a_\mathrm{p}a_\mathrm{n}`$ plane with DAMA’s spin-dependent annual modulation allowed region when DAMA publishes relevant data to calculate the allowed region in the $`a_\mathrm{p}a_\mathrm{n}`$ plane.

## Acknowledgements

We would like to thank all the staffs of Kamioka Observatory of the Institute for Cosmic Ray Research, the University of Tokyo, to whose hospitality we owe a great deal in using the facilities of the Observatory. We would also like to express our gratitude to Prof. Komura at Kanazawa University for providing us with 25kg of low-activity old lead. This research is supported by the Grant-in-Aid for COE Research by the Japanese Ministry of Education, Culture, Sports, Science and Technology.