text:
# 1 INTRODUCTION

## 1 INTRODUCTION

The study of phase changes in macroscopic systems has a long tradition in statistical physics. Its basic assumption is that the dimensions of the macroscopic system are very large when compared with that of the constituting elements: phase transitions are well-defined only for infinite systems. However, there are many phenomena in finite systems which resemble phase transitions, see, for instance, Refs. . Because of their importance for the understanding of the physics of clusters of atoms  or the folding of proteins and other biological molecules , to name only a few examples, these “phase transitions” in small systems have recently attracted renewed interest. The main question is how the observed effects in small systems can be related to true phase transitions in macroscopic (or infinite) systems. A few attempts were recently made in this direction, either through studying for finite and small systems the topology of curvatures of the entropy-density surface $`s(e,n)`$  in the microcanonical ensemble, or by exploring the density of complex zeros of their canonical partition function .

The later approach is closer to the traditional view on phase transitions. For infinite systems the physical features of a phase transition can be obtained from the distribution of the complex zeros of partition functions of finite systems. As the number of complex zeros $`\beta _j=\beta _j(L)`$ grows with system size $`L`$ they will (for a system with no external field) pinch the positive real $`\beta `$-axis, and for large $`L`$ the corresponding value is the inverse of the physical critical temperature $`T_c`$. One example for the extension of these ideas to finite systems is the classification scheme by Borrmann et al.  which explores the linear behavior for the limiting density of zeros . Another description was suggested by Janke and Kenna (JK) who propose a new scaling relation to identify the order and strength of a transition from the behavior of small systems .

In this paper we try to evaluate the usefulness and validity of these two approaches. For this purpose we apply them to systems where a thermodynamic limit exists. In order to be useful the two approaches should be able to identify the macroscopic phase transition from investigations of small system sizes. Our test case is the 2D Potts model with $`q=5`$ states where the order of the phase transition is difficult to distinguish and which is therefore a challenging test for the new approaches. Simulations of the Potts model with $`q=4`$ states add data for the case of a second order phase transition. Our numerical investigations complement earlier work in Ref.  on the $`q=10`$ Potts model which has a pronounced first order phase transition.

Guided by our results for the Potts model we finally consider the helix-coil transition in polyalanine. The purpose of this investigation is to test whether the two approaches allow a characterization of structural transition in biological molecules.

## 2 PARTITION FUNCTION AND DENSITY OF COMPLEX ZEROS

In the canonical ensemble a system is completely described by its partition function

$$Z(\beta )=\underset{E}{}n(E)\mathrm{exp}(\beta E).$$
(1)
Introducing variables $`u=\mathrm{exp}(k\beta )`$ with conveniently defined constant $`k`$ allows to write the partition function for discrete energy models as a polynominal:

$$Z(u)=\underset{E}{}n(E)u^E.$$
(2)
The number of complex zeros $`u_j`$, ($`u_j=\mathrm{exp}(k\beta _j),j=1,2,\mathrm{}`$) of this polynominal will grow with system size. In the case of a phase transition, we expect the zeros (or at least the ones close to the real axis) condense for large enough system size $`L`$ on a single line

$$u_j=u_c+r_j\mathrm{exp}(i\phi ).$$
(3)
As the system size $`L`$ increases, those zeros will move towards the positive real $`u`$-axis and the corresponding value is for large $`L`$ the inverse of the physical critical temperature $`T_c`$. Crucial information on phase transitions can be obtained from the way in which the first zero approaches the real $`\beta `$-axis. However, such an analysis depends on the extrapolation towards the infinite large system and does not allow characterization of the situation in small systems.

One possible extension of the above ideas to “phase transitions” in small systems is the classification scheme by Borrmann et al. . The idea behind this scheme is to use not only the first complex zero $`u_1`$, but also higher zeros $`u_2,u_3`$ and $`u_4`$. Writing the complex zeros as $`u_k=\mathrm{Re}(u_k)+i\tau _k`$, where $`\tau _k`$ stands for $`\mathrm{Im}(u_k)`$, the assumed distribution of zeros on a straight line allows to define two parameters $`\alpha _u`$ and $`\gamma _u`$:

$$\alpha _u=\frac{\mathrm{ln}\varphi (\tau _3)\mathrm{ln}\varphi (\tau _2)}{\mathrm{ln}\tau _3\mathrm{ln}\tau _2},\mathrm{where}\varphi (\tau _k)=\frac{1}{2}\left(\frac{1}{|u_ku_{k1}|}+\frac{1}{|u_{k+1}u_k|}\right),$$
(4)
with $`k`$ labelling the first zeros, and

$$\gamma _u=(\mathrm{Re}(u_2)\mathrm{Re}(u_1))/(\tau _2\tau _1).$$
(5)
Note, that our notation differs from that in Ref.  in that we define the discrete line density $`\varphi `$ in function of the $`u`$-zeros instead of the temperature $`\beta `$. Following the classification scheme by Grossmann and Rosenhauer , phase transitions can now be classified according to the values of these two parameters: for $`\alpha _u0`$ and $`\gamma _u=0`$ one has a phase transition of first order, it is of second order if $`0<\alpha _u<1`$ and arbitrary $`\gamma _u`$, and for $`\alpha _u>1`$ and arbitrary $`\gamma _u`$ one has a higher order transition. Besides the above parameters, $`\tau _1`$ plays also an important role: only for $`\tau _10`$ one obtains in the thermodynamic limit a real temperature for a phase transition.

Another extension of partition function zero analysis to small systems is the approach by Janke and Kenna (JK)  which uses that the average cumulative density of zeros

$$G_L(r_j)=\frac{2j1}{2L^d}$$
(6)
can be written in the thermodynamic limit and for a first order transition as

$$G_{\mathrm{}}(r)=g_{\mathrm{}}(0)r+ar^{w+1}+\mathrm{}.$$
(7)
Here, the slope at the origin is related to the latent heat, $`\mathrm{\Delta }eg_{\mathrm{}}(0)`$. Equations (3) and (7) imply that the distance $`r_j`$ of a zero from its critical point can be written for large enough lattice sizes as $`\mathrm{Im}u_j(L)`$ since $`\mathrm{Re}u_j(L)u_c`$. Hence, in this limit Eq. (6) and (7) lead to the following scaling relation for the cumulative density of zeros as an equation in $`j`$ and $`L`$,

$$\frac{2j1}{2L^d}=a_1(\mathrm{Im}u_j(L))^{a_2}+a_3.$$
(8)
A necessary condition for the existence of a phase transition is that $`a_3`$ is compatible with zero, else it would indicate that the system is in a well-defined phase. The values of the constants $`a_1`$ and $`a_2`$ then characterize the phase transition. For instance, for first order transitions should the constant $`a_2`$ take values $`a_21`$ for small $`r`$, and in this case the slope of this equation is related to the latent heat through the relation

$$\mathrm{\Delta }e=ku_c\mathrm{\hspace{0.17em}2}\pi a_1,$$
(9)
with $`u_c=\mathrm{exp}(k\beta _c)`$. On the other hand, a value of $`a_2`$ larger than $`1`$ indicates a second order transition whose specific heat exponent is given by $`\alpha =2a_2`$.

The above approach was originally developed and tested for systems with well defined first order phase transitions such as the 2D 10-state Potts and 3D 3-state Potts model. The obtained results agree with previous work from numerical simulations and partition function zeros analysis of systems up to $`L=64`$  ($`L=36`$ for the 3D case ).

The classification of phase transitions by Borrmann et al. has been tested for the finite Bose-Einstein condensates in a harmonic trap and for small magnetic clusters and nuclear multifragmentation .

In the following we will show that despite these successful tests both approaches can lead to wrong conclusions if applied blindly. For this purpose we study systems where determination of the order of the phase transition is known to be computational difficult. We will concentrate on the 2D Potts model with $`q=4`$ (second order transition) and $`q=5`$ (weak first order transition) states since a series of exact results  exist for both models. For instance, the critical temperatures are known to be $`\beta _c=\mathrm{ln}(1+\sqrt{q})`$ and the latent heat is given by $`\mathrm{\Delta }e(q=5)=0.0529187`$. Finally, we research whether the two approaches can be used to characterize structural transitions in biological molecules.

## 3 RESULTS FROM POTTS MODEL SIMULATIONS

We start by presenting our results for the two Potts models with $`q=4`$ and $`q=5`$ on small lattices. Our study was performed on lattice sizes $`L=16,32,64`$ and $`128`$ for $`q=4`$ and $`L=8,12,16,20,24,32,64,96`$ and $`128`$ for $`q=5`$. For these models and system sizes we have evaluated the complex partition function zeros from heat-bath simulations with large statistics at temperatures $`\beta _0`$ listed in Table 1. For each lattice size, our results rely on 16 bins of 500000 measurements. For $`L<96`$ each measurement was separated by an additional Monte Carlo sweep which was discarded; for larger sizes measurements were only taken every 5th sweep.

For lattice sizes up to $`L=20`$, we could calculate all complex zeros from the polynominal form of the partition function using MATHEMATICA. For larger lattice sizes such a direct evaluation is no longer possible with standard numerical algorithms, even when double precision was used. Only by using the scan method (see  and references therein) we were able to obtain reliable estimates for the first $`J`$ complex zeros. Note, that $`J`$ is limited by the precision of our data: insufficient statistics of our simulation will lead to presence of so-called fake zeros . For instance, for the $`q=4`$ Potts model we could only get reliable estimates for the first three zeros (Table 2) while for the $`q=5`$ Potts model we were able to obtain estimates for the first four partition function zeros (presented in Table 3). The error bars were calculated from the fluctuation of the averages taken over each of those 16 bins.

### A The Borrmann Approach

From the listed zeros in Tables 2 and 3 we obtain our first result: even for our smallest lattice size ($`L=16`$ for $`q=4`$ and $`L=8`$ for $`q=5`$) the zeros stay on a straight line and the real parts are approximately constants. This is a necessary condition for both the Borrmann and the JK approach. We first calculate for both $`q=4`$ and $`q=5`$ Potts model the parameters $`\alpha _u`$ and $`\gamma _u`$ on which the Borrmann classification scheme relies. Their values are listed in Table 4 for all lattice sizes, and also plotted in Fig. 1a for the $`q=4`$ Potts model and in Fig. 1b for the case of $`q=5`$. For the $`q=4`$ Potts model we were only able to obtain reliable estimates for the first three zeros for each lattice size. However, in order to calculate the parameter $`\alpha _u`$ we need four zeros. For this reason, we have included here the less reliable fourth zero which is not listed in Table 2. As a consequence, our results for the $`q=4`$ Potts model are less reliable than these for the $`q=5`$ model where we have also acceptable estimates for the fourth zeros. Estimates and standard deviations presented in Table 4 were obtained by means of the bootstrap method  based on our statistics of 16 independent bins for each lattice size. Details on this method are presented in Appendix A.

For both models, the obtained values of Borrmann parameters $`\alpha _u`$ and $`\gamma _u`$ show no size dependence, but rather seems to fluctuate around some average value. For the case $`q=4`$, the values of $`\alpha _u`$ and $`\gamma _u`$ are compatible with the well-known fact that this model has a second order phase transition in the thermodynamic limit. However, our approach seems to fail for the case $`q=5`$ which has a weak first order phase transition and where we would expect $`\alpha _u0`$ and $`\gamma _u=0`$. Our results rather indicate a second order transition ($`0<\alpha _u<1`$ and $`\gamma _u`$ arbitrary). However, the $`q=5`$ Potts model is well-known to have a very large correlation length (of order 2000 lattice units ), and hence, its true behavior may only be caught for very large lattice sizes. Fig. 1b does not give any indications that we are even close to lattice sizes where the approach by Borrmann would lead to the correct result since the values for $`\alpha _u`$ and $`\gamma _u`$ show no systematic size dependence for the lattice sizes which we have studied ($`L128`$).

Our results indicate that an uncritical application of Borrmann approach may lead to wrong conclusions on the nature of phase transition in a system. In the case of the $`q=5`$ Potts model this approach fails to identify the nature of a phase transition from the distribution of zeros on small lattices. This seems to limit the usefulness of this method to systems where the order of the transition is clear.

### B The JK Approach

A similar method to study phase transitions in small systems is the one proposed by Janke and Kenna . In this approach one has to calculate the average cumulative density of zeros $`G(r)`$ from the zeros listed in Table 2 and 3. In order to investigate the phase transitions in our two systems one has to fit the cumulative density $`G(r)`$ to

$$G(r)=a_1r^{a_2}+a_3.$$
(10)
The aim of this fit is to obtain an estimate for $`a_3`$ since $`a_3=0`$ indicates the existence of a phase transition. However, a simple evaluation of this fit may be misleading if lattice sizes are too small. In this case it may be necessary to study the FSS dependence of these quantities, i.e. how their estimates are related to ones obtained for larger systems or even in the thermodynamic limit.

In the case of the $`q=4`$ Potts model we have to rely only on the first three zeros for each lattice size which is a too small number for a meaningful three-parameter fit. For this reason, we have combined the zeros of two neighboring lattice sizes and our fit therefore relies on six zeros for each pair. The so obtained estimates for $`a_3(L)`$ are listed in Table 5. For the case of the $`q=5`$ Potts model we have four zeros for each lattice size, but we have also calculated $`a_3(L)`$ from fits where the eight zeros of two neighboring lattice sizes were combined, which leads to a more robust estimate of this parameter. The so obtained sets of $`a_3(L)`$ are listed in Table 6. We note that for both models the values of $`a_3(L)`$ are for all lattice sizes compatible with zero, demonstrating that the two Potts models have indeed a phase transition.

The next question is whether and for what sizes the above approach is able to identify the order of the transition. This requires to calculate an accurate estimate for the quantity $`a_2`$ in Eq. (10). For this purpose, we set $`a_3`$ in Eq. (10) to zero and replace that equation by the simpler two-parameter fit

$$G(r)=a_1r^{a_2}.$$
(11)
Extracting the parameters $`a_1`$ and $`a_2`$ from this fit requires an extremely careful data analysis and error estimation. For this reason we have again used the bootstrap method for estimating averages and standard errors. The obtained values of $`a_1`$ and $`a_2`$ are presented for all lattice sizes and both $`q=4`$ and $`q=5`$ in Tables 5 and 6. In Table 6 we list these parameters as a function of the first four complex zeros.

Let us first consider the case $`q=4`$. Fig. 2 displays the values for $`a_2(L)`$ as a function of lattice size. In this plot we do not observe any dependency of $`a_2`$ from the system size and the possibility of a first order transition ($`a_2=1`$) is clearly excluded. Hence, our analysis reproduces the well-known fact that the $`q=4`$ Potts model has a second order phase transition. However, with our estimate of $`a_2`$ from our largest lattice size, we find as critical exponent $`\alpha =2a_2=0.430(18)`$. This value is far from the thermodynamic limit $`2/3`$ . This discrepancy may be due to the dependence of $`a_2(L)`$ on the higher and less precise zeros in the fits. This dependence is less pronounced if we merge the zeros of all lattice sizes. Evaluating the resulting cumulative density leads to a value of $`a_2=1.466(23)`$ which corresponds to $`\alpha =0.534(23)`$. This value of the critical exponent $`\alpha `$ is closer to, but still far from the theoretical value. It follows that on the simulated lattice sizes the JK scaling relations do allow a qualitative characterization of the transition in the $`q=4`$ Potts model, but not to obtain quantitative results such as the numerical value of critical exponents. This limitation is a somehow surprising result since the density and the cumulative density of zeros seem to be free  from the multiplicative logarithmic corrections to the leading power-law finite-size behavior  which hamper the determination of the critical exponents $`\alpha `$ and $`\nu `$ in other approaches.

At first sight, the situation seems to be worse for the $`q=5`$ Potts model. Even if we discard the fourth zero, which has larger fluctuations, and repeat our analysis only for the more reliable first three ones, our data do not lead to the expected result $`a_21`$. This can also be seen from Fig. 2 where we plot also for this model the parameter $`a_2(L)`$ as a function of lattice size $`L`$. However, unlike in the Borrmann analysis of the $`q=5`$ Potts model, our data show a clear trend towards the expected value for a first order phase transition ($`a_21`$). Our data suggest, that through by an extrapolation towards large (ideally infinite) lattices the true value of $`a_2`$ can be determined with the JK approach even for the $`q=5`$ Potts model. This is consistent with the fact that the $`q=5`$ Potts model has a very large correlation length (of order 2000 lattice units), and hence, its true behavior can only be caught for very large lattice sizes. The log-log plot of $`a_2(L)`$ in Fig. 3 indicates that one should use a polynominal fit for the extrapolation towards the thermodynamic limit:

$$a_2(L)=a_2(\mathrm{})+bL^c.$$
(12)
In order to see how acceptable the expected limit $`a_2=1`$ is, we replace the above equation by a two-parameter fit $`\mathrm{ln}(a_2(L)1)=\mathrm{ln}bc\mathrm{ln}L`$ (i.e. we set $`a_2=1`$). In this way, we obtain $`b=0.717(16)`$ and $`c=0.1765(69)`$, with the goodness of fit  $`Q=0.72`$ for four zeros and $`b=0.7910(88),c=0.2050(40)`$, with a smaller $`Q=0.09`$ for three zeros. The corresponding fit for the case of four zeros is shown in Fig. 4 and demonstrates that our data are indeed compatible with the expected value $`a_21`$ for a first order phase transition. However, we still can obtain acceptable fits ($`Q>0.7`$) for the whole range of $`0a_21.013`$ (but note that values of $`a_2<1`$ are difficult to interprete in the JK theory, and rather indicate numerical instabilities) which shows the difficulty in determining the order of the phase transition for the $`q=5`$ Potts model. Only if we restrict our analysis to the first three zeros the goodness of fit shows a maximum for $`a_21`$.

The JK approach allows us also to calculate the latent heat for the case of a first order phase transition by means of Eq. (9). However, that equation is valid only when $`a_2(L)=1`$. In the case of the $`q=5`$ Potts model we have $`a_2(L)>1`$ for all lattice sizes. Hence, we can not use Eq. (9) to calculate the latent heat from our values of $`a_1(L)`$. Instead we have replaced Eq. (11) by

$$G_L(r)=A_1(L)r,$$
(13)
i.e. $`a_2(L)`$ is substituted in Eq. (11) by $`a_2(\mathrm{})`$, and have studied the finite size scaling of the new quantity $`A_1(L)`$ which has the same limiting value than $`a_1(L)`$. Values for $`A_1`$ are listed in Table 6. In order to calculate the latent heat of the $`q=5`$ Potts-model we evaluate firstly $`a_1=A_1(\mathrm{})`$ from the finite-size-scaling fit

$$A_1(L)=A_1(\mathrm{})+BL^C,$$
(14)
which leads with a goodness-of-fit $`Q=0.68`$ to the values $`a_1=0.0284(18),B=1.381(21)`$ and $`C=0.6187(80)`$ when the first four zeros are used for each lattice size. Restricting the analysis to $`A_1(L)`$ calculated only from the first three zeros at each lattice size leads to $`a_1=0.03085(78)`$ (with a goodness of fit $`Q=0.75`$). Applying Eq. (9) we therefore find as latent heat $`\mathrm{\Delta }e=0.0551(35)`$ ($`\mathrm{\Delta }e=0.0599(15)`$ in case of three zeros) which is very close to the theoretical value $`\mathrm{\Delta }e(q=5)=0.0529\mathrm{}`$. This result is surprisingly good when compared with other numerical estimates. For instance, a recent study using lattices of up to $`L=4000`$ led to a latent heat value of $`\mathrm{\Delta }e=0.054`$ .

Note that Eq. (14) corresponds to the finite size scaling behavior of the specific heat if we identify $`C=(1\alpha )/\nu `$ since the latent heat scales with $`C_v\mathrm{\Delta }T_c(L)`$, and at the critical temperature we have $`\mathrm{\Delta }T_c(L)L^{1/\nu }`$ and

$$C_v^{\mathrm{max}}(L)=c_1+c_2L^{\alpha /\nu }.$$
(15)
This allows us to identify theoretically a second correction term in Eq. (14):

$$A_1(L)=A_1(\mathrm{})+B_1L^{(1\alpha )/\nu }+B_2L^{1/\nu }.$$
(16)
A check shows that our value of $`C=0.6187(80)`$ is indeed close to $`(1\alpha )/\nu =0.627`$ where the values of the called pseudocritical exponents $`\alpha =0.63(5)`$ and $`\nu =0.59(3)`$ were taken from finite estimates in Ref. .

We also see from Eq. (16) that for a strong first order phase transition (as for instance in the earlier studied $`q=10`$ Potts model ) where we have $`\nu =1/d`$ and $`\alpha =1`$, we do not have the first correction term in Eq. (16). This explains why even from small lattice sizes very good estimates of the latent heat ($`\mathrm{\Delta }e=0.698(2)`$ which one has to compare with the exact value for the $`q=10`$ Potts model: $`\mathrm{\Delta }e(q=10)=0.6960494`$ ) could be obtained for the $`q=10`$ Potts model in Ref. . In fact, these estimates can be easily improved by including a second correction term, which goes as $`L^2`$ for this model towards the new value $`\mathrm{\Delta }e=0.696(92)`$. The large error is due to the fact that only the largest sizes $`L=32,38,48\mathrm{and}\mathrm{\hspace{0.17em}64}`$, for which the first correction term can be discarded, were considered in our calculation. Here, our data relies on the values quoted in Ref. .

Our results indicate the JK scaling relations are more suitable than the Borrmann approach for studying phase transitions from the behavior of small systems. Unlike the Borrmann approach the order of the phase transition could be determined for both the $`q=4`$ and the $`q=5`$ Potts model. In the latter case (and for the $`q=10`$ Potts model which has also a first order phase transition) it was possible to calculate the latent heat with good accuracy. However, the approach failed to give the correct value for the critical exponent $`\alpha `$ in the case of the $`q=4`$ Potts model which has a second order transition.

## 4 Helix-coil transition in polyalanine

A common, ordered structure in proteins is the $`\alpha `$-helix and it is conjectured that formation of $`\alpha `$-helices is a key factor in the early stages of protein-folding . It is long known that $`\alpha `$-helices undergo a sharp transition towards a random coil state when the temperature is increased. The characteristics of this so-called helix-coil transition have been studied extensively . In previous work  evidence was presented that polyalanine exhibits a phase transition between the ordered helical state and the disordered random-coil state when interactions between all atoms in the molecule are taken into account. Here, we revisit polyalanine and investigate the helix-coil transition by means of the partition function analysis with the classification scheme of Borrmann and Janke/Kenna.

Our investigation of the helix-coil transition for polyalanine is based on a detailed, all-atom representation of that homopolymer. Since one can avoid the complications of electrostatic and hydrogen-bond interactions of side chains with the solvent for alanine (a nonpolar amino acid), explicit solvent molecules were neglected. The interaction between the atoms was described by a standard force field, ECEPP/2, (as implemented in the KONF90 program ). Chains of up to $`N=30`$ monomers were considered, and our results rely on multicanonical simulations  of $`N_{sw}`$ Monte Carlo sweeps starting from a random initial conformation, i.e. without introducing any bias. We chose $`N_{sw}`$=400,000, 500,000, 1,000,000, and 3,000,000 sweeps for $`N=10`$, 15, 20, and 30, respectively. Measurements were taken every fourth Monte Carlo sweep. Additional 40,000 sweeps ($`N=10`$) to 500,000 sweeps ($`N=30`$) were needed for the weight factor calculations by the iterative procedure described first in Refs. . In contrast to our first calculation of complex zeros presented in Ref. , where we divided the energy range into intervals of lengths 0.5 kcal/mol in order to make Eq. (2) a polynomial in the variable $`u=e^{\beta /2}`$, we avoided any approximation scheme in the present work. This because the above approximation works very well for the first zero, but not for the next ones. Since we need high precision estimates also for the next zeros we applied again the scan method.

In Table 7 we present our first four partition function zeros for 7 bins, although the fourth one is less reliable due to the presence of fake zeros. It is hardly possible to divide our production data in larger number of bins due to the limited statistics of our runs. Using again the bootstrap method we first calculated out of the zeros for each bin the parameters $`\alpha _u`$ and $`\gamma _u`$ which characterize in the Borrmann approach phase transitions in small systems. As one can see in Table 8 the so obtained values for polyalanine are characterized by large error bars and show no clear trend with chain length. It seems that the median of the $`\alpha _u`$ values is $`\alpha _u=0`$ which would indicate a first order transition. However, our data have too large errors to draw such a conclusion on the nature of the helix-coil transition.

For this reason, we tried instead the JK scaling relations. Table 9 lists the parameter $`a_3(N)`$ of Eq. (8). Here, the average cumulative density of zeros is replaced by

$$G_N(r_j)=\frac{2j1}{2N},$$
(17)
where we have translated the linear length $`L`$ as $`N^{1/d}`$ . Therefore all finite size scaling relations can be written is terms of the number of monomers $`N`$.

The values $`a_3`$ are compatible with zero for chains of all length indicating that we have indeed a phase transition. In order to evaluate the kind of transition we also calculate the parameters $`a_1(N)`$ and $`a_2(N)`$ which we also summarize in Table 9. Similar to the case of the $`q=5`$ Potts-model, the parameter $`a_2(N)`$ decreases with system size. The log-log plot of this quantity as a function of chain length in Fig. 5 suggests again a scaling relation

$$a_2(N)=a_2+bN^c.$$
(18)
A numerical fit of our data to this function leads to a value of $`a_2=1.31(4)`$ with $`Q=0.95`$. Using $`\alpha =2a_2`$ we find $`\alpha =0.70(4)`$ which is barely compatible with our previous value of $`\alpha =0.86(10)`$ in Ref. , obtained from the maximum of specific heat. A fit of all four chain lengths can also not exclude a value $`a_2=1`$ since we can find acceptable fits with $`Q>0.55`$ in the range $`0.92<a_2<1.44`$. However, a close examination of Fig. 5 shows that the $`N=30`$ data point shows a considerable deviation from the trend suggested by the smaller chain length. Since the $`N=30`$ data are the least reliable, we also evaluated Eq. (18) omitting the $`N=30`$ chain. This leads to a value of $`a_2=1.16(1)`$ and a critical exponent $`\alpha =0.84(1)`$ which is now compatible with our previous result $`\alpha =0.86(10)`$.

Let us summarize our results for the helix-coil transition studies of polyalanine. The JK approach is able to reproduce for polyalanine results obtained in previous work , but does not lead here to an improvement over other finite size scaling techniques. Especially, the JK approach does not allow to establish the order of the helix-coil transition from simulation of small chains. Our results for the parameter $`a_2`$ seem to favor a second order transition, but are due to large errors and disputed by Ref.  were indications for a finite latent heat were found. In the present study we considered only a special kind of biological molecules, homopolymers of amino acids, where in principle the thermodynamic limit can be considered. This allows finite size scaling which was an essential tool for obtaining the correct results with the JK approach. However, in general such finite size scaling is not possible for biomolecules that have a distinct size and composition. In these cases we have to rely on Borrmann’s approach. However, in our example of polyalanine, Borrmann’s approach led to even less conclusive results since the error bars were large. Hence, it seems that both approaches are limited in their usefulness for the study of phase-transitions in biomolecules.

## 5 Conclusion

We have evaluated two recently proposed schemes for characterizing phase transitions in small systems. Simulating the $`q=4`$ and $`q=5`$ Potts model, where for the thermodynamic limit we can compare our data with exact results, we found that both Borrmann and JK approach work well when the order of the phase transition is not in question (such as in the case of the $`q=4`$ Potts model (second order) or the earlier studied first order transition in the $`q=10`$ Potts model). The situation is different for cases such as the $`q=5`$ Potts model where it is difficult to distinguish between a weak first order or a strong second order transition. A careful application of the JK approach led to the correct result of a weak first order transition for the $`q=5`$ Potts model while the Borrmann approach did not allow a correct identification of this transition with the lattice sizes studied by us. Our results from Potts model simulations indicate that both approaches have to be applied with great care if one wants to avoid wrong conclusions on the nature of phase transition in a system. This may limit their usefulness to systems where the order of the transition is clear and to systems which are not too small. Especially, a first application of both approaches to helix-coil transitions in polyalanine demonstrates the difficulties appearing when they are applied to the study of phase transitions in biomolecules.

## 6 Appendix: Bootstrap method

Bootstrap is a simulation method based on a given data sample, to produce statistical inferences like standard error and bias-corrected estimator for the mean sample . The bootstrap algorithm consider that our data can be obtained from an unknown probability distribution $`F`$ by random sampling

$$F𝐗=(x_1,x_2,\mathrm{},x_n).$$
(19)
Here the points $`x_i`$ refer to our $`n=16`$ bins for each lattice size in the Potts model or to $`n=6`$ bins for polyalanine. For each bin, the points $`x_i`$ contain the first four complex zeros.

Our statistics of interest are estimates for mean values $`\widehat{\theta }=\theta (𝐗)`$ for the parameters $`\alpha `$ and $`\gamma `$ in the Borrmann approach

$`\widehat{\alpha }=(𝐗)`$ (20)
$`\widehat{\gamma }=𝒢(𝐗)`$ (21)
and their respectives standard errors $`\mathrm{se}(\widehat{\alpha })`$ and $`\mathrm{se}(\widehat{\gamma })`$. Here $``$ and $`𝒢`$ stand for the application of functions in Eqs. (4) and (5).

The bootstrap algorithm follows by considering our sample as an empirical distribution $`\widehat{F}`$, where each data has the probability $`1/n`$. The way bootstrap assigns an accuracy to our parameters does not depend on any theoretical calculation but random samples of size $`n`$ drawn with replacement from $`\widehat{F}`$, called a bootstrap sample,

$$\widehat{F}𝐗^{}=(x_1^{},x_2^{},\mathrm{},x_n^{}),$$
(22)
where each value $`x_i^{}`$ equals any one of the $`n`$ values $`x_i`$.

Therefore for each bootstrap sample $`𝐗^{}`$ there is a bootstrap replication of $`\widehat{\theta }`$ which we denote as $`\widehat{\theta }^{}`$. If we repeat this process, the sample standard deviation can be obtained from $`B`$ replications

$$\widehat{se}_B=\{\underset{b=1}{\overset{B}{}}[\widehat{\theta }^{}(b)\widehat{\theta }^{}()]^2/(B1)\}^{1/2},$$
(23)
where the bootstrap mean $`\widehat{\theta }^{}()=_{b=1}^B\widehat{\theta }^{}(b)/B`$. The limit of $`\widehat{se}_B`$ as $`B`$ goes to infinity gives the ideal bootstrap estimate of $`se(\widehat{\theta })`$.

The histogram in Fig. 6 represents our $`\widehat{\theta }^{}`$ distribution for the parameter $`\alpha `$ in our simulation of the $`q=4`$ Potts model for $`L=128`$ with $`B=400`$ replications. We obtain the bootstrap mean $`\widehat{\alpha }^{}=0.46(17)`$. On the other hand, if calculating simple averages from our 16 measurements leads to $`\widehat{\alpha }=0.41(18)`$.

The bootstrap estimate of bias based on $`B`$ replications is given by

$$\widehat{bias}_B=\widehat{\theta }^{}()\widehat{\theta }.$$
(24)
Therefore, the bias-corrected estimator is $`\overline{\theta }=\widehat{\theta }\widehat{bias}`$. It is also convenient to evaluate the ratio of estimated bias to standard deviation $`\widehat{bias}_{400}/\widehat{se}_{400}`$. A small number indicates an unneeded bias correction in face of the standard error. The called root mean square error of estimator $`\widehat{\theta }`$ for $`\theta `$ can be defined to take into account both bias and standard deviation ,

$$\sqrt{E_F[(\widehat{\theta }\theta )^2]}=\sqrt{se_F(\widehat{\theta })^2+bias_F(\widehat{\theta },\theta )^2}se_F(\widehat{\theta })[1+\frac{1}{2}(\frac{bias_F}{se_F})].$$
(25)
If we take $`\widehat{bias}=\widehat{bias}_{400}`$, our final result is $`\overline{\alpha }=0.36(18)`$ for $`q=4`$ and $`L=128`$. In the tables we quote our final estimates $`\overline{\theta }=\widehat{\theta }\widehat{bias}`$ and the above bias-corrected standard errors.

Acknowledgements
U.H.E. Hansmann gratefully acknowledges support by research grant of the National Science Foundation (CHE-9981874), and N.A. Alves support by CNPq (Brazil). Part of this article was written when U.H. was visiting the University of Central Florida. He thanks Brian Tonner, Alfons Schulte and the other faculty at the Department of Physics for the kind hospitality extended to him.

Figure Captions:

$`\alpha _u(L)`$ and $`\gamma _u(L)`$ estimates for $`q=4`$ Potts model data (a) and for $`q=5`$ Potts model data (b).

The parameter $`a_2(L)`$ as function of system size $`L`$ for the $`q=4`$ Potts model and the $`q=5`$ Potts model.

The parameter $`a_2(L)`$ as a function of system size $`L`$ for the $`q=5`$ Potts model in a log-log plot.

$`a_2(L)1(=a_2(\mathrm{}))`$ as a function of system size $`L`$ in a log-log plot. The straight line through the data points is from our fit.

The parameter $`a_2(N)`$ for polyalanine molecules of length $`N`$ in a log-log plot.

Histogram of 400 bootstrap replications for the parameter $`\alpha _u`$ from four-state Potts model data.

Table Captions:

Heat-bath simulations at $`\beta _0`$ for four and five-state Potts model in two dimensions.

Complex partition function zeros $`u_j,(j=1,2\mathrm{and}3)`$ for four-state Potts model.

Complex partition function zeros $`u_j,(j=1,2,3\mathrm{and}4)`$ for five-state Potts model.

Bootstrap bias-corrected estimates and bias-corrected standard error for the parameters $`\alpha _u`$ and $`\gamma _u`$ for $`q=4`$ and $`q=5`$.

Bootstrap bias-corrected estimates and bias-corrected standard error for the JK parameters from first three zeros for $`q=4`$.

Bootstrap bias-corrected estimates and bias-corrected standard error for the JK parameters from first four zeros for $`q=5`$.

Partition function zeros for polyalanine.

Bootstrap bias-corrected estimates and bias-corrected standard error for the parameters $`\alpha _u`$ and $`\gamma _u`$ for polyalanine model.

Bootstrap bias-corrected estimates and bias-corrected standard error for the JK parameters for polyalanine.


format: text

text:
# Ultra high energy cosmic rays: clustering, GUT scale and neutrino masses

## 1 Introduction

The interaction of protons with the microwave background predicts a drop in the cosmic ray flux above the GZK  cutoff $``$5$``$1019 eV. The data shows no such drop. Section 2 studies the propagation and determines the probability $`P(r,E,E_c)`$ that protons created at distance $`r`$ with energy $`E`$ reach earth above a threshold $`E_c`$. Using this $`P`$ one can give the observed spectrum by one numerical integration for any injection spectrum.

It is an interesting phenomenon that the UHECR events are clustered. Recently, a statistical analysis  based on the multiplicities of the clustered events estimated the source density. In Section 3 we extend the above analysis.

In Section 4 we study the scenario that the UHECRs are coming from decaying superheavy particles (SP) and we determine their masses $`m_X`$ by an analysis of the observed UHECR spectrum.

Ultrahigh energy neutrinos (UHE$`\nu `$) scatter on relic neutrinos (R$`\nu `$) producing Z bosons, which can decay hadronically (Z-burst) . In Section 5 we compare the predicted proton spectrum with observations and determine the mass of the heaviest R$`\nu `$ via a maximum likelihood analysis.

The details of the presented results and a more complete reference list can be found in .

## 2 Propagation of UHECR protons

Using pion production as the dominant effect of energy loss for protons at energies $`>`$1019 eV ref.  calculated $`P(r,E,E_c)`$ for three threshold energies. We extended the results of . In our Monte-Carlo approach protons are propagated in small steps ($`10`$ kpc), and after each step the probabilities of pion production and the energy losses due to pair production, pion production and the adiabatic expansion are calculated. We used the following type of parametrization $`P(r,E,E_c)=\mathrm{exp}\left[a(r/1\mathrm{Mpc})^b\right]`$. Fig. 1 shows $`a(E/E_c)`$ and $`b(E/E_c)`$ for a range of three orders of magnitude and for five different $`E_c`$.

## 3 Density of sources

The arrival directions of the UHECRs measured by experiments show some peculiar clustering: some events are grouped within $`3^o`$, the typical angular resolution of an experiment. Above $`410^{19}`$ eV 92 cosmic ray events were detected, including 7 doublets and 2 triplets. Above $`10^{20}`$ eV one doublet out of 14 events were found . The chance probability of such a clustering from uniform distribution is rather small .

The clustered features of the events initiated an interesting statistical analysis assuming compact UHECR sources . The authors found a large number, $`400`$ for the number of sources within the GZK sphere. We generalize their analysis. The most probable value for the source density is really large; however, the statistical significance of this result is rather weak.

Fig. 2 shows the CL regions for one of our models (injected energy distribution $`c(E)E^3`$; Schechter’s luminosity distribution: $`h(j)dj(j/j_{})^{1.25}\mathrm{exp}(j/j_{})d(j/j_{})`$). The regions are deformed, thin ellipse-like objects. For this model our final answer for the density is $`180_{165(174)}^{+2730(8817)}10^3`$ Mpc-3, where the first errors indicate the 68%, the second ones in the parenthesis the 95% CLs, respectively. The choice of  –h(j)$`\delta `$(j)– and, e.g. $`E^2`$ energy distribution gives much smaller value: $`2.77_{2.53(2.70)}^{+96.1(916)}10^3`$ Mpc-3, which is in a quite good agreement with their result.

## 4 Decay of GUT scale particles

An interesting idea discussed by refs. is that SPs could be the source of UHECRs. The hadronic decay of SPs yields protons. They are characterized by the fragmentation function (FF) $`D(x,Q^2)`$ which gives the number of produced protons with momentum fraction $`x`$ at energy scale $`Q`$. For the proton’s FF at present accelerator energies we use ref. . We evolve the FFs in ordinary and in supersymmetric QCD to the energies of the SPs. This result can be combined with the prediction of the MLLA technique, which gives the initial spectrum of UHECRs at the energy $`m_X`$ (cf. Fig.. 3). Similar results are obtained by .

Depending on the location of the source –halo or extragalactic (EG)– and the model –SM or MSSM– we study four different scenarios. In the EG case protons loose some fraction of their energies, described by $`P(r,E,E_c)`$. We compare the predicted and the observed spectrums by a maximum likelihood analysis. This analysis gives the mass of the SP and the error on it.

Fig. 4 shows the measured UHECR spectrum and the best fit, which is obtained in the EG-MSSM scenario. The goodnesses of the fits for the halo models are far worse. The SM and MSSM cases do not differ significantly. The most important message is that the masses of the best fits (EG cases) are compatible within the error bars with the MSSM gauge coupling unification GUT scale: $`m_X=10^b`$ GeV, where $`b=14.6_{1.7}^{+1.6}`$.

## 5 Z-burst scenario

Already in the early 80’s there were discussions that the UHE$`\nu `$ spectrum could have absorption dips at energies around $`E_{\nu _i}^{\mathrm{res}}`$=$`M_Z^2/(2m_{\nu _i})`$=$`4.210^{21}`$ (1 eV/$`m_{\nu _i}`$) eV due to resonant annihilation with R$`\nu `$s of mass $`m_\nu `$, predicted by the hot Big Bang, into Z bosons of mass $`M_Z`$ . Recently it was realized that the same annihilation mechanism gives a possible solution to the GZK problem . It was argued that the UHECRs above the GZK cutoff are from these Z-bursts.

We compare this scenario with observations.

The density distribution of R$`\nu `$s as hot DM follows the total mass distribution; however, it is less clustered. Thus we, as opposed to practically all previous authors , do not follow the unnatural assumption of having a relative overdensity of $`10^2÷10^4`$ in our neighborhood (for an approach with lepton asymmetry see ).

We give the energy distribution of the produced protons in our lab system, which is obtained by Lorentz transforming the CM collider results.

The next ingredient is the propagation of the protons, which can be described by $`P(r,E_p,E)`$.

Finally, we compare the predicted and observed spectrum and extract the mass of the R$`\nu `$ and the necessary UHE$`\nu `$ flux by a maximum likelihood analysis. In the Z-burst scenario small R$`\nu `$ mass needs large $`E_\nu ^{\mathrm{res}}`$ to produce a Z. Large $`E_\nu ^{\mathrm{res}}`$ results in a large Lorentz boost, thus large proton energy. In this way the detected energy determines the mass of the R$`\nu `$. The analysis is completely analogous to that of the previous section.

Our best fits to the observed data gives for the neutrino mass $`2.75_{0.97(1.89)}^{+1.28(3.15)}`$ eV for the “halo”- and $`0.26_{0.14(0.22)}^{+0.20(0.50)}`$ eV for the “EG”-case, respectively (in the halo/EG cases the power part of the UHECR spectrum is generated in the halo/EG). This gives an absolute lower bound on the mass of the heaviest $`\nu `$ of $`0.06`$ eV at the 95% CL.

The most attractive pattern for $`\nu `$ masses is hierarchical. Using the mass difference of the atmospheric $`\nu `$ oscillation for the heaviest mass , one obtains values between 0.03 and 0.09 eV. It is an intriguing feature of our result that the smaller one of the predicted masses is compatible on the $``$ 1.3 $`\sigma `$ level with this scenario.

The necessary UHE$`\nu `$ flux can be obtained from our fits. We summarized them in Fig. 5, together with some upper limits. We also compared our $`\gamma `$-flux with the EGRET bound . Our $`\gamma `$ flux is somewhat smaller than that of EGRET.

The nice collaboration with S.D. Katz and A. Ringwald and the careful reading of the manuscript is acknowledged. This work was partially supported by Hung. Sci. grants No. OTKA-T34980/T29803/T22929/M28413/OM-MU-708/IKTA111/NIIF.


text:
# Differential Elliptic Flow in 2 - 6 AGeV Au + Au Collisions: A New Constraint for the Nuclear Equation of State

## Abstract

Proton elliptic flow is studied as a function of impact-parameter $`b`$, for two transverse momentum cuts in 2 - 6 AGeV Au + Au collisions. The elliptic flow shows an essentially linear dependence on b (for $`1.5<b<8`$ fm) with a negative slope at 2 AGeV, a positive slope at 6 AGeV and a near zero slope at 4 AGeV. These dependencies serve as an important constraint for discriminating between various equations of state (EOS) for high density nuclear matter, and they provide important insights on the interplay between collision geometry and the expansion dynamics. Extensive comparisons of the measured and calculated differential flows provide further evidence for a softening of the EOS between 2 and 6 GeV/nucleon.

For several years now, the study of nuclear matter at high energy density has held the promise of providing valuable insights on the nuclear equation of state (EOS) and on the predicted phase transition to a quark-gluon plasma (QGP). At AGS energies of $`114`$ AGeV, elliptic flow has emerged as an invaluable probe of high density nuclear matter. This flow has been attributed to a delicate balance between (i) the ability of compressional pressure to effect a rapid transverse expansion of nuclear matter and (ii) the passage time for removal of the shadowing of participant hadrons by the projectile and target spectators. If the passage time is long compared to the expansion time, spectator nucleons serve to block the path of participant hadrons emitted toward the reaction plane, and nuclear matter is squeezed-out perpendicular to this plane giving rise to negative elliptic flow. For shorter passage times, the blocking of participant matter is significantly reduced and preferential in-plane emission or positive elliptic flow is favored because the geometry of the participant region exposes a larger surface area in the direction of the reaction plane. Thus, elliptic flow is predicted and found to be negative for beam energies $`4`$ AGeV and positive for higher beam energies.

Recent theoretical studies of elliptic flow have suggested a sensitivity to the pressure at maximum compression  and thus to the stiffness of the EOS, and to possible QGP formation. Despite this sensitivity, the commonly calculated patterns for elliptic flow very often do not constrain the EOS uniquely. This being the case, it is important to investigate additional experimental observables which may provide more stringent constraints for the EOS. Here, we investigate the utility of differential elliptic flow measurements $`v_2(b)`$ and $`v_2(b,p_T)`$, as possible constraints.

The measurements were performed at the Alternating Gradient Synchrotron (AGS) at the Brookhaven National Laboratory. Beams of 197Au at $`E_{Beam}=2`$, 4, and 6 AGeV were used to bombard a 197Au target of thickness calculated for a 3% interaction probability. Typical beam intensities resulted in $`10`$ spills/min with $`10^3`$ particles per spill. Charged reaction products were detected in the Time Projection Chamber (TPC) of the E895 experimental setup. The TPC located in the MPS magnet (typically at 1.0 Tesla) provided good acceptance and charge resolution for charged particles $`1<Z<6`$ at all three beam energies. However, a unique mass resolution for $`Z=1`$ particles was not achieved for all rigidities. Data were taken in two experimental runs with a trigger which allowed for a wide range of impact-parameter selections as presented below.

Our flow analysis follows the now standard procedure  of using the second Fourier coefficient, $`v_2=\mathrm{cos}2\varphi `$, to measure the elliptic flow or asymmetry of the proton azimuthal distributions at mid rapidity ($`|y_{cm}|<0.1`$). This distribution can be expanded as

$$\frac{dN}{d\varphi }\left[1+2v_1\mathrm{cos}(\varphi )+2v_2\mathrm{cos}(2\varphi )\right],$$
(1)
where $`\varphi `$ represents the azimuthal angle of an emitted proton relative to the reaction plane. Near mid rapidity in a symmetric system $`v_10`$. The reference azimuthal angle $`\mathrm{\Phi }_{plane}`$ of the reaction plane is determined using the vector Qi = $`_{ji}^nw(y_j)𝐩_{𝐓}^{}{}_{j}{}^{}/p_{T}^{}{}_{j}{}^{}`$. Here, $`𝐩_{𝐓}^{}{}_{j}{}^{}`$ and $`y_j`$ represent, respectively, the transverse momentum and the rapidity of baryon j (Z$`2`$) in an event. The weight $`w(y_j)`$ is assigned the value $`<p_x>/<p_T>`$, where $`p_x`$ is the transverse momentum in the reaction plane. The average $`<p_x>`$ is obtained from an earlier pass of an iterative procedure employed for each energy and impact-parameter selection.

The orientation of the impact-parameter vector follows azimuthal symmetry about the beam axis. Therefore, the azimuthal distribution of the determined reaction plane should be uniform or flat. We have established that deviations from this uniformity can be attributed to deficiencies in the acceptance of the TPC and have applied rapidity and multiplicity dependent corrections following Ref. . The corrections were applied for each of several impact-parameter selections at each beam energy; they ensure the absence of spurious elliptic flow signals which might result from distortions in the reaction plane distribution. The dispersion of the reaction plane $`<|\varphi _{12}|>/2`$ was estimated for each impact-parameter $`b`$ via the sub-event method. Suffice to say, a reasonable resolution was observed over the entire range of energy and $`b`$ studied. These estimates for the reaction plane dispersion serve as a basis for evaluating the dispersion corrections summarized in Table 1; these corrections have been applied to the extracted flow values discussed below.

The event multiplicity of identified charged particles $`M_{filt}`$ was used for centrality selection. That is, several multiplicity bins were selected in the range from 0.4 to 1.0 $`M_{max}`$ where $`M_{max}`$ is the point in the charged particle multiplicity distribution where the height of the distribution has fallen to half its plateau value. Impact-parameter estimates have also been made for these centrality selections, at each beam energy, via their respective fraction of the minimum bias cross section.

Figure 1 shows representative distributions in the azimuthal angle $`\varphi `$ obtained at the energies of 2, 4 and 6 AGeV for mid-rapidity ($`|y_{cm}^{(o)}|<0.1`$) protons. The panels from left to right represent the three beam energies, respectively, and from top to bottom the three impact-parameter ranges of $`0b3`$, $`4b6`$ and $`7b8`$ fm. For visual clarity, a $`p_T`$ cut has been applied to the distributions shown for both the 4 and 6 AGeV data, as indicated. Within each $`b`$-range in Fig. 1, the previously reported transition from negative to positive elliptic flow at $`4`$ AGeV is clearly seen. That is, the elliptic flow is negative at 2 AGeV, positive at 6 AGeV and essentially zero at 4 AGeV. An apparent increase of the anisotropy of the distributions with increasing $`b`$ can also be discerned for the 2 and 6 AGeV data shown in Fig. 1. We attribute this trend to an interplay of the changing geometry with the expansion of excited participant matter as discussed below.

Figure 2 shows the $`v_2`$ coefficients for the full $`p_T`$ range, as a function of $`b`$ for data (filled stars) obtained at 2, 4, and 6 AGeV in the three panels, respectively. These coefficients have been obtained by evaluating the $`\mathrm{cos}2\varphi `$ for each azimuthal distribution obtained for a given impact-parameter at each beam energy. A correction has been applied to some of these coefficients to account for biases resulting from (i) low $`p_T`$ acceptance losses in the TPC for the 2, 4, and 6 AGeV beams, (ii) high $`p_T`$ acceptance losses in the TPC for the 2 AGeV beam, and (iii) $`\pi ^+`$ contamination of the proton sample at 4 and 6 AGeV. A procedure for effecting these corrections has been detailed in Ref.. That is, we first plotted the observed Fourier coefficient $`\mathrm{cos}2\varphi ^{}`$ vs. $`p_T`$ with $`p_T`$ thresholds which allowed clean particle separation ($`p_T1`$ GeV/c). We then extracted the coefficients for the quadratic dependence of $`\mathrm{cos}2\varphi ^{}`$ on $`p_T`$. These quadratic fits are restricted by the requirement that $`\mathrm{cos}2\varphi ^{}=0`$ for $`p_T=0`$. Next, we corrected the proton $`p_T`$ distributions for possible high and low $`p_T`$ losses. A weighted average (relative number of protons in a $`p_T`$ bin times the $`\mathrm{cos}2\varphi ^{}`$ for that bin) was then performed to obtain $`\mathrm{cos}2\varphi ^{}`$ for each beam energy. The corrections which result from this procedure are $``$ 5% for the 4 and 6 AGeV beams and $``$ 15% for the 2 AGeV beam. Subsequent to these evaluations, the $`v_2`$ values were corrected for reaction plane dispersion using the procedures detailed in Refs. .

The $`v_2`$ values represented by filled stars in Fig. 2 indicate an essentially linear dependence on impact parameter. The slope of this dependence is clearly negative and positive for the 2 and 6 AGeV data, respectively. By contrast, an essentially flat dependence is observed for the 4 AGeV data suggesting that the beam energy at which the elliptic flow changes sign is not very sensitive to $`b`$ for $`0b8`$fm. The approximately linear dependence exhibited by the data can be understood in terms of the collision geometry and the development of transverse expansion within the participant matter. At 2 AGeV an expansion perpendicular to the reaction plane develops over the characteristic time of $`d/c_s`$ while the spectators are present. Here, $`c_s=\sqrt{p/e}`$ represents the speed of sound for a given pressure $`p`$, and energy density $`e`$, and $`d`$ is the perpendicular distance from the center of the participant region to the surface. The spectator passage time (estimated in sharp cut-off geometry) first increases and then remains essentially constant as $`b`$ increases over the range of interest. On the other hand, the expansion time decreases with increasing $`b`$ due to a decrease in $`d`$. It is this decrease in the expansion time coupled with an essentially constant passage time, which provides the driving force for more matter to escape the interaction region as $`b`$ is increased i.e. an increase in “squeeze-out” with $`b`$. The magnitude of the “squeeze-out” follows an approximately linear dependence because $`d`$ is roughly proportional to $`1/b`$ for the Au + Au impact parameter range 1 - 8 fm.

At 6 AGeV the spectator passage time is very short compared to the expansion time and preferential in-plane emission dominates. In this case, however, the linear increase of $`v_2`$ with increasing impact parameter is strongly influenced by the initial spatial asymmetry of the nuclei overlap region or participant matter. This asymmetry is commonly characterized in terms of the width $`L_x`$ and height $`L_y`$ of the overlapping region via $`\alpha _s=(L_yL_x)/(L_y+L_x)`$ and can be shown to be nearly linearly proportional to the impact parameter for medium b values. The essentially flat dependence of $`v_2`$ observed at 4 AGeV suggests that, at the transition energy, the reduction in the expansion time \[in competition with the spectator passage time\] with increasing $`b`$, is compensated for by the (later) increased in-plane-emission from the preserved initial spatial asymmetry.

The solid circles, open squares and solid triangles shown in Figs. 2, represent results from calculations with a recent version of the Boltzmann Equation Model BEM which assumes a soft ($`K=210`$ MeV), a stiff ($`K=380`$ MeV) and an intermediate ($`K=300`$ MeV) EOS respectively. The calculations include momentum dependent forces. A comparison of the calculated $`v_2`$ values indicate sizeable differences between the predictions for a stiff and a soft EOS for all three beam energies. For both 2 and 4 AGeV this distinction increases with increasing impact parameter indicating that the impact parameter dependence of elliptic flow lends a new and important constraint for the EOS. At 2 AGeV, the $`v_2`$ values for the stiff EOS show good agreement, both in magnitude and trend, with the experimental data. At 4 AGeV the measured $`v_2`$ values lie between the calculated result for a stiff and a soft EOS, and appear to be in better overall agreement with an intermediate form of the EOS. At 6 AGeV the data is less compatible with a stiff EOS, but does not allow a clear distinction between the soft and the intermediate ($`K=300`$ MeV) EOS. The results of these comparisons can be taken as being suggestive of a softening of the EOS as previously reported in Ref.. However, it is interesting to investigate whether or not the differential flow measurement $`v_2(b,p_T)`$ provides further constraints. This line is pursued below.

Figs. 3 compares experimental (stars) and calculated (circles, triangles and squares) differential elliptic flow $`v_2(b,p_T)`$ for 2, 4 and 6 AGeV as indicated. At each beam energy, the BEM calculations have been carried out for the same $`p_T`$ and $`b`$ selections applied to the data. Fig. 3 indicates good agreement between the data and the calculated results for a stiff EOS at 2 AGeV. At 4 AGeV the data again shows better overall agreement with the intermediate and soft EOS. At 6 AGeV the comparison indicates quite good agreement (both in magnitude and trend) between the data and the results from the calculations which assume a soft EOS. The latter agreement is in contrast to the results obtained from the comparison made in Fig. 2, and clearly indicates that the differential flow $`v_2(b,p_T)`$, does indeed provide an additional constraint for making a relatively clear distinction between the different EOS’s at 6 AGeV.

To summarize, we have studied differential proton elliptic flow in 2 - 6 AGeV Au + Au collisions. The elliptic flow shows an essentially linear dependence on $`b`$, in the range $`1.5b8`$ fm, with a negative slope at 2 AGeV, approximately zero slope at 4 AGeV, and a positive slope at 6 AGeV. This dependence can be understood in terms of (a) the relationship between the collision geometry, (b) the relative magnitude of the time for development of the transverse expansion, and (c) the passage time for removal of the shadowing of participant hadrons by the projectile and target spectators. Detailed comparisons between the measured differential elliptic flow $`v_2(b,p_T)`$, and $`v_2(b)`$, and the results obtained from a relativistic Boltzmann-equation calculation not only suggest a softening of the EOS but also indicate that differential flow measurements provide very important constraints for the determination of the EOS of high density nuclear matter.

###### Acknowledgements.

This work was supported in part by the U.S. Department of Energy under grants DE-FG02-87ER40331.A008, DE-FG02-89ER40531, DE-FG02-88ER40408, DE-FG02-87ER40324, and contract DE-AC03-76SF00098; by the US National Science Foundation under Grants No. PHY-98-04672, PHY-9722653, PHY-0070818, PHY-9601271, and PHY-9225096; and by the University of Auckland Research Committee, NZ/USA Cooperative Science Programme CSP 95/33.

|  |  | Dispersion Correction Factor |  |
| --- | --- | --- | --- |
| b range (fm) | 2 AGeV | 4 AGeV | 6 AGeV |
| $`0<b<3`$ | 1.71 | 2.64 | 4.65 |
| $`4<b<6`$ | 1.22 | 1.59 | 2.47 |
| $`7<b<8`$ | 1.26 | 1.99 | 2.86 |

Table 1: Correction factors for reaction plane dispersion for several impact parameter ranges for the 2, 4 and 6 AGeV beam energies.


text:
# Periodic Oscillations of Josephson-Vortex Flow Resistance in Bi2Sr2CaCu2O8+y

## Abstract

To study the Josephson-vortex system in the intrinsic Josephson junctions, we have measured the flow resistance as a function of magnetic field parallel to the $`ab`$ plane in Bi2Sr2CaCu2O8+y single crystals. Although it was known that the flow resistance increases smoothly with increasing field, we have found novel oscillations of vortex flow resistance in the small current regime. The period of the oscillations corresponds to the field which is needed to add “one” vortex quantum per “two” Josephson junctions. Commensurability between the lattice spacing of Josephson vortices along the $`ab`$ plane and the size of the junction is related to the flow velocity of vortices. The results show that Josephson vortices form triangular lattice in the state where the oscillations occur.

When an external magnetic field is applied along the plane of Josephson junctions in the superconducting state, the field penetrates into the junctions as quantized Josephson vortices. Since high-$`T_\mathrm{c}`$ cuprates such as Bi2Sr2CaCu2O8+y (BSCCO) naturally consist of atomic-scale Josephson junctions (intrinsic Josephson junctions; IJJ’s) based on the layered crystal structures and large anisotropy , the magnetic field penetrated along superconducting planes forms a Josephson-vortex system. Recently, the Josephson-vortex system in high-$`T_c`$ cuprates has attracted much attention because of fundamental interests and possibilities for applications.

When the intrinsic pinning by the superconducting layers  becomes strong enough, the commensuration of the vortex lattice with the underlying layered structure becomes important, so that the drastic changes in the vortex phase diagram is expected. To clarify thermodynamic phases of the Josephson-vortex system, both theoretical and experimental studies have been done extensively. Theoretically, several new phases , a smectic phase , and a tricritical point where the order of the melting transition changes  have been predicted. Experimentally, the important facts such as the melting transition of the Josephson-vortex systems , the oscillatory melting temperature , and the existence of the vortex smectic phase  have been found in YBa2Cu3O7-δ (YBCO).

Interesting phenomena appear when the vortex dynamics is influenced by a periodic one-dimensional pinning potential. Periodic oscillations of pinning strength as a function of the magnetic field $`H`$ have been observed in an artificially structured superconducting film with channels of weak pinning in a strong-pinning environment . The commensurability between vortex lattice and the channel width is related to this phenomenon. Additionally, Gurevich et al. have reported the periodic oscillations of in-plane resistance $`R(H)`$ in superconducting multilayers of NbTi/Cu when the Josephson vortices move perpendicular to the plane . They concluded that the origin is a dynamic matching effect in terms of the generation of standing electromagnetic waves in moving vortex lattice along the $`c`$-axis.

Although the influence of the periodic pinning potential is expected to be more important in highly anisotropic materials such as BSCCO than in YBCO and others, there is little information on the Josephson vortices in BSCCO. In this Letter, we report on novel periodic oscillations of the flow resistance of Josephson vortices as a function of the in-plane magnetic field. The period of the oscillations suggests that the ground state of Josephson vortices is in triangular lattice when the oscillations occur.

High quality single crystals of BSCCO were grown by traveling-solvent floating-zone technique . A platelet of single crystals was carefully cut into narrow strips. After four contacts with electrodes were formed by using a silver paste, the center of strips of the single crystals was cut by focused ion beam (FIB). The same process to fabricate IJJ’s has been used in the case of BSCCO whiskers . An FIB image of a typical sample is shown in Fig. 1. The inset is the schematic illustration of the junction. The dimensions (width ($`w`$)$`\times `$length ($`l`$)$`\times `$thickness ($`t`$)) of the measured samples A, B and C are 7.3$`\times `$13$`\times `$1.3, 18$`\times `$16$`\times `$1.6 and 31$`\times `$17$`\times `$2.0 $`\mu `$m3, respectively. The superconducting transition temperature $`T_\mathrm{c}`$ is 86 K in all samples.

While resistance of strips was smaller than 1$`\mathrm{\Omega }`$ before the fabrication of junctions, the resistance of sample A, B and C increased to 1400, 560 and 315 $`\mathrm{\Omega }`$ at room temperature after the fabrications, respectively. Therefore, the measured resistance is almost $`c`$-axis resistance of the junction parts and the contribution of in-plane resistivity can be ignored. In the small mesa crystals of BSCCO, the flow of Josephson vortices, which is not affected very much by pancake vortices, has been realized in the so-called lock-in state . The $`ab`$ plane of the sample was possible to orient parallel to the magnetic field to an accuracy of 0.005 degrees by adjusting the flux-flow resistance as a function of the misorientation angle to its maximum.

When the field is applied parallel to the $`ab`$ planes in the superconducting state, the current along the $`c`$-axis drives Josephson vortices to the direction perpendicular to both current and field. The motion of the vortices causes the flow resistance. Up to the present it has been reported that the flow resistance increases smoothly with increasing field because the number of vortices increases . However, we have discovered novel oscillations as a function of field as shown in Fig. 2. When the magnetic field is increased, finite flow resistance appears in 2 $``$ 6 kOe depending on the sample and the current density. This field may be related to the pinning effect of Josephson vortices due to the quality of the surface boundary of samples. The oscillations start from a constant field $`H_{\mathrm{start}}`$ 7-8kOe in all samples. The period $`H_\mathrm{p}`$ of the oscillations is quite constant in the wide range of magnetic fields. As the magnetic field is increased, the flow resistance increases accompanying with the oscillations and suddenly drops above a field $`H_{\mathrm{stop}}`$.

$`H_{\mathrm{stop}}`$ decreases with increasing the absolute value of the angle ($`\theta `$) between the magnetic field and the $`ab`$ plane. This may be due to the breakdown of the lock-in state because the pancake vortices penetrate into the junction and are pinned to stop the flow of Josephson vortices when the $`c`$-axis component of the magnetic field exceeds the $`c`$-axis lower critical field $`H_{\mathrm{c1}}`$. Actually, $`H_{\mathrm{stop}}`$ is roughly proportional to $`1/\mathrm{sin}\theta `$ and the coefficient is of the order of $`H_{\mathrm{c1}}`$. The maximum of $`H_{\mathrm{stop}}`$ could be reached to 40 kOe with our alignment set-up. On the other hand, $`H_{\mathrm{start}}`$ and $`H_\mathrm{p}`$ does not almost change even if the field is intentionally shifted from the optimal position.

The temperature dependence of flow resistance has also been measured from 4.2 K to $`T_\mathrm{c}`$. The periodic oscillations can be observed at all temperatures below $`(T_\mathrm{c}5)`$ K. $`H_\mathrm{p}`$ is independent of temperatures. These facts show that the oscillations relate only to the field component parallel to the $`ab`$ plane, namely, the dynamics of the Josephson vortices.

It is interesting that periodic oscillations can be observed only in sufficiently small current regime. $`I`$-$`V`$ characteristics of the sample B are shown in Fig. 3. In the fields where the resistance is minimum, the nonlinearity of $`I`$-$`V`$ curves becomes large and a kink structure is observed around 50$`\mu `$A. The periodic appearance of this structure as a function of fields corresponds to the oscillations of flow resistance. In the larger currents $`>`$ 200$`\mu `$A the oscillatory behavior disappears thoroughly.

$`H_\mathrm{p}`$ depends on the dimensions of the junctions. To extract $`H_\mathrm{p}`$, the power spectral densities for three samples which are obtained from the data in Fig. 2 by FFT are demonstrated in Fig. 4. The sharp fundamental peaks indicate that $`H_\mathrm{p}`$ is uniform over the wide fields. The small peaks corresponding to the second harmonic waves are observed because the oscillations are distorted slightly from sinusoidal. The positions of the peaks, which are inversely proportional to $`H_\mathrm{p}`$, are plotted as a function of the width ($`w`$) of the junctions in the inset of Fig. 4. We see that the $`1/H_\mathrm{p}`$ is proportional to the width ($`w`$).

What is the origin of these oscillations? It should be checked on the relation between $`H_\mathrm{p}`$ and the number of vortices in the junctions. Supposing that the Josephson vortices uniformly penetrate into every junction, the increment of the field ($`H_0`$) which is needed to add one vortex per a junction is $`\varphi _0/ws`$, where $`\varphi _0`$ ($`2.07\times 10^7`$ G cm2) is flux quantum and $`s`$ (=15Å) the junction thickness, namely, the distance between superconducting CuO2 layers. $`1/H_0(w)`$ is shown as a broken line in the inset of Fig.4. Interestingly, the experimental values of $`1/H_\mathrm{p}`$ are close to the twice of $`1/H_0(w)`$ as shown by a solid line, which means that $`H_\mathrm{p}`$ corresponds to adding “one” flux quantum per “two” IJJ’s. We can check whether the resistance takes a maximum or a minimum when the field is $`nH_\mathrm{p}`$ ($`n`$ : integer). Experimentally, we have shown that the oscillations have minima at $`nH_\mathrm{p}`$.

The above experimental results can be explained by assuming that; (1) the lattice structure of the Josephson-vortex system is the triangular lattice, and (2) at both right and left surfaces of the junction, a surface barrier, i.e., Bean-Livingston barrier  works on vortices, although both attractive image force and repulsive force by Meissner current are quite small because of a little inhomogeneity of the magnetic field around vortices in junctions. For example, let’s consider that the field is $`n(H_0/2)nH_\mathrm{p}`$. The Josephson vortices are distributed as schematically shown in Fig.5(a). This ordered state fits to the width of junctions. When the vortex lattice flow from left to right, vortices at right side exit beyond small surface potentials and vortices would enter at left side simultaneously. Now, the total potential which the lattice feels at both sides becomes maximum. At this situation, the averaged velocity and hence resistance of vortex flow has minimum. In a short time scale, the velocity of the vortex lattice would oscillate with time whenever vortices go over the boundaries. As the field is increased, additional vortices are forced to enter into the ordered lattice. The Josephson vortex lattice could move more easily because of the incommensurability between the width of the junctions and the lattice spacing of vortices. Then, the flow resistance becomes larger than that in the above fitted case. When the field goes over $`(n+1/2)(H_0/2)`$, the lattice starts to return to the fitted state. At $`(n+1)(H_0/2)`$, matching between the Josephson vortex lattice and the size of junction recovers as shown in Fig.5(b), where the Josephson vortices enter in alternate layers comparing with Fig.5(a). Hence, it can be plausible that the oscillation period of the flow resistance is a half of $`H_0`$. When the Lorentz force from the applied current becomes large enough comparing with the surface potential, it is expected that the oscillations smear out in this scenario, which can explain the $`I`$-$`V`$ characteristics qualitatively. More complicated models are also possible, which may explain the details of flow dynamics.

There are several theoretical studies with respect to the thermodynamic phase diagram of the Josephson vortex system in anisotropic layered superconductors . Ivlev et al.  have derived the structure of Josephson-vortex lattice as the compressed hexagons of triangular lattice pointing along the $`c`$-axis in high enough fields by using London theory. However, in the case of strong magnetic fields of order $`\varphi _0/\gamma s^2`$, where $`\gamma `$ is anisotropy ratio, the three-dimensional(3D) anisotropic London or Ginzburg-Landau approach is not adequate . This field is estimated as 46 kOe in $`\gamma 200`$ . It is delicate whether the London theory can be applied at the field where the oscillations occur. On the other hand, the extensive Monte Carlo simulations with the 3D anisotropic frustrated XY model by Hu et al. have also shown that a triangular lattice is the ground state where the Josephson vortices distribute in every layer for the anisotropy parameter $`\gamma `$8 when the filling factor is 1/32 . For $`\gamma =`$150 as is the case of BSCCO, it is known that the compressed triangular lattice is achieved even below the tricritical point value ($``$ 15 kOe) by rescaling . These theoretical results are consistent with our experimental ones. For the magnetic fields lower than 7kOe, the oscillations cannot be observed in all samples in spite of the existence of finite flow resistance. A reason why the oscillations smear out at the low magnetic fields may be that the distance between the vortices is not close enough to form a triangular lattice. Ivlev et al. have predicted that a triangular array becomes unstable with respect to the formation of several new types of vortex lattices in the low fields . In addition to the ground state, the existence of several metastable states has been indicated . In this situation, the periodic oscillations could not be expected.

Similar phenomena have been reported in artificially structured materials as described above , although the origins are different from our case. Recently, in YBCO single crystals of 60 K phase which is more anisotropic than 90 K phase, an oscillatory behavior of the in-plane resistivity in the field parallel to the $`ab`$ plane have been reported . In their experiments, the oscillations have minima when the period of the layered structure of YBCO matches the intervortex spacing in the direction normal to the layers. The period of the oscillation is proportional to $`B^{0.5}`$ and is independent of the size of samples. This phenomenon originates from the commensurability between $`s`$ and vortex-lattice space along the $`c`$-axis, while the two lengths, the junction width $`w`$ and vortex-lattice spacing along the $`ab`$ plane are important in present case.

In summary, we have measured the flow resistance of Josephson vortices in IJJ’s of BSCCO single crystals. The novel oscillations of the resistance as a function of field by the Josephson-vortex flow have been found. The period of the oscillations depends only on the width of junctions. The fact that the period is $`H_0/2`$, which is needed to add “one” Josephson vortex per“two” layers, suggests that the lattice structure of the Josephson vortices is triangular. The oscillations originate from commensuration between the triangular lattice and the width of the junction perpendicular to the field. At present, our measurements do not access the boundary between the lattice and liquid phase of Josephson-vortex system. In the future, further improvements of measurements would make it possible to investigate the phase diagram of Josephson-vortex system including a liquid state by using this phenomenon.

The authors thank X. Hu and M. Machida for useful discussions, N. Ishikawa for technical supports of FIB.


text:
# Critical Dynamics of Burst Instabilities in the Portevin-Le Châtelier Effect

\[

## Abstract

We investigate the Portevin-Le Châtelier effect (PLC), by compressing Al-Mg alloys in a very large deformation range, and interpret the results from the viewpoint of phase transitions and critical phenomena. The system undergoes two dynamical phase transitions between intermittent (or “jerky”) and “laminar” plastic dynamic phases. Near these two dynamic critical points, the order parameter $`1/\tau `$ of the PLC effect exhibits large fluctuations, and “critical slowing down” (i.e., the number $`\tau `$ of bursts, or plastic instabilities, per unit time slows down considerably).

\]

The motion of lattice dislocations is the main source of macroscopic plastic deformation in crystalline solids. In principle, the motion of a single dislocation in an otherwise perfect crystal is a simple phenomenon. However, in a real crystal, many cooperative and core dislocation effects may give rise to unusual and complex plastic properties. An intriguing example of the latter is provided by the Portevin-Le Châtelier (PLC) effect , as observed in physical metallurgy (see, e.g., ). In stress-strain experiments in which the deformation rate is held constant, the stress usually increases monotonically with the strain, and the solid deforms homogeneously in a continuously flowing or “laminar” plastic regime. When the PLC effect arises, notably in some alloys and at high deformation, intermittent yielding points, or plastic instabilities are observed, and the deformation is spatially inhomogeneous and apparently confined to mesoscopic channels in the sample. This intermittently flowing or “jerky” plastic regime has been analyzed in terms of microscopic ageing models, which ascribe it to dislocation-avalanches triggered by the interaction of dislocations and clouds of mobile impurities, or to closely related mechanisms (see, e.g., ). In these models, mobile point defects are attracted by the stress field of an immobile dislocation, condensing in a cloud around the dislocation. Thus, the pinning force necessary to move the dislocation increases in time and the system ages. By increasing the applied force, and as soon as the stress becomes large enough, the pinning breaks down and the dislocation flashes over a large distance to stop again. This distance is determined by the arrangement of the other dislocations sliding along the same channel and by strong pinning defects. This cooperative process repeats itself again and again, producing a succession of plastic instabilities.

Beyond its importance in metallurgy, the PLC effect is a paradigm for a general class of nonlinear complex systems with intermittent bursts. Indeed, the succession of plastic instabilities shares both physical and statistical properties with many other systems exhibiting loading-unloading cycles. This is because the macroscopic statistical properties are a consequence of similar underlying microscopic physical mechanisms (e.g., log-interacting dislocation lines or magnetic flux lines, slowly driven in a disordered landscape with pinning traps, and obeying glassy dissipative dynamics). Weertman, and more recently Lebyodkin et al., have analyzed the analogy of the PLC effect with earthquakes. But similarities also exist to stick-slip phenomena in sheared granular media, avalanches of magnetic vortices in superconductors, “starquakes” with $`\gamma `$–ray burst activity in magnetars, and even bursts of economic activity and stock market crashes. The dynamic behavior of some of these systems can be viewed in terms of dynamic phases separated by dynamical critical points (see, e.g.,  and refs. therein). This is achieved by adapting concepts of equilibrium phase transitions  to non-equilibrium cases. It is tempting to extend the analogy to the PLC effect by further studying the dynamics of the jerky regime. Here we present results of compression experiments on $`34`$ polycrystalline Al-Mg samples with nominal 4%Mg. In contrast to the more common traction technique, compression experiments allow to access very large deformations without sample failure, and to collect very long time-series. We interpret our results from the viewpoint of phase transitions and critical phenomena.

Experiments were conducted on $`34`$ samples of the alloy Al-4at.%Mg; a classical system exhibiting the PLC effect. Data were collected at room temperature, about 18$`{}_{}{}^{}C`$. Samples were slightly oblique cylinders of about 6 mm length and 4 mm diameter, inclined 2 degrees from the vertical; Grain size in the virgin samples was about 100 $`\mu `$m. In Fig. 1 we show a typical segment of a stress versus time curve, for an imposed compression speed $`v=0.18\mu `$m/s, and when the sample is in the range of a 17% deformation. (The compression velocity is $`v=dl/dt`$, where $`l`$ is the actual length of the sample. Deformation is calculated assuming constant volume by $`ϵ=\mathrm{ln}\left(1\mathrm{\Delta }l/l_o\right)`$, where $`\mathrm{\Delta }l`$ is the elongation and $`l_o`$ the initial length of the sample ). One sees various sequences of elastic reloading of the sample until a yield point is reached, followed by the sharp stress drop typical of the PLC effect. In Fig. 1 we show the negative derivative of the stress versus time, $`d\sigma /dt`$. The latter evidences the stress drops as sharp “bursts” of energy and permits us to clearly define the waiting time $`\tau _n`$ between successive events, and the size $`s_n`$ of the nth burst. In this way, we also eliminate difficulties related to the small drift of the average stress with the strain hardening, which results in an almost constant offset in $`d\sigma /dt`$.

FIG. 1. Short time-sequence, taken in the range of $`17\%`$ deformation, of the stress $`\sigma `$ (and also $`d\sigma /dt`$) versus time for an Al-4at.%Mg polycrystalline sample. For all figures, the imposed compression speed $`v`$ is $`0.18\mu `$m/s, unless otherwise indicated.

FIG. 2. Various $`d\sigma /dt`$ versus time for short time-sequences taken from a single measurement over a long time. (a) The first plastic instability appears at about 16% deformation. (b)-(c) The average time $`\tau `$ between successive plastic instabilities diminishes when increasing the deformation. (d)-(e) For large deformations, $`\tau `$ increases and becomes more periodic.

The burst history of a sample is summarized in Fig. 2, which shows five segments of $`d\sigma /dt`$ versus time, taken along a single stress versus strain curve measured at constant compression speed $`v=0.18\mu `$m/s. The first burst appears at about 16% deformation (about 6000 sec after the beginning of the experiment). Above this lower critical deformation, in Fig. 2(a), we observe intense bursts which come in groups of one, two, two, three, and so on. We follow this jerky dynamic regime up to very large deformations. When increasing the deformation, $`\tau _n`$ initially tends to diminish, as the number of events per unit time increases (Fig. 2(b,c)); but beyond about 35% deformation, $`\tau _n`$ tends to increase (Fig. 2(d,e)). At low deformations, below 35%, the sequence of bursts seems random, possibly chaotic . For large compressions, the sequence of bursts develops a periodic component (Fig. 2(d,e)). We do not observe here the “Gutenberg-Richter” power law in the distribution of burst sizes, as in ref. , but a bell-shaped distribution instead.

The variations of the waiting time $`\tau _n`$ between successive events with the deformation, $`\epsilon _n,`$ is shown in Fig. 3. These data (as well as other measurements in different samples at similar compression speeds) show that $`\tau _n`$ diverges, at both low and high deformations, according to power-laws of the form $`\tau _n=A_1\left(ϵ_nϵ_{c1}\right)^{\beta _1}`$, and $`\tau _n=A_2\left(ϵ_{c2}ϵ_n\right)^{\beta _2}`$. For the divergence at small deformations, we find $`ϵ_{c1}`$=0.13 and $`\beta _1`$=0.58. The second divergence occurs at very large deformations, with $`ϵ_{c2}`$=1.1 and $`\beta _2`$=1.8. In Fig. 3 we show that the system undergoes transitions between jerky and laminar dynamic phases, and exhibits large fluctuations and critical slowing down in the number of bursts per unit time in the jerky phase close to two non-equilibrium dynamic critical points. The divergences are of the form $`\tau \left|ϵϵ_{ci}\right|^{\beta _i}`$; where $`t_i=\left|ϵϵ_{ci}\right|`$, ($`i=1,2`$), behave as the analog of the “reduced temperature” for each one of the two phase boundaries separating the different dynamic phases. Surprisingly, the system recovers the laminar plastic flow in the very large deformation regime.

The effect of the compression speed $`v`$ on the jerky regime is analyzed in Fig. 4, which shows four segments of $`d\sigma /dt`$ versus time, taken at four different compression speeds (in four different samples) at about 25% deformation. By increasing $`v`$, the time $`\tau _n`$ between events decreases roughly following a power law. This is shown in Fig. 5, where the local average $`\tau _n\tau `$ is plotted versus $`v`$. This average for $`\tau _n`$ is local, in the sense that it is taken over about 30 events around a given deformation. In Fig. 5 we show the local average $`\tau _n`$ for 25% and 60% deformation, respectively. Only as a guide to the eye, a power-law ($`\tau =Bv^\alpha `$) fit to the data is also shown. We find $`\alpha 0.8`$ at both low and high deformation. Notice that the parameters describing the critical slowing down close to the two dynamic critical points do not seem to be affected by $`v`$. For example, for twice the compression velocity $`v`$ used in Fig. 3, we still obtain two dynamical critical points. One for small deformations, with $`ϵ_{c1}`$=0.17 and $`\beta _1`$=0.51, and a second divergence for very large deformations, with $`ϵ_{c2}`$=1.07 and $`\beta _2`$=1.6.

Fast driving rates and large deformations are almost equivalent in their effect on the PLC dynamics. Fast driving rates push the system into a more periodic-response regime, as do large deformations, while small driving rates favor an aperiodic (possibly random or chaotic) response. A qualitative microscopic explanation can be given. For very small driving rates, the system has time to undergo complex loading-unloading cycles in which a large number of plastic channels are activated and deactivated. At fast driving rates, few percolation-like channels quickly dominate the dynamics and characteristic times appear. Similarly, as the system is squeezed at large deformations, only few channels remain active and the system acquires a periodic response. Similar behavior has been observed in the plastic instabilities and burst-like cooperative motion of “magnetic flux lines” (instead of our “dislocation lines” here) over pinning centers.

The revolution in the general understanding of equilibrium critical phenomena was preceded by the gradual realization that apparently dissimilar and unrelated phenomena (e.g., chemical, magnetic, and superfluid transitions) shared some commonalities near critical points. More recently, the focus has shifted to systems far from equilibrium or in metastable or steady states, and to the search for common behaviors and trends near their phase transitions. In our samples, we apply a generalized force, the externally imposed stress, and the samples respond by generating intermittent bursts. This reponse is quantified by the rate of burst generation, $`1/\tau `$, which can be seen as the order parameter for the PLC effect; the analog, e.g., of the magnetization $`M`$ in a magnet, or the density of paired electrons in superconductors. In the PLC effect, the order parameter is a temporal average or “current”, the number of bursts per unit time; instead of a spatial average typical for equilibrium systems. Thus, the order parameter for the PLC dynamic phase behaves as $`1/\tau \left|ϵϵ_{ci}\right|^\beta `$, with critical exponent $`\beta `$, as obtained in Fig. 3. The relative deformation, $`\left|ϵϵ_{ci}\right|`$, acts as the reduced temperature. In other words, in our samples: $`1/\tau 0`$ in the jerky PLC phase, and $`1/\tau =0`$ in the laminar phase. For magnets: $`M0`$ in the ordered ferromagnetic phase, and $`M=0`$ in the disordered phase. Microscopic ordering of a spin is induced by the cooperative alignment of the neighboring spins. Microscopic motion of a dislocation is induced by the cooperative motion of the neighboring dislocations, in the PLC phase. As clearly seen in Fig. 3, fluctuations in the order parameter $`1/\tau `$ increase near the two dynamic critical points, in analogy with equilibrium phase transitions.

That a far from equilibrium transition can be described by an order parameter, is already known for lasers. We add the PLC effect. Like in the PLC effect, a laser is a driven system kept in a stationary state out of equilibrium. When the optical pump power crosses a threshold, a coherent electromagnetic (EM) wave is generated. In a laser, the internal EM field induces order in the form of temporally correlated emission events, while random decay events destroy correlations. In our samples, temporally and spatially correlated cooperative jump events occur in the PLC phase (e.g., some nearby dislocations jump, while the rest do not move); while many small uncorrelated dislocations move in the laminar phase. Near the critical point, the laser’s relaxation time increases. In the PLC effect, also the relaxation time increases—in the sense that the average time $`\tau `$ between bursts increases.

In summary, we have studied the PLC effect in alloys, focusing on its critical behavior near its two dynamical phase boundaries with enhanced fluctuations and diverging relaxation times (i.e., “critical slowing down” of the jerky activity). In the large compression regime, the time between successive plastic instabilities diverges, until finally the jerky regime stops. The discovery of critical behavior in the PLC effect is relevant for the general study of non-equilibrium driven systems with loading-unloading cycles, which produce plastic nonlinear instabilities and complex stick-slip dynamics. We do not know to which extent this dynamic behavior is universal, and if it provides a key for predicting the future evolution of jerky nonlinear systems. However, our findings show an up to now unexplored approach for the PLC, possibly valid for other systems which share similar dynamics, including the seismic activity of crustal faults and the starquakes in the solid crust of superdense stars.

We thank J. Martin and J. Bonneville for providing the testing machines, E. Giacometti for his help in using them, and M. Bretz and C. Kurdak for comments on the text. This work is supported by the Swiss National Science Foundation. FN acknowledges partial support from the UM Center for the Study of Complex Systems, and DOE contract No. W-31-109-ENG-38.

FIG. 1. Short time-sequence, taken in the range of $`17\%`$ deformation, of the stress $`\sigma `$ (and also $`d\sigma /dt`$) versus time for an Al-4at.%Mg polycrystalline sample. For all figures, the imposed compression speed $`v`$ is $`0.18\mu `$m/s, unless otherwise indicated.

FIG. 2. Various $`d\sigma /dt`$ versus time for short time-sequences taken from a single measurement over a long time. (a) The first plastic instability appears at about 16% deformation. (b)-(c) The average time $`\tau `$ between successive plastic instabilities diminishes when increasing the deformation. (d)-(e) For large deformations, $`\tau `$ increases and becomes more periodic.

FIG. 3. The time $`\tau _n`$ between successive instabilities, versus the deformation, $`\epsilon _n,`$ and power-law fits shown by dashed lines. Large fluctuations and critical slowing down are evident near the two dynamic critical points.

FIG. 4. Various $`d\sigma /dt`$ versus time for short segments at four different compression speeds $`v`$ in different samples, at about the same deformation of $`25\%`$.

FIG. 5. The average time between bursts versus compression speed $`v`$, and power-law fits. These fits have no bearing into the main results (e.g., large fluctuations and diverging relaxation times, critical slowing down, in the order parameter $`1/\tau `$ near the two dynamical critical points) because the compressions shown, $`25\%`$ and $`60\%`$, are far from the critical points.