text:
# Phases of a conserved mass model of aggregation with fragmentation at fixed sites

## I Introduction

Reaction-Diffusion systems form an important class of nonequilibrium systems whose dynamics and steady state depend on various factors such as the nature of the reaction (aggregation, annihilation, birth, fragmentation), number or type of reactants involved (single or multi-species), velocity of reactants (ballistic or diffusion controlled) and presence of external input (injection)  . An interesting class with wide ranging applications involves the elementary moves of aggregation (coalescence on contact) and fragmentation (break-up of clusters of masses), besides diffusion. A number of analytical results including the occurrence of nonequilibrium phase transitions have been obtained for such systems with translationally invariant geometries . A natural question arises – what is the effect of quenched disorder on the possible phases of such systems and the transitions between them? We may anticipate interesting effects, as quenched disorder is known to strongly influence the character of the steady state in other nonequilibrium systems . For instance, the steady state of driven diffusive systems on a one-dimensional lattice with bias shows phase separation in the presence of even a single defect .

To see what effect quenched disorder might have in a system with diffusion and aggregation, consider the process of polymerisation in a random medium with traps at certain fixed sites in which the polymer can get stuck. Aggregation occurs when two diffusing chains meet and coalesce; the reduced mobility of the aggregates at the trap sites promotes the formation of large, *localised* aggregates at such sites. If these traps are not perfect and allow monomers to detach and leave, there is also a possibility of formation of large *mobile* aggregates in the bulk.

To elucidate under what circumstances which types of aggregates may form as a result of these physical effects, namely, diffusion-aggregation, chipping and trapping, we consider a simple, reduced model in which the mass leaves as a single aggregate to a nearest neighbour from all sites except at certain, fixed sites from which it is allowed to leave only by single monomer dissociation. We refer to the loss of unit mass as chipping and these special sites as chippers. The total mass in the system is conserved as is evident from the dynamical rules described above. The quenched character of the chipper sites is important in determining the steady state of this model which is quite different from that in the translationally invariant uniform-chipping model where chipping can occur at every site .

As a first step towards understanding such a spatially inhomogeneous system, we study in detail the case when only a single chipper is present. Despite the simplicity of the model, we find that that the presence of a chipper gives rise to interesting steady states. We study the model both analytically and numerically in arbitrary dimensions, both with and without a global bias which sets up an overall mass current. The constraint of conservation of total mass in the model plays an important role in determining the steady state. We find that the system may exist in three possible phases. Depending on the dimension and the presence of drive, the system either stays in one of the phases or else make transitions from one phase to the another as the parameters are varied. These phases are characterised by the presence or absence of aggregates with mass proportional to the volume of the system. In the thermodynamic limit, the mass of such an aggregate diverges, so we refer to it as an infinite aggregate. The three phases are described below:

Pinned Aggregate phase (PA phase): In this phase, with probability one, an infinite aggregate occurs only at the chipper site and not in the bulk. This infinite aggregate acts as a spatially localised particle bath for rest of the system.

Unpinned Aggregate phase (UA phase): In this phase, a localised infinite aggregate can exist at the chipper site together with mobile infinite aggregates in the bulk. Interestingly, both types of infinite aggregate can exist simultaneously. In this sense, this state is different from the aggregate phases in translationally invariant systems such as the no-chipper limit of this model or the uniform chipping model, in which only one mobile infinite aggregate occurs at a time.

Non Aggregate phase (NA phase): This phase is characterised by the absence of an infinite aggregate anywhere in the system. The mass is spread out all over the system in clusters, each of which has a vanishing fraction of the total mass in the thermodynamic limit.

We analysed the system within a Mean Field Theory (MFT) allowing for the spatial dependence in the mass distributions. We find that the MFT predicts the existence of only the UA phase. On comparing this result with Monte Carlo simulations in $`d=1`$ and $`d=2`$, we find qualitative agreement except in the $`1d`$, biased case. In this exceptional case, there is a phase transition from the NA phase to the PA phase as the density is increased. This case is different from the rest due to an interesting interplay of the ballistic scale of motion and the diffusive scale of coalescence in one dimension.

In the model described above, an infinite aggregate is formed due to both interactions (coalescence of diffusing particles) and the presence of an inhomogeneity in the form of the chipper. It is useful to contrast the behaviour of this model with that of a model of noninteracting particles in which the aggregate is formed solely due to the presence of disorder. We solve this model in the presence of a single defect by mapping it to a Zero Range Process  and show that it exhibits a phase transition from the NA phase to the PA phase in all dimensions, irrespective of bias. This is in contrast to the aggregation model described earlier which does not show a phase transition in higher dimensions in the presence of single chipper.

Although we mainly discuss the case of a single chipper in this paper, we also discuss what happens in the presence of extensive disorder in both the models. We argue that in the presence of an extensive number of chippers, the interaction effects (i.e. coalescence) can be ignored in the aggregation model on large enough length and time scales and it behaves like the free-particle model. The latter model can be solved in the presence of extensive disorder and shows a phase transition in all dimensions for all bias from the PA phase to the NA phase as the density is decreased.

The remainder of the paper is organised as follows. We define the single-chipper aggregation model in Section II and discuss the possible phases on the basis of the conservation law. We analyse the system within a MFT in Section III and show that it predicts the occurrence of the UA phase. We also discuss the numerical results which support this broad conclusion in several cases. In Section IV we turn to the exceptional case namely the one dimensional, biased case. We present analytical results in the PA phase and numerical results in the NA phase and at the critical point. In Section V, we present the solution of the single-chipper free particle model. We also discuss the likely behaviour of the system with extensive disorder in Section VI. Finally we summarise our results in Section VII.

## II Single-Chipper Aggregation model and its phases

### A The Model

Our model is defined on a $`d`$-dimensional hypercubic lattice of length $`L`$ with periodic boundary conditions, and is studied both in the presence and absence of bias. The system evolves via following rules: at any site except the origin,

all the mass at that site moves as a whole to a nearest neighbour at a rate $`1`$, while at the origin only a single monomer chips off at a rate $`w`$ leaving rest of the mass behind. Once the mass moves from any site to its neighbouring site, it simply coalesces with the mass present on that site (see Fig. 1). Evidently, the total mass $`M=\rho L^d`$ of the system is conserved, where $`\rho `$ is the mass density. The evolution of the system can be described by the time evolution equation for the probability $`P_k(m,t)`$ that there is mass $`m`$ at site $`k`$ at time $`t`$. As the moves at the origin are different from that at the rest of the system, the mass distribution at the origin (denoted by $`0`$) and its nearest neighbours which receive mass from the origin (denoted by set $`A`$) obey equations different from those obeyed by rest of the system (set $`B`$).

The set of time evolution equations obeyed by the sites in the bulk (set $`B`$) are given by

$`{\displaystyle \frac{P_k(m,t)}{t}}`$ $`=`$ $`{\displaystyle \underset{\delta }{}}[{\displaystyle \underset{m^{}=1}{\overset{m}{}}}P_{k,k+\delta }(mm^{},m^{},t){\displaystyle \underset{m^{}0}{}}P_{k,k+\delta }(m,m^{},t)P_k(m,t)],kB`$ (2)
$`{\displaystyle \frac{P_k(0,t)}{t}}`$ $`=`$ $`{\displaystyle \underset{\delta }{}}[{\displaystyle \underset{m0}{}}P_k(m,t){\displaystyle \underset{m^{}0}{}}P_{k,k+\delta }(0,m^{},t)],kB,`$ (3)
where $`P_{k,k+\delta }(m,m^{},t)`$ is the joint probability that site $`k`$ and its nearest neighbour $`k+\delta `$ have mass $`m`$ and $`m^{}`$ respectively. The index $`\delta `$ runs over the $`2d`$ nearest neighbours in the unbiased case and $`d`$ forward neighbours for totally biased case. The convolution term in Eq.(2) is the gain term through which the deficit mass is supplied to site $`k`$ via diffusion-aggregation moves from its nearest neighbours. The gain term for zero mass in Eq.(3) results from hopping the mass from site $`k`$ to one of its nearest neighbours. The system can get out of the configuration in which the site $`k`$ has mass $`m`$ (including zero) when (a) the mass hops into site $`k`$ from its nearest neighbours and (b) the nonzero mass at site $`k`$ hops to its nearest neighbours.

For the sites in the immediate neighbourhood of the origin (set $`A`$), the time evolution equation for $`P_k(m,t)`$ is similar to that for set $`B`$ except that the contribution of the origin needs to be taken into account separately. We have

$`{\displaystyle \frac{P_k(m,t)}{t}}`$ $`=`$ $`{\displaystyle \underset{\delta }{}}{\displaystyle {}_{}{}^{}\underset{m^{}=1}{\overset{m}{}}}P_{k,k+\delta }(mm^{},m^{},t)+w{\displaystyle \underset{m^{}0}{}}P_{k,0}(m1,m^{},t){\displaystyle \underset{\delta }{}}{\displaystyle {}_{}{}^{}\underset{m^{}0}{}}P_{k,k+\delta }(m,m^{},t)`$ (6)
$`w{\displaystyle \underset{m^{}0}{}}P_{k,0}(m,m^{},t){\displaystyle \underset{\delta }{}}P_k(m,t),kA`$
$`{\displaystyle \frac{P_k(0,t)}{t}}`$ $`=`$ $`{\displaystyle \underset{\delta }{}}{\displaystyle \underset{m0}{}}P_k(m,t){\displaystyle \underset{\delta }{}}{\displaystyle {}_{}{}^{}\underset{m{}_{}{}^{}0}{}}P_{k,k+\delta }(0,m^{},t)w{\displaystyle \underset{m^{}0}{}}P_{k,0}(0,m^{},t),kA,`$ (7)
where the primed sum denotes the sum over the nearest neighbours, excluding the origin. The contribution of the origin is taken care of by the terms with a coefficient $`w`$ which accounts for the gain or loss in mass via chipping from the origin.

At the origin itself, the set of evolution equations obeyed are given by

$`{\displaystyle \frac{P_0(m,t)}{t}}`$ $`=`$ $`{\displaystyle \underset{\delta }{}}[{\displaystyle \underset{m^{}=1}{\overset{m}{}}}P_{0,0+\delta }(mm^{},m^{},t)+wP_0(m+1,t){\displaystyle \underset{m^{}0}{}}P_{0,0+\delta }(m,m^{},t)wP_0(m,t)]`$ (9)
$`{\displaystyle \frac{P_0(0,t)}{t}}`$ $`=`$ $`{\displaystyle \underset{\delta }{}}[wP_0(1,t){\displaystyle \underset{m^{}0}{}}P_{0,0+\delta }(0,m^{},t)].`$ (10)
In the gain term for nonzero mass, besides the convolution term, there is an extra term due to the possibility of chipping off one extra monomer with a rate $`w`$. The loss terms are similar to those discussed for the sites in set $`B`$.

### B The Phases

One can obtain useful information about the nature of the steady state from the following simple analysis. In the steady state, the mass current into and out of any site must be equal, i.e.

$`{\displaystyle \underset{\delta }{}}m_k`$ $`=`$ $`{\displaystyle \underset{\delta }{}}m_{k+\delta },kB`$ (12)
$`{\displaystyle \underset{\delta }{}}m_k`$ $`=`$ $`ws_0+{\displaystyle \underset{\delta }{}}{}_{}{}^{}m_{k+\delta },kA`$ (13)
$`{\displaystyle \underset{\delta }{}}ws_0`$ $`=`$ $`{\displaystyle \underset{\delta }{}}m_{0+\delta },`$ (14)
where $`m_k`$ is the average mass at site $`k0`$ and $`s_0=1P_0(0)`$ is the probability that the origin is occupied. The solution of the above equations gives $`m_k=ws_0`$ for all $`k0`$. Thus the average mass on a site in the bulk is uniform in spite of broken translational invariance. Since the total mass of the system is conserved, the average mass at the origin $`M_0`$ must satisfy

$$M_0=\rho L^dws_0(L^d1).$$
(15)
Using this constraint equation, one can deduce some characteristics of the possible phases in the system as follows. Since the LHS in Eq.(15) is non-negative and $`s_0`$ is bounded above by one, the allowed domain for $`s_0`$ is constrained as shown in Fig. 2. Depending on the value of $`s_0`$ (and hence $`M_0`$), the system can exist in one of the three distinct phases described below.

PA phase: The system exists in this phase when $`s_0`$ is pinned to its maximum value $`1`$ for all $`\rho >w`$ (shown by the bold line in Fig. 2). From the constraint equation Eq.(15), one deduces that the average mass at the origin grows as $`L^d`$ giving rise to an infinite aggregate in the thermodynamic limit. Thus the Pinned Aggregate phase is defined as that in which an infinite aggregate is present at the origin with a probability one.

UA phase: The system exists in this phase for all values of $`s_0`$ lying in the shaded region in Fig. 2. This phase is characterised by $`s_0`$ strictly less than $`1`$, and not equal to $`\rho /w`$. As in the PA phase, the average mass at the origin grows as $`L^d`$ but $`s_0`$ is not fixed and varies with $`\rho `$. Since $`s_0<1`$, the infinite aggregate is present at the origin only for a finite fraction of the time. For this reason, we refer to this phase as the Unpinned Aggregate phase.

NA phase: This phase is characterised by $`s_0\rho /w`$ as $`L\mathrm{}`$ (shown by the solid line in Fig. 2). Using Eq.(15), one immediately obtains $`M_0/L^d0`$ in the thermodynamic limit. Due to the absence of an aggregate with mass of $`O(L^d)`$ at the origin, we call this the Non Aggregate phase.

A more detailed study of the system shows that depending on the dimension and drive, the system can either stay in one of the three phases or else make transitions from one phase to another as $`\rho `$ is varied, keeping $`w`$ fixed. We have solved for $`M_0`$ and $`s_0`$ in $`d`$ dimensions within a sitewise inhomogeneous Mean Field Theory (Fig. 2) which predicts that the system exists only in the UA phase in all dimensions, regardless of the bias. This prediction was tested numerically in several cases and seen to be qualitatively correct except in the $`1d`$, biased case.

In the latter case, there is a phase transition from the NA phase to the PA phase as one crosses the critical line $`\rho =w`$ either by increasing $`\rho `$ with $`w`$ held fixed or by decreasing $`w`$ for fixed $`\rho `$.

## III Unpinned Aggregate Phase: Mean Field Theory and Monte-Carlo simulations

In this Section, we analyse the aggregation model within the mean field approximation. We implement this approximation by replacing the joint probability distribution $`P_{k,k+\delta }(m,m^{},t)`$ by the product $`P_k(m,t)`$ $`P_{k+\delta }(m^{},t)`$ of the single site probability distributions. We begin with the MFT analysis in arbitrary dimensions in Section III A. We will see that MFT predicts only the existence of the UA phase in all dimensions, irrespective of bias. In Section III B, we study the UA phase in one dimension in more detail. We discuss some limitations of this MFT and also present some numerical results. We close this Section with a qualitative discussion of the UA phase.

### A MFT in arbitrary dimensions

In this Section, we solve for $`M_0`$ and $`s_0`$ in arbitrary dimensions and for all bias within MFT and show that it predicts the occurrence of only the UA phase in all cases.

In the steady state, the mass distribution at any site $`k`$ is independent of time and is determined by setting the LHS of Eqs.(2)-(10) equal to zero. Further we replace $`P_{k,k+\delta }(m,m^{})`$ by $`P_k(m)`$ $`P_{k+\delta }(m^{})`$ in Eqs.(2)-(10) to obtain the mean field equations for $`P_k(m)`$. To study these equations, we define $`Q_k(z)=_{m=0}P_k(m)z^m`$, the Laplace transform of $`P_k(m)`$ with respect to mass, and find that it obeys the following set of equations,

$`{\displaystyle \underset{\delta }{}}Q_k`$ $`=`$ $`{\displaystyle \underset{\delta }{}}\left[1+Q_k(Q_{k+\delta }1)\right],kB`$ (17)
$`Q_k`$ $`=`$ $`{\displaystyle \frac{_\delta 1}{_\delta 1+ws_0(1z)_\delta {}_{}{}^{}(Q_{k+\delta }1)}},kA`$ (18)
$`Q_0`$ $`=`$ $`{\displaystyle \frac{_\delta w(1z)(1s_0)}{_\delta [w(1z)+z(Q_{0+\delta }1)]}}.`$ (19)
Since $`P_k(m)`$ is normalised to one, $`Q_k(1)=1`$ for all $`k`$ in the above equations. The average mass in the bulk $`m_k=Q_k^{}(z)|_{z=1}`$ still obeys Eqs.(12)-(14) which gives uniform mass in the bulk $`m=ws_0`$, and Eq.(15) still holds.

Although we could not solve the above set of nonlinear coupled equations, we were able to compute two quantities of primary interest, namely, $`M_0`$ and $`s_0`$. We begin by observing that the average mass at the origin can be written in terms of the mean-squared mass at its nearest neighbours,

$$M_0=\frac{_\delta \left(2ws_0+m_{0+\delta }^2m\right)}{_\delta 2w(1s_0)},$$
(20)
where we have used that $`M_0=Q_0^{}(z)|_{z=1}`$ and $`m_k^2=Q_k^{\prime \prime }(z)|_{z=1}+m`$, $`k0`$.

The mean-squared mass at sites other than the origin obeys the following set of linear equations,

$`{\displaystyle \underset{\delta }{}}m_k^2`$ $`=`$ $`{\displaystyle \underset{\delta }{}}\left(m_{k+\delta }^2+2m^2\right),kB`$ (22)
$`{\displaystyle \underset{\delta }{}}m_k^2`$ $`=`$ $`{\displaystyle \underset{\delta }{}}{}_{}{}^{}(m_{k+\delta }^2+2m^2)+2m^2+m,kA.`$ (23)
It suffices to calculate $`_\delta m_{0+\delta }^2`$ in order to obtain $`M_0`$. Adding the above $`L^d1`$ equations , one obtains

$$\underset{\delta }{}m_{0+\delta }^2=\underset{\delta }{}\left(2m^2(L^d1)+m\right),$$
(24)
which further yields

$`M_0`$ $`=`$ $`{\displaystyle \frac{ws_0+m^2(L^d1)}{w(1s_0)}}=\rho L^dws_0(L^d1),`$ (25)
where the mass conservation equation Eq.(15) has been used in the last identity. Solving for $`s_0`$ and $`M_0`$ in terms of $`\rho ,w`$ and $`L`$, we find

$`M_0`$ $`=`$ $`{\displaystyle \frac{\rho L^d(1+\rho L^d)}{(1w)+(\rho +w)L^d}}`$ (27)
$`s_0`$ $`=`$ $`{\displaystyle \frac{\rho L^d}{(1w)+(\rho +w)L^d}}.`$ (28)
To lowest nonvanishing order in $`L`$, the above equations give $`M_0L^d`$ and $`s_0<1`$ which are the signatures of the UA phase. Thus MFT predicts that the system exists only in the UA phase for all $`\rho `$ and $`w`$ (see Fig. 2). Our numerical simulations show that this prediction is true at least qualitatively in all cases except in the $`1d`$, biased case. In Fig. 3, we show our simulation results for $`M_0`$ in the absence of bias in $`1d`$ and $`2d`$ plotted along with the MFT results, and find qualitative agreement. We also measured $`s_0`$ for various system sizes and densities, fixing $`w=1`$ and find that it is independent of $`L`$ in accordance with Eq.(28). For $`\rho =0.5`$, we found that $`s_0^{num}0.42`$ to be compared with $`s_0^{MFT}0.33`$; for $`\rho =1.5`$, numerically $`s_0^{num}0.84`$ whereas $`s_0^{MFT}=0.6`$.

### B UA phase in one dimension

In this Section, we analyse the UA phase in more detail in the $`1d`$, unbiased case. We solve Eq.(22) and Eq.(23) for the mean-squared mass in the bulk in this case and find that the MFT seems to violate the conservation law. We comment on this limitation of MFT. We also present numerical evidence which shows the simultaneous presence of more than one infinite aggregate in the system.

For convenience, we will choose $`L`$ to be odd. The mean-squared mass at sites other than the origin obeys the following equations,

$`2m_k^2`$ $`=`$ $`m_{k1}^2+m_{k+1}^2+4m^2,|k|2`$ (30)
$`2m_1^2`$ $`=`$ $`m_2^2+m+4m^2`$ (31)
$`2m_1^2`$ $`=`$ $`m_2^2+m+4m^2.`$ (32)
This set of equations can be solved and one obtains

$`m_k^2`$ $`=`$ $`2m^2(L1)+2m^2(|k|1)(L|k|1)+m,k0.`$ (33)
To leading order in $`L`$, the above result can be written in the scaling form $`m_k^2=L^2f(|k|/L)`$ where $`f(x)=x(1x)`$ .

The mean squared mass $`m_k^2`$ is actually a measure of the typical mass at site $`k`$. To see this, note that we may define two types of average over mass distributions at site $`k`$, namely, an average $`\mathrm{}`$ over all mass occupations including $`m_k=0`$, and an average $`\mathrm{}`$ over mass occupations excluding $`m_k=0`$. It is straightforward to see that $`\mathrm{}=\mathrm{}/s_k`$ where $`s_k=1P_k(0)`$ is the probability that site $`k`$ is occupied. For instance, the average mass $`m_k`$

is same on all sites while the typical mass $`m_k=m_k/s_k`$ will depend on $`k`$. Similarly, $`m_k^2`$ is the square of the typical mass at site $`k`$, implying that $`m_k^2=m_k^2s_k`$ may be interpreted as a measure of the typical mass at site $`k`$.

We performed numerical simulations to test the mean field prediction for the scaling form $`m_k^2=L^2f(|k|/L)`$ where $`f(x)=x(1x)`$. As shown in Fig. 4, we obtain a data collapse with the scaling form $`m_k^2=Lg(|k|/L)`$ with $`g(x)\sqrt{x(1x)}`$ as the numerical fit for the scaling function. Since $`g(x)\sqrt{x}`$ for $`x`$ close to zero and is constant for $`x1/2`$, the typical mass scales as $`\sqrt{L}`$ for $`k`$ close to the origin but as $`L`$ at sites diametrically opposite to the origin. This points to the existence of an aggregate in the bulk as well, consistent with the nomenclature Unpinned Aggregate Phase.

One notes that the MFT solution seems to violate the conservation law since it predicts typical masses of $`O(L^2)`$ at sites situated $`O(L)`$ away from the chipper. On the other hand, one can also check that Eq.(15) is true using the MFT results for $`M_0`$ and $`s_0`$. Thus MFT seems to be able to describe the vicinity of the chipper more correctly than the bulk. The reason for this could be that the mass fluctuations about the mean near the chipper are small due to fragmentation unlike those in the bulk where only the aggregation move operates. At a distance of $`O(L)`$ away from the chipper, to a good approximation, one can neglect the presence of the chipper. Then as all the mass resides only on one site, the mass-mass correlation function at two different sites $`m_im_j`$ is exactly zero in the steady state, in strong contrast to the mean field approximation $`m_im_j=m^2`$. Thus a more refined approximation is required in the regions where the aggregation move dominates.

In Fig. 5, we present the numerical evidence which indicates the simultaneous presence of more than one infinite aggregate in the system. We measured the conditional average masses $`M_{max}^{(1)}`$ and $`M_{max}^{(2)}`$ of the clusters of the largest and second largest mass in the bulk respectively, given there is a mass $`M_0`$ at the chipper. We find that data collapse is obtained with the scaling forms $`M_{max}^{(1)}=Lf_1(M_0/L)`$ and $`M_{max}^{(2)}=Lf_2(M_0/L)`$, where the scaling functions $`f_1(x)`$ and $`f_2(x)`$ decay almost linearly. Thus a localised infinite aggregate at the chipper and one or more mobile infinite aggregates in the bulk can be present at the same time; we would expect more than one infinite aggregate to be present in higher dimensions as well.

In any dimension $`d`$, far away from the chipper site, the state resembles that in the absence of the chipper. In the latter case, there is a single mobile infinite aggregate which is equally likely to be present at any site with a probability $`1/L^d`$. Although this mobile aggregate with mass of order $`L^d`$ arrives at the chipper infrequently, the probability of occupation $`s_0`$ of the chipper is of $`O(1)`$ – this enhancement occurs because mass can leave only one unit at a time, so that it stays for a time of at least order $`L^d`$. It is implicit in the above argument that an infinite aggregate can be formed in the bulk before it hits the chipper. However, this fails to be true in the $`1d`$, biased case which explains the absence of the UA phase in that case.

## IV ASYMMETRIC SINGLE-CHIPPER AGGREGATION MODEL IN ONE DIMENSION

The mean field prediction that the system exists in the UA phase fails in $`1d`$, biased case. In the presence of a drift velocity, the one-dimensional system undergoes a phase transition from the NA phase to the PA phase as $`\rho `$ is increased, keeping $`w`$ fixed. This exceptional case is the subject of this Section. Although the phase transition survives for all nonzero bias, we will only discuss the extreme case when the mass moves only forward (i.e. infinite bias).

The time required to form an aggregate with mass of $`O(L)`$ in the bulk is $`O(L^2)`$. But in this case, due to ballistic motion, the mass clusters return to the chipper in time of $`O(L)`$ ruling out the formation of an infinite aggregate in the bulk. For small $`w`$, the sublinear mass arriving at the chipper cannot leave it easily and is temporarily trapped giving rise to a localised infinite aggregate (PA phase). As $`w`$ is increased, the mass leaves more frequently rendering the trapping less effective; also this chipped off mass cannot return before the chipper gets empty so that an infinite aggregate cannot be sustained at the chipper for large $`w`$ (NA phase). Thus there is a phase transition in the $`1d`$, biased case as $`w`$ (or alternatively $`\rho `$) is varied.

The critical point is determined exactly to be at $`\rho _c=w`$ in the thermodynamic limit by setting the LHS to be zero and $`s_0=1`$ in Eq.(15). For $`\rho <\rho _c`$, the system exists in the NA phase which is characterised by $`s_0=\rho /w`$ and $`M_0`$ growing sublinearly with $`L`$. For $`\rho >\rho _c`$, the system exists in the PA phase in which $`s_0`$ is pinned to its maximum value $`1`$ and $`M_0`$ grows linearly with $`L`$ (see Fig. 2). The probability $`P_0(0)=1s_0`$ that the origin is empty serves as an order parameter. Our numerical simulations indicate that to lowest order in $`L`$,

$`1s_0=P_0(0)\{\begin{array}{cc}O(e^{L/L_0})\hfill & ,\rho >w\hfill \\ O(L^{\beta 1})\hfill & ,\rho =w\hfill \\ \text{constant}\hfill & ,\rho <w.\hfill \end{array}`$ (37)
Since $`P_0(0)`$ varies continously from the NA phase to PA phase as $`\rho `$ increases, this is a second order phase transition.

Before delving into details, we first present a pictorial representation of the phases. Let us imagine monitoring the instantaneous mass at the chipper site in the steady state. Figure 6 shows the time series for the mass at the origin as a function of time in the two phases and at the critical point. The mass at the chipper increases due to the mass input from the site behind the chipper, and decreases due to fragmentation. Let $`t_0`$ and $`t_1`$ respectively denote the consecutive number of time steps during which the chipper is empty and occupied. Then $`s_0`$ can be related to $`t_0`$ and $`t_1`$ through $`s_0^11=t_0/t_1`$. Depending on whether the ratio $`t_0/t_1`$ is zero or not in the thermodynamic limit, $`s_0`$ is either pinned to $`1`$ (PA phase) or strictly less than $`1`$ (NA phase). In the PA phase, $`t_1/t_01`$ holds due to long cascades of successive mass inputs, thus enabling the origin to maintain $`s_0`$ equal to $`1`$. The long $`t_1`$ stretches (see Fig. 6(a)) enables the chipper to build an aggregate with mass of order of system size. In contrast, in the NA phase, $`t_0t_1`$ holds as depicted in Fig. 6(b) by the stretches of time during which the origin is empty (due to absence of multiple inputs) which reduces $`s_0`$ from its maximum value $`1`$. Since $`t_1`$ is typically not very long, the chipper cannot sustain an infinite aggregate in NA phase. Finally, at the critical point, there are multiple mass inputs to the chipper but these cascades are not very long as shown in Fig. 6(c). We now turn to a systematic discussion of each phase.

### A The Pinned Aggregate Phase

Since the total mass in the system is conserved and the $`L1`$ sites in the bulk have nonzero average mass $`m`$ and sublinear fluctuations about the mean mass (due to absence of a cluster of mass of order $`L`$ in the bulk as argued above), the probability $`P_0(m)`$ that the origin has mass $`m`$ is $`1/\sqrt{L}`$ exp$`((mm_0)^2/L)`$ where $`m_0MmL`$ which gives $`P_0(0)`$ exp$`(L)`$. Thus in this phase, the chipper site is occupied with probability one in the thermodynamic limit so that it acts as a reservoir of particles for the bulk. One can think of the system as a semi-infinite, one dimensional lattice with a perfect, localised source at the origin injecting monomers into it at the rate $`w`$.

The problem of aggregation in the presence of such a source has been considered previously as well. Some properties were studied in  using the technique of interparticle distribution function (IPDF) introduced in , and in  by mapping this $`1d`$ problem to a bounded random walk in $`2d`$. Here we calculate the steady state mass distribution at site $`k`$ denoted by $`P_k(m)`$ by the generating function method. We define the $`r`$-point characteristic function for site $`k`$ at time $`t`$, $`Z_r^{(k)}(\lambda ,t)=e^{\lambda _{j=k}^{k+r1}m_j(t)}`$. The time evolution equations obeyed by $`Z_r^{(k)}(\lambda ,t)`$ with a perfect, localised source at $`k=0`$ are given by

$`{\displaystyle \frac{Z_r^{(k)}}{t}}`$ $`=`$ $`Z_{r+1}^{(k1)}+Z_{r1}^{(k)}2Z_r^{(k)},k1,r0`$ (39)
$`{\displaystyle \frac{Z_r^{(1)}}{t}}`$ $`=`$ $`Z_{r1}^{(1)}+(we^\lambda 1w)Z_r^{(1)},r0`$ (40)
with the boundary condition $`Z_0^{(k)}(\lambda ,t)=1`$ for all $`k>0`$. In the steady state, $`Z_r^{(k)}(\lambda ,t)`$ is independent of $`t`$.

We need to solve for $`Z_1^{(k)}(\lambda )`$ which is the Laplace transform of $`P_k(m)`$ with respect to $`m`$. One can easily solve for $`Z_1^{(1)}(\lambda )`$ using Eq.(40) and the boundary condition $`Z_0^{(1)}=1`$ which on taking the inverse Laplace transform for large $`m`$ gives $`P_1(m)=e^{m/w}/w`$. Thus the mass distribution at the first site decays exponentially. Now to find the probability distribution in the bulk, we define

$$H(x,y,\lambda )=\underset{k=2,\mathrm{}}{}\underset{r=1,\mathrm{}}{}Z_r^{(k)}(\lambda )x^ky^r.$$
(41)
Using equations Eq.(39) and Eq.(40), one obtains

$`H(x,y,\lambda )`$ $`=`$ $`{\displaystyle \frac{xy}{y^22y+x}}\left[G(x,\lambda )xy\left({\displaystyle \frac{1}{1x}}+{\displaystyle \frac{f^2}{1fy}}\right)\right],`$ (42)
where $`1/f=1+wwe^\lambda `$ and $`G(x,\lambda )=_{k=2,\mathrm{}}Z_1^{(k)}(\lambda )x^k`$. Clearly, the inverse Laplace transform of $`G(x,\lambda )`$ w.r.t. $`\lambda `$ and $`x`$ gives $`P_k(m)`$ for $`k>1`$. To calculate $`G(x,\lambda )`$, we use the same method as in . We first note that the denominator of $`H(x,y,\lambda )`$ has roots at $`y_\pm =1\pm \sqrt{1x}`$. The root at $`y=y_{}`$ lies inside the unit circle which in the real space (i.e. $`r`$-space) gives an exponentially increasing solution which is disallowed since the inverse Laplace transform of $`H(x,y,\lambda )`$ gives a probability, which is always bounded above. To avoid this pathological solution, we demand that $`y=y_{}`$ is a zero of the numerator as well, i.e.

$$G(x,\lambda )=xy_{}\left(\frac{1}{1x}+\frac{f^2}{1fy_{}}\right).$$
(43)
For $`\lambda 0`$ and $`x1`$, the inverse Laplace transform can be easily found. We find that $`P_k(m)`$ has a scaling form $`\varphi (u)/m^2`$ where $`u=m/2w\sqrt{k}`$ and the scaling function is

$$\varphi (u)=\frac{4w}{\sqrt{\pi }}u^3e^{u^2}.$$
(44)
One can also obtain the probability distribution $`P_k(0)`$ of zero mass at site $`k`$ by solving for $`G(x,\lambda )`$ when $`\lambda \mathrm{}`$ and $`x1`$. To leading order in $`k`$, one obtains $`P_k(0)=1s_k=11/\sqrt{\pi k}`$. We tested our calculations against the numerical simulations performed on the single-chipper model and found reasonably good agreement (see Fig. 7). Thus the typical mass at site $`k`$ grows as $`\sqrt{k}`$.

We now calculate the typical spacing between the masses using the method of IPDF . We define $`E_k(r,t)`$ as the probability that the sites $`k`$ to $`k+r`$ are empty (including $`k`$ and $`k+r`$). Then $`E_k(r,t)`$ satisfies the following equations,

$`{\displaystyle \frac{E_k(r,t)}{t}}`$ $`=`$ $`E_k(r1,t)2E_k(r,t)+E_{k1}(r+1,t),k1`$ (46)
$`{\displaystyle \frac{E_1(r,t)}{t}}`$ $`=`$ $`E_1(r1,t)(1+w)E_1(r,t),`$ (47)
with the boundary condition $`E_k(1,t)=1`$ for all $`k>0`$. In the steady state, $`E_k(r,t)`$ is independent of $`t`$ and one can write the above equation in terms of $`F_1(z,y)=_{k=2,\mathrm{}}_{r=0,\mathrm{}}E_k(r)z^ky^r`$ and $`R(z,r)=_{k=2,\mathrm{}}E_k(r)z^k`$ , $`r=0,1,\mathrm{}`$. We find

$`F_1(z,y)`$ $`=`$ $`{\displaystyle \frac{y}{y^22y+z}}\left[{\displaystyle \frac{zR(z,0)}{y}}{\displaystyle \frac{z^2}{(1+w)(1+wy)}}{\displaystyle \frac{z^2}{1z}}\right].`$ (48)
With the same reasoning as in the previous calculation, we demand that the numerator of $`F_1(z,y)`$ evaluated at $`y_{}=1\sqrt{1z}`$ be zero. This condition gives

$$R(z,0)=\left[\frac{1}{1z}+\frac{1}{(1+w)(1+wy_{})}\right]zy_{},$$
(49)
which on inverting the Laplace transform gives the probability $`E_k(0)`$ which is same as $`P_k(0)`$. One can check that the result for $`P_k(0)`$ quoted above is reproduced. Substituting $`R(z,0)`$ in $`F_1(z,y)`$, we obtain

$$F_1(z,y)=\frac{z^2}{yy_+}\left[\frac{1}{1z}+\frac{1}{(1+wy)(1+wy_{})}\right],$$
(50)
where $`y_+=1+\sqrt{1+z}`$.

We further define $`D_k(r)`$ as the probability that both sites $`k`$ and $`k+r`$ are occupied but no sites in between. Then

$`D_k(r)`$ $`=`$ $`E_{k+1}(r2)E_{k+1}(r1)E_k(r1)+E_k(r),k1,r0.`$ (51)
Defining the Laplace transform of $`D_k(r)`$ with respect to $`k`$ and $`r`$ as $`F_2(z,y)=_{k=2,\mathrm{}}_{r=0,\mathrm{}}D_k(r)z^ky^r`$, we obtain

$`F_2(z,y)`$ $`=`$ $`{\displaystyle \frac{1}{z}}\left(R(z,2)+{\displaystyle \frac{yz^2}{1z}}\right){\displaystyle \frac{z(1+z)}{1z}}+F_1(z,y)\left({\displaystyle \frac{y^2}{z}}+(1y)(1+{\displaystyle \frac{1}{z}})\right)`$ (53)
$`+z(1y)\left(1+{\displaystyle \frac{y}{2y}}\left\{1+{\displaystyle \frac{1}{(1+w)(1+wy)}}\right\}\right),`$
where $`F_1(z,y)`$ was calculated above.

Then the typical empty space $`r_k`$ in front of site $`k`$ can be obtained by taking the inverse Laplace transform of $`(F_2/y)|_{y=1}`$ with respect to $`z`$. One finds

$$r_k=1+\frac{2}{\sqrt{\pi k}}we^{w^2k}\text{erfc}(w\sqrt{k}),$$
(54)
which for large $`k`$ implies that the typical empty space $`r_k=r_k/s_k`$ in front of occupied site $`k`$ varies as $`\sqrt{k}`$.

Thus, in the PA phase, the mass in the bulk is distributed in $`\sqrt{L}`$ clusters with typical mass $`\sqrt{k}`$ at site $`k`$, and an empty stretch of length $`\sqrt{k}`$ in front of it.

### B The Non Aggregate Phase

As discussed earlier, in the NA phase, $`s_0`$ is strictly less than $`1`$ due to substantial time stretches of typical length $`t_0`$ during which the origin is not occupied (Fig. 6). Thus the origin does not act as a perfect source unlike in the PA phase and we could not calculate the mass distribution in this phase. However, one can obtain useful information about the nature of this phase by simple scaling arguments. We begin by arguing that $`M_0`$ grows sublinearly with $`L`$ and then provide numerical evidence for the absence of an infinite aggregate (i.e. mass proportional to $`L`$) in the bulk, thus concluding that there is no infinite aggregate anywhere.

Monte Carlo simulations indicate that $`M_0`$, $`t_0`$ and $`t_1`$ vary as a power law in $`L`$ with a constant as a next order correction (see Fig. 9). Therefore, we assume $`M_0=aL^\beta +b`$, $`t_0=c_0L^\gamma +d_0`$ and $`t_1=c_1L^\eta +d_1`$. Solving for $`s_0`$ using Eq.(15), we obtain

$$s_0=\frac{\rho }{w}+O\left(\frac{1}{L^{1\beta }}\right).$$
(55)
Then since the lowest order term for $`s_0`$ is a constant and $`s_0^11=t_0/t_1`$, one obtains $`\gamma =\eta `$. Using this identity and retaining terms to lowest order in $`L`$ we further obtain

$$s_0=\frac{c_1}{c_0+c_1}+O\left(\frac{1}{L^\gamma }\right).$$
(56)
Comparison of Eq.(55) and Eq.(56) yields $`\beta +\gamma =1`$.

Now we consider the typical mass $`m_k`$ in the bulk where the site index $`k`$ is labelled by positive integers in the clockwise direction with the chipper at the origin. Numerically, $`m_k`$ is observed to obey the scaling form, $`m_k=\sqrt{k}f\left(k/L\right)`$ as shown in the data collapse in Fig. 8. The scaling function is a slowly varying function, which gives sublinear mass everywhere in the bulk. Thus in the NA phase, the mass is distributed in clusters of typical mass $`O(\sqrt{L})`$ for $`kO(L)`$.

Since the mass at the site just behind the chipper is of $`O(\sqrt{L})`$ and as one can see in Fig. 6, the number of mass inputs to the chipper is typically of $`O(1)`$, it follows that $`t_1O(\sqrt{L})`$ so that $`\eta =0.5`$, which further gives $`\beta =\gamma =\eta =0.5`$ due to the two scaling relations above. These exponent values are checked numerically, as shown in Fig. 9.

### C The Critical Point

We studied the critical point mainly numerically by studying the $`L`$ dependence of $`M_0`$, $`t_0`$ and $`t_1`$. Assuming a power law dependence for $`M_0L^\beta `$, $`t_0L^\gamma `$ and $`t_1L^\eta `$ , a naive best fit in the range $`L=32`$ to $`2048`$ gives $`\beta 0.62`$, $`\gamma 0.38`$ and $`\eta 0.80`$. However, the effective exponent calculated using the successive ratios of $`L`$ used in the simulations shows systematic trends: $`\beta `$ decreases as $`L`$ increases while both $`\gamma `$ and $`\eta `$ increase with $`L`$.

A good fit is obtained if one allows for logarithmic corrections in the power laws for $`M_0`$ and $`t_0`$ with $`\beta =\gamma =1/2`$ (see Fig. 10 and Fig. 11). This indicates that both $`M_0`$ and $`t_0`$ vary as $`\sqrt{L}`$ with multiplicative logarithmic corrections. Using Eq.(15) and the condition of criticality, one can easily show that $`\eta +\beta \gamma =1`$. Using $`\beta =\gamma =1/2`$ as suggested by the discussion above, we obtain $`\eta =1`$. The above quoted values for the exponents $`\beta `$ and $`\gamma `$ are consistent with the scaling relation $`\beta +\gamma =1`$ which seems to be true in the NA and the PA phase as well. In the NA phase, we have seen in Section IV B that this scaling relation holds. In the PA phase, we found numerically that $`t_0`$ is independent of $`L`$ in the thermodynamic limit implying that $`\gamma =0`$; further since $`\beta =1`$, it follows that the scaling relation $`\beta +\gamma =1`$ holds in this phase also.

We have also numerically studied $`P(t_1)`$ which is the probability that the origin is occupied for $`t_1`$ consecutive time steps. At the critical point, this distribution shows two peaks – a broad peak occurs at the typical time scale $`t_1L^\beta `$ while there is a narrow peak at $`t_1=N_p/w`$. We have not been able to reliably separate out these two contributions to $`P(t_1)`$ but since the scaling relation yields $`\eta =1`$, one may expect that the second peak dominates the asymptotic value of $`\eta `$.

## V Single-Chipper Free particle model

### A The Model

In this Section, we study the steady state of a model of noninteracting particles in which a localised infinite aggregate is formed solely due to disorder. In this model the particles diffuse freely in the bulk except at certain, quenched sites referred to as the chipper sites. It can be shown that this system exhibits a phase transition from the NA phase to the PA phase (as defined in Section II B) as the particle density is increased in all dimensions for all bias. In this Section, we will demonstrate this result in the presence of a single chipper. We will see that even a single defect is capable of inducing a phase transition in all dimensions, unlike in the aggregation model. We will discuss the solution with an extensive number of chipper sites in Section VI.

Our model is defined on a $`d`$-dimensional hypercubic lattice with periodic boundary conditions on which we consider the biased and unbiased diffusion of a conserved number of particles. At any site except the chipper, each particle attempts to hop out at rate $`1`$. Since the particles do not interact, the hopping rate out of site $`k`$ occupied by $`m`$ particles is $`m`$. At the chipper site, assumed to be at the origin, the hopping rate is a constant $`w`$, independent of the number of particles. Thus, the rate $`u_k(m)`$ at which a particle leaves site $`k`$ which has $`m`$ particles is given by

$$u_k(m)=w\delta _{k,0}+m(1\delta _{k,0}),m0.$$
(57)
In the fully asymmetric, $`1d`$ case, the model described above can be mapped onto a traffic model on a one lane road with no overtaking allowed. We represent each site in this model as a car and each particle as a vacant site. Then a system of $`M`$ free particles on a lattice of size $`L`$ maps onto a system of $`L`$ particles with hard core interactions on a lattice of size $`L+M`$. In this new representation, the special car (corresponding to the chipper site) moves with a constant rate, irrespective of the headway in front of it and rest of the cars (sites other than the chipper) move with a rate proportional to the headway in front of it. Note that the sitewise disorder in the free particle model corresponds to particlewise disorder in the traffic model.

### B Phase Transition in arbitrary dimensions

The steady state of this system can be found exactly by noting that the hopping rates Eq.(57) in this model correspond to a special choice of rates in the Zero Range Process . In this process, a particle at site $`k`$ hops to its nearest neighbour independent of the state at the target site, so that the interaction has zero range. The steady state measure of this process can be found in any dimension with or without bias. A convenient way to find it is by using the condition of pairwise balance  which states that for any given configuration $`C`$, one can find a configuration $`C^{}`$ in one-to-one correspondence with $`C^{\prime \prime }`$ such that

$$W(C^{}C)P(C^{})=W(CC^{\prime \prime })P(C),$$
(58)
where $`P(C)`$ is the probability of a configuration $`C`$ in a system with $`N`$ sites and $`W(CC^{})`$ is the transition rate from $`C`$ to $`C^{}`$. The above condition is satisfied by

$$P(C)=\frac{1}{𝒩}\underset{k=1}{\overset{N}{}}f_k(m),$$
(59)
where

$`f_k(m)`$ $`=`$ $`{\displaystyle \underset{i=1}{\overset{m}{}}}{\displaystyle \frac{1}{u_k(i)}}\text{for}m>0`$ (60)
$`=`$ $`1\text{for}m=0,`$ (61)
and $`u_k(m)`$ is the rate at which a particle hops out of a site $`k`$ having $`m`$ particles and $`𝒩`$ is the normalisation constant. $`f_k(m)`$ is defined upto a multiplicative factor $`v^m`$ where $`v`$ can be interpreted as the fugacity.

For our model, using the above solution for the particle distribution $`f_k(m)`$ and Eq.(57), we find that the normalised probability distribution $`P_k(m)`$ that the site $`k`$ in the bulk has $`m`$ particles has the Poisson form,

$$P_k(m)=P(m)=\frac{e^vv^m}{m!},k0.$$
(62)
The fugacity $`v`$ will be determined using the conservation law. At the chipper site, one can obtain solution for $`P_0(m)`$ when $`v<w`$,

$$P_0^{NA}(m)=\left(1\frac{v}{w}\right)\left(\frac{v}{w}\right)^m,\frac{v}{w}<1,$$
(63)
where NA in the superscript stands for the Non-Aggregate phase which we have added in the anticipation of a phase transition.

Using equations Eq.(62) and Eq.(63), one can easily see that the average particle number at site $`k0`$ is $`m_k=v=ws_0`$ which further leads to constraint equation Eq.(15). Since this constraint equation is identical to that in the aggregation model, then as discussed in Section II B, the steady state of this model also has three possible phases (see Fig. 2). But, as we will see, this system never exists in the UA phase and there is a phase transition from the NA phase to the PA phase as $`\rho `$ is increased, keeping $`w`$ fixed, in all dimensions, irrespective of bias.

Solving for the average number of particles at the origin $`M_0`$ using Eq.(63) and substituting in Eq.(15), we obtain

$$\frac{s_0}{1s_0}=\rho L^dws_0(L^d1),\mathrm{\hspace{0.33em}\hspace{0.33em}0}s_0<1.$$
(64)
To leading order in $`L`$, one obtains two solutions, namely, $`s_0=\rho /w`$ and $`s_0=1`$. Since the above analysis holds only when $`s_0`$ is strictly less than $`1`$, the only valid solution is $`s_0=\rho /w`$ for $`\rho <w`$. For $`\rho w`$, $`s_0`$ is pinned to its maximum value $`1`$. Thus there is a phase transition from the NA phase to the PA phase at $`\rho _c=w`$, as $`\rho `$ is increased, keeping $`w`$ fixed.

One can solve for $`P_0(m)`$ in the PA phase and at the critical point using the conservation law Eq.(15) and the fact that Eq.(62) is valid for all $`s_0`$ in the range . Since the total number of particles is conserved, $`P_0(m)=_CP(C)\delta (_{k0}m_k=Mm)`$. The quantity on the RHS can be calculated strightforwardly. The result is

$`P_0^{PA}(m)`$ $`=`$ $`\sqrt{{\displaystyle \frac{1}{\pi m^2(L^d1)}}}\text{exp}\left({\displaystyle \frac{(Mmm(L^d1))^2}{(L^d1)m^2}}\right),`$ (65)
where $`m`$=$`v`$ and $`m^2`$=$`v(1+v)`$ from Eq.(62). One can see from the above equation that $`M_0=Mw(L^d1)L^d`$ and $`P_0^{PA}(0)=1s_0e^{L^d}`$ in the PA phase. Further, at the critical point, we have $`\rho L^d=w(L^d1)`$ using which in Eq.(65), one obtains $`P_0(m)`$ at the critical point,

$$P_0^{cp}(m)=\sqrt{\frac{4}{\pi m^2(L^d1)}}\text{exp}\left(\frac{m^2}{(L^d1)m^2}\right),$$
(66)
which gives a power law decay in $`L`$ for $`P_0^{cp}(0)`$.

The phase transition is brought about in this model in the following way. In the NA phase which is a low density phase, the typical number of particles at all sites including the chipper is of $`O(1)`$. As the density is increased, there is a phase transition to the PA phase at $`\rho _c=w`$. In this high density phase, each site in the bulk still supports only $`O(1)`$ number of particles in accordance with Eq.(62) so that the extra particles condense on the chipper giving rise to $`M_0L^d`$. The mechanism of phase transition in this model is similar to that in Bose-Einstein condensation as was pointed out in  in a similar $`1d`$ model which shows a phase transition with $`u_k(m)=w\delta _{k,0}+(1\delta _{k,0})`$.

## VI Extensive disorder

In this Section, we describe a possible scenario in the more interesting and physically relevant situation when there is an extensive number of chipper sites. These are assumed to be placed randomly, with quenched random chipping rates $`w_k`$ at site $`k`$ distributed according to a distribution $`\text{Prob}(w_k)`$. We are interested primarily in the aggregation model but as argued below, on large space and time scales the behaviour of this model with extensive disorder resembles that of the corresponding free particle model (a generalization of the model in Section V with an extensive number of chippers) which is exactly solvable.

Let us consider the diffusion-aggregation process with no bias and take the initial condition to have a random distribution of masses. Then the finite concentration $`c`$ of chipper sites brings in new length and time scales into the problem, namely the mean spacing $`l_cc^{1/d}`$ between the chippers and the associated diffusion time $`t_c=l_c^2`$. Evidently, $`l_c`$ and $`t_c`$ define respectively the relevant length and time scales over which the diffusing mass clusters sense the presence of an extensive number of defects.

Let us consider a finite but low concentration of chipper sites so that $`1l_cL`$. On time scales $`tt_c`$, we would expect the system to behave roughly as independent, finite systems of size $`l_c`$ with typically a single chipper each. The typical state on these short time scales thus resembles the UA-like steady state discussed in Section III. The typical mass of clusters, both mobile and localised is limited by the time and grows proportional to $`t^{d/2}`$ as long as $`tt_c`$.

As the time crosses $`t_c`$, the finite spacing between the chippers becomes relevant. The primary effect is to limit the size of the aggregates formed in the non-chipper region by the diffusion-aggregation process to $``$ $`l_c^d`$, since on a time scale $`t_c`$ the mass cluster is likely to encounter a chipper and get trapped. For $`tt_c`$, the coarse grained view of the system is that of mass exchanged between close-by random rate chippers with each exchange taking a time $``$ $`t_c`$. To the extent that only finite aggregates (with mass $``$ $`l_c^d`$ $`L^d`$) are formed in the transit between chippers, it is plausible that on large length and time scales, we may ignore interaction effects (i.e. coalescence) and think of the system as effectively free-particle like in the non-chipper region.

The free particle model with extensive disorder is solvable along the same lines as the single-chipper problem described in Section V. It defines a Zero Range Process with the hopping rates,

$`u_k(m)`$ $`=`$ $`w_k,\text{if }k\text{ is a chipper site}`$ (67)
$`=`$ $`m,\text{if }k\text{ is not a chipper site}`$ (68)
Let $`x`$ denote the fraction of chipper sites. We recover the single-chipper model as $`x0`$, while $`x1`$ corresponds to every site being a chipper site, and is the model considered in . For all $`x`$, in the low-density phase, the mass distribution on the non-chipper sites follows the probability distribution of Eq.(62), while at chipper sites, Eq.(63) is valid with $`w`$ replaced by $`w_k`$. Let $`s_k`$ be the occupation probability of site $`k`$, and $`s_0`$ refer to the occupation probability of the site with the lowest chipping rate, $`w_0=Min\{w_k\}`$. Then in the steady state, the spatial uniformity of the current leads to (i) $`w_ks_k=w_0s_0`$ if $`k`$ is a chipper site, and (ii) $`m_k=w_0s_0`$ if $`k`$ is not a chipper site. Using these relations, one can write the mass conservation equation analogous to Eq.(64) as

$$\frac{1}{L^d}\frac{s_0}{1s_0}+x𝑑w\frac{\text{Prob}(w)s_0}{w/w_0s_0}+(1x)s_0w_0=\rho .$$
(69)
where we have separated out the first term corresponding to the slowest chipper.

One can analyse Eq.(69) for various $`x`$ as follows: (a) In the limit of a single chipper (reached as $`x`$ approaches $`0`$), we know from Section V that as the density is increased, there is a phase transition from the NA phase to the PA phase with an infinite aggregate at the chipper site. (b) In the limit $`x1`$, the model reduces to that considered in , where it is shown that for Prob$`(w_k)(w_kw_0)^\nu `$ as $`w_kw_0`$, the system stays in the NA phase for all $`\rho `$ if $`\nu 0`$. On the other hand, if $`\nu >0`$, then there is a transition to the PA phase with an infinite aggregate at the site with chipping rate $`w_0`$ when the density crosses the critical density given by

$$\rho _c(x=1)=𝑑w\frac{\text{Prob}(w)}{w/w_01}.$$
(70)
(c) For $`0<x<1`$, there is a transition from the NA to the PA phase if Prob$`(w)`$ is chosen as in (b) above. The critical density can be determined by taking $`s_01`$ and $`L^d(1s_0)\mathrm{}`$ in Eq.(69), with the result

$$\rho _c(x)=x\rho _c(x=1)+(1x)w_0.$$
(71)
Thus the critical density interpolates linearly between its values in the limits $`x=0`$ and $`x=1`$.

In view of the correspondence discussed at the beginning of this Section, we would expect the aggregation model with a similar distribution of chipper sites to show a phase transition from the NA to the PA phase. In the former phase, there are both localized and mobile aggregates with typical mass $`l_c^d`$. In the PA phase, the distribution of masses is similar to that in the NA phase at all sites except the slowest chipper with $`w_k=w_0`$; at this slowest site, there is an aggregate with mass of order volume.

## VII Summary

In this paper, we introduced a minimal model to study the effect of quenched, sitewise disorder in a aggregation-fragmentation system. Our model had some simplifications: the fragmentation was allowed to occur only at the trapping sites, and mass-independent kernels for aggregation and fragmentation were considered. Despite these simplifications, it retains the important physical effects of diffusion, aggregation, fragmentation and trapping.

We studied the case of a single-chipper aggregation model in detail. In all cases except the $`1d`$, biased case, the system exists in the UA phase in which an infinite aggregate localised at the chipper site can coexist with mobile infinite aggregates in the bulk. The simultaneous existence of more than one infinite aggregate is a new feature absent in previous studies of translationally invariant systems. In the $`1d`$, biased case, there is a phase transition from a phase in which a localised aggregate is formed at the chipper site (PA phase) to the one in which no aggregate is formed anywhere in the system (NA phase) as the density is decreased.

We also studied a variant of the above aggregation model in which particles chip off at a single site but diffuse freely in the bulk. This model can be solved exactly and shows a phase transition from the PA phase to the NA phase as density is increased in all dimensions for all bias.

Finally we discussed a likely scenario for the aggregation model in the presence of extensive disorder and argued that interaction effects arising due to coalescence can be ignored on large enough time scales. It would be interesting to check this expectation by a more detailed study of the system with extensive disorder.

Acknowledgements: We thank D. Dhar and Rajesh R. for useful discussions and comments on an earlier version of the manuscript.


format: text

text:
# Strange Cepheids and RR Lyrae

## 1. Introduction

Strange modes were discovered and named by Saio, Wheeler & Cox (198) in highly nonadiabatic stars. (For a general reference cf. Gautschy & Saio 1996). In a recent paper (Buchler, Yecko & Kolláth 1997; hereafter Paper I) it was demonstrated that strange, or surface modes not only exist also in classical Cepheids, thrive there. Paper I showed that the strange modes are vibrational modes that are predominantly surface modes that exist quite naturally even in the adiabatic limit.

This was exhibited with the help of simple transformations of the radial coordinate to an acoustic depth $`z`$, defined through $`dz=dr/c_s`$ and of the radial displacement vector $`\psi =\sqrt{\rho c_s}r\delta r`$, where $`c_s`$ is the sound speed, with which the adiabatic radial pulsation problem for radiative models can be reduced to a time-independent Schrödinger equation

$$\frac{d^2\psi }{dz^2}+V(z)\psi =\omega ^2\psi $$
(1)
The potential

$$V(z)=\frac{1}{\sqrt{\rho c_s/r^2}}\frac{d^2}{dz^2}\sqrt{\rho c_s/r^2}\frac{4Gm}{r^3}$$
(2)
incorporates thus the spatial variations of density $`\rho `$, sound speed $`c_s`$ and spherical effects, as well as gravity. This problem is very similar to that encountered in the study of wind instruments, such as the trumpet, where Eq. 1 is known as the Bernoulli-Webster equation.

In Cepheids and RR Lyrae the hydrogen partial ionization, in which $`c_s`$ varies extremely rapidly, produces a very sharp and enormously high potential barrier, that splits the pulsating envelope into two almost disjoint parts (cf. Paper I). Using the scaled eigenvectors $`\mathrm{\Psi }`$ one readily sees that the ’normal’ vibrational modes have a large amplitude in the interior region. It is however also possible for some modes to be predominantly trapped in the outer region – these surface modes are the ’strange’ modes.

When, as is customary, the modes are classified with increasing number of nodes (and thus increasing frequency) the strange modes appear when a quarter wavelength fits approximately into the surface region. For a typical 5$`M_{}`$ Cepheid on the left of the instability strip one finds that the strange modes have $`n`$ = 10 – 12, and somewhat smaller $`n`$ on the cold side of the instability strip because there is more matter outside the partial ionization region of H.

When nonadiabatic effects, which are weak in Cepheids, are properly taken into account, the modal frequencies acquire of course an imaginary component. Most of the stellar envelope is damping, but the partial ionization zones cause strong pulsational driving. Since the strange mode has at best a very small amplitude in the He partial ionization region it thus experiences no driving or damping there. All its driving must occur in the H partial ionization region. It is not possible to extend the simplified description of driving found in textbooks (e.g., Cox & Giuli 1965) to modes higher than the fundamental or the first overtone. For the higher modes the driving or damping are a delicate function of the relative phase of the pressure and specific volume eigenfunctions which exhibit an increasing number of oscillations (related to the nodes in the adiabatic limit) (e.g., Glasner & Buchler 1993, Buchler, Yecko, Kolláth 1997) as the work-integral shows

$`W_k`$ $`=`$ $`{\displaystyle \frac{2\pi }{I_k}}\mathrm{Im}{\displaystyle \delta p_k^{}\delta v_k𝑑m}`$ (3)
$`=`$ $`{\displaystyle \frac{2\pi }{I_k}}{\displaystyle }|\delta p_k||\delta v_k|\mathrm{sin}(\varphi _v\varphi _p)_kdm`$ (4)
$`I_k`$ $`=`$ $`\omega _k^2{\displaystyle |\delta r_k|^2𝑑m}`$ (5)
Fig. 1. – Cepheid model, M=8, L=10000, $`T_{\mathrm{ef}\mathrm{f}}`$=6700 K model; Top: work-integrand w vs. Log10(pressure) \[cgs\] for modes 0 (fundamental), 1 (first overtone) and 9 (strange mode); middle: radial displacement eigenvectors for modes 0, 1, 8 and 9; bottom: transfored radial displacement eigenvector ($`\mathrm{\Psi }`$) as a function of acoustic depth z (for clarity mode 1 has been scaled up by a factor of 2 and modes 8 –10 by 12).

In Fig. 1 (middle) we exhibit the linear displacement eigenvectors (modulus) as a function of Log p of the fundamental, first overtone and strange 9th overtone for a Cepheid model with M=8$`M_{}`$, L=10 000$`L_{}`$, $`T_{\mathrm{ef}\mathrm{f}}`$=6700 K, X=0.70, Z=0.02. Clearly that the strange mode is very strongly restricted to the surface region.

This even more apparent in the bottom of Fig. 1 which displays the scaled displacement eigenvectors $`\psi `$ (cf. Eq. 1) as a function of acoustical depth $`z`$ for the fundamental, first overtone, and for overtones 8 through 10. As pointed out in Paper I, this representation of the modes as $`\psi (z)`$ brings out the vibrating string nature of the modes, and shows that the strange mode (9th overtone) is mainly a surface mode in contrast to its neighbors.

In the top Fig. 1 which displays the work-integrand as a function of Log p, the integrand is normalized so that the area under the curve represents the relative growth-rate $`\eta =2\kappa _kP_k`$. The solid line represents the strange mode (9th overtone). For comparison we also show the work-integrands (multiplied by 5 for improved visibility) for the fundamental (dotted) and first overtone (dashed). The latter modes are very stable for this model which is on the left of the instability strip (cf. Fig. 1). The strange mode is restricted to the surface regions and its driving comes entirely from the H partial ionization region around Log p $`3`$.

Some simple a posteriori explanation of the odd stability properties of the strange modes can be given. The strange modes penetrate much less deeply than the normal modes. For the strange modes $`W_k`$ is thus (a) totally dominated by the outer regions, in particular the hydrogen partial ionization region, and (b) $`W_k`$ is much smaller. However, its the moment of inertia $`I_k`$ is also very small because the mode is localized in the surface. Since the smaller work-integrals $`W_k`$ are normalized by the smaller $`I_k`$, the strange modes still have relative growth-rates that are comparable in magnitude to those of the normal modes. However, because only hydrogen driving contributes, the growth-rates of the strange modes are egregious as was already shown in Paper I. In that paper strange modes were found to be vibrationally unstable to the left (blue) of the usual fundamental and first overtone instability strips.

Furthermore nonlinear hydrodynamical modelling indicated that the strange modes can give limit cycles, albeit with rather small pulsation amplitudes, typically in milli-magnitudes. This has opened the exciting possibility of observing ’strange Cepheids’, and of using such additional information to put constraints on theoretical Cepheid modelling.

Two questions were left unanswered in Paper I, however. The first concerned the surface boundary condition, namely the disregard of the surface impedance, due to the reaction of the isothermal atmosphere and of running wave or acoustic damping. The second resulted from the neglect of turbulence and convection which can also affect driving and damping. Here, in §2 we reexamine and improve on the outer boundary condition. In §3 and §4 we examine turbulent convective Cepheid and RR Lyrae models. In §5 we present nonlinear pulsations of turbulent convective models and we conclude in §6.

## 2. Isothermal Atmosphere and Boundary Conditions

In paper I the models were integrated up to a point where the gas pressure was small, typically a few times the radiation pressure. In the linearization of the discretized equations we imposed a standard condition of a zero Lagrangean pressure variation, $`\delta p_{}=0`$, i.e., to an open and perfectly reflecting boundary. This boundary condition disregards two effects, first, running, escaping waves, and second, inertial effects due to the isothermal atmosphere. In acoustic parlance, these effects give rise to resistive and reactive surface loads (Morse & Ingard 1968).

It is interesting to note that the change of variables from $`\delta r`$ to $`\mathrm{\Psi }`$, because of the density dependence, effectively changes the problem from an ’open quarter-wave length tube’ to a ’closed half-wave length tube’. (Here we have taken the stellar surface to be defined as the point where the gas pressure equals two times the radiation pressure, i.e., at a very low density.)

Another way of looking at this is to say that when $`\omega `$ is sufficiently large, the concomitant wave length gets sufficiently short compared to the the isothermal scale height $`1/\alpha `$ so that reflection at the surface becomes incomplete and acoustic energy is lost into the atmosphere.

Fig. 2. – Cepheid model, M=8$`M_{}`$, L=10000$`L_{}`$. Left: growth-rates as a function of $`T_{\mathrm{ef}\mathrm{f}}`$, positive growth-rates correspond to linear unstable modes. The modes are labelled with the modal number. Right: period ratios $`P_k/P_0`$ as a function of $`T_{\mathrm{ef}\mathrm{f}}`$. The fundamental period $`P_0`$ \[d\] is indicated on top. The triangles denote the unstable overtone modes. The solid line at the bottom is the isothermal cut-off period ratio.

In the optically thin outer region of the Cepheid model is essentially isothermal $`T(1/2)^{1/4}`$ $`T_{\mathrm{ef}\mathrm{f}}`$ . To a good approximation the density decreases exponentially with scale height $`\alpha ^1=T/(\mu g)`$, where $`g=GM/r^2`$ is essentially constant. This gives rise to a potential barrier which rapidly flattens to a constant value of $`V_{}=(c_s\alpha /2)^2`$. In Buchler, Whiting & Kollath (2001) it is shown that the presence of the isothermal atmosphere can be approximately taken into account through a boundary condition that is imposed at the base of this atmosphere for the modes below the barrier ($`\omega <\omega _c`$). The modes above the barrier disappear into the continuum and are no longer astrophysically relevant.

$$𝒮(\mathrm{\Psi }^{}/\mathrm{\Psi })_o=ik\frac{(\omega +k)+(\omega k)\mathrm{exp}(2ik\mathrm{\Delta })}{(\omega +k)(\omega k)\mathrm{exp}(2ik\mathrm{\Delta })}$$
(6)
where $`kk_r+ik_i=\sqrt{\omega ^2V_{}}`$. When $`|k_i\mathrm{\Delta }|`$ is large then $`𝒮ik\mathrm{sgn}k_i`$. The approximation that is inherent in the expression for the boundary condition is that the potential is a square barrier of width $`\mathrm{\Delta }`$. This translates into a boundary condition

$$\delta p_{}=\mathrm{\Gamma }_1p_{}\left(1+(\alpha /2+𝒮)R_{}\right)\delta r/R_{}$$
(7)
With this new boundary condition our eigenvalue problem becomes more complicated and delicate. We compute the eigenvalues with a modified Castor method in which the $`\delta p_{}`$ of Eq. 7 replaces the usual $`\delta p_{}=0`$ boundary condition. The presence of strange modes and their frequent proximity to the regular modes makes the Castor search difficult. We therefore first compute all the nonadiabatic eigenvalues for the $`\delta p_{}=0`$ case with a general eigenvalue solver (e.g., Glasner & Buchler 1993) and then switch on the new boundary condition gradually.

Fig. 3. – Cepheid model, M=6$`M_{}`$, L=3000$`L_{}`$; cf. Fig. 2 for legend.

## 3. Turbulent Convective Cepheid Models

Turbulent convection changes the structure of the outer envelope, as well as the pulsational driving. To study its effects on the strange modes and to ascertain that they survive convection, we have used the Florida code (Yecko, Kolláth & Buchler 1998) with the $`\alpha `$’s defined as in Feuchtinger, Buchler & Kollath (2000) with values that are very similar to those used in our study of first overtone Cepheids, viz. $`\alpha _d`$= 2.177, $`\alpha _c`$= 0.4, $`\alpha _s`$= 0.433, $`\alpha _\nu `$= 0.12, $`\alpha _t`$= 0.001, $`\alpha _r`$= 0.4, $`\alpha _p`$= 0 and $`\alpha _l`$= 1.5. We use the OPAL opacities of Iglesias & Rogers (1996) and of Alexander & Ferguson (1994).

In Figs. 2, 3 and 4 we show the results for 3 sequences of Cepheid models of fixed mass and luminosity over a broad range of $`T_{\mathrm{ef}\mathrm{f}}`$. The left subfigures show the growth-rates as a function of the $`T_{\mathrm{ef}\mathrm{f}}`$ of the equilibrium models. The growth-rates for the unstable modes are positive. Some modes are very stable and lie below the graph.

As a function of $`T_{\mathrm{ef}\mathrm{f}}`$ the modal number (ordering with period) of the strange mode changes since it switches its properties with the regular modes through avoided crossings (cf. Paper I). It is therefore ambiguous how we should label the modes along a sequence, but this labelling is not important for us since we are interested in finding out whether there is an unstable strange mode and what its period is, both of which are unambiguous.

The instability strip of the 8$`M_{}`$ in Fig. 2 extends from the fundamental blue edge at $`T_{\mathrm{ef}\mathrm{f}}`$= 5700 K to the fundamental red edge at $`T_{\mathrm{ef}\mathrm{f}}`$= 4600 K. An unstable strange mode exists from 7700 K through 6100 K which is the 11th overtone at the lower $`T_{\mathrm{ef}\mathrm{f}}`$, but which gradually morphs into the 8th overtone by the time it gets to 8000 K. A correspondingly erratic behavior shows up in the period ratios. It is not a numerical artifact, but arises from avoided level crossings (cf. Paper I) that are also familiar from the qunatum mechanical behavior of an atom in a B field and from the nonradial $`p+g`$ spectrum in unevolved stars.

Fig. 4. – Cepheid model, M=4$`M_{}`$, L=860$`L_{}`$; cf. Fig. 2 for legend.

The instability strip of the 6$`M_{}`$ model in Fig. 3 extends from the first overtone blue edge at $`T_{\mathrm{ef}\mathrm{f}}`$= 6200 K to the fundamental red edge at $`T_{\mathrm{ef}\mathrm{f}}`$= 5150 K. From 7000 K to 6450 K the 10th overtone is unstable. Around 5000 K the 8th overtone just misses being unstable.

The instability strip of the 4$`M_{}`$ model Fig. 4 extends from the first overtone blue edge at $`T_{\mathrm{ef}\mathrm{f}}`$= 6550 K to the fundamental red edge at $`T_{\mathrm{ef}\mathrm{f}}`$= 5500 K, and it includes a small region with second overtone instability. From 7700 K to 6950 K the 12th overtone is unstable, and from 5200 K to 5100 K it is the 9th overtone.

## 4. Turbulent Convective RR Lyrae Models

We consider an RR Lyrae model with M=0.65$`M_{}`$, L=45$`L_{}`$ with X=0.70, Z=0.02. Fig. 5 shows that the instability strip of the model extends from the first overtone blue edge at $`T_{\mathrm{ef}\mathrm{f}}`$= 7250 K to the fundamental red edge at $`T_{\mathrm{ef}\mathrm{f}}`$= 6000 K. From 8300 K to 7200 K the 10th overtone is unstable, and from 5900 K to 5800 K the 8th overtone is unstable.

Again, these results show that an excited strange mode can exist on either side of the classical instability strip.

## 5. Nonlinear Behavior

The full amplitude pulsations have been computed with the Florida hydrodynamics code (Kolláth, Beaulieu, Buchler & Yecko, 1998) with the same mesh and parameters as used for the linear properties. The pseudo-viscosity was taken to be zero in our calculations.

Fig. 5. – RR Lyrae model, M=0.65$`M_{}`$, L=45$`L_{}`$; cf. Fig. 2 for legend.

In Fig. 6 we exhibit the light and radial velocity curves for three selected Cepheid models and on RR Lyrae model. The lightcurves are measured in millimags, and the radial velocities in units of 100 m/s. Note that five of the models (1, 3 – 6) are blue strange variables, i.e., they lie to left of the blue edge, and two (2 and 7) are red strange variables, i.e., they lie to right of the blue edge. Interestingly they all have about the same pulsation amplitudes. Because the oscillation amplitude is so small and it is so confined to the surface region, there is not much matter below the photosphere that can modulate the heat flow. The light-curve amplitudes are therefore also small, in the millimagnitude regime, but they are observable.

The temporal structure of the photospheric velocities appears strange at first. In Fig. 7 we therefore show the behavior of the velocities throughout most of the model, in fact Lagrangean zones 65 through 145 of a 150 mass-shell model. The interior zones are essentially quiescent because of the high modal number of this mode. On the other hand the amplitude of the very outer zones increases considerably (but is not relevant because it contains so little mass and is considerably above the photosphere. The photosphere itself, which is shown as a single dashed line, is relatively quiescent. For reference we also show the hydrogen partial ionization region with three adjacent dashed lines. The helium ionization lies below in zones that are not shown because the velocities are small. The photosphere is seen to lie in a relatively quiescent region and is subjected to inward and outwards moving ripples that account for its sometimes strange looking behavior in Fig. 6.

Fig. 6. – Left: Light curves (in millimag): Right: radial velocity curves (in 100m/s). From top to bottom, Cepheids: M=4$`M_{}`$\[A\], $`T_{\mathrm{ef}\mathrm{f}}`$=7300 and 5150 K, M=6$`M_{}`$\[C\], $`T_{\mathrm{ef}\mathrm{f}}`$=6700 K, M=8 \[E\], $`T_{\mathrm{ef}\mathrm{f}}`$=7100 and 6700 K, RR Lyrae: M=0.65$`M_{}`$\[R\], $`T_{\mathrm{ef}\mathrm{f}}`$=7800 and 5900 K.

There are two shortcomings to these calculations. First, the (linear) special boundary condition at the outer surface (Eq. 7) has not been incorporated in the nonlinear code, but we do not expect that to cause a major distortion, for the pressure is so small at the boundary. Second, what is perhaps more severe, there is no provision for acoustic leakage in the nonlinear hydrodynamics code for waves with frequencies above the isothermal threshold $`\omega _c`$.

The presence of higher frequency reflected waves is perhaps the reason for the pronounced harmonic wiggles on the light curves and especially in the radial velocity curves. Indeed, because spherical effects are negligible, the higher order strange modes are all in resonance. However, if they are excited along with the lowest strange mode, then the even harmonics should be absent, but they are not.

## 6. Discussion and Conclusion

The improved boundary conditions have a small effect on the behavior of the deep modes, mostly because the pressure variation at the stellar surface $`\delta p_{}`$ is always small. However, they do affect the higher overtones as they approach the cutoff frequency by shifting their periods, sometimes causing inversions in the modal ordering and affecting the stability of the modes.

Turbulent convection does not kill the linear instability of the strange modes, as one might have feared. On the contrary, they are as unstable as with the radiative models. One even finds unstable strange modes to the right of the fundamental red edge.

The growth-rates of the strange modes on the hot side of the instability strip are not very sensitive to turbulent convection. On the other hand convection is very strong in the cooler models, and the values of the stability coefficients of the strange modes on the cool side are consequently less certain than they are on the hot side.

Fig. 7. – Radial velocity profile for the Cepheid model M=6$`M_{}`$, $`T_{\mathrm{ef}\mathrm{f}}`$=6700 K; the individual Lagrangean zones (zones 65 through 145 for a 150 zone model) are shifted vertically for better visualization; the dashed line corresponds to the photospheric velocity. the three adjacent dashed lines denote the hydrogen partial ionization region

So, do strange Cepheids really exist in nature? We find the occurrence of self-excited strange modes to be quite robust, i.e., largely independent of any modelling parameters. For an individual model the change in modelling parameters may have a large effect, but when one considers a whole sequence of models on finds that modelling parameters (e.g., the fineness of the grid) merely cause a small shift (e.g., in $`T_{\mathrm{ef}\mathrm{f}}`$). Especially in the high luminosity Cepheid models, the strange mode has a period well above the isothermal cutoff period as the right hand side of Fig. 2 shows. For the low mass Cepheids and the RR Lyrae some of the unstable strange modes occur very close to the isothermal cutoff frequency $`\omega _c`$. This cutoff frequency however is somewhat uncertain because it depends on the structure of the dilute atmosphere which itself is sensitive to the treatment of radiative transport (treated here in the customary diffusion approximation).

We find furthermore that the full amplitude (nonlinear) pulsations are quite robust in all cases. However, the values of the pulsation amplitudes are somewhat uncertain because they depend on modelling parameters as they did for the classical Cepheid pulsations. In particular the turbulent eddy viscosity has a large effect. (We recall that our convective parameters $`\alpha `$’s were calibrated to give good agreement with the observations of classical Cepheids, Feuchtinger et al. 2000).

The detection of this type of pulsation would be very important because it would add very strong additional constraints on the Cepheid and RR Lyrae modelling. The search for millimagnitude variability on either side of the classical instability strip with periods approximately 1/4 to 1/5th the fundamental period (cf. Figs. 1 –3, 5) would be required. Hopefully this will become possible with COROT.

## 7. Acknowledgements

It is a pleasure to thank Phil Yecko and Bernard Whiting for fruitful discussions. This work has been supported by NSF (AST 98-19608).