text:
# Proposition 1 ([]).

Pattern avoidance in circular permutations
DAVID CALLAN

Department of Statistics

University of Wisconsin-Madison

1210 W. Dayton St

Madison, WI 53706-1693

callan@stat.wisc.edu

In the present context, a *circular permutation* is an arrangement of distinct positive integers (called letters) clockwise around a circle, two arrangements being considered the same if they differ only by a rotation. A (linear) permutation is a list of distinct positive integers and its *reduced form* is obtained by replacing its smallest entry by 1, its next smallest by 2, and so on. A *pattern* is any such reduced form. An occurrence of a pattern $`\tau `$ in a cyclic permutation $`\pi `$ is a sequence of letters in clockwise order in $`\pi `$ (not extending more than one revolution) whose reduced form is $`\tau `$. We will always represent a circular permutations on $`[n]:=\{1,2,\mathrm{},n\}`$ as a linear permutation that ends in $`n`$. For example, $`\mathrm{5\hspace{0.17em}6\hspace{0.17em}4\hspace{0.17em}2\hspace{0.17em}3\hspace{0.17em}1\hspace{0.17em}7}`$ has just one occurrence of the pattern 1234, the letters being $`\mathrm{2\hspace{0.17em}3\hspace{0.17em}5\hspace{0.17em}6}`$ (keep in mind we may wrap around from the end to the start).

Much work has been done counting permutations with restrictions on the patterns they contain and possibly other restrictions as well; see  for a survey. There are two 3-letter patterns up to cyclic rotation: 123 and 132. Only the identity circular permutation avoids 132 and only the reverse identity avoids 123. Here we find simple explicit formulas for the number of circular permutations avoiding a single four letter pattern. There are 6 such patterns up to cyclic rotation and “reversing” cuts the number to 3: 1234, 1324 and 1342 are representative. In each case, an analysis of the structure according to the position of $`n1`$ (recall $`n`$ is at the end) yields a recurrence relation. For completeness we include a proof of the following result.

###### Proposition 1 ().

There are $`2^{n1}`$ $`(`$linear$`)`$ permutations on $`[n]`$ that avoid both $`213`$ and $`231`$.

###### Proof.

Consider such a permutation $`\pi `$ and any letter $`a`$ in $`\pi `$. The letters following $`a`$ in $`\pi `$ must either all be $`>a`$ or all be $`<a`$, else $`a`$ would be the first letter of a proscribed pattern. So the the first letter of $`\pi `$ is 1 or $`n`$ and we have two choices for each successive letter except the last. This observation translates into a bijection from 0-1 sequences of length $`n1`$ to permutations on $`[n]`$ that avoid 213 and 231: given $`(b_i)_{i=1}^{n1},b_i=0`$ or 1, form $`\pi `$ left to right as follows. For $`i=1`$ to $`n1`$ in turn, if $`b_i=1`$ (resp. 0), set $`\pi _i=`$ the largest (resp. smallest) letter in $`[n]`$ that has not already been used in $`\pi `$. For example, with $`n=8,(0,1,1,1,0,1,0)(1,8,7,6,2,5,3,4)`$. ∎

###### Theorem 1.

The number of $`1324`$-avoiding circular permutations on $`[`$n$`]`$ is the Fibonacci number $`F_{2n3}`$

###### Proof.

Let $`u_n`$ denote the desired number. Then $`u_1=u_2=1,u_3=2`$. For $`n4`$, we count by the location $`k[n1]`$ of $`n1`$ in $`\pi `$. First suppose $`kn2`$. so that $`n1`$ and $`n`$ are not adjacent. The first $`k1`$ letters in $`\pi `$ must all be $`>`$ than all the letters between $`n1`$ and $`n`$, else there would be a 1324 pattern of the form $`an1bn`$. So the first $`k`$ entries are $`\{nj\}_{j=1}^k`$. Furthermore, they must occur in increasing order, else there would be a $`cbn1a`$ pattern, rotationally equivalent to a $`1324`$. We leave the reader to verify that any such circular permutation is in fact 1324-avoiding. It is now easy to check that for $`kn2`$ (including $`k=1`$), deleting the first $`k`$ entries and reducing (in this case, reducing means replacing $`n`$ by $`nk`$) is a bijection to 1324-avoiding circular permutations on $`[nk]`$, counted by $`u_{nk}`$. Now suppose $`k=n1`$. Deleting the $`n`$ gives a bijection to the $`n1`$ case, counted by $`u_{n1}`$. Summing over $`k`$, we get a recurrence relation: $`u_n=u_{n1}+_{k=1}^{n1}u_{nk}n3`$ with initial condition $`u_1=u_2=1`$, whose solution is $`u_n=F_{2n3}`$.

Here is a short combinatorial proof that $`u_n:=F_{2n3}`$ satisfies the recurrence. Let $`_n`$ denote the set of 0-1 sequences of length $`n1`$ that start with 0 (unless they’re empty) and do not contain two consecutive 1’s, so that $`|_n|=F_n`$. Set $`𝒢_n=_{2n3}`$ so that $`|𝒢_n|=u_n`$. If $`w𝒢_n`$ begins 01…, deleting this 01 is a bijection to $`𝒢_{n1}`$. Otherwise $`w`$ begins $`00\mathrm{}`$ and we count by number $`k`$ of consecutive 10 pairs immediately following the initial 00. Clearly, $`k[0,n3]`$ and deleting the initial 00 and all $`k`$ of these 10’s is a bijection to $`𝒢_{n1k}`$, counted by $`u_{n1k}`$. ∎

###### Theorem 2.

There are $`2^{n1}(n1)`$ $`1342`$-avoiding circular permutations on $`[n]`$.

###### Proof.

Let $`u_n`$ denote the number of 1342-avoiding circular permutations on $`[n]`$. Once again $`u_1=u_2=1,u_3=2`$, so let $`n4`$. If $`n1`$ is in first position, deleting $`n`$ (and rotating) gives a bijection to 1342-avoiding circular permutations on $`[n1]`$, counted by $`u_{n1}`$.

Now suppose $`n1`$ is in position $`k`$ with $`k[2,n2].`$ Write $`\pi =w_1n1\mathrm{}n`$ with $`w_1`$ the initial subword of length $`k1`$. Then $`w_1`$ is nonempty and its letters are increasing, else $`n1,n`$ would be the middle letters of a 1342 pattern. Also, the letters in $`w_1`$ form an interval of consecutive integers. To see this, suppose $`a<c`$ are letters of $`w_1`$ and $`b[a+1,c1]`$. Then $`b`$ can’t precede $`a`$ in $`\pi `$ or fall between $`c`$ and $`n1`$ because the letters in $`\pi `$ increase. Neither can $`c`$ fall between $`n1`$ and $`n`$, else $`acn1b`$ would be a $`1342`$. So $`b`$ must fall between $`a`$ and $`c`$ and $`w_1`$ is a list of consecutive integers.

It follows that $`w_1`$ splits the remaining letters in $`[n2]`$ into two sets $`A,B`$ (either or both may be empty, $`A`$ denotes the set of smaller letters). Write $`\pi `$ as $`w_1n1w_0n`$ so that $`w_0`$ is comprised of the letters in $`AB`$. The letters in $`A`$ are in increasing order in $`\pi `$ because if $`w_0=\mathrm{}a_2\mathrm{}a_1\mathrm{}`$ with $`a_1<a_2`$, then $`a_1bn1a_2`$ would be a 1342 for any $`bw_1`$. All letters of $`C`$ precede all letters of $`A`$ because if $`w_0=\mathrm{}a\mathrm{}c\mathrm{}`$ with $`aA,cC`$, then $`acnb`$ would be a 1342 for any $`bw_1`$.

So far, then, we have shown that $`\pi `$ has the form

$$w_1n1w_2w_3n$$

with $`w_1`$ of length $`k1`$ (and so nonempty), all letters in $`w_2>`$ those in $`w_1`$, all letters in $`w_3`$ increasing and $`<`$ those in $`w_1`$. Now we claim further that the word $`w_2`$ must avoid the patterns 213 and 231. Indeed, if $`w_2=\mathrm{}b\mathrm{}a\mathrm{}c\mathrm{}`$ with $`a<b<c`$, then $`acnb`$ is a 1342 in $`\pi `$. And if $`w_2=\mathrm{}b\mathrm{}c\mathrm{}a\mathrm{}`$, then $`xbca`$ is a 1342 for any $`x`$ in $`w_1`$. We leave the reader to verify that any circular permutation with $`n1`$ in position $`k[2,n2]`$ that meets all these specifications is in fact 1342-avoiding. Let us count these latter permutations. The initial segment $`w_1`$ is determined by its first letter $`j`$ which must lie in $`[nk]`$. Then $`w_3`$ is fully determined as are the letters of $`w_2`$. Since $`w_2`$ avoids 213 and 231 (Prop. 1), there are $`2^{|w_2|1}=2^{nkj1}`$ possibilities for $`w_2`$ unless $`j=nk`$ in which case $`w_2`$ is empty and there is just one possibility. Summing over $`j`$, we have $`1+_{j=1}^{nk1}2^{nkj1}=2^{nk1}`$ possibilities in all. Now summing over $`k`$, we get the defining recurrence $`u_1=1,u_2=1,u_n=u_{n1}+2^{n2}1`$ for $`n3`$, with solution $`u_n=2^{n1}(n1)`$ valid for all $`n1`$.

The preceding analysis could be worked up into a bijection. Briefly, let $`𝒜_n`$ denote the set of $`n`$-bit 0-1 sequences that do not contain precisely one 1 so that $`|𝒜_n|=2^nn`$. Given $`u𝒜_{n1}`$ build up $`\pi =w_1n1w_2w_3n`$ as follows. If $`u_1=0`$, put $`n1`$ in position 1 and proceed recursively. Otherwise, $`u_1=1`$ and put $`n1`$ in the position $`k2`$ of the next 1 in $`u`$ (there will be one!). Then let $`j[nk]`$ be the position of the next 1 after that in $`u_{k+1}\mathrm{}u_{n1}`$ (with $`j=nk`$ if there are no more 1’s). Use $`j`$ as the initial letter that determines $`w_1`$ and use $`u_{k+j+1}\mathrm{}u_{n1}`$ to construct $`w_2`$ (Prop. 1). ∎

###### Theorem 3.

There are $`2^n+12n\left(\genfrac{}{}{0pt}{}{n}{3}\right)`$ $`1234`$-avoiding circular permutations on $`[n]`$.

###### Proof.

Again let $`u_n`$ denote the desired number and $`k`$ the position of $`n1`$. If $`k=1`$, deleting the $`n1`$ and reducing is a bijection to the $`n1`$ case. If $`k=2`$, write $`\pi `$ as $`an1w_1n`$ with $`a[n2]`$. All letters $`>a`$ in $`w_1`$ clearly occur in decreasing order, as must all letters $`>a`$ (else $`n1`$ would terminate a 1234). So the only freedom is in placement of the letters $`<a`$, giving $`_{a=1}^{n2}\left(\genfrac{}{}{0pt}{}{n3}{a1}\right)=2^{n3}`$ possibilities with $`k=2`$.

Now suppose $`k3`$ so that $`\pi =w_1n1w_2n`$ with $`|w_1|=k12`$. Let $`m`$ and $`M`$ denote the minimum and maximum letters in $`w_1`$ and let $`A`$ consist of the remaining $`k3`$ letters in $`w_1`$. Then the letters of $`w_2`$ can be partitioned into three sets: $`B=[1,m1],C=[m+1,M1]\backslash A,D=[M+1,n2]`$ of sizes $`|B|=m1,|C|=Mm1K,|D|=nM2`$ respectively, where $`K=|A|=k3`$. All letters $`<M`$ in $`w_2`$ must occur in decreasing order, else there would be a 1234 terminating in $`Mn1`$; likewise for all letters $`>m`$, else $`muvn`$ would be a 1234 for some $`u,vw_2`$. This means that unless $`M=m+1`$, *all* letters in $`w_2`$ must decrease and in case $`M=m+1`$ (and hence $`C`$ is empty) the letters of $`B`$ and $`D`$ separately must occur in decreasing order in $`w_2`$ giving $`\left(\genfrac{}{}{0pt}{}{n+mM3}{m1}\right)`$ choices to place the letters of $`B`$ and $`D`$ in $`w_2`$. Again we leave the reader to verify that if all these conditions are met, the proscribed pattern is avoided.

Using standard binomial coefficient identities, we can now count possibilities. Fix $`k[3,n1]`$ (so $`K=k3`$ is also fixed) and sum over $`m`$ and $`M`$. There are $`\left(\genfrac{}{}{0pt}{}{Mm1}{K}\right)`$ ways to select the set $`A`$. For $`C=\mathrm{}`$ (in which case $`M=m+1+K`$), the count is $`_{m=1}^{nK3}\left(\genfrac{}{}{0pt}{}{nK4}{m1}\right)=2^{nK4}=2^{n1k}`$. For $`C\mathrm{}`$, the count is

$`{\displaystyle \underset{m=1}{\overset{nK4}{}}}{\displaystyle \underset{M=m+K+2}{\overset{n2}{}}}\left({\displaystyle \genfrac{}{}{0pt}{}{Mm1}{K}}\right)`$ $`=`$ $`{\displaystyle \underset{m=1}{\overset{nK4}{}}}\left(\left({\displaystyle \genfrac{}{}{0pt}{}{nm2}{K+1}}\right)1\right)`$
$`=`$ $`\left({\displaystyle \genfrac{}{}{0pt}{}{n2}{K+2}}\right)n+K+3`$
$`=`$ $`\left({\displaystyle \genfrac{}{}{0pt}{}{n2}{k1}}\right)n+k`$
Taking the cases $`C=\mathrm{}`$ and $`C\mathrm{}`$ together, the total count for fixed $`k[3,n1]`$ is $`2^{n1k}+\left(\genfrac{}{}{0pt}{}{n2}{k1}\right)n+k`$, while, as shown above, for $`k=2`$ it is $`2^{n3}`$ and for $`k=1`$ it is $`u_{n1}`$. Summing over $`k`$ now yields the recurrence relation $`u_1=u_2=1,u_n=u_{n1}+2^{n1}n\left(\genfrac{}{}{0pt}{}{n2}{2}\right)`$ for $`n3`$, with solution $`u_n=2^n+12n\left(\genfrac{}{}{0pt}{}{n}{3}\right)`$. ∎

We conclude by noting that the preceding exact enumerations yield growth constants ($`lim_n\mathrm{}u_n^{1/n}`$) for the patterns 1324, 1342 and 1234 of $`\varphi ^22.618,2`$ and 2 respectively, where $`\varphi `$ denotes the golden ratio.


format: text

text:
# 1. Introduction.

## 1. Introduction.

Dielectric relaxation in solids, a process commonly defined as an approach to equilibrium of a dipolar system driven out of equilibrium by a step or alternating external electric field, represents one of the most intensively researched topics in experimental and theoretical physics, see e.g. \[1-21\]. Wide-ranging experimental information resulting from the studies of the basic dielectric characteristics (i.e., the time decay of the depolarization current $`i(t)`$ and the frequency-dependent complex dielectric permittivity $`ϵ(\omega )`$ or susceptibility $`\chi (\omega )ϵ(\omega )ϵ_{\mathrm{}}`$, where $`ϵ_{\mathrm{}}`$ is the asymptotic value of the dielectric permittivity $`ϵ(\omega )`$ at high frequencies) has led to the conclusion that the classical phenomenology of relaxation breaks down in complex materials. It has been found that the Debye behaviour, represented by the exponentially decaying relaxation function

$$\varphi (t)=\mathrm{exp}(\omega _pt)$$
(1)
yielding

$$i(t)f(t)=\frac{d\varphi (t)}{dt}=\omega _p\mathrm{exp}\left(\omega _pt\right)$$

and

$$\chi (\omega )\varphi ^{}(\omega )=\underset{0}{\overset{\mathrm{}}{}}e^{i\omega t}f(t)𝑑t=\frac{1}{1+i\omega /\omega _p}$$
(2)
(where the constant $`\omega _p`$ denotes the loss peak frequency) is hardly ever found in nature and that the deviations from it for many dielectrics may be relatively large . A class of systems exhibiting the non-Debye relaxation patterns includes various complex materials such as supercooled liquids, amorphous semiconductors and insulators, polymers, disordered crystals, molecular solid solutions, glasses, etc.

For a long time effort was being diverted to a purely qualitative representation of the shape of the non-Debye dielectric functions $`\varphi (t)`$ or $`\varphi ^{}(\omega )`$ in terms of certain mathematical expressions without in any way going into the physical significance of these representations. It turns out that all dielectric data are characterised well enough by a few empirical functions . The time-domain relaxation data usually are fitted by means of the stretched exponential relaxation function

$$\varphi (t)=\mathrm{exp}\left((\omega _pt)^\alpha \right)\mathrm{with}\mathrm{\hspace{0.33em}\hspace{0.33em}0}<\alpha <1$$
(3)
known as the Kohlrausch-Williams-Watts (KWW) function and coincident with the Debye case (1) for $`\alpha =1`$. The most popular analytical expression applied to the complex susceptibility data is given by the Havriliak-Negami (HN) function

$$\varphi ^{}(\omega )=\frac{1}{\left(1+(i\omega /\omega _p)^\alpha \right)^\gamma }$$
(4)
where $`0<\alpha ,\gamma <1`$. For $`\alpha =1`$ and $`\gamma <1`$, formula (4) takes the form known as the Cole-Davidson (CD) function; for $`\gamma =1`$ and $`\alpha <1`$ it takes the form of the Cole-Cole (CC) function, and for $`\alpha =1`$ and $`\gamma =1`$ one obtains the classical Debye form (2).

We repeat that these various model functions are essentially attempts to characterise the observed behaviour without in any way indicating the physical mechanisms involved. To that extent, we note that the problem of understanding the nature of dielectric relaxation is as yet largely open since we do not have a sufficiently general method of approach to the relaxation processes. At the same time we note certain fundamental features of relaxation response which make it clear that we are faced with some very general processes into which we require a much deeper insight.

The HN, CC, CD, and KWW model functions all show the high-frequency limit of the response characterised by the fractional power law :

$$\chi ^{}(\omega )=\mathrm{tan}\left(\frac{n\pi }{2}\right)\chi ^{\prime \prime }(\omega )\omega ^{n1}\text{for some}\mathrm{\hspace{0.33em}\hspace{0.33em}0}<n<1$$
(5)
which has the straightforward consequence of the constant ratio of the imaginary to the real components of the dielectric susceptibility $`\chi (\omega )=\chi ^{}(\omega )i\chi ^{\prime \prime }(\omega )`$:

$$\frac{\chi ^{\prime \prime }(\omega )}{\chi ^{}(\omega )}=\mathrm{cot}\left(n\frac{\pi }{2}\right)\mathrm{for}\omega \omega _p.$$
(6)
Let us note that this relation implies that the ratio of energy lost per radian to the energy stored at the peak is independent of frequency. This very simple and intuitively evident property is characteristic of the prevailing universal response of most dielectric materials. We call this the “energy criterion” and it is neither known nor predicted by any of the accepted theories of relaxation.

Important for our discussion of relaxation phenomena is the concept of residual loss  which remains after the removal of such strongly frequency-dependent processes as the direct-current contribution going as $`\omega ^1`$, or low-frequency dispersion where loss is proportional to $`\omega ^{n1}`$ with $`n0`$, or dipolar loss which shows pronounced peaks. The residue left is broadly frequency-independent or only weakly dependent over extended ranges of frequency which may be approximated by the limit $`n1`$. Such a limit corresponds to frequency-independent components of the dielectric susceptibility, $`\chi ^{}(\omega )`$ and $`\chi ^{\prime \prime }(\omega )`$, with their ratio tending to infinity which is impossible practically but can be seen experimentally in various approximations. Figure 1 shows a compilation of data for different materials taken from much wider sets for different temperatures. These data are representative of a very wide range available in the published literature and are quoted here purely by way of examples. The data shown correspond to low or relatively low losses so that the real part of the permittivity $`ϵ^{}(\omega )`$ does not vary by a large amount over the frequency range in question; and it is therefore possible to normalise the data by dividing $`\chi ^{\prime \prime }(\omega )`$ by the corresponding $`ϵ^{}(\omega )`$, thus defining the loss tangent $`\mathrm{tan}\delta =\chi ^{\prime \prime }(\omega )/ϵ^{}(\omega )`$, which gives absolute values of loss for the various samples while retaining the frequency dependence of $`\chi ^{\prime \prime }(\omega )`$. The essential conclusion, coming from the analysis of the data presented in Figure 1, is that the residual loss follows in some cases law (5) while in other cases it is essentially “flat” in frequency with minor perturbations. It is interesting to note that the addition of 17% of carbon black to polyethylene causes a rise of $`\mathrm{tan}\delta `$ by one order of magnitude but hardly any change of the frequency dependence. It is also noteworthy that the total range of $`\mathrm{tan}\delta `$ covers at most three decades for materials for which the direct-current conductivity varies by many orders of magnitude. It is not profitable to fit these data to any particular law since they are the result of the interplay of minor accidental factors superimposed on the prevailing broad “flatness”. A theoretical justification for the existence of the fractional power law of type (5) with the exponent $`n`$ arbitrarily close to unity would be sufficient for the understanding of all flat losses, any minor deviations being purely accidental.

In our theoretical approach to the analysis of the universal-response characteristics (5) we make use of the concept of dipolar screening introduced by Jonscher . In the framework of this concept, the effectiveness of screening depends on the relative magnitudes of the thermally-activated-dipole density $`N_d`$ with an energy $`W`$, $`N_d\mathrm{exp}(W/kT)`$, and the “critical” density $`N_1kT/\mu ^2`$ where $`\mu `$ is the dipole moment of the dipole being screened. The theory of dipolar screening predicts that for low dipolar density, $`N_dN_1`$, the number of dipoles within the field of any one dipole increases almost exponentially, so that their behaviour is strongly collective and the etire system behaves in the universal manner (5). On the other hand, for a high dipolar density, $`N_dN_1`$, screening is effective, meaning that dipoles do not “see” their neighbours and behave individually in a Debye-like manner. This is summarized in Table 1. The situation is in many respects similar to the classical screening by charged particles but the onset of screening is much more rapid with a fall in the dipolar density $`N_d`$.

A considerable effort that has been devoted in the past to finding a theoretical explanation of the empirically observed results points on the two most widespread and at the same time least understood properties of the relaxation responses:

* the existence of the characteristic property (5) with its limiting form of virtually frequency-independent loss;
* the fact that this property of the relaxation responses is common to a very wide range of materials with very different physical and chemical interactions.

As a consequence, in theoretical attempts to model relaxation it has been commonly assumed that the empirical relaxation laws reflect a kind of general behaviour which is independent of the details of examined systems. This idea has stimulated the proposal of several relaxation mechanisms that differ mainly in the interpretation of the relaxation function. In recent attempts to find the origins of the non-Debye relaxation patterns the idea of complex systems as the “structures with variations”  that are characterised through a large diversity of elementary units and strong interactions between them is of special importance. The evolution in course of time of physical properties of a complex system is nonpredictable or anomalous , and the main feature of all the dynamical processes in such a system is their stochastic background. In the framework of statistical models (see e.g. \[4-6, 8, 9, 11-14, 17-21, 26-31\] the fact that the large scale behaviour of complex systems shows universality, i.e., that it is to some extent independent of the precise local nature of the considered system, should come as no surprise. Intuitively, one expects “averaging principles” like the law of large numbers to be in force. However, it turns out to be very hard to make this intuition precise in concrete examples of stochastic systems with a large number of locally interacting components.

The empirical facts stress the need for a completely novel approach to the modelling of the dielectric relaxation and, to a certain extent, also mechanical relaxation, photoconduction, photoluminescence and chemical reaction kinetics (sharing some common features, see ). The need to understand the connections between the macroscopic property (5) of the relaxing complex system and the statistical properties of individual molecular or dipolar species requires the introduction into the relaxation theory of advanced methods of stochastic analysis. As shown by us \[8, 32-34\], a general formalism of limit theorems of probability theory plays an important role in constructing tools to relate the local random characteristics of the complex system to the empirical, deterministic relaxation laws, regardless of the specific nature of the system considered. The significance of the present paper lies in the fact that no one has clarified in a simple and plausible way, let alone one based on a stochastic argument, why the universal relation should exist at all and why the residual loss is such slowly variable function of frequency. Our approach provides a rigorous explanation for the most widely observed form of frequency dependence of the permittivity (or susceptibility) and thereby opens up a new and very powerful way of interpreting relaxation phenomena not only in the dielectric context.

The main objective of this paper is to focus on the approach to relaxation in the framework proposed in  as consistent with Jonscher’s screening and energy criterion ideas and providing their strict mathematical formulation. We also bring to light the spatio-temporal scaling conditions hidden behind all the well-known empirical responses. For the reader’s convenience, the mathematical details neccessary for the stochastic construction of the effective representation of a relaxing complex system are followed by extended comments.

In Section 2 we introduce the basic mathematical concept underlying the stochastic transition of a complex system from its initially imposed state. We show that the relaxational properties of the entire system can be represented by means of a random effective relaxation rate which contains information on the internal stochastic structure of the investigated system. In Section 3 we point to the origins of the statistical properties of the effective relaxation rate. We discuss the role of limit theorems of probability theory as the mathematical technique which allows us to derive the explicit relaxation formulas even with rather restricted information on local random properties of the system. In Section 4 we relate the statistical conditions yielding the well-known empirical responses to the spatio-temporal scaling properties of the relaxing complex system, and we show that they underlie Jonscher’s energy-criterion hypothesis. The last section contains conclusions.

## 2. Transition and survival probabilities

(i) Let us consider a complex physical system containing identical objects undergoing irreversible transitions from state $`A`$, imposed at time $`t=0`$, to state $`B`$ at random instants of time. States $`A`$ and $`B`$ differ in some physical parameter, so that the transition $`AB`$ is defined as the change of this particular parameter (changes in all other parameters may also have an influence on the transition). Let us choose one of the objects. Consider the conditional probability $`p(t,dt)`$ that this object will undergo the transition during the time interval $`(t,t+dt)`$ if the transition has not occured before time $`t`$, i.e.

$$p(t,dt)=\mathrm{Pr}(t\theta t+dt|\theta t)$$
(7)
where $`\theta `$ is the random waiting time for the transition of the chosen object. The conditional probability $`p(t,dt)`$ defined in (7) can be expressed in a form

$$p(t,dt)=\frac{\mathrm{Pr}(\theta t+dt)\mathrm{Pr}(\theta t)}{\mathrm{Pr}(\theta t)}$$

where $`\mathrm{Pr}(\theta t+dt)`$ and $`\mathrm{Pr}(\theta t)`$ are the survival probabilities, i.e. the probabilities that the considered object will remain in state $`A`$ until time $`t+dt`$ and $`t`$, respectively. The survival probability of the object can be expressed as

$$\mathrm{Pr}(\theta t)=1p(t)$$

where

$$p(t)=\mathrm{Pr}(\theta <t)$$

is the waiting time distribution of the object, i.e. the total probability of its transition $`AB`$ until time $`t`$. One can rewrite the survival probability in the form

$$\mathrm{Pr}(\theta t)=\mathrm{exp}(\underset{0}{\overset{t}{}}r(s)𝑑s)$$
(8)
which is dependent on a non–negative quantity $`r(s)`$ called the intensity of transition . This quantity is time–dependent, in general; and moreover, because of random impacts affecting each object, for different objects in the system it can take various values at the same instant of time .

The survival probability $`\mathrm{Pr}(\theta t)`$ can be derived if one knows the explicit form of the intensity $`r(s)`$. For the time independent intensity $`r(s)=b_0=const`$ one gets

$$\mathrm{Pr}(\theta t)=\mathrm{exp}(b_0t)$$

what recovers the classical exponential evolutionary law for the object; the value $`b_0`$ determines the relaxation rate of the transition process. If, on the contrary, the intensity of transition is essentially time–dependent, then the evolutionary law for the object is of the nonexponential form. In this case the resulting relaxation rate is not directly given by the intensity of transition as it is in the exponential case. In fact, taking into account the peculiarity of local random environment of an object, one should assume that, in general, the survival probability $`\mathrm{Pr}(\theta t)`$ has the form of the weighted average of an exponential decay with respect to the probability distribution $`F_\beta (b)`$ of the relaxation rate of the object

$$\mathrm{Pr}(\theta t)=\underset{0}{\overset{\mathrm{}}{}}e^{bt}𝑑F_\beta (b).$$
(9)
In other words, the relaxation rate of the object is the random variable $`\beta `$ such that the total survival probability has the form

$$\mathrm{Pr}(\theta t)=\mathrm{exp}(\beta t)$$
(10)
where the mean value $``$ is taken with respect to the relaxation-rate probability distribution $`F_\beta (b)`$. Expressions (9) and (10) are a generalisation of the Debye relaxation for which the relaxation rate $`\beta `$ takes a constant value, say $`b_0`$, with probability 1, $`\mathrm{Pr}(\beta =b_0)=1`$. Then $`\beta `$ has a degenerate probability distribution, $`dF_\beta (b)=\delta (bb_0)db`$ (where $`\delta ()`$ denotes the Dirac delta function), and the integral on the right–hand side of (9) equals

$$\underset{0}{\overset{\mathrm{}}{}}e^{bt}\delta (bb_0)𝑑b=\mathrm{exp}(b_0t).$$

In general, in a complex system the probability distribution $`F_\beta (b)`$ and, consequently, the survival probability $`\mathrm{Pr}(\theta t)`$ of an object are of unknown forms.

From equation (8) the intensity of transition from an initial state, $`r(t)`$, is related to the total survival probability $`\mathrm{Pr}(\theta t)`$ of the object as follows

$$r(t)=\frac{d}{dt}\mathrm{ln}\mathrm{Pr}(\theta t).$$
(11)
Comparing (11) and (9) we obtain the relationship between the time–dependent intensity of transition and the relaxation rate distribution

$$r(t)=\frac{d}{dt}\mathrm{ln}\underset{0}{\overset{\mathrm{}}{}}e^{bt}𝑑F_\beta (b).$$

The above mentioned lack of information about the relaxation rate distribution $`F_\beta (b)`$ yields that the corresponding intensity $`r(t)`$ is of unknown form. The explicit form of $`r(t)`$ depends on the characteristics of the random environment around the examined object and the rules needed to specify the sets of deterministic and stochastic parameters.

(ii) Consider (as before) a system of $`N`$ objects, each waiting for transition $`AB`$ for some random time. For the $`i`$th object, $`1iN`$, let the waiting time be denoted by $`\theta _{iN}`$. The notation here (i.e. index “$`iN`$”) emphasizes the impact of the system size $`N`$ on the behaviour of each individual object that has to be taken into account. The non–negative waiting times $`\theta _{1N},\mathrm{},\theta _{NN}`$ form a sequence of independent, identically distributed random variables. The waiting time distribution $`p_{iN}(t)=\mathrm{Pr}(\theta _{iN}<t)`$ is hence the same for each $`i`$ and is equal to some function $`F_\theta (t)`$ that may depend on $`N`$.

The objects undergo transition in a certain order that can be reflected in the notion of order statistics  $`\theta _{(1)}\mathrm{}\theta _{(N)}`$, which is simply a nondecreasing rearangement of times $`\theta _{1N},\mathrm{},\theta _{NN}`$. Traditionally, $`\theta _{(l)}`$ is called the $`l`$–th order statistics of sample $`\theta _{1N},\mathrm{},\theta _{NN}`$. Note that $`\theta _{(1)}=\mathrm{min}(\theta _{1N},\mathrm{},\theta _{NN})`$ and $`\theta _{(N)}=\mathrm{max}(\theta _{1N},\mathrm{},\theta _{NN})`$.

For a fixed size $`N`$ of the system, the ratio of objects not transformed up to time $`t`$ is equal to

$$1\frac{\eta _N(t)}{N}$$
(12)
where $`\eta _N(t)`$ denotes an unknown (random!) number of objects already transformed at time $`t>0`$. For $`l=0,1,\mathrm{}N`$ the events $`\{\eta _N(t)=l\}`$ that up to time $`t`$ exactly $`l`$ transitions occured in the system can be expressed via order statistics in the following way:

$$\begin{array}{ccc}\{\eta _N(t)=0\}\hfill & =\hfill & \{\theta _{(1)}>t\},\hfill \\ \{\eta _N(t)=l\}\hfill & =\hfill & \{\theta _{(l)}t,\theta _{(l+1)}>t\}\mathrm{for}l=1,\mathrm{}N1,\hfill \\ \{\eta _N(t)=N\}\hfill & =\hfill & \{\theta _{(N)}t\},\hfill \end{array}$$
(13)
and the probability of such events is given by the Bernoulli model

$$\mathrm{Pr}(\eta _N(t)=l)=\left(\genfrac{}{}{0pt}{}{N}{l}\right)\left[p_{iN}(t)\right]^l\left[1p_{iN}(t)\right]^{Nl},l=0,1,\mathrm{}N.$$
(14)
Thus the random number $`\eta _N(t)`$ has the Bernoulli (binomial) distribution $`(N,p)`$ with parameter $`p=p_{iN}(t)`$. Consequently, the fraction in (12) is a random variable. However, it follows from the strong law of large numbers  that for $`N`$ large enough

$$\frac{\eta _N(t)}{N}\frac{\eta _N(t)}{N}=p_{iN}(t)=F_\theta (t)$$
(15)
i.e. the random nature of the fraction in (12) vanishes for large $`N`$. Hence (12) is asymptotically equal to $`(1p_{iN}(t))`$ which is the survival probability $`\mathrm{Pr}(\theta _{iN}t)`$ of any single object:

$$1\frac{\eta _N(t)}{N}1p_{iN}(t)=\mathrm{Pr}(\theta iNt).$$
(16)
The basic property of the relaxation function $`\varphi (t)`$ is its monotonic decrease from 1 at $`t=0`$ to 0 as $`t\mathrm{}`$. This is in fact a property of the survival probability and therefore (12) has been used as a definition of the relaxation function in several models. Yet (12) can be used to describe the evolution of the entire system only if the behaviour of the system is represented by any individual object (from those forming the system). Unfortunately, this strict condition seems to contradict the idea of complexity of the investigated systems, and hence models defining the survival probability of a complex system as in (16) do not capture the nature of relaxation phenomena. Moreover, the relation (15) holds for any arbitrarily chosen form of the waiting time distribution (e.g. $`F_\theta (t)=0`$ for $`t<0`$, $`=t`$ for $`0t1`$, and $`=1`$ for $`t>1`$) so that the ratio in (15) is not uniquely determined. Therefore the attempts leaving the probabilistic analysis of the irreversible stochastic transitions in complex systems at this stage can only propose the form of distribution $`F_\theta (t)`$ fitting the data most exactly, see e.g. . They do not explain the observed fractional power law (5) indicating strictly limiting properties of the survival probability of the initial state of a complex system.
(iii) The considerations of irreversible stochastic transitions in complex systems show that the behaviour of the system as a whole, in general, cannot be attributed to any chosen object forming the system . Crucially relevant to this statement is the question of a proper mathematical construction of an “averaged” imaginary object respresenting the entire system. As it follows from (8), the intensity of transition from an initial state for any (real or imaginary) object depends on its survival probability. Let us denote the survival probability of an imaginary object representing the system as a whole by $`\mathrm{Pr}(\stackrel{~}{\theta }_Nt)`$ and by $`\stackrel{~}{\theta }_N`$ the effective waiting time for the entire system. Since the considered imaginary object represents the whole system, the survival probability $`\mathrm{Pr}(\stackrel{~}{\theta }_Nt)`$ is the probability that the transition of the system as a whole from its initial state (imposed by external constraints at $`t=0`$) has not happened prior to a time instant $`t`$. This mathematical quantity is defined as the probability that there is no transition occuring in the system up to time $`t`$ so that

$$\mathrm{Pr}(\stackrel{~}{\theta }_Nt)=\mathrm{Pr}(\eta _N(t)=0).$$

By means of the Bernoulli scheme (14) with $`l=0`$, the survival probability of the system is the product of $`N`$ factors, each asymptotically equal to (12), see equation (16):

$$\mathrm{Pr}(\stackrel{~}{\theta }_Nt)=\mathrm{Pr}(\eta _N(t)=0)=(1p_{iN}(t))^N.$$

On the other hand, by means of the order statistics $`\{\eta _N(t)=0\}=\{\theta _{(1)}t\}`$, see equation (13); and the survival probability of the entire system

$$\mathrm{Pr}(\stackrel{~}{\theta }_Nt)=\mathrm{Pr}(\eta _N(t)=0)=\mathrm{Pr}(\theta _{(1)}t)=\mathrm{Pr}(\mathrm{min}(\theta _{1N},\mathrm{},\theta _{NN})t)$$

is just the probability that the first passage of the system from the initial state has not happened before time $`t`$.

Since the waiting times $`\theta _{1N},\mathrm{},\theta _{NN}`$ are independent and identically distributed, we have

$$\mathrm{Pr}(\stackrel{~}{\theta }_Nt)=\mathrm{Pr}(\theta _{1N}t)\mathrm{}\mathrm{Pr}(\theta _{NN}t).$$

Assuming (10) for each $`\theta _{iN}`$, i.e.

$$\mathrm{Pr}(\theta _{iN}t)=\mathrm{exp}(\beta _{iN}t)$$

where random variables $`\beta _{1N},\mathrm{},\beta _{NN}`$ are also independent and identically distributed, we obtain that the survival probability of the system is of the form analogous to (10):

$$\mathrm{Pr}(\stackrel{~}{\theta }_Nt)=\mathrm{exp}(t\stackrel{~}{\beta }_N)$$
(17)
where

$$\stackrel{~}{\beta }_N=\underset{i=1}{\overset{N}{}}\beta _{iN}.$$
(18)
The random variable $`\stackrel{~}{\beta }_N`$ can be considered as the effective relaxation rate representing the entire system; and, in accordance with the rate-theory concept (see e.g. ), individual relaxation rates $`\beta _{iN}`$’s are hence some contributions to the total rate $`\stackrel{~}{\beta }_N`$.

In fact, in any complex dielectric system under a weak external electric field only a part of the total number $`N`$ of dipoles, referred to as “active” dipoles, is directly governed by changes of the field. The exact number of such dipoles depends on temperature, density of the dipolar species in the system, and interactions between them, and it is usually unknown. (In some cases it may take the value of the total number of the dipoles in the system.) Therefore, there can be only a part of dipoles contributing to the effective relaxation rate, and it is reasonable to modify formula (18) taking some unknown number $`\nu _N`$ instead of exactly $`N`$ components:

$$\stackrel{~}{\beta }_N=\underset{i=1}{\overset{\nu _N}{}}\beta _{iN}.$$
(19)
## 3. Origins of the random effective relaxation rate

The relaxation function $`\varphi (t)`$, which is the probability that the initially imposed state of a macroscopic system survives by time $`t`$, is well known to be expressed as the weighted average of an exponential decay with respect to the distribution of the effective relaxation rate $`\stackrel{~}{\beta }_N`$ , and its explicit form depends on statistical properties of the rate. In view of (17), we can use the survival probability $`\mathrm{Pr}(\stackrel{~}{\theta }_Nt)`$ as a definition of the relaxation function and we obtain in this case that

$$\varphi (t)=\mathrm{exp}(\stackrel{~}{\beta }_Nt)$$
(20)
where $`\stackrel{~}{\beta }_N`$ is given by (19).

The evaluation of $`\nu _N`$, the number of components in (19), which is in agreement with the nature of the relaxation phenomenon, has been presented in . It is based on the following idea:

Depending on the strength of screening, the active dipoles in the complex system may “see” each other to some extent. The collective behaviour of active dipoles results in the appearance of mesoscopic cooperative regions contributing to the macroscopic relaxation process . The number $`L_N`$ of such mesoscopic regions in the system is determined by their sizes $`M_1,M_2,\mathrm{}`$ and by the number of the active dipoles in the system .

The effective relaxation rate $`\stackrel{~}{\beta }_N`$ consists of the contributions $`\overline{\beta _{jN}}`$ of all $`L_N`$ mesoscopic cooperative regions

$$\stackrel{~}{\beta }_N=\underset{j=1}{\overset{L_N}{}}\overline{\beta _{jN}}.$$
(21)
Similarily, the relaxation rate $`\overline{\beta _{jN}}`$ of the $`j`$th mesoscopic region is the sum of the contributions $`\beta _{iN}`$ of all active dipoles over the region \[13, 32-34\] and hence it is equal to

$$\overline{\beta _{jN}}=\underset{i=M_1+\mathrm{}+M_{j1}+1}{\overset{M_1+\mathrm{}+M_j}{}}\beta _{iN}.$$
(22)
As mentioned in section 2 (i), the individual relaxation rate, denoted here by $`\beta _{iN}`$, reflects the random intra-cluster dynamics, i.e. statistical properties of the interactions of the $`i`$th active dipole with some inactive neighbours forming around it a cluster of size $`N_i`$.

Comparing (19) and (21) with $`\overline{\beta _{jN}}`$ of the form (22) one obtains the explicit formula for $`\nu _N`$, the number of the contributions to the total relaxation rate, in the framework of the proposed approach. For the sake of this paper it is not neccessary to cite the entire complicated formula derived (for details see ). It is only important to note that $`\nu _N`$ is fully determined by $`N`$, the system size; by the cluster sizes $`N_1,N_2,\mathrm{}`$; and by the cooperative-region sizes $`M_1,M_2,\mathrm{}`$. The number of the dipoles directly engaged in the relaxation process, their location, interaction range, relaxation rates, and all the quantities defined by them are random. Their precise stochastic characteristics determining the relaxation response of the system are unknown, unfortunately. However, as a rule, the relaxing systems consist of a large number of dipoles so that the distribution of $`\stackrel{~}{\beta }_N`$ can be satisfactorily approximated by the weak limit

$$\stackrel{~}{\beta }=\underset{N\mathrm{}}{lim}\stackrel{~}{\beta }_N.$$
(23)
In practice, even $`N10^6`$ can suffice to replace adequately $`\stackrel{~}{\beta }_N`$ in (20) by the limit $`\stackrel{~}{\beta }`$.

On the basis of the limit theorems of probability theory the distribution of the effective relaxation rate $`\stackrel{~}{\beta }`$ in (23) and the resulting form of the relaxation response (20) in time or the corresponding susceptibility in frequency domain

$$\chi ^{}(\omega )i\chi ^{\prime \prime }(\omega )\varphi ^{}(\omega )=\underset{0}{\overset{\mathrm{}}{}}e^{i\omega t}\left(\frac{d\varphi }{dt}(t)\right)𝑑t$$

can be derived even with rather restricted information on properties of the micro/mesoscopic statistical levels in the system. It appears  that each of the well-known dielectric responses can be obtained in this way if appropriate conditions are imposed on the basic random quantities: the cluster sizes $`N_i`$, the sizes $`M_j`$ of the cooperative regions of correlated clusters, and the individual active-dipole relaxation rates $`\beta _{iN}`$.

## 4. Foundations of the universal response

Table 2 collects the statistical properties shown  to underlie the KWW, HN, CC, CD, and Debye relaxation responses, given by formulas (3) and (4). As it follows from the table, the empirical responses are realized by various stochastic schemes. The conditions imposed on $`N_i`$, $`M_j`$, and $`\beta _{iN}`$ here take one of two forms; either the distribution of the quantity considered is of finite mean value or it has a heavy tail. Let us explain that the distribution of a nonnegative random variable, say $`X`$, has a heavy tail if for some $`0<a<1`$ the tail function $`\mathrm{Pr}(X>x)x^a`$ for large $`x`$, i.e.,

$$\underset{x\mathrm{}}{lim}\frac{\mathrm{Pr}(X>x)}{x^a}=\mathrm{const}>0.$$
(24)
Classical example of the heavy-tailed distribution is the Pareto distribution

$$\mathrm{Pr}(X>x)=\frac{1}{1+(Ax)^c},x>0$$

with the shape parameter $`0<c<1`$ and any positive scale constant $`A`$. (In this case the heavy-tail exponent $`a=c`$.) Another distribution possesing this property is a generalisation of the Pareto example; namely, the Burr law

$$\mathrm{Pr}(X>x)=\frac{1}{(1+(Ax)^c)^d},x>0$$

with the shape parameters $`c,d>0`$ such that $`0<cd<1`$ and any positive scale constant $`A`$. (Here the heavy-tail exponent $`a=cd`$.) More difficult but interesting heavy-tailed distributions are the completely asymmetric Lévy-stable laws  defined by their Laplace transforms $`(X;t)=\mathrm{exp}((At)^c)`$ with $`0<c<1`$ and $`A>0`$. (In this case the heavy-tail exponent $`a=c`$.)

Condition (24) applied to any random variable $`X`$ expresses the following scaling property of the magnitude represented by $`X`$:

$$\mathrm{Pr}(XCx)C^a\mathrm{Pr}(Xx),x\mathrm{},$$
(25)
for any fixed constant $`C>0`$. (If there exists the corresponding probability density $`w(x)`$, the condition (25) yields $`w(Cx)C^{a1}w(x)`$ for large $`x`$.) Hence, in the presented approach the heavy-tail property is directly related to the spatial (if referred to $`N_i`$ and $`M_j`$, the cluster and the cooperative-region sizes) or temporal (if referred to the relaxation rate $`\beta _{iN}`$) scaling properties of the system. Comparing the non-Debye responses with the Debye one, as done in Table 3, one can see that the difference between them lies in the presence of such a scaling in the system. Namely, in the Debye case, that of the individual behaviour of the active dipoles, there is no property of this form. By contrast, spatial or/and temporal scaling can be found for other responses. Moreover, at least one of the two variates $`M_j`$ and $`\beta _{iN}`$ has a heavy-tailed distribution in the non-Debye cases. Since, as equation (22) shows, the correlated-cluster relaxation rate $`\overline{\beta _{jN}}`$ includes both those quantities; its distribution has also a heavy tail, “producing” a hierarchy of mesoscopic relaxation rates; i.e. the same proportion of smaller or larger mesoscopic contributions to the effective response no matter the scale at which one is looking at the relaxation rate distribution. This property follows from both the stochastic properties of the active dipoles generated by the intra-cluster dynamics and the strength of screening yielding the range of interactions between them.

As a consequence of the proposed approach, the origins of the macroscopic energy criterion coefficient can be pointed out. As we have already explained, the heavy-tailed distribution of any (micro/meso/macroscopic) relaxation rate $`\beta `$ reflects the scaling property of the rate

$$\mathrm{Pr}(\beta Cb)C^a\mathrm{Pr}(\beta b),b\mathrm{},$$
(26)
for any fixed constant $`C>0`$ (where $`a`$ is the heavy-tail exponent). As in (10) and (17), we have $`\mathrm{Pr}(\theta t)=\mathrm{exp}(\beta t)`$ so that the properties of the micro/meso/macroscopic relaxation rates determine the survival probabilities of the microscopic (active dipole), mesoscopic (cooperative region), and macroscopic (entire system) objects, respectively; thus, by analogy to (20), defining the function $`\varphi (t)=\mathrm{exp}(\beta t)`$, for simplicity referred to as the micro/meso/macroscopic relaxation function.

The asymptotic behavior (26) of the relaxation-rate distribution at large $`b`$ is connected with the short-time behavior of the associated relaxation function $`\varphi (t)`$ . Namely, at the origin $`t0`$ the response function $`f(t)={\displaystyle \frac{d\varphi (t)}{dt}}`$ takes the form

$$f(t)t^{a1}L(t)$$
(27)
where $`L(t)`$ is a function slowly varying at $`t=0`$ (i.e., $`L(Ct)/L(t)1`$ as $`t0`$ for any fixed positive constant $`C`$). Then, it can be shown  that the short–time property (27) of the response function $`f(t)`$ corresponds to the following asymptotic behavior of $`\chi (\omega )\varphi ^{}(\omega )=\underset{0}{\overset{\mathrm{}}{}}e^{i\omega t}f(t)𝑑t`$ for the high-frequency region:

$$\chi (\omega )=\chi ^{}(\omega )i\chi ^{\prime \prime }(\omega )(i\omega )^aL(1/\omega )\mathrm{for}\omega \mathrm{}.$$
(28)
Property (28) leads straightforwardly to the constant ratio

$$\frac{\chi ^{\prime \prime }(\omega )}{\chi ^{}(\omega )}=\mathrm{cot}\left(\eta \frac{\pi }{2}\right)\mathrm{for}\omega \omega _p$$

with $`\eta =1a`$. This result, interpreted as the (micro/meso/macroscopic) energy criterion, is consistent with the energy criterion (6) if only $`\eta `$ is equal to the macroscopic high-frequency power-law exponent $`n`$. Summing up, the heavy-tail property of the relaxation rate with the heavy-tail exponent $`a`$ leads to the energy criterion with the characteristic constant $`\eta =1a`$.

Table 4 contains the detailed discussion of the energy criterion at the micro/meso/ macroscopic levels for the cases related to the KWW, HN, CC and CD empirical relaxation responses (3) and (4). Observe that in each particular response the parameter $`\eta `$ followed from the mesoscopic-relaxation-rate tail is equal to the macroscopic energy-criterion coefficient $`n`$, and hence, that the high-frequency universal laws (5) and (6) and its parameter $`n`$ are determined by the heavy-tail properties of the mesoscopic (correlated-cluster) relaxation rate distributions. The above proves Jonscher’s hypothesis introduced in his screening theory of relaxation and stating that in order to observe property (6) at the macroscopic level, the same property must characterise the behaviour of the relaxing dipoles. As shown in Table 4, in general, the energy-criterion parameter agrees with $`n`$ from equation (6) on the mesoscopic level, only.

## 5. Conclusions

To our knowledge, this is the first attempt to relate rigorously the widely observed universal macroscopic laws to the micro/mesoscopic statistical characteristics of the complex system, thereby establishing a basis for the understanding of the stochastic origins of the relaxation phenomenon. The essential element which was required for the derivation of these results is the strict probabilistic formalism in terms of which the dipolar screening and energy criterion concepts can be expressed.

In terms of our analysis, the conditions for the universal relaxation response may be stated as follows:

1. at the microscopic level, a random number of the active dipoles, those that follow changes of the external field, is selected; their individual relaxation rates are determined by the interactions of the active dipoles with inactive neighbours forming random-sized clusters around them;
2. at the mesoscopic level, correlated-cluster regions of sizes depending on the strength of screening appear; the collective rate of relaxation of the active dipoles in such a mesoscopic cooperative region becomes correlated with the number of the active dipoles in the region and with the stochastic properties of their individual relaxation rates, see equation (22); the mesoscopic correlated-cluster relaxation rate combines both the spatial and the temporal scaling properties of the system;
3. at the macroscopic level, averaging over the number of effective contributions and their rates leads to the universal relaxation for the entire system, giving law (5); the sufficient condition of this is that the high-frequency energy criterion is satisfied at the mesoscopic level; the mesoscopic energy criterion is provided by the heavy-tail property of the correlated-cluster relaxation rates.

Table 1.
Dipolar screening schemes.

Table 2.
Stochastic origins of the empirical relaxation responses.

Table 3.
Comparison of the non-Debye relaxation responses to the Debye one.

Table 4.
Discussion of the energy criterion at the micro/meso/macroscopic levels.

Figure captions

Figure 1.
Typical examples of the frequency dependence of the residual loss, after substraction of dc component, low-frequency dispersion and strong dipolar peaks, plotted logarithmically against frequency, for a range of dielectric materials. The data cover more than ten decades of frequency, and the prevailing trend is a “flat” or frequency-independent loss with only minor deviations, for instance the slope $`0.2`$ over six decades for pure polyethylene. The $`\mathrm{tan}\delta `$ values cover a range of barely three decades between the lowest and highest losses; and the surprising feature is the very slight effect of the addition of conducting species (like 12% graphite to polyethylene raising the loss by a factor of 3, and 17% of graphite by a further factor of 3) with very little change of the frequency dependence. Likewise, the addition of 55% of graphite to boron nitride does not raise the extremely low loss of the system. Most of the data quoted here are taken from much wider ranges for variable temperature and for ranges of composition. The evident trends are the general “flatness” of the losses in frequency which is not compatible with any known mechanism and that irrespective of the absolute level of loss. Detailed information: $`\mathrm{Pb}_2\mathrm{Sc}_{0.9}\mathrm{Ta}_{0.9}\mathrm{Ti}_{0.1}\mathrm{O}_6`$ and $`\mathrm{Pb}_2\mathrm{ScTaO}_6`$ are ferroelectric ceramics, from A. Isnin and A.K. Jonscher, Ferroelectrics, 210 (1998) 47. Pure polyethylene (PE) showing nine decades of frequency with a very slowly varying $`\mathrm{tan}\delta `$ in the range $`3\times 10^4`$ \- $`2\times 10^3`$ with a frequency dependence proportional to $`f^{0.2}`$ over most of the middle range; also shown are data for polyethylene with admixtures of 12% and 17% of graphite; from D.S. McLachlan and M.B. Heaney, Phys. Rev. B 60 (1999) 12746. Single crystal $`\mathrm{CaTiO}_3`$ over a range of temperatures with $`\mathrm{tan}\delta =59\times 10^4`$, show extremely small values of $`1n=0.000325`$ to $`0.000625`$; from B.S. Lim, A.V. Vaysleyb, and A.S.Nowick, Applied Physics A 56 (1993) 8. Alpha-poly-vinylidene di-fluoride (PVDF$`\alpha `$) data from A.K. Jonscher and G. Menegotto, IEEE Trans. Dielectrics and EI 7 (2000) 303. Percolation system graphite-boron nitride (45%BN+55%C) from J.J. Wu and D.S. McLachlan, Phys. Rev. B 56 (1997) 1236; Phys. Rev. B 58 (1998) 14880. Zeolite plus NLO-active para-aniline from M. Wuebbenhorst, private communication.

Prof. Andrew K. Jonscher: Born 1922 Warsaw. 1944 Joined the Polish Army in Italy. Engineering studies in London 1945; B.Sc. (Eng) 1st class Hons 1948; Ph.D. in semiconductors from Queen Mary College, University of London, 1951. 1952-1963 on the Staff of the Hirst Research Center of the General Electric Company in Wembley, England, finishing as Leading Scientific Staff. 1963-1965 Reader in Solid State Physics at Chelsea College, University of London; 1965-1985 Professor of Solid State Electronics there. 1970 established Chelsea Dielectrics Group, one of the leading groups working on dielectric relaxation in solids, formulated the “universal” dielectric response, known in the literature as Jonscher’s law. 1983 set up Chelsea Dielectrics Press to publish his first monograph “Dielectric Relaxation in Solids” which was too “hot” for commercial publishers to touch. 1987 reached retirement age and transferred as Visiting Professor to Royal Holloway, University of London, where he continued experimental and theoretical work on relaxation. 1990 nominated Whitehead Memorial Lecturer of IEEE Electrical Insulation Society. 1996 Chelsea Dielectrics Press publishes his second monograph “Universal Relaxation Law”. 1998 Honorary Doctorate of the Technical University in Łódź, Poland, for “pioneering work on dielectric relaxation”. Most recently theoretical studies of relaxation jointly with the Wrocław group.
Dr Agnieszka Jurlewicz: Ph.D. in Mathematics (1994), Institute of Mathematics, Faculty of Basic Problems of Technology, Wrocław University of Technology. Since October 1995 Assistant Professor at the Institute of Mathematics, Wrocław University of Technology. Research topics: Theoretical investigations of a new type of continuous-time random walk; application of the obtained results in the probabilistic model of relaxation phenomena in disordered systems; studies on the relationship between the model and the Adam-Gibbs theory.
Prof. Karina Weron: Ph.D. in Theoretical Condensed Matter Physics (1978), Faculty of Basic Problems of Technology, Wrocław University of Technology. Since May 2001 Full Professor at the Institute of Physics, Wrocław University of Technology. Research topics: Statistical and stochastic methods in physics (random matrices, random walks, self-similar stochastic processes, time series analysis); modelling of relaxation phenomena; dynamical properties of complex systems; ionic transport through biological membranes.


text:
# Discovery of A New Class of Pulsating Stars: Gravity-mode Pulsators among Subdwarf B Stars

## 1 Introduction

Subdwarf B (sdB) stars are evolved hot stars (22,000 to 40,000 K) commonly found in the disk of our Galaxy (e.g. Saffer et al. 1994). They are widely believed to be core He-burning stars of $`0.5M_{\mathrm{}}`$ on the extended horizontal branch (EHB) (Heber 1986; Saffer et al. 1994), having lost nearly all of their envelopes near the first red giant tip. The existence of multimode pulsators among sdB stars presents a powerful opportunity for further investigation of their structure using asteroseismological techniques. In addition to clarifying the uncertain evolutionary histories of the sdB stars themselves, the results will be relevant to the later evolution of the vast majority of normal stars. Buried under their huge envelopes, all low mass He-burning red giants have cores that are basically identical to sdB stars. It is our good fortune that the latter have shed all but a few wisps of their obscuring envelopes, and that some appear to pulsate with exactly the right complexity required for asteroseismology.

The first pulsating sdB stars were discovered only recently (Kilkenny et al. 1997). Commonly called EC 14026 stars, after the prototype111EC 14026 stars are now officially V361 Hya stars., they are multimode pulsators with typical periods of 100–250 s, and a total range of about 80–600 s. Independently and nearly simultaneously, their existence was predicted by Charpinet et al. (1996, 1997), due to the presence of low-order and low-degree radial and nonradial pressure modes (p-modes) excited by an iron opacity $`\kappa `$ mechanism in the thin diffusion-dominated envelopes. Additional EC 14026 stars were soon found by teams at the South African Astronomical Observatory (Kilkenny 2002, and references therein), the University of Montreal (Billères et al. 2002), and the Nordic Optical Telescope (NOT) (Silvotti et al. 2002, and references therein), bringing the current total to about 30. The justifiable excitement over these stars (Koen et al. 1998a; Fontaine et al. 1998; O’Donoghue et al. 1999; Charpinet, Fontaine, & Brassard 2001; Kilkenny 2002) stems from the fact that stellar and atmospheric models for sdB stars appear to be sufficiently realistic for nonadiabatic calculations to be extremely successful in explaining their properties (Charpinet et al. 2001). In agreement with the theoretical results, EC 14026 variables are found primarily among the hotter sdB stars, clustering around $`T_{\mathrm{eff}}`$ $``$ 33,500 K and log g $``$ 5.8, with only a few low gravity examples at somewhat cooler temperatures.222PG 1605+072 (Koen et al. 1998b); Feige 48 (Koen et al. 1998c); HS 2201+2610 (Østensen et al. 2001a); HS 0702+6043 (Dreizler et al. 2001) The pulsation amplitudes are generally less than a few hundredths of a magnitude.

In this letter, we introduce another chapter in sdB asteroseismology: the discovery of a new class of multimode pulsators. The newly identified stars have apparent periods at least a factor of ten longer than those observed in EC 14026 stars, implying that they must be due to gravity modes (g-modes). The longer periods present both advantages and disadvantages for further observations. The most important gain may be the greatly increased range and number of suitable EHB candidates for asteroseismology. The following sections describe the observations and apparent trends identified in the data, discuss their implications for the unknown excitation mechanism, and summarize our conclusions. Subsequent papers will examine individual pulsators in more detail.

## 2 Observational Characteristics

The prototype long period sdB variable, PG 1716+426, is an unremarkable faint blue star less than 1$`\mathrm{°}`$ south of M92 in Hercules ($`V`$ 13.97, Wesemael et al. 1992). Its variability was first detected in 1999, though its nature wasn’t fully appreciated for nearly two years. The discovery light curve (shown at the top of Figure 1) was taken amidst thin cirrus at the Steward 2.3 m telescope on Kitt Peak, during initial tests of a pilot program to monitor sdB binaries with large radial velocity variations. We were searching for eclipses by white dwarfs, ellipsoidal effects, and reflection effects, all of which have amplitudes less than a few hundredths of a magnitude and timescales longer than 250 s. Except where noted, all data for this paper were obtained at the Arizona 2.3 m and 1.6 m telescopes on Kitt Peak and Mt. Bigelow using conventional 2K CCD’s binned 3x3 (0.45″/pixel) to reduce the overhead time to 24 s. Depending on the target, integration times of 15–60 s were chosen to maximize the counts and minimize photometric errors. Differential magnitudes were derived from aperture photometry of the sdB relative to reference stars of comparable magnitude within each frame.

Six additional light curves for PG 1716+426, 2–4 hours each during April–June 2000, showed persistently irregular variations in both amplitude and phase, that were not present in any of several nearby reference stars. Both low resolution (Saffer et al. 1994) and very high S/N 1Å resolution MMT spectra (Green et al., in preparation), and multicolor photometry (Wesemael et al. 1992) proved that PG 1716+426 is a normal, moderately cool sdB star, but the timescales seemed much too long for an EC 14026 star. Eventually we detected similar variations in a second sdB, PG 0101+039, although at lower amplitude. We were allocated 18 nights at Mt. Bigelow to monitor PG 0101+039, but were able to observe for only one night due to unseasonable bad weather. We finally acquired 81 hours of CCD photometry for PG 1716+426 during 14 nights in April–June 2001: 12 clear nights at the Mt. Bigelow 1.6 m telescope, several adjoining hours at the end of a partially cloudy night on the Calar Alto 2.2 m, a separate night on the NOT, and another with the Fick Observatory 0.6 m telescope. At Mt. Bigelow and at Fick, all data were taken with an $`R`$ filter. The NOT observations cycled rapidly between $`V`$, $`R`$, and $`I`$ with ALFOSC (Østensen et al. 2001b), while the BUSCA camera at Calar Alto allowed simultaneous photometry in four wide passbands from $`U`$ to $`I`$ (Falter et al. 2003). A first analysis of all the $`R`$ data found at least 3–5 pulsation modes with periods between 0.8 and 1.4 hours. Further details will be presented by Reed et al. (in preparation).

Once pulsations were confirmed in PG 1716+426, we broadened the 1.6 m monitoring program to include a nearly random sample of bright sdB stars ($`V\mathrm{\hspace{0.17em}14}`$). We have now identified 20 stars (including PG 1716+426) that appear to be multimode pulsators with periods of the order of an hour, from a total sample of 72. Each star was monitored for 2 to 7 hours at a time, and most were observed on multiple nights to allow for the possibility of occasional cancellation of modes. (PG 1716+426 exhibited a nearly flat light curve for nearly four hours once, in 14 nights of observation.) The peak-to-peak amplitudes in the light curves are comparable with those of most EC 14026 stars. The strongest new pulsators clearly show the characteristic beating effects of multiple modes with close periods (Fig. 1). However, about half of our variables pulsate so weakly that relative light variations are detectable only because of the stability, precision, and good sky coverage of modern CCD’s. We classify the weakest pulsators as bonafide members of the class only when light variations are clearly detected on more than one night, but even in their “low” states, the light curves are usually not as flat as those of non-pulsators. Daily aliasing is obviously a major problem, since fewer than 8 to 10 cycles per night can be observed from a single site. Extensive, well-coordinated multi-site campaigns or space observations over fairly long time baselines are required for good mode characterization.

The only selection biases arise from a preference for brighter targets, mostly from the PG catalog (Green, Schmidt, & Liebert 1986), and our avoidance of sdB stars with composite spectra (for better temperature and gravity determinations). Otherwise, we selected stars having the best sky positions at the time of observation. We therefore expect our sample to be representative of all sdB stars for properties such as temperature, gravity, presence of close companions, etc. For example, PG 1716+426333We derive an orbital period of 1.77739d for PG 1716+426; Morales-Rueda et al. (2002) found 1.77732d. and the next few members of the class were all found in single-lined spectroscopic binaries with orbital periods of less than a few days, but this was due only to our pre-discovery sample selection. Now that we have obtained light curves for most of a representative subsample of 52 sdB stars for which we have precise radial velocities from multiple spectra, we find no binary correlation. Specifically, 29 of the 52 sdB stars (56%) are short period binaries, having the same apparent distribution in temperature and gravity as the rest of the sample. 13 of the 52 are also long period pulsators. If the pulsations are independent of close binary status, we expect 7.25 $`\pm `$ 2.4 binary pulsators. The observed number is 7, indicating that the fraction of short period binaries among PG 1716-like stars is no different from that of sdB stars in general.

The present data suggest an upper limit of 90–120 minutes for the major periods of all the new pulsators. We have not found any sdB pulsators with obvious periods as long as the 7.8 hour sawtooth modulation detected by Orosz & Wade (1999) in KPD 0422+5421, although we did find three low amplitude variables with periods of several hours or more. Two are sdB stars exhibiting apparent reflection effects, having appropriately phased sinusoidal light curves the same length as their orbital periods ($``$ 6h and 12h). The third was a reference star with an early F spectrum and a 6h period, a likely $`\delta `$ Scuti star.

A general characteristic of nonradial pulsators is the expected amplitude dependence of the pulsations as a function of wavelength, although the largest effects are in the ultraviolet (e.g. Robinson et al. 1995; Brassard, Fontaine, & Wesemael 1995). Multicolor amplitudes have been reported for two short period (p-mode) pulsators. Koen (1998) observed somewhat larger $`U`$ amplitudes in KPD 2109+4401, but no significant trends in $`BVR`$. Falter et al. (2003) found similar behavior in PG 1605+072. Following Brassard et al. (1995), we computed the wavelength dependence of the amplitude of a g-mode in a 30,000 K, log g = 5.5 sdB model. We find that the theoretical curve is rather flat in the red and increases somewhat toward the blue: predicted amplitude ratios are 1.09 for $`B/V`$, 1.38 for $`U/V`$. Our observations appear mostly consistent. The NOT data for PG 1716+426 show no obvious systematic trends in $`VRI`$, and (non-simultaneous) 2.3 m $`BVR`$ amplitudes are fairly comparable. At first glance, the largest amplitude in the short and somewhat noisier BUSCA data appears to be in the $`I`$ filter, but closer inspection reveals that the amplitudes cannot be disentangled from surprisingly large slopes in the baselines (mainly in the wide $`U`$ and $`I`$ filters). Subsequent experiments with simultaneous $`B`$ and $`R`$ photometry in Arizona demonstrated that the baseline trends are almost certainly due to differential extinction between the sdB and the much redder comparison stars. With hour long timescales, there are serious problems at bluer wavelengths or with wider bandpasses, because the effect of differential extinction is comparable to the pulsational amplitudes. At Mt. Bigelow and Kitt Peak, we generally see significant differential extinction in $`B`$, a less noticeable effect in $`V`$, and no detectable effect in the (preferred) $`R`$ filter ($`I`$ is usually avoided due to fringing).

Figure 2 compares the positions of the long period and short period sdB pulsators in the log g vs $`T_{\mathrm{eff}}`$ diagram. The 13 newly identified pulsators (red and green error bars) for which we derived homogeneous temperatures and gravities are all cooler than 30,000 K, in marked contrast with the hotter EC 14026 stars (blue filled circles). We lack spectra for seven (unplotted) new pulsators, but six of those with temperatures from other sources also appear to be relatively cool (Heber 1986; Wesemael et al. 1992; Saffer et al. 1994; Morales-Rueda et al. 2002). At first glance, there seems to be little overlap in temperature and gravity between the long and short period pulsators, and no long period variations have been found in any EC 14026 stars444This is a weak constraint, since only a few EC 14026 stars are bright enough for our survey, and several of these have such complicated light curves that it would be very difficult to detect additional long period components, e.g. PG 1605+072 (extreme amplitudes), PG 1336-018 (eclipsing sdB + dM), or KPD 1930+2752 (sdB \+ wd with ellipsoidal variations).. However, we note that the inhomogeneous atmospheric parameters for the latter (Charpinet 2001) are not necessarily on the same scale as our data points. Our temperatures are systematically somewhat hotter than those derived by some other investigators for $`T_{\mathrm{eff}}`$ $`<`$ 32,000 K, and our gravities consequently slightly higher. (We are continuing to pursue the longstanding disagreement between various atmospheric codes for sdB stars (Wesemael et al. 1997), but it is beyond the scope of this paper.) Thus, it is probably premature to speculate about the (small) degree of separation or overlap between the two types of pulsators in Figure 2. The important point is that long period sdB pulsators clearly populate a separate temperature and gravity regime than do short period pulsators, in addition to their very different timescales.

The fraction of sdB stars exhibiting long period light variations is surprisingly high, in further contrast with the EC 14026 stars: roughly 75% of sdB stars cooler than 30,000 K, or 25–30% of all sdB stars.

## 3 Discussion

Observed periods of the order of an hour require the newly discovered sdB variables to be gravity mode pulsators. Preliminary calculations show that high radial order g-modes are needed to produce the range of periods observed in PG 1716+426. The $`\kappa `$ mechanism that has been so successful for understanding p-mode pulsations in EC 14026 stars (Charpinet et al. 1997; Brassard et al. 2001) does not appear to be capable of driving such deep pulsational modes. The g-mode driving mechanism is therefore unknown.

The limited temperature range of the long period variables in Figure 2 implies the pulsation mechanism is intrinsic to the star, rather than an external driver such as tidal excitation (Fontaine et al. 2003). The lower temperatures result from increased hydrogen envelope masses, up to 0.003–0.004 $`M_{\mathrm{}}`$ in the coolest sdB stars. However, it is not clear how such apparently negligible envelopes could cause longer period pulsations.

Canonical EHB evolutionary sequences (e.g. Dorman’s tracks, from Charpinet et al. (2000), shown in Figure 2) indicate that a star initially expands away from the zero age EHB, and then begins contracting toward higher temperatures as helium is consumed in the core. Following central helium exhaustion (where the slope of the tracks changes sign for the second time, very abruptly), helium shell burning produces a feeble attempt at further expansion, which is quickly succeeded by rapid contraction to the white dwarf stage. Our data suggest that larger amplitude g-mode pulsators (in red) occur preferentially (only?) toward the end of core helium burning (after correcting for a small apparent offset by shifting the theoretical tracks down until their starting points form a lower envelope to all the cooler stars). If strong pulsators do not exist near the zero age EHB, it is tempting to think that they might be driven by a slowly reviving hydrogen burning shell. However, in standard EHB models, the miniscule envelope does not allow the hydrogen shell to “turn on” until after the star has contracted to a very hot white dwarf. Thus, if the hydrogen envelope is really of key importance in the driving mechanism, its role could be elucidated only with nonstandard models of sdB stars. It would be interesting to know if similar g-mode pulsations can occur in sdB stars whose spectra show little if any hydrogen: the He-sdB, or sdB:He4 stars (Jeffery et al. 1997, and references therein), if any exist with comparable temperatures and gravities.

P. Moskalik (private communication) pointed out that the distribution of short- and long-period sdB pulsators in the log g vs $`T_{\mathrm{eff}}`$ diagram resembles the situation in main sequence B stars. In both cases, the hotter objects are low-order p-mode pulsators, while the cooler stars are high-order g-mode pulsators. Dziembowski, Moskalik, & Pamyatnykh (1993) found that the same mechanism could drive both types of pulsation on the main sequence. Perhaps a variant of Charpinet et al.’s (1997) $`\kappa `$ mechanism could be found to work in the cooler sdB stars, though this remains to be seen.

## 4 Summary and Conclusions

We have so far identified 20 subdwarf B stars that appear to be multimode pulsators on timescales of about an hour, more than a factor of ten longer than typical periods of the recently discovered EC 14026 pulsators. Observed peak-to-peak amplitudes are $``$ 0.05 magnitudes. Due to the longer timescales, only the prototype, PG 1716+426, has been sufficiently well observed to detect individual modes. Our first attempt, with 81 hours monitoring spread over 2 months, indicates the presence of at least 3–5 probable modes with periods between 0.8 and 1.4 hours. All long period sdB pulsators with known, spectroscopically determined temperatures are relatively cool (Teff $`<`$ 30,000 K), in contrast with the hotter EC 14026 stars. PG 1716 stars are also about six times more common.

Many of the newly identified variables are easily observed with small telescopes capable of precise differential photometry. Fast frame transfer CCD’s are not required; conventional CCD’s may have an advantage if their larger field of view permits a better selection of reference stars. Good mode characterization will require photometry over a much larger fraction of a day than can be obtained from a single site, in addition to long time baselines. With ground-based campaigns, differential extinction is a much worse problem for long period sdB pulsators than for the short period ones, but it can be minimized by using an R filter. PG 1716 stars are excellent candidates for observation from space, to eliminate both aliasing and atmospheric problems.

The long periods found in PG 1716+426, and strongly indicated in the light curves of nineteen other cool sdB variables, imply the detection of high radial order g-modes. The driving mechanism for g-modes in sdB stars is currently unknown. The distribution of the stars in the log g – $`T_{\mathrm{eff}}`$ plane suggests that their longer pulsations might be related to thicker hydrogen envelopes, but this is not easily understood from standard EHB models. Alternatively, some variant of the EC 14026 iron opacity mechanism might also be at work in cooler sdB stars, analogous to the situation in main sequence B stars (Dziembowski et al. 1993).

The relatively long pulsation timescales, the distinctly different temperature range, and the inability of the existing iron opacity mechanism to drive g-modes, all support the identification of the long period sdB stars as a separate class of pulsating variables, one with great potential to expand our understanding of evolved stars.

This investigation was made possible by NSF grants AST-9731655 and AST-0098699. We thank the Steward TAC for their generous time allocations, and P. Strittmatter for his continuing support of undergraduate research. We are grateful to the SO staff, especially E. Christensen, B. Wood, and G. Rosenbaum, for their patience with our lengthy calibrations, and all their other assistance. EMG thanks Bill Peters for many tangible and intangible contributions too numerous to mention.


text:
# Collapsed 2-Dimensional Polymers on a Cylinder

## Abstract

Single partially confined collapsed polymers are studied in two dimensions. They are described by self-avoiding random walks with nearest-neighbour attractions below the $`\mathrm{\Theta }`$-point, on the surface of an infinitely long cylinder. For the simulations we employ the pruned-enriched-Rosenbluth method (PERM). The same model had previously been studied for free polymers (infinite lattice, no boundaries) and for polymers on finite lattices with periodic boundary conditions. We verify the previous estimates of bulk densities, bulk free energies, and surface tensions. We find that the free energy of a polymer with fixed length $`N`$ has, for $`N\mathrm{}`$, a minimum at a finite cylinder radius $`R^{}`$ which diverges as $`TT_\theta `$. Furthermore, the surface tension vanishes roughly as $`(T_\theta T)^\alpha `$ for $`TT_\theta `$ with $`\alpha 1.7`$. The density in the interior of a globule scales as $`(T_\theta T)^\beta `$ with $`\beta 0.32`$.

Although the behaviour of polymers has been studied for many years, there are still a number of open questions. At high temperatures or in good solvents repulsive interactions (the excluded volume effect) dominate the conformation and the polymer chain tends to swell to a random coil. At low temperatures or in poor solvents attractive interactions between monomers dominate the conformation and the polymer chain tends to collapse and forms a compact dense globule. The coil-globule transition point is called the $`\mathrm{\Theta }`$ point.

The behaviour in good solvents is accessible by field theory and thus well understood deGen , e.g. the partition function scales as $`Z\mu ^NN^{\gamma 1}`$ at $`T>T_\theta `$. In poor solvents, much less is known exactly. For a long time it was e.g. widely believed that the partition sum scales in the same way, although already Lifshitz mean-field theory Grosberg ; Lifshitz  had predicted for $`T<T_\theta `$ that there should be also a surface term,

$$Z(N,T)a^Nb^{N^s}N^{\gamma 1},$$
(1)
with $`s=(d1)/d`$ and $`b>1`$. First numerical indications for this were given by Owczarek ; hegger ; nidras , and theoretical arguments in favour of Eq. (1) were given in duplantier . More recently, Eq. (1) was checked numerically in g02 . But it seems that there exists still no proof for Eq. (1), and it is not known how the surface tension (related to $`b`$) scales for $`TT_\theta `$. Moreover, the exponent $`\gamma `$ is not known Owczarek ; duplantier ; nidras , and it is not even clear whether it is universal in the collapsed phase.

In g02  we had studied two different geometries – globules on infinite lattices and polymers filling finite lattices with periodic boundary conditions – where surface tension should make very different contributions. In the present work we study a third geometry, namely polymers on the surface of an infinitely long cylinder, where the surface tension should make yet another contribution. Our main motivation was originally to have a simple geometry with straight “surfaces” so that the surface tension and its $`T`$-dependence can be more easily measured. But we found also another interesting effect: the bulk fugacity $`\mu `$ depends non-monotonically on the cylinder radius. It has a minimum at a finite radius $`R^{}`$ which diverges for $`TT_\theta `$. Thus a collapsed polymer, when wrapped around an elastic cylinder, will tend to squeeze it when $`R>R^{}`$ and widen it when $`R<R^{}`$.

A collapsed polymer on the surface of an infinitely long cylinder is modeled as an interacting self-avoiding random walk (ISAW) on a square lattice with periodic boundary conditions in the vertical direction, as shown in Fig. 1. The vertical lattice size is $`h=2\pi R`$. The collapse is mediated by nearest neighbour attraction between monomers, described by a Boltzmann factor $`q=e^{\beta ϵ}`$ for each pair of non-bonded monomers occupying nearest neighbour sites on the lattice. Here $`\beta =1/kT`$ and $`ϵ<0`$ is the attractive potential. The partition sum is therefore

$$Z=\underset{walks}{}q^m,$$
(2)
where $`m`$ is the number of non-bonded occupied nearest-neighbour pairs. We are interested in polymers which form collapsed globules, therefore $`q>q_\theta `$ with $`q_\theta =e^{ϵ/kT_\theta }1.95`$ barkema .

As argued in Grosberg ; Lifshitz  and verified in Owczarek ; g02 , below the $`\mathrm{\Theta }`$ point, the polymer forms a globule with uniform monomer density $`\rho `$ for large $`N`$ and $`\rho 0`$ as $`TT_\theta `$. Thus, the free energy for a collapsed polymer with chain length $`N`$ in infinite volume consists (except for a logarithmic term resulting from the last factor in Eq. (1)) of two parts: An extensive bulk contribution $`N`$, and a surface contribution $`N^{(d1)/d}`$. For $`d=2`$, this reads (up to a factor $`1/T`$ which will be suppressed in the following)

$$\mathrm{ln}Z_N(q)\mu _{\mathrm{}}(q)N+\stackrel{~}{\sigma }(q)\sqrt{N}+O(\mathrm{ln}N)$$
(3)
where $`\mu _{\mathrm{}}`$ is the chemical potential per monomer in an infinite chain, $`\stackrel{~}{\sigma }`$ is related to the surface tension (=free energy per unit of perimeter length) $`\sigma `$ and to the monomer density $`\rho `$ by

$$\sigma =\frac{1}{2}\stackrel{~}{\sigma }\sqrt{\frac{\rho }{\pi }}.$$
(4)
The latter follows by assuming the globule to be a disk with radius $`R_g`$ fixed by $`N=\pi R_g^2\rho `$.

Surface terms are absent from the free energy if one has a finite $`L\times L`$ lattice with periodic boundary conditions, and if one chooses the polymer length such that the entire lattice is filled with finite density. More precisely, consider the quantity $`\mathrm{ln}Z_N(q)+\mu N`$. For $`T<T_\theta `$ this is non-convex, and for $`\mu =\mu _L(q)`$ it has two peaks of equal height, corresponding to the coexistence of two phases: One phase with very short chains and density $`0`$, and another phase with density $`\rho _L(q)`$. For $`L\mathrm{}`$, $`\mu _L(q)`$ and $`\rho _L(q)`$ tend monotonically g02  to $`\mu _{\mathrm{}}(q)`$ and $`\rho `$.

Finally, on the surface of an infinitely long cylinder with perimeter $`h`$, a collapsed polymer forms a cylindrical blob with length $`L=\rho ^1N/h`$ and surface free energy $`2\sigma h`$,

$$\mathrm{ln}Z_N(q,h)\mu _h(q)N+2\sigma (q)h.$$
(5)
Notice that the chemical potential per monomer will in general depend on $`h`$, even in the limit $`N\mathrm{}`$. For $`h\mathrm{}`$ we expect again to obtain the bulk limit,

$$\mu _h(q)\mu _{\mathrm{}}(q)h\mathrm{}.$$
(6)
In principle we should also expect $`\sigma (q)`$ and $`\rho (q)`$ to depend on $`h`$. But we shall not study the detailed behaviour of $`\rho `$, and we shall neglect the $`h`$-dependence of $`\sigma (q)`$, since the surface term is a small correction anyhow. We also might expect logarithmic terms in $`N`$ to be present (just as also in the previous geometries), but they will also be neglected.

Our aim is to verify Eq. (5) numerically, to estimate the functions $`\mu _h(q)`$ and $`\sigma (q)`$, and to see whether they are compatible with the results for the other two geometries.

The pruned-enriched-Rosenbluth method (PERM) has been employed to study collapsed 2-d polymers in barkema ; g02 . In g02  the predicted behaviour was verified numerically for globules on infinite lattices and for polymers filling finite lattices with periodic boundary conditions. More precisely, finite lattices were studied first, and from them estimates were obtained for $`\mu _{\mathrm{}}`$ and $`\rho `$ by extrapolation. Given these, subsequent simulations on infinite lattices gave rather precise estimates for $`\sigma (q)`$. But the latter was not done systematically for a wide range of temperatures, in particular the behaviour of $`\sigma (q)`$ for $`qq_\theta `$ was not studied. Thus we want to obtain also the latter in the present paper from cylinder geometry simulations.

For the first simulations, we choose $`q=2.4`$ which is deep inside the collapsed regime. From Eq. (5) we expect that plots of $`\mathrm{ln}Z_N(q,h)+\mu N`$ against $`N`$ should give horizontal lines for $`N\mathrm{}`$, if and only if $`\mu =\mu _h(q)`$. This gives estimates of $`\mu _h(q)`$. Using them, we plot in Fig. 2 $`\mathrm{ln}Z_N(q,h)+\mu _h(q)N`$ versus $`N`$ for various values of the cylinder circumference $`h`$. For small $`h`$ we find perfect agreement with our expectations. For larger $`h`$ we see that PERM has problems in sampling long chains correctly, as should be expected. Fluctuations similar to those seen in the curve for $`h=16`$ prevent us from going to much larger $`h`$, except for $`q`$ close to $`q_\theta `$ where the efficiency of PERM increases.

Values of $`\mu _h(q)`$, again for $`q=2.4`$, are plotted in Fig. 3. For large $`h`$, we obtain $`\mu _{\mathrm{}}(q)=1.3211(1)`$ in perfect agreement with g02 . But, surprisingly, the dependence of $`\mu _h(q)`$ on $`h`$ is not monotonic. Instead, we see a minimum at a finite value of $`h`$. In our simulations $`h`$ is an integer, but we can interpolate to continuous values by fitting a parabola, $`\mu _h(q)=\mu _{\mathrm{min}}+b(hh^{})^2`$, through the three lowest points, as shown in Fig. 3. Such a fit gives $`h^{}(q=2.4)=6.13(1)`$ and $`\mu _{\mathrm{min}}(q=2.4)=1.3243(1)`$.

Finally, in order to obtain the surface tension, we plot in Fig. 4

$$A(q,h)=\underset{N\mathrm{}}{lim}[\mathrm{ln}Z_N(q,h)+\mu _h(q)N]$$
(7)
against $`2h`$. The limit is of course obtained by using the parts of the curves in Fig. 2 where they are horizontal. From Eq. (5) we expect this to give a straight line with slope $`\sigma `$. The result shows that $`\sigma =0.105(5)`$ at $`q=2.40`$. Actually, the line is slightly curved, showing that systematic errors are not negligible. For large $`h`$ they seem to arise from insufficient sampling. For small $`h`$ they result from the neglected terms discussed above. The present estimate should be compared to the estimate $`\sigma =0.119(5)`$ from g02 . Although the discrepancy is about 2 standard deviations, we consider this as reasonable agreement. Results for $`q=2.0`$ are also shown in Fig. 4, and illustrate the large corrections near the $`\mathrm{\Theta }`$-point.

Analogous simulations as shown in Fig. 2-4 were also done at other values of $`q`$. The lengths $`N`$ depended strongly on $`q`$. Close to the $`\mathrm{\Theta }`$-point the algorithm is very efficient and we obtained reliable high statistics data for $`N`$ up to 15,000. The more the polymer is collapsed, the more difficult are the simulations. We stopped at $`q=2.4`$. As seen from Fig. 2, the errors decrease when $`h`$ becomes small, most likely because the walk is effectively Markovian with memory $`h^2`$, and PERM is most efficient when memory is short. In all those simulations the results were similar. $`\mu _h(q)`$ had a minimum at a finite $`h^{}(q)`$, and $`A(q,h)`$ versus $`h`$ was roughly linear with a negative slope. Results of $`h^{}(q)`$, $`\mu _{\mathrm{min}}(q)`$, $`\rho (q)`$ and $`\sigma (q)`$ for a wide range of $`q`$ are shown in Fig. 5-8.

In Fig. 5 we can see that $`h^{}(q)`$ increases as $`q`$ decreases. This is to be expected since the $`\mathrm{\Theta }`$ transition is continuous, with a diverging correlation length for $`TT_\theta `$. Apart from the lattice constant (which is set to 1 here), the correlation length is the only relevant length scale. The divergence of $`h^{}(q)`$ when approaching the $`\mathrm{\Theta }`$ point thus indicates that the existence of the minimum is not a lattice artifact and is presumably a genuine feature of the continuum theory. The simulations are not precise enough to allow a meaningful fit to a power law and a determination of a critical exponent. This exponent would also be effected strongly by the uncertainty of $`T_\theta `$. To reduce the latter to a minimum, we performed extensive simulations of free (i.e. infinite lattice) ISAWs near the estimated $`\mathrm{\Theta }`$ point, with the result $`q_\theta =1.9487(5)`$ which is in good agreement with the best previous estimate barkema .

In Fig. 6 we plot $`\mu _{\mathrm{min}}(q)`$ and $`\mu _{\mathrm{}}(q)`$ versus $`(qq_\theta )`$. We see that both decrease linearly with temperature. The difference between them depends weakly on the temperature. It seems to vanish for $`TT_\theta `$.

The monomer density $`\rho `$ inside the “globule” is constant for $`Nh^21`$ and fixed temperature. It can be measured by measuring the $`x`$-coordinate (parallel to the cylinder) of the end point of the chain. Since the chain has to grow either into the positive or into the negative direction, we have

$$\frac{d|x_N|}{dN}(h\rho )^1,N\mathrm{}.$$
(8)
For small values of $`h`$ there are finite size corrections to the density, but it is easy to extrapolate to $`h\mathrm{}`$. The resulting densities are plotted versus $`q`$ in Fig. 7. They show a power law

$$\rho (qq_\theta )^\beta $$
(9)
with $`\beta =0.32\pm 0.02`$. We are not aware of previous estimates of $`\beta `$, either theoretically or by simulations.

The values of the surface tension $`\sigma (q)`$ obtained by the same way as shown in Fig. 4 are shown in Fig. 8. They follow a power law, $`\sigma (qq_\theta )^\alpha `$, with $`\alpha 1.7\pm 0.1`$. For the parameter $`\stackrel{~}{\sigma }`$ this implies that

$$\stackrel{~}{\sigma }(qq_\theta )^{\alpha \beta /2}(qq_\theta )^{1.54\pm 0.1}.$$
(10)
This is incompatible with the theoretical prediction $`\stackrel{~}{\sigma }(T_\theta T)^{7/6}`$ quoted in nidras .

In summary, we have studied collapsed polymers on the surface of a cylinder. Using the PERM algorithm we could simulate chains of lengths between $`3000`$ (far below $`T_\theta `$) and $`10^4`$ (near $`T_\theta `$) for a wide range of temperatures. We have verified that the concept of a surface tension applies to 2D collapsed polymers, and have determined its scaling near the $`\mathrm{\Theta }`$ point. We also determined the scaling of the monomer density in the collapsed phase, as $`TT_\theta `$. In addition, we found a surprising non-monotonic dependence of the chemical potential on the perimeter of the cylinder. It shows that for any temperature there is a special value of the perimeter where the free energy of the polymer is minimal. The fact that this value diverges when the $`\mathrm{\Theta }`$-temperature is approached suggests that this is not an artifact of the square lattice but a generic feature.

Acknowledgements: We thank Dr. Walter Nadler for valuable discussions. Parts of the simulations were done on a PC-farm in Taipei. H. P. thanks the Computing Centre of Academia Sinica in Taipei, Taiwan for providing the computing facilities.