text:
# Effects of point defects on the phase diagram of vortex states in high-𝑇_c superconductors in 𝐵⃗∥𝑐 axis

\[

## Abstract

The phase diagram for the vortex states of high-$`T_\mathrm{c}`$ superconductors with point defects in $`\stackrel{}{B}c`$ axis is drawn by large-scale Monte Carlo simulations. The vortex slush (VS) phase is found between the vortex glass (VG) and vortex liquid (VL) phases. The first-order transition between this novel normal phase and the VL phase is characterized by a sharp jump of the density of dislocations. The first-order transition between the Bragg glass (BG) and VG or VS phases is also clarified. These two transitions are compared with the melting transition between the BG and VL phases.

\] Vortex states in high-$`T_\mathrm{c}`$ superconductors in $`\stackrel{}{B}c`$ axis have been intensively studied. Although the melting transition in pure systems has now been understood very well, experimental phase diagrams are more complicated owing to effects of impurities. In the present Letter, point defects are taken into account as the simplest impurity.

A schematic phase diagram of high-$`T_\mathrm{c}`$ vortex states with point defects is given in the inset of Fig. 1. These three phases have been observed in YBa2Cu3O7-δ (YBCO)  and Bi2Sr2CaCu2O8+y (BSCCO).  The Bragg glass (BG) phase  is characterized by the power-law decay of correlation functions of vortex positions  and the triangular Bragg pattern of the structure factor. The vortex glass (VG) phase was first defined on the basis of phase variables,  and can be defined alternatively on the basis of vortex positions.  Recent simulations including screening effects  suggest the absence of the phase-coherent VG. However, the positional VG is expected to be more stable than the phase-coherent VG.  Therefore, whether the VG exists as the thermodynamic phase or not is still unsettled.

The BG–VG phase boundary was studied analytically,  essentially based on the Lindemann criterion. Numerically, difference between the BG and VG phases was discussed,  and the overall phase diagram was obtained recently.  However, studies based on thermodynamic quantities are still lacking.

Quite recently, another phase has been observed experimentally. Nishizaki et al. have found a sharp jump of the resistivity  and the magnetization  in optimally-doped YBCO in the vortex liquid (VL) region. They pointed out that this might be the transition to the vortex slush (VS) phase  originally found in irradiated YBCO. The VS phase is characterized by evolution of a short-range order, and the boundary between the VS and VL phases terminates at a critical point. Similar anomalies were also observed in BSCCO.

In the present Letter, we compose the phase diagram in the pinning-strength ($`ϵ`$)–temperature ($`T`$) plane (Fig. 1) on the basis of thermodynamic quantities. The first-order melting transition occurs between the BG and VL phases as in pure systems.  A first-order transition line stretches from the melting line into the VL region dividing the VS and VL phases. The VG phase exists at much lower temperatures. The boundary between the BG and VG or VS phases is almost independent of temperature, and the phase transition is of first order. This phase diagram is quite similar to the $`B`$$`T`$ phase diagram observed experimentally.

Our formulation is based on the three-dimensional frustrated XY model on a simple cubic lattice,

$``$ $`=`$ $`{\displaystyle \underset{i,jab\mathrm{plane}}{}}J_{ij}\mathrm{cos}\left(\varphi _i\varphi _jA_{ij}\right)`$ (2)
$`{\displaystyle \frac{J}{\mathrm{\Gamma }^2}}{\displaystyle \underset{m,nc\mathrm{axis}}{}}\mathrm{cos}\left(\varphi _m\varphi _n\right),`$
$`A_{ij}`$ $`=`$ $`{\displaystyle \frac{2\pi }{\mathrm{\Phi }_0}}{\displaystyle _i^j}𝐀^{(2)}d𝐫^{(2)},\stackrel{}{B}=\stackrel{}{}\times \stackrel{}{A},`$ (3)
with the periodic boundary condition. Screening effects are not included in this model. A uniform magnetic field $`\stackrel{}{B}`$ is applied along the $`c`$ axis, and the averaged number of flux lines per plaquette is given by $`f=l^2B/\mathrm{\Phi }_0`$. Here $`l`$ stands for the grid spacing in the $`ab`$ plane. Point defects are introduced as the plaquettes which consist of four weaker couplings and are randomly distributed in the $`ab`$ plane with probability $`p`$. Couplings are given by $`J_{ij}=(1ϵ)J`$ ($`0<ϵ<1`$) on the point defects, and $`J_{ij}=J`$ elsewhere. Here we concentrate on the model with $`\mathrm{\Gamma }=20`$, $`f=1/25`$ and $`p=0.003`$. We do not consider the lower critical point  at present. The system size is $`L_x=L_y=50`$ and $`L_c=40`$, which is large enough to analyze the pure system ($`ϵ=0`$).

For each $`ϵ`$, Monte Carlo simulations are started from a very high temperature, and systems are gradually cooled down. After such annealing, further equilibration and measurement are performed at each temperature. Typical simulation times are $`412\times 10^7`$ Monte Carlo steps (MCS) for equilibration, and $`25\times 10^7`$ MCS for measurement. The present simulations are based on one sample. Since configurations of point defects are independent in different $`ab`$ planes, it is reasonable to expect that there exists a self-averaging effect. We calculate the internal energy $`e`$, the specific heat $`C`$, the helicity modulus along the $`c`$ axis $`\mathrm{{\rm Y}}_\mathrm{c}`$,  and the phase difference between the nearest-neighbor $`ab`$ planes $`\mathrm{cos}(\varphi _n\varphi _{n+1})`$, together with the ratio of entangled flux lines to total flux lines $`N_{\mathrm{ent}}/N_{\mathrm{flux}}`$,  the density of dislocations in the $`ab`$ plane $`\rho _\mathrm{d}`$, and the structure factor of flux lines in the $`ab`$ plane. The helicity modulus is proportional to the superfluid density, and is the order parameter of superconductivity. The inter-layer phase difference is related  to the frequency of the Josephson plasma resonance (JPR).

BG–VL and VS–VL phase transitions.—Temperature dependence of $`e`$, $`C`$, $`\mathrm{{\rm Y}}_c`$ and $`\mathrm{cos}(\varphi _n\varphi _{n+1})`$ at $`ϵ=0.05`$, $`0.07`$ and $`0.11`$ is displayed in Figs. 2(a)–(d). The BG–VL transition is observed at $`ϵ=0.05`$, and the VS–VL one at $`ϵ=0.07`$ and $`0.11`$. Jumps of $`e`$ and $`\mathrm{cos}(\varphi _n\varphi _{n+1})`$ and the $`\delta `$-function peak of $`C`$ are observed at the melting temperature $`T_\mathrm{m}`$ or the slush temperature $`T_{\mathrm{sl}}`$. Finite latent heats $`Q=\mathrm{\Delta }e`$ indicate that these phase transitions are of first order. The anomalies of the VS–VL transition at $`ϵ=0.07`$ are as large as those of the BG–VL transition at $`ϵ=0.05`$, and gradually lose sharpness as $`ϵ`$ is increased (see Figs. 2(c) and 2(d)). At $`ϵ=0.14`$, no anomalies are observed other than a small jump of $`\mathrm{cos}(\varphi _n\varphi _{n+1})`$. These properties can be explained by the existence of the critical point  which terminates the first-order VS–VL transition line around $`ϵ=0.14`$.

The most important difference between the BG–VL and VS–VL transitions is seen in $`\mathrm{{\rm Y}}_c`$. In the BG–VL transition, this quantity appears discontinuously at $`T_\mathrm{m}`$. On the other hand, $`\mathrm{{\rm Y}}_c`$ remains vanishing for $`T<T_{\mathrm{sl}}`$ in the VS–VL transition. In other words, the BG–VL phase transition is the superconducting–normal one, while the VS–VL transition occurs between two normal phases. The latter transition does not contradict the existence of the critical point. The proliferation of $`\mathrm{{\rm Y}}_c`$ at lower temperatures signals the phase transition to the VG phase, though error bars of $`\mathrm{{\rm Y}}_c`$ are fairly large owing to very large correlation time in the VG phase.

Direct observation of flux lines also reveals the difference between these two first-order transitions. Temperature dependence of $`N_{\mathrm{ent}}/N_{\mathrm{flux}}`$ and $`\rho _\mathrm{d}`$ at $`ϵ=0.05`$ and $`ϵ=0.11`$ is displayed in Figs. 3(a) and 3(b), respectively. $`N_{\mathrm{ent}}/N_{\mathrm{flux}}`$ shows a sharp jump at $`T_\mathrm{m}`$ (Fig. 3(a)) as in pure systems,  while its temperature dependence is continuous around $`T_{\mathrm{sl}}`$ (Fig. 3(b)). The quantity $`\rho _\mathrm{d}`$ shows sharp jumps both at $`T_\mathrm{m}`$ and $`T_{\mathrm{sl}}`$ as predicted by Kierfeld and Vinokur.  These properties are consistent with the structure factors shown in the same figures. A ring-like pattern is seen in the VL phase as in pure systems  both at $`ϵ=0.05`$ and $`0.11`$. The clear triangular Bragg pattern for $`T<T_\mathrm{m}`$ at $`ϵ=0.05`$ represents the formation of a hexatic quasi long-range order. The obscure Bragg pattern with a $`6`$-fold symmetry for $`T<T_{\mathrm{sl}}`$ at $`ϵ=0.11`$ might stand for domains of a short-range hexatic order in the $`ab`$ plane divided by dislocations.

BG–VG and BG–VS phase transitions.—Pinning-strength dependence of $`e`$ and $`\mathrm{cos}(\varphi _n\varphi _{n+1})`$ is shown in Figs. 4(a) and 4(b) for $`T=0.06`$, $`0.07`$ and $`0.08J/k_\mathrm{B}`$. Finite latent heats $`Q`$ indicate that these phase transitions are of first order. The quantity $`\mathrm{cos}(\varphi _n\varphi _{n+1})`$ jumps simultaneously, which is consistent with Gaifullin et al.’s JPR experiment of BSCCO.  Since the anomalies are always observed between $`ϵ=0.065`$ and $`0.070`$, the phase boundary is almost independent of temperature as shown in Fig. 1. We also find a sudden jump of $`N_{\mathrm{ent}}/N_{\mathrm{flux}}`$ on this phase boundary as predicted theoretically.

The jumps of $`e`$ and $`\mathrm{cos}(\varphi _n\varphi _{n+1})`$ (abbreviated as $`\mathrm{\Delta }\mathrm{cos}`$ from now on) on this almost-flat phase boundary are about one order smaller than those on the melting line (see Figs. 2(b) and 4(b)). The ratio of the jump of the Josephson energy $`\mathrm{\Delta }e_\mathrm{J}=(J/\mathrm{\Gamma }^2)\mathrm{\Delta }\mathrm{cos}`$ to $`Q`$ reveals the difference between these two transitions clearly. As in extremely anisotropic pure systems,  this ratio is about one half on the melting line, e.g. $`\mathrm{\Delta }e_\mathrm{J}5.5\times 10^4J`$ and $`Q1.1\times 10^3J`$ at $`ϵ=0.05`$ (see Figs. 2(a) and 2(b)). On the other hand, $`\mathrm{\Delta }e_\mathrm{J}/Q`$ is about unity on this almost-flat phase boundary, e.g. $`Q\mathrm{\Delta }e_\mathrm{J}9.0\times 10^5J`$ at $`T=0.07J/k_\mathrm{B}`$ (see Figs. 4(a) and 4(b)). The latter result shows a sharp contrast to the JPR experiment of BSCCO, where $`\mathrm{\Delta }e_\mathrm{J}/Q1`$ seems to be satisfied.

Discussions.—On the melting line of pure systems, $`\mathrm{\Delta }\mathrm{cos}`$ is proportional to $`\mathrm{\Gamma }^2f`$ and gradually approaches a saturated value,  $`\mathrm{\Delta }\mathrm{cos}0.3`$.  When the anisotropy is as small as that of YBCO ($`\mathrm{\Gamma }78`$), $`\mathrm{\Delta }\mathrm{cos}`$ is small both on the melting line and on the BG–VG phase boundary. In the present system ($`\mathrm{\Gamma }=20`$), $`\mathrm{\Delta }\mathrm{cos}`$ ($`0.22`$ at $`ϵ=0.05`$) is as large as the saturated value on the melting line, while it is small on the BG–VG/VS phase boundary. Therefore, a jump in $`\mathrm{cos}(\varphi _n\varphi _{n+1})`$ inevitably occurs in the VL region, which results in the VS–VL phase transition. On the other hand, when the anisotropy is as large as that of BSCCO ($`\mathrm{\Gamma }150`$), $`\mathrm{\Delta }\mathrm{cos}`$ has reached the saturated value both on the melting line and on the BG–VG phase boundary as shown experimentally.  In such a case, it might be difficult to observe a jump of $`\mathrm{cos}(\varphi _n\varphi _{n+1})`$ outside of the BG phase by the JPR.

Finally, the present results summarized in Fig. 1 are compared with theoretical studies related to the VS phase in literature. When Worthington et al.  proposed the VS–VL transition line as a reminiscent of the melting line in pure systems, the BG phase was out of the scope. Ikeda  derived a phase diagram consisting of the VG, VS and VL phases, and argued that the BG phase and the VS phase cannot coexist. Quite recently, he modified his argument and proposed a possible phase diagram including both the BG and VS phases.  However, he simply assumed the existence of the BG phase in this article. Kierfeld and Vinokur  obtained a phase diagram consisting of the BG phase and a first-order transition line stretching from the melting line into the VL region with a critical point. Although they interpreted this transition line as the VG–VL one, it turns out to correspond to the VS–VL one as shown in the present study. Reichhardt et al.  numerically found a window-glass-like region with diverging time scales and a finite correlation length, which might also correspond to the VS phase. We believe that the present simulations are the first theoretical derivation of the phase diagram consisting of four phases (Fig. 1) based on thermodynamic quantities and from a single model. Although finite-size analyses of $`\mathrm{{\rm Y}}_c`$  have not been performed at present, our system size, $`L_c=40`$ for $`\mathrm{\Gamma }=20`$ and $`f=1/25`$, is already very large. According to the anisotropy scaling,  this size corresponds to $`L_c113`$ for the parameters in Ref. , $`\mathrm{\Gamma }=\sqrt{40}`$ and $`f=1/20`$.

We would like to thank Y. Matsuda and T. Nishizaki for communications, and A. Tanaka for comments. One of us (Y. N.) also thanks K. Kadowaki and R. Ikeda for sending their preprints. Numerical calculations were performed on Numerical Materials Simulator (NEC SX-5) at National Research Institute for Metals, Japan.


format: text

text:
# Critical behavior of the long-range Ising chain from the largest-cluster probability distribution

\[

## Abstract

Monte Carlo simulations of the 1D Ising model with ferromagnetic interactions decaying with distance $`r`$ as $`1/r^{1+\sigma }`$ are performed by applying the Swendsen-Wang cluster algorithm with cumulative probabilities. The critical behavior in the non-classical critical regime corresponding to $`0.5<\sigma <1`$ is derived from the finite-size scaling analysis of the largest cluster.

PACS numbers: 05.50.+q, 64.60.Cn

\]

The calculation of the critical behavior of the Ising model with long-range (LR) interactions decaying with distance with a power law $`1/r^{d+\sigma }`$ is not an easy task, even in one dimension, where phase transition at finite temperature occurs for $`0<\sigma 1`$ . Exact analytical expressions for the critical exponents may be derived for $`\sigma 0.5`$ , corresponding to the classical regime, while for $`\sigma >0.5`$ only the approximate results exist. Various analytical and numerical approaches have been applied, from direct numerical calculations on finite chains, to several approaches based on the renormalization group (RG) and scaling. The nonlocal character of interactions reduces the efficiency of most of the standard approaches so that the values of critical exponents obtained by all these methods differ considerably.

In numerous cases of phase transitions in systems with short-range interactions, a useful complementary tool for obtaining both qualitative and quantitative results is provided by Monte Carlo (MC) simulations in combination with finite-size scaling.

Such systematic studies were lacking for the LR models until recently. Namely, when applied to models with LR interactions, the standard MC approaches based either on Metropolis or on various cluster algorithms are particularly time consuming, since the number of operations per spin-flip is proportional there to the size of the system. Recently, this problem was successfully resolved by Luijten and Blöte  who used the cumulative probabilities within the Wolff cluster algorithm, which they applied to the Ising and similar models reducing the computing time by several orders of magnitude.

Their very exhaustive studies concentrate on questions related to the mean-field (MF) regime, while very little or no interest has been dedicated yet to the regime of the non-classical critical behavior corresponding to $`\sigma >0.5`$.

The purpose of this work is to extend the MC studies of the critical behavior of the LR Ising model to the non-classical regime. At the same time this is a suitable example to examine the efficiency of using only the cluster statistics in deriving the critical properties of the LR model.

The 1D Ising model with LR interactions written in form of a special case ($`q=2`$) of the Potts model is described by the Hamiltonian

$$H=\underset{i<j}{}\frac{J}{|ij|^{1+\sigma }}\delta _{s_i,s_j},$$
(1)
where $`J>0`$, $`s_i`$ is a two-state Potts variable at the site $`i`$, $`\delta `$ is the Kronecker symbol and the summation is over all pairs of the system. By the substitution $`\delta _{s_i,s_j}=(S_iS_j+1)/2`$, and $`J_I=J/2`$, where $`S_i=\pm 1`$ and $`J_I`$ denote the Ising spins and the interaction constant, respectively, one recovers the standard definition of the Ising model. In the mean-field regime, $`\sigma 0.5`$, the critical exponents have classical values ($`\nu =1/\sigma `$, $`\eta =2\sigma `$). We focus here on the region $`0.5<\sigma <1`$, where the critical exponents are non-trivial and known only approximately.

When the shortcut using the cumulative probabilities is implemented in the MC simulations, the full extent of the numerical advantages is achieved when only the distribution of spins has to be calculated, while for the energy sampling the reduction in the CPU-time is less important. Recently the calculation by energy sampling was improved by calculating the energy in the momentum space and applying the fast Fourier transform. Our intent here is to avoid the energy sampling by making a more complete study of the cluster statistics. To this purpose we use the Swendsen-Wang cluster algorithm. The implementation of cumulative probabilities to the Swendsen-Wang cluster algorithm is straightforward. Each step of this iterative procedure consists in identifying all the clusters in a given spin configuration of the system following the rule that two particles belong to the same cluster with the probability $`p_{ij}=(1exp[J(i,j)/T])\delta _{s_i,s_j}`$, and then flipping all the clusters randomly. The cumulative probabilities are applied at the point of identifying the individual clusters along the same lines as was done for the Wolf single-cluster algorithm, and this reduces the number of required operations per a single spin-flip by a factor of system size. Like the Wolf algorithm, the Swendsen-Wang algorithm suppresses the slowing-down at criticality, and though it might be somewhat more costly regarding CPU time, it gives a more complete insight into the cluster statistics and the related probability distributions, such as the distribution of cluster sizes or that of the largest cluster.

The probability distributions related to clusters are mostly used in the description of critical properties of geometrical transitions, like the percolation. As it follows from the graph expansion by Kasteleyn and Fortuin , the very basis of the Swendsen-Wang algorithm, this cluster statistics may be related to the thermodynamic quantities of the thermal transitions as well. The only important difference is that in percolation one deals with the simple geometrical clusters, while in thermal transitions the clusters formed by the active bonds (in terms of the cluster algorithm) are considered. In systems with LR interactions, such as the present model, this difference is obvious.

The cluster-related distributions were much less exploited in the case of thermal transitions, although some detailed studies exist (see e.g. refs ) for models with short-range interactions. More recently, there has been renewed interest in this subject within a somewhat different context .

By performing $`10^5`$ MC iterations, which is sufficient for most of data presented here, the chains of sizes up to a few tens of thousands of sites could be simulated without requiring too much numerical effort on a modest workstation.

We limit our analysis here to the size distribution of the largest cluster, the quantity that was left out of reach within the earlier single-cluster approach.

The probability that the largest cluster will be of the size $`l`$ is defined as

$$P_{max}(l)=\frac{1}{n_{MCS}}N_{max}(l),$$
(2)
where $`N_{max}(l)`$ is the total number of occurrences of the largest cluster of size $`l`$ during $`n_{MCS}`$ MC swaps of the system.

Several thermodynamic quantities of interest may be expressed through the corresponding moments defined as

$$<l^k>=\underset{l}{}l^kP_{max}(l).$$
(3)
The first moment

$$<l>=\underset{l}{}lP_{max}(l)$$
(4)
gives the mean size of the largest cluster. Below the critical temperature it describes the order parameter, which is usually defined as the average of the largest component, i.e. $`M=(q<max\{m_\alpha \}>1)/(q1)`$, where $`m_\alpha `$ denotes the fraction of spins of the system in the state $`\alpha `$.

A quantity of interest to consider is the Binder’s fourth-order cumulant  ratio for the largest cluster,

$$R_L(T)=\frac{<l^4>}{<l^2>^2},$$
(5)
defined in the same way as for the standard order parameter. Differences between the two corresponding distribution functions, however, result in different shapes of the same ratio defined for the two quantities. As may be seen in Fig. 1, the ratio $`R_L(T)`$ at high temperatures and in the thermodynamic limit tends to 1 and not to 3, while the corresponding probability distribution is expected to be non-Gaussian, as found already in the case of percolation  and the short-range Ising model . However, as well as in the mentioned cases, these differences should not change the basic scaling properties in the critical regime. At low temperatures and around $`T_c`$ the behavior of the two ratios coincides. Like in the case of the standard order parameter, there is a common crossing point of the curves $`R_L(T)`$ of different sizes where the ratio (5) almost does not depend on $`L`$ and which can be identified as the critical temperature. In Fig. 1, $`R_L(T)`$ is shown in the case $`\sigma =0.9`$ and for sizes varying from 1000 to 20000.

By taking into account sizes within the range $`1000L20000`$, the determination of this crossing point can be done with precision of 4 digits, which gives the estimated error of $`0.2\%`$.

In table I are presented the critical temperatures for several values of $`\sigma `$ in the non-classical regime $`0.5<\sigma <1`$ incremented by $`0.1`$, obtained as common crossings of $`R_L(T)`$. For comparison, we quote the earlier results obtained by numerical calculations , coherent anomaly approach , finite-range scaling (FRS) , and real-space RG . When quoting ref.  we mention only one of the two given sets of results for $`T_c`$, which fits better to the later work  performed with a significantly improved precision, but, unfortunately, available only for two values of the range considered here.

Near $`T_c`$ (i.e. in the regime $`L<\xi `$), the moments (3) and their ratio (5) should obey finite-size scaling. Thus

$$R_L(T)=L^xf(L^{1/\nu }\tau ),$$
(6)
where $`\tau =(TT_c)/T_c`$ is the reduced temperature. Since $`R_{\mathrm{}}(T_c)`$ is finite, $`x=0`$.

It turns out that the critical exponent $`\nu `$ may be estimated with a rather good precision by a simple tuning of the unknown exponent $`1/\nu `$ in the Eq. (6) until the curves collapse in the vicinity of the crossing point.

In Fig. 2 we present the curves from Fig. 1, collapsed after taking $`1/\nu =0.40`$.

The results for $`\nu ^1`$ obtained by collapsing fits for other considered values in the interval $`0.5<\sigma <1`$ are given in table II.

In order to reduce the degree of arbitrariness of such fits we require the curves to be indistinguishable in the region of temperatures around $`T_c`$, up to the temperature of the maximum of $`R_L(T)`$, $`T_{max}`$. This temperature may be related to the finite-size crossover temperature up to which the finite-size scaling equation (6) should hold. In Fig. 1 one can observe that the maximum of $`R_L(T)`$ shrinks and shifts towards $`T_c`$ with increasing L. The rough fit of the difference $`T_{max}T_c`$ to the power-law form $`(1/L)^\varphi `$ gives the values of $`\varphi `$ which roughly correspond to $`1/\nu `$ for a given $`\sigma `$, as expected for the finite-size crossover exponent . The interval of tuning in which the curves remain indistinguishable determines the error margins of the obtained values for the exponent $`\nu ^1`$, which increase with decreasing $`\sigma `$ but do not exceed 4% in the considered range of $`\sigma `$’s.

Alternative calculations of the exponent, e.g. by using the derivatives of moments, turn out to be less advantageous. (The purely numerical derivative would require a much higher numerical precision and consequently far more extensive runs; the derivative obtained by histogram interpolations or by mixed moments would require an energy sampling and an increased numerical precision in addition.) We also mention at this point that the calculations of the Binder’s cumulant ratio may equally be performed for the standard order parameter using the Wolf cluster algorithm. One advantage of the approach which uses the Swendsen-Wang cluster algorithm is that it has similar precision and efficiency both in low- and high-temperature regimes, which improves the scaling fit.

The degree of precision of the obtained data permits the comparison to the results by other approximate methods. These results, which are quite few, are quoted in table II in chronological order and correspond to numerical calculations on finite chains using Padé series extrapolations , $`ϵ`$-expansion near $`2\sigma =d`$ , the coherent anomaly method , finite-range scaling , and real-space RG . In refs  and  only the numerical data for the exponents $`\gamma `$ and $`\beta `$ are available. The exponent $`\nu `$ may be expressed in two ways, which, due to the approximate character of the results, do not give identical values. By ”a” are labelled the values obtained from the scaling relation $`\nu =\gamma /(2\eta )`$ and the exact analytical expression $`\eta =2\sigma `$ expected to hold also in the whole non-classical regime. By ”b” are labelled the values obtained from the scaling relation involving both exponents, $`\nu =\gamma +2\beta `$. The quoted results of ref.  are the ones obtained by the expansion in $`\mathrm{\Delta }\sigma =\sigma 0.5`$ up to the second order. The expansion in $`(1\sigma )`$  is not quoted in the table. It holds very close to $`\sigma =1`$ and we may quote only the value for $`\sigma =0.9`$ which is equal to $`\nu ^1=0.447`$. The results of ref.  are also expected to be more reliable closer to $`\sigma =0.5`$. Other approaches cover the whole range of $`\sigma `$ evenly, but their results are rather different and their accuracy as well. Our MC results agree the best with the FRS results.

In ref. , which deals with the Potts model with arbitrary number of states q, we found intriguing the conjecture that the exponent $`\nu `$ would remain the same for other values of $`q`$ (when the transition is of the second-order), which is in contrast to earlier results of FRS . A MC approach can thus be used to provide an independent calculation. (It also has an advantage over both used approaches that it is able to distinguish the first- from the second-order transition). Our first calculations , performed along the same lines as above, give for the case $`q=3`$, $`\sigma =0.9`$ the value $`1/\nu =0.48\pm 0.01`$, which would rather be in favor of the FRS results.

In conclusion, we have applied the MC simulations to the two-state Potts model with LR interactions on chains up to 20000 sites by using the Swendsen-Wang algorithm to calculate the cluster statistics, in particular the size distribution of the largest cluster. The approach was used as an alternative calculation of the critical behavior in this model in the non-classical regime $`0.5<\sigma <1`$. We have shown that the scaling analysis of the Binders 4th-order cumulant ratio of the largest-cluster size gives the critical temperature and the critical exponent $`\nu `$ with reasonable accuracy (better then 0.2% and 4%, respectively). Although the overlapping fit is not expected to be very precise, nor were the simulations pushed to their extreme, the estimated errors are significantly smaller than the differences between the results obtained by different approaches. The obtained results, presented in tables I and II, are closest to the values obtained earlier by the FRS. While the present work deals only with the Binder’s 4th-order cumulant ratio, a wider analysis, which includes other quantities that may be derived from the cluster statistics, should be performed in future. We believe that the approach through cluster statistics might be useful to study other cases of the Potts model and other models with discrete symmetries, where in the case of LR interactions, the alternative methods are very restricted.


text:
# Waves attractors in rotating fluids: a paradigm for ill-posed Cauchy problems

\[

## Abstract

In the limit of low viscosity, we show that the amplitude of the modes of oscillation of a rotating fluid, namely inertial modes, concentrate along an attractor formed by a periodic orbit of characteristics of the underlying hyperbolic Poincaré equation. The dynamics of characteristics is used to elaborate a scenario for the asymptotic behaviour of the eigenmodes and eigenspectrum in the physically relevant régime of very low viscosities which are out of reach numerically. This problem offers a canonical ill-posed Cauchy problem which has applications in other fields.

\]

Rotating fluids encompass all fluids whose motions are dominated by the Coriolis force. These flows play an important role in astrophysics or geophysics where the large size of the bodies makes the Coriolis force a prominent force. Some engineering problems like the stability of artificial satellites also require the study of rotating fluids because of their liquid-filled tanks . This latter problem is related to the existence of waves specific to rotating fluids, namely inertial waves, which easily resonate. These waves play also an important part in the oscillation properties of large bodies like the atmosphere, the oceans, the liquid core of the Earth, rapidly rotating stars or neutron stars. As such, they have been considered since the work of Poincaré on the stability of figures of equilibrium of rotating masses . Pressure perturbations of inertial modes for inviscid fluids obey the Poincaré equation (PE) (christened by Cartan) which reads $`\mathrm{\Delta }P(2\mathrm{\Omega }/\omega )^2^2P/z^2=0`$ where $`\mathrm{\Omega }\stackrel{}{e}_z`$ is the angular velocity of the fluid and $`\omega `$ is the frequency of the oscillation. Since $`\omega <2\mathrm{\Omega }`$ , the PE is hyperbolic (energy propagates along characteristics) and since its solutions must meet boundary conditions, the problem is ill-posed mathematically. Although some smooth solutions exist (for instance for a fluid contained in a full sphere or a cylinder), one should expect singular solutions in the general case. These latter solutions have been made explicit only recently thanks to numerical simulations which include viscosity to regularize the singularities and let this parameter be very small as in real systems.

In this letter we wish to present a scenario, based on analytical and numerical results, for the asymptotic behaviour of inertial modes at small viscosities. We use the case of a spherical shell as a container, which is relevant for astrophysical or geophysical problems, but it will be clear that this case is general. We will only sketch the main results, more details can be found in . While the fluid mechanical problem is of much interest by itself, it opens new perspectives in the theory of Partial Differential Equations (PDE) and also offers a toy model for some (very involved) problems of General Relativity which we shall present briefly.

The model we use is a spherical shell whose inner radius is $`\eta R`$ and outer radius $`R`$ ($`\eta <1`$). The fluid is assumed incompressible with a kinematic viscosity $`\nu `$. We write the linearized equations of motion for small amplitude perturbations for the velocity $`\stackrel{}{u}`$ in a frame corotating with the fluid; momentum and mass conservation imply:

$$\frac{\stackrel{}{u}}{t}+\stackrel{}{e}_z\times \stackrel{}{u}=\stackrel{}{}p+E\mathrm{\Delta }\stackrel{}{u},\stackrel{}{}\stackrel{}{u}=0$$
(1)
when dimensionless variables are used; $`(2\mathrm{\Omega })^1`$ is the time scale and $`E=\nu /2\mathrm{\Omega }R^2`$ the Ekman number. When $`E`$ is set to zero and $`\stackrel{}{u}`$ is eliminated, one obtains the Poincaré equation. In nature $`E1`$ and one is tempted to use boundary layer theory and singular perturbations to solve (1). However, this is feasible only when regular solutions exist for $`E=0`$; this is the case when the container is a full sphere but not when the container is a spherical shell. Indeed, numerical solutions of the eigenvalue problem issued from (1), where solutions of the form $`\stackrel{}{u}(\stackrel{}{r})e^{\lambda t}`$ are searched for (with $`1\omega =Im(\lambda )1`$), yield eigenmodes of the kind shown in Fig. 1. In this figure we see that the amplitude of the mode is all concentrated along a periodic orbit of characteristics of the PE; we found this property to be quite general, after extensive numerical exploration of least-damped modes of (1) , and will now explain its origin and consequences on the asymptotic spectrum of inertial modes. For this purpose we will use axisymmetric modes since the azimuthal dependence of solutions can always be separated out because of the axial symmetry of the problem.

For understanding the concentration of kinetic energy along a periodic orbit of characteristics, it is necessary to consider in some details the dynamics of these lines. Characteristics of PE are, in a meridional plane, straight lines making the angle $`\mathrm{arcsin}\omega `$ with the rotation axis. A numerical calculation of their trajectories shows that they generally converge towards a periodic orbit which we call, after, an attractor. The periodic orbit of Fig. 1 is one example of such an attractor.

The Lyapunov exponent (LE) of a trajectory, defined by $`\mathrm{\Lambda }=lim_N\mathrm{}\frac{1}{N}_{n=1}^N\mathrm{ln}\left|\frac{d\varphi _{n+1}}{d\varphi _n}\right|`$ ($`\varphi _n`$ is the latitude of the nth reflection point), describes how fast characteristics are attracted or repelled. Its computation as a function of frequency shows that attractors ($`\mathrm{\Lambda }<0`$) are ubiquitous in frequency space (see Fig. 2). Their existence shows that the dynamical system described by the characteristics is not hamiltonian; the “dissipation” is purely geometrical and is due to the fact that, unlike billiards, the reflection on boundaries is not specular but conserves the angle with the rotation axis. In fact, the dynamics of rays is a one-to-one one-dimensional map (from the outer boundary to itself), piecewise smooth, but with a finite number (twelve) of discontinuities. This kind of map has not been studied in the literature of dynamical systems, perhaps because it does not produce chaos because of its invertibility. Iterations of such a map generate fixed points which either correspond to attractors or to some neutral periodic orbits. Indeed, if $`\eta =0`$ (i.e. the sphere is full), all orbits such that $`\omega =\mathrm{sin}(p\pi /2q)`$ with $`(p,q)\mathrm{I}\mathrm{N}^2`$, are neutral and periodic while those such that $`\omega =\mathrm{sin}(r\pi )`$, $`r`$ being irrational, are neutral ergodic (quasiperiodic). When $`\eta `$ is non-zero only a finite number of such neutral periodic orbits subsist; for instance, if $`\eta =0.35`$ which is the aspect ratio of the Earth’s liquid core, $`q=1,2,3,4`$ are the only possibilities. Interestingly, we face here a situation which is just the opposite of the one described by the KAM theorem in Hamiltonian systems: when the full sphere is perturbed by the introduction of an inner sphere, all ergodic orbits are instantaneously destroyed while the longer periodic orbits survive the smaller the denominator $`q`$ is.

Apart from these isolated frequencies which become rarer and rarer as $`\eta `$ increases, generic trajectories are in the basin of attraction of attractors. We were able to show  that the number of attractors at a given frequency is finite. The inset of Fig. 2 shows the typical case where an attractor exists in a frequency band $`[\omega _1,\omega _2]`$ with $`\mathrm{\Lambda }(\omega _1)=0`$, $`\mathrm{\Lambda }(\omega _2)=\mathrm{}`$ and $`\omega _2\omega _11/N^2`$ where $`N`$ is the length of the attractor defined as its number of reflection points. Near $`\omega _1`$, $`\mathrm{\Lambda }\sqrt{\omega \omega _1}`$ and near $`\omega _2`$, $`\mathrm{\Lambda }(\omega )\frac{1}{N}\mathrm{ln}N(\omega \omega _2)`$. The latter implies that long attractors have small LE in a large fraction of $`[\omega _1,\omega _2]`$ (all these results are shown in ).

The existence of attractors for characteristics implies that solutions of the inviscid problem (i.e. of PE) are singular. This property can be made explicit in the simplified case of a 2D problem. Indeed, in this case the PE may be written $`^2P/u_+u_{}=0`$ using characteristics coordinates; solutions may be constructed explicitly from an arbitrary function but, as shown in , regular eigenmodes exist only when neutral periodic orbits exist and eigenvalues are infinitely degenerate. When attractors are present, the scale of variations of the pressure vanishes on the attractors while its amplitude remains constant. As velocity depends on the pressure gradient, it diverges on the attractor; this divergence is like the inverse of the distance to the attractor which makes the velocity field not square integrable. This result seems to be valid also in 3D .

We therefore understand why solutions of (1) look like Fig. 1: the inviscid part of the operator focuses energy of the modes thanks to the action of the mapping made by characteristics while viscosity opposes to this action via diffusion. The resulting picture of Fig. 1 therefore comes from a balance between inviscid terms and viscous ones; let us make this more quantitative.

For this purpose we first observe that the patterns drawn by the kinetic energy of the mode in Fig. 1 is in fact a shear layer whose width scale with $`E^\sigma `$ and $`\sigma 1/4`$. Such a scaling is observed numerically and seems generic ; it implies that the damping rate of such modes scales like $`E^{1/2}`$ as clearly shown in Fig. 3. Now we may consider a wave packet travelling around an attractor in a slightly viscous fluid. The above mentioned balance, when applied to both the width and the amplitude of the packet, leads to a relation between the LE and the Ekman number such as $`\mathrm{\Lambda }E^{13\sigma }`$ with $`\sigma <1/3`$ for an eigenmode of the viscous problem. We see that the constraint $`\sigma <1/3`$ is met by actual shear layers. It therefore turns out that frequencies of eigenmodes of the viscous problem are such that $`\mathrm{\Lambda }0`$ when $`E0`$ which means that they will gather around the roots of the equation $`\mathrm{\Lambda }(\omega )=0`$.

The above result shows the importance of the scaling verified by shear layers. A boundary layer analysis reveals that these shear layers are in fact nested layers which consist of an inner $`\sigma =1/3`$-layer surrounded by a thicker layer. The inner 1/3-layer can be fully explicited. Using coordinates along the shear layers ($`x`$) and perpendicular to it ($`y`$), we find that the $`\phi `$-component of the velocity verifies $`\frac{^3u_\phi }{Y^3}=i\frac{u_\phi }{q}`$, with $`Y=y/E^{1/3}`$ and $`q=x/\sqrt{1\omega ^2}`$ which is also the equation verified by the stream function in a steady shear layer of a rotating fluid. Solutions which vanish in $`Y=\pm \mathrm{}`$ are self-similar and of the form $`u_\phi =q^\alpha H_\alpha \left(Y/q^{1/3}\right)`$ with $`H_\alpha (t)=_0^{\mathrm{}}e^{ipt}e^{p^3}p^{3\alpha 1}𝑑p`$. Besides, $`\alpha =\frac{1}{3}`$ is the only admissible value to ensure a coherent evolution of the width and amplitudes after reflection on a boundary.

We are now in position to propose a scenario for the asymptotic behaviour of inertial modes when the viscosity vanishes. Eigenfunctions reduce to nested shear layers concentrated along attractors while the associated eigenvalues converge toward the frequency $`\omega _i`$ such that $`\mathrm{\Lambda }(\omega _i)=0`$ for the associated attractor. Furthermore, we can constrain this convergence of eigenfrequencies; indeed, since $`\mathrm{\Lambda }\sqrt{\omega \omega _i}`$, one finds that $`\omega =\omega _i+aE^{26\sigma }+\mathrm{}`$ and $`\tau =Re(\lambda )=bE^{12\sigma }`$ when $`E0`$; Fig. 3 shows that this law agrees well with the numerical results, in the case shown, with $`\sigma =1/4`$.

In addition, we noticed earlier that for a finite number of $`\omega `$ such that $`\omega =\mathrm{sin}(p\pi /2q)`$ all orbits of characteristics are periodic; this implies that in the vicinity of these frequencies very long attractors with very small average LE accumulate as shown by Fig. 4; therefore, these frequencies will be accumulation points of the asymptotic spectrum. Moreover, around these frequencies eigenmodes are weakly damped. On the contrary, modes whose frequency is in the frequency band of short attractors (like the one of Fig. 1) are more strongly damped. It therefore turns out that the LE curve in Fig. 2 will strongly constrain the distribution of least-damped modes in the complex plane at finite viscosities: such modes will avoid the large frequency bands of short-period attractors and concentrate around frequencies where $`\mathrm{\Lambda }(\omega )=0`$ especially those with $`\omega =\mathrm{sin}(p\pi /2q)`$.

This general evolution of the spectrum is well illustrated in Fig. 5. Here, the least-damped eigenvalues have been computed for $`E=10^8`$. We clearly see frequency bands of attractors avoided by weakly damped eigenvalues but see the gathering of these eigenvalues around $`\mathrm{sin}(\pi /4)`$ and, but less conspicuously, around $`\mathrm{sin}(\pi /6)`$.

To complete the picture, we need now mentioning that a few regular modes survive among all these singularities; such modes are purely toroidal modes or r-modes  which are non-axisymmetric. They avoid the constraint of characteristics for their velocity field has no radial component; this property makes their characteristics independent of frequency (they are circles and vertical straight lines) and authorizes smooth solution at zero viscosity. The associated eigenvalues $`\omega =1/(m+1),m\mathrm{I}\mathrm{N}^{}`$ seem to be the only eigenvalues of the Poincaré operator in a spherical shell.

Ending this letter, it is worth emphasizing the role of the geometrical approach allowed by the dynamics of characteristics, for describing the asymptotic properties of inertial modes; in the domain of very low Ekman numbers ($`10^{10}10^{20}`$), typical of astrophysical or geophysical fluids, these modes are out of reach numerically.

The foregoing presentation shows that inertial modes display a very rich dynamical behavior which comes from the ill-posedness of the underlying inviscid problem. Here we discussed the case of the spherical shell, but our results are general and can be extended to any container; this is important since natural containers are usually not perfect geometrical objects. Hence, fortunately, a curve like Fig. 2 is structurally stable (see our discussion relative to the core of the Earth in ).

We note that the relevance of attractors has also been shown experimentally in stratified fluids. Some configurations of conducting fluids bathed by a magnetic field, obeying the PE, will also display attractors. These properties are in fact very general and extend to mixed-type PDE as illustrated by the case of gravito-inertial modes. We think that similar results should hold for systems which are solutions of PDE of hyperbolic or mixed type meeting boundary conditions. As an example, our results may have applications in General Relativity and the problem of “closed timelike curves” (CTC), that is the problem of the existence of physical systems which permit travels backward in time. Such systems like wormholes have been studied by various authors; they set many problems among which that of causality. Such a problem is also at the origin of the ill-posedness of the Poincaré problem and we showed that it leads to many kinds of singularities.

We therefore see that inertial oscillations of a fluid inside a container offers a paradigm which may guide our intuition for problems in other fields of physics which are also ill-posed Cauchy problems.

We would like to thank Boris Dintrans and Leo Maas for very helpful discussions. Part of the calculations have been carried out on the Cray C98 of IDRIS at Orsay and on the CalMip machine of CICT in Toulouse which are gratefully acknowledged.


text:
# New access to very weak interactions in molecules

## I Introduction

Many symmetrical molecules exist in Nature only in the form of nuclear spin isomers . These isomers differ by symmetry of nuclear spin wave function and, consequently, by symmetry of molecular spatial wave function. Relaxation between different spin states (spin conversion) is extremely slow process if it is compared with other gas kinetic rates, e.g., vibrational relaxation. It makes nuclear spin isomers unique objects with many potential applications.

It was shown in a recent paper  that external radiation can influence conversion of spin isomers both through level populations and through optically induced coherences. Purpose of the paper  was to investigate main features of the coherent control of isomer enrichment and conversion. In order to achieve the goal, the process was considered in a simplest arrangement in which microwave radiation excited molecular rotational transition.

In the present paper we will consider the process in more complicated arrangement which promises significantly better control of spin conversion. One outcome of this high efficiency is that the coherent control in new arrangement can be used for detection of very weak ($`10100`$ Hz) interactions in molecules.

## II Equation of change

First, we give qualitative picture of the process. Let us assume that a test molecule has two nuclear spin states, ortho and para, and that there is a laser radiation resonant to the rovibrational transition, $`mn`$, in the ortho subspace, Fig. 1. The low state $`n`$ is not mixed with para states, but the upper state $`m`$ is mixed by the intramolecular perturbation $`\widehat{V}`$ with the para state $`k`$. In addition, there is the ortho-para level pair, $`m^{}k^{}`$, in the ground vibrational state mixed by another intramolecular perturbation $`\widehat{V}^{}`$. We assume that collisions of the test molecules cannot alter their spin state. This arrangement corresponds to the general formulation of quantum relaxation in which one has collisionally isolated subspaces of states mixed together by internal perturbation . Suppose that the molecule is placed initially in the ortho subspace of the ground vibrational state. Due to collisions the test molecule will undergo fast rotational relaxation inside the ortho subspace. This will proceed until molecule jumps to the state $`m^{}`$, which is directly mixed with the para state $`k^{}`$, or to the state $`n`$ which is mixed by combined action of the external field and intramolecular perturbation $`\widehat{V}`$ with the para state $`k`$. Admixture of a para state implies that the next collision can move the molecule to another para states and thus localize it inside the para subspace. It is clear that resonant radiation can significantly modify the spin conversion. One knows from the literature on resonant interaction of strong laser radiation with matter (see, e.g., ) that radiation can change population of states, split states and create coherences in the system. All these can affect the mixing of ortho and para states and thus the spin conversion process.

Quantitative description of the problem can be performed with the help of kinetic equation for density matrix, $`\widehat{\rho }`$. The molecular Hamiltonian consists of four terms,

$$\widehat{H}=\widehat{H}_0+\mathrm{}\widehat{G}+\mathrm{}\widehat{V}+\mathrm{}\widehat{V}^{}.$$
(1)
The main part, $`\widehat{H}_0`$, has the eigen ortho and para states shown in Fig. 1. $`\mathrm{}\widehat{G}`$ describes the molecular interactions with the external radiation,

$$\widehat{G}=(𝐄_0\widehat{𝐝}/\mathrm{})\mathrm{cos}\omega _Lt,$$
(2)
where $`𝐄_0`$ and $`\omega _L`$ are the amplitude and frequency of the electromagnetic wave; $`\widehat{𝐝}`$ is the operator of the molecular electric dipole moment. We have neglected molecular motion in the operator $`\widehat{G}`$ in order to simplify the theory.

In the representation of the eigen states of $`\widehat{H}_0`$ kinetic equation reads,

$$d\widehat{\rho }/dt=(d\widehat{\rho }/dt)_{coll}i[\widehat{G}+\widehat{V}+\widehat{V}^{},\widehat{\rho }],$$
(3)
where $`(d\widehat{\rho }/dt)_{coll}`$ is the collision integral. Further, collisions in our system will be described by the model standard in the theory of molecular interaction with laser radiation. The off-diagonal elements of $`(d\widehat{\rho }/dt)_{coll}`$ will be assumed to have only decay terms,

$$(d\rho _{j,j^{}}/dt)_{coll}=\mathrm{\Gamma }\rho _{jj^{}};jj^{}.$$
(4)
Here, $`j`$ and $`j^{}`$ indicate rovibrational states of the molecule. These states are assumed to have no degeneracy. The decoherence rates were taken equal for all off-diagonal elements of collision integral. Diagonal terms of the collision integral will be described in the framework of the strong collision model.

Our goal is to determine time dependence of the total concentration of molecules in one spin state. For example, for the total concentration of ortho molecules, $`\rho _o`$, one can get from Eq. (3) the following equation of change,

$$d\rho _o/dt=2Rei(\rho _{mk}V_{km}+\rho _{m^{}k^{}}V_{k^{}m^{}}^{}).$$
(5)
In fact, this result is valid for any model of collision integral, as long as collisions do not alter the molecular spin state which implies, $`_j(d\rho _{jj}/dt)_{coll}=0`$, if $`j`$ ortho, or $`j`$ para.

One has to make a few simplifications in order to find the off-diagonal density matrix elements, $`\rho _{mk}`$ and $`\rho _{m^{}k^{}}`$. We assume $`\widehat{V}`$ and $`\widehat{V}^{}`$ being small and consider zero and first order terms of the density matrix,

$$\widehat{\rho }=\widehat{\rho }^{(0)}+\widehat{\rho }^{(1)}.$$
(6)
We start with zero order perturbation theory. $`\rho ^{(0)}`$ is determined by the kinetic equation,

$$d\widehat{\rho }^{(0)}/dt=(d\widehat{\rho }^{(0)}/dt)_{coll}i[\widehat{G},\widehat{\rho }^{(0)}].$$
(7)
In zero order perturbation theory, para molecules are at equilibrium,

$$\rho _p^{(0)}(g,j)=(n\rho _o^{(0)})w(j);\rho _p^{(0)}(e,j)=0,$$
(8)
where $`n`$ is the total concentration of the test molecules; $`w(j)`$ is the Boltzmann distribution over rotational states. We will assume the same function $`w(j)`$ for each of four vibrational states. We have neglected in Eq. (8) vibrational resonance exchange between ortho and para molecules which would populate the upper vibrational state of para molecules. It was taken into account also that ortho-para exchange is on many orders of magnitude slower than vibrational and rotational relaxations.

Equations for stationary level populations of ortho molecules are obtained from Eqs. (4) and (7). Thus one has,

$`(\nu _V+\nu _R)\rho _o^{(0)}(e,j)`$ $`=`$ $`\nu _Rw(j)\rho _o^{(0)}(e)+\rho _o^{(0)}p\delta _{jm};`$ (9)
$`\nu _R\rho _o^{(0)}(g,j)`$ $`=`$ $`\nu _Rw(j)\rho _o^{(0)}(g)+\nu _Vw(j)\rho _o^{(0)}(e)\rho _o^{(0)}p\delta _{jn};`$ (10)
$`\rho _o^{(0)}p`$ $`=`$ $`{\displaystyle \frac{2\mathrm{\Gamma }G^2}{\mathrm{\Gamma }^2+\mathrm{\Omega }^2}}\left[\rho _o^{(0)}(g,n)\rho _o^{(0)}(e,m)\right],`$ (11)
where $`\rho _o^{(0)}(e)`$ and $`\rho _o^{(0)}(g)`$ are the concentrations of ortho molecules in excited and ground vibrational states; $`\nu _V`$ and $`\nu _R`$ are the rates of vibrational and rotational relaxations; $`p`$ is the probability of optical excitation of ortho molecules. Introduction of different relaxation rates for different degrees of freedom makes the model of strong collisions more accurate. Eqs. (11) correspond to rotational wave approximation. Matrix element of $`\widehat{G}`$ is given by

$$G_{mn}=Ge^{i\mathrm{\Omega }t};GE_0\overline{d}_{mn}/2\mathrm{},$$
(12)
where $`\mathrm{\Omega }=\omega _L\omega _{mn}`$ is the radiation frequency detuning from the absorption line center, $`\omega _{mn}`$; the line over symbol indicate a time-independent factor. Rabi frequency, $`G`$, is assumed to be real.

Solution of Eqs. (11) has no difficulty. Concentration of ortho molecules in excited vibrational state, $`\rho ^{(0)}(e)`$, and in the state $`m`$ which one needs for further calculations read,

$`\rho _o^{(0)}(e)`$ $`=`$ $`\rho _o^{(0)}{\displaystyle \frac{p}{\nu _V}};\rho _o^{(0)}(e,m)=\rho _o^{(0)}{\displaystyle \frac{p}{\nu _m}};p={\displaystyle \frac{2\mathrm{\Gamma }G^2w(n)}{\mathrm{\Gamma }_B^2+\mathrm{\Omega }^2}};`$ (13)
$`\nu _m^1`$ $`=`$ $`w(m)\nu _V^1+(1w(m))(\nu _V+\nu _R)^1.`$ (14)
Here, $`\mathrm{\Gamma }_B`$ is the homogeneous linewidth of the absorption spectrum profile, $`\mathrm{\Gamma }_B^2=\mathrm{\Gamma }^2+2\mathrm{\Gamma }\tau G^2;\tau =\nu _m^1+\nu _n^1;\nu _n^1=w(n)\nu _V^1+(1w(n))\nu _R^1`$. $`\nu _m`$ and $`\nu _n`$ are the effective population decay rates of the corresponding states. In a similar way, one can calculate from Eq. (7) the off-diagonal density matrix element, which amplitude is equal to,

$$\overline{\rho }_o^{(0)}(m|n)=iG\rho _o^{(0)}w(n)\frac{\mathrm{\Gamma }+i\mathrm{\Omega }}{\mathrm{\Gamma }_B^2+\mathrm{\Omega }^2}.$$
(15)
In zero order perturbation theory, one neglects perturbations $`\widehat{V}`$ and $`\widehat{V}^{}`$. It implies that there are no coherences between ortho and para states, $`\rho _{mk}^{(0)}=0`$; $`\rho _{m^{}k^{}}^{(0)}=0`$. Consequently, one has,

$$d\rho _o/dt=2Rei(\rho _{mk}^{(1)}V_{km}+\rho _{m^{}k^{}}^{(1)}V_{k^{}m^{}}^{}),$$
(16)
instead of Eq. (5). Note, that the spin conversion appears in the second order approximation. The first order correction to the density matrix, $`\rho ^{(1)}`$, is determined by the equation,

$$d\widehat{\rho }^{(1)}/dt=(d\widehat{\rho }^{(1)}/dt)_{coll}i[\widehat{G},\widehat{\rho }^{(1)}]i[\widehat{V}+\widehat{V}^{},\widehat{\rho }^{(0)}].$$
(17)
For $`\rho _{m^{}k^{}}^{(1)}`$ one has from this equation,

$$\rho _{m^{}k^{}}^{(1)}=\frac{iV_{m^{}k^{}}^{}}{\mathrm{\Gamma }+i\omega ^{}}[\rho _p^{(0)}(g,k^{})\rho _o^{(0)}(g,m^{})],$$
(18)
where $`\omega ^{}\omega _{m^{}k^{}}`$. $`\rho _{mk}^{(1)}`$ can be obtained from equations which are deduced from Eq. (17),

$`d\rho _{mk}^{(1)}/dt+\mathrm{\Gamma }\rho _{mk}^{(1)}+iG_{mn}\rho _{nk}^{(1)}`$ $`=`$ $`iV_{mk}\rho _o^{(0)}(m);`$ (19)
$`d\rho _{nk}^{(1)}/dt+\mathrm{\Gamma }\rho _{nk}^{(1)}+iG_{nm}\rho _{mk}^{(1)}`$ $`=`$ $`iV_{mk}\rho _o^{(0)}(n|m).`$ (20)
Substitutions, $`V_{mk}=\overline{V}e^{i\omega t},(\omega \omega _{mk});\rho _{mk}^{(1)}=\overline{\rho }_{mk}^{(1)}e^{i\omega t};\rho _{nk}^{(1)}=\overline{\rho }_{nk}^{(1)}e^{i(\omega _L\omega _{kn}t)}`$, transform Eqs. (20) to algebraic equations from which one finds $`\rho _{mk}^{(1)}`$. Using $`\rho _{mk}^{(1)}`$ and $`\rho _{m^{}k^{}}^{(1)}`$ from Eq. (18) one has an equation of change (16) in the form,

$`{\displaystyle \frac{d\rho _o}{dt}}`$ $`=`$ $`{\displaystyle \frac{2\mathrm{\Gamma }|V^{}|^2}{\mathrm{\Gamma }^2+\omega ^2}}\left[\rho _p^{(0)}(g,k^{})\rho _o^{(0)}(g,m^{})\right]`$ (22)
$`2|V|^2Re{\displaystyle \frac{[\mathrm{\Gamma }+i(\mathrm{\Omega }+\omega )]\rho _o^{(0)}(m)+iG\overline{\rho }_o^{(0)}(n|m)}{(\mathrm{\Gamma }+i\omega )[\mathrm{\Gamma }+i(\mathrm{\Omega }+\omega )]+G^2}}.`$
## III Enrichment and conversion

The denominator of the second term in the right-hand side of Eq. (22) is convenient to present as, $`(\mathrm{\Gamma }+i\omega _1)(\mathrm{\Gamma }+i\omega _2)`$, where

$$\omega _{1,2}=\omega +\frac{\mathrm{\Omega }}{2}\pm \sqrt{\left(\frac{\mathrm{\Omega }}{2}\right)^2+G^2}.$$
(23)
New parameters, $`\omega _1`$ and $`\omega _2`$, can be interpreted as the gaps between the two components, $`|m_1>`$ and $`|m_2>`$, of the ortho state $`|m>`$, split by the optical field, and the para state $`|k>`$. The ortho state, $`|m_2>`$, crosses the para state, $`|k>`$, at $`\mathrm{\Omega }=\omega (1G^2/\omega ^2)`$, see Fig. 2. The splitting of states by resonant laser radiation is well-known phenomenon in nonlinear spectroscopy .

Using Eq. (22) one can present equation of change in the final form,

$$d\rho _o/dt=n\gamma _{op}^{}\rho _o\gamma ;\gamma \gamma _{op}^{}+\gamma _{po}^{}\gamma _n^{}+\gamma _n+\gamma _{coh}.$$
(24)
In writing this equation we have neglected in the right-hand side of Eq. (22) small difference between $`\rho _o^{(0)}`$ and the total concentration of ortho molecules, $`\rho _o`$. In Eq. (24) the following partial conversion rates have been introduced. The field independent rates,

$$\gamma _{op}^{}=\frac{2\mathrm{\Gamma }|V^{}|^2}{\mathrm{\Gamma }^2+\omega ^2}w(k^{});\gamma _{po}^{}=\frac{2\mathrm{\Gamma }|V^{}|^2}{\mathrm{\Gamma }^2+\omega ^2}w(m^{}).$$
(25)
The rate $`\gamma _{free}\gamma _{op}^{}+\gamma _{po}^{}`$ determines the equilibration rate in the system without an external field. The field dependent term,

$$\gamma _n^{}=\gamma _{po}^{}p/\nu _V,$$
(26)
appears because of depletion of the ground vibrational state of ortho molecules by optical excitation, Index $`n`$ in $`\gamma _n^{}`$ comes from “noncoherent”, i.e., induced by level populations. Another term of similar “noncoherent” origin appears due to the level population, $`\rho _o^{(0)}(m)`$, in Eq. (22),

$$\gamma _n=2|V|^2\frac{p}{\nu _m}Re\frac{\mathrm{\Gamma }+i(\mathrm{\Omega }+\omega )}{(\mathrm{\Gamma }+i\omega _1)(\mathrm{\Gamma }+i\omega _2)}.$$
(27)
And finally the “coherent” term, $`\gamma _{coh}`$, originated from $`\overline{\rho }_o^{(0)}(n|m)`$, in Eq. (22),

$$\gamma _{coh}=2|V|^2\frac{p}{2\mathrm{\Gamma }}Re\frac{\mathrm{\Gamma }i\mathrm{\Omega }}{(\mathrm{\Gamma }+i\omega _1)(\mathrm{\Gamma }+i\omega _2)}.$$
(28)
Solution to Eq. (24) can be presented as,

$$\rho _o=\overline{\rho }_o+(\rho _o(0)\overline{\rho }_o)\mathrm{exp}(\gamma t);\overline{\rho }_o=n\gamma _{op}^{}/\gamma .$$
(29)
Here $`\gamma `$ is the equilibration rate in the system in the presence of external field; $`\overline{\rho }_o`$ is the stationary concentration of ortho molecules. Without an external radiation (at the instant $`t=0`$), the equilibrium concentration of para molecules is equal to,

$$\rho _p(0)=n\rho _o(0)=n\gamma _{po}^{}/\gamma _{free},$$
(30)
if the Boltzmann factors are assumed to be equal, $`w(k^{})=w(m^{})`$. This implies $`\gamma _{op}^{}=\gamma _{po}^{}`$ (see Eq. (25)), the laser field produces a stationary enrichment of para molecules,

$$\beta \frac{\overline{\rho }_p}{\rho _p(0)}1=12\frac{\gamma _{op}^{}}{\gamma }.$$
(31)
One can see from this equation that external field changes concentration of para isomers if $`\gamma \gamma _{free}`$.

We assume in further analysis the following parameters, $`\omega =100`$ MHz, $`\omega ^{}=130`$ MHz, $`V_{m^{}k^{}}^{}=5`$ kHz, $`\mathrm{\Gamma }=210^8`$ s-1/Torr and the Boltzmann factors of the states $`m^{}`$, $`k^{}`$, $`m`$, and $`k`$ all equal $`10^2`$. This set of parameters gives the field free conversion rate, $`\gamma _{free}=10^2`$ s-1/Torr, which coincides with the conversion rate in 13CH3F. Nuclear spin conversion in these molecules is governed by quantum relaxation (see the review ). The rotational and vibrational relaxation rates will be taken equal, $`\nu _R=0.1\mathrm{\Gamma }`$ and $`\nu _V=0.01\mathrm{\Gamma }`$, respectively.

First, we consider relatively low optical fields, thus small $`G`$. In this case one has two peaks in enrichment at frequencies $`\mathrm{\Omega }\omega `$ and $`\mathrm{\Omega }=0`$, see Fig. 3. The data shown in this figure correspond to $`V_{mk}=3`$ kHz, and $`\mathrm{\Gamma }=2`$ MHz. The peak at $`\mathrm{\Omega }=0`$ appears because the excitation probability, $`p`$, has maximum at this frequency. Amplitude of this peak is determined mainly by the rate $`\gamma _n`$. As $`G`$ grows, the amplitude of the peak 2 reaches the value $`\gamma _n/\gamma _{free}(V\omega ^{}/2V^{}\omega )^2`$ which constitutes $`15\%`$. The peak at $`\mathrm{\Omega }=0`$ in isomer enrichment was predicted in  by considering only the level population effects.

Peak at $`\mathrm{\Omega }\omega `$ appears because the ortho state $`|m_2>`$ crosses the para state $`|k>`$ at this frequency of the external field (see Fig. 2). This peak is determined mainly by $`\gamma _{coh}`$. When $`G`$ increases its amplitude grows to much bigger values than the amplitude of the peak at $`\mathrm{\Omega }=0`$. At resonant frequency $`\mathrm{\Omega }\omega `$ the rate $`\gamma _{coh}`$ is enhanced by large factor $`(\omega /\mathrm{\Gamma })^2`$. This explains much larger enrichment at $`\mathrm{\Omega }\omega `$. Note, that large enrichment occurs only if the excitation probability at this frequency, $`p(\omega )`$, is not very low.

The data shown in Fig. 4 correspond to strong optical field, $`G=50`$ MHz, and three values of $`V_{mk}`$. $`\mathrm{\Gamma }`$ was taken equal 2 MHz. One can see, that strong optical field is able to convert almost all molecules to the para state if $`V_{mk}V_{m^{}k^{}}`$. Thus relatively weak (3 kHz) coupling in upper state is able to produce macroscopic effect, viz., almost complete enrichment of spin isomers. It is of fundamental importance, that even for much weaker coupling in upper state, enrichment is still significant. For example, if the perturbation in upper state, $`V_{mk}=30`$ Hz, one has the enrichment, $`\beta 1`$%. Enrichment at this level can easily be measured. It is important that the enrichment peak at $`\mathrm{\Omega }\omega `$ is narrow (the width $`\mathrm{\Gamma }`$) and thus can be distinguished from much wider structures (the width $`\mathrm{\Gamma }_B`$) induced by population effects.

Equilibration rate in the system is given by $`\gamma `$, see Eq. (29). It is convenient to characterize the conversion rate in relative units,

$$\gamma _{rel}=\gamma /\gamma _{free}1,$$
(32)
Conversion rate, like enrichment, has two peaks in its frequency dependence at low $`G`$. If Rabi frequency, $`G`$, is large and the ortho-para couplings in upper and low states have the same order of magnitude, conversion can be significantly enhanced (Fig. 4, upper panel). Again, this enhancement appears because of the crossing of ortho and para states in upper vibrational state by external field.

## IV Discussion

The phenomenon considered in the paper is based on the level splitting produced by resonant electromagnetic radiation. Sometimes, this splitting is called in optics the dynamic Stark effect. The essence of the effect can be understood as follows. Mixing of ortho and para states depends on magnitude of the perturbation $`\widehat{V}`$ but also on the ortho-para level gap. Optical field splits the molecular state and thus change the gaps between the ortho and para states. Conversion rate is significantly enhanced when the ortho and para states cross. Similar enhancement occurs when ortho and para states are crossed by ordinary Stark effect in an external DC electric field .

In the same way, one can understand high sensitivity of the phenomenon to weak ortho-para couplings in excited state. Equilibrium concentrations of ortho and para molecules are achieved when one has the ortho-to-para flux in excited state equal to the back flux in the ground vibrational state. Back flux is slow because it is determined by the non-degenerate ortho-para level pair. On the other hand, the flux in excited state can be significantly enhanced by proper choice of radiation parameters which allows to cross the ortho and para states in upper vibrational state.

High efficiency of the proposed enrichment method can be used to detect weak perturbations in excited vibrational state. First of all, it can be the hyperfine perturbations of the same origin and similar magnitude as the perturbations in the ground vibrational state. We have seen that hyperfine coupling of the order of $`10^3`$ Hz is able to convert almost all molecules in one spin state. One can also detect much weaker interactions in molecules. In this case one should select the ortho-para level pair in excited vibrational state which is not mixed by “ordinary” hyperfine interactions in order to avoid the weak interaction to be hidden by stronger, ordinary hyperfine interactions. An interesting case is the crossings of states having opposite parity. Mixing of such ortho and para states can be performed only by spin-dependent, parity-odd interactions which are not observed in molecules yet.

## V Conclusions

We have performed analysis of the spin isomer enrichment and conversion governed by molecular rovibrational excitation. This analysis was done using a few simplifications. We have neglected the Doppler broadening of the absorbing transition, degeneracy of molecular states, and resonant vibrational exchange between excited and unexcited molecules. These simplifications are not crucial for the existence of the phenomenon. More detailed analysis will be done elsewhere.

We have shown that coherent control of nuclear spin conversion in molecules can be efficiently performed by strong radiation resonant to rovibrational molecular transition. A possible applications of this phenomenon is the enrichment of molecular spin isomers. Another application is the detection of very weak ($`10100`$ Hz) interactions in molecules, which can be, e.g., parity-odd interactions.

## Acknowledgments

This work was made possible by financial support from the Russian Foundation for Basic Research (RFBR), grant No. 98-03-33124a.


text:
# Flow at the SPS and RHIC as a Quark-Gluon Plasma Signature

\[

## Abstract

Radial and elliptic flow in non-central heavy ion collisions can constrain the effective Equation of State(EoS) of the excited nuclear matter. To this end, a model combining relativistic hydrodynamics and a hadronic transport code(RQMD) is developed. For an EoS with a first order phase transition, the model reproduces both the radial and elliptic flow data at the SPS. With the EoS fixed from SPS data, we quantify predictions at RHIC where the Quark Gluon Plasma(QGP) pressure is expected to drive additional radial and elliptic flow. Currently, the strong elliptic flow observed in the first RHIC measurements does not conclusively signal this nascent QGP pressure. Additional measurements are suggested to pin down the EoS.

\]

1. By colliding heavy nuclei at the SPS and RHIC accelerating facilities, physicists hope to excite hadronic matter into a new phase consisting of deconfined quarks and gluons – the Quark Gluon Plasma(QGP). After the collision, the produced particles move collectively or $`flow`$ and this flow may quantify the effective Equation of State(EoS) of the matter. In central PbPb collisions at the SPS, a strong radial flow is observed. The matter develops a collective transverse velocity approaching (1/2)c. In non-central collisions, a radial and an $`elliptic`$ flow are observed. Since in non-central collisions the initial nucleus-nucleus overlap region has an elliptic shape, the initial pressure gradient is larger along the impact parameter and the matter moves preferentially in this direction.

The phase transition to the QGP influences both the radial and elliptic flows. QCD lattice simulations show an approximately 1st order phase transition. Over a wide range of energy densities $`e=.51.4GeV/fm^3`$, the temperature and pressure are nearly constant. Over this range then, the ratio of pressure to energy density $`p/e`$, decreases and reaches a minimum at a particular energy density known as the softest point, $`e_{sp}1.4GeV/fm^3`$. When the initial energy density is close to $`e_{sp}`$, the small pressure (relative to $`e`$) cannot effectively accelerate the matter. However, when the initial energy density is well above $`e_{sp}`$, $`p/e`$ approaches 1/3, and the larger pressure drives collective motion . At a time of $`1fm/c`$, the energy densities at the SPS($`\sqrt{s}_{NN}=17GeV`$) and RHIC ($`\sqrt{s}_{NN}=130GeV`$) are very approximately 4 and $`7GeV/fm^3`$ respectively. Based on these experimental estimates, the hard QGP phase is expected to live significantly longer at RHIC than at the SPS. The final flows of the produced particles should reflect this difference. In this paper we pose the question: Can both the radial and elliptic flow at the SPS and RHIC be described by a single effective EoS?

Since the various hadron species have different elastic cross sections, they freezeout (or decouple) from the hot fireball at different times. Because flow builds up over time, it is essential to model this differential freezeout. It was ignored in previous hydrodynamic simulations of non-central heavy ion collisions and elliptic flow was over-predicted flow by a factor of two.

2. The Hydro to Hadrons(H2H) model will be described in detail elsewhere. Other authors have previously constructed a similar model for central collisions. The model evolves the QGP and mixed phases as a relativistic fluid, but switches to a hadronic cascade (RQMDv2.4) at the beginning of the hadronic phase to model differential freezeout. The computer code consists of three distinct components. Assuming Bjorken scaling, the first component solves the equations of relativistic hydrodynamics in the transverse plane and constructs a switching surface at a temperature, $`T_{switch}=160MeV`$. The second component generates hadron on the switching surface using the Cooper-Frye formula with a theta function rejecting backward going particles . Finally, the third component (RQMD) sequentially re-scatters the generated hadrons until freezeout.

For the hydrodynamic evolution, a family of EoSs was constructed with an adjustable Latent Heat(LH)

(see Fig. 1). LH$`\mathrm{}`$ is considered as a limiting case, mimicking non-equilibrium phenomena. The hadron phase exists up to a critical temperature of $`T_c=165MeV`$, and consists of ideal gas mixture of the meson pseudoscalar and vector nonets and the baryon octet and decuplet. The hadron phase is followed by a mixed phase with a specified LH, which is finally followed by a QGP phase with $`C_s^2=1/3`$. In addition, a Resonance Gas(RG) EoS was constructed with a constant speed of sound above the hadron phase.

3. Radial flow is quantified experimentally by slope parameters, $`T_{slope}`$; the momentum spectrum of each particle is fit to the form $`dN/dM_T^2dy|_{y=0}=Ce^{M_T/T_{slope}}`$ where $`M_T^2=P_T^2+m^2`$. $`T_{slope}`$ incorporates random thermal motion and the collective transverse velocity.

In Fig. 2, the pion and proton slope parameters are plotted as functions of the total charged particle multiplicity in the collision. Look first at the leftmost points at SPS multiplicities and compare the model and experimental slopes: The proton slope data favor a relatively hard EoS – LH8 or harder. A direct comparison of the model to published spectra  supports this claim. A RG EoS can also reproduce the proton flow. A similar analysis of elliptic flow (shown and quantified below) favors a relatively soft EoS – LH8 or softer. With some caveats, LH8 represents a happy middle which can reproduce both the radial and elliptic flow at the SPS.

Look now at the energy/multiplicity dependence of the slopes. For all EoSs, $`T_{slope}`$ increases with the collision energy. For a soft EoS (e.g. LH$`\mathrm{}`$) the increase is small and for a hard EoS (e.g. LH8) the increase is large. At RHIC multiplicities, the difference between the slope parameters is large and easily experimentally observable.

4. Elliptic flow is quantified experimentally by the elliptic flow parameter, $`v_2=\mathrm{cos}(2\mathrm{\Phi })`$; here $`\mathrm{\Phi }`$ is the angle around the beam measured relative to the impact parameter and $``$ denotes an average over the single particle distribution, $`\frac{dN}{dP_Td\mathrm{\Phi }}`$. $`v_2(P_T)`$ is found by holding $`P_T`$ constant in while averaging $`\mathrm{cos}(2\mathrm{\Phi })`$ over $`\frac{dN}{dP_Td\mathrm{\Phi }}`$. $`v_2`$ measures the response of the fireball to the the spatial deformation of the overlap region, which is usually quantified in a Glauber model by the eccentricity $`ϵ=y^2x^2/x^2+y^2`$. Since the response($`v_2`$) is proportional to the driving force($`ϵ`$), the ratio $`v_2/ϵ`$ is used to compare different impact parameters and nuclei.

In Fig. 3(a), the number elliptic flow ($`v_2`$) is plotted as a function of charged particle multiplicity at an impact parameter of 6 fm. Before studying the energy dependence, look at the magnitude of the elliptic flow at the SPS. For LH8, the stars show the pion $`v_2`$ when the matter is evolved as a fluid until a decoupling temperature of $`T_f=120MeV`$; they illustrate the excessive elliptic flow typical of pure hydrodynamics. Once a cascade is included, LH8 (the squares) is only $`20\%`$ above the data – a substantial improvement. Typically in hydrodynamic calculations, the freezeout temperature $`T_f`$ is adjusted to fit the proton $`P_T`$ spectrum. However, protons are driven by a pion “wind” and decouple from the fireball $`5fm/c`$ after the pions on average. This pion wind accounts for the strong proton flow at the SPS and is not described by ideal hydrodynamics. In order to match the observed proton flow, hydrodynamic calculations must decouple at low freezeout temperatures, $`T_{frz}120MeV/c`$. This low temperature has two consequences for elliptic flow: First the reduction of elliptic flow due to resonance decays is small $`15\%`$, compared to $`30\%`$ in the H2H model. Second, compared to a cascade, the hydrodynamics generates twice as much elliptic flow during the late cool hadronic stages of the evolution. By including the pion “wind”, and more generally by decoupling differentially, we can simultaneously describe the radial and elliptic flow data at the SPS.

The energy dependence of $`v_2`$ is the central issue. As seen in Fig. 3, the H2H model predicts an increase in elliptic flow by a factor $`1.4`$ and is in reasonable agreement with SPS and RHIC flow data. This result was presented prior to the publication of RHIC data . In contrast, UrQMD, a hadronic cascade based on string dynamics, predicts a decrease by a factor of $`2`$ . This is because the UrQMD string model has a super-soft EoS at high energies. For pure hydrodynamics as illustrated by the stars, $`v_2`$ is approximately constant (but see ). For HIJING, a model which considers only the initial parton collisions, $`v_2`$ is $`0`$. The first RHIC data clearly contradict these models.

The increase in $`v_2`$ is now used to constrain the EoS of the excited matter. The QCD phase diagram has two

distinguishing features. It is soft at low energy densities and subsequently hard at high energies. A RG EoS(the open squares) has no softness and the elliptic flow is clearly too strong both at the SPS and RHIC. The entire family of EoSs, LH8 through LH$`\mathrm{}`$, reproduces the elliptic flow data in both energy regimes. Counter-intuitively, as the latent heat is increased, $`v_2`$ first decreases and then increases. In the final count, LH8 and LH$`\mathrm{}`$ have roughly the same $`v_2`$. However, they develop the $`v_2`$ in different ways. For LH8, the EoS shifts from hard to soft and the early pressure starts an early elliptic expansion. For LH$`\mathrm{}`$, the EoS is just soft and the elliptic expansion stalls. However because the expansion is stalled, the LH$`\mathrm{}`$ collision lifetime($`13fm/c`$ at RHIC) is significantly longer than the LH8 lifetime($`9fm/c`$ at RHIC). Over the long LH$`\mathrm{}`$ lifetime, $`v_2^{LH\mathrm{}}`$ steadily grows and is finally comparable to $`v_2^{LH8}`$. As the latent heat is increased from LH8 to LH16, the EoS becomes becomes softer and $`v_2`$ at first decreases. However, as the EoS is made softer still, the lifetime increases and $`v_2`$ rises again.

5.Impact Parameter Dependence. In Fig. 4, $`v_2`$ for LH8 as a function of the number participants($`N_p`$) is compared to data. Different EoSs show a similar participant (or b) dependence. The agreement is good at RHIC where the multiplicity is high. For ideal hydrodynamics, $`v_2ϵ(N_p^{max}N_p)`$. In the low density limit, since the

response is proportional to the number of collisions, $`v_2ϵ\frac{dN}{dy}(N_p^{max}N_p)N_p`$. Therefore, $`v_2`$ has a different $`N_p`$(or b) dependence in the hydrodynamic and low density limits. At RHIC, except in very peripheral collisions, the $`N_p`$ dependence is clearly linear and strongly supports the hydrodynamic limit. At the SPS, the $`N_p`$ dependence may not be clearly linear, but it also does not follow the low density limit. Two-pion correlations, may change the data analysis, reduce $`v_2`$ in the periphery and improve the low density agreement.

Finally in Fig. 5, $`v_2`$ is studied both as a function of transverse momentum and impact parameter. For both LH8 and LH$`\mathrm{}`$, the calculation produces too much elliptic flow in peripheral collisions(45-85%), and too $`little`$ elliptic flow in the most central collisions(0-11%). The $`P_T`$ dependence of $`v_2`$ also clarifies the difference between and LH8 and LH$`\mathrm{}`$: LH$`\mathrm{}`$, a super soft EoS, generates elliptic flow only at low momentum while LH8, a hard EoS, generates elliptic at high momentum.

6.Summary and Discussion. By incorporating differential freezeout, the Hydro to Hadrons(H2H) model simultaneously reproduces the radial and elliptic flow at the SPS and RHIC. At the SPS, the radial flow demands an EoS with a latent heat $`LH0.8GeV/fm^3`$, while elliptic flow demands an EoS with a latent heat $`LH0.8GeV/fm^3`$. Further, in contrast to string and collision-less parton models, the increase in $`v_2`$ is naturally explained using hydrodynamics. This challenges the prevailing view that the SPS is in the low density regime and that the increase in $`v_2`$ represents a transition

to the hydrodynamic regime. However, the increase in $`v_2`$ does not uniquely signal the asymptotic QGP pressure. Indeed, at RHIC collision energies, a very soft EoS can have the same $`v_2`$ as an EoS with a well developed QGP phase. This EoS is not academic since softness can mimic non-equilibrium phenomena. To reveal the underlying EoS and the burgeoning QGP pressure, the collision energy should be scanned from the SPS to RHIC. If the prevailing low density view of the SPS is correct, a transition in the b dependence of elliptic flow should be observed over the energy range. In addition, for different EoSs, $`v_2`$ depends differently on collision energy and transverse momentum (Fig. 3 and Fig. 5). Taken with the radial flow (Fig. 2), this experimental information would help settle the EoS of hot hadronic matter.

Acknowledgments. The work is partly supported by the US DOE grant No. DE-FG02-88ER40388 and grant No. DE-FGO2-87ER 40331.

Note Added. After the submission of this work, the STAR and PHENIX collaborations reported proton at anti-proton spectra. The preliminary spectra favor LH8-LH16 and disfavor LH$`\mathrm{}`$.