text:
# Crossover effects in a discrete deposition model with Kardar-Parisi-Zhang scaling

## I Introduction

Surface growth processes and deposition of thin films are of great interest due to potential technological applications (such as production of nanostructures for microelectronic devices, the possibility of growth of designed quantum objects, magnetic storage devices, among others) and due to the fundamental role these systems play in non-equilibrium statistical physics . Several models have been investigated in the last decade, most of them involving one kind of particle and a simple microscopic aggregation rule. The competition between different growth mechanisms have received less attention, but is essential to describe some practical situations, such as growth of materials designed to have specific electronic, mechanical or magnetic properties, which involves deposition of two or more chemical species. In this framework, some authors considered growth models with two kinds of particles and different aggregation rules. Other situations involving competition between two growth mechanisms have also been considered .

These models usually show crossover effects from one dynamics at small times $`t`$ or short length scales $`L`$ to another dynamics at long $`t`$ and large $`L`$. One typical example is Kardar-Parisi-Zhang ($`KPZ`$) growth at small nonlinearities . The Langevin-type equation

$$\frac{h}{t}=\nu ^2h+\frac{\lambda }{2}\left(h\right)^2+\eta (\stackrel{}{x},t),$$
(1)
known as $`KPZ`$ equation, was proposed as a hydrodynamic description of kinetic surface roughening. Here $`h`$ is the height at the position $`\stackrel{}{x}`$ in a $`d`$-dimensional substrate at time $`t`$, $`\nu `$ represents a surface tension, $`\lambda `$ represents the excess velocity and $`\eta `$ is a Gaussian noise  with zero mean and variance $`\eta (\stackrel{}{x},t)\eta (\stackrel{}{x^{}},t^{})=D\delta ^d\left(\stackrel{}{x}\stackrel{}{x^{}}\right)\delta \left(tt^{}\right)`$. When the coefficient $`\lambda `$ of the nonlinear term is small, a crossover is observed from linear growth ($`\lambda =0`$, known as Edwards-Wilkinson theory - $`EW`$ to $`KPZ`$ behavior.

In discrete models, the interface width, which characterizes the roughness of the interface, is defined as

$$W(L,t)=\left[\frac{1}{L^d}\underset{i}{}\left(h_i\overline{h}\right)^2\right]^{1/2}$$
(2)
for deposition in a $`d`$-dimensional substrate of length $`L`$ ($`h_i`$ is the height of column $`i`$ at time $`t`$, the bar in $`\overline{h}`$ denotes a spatial average and the angular brackets denote a configurational average). For short times it scales as $`Wt^\beta `$ and for long times, in the steady state regime, it saturates at $`W_{sat}L^\alpha `$. The dynamical exponent $`z=\alpha /\beta `$ characterizes the crossover from the growth regime to the steady state regime. For systems belonging to the $`EW`$ universality class, we have $`\alpha _0=1/2`$, $`\beta _0=1/4`$ and $`z_0=2`$ in $`d=1`$ (in this paper, the subscript $`0`$ will refer to exponents of the $`EW`$ theory). For systems in the $`KPZ`$ class, in $`d=1`$, we have $`\alpha =1/2`$, $`\beta =1/3`$, $`z=3/2`$ .

Considering the crossover from $`EW`$ to $`KPZ`$ scaling in $`d=1`$, Grossmmann, Guo and Grant ($`GGG`$)  and Nattermann and Tang ($`NT`$)  (see also the review by Forrest and Toral ) proposed multiscaling relations that are equivalent to

$$W(L,t)=L^\alpha f(\frac{t}{t_c},\frac{L}{\xi _c}),$$
(3)
in which $`\xi _ct_c^{1/z_0}`$. $`GGG`$ also proposed that the characteristic time of crossover from $`EW`$ to $`KPZ`$ dynamics was

$$t_c\lambda ^\varphi ,$$
(4)
with $`\varphi >0`$, since the $`EW`$-$`KPZ`$ crossover disappears for $`\lambda =0`$. Through scaling arguments, those authors obtained $`\varphi =z_0/(\alpha _0+z_02)`$, which gives $`\varphi =4`$ in $`d=1`$. This was confirmed through one-loop renormalization group calculations by $`NT`$. The scaling analysis of the $`KPZ`$ equation by Amar and Family ($`AF`$) and the assumption of Family-Vicsek scaling were used to show that, in the nonlinear and saturation regimes,

$$W(L,t)L^{1/2}g(|\lambda |\frac{t}{L^{3/2}}),$$
(5)
in which $`g`$ is a scaling function and the dependence of $`W`$ on the parameters $`\nu `$ and $`D`$ of Eq. (1) was omitted. A generalized scaling relation equivalent to Eq. (3), which is a more general result than Eq. (5), was also obtained by Derrida and Mallick in the context of the connection to the one-dimensional asymmetric exclusion model . Amar and Family  have shown that the scaling form (5) also predicts a crossover exponent $`\varphi =4`$.

On the other hand, all previous numerical results suggested $`\varphi 3`$; for instance, $`GGG`$ obtained this value using data collapse methods . Thus, it would be desirable to confirm numerically the scaling properties predicted for a $`KPZ`$ system in order to solve this controversy.

The purpose of this work is to study a competitive growth process with $`EW`$ to $`KPZ`$ crossover, involving ballistic deposition ($`BD`$ and random deposition with surface relaxation (Family model)  in $`d=1`$. In this model, incident particles aggregate to the deposit according to the rules of $`BD`$ with probability $`p`$ and according to the rules of the Family model with probability $`1p`$. It is known that the Family model is in the $`EW`$ universality class, while $`BD`$ is in the $`KPZ`$ class. This competitive model was introduced by Pellegrini and Jullien , whose main interest was the connection to the roughening transition present in higher dimensions. Although it is expected that this model is in the $`KPZ`$ class for any $`p>0`$, the crossover in $`d=1`$ was not studied in detail in their original work and, for $`p0.3`$, effective exponents very near the $`EW`$ values were obtained .

Here we will simulate that model in order to analyze the interface width scaling in the nonlinear regime, the crossover to the saturation regime and to connect the parameter $`p`$ and the coefficient $`\lambda `$ of the $`KPZ`$ equation in the corresponding continuum limit. The amplitudes of typical saturation times and of interface width scaling in the growth regime are consistent with multiscaling concepts  and refine previous numerical estimates for related systems. The crossover exponent $`\varphi =4`$ follows directly from our numerical results and, together with the observed relation $`\lambda p^{2.1}`$, indicate that the crossover at small $`p`$ is very slow. The analysis of this apparently simple problem shows that, in order to obtain reliable asymptotic exponents governing various quantities, it is essential to account for corrections to the leading terms in the scaling relations. Thus, this work may also be relevant to the analysis of other systems with slow crossover to $`KPZ`$ scaling, whose interest increased after the recent debate on the problem of Fisher waves and their discrete realizations in $`d=1`$ dimensions . For that reason, the crossover effects identified in our simulations’ data will be discussed in detail.

It is also relevant to point out that a related competitive model was recently studied in $`d=1`$ and $`d=2`$ , showing evidence of the asymptotic $`KPZ`$ behavior. However, that work did not study the relation between the parameters of the discrete and the continuous ($`KPZ`$) model nor the scaling amplitudes that will be considered here.

The rest of this work is organized as follows. In Sec. II, we will define precisely the discrete model and connect it to the $`KPZ`$ equation using the scaling properties of the growth velocity. In Secs. III and IV we will present results for the interface width scaling in the discrete model at the nonlinear growth regime and at the steady state regime, respectively. In Sec. V we summarize our results and present our conclusions.

## II The discrete model and its connection to the $`KPZ`$ theory

We considered a model in which particles are aggregated following the rules of $`BD`$ with probability $`p`$ or the rules of random deposition with surface relaxation (Family model) with probability $`1p`$. In $`BD`$ (Fig. 1a), the incident particle follows a straight trajectory perpendicular to the surface and sticks upon first contact with a nearest neighbor occupied site. It leads to the formation of a porous deposit. In the Family model (Fig. 1b), the particle falls towards the surface along the incident column and sticks at the top of that column if its height is lower than or equal to the heights of the neighboring columns. Otherwise, the particle diffuses to the neighboring column which has the lowest height and, if two or more neighbors have the same height, it chooses one of them randomly.

For $`p=0`$, we have the Family model, which is in the $`EW`$ universality class. For any $`p0`$, in $`d=1`$, we expect the $`BD`$ process to change the universality class to $`KPZ`$ in the continuum limit (see the analysis in Ref.  for a related model). Then the coefficient $`\lambda `$ of the nonlinear term vanishes with $`p`$ in the form

$$\lambda p^\gamma ,$$
(6)
with $`\gamma >0`$ (to be estimated below). For small $`p`$ and sufficiently large $`L`$, the interface width $`W(L,t,p)`$ must scale analogously to the weak coupling regime of the KPZ theory , in which three regimes were identified: a linear ($`EW`$) growth regime at early times ($`tt_c`$), a nonlinear ($`KPZ`$) growth regime for $`t_ct\tau `$ and the saturation regime for $`t\tau `$, as illustrated in Fig. 2 ($`\tau `$ is the characteristic time for the interface width saturation).

In order to calculate the exponent $`\gamma `$, we considered the scaling of the interface growth velocity. The difference between the growth velocity in an infinitely large substrate, $`v_{\mathrm{}}`$, and the velocity in the steady state of a finite lattice (thick films), $`v(L)`$, scales as

$$\mathrm{\Delta }v(L)v_{\mathrm{}}v(L)\lambda L^\alpha _{},$$
(7)
with $`\alpha _{}=1`$ in $`d=1`$ . Defining

$$b_v(L)\mathrm{\Delta }v(L)\times L,$$
(8)
we expect that, as $`L\mathrm{}`$,

$$b_v(L)B_v=B\lambda ,$$
(9)
where $`B`$ is a constant.

In the discrete model, $`\lambda `$ varies with $`p`$, consequently $`b_v`$ is a function of $`L`$ and $`p`$ which has a finite limiting value $`B_v(p)`$ as $`L\mathrm{}`$. For very large $`L`$, Eqs. (6) and (9) shows that $`B_v`$ scales with $`p`$ with exponent $`\gamma `$.

Simulations of the model were performed in lattice of lengths from $`L=16`$ to $`L=4096`$ until the saturation regime, and in lattices with $`L=2^{16}=65536`$ during the growth regime (linear and nonlinear), for several values of the probability $`p`$ between $`p=0.15`$ and $`p=0.5`$. The results presented in this paper are averages typically over $`10^5`$ realizations for the smallest lattices ($`L256`$), $`10^4`$ realizations for $`256L4096`$ and $`10^2`$ realizations for $`L=65536`$. The growth velocities were calculated from numerical derivatives of the average heights of the deposits, with accuracies from $`5`$ to $`6`$ decimal places, in lattices of lengths $`L128`$ ($`L512`$ for $`p=0.15`$). We considered the data for $`L=65536`$ as representative of an infinite lattice in the growth regime (some simulations in $`L=131072`$ supported this assumption), and also obtained $`v_{\mathrm{}}`$ with high accuracy. These data provided estimates of $`b_v(p,L)`$ with accuracy from $`0.5\%`$ to $`5\%`$. For larger lengths, poorer results were obtained due to the much smaller number of realizations.

In Fig. 3 we show $`b_v(p,L)`$ versus $`1/L`$ for the three smallest values of $`p`$ considered in this work. The variable $`1/L`$ in the abscissa was the best choice to represent finite-size corrections in $`b_v`$ as $`L\mathrm{}`$, and is related to higher order terms ($`1/L^2`$) in Eq. (7). Such scaling corrections have been previously observed in the analysis of small $`L`$ data for $`BD`$ and for the restricted solid-on-solid model by Krug and Meakin . The corrections were considered in the extrapolation of the data in Fig. 3, which provided estimates of $`B_v(p)`$ for several values of $`p`$ (intercepts with the vertical axis in Fig. 3).

Crossover effects may be crucial in the extrapolation procedure discussed above, and may severely affect the estimates of $`B_v`$ for small $`p`$. For $`0.25p0.5`$, four values of $`b_v`$ ($`16L128`$) were well fitted by straight lines in the $`b_v(p,L)\times 1/L`$ plots (these results were not shown in Fig. 3, except for $`p=0.25`$). For $`p=0.2`$, the data for $`32L256`$ confirm the presence of the $`1/L`$ correction and was also used to estimate $`B_v`$ (the estimate for $`L=16`$ deviates from this trend). On the other hand, for $`p=0.15`$, the result for $`L=256`$ showed a crossover in $`b_v`$, which suggested calculations for $`L=512`$. Fig. 3 shows that $`b_v(0.15,L)`$ slowly increases for $`16L128`$, but decreases for $`128L512`$. Consequently, the extrapolation considered only the three last points (see Fig. 3) and gave $`B_v0.057`$. However, if the extrapolation to $`L\mathrm{}`$ was performed only with results for $`L128`$, then a $`7\%`$ larger value of $`B_v`$ would be obtained. Smaller values of $`p`$ were not studied here because such crossover would appear for much larger $`L`$ and, consequently, the extrapolations based on small systems’ data would provide unreliable estimates of $`b_v(p,\mathrm{})`$.

In Fig. 4 we show $`\mathrm{ln}\left[B_v(p)\right]`$ versus $`\mathrm{ln}p`$ using the extrapolated values of $`B_v`$, as discussed above. The linear fit in Fig. 4 gives $`B_vp^{2.1}`$. Considering the error bars in $`B_v`$, we obtain an exponent $`\gamma =2.1\pm 0.2`$ (Eq. 6).

The large value of the exponent $`\gamma `$ explains the crossover effect discussed above. Since $`\lambda `$ decreases rapidly with $`p`$, the coefficient of the leading term in $`\mathrm{\Delta }v`$ (Eq. 7) is small compared to higher order corrections ($`1/L^2`$, $`1/L^3`$ etc) for small $`p`$. Thus, very large values of $`L`$ are needed to provide reliable extrapolations with a single correction term, which prevented us to use values of $`p<0.15`$ in our study.

## III Interface width at the nonlinear growth regime

In the nonlinear growth regime ($`t_ct\lambda ^1L^z`$) for sufficiently large substrates in $`d=1`$, the interface width does not depend on $`L`$ (weak finite-size effects). Then the scaling function of Eq. (5) behaves as

$$g(x)Cx^\beta ,\beta =1/3,$$
(10)
with constant $`C`$, so that $`W`$ does not depend on $`L`$, except for vanishing corrections to scaling. Consequently, the $`\lambda `$-dependent scaling of $`W`$ in this regime is

$$WC\lambda ^\beta t^\beta .$$
(11)
In this section, we will verify this $`\lambda `$-dependence through a careful analysis of simulations’ data of our discrete model.

However, first we will show that Eq. (11) gives $`\varphi =4`$ (Eq. 4) in a simple way, as follows. The crossover $`EWKPZ`$ (at $`tt_c`$) occurs when the scaling relation (11) matches the $`EW`$ scaling

$$W(t,L)C^{}t^{\beta _0},$$
(12)
with $`C^{}`$ constant. Thus

$$t_{c}^{}{}_{}{}^{\beta _0}\lambda ^\beta t_{c}^{}{}_{}{}^{\beta },$$
(13)
then we obtain $`\varphi =\frac{\beta }{(\beta \beta _0)}=\frac{z_0}{(z_0z)}=4`$.

We conclude that the numerical test of Eq. (11), in particular of the dependence on the parameter $`\lambda `$, may be used to test the proposal $`\varphi =4`$. The first step is to extract the amplitude of $`t^\beta `$ scaling in Eq. (11), which motivates the definition of the amplitude $`a(p,t)`$ as

$$a(p,t)W(L\mathrm{},t)/t^{1/3}.$$
(14)
In Fig. 5a we show $`a(p,t)`$ versus $`1/t^{1/3}`$ for several values of $`p`$, using the data for $`L=2^{16}=65536`$. Different variables in the form $`1/t^x`$ ($`x>0`$) were tested in the abscissa, but the variable $`1/t^{1/3}`$ of Fig. 5a provided the best linear fits for most values of $`p`$. The fact that in Fig. 5a $`a(p,t)`$ is still decreasing for large $`t`$ indicates the presence of a constant (independent of $`t`$) correction to the leading behavior in Eq. (11). It proves again the relevance of accounting for scaling corrections in this problem, although we are not able to justify these corrections on theoretical grounds.

As $`t\mathrm{}`$, $`a(p,t)`$ converges to a finite limiting value

$$A(p)a(p,\mathrm{}).$$
(15)
$`A(p)`$ is the complete amplitude of $`t^\beta `$ scaling of the interface width in the nonlinear growth regime (Eq. 11). Our estimates of $`A(p)`$ were obtained from linear extrapolations of $`a(p,t)\times 1/t^{1/3}`$ plots to $`t\mathrm{}`$ (intercepts with the vertical axis in Fig. 5a).

From Eq. (11), it is expected that the amplitude $`A(p)`$ scales as $`\lambda ^\beta `$. From the connection relation (6), it is expected that

$$A(p)p^\delta ,$$
(16)
with

$$\delta =\gamma \beta .$$
(17)
Then, the test of Eq. (11) reduces to the test of Eq. (17) for the amplitude exponent $`\delta `$.

In order to calculate the exponent $`\delta `$ in the relation (16), our first step was to plot $`\mathrm{log}\left[A(p)\right]`$ versus $`\mathrm{log}p`$, but we noticed that it showed decreasing slopes as $`p`$ decreased. We analyzed the evolution of the slopes of $`\mathrm{log}\left[A(p)\right]\times \mathrm{log}p`$ plots by calculating the following effective exponents for consecutive values $`p=p^{}`$ and $`p=p^{\prime \prime }`$:

$$\delta _p=\frac{\mathrm{ln}\left[a(p^{},\mathrm{})/a(p^{\prime \prime },\mathrm{})\right]}{\mathrm{ln}\left(p^{}/p^{\prime \prime }\right)},p=\sqrt{p^{}p^{\prime \prime }},$$
(18)
so that, as $`p0`$ ($`\lambda 0`$), we expect that $`\delta _p\delta `$.

In Fig. 5b we show $`\delta _p`$ versus $`p^2`$, which gives a reasonable linear fit and indicates that $`\delta =0.7\pm 0.2`$. Again the variable $`p^2`$ in the abscissa is the one that provides the best linear fit of the central estimates of $`\delta _p`$, chosen among other variables in the form $`p^y`$ ($`y>0`$). In Fig. 5b, the effective exponents systematically decrease as $`p`$ decreases, which reflects our previous observation of decreasing slopes in $`\mathrm{log}\left[A(p)\right]\times \mathrm{log}p`$ plots.

Our estimates $`\delta =0.7\pm 0.2`$ and $`\gamma =2.1\pm 0.2`$ (Sec. II) are consistent with relation (17) with $`\beta =1/3`$. Even considering that the error bars are large, it is relevant to notice that the central estimates confirm that relation exactly, which gives additional support to our analysis.

## IV Interface width near and at the steady state regime

Our numerical results in the steady state regime provide additional support for the scaling picture proposed for the problem.

From Eq. (5), we expect that the crossover from the nonlinear to the steady state regime takes place at a characteristic time $`\tau `$ that scales as

$$\tau \lambda ^1L^z,$$
(19)
with $`z=3/2`$ in $`d=1`$. In this section, our main purpose is to test the $`\lambda `$-dependence of this characteristic time in our discrete model.

The saturation time $`\tau `$ is usually estimated using some arbitrary recipe. Here, instead of estimating the saturation time $`\tau `$ (which may be defined from the time dependence of the interface width as it converges to the saturation value), we calculated a characteristic time $`\tau _0`$ which is proportional to $`\tau `$, according to a recently proposed method . That method provided accurate estimates of dynamic exponents for several growth models in $`d=1`$ and $`d=2`$, including the Family and the $`BD`$ models.

First, the saturation width $`W_s`$ is estimated, for fixed $`p`$ and $`L`$. Then we define $`\tau _0`$ through

$$W\left(\tau _0\right)=kW_s,$$
(20)
with fixed $`k`$ ($`k1`$) . Using the Family-Vicsek relation $`W(L,t)=L^\alpha f(tL^z)`$ and considering that $`W_sL^\alpha `$, we conclude that $`\tau _0L^z`$, i.e., $`\tau _0`$ is proportional to the saturation time $`\tau `$. For the particular case of a $`KPZ`$ system, Eq. (5) gives

$$\tau _0\lambda ^1L^{3/2}.$$
(21)
Extending the procedure of previous work , we considered $`k=11/e=0.6321\mathrm{}`$ in Eq. (20) to estimate $`\tau _0`$. This value of $`k`$ gave $`\tau _0\tau `$ for $`BD`$, where $`\tau `$ was estimated from the decay of $`W_sW`$ . In the present model, for fixed $`p`$, we calculated the ratios $`\tau _0/L^{3/2}`$ for several lengths $`L`$ and obtained the asymptotic amplitude

$$D(p)=\frac{\tau _0}{L^{3/2}},L\mathrm{}.$$
(22)
The extrapolation procedure follows the same lines of the calculation of $`A(p)`$ from $`a(p,t)`$ in Sec. III. However, only results for $`p0.2`$ could be obtained using data for lattice sizes $`L4096`$, since the saturation for smaller values of $`p`$ is typically of $`EW`$ type ($`\tau L^2`$) in this range of $`L`$.

From Eqs. (22) and (21), we expect that $`D(p)\lambda ^1`$. Consequently, it must scale with $`p`$ as

$$D(p)p^\gamma .$$
(23)
In Fig. 6 we show $`\mathrm{ln}D(p)`$ versus $`\mathrm{ln}p`$, with a linear fit that gives $`D(p)p^{2.1}`$. This result is consistent with the independent estimate of $`\gamma `$ from Eq. (6) (Sec. II).

We also analyzed the scaling of the saturation width $`W_s`$. For lattice sizes $`L1024`$, we obtained $`W_sL^\alpha `$ with $`\alpha =1/2`$ and weak corrections to scaling. Using the data for $`L=1024`$, we defined

$$\mathrm{\Delta }W_sW_s(p)W_s(0)$$
(24)
as the difference between the saturation width for a given probability $`p`$ and the saturation width for the Family model ($`p=0`$).

In Fig. 7 we show $`\mathrm{ln}\mathrm{\Delta }W_s`$ versus $`\mathrm{ln}p`$. The linear fit suggests $`\mathrm{\Delta }W_sp^{3/2}`$, thus we obtain the complete form for the saturation width as

$$W_s(C_1+C_2p^{3/2})L^\alpha $$
(25)
with $`\alpha =1/2`$, $`C_1`$ and $`C_2`$ constants. The amplitude of $`W_s`$ scaling is $`\left(\frac{D}{24\nu }\right)^{1/2}`$ , i.e., the heights’ fluctuations depend only on the parameters $`\nu `$ and $`D`$ of the $`KPZ`$ equation (1), but not on the nonlinearity parameter $`\lambda `$. Thus we conclude that the dependence on $`p`$ in Eq. (25) is related to the dependence on $`p`$ of the surface tension parameter $`\nu `$: when $`p`$ decreases, the amplitude in Eq. (25) decreases, then the parameter $`\nu `$ increases. Indeed, this term is physically expected to increase in the crossover from $`BD`$ (low $`\nu `$) to the Family model (high $`\nu `$).

## V Summary and conclusions

We studied a competitive growth model in $`1+1`$ dimensions involving two dynamics: ballistic deposition with probability $`p`$ and random deposition with surface relaxation (Family model) with probability $`1p`$. This model is a discrete realization of the continuum $`KPZ`$ equation with an adjustable nonlinear coupling $`\lambda `$ related to $`p`$. At the critical probability $`p_c=0`$, the process belongs to the $`EW`$ universality class, while any finite value of $`p`$ drives the system to $`KPZ`$ class.

We established the connection between the parameters $`p`$ and $`\lambda `$ as $`\lambda p^{2.1}`$ and showed that $`Wp^{0.7}t^{1/3}`$ in the growth regime. This indicates that the discrete model presents a very slow crossover from $`EW`$ to $`KPZ`$ scaling at small values of $`p`$, since the crossover time is $`t_c\lambda ^4p^{8.4}`$. This slow crossover explains the discrepancies in the effective exponents $`\beta `$ measured in that regime in previous works .

We also obtained the saturation time $`\tau p^{2.1}L^{3/2}`$. The condition $`\tau t_c`$ is necessary to observe the crossover to $`KPZ`$ scaling, while the opposite condition leads to $`EW`$ saturation without an intermediate $`KPZ`$ growth of the interface width. A critical system size $`\xi _c`$ separates systems which present $`EW`$ or $`KPZ`$ saturation, and $`\xi _c`$ can be estimated from the condition $`\tau t_c`$, which gives $`\xi _c\lambda ^2p^{4.2}`$. This large exponent proves that simulations in very large system sizes are necessary in order to observe all features of $`KPZ`$ scaling for small $`p`$.

Our results are consistent with the scaling theories for the weak coupling regime of the $`KPZ`$ equation proposed by several authors and refine previous numerical analysis. Then we expect that the methods presented here may be helpful to analyze other growth models with slow crossovers to $`KPZ`$ scaling, in which scaling theories cannot be easily developed.


format: text

text:
# Thermodynamical description of a family of partially relaxed stellar systems

## 1 Introduction

The possibility of providing a thermodynamical description of self-gravitating stellar systems has motivated a number of investigations in galactic dynamics, starting with the pioneering work of Antonov and Lynden-Bell in the 60s. After the realization that violent relaxation is likely to lead to partially relaxed configurations in dynamical equilibrium (Lynden-Bell, 1967), a re-examination of the problem of the isothermal sphere, studied earlier by Bonnor (1956) for a self-gravitating gas, led to the interesting possibility that stellar systems may undergo the process of gravothermal catastrophe (Lynden-Bell and Wood 1968; see also Antonov 1962). However, a rigorous derivation of the onset of the gravothermal catastrophe from a study of the Boltzmann entropy appeared to be available only for the case of ideal systems confined by a spherical reflecting wall. A number of convincing qualitative arguments made it clear that also unbound stellar systems with finite mass, such as those described by the King sequence (King 1966; these spherical models have a finite radius, but do not require an external wall), should probably fall into the same physical framework and indeed the paradigm received a lot of attention, especially in the context of the dynamics of globular clusters (see Spitzer 1987), which are known to possess, at least to some extent, the desired internal collisionality (see also Lynden-Bell and Eggleton 1980). \[An indirect indication that the general physical picture of the gravothermal catastrophe is likely to be robust comes also from the proof that the behaviour of the classical gas case is basically independent of the assumption of spherical geometry (Lombardi and Bertin, 2001).\] Several investigations have aimed at producing a rigorous derivation of the gravothermal catastrophe for unbound stellar systems, focusing on the underlying argument that refers to the Poincaré stability of linear series of equilibria (Katz, 1978, 1979; Padmanabhan, 1989), but the proof has always been centered on an unjustified ansatz in order to connect the underlying entropy $`S`$ with the global temperature $`T=1/(S/E_{tot})`$ (see Appendix V in the article by Lynden-Bell and Wood 1968; Katz 1980; Magliocchetti et al. 1998). Other investigations have explored the possibility of setting the discussion in the context of non-standard entropies (e.g., the Tsallis entropies; see the study of the polytropic spheres by Chavanis 2002 and Taruya and Sakagami 2002). Note that the concept of entropy for collisionless systems is quite subtle (e.g., see Stiavelli and Bertin 1987 and references therein). One might even argue whether it is actually compatible with the long-range nature of gravity, given the fact that self-gravitating systems lack additivity, a key ingredient in thermodynamics.

In the meantime, inspired by N-body simulations of collisionless collapse (van Albada, 1982), which confirmed the general picture of incomplete violent relaxation and showed that it can lead to systems with realistic density profiles without ad hoc tuning of the initial conditions, some families of models were constructed able to reproduce, for quasi-spherical configurations, the characteristic feature of the anisotropy profile with an inner isotropic core and an outer radially biased envelope (Bertin and Stiavelli 1984; see Bertin and Stiavelli 1993 and references therein): these families turned out to exhibit the characteristic $`R^{1/4}`$ projected density profile and indeed were shown to match nicely the observed photometric and kinematical characteristics of bright ellipticals. In an attempt at providing a justification of these models (in particular, of the so-called $`f_{\mathrm{}}`$ models, constructed initially only from dynamical arguments) from statistical mechanics, two routes were pursued (Stiavelli and Bertin, 1987).

The first combines an explicit statement of partial relaxation, i.e. of a relaxation process that is expected to be inefficient in the outer regions, and the existence of a suitable weight, related to the orbital period, for the cells that make the relevant partition of phase space; it follows qualitative arguments proposed by Lynden-Bell (1967) and is physically appealing (see also Tremaine 1986). It was indeed shown to lead naturally to the $`f_{\mathrm{}}`$ models. However, this route is not fully satisfactory from the mathematical point of view, especially since it involves an explicit approximation for the orbital period that is applicable only to the low binding energy limit of quasi-Keplerian orbits. The second route is straightforward from the mathematical point of view, being based on the classical Boltzmann entropy and on the assumed explicit conservation of a third quantity $`Q`$, in addition to the total mass $`M`$ and to the total energy $`E_{tot}`$. It was shown to lead to an analytically different family of models (that we may call the $`f^{(\nu )}`$ models; see definition in Sect. 3 below), with qualitative properties very similar to those of the $`f_{\mathrm{}}`$ models. Those models were not studied much further and did not receive great attention, not only because the relevant distribution function is not as simple as that of the $`f_{\mathrm{}}`$ models, but especially because the conservation of $`Q`$ could only be justified approximately by inspection of a number of N-body simulations, without a clear-cut physical justification (see Stiavelli and Bertin 1987; in contrast, the conservation of the additional $`𝐀𝐁`$ invariant sometimes invoked in plasma physics is rather transparent; see Chandrasekhar and Woltjer 1958).

In this paper we take advantage of the simple statistical mechanics foundation of the $`f^{(\nu )}`$ family of unbound partially relaxed stellar systems to explore the possibility of a thermodynamical description of stellar-dynamical models that are endowed with realistic properties.

## 2 Comment on the evolution of elliptical galaxies

Before proceeding to illustrate the results of this paper, we make a short digression in order to bring out the connections between the present analysis and the evolution of elliptical galaxies. We start by recalling that, formally, the sequence of King (1966) models is one special family of solutions of the collisionless Boltzmann equation. Yet, it is recognized to provide a reasonable description of the current properties of globular clusters (see Djorgovski and Meylan 1994), within a framework where these stellar systems continually evolve as a result of a variety of mechanisms (among which star evaporation and disk shocking; see Vesperini 1997 and references therein) and where the paradigm of the gravothermal catastrophe can be applied (see Spitzer 1987). Of course, it is well known that the level of internal collisionality in globular clusters is relatively high, so that the above approach is quite natural.

In contrast, one might at first think of dismissing the possibility that the paradigm of the gravothermal catastrophe should be of interest for the study of elliptical galaxies, because these large stellar systems lack the desired level of collisionality, judging from the estimate of the relevant star-star relaxation times. Here, following the spirit of earlier investigations (starting with Lynden-Bell and Wood 1968), we note that real elliptical galaxies are actually complex systems the evolution of which goes well beyond the idealized framework of the collisionless Boltzmann equation. In other words, splitting their description into past (formation) processes and present (mostly collisionless equilibrium) conditions should be considered only as an idealization introduced in order to assess the properties that define their current basic state.

In practice, elliptical galaxies are expected to be in a state of continuous evolution, for which we can list several specific dynamical causes: (1) Left-over granularity of the stellar system itself from initial collapse. Clumps of stars are likely to continue to form and dissolve in phase space even after the system has reached an approximate steady state. This acts as internal collisionality thus making some relaxation proceed even at current epochs. Indeed, numerical simulations of violent (partial) relaxation show that some evolution continues well after the initial collapse has taken place. (2) Drag of a system of globular clusters or other heavier objects towards the galaxy center. A globular cluster system or the frequent capture of small satellites (mini-mergers) may provide an internal heating mechanism associated with the process of dynamical friction by the stars on the heavier objects (Bertin at al., 2002). (3) Long-term action of tidal interactions of the galaxy with external objects. (4) Presence of gas in various phases (cold, warm, and hot). Significant cooling flows have been observed in bright ellipticals. Traditionally, studies of processes of this kind focus on the dynamics of the cooling gas and keep the background stellar system as ‘frozen’. In reality, energy and mass exchanges take place between the stellar system and the interstellar medium. (5) Interaction between the galactic nucleus and the galaxy. A number of interesting correlations have been found between the properties of galaxy nuclei and global properties of the hosting galaxies (e.g., see Pellegrini 1999). These correlations suggest that significant energy exchanges are taking place between the galaxy and its nucleus. Eventually, if a sufficiently concentrated nucleus is generated, then star-star relaxation in the central regions may also become a significant cause of dynamical evolution.

All of the above are specific mechanisms that are expected to make elliptical galaxies evolve in spite of their very long typical star-star relaxation time. Most of these processes are hard to model and to calculate in detail. As for the evolution of other complex many-body systems, it is hoped that thermodynamical arguments may help us identify general trends characterizing such evolution. This is the basic physical scenario in which the calculations presented in this paper are expected to be of interest for real elliptical galaxies.

## 3 Partially relaxed, unbound, finite mass systems from the Boltzmann entropy

Let us consider the standard Boltzmann entropy $`S=f\mathrm{ln}fd^3xd^3v`$ and look for functions that extremize its value under the constraint that the total energy $`E_{tot}=(1/3)Efd^3xd^3v`$, the total mass $`M=fd^3xd^3v`$, and the additional quantity

$$Q=J^\nu |E|^{3\nu /4}fd^3xd^3v$$
(1)
are taken to be constant. Here the functions $`E`$ and $`J^2`$ represent specific energy and specific angular momentum square of a single star subject to a spherically symmetric mean potential $`\mathrm{\Phi }(r)`$. As shown elsewhere (Stiavelli and Bertin, 1987), this extremization process leads to the following family of distribution functions

$$f^{(\nu )}=A\mathrm{exp}\left[aEd\left(\frac{J^2}{|E|^{3/2}}\right)^{\nu /2}\right],$$
(2)
where $`a`$, $`A`$, and $`d`$ are positive real constants. One may think of these constants as providing two dimensional scales (for example, $`M`$ and $`Q`$) and one dimensionless parameter; the dimensionless parameter can be taken to be $`\gamma =ad^{2/\nu }/(4\pi GA)`$. In principle, $`\nu `$ is any positive real number; in practice, we will focus on values of $`\nu 1`$. The $`f^{(\nu )}`$ non-truncated models are constructed by taking this form of the distribution function for $`E0`$, a vanishing distribution function for $`E>0`$, and by integrating the relevant Poisson equation under the condition that the potential $`\mathrm{\Phi }`$ be regular at the origin and behaves like $`GM/r`$ at large radii. This integration leads to an eigenvalue problem (see Appendix) for which a value of $`\gamma `$ is determined by the choice of the central dimensionless potential, $`\gamma =\gamma (\mathrm{\Psi })`$, with $`\mathrm{\Psi }=a\mathrm{\Phi }(r=0)`$.

The main point of the following analysis is the determination of the Boltzmann entropy $`S(M,Q,\mathrm{\Psi })`$ and of the total energy $`E_{tot}(M,Q,\mathrm{\Psi })`$ along the sequence of models, i.e. as a function of the concentration parameter $`\mathrm{\Psi }`$ defined above. These functions, at constant $`M`$ and $`Q`$, are illustrated in Fig. 1. They have been obtained by noting that, from the definitions of $`S`$ and $`f^{(\nu )}`$,

$$S=M\mathrm{ln}A+3aE_{tot}+dQ.$$
(3)
From the definitions $`Q=Aa^{9/4}d^{13/\nu }\widehat{Q}(\mathrm{\Psi })`$ and $`M=Aa^{9/4}d^{3/\nu }\widehat{M}(\mathrm{\Psi })`$ and the definition of $`\gamma `$, we can express the variables $`(A,a,d)`$ in terms of the variables $`(M,Q,\mathrm{\Psi })`$ and thus find that the entropy per unit mass can be written as $`S/M=S_0(M,Q)+\sigma (\mathrm{\Psi })`$, where $`S_0`$ is constant when the values of $`M`$ and $`Q`$ are fixed, with

$$\sigma =\mathrm{ln}\left(\widehat{M}^{\frac{4\nu 6}{5\nu }}\widehat{Q}^{\frac{6}{5\nu }}\gamma ^{\frac{9}{5}}\right)+\frac{3\widehat{E}}{\widehat{M}}+\frac{\widehat{Q}}{\widehat{M}}.$$
(4)
Here $`\widehat{E}=\widehat{E}(\mathrm{\Psi })`$ is the dimensionless total energy defined from $`E_{tot}=Aa^{13/4}d^{3/\nu }\widehat{E}`$. From the identity $`aE_{tot}/M=\widehat{E}/\widehat{M}`$ and the expression of $`a=a(M,Q,\mathrm{\Psi })`$ obtained previously, we find $`E_{tot}/M=H(M,Q)ϵ(\mathrm{\Psi })`$, with:

$$ϵ=\gamma ^{\frac{4}{5}}\widehat{M}^{\frac{9\nu +4}{5\nu }}\widehat{Q}^{\frac{4}{5\nu }}\widehat{E}.$$
(5)
The factor $`H(M,Q)`$ is a constant when $`M`$ and $`Q`$ are taken to be constant. The quantities $`\gamma (\mathrm{\Psi })`$, $`\widehat{M}(\mathrm{\Psi })`$, $`\widehat{Q}(\mathrm{\Psi })`$, and $`\widehat{E}(\mathrm{\Psi })`$ that enter the expression of $`\sigma `$ and $`ϵ`$ depend only on $`\mathrm{\Psi }`$ and are evaluated numerically on the equilibrium sequence.

This completes the derivation that allows us to draw the analogy with the classical paper of Lynden-Bell and Wood (1968). This step, straightforward for the $`f^{(\nu )}`$ models, is by itself interesting and new. In fact, other attempts at applying the paradigm of the gravothermal catastrophe to stellar dynamical equilibrium sequences were either based on an unjustified ansatz for the identification of the relevant temperature (e.g., see Appendix V in the article by Lynden-Bell and Wood 1968; Katz 1980; Magliocchetti et al. 1998) or on the use of non-standard entropies (for less realistic models; Chavanis 2002).

## 4 The high concentration regime: gravothermal catastrophe

When the $`f^{(\nu )}`$ models were constructed (Stiavelli and Bertin, 1987), it was immediately realized that they have general properties similar to those of the $`f_{\mathrm{}}`$ models (Bertin and Stiavelli, 1984); in particular, for values of $`\nu 1`$, sufficiently concentrated models along the sequence tend to settle into a “stable” overall structure, except for the development of a more and more compact nucleus, as the value of $`\mathrm{\Psi }`$ increases, and are characterized by a projected density profile very well fitted by the $`R^{1/4}`$ law characteristic of the surface brightness profile of bright elliptical galaxies. This property is illustrated in Fig. 2.

Now, by inspection of Fig. 1 and by analogy with the study of the isothermal sphere (Lynden-Bell and Wood, 1968), we can identify the location at $`\mathrm{\Psi }9`$ as the location for the onset of the gravothermal catastrophe. This sequence of models thus has the surprising result that the value of $`\mathrm{\Psi }`$ that defines the onset of the gravothermal catastrophe is precisely that around which the models appear to become realistic representations of bright elliptical galaxies. We leave to other papers (see Bertin and Stiavelli 1993) the detailed discussion of the issues that have to be addressed when comparison is made with the observations.

We note that in this regime of high concentration the general properties of the gravothermal catastrophe are reasonably well recovered by the use of the ansatz that the temperature parameter conjugate to the total energy is $`a`$, a quantity directly related to the velocity dispersion in the central regions. Basically, this was the ansatz made in the discussion of the possible occurrence of the gravothermal catastrophe for the King models or for other sequences of models (e.g., see Lynden-Bell and Wood 1968, Katz 1980, Magliocchetti et al. 1998). Here we have proved that the application of a rigorous derivation, which is available in our case, gives rise to relatively modest quantitative changes in the $`(E_{tot},1/T)`$ diagram for values of $`\mathrm{\Psi }`$ close to and beyond the onset of the catastrophe (see Fig. 3). However, in Sect. 6 we will draw the attention to an interesting, qualitatively new phenomenon missed in the previous derivations based on the use of the $`a`$-ansatz.

In passing, we note that in this regime of relatively high concentrations, the $`f^{(\nu )}`$ models possess one intrinsic property that makes them more appealing than the widely studied $`f_{\mathrm{}}`$ models. This is related to the way the models compare to the phase space properties of the products of collisionless collapse, as observed in N-body simulations (van Albada, 1982). In fact, one noted unsatisfactory property of the concentrated $`f_{\mathrm{}}`$ models was their excessive degree of isotropy with respect to the models produced in the simulations.111 Merritt et al. (1989) stressed this point and thus argued that a better representation of N-body simulations would be obtained by considering the $`f_{\mathrm{}}`$ family of models extended to the case of negative values of the coefficient $`a`$ multiplying the energy in the exponent. Unfortunately, their proposed solution, in terms of models characterized by such a peculiar phase-space structure, is unable to reproduce both the core velocity distribution observed in numerical experiments and the modest amount of radial anisotropy revealed by observed line profiles. In addition, their proposed solution is not viable because for negative values of $`a`$ the radial anisotropy level is so high that the models are violently unstable, on an extremely short timescale, with the result that their structural properties would be drastically changed by rapid evolution (Stiavelli and Sparke, 1991; Bertin et al., 1994). Here we can easily check that the anisotropy level of the concentrated $`f^{(\nu )}`$ models, while still within the desired (radial orbit) stability boundary and still consistent with the modest amount of radial anisotropy revealed by the observations, seems much closer to that resulting from N-body simulations of collisionless collapse; in particular, the anisotropy radius $`r_\alpha `$, defined from the relation $`\alpha (r_\alpha )=1`$, with $`\alpha =2(v_\theta ^2+v_\varphi ^2)/v_r^2`$, is close to the half-mass radius $`r_M`$ (while for the $`f_{\mathrm{}}`$ models it is about three times as large). This is illustrated in Fig. 4. In any case, we would like to emphasize that, while we have been clearly taking inspiration from simulations of collisionless collapse, our main interest is in comparing the structure of our models with that of observed objects rather than in providing a detailed fit to the results of N-body simulations.

## 5 The intermediate concentration regime: the $`R^{1/4}`$ law and deviations from it

The intermediate concentration regime (the precise point that marks the low concentration regime will be identified in the next Section) is a regime where the models appear to be stable, with respect not only to the gravothermal catastrophe (following the arguments provided earlier; but we should recall that the catastrophe is expected to require a sufficiently high level of effective collisionality in order to take off) but also to other instabilities (see the discussion given by Bertin et al. 1994 and references therein). The relatively wide variations, between $`\mathrm{\Psi }=3.5`$ and $`\mathrm{\Psi }=9`$, in all the representative quantities that characterize the equilibrium models suggest that this part of the sequence could be used to model the weak homology of bright elliptical galaxies (see Bertin et al. 2002), much like the sequence of King models is able to capture observed systematic variations in the structure of globular clusters (see Djorgovski and Meylan 1994).

## 6 The low concentration regime: negative global temperature and radial orbit instability

The low concentration regime is marked by an unexpected and significant difference with respect to the low concentration limit of the classical isothermal sphere (Bonnor, 1956; Lynden-Bell and Wood, 1968). In fact, while the classical case reduces to the ideal non-gravitating gas, to which Boyle’s law applies, for the $`f^{(\nu )}`$ models the system remains self-gravitating, although it develops a wide core in the density distribution. A clear-cut proof of this difference is given by inspecting the behavior of the global temperature $`T`$, identified from the thermodynamical definition $`T=1/(S/E_{tot})`$. While the temperature defined by the a-ansatz remains obviously positive definite, by definition, if we look at Fig. 1 we see that the global temperature $`T`$ changes sign at $`\mathrm{\Psi }3.5`$. This marks a drastic qualitative deviation from the classical studies.

Here we note a curious coincidence of this transition value of $`\mathrm{\Psi }`$ with the value around which the sequence is bound to change its stability properties with respect to the radial orbit instability (Polyachenko and Shukhman, 1981). Indeed, the location where the sequence is expected to become unstable in this regard is precisely that defined by $`\mathrm{\Psi }3.5`$, as can be judged from inspection of Fig. 4; around those values of $`\mathrm{\Psi }`$ the level of radial anisotropy, as measured by $`2K_r/K_T`$ reaches the threshold value of $`1.82`$, known to be sufficient for the excitation of the instability (the precise value of $`2K_r/K_T`$ corresponding to marginal stability is model dependent; for some sequences the reported value is below the range suggested by Polyachenko and Shukhman 1981).

These clues appear to be interesting and important, but more work is required before a final claim can be made that there is indeed a direct relation between the dynamical radial orbit instability and the fact that the system possesses a negative global temperature, as we found based on the simple work presented here.

## 7 Conclusions

A relatively straightforward and simple thermodynamical description of an equilibrium sequence constructed earlier and known to possess realistic characteristics with respect to bright elliptical galaxies shows that, on the one side, the paradigm of the gravothermal catastrophe may be adequate to explain the occurrence of realistic properties in models of collisionless stellar systems and, on the other side, a long-known dynamical instability might turn out to be interpreted in terms of a thermodynamical argument.

Probably the main open question regarding the models discussed in this paper, partially addressed in previous studies (see Stiavelli and Bertin 1987), is to what extent the quantity $`Q`$ is actually reasonably well conserved during violent (partial) relaxation. It is likely that a thorough investigation of this issue may give indications that the quantity is best conserved only in certain ranges of $`\nu `$. Studies of this type, combined with other dynamical and thermodynamical considerations, may turn out to lead to the identification of a family of equilibrium models with optimal behavior with respect to statistical mechanics, with respect to what we know about collisionless collapse, and with respect to the problem of providing a realistic representation of bright elliptical galaxies.

We would like to thank L. Ciotti, M. Stiavelli, and T. van Albada for interesting discussions and suggestions. This work has been partially supported by MIUR of Italy.

## Appendix A Numerical integrations

###

The $`f^{(\nu )}`$ models are constructed by solving the Poisson equation

$$\frac{1}{\widehat{r}^2}\frac{d}{d\widehat{r}}\widehat{r}^2\frac{d}{d\widehat{r}}\widehat{\mathrm{\Phi }}(\widehat{r})=\frac{1}{\gamma }\widehat{\rho }(\widehat{r},\widehat{\mathrm{\Phi }}),$$
(A1)
in which $`\gamma `$ is considered as an eigenvalue to be determined by imposing the two natural boundary conditions $`\widehat{\mathrm{\Phi }}(0)=\mathrm{\Psi }`$ and $`\widehat{\mathrm{\Phi }}(\widehat{r})\widehat{M}(\widehat{r})/(4\pi \gamma \widehat{r})`$ as $`\widehat{r}\mathrm{}`$. Here the symbol $`\widehat{}`$ indicates that the quantity is suitably expressed in dimensionless form.

###

We have computed the two-dimensional integral for the density with an adaptive seven-point scheme (Berntsen et al., 1991) in order to properly handle the presence of a peaked integrand for certain values of the pair $`(\widehat{r},\widehat{\mathrm{\Phi }})`$. The Poisson equation has then been solved with a fourth order Runge-Kutta code by starting from $`\widehat{r}=0`$ with a seed value for $`\gamma `$ and iterating the procedure until the boundary condition at large radii is matched within a certain accuracy. Finally, we have proceeded to calculate the global quantities $`\widehat{M}(\mathrm{\Psi })`$, $`\widehat{Q}(\mathrm{\Psi })`$, and $`\widehat{E}(\mathrm{\Psi })`$ from their definitions.

###

In order to check the accuracy of the numerical integration we have performed the following tests: (1) The virial theorem is satisfied with accuracy of the order $`10^6`$ or better; (2) The integrated mass (from its definition) and the mass derived from the asymptotic behaviour of the potential at large radii are the same with accuracy of the order $`10^4`$; (3) The expression for $`\widehat{\mathrm{\Phi }}(\widehat{r})`$ at large radii to two significant orders in the relevant asymptotic expansion has been checked to be correct with an accuracy from $`10^3`$ to $`10^4`$; (4) The asymptotic analysis allows us to estimate the contributions to $`\widehat{M}`$, $`\widehat{Q}`$, and $`\widehat{E}`$ external to a sphere of large radius $`R`$; this has been checked to help improve the numerical determination of the relevant global quantities.

###

We estimate that the final relative error in the quantities along the equilibrium sequence is of the order of some parts times $`10^4`$ for $`\widehat{M}`$ and $`\widehat{Q}`$ and some parts times $`10^5`$ for $`\widehat{E}`$. The total energy is less sensitive to the finite radius truncation error, due to its $`1/R^2`$ convergence. The propagation of these errors leads to the error bars plotted in Fig. 1, where the entropy $`\sigma `$ is the more difficult to determine with good accuracy.


text:
# Model-independent dark energy differentiation with the ISW effect

## Abstract

We study the integrated Sachs-Wolfe effect using a model-independent parameterization of the dark energy equation of state, w(z). Cosmic variance severely restricts the class of models distinguishable from $`\mathrm{\Lambda }`$CDM unless w(z) currently satisfies $`w_Q^o>0.8`$, or exhibits a rapid, late-time, transition at redshifts $`z<3`$. Due to the degeneracy with other cosmological parameters, models with a slowly varying w(z) cannot be differentiated from each other or from a cosmological constant. This may place a fundamental limit on our understanding of the origin of the presently observed acceleration.

cosmology: CMB anisotropies

1) Introduction – Cosmological distance measurements of high redshift type Ia supernova, combined with measurements of the cosmological volume via galaxy cluster surveys, weak lensing tomography and quasar clustering can provide a new insight into the nature of the dark energy sourcing the present acceleration of the universe. In spite of their potential for breaking degeneracies between cosmological parameters such tests are limited by our ignorance of effects such as possible supernova evolution or nonlinear galaxy cluster physics. In this sense the cosmic microwave background (CMB) observations are uniquely pure, limited only by systematic experimental effects and cosmic variance (an unavoidable theoretical error). This purity motivates us to ask the following question: how much model-independent information can the CMB give us about the dark energy ? We will show that if the dark energy only clusters on very large scales then only models with equations of state, $`w(z)`$, which vary rapidly at low redshifts $`z<3`$ can be distinguished from the ‘concordance’ $`\mathrm{\Lambda }`$CDM model using the CMB alone. The importance of differentiating between a cosmological constant ($`\mathrm{\Lambda }`$) and dynamical models of dark energy such as one (or more) evolving scalar fields (‘Quintessence’ - Q) (WETT, ; RATRA, ; STEIN, ) can hardly be overstated. $`\mathrm{\Lambda }`$ is essentially anti-gravity, the other is a new long-range force. We tackle this differentiation issue using a very general parameterization of $`w(z)p/\rho `$ described in (CORAS2, ) and used in simplified form in (BASS, ). This form accurately encompasses most quintessence models, independent of the number of fields involved, and many other models for the acceleration. In principle cosmological distance measurements can distinguish between $`\mathrm{\Lambda }`$CDM and QCDM models by constraining the present value of the dark energy equation of state $`w_Q^o`$. This does not work however, since current CMB data suggest that $`w_Q^o`$ is close to the cosmological constant value of $`1`$ (CORASA, ; BACCI, ; BEAN, ; BASS, ; HANN, ). However, the situation is not hopeless since dark energy clusters on very large scales and hence can leave a distinctive contribution in the anisotropies of the Cosmic Microwave Background radiation (CMB). This is particularly clear in the case of a rapid late-time transition in $`w(z)`$. Indeed current data does (weakly) prefer exactly such a situation (BASS, ).
In this letter we focus on the relation between the dark energy properties and the integrated Sachs-Wolfe effect (ISW) in the CMB power spectrum. Although this is model dependent we have recently shown that an appropriate parameterization of the equation of state, $`w(z)`$, accurately describes a large class of models given by eq. (4) in ref. (CORAS2, ). Assuming that the dark energy contribution during the radiation era is negligible, then different dark energy models are specified by the vector $`\overline{W}=(w_Q^o,w_Q^m,a_c^m,\mathrm{\Delta })`$ which respectively specifies the equation of state today (o) and during the matter era (m); the scale factor where the equation of state changes from $`w_Q^m`$ to $`w_Q^o`$ and the width of the transition. We now discuss how these parameters affect the CMB and show that only a small range of $`\overline{W}`$s leave an imprint on the CMB distinguishable from $`\mathrm{\Lambda }`$CDM.
2) ISW effect in dark energy models – A late-time mechanism to generate anisotropies is due to CMB photons climbing in and out of evolving gravitational potentials (REES, ). During the matter dominated era the gravitational potential $`\mathrm{\Phi }`$ associated with the density perturbations is constant and there is no ISW effect. However in $`\mathrm{\Lambda }`$CDM models $`\mathrm{\Phi }`$ starts decaying at redshifts when $`\mathrm{\Lambda }`$ starts to dominate, producing large angular scale anisotropies (CRI, ). In dark energy scenarios the cosmic acceleration is not the only contribution to the decay of $`\mathrm{\Phi }`$: on large scales the clustering properties alter the growth rate of matter perturbations (CALD, ; MA2, ). It is the signal of this clustering HU0  that we are hunting in as model-indepenent a way as possible. We assume a flat spatial geometry and fix the value of the Hubble constant $`H_o=70`$ Km$`s^1`$$`Mpc^1`$, the scalar spectral index $`n=1`$, the baryon density $`\mathrm{\Omega }_bh^2=0.021`$ and the amount of matter (CDM) $`\mathrm{\Omega }_m=0.3`$. We can usefully distinguish two classes of models,

* first those with a slowly varying equation of state for which $`0<a_c^m/\mathrm{\Delta }<1`$, as in the case of the inverse power law potential (ZLATEV, ), and second:
* a rapidly varying $`w(z)`$, such as the ’Albrecht-Skordis‘ model (ALB, ) and the two exponential potential (BARRE, ), with $`a_c^m/\mathrm{\Delta }>1`$. This class also includes many interesting radical models such as vacuum metamorphosis (PARK, ), late-time phase transitions (HILL, ), and backreaction-induced acceleration WETT2 .

We show these two classes in fig.1. The red solid line corresponds to dark energy that tracks the dust during the matter era ($`w_Q^m=0.0`$) and evolves slowly toward $`w_Q^o=1`$, and the blue dotted line corresponds to a model with a rapid transition in its equation of state at $`a_c^m=0.1`$ ($`z=9`$). Given current data it is worth studying the case with $`w_Q^o=1`$, (since it is also the most difficult to di stinguish from $`\mathrm{\Lambda }`$CDM) whilst allowing the other parameters $`w_Q^m`$ and $`a_c^m`$, to vary.

Fig.2 shows the anisotropy power spectrum, $`C_l^{ISW}`$, produced through the integrated Sachs-Wolfe effect by a rapidly evolving ($`top`$ $`panels`$) and a slowly evolving ($`bottom`$ $`panels`$) equation of state; the red (solid) line corresponds to the $`\mathrm{\Lambda }`$CDM model. As we can see in the top left panel (fig.2a), varying $`a_c^m`$ can produce a strong ISW. The effect is larger if the transition in the equation of state occurs at redshifts $`z<3`$. On the other hand the $`C_l^{ISW}`$ is the same as in cosmological constant regime if $`a_c^m<0.2`$ ($`z>4`$). In the top right panel (fig.2b) we plot the ISW for two different values of $`w_Q^m`$, corresponding to $`w_Q^m=0.0`$ (dashed line) and $`w_Q^m=0.1`$ (dot-dash line). We note that the signal is larger if the quintessence field is perfectly tracking the background component. But as $`w_Q^m`$ diverges from the dust value the ISW effect becomes the same as in $`\mathrm{\Lambda }`$CDM. This means that even for rapidly varying $`w(z)`$ (small $`\mathrm{\Delta }`$), the ISW is distinguishable from that in the $`\mathrm{\Lambda }`$CDM scenario only if $`w(z)`$ during matter domination closely mimics the dust value and the transition occurs at low redshifts, $`z<3`$. We can see the amplitude of the integrated Sachs-Wolfe effect is smaller in slowly varying models ($`bottom`$ $`panels`$). As we expect the $`C_l^{ISW}`$ is independent of $`a_c^m`$ (fig.2c), since for these models a different value in the transition redshift does not produce a large effect on the evolution of the dark energy density. In fig.2d, the ISW power spectrum is large for $`w_Q^m=0.0`$ (dash line) and becomes smaller than the cosmological constant on horizons scales as $`w_Q^m`$ has negative values (dot-dashed line), and increases toward $`\mathrm{\Lambda }`$ for $`w_Q^m`$ approaching $`1`$. This class of models is then more difficult to distinguish from the $`\mathrm{\Lambda }`$CDM if the equation of state today is close to $`w_\mathrm{\Lambda }=1`$. This can be qualitatively explained noting that perfect tracking between dark energy and CDM causes a delay in the time when the gravitational potential starts to decay, compared to the case of $`\mathrm{\Lambda }`$CDM. This effect is stronger for models with rapidly varying equation of state since the rapid change in $`w_Q`$ produces a stronger variation in the gravitational potential.

3) CMB power spectrum – The imprint of the ISW effect in the CMB spectrum results in a boost in power at low multipoles. This affects the position of the first acoustic peak and the COBE normalization. In particular, since some of the anisotropies at large scales are produced at late times, as we have seen in the previous section, normalizing the spectrum relative to the COBE measurements produces a suppression of power at smaller scales. A simple way to characterize the amplitude of the CMB spectrum is to consider the height of the first three acoustic peaks relative to the power at $`l=10`$, i.e. $`H_i=C_{l_i}/C_{10}`$, and their multipole positions $`\mathrm{}_i`$, $`i=1,2,3`$ (HU, ). These numbers allows us to quantify the discrepancy between the dark energy models and $`\mathrm{\Lambda }`$CDM. In particular since $`H_1`$ depends on the amplitude of the power spectrum at the multipoles characteristic of the ISW effect, we expect $`H_1`$ to be more sensitive to Quintessence signatures. We have computed the CMB spectra for the class of models with $`w_Q^o=1`$ described in the previous section, determined the parameters, $`H_i`$, (plotted in fig.3) and compared them with those of the $`\mathrm{\Lambda }`$CDM spectrum. The rapidly varying models are shown in the top panels. We can see the strong ISW effects produced by changing $`a_c^m`$ are now evident in the large discrepancy between $`H_1`$ and $`H_1^\mathrm{\Lambda }`$ (blue line) (fig.3a): it can be larger then 20 per cent for $`a_c^m>0.6`$. The effect on $`H_2`$, $`H_3`$ is smaller. However varying $`w_Q^m`$ (fig.3b) produces a discrepancy of only order $`4`$ per cent on $`H_1`$, while $`H_2`$ and $`H_3`$ remain the same as in $`\mathrm{\Lambda }`$CDM. For a slowly varying equation of state, $`H_1`$, $`H_2`$ and $`H_3`$ are independent of $`a_c^m`$ (fig.3c). The dark energy imprint is only on $`H_1`$ for which the discrepancy to the $`\mathrm{\Lambda }`$ case is about 10 per cent. Such discrepancy decreases when changing the value of the quintessence equation of state during matter from $`w_Q^m=0.0`$ to $`w_Q^m=1`$ (fig.3d). Values of the equation of state today $`w_Q^o>1`$ imply a stronger ISW effect. Consequently the curves of fig.3 are shifted upwards. For instance in fig.4 we plot the class of models previously analysed, with $`w_Q^o=0.88`$. We note the same behavior as we vary the dark energy parameters, but the discrepancy with the $`\mathrm{\Lambda }`$CDM model is now larger. In fig.4d it is worth noticing the case $`w_Q^m=1`$, that corresponds to a model very similar to a ‘k-essence’ model MUKA . We can see that the relative difference with the $`\mathrm{\Lambda }`$CDM case is of the order of a few percent, in agreement with ERIC  for the same value of $`w_Q^o=0.88`$. At this point we ask the key question whether such differences are observable. We have shown that $`H_1`$ is a good estimator of the ISW effect, and that it is a tracer of the dark energy imprint on the CMB. However its estimation from the data will be affected by cosmic variance at $`l=10`$. Hence with even perfect measurements of the first acoustic peak the uncertainty on $`H_1`$ will be dominated by the $`30`$ per cent uncertainty due to cosmic variance. With the plots of fig.3 and fig.4 in mind, this means that if the present value of the equation of state is close to $`1`$, slowly varying dark energy models are hardly distinguishable from $`\mathrm{\Lambda }`$CDM, while rapidly varying ones can produce a detectable signature only if the transition in the equation occurred at $`a_c^m>0.7`$, but in any case it will be difficult to constrain $`w_Q^m`$. The degeneracy of $`H_1`$ with the baryon density is marginal since varying $`\mathrm{\Omega }_bh^2`$ mainly affects $`H_2`$ which is insensitive to dark energy effects. Hence only an accurate determination of the angular diameter distance, inferred from the location of the acoustic peaks, would allow detection of such deviations from the cosmological constant model. The relation between the position of the CMB peaks and the dark energy has been widely discussed ((DORAN2, ) and references therein). The shift of the multipole positions ($`\mathrm{}_i`$) of the acoustic peaks caused by the evolution of the dark energy in the class of models analysed in fig.3 can be seen in fig.5, where we plot the relative difference of $`l_1,l_2`$ and $`l_3`$ to the $`\mathrm{\Lambda }`$ case. We note that due to the additional shift induced on the first acoustic peak by the ISW effect the difference with the $`\mathrm{\Lambda }`$CDM model for the first peak is generally larger than for the second and third peaks. As with the comparison of the amplitude of the CMB spectrum, the largest effect is produced by models with a rapid transition occurring at small redshifts. However the degeneracy of the angular diameter distance, in particular with the value of Hubble constant and the amount of dark energy density, will limit our ability to put tight constraints on the dark energy parameters. Fortunately there are alternative ways in which these problems can be alleviated, for instance, cross-correlating the ISW effect with the large scale structure of the local universe (HUW, ; CO, ; BOUGHN, ). An efficient approach would also be to combine different observations in order to break the degeneracies with the cosmological parameters (WAS, ; FRIEMAN, ).

4) Conclusions – The next generation of high resolution CMB experiments will measure the anisotropy power spectrum with an accuracy close to theoretical limits. It is therefore of particular interest to study the sensitivity of CMB observations to the effects produced by a dark energy component in the CMB. A quintessential contribution leaves a distinctive signature in the ISW effect. On the other hand such an imprint occurs at low multipoles, consequently cosmic variance strongly limits the possibility of differentiating dark energy models from a cosmological constant. In particular we find that for values of $`1w_Q^o<0.8`$, the only distinguishable cases are those with a rapidly varying equation of state. In fact for the slowly varying models, the ISW is the same as in $`\mathrm{\Lambda }`$CDM and a deviation from $`w_Q^o=1`$ can be inferred only from an accurate determination of the location of the acoustic peaks. However such measurements will be affected by the degeneracy of the angular diameter distance with the value of the Hubble constant and the amount of dark energy density. Therefore CMB observations are insensitive to this class of models and we will be left with a fundamental uncertainty as to whether the cosmic acceleration is due to an evolving field or a cosmological constant, an issue of great theoretical importance for our understanding of the foundations of quantum gravity and string theory.

###### Acknowledgements.

We thank Rob Lopez for discussions and supplying part of the code upon which this work is based.
PSC is supported by a University of Sussex bursary. The research of BB and CU was supported by PPARC grant PPA/G/S/2000/00115.


text:
# Toy models for the falling chimney

## I Introduction

One of the most interesting demonstrations for an introductory mechanics course is the “Falling Chimney - Free Fall Paradox,” as it was named by Sutton in his classical book Demonstration Experiments in Physics.Sutton3  In the original version of this demonstration a ball is placed at one end of a uniform stick, which is pivoted at the other end and makes initially an angle of about $`30^{}`$ with the horizontal. The elevated end of the stick is suddenly dropped, together with the ball, thus showing a very counter-intuitive behavior. The falling end of the stick accelerates at a greater rate than the free-falling ball, proving that its acceleration is greater than $`g`$, the acceleration of gravity.

A simplified version of the experiment can be performed just with a meter stick and a coin. The stick is supported in the horizontal position by two fingers, placed near the two ends. A coin is set on the stick near one end, which is suddenly released. The effect is similar to the previous demonstration: the falling end of the rotating stick eventually acquires an acceleration greater than that of the freely falling coin, which loses contact with the stick surface and lags behind the falling stick.

A description of the first version of the experiment can be found in almost every book of physics demonstrations,Hilton2 ; Freier ; DickRae  sometimes with different names (“Falling Stick”, “Hinged Stick and Falling Ball”, and others). Photographic descriptions or even video-clips of this demo can be found on-line in several web-pages (see our web-page,kaoru  for a collection of related links).

In addition, countless papers exist in the literature; we have traced several of these, from the 1930’s to the present. Some of the earliest discussions can be found in ConstantinidesConstantinides  and LudekeLudeke  (as well in the book by SuttonSutton3 ), followed by many others.Hilton1 ; Young ; Phelps ; Theron ; Hartel  These concentrate mostly on the simple explanation of the effect, which relies on the concept of “center of percussion” of the rotating stick (a simple introduction to this concept can be found in BloomfieldBloomfield ). This particular point of the stick (located at a distance from the hinged end equal to two thirds of the length, for a uniform stick) is moving with the same acceleration as a particle under gravity, constrained to move along the same circular path. Points on the stick beyond the center of percussion descend with accelerations greater than that of particles freely moving under gravity, on their respective circular paths. As a consequence of this, if the initial angle formed by the stick with the horizontal is less than about $`35^{}`$, the end point will possess at all times a vertical component of the acceleration greater than $`g`$, producing the effect described above.

Several variations of the basic demonstration also exist, Sutton2 ; Miller ; Ficken ; Bartlett1 ; Adams ; Haber  the majority of which suggest attaching an additional mass to the rotating stick at different positions. The effect for the student or the viewer is even less intuitive than the original version: an additional mass placed near the end of the stick actually reduces the acceleration of the end point, affecting substantially the outcome of the experiment. In general the addition of a mass at any point on the stick will increase both the total torque on the system (thus increasing the rotational acceleration) and the moment of inertia of the system around the axis of rotation (resulting in a decreased rotational acceleration). The center of percussion of the stick still plays a key role: if the additional mass is placed beyond it, the effect of the increased moment of inertia dominates and the acceleration of the rotational motion will be reduced. If the mass is placed before the center of percussion, the increase in the torque will dominate and the rotational motion will be enhanced. The effect is null if the mass is placed exactly at the center of percussion (a complete discussion of this effect can be found in BartlettBartlett1  and Haber-SchaimHaber ).

The next logical step is to analyze the behavior of a real falling chimney. Almost invariably a tall chimney, falling to the ground like the stick in the previous discussion, will break in mid-air at some characteristic height. This is well documented in several photos reproduced in the literature, such as the one which appeared on the cover of the September 1976 issue of The Physics Teacher (other photos can be found in BundyBundy  and Bartlett,Bartlett2  or also on our web-pagekaoru ). The causes of such breaking, the height of the rupture point and the angle at which the breaking is most likely to occur, are the most natural questions which arise.

The first analysisSutton1  of this problem compared the fall of the real chimney to the fall of the hinged stick, but wrongly identified the center of percussion (at about two thirds of the height) as the probable point of rupture. ReynoldsReynolds  first identified the possible causes of the breaking with the shear forces and the bending moment originating within the structure of the toppling chimney. More detailed analyses were given by BundyBundy  and MadsenMadsen  (the most complete papers we found on the subject) while simplified explanations are also reported.Jones ; Bartlett2 ; Walker  It even appears in graduate student study guides,Cahn ; Cronin  although the chimney is shown bending the wrong way in one of these books.

In this paper we review the theory of the real falling chimney, outlined by Madsen,Madsen  aiming for a complete and clear explanation of this phenomenon in Sects. II and III. Then, in Sect. IV, we propose simple ways of using small scale models (literally toy models \- made with toy blocks and bricks) to test effectively the outlined theory. More information on these toy models can also be found on our web-site,kaoru  together with photos and movie clips of the experiments we have performed.

## II Rotational motion of the falling chimney

The rotational motion of a falling chimney under gravity is equivalent to that of the falling hinged stick of Sect. I. We can describe it as in Fig. 1,

where we use polar coordinates $`r`$ and $`\theta `$ (with $`\widehat{e}_r`$ and $`\widehat{e}_\theta `$ as unit vectors) for the position of an arbitrary point $`A`$ on the longitudinal axis of the chimney, measuring the angle $`\theta `$ from the vertical direction. We treat the chimney as a uniform rigid body of mass $`m`$ and height $`H`$, under the action of its weight $`𝐖=m𝐠`$, applied to the center of gravity (basically the center of mass -CM- of the body), and a force $`𝐅`$ exerted by the ground on the base of the chimney, assumed to act at a single point (we neglect air resistance, or any other applied force). In plane polar coordinates:

$`𝐖`$ $`=W_r\widehat{e}_r+W_\theta \widehat{e}_\theta =mg\mathrm{cos}\theta \widehat{e}_r+mg\mathrm{sin}\theta \widehat{e}_\theta `$ (1a)
$`𝐅`$ $`=F_r\widehat{e}_r+F_\theta \widehat{e}_\theta .`$ (1b)
The moment of inertia of the chimney can be approximated with the one for a uniform thin rod, with rotation axis perpendicular to the length and passing through one end:111This is in general a good approximation. For a right circular hollow cylinder of outer and inner radii $`r_1`$, $`r_2`$ and length $`H`$, hinged at one end, the moment of inertia is $`I=m\left(\frac{r_1^2+r_2^2}{4}+\frac{1}{3}H^2\right)\frac{1}{3}mH^2`$ for $`r_1,r_2H`$. For a typical chimney $`r_1\frac{H}{10}`$, thus $`\frac{r_1^2/4}{H^2/3}0.0075`$, i.e., a correction of less than one percent. The correction due to $`r_2`$ is even smaller.

$$I=\frac{1}{3}mH^2.$$
(2)
Applying the torque equation $`I\stackrel{..}{\theta }=\tau _z`$, for a rotation around the origin, with an external torque given by $`\tau _z=mg\frac{H}{2}\mathrm{sin}\theta `$, we find the angular acceleration

$$\stackrel{..}{\theta }=\frac{\tau _z}{I}=\frac{3}{2}\frac{g}{H}\mathrm{sin}\theta .$$
(3)
A simple integration, using $`\stackrel{..}{\theta }=\frac{d\stackrel{.}{\theta }}{dt}=\frac{d\stackrel{.}{\theta }}{d\theta }\stackrel{.}{\theta }=\frac{3}{2}\frac{g}{H}\mathrm{sin}\theta `$, gives the angular velocity

$$\stackrel{.}{\theta }^2=3\frac{g}{H}\left(1\mathrm{cos}\theta \right),$$
(4)
assuming that the chimney starts moving from rest and is initially in the vertical direction. A further integration of Eq. 4, can lead to $`\theta (t)`$ in terms of elliptic integrals.

We recall that the acceleration in polar coordinates can be written as $`𝐚`$=$`\stackrel{..}{𝐫}`$=$`(\stackrel{..}{r}r\stackrel{.}{\theta }^2)\widehat{e}_r+(r\stackrel{..}{\theta }+2\stackrel{.}{r}\stackrel{.}{\theta })\widehat{e}_\theta `$, so that, for a point $`A`$ at a fixed distance $`r`$ from the origin, it becomes

$$𝐚=a_r\widehat{e}_r+a_\theta \widehat{e}_\theta =r\stackrel{.}{\theta }^2\widehat{e}_r+r\stackrel{..}{\theta }\widehat{e}_\theta .$$
(5)
For a point at two thirds of the height, $`r=\frac{2}{3}H`$, combining Eqs. 3 and 5 we get $`a_\theta (r=\frac{2}{3}H)=\frac{2}{3}H\stackrel{..}{\theta }=g\mathrm{sin}\theta `$, proving that this particular point is the center of percussion of the body, as already mentioned in Sect. I.

The torque equation allowed us to determine the angular acceleration of the motion in Eq. 3. We can use this result and Newton’s second law for the motion of the center of mass (CM) of the whole chimney to determine the unknown force $`𝐅`$ at the base. The vectorial equation is

$$m\stackrel{..}{𝐫}_{CM}=𝐖+𝐅$$
(6)
which, for $`r=\frac{H}{2}`$, splits into radial and angular equations,

$`m{\displaystyle \frac{H}{2}}\stackrel{.}{\theta }^2`$ $`=F_rmg\mathrm{cos}\theta `$ (7a)
$`m{\displaystyle \frac{H}{2}}\stackrel{..}{\theta }`$ $`=F_\theta +mg\mathrm{sin}\theta .`$ (7b)
Using Eqs. 3 and 4, the two components of the force $`𝐅`$ are easily determined:
$`F_r`$ $`={\displaystyle \frac{5}{2}}mg\left(\mathrm{cos}\theta {\displaystyle \frac{3}{5}}\right)`$ (8a)
$`F_\theta `$ $`={\displaystyle \frac{1}{4}}mg\mathrm{sin}\theta .`$ (8b)
## III Internal forces and bending moment

We now move to the analysis of the internal forces which develop inside the structure of the falling chimney. The resulting stresses and bending moment are the causes of the rupture of the toppling chimney. Consider, as in Fig. 2,

an arbitrary lower portion of the chimney of height $`r`$ (as opposed to the total height $`H`$) and the forces acting on this part of the structure due to the upper portion and the base. The weight of the lower portion is now $`W(r)=mg\frac{r}{H}`$ (assuming again a uniform structure, so that the weight is proportional to the height of the considered portion) and it is applied to the center of gravity of this lower portion (at a distance $`\frac{r}{2}`$ from the origin). In polar coordinates:

$$𝐖(r)=W_r(r)\widehat{e}_r+W_\theta (r)\widehat{e}_\theta =mg\frac{r}{H}\mathrm{cos}\theta \widehat{e}_r+mg\frac{r}{H}\mathrm{sin}\theta \widehat{e}_\theta .$$
(9)
The force $`𝐅`$, applied at the base, is still the same as in Eqs. 8a, 8b, but we have to add the action of the upper part on the lower portion. We follow here the general analysis of the internal forces and moments which can be found in every textbook on Statics (see for example Statics ; Statics2 ) and which can be easily adapted to our case.

The distribution of all the internal forces, at the cross section being considered, can be equivalently described by a resultant force and a resultant moment acting at a specific point of the cross section (typically the “centroid” of the sectioned area, in our case simply the central point of the section, on the longitudinal axis). In particular, the resultant force can be decomposed into a transverse shearing force $`𝐒=S_\theta \widehat{e}_\theta `$, and a longitudinal stress force (tension or compression) $`𝐏=P_r\widehat{e}_r`$, applied as in Fig. 2, at the cross section between the upper and lower portions, and assumed positive in the $`\widehat{e}_\theta `$, $`\widehat{e}_r`$ direction respectively.

In addition, we have to consider the resultant moment of the forces at the cross section, which is usually called the “bending moment” $`𝐍_b`$, because its effect will ultimately result in bending the structure. It is represented in the picture by the curved arrow. Since we treat this as a plane problem, the bending moment can only have a component perpendicular to the plane of the figure, i.e., in the $`z`$ direction. No other components are considered here, in particular we assume that no torsional moment exist in the structure, which would tend to twist the chimney around its longitudinal axis.

The bending moment $`𝐍_b=N_b\widehat{e}_z`$ can be thought as originating from a couple of forces, $`𝐟`$ and $`𝐟`$, that can be regarded as applied to the leading and trailing edge of the structure, at the cross section considered. This couple of forces is shown explicitly in the papers by BundyBundy  and Madsen,Madsen  but we prefer to use directly $`𝐍_b`$ in our treatment, because the bending moment is the result of the whole distribution of forces at the cross section considered. The two small insets inside Fig. 2 explain the definition of the bending moment in terms of the couple of forces $`𝐟`$ and $`𝐟`$. We also show the resulting deformation of the structure due to a “clockwise” (diagram a), or a “counter-clockwise” bending moment (diagram b). The latter case will be the actual deformation of the falling chimney. $`N_b`$ will be assumed to be positive if it acts as in the figure, i.e., a positive component of the torque in the $`z`$ direction (we assume here the use of a right-handed system of coordinate axis). In the following we will refer to $`N_b`$ as the bending moment, acting on the lower portion of the chimney.222The bending moment and the stress forces can be thought as applied to the center (centroid) of the cross sectional area we are considering, on the face belonging to the lower portion of the structure. Equal and opposite forces and moments would originate on the face belonging to the upper portion.

Again, we will consider the torque equation and the second law for the motion of the center of mass (located at $`\frac{r}{2}`$) for just the lower portion of the chimney (of mass $`\frac{r}{H}m`$). It is better to analyze the CM motion first. The vector equation

$$m\frac{r}{H}\stackrel{..}{𝐫}_{CM}=𝐖(r)+𝐅+𝐏+𝐒,$$
(10)
will split into the radial and angular directions,

$`m{\displaystyle \frac{r^2}{2H}}\stackrel{.}{\theta }^2`$ $`=`$ $`{\displaystyle \frac{3}{2}}mg(1\mathrm{cos}\theta ){\displaystyle \frac{r^2}{H^2}}=mg{\displaystyle \frac{r}{H}}\mathrm{cos}\theta +{\displaystyle \frac{5}{2}}mg\left(\mathrm{cos}\theta {\displaystyle \frac{3}{5}}\right)+P_r`$ (11a)
$`m{\displaystyle \frac{r^2}{2H}}\stackrel{..}{\theta }`$ $`=`$ $`{\displaystyle \frac{3}{4}}mg\mathrm{sin}\theta {\displaystyle \frac{r^2}{H^2}}=mg{\displaystyle \frac{r}{H}}\mathrm{sin}\theta {\displaystyle \frac{1}{4}}mg\mathrm{sin}\theta +S_\theta ,`$ (11b)
having used Eqs. 3, 4, 8a, 8b, and 9. We can solve for the longitudinal and transverse forces
$`P_r`$ $`=`$ $`{\displaystyle \frac{1}{2}}mg\left(1{\displaystyle \frac{r}{H}}\right)\left[\left(5+3{\displaystyle \frac{r}{H}}\right)\mathrm{cos}\theta 3\left(1+{\displaystyle \frac{r}{H}}\right)\right]`$ (12a)
$`S_\theta `$ $`=`$ $`{\displaystyle \frac{3}{4}}mg\mathrm{sin}\theta \left({\displaystyle \frac{r^2}{H^2}}{\displaystyle \frac{4}{3}}{\displaystyle \frac{r}{H}}+{\displaystyle \frac{1}{3}}\right),`$ (12b)
which depend on the fraction of height $`\frac{r}{H}`$, the angle of rotation $`\theta `$, and also the total weight $`mg`$. Following the analysis by Bundy,Bundy  we plot these two forces in Figs. 3 and 4 respectively, normalized to the total weight $`mg`$, as a function of the height fraction, for several angles.
From Fig. 3 we see that $`P_r`$ is negative (a compression) for smaller angles, but eventually becomes positive (a tension) for angles greater than about $`45^{}`$. $`P_r`$ also depends critically on $`\frac{r}{H}`$ (for $`\theta =0^{}`$, $`P_r`$ represents simply the compression due to the weight of the upper part acting on the lower part). This longitudinal force will be combined later with the bending moment to determine the total stress at the leading and trailing edges, which is the most typical cause of the rupture.

In Fig. 4 we plot the (transverse) shear force $`S_\theta `$, which can be the other leading cause of rupture. It is easily seen that, for any considered angle, the magnitude of the shear force $`\left|S_\theta \right|`$ has an absolute maximum at $`\frac{r}{H}=0`$ (and a positive value), meaning that large shear forces, in the $`\widehat{e}_\theta `$ direction, usually originate near the base. The shear force is always zero at one third of the height, and $`\left|S_\theta \right|`$ also has a (relative) maximum at $`\frac{2}{3}H`$ (with a negative value, therefore $`S_\theta `$ is in the $`\widehat{e}_\theta `$ direction), but this value is smaller than the one near the base.

From this analysis, it is typically concluded that if the structure breaks from shear stress alone, this is usually more likely to happen near the base. This can be seen for example in the already mentioned cover photo of the September 1976 issue of The Physics Teacher,Bartlett2 ; kaoru  showing the fall of a chimney in Detroit. The two ruptures at the bottom are likely due to shear forces, while a third rupture can be seen at about $`r=0.47H`$, and this is due to the combination of bending moment and longitudinal force $`P_r`$, as we will explain in the following. More photos and detailed pictorial descriptions of chimney ruptures can be found in the paper by Bundy.Bundy ; kaoru

The “bending moment” $`N_b`$ can be calculated from the torque equation $`I(r)\stackrel{..}{\theta }=\tau _z`$, where now $`I(r)=\frac{1}{3}m\frac{r}{H}r^2`$ is the moment of inertia of just the lower part. $`\stackrel{..}{\theta }`$ will come from Eq. 3, and the total external torque is now $`\tau _z=\frac{r}{2}W_\theta (r)+rS_\theta +N_b`$. Using Eqs. 9 and 12b, we can solve the torque equation for the bending moment, obtaining

$$N_b=\frac{1}{4}mg\mathrm{sin}\theta r\left(1\frac{r}{H}\right)^2,$$
(13)
or, in a non-dimensional form

$$\frac{N_b}{mgH}=\frac{1}{4}\mathrm{sin}\theta \frac{r}{H}\left(1\frac{r}{H}\right)^2,$$
(14)
which is plotted in Fig. 5, as a function of the height fraction and for various angles.

$`N_b`$ is always negative, showing that it is actually directed in the opposite way of Fig. 2 (or as in diagram b of Fig. 2). This particular direction of the bending moment will induce a tension in the leading edge of the chimney and a compression in the trailing edge. The structure will bend accordingly, with the concavity on the side of the trailing edge, and will eventually break in the way shown by the many existing photos. We can see that, for any angle, the bending moment is obviously zero at the bottom and at the top of the chimney, while it assumes its maximum absolute value at exactly one third of the height $`H`$. The bending moment alone would therefore induce a rupture at one third of the structure, but the total longitudinal stress at the leading edge is also due to the force $`P_r`$, as we will show next.

Another interesting relation between the bending moment $`N_b`$ and the shear force $`S_\theta `$ is that they are in general related by a simple derivative, i.e., $`S_\theta =\frac{N_b}{r}`$, as it is easy to check from Eq. 12b and Eq. 13.333A positive sign in the relation between bending moment and shear is usually reported in textbooks on Statics, due to a different choice of the sign of the moment. This is a well known relationship of the Statics of beams and other structural members (for a complete proof see for example HibbelerStatics2 ). It is a direct consequence of the equilibrium equations applied to an infinitesimal longitudinal portion of the beam: the change in bending moment along the beam is always equal to the shear force applied to that portion of the beam.

Finally, we can combine $`N_b`$ and $`P_r`$ to compute the total longitudinal stress on the cross sectional area between the lower and upper parts. We follow the theory of elasticity and deformations in beams, which can be found in classical treatises such as SommerfeldSommerfeld  and Landau-Lifshitz,Landau  or again in Bundy’s paper.Bundy  The longitudinal stress is maximum at the leading and trailing edges, located at the maximum distance from the longitudinal (centroidal) axis of the chimney which lies within the “neutral surface” of the structure, the surface which is neither stretched nor compressed.

For simplicity, we will only consider from here on, structures with uniform square cross section of side $`a`$, as this is the case of the toy models described in Sect. IV. In this case the stresses at the leading and trailing edge, $`\sigma _L`$ and $`\sigma _T`$ respectively, can be evaluated from

$$\sigma _{L/T}=\frac{P_r}{A}\frac{aN_b}{2J}$$
(15)
(the upper sign is for $`\sigma _L`$, the lower for $`\sigma _T`$) where $`A=a^2`$ is the area of the square cross section of side $`a,`$ with the factor $`\frac{a}{2}`$ representing the distance between the longitudinal axis, considered as the “neutral axis,” and the two edges. $`J=\frac{a^4}{12}`$ is the moment of inertia of the cross sectional area computed about the neutral axis (see SommerfeldSommerfeld  for details). Using also the expressions for $`P_r`$ and $`N_b`$, we obtain

$$\frac{\sigma _{L/T}a^2}{mg}=\frac{1}{2}\left(1\frac{r}{H}\right)\left[\left(5+3\frac{r}{H}\right)\mathrm{cos}\theta 3\left(1+\frac{r}{H}\right)\right]\pm \frac{3}{2}\frac{H}{a}\mathrm{sin}\theta \frac{r}{H}\left(1\frac{r}{H}\right)^2,$$
(16)
where we normalized $`\sigma _{L/T}`$, dividing by $`\frac{mg}{a^2}`$, in order to obtain a dimensionless quantity which is plotted in Fig. 6, as a function of the height fraction and for several angles.

This quantity depends also on the ratio $`\frac{H}{a}`$, which for a real chimney is of the order $`\frac{H}{a}10`$. For the toy models described in Sect. IV, the value of this ratio is even bigger: $`\frac{H}{a}2461`$, enhancing the contribution of the second term of Eq. 16, which comes from the bending moment $`N_b`$. In Fig. 6 we show the plot for $`\frac{H}{a}=24`$, but similar figures can be obtained for different values of the ratio.

In Eqs. 15-16, and in Fig. 6 the total stresses are considered positive if they represent tensions, negative if they are compressions. It is easily seen from the figure that the stress at the leading edge $`\sigma _L`$, is initially a compression, but eventually becomes a tension, constantly increasing for larger angles; $`\sigma _T`$ on the contrary is usually a compression. It is therefore the combination of these intense tensile stresses in the leading part of the chimney (and also compressions on the trailing side) that causes the rupture of the chimney.

This type of breaking is more likely to occur at the positive maximum value of $`\sigma _L`$. This maximum value depends critically on the $`\frac{r}{H}`$ ratio, for a certain angle of rupture $`\theta `$, so it is possible, just by looking at the maxima of the solid curves in Fig. 6 (marked by solid points), to roughly match the height of the rupture point to the angle at which the breaking started to occur. As far as the actual prediction of the point of rupture, this would obviously depend on the building materials and the construction of the chimney or tower, an analysis of which goes beyond the scope of this work.

It is interesting to note, from Fig. 6 again, that the stress $`\sigma _L`$ is not always maximum at one third of the height (as for the bending moment of Fig. 5). For small angles of about $`\theta 5^{}20^{}`$ it reaches a maximum for $`\frac{r}{H}0.40.5`$, while for larger angles it approaches the typical ratio $`\frac{r}{H}\frac{1}{3}`$. This means that if the chimney breaks early in its fall, for small angles, it is more likely to break near the center; on the contrary if the rupture occurs at larger angles, the breaking point is usually shifted toward one third of the height. This is the effect of the $`P_r`$ term in Eq. 15, which modifies the position of the maximum $`\sigma _L`$ (or $`\sigma _T`$).

As already mentioned, another factor to be considered is the ratio $`\frac{H}{a}`$, which is usually in the range $`\frac{H}{a}520`$ for real chimneys, but can be increased up to about $`\frac{H}{a}100`$ with our toy models. In Fig. 7 we plot the maxima of the leading edge stress curves, for several values of the ratio $`\frac{H}{a}=5,10,20,30,\mathrm{},100`$ (for the group of curves between the values of $`30`$ and $`100`$, the parameter is increased by $`10`$ units for each curve).

These maxima are plotted as continuous curves, showing the corresponding values of the ratio $`\frac{r}{H}`$ and the angle $`\theta `$, at which the structure is more likely to break, for a given value of $`\frac{H}{a}`$. In other words, connecting for example the solid points of Fig. 6, we would obtain a corresponding continuous curve in Fig. 7, with the angle of probable rupture on the vertical axis, instead of the normalized stress.

The dependence of the rupture point on these quantities can be noted in several of the existing photos of falling chimneys, when the breaking occurs due to the bending of the structure and not for the transverse shear stress near the base mentioned at the beginning of this section. The angle of rupture can be roughly estimated by measuring the angle the upper part of the chimney forms with the vertical in the photos. This angle tends in fact not to change much after the rupture, since the upper part falls without much additional rotation. For example the photos in the paper by BartlettBartlett2 ; kaoru  refer to chimneys with $`\frac{H}{a}10`$, breaking at about $`\theta 20^{}25^{}`$ and for $`\frac{r}{H}0.47`$, consistent with the values of Fig. 7, for the $`\frac{H}{a}=10`$ curve. Similar behavior can be seen in other photos,Bundy ; kaoru  since real chimneys tend to break very early in their fall, due to the intense stresses originating within their structure.

## IV Toy models

In this section we discuss our efforts to reproduce the effects described above, with the help of toy models of the falling chimney. These models were constructed with simple toy blocks of different type and size, and their fall was filmed with a digital camera, so that we could analyze the events frame by frame, to test the theory. Complete details on the type of blocks used, experimental settings and video-capture techniques, as well as the complete set of our video-recordings and still photos can be found on our web-site,kaoru  and in an upcoming publication.kaoru2

Bundy noted in his workBundy  that the use of a model to test a real chimney would be useless due to a “scale effect.” The stresses inside the chimney depend roughly on the scale of the object, so that real chimneys would develop bigger stresses than equivalent small-scale models, therefore breaking earlier in their fall. Nevertheless we found it interesting to reproduce these effects in small scale models to test especially the discussion based on Figs. 6 and 7. It was also supposed to be difficultBundy  to show these effects with toy models.444We did not find any use of toy models for this problem, either in the literature or in web-pages devoted to physics demos. We will be grateful to receive information about similar experiments, if any.

Fig. 8 is the first example of one of our toy models. We made a tower by simply stacking $`24`$ wooden toy blocks of cubic shape, for a total height $`H=0.76m`$, mass $`m=0.32kg`$ and a ratio $`\frac{H}{a}=24`$, the value used in Fig. 6. The tower was set into the falling motion by removing a support at the bottom, inducing a rotation without slipping at the bottom point. The picture clearly shows the “rupture” due to bending of the structure at exactly half the height, $`\frac{r}{H}=0.50`$, and for a small angle $`\theta 10^{}`$, which again can be estimated by measuring the angle the upper part of the chimney forms with the vertical direction. This is in good agreement with the position of the maximum for the solid-$`10^{}`$ curve in Fig. 6 and also with the data of Fig. 7, using the $`\frac{H}{a}=20`$ curve, showing that the theory is applicable also to these small-scale models.

Our second example is a taller tower ($`H=1.9m`$, $`m=0.65kg`$, and $`\frac{H}{a}=61`$) made with 100 plastic blocks of a very popular brand of toy bricks. The blocks are inserted on top of each other so that bending of the structure is allowed, but shear stress cannot possibly break the tower. The 100 toy blocks are arranged by color to subdivide the structure into three equal parts, and also the position of the center is marked. This time the rupture occurs for an angle around $`\theta 30^{}35^{}`$, and at the height ratio $`\frac{r}{H}=0.40`$. This is consistent with the $`30^{}`$ solid curve in Fig. 6, while the data from Fig. 7 (for $`\frac{H}{a}=60`$) would suggest a smaller angle of rupture. It is actually difficult, with this type of toy bricks, to estimate the angle at which the structure begins to bend.

We performed several other experiments, varying the dimensions of the towers, the type of blocks, always obtaining results consistent with the theory. We can conclude that it is actually easy to reproduce the bending and breaking of chimneys with small scale towers, and this type of experiment could be made part of an undergraduate laboratory class for rotational mechanics, with some minor adaptations and changes.kaoru2

## V Conclusion

In this paper we reviewed the theory of the falling chimney, showing that the rupture can be caused by either shear forces typically near the base, or by the bending of the structure which is caused primarily by the bending moment, but is also affected by the longitudinal stress force. In the latter case the breaking is more likely to occur between one third and one half of the height of the chimney.

This point of rupture is also related to the angle of rupture and this relationship can be verified in the many existing photographic reports of falling chimneys. We also constructed several small scale toy models, showing that it is possible with them to reproduce the dynamics of the fall. By examining photos taken during the fall of these models we were able to confirm the theoretical model outlined in this paper.

###### Acknowledgements.

This research was supported by an award from Research Corporation. The authors would like to thank Dr. V. Coletta for the many useful discussions regarding this paper. G.V. would like to acknowledge and thank his friend and physics teacher, Prof. G. Tonzig, whose excellent bookTonzig  “Cento Errori di Fisica” inspired the original idea for this paper.