text:
# LU TP 01-13LUNFD6/(NFFL-7194)2001 Hadronization corrections to helicity components of the fragmentation function

## 1 Introduction

Studies of fragmentation functions have always been important, since these distributions can not be predicted theoretically, but only be measured experimentally, and consecutively be described by phenomenological models. Hadronic decays of $`\gamma ^{}/Z^0`$ provide a particularly convenient set of events for analysis and interpretation. Helicity components of the fragmentation function, measured in such events, can be used in various QCD studies, e.g., the extraction of the gluon fragmentation function, and the evaluation of $`\alpha _s`$. However, existing theoretical calculation being restricted to the perturbative region, hadronization corrections must be taken into account. In what follows, methods for applying such corrections using the Jetset event generator  as an example, will be discussed.

Consider the angular distribution in the process $`e^+e^{}\gamma ^{}/Z^0q\overline{q}`$ in its rest frame. Assuming that the final quark and antiquark are not charge-tagged, i.e. that the forward–backward asymmetry is not accessed, the cross section can be written as

$$\frac{\mathrm{d}\sigma }{\mathrm{d}(\mathrm{cos}\theta )}=\frac{3}{8}(1+\mathrm{cos}^2\theta )\sigma _\mathrm{T}+\frac{3}{4}\mathrm{sin}^2\theta \sigma _\mathrm{L}.$$
(1)
Here $`\sigma _\mathrm{T}`$ ($`\sigma _\mathrm{L}`$) is the cross section associated with a transverse (longitudinal) gauge boson polarization state with respect to the $`q\overline{q}`$ axis, and $`\theta `$ is the polar angle of a particle with respect to the incoming lepton axis. To lowest order, only mass effects contribute to a non-vanishing $`\sigma _\mathrm{L}`$, but only for the vector part of the cross section, and there only with a coefficient $`\sigma _\mathrm{L}/\sigma _\mathrm{T}=2m_q^2/E_{\mathrm{cm}}^2`$. Even for the $`b`$ quarks this gives a negligible $`\sigma _\mathrm{L}`$ contribution at the energies around the $`Z^0`$ peak. Therefore, $`\sigma _\mathrm{L}`$ effectively starts in $`𝒪(\alpha _\mathrm{s})`$ of perturbation theory, associated with the emission of gluons.

Since partons are not directly observable, one may define a hadron-level analogue of Eq.(1,

$$\frac{\mathrm{d}^2\sigma ^\mathrm{h}}{\mathrm{d}x\mathrm{d}(\mathrm{cos}\theta )}=\frac{3}{8}(1+\mathrm{cos}^2\theta )\frac{\mathrm{d}\sigma _\mathrm{T}}{\mathrm{d}x}+\frac{3}{4}\mathrm{sin}^2\theta \frac{\mathrm{d}\sigma _\mathrm{L}}{\mathrm{d}x}.$$
(2)
Here $`x`$ would preferably be associated with the energy fraction taken by a hadron, $`x_E=2E/E_{\mathrm{cm}}`$, so that $`x_E=2`$ in each event. Experimentally it is more convenient to use the momentum fraction $`x_p`$. The transverse and longitudinal fragmentation functions are defined by a normalization to the total cross section $`\sigma _{\mathrm{tot}}=\sigma _\mathrm{T}+\sigma _\mathrm{L}`$ ,

$$F_\mathrm{T}(x)=\frac{1}{\sigma _{\mathrm{tot}}}\frac{\mathrm{d}\sigma _\mathrm{T}}{\mathrm{d}x},F_\mathrm{L}(x)=\frac{1}{\sigma _{\mathrm{tot}}}\frac{\mathrm{d}\sigma _\mathrm{L}}{\mathrm{d}x}.$$
(3)
The former is dominated by the fragmentation of quark jets, whereas the latter receives a major contribution from gluon fragmentation. Therefore an experimental determination of $`F_\mathrm{L}(x)`$ is a first step towards an extraction of the gluon fragmentation function, alternative to what is offered by more direct methods in three-jet events . Several experimental $`F_\mathrm{L}(x)`$ studies have also been presented .

A complication is that hadrons are not moving in the direction of their imagined mother parton. Already in lowest order of perturbation theory, for $`q\overline{q}`$ two-jet events, nonperturbative hadronization gives an effective $`p_{}`$ smearing that induces a nonvanishing $`F_\mathrm{L}(x)`$ even where none is expected. Furthermore, the association of a hadron to a single mother parton is not in agreement with our current best understanding of the hadronization process, where it is rather the colour field between a colour-connected pair of partons (a string piece , or a cluster ) that mediates the hadron production. Therefore the structure of smearing effects may become rather nontrivial. Obviously, the effects are especially important at small $`x`$, which is also the region where one would hope to have some sensitivity to the gluon fragmentation function. The string picture also casts in doubt the concept of a gluon fragmentation function defined from inclusive quantities, since the string fragmentation of a parton depends on the angles to other colour-connected partons.

The issue of hadronization corrections to fragmentation functions was addressed in . The emphasis was on the $`\sigma _\mathrm{L}`$ that can be extracted from $`_0^1F_\mathrm{L}(x)xdx`$ rather than on $`F_\mathrm{L}(x)`$ itself, however. Therefore we here address how hadronization affects $`F_\mathrm{L}(x)`$ (and $`F_\mathrm{T}(x)`$). One main conclusion is that a simple smearing approach is not sufficient to describe hadronization effects. Thus it appears impossible to define a completely model-independent, hadronization-smearing-corrected $`F_\mathrm{L}(x)`$, that could be used to extract a gluon fragmentation function. We further suggest a correction procedure, based on a cluster search strategy, that should give a less model dependent $`F_\mathrm{L}(x)`$, but at the price of introducing the cluster resolution scale $`y`$ as a new parameter in the problem.

## 2 The simple smearing

The string model description of $`q\overline{q}`$ events introduces a Gaussian transverse momentum smearing of primary hadrons, $`\mathrm{exp}(p_{}^2/2\sigma ^2)\mathrm{d}^2p_{}`$, where $`\sigma 0.36`$ GeV . Many primary hadrons are unstable and decay further; this distorts the original Gaussian spectrum and reduces the average $`p_{}`$. Since decay products have smaller $`p_{}`$, momentum parallel to the jet axis, some correlation is also introduced between $`p_{}`$ and $`p_{}`$. Therefore no simple parameterization is proposed, but instead a Monte Carlo simulation with Jetset 7.4  is used to histogram the amount of angular smearing for different $`x_p`$ bins.

It is now assumed that this smearing should be applied both to quark and gluon jets, so that a hadron will not move exactly in the direction of its mother parton. There are obvious shortcomings to equating different kinds of jets, like that gluon jets have a lower energy and do not contain decays of charm and bottom hadrons, but those particular issues only introduce moderate corrections. More severe objections can be raised to the association of hadrons to individual partons, as we will discuss further in the next section, but forget for the moment.

Had particles not been smeared in $`p_{}`$, but parallel with their parton of origin, the shape of the angular distribution

$$F(x_p,\mathrm{cos}\theta )=\frac{3}{8}(1+\mathrm{cos}^2\theta )F_\mathrm{T}(x_p)+\frac{3}{4}\mathrm{sin}^2\theta F_\mathrm{L}(x_p)$$
(4)
in a bin of $`x_p`$ could be used to extract $`F_\mathrm{T}(x_p)`$ and $`F_\mathrm{L}(x_p)`$ in that bin. The abovementioned smearing will now modify this. The two angular shapes, $`(3/8)(1+\mathrm{cos}^2\theta )`$ and $`(3/4)\mathrm{sin}^2\theta `$, both normalized to unity, are therefore convoluted with the $`x_p`$-dependent smearing distributions, characterized by a distribution in the smearing angle $`\theta ^{\mathrm{sm}}`$ and an isotropic azimuthal distribution $`\phi ^{\mathrm{sm}}`$. That is, a parton at an angle $`\theta ^\mathrm{p}`$ will produce a hadron at an angle $`\theta ^\mathrm{h}`$, where

$$\mathrm{cos}\theta ^\mathrm{h}=\mathrm{cos}\theta ^\mathrm{p}\mathrm{cos}\theta ^{\mathrm{sm}}\mathrm{sin}\theta ^\mathrm{p}\mathrm{sin}\theta ^{\mathrm{sm}}\mathrm{cos}\phi ^{\mathrm{sm}}.$$
(5)
Data can now be fitted both to the “non-smeared” angular distribution form, Eq.(4), and to the convoluted (“smeared”) analogue. The resulting distributions for $`F_T`$ and $`F_L`$, obtained using the Jetset generated events, are shown in Fig. 1. The effect on $`F_T`$ and $`F_L`$ is visible for hadron momenta below 10% of the beam energy. The low momentum region is affected the most, giving $`F_L`$ values reduced up to one order of magnitude.

## 3 Objections to the simple smearing

The above smearing procedure is correct to lowest order in $`\alpha _\mathrm{s}`$, i.e. it describes how two-jet events can induce a nonvanishing $`F_\mathrm{L}(x_p)`$. We know, however, that hadronization of three-jet events cannot be described in terms of a simple incoherent sum of three $`q`$, $`\overline{q}`$ and $`g`$ jets. One example is the string/drag effect , i.e. that particle production is suppressed in the angular region between the $`q`$ and $`\overline{q}`$ and enhanced in the other two regions, well confirmed experimentally . High-momentum hadrons still essentially follow the separate parton directions, but low-momentum ones are significantly affected. These are the ones where the angular smearing effects are large to begin with. It is well-known that the string effect leads to more two-jetlike events, e.g. in terms of thrust $`T`$, than implied by symmetric smearing . The reason is to be found in the enhanced production of particles between two colour-connected partons that are close in angle, leading to them seemingly being even closer, e.g. that the opening angle between the reconstructed jets typically is smaller than that between the original partons.

There is a nontrivial topology dependence on string effects, especially when multiple gluon emission is considered. The issue is therefore best studied in an event generator, by comparing angular distributions on the parton ($`i=`$p below) and on the hadron ($`i=`$h) level. As a simple measure of the jettiness of events we use $`\overline{T}=1T`$, defined on the parton level of each event. The integrated $`x_E`$ spectrum is decomposed as

$`\left({\displaystyle \frac{\mathrm{d}\sigma _{\mathrm{tot}}}{\mathrm{d}\overline{T}}}\right)^1{\displaystyle _0^1}{\displaystyle \frac{\mathrm{d}^3\sigma ^i}{\mathrm{d}\overline{T}\mathrm{d}x_E\mathrm{d}(\mathrm{cos}\theta )}}x_Edx_E`$ $`=`$
$`={\displaystyle \frac{3}{8}}(1+\mathrm{cos}^2\theta )F_\mathrm{T}^i(\overline{T})+{\displaystyle \frac{3}{4}}\mathrm{sin}^2\theta F_\mathrm{L}^i(\overline{T})`$ $`,`$ (6)
where the $`x_E`$-weighting ensures a common normalization $`F_\mathrm{T}^i(\overline{T})+F_\mathrm{L}^i(\overline{T})=2`$ at parton and hadron level (values of $`F_\mathrm{T}^i`$ and $`F_\mathrm{L}^i`$ most conveniently are obtained by weighting each particle with an appropriate angular factor ). Then $`\mathrm{\Delta }F_\mathrm{L}(\overline{T})=F_\mathrm{L}^\mathrm{h}(\overline{T})F_\mathrm{L}^\mathrm{p}(\overline{T})`$ is a simple measure of the hadronization impact on $`F_\mathrm{L}`$. This quantity is shown in Fig. 2, for one realistic simulation and two toy ones, for $`u\overline{u}`$ events at 91.2 GeV. In the realistic case, a parton shower is used to generate multiparton configurations, followed by string fragmentation. The shower develops down to a cut-off scale $`Q_01`$ GeV, so that also events in the first bin, $`\overline{T}<0.01`$, can contain some gluons. The other two histograms are based on $`𝒪(\alpha _\mathrm{s})`$ matrix elements, where only 2- and 3-parton configurations are generated, with a cut $`\overline{T}>0.01`$ on the latter to avoid the singularities of the 3-parton matrix element. Thus the first bin here represents pure 2-parton events. While one simulation is again based on string fragmentation, the other assumes isotropic smearing around the jet axes, basically the independent fragmentation scheme of Hoyer et al. . (The same fragmentation parameters, tuned to the shower model, have been used in all three cases. A retuning of parameters for the $`𝒪(\alpha _\mathrm{s})`$ simulations would have given a larger nonperturbative $`p_{}`$ width $`\sigma `$ to cover for the lack of perturbative gluons, and so would have implied even larger fragmentation smearing.)

The isotropic smearing is, as expected, giving a rather constant hadronization correction $`\mathrm{\Delta }F_\mathrm{L}(\overline{T})`$. There is some jump up in going from two to three jets that are smeared, followed by a slow but steady drop with $`\overline{T}`$, since the longitudinal component itself is increasing in importance with $`\overline{T}`$ and therefore gives an increasing hadronization smearing of the longitudinal component onto the transverse one rather than only the other way around. By contrast, the string fragmentation provides a much steeper drop of $`\mathrm{\Delta }F_\mathrm{L}`$ with $`\overline{T}`$, kicking in immediately when going from two to three partons, and enhanced in the shower simulation relative to the simpler $`𝒪(\alpha _\mathrm{s})`$ one. At large $`\overline{T}`$ the overall hadronization correction can even turn negative. Averaging over the $`\overline{T}`$ spectrum (with mean value $`\overline{T}0.05`$), we conclude that the typical hadronization smearing contribution is only about a third of the naively expected one, as obtained from two-parton results. (Qualitatively this agrees with and probably explains a similar observation in  of smaller-than-expected hadronization corrections when using herwig .) That is, if hadronization corrections are viewed as a power series in $`\alpha _\mathrm{s}`$, the $`𝒪(\alpha _\mathrm{s})`$ term is of opposite sign and almost as big as the $`𝒪(1)`$ one.

It should be remembered, however, that this is integrated over all $`x_E`$, and that we have no similar way of addressing results in specific $`x`$ bins, since the parton and hadron $`x`$ spectra are quite different. Thus it is likely that the $`F_\mathrm{L}(x_p)`$ derived in the previous section is an underestimation, just like an $`F_\mathrm{L}(x_p)`$ found without any smearing corrections is likely to be an overestimation, but it appears impossible to find the “correct” $`F_\mathrm{L}(x_p)`$ without making detailed assumptions about the hadronization process.

## 4 Clustering

Given the problems with the above smearing recipe, we introduce a new strategy, based on the clustering approach. In a nutshell, we propose to rotate all hadrons to the direction of the cluster they belong to, as an approximate way of removing hadronization smearing effects. Only thereafter is $`F_\mathrm{L}(x)`$ extracted from this modified $`\mathrm{cos}\theta `$ distribution. The strategy is explained further in the following.

In clustering algorithms, nearby hadrons are combined to form clusters/jets, in a way that should reflect the underlying partonic state, to some approximation. The combination process is controlled by (at least) one separation parameter, call it $`y_{\mathrm{min}}`$, such that the final state contains no pair of clusters closer to each other than that. Clustering algorithms can be applied also to a partonic state, and here $`y_{\mathrm{min}}`$ provides a regularization of soft and collinear divergences in the perturbative cross sections. It is then meaningful to calculate the distribution of partons at a factorization scale $`\mu ^2=y_{\mathrm{min}}E_{\mathrm{cm}}^2`$, and define scale-dependent fragmentation functions parameterizing the subsequent soft-perturbative and nonperturbative hadronization. The latter should obey standard QCD evolution equations, starting from some unknown nonperturbative form at a low reference scale.

Over the years many cluster algorithms have been proposed , each with its strengths and weaknesses. In this article we adopt the Durham one , which is a standard for many perturbative calculations. The distance measure between two clusters $`i`$ and $`j`$ is

$$y_{ij}=\frac{2\mathrm{min}(E_i^2,E_j^2)(1\mathrm{cos}\theta _{ij})}{E_{\mathrm{vis}}^2},$$
(7)
so that $`\sqrt{y_{ij}}`$ roughly corresponds to the relative transverse momentum, scaled to the total visible energy $`E_{\mathrm{vis}}`$ ($`=E_{\mathrm{cm}}`$ for an ideal detector).

If we begin by considering a simple $`q\overline{q}`$ event, it should reconstruct back to two clusters, unless $`y_{\mathrm{min}}`$ has been chosen very small. Since the momentum of a cluster is given by the vector sum of its constituent hadrons, it would resum opposite and compensating $`p_{}`$ kicks imparted to hadrons in the fragmentation process. The cluster direction should therefore be a better measure of the $`q\overline{q}`$ axis than that provided by the individual hadron momenta. It is the angular distribution of this axis that relates back to the polarization character of the $`\gamma ^{}/Z^0q\overline{q}`$ decay, and that we want to be reflected in our extracted $`F_\mathrm{L}(x)`$ and $`F_\mathrm{T}(x)`$. Therefore it would be an improvement to rotate all hadrons in a cluster to sit along the cluster direction. That is, the $`\theta `$ of a hadron is redefined while its $`x`$ value is unchanged.

At this level there is no contradiction with the smearing approach studied earlier. Then we smeared the simple partonic angular shapes to arrive at realistic hadronic ones to compare with data, now we un-smear the hadronic angles to approach the simple partonic distributions. There is one advantage, however: the clustering approach is not sensitive to the width of the $`p_{}`$ distribution, i.e. the $`\sigma `$ parameter, unlike the smearing procedure. Of course, the $`p_{}`$ width still affects the typical error between the $`q\overline{q}`$ and cluster axes.

When considering multijet production, the $`y_{\mathrm{min}}`$ choice does become relevant, with $`\mu ^2=y_{\mathrm{min}}E_{\mathrm{cm}}^2`$ acting as a factorization scale, as noted above. For a large $`y_{\mathrm{min}}`$ all activity is clustered into the two quark jets, and neither the gluon structure nor $`F_\mathrm{L}(x)`$ would be probed. For $`y_{\mathrm{min}}0`$ each hadron or parton is a cluster unto itself, and we are back at the starting point. So obviously some intermediate scale is to be preferred. Given that the typical hadronization $`p_{}`$ width is $`0.4`$ GeV, with a tail to larger values, one would conclude that clustering up to $`p_{}1`$ GeV would be a sensible minimum to eliminate the bulk of the hadronization $`p_{}`$ smearing. At the $`Z^0`$ peak this translates into $`y_{\mathrm{min}}\text{ }\stackrel{>}{}\text{ }0.0001`$. In the upper end, we want to stay with a picture of multiple gluon emission as the norm, i.e. retain $`F_\mathrm{L}(x)`$ as an inclusive quantity, in order not to overlap with traditional studies of gluon jets in exclusive three-jet events. Since the average number of clusters per event is three for $`y_{\mathrm{min}}0.0025`$, we conclude that $`0.0001\text{ }\stackrel{<}{}\text{ }y_{\mathrm{min}}\text{ }\stackrel{<}{}\text{ }0.0025`$ is a reasonable range, over which to study a scale-dependent $`F_\mathrm{L}(x,\mu )`$.

In Fig. 3, the distribution of event-by-event $`x_E`$-weighted and -integrated hadronization corrections $`\mathrm{\Delta }\sigma _\mathrm{L}/\sigma _{\mathrm{tot}}=(\sigma _\mathrm{L}^h\sigma _\mathrm{L}^p)/\sigma _{\mathrm{tot}}`$ is shown for some different $`y_{\mathrm{min}}`$ scales, for events generated with parton showers and string fragmentation (superscripts $`h`$ and $`p`$ stand for hadron and parton level, respectively). We note the significant width of these distributions, showing that event-by-event fluctuations in the hadronization process are important and can be of either sign. Even if small by comparison, the mean $`\mathrm{\Delta }\sigma _\mathrm{L}/\sigma _{\mathrm{tot}}`$ does show a systematic bias, positive for small $`y_{\mathrm{min}}`$ and negative for large $`y_{\mathrm{min}}`$. That is, at small $`y_{\mathrm{min}}`$ the hadronization smearing wins over the string effects, while it is the other way around for large $`y_{\mathrm{min}}`$ — but remember that this is only true when averaging over many events. Nevertheless, one possible criterion for a good choice of $`y_{\mathrm{min}}`$ would be where the two effects cancel, which then gives $`y_{\mathrm{min}}0.0002`$, i.e. $`\mu 1.3`$ GeV. While a sensible reference value, one should not take this particular value too seriously, since it is for one specific model, and for one specific set of model parameters. Somewhat different parameter values, like for the parton shower cut-off $`Q_01`$ GeV, defining the parton level of the events studied, could lead to slightly different “preferred” $`\mu `$ values.

It is important to note that we here have been considering the $`x`$-integrated quantity. This is of relevance if one e.g. would like to extract an $`\alpha _\mathrm{s}`$ from an $`\sigma _\mathrm{L}`$ measurement, and so not uninteresting. For the purpose of determining the differential $`x`$ distribution, $`F_\mathrm{L}(x)`$, however, one would have reason to fear that any bias could have an $`x`$ dependence that would not be caught. In the string model, a string piece connecting two partons is boosted by an increasing velocity vector as the relative opening angle between the partons is decreased, and so the string effects spread upwards to larger $`x`$ values. A warning signal is then that $`\mathrm{\Delta }F_\mathrm{L}`$ does depend quite significantly on $`\overline{T}`$, Fig. 4, i.e. clustering does not reduce the $`\overline{T}`$ dependence noted in Fig. 2, but mainly shifts the overall level. Since $`\overline{T}`$ probes the topology of events, we also do expect this topology to reflect itself in an $`x`$ dependence of hadronization corrections. As in the previous studies, this dependence is then likely to show up mainly in the lower end of the $`x`$ range. At larger $`x`$, hadrons are rather well aligned with the jet axes, so, even with $`x`$-weighting, the few particles out there give a small contribution to the $`\mathrm{\Delta }F_\mathrm{L}`$. In Fig 5, the relative difference between the inclusive $`F_L(x_p)`$ and $`F_L^{cluster}(x_p)`$, obtained by replacing hadron angles with cluster angles for different $`y_{min}`$ values, is shown. It is clearly seen, indeed, that the hadronization corrections are only important at low $`x_p`$, unless $`y_{\mathrm{min}}`$ is chosen too high. For the transverse component of the fragmentation function, $`F_T(x_p)`$, corrections have the same absolute amplitude but the opposite sign.

On the up side, the string effect has its perturbative equivalent in the colour dipole . That is, asymmetries also exist in the production of soft gluons around the direction of the harder partons of an event. Such soft parton emission, below the cut-off scale $`Q_01`$ GeV we have used, would thus largely fill in the same regions as the nonperturbative hadron production, and with the same topology dependence. If one takes Local Parton-Hadron Duality  seriously, this equivalence should come very close. Our proposed strategy, to reset the $`\theta `$ angle of particles to that of the cluster they belong to, would be applicable also to such perturbatively calculated parton topologies.

## 5 Summary

The coherence phenomenon  kills the concept of gluon fragmentation functions that can be defined independently of the environment they are found in. The “hump-backed” shape of inclusive $`x`$ spectra  is an excellent illustration: by coherence the multiplication of partons/hadrons at small $`x`$ is much less than if the hard partons could radiate/hadronize independently. The immediate consequence is that the expected “softer gluon than quark jets” picture is difficult to test. This impacts both on studies of gluon jets directly in identified three-jet events and indirectly via $`F_\mathrm{L}(x)`$. In this article we have illustrated some of these issues for the latter observable.

It appears safe to conclude that a straightforward extraction of $`F_\mathrm{L}(x)`$ from hadron angular distributions exaggerates the rate of particles at small $`x`$ that should be attributed to gluon jets, since even the hadronization of pure $`q\overline{q}`$ events induces a ‘false’ $`F_\mathrm{L}(x)`$ by $`p_{}`$ smearing. We have also here shown that a symmetric smearing around jet axes introduces a bias in the other direction, since it misses important string/drag effects that tend to make three-jet events more two-jetlike. In summary, there is no model-independent extraction of a unique $`F_\mathrm{L}(x)`$, especially not at small $`x`$ values.

We therefore propose to introduce a scale-dependent quantity $`F_\mathrm{L}(x,\mu ^2)`$. Particles are clustered, e.g. with the Durham algorithm, and thereafter assigned the $`\theta `$ angle of the cluster they belong to, while retaining their $`x`$ value. Thus $`\mu ^2`$ sets an ‘un-smearing’ scale, below which $`p_{}`$ fluctuations are killed. We find that a $`\mu 1.3`$ GeV gives opposite and compensating $`p_{}`$ smearing and string effects in Jetset simulations. While the exact number certainly is model-dependent, the order is a sensible one, given that the average hadronization $`p_{}`$ is of the order of 0.4 GeV. If this then sets a reasonably lower limit, an upper one is related to the desire to stay away from the region of exclusive two- or three-jet events. Over an intermediate range, one could imagine several experimental determinations providing the scale dependence. We also remind that, so far, our studies have only been intended for LEP1 energies. Coverage of a wider energy range, e.g. at LEP2, introduces $`s`$ as a further scale of the process and allows more differential tests.


format: text

text:
# Commensurate and Incommensurate Vortex Lattice Melting in Periodic Pinning Arrays

## I Introduction

Vortices interacting with periodic pinning arrays have attracted growing interest with the advances of nano-lithographic techniques in which the size and shape of individual pinning sites, as well as the geometry of the pinning array, can be carefully controlled. These pinning arrays can be composed of 1D periodic modulations , micro-holes  or magnetic dots  in which the holes or dots have a diameter $`d`$ significantly less than the dot spacing $`a`$. The holes and dots suppress the order parameter, allowing the flux lines in the mixed state to lower their energy, and thereby pinning the vortices. Strong commensurability effects in these systems are observed at the matching fields, when the number of particles equals an integer multiple of the number of divots . Transport measurements near $`T_c`$ above the first matching field indicate that some of the vortices are only weakly pinned in the interstitial regions. These vortices become mobile at a lower applied current than the vortices at the pinning sites . Direct evidence for commensurate states and interstitial vortices has been obtained using Lorentz force microscopy  and scanning Hall probe microscopy . Magnetization measurements by Baert et al.  interpreted the loss of pinning above the first matching field near $`T_c`$ as due to a coexistence of an interstitial vortex liquid with a pinned vortex solid at the pinning sites.

In addition to the technological applications of superconductors which require strong pinning, vortices interacting with periodic pinning arrays are also an ideal realization of an elastic lattice interacting with a periodic substrate which can be generalized to a wide variety of condensed matter systems including atoms adsorbed on a surface , colloidal crystals interacting with interfering laser beams, and vortices in Josephson-junction arrays . A very similar system to vortices in periodic arrays that has recently been experimentally realized is colloidal particles interacting with optical trap arrays . Recent experiments in these systems for colloids interacting with a square array observed an interesting transition from a commensurate square colloidal crystal to a triangular colloidal crystal for increasing colloidal densities when the colloidal interactions began to dominate over the substrate .

Although studies of the effects of temperature have been conducted on the dynamics of vortices interacting with periodic pinning arrays , the static melting has not been investigated. Here one could expect that vortex lattice melting behavior may be significantly different at commensurate vs incommensurate fields. Melting of vortices in wire networks has been studied  for systems for filling fractions far less than one with Monte-Carlo simulations. The Monte-Carlo simulations of vortices in wire networks at low filling factors by Franz and Teitel  have shown a transition from a pinned solid to a floating solid phase, both of which have the same symmetry as the underlying pinning lattice. For higher temperatures this floating solid phase was found to melt into a liquid. In the wire network case there are no interstitial regions, and at the matching field the floating solid phase is absent . The wire network system has several differences from vortices in periodic arrays in that in the former the vortices can be considered to move in an egg-carton potential. For vortices in periodic arrays the vortices would move in a muffin-tin substrate in which there are flat potential regions between the minima. In this system at low temperatures at commensurate matching, particles in the overlying lattice will sit in the divots or energy minima, while at the incommensurate matching, some particles will sit in the interstitial regions between the divots . Even at the matching field it can be possible for vortices to move into the interstitial regions when the vortex-vortex interactions become dominant or thermal fluctuations are strong enough.

Another difference between vortices in wire networks compared to vortices in periodic pinning arrays appears for fields greater than the first matching field. The additional vortices in the imaging experiments by Harada et al.  and simulations by Reichhardt et al  indicate that in periodic pinning arrays a wide variety of interstitial vortex crystals can be stabilized at different matching fields.

The vortices located in the interstitial regions can have significantly different static, dynamic and thermal properties than the commensurate particles, and the coupling between the two different species can give rise to interesting behavior. Recent imaging experiments by Grigorenko et al  using high resolution scanning Hall probe microscopy near the matching field $`B_\varphi `$ indicate that, for fields slightly below and above the first matching field, over time the vacancies and interstitials appear to show site to site hopping. The nature of the hopping dynamics is however not known, such as whether the interstitial vortices jump from one interstitial site to another or if the interstitials undergo an exchange process with the vortices at the pinning sites. These experiments were done for square arrays and it is not known how the vortex hopping might occur in different geometries such as triangular pinning arrays.

The melting of interstitial vortices which coexist with strongly pinned vortices has also been proposed by Radzihovsky  for vortices interacting with randomly placed columnar defects where there are more vortices than pinning sites.

In this work we report the study of two-dimensional melting of vortices interacting with a square muffin tin pinning potential with the use of molecular dynamics simulations. For the commensurate case at low substrate strengths we observe a novel order-order transition from a square pinned vortex lattice with long-range order to a triangular ordered floating solid phase with quasi-long range order. Our simulations also allow us to examine the vortex dynamics at the transition into the floating solid where the vortices show a collective 1D avalanche motion along the pinning rows. For increased substrate strength the floating solid phase is lost and the pinned solid melts directly to a liquid. In the incommensurate case for vortex densities slightly above the first matching field for strong substrates we observe a multi-stage melting in which the initial vortex motion occurs by an exchange process of interstitial and pinned vortices.

We have also examined the effects of temperature for higher matching and non-matching fields greater than one. At the higher matching fields we find that the vortices form a solid at low temperatures. For higher temperatures we observe a coexistence of an interstitial vortex liquid with a pinned vortex solid in agreement with the experiments of Baert et al. At high incommensurate fields the extra vortices or vacancies in the interstitial vortex liquid become mobile. We observe the same general behaviors for triangular pinning arrays; however, at the incommensurate fields the motion of extra vortices or vacancies can occur by novel collective ring excitations where a group of vortices rotate around a vortex in a pinning site. These excitations occur when an incommensurate number of vortices cage a vortex located at a pinning site.

Our system should not only be of relevance for vortices and colloids interacting with periodic substrates but also for the square to hexatic vortex transitions observed in certain superconductors , where little is known about the dynamics at the transitions.

## II Simulation

We model a 2D system with periodic boundary conditions in $`x`$ and $`y`$ with $`N_v`$ vortices interacting with a square array of $`N_p`$ pinning sites. We use a Langevin simulation for vortices at finite temperature using the model of Brass and Jensen . We consider the vortex motion to be overdamped and integrate the following equation of motion.

$$𝐟_i=\eta 𝐯_i=𝐟_i^{vv}+𝐟_i^{vp}+𝐟_i^T$$
(1)
Here $`𝐟_i`$ is the total force acting on vortex $`i`$ and $`\eta `$ is the damping coefficient. The repulsive interaction of vortex $`i`$ with the other vortices is $`𝐟_i^{vv}=_{j=1}^{N_v}f_0K_1(|𝐫_i𝐫_j|/\lambda )\widehat{𝐫}_{ij}`$. Here, $`\lambda `$ is the penetration depth and $`K_1(r/\lambda )`$ is a modified Bessel function which falls off exponentially for $`r>\lambda `$, allowing us to place a cutoff in the interactions at $`r=6\lambda `$ for computational efficiency. The pinning is modeled as $`N_p`$ parabolic traps with $`𝐟_i^{vp}=_{k=1}^{N_p}(f_p/r_p)(|𝐫_i𝐫_k^{(p)}|)\mathrm{\Theta }(r_p|𝐫_i𝐫_k^{(p)}|)\widehat{𝐫}_{ik}^{(p)}`$. Here, $`\mathrm{\Theta }`$ is the Heaviside step function, $`r_p`$ is the range of the pinning potential, $`f_p`$ is the maximum pinning force of the wells, $`\widehat{𝐫}_{ij}=(𝐫_i𝐫_j)/|𝐫_i𝐫_j|`$ and $`\widehat{𝐫}_{ik}^{(p)}=(𝐫_i𝐫_k^{(p)})/|𝐫_i𝐫_k^{(p)}|`$. The thermal force has the properties $`<f_i^T(t)>=0`$ and $`<f_i^T(t)f_j^T(t^{^{}})>=2\eta k_BT\delta _{ij}\delta (tt^{^{}})`$. In this work we set $`\eta =k_B=1`$. All lengths, fields and forces are given in units of $`\lambda `$, $`\mathrm{\Phi }_0/\lambda ^2`$ and $`f_0`$ respectively. For most of the results shown here the density of vortices is $`0.25\mathrm{\Phi }_0/\lambda ^2`$. The initial ground state of the system is found using simulated annealing. We then slowly increase $`T`$ and analyze the coordination number, $`P_n`$ (from the Voronoi or Wigner Seitz cell construction), and the structure factor,

$$S(𝐤)=(1/L^2)\underset{i,j}{}e^{i𝐤[𝐫_i𝐫_j]}$$

where $`L`$ is the system size length.

Another measure we employ to characterize the thermally induced transition is the sum of the displacements of the vortices from their initial positions at $`T=0`$ to their positions at a higher $`T`$:

$$d_w(T)=\underset{i=1}{\overset{N_v}{}}|𝐫_i(T)𝐫_i(0)|^2$$

If the vortex lattice configuration is the same as the initial state than $`d_w=0`$. If a portion of the vortices move but then become frozen again $`d_w`$ will jump up. If the vortex lattice becomes molten $`d_w`$ will increase to a saturation value of $`L/2`$.

To visualize vortex motion we plot the trajectories that the vortices follow for a period of time. To examine the dynamic response in the phases we can add a driving term $`𝐟_i^d`$ representing a Lorentz force from an applied current to Eq. 1 and sum over the net vortex velocities $`V_x=_{j=1}^{N_v}𝐯_i\widehat{𝐱}`$. To analyze finite size scaling in $`S(𝐤)`$ we conduct simulations with the same vortex and pin densities for system sizes with a length $`L`$ for systems ranging from $`L=24\lambda `$ to $`L=60\lambda `$ with $`N_v=144`$ to $`N_v=900`$.

## III Commensurate Melting

In Fig. 1 we present the evolution of the coordination numbers $`P_6`$ and $`P_4`$, the average power of the secondary peaks in the structure factor $`<S(𝐤)>`$, and the vortex displacements $`d_w`$ for increasing $`T`$ for a system at the matching field with $`L=48\lambda `$. Here $`f_p=0.03f_0`$, $`r_p=0.17\lambda `$ and $`a=2.0\lambda `$. For $`0<T<0.0012`$ the pinned vortex lattice has the same square symmetry as the pinning lattice, as indicated by $`P_41.0`$. The initial increase in $`d_w`$ for $`T<0.0012`$ is due to the vortices moving randomly within pinning sites but remaining confined to a distance $`r_p`$. Near $`T=0.0012`$ there is a structural transition in the vortex lattice from square lattice to triangular, indicated by the sudden drop in $`P_4`$ to $`0`$ and the large increase in $`P_6`$ to $`0.98`$. In Fig. 1(b) we plot the average $`<S(𝐤)>`$ of the appropriate secondary peaks in $`S(𝐤)`$. For $`T<0.0012`$ we observe only four secondary peaks in $`S(𝐤)`$ as indicated in the left inset in Fig. 1(a). At the square-triangular transition there is initially a large drop in $`<S(𝐤)>`$ from $`0.95`$ to $`0.2`$ and then a recovery to $`0.55`$. The position of the secondary peaks shifts at the transition from four-fold to six-fold order. During the transition a mixture of both four and six-fold peaks appear, reducing the overall power at each peak. In the triangular phase there is almost no power for the four-fold order. The square to triangular transition can also be observed as a jump in $`d_w`$ (Fig. 1(c)), coinciding with the drops in $`P_4`$ and $`<S(𝐤)>`$, which is due to the shifting of the vortex positions during the transition. In the triangular lattice phase, $`d_w`$ remains roughly constant, indicating that diffusion is not occurring and that the vortex lattice is still a solid. Near $`T=0.0049`$ the triangular vortex lattice melts into a liquid, indicated by the increase in $`d_w`$ as well as the drops in $`<S(𝐤)>`$ and $`P_6`$. We note that $`P_6`$ will not drop to zero for a liquid state since even for random particle distributions a portion of the vortices will always have six-fold coordination numbers. We have also conducted simulations where we cool down in temperature starting at $`T=0.008`$. In this case we do not observe any hysteresis in the liquid to triangular phase; however, hysteresis is present in the triangular to square transition, suggesting that the transition is first order in nature. We have also examined the response of the phase to an applied drive and find that only the square pinned phase has a finite critical depinning force.

The finite size scaling behavior of the structure factor

$$S(𝐤)/L^2L^\theta $$

which is described more fully in  gives us further information about the nature of the phases. For long range order, $`\theta =0`$, for floating solids or hexatic phases $`0<\theta <2.0`$, and for the liquid phase $`\theta =2.0`$. For different system sizes we average the power of the secondary peaks in $`<S(𝐤)>`$ for 150 frames. In the inset of Fig. 1(b) we plot $`<S(𝐤)>/L^2`$ (where $`N_vL^2`$) versus $`N_v`$ for the three different phases. In the commensurate phase ($`T=0.0005`$) we find $`\theta 0.0`$ corresponding to the pinned lattice. For $`T=0.003`$, in the triangular phase, $`\theta =0.33\pm 0.04`$ which is consistent with a floating solid phase. In the floating solid phase we find a small variation of $`\theta `$ with $`T`$ which we will address elsewhere. In the liquid phase ($`T=0.007`$) we find $`\theta =1.98\pm 0.04`$.

### A Dynamics of the Order-Order Transition

In order to gain insight into the dynamics of the square to triangular transition, in Fig. 2 we show the real space images as well as the individual vortex trajectories in the pinned solid, the pinned-solid to floating solid transition, and the liquid phase for the system in Fig. 1. In the pinned solid phase shown in Fig. 2(a) ($`T=0.0005`$) we find a square vortex lattice with little or no detectable vortex motion. Fig. 2(b) ($`T=0.00135`$) shows the dynamics at the beginning of the pinned to floating solid transition. Here alternate rows of vortices move in a sudden 1D avalanche where the vortices in one row move in a collective manner while the other rows remain static. All the vortices in a given row shift together as the vortices jump out of the pinning sites and move to the interstitial regions. The motionless rows shown in Fig. 2(b) also begin to move at a slightly later time. This shifting of the vortex positions produces the jump in $`d_w`$ at $`T=0.0012`$ seen in Fig. 1(c) at the square to floating solid transition. In Fig. 2(c) ($`T=0.0015`$), in the floating solid phase, the vortex lattice has triangular order and shows increased thermal wandering around the equilibrium positions as can be seen in the vortex trajectories. This increased thermal wandering also smears the secondary peaks in $`S(𝐤)`$, reducing the maximum value to $`0.55`$, as in Fig. 1(b). In Fig. 2(d) ($`T=0.0065`$), in the molten state the vortices are disordered and diffuse in a random manner corresponding to the large value of $`d_w`$ in Fig. 1(c).

### B Phase Diagram for the Commensurate Case

By conducting a series of simulations with different substrate strengths and using the criteria from Fig. 1 to determine the onset of the different phases, we construct the phase diagram in Fig. 3. The range of temperatures over which the commensurate solid appears increases linearly for increasing $`f_p`$. The floating solid phase only appears for $`f_p<0.125f_0`$. The line separating the floating solid to liquid phase is roughly constant with temperature for lower substrate strengths at a value that corresponds to roughly the zero pinning melting temperature, $`T_m=0.0048`$ (as also obtained with the diffusion and structure factor measurements). All three phase boundaries can be understood from simple considerations. The floating solid to liquid line is independent of $`f_p`$ because it is determined by the vortex-vortex interactions and not the pinning force. For $`f_p0.125f_0`$ the commensurate solid melts directly into a liquid and $`<S(𝐤)>`$ has the behavior shown in the inset of Fig. 3. The melting transition is seen to shift to higher temperatures with increasing pinning strength, as the melting occurs when thermal fluctuations are able to overcome the increasing pinning strength. Finally the commensurate-to-floating solid follows $`f_pT`$ with the low $`f_p`$ being washed out before the vortex-vortex interactions are washed out. The nature of the floating solid to liquid transition, such as whether a dislocation-mediated melting transition occurs , is beyond the scope of this work.

To further compare our results with those of vortices in wire networks , in which the floating solid phase with the same symmetry as the pinned phase is observed as a function of commensurability, in our case the pinned to floating solid transition is seen as a function of substrate strength. Further, we observe a square to triangular transition due to the fact that a portion of the vortices shift into pin free regions, which is not possible in wire networks.

We note that in our simulation we imposed nearly square boundary conditions which can create defects or stress in the triangular lattice. We do not find any defects in our triangular lattice and our finite size scaling is consistent with that in . It is also possible that the square pinned solid phase may be enhanced for a higher temperature range by the square boundary conditions; however, the general trends of the phase diagram should still hold.

## IV Incommensurate Melting Near the First Matching Field

We now consider the vortex behavior above the first matching field when there are more vortices than pins. We find that for large substrate strengths, $`f_p>0.25f_0`$, a commensurate sub-lattice forms, with the additional vortices sitting in the interstitial regions. In Fig. 4 we show the melting behavior for a system with $`N_v=1.063N_p`$ and $`f_p=0.6f_0`$. In the imaging experiments of Grigorenko et al. , the interstitial motion was observed for a field of $`B/B_\varphi =1.04`$. Here, for $`T<0.0016`$ both the commensurate vortices and the interstitial vortices remain pinned as seen in Fig. 4(a). At these low temperatures, the interstitial vortices are pinned by the potential cage created by vortices at the pinning sites. At $`T>0.0016`$ vortex diffusion begins to occur as indicated by the increase in $`d_w`$ in Fig. 5(b). This is due to the onset of motion of the incommensurations as seen in the vortex trajectories in Fig.4(b). The motion is not restricted to interstitial vortices. Instead, interstitial and pinned vortices exchange places, and over a period of time all the vortices take part in the motion. The exchange process can be understood by considering that the interstitial vortex will produce a force $`f_{in}`$ on the nearby commensurate vortices which effectively reduces the pinning force to $`|f_p||f_{in}|`$. Since the melting transition is a function of the pinning force the melting temperature will be reduced.

In Fig. 5(a) the onset of the incommensuration motion is also marked by a drop in $`<S(𝐤)>`$ from $`0.96`$ to approximately $`0.7`$. In contrast to the sharp drop in $`<S(𝐤)>`$ seen in the commensurate case, here $`<S(𝐤)>`$ still retains strong four-fold order due to the presence of the background lattice of pinned vortices as seen in Fig. 4(b). For higher temperatures $`<S(𝐤)>`$ gradually decreases with most order being lost for $`T>0.0175`$, corresponding to the melting temperature for the commensurate case with the same substrate strength. For simulations where we cool down from $`T>0.0175`$ we do not observe any hysteresis. We have also conducted simulations in which we apply a constant drive of $`f_d=0.012`$, gradually increase $`T`$ and measure the average vortex velocity $`<V_x>`$. For $`T<0.0016`$, in the pinned interstitial phase $`<V_x>=0`$. At the interstitial transition there is a jump in $`<V_x>`$ indicating the onset of vortex motion. Finally at the overall vortex melting transition, $`<V_x>`$ jumps to a value corresponding to the entire lattice flowing.

## V Melting for Higher Commensurate and Incommensurate Fields

### A Commensurate Melting

We have also investigated the melting in square and triangular pinning arrays for the higher matching fields $`B/B_\varphi =2,3,4`$, $`5`$ and $`8`$. At these fields ordered vortex crystals are formed . In Fig. 6 we show the vortex positions and trajectories for the square pinning array with $`B/B_\varphi =4.0`$ for increasing temperatures. For $`T<0.004`$ the vortex lattice is ordered as seen in experiments and simulations. For $`T>0.004`$ \[Fig. 6(b)\] the interstitial vortex lattice melts; however, the vortices at the pinning sites remain immobile so that a coexistence of a solid and a liquid phase occurs. The interstitial vortices are constrained to move in a square grid and there is a region around the pinning sites which the interstitial vortices do not enter due to the repulsion of the vortices located at the pinning sites. Under an applied drive the depinning force is finite for low temperatures but goes to zero at the onset of the interstitial vortex liquid transition. This result is in agreement with the interpretations of the experiments by Baert et al. in which the pinning force goes to zero near $`T_c`$ when interstitial vortices are present, but is still finite for fields at which vortices are only located at the pinning sites. For increased $`T`$, as in Fig. 6(c), the vortices at the pinning sites become depinned as well and the entire lattice is in the molten state. The general behavior of the melting from Fig. 6 is also seen at the other matching fields we have investigated with transitions from the solid, to solid-liquid coexistence, to the liquid state.

### B Incommensurate Melting

In Fig. 7 we show the melting behavior for the incommensurate field $`B/B_\varphi =4.1`$. Here the onset of vortex diffusion occurs for a much lower temperature than for $`B/B_\varphi =4.0`$. The vortex motion occurs at defect sites where additional interstitial vortices are present in the ordered $`B/B_\varphi =4.0`$ interstitial vortex crystal. These extra interstitial vortices are less strongly pinned than the other interstitial vortices. The vortex motion can occur by the continuous motion of a single vortex over a certain distance, or by a pulse like motion in which a series of individual vortices move by one lattice constant as the pulse moves through. If the trajectories are drawn for a longer time period the interstitial vortex diffusion can be seen to occur through the entire lattice. As the temperature is increased the vortices at the pinning sites will become depinned.

### C Collective Ring Excitations for Incommensurate Fields with Triangular Pinning

We have found that the melting behavior for square and triangular pinning arrays is similar at the matching fields; however, at the incommensurate fields we observe clear differences between the two pinning geometries. In Fig. 8 we show the melting behavior for a system with triangular pinning for $`B/B_\varphi =3.08`$. At low temperatures the system is frozen but as $`T`$ is increased motion begins to occur in the form of ring excitations. Here a portion of the vortices rotate around a vortex at a pinning site; however, no net diffusion of interstitial vortices is occurring. For $`B/B_\varphi =3.0`$ each vortex at a pinning site is surrounded by six interstitial vortices. For fields above or below $`B/B_\varphi =3.0`$ a portion of the pinned vortices will be surrounded by 5 or 7 vortices, and it is at these sites where the ring excitations occur. Due to the symmetry of the pinning array the vortices at the pinning sites create a potential around each pinning site which has six minima, so that at $`B/B_\varphi =3.0`$ the vortex configuration is particularly stable. When there are seven or five vortices around a pinning sites, they are incommensurate with the six mimima, and therefore the thermal kicks can cause the interstitial vortices to slowly rotate.

As the applied field is moved further away from the matching field, the number of incommensurations increases and the number of ring excitations increases as well. Such ring excitations are observed in triangular pinning arrays for most incommensurate fields but with different numbers of vortices forming the ring. As the temperature increases, the initial diffusion of the vortices occurs via a hopping mechanism in which the extra vortex in the rotating ring hops out of the ring and starts the rotation of a different ring. The rate of this hopping increases with the temperature.

## VI Summary

In summary, we have investigated the melting of vortex lattices interacting with periodic pinning arrays for both the commensurate and incommensurate case. For the first matching field with weak substrates we have observed a transition in the vortex lattice from a square commensurate solid to a triangular floating solid which melts into a liquid at higher temperatures. This transition occurs by 1D collective motion of vortices shifting from the pinning sites into the interstitial regions. For strong substrates the floating solid phase is lost and the vortex lattice melts directly into a liquid. For fields slightly above the first matching field for strong disorder we observe a multi-stage melting where the interstitial vortices are highly mobile consistent with recent imaging experiments. The motion of the interstitial vortices occurs by the exchange of interstitial with pinned vortices rather than with the interstitials hopping from one interstitial site to another. As the temperature increases more vortices become mobile.

We have also examined the melting behavior for matching fields greater than one. At low temperature the interstitial vortices form a crystalline state. At higher temperatures we observe a coexistence of an interstitial vortex liquid with a commensurate pinned vortex lattice. The interstitial vortices can diffuse in the regions between the pinned vortices. We also observe a depletetion zone around the pinned vortices which interstitial vortices do not enter due to the repulsion from the pinned vortex. At higher temperatures the entire vortex lattice melts. At higher incommensurate fields the additional defects (extra interstitials or vacancies) in the interstitial vortex crystal become mobile and diffuse in the interstitial regions. With triangular pinning arrays at incommensurate fields the initial vortex motion occurs through a novel collective ring excitation where interstitial vortices can rotate around a pinned vortex. These rotations occur when an incommensurate number of interstitial vortices surround a pinned vortex.

We briefly discuss some experimental systems in which these effects could be observed. For vortex lattices in periodic pinning arrays, vortex lattice melting is most relevant to high Tc superconductors. Our results are only for a 2D system whereas 3D effects can be relevant to melting. Our model would thus best describe high-temperature superconductors with periodic columnar defects where the vortices would have line like behavior. However, several of the results here should still be relevant for low temperature superconductors. For the low temperature superconductors melting would only occur very near $`T_c`$. For strong pinning such as arrays of holes, a floating solid phase at the first matching field would not be observable and the system should go directly from the pinned solid to the normal phase with a very small region of a disordered vortex lattice. For weak pinning such as small magnetic dots or weak defect arrays the pinned solid to floating solid transition should occur well below $`T_c`$ and should be observable for the low temperature superconductors. For a low temperature superconductor the phase diagram in Fig. 3 would be modified, with the floating solid to liquid transition replaced by the floating solid to normal transition, although again a small vortex liquid state may occur just below $`T_c`$. Also the pinned solid to vortex liquid transition would remain flat rather than increasing with $`f_p`$. The solid to floating solid transition could be imaged with scanning Hall probe measurements, neutron scattering, or Lorentz microscopy. In addition transport measurements would also be able to reveal the loss of pinning at the transition.

The vortex behavior at the incommensurate fields should also be visible in the low-temperature superconductors. The recent imaging experiments of Grigorenko et al.  have already found evidence for the motion of the highly mobile interstitial vortices in these types of systems. These same imaging techniques should be able to image the pinned vortex lattice coexisting with the interstitial vortex liquid. In triangular pinning arrays collective ring excitations would appear as a smeared ring around pinning sites with the Hall probe arrays and the dynamics could be directly imaged with Lorentz microscopy.

Our work is also directly relevant to colloidal particles interacting with periodic pinning. A particularly promising realization of this system is optical trap arrays in which the pinning strength can be easily tuned. In such a system our phase diagram can be directly tested. Recent experiments have already seen evidence for a square to triangular transition as the colloid density is increased . In addition the dynamics of the interstitial colloids can be directly imaged with video microscopy to determine if the interstitials hop directly from one interstitial site to another, or if they move through an exchange process with a pinned colloid. In addition it should be also be possible to observe the melting behavior for higher matching colloidal densities and the collective ring excitations for colloids interacting with triangular optical trap arrays.

We thank S. Bending, L. DeLong, S. Field, D. Grier, N. Grønbech-Jensen, P. Korda, V. Moshchalkov, and G. Spalding for useful discussions. Funding provided by NSF-DMR-9985978, CLC, CULAR, and DOE grant W-7405-EBG-36.


text:
# Software Corrections of Vocal Disorders

## Phase Space Reconstruction

The first simple difference between dysphonic and healthy voices is shown in Fig.1, where the time evolution of the amplitude of a microphone-registered sound is represented. The upper panel could be interpreted as a highly noisy time series, but careful investigations reveal that this is not the case: Applying a simple low-pass filter would only introduce a distortion bigger than the original noise level. Some of the noise-like structures belong to the time series and one has to be able to correctly identify what is worth keeping and what has to be eliminated during the correction procedure.

From a theoretical point of view, this paper relies on the theory of dynamical systems and deterministic chaos; the former implies that the time evolution is defined in some phase space, the latter offers a striking explanation for irregular behaviour and anomalies in systems which do not seem to be inherently stochastic. Even very simple chaotic dynamical systems can exhibit strongly irregular time evolution without random inputs.

Consider for a moment a purely deterministic system. Once its present state is fixed, the states at all future times are determined as well. Thus it will be important to establish a vector space, the so-called *phase space*, for the system such that specifying a point in this space specifies the state of the system, and vice versa. This implies that we can study the dynamics of the system by studying the dynamics of the corresponding phase space points.

Unfortunately what we observe in an experiment is not a phase space object but a time series, most likely only a scalar sequence of measurements. We therefore have to convert the observations into state vectors: This is the problem of *phase space reconstruction* which is technically solved by the method of delays or related constructions takens ; sauer ; kantz1 . Most commonly, the time series is a sequence of scalar measurements of some quantity which depends on the current state of the system, considered at multiples of a fixed sampling time:

$$s_n=s(𝐱(n\mathrm{\Delta }t))+\eta _n,$$
(1)
namely we look at the system through some measurement function $`s()`$ and make observations only up to some random fluctuations $`\eta _n`$, the measurement noise. A *delay reonstruction* in $`m`$ dimensions is then formed by the vectors $`\beta _n`$, given as

$$\beta _n=(s_{n(m1)\nu },s_{n(m2)\nu },\mathrm{},s_{n\nu },s_n).$$
(2)
The time difference in number of samples $`\nu `$ or in time units $`\nu \mathrm{\Delta }t`$ between adjacent components of the delay vectors is referred to as the *lag* or *delay time*.

A number of embedding theorems are concerned with the question under which circumstances and to what extent the geometrical object formed by the vectors $`\beta _n`$ is equivalent to the original trajectory $`𝐱_n`$. Here equivalent means that they can be mapped onto each other by a uniquely invertible smooth map and, under quite general circumstances the attractor formed by $`\beta _n`$ is equivalent to the attractor in the unknown space in which the original system is living if the dimension $`m`$ of the delay coordinate space is sufficiently large.

## Recurrence plots of the voice

Human voices form an aperiodic and highly nonstationary signal. A sentence can be decomposed in subunits, called phonemes, which can be considered as different types of dynamics. Careful investigation of time and length scales shows that the sound wave characterizing a single phoneme (duration between 50 and 200 ms) has a characteristic profile (*pitch*) of about 5-15 ms. A kind of phase angle on this highly nontrivial oscillation will then identify the instantaneous amplitude. Thus a delay reconstruction should allow us to identify the actual phoneme and the phase inside it.

Looking at the Eq.2 one realizes that at least two parameters are involved in the delay reconstruction from a scalar time series, namely $`m`$ and $`\nu `$. Some recipes for an optimal tuning of them are available at the moment, but an adequate theory is still missing; in any case, a good tool is represented by the *recurrence plot* (see recurrenceplot ; recurrenceplot2 ). This method was used for the first time to study recurrencies and nonstationary behaviour occurring in dynamical systems. It allows to identify system properties that cannot be observed using other linear and nonlinear approaches and it is especially useful for analysis of nonstationary systems with high dimensional and/or noisy dynamics.

Recurrence plots are constructed on the basis of mutual distances between point belonging to the same trajectory. In the plane of indices $`i`$ and $`j`$ a dot is printed whenever the delay vectors $`\beta _i`$ and $`\beta _j`$ fulfill the relation $`|\beta _i\beta _j|<ϵ`$. So a recurrence plot depends also on the parameter $`ϵ`$. The Fig.2 proves that our delay vectors really represent meaningful states, where the line structure shows the approximate periodicity inside the phonemes and the number of intraphoneme neighbours.

The distance between two consecutive lines indicates the duration of the basic structure (*pitch*) inside a phoneme. We can see from Fig.2 that this involves approximately 200 points. Considering that the sentence has been sampled at a 22,05 kHz rate, this implies a profile of about 9 ms. In order to correctly identify such a structure, in the delay reconstruction of Eq.2, the product $`\nu m`$ (the so-called *time window*) has to be bigger than 200. A more detailed insight into the role of this two parameters is given in Ref. dimitris .

## Nonlinear noise reduction

Noise reduction means that one tries to decompose a time series into two components, one containing the signal, the other random fluctuations. Implicitly we always assume that the data represents an additive superposition of two different components which have to be distinguishable through some objective criterion.

The classical statistical tool for obtaining this distinction is the power spectrum. Random noise has a flat, or at least a broad, spectrum, whereas periodic or quasi-periodic signals have sharp spectral lines. After both components have been identified in the spectrum, a Wiener filter can be used to separate the time series accordingly. This approach fails for our purposes here because the undesirable part of the signal is not what is usually considered to be noise. It is very strongly correlated to the clean part of the signal, indeed it is part of the signal. Even if parts of the spectrum can be clearly associated with the signal, a separation into signal and noise fails for most parts of the frequency domain.

The filter we use has been proposed in lms  and arises from the chaotic deterministic systems field, where the determinism yields a criterion to distinguish the signal and the noise (which is supposed not to be deterministic). Let the time evolution of the signal be deterministic with an unknown map $`f`$. All that we have knowledge of are noisy measurements of this signal:

$$s_n=x_n+\eta _n,x_n=f(x_{nm},\mathrm{},x_{n1}).$$
(3)
Here $`s_n`$ is the scalar time series we can measure, $`x_n`$ the clean signal and $`\eta _n`$ the superimposed noise. In order to obtain an estimate $`\widehat{x}_{n_0m/2}`$ for the value of $`x_{n_0m/2}`$ we form delay vectors

$$𝐬_n=(s_{nm+1},\mathrm{},s_n)$$
(4)
and determine those which are close to $`𝐬_{n_0}`$. The average value of $`s_{nm/2}`$ is then used as a cleaned value $`\widehat{x}_{n_0m/2}`$:

$$\widehat{x}_{n_0m/2}=\frac{1}{|U_ϵ(𝐬_{n_0})|}\underset{𝐬_nU_ϵ(𝐬_{n_0})}{}s_{nm/2}.$$
(5)
Here $`|U_ϵ(𝐬_{n_0})|`$ denotes the number of elements of the neighbourhood $`U_ϵ(𝐬_{n_0})`$ of radius $`ϵ`$ around the point $`𝐬_{n_0}`$, which is never empty, no matter how small a value of $`ϵ`$ we choose: It always contains at least $`𝐬_{n_0}`$. This is good to know since we have to make some choice of $`ϵ`$ when we use this algorithm. It is guaranteed that if we choose $`ϵ`$ too small the worst thing that can happen is that the only neighbour found is $`𝐬_{n_0}`$ itself. This, however, yields the estimate $`\widehat{x}_{n_0m/2}=s_{n_0m/2}`$ which just means that no correction is made at all.

To get an impression of how the local projective noise reduction scheme works, assume that one has to eliminate noise from music stored on an old-fashioned long playing record, induced by scratches on the black disc. The task becomes almost trivial if one can make use of several samples of this LP. When playing them synchronously, the signal part of the different tracks id identical, whereas the noise part is independent: As a consequence of that, already a simple averaging would enhance the sound quality. In deterministic chaotic signals, this redundancy is stored in the past: Similar initial conditions will behave in a similar way, at least for short periods. In human voice signals there is no need to suppose a chaotic behavior, since every phoneme is made up of pitches in an almost periodic fashion. This means that every logical unit provides all the redundancy required for its filtering.

How can we make sure that all we do is to reduce noise without distorting the signal? In Fig.1 two typical time series, one related to a healthy voice, the other to a sick one, are compared. In the lower panel we can see the recurrence of structures of length approximately 9 ms inside a phoneme, in total agreement with the recurrence plot analysis. The person who has spoken this sentence is absolutely healthy and the corresponding signal looks very clean. The same repetition is not that clear in the upper panel, where a patient affected by a T1A glottis cancer was asked to say the same group of words in the same environmental conditions, namely in a professional quiet room. The signal looks very noisy, but listening to it one perceives only the sick voice and no indication of noise, at least of what one usually means saying noise.

Spectral analyses show the same result: The vocal disorder reveals itself as a noisy signal, but the nature of this noise is not the conventional one. It is not addivite and it is correlated with the other sub-part of the sentence one would like to isolate. In herzel2 ; ishizaka ; manfredi ; reuter ; pinto ; kumar  it is shown how some of the complexities observed in disordered voices are not caused by random external input to the vocal apparatus, but by the intrinsic nonlinear dynamics of vocal fold movement. Normal phonation corresponds to an essentially synchronized motion of all vibratory modes. A change of parameters such as muscle tension or localised vocal fold lesions may lead to a desynchronization of certain modes resulting in the appearance of these new features that look like noise.

The distinction of the signal in a deterministic part and in a stochastic component is then not possible. Nevertheless one would still like to extract the main structure from the signal, suppressing the noise-like features. For this purpose the previously described filter is worth applying to correct the disorders: If we choose suitably the involved parameters we are able to identify the structures inside a phoneme and to perform an average of them via Eq.5. Every point belonging to a structure is replaced by a local average of similar points coming from different structures (this is done searching for neighbours in the embedding space). The resulting signal does not show just an exact repetition of these parts because the averaging procedure acts only locally and the points involved in the computation varies from structure to structure. Also in a normal signal, in fact, the repetitions are not perfect.

One has then to be careful when choosing the $`ϵ`$: A too big value would result in a drastic averaging and the resulting signal would sound too artificial. On the other hand a small value of $`ϵ`$ is not able to perform any correction. In Fig.3 we show what happens when using the filter with a *good* set of parameters: The main shape of the sub-structures is preserved, but the noise-like features are attenuated. The meaning of good here is the following: The recurrence plot of Fig. 2 shows that the system generating the time series we are analyzing really possesses dynamical regimes. Since we want to be able to distinguish between different states, the vectors $`\beta _n`$ of the delay reconstruction has to comprise all the points belonging to a pitch. So the product $`\nu m`$ of the lag and the embedding dimension has to be approximately 200 in the example depicted on Fig. 2. A choice of a big $`m`$ results on big computational efforts; on the other hand a big $`\nu `$ produces worse results. We refer again to Ref. dimitris  for more details. The parameter $`ϵ`$ is used for the neighbourood relations between vectors in the embedding space. We draw a point in $`(i,j)`$ on the recurrence plot if $`|\beta _i\beta _j|<ϵ`$. The value of $`ϵ`$ has to be bigger than the noise level. A good way to tune it is through the recurrence plot: Whit a small $`ϵ`$ we get almost no recurrences, a very big $`ϵ`$ is such that the relation $`|\beta _i\beta _j|<ϵ`$ is fulfilled from almost all the vector pairs. The optimal value of $`ϵ`$ is mirrored on a recurrence plot where the lines are as long as possible, as thin as possible and all the recurrences belong to lines.

If we tune the parameters according to this recipe, the filtered signal sounds more *normal* than the original, even if some characteristic aspect of the voice has been lost. In the horizontal axis of Fig.3 we report the index of the time series in order to show the number of points involved in the filtering procedure. Every structure lasts 80 points; in a space of such a dimensionality, the vector having these elements as components is represented as a single point. The four structures of Fig.3 are neighbours because the involved sequences of points are approximately the same. As shown in dimitris , it is not necessary to consider such a big space, but one can skip some intermediate point, taking care to cover the full structure in any case.

## Results

We are dealing with three categories of subjects: (i) patients affected by T1A glottis cancer, a tumour confined to the glottis region with mobility of the vocal cords (ii) healthy patients and (iii) patients under medical treatment, operated via endoscopic laser or traditional lancet technique. The sentence spoken by all of them is the italian word *aiuole* (flower-beds), since it contains a lot of vowels: Patients suffering from dysphonia need a great effort in saying this word.

To classify the sentences we make use of the feature space introduced in fspace , whose components are special quantities extracted from the time series: Spectral factor, pseudo-entropy, pseudo-correlation dimension, first zero-crossing of the autocorrelation function, first Lyapunov exponent, prediction error, jitter, shimmer, peaks in the phoneme transition, residual noise. It is not necessary here to go into details, one has just to know that this space represents the proper object where a good classification of vocal pathologies is feasible.

In Fig.4 we see a projection of the feature space onto the spectral factor and the pseudo-entropy dimensions; we can see how healthy and sick patients populate different regions and which are the results of the medical treatments performed on four patients (indicated by a bold cross). Unfortunately we don’t have samples from the same treated people before the surgical operation, so that it is not possible to link the four cases directly to the pre-operatory phase. The dotted lines link the points before the filtering to the points after the attenuation of the noise-like features. We can force the algorithm to perform a stronger filtering of the dysphonic time series (as explained in the Appendix, this is possible increasing the number of iteration of the correction routine), but one has always to be careful: In order to get a better value of some components of the feature space, one has to come to a compromise with some other quantities and the best answer comes, as usual in these cases, from a direct listening of the time series (after the convertion in a *.wav* format).

The results shown in Fig.4 have been obtained with the following set of parameters: $`m=30`$, $`\nu =4`$ and $`ϵ=0.3\sigma `$, with $`\sigma `$ being the variance of the original data. Dysphonic patients, in the (spectral factor, pseudo-entropy) plane, are spread around the average point (5,1.8); the centre of mass of the healthy cluster is (35,0.5). The surgical treatment produces an average correction located in (20,0.6). Our software moves the centre of mass of sick patients to (15,0.7). In getting these results the value of $`m`$ and $`\nu `$ is not so crucial, provided that $`m\nu >l_p`$, being $`l_p`$ the extension of a pitch. The program is more sensitive to the choice of $`ϵ`$, in the following way: The length of the dotted line of Fig.4 is somehow proportional to the value of $`ϵ`$. Unfortunately the angle is not constant: This means that up to a certain value of $`ϵ`$ all the corrections act along the same direction, beyond that threshold they start to deteriorate the voice. In the extreme case of $`ϵ`$ as big as the full embedding space, the averaging procedure performed by Eq. 5 would destroy almost completely the signal (producing a pseudo-entropy and a spectral factor close to zero).

The sensitivity of the program to the choice of $`ϵ`$ is illustrated in Fig. 5. There we have filtered one sample with 8 different values of the neighbourood size. The original position is the one with the biggest value of the pseudo-entropy. $`ϵ`$ is increasing along the direction indicated by the arrow and only for two values of it the corrected point lies in the healthy region. For the last three corrections the neighbourood size was absolutely too big. The other samples behave in a qualitatively similar way.

## Conclusions

To summarize, we have proposed a procedure to correct numerically some kind of vocal disorders. The main ingredients are the theory of dynamical systems, the delay reconstruction from a scalar measurement, the method of recurrence plot and a nonlinear noise reduction scheme. We also make use of a feature space in order to visualize the results of the filtering procedure. A detailed description of how the idea has been implemented is also provided in Appendix. We have discussed the meaning, implications, methodology and mathematical background of the most relevant parameters and the sensitivity of the program to their choice. The software corrections perform an improvement on the voice quality that is comparable to what a surgical operation is able to do. This suggests the idea of implementing the procedure in a physical device able to help people correcting their voice, without having to undergo a medical treatment.

## Acknowledgements

The authors would like to thank for fruitful collaborations Prof. Hanspeter Herzel about modeling of the voice, Priv. Doz. Holger Kantz, Priv. Doz. Thomas Schreiber and Dr. Rainer Hegger about Chaos and Noise Reduction, Prof. Cecilia Salimbeni and the Phoniatric Section of the Careggi Hospital in Florence about Vocal Disorders.

## Appendix: Implementation outline

We present here the sketch of the subroutine *filter* written in a *pseudo-C* code. We skip all the details related to the syntax of the programming language for a better readability.

```

void show_options(char *progname)
{
  fprintf(stderr,"\t-m embedding dimension [Default: 5]\n");
  fprintf(stderr,"\t-T max distance in time [Default: no limit]\n");
  fprintf(stderr,"\t-d delay [Default: 1]\n");
  fprintf(stderr,"\t-r minimal neighbourhood size \n\t\t"
Ψ  "[Default: (interval of data)/1000]\n");
  fprintf(stderr,"\t-R maximal neighborhoodsize[Default: not set]\n");
  fprintf(stderr,"\t-i # of iterations [Default: 1]\n");
}

```

The input of the program has to be a scalar time series corresponding to the amplitude of the voice signal. Every conventional ASCII file is accepted with the data in one column. The parameter $`m`$ is the embedding dimension (typical values in the range $`[20,40]`$), $`T`$ is the portion of the time series where to look for neighbours: Since every pitch gets neighbours only within the same phoneme, the typical value of $`T`$ is 2000 (with a sampling rate of 22,05 kHz, 2000 points are almost 100 ms, approximately the length of a phoneme). The parameter $`d`$ is the time delay, usually in the range $`[3,10]`$ such that the product $`md`$ is as big as the length of a pitch. The size of the neighbourood $`ϵ`$ is specified through $`r`$ and $`R`$, its minimal and maximal value respectively: $`r`$ should equal the noise level and $`R`$ is usually $`1.5r`$. With the parameter $`i`$, number of iterations, one can set the strength of the filtering. Tipically we use $`i=1`$.

```

unsigned long lfind_neighbors(long act,double eps)
{
  k=(int)((searchdim-1)*delay);
  k1=(int)((searchdim/2)*delay);
  i=(int)(series[act-k]/eps)&ib;
  j=(int)(series[act-k1]/eps)&ib;
  n=(int)(series[act]/eps)&ib;

  for (i1=i-1;i1<=i+1;i1++) {
    i2=i1&ib;
    for (j1=j-1;j1<=j+1;j1++) {
      j2=j1&ib;
      for (n1=n-1;n1<=n+1;n1++) {
Ψelement=box[i2][j2][n1&ib];
Ψwhile (element != -1) {
Ψ  if (labs(act-element) < maxdist) {
Ψ    dx=0.0;
Ψ    for (k=0;k<searchdim;k++) {
Ψ      k1= k*(int)delay;
Ψ      dx += fabs(series[act-k1]-series[element-k1]);
Ψ    }
Ψ    if (dx/(double)searchdim <= eps) {
Ψ      dist[nf]=dx;
Ψ      flist[nf++]=element;
Ψ    }
Ψ  }
Ψ  element=list[element];
Ψ}
      }
    }
  }
  return nf;
}

```

Two points in the embedding space are considered neighbours if *(labs(act-element) $`<`$ maxdist)*. For every point we build here a list containing all its neighbours. We do it firstly in a two dimensional space to speed up the process. If two points are neighbours in a five dimensional space, this relation holds also in two dimension. Of course the opposite is not true.

```

void correct(unsigned long n)
{
  epsinv=1./eps;
  for (i=0;i<dim*delay;i++)
    hcor[i]=0.0;

  i=(int)(series[n-(dim-1)*delay]*epsinv)&ibox;
  j=(int)(series[n]*epsinv)&ibox;

  for (i1=i-1;i1<=i+1;i1++) {
    i2=i1&ibox;
    for (j1=j-1;j1<=j+1;j1++) {
      element=box[i2][j1&ibox];
      while (element != -1) {
Ψif (labs(n-element) < maxdist) {
Ψ  for (k=0;k<dim;k++) {
Ψ    k1=k*delay;
Ψ    dx=fabs(series[n-k1]-series[element-k1]);
Ψ    if (dx > eps)
Ψ      break;
Ψ  }
Ψ  if (k == dim) {
Ψ    flist[nfound++]=element;
Ψ    for (k=0;k<dim*delay;k++)
Ψ      hcor[k] += series[element-k];
Ψ  }
Ψ}
Ψelement=list[element];
      }
    }
  }
  for (k=0;k<dim*delay;k++) {
    corr[n-k] += (hcor[k]=series[n-k]-hcor[k]/nfound);
    nf[n-k]++;
  }
  for (i=0;i<nfound;i++) {
    j=flist[i];
    for (k=0;k<dim*delay;k++) {
      trend[j-k] += hcor[k];
      tcount[j-k]++;
    }
  }
}

```

The second to last *for* iteration performs exactly what indicated in Eq.5. The variable *nfound* is $`|U_ϵ(𝐬_{n_0})|`$.

```

int main(int argc,char **argv)
{
  series=(double*)get_series(infile,&length,exclude,column,1);

  rescale_data(series,length,&d_min,&d_max);

  resize_eps=0;
  for (iter=1;iter<=iterations;iter++) {
    epsilon=mineps;
    all_done=0;
    epscount=1;
    allfound=0;
    fprintf(stderr,"Starting iteration %d\n",iter);
    while(!all_done) {
      put_in_box(epsilon);
      all_done=1;
      for (n=(searchdim)*delay-1;n<length;n++)
Ψif (!ok[n]) {
Ψ  nfound=lfind_neighbors(n,epsilon);
Ψ  if (nfound >= minn) {
Ψ    correct(n);
Ψ    ok[n]=epscount;
Ψ    if (epscount == 1)
Ψ      resize_eps=1;
Ψ    allfound++;
Ψ  }
Ψ  else
Ψ    all_done=0;
Ψ}
      fprintf(stderr,"Corrected %ld points with epsilon= %e\n",allfound,
Ψ      epsilon*d_max);
      epsilon *= epsfac;
      epscount++;
      if (epsilon > maxeps)
Ψbreak;
    }

    sprintf(ofname,"%s.%d",outfile,iter);
    file=fopen(ofname,"w");
    fprintf(stderr,"Opened %s for writing\n\n",ofname);
    for (i=0;i<length;i++) {
      fprintf(file,"%e\n",series[i]*d_max+d_min);
      if (stdo && (iter == iterations))
Ψfprintf(stdout,"%e\n",series[i]*d_max+d_min);
    }
    fclose(file);
  }
  return 0;
}

```

The flow is the following: Get the time series, rescale it to a proper interval, apply the filter (call to the function *correct*) a *iter* number of times, print the corrected time series to a new data file.


text:
# 1 Introduction

## 1 Introduction

Gauge hierarchy problem is one of the deep problems in the Standard Model. About three years ago, it was suggested that the large compactified extra dimensions may also be the solution to the gauge hierarchy problem , because a low ($`4+n`$)-dimensional Planck scale ($`M_X`$) may result in the large 4-dimensional Planck scale ($`M_{pl}`$) due to the large physical volume ($`V_p^n`$) of extra dimensions: $`M_{pl}^2=M_X^{2+n}V_p^n`$. In addition, about two years ago, Randall and Sundrum  proposed another scenario that the extra dimension is an orbifold, and the size of extra dimension is not large but the 4-dimensional mass scale in the Standard Model is suppressed by an exponential factor from 5-dimensional mass scale because of the exponential warp factor. Furthermore, they suggested that the fifth dimension might be coordinate non-compact , and there may exist only one brane with positive tension at origin, but, there exists the gauge hierarchy problem and the physical size of the fifth dimension is indeed compact. The remarkable aspect of the second scenario is that it gives rise to a localized graviton field. Recently, a lot of 5-dimensional models with 3-branes were built \[4-5\]. We constructed the general models with parallel 3-branes on the five-dimensional space-time, and obtained that the 5-dimensional GUT scale on each brane can be identified as the 5-dimensional Planck scale, but the 4-dimensional Planck scale is generated from the low 4-dimensional GUT scale exponentially in our world. Furthermore, the models with codimension-1 brane(s) were constructed on the six-dimensional and higher dimensional space-time \[6-9\].

On the other hand, it is well known that five kinds of superstring theories live in the 10-dimensional space-time manifold and M-theory lives in the 11-dimensional space-time manifold due to the anomaly cancellation. In short, it seems that the low energy phenomenology and the current fundamental theories favour the existence of extra dimensions. Therefore, we might want to ask the question: why are the extra dimensions physical compact at low energy scale? From our philosophy, in the early Universe or when the energy scale or temperature is very high, all the space-time dimensions might be phyical non-compact, furthermore, all the space dimensions might be equivalent, i. e., there may be no difference between the 3-dimensional space we observe and the extra space dimensions. If this philosophy was right, then, we should understand how the extra dimensions become physical compact at low energy scale.

We conjecture that at high energy scale or high temperature, the extra dimensions are physical non-compact. After the cosmological phase transition, or the global or gauge symmetry breaking, or the supersymmetry breaking in supergravity theory, the bulk cosmological constant becomes negative, and then, the extra dimensions become physical compact.

In this paper, in order to support our conjecture, we give a 5-dimensional toy brane model with three paprallel 3-branes and one real bulk scalar which has temperature dependent potential. We assume that at high energy scale or temperature, the bulk cosmological constant is zero, because in string theory, it is natural to take the bulk cosmological constant to be zero since the tree-level vacuum energy in the generic critical closed string compactifications (supersymmetric or not) vanishes, and the zero bulk cosmological constant is natural in the scenario in which the bulk is supersymmetric (though the brane need not be), or the quantum corrections to the bulk are small enough to be neglect in a controlled expansion. For simplicity, we also assume that at high energy scale, the metric is flat and diagonal. When the energy scale or temperature goes down, there is a second order cosmological phase transition. After the cosmological phase transition, the bulk cosmological constant becomes negative, and then, the metric has warp factor, which makes the fifth dimension physical compact. In fact, at low energy scale, this model is the same model discussed previous in . Of course, this idea can be generalized to all the previous brane models and brane networks with constant bulk cosmological constant \[2-9\].

The key point to have the spontaneous physical compactification of the extra dimensions is that, after the cosmological phase transition or symmetry breaking, the bulk cosmological constant becomes negative, and then, the extra dimensions become physical compact due to the warp factor. Therefore, after the global or gauge symmetry breaking, or the supersymmetry breaking in supergravity theory, the spontaneous physical compactification of extra dimensions might be happened for the cosmological constant might become negative. And it is interesting to discuss the spontaneous physical compactification of the extra dimensions in the model buildings, and in the string theory or M-theory compactifications.

## 2 One Toy Model

In this section, we would like to present a toy brane model with spontaneous physical compactification of the extra dimension.

We consider the space-time is five-dimensional, and the fifth dimension is $`R^1`$. In addition, in order to obtain the cosmological phase transition, we introduce a real scalar $`\varphi `$ whose bulk potential is temperature dependent. Assuming that we have three parallel 3-branes along the fifth dimension, and their fifth coordinates are: $`\mathrm{}<y_1<y_2<y_3<+\mathrm{}`$. We obtain the metric in each brane from the five-dimensional metric $`g_{AB}`$ where $`A,B=\mu ,y`$ by restriction

$$g_{\mu \nu }^{(i)}(x^\mu )g_{\mu \nu }(x^\mu ,y=y_i).$$
(1)
In this paper, for simplicity, we assume that the metric is diagonal, and at high energy scale or temperature, the metric is $`diag(1,1,1,1,1)`$.

The classical action with a bulk scalar is given by

$`S`$ $`=`$ $`S_{\mathrm{Bulk}}+S_{\mathrm{Brane}},`$ (2)
$`S_{\mathrm{Bulk}}`$ $`=`$ $`{\displaystyle d^4x𝑑y\sqrt{g}\{\frac{1}{2}M_X^3R\frac{1}{2}_A\varphi ^A\varphi \mathrm{\Lambda }(\varphi ,T)\}},`$ (3)
$`S_{\mathrm{Brane}}`$ $`=`$ $`{\displaystyle \underset{i=1}{\overset{3}{}}}{\displaystyle d^4x𝑑y\sqrt{g^{(i)}}\{_iV_i(\varphi ,T)\}\delta (yy_i)},`$ (4)
where $`M_X`$ is the 5-dimensional Planck scale, $`T`$ is the temperature, $`\mathrm{\Lambda }(\varphi ,T)`$ is the temperature dependent bulk potential for $`\varphi `$ , and $`V_i(\varphi ,T)`$ where $`i=1,2,3`$ is the brane tension which is dependent on $`\varphi `$ and $`T`$. Our ansatz for $`\mathrm{\Lambda }(\varphi ,T)`$ and $`V_i(\varphi ,T)`$ is222There might exist the terms proportional to $`T^5`$ and $`T^4`$ in the $`\mathrm{\Lambda }(\varphi ,T)`$ and $`V_i(\varphi ,T)`$, respectively, which will vanish at low temperature. Those additional terms might not change the compactness of the space, but, will make the discussions of the spontaneous physical compactification of extra dimension a little bit complicated. So, as a toy model, we do not consider those terms and do not consider how to derive the ansatz Eqs. (5) and (6) here.

$`\mathrm{\Lambda }(\varphi ,T)`$ $`=`$ $`D(T^2T_o^2)\varphi ^2+{\displaystyle \frac{\lambda (T)}{4M_X}}\varphi ^4,`$ (5)
$`V_i(\varphi ,T)`$ $`=`$ $`C_i(T)M_X\varphi ^2,`$ (6)
where $`D`$ and $`T_o`$ are constant terms, $`\lambda (T)`$ is a slowly varying function of $`T`$. The $`C_i(T)`$ is also a slowly varying function and will be determined later. The temperature dependent potential for $`\varphi `$ has $`Z_2`$ symmery, which is invariant under $`\varphi \varphi `$.

The solutions to the minimum or maximum of the potential $`\mathrm{\Lambda }(\varphi ,T)`$ are given by

$$<\varphi >|_T=0,$$
(7)
$$<\varphi >|_T=\sqrt{\frac{2D(T_o^2T^2)M_X}{\lambda (T)}}.$$
(8)
Therefore, the critical temperatue is given by $`T_o`$. At $`T>T_o`$, the origin $`<\varphi >=0`$ is the minimum. At the same time, only the solution $`<\varphi >=0`$ does exist. In this case, the bulk cosmological constant is zero, and all the brane tensions are zero or one can think that the brane tensions are very small. At $`T=T_o`$, both solutions collapse at $`<\varphi >=0`$. At $`T<T_o`$, the origin $`<\varphi >=0`$ is the maximum, and the solution $`<\varphi >=\sqrt{\frac{2D(T_o^2T^2)M_X}{\lambda \left(T\right)}}`$ gives the minimum, so, the $`Z_2`$ symmetry is broken. This cosmological phase transition is called of second order because there is no barrier during the cosmological phase transition.

When $`T<T_o`$, we can expand the $`\varphi `$ around its minimum, and write $`\varphi `$ as

$`\varphi `$ $`=`$ $`\varphi ^{}+\sqrt{{\displaystyle \frac{2D(T_o^2T^2)M_X}{\lambda (T)}}},`$ (9)
so, we can express the $`\mathrm{\Lambda }(\varphi ,T)`$ and $`V_i(\varphi ,T)`$ as $`\mathrm{\Lambda }(\varphi ^{},T)`$ and $`V_i(\varphi ^{},T)`$, respectively

$`\mathrm{\Lambda }(\varphi ^{},T)`$ $`=`$ $`{\displaystyle \frac{\lambda (T)}{4M_X}}\left(\varphi ^2+2\varphi ^{}\sqrt{{\displaystyle \frac{2D(T_o^2T^2)M_X}{\lambda (T)}}}\right)^2`$ (10)
$`{\displaystyle \frac{D^2M_X}{\lambda (T)}}(T_o^2T^2)^2,`$
$`V_i(\varphi ^{},T)=C_i(T)M_X\left(\varphi ^2+2\varphi ^{}\sqrt{{\displaystyle \frac{2D(T_o^2T^2)M_X}{\lambda (T)}}}+{\displaystyle \frac{2D(T_o^2T^2)M_X}{\lambda (T)}}\right).`$ (11)
Neglecting the field $`\varphi ^{}`$ for $`<\varphi ^{}>=0`$ , during the cosmological phase transition, we obtain the dominant classical action

$`S_{\mathrm{Bulk}}={\displaystyle d^4x𝑑y\sqrt{g}\left(\frac{1}{2}M_X^3R+\frac{D^2M_X}{\lambda (T)}(T_o^2T^2)^2\right)},`$ (12)
$`S_{\mathrm{Brane}}={\displaystyle \underset{i=1}{\overset{3}{}}}{\displaystyle d^4x𝑑y\sqrt{g^{(i)}}\left(_i\frac{2C_i(T)D(T_o^2T^2)M_X^2}{\lambda (T)}\right)\delta (yy_i)}.`$ (13)
Because the calculations are similar to those in Ref. , we just give the result here. Assuming that there exists a solution that respects 4-dimensional Poincare invariance in the $`x^\mu `$-directions, one obtains the 5-dimensional metric:

$`ds^2=e^{2\sigma (y,T)}\eta _{\mu \nu }dx^\mu dx^\nu +dy^2.`$ (14)
With this metric, the Einstein equations reduce to

$`\left({\displaystyle \frac{d\sigma }{dy}}\right)^2`$ $`=`$ $`{\displaystyle \frac{D^2(T_o^2T^2)^2}{6M_X^2\lambda (T)}},`$ (15)
$`{\displaystyle \frac{d^2\sigma }{d^2y}}`$ $`=`$ $`{\displaystyle \underset{i=1}{\overset{3}{}}}{\displaystyle \frac{2C_i(T)D(T_o^2T^2)}{3M_X\lambda (T)}}\delta (yy_i).`$ (16)
The general solution to above differential equations is

$$\sigma (y,T)=\underset{i=1}{\overset{3}{}}(1)^{i+1}k(T)|yy_i|+f[T_o^2T^2],$$
(17)
where $`f[T_o^2T^2]`$ is the polynomial with variable $`T_o^2T^2`$ and without constant term, and $`f[T_o^2T^2]`$ is similar to the constant $`c`$ in .

From Eqs. (15-17), we obtain

$`k(T)=\sqrt{{\displaystyle \frac{1}{6\lambda (T)}}}{\displaystyle \frac{D(T_o^2T^2)}{M_X}},`$ (18)
$`C_i(T)=(1)^{i+1}\sqrt{{\displaystyle \frac{3\lambda (T)}{2}}}.`$ (19)
In addition, for simplicity, we define that $`\sigma (y)\sigma (y,T=0)`$, $`\lambda _o\lambda (T=0)`$, $`C_i^oC_i(T=0)`$, and $`k_ok(T=0)`$. After the cosmological phase transition, i. e., at very low temperature or energy scale, we obtain the metric

$$\sigma (y)=\underset{i=1}{\overset{3}{}}(1)^{i+1}k_o|yy_i|+f[T_o^2],$$
(20)
where

$`k_o=\sqrt{{\displaystyle \frac{1}{6\lambda _o}}}{\displaystyle \frac{DT_o^2}{M_X}}.`$ (21)
And we have

$`C_i^o=(1)^{i+1}\sqrt{{\displaystyle \frac{3\lambda _o}{2}}}.`$ (22)
From above metric, we obtain that, after the cosmological phase transition, the fifth dimension becomes physical compact because of the warp factor. And the gauge hierarchy problem can be solved in this model, as discussed in Ref. .

One might notice that in $`\mathrm{\Lambda }(\varphi ^{},T)`$, the lowest order $`\varphi ^{}`$ term is order of $`\varphi ^2`$, but in $`V_i(\varphi ^{},T)`$, the lowest order $`\varphi ^{}`$ term is order of $`\varphi ^{}`$. In fact, the lowest order $`\varphi ^{}`$ term in $`V_i(\varphi ^{},T)`$ can be order of $`\varphi ^2`$, for example, if we defined the original $`V_i(\varphi ,T)`$ as

$`V_i(\varphi ,T)`$ $`=`$ $`C_i(T)M_X\left(\varphi ^22\theta (T_oT)\sqrt{{\displaystyle \frac{2D(T_o^2T^2)M_X}{\lambda (T)}}}\varphi \right),`$ (23)
or

$`V_i(\varphi ,T)`$ $`=`$ $`C_i(T)M_X\left(\varphi ^22\sqrt{{\displaystyle \frac{2D|T_o^2T^2|M_X}{\lambda (T)}}}\varphi \right),`$ (24)
where $`\theta (x)`$ is the step function. The definition in Eq. (23) preserves the $`Z_2`$ symmetry before the cosmological phase transition. With definition of $`V_i(\varphi ,T)`$ in Eq. (23) or Eq. (24), the discussions of the spontaneous physical compactification of extra dimension are similar to the above discussions if one made the following transformation: $`C_i(T)C_i(T)`$, and $`C_i^oC_i^o`$, so, we will not repeat it here.

One can also discuss the first order cosmological phase transition if one introduced the term $`E(T)\sqrt{T}\varphi ^3`$ in the bulk potential, however, in this case, there exists barrier during the cosmological phase transition .

In addition, this idea can be generalized to all the five-dimensional brane models with parallel 3-branes and the brane networks discussed in \[2-9\].

Furthermore, if one considered a compex scalar in the bulk, the global symmetry will be $`U(1)`$. Global symmetries are argued to be broken by quantum gravity effects . So, one might need to consider it as a gauge $`U(1)`$ symmetry and put additional particles in the bulk. Similarly, the spontaneous physical compactification of extra dimension might also be happened in this scenario.

## 3 Discussion and Conclusion

The key point in above model to have the spontaneous physical compactification of extra dimension is that, after the cosmological phase transition or symmetry breaking, the bulk cosmological constant becomes negative, and then, the extra dimension becomes physical compact due to the warp factor. Therefore, in order to search other scenarios where the spontaneous physical compactification of extra dimensions may be realized, we need to explore other possibilities to have negative vacuum energy after the symmetry breaking.

Obviously, the global or gauge symmetry breaking (Higgs mechanism) is one candidate, because after the symmetry breaking, the cosmological constant will be negative if the cosmological constant is zero in original theory. For example, we can put the gauge group and some particles in the bulk, the gauge symmetry might be broken by radiative corrections at low energy scale (similar to the electroweak symmetry breaking in the mSUGRA model), then, the bulk cosmological constant might become negative and the spontaneous physical compactification of extra dimensions might be happened.

The other possibility is the supersymmetry breaking. As we know, in the supergravity theory, the gravitino mass term is negative in the potential. And if the supersymmetry was broken, the gravitino mass will give a negative contribution to the cosmological constant. So, the cosmological constant might be negative after the supersymmetry breaking, and then, we might have the spontaneous physical compactification of extra dimensions . Of course, at high energy scale, the supersymmetry is preserved, so, the cosmological constant is zero and the extra dimensions are physical non-compact.

How to realize the spontaneous physical compactification of extra dimensions due to the global or gauge symmetry breaking, or the supersymmetry breaking in the model buildings and in the string theory or M-theory compactifications is an interesting subject and deserves further study.

In short, we conjecture that at high energy scale or high temperature, the extra dimensions are physical non-compact. After the cosmological phase transition, or the global or gauge symmetry breaking, or the supersymmetry breaking in supergravity theory, the bulk cosmological constant becomes negative, and then, the extra dimensions become physical compact. We show this in a five-dimensional toy brane model with three parallel 3-branes and a real scalar whose bulk potential is temperature dependent. In the mean time, the gauge hierarchy problem can be solved in the toy model.

## Acknowledgments

This research was supported in part by the U.S. Department of Energy under Grant No. DOE-EY-76-02-3071.


text:
# Acknowledgements

## Acknowledgements

The work is supported in part by the Russian Fund for Basic Research under contracts 99-01-00091 and 01-02-16171. A.A. Pivovarov is Alexander von Humboldt fellow.