text:
# 1 Introduction

## 1 Introduction

The ongoing and planned experiments involving heavy ion collisions at RHIC, CERN’s SPS and at LHC are aimed at studying the behavior of quantum chromodynamics (QCD) at high energies. It is believed that at these high energies  (100 GeV/nucleon at RHIC and 3 TeV/nucleon at LHC in the c.m. frame) a hot baryon free plasma of quarks and gluons (QGP) is expected to be created even as at lesser energies, nuclear collisions at Brookhaven’s AGS and at SPS show large stopping of the nuclei and hint at the baryon rich matter. The QGP plasma will soon undergo a phase transition to hadronic matter. The phase transition from QGP to hadron resonance gas (HRG) has been extensively studied in the context of the early Universe where the time scale is of the order of $`10^610^7`$ s \[2-4\]. Perturbative finite temperature QCD methods along with relativistic transport theory have been used to study the early stages of the nuclear collisions by computer simulations . The initial quantitative picture that emerges is that of a QGP fireball of $`150fm^3`$ created at the initial temperature of $`300350`$ MeV, which is about twice the expected critical temperature $`T_c`$. Studies of quark-hadron phase transition in the early Universe, in heavy ion collisions and in high density nuclear matter in the core of compact stars  point to the importance of interactions in the two phases. It has been shown that unless hadronic interactions are taken into account, the hadronic phase again becomes thermodynamically the preferred phase at high temperature and that QCD interactions are important in the QGP phase  making the critical parameters depend on the interactions. The lattice results also point to the departure of the QGP equation of state (EOS) from the ideal gas EOS even at temperatures well above $`T_c`$. However there still remain uncertainties in the equation of state in both the confined and the deconfined phases. It has also been pointed out that in situations in which nucleation takes place, the curvature energy $`8\pi \gamma r`$, a term in the free energy in addition to the surface energy term $`4\pi r^2\sigma `$ (where $`\sigma `$ and $`\gamma `$ are the surface and curvature energy densities respectively) also plays an important role and should be kept in the calculation of the nucleation rate \[9-12\]. It is in this context that the above experiments will go a long way in establishing the dynamics of the phase transition and the final picture that will emerge would be of great interest in astrophysics and cosmology. If the phase transition from quark gluon plasma to confined hadronic matter is a first order thermodynamic phase transition, the transition would proceed through the nucleation of hadronic bubbles. Conventionally, the QGP would supercool to temperature below $`T_c`$ only then the nucleation starts. The bubbles greater than the critical size would expand driving the QGP phase into smaller and smaller regions thereby converting more and more matter into the hadronic phase. The expansion of the bubble would be accompanied by the release of latent heat which will reheat the plasma shutting off further nucleation. If the phase transition has to be completed, the latent heat must be gotten rid of. Such studies have been made in the context of the early Universe where it has been shown that phase transition proceeds in the expanding Universe after the Universe supercools to a temperature below $`T_c`$, the expanding bubbles then reheat the Universe expelling the QGP phase till the entire Universe is converted to the hadronic phase. The time of completion and the degree of supercooling depends on the EOS and the parameters like the Bag pressure and the surface tension .

Csernai and Kapusta  have applied a recently computed nucleation rate to a first order phase transition in a set of rate equations to study the time evolution of QGP in heavy ion collisions. Based on Bjorken hydrodynamics and on current parameter values, they find the transition generates about 30 percent extra entropy and also a time delay of $`11fm/c`$ in completion of the transition. They also estimate that the system should supercool about 20 percent below $`T_c`$ before nucleation of hadronic bubbles is rapid enough to start reheating the system. Csernai and Kapusta for the purposes of illustration used noninteacting system of hadrons and QGP. They modelled the hadronic phase by a massless gas of pions and the plasma phase by a gas of gluons and massless quarks of two flavors. In the light of the importance of taking interactions both in the hadronic as well as in the plasma phase for purposes of investigating phase transitions at high energies, temperature and densities as discussed above, we calculate the rate of formation of the hadron bubbles for interacting QGP and hadron resonance gas. It is our endeavour here to examine closely the effects of interactions on the dynamics of phase transition as the hot plasma of quarks and gluons is created and the temperature drops to $`T_c`$ where a phase mixture of QGP and HRG develops. For this purpose QGP is treated as a gas of massless $`u`$,$`d`$ quarks, massive $`s`$ quarks and massless gluons. QCD interactions are treated perturbatively to order $`g^3`$ and the long range confinement effects are parameterized by the bag pressure $`B`$. For the hadron gas we use the known spectrum of low lying baryons and meson resonances and the repulsive interactions between them by incorporating the hard core excluded volume effects. In the literature hadronic interactions have been taken in many different ways (see Ref.7 for example), hadronic interactions can be accounted for through an average mean field repulsive potential arising from the Reid potential and $`\pi \pi `$ effective interactions or alternatively, repulsive effects can also be included by modifying the chemical potential as done by Kapusta et al. in a simple method of implementing the mean field theory. One could also use the full fledged framework of relativistic mean field theory itself. It has however been shown by Goyal et al. that in the baryon free case , all the interaction schemes mentioned above give roughly the same value of the hadronic pressure, there are however, marked deviations at very high $`T`$ and/or at high $`\mu _B`$ regime.

The effect of curvature energy term on the nucleation of hadronic bubbles has been studied earlier . It has been found that it leads to a minima in the free energy in addition to a maxima below the critical temperature and also a minima above $`T_c`$, giving rise to spontaneous nucleation of equilibrium sized bubbles. This may contribute to the completion of the phase transition and has interesting consequences. We study here the effect of curvature on the dynamics of the phase transition for several values of the parameters and discuss their consequences.

The plan of the paper is as follows. In section 2 we discuss the bubble nucleation along with the importance of the prefactor which plays a crucial role in guiding the behavior of the nucleation rate. It turns out that the exponential factor is of less relevance here unlike the case of early Universe where exponent dominates the behavior of the nucleation rate. We also study here the effect of curvature on the nucleation rate and give the two rate equations to study the time evolution of the transition for different parameter values. In section 3 we discuss our result and compare them with the case when interactions and the curvature energy are not present. We also look at the supercooling that takes place and evaluate the behavior of the system with different parameter values. Finally in section 4 we give our conclusions.

## 2 Bubble nucleation

When the QGP cools through the critical temperature $`T_c`$, the new hadron phase becomes preferred. Energetically the new phase remains unfavourable as there is surface free energy associated with the QGP and the hadron gas interface. Creation of bubbles of the new phase are thus unfavourable and all nucleated bubbles with radius less than the critical radius die out but those with radii grater than the critical radius expand untill they coalesce with each other. Thus supercooling occurs before the new hadron phase actually appears and takes over. Consequently reheating takes place due to the release of latent heat. The critical radius is obtained by extremizing the thermodynamic work expanded to create a bubble, i.e.

$$W\frac{4\pi }{3}r^3(P_hP_q)+4\pi r^2\sigma 8\pi (\gamma _q\gamma _h)r,$$
(1)
where $`P_h`$ and $`P_q`$ are the pressures in the hadronic and quark phase respectively, $`\sigma `$ is the surface tension and $`\gamma =(\gamma _q\gamma _h)`$ is the curvature coefficient and has been estimated in . Following Mardor and Svetitsky the curvature energy in the MIT model for massless quarks is given by

$$E_c=\frac{gr}{3\pi }_0^{\mathrm{}}𝑑kk\{1+exp(k\mu _B)/T\}^1$$
(2)
where $`g`$ is the statistical weight and $`r`$ the radius. We get $`\gamma \gamma _q\frac{T^2}{16}`$, where we have assumed $`\gamma _h<<\gamma _q`$ following  and $`\mu _B`$ to be negligible compared to $`T`$. According to Madsen  , the assumption of zero strange-quark mass in the curvature term should be relaxed but presently no calculation of the curvature term for massive quarks exists. The surface tension is mainly due to the finite strange quark mass as the massless particles do not contribute to the surface term in the free energy. The critical radius is obtained by putting $`\frac{W}{r}=0`$ and we get an extremum below $`T_c`$ given as

$$r_{c+}=\frac{\sigma }{\mathrm{\Delta }P}(1+\sqrt{1\beta })$$
(3)
where $`\mathrm{\Delta }P=P_hP_q`$ and $`\beta =\frac{2\mathrm{\Delta }P}{\sigma ^2}\gamma `$. Since $`\beta `$ is positive definite, a real solution exists only if $`\beta 1`$. This critical radius is actually a maximum. There is another extremum solution (a minima) given by the critical radius

$$r_c=\frac{\sigma }{\mathrm{\Delta }P}(1\sqrt{1\beta })$$
(4)
These bubbles correspond to equilibrium bubbles which stabilse to their critical size $`r_c`$ and do not grow. The free energies of the critical sized bubbles are given by

$$W_{c+}=\frac{4\pi \sigma ^3}{3\mathrm{\Delta }P^2}[2+2(1\beta )^{3/2}3\beta ]$$
(5)
and

$$W_c=\frac{4\pi \sigma ^3}{3\mathrm{\Delta }P^2}[22(1\beta )^{3/2}3\beta ]$$
(6)
respectively for the two critical radii. $`\beta `$ decides the effect of the curvature term and decreases the critical radius for the expanding bubbles. The existence of a critical radius corresponding to the minimum of the free energy is a consequence of the curvature term. In fact there is a extremum in the free energy even above the critical temperature for all starting temperatures of the plasma. This extremum is also a minimum and the corresponding critical radius and the free energy are given by

$$r_{c+>}=\frac{\sigma }{\mathrm{\Delta }P}(1+\sqrt{1+\beta })$$
(7)
$$W_{c+>}=\frac{4\pi \sigma ^3}{3\mathrm{\Delta }P^2}[22(1+\beta )^{3/2}+3\beta ]$$
(8)
Significantly, the restriction of $`\beta 1`$ is not present for these bubbles. This feature of the curvature term ensures that the phase transition actually begins well above the critical temperature (as soon as the plasma is formed) by the nucleation of equilibrium sized bubbles of the hadron phase, and this continues till below the critical temperature when the expanding bubbles also join them and start the rapid nucleation and expansion of the hadron phase. As pointed out earlier, below $`T_c`$ if $`\beta `$ becomes greater than one, nucleation of both types of critical bubbles dissapear. The phase transition can now take place by the expansion of all the already nucleated critical bubbles $`r_{c+}`$, $`r_c`$ and $`r_{c+>}`$. This is because the critical radius $`r_c`$ and the critical free energy $`W_c`$ necessary for creation of growing bubbles become zero. This essentially means that the phase transition now gets rolling and completes fast. Sooner $`\beta `$ becomes greater than $`1`$, faster will this process dominate the phase tranisition. As we will see , for certain reasonable values of the parameters ($`B`$ and $`\sigma `$) $`\beta `$ becomes greater than $`1`$ almost near $`T_c`$ causing unrestrained expansion of all the bubbles (ones above $`T_c`$ as well as whatever exists below it) and completion of the phase transition.

The bubble nucleation rate(number of bubbles formed per unit time per unit volume) at temperature $`T`$ is given by

$$I=I_0e^{W_c/T}$$
(9)
where $`I_0`$ is the prefactor having dimensions of $`T^4`$. The prefactor used traditionally in the early Universe studies is given by $`I_0=(W_c/2\pi T)^{3/2}T^4`$. Csernai and Kapusta  have recently calculated this in a course grain effective field theoretic approximation to QCD and give

$$I_0=\frac{16}{3\pi }\left(\frac{\sigma }{3T}\right)^{3/2}\frac{\sigma \eta _qr_c}{\xi _q^4(\mathrm{\Delta }w)^2}$$

where $`\eta _q=14.4T^3`$ is the shear viscosity in the plasma phase, $`\xi _q`$ is a correlation length of order 0.7 fm in the plasma phase, and $`\mathrm{\Delta }w`$ is the difference in the enthalpy densities of the two phases. This nucleation rate is limited by the ability of dissipative process to carry latent heat away from the bubbles surface as indicated by the dependence on the viscosity. The above expression however gets modifed in the present case because of the presence of the equilibrium bubbles. We now have different $`I`$’s depending on which kind of bubbles are being discussed, namely $`I_+`$, $`I_{}`$ and $`I_{+>}`$ corresponding to critical radiuses $`r_{c+}`$, $`r_c`$ and $`r_{c+>}`$ respectively.

The pressure in the QGP phase can be calculated by using the thermodynamic potential given in Kapusta and Shuryak  . For the hadron phase we use the known masses of the low lying 33 baryons and 45 mesons whose masses and degeneracy factors are taken from the Particle Data Group Summary . The hadronic pressure and number densities for non-interacting point baryons and mesons are given by the usual thermodynamic relations. One of the simplest ways to account for the short range repulsive forces as discussed above is by considering the finite volume of the baryons that modify the space available for occupation in a manner akin to that of a Van der Waal’s equation. We take baryons and antibaryons to have the same size as protons, given by $`V_p=m_p/4B`$ where $`B`$ is the bag pressure. A radius $`r_p`$ lying between 0.6 fm and 0.8 fm is often used  . The hadronic pressure and baryon densities corrected for finite volume effects are now given by

$$P_h=\frac{_bP_b^{pt}}{1+_bn_b^{pt}V_p}+\frac{_{\overline{b}}P_{\overline{b}}^{pt}}{1+_{\overline{b}}n_{\overline{b}}^{pt}V_p}+\underset{m}{}P_m^{pt}$$
(10)
$$n_B=\frac{_bn_b^{pt}}{1+_bn_b^{pt}V_p}+\frac{_{\overline{b}}n_{\overline{b}}^{pt}}{1+_{\overline{b}}n_{\overline{b}}^{pt}V_p}$$
(11)
where $`b`$, $`\overline{b}`$ and m stand for baryons, anti-baryons and mesons respectively. The pressure equilibrium between the two phases sets $`T_c`$ which is decided by the parameters in the theory, essentially by the parameters $`B`$ and $`\sigma `$. With nucleation rate given by eq.(9), one would like to know the fraction of space $`h(t)`$ which is in the hadronic phase at the given time t as measured in a local comoving frame of the expanding system. Guth and Weinberg  have proposed a formula for $`h(t)`$ in the cosmological first order phase transition. Csernai and Kapusta  have given a more accurate kinetic equation than the one used by Guth and Weinberg and it allows for the transition to complete. If the system cools to $`T_c`$ at time $`t_c`$, then at some later time $`t`$, the fraction of space which has been converted to hadronic gas is given by them as

$$h(t)=_{t_c}^t𝑑t^{}I(T(t^{}))\{1h(t^{})\}V(t^{},t)$$
(12)
where $`V(t,t^{})`$ is the volume of a bubble at time $`t`$ which had been nucleated at the earlier time $`t^{}`$; taking bubble growth into account. Here we have neglected bubble collisions and fusion. The above expression however gets modifed in the present case because of the presence of the equilibrium bubbles. Now we have an additional contribution to the hadron fraction from the bubbles that are nucleated above the critical temperature plus those that are nucleated below the critical temperature and do not grow in size (unless $`\beta >1`$). So, the new expression is now given by

$$h(t)=_{t_c}^t𝑑t^{}I_+(T(t^{}))\{1h(t^{})\}V(t^{},t)$$

$$+_{t_c}^t𝑑t^{}I_{}(T(t^{}))\{1h(t^{})\}\frac{4\pi }{3}r_c^3(T(t^{}))$$

$$+_{t_0}^{t_c}𝑑t^{}I_{+>}(T(t^{}))\{1h(t^{})\}\frac{4\pi }{3}r_{c+>}^3(T(t^{}))$$
(13)
The dynamical equation which couples time evolution of temperature to $`h(t)`$ is obtained by using the longitudinal scaling hydrodynamics of Bjorken. It is given by

$$\frac{de}{dt}=\frac{w}{t}$$
(14)
where $`e(T)`$ and $`w(T)`$ are the energy and enthalpy densities respectively being given by $`e(T)=h(t)e_h(t)+(1h(t))e_q(t)`$, $`e_h`$ and $`e_q`$ being the hadron and quark energy densities at time t respectively, likewise for the enthalpy density. Following  the growth velocity of bubbles is taken to be $`v(T)=v_0(1T/T_c)^{3/2}`$ where $`v_0=3c`$ for $`T>2/3T_c`$. This gives $`v<c/\sqrt{3}`$, the speed of sound for a massless gas. We have assumed that if supercooling takes $`T`$ below $`2/3T_c`$ then $`v`$ is given by $`c/\sqrt{3}`$. This ensures that closer $`T`$ is to $`T_c`$ slower the bubbles grow. The volume of an expanding bubble at time $`t`$ is then given by

$$V(t^{},t)=\frac{4\pi }{3}\{r_{c+}(T(t^{}))+_t^{}^t𝑑t^{\prime \prime }v(T(t^{\prime \prime }))\}^3$$
(15)
In calculations we use the parameters $`\xi _q=0.7fm`$ and $`\eta _q=14.4T^3`$ as in . $`\sigma `$ and $`B`$ are varied to study the behavior of the phase transition as the QGP cools and hadron bubbles are formed. We solve the two coupled equations (13) and (14) for different values of the parameters $`B`$ and $`\sigma `$. The nucleation time is defined as

$$\tau ^1=\frac{4\pi }{3}r_c^3I$$

neglecting bubble growth. In our case we have a contribution to $`\tau `$ from $`I_+`$, $`I_{}`$ below $`T_c`$ and $`I_{+>}`$ above $`T_c`$.

## 3 Results and Discussion

We find that the inclusion of the interactions in both the plasma phase and the hadron gas phase lowers the critical temperature, so higher values of the bag pressure are required for having reasonable values of $`T_c`$. For example, the value $`B^{1/4}=300MeV`$ gives a $`T_c188MeV`$ whereas $`B^{1/4}=235MeV`$ gives a $`T_c144MeV`$ . We show below plots of free energy vs radius of a hadron bubble, nucleation time vs temperature, temperature vs time and hadron fraction vs time etc. In fig.1 we plot the free energy of a hadron bubble against it’s radius in the QGP plasma. For temperatures just above the critical temperature we see that the effect of the curvature term is to produce a minima in the free energy which is to be contrasted with no extremum without a curvature term. Fig.1a pertains to $`T=1.01T_c`$ and the dashed and the dotted curves are without the curvature term. It is clearly seen that equilibrium sized hadron bubbles appear even above the critical temperature (solid and long dashed curves). Fig.1b is plotted for a different set of parameters. The lower panel is plotted just below the critical temperature at $`T=0.99T_c`$. Here we find that below the critical temperature the curvature effect produces a minima as well as a maxima in the free energy. The maxima bubbles grow whereas the minima bubbles reach an equilibrium size. The dashed and the dotted curves are without the curvature effect showing the presence of only the maximum sized bubbles. Figures 1c and 1d plot the same curves but for a different set of the parameters. It is interesting to note that for $`B^{1/4}=300MeV`$, $`\sigma 7MeV/fm^2`$ ($`\sigma ^{1/3}=65MeV`$) and $`T=0.99T_c`$ the free energy has no extremums, only a maximum value at zero (long dashed), indicating that critical radius is zero and all existing bubbles can now expand to reduce the free energy. In such cases $`\beta `$ becomes greater than 1 soon after $`T_c`$ is crossed and there cannot be nucleation of any fresh bubbles. The expansion of all the existing bubbles ensures completion of the transition. In Fig.2 we have shown the plot of $`T/T_c`$ agains time. The temperature of the QGP plasma decreases with time till $`T_c`$, when the phase transition conventionally begins and supercooling starts. The degree of supercooling depends on the parameters and can be a significant fraction of $`T_c`$. The important feature to see here is that the prefactor $`I_0`$ is the deciding factor for supercooling and the duration of the transition rather than the exponent. On the other hand in the early Universe studies, it is the exponent which decides the behavior of the transition. In fig.2a the dashed curve (without curvature) and the solid curve (with curvature) are for parameters $`B^{1/4}=235MeV`$ and $`\sigma =50MeV/fm^2`$ ($`\sigma ^{1/3}125MeV`$). As can be seen they overlap, indicating that in this case the curvature term does not make any difference. This value of the surface tension is rather large and in fact exceeds the upper limit set in . However, for the parameters $`B^{1/4}=235MeV`$ and $`\sigma =20MeV/fm^2`$ ($`\sigma ^{1/3}92MeV`$) there is a delay in the completion of the phase transition when the $`\gamma `$ term is included (long dashed) compared to the no $`\gamma `$ (dotted) case. The delay is of the order of 13 fm and is because of $`\beta `$ becoming greater than one a little below $`T_c`$ with consequent stoppage of all bubble nucleation. When the $`\gamma `$ term is not taken, there is no such stoppage and continous nucleation ensures faster completion. Figure 2b shows the same curves for the parameters $`B^{1/4}=235MeV`$, $`\sigma 7MeV/fm^2`$ and $`B^{1/4}=300MeV`$, $`\sigma 7MeV/fm^2`$. In the latter case, the curvature term completly inhibits bubble nucleation below $`T_c`$ and only the equilibrium bubbles formed above $`T_c`$ expand to complete the phase transition (long dashed). This corresponds to the case of no extremum in free energy (fig.2b). The inclusion of $`\gamma `$ term in this case actually reduces the time to complete the phase transition. This can be contrasted with the case corresponding to the other set of parameters in fig.1d and the cases shown in fig.2a. Comparing with the no curvature case, It is seen that increasing $`B`$ speeds up the phase transition whearas decreasing it slows it down if the $`\gamma `$ term is included. We also see that the decrease in the surface tension increases the supercooling of the plasma by significant amount (can be close to 50 percent). However, increasing $`B`$ reduces supercooling when $`\gamma `$ is included.

Next we look at the nucleation rate vs $`T/T_c`$. One can see from these graphs (figures 3a and 3b) that nucleation of hadronic bubbles can start even above the critical temperature when we include the curvature energy. The effect of decreasing $`\sigma `$ is to suppress nucleation. We may note that in the early Universe reverse is the case, and it happens because there the exponent drives the nucleation rate whearas in QGP the prefactor drives the nucleation. The point at which $`\beta `$ becomes greater than one is clearly visible when curvature term is included. In figure 3a the long dashed curve terminates at the point where further nucleation stops. It is also clear that nucleation starts the moment the plasma is created well above the critical temperature if the curvature term is included. Figure 3b shows that for smaller surface tension the $`\beta `$ term becomes greater than one soon after crossing through $`T_c`$ (solid and the long dashed curves). The other two curves shown are without the curvature terms and so begin below the critical temperature.

In figure 4 upper panel the critical radius of a bubble is plotted as a function of $`T/T_c`$ for those bubbles which are nucleated at the minima in the free energy. These bubbles in general are nucleated above $`T_c(r_{c+>})`$ as well as below $`T_c(r_c)`$. In fig.4a we see that for the parameters $`B^{1/4}=235MeV`$, $`\sigma =20MeV/fm^2`$ the bubbles keep nucleating from above the critical temperature to well below it. However for another set ($`B^{1/4}=300MeV`$, $`\sigma 7MeV/fm^2`$) these bubbles are no longer nucleated soon after crossing $`T_c`$. In figure 4b we find this feature in both the parameter sets. In all such situations (where $`\beta `$ becomes greater than 1) the completion of the phase transition is ensured by expansion of already nucleated bubbles. The lower panel plots the critical bubble radius $`r_{c+}`$ corresponding to the maxima in the free energy as a function of $`T/T_c`$ for temperatures below the critical temperature ($`T<T_c`$). We see in both the graphs 4c and 4d that the critical radius goes down with temperature. In figure 4c we see that the critical bubble radius goes to zero when $`\beta >1`$ (long dashed curve having $`\gamma 0`$). If no curvature term is present then the critical size does not go to zero. Even with curvature term, large $`\sigma `$ does not allow critical radius to go to zero (solid curve). Figure 4d plotted for a different set of parameters has similar features except that for $`B^{1/4}=300MeV`$, $`\sigma 7MeV/fm^2`$ the critical radius below $`T_c`$ is zero. The only bubbles that can expand are the ones already nucleated above $`T_c`$.

Finally, figures 5a and 5b show the hadron fraction $`h`$ as a function of time $`t`$ as the phase transition proceeds to it’s completion. Figure 5a shows the overlap for one parameter set given there (solid and dashed curves) and clearly shows that if we decrease the surface tension the phase transition proceeds slowly (long dashed curve) unlike the case in early Universe studies. The effect of $`\gamma `$ has a mixed behavior depending on the value of the bag pressure. Figure 5b shows that for low bag pressures, the effect of $`\gamma `$ is to delay the phase transition (solid curve) but for larger values of $`B`$ we see that the curvature correction actually speeds up the process of hadronization (long dashed curve).

The idealized Maxwell construction with the correction to the QGP as well as the the hadron gas that we have studied here shows that the time range of the phase transition is actually less than what has been given in . For the same parameter values as given in  we find that in the ideal case with corrections, the phase transition ranges from 2.5 fm to 22 fm whereas without the corrections it ranges from 3 fm to 37 fm. If we compare our study with this corrected ideal Maxwell construction, we find that for $`B^{1/4}=235MeV`$ the presence of the curvature term delays the transition to 31 fm for $`\sigma =50MeVfm^2`$, to 53 fm for $`\sigma =20MeVfm^2`$ and to more than 60 fm for $`\sigma 7MeVfm^2`$ . For $`B^{1/4}=300MeV`$ the ideal case streches from 1.9 fm to 6.8 fm but curvature introduces a delay and the transition extends to 12 fm for $`\sigma =50MeVfm^2`$, to 16 fm for $`\sigma =20MeVfm^2`$ and to 27 fm for $`\sigma 7MeVfm^2`$. It is interesting to note that the ideal Maxwell case range is 1.5 fm to 2.8 fm for $`B^{1/4}=350MeV`$, but $`\sigma 7MeVfm^2`$ with the curvature term produces a range from 1.5 fm to 16 fm.

## 4 Conclusions

In conclusion, we have presented in detail the mechanism of the phase transition of the QGP fireball from its creation to it’s subsequent hadronization in an ultrarelativistic nuclear collision. The inclusion of the interactions in both phases along with the the curvature energy term introduces subtlities that are both interesting and quantitative in nature. Broadly, we can summarize that whereas interactions lower the critical temperature, the combined effect of interaction and curvature energy is to affect both the degree of supercooling and the extent of the delay in the phase transition. Smaller values of the surface tension can enhance supercooling to almost 50 percent and also increase the delay in the hadronization process whereas larger values of the bag pressures speed up the hadronization process. $`\sigma ^{1/3}=65MeV`$ ($`\sigma 7MeVfm^2`$) is a reasonable value according to the literature and anything above $`\sigma ^{1/3}=105MeV`$ ($`\sigma 30MeVfm^2`$) seems to be on the higher side . The critical temperature is decided by the bag pressure and the inclusion of interactions in the two phases suppresses it. $`T_c170180MeV`$ is favored in the literature, so a bag pressure in excess of $`300MeV`$ seems to be on the higher side. Inclusion of the curvature effect has a mixed behavior in that, for larger values of $`B`$ the transition is speeded up whereas for smaller values of $`B`$ delay is enhanced. Another interesting observation is that the hadronization actually begins the moment the QGP is created, through equilibrium sized bubbles and is later joined by the expanding bubbles below the critical temperature. For reasonable parameter values the nucleation of bubbles below $`T_c`$ can get stopped when $`\beta >1`$ and now even the equilibrium sized bubbles join the expansion to complete the phase transition. Finally, for the ideal Maxwell construction, inclusion of interactions reduces the range of the transition but the presence of the curvature term introduces a delay in the completion of the transition compared to the ideal case.

Acknowledgements

This work was partially supported under the SERC Scheme of The Department of Science and Technology (DST), India. One of us (A.G.) is thankful to U.G.C. for partial financial support.

Figure Captions

Fig. 1. Free energy $`W`$ in units of MeV as a function of bubble radius $`r`$ in fm. In fig.1a Solid (with $`\gamma `$) and dashed (without $`\gamma `$) curves are for $`B^{1/4}=235MeV`$ and $`\sigma =50MeVfm^2`$. Long dashed (with $`\gamma `$) and dotted (without $`\gamma `$) curves are for $`B^{1/4}=235MeV`$ and $`\sigma =20MeVfm^2`$. In fig.1b Solid (with $`\gamma `$) and dashed (without $`\gamma `$) curves are for $`B^{1/4}=235MeV`$ and $`\sigma 7MeVfm^2`$. Long dashed (with $`\gamma `$) and dotted (without $`\gamma `$) curves are for $`B^{1/4}=300MeV`$ and $`\sigma 7MeVfm^2`$. Both the curves 1a and 1b are at a temperature of $`1.01T_c`$. In the lower panel the curves 1c and 1d are labelled as in Fig.1a and 1b respectively except that both are at a temperature of $`.99T_c`$.

Fig. 2. Temperature $`T/T_c`$ as a function of time $`t`$ in fm. In fig.2a and 2b Solid, dashed, dotted and long dashed curves are labelled as in fig. 1a and 1b respectively.

Fig. 3. Log of the nucleation time $`\tau `$(Tau) in units of fm/c as a function of temperature. The figures 3a and 3b are labelled as in fig.1a and 1b respectively.

Fig. 4. Upper panel shows critical bubble radius $`r_c`$ in fm as a function of $`T/T_c`$. They are the minimum free energy equilibrium critcal bubbles. The solid and long dashed curves are both with $`\gamma `$ correction. These bubbles do not exist without the curvature correction. Figures 4a and 4b are labelled as in fig.1a and 1b respectively. Lower panel shows critical bubble radius $`r_c`$ in fm as a function of $`T/T_c`$. They are the maximum free energy expanding critcal bubbles. Figures 4c and 4d are labelled as in fig.1a and 1b respectively.

Fig. 5. The hadron fraction $`h`$ as a function of time $`t`$ in fm. Curves in fig.5a and 5b are labelled as in fig. 1a and 1b respectively.


format: text

text:
# Crossover behaviors in liquid region of vortex states above a critical point caused by point defects

Also at \]Lyman Lab. of Physics, Harvard University, Cambridge, MA 02138

## Abstract

Vortex states in high-$`T_\mathrm{c}`$ superconductors with point defects are studied by large-scale Monte Carlo simulations of the three-dimensional frustrated XY model. A critical point is observed on the first-order phase boundary between the vortex slush and vortex liquid phases. A step-like anomaly of the specific heat is detected in simulations of finite systems, similar to an experimental observation \[F. Bouquet et al., Nature (London) 411, 448 (2001)\]. However, it diminishes with increasing system size, and the number and size distribution of thermally-excited vortex loops show continuous behaviors around this anomaly. Therefore, the present study suggests a crossover rather than a thermodynamic phase transition above the critical point.

Intensive studies on vortex states of high-$`T_\mathrm{c}`$ superconductors have revealed complex structures of the phase diagram in the magnetic field along the $`c`$ axis. Point defects such as oxygen deficiencies are essential for such complexity. The Bragg glass (BrG) Nattermann ; Giamarchi  phase is stable at low temperatures and low magnetic fields. This phase melts into the vortex liquid (VL) or vortex glass (VG) Fisher  phases as temperature $`T`$ or magnetic field $`B`$ increases, respectively. The BrG phase has a closed first-order boundary in the $`B`$-$`T`$ phase diagram. Ertas ; Giamarchi2  The vortex slush (VS) Worthington  phase locates between the VG and VL phases, and another first-order boundary between the VS and VL phases terminates at a critical point. Worthington ; Kierfeld

The BrG-VL, Safar92 ; Kwok94 ; Zeldov95 ; Welp96 ; Schilling96  BrG-VG Cubitt93 ; Safar95 ; Khaykovich97 ; Nishizaki98  and VS-VL Nishizaki00 ; Wen01  phase boundaries have been observed experimentally. Similar phase diagrams have been obtained by numerical simulations including point defects, Gingras ; Ryu ; Otterlo ; Sugano ; Nonomura  though the stability of the VG phase is still controversial. Olsson ; Vestergren

Recently Bouquet et al. Bouquet  observed a $`\delta `$-function peak and a sharp step of the specific heat below and above a critical point, respectively. They argued that this step-like anomaly corresponds to a second-order phase transition related to a vortex-loop blowout. Tesanovic ; Nguyen

In the present article, we investigate the phase diagram in the vicinity of the critical point of the VS-VL phase boundary. As displayed in Fig. 1, we find that the critical point is located at $`ϵ0.15`$, which is just above the maximum value of $`ϵ`$ in our previous work. Nonomura  A step-like anomaly of the specific heat is observed above the critical point in finite systems. However, the size dependence of this anomaly reveals that this step does not indicate a thermodynamic phase transition but a crossover.

In the present study we use the same model as investigated in our previous article on vortex states with point defects, Nonomura  namely the three-dimensional anisotropic, frustrated XY model on a simple cubic lattice, Li ; Hu

$``$ $`=`$ $`{\displaystyle \underset{i,jab\mathrm{plane}}{}}J_{ij}\mathrm{cos}\left(\varphi _i\varphi _jA_{ij}\right)`$ (1)
$`{\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)},𝑩=\mathbf{}\times 𝑨,`$ (2)
with the periodic boundary condition in all directions. A uniform magnetic field $`𝑩`$ is applied along the $`c`$ axis, and the averaged number of flux lines per plaquette is $`f`$. A point defect is introduced as the plaquette surrounded by four weaker couplings and distributed randomly in the $`ab`$ plane with probability $`p`$. Couplings are $`J_{ij}=(1ϵ)J`$ ($`0<ϵ<1`$) on the point defects, and $`J_{ij}=J`$ elsewhere. As in our previous article, we fix the parameters $`\mathrm{\Gamma }=20`$, $`f=1/25`$ and $`p=0.003`$, and vary the strength of point defects $`ϵ`$ and temperature $`T`$. Configurations of defects are also the same as in our previous work. Monte Carlo simulations are started from a very high temperature, and systems are gradually cooled down. Typical Monte Carlo steps are $`414\times 10^7`$ for equilibration, and $`216\times 10^7`$ for measurement.

First, we concentrate on the system with $`L_a=L_b=50`$ and $`L_c=40`$. Temperature dependence of the specific heat $`C`$ for $`0.14ϵ0.2`$ is displayed in Fig. 2(a). At $`ϵ=0.14`$, $`C`$ has a $`\delta `$-function peak at $`T_{\mathrm{sl}}=0.0760J/k_\mathrm{B}`$. This peak becomes very small at $`ϵ=0.15`$, and is replaced by a step for $`ϵ0.16`$. This result suggests that the critical point locates at $`ϵ0.15`$. The phase diagram obtained from this finite system is summarized in Fig. 1.

Consistent behaviors are seen in other quantities. Temperature dependence of the density of dislocations in the $`ab`$ plane $`\rho _\mathrm{d}`$ and the ratio of point defects which trap flux lines $`r_{\mathrm{tr}}`$ are given in Figs. 2(b) and 2(c), respectively. The quantity $`\rho _\mathrm{d}`$ characterizes Kierfeld ; Nonomura  the VS-VL phase transition. At $`ϵ=0.14`$, $`\rho _\mathrm{d}`$ jumps between $`T=0.072J/k_\mathrm{B}`$ and $`0.076J/k_\mathrm{B}`$, where $`C`$ sharply increases. This jump becomes much smaller at $`ϵ=0.15`$ and invisible for $`ϵ0.16`$. Keeping possible statistical errors in mind, a tiny hump and dip seem to exist at the temperature region $`T^{}`$ where the specific heat has a step-like anomaly for $`ϵ0.15`$, but the hump and dip are one order smaller than the jump of $`\rho _\mathrm{d}`$ related to the VS-VL phase transition, and therefore considered to be finite-size effects.

At $`ϵ=0.14`$, the quantity $`r_{\mathrm{tr}}`$ gradually increases from $`T=0.070J/k_\mathrm{B}`$ as approaching the transition temperature $`T_{\mathrm{sl}}`$, because the short-range order of flux lines in the VS phase becomes fragile as temperature increases. Then, it rapidly decreases just below $`T_{\mathrm{sl}}`$, because the mobility of flux lines becomes larger in the coexisting region of the first-order phase transition. This anomaly becomes obscure at $`ϵ=0.15`$ and disappears for $`ϵ0.16`$. Although a tiny hump and dip also seem to exist at $`T^{}`$ for $`ϵ0.15`$, they have an opposite behavior to the cusp at $`ϵ=0.14`$; $`r_{\mathrm{tr}}`$ increases at $`T^{}`$ and drops just above $`T^{}`$, and therefore have a different origin from the cusp related to the VS-VL phase transition.

Next, in order to explore whether this step-like anomaly of the specific heat is a thermodynamic phase transition or not, its size dependence is investigated. Since flux lines in $`𝑩\widehat{𝒄}`$ are one-dimensional objects along the $`c`$ axis, physical quantities mainly depend on $`L_c`$. Nonomura2  Even in the vortex liquid region in pure systems, the displacement of flux lines along the $`ab`$ plane is smaller than the random-walk value $`L_c^{1/2}`$ due to the blocking by other flux lines. In the present case, this displacement would be even smaller due to the pinning by defects. We thus fix $`L_a=L_b=50`$ and change $`L_c`$. Temperature dependence of the specific heat for $`L_c=40`$, $`56`$ and $`80`$ at $`ϵ=0.2`$ is displayed in Fig. 3(a). As $`L_c`$ increases, the height of the step of the specific heat becomes smaller. Since it is clearly proportional to $`1/L_c`$ (Fig. 3(b)), this step is expected to vanish in the thermodynamic limit. The anomaly of the specific heat above the critical point of the VS-VL phase boundary is therefore a crossover rather than a thermodynamic phase transition.

Finally, the number and size distribution of thermally-excited vortex loops are investigated in order to check whether this anomaly of the specific heat is related to the vortex-loop blowout Tesanovic ; Nguyen  or not. Temperature dependence of the number of vortex loops for $`0.14ϵ0.2`$ is displayed in Fig. 4(a). There exist no vortex loops including pancake vortices at these parameter regions. Size distributions of Josephson loops for $`ϵ=0.14`$ and $`0.2`$ are shown in Figs. 4(b) and 4(c), respectively. Note

Although the number of vortex loops monotonically increases as temperature increases, no symptom of the vortex-loop blowout is seen around $`T_{\mathrm{sl}}`$ or $`T^{}`$ for each $`ϵ`$. This is also the case for the size distribution of Josephson loops. In Figs. 4(b) and 4(c), the number of Josephson loops $`N_\mathrm{J}(L_\mathrm{J})`$ with a given perimeter $`L_\mathrm{J}`$ is plotted versus $`L_\mathrm{J}`$ in a semi-log scale, and the exponential size distribution of $`N_\mathrm{J}(L_\mathrm{J})`$ is clearly observed. There exists no qualitative difference in the behaviors below and above the critical point, and the size distribution changes continuously as temperature is varied around $`T_{\mathrm{sl}}`$ or $`T^{}`$.

The vortex-loop blowout is characterized by the abrupt increase of the number of vortex loops and the change of the size distribution from an exponential decay to a power-law one at the transition temperature, and the proliferation of vortex loops as large as the system size above the transition temperature. Tesanovic ; Nguyen  None of these behaviors are observed around $`T^{}`$ in the present simulations, which casts serious doubt on the argument that the step of the specific heat is related to the vortex-loop blowout. Bouquet  In this argument, the vortex-loop blowout scenario was borrowed in order to explain the observation that the value of the specific heat above the step-like anomaly $`C_+`$ is larger than that below the anomaly $`C_{}`$. However, this property is quite general in vortex states, and actually occurs even at the melting transition in pure systems. Hu  In vortex states, the density of vortices increases as temperature increases owing to excitation of vortex loops and wandering of flux lines, and consequently the relation $`C_+>C_{}`$ holds. The vortex-loop blowout is merely a special case of this general property.

It remains unsolved why a step of the specific heat is clearly observed in experiments even though the sample size is much larger than the system size in our simulations. Several possibilities are: (a) The height of the step of $`C`$ might saturate at larger $`L_c`$, and our maximum system size $`L_c=80`$ would be still small. (b) Our model (1) might be still too simple to describe the step of $`C`$, e.g. the inhomogeneity of the strength of point defects is not taken into account. (c) The step of $`C`$ observed experimentally might be a non-equilibrium effect.

In summary, vortex states of high-$`T_\mathrm{c}`$ superconductors with point defects in $`𝑩\widehat{𝒄}`$ has been investigated by large-scale Monte Carlo simulations. The critical point of the first-order phase boundary between the vortex slush and vortex liquid phases is identified as the vanishing point of the $`\delta `$-function peak of the specific heat. The shape of the step-like anomaly of the specific heat above the critical point is similar to that of the experiment by Bouquet et al. As the system size along the $`c`$ axis $`L_c`$ increases, the height of this step seems proportional to $`1/L_c`$, which suggests that this anomaly is a crossover rather than a thermodynamic phase transition. Consistent behaviors are observed in the temperature dependence of the density of dislocations in the $`ab`$ plane and the ratio of point defects which trap flux lines. There is no discontinuity in the number and size distribution of vortex loops around the step-like anomaly of the specific heat, which clearly indicates that this anomaly is not related to the vortex-loop blowout.

Numerical calculations were performed on Numerical Materials Simulator (NEC SX-5) at Computational Materials Science Center, National Institute for Materials Science, Japan. Y. N. was supported by the Atomic Energy Research Fund from Ministry of Education, Culture, Sports, Science and Technology, Japan. This study is partially supported by the Priority Grant No. 14038240 from the same Ministry.


text:
# Filament bifurcations in a one-dimensional model of reacting excitable fluid flow

## 1 Introduction

Chemically reacting substances in fluid flows are complex systems important both from a fundamental point of view and for its relevance in industrial and environmental contexts . Typical incompressible chaotic flows stretch fluid parcels along particular directions, with the consequent contraction along the others. Thus, passive markers draw filamental or lamellar structures that have been investigated both from theory as from experiment. Chemical reactants are also stretched in this way, producing a great increase in the surface of contact between different species with deep effects on the global chemical kinetics. In addition to the strictly chemical interactions, biological interactions in ecosystems (predation, grazing, consumption, competition, …) can also be described in the same framework, so that stability and dynamics of aquatic ecosystems are also strongly influenced by this filamentation process.

In this Paper we analyze a simplified model describing the transverse chemical structure of these filaments, for the case in which the chemical or biological activity is of the excitable type. This includes in particular well known laboratory reactions, such as the Belousov-Zhabotinsky reaction, but also reactions of environmental importance, such as phytoplankton-zooplankton competition. Such kind of dynamics in two-dimensional chaotic flows was analyzed in . There it was found that some of the dynamic regimes and transitions between them can be understood in terms of a simple one-dimensional model, of the kind already considered in  that focusses in the transverse structure of the filaments. It was mentioned in  that, in addition to the bifurcations considered there in detail, there was a range of parameters for which a complex coexistence of solutions and bifurcations occurred, that may be of relevance to understand qualitative behaviors of the two-dimensional chaotic flows. The present paper focusses in that regime. We use the FitzHug-Nagumo model as a prototypical excitable dynamics, but we expect that the same qualitative results will be also obtained for other chemical or biological excitable schemes. Similar behavior has been observed in an excitable model of oceanic plankton population in .

## 2 The one-dimensional filament model

Given a set of kinetic laws describing the time evolution of a number of chemical or biological concentrations $`\{C_i(𝐱,t)\}_{i=1,\mathrm{},N}`$ in a homogeneous system,

$$\frac{dC_i}{dt}=kF_i(C_1,\mathrm{},C_N)$$
(1)
($`k`$ is a global reaction rate), a general partial differential equation model describing their spatiotemporal evolution under the simultaneous additional effects of advection in an incompressible velocity field $`𝐯(𝐱,t)`$, and diffusion with diffusion coefficients $`\{D_i\}`$, is the one given by the following advection-reaction-diffusion equations:

$$\frac{C_i(𝐱,t)}{t}+𝐯(𝐱,t)C_i(𝐱,t)=kF_i(C_1,\mathrm{},C_N)+D_i^2C_i(𝐱,t).$$
(2)
In the following we restrict for simplicity to the case of equal diffusion coefficients $`D_i=D,i`$. In the immediate vicinity of filaments or lamelae, one can approximate the flow by its linearization around a point co-moving with a fluid element. In areas where the advection dynamics is hyperbolic this linearization identifies principal directions $`\{x_j\}`$ along which the corresponding velocity components read $`v_j=\lambda _jx_j`$. Positive and negative values of the strain rate $`\lambda _j`$ identify expanding and contracting directions, respectively. Along the expanding directions, diffusion and advection cooperate and the chemical concentrations are fastly homogenized, whereas strong gradients build-up in the contracting directions. This is the origin of filamental (one expanding direction) or lamellar (two expanding directions) structures. Gradients in Eq. (2) can be safely neglected except along the contracting directions. We consider here the case of filaments in two-dimensional flows, or lamellae in three dimensions, so that there is only one contracting direction, of strain rate $`\lambda _i\lambda `$, with $`\lambda >0`$. In this case, (2) reduces to an effective one-dimensional model for the transverse profile of the chemical distributions. By measuring time in units of $`\lambda ^1`$, and space in units of $`\sqrt{D/\lambda }`$, we arrive at:

$$\frac{C_i(x,t)}{t}x\frac{}{x}C_i(x,t)=\mathrm{Da}F_i(C_1,\mathrm{},C_N)+\frac{^2}{x^2}C_i(x,t).$$
(3)
We have introduced a Damköhler number $`\mathrm{Da}=k/\lambda `$ as the ratio between the chemical and the strain rate.

Several limitations are inherent to (3). It is a strictly local model valid close to (moving) hyperbolic trajectories. The linearization assumption will only be justified if the full filament width is smaller than any structure in the velocity field. This implies in particular a very small diffusion coefficient so that we are deeply in the Batchelor regime. Additional qualitative changes of behavior, not contained in (3), are expected when increasing the diffusion coefficient in (2). Curvature effects or, more importantly, any interaction between filaments, are completely neglected. Finally, the use of a constant $`\lambda `$ neglects any fluctuations of the stretching rate.

The FitzHugh-Nagumo (FN) model consists in a dynamics of the type (1) for two interacting species, of concentrations $`C_1`$ and $`C_2`$, and reaction terms:

$`F_1`$ $`=`$ $`C_1(aC_1)(C_11)C_2,`$ (4)
$`F_2`$ $`=`$ $`ϵ(C_1\gamma C_2).`$ (5)
The FN model behaves excitably when $`ϵ1`$ so that there is a separation between the fast evolution of the active component or activator, $`C_1`$, and the slow evolution of $`C_2`$, the passive one or inhibitor. We concentrate in the parameter values $`ϵ=10^3`$, $`\gamma =3.0`$, and $`a=0.25`$, for which robust excitable behavior is obtained.

## 3 Filament solutions

The unexcited solution, $`C_1=C_2=0`$, is a linearly stable exact solution of system (3)-(5). In addition, the most notable solution is a pulse like steady solution, centered on $`x=0`$, in which the activator is fully excited in a central region, whereas the inhibitor remains at reduced concentrations. This chemical distribution can be thought as the transverse cut of one of the excited filaments seen in two-dimensional excitable fluid flows. Pulse solutions of this type are shown in Fig. 1.

The width, $`w_s`$, of this solution is determined by the competition between the tendency to expand, due to the joint effects of the excitable chemistry and diffusion, and the compression by the flow, and can be estimated  to be, for small $`ϵ`$, $`w_s(12a)\sqrt{\mathrm{Da}/2}`$. It is very similar to what is obtained in one-component bistable chemical models, so that the presence of the inhibitor plays here only a secondary rôle.

As in bistable chemical models, this stable pulse solution coexists in a range of parameters with a unstable pulse, of width given at small $`ϵ`$ approximately by  $`w_u2/\sqrt{a\mathrm{Da}}`$. When the value of $`\mathrm{Da}`$ decreases (meaning that the chemistry becomes slower, or the stretching faster) the two solutions approach, collide and disappear from phase space in a saddle-node bifurcation . This happens at a value of $`\mathrm{Da}_c2\sqrt{2/a}(12a)^1`$. Thus, the excited pulse solution does not exists at small $`\mathrm{Da}`$, and the existence of a critical $`\mathrm{Da}_c`$ explains  the dynamical transition occurring in closed two-dimensional chaotic flows between a situation in which local perturbations of the unexcited state have limited impact on the system, and a state in which they give rise to a global excitation of the whole fluid. A related transition occurs in open flows .

At larger values of $`\mathrm{Da}`$, the influence of the inhibitor becomes more noticeable and, as a result, steady filament solutions have a ‘two-humped’ shape (Fig.2a). It was noted in  that this kind of ‘double-filament’, also seen in two-dimensional flow simulations, can be regarded as a bound state of two asymmetric counter-propagating excitable pulses, which are also steady solutions of model (3)-(5) coexisting with the symmetric ones at large $`\mathrm{Da}`$ (Fig. 2b). Of course, for each value of $`\mathrm{Da}`$ there are two asymmetric solutions of this kind, each one being the specular image of the other.

## 4 Bifurcation behavior at intermediate Da

It was noted in  that a complex bifurcation scenario occurred between the small and the large values of $`\mathrm{Da}`$ for which the solutions in Figs. (1) and (2) were respectively obtained. It happens that these three solution branches are not directly connected, at least for the values of $`ϵ`$ and $`\gamma `$ considered here. We have followed the three branches of steady stable pulse solutions by starting with well developed solutions at large and small $`\mathrm{Da}`$, and then performing slow changes in $`\mathrm{Da}`$. In this way we can only follow stable solutions. A more complete characterization would require also the continuation of unstable branches, that will be performed elsewhere. We monitor the height of the activator at $`x=0`$, and plot it in Fig. 3.

When increasing the value of $`\mathrm{Da}`$ on the symmetric one-humped solution branch, the inhibitor value increases in the center, so that for $`\mathrm{Da}45`$, the central concentration of the activator becomes a shallow minimum (Fig. 1b). So, above this parameter value, this branch can be called ‘slightly two-humped’, and it remains so until $`Da77.1`$. It is however clearly different from the ‘strongly two-humped’ solution branch to which Fig. 2a pertains. At $`Da77.1`$, a Hopf bifurcation occurs: the filament (both the $`C_1`$ and the $`C_2`$ concentrations) pulsates in height, width and shape. We show in Fig. (4) this behavior. For $`\mathrm{Da}`$ above this Hopf bifurcation, we plot in Fig. 3 the maximum and minimum central height values attained during the oscillation. At $`\mathrm{Da}81.0`$ the width becomes too narrow at some moment of the oscillation and the filament collapses to the unexcited ($`C_1=C_2=0`$) solution, so that the limit cycle solution disappears. This collapse probably reveals the collision of the filament limit cycle with an unstable pulse solution that we have not characterized. The basin of attraction of the oscillating solution is rather narrow: at these values of $`\mathrm{Da}`$ it is easier to be attracted by the two-hump symmetric steady solution branch or the unexcited state.

On the other hand, if we start with a well developed two-hump filament such as the one in Fig. 2a at high $`\mathrm{Da}`$, the two humps approach each other when decreasing $`\mathrm{Da}`$. Before they fuse together, there is a discontinuous jump (at $`\mathrm{Da}75.4`$) to the steady ‘slightly two-humped filament’ followed before. This probably reveals a saddle-node bifurcation in which the stable two-hump branch collides with an unstable one, that we have not followed.

Finally, starting from the asymmetric filament of Fig. 2b, it approaches the $`x=0`$ axis when decreasing $`\mathrm{Da}`$, becoming more and more symmetric. For the values of $`ϵ`$ and $`\gamma `$ used here, however, it does not join smoothly with the symmetric branch in a forward pitchfork bifurcation in which it would also join the other asymmetric filament of opposite symmetry. Rather, it performs a (small) discontinuous jump to the symmetric filament branch at $`\mathrm{Da}70.48`$. This is probably the signature of a backward pitchfork bifurcation, or of some other complex coexistence with unstable solutions, that would also explain the observed range of bistability between the symmetric and asymmetric filament solutions.

## 5 Summary

We have investigated several pulse-like solutions of a one-dimensional model intended to characterize the chemical structure transverse to filaments or lamelae generated by chaotic flows in reacting systems. The most striking feature has been the finding of an oscillating solution, in which an activator-inhibitor competition is sustained forever by the compressing strain. Moreover, parameter ranges have been found in which there is coexistence of up to three stable pulse solutions, in addition to the unexcited one. The transition behavior observed in two-dimensional open and closed chaotic flows at intermediate $`\mathrm{Da}`$ numbers  is probably related to this complex phase space. Further work is needed to fully characterize this relationship.

## Acknowledgments

We acknowledge support from MCyT (Spain) projects BFM2000-1108 (CONOCE) and REN2001-0802-C02-01/MAR (IMAGEN).


text:
# Precise estimation of shell model energy by second order extrapolation method

## Abstract

A second order extrapolation method is presented for shell model calculations, where shell model energies of truncated spaces are well described as a function of energy variance by quadratic curves and exact shell model energies can be obtained by the extrapolation. This new extrapolation can give more precise energy than those of first order extrapolation method. It is also clarified that first order extrapolation gives a lower limit of shell model energy. In addition to the energy, we derive the second order extrapolation formula for expectation values of other observables.

A pursuit for new numerical algorithms to solve nuclear shell model problems has been one of the most intriguing issues in nuclear structure physics. In the last decade, several new methods smmc ; dmrg ; vampir ; qmcd1 ; qmcd-rev ; txp ; ecm ; xp  were proposed and have been developed. Consequently shell model calculation remarkably widens its capability beyond the limitation of conventional approaches such as Lanczos shell model diagonalization. Among these new methods, several extrapolation techniquestxp ; ecm ; xp  were proposed in order to reinforce a restricted power of Lanczos shell model calculations. We advocated the Lanczos shell model diagonalization with the help of the recently proposed non-empirical extrapolation technique imada1 ; sorella  for correlated electron systems and succeeded in enlarging the feasibility of large-scale shell model calculations by the extrapolation xp .

In general, interpolation can be reliable if the considered quantity behaves smoothly, while extrapolation cannot be guaranteed. A reliable extrapolation method needs a theory predicting how the considered quantity behaves in the extrapolated region. In the recently proposed extrapolation imada1 ; sorella , it has been shown that, exact energies defined in the whole space can be extrapolated from series of energies in truncated Hilbert spaces. If these truncated Hilbert spaces give fairly good approximations, energies of truncated spaces can be shown to behave linearly as a function of the energy variance. Therefore exact energies can be estimated by linearly extrapolated energies of truncated spaces into zero energy variance. In condensed matter physics there are good examples of this extrapolation in analyses of strongly correlated electrons on a lattice imada2 ; imada3 . In Ref. xp  it has been shown that this extrapolation works also when applied to the nuclear shell model. However, for detailed spectroscopic studies, precise excitation energies are important, and this requires a more precise evaluation of the energy eigenvalues. Another problem is, as it will be discussed later, that the extrapolated shell model energies are systematically lower than the true energy eigenvalues. To solve this problem, in the present paper, we successfully improve the extrapolation method by including the second order term. Hereafter this extension is called second order extrapolation, while the previous one is called first order extrapolation.

In order to present our new method, we consider the relation between energy difference and modified energy variance as below. We define the difference $`\delta E`$ between the lowest energy eigenvalue $`\widehat{H}`$ in a given truncated space and true energy eigenvalue $`\widehat{H}_0`$ as

$$\delta E=\widehat{H}\widehat{H}_0.$$
(1)
The modified energy variance $`\mathrm{\Delta }E`$ in the truncated space is also defined as

$$\mathrm{\Delta }E=\widehat{H}^2\widehat{H}^2.$$
(2)
Note that this definition of energy variance is different from the previous one xp . Throughout the present paper, this new definition is used.

An approximate ground state $`|\psi `$ can be decomposed as the sum of the true eigenstate $`|\psi _0`$ and the rest component $`|\psi _r`$ as

$$|\psi =c|\psi _0+d|\psi _r,$$
(3)
$$|\psi _r=\underset{n0}{}c_n|\psi _n,$$
(4)
where $`c`$, $`d`$ and $`c_n`$’s are coefficients, whose normalization condition is taken as $`c^2+d^2=1`$ and $`\underset{n0}{}c_n^2=1`$. The $`|\psi _0`$ and $`|\psi _n`$’s form an orthonormarized complete set spanning the whole space and its energy eigenvalue is denoted by $`E_n`$. By defining the moments $`D_j`$ ($`j=1,2,`$…) as

$$D_j=\underset{n0}{}\frac{c_n^2\left(E_nE_0\right)^j}{E_0^j},$$
(5)
we can rewrite $`\delta E`$ and $`\mathrm{\Delta }E`$ as

$$\delta E=d^2D_1E_0,$$
(6)
and

$$\mathrm{\Delta }E=d^2D_2E_0^2\left(d^2\right)^2\left(D_1E_0\right)^2.$$
(7)
By eliminating $`d^2`$, we obtain

$$\delta E=\frac{1}{2A}\left[1\sqrt{14A^2\mathrm{\Delta }E}\right],$$
(8)
where

$$A=\frac{D_1}{D_2E_0}.$$
(9)
If $`4A^2\mathrm{\Delta }E`$ is small enough compared to the unity, we can expand $`\delta E`$ as

$$\delta E=A\mathrm{\Delta }E+A^3\mathrm{\Delta }E^2+\mathrm{}.$$
(10)
As the ratio of the moments $`A`$ embodies the structure of the rest component, $`A`$ can be a function of $`d^2`$ which is a small parameter. Here we assume that $`A`$ can be expanded by a function of $`d^2`$ as

$$A=A_0+A_1^{}d^2+O(d^4),$$
(11)
where $`A_0`$ and $`A_1^{}`$ are coefficients. This expansion is shown to be reasonable in the later discussion. By eliminating $`d^2`$ by Eq. (7), $`A`$ can be rewritten as

$$A=A_0+A_1\mathrm{\Delta }E+\mathrm{},$$
(12)
where $`A_1`$ is a new coefficient for $`A_1^{}`$ with some additional factors. Therefore, we obtain the following second order extrapolation formula,

$$\delta E=A_0\mathrm{\Delta }E+\left(A_0^3+A_1\right)\mathrm{\Delta }E^2+\mathrm{}.$$
(13)
The expansion parameter of Eqs. (12) and (13) is $`\mathrm{\Delta }E/E_0^2`$ but, as $`E_0`$ is unknown, we present a new extrapolation formula by including the $`E_0`$ depedence in the expansion coefficients. We comment that, a dimensionless definition of energy variance as $`\stackrel{~}{\mathrm{\Delta }}E=\mathrm{\Delta }E/\widehat{H}_0^2`$ leads to the same result.

The first term of Eq. (13) shows the effect of the decrease of the norm of rest component. If we consider only the first term, this extrapolation is nothing but the previous first order extrapolation xp . The second term shows the second order effects of the decrease of the norm of rest component ($`A_0^3`$) and the first order effect of the structure change ($`A_1`$). The sign of the second order coefficient will be clarified later in the discussion. This new method allows us to extrapolate shell model energies in the limit $`\mathrm{\Delta }E0`$ after determining the expansion coefficients in Eq. (13) by fitting to the shell model energies of various truncation spaces.

In order to investigate the properties of the present second order extrapolation method, we utilize numerically solvable shell model calculations. One frontier of the state-of-the-art large-scale shell model calculations is the $`fp`$ shell. Here we take ground state energies of 48Cr and 49Cr as $`fp`$ shell nuclei. For a numerical test, we consider two cases of different nature. The $`M`$-scheme dimensions are about 2 and 5 Millions for 48,49Cr, respectively, and we can easily diagonalize the Hamiltonian matrices by a shell model code, for instance, MSHELL mshell . As an effective shell model interaction, the KB3 interaction kb3  is taken.

We consider two kinds of truncation scheme. The first is $`_{st}(f_{7/2})^{A40s}(r)^s`$ where $`r`$ means the set of the $`f_{5/2}`$, $`p_{3/2}`$ and $`p_{1/2}`$ orbits and $`t`$ is the maximum number of particles allowed to be excited into the $`r`$ orbits. This truncation scheme is natural due to a considerable shell gap between the $`f_{7/2}`$ orbit and the others. We call this truncation scheme (I). The other is the truncation scheme (II): $`_{s+s^{}t_0,s^{}t}(f_{7/2})^{A40ss^{}}(p_{3/2})^s(r^{})^s^{}`$ where $`r^{}`$ means the set of the $`f_{5/2}`$ and $`p_{1/2}`$ orbits. Truncation spaces are specified by $`t`$ and $`t_0`$. For 48,49Cr, the $`t_0`$ is set to 5, which means up to 5 nucleons in the $`p_{3/2}`$ and $`r^{}`$ orbits. The parameter $`t`$ controls the size of the shell model space by specifying the maximum number of nucleons in the $`r^{}`$ orbits. These two sets of truncation scheme gradually cover the whole shell model space in different ways as $`t`$ increases.

In Figs. 1 (A) and (C), the energies for the truncated shell model spaces with $`t=25`$ for the truncation scheme (I) and $`t=15`$ for the truncation scheme (II) are shown as a function of the energy variance $`\mathrm{\Delta }E`$ for the ground states of 48Cr and 49Cr. The solid (I) and dashed (II) curves show the second order extrapolations. Independently from the nature of the nuclei and from the adopted truncation scheme, the shell model energies for truncated spaces are well-fitted by quadratic curves as indicated by the second order extrapolation formula. These extrapolated energies are quite accurate, for examples, for 48Cr, they are -32.91 and -32.97 MeV for (I) and (II), respectively, while exact energy is -32.95 MeV. For 49Cr, extrapolated energies are -35.59 and -35.60 MeV for (I) and (II), respectively, while exact one is -35.59 MeV. Different truncation schemes can give almost the same extrapolation energy. This fact is quite useful for confirming the reliability of the extrapolated energy.

In Fig.2, we compare the results of the first order extrapolation and of the second order extrapolation for the ground state of 49Cr, where the truncation scheme (I) is taken. L1 and L2 are the extrapolations from the shell model energies for $`t=2,3,4`$ and $`t=3,4,5`$, respectively. L2 gives higher extrapolation energy than L1, while it is still lower than the exact energy. However, globally the shell model energies for truncated spaces can be better fitted by a quadratic curve than by straight lines. The first order extrapolation is approximately parallel to the average tangential line of this quadratic curve around the fitted shell model energies. If the quadratic curve is upwardly concave, that is $`A_0^3+A_1>0`$, the extrapolated energy becomes a lower limit of shell model energy and approaches the exact energy from the lower side.

We further go into the details of the second order extrapolation formula. As we can obtain exact ground state wavefunctions for 48,49Cr, the moments $`D_1`$ and $`D_2`$ can also be evaluated. In Figs. 1 (B) and (D), the ratio of the moments $`A`$ is plotted as a function of $`\mathrm{\Delta }E`$. These figures show that $`A`$ can be well described by the linear relation of $`\mathrm{\Delta }E`$. The coefficient $`A_0`$ is almost the same for different truncation schemes, while $`A_1`$ is different. For 48Cr, the moment ratios are $`A=0.037+0.0018\mathrm{\Delta }E`$ for (I) and $`A=0.036+0.0010\mathrm{\Delta }E`$ for (II). In the second order coefficient, $`A_1`$ is positive for both cases. Moreover, $`A_1`$ is much larger than $`A_0^3`$. This means that it is important to take into account the structure change of the rest component in second order term.

As the moments $`D_1`$ and $`D_2`$ are positive by definition, the ratio of the moments $`A`$ is positive then $`A_0`$ is also positive. The sign of the second order coeficient of Eq. (13) depends on $`A_1`$. Here we explain the condition $`A_1>0`$. We define the energy of the rest component relative to the true ground state energy, $`E_r=\psi _r\left|H\right|\psi _rE_0`$. The first moment $`D_1`$ is proportional to $`E_r`$, while the second moment $`D_2`$ is approximately proportional to $`E_r^2`$. Therefore their ratio $`A`$ has a $`1/E_r`$ dependence. If $`E_r`$ increases as $`t`$ increases, $`A`$ decreases as a function of $`t`$, which means that $`A_1>0`$. In the following we consider two physically relevant but different truncation schemes that, however, share a common way for expanding shell model spaces as $`t`$ increases. The shell model orbits are divided into higher and lower orbits. As the parameter $`t`$ is the maximum nucleon number in the higher orbits, we expand truncated spaces by including basis states with larger spherical single particle energies on average. By this construnction of the series of truncation spaces, as $`t`$ increases, the rest component is described by basis states with larger spherical single particle energies and $`E_r`$ becomes larger. Therefore upwardly concave feature is generally expected in energy - energy variance plot.

Next we consider the expectation value of an observable $`\widehat{O}`$, where a similar second order extrapolation formula holds. In the same way as the energy, we define the difference as

$$\delta O=\widehat{O}\widehat{O}_0,$$
(14)
that can be rewritten by imada2

$$\delta Od^2D_1^OE_0,$$
(15)
in which $`D_1^O`$ is defined as

$$D_1^O=\underset{n,m0}{}\frac{c_nc_m\left(O_{nm}O_{00}\delta _{nm}\right)}{E_0},$$
(16)
where $`O_{nm}=\psi _m\left|\widehat{O}\right|\psi _n`$. By eliminating $`d^2`$, $`\delta O`$ can be expanded by the ratio $`\mathrm{\Delta }E/E_0^2`$ as

$$\delta O=B\mathrm{\Delta }E+\left(BA^2\right)\mathrm{\Delta }E^2+\mathrm{},$$
(17)
with

$$B=\frac{D_1^O}{D_2E_0}.$$
(18)
We assume that $`B`$ can be expanded as

$$B=B_0+B_1\mathrm{\Delta }E+\mathrm{},$$
(19)
where $`B_0`$ and $`B_1`$ are coefficients. Then, we obtain the second order formula for $`\delta O`$ as

$$\delta O=B_0\mathrm{\Delta }E+\left(B_0A_0^2+B_1\right)\mathrm{\Delta }E^2+\mathrm{}.$$
(20)
As a further large-scale numerical test, we investigate the ground state of 64Zn, whose $`M`$-scheme dimension is about 501 Million and its exact shell model calculation is quite difficult. However, exact results are reported in Ref. caurier-fp .

In Fig. 3, the results of extrapolations of energies (A), $`B(E2)`$ (B) and Gamow-Teller (C) sum rules for the $`0_1^+`$ state are shown for 64Zn. We consider two truncation schemes. For the truncation scheme (I), we plot the results for $`t=1`$ $``$ 4 spaces, whose $`M`$-scheme dimension is 0.2 $``$ 55 Millions. For the truncation scheme (II), we plot the results for $`t=4`$ $``$ 8 with $`t_0=4`$, whose $`M`$-scheme dimension is 0.2 $``$ 40 Millions. Shell model energies of truncated spaces can be well-fitted by the quadratic curves for both the truncation schemes. Extrapolated energies are -106.38 and -106.34 MeV, for the truncation schemes (I) and (II), respectively. The exact ground state energy is -106.3 MeV caurier-fp . The present extrapolated energies are quite close to the exact one. The extrapolated values of the total B(E2) are 1.17 $`\times 10^3`$ and 1.18 $`\times 10^3`$ $`e^2fm^4`$ for the truncation schemes (I) and (II), respectively, while the exact value is 1157 $`e^2fm^4`$ caurier-fp  with the effective charges as $`e_\pi =1.35e`$ and $`e_\nu =0.35e`$. Here we use the second order extrapolation formula for the expectation values. The extrapolated values of the total B(GT+) with unquenched Gamow-Teller operators are 5.3 and 5.5 for the truncation schemes (I) and (II), respectively, while its exact value is 5.54. The present extrapolations with two truncation schemes can give almost exact values.

In summary, we have proposed a second order extrapolation method for shell model calculations. Shell model energies of the truncated spaces as a function of $`\mathrm{\Delta }E`$ are well described by a quadratic curve and the exact shell model energy is extrapolated in the limit $`\mathrm{\Delta }E0`$. To derive this second order extrapolation formula, it is quite important to consider the second order term of the $`d^2`$ expansion and the structure change of rest components in terms of the ratio of the moments $`A`$. By numerical solvable shell model calculations, we have compared the extrapolated energies to the exact one and showed that their accuracy is improved by the present second order extrapolation and that the first order extrapolation gives a lower limit of the shell model energy. In addition, we have extended the derivation of second order extrapolation to any observable. Finally, we have shown that different truncation schemes have different behavior in extrapolation, but, in principle, extrapolated values must be the same. This fact is used to confirm the accuracy of extrapolated values.

We thank Professor N. Yoshinaga and Dr. M. Serra for reading the manuscript and delightful comments. This work was supported from ’Research for the Future Program’ by Japan Society for the Promotion of Science under grant RFTF97P01103, and was supported in part by Grant-in-Aid for Scientific Research (A)(2) (10304019), and one for Specially Promoted Research (13002001) from the Ministry of Education, Science and Culture.