text:
# On the Local Structure of Topological Charge Fluctuations in QCD

## Abstract

We consider the lattice topological charge density introduced by Hasenfratz, Laliena and Niedermayer and propose its eigenmode expansion as a tool to investigate the structure of topological charge fluctuations in QCD. The resulting effective density is built from local chiralities studied previously. At every order of the expansion the density exactly sums up to the global topological charge, and the leading term describes the maximally smooth space-time distribution of charge relevant for propagating light fermions. We use this framework to demonstrate our previous suggestion that the bulk of topological charge in QCD does not effectively appear in the form of quantized unit lumps. Our conclusion implies that it is unlikely that the mixing of “would-be” zeromodes associated with such lumps is the prevalent microscopic mechanism for spontaneous chiral symmetry breaking in QCD. We also present first results quantitatively characterizing the space-time behavior of effective densities. For coherent fluctuations contained in spherical regions we find a continuous distribution of associated charges essentially ending at $`\mathrm{\hspace{0.17em}0.5}`$.

Some intriguing effects in QCD, such as the large $`\eta ^{}`$ mass and the $`\theta `$-dependence, are related to vacuum fluctuations of topological charge WittenUA(1) . Understanding the local structure of topological charge fluctuations is thus of great interest for building a detailed picture of how these phenomena arise in terms of fundamental degrees of freedom. Another important phenomenon possibly related to topological charge fluctuations is spontaneous chiral symmetry breaking (S$`\chi `$SB). This is based on the proposition that the effective low-energy structure of topological charge fluctuations in QCD is such that in a typical configuration most of topological charge is concentrated in $`N_L`$ typically non-overlapping, four-dimensional space-time regions $`_i`$, each containing a sign-coherent lump of approximately unit topological charge (i.e. the generalization of the instanton liquid picture). If true, this would imply an appealing microscopic explanation for the origin of Dirac near-zero modes ins\_mixing  and hence the origin of S$`\chi `$SB Banks\_Casher . The connection to fermions arises by associating with each lump $`_i`$ a localized chiral mode $`\chi ^i`$ which “would be” a zeromode in the absence of other lumps, and a formation of a topological subspace of low-lying eigenmodes $`\psi ^i_ja_{ij}\chi ^j,i=1,\mathrm{}N_L`$. Verifying this “topological mixing” scenario has long been hindered by inherent difficulties in interpreting the local behavior of topological charge density for typical configurations of lattice-regularized theory. We have argued Hor01A ; Hor02A  that meaningful information can be extracted indirectly by studying the local chirality of low-lying modes. The flat behavior of the associated $`X`$-distribution observed in the initial study suggested that topological mixing was not the origin of the near-zeromodes, and also indicated the lack of extended self-dual excitations Hor01A . However, later studies revealed the double-peaked behavior Followup  which is qualitatively consistent with both these aspects, but not sufficient for their demonstration Hor02A . Being agnostic about self-duality, we observed in Ref. Hor02A  that even with the double-peak structure present, the patterns in low eigenmodes of the overlap operator are inconsistent with the presence of a distinctive topological subspace, and hence with the topological mixing scenario. In this work, we propose a framework generalizing the local chirality method, and use it to demonstrate this conclusion directly. The basis of our approach is the topological charge density $`q_x=\frac{1}{2}\text{tr}\gamma _5D_{x,x}`$ Has98A  associated with $`\gamma _5`$-Hermitian $`D`$ satisfying $`\{D,\gamma _5\}=D\gamma _5D`$ (Ginsparg-Wilson (GW) fermions). Note that since $`\mathrm{tr}\gamma _5=0`$, one can also write $`q_x=\mathrm{tr}\gamma _5(C\frac{1}{2}D_{x,x})`$ with arbitrary constant $`C`$. We choose $`C=1`$ (see also Ref. Fuji99 ) which guarantees that the eigenmode expansion of $`q_x`$ is properly normalized and satisfies the index theorem for arbitrary truncation (see below). In this case we have

$$q_x=\text{tr}\gamma _5(1\frac{1}{2}D_{x,x})=\underset{\lambda }{}(1\frac{\lambda }{2})c_x^\lambda $$
(1)
where $`c_x^\lambda =\psi _x^{\lambda +}\gamma _5\psi _x^\lambda `$ is the local chirality of the mode with eigenvalue $`\lambda `$. For low-lying truncation $`(|\lambda |0)`$ this is effectively the sum of individual local chiralities. Very recently, Gattringer offered the evidence for a large degree of self-duality in the vicinity of peaks of an eigenmode-filtered action density Gat02A . If the extended nature of self-dual regions is confirmed, then combined with our conclusions one is led to a picture of inhomogeneous, but essentially continuous topological charge fluctuations at low energy, with a significant level of self-duality in the vicinity of local maxima. The origin of the double-peaked behavior of local chirality would then be the local attraction of spinorial components by (anti)self-dual field as argued in Hor01A , with topological mixing not playing a role. In the last part of this work, we describe the first attempt to characterize the typical amounts of topological charge associated with coherent fluctuations, as well as some interesting quantitative results that emerged.

1. In the usual discussion of topological mixing there is a paradox rooted in the negativity of Euclidean topological charge density correlator at nonzero distance,

$$<q(x)q(0)>\mathrm{\hspace{0.33em}\hspace{0.17em}0},|x|>0.$$
(2)
This follows from reflection positivity and the fact that the operator $`q(x)`$ is reflection-odd Kov ; SeSt . One consequence is that it is impossible that in typical configurations most of topological charge is concentrated in coherent four-dimensional lumps of finite physical size (e.g. instantons). Indeed, if this was the case and the typical size of such lumps was $`r_c`$, then the average correlator would be positive over $`r_c`$, thus contradicting (2). This means that lumps $`_i`$ cannot be identified in $`q(x)`$.

These considerations do not rule out topological mixing as an effective low-energy scenario for S$`\chi `$SB, but emphasize the necessity of some short-distance filtering when studying these issues. This can be achieved by replacing the local operator $`q(x)`$ by a nonlocal one. While the use of nonlocal operators is generally unacceptable, here they will only serve as a filter to smear the singular contact term in (2) into a finite positive core. The dominance of coherent structures is then possible. Following up on our reasoning in Refs. Hor01A ; Hor02A , we propose that the low energy truncation of the Dirac eigenmode expansion for $`q(x)`$ (see below) can serve as the physically motivated nonlocal operator needed. We stress the physical motivation because, for the problem at hand, the underlying issue is whether the light fermion effectively feels a collection of coherent unit lumps as it moves through the vacuum, or there is some other dynamics governing its propagation. This is to be decided by the fermion.

The desired fermion filtering can be realized starting from lattice-regularized theory. Among the lattice topological charge density operators considered, the ones associated with GW fermionic kernels are perhaps theoretically most appealing Has98A . These are constrained so that they sum up to integer global charge defined by counting the exact zeromodes of the GW kernel used. The continuum gauge-fermion correspondence (index theorem) for smooth gauge fields is thus extended to the lattice by construction Has98A . Before defining the filtered densities starting from (1), it is useful to recall that the spectrum of $`D`$ contains $`N_0`$ zeromodes with global chirality $`+1`$ ($`N_0^+`$) or $`1`$ ($`N_0^{}`$), $`N_2`$ modes at $`\lambda =2`$ with global chirality $`+1`$ ($`N_2^+`$) or $`1`$ ($`N_2^{}`$), and $`N_p`$ pairs of complex modes with zero global chirality. The dimension of Dirac space is $`N=N_0+N_2+2N_p`$, and the topological charge is $`Q=_xq_x=N_0^{}N_0^+=N_2^+N_2^{}`$. We now associate with $`q_x`$ the set of related densities

$$q_x^{(k)}=\underset{i=1}{\overset{N_0}{}}c_x^{0,i}\underset{j=1}{\overset{k}{}}(2\text{Re}\lambda _j)c_x^{\lambda _j}$$
(3)
representing the truncated eigenmode sum that includes all zeromodes and the $`k`$ lowest-lying complex pairs. Complex eigenvalues $`\lambda _j`$ enter the sum ordered on the upper (lower) half of the spectral circle. Eigenmode-filtered densities $`q^{(k)}`$ are real and $`_xq_x^{(k)}=Q=N_0^{}N_0^+`$ thus leading to identical global fluctuations satisfying the index theorem. We have $`q_x^{(N_p)}=q_x`$ and since $`q^{(0)}`$ can be identically zero for $`Q=0`$ configurations, we consider $`q^{(1)}`$ to be the leading order in the expansion. The infrared eigenmodes of $`D`$ are significantly smoother than the underlying gauge field Hor02A , and $`q^{(1)}`$ is expected to be maximally smooth. As $`k`$ increases, $`q^{(k)}`$ gradually incorporates more short-distance structure. This does not mean that there is a strict new cutoff present in filtered densities. However, ultraviolet fluctuations irrelevant for propagation of light quarks are filtered out.

We now use the overlap Dirac operator Neu98BA  (see details in Ref. Hor02A ) on Wilson gauge backgrounds (see Table 1) to demonstrate the basic properties of filtered densities discussed above. We have calculated the full density $`q`$ for configuration $`𝒞_2`$ from ensemble $`_4`$. In Fig. 1 we show its correlator $`C_q(r)`$ (normalized at the origin) and compare it to $`C_{q^{(1)}}`$ and $`C_{q^{(16)}}`$ for filtered densities. Note that $`C_q`$ has a very short-ranged positive core as expected from (2). This provides indirect evidence that the locality of the non-ultralocal operator $`q`$ is quite good. As for the correlators $`C_{q^{(k)}}`$, their range is evidently larger. The shape stabilizes at about $`k=6`$ and changes very slowly from then on. To further characterize the roughness of filtered densities we calculate

$$(G^{(k)})^2\underset{x,\mu }{}(q_{x+\mu }^{(k)}q_x^{(k)})^2$$
(4)
and plot $`G^{(k)}`$ as a function of $`k`$. As expected, $`q^{(k)}`$ becomes rougher as $`k`$ increases. While the physical size of the positive core in $`C_q`$ is expected to go to zero in the continuum limit, this is not necessarily so for $`C_{q^{(k)}}`$. To see that, we have calculated $`C_{q^{(2)}}`$ for ensembles $`_1`$-$`_3`$. The size of the positive core in the average correlator was determined as $`<r>`$ over the probability distribution given by $`C_{q^{(2)}}(r)`$ in the range from zero up to the maximal distance where $`C_{q^{(2)}}(r)`$ is manifestly positive (with errors taken into account). As can be seen in Fig. 1, the size scales well, indicating that if the dominance of coherent four-dimensional structures in $`q^{(2)}`$ can be established, these structures can survive the continuum limit.

2. Introduction of fermion-filtered densities allows us to study the relevance of topological mixing consistently. If the topological subspace of low-lying modes $`\psi ^i_ja_{ij}\chi ^j,i=1,\mathrm{},N_L`$ is formed, then

$$\underset{x_i}{}\underset{j=1}{\overset{N_L}{}}\psi _x^{j+}\gamma _5\psi _x^js_i,\gamma _5\chi ^i=s_i\chi ^i,s_i=\pm 1.$$
(5)
This implies that the structure of quantized unit lumps $`_i`$ will be revealed in $`q^{(k)}`$ when all modes belonging to the topological subspace are included. Indeed, from equation (5) we have

$$\underset{x_i}{}q_x^{(k)}Q_i\pm 1,2kN_LN_0$$
(6)
since $`\text{Re}\lambda 0`$ for modes in topological subspace. Moreover, the value of pure gauge topological susceptibility ($`1\text{fm}^4`$) constrains $`N_L`$ in lump-dominated configurations of volume V fm4 to be $`N_LV`$. This implies that the lumpy structure should typically be saturated in $`q^{(k)}`$ with $`k\stackrel{<}{}V/2`$. Consequently, $`q^{(2)}`$ is expected to be well sufficient for all ensembles in Table 1.

We now ask whether the subset $`_i_i`$ of the lattice containing most of topological charge can be identified using $`q^{(k)}`$. A simple way to proceed is to order lattice points by the magnitude of $`q_x^{(k)}`$ and compute the running sum of positive ($`Q^{+,k}(f)`$) and negative ($`Q^{,k}(f)`$) charge as the fraction $`f`$ of the highest points included increases. Functions $`Q^{+,k}`$ ($`Q^{,k}`$) should stabilize to a constant ($``$ integer-valued) plateau at a well-defined value of $`f`$ corresponding to fraction of volume occupied by $`_i_i`$. In Fig. 2 we show the behavior of $`Q^{+,k}`$ for configuration $`𝒞_2`$ ($`Q=1`$) from $`_5`$. We find no sign of plateaus for any $`k`$, and nothing special happens around $`k=2`$. Such behavior is characteristic for all configurations from the ensembles in Table 1. The same conclusion applies to functions $`Q^k(f)`$ monitoring total charge. The smooth monotonic behavior exhibited in Fig. 2 excludes the possibility that the bulk of topological charge is effectively concentrated in a small subvolume $`_i_i`$ of typically isolated lumps.

For another quantitative test, consider a configuration with topological charge $`Q`$ and assume it is dominated by $`|Q|+\xi `$ unit lumps (antilumps) and $`\xi `$ antilumps (lumps). If the number of zeromodes is minimal, i.e. $`|Q|=N_0`$ (true for all our configurations), then the dimension of the topological subspace is $`N_L=|Q|+2\xi `$ and we should have $`_x|q_x^{(\xi )}|Q^{\xi ,abs}N_L=|Q|+2\xi `$. In fact, the quantity $`|Q|+2kQ^{k,abs}`$ would be close to zero for all $`k\xi `$ (since each mode in the topological subspace contributes approximately unity to $`Q^{k,abs}`$), and it would start to increase rapidly for $`k>\xi `$. In a theory where the topological mixing scenario is relevant, monitoring the $`k`$-dependence of $`|Q|+2kQ^{k,abs}`$ could thus serve as a procedure for determining the dimension of the topological subspace for a given configuration. We have computed this $`k`$-dependence for all our configurations. The minimal value is always achieved at $`k=1`$ and is of order one rather than close to zero. Moreover, in ensembles $`_4`$, $`_5`$, where a wide range of $`k`$ is available, we observe monotonic (linear) increase for every configuration. This robust behavior is illustrated in Fig. 3, confirming again that there are no signs of a distinctive topological subspace and hence no signs of dominance by unit lumps.

3. The above arguments do not depend on specific properties of $`_i`$, such as their shape, volume or particular field content. The results show that the topological charge is effectively carried by the bulk of the lattice, which is inconsistent with dominance by unit quantized lumps. At the same time, we find inhomogeneous behavior with noticeable peaks in $`q^{(k)}`$, accompanied by visible coherence (see $`C_{q^{(k)}}`$ in Fig. 1). As a first step toward understanding this structure we now provide a simple characteristic of the typical values of charge associated with such coherent fluctuations. Given an arbitrary density $`q`$ we consider the sets $`^l`$ of centers of coherent behavior with lattice resolution $`\sqrt{l}`$, namely

$$^l\{x:|q_y|<|q_x|,q_xq_y>0y;\mathrm{\hspace{0.33em}0}<|xy|^2l\}.$$

The elements of $`^l`$ are local maxima of $`|q|`$ over distance $`\sqrt{l}`$, for which the sign of $`q`$ is coherent over at least the same distance. Obviously, $`^1^2^3\mathrm{}`$. To every $`x^l`$ we assign a radius $`R_x`$ defined as the maximal distance from $`x`$ over which the density is still coherent, i.e. $`R_x\sqrt{l}`$. The charge $`Q_x`$ is also assigned by summing the density over the sphere of radius $`R_x`$ centered at $`x`$. However, to be able to interpret $`Q_x`$ as a charge corresponding to an individual fluctuation, there should typically be no other centers within $`R_x`$. This allows for fixing the resolution in a self-consistent way by choosing the smallest $`l`$ with this condition satisfied.

To carry out the above procedure meaningfully for $`q^{(k)}`$, one should work in the range of $`k`$ where coherent fluctuations and their charges are relatively stable with respect to a change of $`k`$. For our ensembles, this is best satisfied for $`_4`$ where $`q^{(9)}`$ is available. In physical terms, this corresponds to including eigenmodes with imaginary parts up to $``$ 500 MeV. The fraction of overlapping centers drops dramatically (to 5 %) at $`l=3`$, which is the value we have used. The corresponding distribution of charges for non-overlapping centers is shown in Fig. 4. A property which is insensitive to the choice of $`k`$ and $`l`$ is that the distribution effectively ends at about $`0.5`$. This appears to hold also in the case of full density $`q_x`$ (with only two configurations available). Interestingly, this behavior might be compatible with the presence of center vortices in the QCD vacuum as the recent discussion of topology in the field of an idealized vortex suggests Rein01A  (The possible manifestation of center vortices in topological charge fluctuations was also recently discussed in Ref. Cor02 .).

To summarize, we have proposed that the low-energy behavior of topological charge density can be studied by using a suitable non-local realization of $`q(x)`$ as a filter for short-distance fluctuations. If one is interested in the aspects of topological charge affecting the low-momentum propagation of light quarks, then the appropriate non-local realization is naturally available through the low-eigenmode expansion of $`q(x)`$. Such fermion filtering can be explicitly realized starting from lattice-regularized theory. We have proposed the expansion of $`q_x`$ associated with GW fermions as an ideal tool for this purpose. This provided us with a consistent framework to test whether the propagation of light fermion is effectively driven by the dominance of unit topological lumps in the QCD vacuum. We find that this is not the case, implying that the picture of S$`\chi `$SB based on the mixing of corresponding topological “would-be” zeromodes is not accurate. We emphasize that we have not ruled out the logical possibility that some unit quantized structures with well-defined boundaries can occur. However, we have shown that the bulk of topological charge does not come in that form. A first attempt to characterize the inhomogeneous nature of topological charge fluctuations in filtered densities resulted in very interesting results indicating that patterns of local behavior can provide us with detailed information about dynamically important structures in the QCD vacuum.

###### Acknowledgements.

The support by U.S. Department of Energy grants DE-FG05-84ER40154, DE-FG02-95ER40907, DE-FG02-97ER41027 and DE-AC05-84ER40150 is acknowledged.


format: text

text:
# Dynamical properties of constrained drops

## I Introduction

Since the advent of accelerators powerful enough to explore the behavior of nuclear systems at intermediate energies, a new field of research has been opened, the thermodynamics of small systems. This problem has emerged since the first suggestion that a highly excited nuclear system might be undergoing a second order phase transition. Such a hypotesis was brought up by the pioneering work of the Purdue group in which a power law was fit to the mass spectra in collisions of highly energetic protons against heavy nuclei. Since then many experiments have been performed in this energy range (for a recent review see ). In order to gain understanding of the physical phenomena involved in such a process different models have been devised. This models can be roughly classified into two main groups, the statistical ones and the dynamical ones. In the first case it is assumed that the highly excited nuclear system resulting from the collision of two heavy ions is able to equilibrate in a fixed volume (usually referred to as the freeze-out volume) and then fragments. This fragmentation process is entirely driven by the available phase space. This model has been quite successful in describing some aspects of the fragmentation phenomena. In what is relevant for our present analysis we recall that the caloric curve ($`CC`$), the functional relationship between the energy and the temperature of the system, predicted in this approach is of the ”rise-plateau-rise” type and consequently displays a vapor branch. The second approach, the dynamical one, assumes that the system is properly described by a given Hamiltonian and puts no constrain in the evolution of the excited system. In this case no equilibration has been found. It has been shown  that the resulting $`CC`$ is of the ”rise-plateau” type, and no vapor branch is present. In this case the $`CC`$ is defined as the functional relationship between the temperature and the energy of the system at fragmentation time.

This two views of the process of fragmentation give interesting insight on the properties of a fragmenting system, but differ in the relevant issue of the degree of equilibration of the system. It is then of primary importance to have a clear understanding of the thermodynamical and dynamical properties of finite systems enclosed in constraining volumes, thus having the possibility of reaching equilibrium. Some studies have been performed in this direction. For example in  the behavior of lattice gases constrained in fluctuating volumes has been analyzed. Moreover in  a rather extensive analysis of this kind of systems was performed. In  a phase diagram was built.

In this paper we focus on the analysis of drops formed by 147 Lennard Jones particles enclosed in different volumes performing a study in terms of the amount of energy added to the system and its density.

The main reason for choosing this particular system is that in a series of previous works we have performed a detailed, though still incomplete, study of its properties when no constrain is imposed . Moreover we have recently given a first step towards the characterization of the behavior of equilibrated systems .The sizes of the constrained volumes have been chosen such that we go from a dilute case to a very dense one. In order to study such a system we performed extensive numerical simulations of the Molecular Dynamics type.

This paper is structured in the following way: In section I we briefly describe the model used in our numerical simulations.

In section II we show the calculated Caloric Curves ($`CC`$) for our test system composed of 147 Lennard Jones particles as a function of its energy for different values of the radius of the constraining volume.

In section III we describe the methodology we use to calculate the MLE and the $`d_{\mathrm{}}`$.

In section IV we show our results. Finally conclusions are drawn.

### A The Model

Following a series of previous works in this field  we will rely heavily on numerical simulations of classical systems interacting via a Lennard Jones potential, which reads:

$$V(r)=\{\begin{array}{cc}4ϵ\left[\left(\frac{\sigma }{r}\right)^{12}\left(\frac{\sigma }{r}\right)^6\left(\frac{\sigma }{r_c}\right)^{12}+\left(\frac{\sigma }{r_c}\right)^6\right]\hfill & r<r_c\hfill \\ 0\hfill & rr_c\hfill \end{array}$$
(1)
We fix the cut-off radius as $`r_c=3\sigma `$. Energy and distance are measured in units of the potential well ($`ϵ`$) and the distance at which the potential changes sign ($`\sigma `$), respectively. The unit of time used is: $`t_0=\sqrt{\sigma ^2m/48ϵ}`$. In our numerical experiments initial conditions were constructed using the already presented  method of cutting spherical drops composed of 147 particles out of equilibrated, periodic, 512 particles per cell L.J. system. We choose this kind of initialization because we consider that the resulting correlations present in the system at initial time are the least biased ones. It is worth mentioning at this point that the general features of the fragmentation process do not depend on the initial state as has been shown in  were the fragmentation of bidimensional drops via the collision with fast aggregates of three particles was studied. A broad energy range was considered such that the asymptotic mass spectra of the fragmented drops displays from ”U shaped” pattern , to exponentially decaying one. Somewhere in between this two extremes a power law like spectra can be found.

Although our system is purely classical and no direct connection with nuclear systems can be established, one has to take into account that the main features of the nuclear interaction (strongly repulsive at very short range and attractive at a longer range) are present in this interaction potential. Then it is quite plausible that the main features of L.J. systems should appear in nuclear systems. In this respect, and of great importance for this work, both systems present an equation of state of the same type.

Because we will be mainly interested in the behavior of constrained systems we have to define the walls that will contain the excited drop. We have defined the walls of our container using a very strongly short ranged repulsive potential (cut and shifted) defined as:

$$V_w=\{\begin{array}{ccc}\mathrm{exp}\left(1/(rR)\right)\mathrm{exp}\left(1/(r_cR)\right)\hfill & \hfill & r_crR\hfill \\ 0\hfill & \hfill & r<r_c\hfill \end{array}$$
(2)
Where $`r_c`$ defines the skin of the constraining volume.

### B Caloric Curves at constant density of constrained Systems

The Caloric Curve ($`CC`$) is one of the main observables in the analysis of multifragmentation. In fact there is no agreement in the nuclear physics community about its properties. The analysis of experimental data has given different results. Different Thermometers have been used and the resulting $`CC`$s are of two types. On the one hand we have those that present a ”rise-plateau-rise” shape which resemble the standard view, inherited from classical thermodynamics, and have induced to recognize a transition from a liquid-like to a vapor-like state. This same result has been obtained when one adheres to statistical models to analyze the phenomena. Because statistical models impose equilibrations in a given ”freeze-out” volume it is natural to get this kind of behavior. On the other hand classical molecular dynamics calculations indicate  that when the process of fragmentation is properly analyzed in phase space  the resulting caloric curve, in this case defined as the temperature of the system at fragmentation time, is of the type ”rise-plateau” and no vapor branch is present. This behavior has been traced to the presence of a collective motion, expansion, that behaves as a ”heat sink”. Recent experimental results seem to confirm this view .

In a recent work we have shown that the effect of confining the excited drop in a given volume, i.e. allowing the system to reach equilibrium, is the appearance of the vapor branch. In this work we will extend those calculations exploring a broader range of densities and incorporing useful dynamical quantities.

From the analysis of molecular dynamics calculations of constrained systems  we have extracted the $`CC`$s that are displayed in fig.1). It is immediate that a very interesting phenomenon takes place. The same system at the same temperature will change its behavior as a function of the density. In this figure it can be clearly seen that for the less dense case a clear loop in the $`CC`$ is obtained. But as the density is increased this loop disappears and is replaced by a change in the slope. Moreover at even higher densities the $`CC`$ looks essentially straight and all signals of a change in behavior are erased.

The origin of the relationship between the $`CC`$ and the density can be understood quite easily by analyzing the biggest cluster formed in the system. For this purpose we use two algorithms already presented in the literature  (see Appendix for details). Very briefly we can say that the, improperly called, Minimum Spanning Tree (MST) algorithm looks for clusters of interacting particles in configuration space, and completely disregards the relative momentum. On the other hand we use the Early Cluster Recognition Algorithm (ECRA) which seeks for the most bound partition in phase space. This algorithm has allowed us to find that the fragmentation in an expanding system takes place very early in the evolution. In fig.2) we show the obtained biggest fragment from the analysis of configurations according to both techniques in constrained systems (see caption for details).

The emerging picture is the following: In the case that the fragment recognition algorithm is the MST, we see that for low densities there is enough room in the constraining volume to allow the formation of drops i.e. At low densities the biggest MST fragment is a decreasing function of the energy. But as we increase the density the biggest MST fragment comprises most of the mass in the system. On the other hand ECRA analysis shows that for all the densities considered the ECRA biggest fragment is a decreasing function of the energy.

Once the $`CC`$ is known it is easy to calculate the Entropy $`S`$ as a function of the energy and the density.

$$S=\frac{dE}{T(E)}$$
(3)
In fig.3) we show $`S`$ for a dilute system (see caption for details). It is immediate that a convex intruder appears which has been proposed to be a signature of a first order phase transition in non extensive systems. (i.e. the formation of surfaces turns the entropy into a non-extensive function in small systems).

The next step is to calculate the behavior of the Thermal response function of such a system.

$$TRF=\left(\frac{dE}{dT}\right)$$
(4)
In fig.4) the results of such calculation are displayed. We can see that for low densities, as a consequence of the presence of a loop, two poles and negative values are attained by this quantity. This has been signaled as an evidence for a first order phase transition. It is due to the fact that surfaces appear in the system. As the density is increased we reach a threshold above which the $`CC`$ only displays a change in the slope. Then the two poles merge in a single finite maximum (finite size effect). We can then state that the system goes from a first order like to a second order like behavior as a function of the available space in configuration space.

### C Maximum Lyapunov exponent and Asymptotic distance in momentum space

We now proceed with our analysis and focus on the dynamical aspects of the constrained system and its relationship with the above described thermodynamical properties.

One of the main tools to study a chaotic system is the maximum Lyapunov exponent , which is a measure of the sensitivity of the system to initial conditions and also gives an idea of the velocity at which the system explores the available phase space. Given two very close initial conditions in phase space, the MLE, $`\widehat{\lambda }`$, is given by the following relation :

$$\lambda =\underset{t\mathrm{}}{lim}\underset{d(0)0}{lim}\left[\frac{1}{t}ln\frac{d(t)}{d(0)}\right].$$
(5)
where $`d(t)`$ is the distance in phase space between two trajectories ($`1`$ and $`2`$) which initially differs each other in a very small quantity $`d_0`$.

In order to calculate this quantity we must define a metric:

$$d_{12}(t)=\sqrt{\underset{i=1,N}{}\left[a\left(r_1\left(t\right)r_2\left(t\right)\right)^2+b\left(p_1\left(t\right)p_2\left(t\right)\right)^2\right]}$$
(6)
In this equation $`a`$ and $`b`$ are constants which take care of the units. It has been shown that the $`MLE`$ are independent of the metric . In our case we have found it useful to take $`a=0`$ and $`b=1/m.`$ With m the mass of the particles that according to the units defined in the model gives $`b=1.`$

If we calculate the distance in momentum space between nearby trajectories we find an exponential growth, followed by a saturation. This saturation is due to the fact that the available phase space is limited (The energy is a constant of motion). In order to handle this we calculate the $`MLE`$ following a method due to Benettin . In this method, after a time step $`\tau <<\tau _{sat}`$, the distance $`d(\tau )=d_1`$ is re-scaled to $`d_0`$ in the maximum growing direction and the quantity $`ln[d_1/d_0]`$ is saved. Repeating the procedure at every time step $`\tau `$, the logarithmic increments $`ln[\frac{d_i}{d_{i1}}]`$ are collected. The MLE is defined as:

$$\lambda =\underset{n\mathrm{}}{lim}\frac{1}{n\tau }\underset{i=1}{\overset{n}{}}ln\left|\frac{d_i}{d_{i1}}\right|$$
(7)
The ratio $`\frac{d_i}{d_{i1}}`$ is a measure of the exponential divergence between two initially nearby orbits along the maximum growth direction at time $`i\tau `$.

Another quantity that has been recently proposed is the $`d_{\mathrm{}}`$ . This is a measure of the maximum distance in momentum space when two initially close trajectories are followed in time. In order to calculate $`d_{\mathrm{}}`$ we use the same metric but no rescaling as in the Benettin method is used. In this case the molecular dynamics evolutions needed to calculate such a quantity are much shorter due to the fact that $`d_{\mathrm{}}`$ reaches a plateau rather fast.

## II Results

In what follows we show the results of our calculations.

In first place we show the $`MLE`$ for constrained systems. In fig.5) we show this quantity for four densities (the same densities at which the $`CCs`$ are shown in fig.1)). The following features are relevant: For energies below 0, the $`MLE`$ is an increasing function of the energy, following the behavior of the $`CC`$. As energy is increased the behavior of the $`MLEs`$ changes abruptly. In the range $`0E1`$ we find that, for the low density case the $`MLE`$ displays a very pronounced loop which is in correspondence with the loop displayed by the $`CC`$. On the other hand for the next two densities a clear loop is present in the $`MLE`$ while the corresponding $`CCs`$ only show a change in the slope in this region. Finally for the highest density considered in this work the $`CC`$ is featureless in this region while the $`MLE`$ shows a valley.

In order to gain insight into the reasons of this behavior we have found it useful to study the following quantities. In first place we look at the mass distributions. When dealing with $`MST`$ spectra we see that two very different behaviors appear. For $`R15\sigma `$ fragments (drops) of different sizes appear in the volume as energy is added to the system. The mass spectra goes from extreme U shape in the liquid like region to exponentially decaying in the vapor like region. It is clear that as the density is raised there will be less and less room for the system to form non interacting drops. It is then seen that for $`R8\sigma `$ only one $`MST`$ cluster appears in the system. At this point the $`CC`$ also shows a change in behavior going from a curve that displays a loop to one that only shows a change in slope. On the other hand when the same configurations are analyzed using the $`ECFM`$ , via its algorithmic implementation : the $`ECRA`$ method, a different picture appears: Regardless the size of the constraining walls, the system behaves as undergoing a phase transition, i.e. drops are formed in phase space.

We have also calculated the Normalized Variance of the size of the Maximum fragment ($`NVM`$). It has been shown  that this magnitude will display a maximum for a system whose mass spectra is well described by a scaling law of the type $`n_ss^\tau f(z)`$ at $`f(z)=1`$. In other words at the point in which the mass distribution is a power law, i.e. a distribution function free of scales. $`NVM`$ is defined as

$$NVM=\frac{\sigma ^2\left(BF\right)}{<BF>}$$
(8)
Where $`BF`$ is the normalized mass of the biggest fragment and the brackets indicate an ensemble averaging. In fig.6) we show the $`NVM`$ for the four densities considered. It can be readily seen that this quantity displays a maximum in the energy range in which the $`MLE`$ displays a loop (or a valley for the highest density considered).

Now we consider the behavior of the other relevant magnitude presented in this work, i.e. $`d_{\mathrm{}}.`$ In order to explore its scaling properties we define $`d_{\mathrm{}}^n`$ which stands for the normalized asymptotic distance

$$d_{\mathrm{}}^n(t)=\left(\frac{_{i=1}^N\left[\left(p_1(t)p_2(t)\right)^2\right]_i}{\frac{1}{2}_{i=1}^N\left[p_1(t)\right]_i^2}\right)^{\frac{1}{2}}$$
(9)
In other words

$$d_{\mathrm{}}^n(t)=\frac{d(t)}{\sqrt{K}}$$
(10)
i.e., the distance in momentum space normalized to the square root of the kinetic energy.

In fig.7) we show the result of such an analysis for $`R=8\sigma `$ and in fig.8) the same but at $`R=15\sigma `$ . In both cases it is seen that the asymptotic distance in momentum space scales as the square root of the kinetic energy (and then as the square root of the temperature) of the system as it was conjectured in .

Finally in fig.9) we show the $`d_{\mathrm{}}`$ as a function of the energy for two densities. It can be seen that its behavior clearly resembles the corresponding one for the $`CC`$, thus this quantity could be used instead of T if the latter is not reliably measurable.

## III Conclusions

The main results of our calculations can be summarized in the following:

a) The $`CC`$ is strongly dependent on the density of our system. As the density is increased the behavior of this function goes from displaying a loop (which we refer to as first order like behavior) to one that only presents a change in the slope (in this case we talk about a second order like behavior). In the first case the $`TRF`$ displays two poles, which converge into a single maximum as the density surpasses a given threshold $`(R8.7\sigma \rho 0.05\sigma ^3).`$ Three regions can then be recognized, The liquid like, the transition region and the vapor like region.

b) The $`MLE`$ is quite sensitive to the transition from liquid like to vapor like states of the finite constrained system. In all cases, even when the $`CC`$ is almost featureless a signal is detected in the $`MLE.`$ This signal has been found to be significant by studying other quantities like the $`NVM.`$

c)The $`d_{\mathrm{}}`$ has been shown to scale like $`\sqrt{K}`$ and then as the $`CC`$, giving useful insight into the amount of Phase Space visited by the system.

## IV Acknowledgment

This work was partially supported by the University of Buenos Aires (UBA) via grant No tw98, and CONICET via grant No 4436/96.

C.O.Dorso is a member of the carrera del investigador (CONICET), P.Balenzuela is a fellow of the Conicet, M.Ison is a fellow of UBA.

## V Appendix

In previous papers the main fragment recognition algorithms currently in use have been fully analyzed . The simplest definition of cluster is basically: a group of particles that are close to each other and far away from the rest. The fragment recognition method known as minimum spanning tree (MST) is based on the last idea (I). In this approach a cluster is defined in the following way: given a set of particles $`i,j,k,\mathrm{}`$, they belong to a cluster $`C`$ if :

$$iC,jC/\left|𝐫_i𝐫_j\right|R_{cl}$$
(11)
where $`𝐫_i`$ and $`𝐫_j`$ denote the positions of the particles and $`R_{cl}`$ is a parameter usually referred to as clusterization radius, and is usually related to the range of the interaction potential. In our calculations we took $`R_{cl}=3\sigma `$.

On the other hand, the early cluster formation model (ECFM) , is based on the next definition: clusters are those that define the most bound partition of the system, i.e. the partition (defined by the set of clusters $`\{C_i\}`$) that minimizes the sum of the energies of each fragment according to:

$$E_{\left\{C_i\right\}}=\underset{i}{}\left[\underset{jC_i}{}K_j^{cm}+\underset{j,kC_i}{}V_{j,k}\right]$$
(12)
where the first sum is over the clusters of the partition, and $`K_j^{cm}`$ is the kinetic energy of particle $`j`$ measured in the center of mass frame of the cluster which contains particle $`j`$. The algorithm (early cluster recognition algorithm, ECRA) devised to achieve this goal is based on an optimization procedure in the spirit of simulated annealing .


text:
# Oscillations of the superconducting critical current in Nb-Cu-Ni-Cu-Nb junctions

## Abstract

We report on experimental studies of superconductor-ferromagnet layered structures. Strong oscillations of the critical supercurrent were observed with the thickness variation of the ferromagnet. Using known microscopic parameters of Ni, we found reasonable agreement between the period of oscillations and the decay of the measured critical current, and theoretical calculations.

The interplay between superconductivity and ferromagnetism is an old subject which was studied extensively over decades . The most striking effect in such systems is the formation of the so called $`\pi `$ phase junction in a superconductor-ferromagnet-superconductor (SFS) structure . One of the manifestations of the $`\pi `$ phase is a non-monotonic variation of the critical temperature , $`T_c`$, with the variation of the ferromagnetic layer thickness, $`d`$. Oscillatory dependence of $`T_c`$ vs. $`d`$ of SF multilayers has been observed by few groups , however, other groups  reported monotonic behavior in similar structures. Interpretations of these experiments suggested the $`\pi `$ junction mechanism as well as some other origins of the effect . The recent observation  of non-monotonic behavior of the critical current as a function of temperature in weak-ferromagnetic layer of CuxNi1-x between two Nb layers is considered as an unambiguous proof of the $`\pi `$ phase formation. Another interesting theoretical prediction  concerns non-monotonic behavior of the critical current as a function of the thickness of the ferromagnetic layer $`d`$. According to the above predictions, the critical current $`I_c`$ is expected to oscillate and decay as $`d`$ is increased. To the best of our knowledge, such a behavior has not been reported so far.

In this paper, we present the experimental evidence of oscillatory behavior of the critical current vs. thickness variation of ferromagnetic Ni layer. We also show a reasonable agreement between our data and the theories  in the appropriate limit of $`E_{ex}>>\mathrm{}/\tau >>k_BT_c`$. Here, $`\tau `$ and $`E_{ex}`$ are the electron relaxation time and the exchange energy of the ferromagnet, and $`T_c`$ is the critical temperature of the superconductor.

We have studied temperature and thickness dependence of the critical current in Nb-Cu-Ni-Cu-Nb junctions. Additionally, we have studied the thickness dependence of critical current of similar junctions without a Ni layer. The junctions with $`10`$x$`10\mu `$m2 area were fabricated with the standard photolithography technique. The process contained three stages of lithography: liftoff of the bottom Nb-Cu layer, liftoff of the variable thickness Ni or Cu layers, and liftoff of the top Cu-Nb layer. Nb films were sputtered using a magnetron gun and in situ covered with the Cu layer by thermal evaporation, for preventing the Nb oxidation. The ferromagnet layers of Ni were e-gun evaporated in a separate vacuum chamber, and subsequently covered in situ by Cu. It is important to emphasize that all samples were prepared simultaneously. The variation of Ni thickness for Nb-Cu-Ni-Cu-Nb junction and variation of Cu for Nb-Cu-Nb junctions was obtained by a specially designed shutter, which exposed the samples in sequence, so that every sample was exposed to the evaporating Ni or Cu for additional fragments of time. This method guaranteed that all the interfaces between each layer in our multilayer structure are identical, and the only difference between the samples is their Ni or Cu thickness. The thickness of each Nb layer was 2000Å. The total thickness of the Cu in the Nb-Cu-Ni-Cu-Nb junctions was 2400Å and the Ni thickness varied from 10Å to 90Å.

In the other set of junctions, namely Nb-Cu-Nb, the Cu thickness varied from 5000Å to 10000Å. The structure of these junctions is shown schematically in the inset of Fig.1.

Typical resistance of both types of junctions was around $`100\mu \mathrm{\Omega }`$, independent of both Ni and Cu thicknesses. The critical temperature of the Nb was about 8.5K. The measured resistivity of Ni layer $`\rho =31\mu \mathrm{\Omega }cm`$, together with the known value $`\rho l=1.510^{11}\mathrm{\Omega }cm^2`$, allowed us to establish the mean free path $`l=48`$Å.

The measurements were performed in 4He cryostat in the range from 6.5K down to 1.5K. The critical current was measured by passing a DC current with a small AC modulation through the sample. The AC voltage which appeared above the critical DC current was picked up by lock-in amplifier operated in a transformer mode.

A typical differential resistance of both types of junctions vs. DC current is plotted in Fig.1, which clearly indicates on the sharp onset of dissipation at $`I=I_c`$. First, we show results of the Nb-Cu-Nb set (Fig.2). As expected , the critical current $`I_c`$ normalized to the length of the junction $`L`$, decreases exponentially with $`L`$. From the decay of the critical current we found the thermal length of the the Cu, $`L_T=\sqrt{\mathrm{}D/2\pi k_BT}=1600`$Å at 4.2K, and the diffusion constant $`D=860cm^2/s`$.

Very different thickness dependence was found in the Nb-Cu-Ni-Cu-Nb junctions. Fig.3 shows the thickness dependence of the critical current in these junctions at $`T=4.2K`$. In spite of the large error bars, the non-monotonic variation of the critical current is quite evident, namely, the deviations of the data from the exponential decay surmounts by far the uncertainty of each measured point.

A further, and even a stronger evidence for the oscillatory behavior is provided by temperature dependence of the critical current. Fig.4 shows a family of $`I_c`$ vs. $`T`$ curves for all Ni thicknesses, which are normalized to their values at $`4.2K`$. We define the slope $`\alpha _T`$ of temperature variation, namely $`\alpha _Td(lnI_c(T))/dT`$, and plot these values as a function of $`d`$ in Fig.3 (squares). Oscillations of $`\alpha _T`$ are very prominent, and are in anti-phase with the oscillations of the critical current. Unlike the critical current, which had quite large experimental error bars, the slope of the temperature dependence had error bars of only few percents. Such a behavior of $`\alpha _T`$ as a function of $`d`$ is intimately related to the variation of the critical current oscillations amplitude with temperature. Although the amplitude of the oscillations decreases at high temperatures, their relative value increases, and on a semi-logarithmic plot of Fig.5, the oscillatory behavior of $`I_c`$ is more pronounced at higher temperatures.

We would like to start our discussion with what is well known for superconductor-normal-superconductor (SNS) junctions. The dependence of critical current of such junctions on the thickness of N layer is exponential, $`I_c\alpha e^{L/L_T}`$ for $`L>L_T`$. In the opposite limit, where $`L<L_T`$, $`I_c`$ approaches a constant value. This behavior has been observed in few experiments for different material systems , including our sample presented in Fig.2.

Non-monotonic behavior of the critical current is not expected theoretically, and had never been observed experimentally in SNS structures. Since the thickness of Ni is the only parameter which is varied in our Nb-Cu-Ni-Cu-Nb samples, it is reasonable to assume that the oscillations observed in the data are due to the presence of Ni in the structure.

As mentioned in the introduction, oscillatory behavior of the critical current vs. the thickness of the ferromagnetic layer is predicted theoretically . The origin of these oscillations is the phase shift acquired by electron-hole Andreev particles upon entrance into the ferromagnet, due to their different spin orientations. Several expressions have been derived for the critical current in various limits of the strength of $`E_{ex}`$, thickness of the ferromagnet $`d`$ and disorder. Since in our experiment we have determined only the magnitude of the critical current $`I_c`$, the formulae below will be written for the absolute value $`|I_c|`$. For the clean, thin and strong ferromagnetic layers $`l>\mathrm{}v_f/E_{ex},d`$, the critical current should vary with the thickness as:

$$I_c|\mathrm{sin}(2E_{ex}d/\mathrm{}v_f)|/(2E_{ex}d/\mathrm{}v_f)$$
(1)
The dependence of absolute value of the critical current is shown on Fig.6 by dashed line. In the opposite limit, namely $`l<d`$, the following dependence was derived (dotted line in Fig.6):

$$I_cye^y|\mathrm{sin}(y+\pi /4)|y=\sqrt{\mathrm{}D/4E_{ex}}$$
(2)
Another expression for the critical current in the limit $`l>\mathrm{}v_f/E_{ex}`$, is given by:

$$I_c\left|Re\underset{\omega _n>0}{}\frac{\mathrm{\Delta }^2}{\mathrm{\Delta }^2+\omega _n^2}_1^1\frac{\mu d\mu }{\mathrm{sinh}(k_\omega d/\mu l)}\right|$$
(3)
where $`\omega _n=\pi Tk_B(2n+1)`$ is the Matsubara frequency, $`n`$ is an integer number, $`k_\omega =(1+2|\omega _n|\tau /\mathrm{})2iE_{ex}\tau /\mathrm{}`$, $`\mu =cos\theta `$, $`\theta `$ is the angle between the momentum and the normal to the SF interface, and $`\mathrm{\Delta }`$ is the order parameter in the superconductor. The solid line of Fig. 6 represents the above expression. In fitting our data to the expressions for all limits, we have used $`v_f=2.8`$x$`10^5m/sec`$ and $`l=48`$Å (based on the resistivity measurements of Ni). Therefore, the only fitting parameter apart from the numerical prefactor was the strength of the exchange interaction $`E_{ex}`$.

The periodicity of oscillations in Eq. 1, $`L_{osc}\pi \mathrm{}v_f/E_{ex}54`$Å fits the best our data when $`E_{ex}=107\pm 3meV`$. Note, that $`L_{osc}`$ is twice larger than the periodicity observed in Fig.6, due to the absolute value taken for $`I_c`$. In order to have similar periodicity in Eq. 2, we have used weaker exchange field $`E_{ex}=86meV`$. All values are close to the recently reported value  $`E_{ex}=115meV`$.

The curve given by Eq. 3 in Fig.6 follows the decay and the phase of the oscillations of the measured $`I_c`$ closer than the other two curves for large $`d`$. Furthermore, we estimate $`\mathrm{}v_f/E_{ex}=17`$Å which is smaller than the measured mean free path $`l=48`$Å, and therefore, the validity of using Eq. 3, $`l>\mathrm{}v_f/E_{ex}`$, is justified. However, since Eq. 3 is valid only for $`d>l`$, the data points should follow Eq. 1 for $`d<40`$Å. We thus conclude that our data is consistent with Eq. 3 and Eq. 1 in their appropriate limits. Therefor, we give the fit of Eq. 1 and Eq. 3 to our data for four temperatures in Fig. 7.

The oscillations of $`\alpha _T`$ (Fig. 3), which slightly smear the oscillations of $`I_c`$ (Fig. 5), are not accounted by any of the mentioned theories and have only a minor effect on the overall quality of the fit presented in Fig. 7. A possible origin of this effect could be related to the presence of normal metal on both sides of the ferromagnetic layer. The structure we have used contains layers of Cu of a combined thickness $`L=2400`$Å, which is much larger than the thickness of the Ni layer. It is obvious that at high temperatures ($`L_T<L`$) the contribution of Andreev particles propagating via short diffusion path is predominant, since paths which are much longer than $`L`$ do not contribute to the supercurrent. At low temperatures, where $`L_T>L`$, the contribution of longer paths is not negligible. The latter implies that the impact of the trajectories crossing Ni layer several times is higher at low temperatures. The multiple crossing of the ferromagnetic layer could produce additional harmonics of the basic frequency of the oscillations, which will smear the non-monotonic behavior. This multiple crossing has not been considered in the existing theories for S-F-S junctions, and therefore might be responsible for the observed temperature dependence of the oscillations in S-N-F-N-S junctions (Fig. 4).

In summary, we have observed the oscillations and the decay of the critical current in Nb-Cu-Ni-Cu-Nb junctions upon the increase of the thickness of the Ni layer. We found a reasonable agreement with the recent theoretical calculations in the appropriate limit. The puzzling temperature dependence of the oscillations’ amplitude remains unexplained.

We would like to thank A. Aharony, K. B. Efetov, O. Entin-Wohlman, V. Fleurov, Y. Imry, K. Kikoin, Z. Ovadiyahu and A. Schiller for fruitful discussions. Partial support by the Israel Science Foundation founded by the Israel Academy of Sciences and Humanities — centers of Excellence Program is gratefully acknowledged.