text:
# A globally accurate theory for a class of binary mixture models

## Abstract

Using the self-consistent Ornstein-Zernike approximation (SCOZA) results for the 3D Ising model, we obtain phase diagrams for binary mixtures described by decorated models. We obtain the plait point, binodals, and closed-loop coexistence curves for the models proposed by Widom, Clark, Neece, and Wheeler. The results are in good agreement with series expansions and experiments.

aelectronic address: adickman@pucminas.br
phone number: 55-31-3319-4176
belectronic address: gstell@mail.chem.sunysb.edu
phone number: +1-631-632-7898
fax number: +1-631-632-7960

1. Introduction

Critical phenomena of mixtures are considerably more complicated to investigate than such phenomena in pure fluids. Purification of components, precise determination of composition, gravitational effects and long equilibration times make experiments on mixtures much more difficult to realize. Similar difficulties are found in theoretical analyses, where the addition of one component greatly complicates the generalization of the lattice-gas model to mixtures, for example the derivation and analysis of series expansions .

The use of decorated lattice-gas models to describe fluid mixtures represents an important gain in precision and information about these models. This is a result of the exact mapping between these models and the single-species lattice-gas or Ising model. The exact results in two dimensions, and exhaustive series analysis in three, allow a detailed description of the phase diagrams, coexistence surfaces and critical exponents of the much more complicated models for mixtures.

In this paper we obtain information related to critical phenomena of binary mixtures for various decorated lattice models using Ising-model results obtained via the SCOZA . This approximation provides very precise results for the lattice-gas Ising model in three dimensions as shown elsewhere .

The next sections define the models we work with and show how to derive the grand partition function for these models. Then we present a brief discussion about SCOZA. Finally, we present the phase diagrams and coexistence curves we obtain for the binary mixture models.

2. Decorated model

Consider a binary mixture of $`N_1`$ molecules of type 1 and $`N_2`$ molecules of type 2 contained in a volume $`V`$. A decorated lattice model for a mixture is constructed as follows . We divide the volume into C primary cells, which are centered on sites of a regular lattice, and secondary cells centered on the bond midpoints of this lattice. If the primary cells have volume $`v_0`$, then the secondary cells occupy a volume $`2v_0/q`$, where $`q`$ is the lattice coordination number. Each primary cell overlaps $`q`$ secondary cells and each secondary cell overlaps two primary cells, as shown in figure 1. A primary cell is occupied if a type-1 molecule lies within its boundaries; secondary cells are occupied by type-2 molecules. The cells are not meant to restrict the molecules to lattice sites. They merely provide a discrete coordinate system for measuring the positions of particles of type 1 and 2 and thus defining intermolecular separations.

\[Insert figure 1 about here\]

All interactions are pairwise. If two molecules of type 1 are in the same primary cell, the interaction energy is $`+\mathrm{}`$; this energy is $`ϵ`$ if they occupy adjacent primary cells, and zero otherwise. The interaction energy for a pair of type-2 molecules is $`+\mathrm{}`$ if they occupy the same secondary cell and zero otherwise. If a primary cell containing a type-1 molecule overlaps a secondary cell occupied by a type-2 molecule, the interaction energy is $`\varphi `$. Particles of different species have no interaction if their cells do not overlap.

The condition that there is no interaction between type-2 molecules in different secondary cells is essential to the existence of an exact mapping between the decorated model and the spin-$`\frac{1}{2}`$ Ising model. As a consequence, the pure component 2 cannot undergo a phase transition .

So far we have described the decorated model in general terms. The specific models for binary mixtures by Widom, Clark, Neece, and Wheeler are obtained by adjusting the parameters for the interaction energy.

Widom  defined the simplest model by setting a hard-core repulsion between molecules of types 1 and 2, the only interaction beyond exclusion of multiple occupancy. The model undergoes a first-order phase transition even though the intermolecular forces are all assumed to be infinitely strong repulsions.

Clark  generalized Widom’s model by allowing a finite repulsion between unlike nearest-neighbour molecules. All other interactions between like molecules are zero and multiple occupancy of cells is not allowed. This model has a coexistence dome with a maximum critical solution temperature.

Neece  extended Widom’s model including the effects of interactions between molecules on the primary lattice. This model exhibits a variety of behaviours with changing temperature, especially if the molecules on the primary lattice interact with an attractive potential. It provides a good qualitative description of certain types of ‘gas-gas equilibrium’ such as that found in the He-Xe system.

Wheeler  studied the binary case with an attractive interaction between molecules on the primary lattice, and a finite interaction (attractive or repulsive) between type-1 and -2 molecules. This version proved useful in understanding critical phenomena in very dilute solutions.

Another model studied by Wheeler  considers a highly directional, short-ranged interaction favoring special orientations of the molecules at low temperatures. Unlike the other models, a primary or secondary lattice site can be occupied either by a type-1 or a type-2 molecule. Interactions between particles of the same species (1-1 or 2-2) are zero regardless of their orientations. The 1-2 interaction depends upon the orientation of the particle in the secondary cell. If the 1-2 bond involves the special contact point of this particle then the energy is negative (attractive) otherwise it is positive (repulsive). Each molecule (1 or 2) has one special contact point and $`(\omega 1)`$ other identical contact points regardless of whether it is on a primary or secondary cell. We consider the case with $`\omega =6`$. Thus, if a molecule of either type is in a secondary cell which has one 1-filled and one 2-filled primary neighbour, then it has one orientation with energy $`U_2<0`$ and five orientations with energy $`U_1>0`$. If both their primary neighbours are of the opposite type, then there are two orientations in which it makes one attractive and one repulsive contact with resulting energy $`U_2+U_1`$, and four orientations in which it makes two repulsive contacts with resulting energy $`2U_1`$. The 1-2 interaction energy is independent of the orientation of the particle in the primary cell. This is required to map the model exactly onto the spin-1/2 Ising model. This model presents closed-loop coexistence curves with both upper and lower critical solution temperatures, analogous to the behaviour found in the nicotine + water and m-toluidine + glycerol systems.

Table 1 shows a summary of the models presented above.

\[Insert table 1 about here\]

3. Derivation of the Grand Partition Function

In this section we derive the grand partition function for the decorated models and show how to map these models into the lattice gas Ising model. We follow the derivation given by Wheeler (for details see ). The grand partition function can be written as

$`𝒵={\displaystyle \underset{𝒞_1}{}}{\displaystyle \underset{𝒞_2}{}}z_1^{N_1}x^{N_{11}}z_2^{N_2}e^{\frac{E^{(2)}}{kT}}`$ (1)
where the sums are over all assignments of filled or empty primary and secondary cells. $`E^{(2)}`$ is the sum of all primary-secondary interactions, $`N_{11}`$ is the number of filled-filled nearest neighbour pairs of primary cells and $`x=\mathrm{exp}(ϵ/kT)`$. $`z_1`$ and $`z_2`$ are dimensionless activities defined as,

$`z_1`$ $`=`$ $`v_0\left({\displaystyle \frac{2\pi m_1kT}{h^2}}\right)^{\frac{d}{2}}e^{\frac{\mu _1}{kT}},`$ (2)
$`z_2`$ $`=`$ $`\left({\displaystyle \frac{qv_0}{2}}\right)\left({\displaystyle \frac{2\pi m_2kT}{h^2}}\right)^{\frac{d}{2}}e^{\frac{\mu _2}{kT}},`$ (3)
where $`d`$ is the dimensionality.

Let $`\eta _{ij}=e^{\psi _{ij}/kT},i,j=0,1`$, where $`\psi _{00},\psi _{01},\psi _{11}`$ are the energies for a molecule of type 2 that lies between two empty, one empty and one filled, or two filled primary cells, respectively. The idea is to sum over filled and empty states for each secondary cell, once a specific configuration of the primary cells is chosen. Thus if a secondary cell lies between empty-empty primary cell pairs, it contributes to the grand partition function a factor $`1+z_2\eta _{00}`$, if it lies between empty-filled primary cell pairs, it contributes a factor $`1+z_2\eta _{01}`$, and if it lies between filled-filled primary cell pairs, it contributes a factor $`1+z_2\eta _{11}`$. Then,

$`𝒵={\displaystyle \underset{𝒞_1}{}}z_1^{N1}x^{N_{11}}(1+z_2\eta _{00})^{N_{00}}(1+z_2\eta _{01})^{N_{01}}(1+z_2\eta _{11})^{N_{11}}`$ (4)
where $`N_{00}`$ and $`N_{01}`$ are the number of empty-empty and empty-filled nearest neighbour primary cells. Using the lattice identities, $`qN_1=2N_{11}+N_{01},qN_0=2N_{00}+N_{01},N_1N`$ and $`N_0=CN`$ we get

$`𝒵=(1+z_2\eta _{00})^{\frac{qC}{2}}{\displaystyle \underset{𝒞_1}{}}z^N\mathrm{exp}\left({\displaystyle \frac{ϵN_{11}}{kT}}\right)=(1+z_2\eta _{00})^{\frac{qC}{2}}𝒵_R.`$ (5)
Comparing this with the partition function for the lattice gas,

$`𝒵_R={\displaystyle \underset{𝒞}{}}(z_R)^N\mathrm{exp}(\beta _Rϵ_RN_{11}),`$
we obtain the relations that map the decorated model exactly to the lattice gas Ising model:

$`𝒵`$ $`=`$ $`(1+z_2\eta _{00})^{\frac{qC}{2}}𝒵_R`$
$`z_R`$ $`=`$ $`z_1(1+z_2\eta _{00})^q(1+z_2\eta _{01})^q`$
$`\mathrm{exp}\left({\displaystyle \frac{ϵ_R}{kT_R}}\right)`$ $`=`$ $`x(1+z_2\eta _{00})(1+z_2\eta _{11})(1+z_2\eta _{01})^2.`$ (6)
The subscript $`R`$ denotes that the quantity refers to the reference system. This is an exact correspondence between the decorated model and the spin-$`1/2`$ Ising model. All thermodynamic properties can be obtained from the grand canonical potential $`\pi (z_1,z_2,x)=lim_c\mathrm{}1/C\mathrm{ln}𝒵`$ by differentiation.

The thermodynamic pressure $`p(\mu ,T)`$ is given by

$`\beta p=\underset{V\mathrm{}}{lim}{\displaystyle \frac{1}{V}}\mathrm{ln}𝒵(T,\mu ,V)={\displaystyle \frac{q}{2}}\mathrm{ln}(1+z_2\eta _{00})+\beta _Rp_R`$ (7)
for the decorated model. The number densities are given by:

$`\rho _1`$ $`=`$ $`\left({\displaystyle \frac{\pi }{\mathrm{ln}z_1}}\right)_{z_2,x}=\rho _R`$
$`\rho _2`$ $`=`$ $`\left({\displaystyle \frac{\pi }{\mathrm{ln}z_2}}\right)_{z_1,x}={\displaystyle \frac{q}{2}}A+q(CA)\rho _R+(A+B2C){\displaystyle \frac{\mathrm{\Phi }_R}{ϵ_R}}`$ (8)
where $`\mathrm{\Phi }_R`$ is the energy density in the reference system and,

$`A={\displaystyle \frac{z_2\eta _{00}}{(1+z_2\eta _{00})}}B={\displaystyle \frac{z_2\eta _{11}}{(1+z_2\eta _{11})}}C={\displaystyle \frac{z_2\eta _{01}}{(1+z_2\eta _{01})}}.`$ (9)
Other quantities of interest are the internal energy, given by $`U=(\pi /\beta )_{z_1,z_2}`$, and the specific heat at constant volume, given by $`C_{v,x}=(U/T)_{\rho _1,\rho _2}`$.

4. Self-Consistent Ornstein-Zernike Approximation(SCOZA)

The SCOZA closure of the Ornstein-Zernike equation yields critical parameters of remarkable accuracy for three-dimensional lattice gases with nearest-neighbour interactions .

In the SCOZA, the direct correlation function $`c(r)`$ is determined in a way that results in a unique free energy function. The starting point is the Ornstein-Zernike relation

$`h_i=c_i+\rho {\displaystyle \underset{Lattice}{}}c_{ij}h_j`$ (10)
where $`h`$ is the total correlation function. Thermodynamic consistency is embodied in the exact relation that must be satisfied by the derivation of the free energy

$`\rho {\displaystyle \frac{^2(\rho u)}{\rho ^2}}={\displaystyle \frac{^2(\beta p)}{\beta \rho }}`$ (11)
where $`u`$ is the internal energy per particle, $`p`$ is the pressure and $`\rho `$ the density. $`c_0c(0)`$ is fixed by the core condition on the total correlation function, $`h_0=1`$. The single approximation of the method is the truncation of the direct correlation function $`c(r)`$ at the nearest-neighbour separation: $`c(r)=0`$ for $`r>1`$. From the relations between $`c_0,c_1`$ and $`h_1`$ we find a partial differential equation for $`c_1(\beta ,\rho )`$ that can be solved numerically. For details see references .

5. Results

In this section we specify the mapping equations for each model of interest and we also show the corresponding phase diagrams. We study binary mixtures in three dimensions. Thus the reference one-component lattice gas has a range of thermodynamic states corresponding to liquid-vapor coexistence terminating in a critical point. At the boundary of the two-phase region in the ($`\mathrm{\Phi }/\mathrm{\Omega },N/\mathrm{\Omega }`$) plane of the reference system the activity $`z_R`$, the temperature $`T_R`$, and the pressure-temperature ratio $`p_R/T_R`$ are all singular, having discontinuous derivatives there. Therefore, we see from equations ( 6) and (7), that the functions $`z_1,z_2,p/T`$ of the two-component models are singular at the corresponding point in the ($`\rho _1,\rho _2`$) plane.

As a reminder, in our notation $`ϵ`$ is the interaction energy between two type 1 molecules when sitting at neighbouring primary cells, and $`\varphi `$ denotes the interaction energy between a type-1 and a type-2 molecule when they occupy overlapping primary and secondary cells.

5.1. Reference system

Our reference system is the lattice gas transcription of the Ising model, which has two-phase behaviour below a critical temperature $`\theta _c`$. The coexistence curve is shown in figure 2. In the two-phase region there are two values of the density which are in equilibrium for each temperature $`\theta `$. The same is true for the binary system after carrying out the transformations mapping it to the reference system. For each temperature below the critical temperature $`T_c`$ in the binary system, the binodal curve corresponds to the coexistence curve in the reference system. By varying the temperature the binodals will then form the coexistence surface.

\[Insert figure 2 about here\]

Before investigating the nature of the coexistence surface it is important to review some of the properties of the one-component lattice gas in the two-phase region. Nearest neighbour lattice gases of dimensionality greater than one have a critical density $`\rho _c=1/2`$, with the critical isochore corresponding to this density. The properties along the critical isochore are known exactly only in the two-dimensional case. However, the behaviour in three dimensional systems is known with sufficient accuracy for our purposes. In this paper we use the results obtained via SCOZA.

The equation of the coexistence curve can be written as $`\theta =\tau (\rho )`$, for which there are two values of $`\rho `$ for each $`\theta `$ below the critical temperature $`\theta _c`$. The activity along the critical isochore is

$`\lambda =\mathrm{exp}\left({\displaystyle \frac{Zϵ}{2k\theta }}\right).`$ (12)
As $`\theta \theta _c`$ in the two-phase region, the usual power law expressions are assumed to hold. The coexistence curve in the neighbourhood of the critical density can be written as

$`\theta _c\tau (\rho )|\rho \rho _c|^{\frac{1}{\beta }},`$ (13)
where the exponent $`1/\beta `$ is very close to 3 in three dimensions, so the coexistence curve is nearly cubic. In SCOZA, $`1/\beta =20/7`$ . The heat capacity (constant volume) diverges as the critical point is approached asymptotically along the critical isochore, following

$`\mathrm{\Gamma }(\theta _c\theta )^\alpha ^{},`$ (14)
where $`\alpha ^{}>0`$ and nearly equal to zero in three dimensions. In SCOZA, $`\alpha ^{}=1/10`$ . Current best estimates for $`\beta `$ and $`\alpha ^{}`$ are $`0.3257\pm 0.0025`$ and $`0.1130\pm 0.0075`$ , respectively.

5.2. Behaviour of a general binary model

The complete three-dimensional phase diagram for a binary system at constant pressure is shown in figure 3.

\[Insert figure 3 about here\]

The coexistence surface separates the one- and two-phase regions of the system. In a constant temperature plane the variables are the number densities of the two components. The two-phase region is bounded by a binodal curve APB, containing tie lines between points representing the compositions of the two coexisting phases. The tie lines vanish at point P, called the plait point (isothermal critical mixing point), where the coexistent phases lose their separate identities to merge into a single homogeneous phase. In the full three-dimensional diagram the plait points form a curve on the coexistence surface. This curve ends up reaching a maximum at the top C, above which only one phase is found. C is then a critical point of the system. We will see that in some cases the binodals are closed loops involving two plait points.

5.3. Widom’s model

The number densities for Widom’s model, which is defined by $`ϵ=0`$ and $`\varphi =+\mathrm{}`$, are

$`\rho _1`$ $`=`$ $`\rho _R`$
$`\rho _2`$ $`=`$ $`{\displaystyle \frac{1}{3}}\left({\displaystyle \frac{q}{2}}q\rho _R+{\displaystyle \frac{\mathrm{\Phi }_R}{ϵ_R}}\right)[\mathrm{exp}(\beta _Rϵ_R)1].`$ (15)
In order to obtain the binodal curve, shown in figure 4, we use the reference system coexistence curve values for the energy density, temperature and number density in equation 15.

\[Insert figure 4 about here\]

The binodal curve APB and the line AB for $`T=0`$ separates the two-phase region from the remainder of the composition triangle corresponding to the homogeneous states of the system. In order to determine the tie line in the binary system we identify the double-valued density in the reference system for a given temperature. The corresponding densities in the binodal are coexistent. The tie line CD connects pairs of points on the binodal curve which represent the compositions of coexistent phases. The longest tie line connects the coexistent densities ($`\rho _1=0,\rho _2=1`$) and ($`\rho _1=1,\rho _2=0`$), corresponding to $`\theta =0`$ in the reference system. The tie lines vanish at the point P, where the conjugate phases merge into a single homogeneous phase. Thus P is a critical point achieved by changing composition at fixed temperature. P is frequently called a plait point in order to distinguish it from the usual critical point achieved by changing temperature. The coordinates $`\rho _1^{},\rho _2^{}`$ of the plait point are

$`\rho _1^{}={\displaystyle \frac{1}{2}},\rho _2^{}={\displaystyle \frac{1}{3}}\left({\displaystyle \frac{\mathrm{\Phi }_c}{ϵ_R}}\right)[\mathrm{exp}(\beta _cϵ_R)1]=0.1968242`$ (16)
The binodal curve is an isothermal section of a more general three-dimensional diagram. For this model the number density and all the other quantities derived are independent of the temperature $`T`$. This happens because the molecular interaction potential is either $`0`$ or $`+\mathrm{}`$, so the total potential energy of this system never deviates from $`0`$, i.e. there is no energy scale. Nevertheless, the model exhibits phase coexistence and an associated plait point.

5.4. Clark’s model

For this model, $`ϵ=0`$ and $`\varphi `$ is finite. The number densities $`\rho _1`$ and $`\rho _2`$ are

$`\rho _1`$ $`=`$ $`\rho _R`$
$`\rho _2`$ $`=`$ $`{\displaystyle \frac{1}{3}}\left[{\displaystyle \frac{q}{2}}A+N{\displaystyle \frac{q}{2}}(CA)+N_{11}(A+B2C)\right]`$ (17)
where $`A,B`$ and $`C`$ were defined by equation 9 when deriving the grand partition function.

Equations 6 and 17 are the mapping between Clark’s model and the lattice gas reference system. In order to obtain the binodal curve, we determine the activities and the number densities for the binary model by using the energy density, number density and temperature values of the reference system coexistence curve.

\[Insert figure 5 about here\]

Figure 5 shows the isothermal binodal curves in the $`(\rho _1,\rho _2)`$ plane, corresponding to $`\varphi =0.5`$ and temperatures $`T=0.1;0.2;0.247162`$. The coexistence surface, figure 3, is the set of all such binodal curves, which close to form a dome whose apex marks the critical point ($`T_c=0.247173`$). Each isotherm in the $`(\rho _1,\rho _2)`$ plane, has an upper and a lower plait point. These plait points form a plait point curve that ends at the critical point $`T_c`$.

5.5. Neece’s model

For this model $`ϵ`$ is nonzero and $`\varphi =+\mathrm{}`$. The number densities are

$`\rho _1`$ $`=`$ $`\rho _R`$
$`\rho _2`$ $`=`$ $`{\displaystyle \frac{1}{3}}\left({\displaystyle \frac{q}{2}}+Nq+N_{11}\right)\left[1\mathrm{exp}\left({\displaystyle \frac{ϵ_R}{k\theta }}{\displaystyle \frac{ϵ}{kT}}\right)\right].`$ (18)
If $`ϵ<0`$ there is an attraction between nearest neighbours on the primary lattice. In our model $`|ϵ_R|=|ϵ|=1`$. The second term in brackets in the expression for $`\rho _2`$ is never zero for $`T>\theta _c`$. If $`T=\theta _c`$, then this factor is zero and $`\rho _2=0`$ at $`\rho _1=1/2`$. If $`T<\theta _c`$, the second factor is zero for some value of $`\theta `$ which in turn corresponds to two values of $`\rho _1`$. The binodal curve for $`ϵ<0`$ is shown in figure 6 at different temperatures. The critical temperature is $`T_c=1.1299`$, represented by curve C. From the inset of figure 6, which is a zoom of curve D, we can see that there are two separate regions corresponding to the homogeneous states of the system for $`T<\theta _c`$. Also, in this case there is no plait point since there is no tie line of vanishing length.

\[Insert figure 6 about here\]

For $`ϵ>0`$, the binodals as a function of temperature are shown in figure 7. The plait points lie along $`\rho _1=1/2`$.

\[Insert figure 7 about here\]

5.6. Wheeler’s model

For this model, $`ϵ<0`$ and $`\varphi `$ is finite (attractive or repulsive). The number densities $`\rho _1`$ and $`\rho _2`$ are given by

$`\rho _1`$ $`=`$ $`\rho _R`$
$`\rho _2`$ $`=`$ $`({\displaystyle \frac{q}{2}})A+q(CA)\rho _R+(A+B2C){\displaystyle \frac{\mathrm{\Phi }_{11}}{ϵ_R}}`$ (19)
where $`A,B`$ and $`C`$ are defined by equation 9.

\[Insert figure 8 about here\]

Figure 8 shows the binodal curves for $`\mathrm{\Phi }<0`$ and $`\mathrm{\Phi }>0`$ at different temperatures. Curve A is for $`T=1.11`$, curve B is for $`T=1.25`$ and curve C for $`T=1.336`$ for attractive and repulsive interactions. The slope of the tie lines is positive for $`\mathrm{\Phi }<0`$ and negative for $`\mathrm{\Phi }>0`$. At $`T=T_c=1.336498`$ the binodal shrinks to a point. In a 3-D diagram, including the temperature, it corresponds to a dome that separates the two-phase region from the homogeneous one.

5.7. Wheeler’s oriented model

In this model particles of the same species (1-1 or 2-2) do not interact. The 1-2 interaction depends upon the orientation of the particle in the secondary cell. If the 1-2 bond involves the special contact point of this particle then the energy is $`U_2<0`$ (attractive), otherwise it is $`U_1>0`$ (repulsive).

\[Insert figure 9 about here\]

Figure 9 shows the coexistence curve for $`r=|U_2|/U_1=0.189`$ (curve A), $`0.18`$ (curve B), and $`0.17`$ (curve C). Each curve forms a closed loop with a lower and upper critical temperature. The ratios $`T_u/T_l`$ of upper to lower critical solution temperature are respectively 1.161, 1.467 and 1.721 for these curves. Our results are in close agreement with those of Wheeler , using low temperature series expansions and with the experimental data.

In figure 10 we show the experimental coexistence curves for solutions of glycerol and guaiacol (circles)  and m-toluidine (squares) . These curves are qualitatively similar to the closed loops obtained from Wheeler’s oriented model.

\[Insert figure 10 about here\]

6. Concluding remarks

We study binary systems through various versions of the decorated lattice model using the self-consistent Ornstein-Zernike approximation. We obtain the binodals, plait points, and the critical temperatures for the models introduced by Widom, Neece, Clark and Wheeler. SCOZA provides a very useful way of getting accurate results for binary systems that can be mapped to the lattice-gas Ising model.

We also study Wheeler’s oriented model, obtaining its coexistence curve for various interaction energies. The results are consistent with experimental data and low temperature series expansions.

In order to analyse the algebraic degree of the binodal at the plait point we can use the idea proposed by Widom . It consists in measuring the vertical distance $`\mathrm{\Delta }\rho _2`$ from the binodal curve to a line drawn tangent to the binodal at the plait point. By combining the number densities with the energy density and the coordinates of the plait point, we find that

$`\mathrm{\Delta }\rho _2|\rho _1{\displaystyle \frac{1}{2}}|^{\frac{(1\alpha ^{})}{\beta }}.`$ (20)
The modification of the critical exponent by the factor $`(1\alpha ^{})`$ is an example of ‘renormalization’ discussed by Fisher . The SCOZA values for $`\beta `$ and $`\alpha ^{}`$ are 7/20 and 1/10, respectively, providing 18/7 for the algebraic degree of the binodal.

These results can be extended to study ternary systems as well as pure fluids. It is not difficult to derive the energy, heat capacity, pressure for all the models studied in this paper.

A.G.D. is grateful to CAPES and Fapemig (Brazil) for financial support. A.G.D. is also grateful to Ann and the staff of SUNY at Stony Brook for their hospitality and kindness during her stay. G.S. gratefully acknowledges the support of the Division of Chemical Sciences, Office of Basic Energy Sciences, Office of Energy Research, US Department of Energy. The authors thank F. Raineri for his assistance.

Figure Captions

Figure 1. Example of a decorated lattice in two dimensions. The solid squares represent the primary cells while the dashed diamonds represent the secondary cells.

Figure 2. Coexistence curve for the lattice gas Ising model using SCOZA results.

Figure 3. Coexistence surface for a binary mixture using the results for Clark’s model. The top C marks a critical point.

Figure 4. Binodal curve for Widom’s model. CD is a tie line and P is the plait point. The composition variables $`\rho _1`$ and $`\rho _2`$ are the number densities of the two components.

Figure 5. Binodal curve for Clark’s model for attractive ($`\mathrm{\Phi }>0`$) and repulsive ($`\mathrm{\Phi }<0`$) interactions. The outer binodal corresponds to $`T=0.1`$; the intermediate one corresponds to $`T=0.2`$ and the inner curve is for $`T=0.247`$, very close to $`T_c=0.247173`$.

Figure 6. Binodal curve for Neece’s model with $`ϵ<0`$. Curve A corresponds to $`T\mathrm{}`$; curve B is for $`T=2.0>T_c`$; the critical temperature $`T_c=1.12989`$ is represented by curve C; and $`T=1.0<T_c`$ at curve D. The inset is a zoom of curve D.

Figure 7. Binodal curve for Neece’s model with $`ϵ>0`$. Curve A corresponds to $`T=0`$; curve B is for $`T>0`$; $`T\mathrm{}`$ is represented by curve C.

Figure 8. Binodal curve for Wheeler’s model for repulsive ($`\mathrm{\Phi }<0`$) and attractive ($`\mathrm{\Phi }>0`$) interactions. Curve A is for $`T=1.11`$, curve B is for $`T=1.25`$ and curve C for $`T=1.336`$ for both diagrams. The critical temperature is $`T=T_c=1.336498`$.

Figure 9. Coexistence curves for Wheeler’s oriented model in the temperature-mole fraction diagram. Curves A, B and C correspond to $`r=0.189,0.18,0.17`$, respectively.

Figure 10. Temperature-mass fraction coexistence curves for solutions of glycerol with guaiacol (circles) and m-toluidine (squares). Temperature is in Celsius.


format: text

text:
# Hypercentral constituent quark model and isospin dependence

## 1 Introduction

Constituent Quark Models have been recently widely used for the description of the internal structure of baryons . The baryon spectrum is usually described well, although the various models are quite different. However the study of hadron spectroscopy is not sufficient to distinguish among the various forms of quark dynamics. To this end one has to study in a consistent way all the physical observables of interest, in particular, besides the spectrum, the photocouplings, the electromagnetic form factors and the strong decay amplitudes. Such a systematic study of baryon properties is better performed within a general framework, and in this respect a hypercentral approach to quark dynamics can be used . The model consists of a hypercentral quark interaction containing a linear plus coulomb-like term, as suggested by lattice QCD calculations . A hyperfine term of the standard form  is added and treated as a perturbation. The few free parameters of the model are fitted to the spectrum, the resulting baryon states are then used in order to calculate the various properties of interest, in particular the photocouplings , the transition form factors  and the elastic electromagnetic form factors of the nucleon . The electromagnetic properties are evaluated using a non relativistic current for pointlike quarks, also taking into account the effects of relativistic corrections . In particular this parameter-free calculation predicts that the ratio of the electric and magnetic proton factors decreases with $`Q^2`$ , as shown by the recent TJNAF experiment .

The description of the non strange baryon spectrum obtained by the hypercentral Constituent Quark Model (hCQM)  is fairly good and comparable to the results of other approaches. In particular, the SU(6)-structure of the levels is accounted for thanks to the spin-independent hypercentral interaction; the $`\mathrm{\Delta }N`$ mass difference is correctly described by the hyperfine splitting and the theoretical energies of the negative parity resonances are in good agreement with data. However, notwithstanding such overall fair description of the spectrum, in some cases the splittings within the various SU(6)- multiplets are too low. This is particularly true for the Roper resonances and for the higher states. A possible origin of these problems could be the (widespread) use of a $`\delta `$like hyperfine interaction. To this end we have introduced different kinds of space smearings; the resulting hyperfine term becomes acceptable from the theoretical point of view, but, as we shall show below, it does not improve the description of the spectrum.

A more important issue is the flavour dependence of the quark interaction. Actually, within the algebraic approach, the quark energy is written in terms of Casimir operators of symmetry groups which are relevant for the three-quark dynamics; in this respect it is quite natural to introduce an isospin dependent term, which turns out to be important for the description of the spectrum . On the other hand, in the chiral constituent quark model recently proposed , the splittings are produced by the one-boson-exchange interaction between quarks and therefore a flavour-dependent potential arises, which seems to be important in order to describe the baryon spectrum, at least below $`1.7GeV`$.

In the following, we shall show that in the hCQM a flavour dependent potential can be introduced  as a perturbative term leading to improved splittings within the $`SU(6)`$-multiplets. In particular, in this way, the Roper resonance is reproduced quite well and the higher states acquire a much larger splitting, in agreement with data.

In Section $`2`$ we remind briefly the model and the main results in the description of the spectrum and the electromagnetic excitation of the baryon resonances. In Section $`3`$ we introduce in the hCQM a generalized SU(6)-breaking interaction treated as a perturbation and we show the results of the model compared with the experimental spectrum. Finally, in Section $`4`$ there are some discussions and conclusions.

## 2 The hypercentral model

The internal quark motion is described by the Jacobi coordinates $`\rho `$ and $`\lambda `$:

$`\rho ={\displaystyle \frac{1}{\sqrt{2}}}(\text{r}_1\text{r}_2),\lambda ={\displaystyle \frac{1}{\sqrt{6}}}(\text{r}_1+\text{r}_22\text{r}_3),`$ (1)
or equivalently, $`\rho `$, $`\mathrm{\Omega }_\rho `$, $`\lambda `$, $`\mathrm{\Omega }_\lambda `$. In order to describe the three-quark dynamics it is convenient to introduce the hyperspherical coordinates, which are obtained substituting the absolute values $`\rho `$ and $`\lambda `$ by

$$x=\sqrt{\rho ^2+\lambda ^2},\xi =arctg(\frac{\rho }{\lambda }),$$
(2)
where $`x`$ is the hyperradius and $`\xi `$ the hyperangle. In this way one can use the hyperspherical harmonic formalism .

In the hypercentral constituent quark model (hCQM), the quark potential, $`V`$, is assumed to depend on the hyperradius $`x`$ only, that is to be hypercentral. Therefore, $`V=V(x)`$ is in general a three-body potential, since the hyperradius $`x`$ depends on the coordinates of all the three quarks. Since the potential depends on $`x`$ only, in the three-quark wave function one can factor out the hyperangular part, which is given by the known hyperspherical harmonics . The remaining hyperradial part of the wave function is determined by the hypercentral Schrödinger equation:

$$[\frac{d^2}{dx^2}+\frac{5}{x}\frac{d}{dx}\frac{\gamma (\gamma +4)}{x^2}]\psi _{[\gamma ]}(x)=2m[EV(x)]\psi _{[\gamma ]}(x),$$
(3)
where $`\psi _{[\gamma ]}(x)`$ is the hypercentral wave function and $`\gamma `$ is the grand angular quantum number given by $`\gamma =2v+l_\rho +l_\lambda `$; $`l_\rho `$ and $`l_\lambda `$ are the angular momenta associated with the $`\rho `$ and $`\lambda `$ variables and $`v`$ is a non negative integer number.

There are at least two hypercentral potentials which lead to analytical solutions. First, the h.o. potential, which has a two-body character, turns out to be exactly hypercentral, since

$$\underset{i<j}{}\frac{1}{2}k(\text{r}_i\text{r}_j)^2=\frac{3}{2}kx^2=V_{h.o.}(x).$$
(4)
The second one is the ’hypercoulomb’ potential

$$V_{hyc}(x)=\frac{\tau }{x}.$$
(5)
This potential is not confining, however it has interesting properties. In fact it leads to a power-law behaviour of the proton form factor  and of all the transition form factors . Moreover it has an exact degeneracy between the first $`0^+`$ excited state and the first $`1^{}`$ states , which can be respectively identified with the Roper resonance and the negative parity resonances. This degeneracy seems to be in agreement with phenomenology and is typical of an underlying O(7) symmetry . This feature cannot be reproduced in models with only two-body forces and/or harmonic oscillator bases since the excited $`L=0`$ state, having one more node, lies above the $`L=1`$ state .

The dynamic symmetry $`O(7)`$ of the hyperCoulomb problem can be used to obtain the eigenvalues using purely algebraic methods, similarly to what is done in the hydrogen atom case with the $`O(4)`$ symmetry. In fact, the hyperCoulomb Hamiltonian can be rewritten as

$$H=\frac{\tau ^2m}{2[C_2(O(7))+\frac{25}{4}]},$$
(6)
where $`C_2(O(7))`$ is the quadratic Casimir invariant of $`O(7)`$, and the energy eigenvalues became

$$E=\frac{\tau ^2m}{2(n+5/2)^2},$$
(7)
where n is a non negative integer .

As a confining hypercentral potential in our model we have assumed a form

$$V(x)=\frac{\tau }{x}+\alpha x,$$
(8)
that means a coulomb-like term plus a linear confining term as suggested by lattice QCD calculations . In order to describe the splittings within the $`SU(6)`$-multiplets we introduce a hyperfine interaction of the standard form  and we treat it as a perturbation. Having fixed the quark mass $`m`$ to $`1/3`$ of the nucleon mass, the remaining three free parameters ($`\tau `$, $`\alpha `$ and the strength of the hyperfine interaction) are fitted to the spectrum. The strength of the hyperfine interaction is determined by the $`\mathrm{\Delta }`$ \- Nucleon mass difference and the spectrum is described with $`\tau =4.59`$ and $`\alpha =1.61fm^2`$. Having fixed the parameters of the potential, the wave functions of the various resonances are completely determined and have been used for the calculation of the photocouplings , the transition form factors to the negative parity resonances , the elastic form factors  and the ratio between the electric and magnetic form factors of the proton .

The resulting overall description of the experimental data is quite good, however there are some problems. As far as the spectrum is concerned, the model fails to reproduce the Roper resonance and the higher levels. The evaluated size of the nucleon is too low, since the resulting r.m.s. radius is $`0.46fm`$ and this affects the description of the elastic form factors. The introduction of relativistic corrections substantially improves the calculated form factors  but it is not enough. The strength of the helicity amplitudes are underestimated, as it happens in all constituent quark models.

There may be different reasons for such discrepancies. The splittings within the multiplets are not all adequately described by the hyperfine interaction. One should remind that the form we have assumed for the hyperfine interaction contains a $`\delta `$-like term, which is troublesome from the theoretical point of view. Furthermore, the splittings can be originated also by other terms, for instance isospin-dependent ones .

Another possible reason of the discrepancies is the fact that the model does not contain any explicit quark-pair creation mechanism, which is expected to be particularly important for the description of the eletromagnetic excitation  but may lead to some residual effects also in the spectrum . In the latter case, the creation of quark-antiquark pairs could be the microscopic origin of an isospin dependent part of the potential. In this respect the constituent quark potential has to be considered as an effective potential in the three-quark subspace, taking into account implicitly the missing Fock-space configurations.

## 3 The spin and isospin splittings

The standard hyperfine interaction is used in order to reproduce the splittings within the $`SU(6)`$multiplets. As mentioned above, it contains a $`\delta `$-like term which is an illegal operator. For this reason we have modified it by introducing a smearing factor given by a gaussian function of the quark pair relative distance :

$$H_S=A_S\underset{i<j}{}\frac{1}{(\sqrt{\pi }\sigma _S)^3}e^{\frac{r_{ij}^2}{\sigma _S^2}}(\text{s}_i\text{s}_j)\text{ ,}$$
(9)
where $`\text{s}_i`$ is the spin operator of the i-th quark and $`r_{ij}`$ is the relative quark pair coordinate. The results for the spectrum are shown in Fig. 1. The fitted parameters are $`\alpha =1.58fm^2`$, $`\tau =4.98`$,$`A_S=38.4fm^2`$,$`\sigma _S=0.8fm`$. The correct limit for vanishing smearing is obtained, in the sense that the spectrum of Ref.  is reproduced.

We have also tried a smearing depending on the hyperradius $`x`$ only of the form :

$$H_{Spin}=B_S\left(\frac{1}{\mathrm{\Lambda }_S}e^{\frac{x}{\mathrm{\Lambda }_S}}\right)\underset{i<j}{}(\text{s}_i\text{s}_j)\text{ ,}$$
(10)
where $`\text{s}_i`$ is the spin operator of the i-th quark and $`x`$ is the hyperradius. The fitted parameters are $`\alpha =1.49fm^2`$, $`\tau =5.01`$,$`B_S=196.4fm^2`$,$`\mathrm{\Lambda }_S=1.6fm`$. There is an improvement for the higher states, however there is a too strong degeneracy of the levels and the good agreement for the negative parity states obtained in the hCQM  is lost. Therefore in the following we shall use only the hyperfine interaction with a two-body smearing.

As quoted in the previous sections, there are different motivations for the introduction of a flavour dependent term in the three-quark interaction. The well known Guersey-Radicati mass formula  contains a flavour dependent term, which is essential for the description of the strange baryon spectrum:

$`\widehat{M}=M_0+a\widehat{C}_2(SU_{SF}(6))+b\widehat{C}_2(SU_F(3))`$
$`+b^{}\widehat{C}_2(SU_I(2))+b^{\prime \prime }\widehat{C}_1(U_Y(1))`$
$`+b^{\prime \prime \prime }\widehat{C}_2(U_Y(1))+\widehat{C}_2(SU_S(2)),`$ (11)
where $`M_0`$ is fixed for any $`SU(6)`$-multiplet, $`\widehat{C}_2(SU(N))`$ is the Casimir of SU(N), $`Y`$ is the hypercharge, $`I`$ is the total isospin and $`S`$ the total spin. For non strange baryons this formula implies an isospin dependence. In the algebraic description of baryon properties , the space part of the mass operator is written in terms of the generators of the U(7) group, while for the internal degrees of freedom the Guersey-Radicati mass formula  is used. In the chiral Constituent Quark Model , the non confining part of the potential is provided by the interaction with the Goldstone bosons, giving rise to a spin- and isospin-dependent part, which is crucial for the description of the spectrum for energies lower than $`1.7GeV`$. It has been also pointed out quite recently that an isospin dependence of the quark potential can be obtained by means of quark exchange . More generally, one can expect that the quark-antiquark pair production can lead to an effective quark interaction containing an isospin (or flavour) dependent term. On the other hand, the fact that the constituent quark model does not contain explicitly this mechanism is may be the reason why the low $`Q^2`$ -behaviour of the electromagnetic transition form factors is not reproduced .

With these motivations in mind, we have introduced isospin dependent terms in the hCQM hamiltonian.

To this end we have added two terms in the three-quark hamiltonian with the hyperfine interaction of Eq.(9). The first one depends on the isospin only and has the form:

$$H_\mathrm{I}=A_I\underset{i<j}{}\frac{1}{(\sqrt{\pi }\sigma _I)^3}e^{\frac{𝐫_{ij}^2}{\sigma _I^2}}(\text{t}_i\text{t}_j)\text{ ,}$$
(12)
where $`\text{t}_i`$ is the isospin operator of the i-th quark and $`r_{ij}`$ is the relative quark pair coordinate. The second one is a spin-isospin interaction, given by

$$H_{\mathrm{SI}}=A_{SI}\underset{i<j}{}\frac{1}{(\sqrt{\pi }\sigma _{SI})^3}e^{\frac{r_{ij}^2}{\sigma _{SI}^2}}(\text{s}_i\text{s}_j)(\text{t}_i\text{t}_j)\text{ ,}$$
(13)
where $`\text{s}_i`$ and $`\text{t}_i`$ are respectively the spin and isospin operators of the i-th quark and $`r_{ij}`$ is the relative quark pair coordinate. The complete interaction is then given by

$$H_{int}=V(x)+H_\mathrm{S}+H_\mathrm{I}+H_{\mathrm{SI}}.$$
(14)
The resulting spectrum for the 3\*- and 4\*- resonances is shown in Fig.2 and in Table I. The $`N\mathrm{\Delta }`$ mass difference is no more due only to the hyperfine interaction. In fact, in this model its contribution is only about $`35\%`$, the remaining splitting comes from the spin-isospin term, $`(50\%)`$, and from the isospin one, $`(15\%)`$.

It should also be noted that the inversion between the Roper and the negative parity resonances is almost entirely due to the spin-isospin interaction, as stated in Ref. .

The tensor term coming from the hyperfine interaction has been kept as well, however its contribution to the spectrum is negligible.

In Table II and Table III we list all the remaining states predicted by the model and when possible we give a temptative assignation to one- and two-star states. We see that the number of predicted states is higher than the presently observed ones, that is, as in other models, we have the problem of missing resonances. It is interesting to observe that while the quality of the reproduction of the spectrum for the 3\*- and 4\*- resonances is quite the same for the various Constituent Quark Models , there are different predictions for the other states, that means for 1\*, 2\* and missing resonances.

Recently in a three-channel multi-resonance amplitude analysis it has been found evidence for a third low-lying $`P_{11}`$ state at $`1740\pm 11MeV`$ . The first two $`P_{11}`$ states at $`1439\pm 19MeV`$ and $`1729\pm 16MeV`$ correspond to the N(1440) and N(1710) of the PDG . In the hCQM the first three $`P_{11}`$ states are at $`1463MeV`$, $`1752MeV`$ and $`1828MeV`$ respectively. A new analysis of kaon photoproduction data  has shown evidence for a third $`D_{13}`$ resonance at $`1895MeV`$, which can be described by one of the states predicted by the present model (see Table II).

## 4 Discussions and conclusions

A considerable improvement in the description of the spectrum is obtained with an isospin dependent potential. As quoted in the previous section, a possible motivation of the isospin-dependent terms of the quark interaction is given by quark-antiquark pair production mechanisms.

This kind of mechanisms have been studied within the string model in connection with the quenching problem of quark models in the meson spectrum . Their contributions turn out to be mainly spin-independent and can be reabsorbed into the string parameter, so the residual part can be considered as a perturbative correction. This can be considered a suggestion also for the baryon case, even if there are differences. Therefore, any interaction leading to unperturbed states with reasonable spin-averaged values of the energy levels can be considered to contain implicitly this kind of contribution. The splittings, which are in general spin and isospin dependent, can be treated perturbatively.

In this article we have shown that the complete interaction including spin and isospin terms (see Eq.(14) ), reproduces the position of the two Roper resonances of the nucleon, while keeping the good description of the negative parity resonances. It should be noted that also the higher states are accounted for. The hypercentral potential is a good starting point for the construction of an unperturbed spectrum and leads to realistic quark states, as it is shown also by the reproduction of the e.m. form factors and transition form factors, which are sensitive to the wave functions. The configuration mixing is usually a higher order correction, apart from some special cases as the neutron form factor or the $`S_{11}(1650)`$ helicity amplitude etc. , for which the $`SU(6)`$-contribution is vanishing. In this respect, the forthcoming more precise data that will be soon available will supply valuable information concerning the possible $`SU(6)`$-breaking terms in the quark interaction. Finally, one can observe that the various Constituent Quark Models give different results concerning the 1\* and 2\* states and the position and number of the missing resonances . Therefore, the expected new data coming from the TJNAF will be very helpful in order to discriminate among thems.


text:
# Colloidal Dynamics on Disordered Substrates

\[

## Abstract

Using Langevin simulations we examine driven colloids interacting with quenched disorder. For weak substrates the colloids form an ordered state and depin elastically. For increasing substrate strength we find a sharp crossover to inhomogeneous depinning and a substantial increase in the depinning force, analogous to the peak effect in superconductors. The velocity versus driving force curve shows criticality at depinning, with a change in scaling exponent occuring at the order to disorder crossover. Upon application of a sudden pulse of driving force, pronounced transients appear in the disordered regime which are due to the formation of long-lived colloidal flow channels.

\]

Colloidal crystals are an ideal system in which to study the general problem of ordering and dynamics in 2D , since the particle size permits direct imaging of the particle locations and motion. A considerable amount of work has been conducted on the melting of 2D colloidal crystals in the absence of a substrate . In addition, a number of experimental and theoretical studies have considered colloidal crystallization and melting in 2D systems with periodic 1D  and 2D substrates , where a rich variety of crystalline states can be stabilized.

Colloid crystals are also ideal for studing the ordering and dynamics of an elastic media interacting with random substrates, a problem that is relevant to a wide variety of systems such as superconducting vortices, Wigner crystals, and charge-density waves (CDWs). Open issues include the nature of the dynamical response to applied forces, as well as whether an order to disorder transition occurs as the strength of the random substrate increases. Recently, LeDoussal and Carpentier have theoretically investigated the effects of quenched disorder on the order and melting of 2D lattices and find a sharp crossover from the ordered Bragg glass (where defects are absent) to a disordered or molten state . They predict that the depinning threshold increases at this crossover due to the softening of the lattice, which allows the particles to better adjust to the substrate. A similar mechanism could account for the peak effect observed in vortex matter in superconductors , in which the depinning threshold rises dramatically when the applied magnetic field is increased. In low temperature superconductors, where the vortices are fairly stiff so that their behavior can be considered as effectively 2D, recent small angle neutron scattering experiments have shown that the peak effect is associated with a sharp disordering or melting transition .

In addition to static properties, the dynamics of elastic media interacting with quenched disorder in 2D is a topic of intense study. In the disordered region the driven system may break up into pinned and flowing regions, as observed in experiments  and simulations  of vortices in superconductors. Conversely, for weak substrate disorder the elastic media is defect free and undergoes elastic depinning, in which the particles keep the same neighbors as they move. Fisher predicted that elastic depinning would show criticality , and that the velocity vs. force curves would scale as $`v=(ff_c)^\beta `$, where $`f_c`$ is the depinning threshold force. This scaling has been studied extensively in 2D CDW systems where $`\beta =2/3`$ . It is, however, not known whether this exponent occurs in other systems undergoing elastic flow. Another intriguing dynamical phenomena is the pronounced transient behavior exhibited by vortices under a sudden applied current pulse at magnetic fields near the peak effect regime . Due to surface barrier effects in the vortex sample, it is not clear whether these transient effects relate to the plasticity of the vortex dynamics or to contamination of the vortex lattice by disorder from the sample edges . Recently, Pertsinidis and Ling  have studied colloids in 2D driven by an electric field and interacting with a disordered substrate. They observe plastic depinning with filamentary or river-like flow of colloids and a velocity-force curve scaling with $`\beta =2.2`$. Under a pulsed drive the system shows very long time transients that fit to a stretched exponential.

Motivated by the recent colloidal experiments as well as the pulse drive experiments in vortex matter, we have conducted Langevin simulations of colloidal particles interacting via a Yukawa potential in 2D systems with random disorder. In simulation, the strength of the disorder can be carefully tuned, which is difficult to achieve in experiments. In addition, the initial conditions of the colloidal arrangements are easily controlled, whereas in experiments, defects generated in the colloidal lattice during preparation may become frozen in by the disorder.

Our simulations show that for weak substrates the colloids form an ordered triangular array which depins elastically without the generation of defects. For increased substrate strength, there is a sharp crossover to a disordered phase where the colloids depin plastically into riverlike structures. This crossover is accompanied by a sharp increase in the depinning threshold, analogous to the peak effect phenomena in superconductors. We

find scaling of the velocity vs. applied drive with an exponent of $`\beta =0.67`$ in the elastic regime, in agreement with studies in 2D CDW’s. In the plastic regime we find $`\beta =1.97`$, close to the value observed in experiments . In the disordered region, long time transients that fit to a stretched exponential occur in response to a sudden applied drive pulse, as also observed in the colloidal and vortex experiments.

The colloids are simulated using Langevin dynamics in 2D . The colloids interact via a Yukawa or screened Coulomb interaction potential $`V(r_{ij})=(Q^2/|𝐫_i𝐫_j|)\mathrm{exp}(\kappa |𝐫_i𝐫_j|)`$. Here $`Q`$ is the charge of the particles, $`1/\kappa `$ is the screening length, and $`𝐫_{i(j)}`$ is the position of particle $`i(j)`$. The length of the system is measured in units of the lattice constant $`a_0`$ and a screening length of $`\kappa =2/a_0`$ is used. The quenched disorder on the substrate is modeled as randomly placed parabolic traps with radius $`r_p<a_0`$ and a maximum force $`f_p`$. This same type of model for pinning has been used previously to model quenched disorder in superconducting vortex systems . Other types of pinning potentials used in vortex simulations produce results similar to the parabolic pinning . The equation of motion for colloid $`i`$ is

$$\frac{d𝐫_i}{dt}=𝐟_{ij}+𝐟_p+𝐟_T+𝐟_d$$
(1)
Here $`𝐟_{ij}=_{ji}^{N_c}_iV(r_ij)`$ is the interaction force from the other colloids, $`𝐟_p`$ is the pinning force, $`𝐟_T`$ is a randomly fluctuating force due to thermal kicks, and $`𝐟_d`$ is the force due to an applied driving field. We start the system at a temperature above melting as determined from the diffusion, and gradually cool below $`T_m`$ to $`T/T_m=0.4`$ for most of the data presented here. To measure velocity $`v`$ vs force curves, care must taken to average over substantial amounts of time to avoid transient effects. The driving force is increased from zero by small increments and the velocity is averaged for $`5\times 10^4`$ time steps at each increment, with typical simulations running for $`10^7`$ time steps. In this model we do not take into account hydrodynamic effects or long-range attractions between colloids.

In Fig. 1(a) we show the depinning force $`f_c`$ vs. substrate strength $`f_p`$ from a series of simulations. For $`f_p<0.18`$ the depinning force increases as a power law, $`f_cf_p^{1.9\pm 0.1}`$. To compare the depinning force to the order in the system, in Fig. 1(b) we show the the percentage of defects or non-six fold coordinated particles $`P_d`$ as calculated from a Delaunay triangulation. This measure indicates that the colloidal crystal is in an ordered state ($`P_d=0.0`$) for $`f_p<0.18`$, and that there is a crossover to a disordered state ($`P_d0`$) at $`f_p=0.18`$. In Fig. 1(c) we show a representative Delaunay triangulation for the ordered state where there are no defects but small distortions in the particle positions can be seen, and in Fig 1(d) we show the disordered state where defects are present. The crossover to the disordered state coincides with a rapid increase in the depinning force as seen in Fig. 1(a) and in the inset of Fig. 1(a) which shows a peak in $`df_c/df_p`$ at the crossover. This behavior is consistent with the recent experiments in superconductors which find that at the peak effect there is an increase in the pinning with a simultaneous disordering of the lattice . For $`f_p>0.2`$, the depinning scales as $`f_cf_p`$, as expected for the single particle pinning regime. The sudden increase in the depinning force results from the fact that the defected colloid lattice is much softer than the ordered lattice, allowing the colloids to adjust their positions to accommodate to the optimal pinning sites. We have also investigated this transition for different colloidal densities and disorder strengths. For increasing $`T`$ the order to disorder transition is shifted to lower values of $`f_p`$. We have also investigated finite size effects for increasing system sizes, and find that the order-disorder crossover shifts only a small amount before saturating, while the sharpness of the transition persists with increased system size.

It is beyond the scope of this paper to determine whether the order to disorder crossover is a first order transition. Although the sharpness suggests a possible first order transition, Carpentier and LeDoussal show that for 2D systems with quenched disorder, a sharp disordering crossover, rather than transition, occurs . In addition, a first order transition is not expected since the Bragg glass in 2D has been shown to have dislocations on large scales at all temperatures. The distance between these dislocations can be arbitrarily large .

In Fig. 2 we show that the order-disorder crossover coincides with the onset of plastic flow above depinning. In Fig. 2(b) the elastic colloid flow is shown for $`f_p=0.12`$ above depinning ($`f_d/f_c=1.1`$). Here each colloid keeps the same neighbors as it moves. In Fig. 2(a) the inhomogeneous or plastic colloidal flow is shown for $`f_d/f_c=1.1`$ for $`f_p=0.25`$. Here, only a portion of the colloids are moving at any one time, and the motion occurs in channels or rivers between pinned regions. The colloid velocities show a bimodal distribution in this regime, split between the stationary and moving colloids. In addition, the channels seen in Fig. 2(a) are not static but change over time, so that any one colloid is only temporarily trapped in a pinning site. These features of the plastic flow are in agreement with observations in colloidal experiments  and in vortex simulations of the strongly pinned regime .

In order to correlate the different types of flow observed in Fig. 2 with properties of bulk measurements, we show in Fig. 3 the scaling of the velocity vs driving force. For elastic depinning in the ordered regime \[Fig. 3(a)\], the velocity vs force curves fit well to $`v=(f_df_c)^\beta `$ with $`\beta =0.66\pm 0.02`$, as illustrated in Fig. 3(b). These results are in good agreement with theoretical predictions  and simulation results  for CDW’s in 2D where the depinning is elastic. In contrast, in driven 2D vortex matter, Higgins and Bhattacharaya  found an exponent of $`\beta =1.2`$ below the peak effect where elastic flow is expected to occur. This may be due to the effects of surface barriers disordering the lattice. In Fig. 3(c,d)

the $`vf_d`$ scaling for the plastic regime is presented, with $`\beta =1.94\pm 0.03`$, close to the value of $`2.2`$ found in the colloid experiments . For larger system sizes, we find that the scaling region is expanded but the exponent is unchanged. Other studies in the plastic flow regime found $`\beta =5/3`$ for electron flow in metallic dots  and $`\beta =2.22`$ for vortex flow in Josephson-junction arrays .

In Fig. 4 we show the response of colloids prepared in an ordered state to the application of a sudden pulse of driving force of different strengths in the plastic flow regime. Since the pulse strength is chosen to be below the depinning threshold value $`f_c`$, the initial colloid velocity is high, and then gradually decreases. We find that a simple functional form cannot be fit to the curves. Instead, we use a stretched exponential fit as performed in experiments : $`v(t)=v_0\mathrm{exp}((t/t_0)^\alpha )+v_1`$. The values of $`t_0`$ and $`\alpha `$ depend on the magnitude of the drive. For the parameters investigated here, $`\alpha `$ falls between $`0.08`$ and $`0.4`$, in agreement with the values found in experiment. A similar stretched exponential decay was also found in vortex matter for the transient response to pulses . We find that in the long time limit, the colloid flow occurs only through a few long-lived channels. In the elastic regime, the decay of $`v`$ is much faster and fits to an initial pure exponential with the velocities going to zero. We note that in the elastic regime the colloids move less than a lattice constant after a pulse is applied, whereas in the plastic regime, colloids in the moving channels can move the entire length of the system. For increased system sizes, the transient times are enhanced in the plastic flow regime but are unchanged in the elastic regime. The long

lived transients in the plastic regime are responsible for the very slow velocity-force sweep necessary to measure an accurate depinning threshold. This sweep-rate dependence is also consistent with the experimentally observed sweep rate dependent critical currents in the peak regime , where slow rates produce larger measured critical currents.

To summarize, we investigated the behavior of 2D colloids interacting with random disorder using Langevin simulations. For weak disorder the colloids form an ordered lattice which depins elastically and shows critical scaling in the velocity vs force curves, with $`\beta =0.67`$, in agreement with studies of 2D CDW’s. For increasing disorder strength, we find a sharp crossover to a disordered state, accompanied by a sharp increase in the depinning force, analogous to the peak effect observed for vortex matter in superconductors. In the disordered region, the colloids depin inhomogeneously into fluctuating channels and the $`vf`$ curve scaling gives $`\beta =1.97`$, in agreement with experiments. In the disordered flow regime, pronounced transients occur in response to a sudden pulse, with the late time dynamics determined by a few long lived channels.

Acknowledgments. We thank S. Ling and A. Pertsinidis for sharing their data before publication and for useful discussions. We also thank S. Bhattacharya, A. Bishop, D. Fisher, D.G. Grier, M.B. Hastings and A.A. Middleton for helpful comments. This work was supported by the US Department of Energy under contract W-7405-ENG-36.


text:
# Numerical Study of the Spin-Flop Transition in Anisotropic Spin-1/2 Antiferromagnets

## Abstract

Magnetization processes of the spin-1/2 antiferromagnetic $`XXZ`$ model in two and three spatial dimensions are studied using quantum Monte Carlo method based on stochastic series expansions. Recently developed operator-loop algorithm enables us to show a clear evidence of the first-order phase transition in the presence of an external magnetic field. Phase diagrams of closely related systems, hard core bosons with nearest-neighbor repulsions, are also discussed focusing on possibilities of phase-separated and supersolid phases.

There have been a lot of interests on magnetic properties of anisotropic quantum antiferromagnets since Néel first predicted the first-order phase transition in the presence of an external magnetic field neel . One of the simplest models for anisotropic antiferromagnets in an external magnetic field is described by the following spin-1/2 $`XXZ`$ model Hamiltonian:

$$H=J\underset{𝐢,𝐣}{}\left(S_𝐢^xS_𝐣^x+S_𝐢^yS_𝐣^y+\mathrm{\Delta }S_𝐢^zS_𝐣^z\right)h\underset{𝐢}{}S_𝐢^z,$$
(1)
where $`S_𝐢^\alpha `$ is the $`\alpha `$$`(=x,y,z)`$ component of the spin-1/2 operator at site i, $`h`$ is an external magnetic field applied in the z-direction, and $`𝐢,𝐣`$ runs over all the nearest-neighbor pairs of spins at sites i and j. $`J(>0)`$ is antiferromagnetic coupling and $`\mathrm{\Delta }(0)`$ is anisotropic constant. Mean-field calculations of the spin-1/2 $`XXZ`$ model mf2 , supporting Néel’s prediction, found the first-order phase transition from the Néel ordered state to the spin-flopping state with increasing the magnetic field $`h`$. On the contrary to these studies, it is known from Bethe ansatz solution that the one-dimensional (1D) spin-1/2 $`XXZ`$ model shows a second-order transition in the presence of the external magnetic field yang . The discrepancy is assigned to the inadequacy of treating quantum fluctuations by the mean-field theories. It is thus important to use an unbiased numerical method for understanding the correct nature of the magnetization process even for the simplest systems such as one given by Eq. (1) in higher spatial dimensions since there exist no analitically exact solutions. This is preciely one of our purposes for this paper takahashi .

Another importance of studing the spin-1/2 $`XXZ`$ model defined by Eq. (1) comes from the fact that the model is mapped onto a system of hard core bosons with nearest-neighbor repulsions described by the Hamiltonian

$$H_\mathrm{B}=t\underset{𝐢,𝐣}{}\left(c_𝐢^{}c_𝐣+c_𝐣^{}c_𝐢\right)+V\underset{𝐢,𝐣}{}n_𝐢n_𝐣\mu \underset{𝐢}{}n_𝐢,$$
(2)
with $`t=J/2`$, $`V=J\mathrm{\Delta }`$, and $`\mu =h+zJ\mathrm{\Delta }/2`$ ($`z`$: coordination number)matsubara . Here $`c_𝐢^{}`$ is a creation operator of hard core boson at site $`𝐢`$ and $`n_𝐢=c_𝐢^{}c_𝐢`$. The total magnetization $`M_z=_𝐢S_𝐢^z`$ thus relates to $`_𝐢\left(n_𝐢1/2\right)`$ in the boson model. The boson Hamiltonian $`H_\mathrm{B}`$ is proposed as a model Hamiltonian to study properties of liquid $`{}_{}{}^{4}\mathrm{He}`$ fisher , granular superconducting arrays cha , and flux lines in superconductors nelson . The ground state phase diagram has been studied using mean-field theories mf2 ; mf  and numerical methods bat2 , and shows a Mott insulating phase, a superfluid phase, and a phase with having both the orders simultaneously ($`supersolid`$ phase). The correspondence of these states to those in the spin model is as follows: Mott insulating and Néel states, superfluid and spin-flopping states, and supersolid and “intermediate” spin states mf2 , respectively.

The main purpose of this paper is to show a clear evidence for the first-order phase transition of the 2D and 3D spin-1/2 $`XXZ`$ models in the presence of a magnetic field by using recently developed numerical method and present the ground state phase diagrams. The presence of a phase-separated phase and the absence of the supersolid phase in the closely related system of hard core bosons with nearest-neighbor repulsions are also discussed.

The magnetization process of the 2D and 3D spin-1/2 $`XXZ`$ models defined by Eq. (1) is studied numerically on square (number of spins $`N_s=L\times L`$) and cubic ($`N_s=L\times L\times L`$) lattices using quantum Monte Carlo (QMC) technique based on stochastic series expansions (SSE) sandvik . Very recently an important technical improvement have been achieved by Sandvik loop . He found a new algorithm of cluster-type updates (operator-loop updates) within the SSE QMC scheme which reduces autocorrelation time drastically compared to simulations using only local updates. While this method is very similar to the loop algorithm in the world-line QMC method proposed by Evertz et al. evertz , one major advantage of the SSE method with operator-loop updates is that there is no difficulty in simulating systems with anisotropic couplings in external magnetic fields owing to not having to have “freezing” configurations and “global” weights which make the loop algorithm in the world-line QMC method highly inefficient loop . This reduced autocorrelation time enables us to go down to very low temperatures in very high magnetic fields and therefore the method is suitable for our purpose. Another advantage of this study using the SSE scheme over other earlier numerical studies takahashi ; bat2 ; bat  is that the simulations are performed directly in the ground canonical ensemble, i.e., magnetization per site $`m_z=M_z/N_s`$ is calculated for a given magnetic field $`h`$. In this paper temperatures $`T`$ are set to be $`J/k_\mathrm{B}T=2L`$ which is low enough to study the ground state properties on finite lattices temp , and periodic boundary condition is used. The exchange coupling $`J`$ will be taken as the energy unit. Since there exists long-range Néel-ordered state only for $`\mathrm{\Delta }1.0`$, from which the first-order spin-flop transition can take place by applying a finite external field, our main focus in this paper will be put on this anisotropic regime.

Let us first study the magnetization process in 2D. A typical example of the magnetization curves is shown in Fig. 1 (a). One can see that at certain critical magnetic field $`h_c`$ the magnetization $`m_z`$ changes discontinuously from 0 to a finite value $`m_z^c`$. For this example in the figure with $`\mathrm{\Delta }=1.5`$ the magnetization jumps at $`h_c1.83`$ to $`m_z^c0.11`$ boson . Working on various system sizes (see in the inset of Fig. 1(a)), we confirm that finite size effects are small and conclude that the jump in magnetization is not due to extrinsic factors of working on finite-size lattices. The result already gives a clear evidence for the first-order phase transition.

In order to illustrate that autocorrelation times of our QMC measurements are short enough, the integrated autocorrelation time $`\tau _{\mathrm{int}}`$ is estimated as follows: we first divide a sequence of Monte Carlo data points into bins of length $`l`$, and for each bin length $`l`$, the average of the data in the $`b`$th bin and the variance $`\sigma ^2(l)`$ of the bin averages are calculated. The integrated autocorrelation time $`\tau _{\mathrm{int}}`$ is then estimated by the asymptotic value of $`\tau (l)=l\sigma ^2(l)/2\sigma ^2(l=1)`$ at large $`l`$ above which $`\tau (l)`$ does not depend on $`l`$ kawa . The results of $`\tau _{\mathrm{int}}`$ for $`m_z`$ and spin structure factor $`S(𝐪)`$ at $`𝐪=(\pi ,\pi )`$ are shown in Fig. 1 (b) for $`\mathrm{\Delta }=1.5`$ and $`h=2.0`$. From this figure $`\tau _{\mathrm{int}}`$’s are estimated to be no longer than 4–5 MC steps. Since we use bin length $`l=1000`$ or more to estimate the statistical errors, we can be sure that our error bar is accurate. It should be noted here that $`\tau _{\mathrm{int}}`$ can be controlled by changing how many times closed loops are constructed in the operator-loop updates. In this paper we have to construct closed loops as many as about 1000 times per MC step for larger values of $`\mathrm{\Delta }`$.

Repeating the same procedure for different values of $`\mathrm{\Delta }`$, the ground state phase diagram of the 2D spin-1/2 $`XXZ`$ model is completed. The result is given in Fig. 1 (c). There exist three different phases denoted by (i) $`m_z=0`$, (ii) $`m_z0`$, and (iii) $`m_z=0.5`$ in the Figure. The third phase (iii) corresponds to a fully saturated ferromagnetic phase and is separated from the second phase (ii) by the critical magnetic field $`h_c^{\mathrm{max}}`$. This transition is trivial and is not our interest here. The critical magnetic field $`h_c^{\mathrm{max}}`$ is indeed easily calculated by going to the boson model described by Eq. (2): the critical chemical potential to have just one particle in the $`d`$-dimensional boson system is $`2dt`$ and therefore $`h_c^{\mathrm{max}}/J`$ is found to be $`d(1+\mathrm{\Delta })`$.

To show further evidence for the first-order transition between phase (i) and phase (ii) we calculate the spin structure factor $`S(𝐪)=1/N_s_{𝐢,𝐣}\mathrm{e}^{i𝐪(𝐢𝐣)}𝐒_𝐢𝐒_𝐣`$ at $`𝐪=(\pi ,\pi )`$ and the spin stiffness (helicity modulus) $`\rho _s`$ as a function of $`h`$ for a fixed $`\mathrm{\Delta }`$. The spin stiffness $`\rho _s`$ is calculated by, for example, $`\rho _s=Tw_x^2+w_y^2(+w_z^2)/dL^{d2}`$ in 2D (3D) where $`w_\alpha `$ is the winding number per linear spatial lattice size $`L`$ in the $`\alpha `$-direction ceperley . $`\rho _s`$ corresponds to the superfluid density in the boson model and is used to detect the superfluid phase in the system mf ; bat2 ; ceperley . The results are shown in Fig. 2(a) for $`\mathrm{\Delta }=1.5`$ as a function of $`h`$. When $`h`$ is small $`S(𝐪)`$ has a peak at $`𝐪=(\pi ,\pi )`$ ($`𝐪`$ dependence is not shown) and $`\rho _s=0`$. With increasing $`h`$ these quantities change discontinuously at $`h_c`$ where the magnetization $`m_z`$ jumps, and $`S(\pi ,\pi )`$ becomes zero while $`\rho _s`$ has a finite value com . Apparently these two phases are (i) Néel ordered phase (for $`h<h_c`$) and (ii) spin-flopping phase (for $`h>h_c`$), and are separated through the first-order transition note .

Earlier studies found the supersolid phase in the boson model characterized by having finite values of $`S(\pi ,\pi )`$ and $`\rho _s`$ simultaneously in a region between phase (i) and phase (ii) of the phase diagram mf ; bat2 . One can indeed barely see small but finite values of $`S(\pi ,\pi )`$ for $`h>h_c`$ where $`\rho _s`$ is finite (see Fig.2 (a)). In order to elucidate systematically the possibility of the existence of the supersolid phase in 2D, we carry out the finite size scaling analyses of $`S(\pi ,\pi )`$ and $`\rho _s`$ for fixed magnetizations. The results are presented in Fig. 2(b). One can see that $`S(\pi ,\pi )/N_s`$ ($`\rho _s`$) stays finite in the limit of $`N_s\mathrm{}`$ for $`h<h_c`$ ($`h>h_c`$) while $`S(\pi ,\pi )/N_s`$ approaches to zero at $`N_s\mathrm{}`$ for $`h>h_c`$. Doing the similar analyses for different values of $`\mathrm{\Delta }`$ and $`m_z`$ it is found that whenever $`h>h_c`$ $`S(\pi ,\pi )/N_s`$ goes to zero in the thermodynamic limit. We therefore conclude that the supersolid phase does not exist in the 2D spin-1/2 $`XXZ`$ model. The results are consistent with very recent studies by Batrouni and Scalettar bat .

We now carry out the similar calculations for the 3D spin-1/2 $`XXZ`$ model to study the nature of the phase transition induced by the external magnet field. Our strong evidence for the first-order transition is provided in Fig. 3(a). In the Figure magnetization $`m_z`$, spin structure factor $`S(\pi ,\pi ,\pi )`$, and spin stiffness $`\rho _s`$ for $`\mathrm{\Delta }=1.5`$ are plotted as a function of the external magnetic field $`h`$. It is clearly seen that as in the case of the 2D model those quantities change discontinuously at the critical magnetic field $`h_c`$. Working with different values of $`\mathrm{\Delta }`$, the ground state phase diagram is constructed and the result is shown in Fig. 3(b). As in 2D the phase diagram consists of three different phases: (i) Néel ordered phase, (ii) spin-flopping phase, and (iii) fully saturated ferromagnetic phase note2 . The critical magnetic field $`h_c`$ in 3D is observed to be larger compared to that in 2D for a given $`\mathrm{\Delta }`$. This is simply because the increased coordination number in 3D makes the needed magnetic field larger to destroy the Néel state. The similar finite size scaling analyses which was done for the 2D model do not find the supersolid phase in the 3D model.

Finally we summarize our results in Fig. 4 by showing the ground state phase diagrams of the 2D and 3D spin-1/2 $`XXZ`$ models on the parameter space of anisotropic constant $`\mathrm{\Delta }`$ and magnetization $`m_z`$. The phase diagrams show (i) the Néel ordered state at $`\mathrm{\Delta }1.0`$ and $`m_z=0.0`$, (ii) the spin-flopping state with $`\rho _s0`$, (iii) the fully saturated ferromagnetic state at $`m_z=0.5`$, and (iv) a phase-separated state (denoted by PS in the Figures). The phase-separated region exists because as was seen in Figs. 1(a) and 3(a) there is a region in magnetization $`m_z`$ to which we can not reach in a thermodynamically stable way no matter how fine the magnetic field $`h`$ is tuned. In other words, if one could have a states with $`m_z`$ which is in this magnetization region, the state would be phase-separated between the Néel state with $`m_z=0`$ and the spin-flopping state with $`m_z=m_z^c`$. Some earlier studies predicted two phase coexisting regions working with the canonical ensemble (i.e., fixed $`m_z`$) and their results were interpreted as the supersolid phase mf ; bat2 . However our calculations conclude that these phases are thermodynamically unstable and are phase-separated. Our conclusions are consistent with recent studies by Batrouni and Scalettar for the 2D boson model bat .

In conclusion, we have studied numerically the magnetization process of the spin-1/2 anisotropic $`XXZ`$ model in two and three spatial dimensions using QMC method based on SSE and shown clear evidences of the first-order phase transition in the presence of an external magnetic field. Based on the calculated ground state phase diagrams, the existence of a phase-separated region and the absence of supersolid phases are pointed out in the related systems of hard core bosons with nearest-neighbor repulsions.

It would be of great interest to study effects of random anisotropic constant $`\mathrm{\Delta }`$ on the phase diagrams wallin  since one can elucidate $`quantum`$ effects on recently proposed impurity-induced quantum-critical-point-like behavior near $`first`$-order phase transitions in the Ising models elbio . Another interesting issue to address is to see effects of long-ranged Coulomb repulsions between bosons on the phase-separated region in the phase diagram. There is a general belief that introducing the Coulomb interactions replaces the phase-separated state by thermodynamically stable states such as a droplet-like state and perhaps a stripe state moreo . The model studied here provide an ideal system to examine the possibilities of those exotic states using unbiased numerical methods.

The author thanks E. Dagotto and A. Dorneich for valuable discussions.


text:
# Josephson Plasma Excitation and Vortex Oscillation Mode in Josephson Vortes State

## I Introduction

When the external magnetic field is applied parallel to the $`ab`$-plane of high-$`T_c`$ superconductors such as Bi2Sr2CaCu2O8+δ (BSCCO), Josephson vortices penetrate into block layers. Since the phase difference between two adjacent CuO2 layers which confine a Josephson vortex is $`\pi `$, the Josephson vortex is a good example of non-linear dynamics in superconductors. Dynamical properties of Josephson vortices are strongly affected by the Josephson plasma because the Josephson vortex and the Josephson plasma are solutions of the non-linear and linear part of the sine-Gordon equation . Furthermore, a strong emission of electromagnetic wave is expected when the velocity of the JV lattice matches the Josephson plasma wave velocity . Therefore, the Josephson plasma resonance in the presence of the Josephson vortices is very powerful tool to reveal dynamical nature of the Josephson vortices.

The Josephson plasma resonance (JPR) is known to manifest the existence of weakly coupled Josephson effects between CuO2 layers sandwitched by the non-superconducting or insulating block layers. In particular, in Bi2Sr2CaCu2O8+δ (BSCCO), which is a typical high-$`T_c`$ superconductor and possesses an intrinsic high anisotropy due to extremely weak Josephson coupling, the plasma frequency $`\omega _pc/\sqrt{ϵ}\lambda _c`$ can be as low as 100 GHz, where $`ϵ`$ and $`\lambda _c`$ are the dielectric constant and the $`c`$-axis penetration depth. Absence of the Landau damping of the quasiparticles in addition allows us to observe extremely sharp resonance at microwave frequencies using conventional high sensitive microwave techniques . The plasma frequency in a magnetic field parallel to the $`c`$-axis can be expressed as a function of field and temperature: $`\omega _p^2(H,T)=\omega _p^2\mathrm{cos}\phi _{l,l+1}(H,T)`$, where $`\phi _{l,l+1}`$ is the gauge invariant phase difference between adjacent $`l`$ and $`l+1`$-th layers and $`\mathrm{}`$ denotes spatial and thermal averages. This means that $`\omega _p(H,T)`$ can directly sense the interlayer phase coherence in a dynamical manner, which is severely influenced by the thermal fluctuation of vortices. Thus, $`\omega _p`$ in general should decrease with increasing field at a given temperature. Especially in the vortex liquid state, $`\omega _p`$ is expected to follow $`H^{1/2}`$ as reported previously .

In contrast to this, the Josephson plasma in fields parallel to the $`ab`$ plane is much more complicated because of the non-linear behavior of the Josephson vortex as mentioned above. In the case of single junctions, the excitation spectrum in the presence of Josephson vortices has been studied by Fetter and Stephen  for example. They derived an excitation spectrum with two branches: one is the gapless vortex branch attributed to the dc sliding mode of JVs, and the other is the plasma branch identical to the transverse Josephson plasma, which monotonically increases with increasing field.

In a layered system having a multi-stack of the intrinsic Josephson junctions like in BSCCO, however, physics becomes richer, because it allows us to observe the longitudinal plasma mode (Nambu-Goldstone mode) that the coherence length comparable to the interlayer distance, and the Josephson plasma can propagate not only as the transverse wave but also as the longitudinal wave or their mixtures. Furthermore, the 2D lattice of JVs may be realized and the properties may be strongly affected by the intrinsic pinning. Although much progress of the theories for the dynamical properties of JVs , experimental development has been hindered due partly to the experimental difficulties restricted by precise magnetic field alignment to the $`ab`$-plane of the crystal, which requires its extremely high quality. Concerning the Josephson plasma in parallel fields, systematic experimental study in order to understand whole spectrum of physical properties of the Josephson plasma has not been performed, although some preliminary work of JPR in parallel fields to the $`ab`$-plane were reported by Matsuda et al..

In this paper, we present our results of JPR measurements in parallel fields to the $`ab`$-plane as functions of microwave frequency, field intensity, and temperature in order to reveal the entire nature of the Josephson plasma in the presence of JVs. As a result, we obtained two resonance branches, which are separated by a finite gap in good accordance with the recent theoretical approaches. One observed at higher frequencies is attributed to the Josephson plasma coupled with the JV lattice, while the other lying at lower frequencies is interpreted due to the vortex oscillation mode hardened by residual pinning.

## II Experiments

We used an under-doped BSCCO single crystal grown by the modified traveling solvent floating zone method. The under-doped crystal was obtained after annealing at 550 C for 12 hours in a vacuum. The superconducting transition temperature $`T_c`$ of this crystal was determined to be 70.2 K by magnetization measurements with a SQUID magnetometer. JPR measurements have been performed in a microwave frequency range between 9.8 and 75 GHz with both reflective and transmission type of bridge balance circuits with rectangular TE102 mode cavity resonators. The sample was placed inside the cavity so that the oscillating electric field $`E_{\mathrm{rf}}`$ is exerted parallel to the $`c`$-axis to selectively excite the longitudinal plasma mode . A frequency-stabilized microwave was generated by the signal swept generator (HP83650B) or the Gunn oscillators, and the resonance was detected as changes of cavity impedance (Q-factor) either by field sweep or by temperature sweep. The magnetic field was applied by a split-pair of superconducting magnet, which can generate horizontal magnetic fields up to 8 T and the angle between field direction and crystal axis was changed by rotating the cavity resonator with respect to the magnet by a precision goniometer. The field direction exactly parallel to the $`ab`$-plane was determined by looking at the symmetry of the angular dependence of the resonance within an accuracy of 0.01 degree.

## III Results and Discussion

Figures 1 (a) and (b) show the resonance curves at 18.9 and 25.5 GHz, respectively, obtained by sweeping temperature. Clear two resonance lines without hysteresis at higher and lower temperatures, and a temperature gap in between the two resonances are found at both frequencies. As the field is increased, the resonance at higher temperatures once shifts to lower temperatures, then turns to shift to higher temperatures above 1 kOe. It is finally smeared out above 1.5 kOe, while the resonance at lower temperatures appears at around 1 kOe from lower temperature side and shifts to lower temperatures with little change its line-shape. Hereafter, we call these resonances at higher and lower temperatures HTB (higher temperature branch) and LTB (lower temperature branch), respectively. The fact that two lines are observed and the resonance temperature of HTB increases as the field is increased is peculiar in comparison with the case in fields parallel to the $`c`$ axis. It is noted that the LTB resonance obtained by field sweep was not considered here, because it showed large hysteresis and moreover the line-width was too wide due to the pinning effect.

As the frequency is increased, these two resonance lines shift to lower temperatures. Figure 2 represents the summary of resonance peaks as functions of field and temperature for various frequencies. Both HTB and LTB shift to lower temperatures and the gap grows rapidly as the frequency increases. LTB goes beyond our experimental temperature range above 30 GHz as shown in the upper panel of Fig. 1, whereas HTB is obtained at all frequencies. At a frequency below 60 GHz, HTB is observed at zero field at a finite temperature $`T_0`$, where the microwave frequency corresponds to the temperature dependent zero-field plasma frequency $`\omega _p(T)`$. Extrapolating $`\omega _p(T)`$ to $`T=0`$ K, the inherent plasma frequency $`\omega _p`$ can be estimated to be 56.8 GHz. The turning of HTB is more significant at higher frequencies, the resonance temperature above 3 kOe exceeds $`T_0`$ above 30 GHz, suggesting that the plasma frequency above 3 kOe may exceed $`\omega _p`$, whereas LTB just shifts to lower temperatures with increasing frequency. The resonance peak above $`\omega _p`$ is also plotted in the lower panel of Fig. 2. Similarly to HTB, the resonance line above $`\omega _p`$ shows no hysteresis in field sweep experiments even at low temperatures and peaks by sweeping field and temperature correspond well. Furthermore, the resonance shifts to higher temperatures as the field is increased, and is smeared out at higher fields. Thus, it is natural that this resonance is attributed to HTB observed below $`\omega _p`$.

In order to estimate in-plane field contribution to the Josephson plasma, the plasma frequency as a function of magnetic field is extracted from Fig. 2. In Fig. 3 the plasma frequencies $`\omega _p(H)`$ divided by the zero-field plasma frequencies $`\omega _p(T)`$ for various temperatures are shown. It is clearly seen that two resonance modes well-separated by a frequency gap exist: the higher one extracted from HTB shows upturn at about 1 kOe, then monotonically rises to higher frequencies even above $`\omega _p`$, while the other extracted from LTB shows a gradual decrease with increasing field. It is also found that both are well scaled with $`\omega _p(T)`$ except for the initial decrease of HTB below 1 kOe. Since LTB can be observed only in a finite field, the existence of JVs would be indispensable for LTB. At low fields below 0.5 kOe, the extrapolation of LTB tends to the zero frequency with decreasing field, suggesting that LTB may be a gapless mode at the zero-field limit. This gives a straightforward interpretation that LTB and HTB are the vortex and the plasma branches as suggested by Fetter and Stephen .

Very recently, Koshelev and Machida have developed formulations of layered superconductors in high parallel fields to the $`ab`$-plane when a uniform oscillating electric field was applied over the $`ab`$-plane parallel to the $`c`$-axis . Assuming that all block layers are fully occupied by JVs, which are forced to make a distorted triangular lattice, they calculated the microwave dissipation due to the Josephson plasma modified by oscillations of the JV lattice. In this case, the Josephson plasma wave propagates an oblique direction along the primitive reciprocal lattice vector of the JV lattice $`𝐪`$, which consists of the $`c`$-axis component $`q_c`$ fixed on $`\pi /s`$ due to the intrinsic pinning and the $`ab`$-plane component $`q_{ab}`$ proportional to magnetic field $`H`$ as $`2\pi sH/\mathrm{\Phi }_0`$. The plasma frequency is determined by the dispersion relation along $`𝐪`$ lying between the longitudinal and the transverse plasma modes. (Actually, it lies very close to the longitudinal mode because of $`q_{ab}q_c`$.) As a result, the peak frequency of the dissipation spectrum was obtained to be proportional to the magnetic field parallel to the layers $`H`$,

$$\omega _p(H)=\omega _p\frac{\pi H\gamma s^2}{\mathrm{\Phi }_0}.$$
(1)
Assuming $`\gamma =1070`$, which is reasonable for under-doped BSCCO samples, this equation gives an excellent agreement with the experimental data above 3 kOe at all temperatures as shown in Fig. 3, where we used $`s=15`$ Å. It denotes that the same equation has derived by Bulaeveskii et al. in Ref. . According to Ichioka, all block layers are occupied by JVs above $`H^{}=\mathrm{\Phi }_0/4.6\gamma s^2`$ , which is obtained to be 2.1 kOe for the above-mentioned $`\gamma `$. This is consistent with the fact that the linear field dependence of HTB is violated below 3 kOe as shown in Fig. 3. We therefore conclude that the linear increase of $`\omega _p(H)`$ of HTB is attributed to the increase of $`q_{ab}`$ of the JV lattice, while $`q_c`$ is fixed on $`\pi /s`$ because of the intrinsic pinning.

On the other hand, LTB is considered to originate from the oscillations of Josephson vortices, which would be hardened by residual pinning. Figure 4 represents results of numerical simulations on the microwave absorption of the single junction model in the presence of randomness of the $`c`$-axis critical current density. Introducing the randomness, which represents pinning center of vortices, a faint peak appears below $`\omega _p`$. With increasing field, the faint peak does not change while the higher frequency peak shifts to higher frequencies and its intensity becomes weaker. This feature is consistent with the experimental results in the following points: the peak of the LTB resonance does not shift to higher frequency and the HTB resonance is smeared as field is increased. However, the intensity of the low-lying resonance is much fainter than that we observed, the intensity of LTB is rather higher than the one of HTB.

Finally, we comment on the initial decrease of HTB resonance frequency in low fields in connection with the Josephson vortex state. In the low field region, the effect of the JV lattice is not significant yet, so that the thermal fluctuation of vortices may play an important contribution to suppress the Josephson coupling likewise the case of JPR in perpendicular fields. More detailed analysis of this feature is necessary.

## IV Conclusion

We have investigated the Josephson plasma mode in parallel fields to the $`ab`$-plane of BSCCO. We obtained two resonance modes with a temperature dependent gap in-between. Strong coupling between Josephson plasma and JV lattice are responsible for them. In high fields where all block layers are occupied with JVs, one lying higher frequency comes from the Josephson plasma oscillation modified by the JV lattice, while the other is attributed to the vortex oscillation hardened by pinning of Josephson vortex lattice.

## Acknowledgements

We would like to thank A.E. Koshelev, T. Koyama, S. Takahashi, and M. Tachiki for valuable discussions. This work has been supported by CREST, Japan Science and Technology Cooperation and is partly supported by the Research Project of The University of Tsukuba.