text:
# Correlation Analysis With Scale-local Entropy Measures

## 1 Introduction

Entropy is a measure of the size of a data distribution contained within a bounded region (distribution support) of some space. In a thermodynamic context this distribution size is interpreted as the number of quantum states accessible to a dynamical system given macroscopic constraints. More generally, if a measure space is partitioned, the measure distribution size is estimated by the effective number of partition elements given the distribution weighting.

Definition of the space partition is a central element of entropy calculation. The partition is sometimes defined as a small-scale limiting partition of the space (e.g., thermodynamic limit, limit procedures in classical analysis ), sometimes based on properties of the data distribution itself and/or on the analysis goals (as in wavelet analysis ).

We define entropy as an explicit function of the partition definition, a scaled binning of the measure space. We calculate the entropy and related quantities as functions of the partition scale(s), similar to the multi-scale partition approaches used in fractal dimension calculations  and image deconvolution . Unlike some analysis methods we invoke no scale limits. Quantities are defined on bounded scale intervals explictly excluding asymptotic limits; the analysis system is in this sense scale local.

The end result of this approach is an entropy which represents arbitrary data correlations as a distribution on scale, particularly useful for problems where the detailed scaling behavior of correlations over substantial scale intervals is of interest (e.g., condensation, coalescence, critical phenomena, strange attractors), where instrumental effects may distort scale distributions, and where the correlation structure is not simply expressible as a power law or other elementary function.

In this paper we describe precision binning methods, define the basic scale-local entropy measure and generalize other aspects of information theory to define the scale-dependent entropy difference or information between an object distribution and a model reference as a differential correlation measure. Based on scale derivatives of entropy and information we define scale-local dimension and dimension transport as generalizations of conventional counterparts based on limit concepts. We apply these correlation measures to several simulations and real data analysis problems.

## 2 Scale-local Entropy

The entropy definition employed here is based on $`C_q(e)`$, the rank-$`q`$ correlation integral at scale $`e`$. Given a data distribution in a $`d`$-dimensional primary measure space spanned by variables $`\{x_i\}`$, we consider a set of corresponding correlation spaces containing all possible $`q`$-point clusters of data points ($`q`$-tuples). There is one such $`q`$-point distribution (in a $`q`$-fold Cartesian product space) for each unique $`q`$ value. The $`q`$-point correlation integral is the projection of the $`q`$-point distribution onto its difference subspace spanned by $`\{x_ix_j\}`$, integrating over the sum variable(s).  The integration limit of the correlation integral on the difference variables is in the simplest case (isotropic binning) the single scale of the analysis.

The reciprocal of the correlation integral estimates the effective bin number in the $`d(q1)`$-dimensional difference subspace. It’s counterpart in the primary measure space is the $`(q1)^{th}`$ root, the effective bin number in the primary space. Defining entropy as the logarithm of the effective primary-space bin number is consistent with entropy as a logarithmic size measure and is a generalization of the thermodynamic definition, the logarithm of the number of accessible states. The correlation integral can be approximated by binning the primary measure space, and expressed in terms of normalized bin contents $`p_i(e)`$, in which case $`C_q(e)_{i=1}^{M(e)}p_i(e)^q`$. This results in the rank-$`q`$ Rényi entropy

$$S_q(e)\mathrm{log}\left[C_q(e)^{\frac{1}{q1}}\right]\frac{1}{1q}\mathrm{log}\left[\underset{i=1}{\overset{M(e)}{}}p_i(e)^q\right].$$
(1)
The entropy $`S_q`$, the number of occupied bins $`M`$, and the bin probability $`p_i`$ are explicit functions of the binning scale $`e`$. More generally, a non-isotropic binning (one utilizing bins without unit aspect ratio) would imply a multidimensional scale space.

## 3 Scale-local Information

Given scale-local entropy we define scale-local information as a basis for differential comparisons between object and reference distributions, data and model. There are a number of possible information definitions from information theory and topology, with significant differences in performance. We define information here as the difference between entropies for object and reference distributions. This implies that the effective bin number for an object distribution is compared in ratio to that of a reference. Nonzero information implies multiplicative reduction of effective bin number by increased correlation structure in the object distribution relative to the reference (cumulant analysis is an alternative differential approach emphasizing linear or additive reduction of distribution size). Scale-local information is then defined as

$`I_q(e)`$ $`=`$ $`S_{q,\text{ref}}(e)S_{q,\text{obj}}(e).`$ (2)
Information so defined provides a differential comparison between a data or object distribution and a model or reference. For example, the reference distribution for an arbitrary point set would be a distribution with the same number of points which maximally ‘fills’ the bounded support – a uniform random distribution. Since the uniform reference is useful in many applications we derive explicit analytic forms for its entropy and dimension.

## 4 Binning and Dithering

Although these scale-local analysis methods are completely general as to the nature of the measure distribution, for the purpose of illustration we emphasize point sets in this paper. The measure distribution is then a set of points in a $`d`$-dimensional embedding space (e.g., the distribution of particles from a heavy-ion collision in momentum space). The analysis begins by applying a partition to the embedding space. For algebraic simplicity we consider a grid of $`d`$-dimensional cubes, an isotropic binning. At each scale there is a continuum of possibilities for the relative position of the binning system on the embedding space. Differing partition placement effects the analysis in general, and there is no a priori reason to prefer any single placement. Thus, we average (dither) over all partition placements to calculate the entropy at each scale.

We define a dithering phase $`\varphi `$ for each of the $`d`$ embedding-space dimensions. The relationship between the partition system and the measure distribution is controlled by varying $`\varphi `$. To implement dithering we calculate the correlation integral $`_{i=1}^{M(e)}p_i(e)^q`$ of an event $`J`$ times at each scale, incrementing $`\varphi `$ each time. Finally we average over these results to obtain the entropy. The dithered entropy $`S_q(e)`$ is thus

$`S_q(e)`$ $`=`$ $`{\displaystyle \frac{1}{1q}}\mathrm{log}\left[{\displaystyle \underset{i=1}{\overset{M_\varphi (e)}{}}}p_{i,\varphi }^q(e)_\varphi \right].`$ (3)
where $`p_{i,\varphi }(e)`$ is the bin probability of the $`i`$th bin for binning phase $`\varphi `$ and scale $`e`$. $`p_i`$ is normalized so that $`_ip_i=1`$ at each scale $`e`$ and dithering phase $`\varphi `$, the sum taken over all occupied bins in each partition.

## 5 A Simple Example

A simple application of scale-local entropy and dithered binning illustrates the analysis process. We obtain the scaled entropy of a 2D uniform distribution of $`N`$ randomly generated points on a unit-square support for several values of index $`q`$. A Monte Carlo simulation with 50k points is shown in figure 1 along with analysis results for $`q=0,2,5`$.

To interpret these results we consider small-, intermediate-, and large-scale regimes. In the small-scale limit, well below the point-separation scale ($`eL/\sqrt{N}`$), each distribution point occupies a single bin; there are N occupied bins ($`M=N`$), each with bin probability $`p_i=1/N`$. Thus, at small scale the entropy approaches

$$\underset{e0}{lim}S_q(e)=\frac{1}{1q}\mathrm{log}\left[\underset{i=1}{\overset{N}{}}\left(\frac{1}{N}\right)^q\right]=\mathrm{log}N.$$
(4)
At intermediate scales we idealize the uniform point set to a continuum ($`N\mathrm{}`$). The number of occupied bins is then simply $`M=(L/e)^2`$ and the bin probability is $`p_i=1/M`$ giving $`S_q=2\mathrm{log}(L/e)`$. At scales substantially greater than the boundary scale ($`eL`$) the entire distribution is contained in a single bin: $`M=1`$ and $`p_i=1`$. The rank-$`q`$ entropy thus vanishes at scales much larger than the distribution boundary size. The detailed $`q`$-dependence near the particle-separation and boundary scales is amenable to an analytic treatment.

## 6 Derivation of $`S_q(e)`$ for a Bounded Uniform Distribution

Information is a relative quantity, a matter of definition. It is impossible to make an absolute determination of the information content of an arbitrary distribution. Using a maximum-entropy reference we can measure information relative to a distribution which is minimally correlated (given certain constraints). This motivates us to derive an algebraic form for the scale-local entropy of a bounded uniform distribution (BUD). This distribution represents a maximum-entropy hypothesis within a boundary, a maximum filling of the distribution support. This is a correlation reference from which any object distribution may deviate with reduced entropy.111Boundedness is itself a form of correlation in a self-consistent description. The derivation is presented in two parts: scales below and above the boundary scale.

### 6.1 Below the boundary scale: $`eL`$

For partitions below the boundary scale there is at least one bin in the interior of the embedding space. To derive an analytic form for the entropy of a BUD we consider a two-dimensional distribution (and later generalize to $`d`$ dimensions). The BUD is defined on a square support with side length $`L`$. Because the distribution is uniform the probability of finding a point in any given bin is simply determined by the bin area.

Figure 2 shows a BUD binned with a general rectangular binning (thin dark lines). We calculate the bin contents for each bin (as a function of scale) and integrate over all possible dithering configurations to determine the analytic form of the entropy. For bins that are entirely contained within the embedding space (interior bins) the bin probability is trivial: $`p_i=e_\text{x}e_\text{y}/L^2`$ is the area of the bin divided by the total area of the support, independent of bin dithering. For edge and corner bins the problem is more complicated. We consider each bin type – interior (white), edge (light grey) and corner (dark grey) – separately.

Corner bins contain both an x- and y-axis support edge. For all dithering configurations there are four corner bins. The effective area of these corner bins is $`e_\text{x}e_\text{y}`$ scaled down by the fraction of the bin along the $`i`$th axis ($`i[1,d]`$) that overlaps the data. By definition, the phase along an axis is $`\varphi _i=(1\alpha _i)`$, where $`\alpha _i`$ is the bin overlap fraction for the first bin on axis $`i`$ ($`\varphi [0,1]`$). We define $`\mathrm{\Delta }_i(L/e_i)\text{int}(L/e_i)`$ and express the amount of overlap between the last bin on axis $`i`$ and the edge of the support along that axis as $`\theta _i=\varphi _i+\mathrm{\Delta }_i\text{int}(\varphi _i+\mathrm{\Delta }_i)`$. With these definitions, calculating the contribution to the correlation integral of the corner bins is a matter of integrating over all $`\varphi `$ values using the relevant bin probabilities. Labeling the corner bins from right to left starting with the upper left bin we write down the $`\varphi `$-dependent corner bin probabilities as

$`p_\text{A}^q`$ $`=`$ $`(1\varphi _\text{x})^q(1\varphi _\text{y})^q\left({\displaystyle \frac{e_\text{x}e_\text{y}}{L^2}}\right)^q`$ (5)
$`p_\text{B}^q`$ $`=`$ $`(1\varphi _\text{x})^q[\varphi _\text{y}+\mathrm{\Delta }_\text{y}\text{int}(\varphi _\text{y}+\mathrm{\Delta }_\text{y})]^q\left({\displaystyle \frac{e_\text{x}e_\text{y}}{L^2}}\right)^q`$
$`p_\text{C}^q`$ $`=`$ $`[\varphi _\text{x}+\mathrm{\Delta }_\text{x}\text{int}(\varphi _\text{x}+\mathrm{\Delta }_\text{x})]^q(1\varphi _\text{y})^q\left({\displaystyle \frac{e_\text{x}e_\text{y}}{L^2}}\right)^q`$
$`p_\text{D}^q`$ $`=`$ $`[\varphi _\text{x}+\mathrm{\Delta }_\text{x}\text{int}(\varphi _\text{x}+\mathrm{\Delta }_\text{x})]^q[\varphi _\text{y}+\mathrm{\Delta }_\text{y}\text{int}(\varphi _\text{y}+\mathrm{\Delta }_\text{y})]^q\left({\displaystyle \frac{e_\text{x}e_\text{y}}{L^2}}\right)^q.`$
We calculate the dither-averaged correlation integral $`C_q=<_ip_i^q>_\varphi `$ by integrating over the different dithering configurations for each bin and summing $`<_ip_i^q>_\varphi =_i<p_i^q>_\varphi `$. Following this approach we integrate the $`p_i^q`$ expressions over the two dithering variables and sum results to calculate the correlation integral. The $`\varphi `$ integral over the first corner bin yields

$`<p_\text{A}^q>_\varphi `$ $`=`$ $`\left({\displaystyle \frac{e_\text{x}e_\text{y}}{L^2}}\right)^q{\displaystyle _0^1}{\displaystyle _0^1}(1\varphi _\text{x})^q(1\varphi _\text{y})^q\varphi _\text{x}\varphi _\text{y}`$
$`=`$ $`\left({\displaystyle \frac{e_\text{x}e_\text{y}}{L^2}}\right)^q\left({\displaystyle \frac{1}{1+q}}\right)^2.`$
The second term is

$`<p_\text{B}^q>_\varphi `$ $`=`$ $`\left({\displaystyle \frac{e_\text{x}e_\text{y}}{L^2}}\right)^q{\displaystyle _0^1}{\displaystyle _0^1}(1\varphi _\text{x})^q[\varphi _\text{y}+\mathrm{\Delta }_\text{y}\text{int}(\varphi _\text{y}+\mathrm{\Delta }_\text{y})]^q\varphi _\text{x}\varphi _\text{y}`$
$`=`$ $`\left({\displaystyle \frac{e_\text{x}e_\text{y}}{L^2}}\right)^q\left({\displaystyle \frac{1}{1+q}}\right)\left[{\displaystyle _{\mathrm{\Delta }_\text{y}}^1}v^qv+{\displaystyle _0^{\mathrm{\Delta }_\text{y}}}v^qv\right]`$
$`=`$ $`\left({\displaystyle \frac{e_\text{x}e_\text{y}}{L^2}}\right)^q\left({\displaystyle \frac{1}{1+q}}\right)^2.`$
The third and fourth terms are similar to the first and second; each of the four corner bins contributes an $`(e_\text{x}e_\text{y}L^2)^q(1+q)^2`$ term to the total dither-averaged correlation integral.

The contributions from edge bins are simpler to calculate. The overlap fraction is unity in the direction parallel to the support edge; along that axis we merely count the number of edge bins. The integral along the second axis is similar to the corner bins. Again there are four terms, but symmetry simplifies the problem

$`<p_{\text{edge}}^q>_\varphi `$ $`=`$ $`\left({\displaystyle \frac{e_\text{x}e_\text{y}}{L^2}}\right)^q{\displaystyle _0^1}{\displaystyle _0^1}\left[\text{int}\left({\displaystyle \frac{L}{e_\text{x}}}+1+\varphi _\text{x}\right)2\right](1\varphi _\text{y})^q\varphi _\text{x}\varphi _\text{y}`$
$`=`$ $`\left({\displaystyle \frac{e_\text{x}e_\text{y}}{L^2}}\right)^q{\displaystyle _0^1}{\displaystyle _0^1}\left[\text{int}\left({\displaystyle \frac{L}{e_\text{x}}}+1+\varphi _\text{x}\right)2\right][\varphi _\text{y}+\mathrm{\Delta }_\text{y}\text{int}(\varphi _\text{y}+\mathrm{\Delta }_\text{y})]^q\varphi _\text{x}\varphi _\text{y}`$
$`=`$ $`\left({\displaystyle \frac{e_\text{x}e_\text{y}}{L^2}}\right)^q\left({\displaystyle \frac{L}{e_\text{x}}}1\right)\left({\displaystyle \frac{1}{1+q}}\right).`$
which is the expression for the x-axis border bins. Contributions from the y-axis bins are obtained by switching indices. The contribution from edge bins is thus

$$2\left(\frac{e_\text{x}e_\text{y}}{L^2}\right)^q\left[\left(\frac{L}{e_\text{x}}1\right)\left(\frac{1}{1+q}\right)+\left(\frac{L}{e_\text{y}}1\right)\left(\frac{1}{1+q}\right)\right].$$
(9)
There remains the integral over interior bins, a simple matter of bin counting

$$<p_{\text{interior}}^q>_\varphi =\left(\frac{e_\text{x}e_\text{y}}{L^2}\right)^q\left(\frac{L}{e_\text{x}}1\right)\left(\frac{L}{e_\text{y}}1\right).$$
(10)
The full correlation integral can be assembled from the corner, edge, and interior bin integrals

$`C_q(e)=<{\displaystyle \underset{i}{}}p_i^q>_\varphi ={\displaystyle \underset{i}{}}<p_i^q>_\varphi `$ (11)
$`=\left({\displaystyle \frac{e_\text{x}e_\text{y}}{L^2}}\right)^q\left[4\left({\displaystyle \frac{1}{1+q}}\right)^2+2\left({\displaystyle \frac{1}{1+q}}\right)\left({\displaystyle \frac{L}{e_\text{x}}}1\right)+2\left({\displaystyle \frac{1}{1+q}}\right)\left({\displaystyle \frac{L}{e_\text{y}}}1\right)+\left({\displaystyle \frac{L}{e_\text{x}}}1\right)\left({\displaystyle \frac{L}{e_\text{y}}}1\right)\right]`$
$`=\left({\displaystyle \frac{e_\text{x}e_\text{y}}{L^2}}\right)^{q1}\left[1+\left({\displaystyle \frac{1q}{1+q}}\right)\left({\displaystyle \frac{e_\text{x}}{L}}\right)\right]\left[1+\left({\displaystyle \frac{1q}{1+q}}\right)\left({\displaystyle \frac{e_\text{y}}{L}}\right)\right].`$
Inserting this result into the definition of the rank-$`q`$ entropy we find that

$`S_q(e)={\displaystyle \frac{1}{1q}}\mathrm{log}[C_q(e)]`$ (12)
$`={\displaystyle \frac{q1}{1q}}\mathrm{log}\left({\displaystyle \frac{e_\text{x}e_\text{y}}{L^2}}\right)+{\displaystyle \frac{1}{1q}}\mathrm{log}\left[1+\left({\displaystyle \frac{1q}{1+q}}\right)\left({\displaystyle \frac{e_\text{x}}{L}}\right)\right]+{\displaystyle \frac{1}{1q}}\mathrm{log}\left[1+\left({\displaystyle \frac{1q}{1+q}}\right)\left({\displaystyle \frac{e_\text{y}}{L}}\right)\right]`$
$`=\left\{\mathrm{log}\left({\displaystyle \frac{L}{e_\text{x}}}\right)+{\displaystyle \frac{1}{1q}}\mathrm{log}\left[1+\left({\displaystyle \frac{1q}{1+q}}\right){\displaystyle \frac{e_\text{x}}{L}}\right]\right\}+\left\{\mathrm{log}\left({\displaystyle \frac{L}{e_\text{y}}}\right)+{\displaystyle \frac{1}{1q}}\mathrm{log}\left[1+\left({\displaystyle \frac{1q}{1+q}}\right){\displaystyle \frac{e_\text{y}}{L}}\right]\right\}.`$
Because of the additivity of entropy under product this entropy expression contains two terms, one for each axis, which suggests the $`d`$-dimensional generalization: adding an equivalent term for each additional dimension.

### 6.2 Above the boundary scale: $`eL`$

For the derivation of Eq. (12) we considered contributions from three types of bins: corner, border, and interior. In the 2D derivation this approach is only valid when both $`e_\text{x}L`$ and $`e_\text{y}L`$. If either of the scales is larger than the support size there are no interior bins and Eq. (12) is not valid. Thus, to obtain the entropy expression for a BUD at large scale consider a single axis (we now exploit the additivity of entropy with respect to dimension for uncorrelated systems) with $`eL`$. When there is a bin edge in the support there are exactly two corner bins; we express the bin probability of the second bin in terms of the first

$`p_1^q`$ $`=`$ $`\left[\alpha \left({\displaystyle \frac{e}{L}}\right)\right]^q`$ (13)
$`p_2^q`$ $`=`$ $`[1p_1]^q.`$
When the support fits completely inside a single bin $`p_0=1`$; the distance from the support edge to the nearest bin edge is larger than the size of the support ($`\alpha eL`$). We now evaluate the relevant integrals

$$<p_0^q>_\alpha =_{\frac{L}{e}}^1[1]^q\alpha =1\frac{L}{e},$$
(14)
$$<p_1^q>_\alpha =_0^{\frac{L}{e}}\alpha ^q\left(\frac{e}{L}\right)^q\alpha =\frac{L}{e}\left(\frac{1}{1+q}\right)=<p_2^q>_\alpha ,$$
(15)
Since the definition of corner bins is arbitrary, we could relabel and get the same result; symmetry requires that $`<p_1^q>_\alpha =<p_2^q>_\alpha `$. These results can now be assembled to calculate the 1D scaled entropy for scales larger than the support size

$`S_q(e)`$ $`=`$ $`{\displaystyle \frac{1}{1q}}\mathrm{log}[C_q(e)]`$
$`=`$ $`{\displaystyle \frac{1}{1q}}\mathrm{log}\left[1{\displaystyle \frac{L}{e}}+{\displaystyle \frac{L}{e}}\left({\displaystyle \frac{2}{1+q}}\right)\right]`$
$`=`$ $`{\displaystyle \frac{1}{1q}}\mathrm{log}\left[1+\left({\displaystyle \frac{1q}{1+q}}\right){\displaystyle \frac{L}{e}}\right].`$
Combining results, we obtain the rank-$`q`$ scaled entropy for each degree of freedom of a BUD over all scales

$$S_q(e)=\{\begin{array}{cc}\mathrm{log}\left(\frac{L}{e}\right)+\frac{1}{1q}\mathrm{log}\left[1+\left(\frac{1q}{1+q}\right)\frac{e}{L}\right],\hfill & \text{for }eL\hfill \\ \frac{1}{1q}\mathrm{log}\left[1+\left(\frac{1q}{1+q}\right)\frac{L}{e}\right],\hfill & \text{for }eL\text{.}\hfill \end{array}$$
(17)
We have generated the exact expression for the scale-local entropy in the general case of a $`d`$-dimensional bounded uniform distribution.

## 7 Point Sets as BUDs

Eq. (17) was derived in the continuum limit, it does not describe a uniform random point set. We generalize to the reference for discrete distributions by introducing an additional factor in Eq. (17). The entropy of a bounded random uniform distribution (BRUD) follows the BUD reference at larger scales but eventually approaches $`\mathrm{log}(N)`$ as a limiting value for scales smaller than the typical two-point separation. The BUD reference can be generalized to the discrete case if we account for the appearance of void bins at smaller scales. The fraction of void bins at scale $`e`$ for a Poisson-distributed point set is $`\mathrm{exp}[\mu _0(e)]`$, where $`\mu _0(e)Ne/L`$ is the average bin occupancy for all bins within the boundary. The fraction of occupied bins (the support of the point set) is then $`f(e)=1\mathrm{exp}(Ne/L)`$. Incorporating $`f(e)`$ into the BUD entropy yields the BRUD result

$$S_q(e)=\{\begin{array}{cc}\mathrm{log}\left(\frac{L}{e}\right)+\frac{1}{1q}\mathrm{log}\left[1+\left(\frac{1q}{1+q}\right)\frac{e}{L}\right]+\mathrm{log}\left[1\mathrm{exp}\left(\frac{Ne}{L}\right)\right],eL\hfill & \\ \frac{1}{1q}\mathrm{log}\left[1+\left(\frac{1q}{1+q}\right)\frac{L}{e}\right]+\mathrm{log}\left[1\mathrm{exp}\left(\frac{Ne}{L}\right)\right],eL.\hfill & \end{array}$$
(18)
This entropy is precise for a discrete, random point set and $`q=0`$. For $`q>0`$ the reference entropy differs from the entropy for real point sets in a small scale interval near $`eL/\sqrt{N}`$, the typical two-point separation. This information (entropy difference) derives from the fact that the form of Rényi entropy employed here is based on ordinary moments (averaged powers of $`p_i(e)`$), whereas the point-set data are Poisson distributed. A correlation integral based on factorial moments would give zero information for uniform (uncorrelated) point sets .

We could conclude that a factorial-moment approximation to the correlation integral should be used for point sets at least, if not other applications. However, extending our previous observation that boundedness itself is a form of correlation in a self-consistent system we can also view the point set as a result of increasing correlation (coalescence) of a continuous distribution at a characteristic scale. The point set does have correlations additional to the boundedness of the continuous BUD. The apparent discrepancies in the entropies of BRUD and real points sets (the nonzero $`I_q`$ in Fig. 3) reveal a genuine correlation feature in a discrete point set compared to a continuum. The Rényi entropies based on ordinary moments are preferred as a more general formulation applicable to arbitrary measure distributions. The information corresponding to the difference between factorial moments and ordinary moments for the uniform random point set is meaningful and can be expressed analytically (the subject of a future paper).

## 8 Scale-local Dimension

Having obtained entropy and information as scale-local distributions we can similarly express the dimension of a distribution as a function of scale. We start with a conventional definition of dimension based on asymptotic limits

$$d_q\frac{1}{q1}\underset{e0}{lim}\left[\frac{\mathrm{log}\left(_{i=1}^{M(e)}p_i^q(e)\right)}{\mathrm{log}(e)}\right]=\underset{e0}{lim}\left[\frac{S_q(e)}{\mathrm{log}(e)}\right].$$
(19)
In this approach dimension is a single number defined in the limit of a zero-scale partition. The slope of the $`S_q(e)`$ distribution is evaluated in the limit of zero scale. This definition favors certain specific types of correlation (power laws), and assumes that the limit exists, which may not be true even in principle. For more general cases the results can be misleading. We relax the asymptotic limit restriction as we did with scale-local entropy to obtain a more general scale-local dimension

$`d_q(e)`$ $`=`$ $`\underset{\delta e0}{lim}\left[{\displaystyle \frac{S_q(e+\delta e)S_q(e)}{\mathrm{log}(e+\delta e)\mathrm{log}(e)}}\right]={\displaystyle \frac{[S_q(e)]}{[\mathrm{log}(e)]}}.`$ (20)
Applying this definition to the generalized entropy of a BRUD in Eq. (18) yields

$$d_q(e)=\{\begin{array}{cc}1\frac{1}{1q}\frac{\left(\frac{1q}{1+q}\right)\frac{e}{L}}{1+\left(\frac{1q}{1+q}\right)\frac{e}{L}}\frac{\frac{Ne}{L}\mathrm{exp}(\frac{Ne}{L})}{1\mathrm{exp}(\frac{Ne}{L})},eL\hfill & \\ \frac{1}{1q}\frac{\left(\frac{1q}{1+q}\right)\frac{L}{e}}{1+\left(\frac{1q}{1+q}\right)\frac{L}{e}}\frac{\frac{Ne}{L}\mathrm{exp}(\frac{Ne}{L})}{1\mathrm{exp}(\frac{Ne}{L})},eL.\hfill & \end{array}$$
(21)
Dimension expressed as a scale-local distribution is a novel aspect of this entropy treatment. To understand how scale can affect the inferred dimension of an object consider the apparent dimension of a planet in the solar system from different viewpoints. For an observer on one planet other planets appear to the eye as isolated points (with zero dimension). With a powerful telescope the resolution size (scale) of the observation decreases substantially. Planets appear as 2D disks with 1D border. A radar probe orbiting a planet can determine that the planet at smaller scale has a rich 3D surface and internal structure that could not be supported by a 1D point or 2D plane. Continuing to the atomic scale planets are made of atoms and molecules that appear point-like. At this scale a planet’s dimensionality returns to zero. This general principle is illustrated in figure 4.

## 9 Analysis of Hierarchically Organized Point Sets

As an exercise in precision correlation analysis using scale-local entropy and dimension we model cluster formation via condensation by generating a hierarchical point distribution with correlation features distributed over a range of scales. This model is relevant to phase transitions and complex systems analysis. To create a two-dimensional, two-tier cluster hierarchy on a square region of side $`L`$ we generate a uniform random distribution of $`N_0`$ cluster sites, providing correlations at the characteristic length scale $`L/\sqrt{N_0}`$ – the mean site separation. At each cluster site we throw a randomly generated uniform distribution of $`N_1`$ points with width $`\delta _1`$, giving the distribution a second characteristic length scale $`\delta _1/\sqrt{N_1}`$.

If the two tiers of the hierarchy are sufficiently separated on scale, as in Fig. 5, the sub-structure of the clusters is not evident to the analysis at large scale ($`e\delta _1`$) . The distribution appears to be a BRUD of $`N_0`$ random points. At smaller scale ($`e\delta _1`$) the apparent structure is dominated instead by the internal cluster structure, a BRUD of $`N_1`$ points. This two-tiered hierarchical distribution is a first approximation to a self-similar distribution: it appears as a BRUD at scales $`L`$ and $`\delta _1`$ simultaneously (see Fig. 5). Extending the hierarchy by recursive self-similar cluster generation would converge to the limiting case of a fractal point distribution over an arbitrarily large scale interval (e.g., Cantor set). The ‘fractal’ dimension would depend on the relations among the hierarchy scale separation, the cluster size and the point count.

The lower right panel of figure 5 shows the dimension transport for the two-tier hierarchy. Dimension transport is defined as the scale derivative of information, $`\mathrm{\Delta }d_q(e)[I_q(e)]/[\mathrm{log}(e)]`$, and is a measure of the scale-dependent dimension difference between reference and object distributions. Dimension transport measures increasing correlation as the transport of dimensionality from larger to smaller scale. In the example of the two-tier hierarchy correlation (relative to the BRUD of the same multiplicity) anticorrelation of points at larger scale is achieved when points condense toward the cluster sites. The anti-correlation at larger scale results in a reduction of larger-scale dimensionality (the system appears more point-like) and an increase in the local point density and dimensionality at smaller scale.

The extended two-tier condensation example in figure 6 shows how small-scale correlations increase by condensing points of a BRUD onto cluster sites. At the onset of cluster formation ($``$3000 cluster sites, $``$3 points per cluster) the transport of dimension to smaller scale is barely visible (but still non-statistical). When the size of the clusters becomes significant ($``$1000 cluster sites, $``$10 points per cluster) the analysis indicates what the eye perceives directly, that the distribution of points is in some way correlated. When the cluster size is 10% of the number of clusters ($``$300 cluster sites, $``$30 points per cluster) the dimension transport shows quite dramatically and quantitatively the transport of dimension from larger to smaller scale.

## 10 Summary and Conclusions

We have developed a novel analysis system which is well-suited to the task of precision correlation analysis for general measure distributions and especially for systems which exhibit clustering or other self-similar behavior. By extending the Rényi-entropy concept to a locally-defined function of scale we are able to establish a more complete picture of data correlations and make precision comparisons among data, simulations and model distributions in the context of information theory. Comparison of Monte Carlo results and analytic distributions have led to a detailed understanding of scale-local entropy measures. Analysis of simulated clustering data suggests the power of this method in the quantification of scale-dependent correlation structure.

## 11 Acknowledgments

The authors would like to thank all of the people who have contributed to and supported this work. In particular, we would like to thank Dhammika Weerasundara who was instrumental in the development of scale-local entropy methods, as well as Stephen Bailey, Justin Prosser, and Curtis Reynolds who helped implement several applications of this analysis.

## 12 References


format: text

text:
# Experimental Identification of the Kink Instability as a Poloidal Flux Amplification Mechanism for Coaxial Gun Spheromak Formation

Now at: \]Los Alamos National Laboratory, Los Alamos, NM 87545

(April 28, 2003; accepted by Phys. Rev. Lett.)

## Abstract

The magnetohydrodynamic kink instability is observed and identified experimentally as a poloidal flux amplification mechanism for coaxial gun spheromak formation. Plasmas in this experiment fall into three distinct regimes which depend on the peak gun current to magnetic flux ratio, with (I) low values resulting in a straight plasma column with helical magnetic field, (II) intermediate values leading to kinking of the column axis, and (III) high values leading immediately to a detached plasma. Onset of column kinking agrees quantitatively with the Kruskal-Shafranov limit, and the kink acts as a dynamo which converts toroidal to poloidal flux. Regime II clearly leads to both poloidal flux amplification and the development of a spheromak configuration.

The spheromak Jarboe94 ; Bellan00  is a simply-connected plasma configuration in which the magnetic fields are largely determined by dynamo-driven plasma currents. Because of its topological simplicity and ease of formation, the spheromak is of interest as a magnetic fusion confinement scheme Mclean02 . Spheromak formation has traditionally been explained by Taylor’s hypothesis Taylor86  that a turbulent magnetohydrodynamic (MHD) system relaxes to a state of minimum magnetic energy subject to the constraint of constant magnetic helicity. While this hypothesis has successfully explained the existence and many equilibrium properties of spheromaks, it says nothing about the actual 3D dynamics underlying the relaxation process, e.g., for coaxial gun spheromaks, the dynamo mechanism which converts injected toroidal flux into required poloidal flux al-Karkhy93 . The dynamics must be 3D because Cowling’s theorem Cowling34  shows that purely axisymmetric processes cannot accomplish this. This Letter experimentally identifies the MHD kink instability Freidberg87  as a 3D mechanism which converts toroidal to poloidal flux in a coaxial gun system, thereby leading to spheromak formation. This mechanism should also be of fundamental importance to coaxial helicity injection in spherical tori Raman03 , relaxation in reversed-field pinches Taylor86 , solar coronal plasma instabilities Vrsnak91 , and astrophysical jets Li01 .

Plasmas in this experiment fall into three regimes depending on peak $`\lambda _{\mathrm{gun}}=\mu _0I_{\mathrm{gun}}/\psi _{\mathrm{gun}}`$ (where $`I_{\mathrm{gun}}`$ is the gun current and $`\psi _{\mathrm{gun}}`$ is the bias poloidal magnetic flux intercepting the inner gun electrode), with (I) low values resulting in a straight plasma column with helical magnetic field along the symmetry axis, (II) intermediate values leading to kinking of the column axis, and (III) high values leading immediately to a detached plasma with $`B_{\mathrm{tor}}B_{\mathrm{pol}}`$. Onset of column kinking agrees quantitatively with the Kruskal-Shafranov limit Freidberg87 , and the kink acts as a dynamo which converts toroidal to poloidal flux. Regime II clearly leads to both poloidal flux amplification and magnetic field profiles consistent with spheromak formation. These results are qualitatively consistent with recent time-dependent resistive MHD numerical simulations of electrostatically driven spheromaks Sovinec01 .

Early coaxial gun experiments Lindberg64  demonstrated poloidal flux amplification, and it was postulated that this was due to an observed helical instability. More recently, toroidal mode number $`n=1`$ central column instabilities in coaxial gun experiments were reported by several spheromak research groups Knox86 ; Browning92 ; Nagata93 ; Duck97 ; Brennan99 . Typically, the observations were based on edge magnetic measurements, and the mode was studied in the context of relaxation during sustainment. Two hypotheses were proposed to explain the mode: (1) development of a $`q=1`$ surface in the closed-flux region Knox86 ; Nagata93  or equivalently a magnetic axis kink, and (2) current-driven instability of the central column Duck97 ; Brennan99  or equivalently a geometric axis kink. Additionally, the mode was shown to couple power from the central column to the spheromak Browning92 . The present work offers a significant new result in that a nonlinear $`n=1`$ geometric axis kink is directly observed to produce the dynamo and poloidal flux amplification leading to spheromak formation.

Building on a prior Caltech spheromak formation experiment Yee00 , the present experiment simplifies, and makes more accessible, the spheromak formation process by using a novel coaxial gun Hsu02a . A schematic of the gun is shown in Fig. 1, along with the cylindrical coordinate system. The gun consists of two concentric, co-planar copper electrodes: a 20 cm diameter disk (blue) biased to negative polarity high voltage surrounded by a 50 cm outer diameter annulus (green) connected to vacuum chamber ground. The width of the gap between electrodes is 6 mm. The inner electrode is mounted on the end of a vacuum re-entrant port; both are electrically insulated from the vacuum chamber by a ceramic break. The co-planar arrangement of the electrode surfaces is geometrically simple, allows the spheromak formation process to be diagnosed with no mechanical obstructions, and is feasible to model numerically Tokman02 . Bias field poloidal flux $`\psi `$-contours are also shown in Fig. 1; $`\psi _{\mathrm{gun}}`$ can be adjusted from 0–7 mWb. Both $`I_{\mathrm{gun}}`$ and bias field point toward the inner electrode. The rise time of $`\psi _{\mathrm{gun}}`$ is 10 ms, making it essentially stationary on the much faster time scale of the gun discharge (tens of $`\mu `$s). The gun is installed at one end of a much larger vacuum chamber (about 2 m long and 1.5 m diameter), and thus boundary effects on the spheromak formation process are minimized. The gun is powered by an ignitron-switched 120 $`\mu `$F, 20 kV capacitor bank. Hydrogen gas is injected transiently using fast puff valves at eight equally spaced toroidal positions on each electrode. Due to the Paschen effect, the optimum path for plasma breakdown is along the bias field and not at the gap between electrodes. The capacitor bank is discharged at $`t=0`$, at which time the bias field and gas puff have already been introduced, and breakdown occurs at approximately 4 $`\mu `$s.

The main diagnostics are a multiple-frame fast CCD camera and a 60-channel magnetic probe array. The camera takes a sequence of time-resolved images in one plasma shot. The capability to visually follow the global evolution of a plasma in one plasma shot is a new and unique aspect of this work. The inter-frame time is typically 1 or 1.5 $`\mu `$s, and the exposure time of each frame is 10 or 20 ns. False color is applied to the images for viewing. The magnetic probe array, shown in Fig. 1, measures all three components ($`R`$, $`\varphi `$, $`Z`$) of $`𝐁`$ at 20 radial positions with 2 cm radial spacing. Probe $`\dot{B}`$ signals are acquired using a digital acquisition system and integrated numerically on a computer. For all $`B`$ measurements in this paper, the probe is located at $`Z=22`$ cm from the plane of the gun electrodes. Kink occurrence is independent of the probe (see Figs. 2 and 3). Propagation of the plasma in $`Z`$ past the stationary probe relates temporal information to $`Z`$-spatial information. This relationship improves as the propagation past the probe becomes fast compared to the plasma expansion rate; this was exploited in the prior Caltech spheromak experiment Yee00 . The time dependence thus acts as a proxy for the $`Z`$-dependence. $`I_{\mathrm{gun}}`$ is measured with a Rogowski coil surrounding the ceramic break, and gun voltage $`V_{\mathrm{gun}}`$ is measured using a high-voltage probe. Typical parameters are: $`V_{\mathrm{gun}}=4`$–6 kV (charge voltage) and 2–2.5 kV (after breakdown), peak $`I_{\mathrm{gun}}=70`$–120 kA, $`\psi _{\mathrm{gun}}=0.5`$–2 mWb, $`B0.1`$–1 kG, $`n10^{14}`$ cm-3, and $`T_\mathrm{e}T_\mathrm{i}520`$ eV.

Figure 2 shows the time evolution of a kinked plasma (shot 2472). Gun electrodes are on the right side of each frame, and the $`Z`$-axis is oriented horizontally across the middle of each frame. At 7 $`\mu `$s, bright arches are visible soon after breakdown, showing that breakdown occurs along the vacuum bias field lines. As $`I_{\mathrm{gun}}`$ increases, the arches expand and quickly merge together (8.5 $`\mu `$s), forming a plasma column (10 $`\mu `$s) which begins to kink (11.5 $`\mu `$s) and then becomes strongly kinked (13 $`\mu `$s). It will be shown that this sequence leads to a spheromak. Depending on peak $`\lambda _{\mathrm{gun}}`$, three distinct plasmas result, as shown in Fig. 3. Regime I leads to a stable column, II to a kinked column and then a spheromak, and III to an immediately detached plasma. The transition from regime I to II (II to III) occurs for peak $`\lambda _{\mathrm{gun}}40`$ (60) m-1.

In order to show that the helical perturbation is an MHD kink mode, consider the time evolution of the $`q(R)`$ profile (relative to geometric axis), which is shown in Fig. 4 (shot 2472). The kink mode becomes linearly unstable when $`q(R)=2\pi RB_\mathrm{Z}/LB_\varphi =1`$, where $`L`$ is the column length; this is the Kruskal-Shafranov limit. Since the kink involves a shift of the current channel, this requires $`q=1`$ on axis and in its vicinity. As seen in Fig. 4, $`q`$ near the axis is greater than unity at 9.5 $`\mu `$s but flattens and approaches unity by 11.5 $`\mu `$s, coinciding with onset of column kinking (Fig. 2). The kink is fully developed at 13 $`\mu `$s and has broken apart by 14.5 $`\mu `$s. The $`q`$ profile is determined using local $`B_\mathrm{Z}`$ and $`B_\varphi `$ data from the magnetic probe array, and $`L`$ is determined from the CCD images. In addition, for an ensemble of plasma shots in which a column forms (regimes I and II), the Kruskal-Shafranov limit (recast Hsu02b  as $`\lambda _{\mathrm{gun}}=4\pi /L`$) is a good predictor Hsu02b  of whether a kink actually develops (Fig. 5). Thus, two sets of independent data (magnetic probe measurements of $`B_\mathrm{Z}`$ and $`B_\varphi `$, and Rogowski/flux loop measurements of $`\lambda _{\mathrm{gun}}`$) both give direct evidence that the helical perturbation is an MHD kink mode. It is interesting to note that the observed kinks always have one axial wavelength, even though there is no fixed boundary at the end of the plasma column. One explanation Nakamura01  for this is that there is a strong axial gradient in $`Bn^{1/2}`$, and thus an effective Alfvén speed discontinuity, which acts as a rigid boundary. This conjecture is consistent with the “mushroom cap” on the left side of the kink in Fig. 3.

Next, it is shown that the kink is followed immediately by three signatures of spheromak formation: (1) appearance of closed $`\psi `$-contours \[calculated assuming axisymmetry, $`\psi (R,t)=_0^R2\pi R^{}B_\mathrm{Z}(R^{},t)dR^{}`$\], (2) $`\psi `$-amplification, and (3) magnetic field radial profiles consistent with spheromak formation. It is important to note that the relationship between closed $`\psi `$-contours and closed flux surfaces becomes ambiguous when axisymmetry is broken. Thus, in the presence of a non-axisymmetric rotating kink, closed $`\psi `$-contours indicate closed flux surfaces only in a time-averaged way. Shortly after 13 $`\mu `$s, the kink breaks apart (Fig. 2). This coincides with signatures of spheromak formation as observed in the magnetic probe measurements (Fig. 6). At approximately 13 $`\mu `$s, closed $`\psi `$-contours appear and $`\psi _{\mathrm{max}}`$ is amplified to larger than $`\psi _{\mathrm{gun}}`$ ($`\psi `$-amplification is due mainly to broadening of the $`B_\mathrm{Z}`$ profile after 12 $`\mu `$s). At 15 $`\mu `$s, magnetic field profiles consistent with spheromak formation are observed (Fig. 7). The measured radial profiles of $`B_\mathrm{Z}`$ and $`B_\varphi `$ at 15 $`\mu `$s are compared with Taylor state solutions Taylor86  in cylindrical geometry, i.e. uniform $`\lambda `$ solutions of

$$\times 𝐁=\lambda 𝐁.$$
(1)
The solutions are $`B_\mathrm{Z}J_0(\lambda R)`$ and $`B_\varphi J_1(\lambda R)`$, where $`J_0`$ and $`J_1`$ are Bessel functions of order zero and one, respectively, and the best fit is found for $`\lambda 15`$ m-1 with a radial offset of 4 cm (displacement of spheromak off the geometric axis). The slight disagreement between measured profiles and the Taylor solution is not surprising since the spheromak is expected to be either (1) still undergoing relaxation toward the Taylor state or (2) in a modified relaxed state since it is still being driven by the gun which has peak $`\lambda _{\mathrm{gun}}50`$ m-1.

The kink modifies the direction of current-flow from $`Z`$ (poloidal) to $`\varphi `$ (toroidal), as seen in Fig. 2. Equivalently, it converts toroidal to poloidal flux. This process is paramagnetic, amplifying $`\psi `$ over the initial applied $`\psi _{\mathrm{gun}}=1.7`$ mWb. The paramagnetism is understood by realizing that kinks involve perturbations with dependence $`\mathrm{exp}(i𝐤𝐱)`$ where $`𝐤𝐁=0`$. The latter means $`k_\mathrm{Z}=nB_\varphi /RB_\mathrm{Z}`$. Thus, a locus of constant phase of the kink is given by $`\varphi =(B_\varphi /RB_\mathrm{Z})Z`$, meaning that the kinked current channel is a right (left) handed helix if $`J_\mathrm{Z}B_\mathrm{Z}>0`$ ($`<0`$). This will always lead to amplification of the original $`B_\mathrm{Z}`$. The additional $`\psi `$ introduced by the kink can be estimated by approximating the helix as a solenoid with current $`I`$, turns per length $`1/L`$, and radius $`a`$. The solenoid formula for the field inside the solenoid is $`B_\mathrm{Z}=\mu _0I/L`$; the $`\psi `$ produced by the solenoid is $`\pi a^2B_\mathrm{Z}`$ and thus depends non-linearly on the kink amplitude $`a`$. Using the measured values $`a5`$ cm, $`L20`$ cm, and $`I60`$ kA at 13.5 $`\mu `$s, the $`\psi `$ generated by the kink is predicted to be approximately 1 mWb, which is within a factor of 2 of the observed amplification of $`\psi _{\mathrm{max}}`$ over $`\psi _{\mathrm{gun}}`$. The discrepancy is within the accuracy of $`a`$ and $`L`$ measurements and of the $`\psi `$ calculation assuming axisymmetry in the presence of the rotating kink. Because the coaxial gun can only inject toroidal flux, 3D plasma dynamics must be responsible for $`\psi `$ exceeding $`\psi _{\mathrm{gun}}`$. The dynamics are provided by the kink, which amplifies $`\psi _{\mathrm{gun}}`$ by converting toroidal to poloidal flux; thus, in this case, the kink constitutes the dynamo intrinsic to spheromak formation.

Regime II clearly leads to all three signatures of spheromak formation: closed $`\psi `$-contours, $`\psi `$-amplification, and proper magnetic field profiles. Regime I lacks both closed $`\psi `$-contours and $`\psi `$-amplification. Regime III may have closed $`\psi `$-contours; however, $`B_{\mathrm{tor}}B_{\mathrm{pol}}`$ associated with high $`I_{\mathrm{gun}}`$, which could indicate that the plasma is still relaxing. This work indicates a close relationship between the kink threshold ($`\lambda _{\mathrm{gun}}=4\pi /L`$) and the spheromak formation threshold, which according to a static force-free treatment is also proportional to an inverse length. For example, solutions of Eq. (1) for uniform $`\lambda `$ indicate that closed $`\psi `$-contour equilibria exist for threshold $`\lambda \sqrt{x_{11}^2/a^2+\pi ^2/h^2}`$, where $`a`$ ($`h`$) is the radius (length) of the flux conserver and $`x_{11}`$ is the first root of $`J_1`$ Bellan00 . Geometric differences among different experiments could result in different eigenmode spectra of Eq. (1) relative to the kink stability thresholds.

In summary, the MHD kink instability has been identified experimentally as a poloidal flux amplification mechanism for coaxial gun spheromak formation. An $`n=1`$ central column helical instability was observed during formation using multiple-frame CCD imaging. Onset of the perturbation was shown using two independent sets of data to agree quantitatively with the Kruskal-Shafranov limit. The kink acts as a dynamo which converts toroidal to poloidal flux, and it is followed immediately by three key signatures of spheromak formation: (1) closed $`\psi `$-contours, (2) $`\psi `$-amplification, and (3) magnetic field profiles similar to the Taylor state.

The authors acknowledge S. Pracko, C. Romero-Talamás, and D. Felt for technical assistance and H. Bindslev for asking about kink paramagnetism. Supported by a US-DoE Fusion Energy Postdoctoral Fellowship and US-DoE Grant DE-FG03-98ER544561.


text:
# Switching a spin-valve back and forth by current-induced domain wall motion

on leave from \] the Department of Physics, Faculty of Science, HR-10002 Zagreb, Croatia.

## Abstract

We have studied the current-induced displacement of a domain wall (DW) in the permalloy (Py) layer of a Co/Cu/Py spin valve structure at zero and very small applied field. The displacement is in opposite direction for opposite dc currents, and the current density required to move DW is only of the order of 106 A/cm2. For $`H`$ = 3 Oe, a back and forth DW motion between two stable positions is observed. We also discuss the effect of an applied field on the DW motion.

Switching the magnetic configuration of a micro-device by spin transfer from a spin-polarized current, rather than by applying an external field, is the central idea of a present extensive research Slonczewski ; Berger1 ; Katine ; Sun ; Grollier1 ; Wegrowe ; Tsoi . In 1996, Slonczewski Slonczewski  showed that the magnetic moment of a ferromagnetic layer can be reversed by injecting a spin-polarized current into this layer. This prediction has been convincingly confirmed by series of experiments on pillar-shaped magnetic multilayers Katine ; Sun ; Grollier1 , nanowires Wegrowe  or nano-contacts Tsoi . However, the current density required in the existing experiments is relatively high, of the order of 107 A/cm2, and some reduction of this density is necessary for practical applications.

Another way to change a magnetic configuration is by current-induced motion of a domain wall (DW). DW-drag by a current has been predicted by Berger Berger2  and its theory has been recently revisited by Waintal and Viret Waintal . When a spin-polarized current goes through a DW, the torque, resulting from the interaction of the conduction electron spins with the exchange field in the DW, progressively rotates the spin polarization of the current. Reciprocally, the spin-polarized current exerts an exchange torque on the magnetization within the DW, which is the origin of the DW motion predicted by Berger Berger2 .

Freitas and Berger Freitas  have obtained some experimental evidence of DW-drag by using Kerr microscopy to detect the DW position. In recent similar experiments Gan et al. Gan  have also measured DW displacement due to current pulses by imaging the DW by MFM before and after current pulses. The main features in these two sets of experiments are that the direction of the the DW motion is reversed when the direction of the current pulses is reversed, and that the order of magnitude of the current pulses needed to move the DW is about 107 A/cm2.

In recent experiments on Co/Cu/Py spin valves Grollier2 , we have also found that a dc current can switch the magnetic configuration of the spin valve by moving a DW in the permalloy (Py) free layer. We observed that a DW can be moved away from an artificial pinning center (notch) when the current density exceeds a threshold value of the order of 107 A/cm2. However, the mechanism of the DW displacement was not completely clear. In fact, the DW could not be displaced at zero field, but only by combining current and applied field. Also, the motion direction was determined by the field direction and not reversed when the current was reversed. These results suggested a more complex mechanism than described by Berger Berger2 , with a possible effect of the applied field on the DW distorted by the current. In this paper we present much clearer results obtained on spin valves but with weaker DW pinning. The displacement of the DW is obtained at zero field, in opposite directions for opposite current directions and with definitely lower current densities.

Our samples are 300 nm wide and 20 $`\mu `$m long stripes patterned by e-beam lithography using a lift-off technique. The spin valves, deposited by sputtering, have a final structure CoO(30 $`\AA `$)/Co(70 $`\AA `$)/Cu(100 $`\AA `$)/Py(50 $`\AA `$)/Au(30 $`\AA `$). The top Au electrodes are processed by UV lithography. In contrast with our previous experiments Grollier2 , the only pinning centers for the DW in the Py soft layer are natural defects of the stripe. All the measurements were performed at room temperature.

In Fig. 1 we show a typical GMR minor cycle associated with the reversal of the magnetization in the Py layer, i.e. with the motion of a DW from one end of the Py stripe to the other one. The plateaus are due to the pinning of the DW on defects in the Py stripe. We emphasize that the series of plateaus on the GMR curve is highly reproducible. As shown for the left half of the cycle, a DW remains pinned on the same defect when the field is brought back to zero. We can therefore start an experiment at zero field with the DW pinned in one of the three positions (sketched) corresponding to the resistance levels 1, 2 and 3. The results presented below correspond to experiments performed with a DW initially pinned in the configuration 2. The first series of experiments are performed by varying the current at zero or very low field (parallel to the stripe). In another series of experiments, we study the influence of a larger bias field on the current-induced DW motion.

In Fig. 2 we present results obtained by varying the dc current at constant field close to zero (4 and 3 Oe). As shown in Fig. 2a, starting from the DW in position 2, we can move the DW to position 3 by increasing the current above the positive critical value $`j_{c2}^+`$(4 Oe) = + 0.65 mA and decreasing it back to zero. Alternatively, the DW is moved in the opposite direction (from 2 to 1) with a negative current exceeding $`j_{c2}^{}`$(4 Oe) = $``$ 1.1 mA (in our notation $`j_{cn}^+`$ and $`j_{cn}^{}`$ are the critical currents required to move the DW from position n to positions n+1 and n-1 respectively). The same type of behavior is observed for all applied fields between 0 and 7 Oe. However, even in this very small field range, there is some field dependence of the critical currents: $`j_{cn}^+`$($`H`$) ($`j_{cn}^{}`$($`H`$)) decrease when $`H`$ decreases (increases) and favors a DW motion from n to n+1(n-1).

Fig. 2b presents an example of back and forth DW motion, namely the motion from 2 to 3 with positive dc current and a return to 2 with a negative dc current. The obvious conditions for this back and forth motion are $`j_{c2}^+(H)<j_{c3}^+(H)`$ (required to stop the first motion in configuration 3) and $`\left|j_{c2}^{}(H)\right|>\left|j_{c3}^{}(H)\right|`$ (necessary for the return to the configuration 2). It turns out that these conditions are satisfied for the pinning centers 2 and 3 of our sample, at least for $`H`$ = 3 Oe.

The behavior observed in the field range close to zero (approximately, 0 $`<H<`$ 7 Oe) can be summarized as follows. A DW can be displaced between pinning centers and, in agreement with what is predicted for a displacement by Berger’s mechanism Berger2 , its motion is in opposite directions for opposite currents. The dc current density needed to move the DW is of the order of 106 A/cm2, that is an order of magnitude smaller than the currents required for the magnetization reversal in pillar-shaped multilayers Katine ; Sun ; Grollier1 . There is however some uncertainty in the exact value of the current density in Py. If the electron mean free paths in Py and Co were much smaller than the thicknesses of the Py and Co layers (which is far from being satisfied at room temperature) and if we could also neglect the specular reflections of the electrons at the interfaces, there would be an uniform current density in the multilayer Barthelemy , that is, for example, 8$``$106 A/cm2 for 0.6 mA. In the opposite limit of almost independent conduction by the magnetic and nonmagnetic layers (this would correspond to layer thicknesses larger than the mean free paths, or also to almost complete specular reflections at the interfaces, with, in both cases, a vanishing GMR), a straightforward calculation, based on the resistivity of the different metals at room temperature, leads to a current density of 1.75 $``$106 A/cm2 in Py. In an intermediate situation (mean free path of the order of the Co and Py thickness, consistently with the small but non-zero GMR, and certainly some current channeling in Cu by specular reflections), the real current density in Py is probably in-between, that is of the order of few 106 A/cm2.

Out of the low field range described above, the behavior becomes more complex. An example of experimental result is shown in Fig. 3 for $`H`$ = $``$ 21 Oe favoring an antiparallel (AP) configuration. A positive current moves the DW from position 2 to the end of the stripe (AP resistance level), which is consistent with the motion direction induced by a positive current at low field. On the other hand, in contrast with the low field behavior, the motion is not reversed for negative currents and the final state is still the AP configuration. For positive fields out of the low field range, the same type of behavior is observed, with a motion towards a more parallel configuration. We can therefore conclude that, out of the low field range, the current is still able to unpin the DW, but the direction of the DW motion is now controlled by the applied field direction. This behavior was also observed in our former experiments Grollier2  with stronger (artificial) pinning centers, where the motion could be obtained only by combining current and applied field.

We will now focus on the interpretation of the DW-drag effects at zero or low field. To start with, we can rule out any contribution from Joule heating. From the small quadratic resistance increase with current (the term subtracted in Figs. 2 and 3), the maximum increase of $`T`$ is about 5 K, and we have checked that, at 300 K, this has practically no effect on the GMR minor loop. An even stronger argument is that heating could not explain that opposite currents produce motions in opposite directions. The Oersted field generated by the current ($``$ 20 Oe), in a perfect structure, is in the plane of the DW, and it cannot favor a motion in one or in the other direction. In the presence of defects, the Oersted field might have a component out of the DW plane, but it can be hardly imagined that different defects give always the same direction for this component and the DW motion. Only the spin transfer mechanism, first proposed by Berger, is consistent with the experimental results at zero or very small field and, particularly, can explain the reversal of the motion with opposite currents. Berger Berger2  expresses the spin transfer by a torque corresponding to the field $`H_B=jP\mathrm{}/e\delta M_s`$, perpendicular to the layer. With $`M_\mathrm{s}`$ = 860 $``$ 103 A/m, $`\delta `$ (DW thickness) = 100 nm, P (polarization) = 1 and $`j`$ = 5 $``$ 106 A/cm2, we obtain $`H_\mathrm{B}`$ 3.8 Oe, just in the range of the pinning fields of the DW in the Py layer. Although, at the present stage of the theory, the connection between the torque of $`H_B`$ and the value of the critical current is still unclear, we can conclude that the DW motions at zero or very low field, characterized by a reversal of the motion in opposite currents, can be ascribed to a spin transfer mechanism. On the other hand, the behavior observed at higher fields, with combined influence of current and applied field, is more complex. As suggested by the model of Waintal and Viret Waintal , it could be that, in this regime, the depinning of the DW is induced by the distortion of the wall, while the succeeding motion is predominately driven by the field.

In conclusion, we have presented experiments in which a spin valve is switched by current-induced DW motion. In zero or very low field, the DW displacement is in opposite directions for opposite dc currents, and back and forth motions between two pinning centers can be obtained. Our results are consistent with the spin transfer mechanism introduced by Berger Berger2 . A more complex and unclear behavior is observed when the effect of the current is combined with the effect of an applied field.

This work was supported by the EU through the RTN ”Computational Magnoelectronics” (HPRN-CT-2000-00143) and the Ministère de la Recherche et de la Technologie through the MRT ”Magmem II” (01V0030) and the ACI contract ”BASIC” (27-01).


text:
# Expansion rates at RHIC

## 1 Introduction and Motivation

Experimental data from RHIC , in particular the systematic analysis of anisotropies in the particle spectra, indicate that at high center of mass energies the collision region rapidly thermalizes and for the major part of its evolution subsequently expands according to the laws of ideal hydrodynamics (for a recent review that collects the arguments and references which support this hydrodynamic scenario see ). This is an important finding by itself as apparently the generated system spends some time at temperatures above the anticipated QCD transition temperature and evolves according to the pressure, pressure-gradients and thus equation of state of deconfined strongly interacting matter.

Furthermore the hydrodynamic evolution can be used to shed light even onto some microscopic scattering phenomena, such as signals from hard probes which propagate through the hydrodynamically evolving ’background’ fields and radiate energy, as well as the late freeze-out dynamics of the reaction. Such probes are significantly dependent on the evolution of the macroscopic fields such as temperature and particle density, but also on the collective flow velocities and their gradients.

We will here give a detailed study of the evolution of these fields. First we investigate a simple analytically accessible parametrization of the flow field and its divergence. Then we will discuss the hydrodynamic evolution of thermodynamic fields in a full hydrodynamic calculation to return later to the collective flow fields and compare to the initial analytic example.

## 2 Expansion rate and dilution rate

Ideal hydrodynamics incorporates the continuity equations for conserved charges $`_\mu n^\mu =0`$ where $`n^\mu =nu^\mu `$ is the current associated with the conserved charge whose density distribution is given by $`n(x)`$ and $`u^\mu (x)`$ is the four-velocity of the collective flow field. In ideal hydrodynamics the entropy density is such a conserved quantity. Slightly reshaping the continuity equation leads to the relation

$$_\mu u^\mu =\frac{1}{n}u^\mu _\mu n.$$
(1)
The expression $`_\mu u^\mu `$ is commonly referred to as expansion rate. Considering the trivial condition of vanishing collective flow $`u^\mu =(1,0,0,0)`$ as is the case in the center of the system due to symmetry, the expansion rate is simply given as $`_\mu u^\mu =\dot{n}/n`$, the dilution rate. Assuming furthermore a power-law decay of the density $`n\tau ^\alpha `$ as is the case for the three-dimensional Hubble expansion of the Universe ($`\alpha =3`$) or the one-dimensional longitudinal Bjorken expansion ($`\alpha =1`$) the expansion rate is given by $`_\mu u^\mu =\alpha /\tau `$ ($`\alpha `$ is referred to as expansion parameter).

Let us now investigate a more general situation, namely a longitudinal Bjorken flow field $`v_z=\mathrm{tanh}\eta =z/t`$ with a transverse flow component given as $`\gamma v_T=\mathrm{tanh}(\xi r)`$. We will later see that the hydrodynamic calculations in fact suggest such a flow profile at freeze-out with $`\xi 0.07`$ fm-1. It is quickly shown that this flow field leads to an expansion rate

$$_\mu u^\mu =\tau ^1\mathrm{cosh}(\xi r)+\xi \mathrm{cosh}(\xi r)+r^1\mathrm{sinh}(\xi r),$$
(2)
which in the limit $`(r,z)0`$ leads to $`_\mu u^\mu 1/\tau +2\xi `$ as discussed in . Keeping the radial dependence in the expression one finds that the expansion rate increases with radial distance from 0.207 fm-1 at $`r=0`$ to $`0.247`$ fm-1 at $`r=10`$ fm (assuming $`\tau =15`$ fm/$`c`$, a typical fireball lifetime). This will become important later in this study.

In the case of such a non-zero transverse flow profile the spatial part of the four product on the right hand side of Eq. (1) does not vanish and the expansion rate is no longer identically equal to the dilution rate. This becomes apparent by assuming a boost invariant distribution of the conserved charge (which means that the charge density only depends on the proper time $`\tau =\sqrt{t^2z^2}`$ and not on the space-time rapidity $`\eta `$) and describing its dilution locally in terms of a power law such that $`n(r,\tau )=n_0(r)(\tau _0/\tau )^\alpha `$. It then quickly follows from Eq. (1) that

$$\tau _\mu u^\mu =\gamma \left(\alpha \tau v_r\frac{_rn_0(r)}{n_0(r)}\right).$$
(3)
As the density is expected to drop with increasing radial distance the last term in the brackets is negative and gives a positive correction to $`\alpha `$ as does the preceding $`\gamma `$ factor. It thus follows that in general the expansion rate $`_\mu u^\mu `$ is greater or equal to the local dilution rate $`\alpha /\tau `$.

## 3 Time evolution of thermodynamic fields in central collisions

We now investigate the evolution of the thermodynamic fields using a hydrodynamic model  which explicitly assumes boost-invariance throughout the evolution. The initial field configuration in the transverse plane is determined by an optical Glauber-parametrization  and the parameters of the calculation, namely the equilibration time $`\tau _{\mathrm{equ}}`$ and the constant $`s_{\mathrm{equ}}`$ which relates the geometrical density to the initial entropy density are obtained by fitting to preliminary particle spectra from central Au+Au collisions with $`\sqrt{s_{\mathrm{NN}}}=200`$ GeV. Good agreement is achieved with $`s_{\mathrm{equ}}=110`$ fm-3 at $`\tau _{\mathrm{equ}}=0.6`$ fm/$`c`$ , to be used in the following.

Evolution of entropy and temperature: We first study the time evolution of the entropy density and the temperature at the center of the system and at 3 and 5 fm radial distance (Fig. 1). Initially the entropy decays with a $`\tau ^1`$ power law due to the initial longitudinal Bjorken expansion (see the dashed lines). As time evolves, the transverse pressure gradients lead to a collective radial expansion and thus a transition from the initial one-dimensional to a full three dimensional expansion in the late stage. At late times the dashed line indicates a decay $`\tau ^3`$.

The equation of state supplemented to the hydrodynamic equations of conservation of energy and momentum features a strong first order phase transition at a critical temperature of $`T_c=164`$ MeV. Above that temperature the system acts like an ideal gas with $`c_s^2=p/e=1/3`$ and $`eT^4`$. For a one-dimensional Bjorken expansion the energy density decays as $`e\tau ^{(1+c_s^2)}`$ and thus $`T\tau ^{1/3}`$. The figure shows that this behavior is neatly followed until just before the system hits the mixed phase where the temperature remains unchanged while the system boils off its latent heat. Further cooling appears then in the resonance gas phase below $`T_c`$.

Evolution of radial velocity: Now we consider the collective radial velocity which is generated through the act of the pressure gradients in the system. The left plot of Fig. 2 shows the radial flow velocity at distance 3 and 5 fm from the fireball center. As pressure gradients are larger at 5 fm than at 3 fm, the acceleration is stronger and thus the rise of $`v_r(\tau )`$. When the system reaches the mixed phase, the acceleration ceases and the velocity locally drops as matter continues to move radially outward, thus transporting the slower fluid cells from the interior to larger radial distances. Only as the matter cools into the hadronic phase, pressure gradients pick up again and more radial flow is generated.

The right hand plot of Fig. 2 shows the radial flow field as a function of radial distance at four different times. The generated flow field 1 fm/$`c`$ after impact has an interesting structure, showing a non-monotonous behavior due to the mixed phase where no collective flow develops. As time evolves the field rapidly adopts a linear behavior which is largely preserved for the rest of the evolution, and can very well be described by $`v_r(r)=\xi r`$ ($`\mathrm{tanh}(\xi r)`$) with $`\xi =0.07`$ fm-1. This ’static’ flow field is at the heart of the successful description of a large number of observables in terms of the popular ’blast wave’ parametrization .

Evolution of dilution- and expansion rates: We now return to the dilution- and expansion rates in the full hydrodynamic calculation.

Fig. 3 shows the evolution of the local expansion parameter $`\alpha =\frac{\tau }{s}\frac{s}{\tau }=\frac{\mathrm{log}s}{\mathrm{log}\tau }`$ in the left panel. It is clearly seen how the initial one-dimensional motion with $`\alpha =1`$ becomes more three-dimensional and just reaches $`\alpha =3`$ at the time of decoupling (assuming freeze-out at $`T_d=100`$ MeV). Note the influence of the radial position during late stages of the expansion. At larger radial distances the dilution rate is smaller, as supplies from the interior are constantly streaming outward.

The right panel of the figure shows the evolution of the product of expansion rate and time, $`\tau _\mu u^\mu `$. Note that with increasing the radial distance we get an opposite behavior as exhibited by the dilution rate. As demonstrated by the analytical calculation before, the expansion rate increases with radial distance. The sign of the difference is understood in terms of Eq. (3).

## 4 Time evolution of anisotropies in non-central collisions

Observables from non-central collisions deliver the strongest arguments for rapid thermalization and a subsequent hydrodynamic expansion of the reaction region . The geometric anisotropy in non-central collision $`ϵ_x`$ results in anisotropic pressure gradients which in turn lead to anisotropies in the collective flow field. Those can be characterized in terms of $`ϵ_p`$, the asymmetry of the transverse part of the energy momentum tensor . This flow anisotropy will finally be observable as an azimuthal dependence of particle emission (i.e. elliptic flow) . From the left part of Fig. 4, which shows contours of constant energy density in the transverse plane at different times, it can be seen how the initial spatial eccentricity quickly diminishes.

This becomes more quantitative in the right part of the figure where the evolution of $`ϵ_x`$ and $`ϵ_p`$ for a system with (RHIC 1) and without (EOS I) a first order phase transition is shown. The large initial geometric eccentricity leads to a rapid generation of momentum anisotropy, which in turn reduces the spatial eccentricity, and thus cuts off the further generation of momentum anisotropy. During the soft phase of the transition pressure gradients vanish (and thus also their anisotropies), which leads to an even earlier saturation of the momentum anisotropy. For systems at RHIC most of the anisotropy is thus generated in the plasma stage of the fireball.

## 5 Summary and Relevance

We have given a detailed account on the evolution of the thermodynamic fields in systems created at RHIC within a hydrodynamic framework. This information is relevant for the microscopic dynamics at different levels. At early stages, hard probes are to be transported over the varying thermodynamic background fields to calculate their energy loss. At late stages, the scattering rates in the resonance gas needs to be compared to the macroscopic expansion and dilution rates to investigate the overall reliability of the idealized macroscopic approach which requires a ’slow’ (adiabatic) evolution and strong rescattering. Furthermore there are sensitive probes of the density around decoupling such as light nuclei, whose formation probability favors a large density, but their survival rate increases with decreasing density. Resonances reconstructed from a hadronic decay channel have similar information content . In summary, essentially all experimental observables are influenced by the temporal evolution of the thermodynamic fields presented here either as the signal is sensitive to all stages of the evolution and thus integrates the temporal history of the system or as it results from a particular stage of the collision with densities and their gradients giving the relevant physical background.

Acknowledgments: This research was supported in part by the U.S. Department of Energy under Grant No. DE-FG02-88ER40388. Support from the Alexander von Humboldt Foundation through a Feodor Lynen Fellowship is gratefully acknowledged.


text:
# Freezeout of resonances and nuclear fragments at RHIC

## 1 Introduction

This talks is based on two papers, one with G.Brown  and one ongoing work with P.Kolb , who will report a part of it related with late time evolution here. Their common goal is to understand what are the conditions which determine the timing of production of $`observable`$ resonances, and estimate at what conditions this happens. We have emphasized $`observable`$ above because there are many resonances produced in the system but unobserved in the final invariant mass spectra because of re-scattering of their decay products. The development of simple analytic model of the kinetics of resonance production/absorption, as well as evaluation of the “optical depth” integral.

An old but still important idea is “matter modification” of hadrons, providing the experimental test of the conditions in question. As such we discuss modification of the mass and shape of two classic $`\pi \pi `$ resonances, $`\rho `$ and $`\sigma `$. Recent STAR observations of these effects are reproduced (for $`\rho `$) and predicted (for $`\sigma `$).

The issue of production of nuclear fragments is also an old one, at RHIC reduced to $`d,He^3`$ and their antiparticles. In literature those are studied either by statistical or coalescence models, which left open many important issues. First of all, like resonances the $`observable`$ fragments must escape all interactions with any particles in order to survive at the end. The second point is that this process is production-rate-limited, thus it can lead to non-thermal quasi-equilibria. Furthermore, new element is the consistent evaluation of the 2-to-1 production rate itself in a recent work by Ioffe et al .

## 2 The optical depth factor

Let me start with the simplest pedagogical points about the observability of resonances and fragments, or “composites” as we may call them collectively, for brevity.

In the next section we will discuss rate equations which can be solved and determine the number of composites $`N(t)`$ at time $`t`$. The “observability condition” of a resonance can be written as

$`\nu _{visible}(t)=\mathrm{\Gamma }N(t)exp\left({\displaystyle _t^{\mathrm{}}}\nu (t^{})𝑑t^{}\right)`$ (1)
where the the l.h.s. is the production rate of visible resonances, $`\mathrm{\Gamma }`$ is the resonance decay width and the exponent is the optical depth factor containing integrated $`\nu (t)`$, the combined scattering rate for all decay products. The $`N(t)`$ decreases with time due to expansion and cooling, while the exponent changes from 0 at early time to 1 at late times: so the product naturally has a maximum at the time $`t_m`$ such that

$`{\displaystyle \frac{1}{N(t_m)}}{\displaystyle \frac{dN(t_m)}{dt}}+\nu (t_m)=0`$ (2)
This condition means that for observable resonances the freezeout condition is different from that for stable particles and reads: the rate of their number change is equal to the absorption rate of all the decay products. For example, for $`\rho `$ and $`\sigma `$ we should not know their scattering rates but just that of two pions. Since for short-lived $`\rho `$ and $`\sigma `$ the first factor is close to overall expansion rate of matter at late time which follows from Hubble-like late-time regime $`d\mathrm{log}N(t_m)/dt3/t`$, and the second is the same, we conclude that ‘visible” $`\rho `$ and $`\sigma `$ are produced at the same time.

The formation rate for a fragment made of $`A`$ nucleons is made by some coalescence, and such rate is obviously proportional to nucleon density to that power, $`(n_N(t))^A`$. After it is produced, however, it still has very small probability to survive. Assuming that the destruction rate for a fragment $`\nu _AA\nu _N`$, where $`\nu _N`$ is a scattering rate for one nucleon, one finds that for $`A`$-fragment the time distribution is approximately the $`A`$-th power of the same universal function

$`n_{fragments}(t)\left[n_N(t)exp\left({\displaystyle _t^{\mathrm{}}}\nu _N(t^{})𝑑t^{}\right)\right]^A`$ (3)
So, the maximum of production of any visible fragment happens at the same time for all $`A`$. Furthermore, the width of the distribution over production time decreases as $`A`$ grows, as $`1/\sqrt{A}`$.

## 3 Solving the rate equations

The equations themselves describing dynamics of resonances are well known, generically they contain the sink (the decay) and the source terms

$`{\displaystyle \frac{n(t,\stackrel{}{r})}{t}}=\mathrm{\Gamma }n(t,\stackrel{}{r})+S(n_i)`$ (4)
(where for expanding source the time should be understood as proper time in the rest frame of all volume elements.) In many papers in literature (e.g. ) the source term is ignored citing “instantaneous hadronization”, but (especially for resonances we consider) it is not true: in fact the primary generation of resonances die out long before the “observable” ones are born.

We use an approximate power fit of the source time dependence $`d^3rS=\mathrm{\Gamma }N_0\left(\frac{t_0}{t}\right)^P`$ Its power can be related to fireball expansion. If the volume $`V(t)1/n(t)t^a`$ the integrated source is proportional to $`V(t)[n_\pi (t)]^N`$ where $`N`$ is the multiplicity in resonance decay ($`N=2`$ for $`\sigma ,\rho `$). The pion number is “chemically frozen”, $`N_\pi =V(t)n_\pi (t)=const(t)`$ from which we conclude that the source term power is $`P=(N1)a`$ (for $`\sigma ,\rho `$ and other binary resonances $`P=a`$ ). An example of such equation solved is shown in Fig.1.

We then evaluate the optical factors for pions and nucleons, using realistic re-scattering rates, with chemically frozen composition, using papers by Hung and myself  and by Tomasik and Wiedemann . Example for the final time distributions for visible $`\rho `$ and $`d`$ is shown in Fig.2. Note that both distribution have maxima we discussed above, and that the “visible” $`d`$ are indeed produced very late. This is our main point: we are speaking about a very dilute matter, after the freezeout of all the basic ingredients of the fireball.

## 4 The resonance modification

It has been argued over the years that in matter the resonances should be modified, with shifted mass, increased width and even significantly changed shape. With the very late stages of RHIC collisions relevant, we can now access very dilute matter in which those effects must be calculable in the lowest order of the density, providing a benchmark test to all such discussions. In such case hadron modification is expressed in terms of their forward scattering amplitude $`M_{ij}(t=0,s)`$. Note that the scattering amplitude is complex, and that this approach gives both the real and imaginary part of the dispersion law modification, also known as the optical potential.

There are two major theoretical approaches to the issue discussed in literature, to be called an $`schannel`$ and a $`tchannel`$ one. The former approach assumes that the scattering amplitude is dominated by s-channel resonances which are known to decay into the $`i+j`$ channel. For most mesons such as $`\pi ,\omega ,\rho ,K`$ in a gas made of pions such calculation has been made e.g. in  related with such $`\pi \rho `$ resonances as $`a_1`$ or $`N\rho `$ resonance $`N^{}(1520)`$  for $`\rho `$ at SPS. Note that the signs of the effects are opposite in those examples, as seen from the following table:

The majority of the particles in the matter are Goldstone bosons $`\pi ,K,\eta `$ which do not interact at small momenta. However attraction between other particles is there. Using a simple expression for the mass shift one gets $`\delta m_\rho ^N28MeV`$ due to all $`\overline{B}+B`$. An additional shift $`\delta m_\rho ^v10MeV`$ comes from scalar exchanges between $`\rho `$ and all other vector mesons $`\rho ,\omega ,K^{}`$. The main difference between the two mechanisms of the mass shift discussed above is that the t-channel attraction is not associated with the broadening, while the s-channel resonances increase the width by about 50 MeV. On the other hand, there is a “kinematic” effect working to the opposite direction. The negative mass shift discussed above automatically reduces the width, both because of the reduced phase space and also due to the power of p in the P-wave matrix element. The magnitude of this effect for the predicted mass shift is

$`\delta \mathrm{\Gamma }_\rho =3{\displaystyle \frac{\delta m_\rho }{m_\rho }}\mathrm{\Gamma }_\rho 50MeV`$ (5)
So, inside the accuracy these two effects cancel each other.

The invariant mass distributions in pp and mid-central AuAu of the $`\pi ^+\pi ^{}`$ system, with a transverse momentum cut $`0.2<p_t<0.9GeV`$ have been measured by STAR , see also the C.Markert’s talk here. I would not have time here to discuss the shift in pp (see ), and I only comment that in AuAu the $`\rho `$ peak is found to be shifted by additional $`40MeV`$ in mass, but the width is the same. This agrees well with estimates above.

The same approach should of course be applied to many other resonances. For more narrow resonances, like $`K^{}`$, we expect smaller shifts, while for wider resonances like $`\sigma `$ we predicted a complete change of shape. At small freezeout $`T100MeV`$ sigma was predicted to be deformed into a much more narrow structure at mass of about 400 MeV, see figures in . Exactly such a peak has been seen by STAR111I have seen it first at this workshop, on the day after my talk, shown by Gary Westfall in a plot of the balance function as a function of $`q_{inv}`$ with a peak he said “nobody understands”. It was rather good timing between the prediction and its experimental confirmation.. Another confirmation of very late freezeout and low $`T`$, the $`\sigma /\rho `$ ratio strongly grows toward central collisions. We are waiting for quantitative analysis of these data with great interest.

## 5 Fragment coalescence

The issue of coalescence, such as $`p+nd`$, was discussed in many papers over the years, and authors struggled with the question how to calculate its rate. In particular, it is clear that when the level crosses zero the wave function at the origin vanishes, and so the production should do so too. And if the production rate is small compared to two other relevant rates, $`\nu _{abs}`$ and $`dlogN/dt`$, there is never thermal equilibrium and one should not use statistical models.

Significant progress has been made in recent paper by Ioffe et al  who have pointed out how to use consistently the in-matter widths of all particles and obtain the production rate. We are now incorporating it into the picture of expansion and the optical depths discussed above, and hope to get quantitative results for $`\overline{d},d`$ spectra soon.

## Acknowledgements

My first thanks go to the organizers, especially to the hard core from MSU, who wonderfully managed to collect us together year after year in the right places. I should also thank my collaborators, G.Brown and P.Kolb, and also Patricia Fachini who made analysis of the STAR data for $`\pi \pi `$ resonances and explained it to us. Finally, my research is supported by the DOE grant No. DE-FG02-88ER40388.


text:
# Travelling waves in pipe flow

## Abstract

A family of three-dimensional travelling waves for flow through a pipe of circular cross section is identified. The travelling waves are dominated by pairs of downstream vortices and streaks. They originate in saddle-node bifurcations at Reynolds numbers as low as $`1250`$. All states are immediately unstable. Their dynamical significance is that they provide a skeleton for the formation of a chaotic saddle that can explain the intermittent transition to turbulence and the sensitive dependence on initial conditions in this shear flow.

Based on decades of studies it is consensus that Hagen-Poiseuille flow through a pipe of circular cross section belongs to the class of shear flows that does not become linearly unstable. Nevertheless, it undergoes an intermittent transition to turbulence for sufficiently high Reynolds numbers and sufficiently large initial perturbations, as first documented in the classic experiments by Reynolds reynolds1883 . Since then many studies have analyzed the mechanisms of transition and the properties of the turbulent state wygnanski73 ; eliahou98 ; ma99 ; boberg88 ; brosa99 . Particularly relevant to the present analysis are the investigations by Darbyshire and Mullin darbyshire95  which clearly show a strong sensitivity to perturbations and a broad intermittent range of decaying and turbulent perturbations in an amplitude vs. Reynolds number plane. Experimental and numerical studies eliahou98 ; schmid01  show that the turbulent flow in the transition region is dominated by downstream vortices and streaks. Various models for their dynamics have been analyzed trefethen93 ; holmes96 . Taking the full nonlinearity into account Waleffe developed the concept of a nonlinear turbulence cycle for the regeneration of vortices and streaks hamilton95 . In addition, stationary solutions and travelling waves have been found in the full nonlinear equations for plane Couette, Taylor-Couette and plane Poiseuille flow cherhabili96 ; busse97 ; waleffe98 ; FE00 ; waleffe01 . It has been suggested that these structures provide a skeleton for the transition to turbulence and the observed intermittency bruno97 ; bruno02 . They clearly dominate various observables in low Reynolds number turbulent flows holmes96 , and are also relevant for an understanding of the effects of non-Newtonian additives waleffe02 .

The existence of exact coherent states in pipe flow has been an object of speculation for some time boberg88 ; darbyshire95 ; brosa99 . As we will show here Hagen-Poiseuille flow supports families of travelling waves with structures similar to those observed in other shear flows as well. This underlines the significance of vortex-streak interactions also in this system and opens alternative routes to modelling and controlling pipe flow.

The existence of stationary states in plane Couette flow is connected with an inversion symmetry in the laminar profile. In the absence of such a symmetry in pipe flow the simplest states we can expect are travelling waves (TWs), i.e. coherent structures that move with constant wave speed. The wave speed depends on shape and structure and is not known in advance. We search for TWs using a Newton-Raphson method applied to the full Navier-Stokes equations for an incompressible fluid,

$`_t𝐮+(𝐮)𝐮`$ $`=`$ $`\nu \mathrm{\Delta }𝐮p`$
$`𝐮`$ $`=`$ $`0.`$
The kinematic viscosity $`\nu `$ is kept as a free parameter for the continuation method and the Reynolds number is determined by the mean flow velocity $`U`$, i.e. $`Re=2RU/\nu `$, where $`R`$ is the pipe radius. For the numerical representation of the Navier-Stokes equation we use a Fourier-Legendre collocation method in cylindrical coordinates with Lagrange multipliers to account for no-slip boundary conditions at the wall and the constraints that the flow field is solenoidal, analytical and regular in the center for $`r=0`$. The program was verified by reproducing literature values for the linearized problem schmid94 , for the nonlinear dynamics of optimal modes zikanov96  and for the statistical properties of fully developed turbulent flow up to Reynolds numbers of $`5000`$ eggels94 . The Newton-Raphson method was implemented with a spatial resolution of $`|n/11|+|m/11|<1`$, where $`n`$ and $`m`$ are the azimuthal and downstream wavenumbers, respectively, and 56 Legendre polynomials radially. This gives us about $`5600`$ dynamically active degrees of freedom.

Initial conditions for the Newton method are obtained by embedding pipe flow into a larger, parameter dependent family of flows. This has been successful in other cases busse97 ; FE00 ; waleffe01 . For pipe flow Barnes and Kerswell barnes00  showed that the instabilities that arise from rotation do not extend to the case of pure pipe flow. Guided by the dominance of downstream vortices in the stationary and travelling states in planar shear flows we therefore adopted the following strategy, similar to the one used in waleffe98 : at low Reynolds number a volume force $`𝐅`$ that generates a transverse flow with translation invariant downstream vortices was added. While there is no specific physical realization underlying this volume force, it could be achieved with differential heating and cooling as predicted by the symmetries of the states, or by magnetic forcing in an electrically conducting fluid brosa91 . The amplitude of the driving was then increased beyond the first bifurcation, where a stable 3-d TW without translational symmetry appears. A path following scheme was then used to track this forced state through parameter space: the Reynolds number was increased and the amplitude of the body force decreased. If a situation with zero body force could be reached, a 3-d TW for the original pipe flow has been found.

In this way we have identified four different states that have a discrete azimuthal rotation symmetry by construction. An $`n`$-fold symmetry $`C_n`$ is defined by the invariance under rotation around the pipe axis by an angle $`2\pi /n`$, i.e.

$$(u_r,u_\varphi ,u_z)(r,\varphi ,z)(u_r,u_\varphi ,u_z)(r,\varphi +2\pi /n,z).$$

The resolution in azimuthal direction was adjusted to the fundamental domains of angle $`2\pi /n`$ to avoid carrying vanishing amplitudes.

The state that bifurcates at the lowest Reynolds number has a three-fold rotational symmetry in azimuthal direction (Fig. 1).

The dynamics of the states is governed by vortices which have a predominant downstream orientation and are slightly tilted inside the volume. They resemble near wall coherent structures in turbulent wall flows holmes96 . They transport slow fluid towards the center and fast fluid towards the wall, thus producing high-speed streaks near the wall and low-speed streaks towards the center (all speeds are relative to the laminar profile for that Reynolds number). The steeper gradients near the wall imply higher friction losses and higher pressure gradients for the same flow speed. The position of the fast-speed streaks hardly changes in downstream direction, whereas the low-speed streaks change shape and move in position quite a bit.

In addition to the $`C_n`$ symmetry externally imposed by the forcing the TW’s dynamically develop a ‘shift-and-reflect’ symmetry, in agreement with the streak instability of other model flows waleffe95B : if $`L_z`$ is the wavelength of the state, then

$$(u_r,u_\varphi ,u_z)(r,\varphi ,z)(u_r,u_\varphi ,u_z)(r,\varphi ,z+L_z/2).$$
(1)
The same shift-and-reflect symmetry is present in 3-d coherent states in plane Couette, Taylor-Couette and plane Poiseuille flow busse97 ; FE00 ; waleffe98 ; waleffe01 .

Each TW is part of a two-parameter family of states that exists for a range of Reynolds numbers and a range of downstream wavelengths (Fig. 2).

The lowest critical Reynolds number is obtained for a wavelength of about $`4.2R`$ for the twofold symmetry and about $`2.5R`$ for the three- and fourfold symmetric state. This period seems a bit shorter than that for optimal states in plane Couette flow where the optimal wavelength is about $`2\pi d`$ with $`d`$ the gap width. However, when the wedge shape of the boundary is taken into account and the lengths are compared to the widths, the comparison is more favorable: the ratio widths to lengths is about 2$`\pi `$:4$`\pi `$ = 1:2 in plane Couette flow and about 1:2.5 in pipe flow.

Other states that we could identify in pipe flow include those with twofold, fourfold and fivefold symmetry (Fig. 3).

All these states have a similar topology: they have $`2n`$ streaks of fast fluid close to the wall and streaks of slow fluid towards the center. The high-speed streaks near the wall remain fairly stationary over one period of the wave, and the low speed streaks in the center oscillate vigorously. All states have the shift-and-reflect symmetry (1). The wavelength that gives the lowest critical Reynolds number for these states varies with the symmetry (Fig. 2). Several properties of these states are listed in Table 1. The lowest critical Reynolds number is obtained for the TW with $`C_3`$ symmetry, for which the arrangement of vortices is optimal in the sense of being closest to a hexagonal packing, the preferred pattern in other systems. The states with $`C_4`$ and $`C_5`$ symmetry are not included in the table as they seem to be the ones which are most sensitive to the numerical resolution and not fully converged. A preliminary estimate for their critical Reynolds number is about $`1590`$ for $`C_4`$ and $`2600`$ for $`C_5`$.

Interestingly, the list of TWs does not include a state with the $`C_1`$\- symmetry of the mode that shows the strongest linear transient growth bergstroem93 ; schmid94 . A state with $`C_1`$ symmetry could be found for non-vanishing volume force, but up to Reynolds numbers of $`4000`$ it could not be continued to pure pipe flow without reconnecting to the laminar profile. This should have implications for low dimensional model building bergstroem99 ; brosa99  where linear arguments have led to an emphasis of the $`C_1`$ mode.

The phase speed $`c`$ of the TWs is a function of both parameters, length and Reynolds number. In all cases the phase speed is larger than the mean speed but slower than the maximal speed possible with a laminar fluid (Fig. 4). That is to say, the wave still propagates downstream when viewed from a frame of reference moving with the mean flow velocity.

By advecting the unperturbed profile small transverse components can produce strong downstream streaks. Table 1 shows that for the TWs the transverse components are about an order of magnitude smaller than the differences between the downstream components of the laminar profile and the TW.

The TWs appear as finite amplitude saddle-node bifurcations. Away from the bifurcations there is an upper (nodal) and a lower (saddle) branch. The saddle has one more positive eigenvalue than the node, but in all cases we found that even at the point of bifurcation at least one eigenvalue is positive: all states are thus unstable right from the point of creation. This is in contrast to plane Couette flow, where at least the lowest state has an interval of stability, although a tiny one busse97 .

The friction properties of all states together with experimental and numerical results for turbulent flows are shown in Fig. 5. The friction factor $`f`$ is defined as

$$f=4R\mathrm{\Delta }p/\rho U^2L_z,$$
(2)
with $`\mathrm{\Delta }p`$ the pressure drop over a distance $`L_z`$ along the pipe and $`\rho `$ the density of the fluid. The saddle-node bifurcations take place at Reynolds numbers significantly below the values where typical perturbations induce turbulent dynamics. At the bifurcation the friction factor is higher than the value extrapolated from turbulent states at higher $`Re`$. The upper branches in Fig. 5 correspond to the lower branches in Fig. 4: states with higher friction have lower phase velocity and vice versa.

Besides the TWs presented here we can expect many other TWs. On the assumption that narrower vortices will require stronger driving we anticipate that they will have critical Reynolds number higher than the $`1250`$ found for the $`C_3`$ state. In addition, there are all the states that bifurcate from the TWs, as well as states that bifurcate in secondary bifurcations. Furthermore, there will be homoclinic and heteroclinic connections between these TWs as well as other periodic and aperiodic states. All these states together then form a chaotic saddle in phase space that supports the turbulent dynamics tel . One characteristic of this chaotic saddle is an exponential distribution of life times and this is indeed found in numerical simulations pipe\_trans . It will be interesting to see to which extent these states can be used to approximate properties of the turbulent flow.

Support by the Deutsche Forschungsgemeinschaft is greatfully acknowledged.