text:
# 1 introduction

## 1 introduction

Model considerations are always very important for the concept clarification in theoretical and experimental physics. In particular, an exact solvable model may provide much more detailed informations which are useful for a complete understanding of the nature of the phenomena. Thus, in the study of magnetic phase transitions of inhomogeneous systems, different ways of modeling the systems can be found in the literature. But to implement the inhomogeneity of real systems into a model usually make the model itself hardly to have exact solution. However, inhomogeneous systems are frequently encountered in nature, and they also play an important role in material science. For example, to develope new materials, the inhomogeneous systems formed by doping with other materials in a bulk system are often grown artificially in laboratory. Thence, it is important to construct some exact solvable models so that we can deepen our insights to the nature of such systems. The hierarchical spin models may be viewed as one of the efforts in this direction.

About one and half decades ago, spin models defined on hierarchical lattices received much attention in the literature$`\left[19\right]`$. Hierarchical lattices are defined as the infinite limit of iterative steps of replacing each fundamental bond or cell with a given type of decoration unit, and thence they are self-similar fractal lattices. Since hierarchical lattices are highly inhomogeneous in the local connectivity and lack of translational invariance, spin models defined on hierarchical lattices may provide a very good frame in understanding the magnetic properties of inhomogeneous systems. Moreover, we can apply the Migdal-Kadanoff renormalization$`[10,11]`$ scheme to obtain exact solutions of the partition functions of hierarchical spin models. Among a variety of hierarchical spin models, the zero-field $`q`$-state Potts model defined on diamond-hierarchical lattice has been widely studied. The results revealed from the exact solutions show that the corresponding systems have some nonclassical thermodynamic features. They include that the second-order phase transition occurs at finite temperature with the cusp-behavior appearing in the specific heat at the critical point$`[8,9,12]`$, the distribution of the partition function zeros owns a multifractal structure on the complex temperature plane, and the free energy near the critical point has a spatially modulated structure$`[6,7]`$.

We notice that owing to the infinite iterative steps of decorations in lattice construction, there exist a well-defined thermodynamic limit for a hierarchical spin model. Thence, the degree of the inhomogeneity is completely specified for a hierarchical spin model. If we had performed any finite iterative steps of decorations in lattice constrution, it would lead to the lack of thermodynamic limit for the model. To improve this while remain the model to be exactly solvable, we may study Ising systems defined on decorated lattices.

A decorated lattice is constructed by implementing bond- or cell- decorations to a regular lattice in an iterative way. The geometric complexity of a decorated lattice then is characterized by the decoration level $`n`$. For a decorated lattice of the decoration level $`n`$, a primary bond or cell defined on the original regular lattice is replaced by an $`n`$-bond or cell. An $`n`$-bond or cell owns a complicated inherent structure formed by $`n`$ iterative steps of replacing a bond or cell with some basic unit of decoration. In the limit of infinite $`n`$, an $`n`$-bond or cell essentially possesses a fractal lattice as the inherent structure. For this type of lattices, the degree of inhomogeneity in the local connectivity is proportional to the decoration level $`n`$ and there exists a well-defined thermodynamic limit for any value of $`n`$.

Triangular type decorated Ising models are examples. The analyses on these models have been carried out by Plechko$`\left[13\right]`$ and Plechko and Sobolev$`\left[14\right]`$, and there exists some non-classical results. In particular, the results of $`\left[14\right]`$ about the specific heat of the Ising model on cell-decorated triangular lattices have revealed some unusual features of ferromagnetic phase transitions. In the curve of the specific heat versus the temperature, as $`n`$ increases, the critical region is narrowing and a round peak appears above the critical point. While the critical point marks the set in for the long-range ordering among different $`n`$-cells, the appearence of the round peak corresponds to the occurrence of the short-range ordering inside an $`n`$-cell. As $`n`$ approaches to the limit of infinity, an $`n`$-cell has the Sierpiński gasket as its inherent structure. Then, in approaching the limit, $`n`$-cells are factorized and become independent in the Ising system. Moreover, the critical point approaches the zero temperature, and the round peak temperature remains finite.

Bond-decorated and cell-decorated lattices are different ways of modelling inhomogeneous systems. As the decoration level $`n`$ increases, the order of ramification remains a constant for a cell-decorated triangular type lattice but it increases unboundedly for a bond-decorated square lattice. Furthermore, it was pointed out by Gefen et. al. that there is no phase transition at finite temperature for a fractal lattice with the finite order of ramification$`\left[15\right]`$. Hence, the non-universal features of phase transitions in cell-decorated Ising models may be quite different from that of bond-decorations. To obtain a more complete picture about the possible features of inhomogeneous systems, we study the diamond-type bond-decorated Ising model. Note that the diamond-type bond-decoration was also adopted in a widely studied model, the diamond-hierarchical $`q`$-state Potts model$`\left[17\right]`$. In this work, we are interested in the following questions: $`\left(i\right)`$ How does the degree of the inhomogeneity in a system, characterized by the decoration level $`n`$, affect the nature of the phase transition? $`\left(ii\right)`$ How do the characteristic features of ferromagnetic phase transitions change with the lattice structure varying from cell-decorated triangular type lattices to bond-decorated square lattices? $`\left(iii\right)`$ How do the non-universal features obtained from the diamond-hierarchical Ising model$`\left[19\right]`$ appear? In particular, can we understand better about the cusp-behavior appearing in the specific heat at the critical point instead of the well-known logarithmic divergence?

The analytic expression of the partition function of square Ising model with diamond-type bond-decorations can be obtained by referring to the exact result of two-dimensional Ising model. The exact solution of the Ising model on a square lattice was first given by Onsager$`\left[16\right]`$. Since then, other alternatives in obtaining the exact solution have been developed. Besides of the combinatorial method, which was first developed by Kac and Ward$`\left[17\right]`$ and then rigorously reformulated by Green and Hurst$`\left[18\right]`$, more recently Plechko used a nonstandard method of Grassman multipliers to obtain the analytic expressions of the partition functions for a variety of cases$`[13,14,1921]`$. Based on the frame provided by these results, in Section 2 we use the bond-renormalization scheme to construct the analytic solution of the free energy. In Section 3, we determine the critical point of an $`n`$-lattice and study how the critical temperature vary with $`n`$. The condition of the critical point is given by the occurrence of the non-analyticity in the free energy. In Section 4, we use the solution of the free energy to calculate the internal energy and specific heat in an analytic way. Since the free energy can be expressed as the sum of those contributed by the local interactions among the sub-bonds of an $`n`$-bond and the long-range interactions among different $`n`$-bonds, we also use the same separation in the internal energy and specific heat. Thence, in general, we can obtain some insights about the nature of the phase transition of an inhomogeneous system by comparing the contributions from the local and the long-range interactions. Our results indicate that the cross over from a finite-decorated system to an infinite-decorated system may not be a smooth continuation, we then calculate the critical values of the internal energy and specific heat of an infinite-decorated system in Section 5. Finally, Section 6 is preserved for the summary and discussion.

## 2 free energy

We study the Ising model defined on an arbitrary $`n`$-lattice. A two-dimensional square lattice is used as backbones for the construction of an $`n`$-lattice. The simple square lattice is referred as $`0`$-lattice, and its connecting bonds between lattice sites are named as $`0`$-bonds. For a $`0`$-lattice, the total site-number is denoted as $`n_s`$ and the total bond-number as $`n_b`$ with $`n_b=2n_s`$. Starting from a $`0`$-bond, we obtain the corresponding $`n`$-bond as the result of replacing a bond with a diamond iteratively $`n`$ times. This $`n`$-bond then is specified by the two sites connected by the $`0`$-bond. For an $`n`$-bond, we have the total site-number $`s^{\left(n\right)}=2(4^n+2)/3`$ and bond number $`b^{\left(n\right)}=4^n`$. The construction of an $`n`$-bond is schematized in Fig. 1. From the primary structure of a $`0`$-lattice, we replace each $`0`$-bond with an $`n`$-bond to form an $`n`$-lattice. For an $`n`$-lattice, the average site- and bond-numbers per unit square are $`N_s^{\left(n\right)}=2s^{\left(n\right)}3`$ and $`N_b^{\left(n\right)}=2b^{\left(n\right)}`$ respectively. Examples of $`n`$-lattices are shown in Fig. 2.

The general form of the partition function for the Ising model defined on an $`n`$-lattice reads

$$Z^{\left(n\right)}=\underset{\left\{\sigma \right\}}{}\left[\underset{i,j}{}\mathrm{exp}\left(\eta \sigma _i\sigma _j\right)\right],$$
(1)
where the sum is over all bond-connected pairs $`i,j`$ of the $`n`$-lattice, and the Ising spin takes two possible values $`\sigma _i=\pm 1`$. Here we consider uniform ferromagnetic couplings characterized by the coupling strength $`J`$, and the dimensionless coupling parameter $`\eta `$ is defined as $`\eta =J/k_BT`$.

To implement the bond-renormalization scheme in the evaluation of the partition function, we rewrite Eq. (1) as

$$Z^{\left(n\right)}=2^{n_b\left(s^{\left(n\right)}2\right)}\underset{\left\{\sigma \right\}}{}\left[\underset{\mu ,\upsilon }{}B_{\mu ,\upsilon }^{\left(n\right)}\right],$$
(2)
where the Boltzmann $`n`$-factor, $`B_{\mu ,\upsilon }^{\left(n\right)}`$, is the contribution to the partition function from the Boltzmann factors associated with the $`n`$-bond specified by the two sites, $`\mu `$ and $`\upsilon `$, the sum is over the spins defined only on the lattice sites of the corresponding $`0`$-lattice, and the front factor is added to compensate the normalization factor added in $`B_{\mu ,\upsilon }^{\left(n\right)}`$. Note that the number of the decorated sites of an $`n`$-bond is $`s^{\left(n\right)}2`$, and the Ising spins on the decorated sites are referred as the inner-spins which couple only to the spins on the same $`n`$-bond. Thence, the $`B_{\mu ,\upsilon }^{\left(n\right)}`$ of Eq. (2) is defined as the result of taking the sum over the inner-spins for the product of all Boltzmann factors associated with the $`n`$-bond,

$$B_{\mu ,\upsilon }^{\left(n\right)}=\left(\frac{1}{2}\right)^{\left(s^{\left(n\right)}2\right)}\underset{\left\{\sigma _k,k=1,\mathrm{},s^{\left(n\right)}2\right\}}{}\mathrm{exp}\left[\eta \left(\sigma _\mu \sigma _1+\sigma _1\sigma _2+\mathrm{}+\sigma _{s^{\left(n\right)}2}\sigma _\nu \right)\right],$$
(3)
where the sum is over the $`\left(s^{\left(n\right)}2\right)`$ inner-spins, and the front factor is added for the normalization of the sum.

To obtain the explicit form of $`B_{\mu ,\upsilon }^{\left(n\right)}`$ for an arbitrary $`n`$, first we express the Boltzmann $`0`$-factor, which takes the form of $`\mathrm{exp}\left(\eta \sigma _\mu \sigma _\upsilon \right)`$, as

$$B_{\mu ,\upsilon }^{\left(0\right)}=\left(\mathrm{cosh}\eta \right)+\sigma _\mu \sigma _\upsilon \left(\mathrm{sinh}\eta \right).$$
(4)
According to the definition of Eq. (3), the Boltzmann $`1`$-factor is given as

$$B_{\mu ,\upsilon }^{\left(1\right)}=\left(\frac{1}{2}\right)^2\underset{\sigma _1,\sigma _2}{}\mathrm{exp}\left[\eta \left(\sigma _\mu \sigma _1+\sigma _\mu \sigma _2+\sigma _1\sigma _\nu +\sigma _2\sigma _\nu \right)\right].$$
(5)
Using Eq. (4) for $`\mathrm{exp}\left(\eta \sigma _\mu \sigma _\upsilon \right)`$, after some algebraic manipulations we obtain the result of Eq. (5) as

$$B_{\mu ,\upsilon }^{\left(1\right)}=\left[\mathrm{exp}(\eta ^{\left(1\right)})\right]\left[\mathrm{cosh}\left(\eta ^{\left(1\right)}\right)+\sigma _\mu \sigma _\upsilon \mathrm{sinh}\left(\eta ^{\left(1\right)}\right)\right],$$
(6)
with

$$\mathrm{exp}\left(\eta ^{\left(1\right)}\right)=\mathrm{cosh}\left(2\eta \right).$$
(7)
By continuing such construction, we can express Eq. (3) as

$$B_{\mu ,\upsilon }^{\left(n\right)}=R^{\left(n\right)}\left(\eta \right)\left[\mathrm{cosh}\left(\eta ^{\left(n\right)}\right)+\sigma _\mu \sigma _\upsilon \mathrm{sinh}\left(\eta ^{\left(n\right)}\right)\right]$$
(8)
for $`n1`$, where the function $`R^{\left(n\right)}\left(\eta \right)`$ is given by

$$R^{\left(n\right)}\left(\eta \right)=\underset{k=1}{\overset{n}{}}\left[\mathrm{exp}\left(\eta ^{\left(k\right)}\right)\right]^{4^{nk}},$$
(9)
and $`\eta ^{\left(k\right)}`$ is the image of $`\eta ^{\left(k1\right)}`$ under the renormalization map defined as

$$\mathrm{exp}\left(\eta ^{\left(k\right)}\right)=\mathrm{cosh}\left(2\eta ^{\left(k1\right)}\right)$$
(10)
with the initial condition, $`\eta ^{\left(0\right)}=\eta `$ for $`1kn`$.

For the Ising system defined on a single $`n`$-bond, the partition function of Eq. (2) has $`n_b=1`$ and it reduces to

$$Z_b^{\left(n\right)}=2^{s^{\left(n\right)}}\left(\frac{1}{2}\right)^2\underset{\sigma _\mu ,\sigma _\upsilon }{}B_{\mu ,\upsilon }^{\left(n\right)}.$$
(11)
After the substituition of Eq. (8) for $`B_{\mu ,\upsilon }^{\left(n\right)}`$, this yields

$$Z_b^{\left(n\right)}=2^{s^{\left(n\right)}}R^{\left(n\right)}\left(\eta \right)\mathrm{cosh}\left(\eta ^{\left(n\right)}\right).$$
(12)
Then, the corresponding free energy per bond per $`k_BT`$, is given as

$$f_b^{\left(n\right)}=\frac{1}{b^{\left(n\right)}}\mathrm{ln}Z_b^{\left(n\right)},$$
(13)
where $`b^{\left(n\right)}`$ is the total bond number contained in an $`n`$-bond. This, up to a constant, yields

$$f_b^{\left(n\right)}=\frac{1}{4^n}\mathrm{ln}\left[\mathrm{cosh}\left(\eta ^{\left(n\right)}\right)\right]\underset{k=1}{\overset{n}{}}\frac{\eta ^{\left(k\right)}}{4^k}.$$
(14)
Note that $`f_b^{\left(n\right)}`$ in the limit of infinite $`n`$ is exactly the free energy density of the diamond-hierarchical Ising model.

Using the result of Eq. (12) for $`Z_b^{\left(n\right)}`$, we can express Eq. (2) in terms of $`Z_b^{\left(n\right)}`$ as

$$Z^{\left(n\right)}=2^{n_s}\left(Z_b^{\left(n\right)}\right)^{n_b}Q^{\left(n\right)},$$
(15)
where $`Q^{\left(n\right)}`$, referred as the reduced partition function, takes the form of

$$Q^{\left(n\right)}=(\frac{1}{2})^{n_s}\underset{\left\{\sigma \right\}}{}\left\{\underset{\mu ,\upsilon }{}\left[1+\sigma _\mu \sigma _\upsilon \mathrm{tanh}\eta ^{\left(n\right)}\right]\right\}.$$
(16)
The reduced partition function of an $`n`$-lattice given by Eq. (16) has exactly the same form as that of simple square lattice when the variable $`\eta ^{\left(n\right)}`$ is changed to $`\eta `$. Thence, the corresponding free energy can be easily written down according to the formal expression of exact solution of simple square lattice provided by Refs. $`[13,16,19,21,22]`$. Up to a constant, we can express the free energy per bond per $`k_BT`$ as

$$f^{\left(n\right)}=f_b^{\left(n\right)}+f_{bb}^{\left(n\right)},$$
(17)
where $`f_b^{\left(n\right)}`$ is given by Eq. (14), and $`f_{bb}^{\left(n\right)}`$ is the contribution from $`Q^{\left(n\right)}`$ and it is given as

$$f_{bb}^{\left(n\right)}=\frac{1}{2N_b^{\left(n\right)}}_0^{2\pi }\frac{d\varphi }{2\pi }_0^{2\pi }\frac{d\theta }{2\pi }\mathrm{ln}\left[A_0^{(n)}A_1^{(n)}\left(\mathrm{cos}\theta +\mathrm{cos}\varphi \right)\right],$$
(18)
with $`A_0^{(n)}`$ and $`A_1^{(n)}`$ defined as

$$A_0^{(n)}=\frac{\mathrm{cosh}^2\left(2\eta ^{\left(n\right)}\right)}{\mathrm{cosh}^4\left(\eta ^{\left(n\right)}\right)},$$
(19)
$$A_1^{(n)}=\frac{2\mathrm{sinh}\left(\eta ^{\left(n\right)}\right)}{\mathrm{cosh}^3\left(\eta ^{\left(n\right)}\right)}.$$
(20)
Since $`f_b^{\left(n\right)}`$ is the free energy density of independent $`n`$-bonds, $`f_{bb}^{\left(n\right)}`$ can be identified as the contribution from the interactions among different $`n`$-bonds. Thus, we refer $`f_b^{\left(n\right)}`$ as the contribution from the local interactions and $`f_{bb}^{\left(n\right)}`$ as that from the long-range interactions.

To simplify the expression of $`f_{bb}^{\left(n\right)}`$, we change the variables, $`\omega _1=\left(\theta +\varphi \right)/2`$ and $`\omega _2=\left(\theta \varphi \right)/2`$, to rewrite Eq. (18) as

$$f_{bb}^{\left(n\right)}=\frac{1}{2N_b^{\left(n\right)}}\mathrm{ln}A_0^{(n)}\frac{1}{N_b^{\left(n\right)}}_0^{\pi /2}\frac{d\omega _2}{\pi }_0^\pi \frac{d\omega _1}{\pi }\mathrm{ln}\left(1\kappa ^{\left(n\right)}\mathrm{cos}\omega _1\mathrm{cos}\omega _2\right),$$
(21)
with

$$\kappa ^{\left(n\right)}=\frac{2A_1^{(n)}}{A_0^{(n)}}.$$
(22)
Then, the integral over $`\omega _1`$ can be easily performed and the result is

$$f_{bb}^{\left(n\right)}=\frac{1}{2N_b^{\left(n\right)}}\mathrm{ln}A_0^{(n)}\frac{1}{N_b^{\left(n\right)}}_0^{\pi /2}\frac{d\omega }{\pi }\mathrm{ln}\left[\frac{1}{2}\left(1+\sqrt{1\kappa ^{\left(n\right)2}\mathrm{sin}^2\omega }\right)\right].$$
(23)
## 3 critical point

The free energy density, $`f_{bb}^{\left(n\right)}`$, is not completely analytic for the physical region of temperature. The non-analyticity arises from the second derivative of $`f_{bb}^{\left(n\right)}`$ with respect to temperature at the point $`\kappa ^{\left(n\right)}=1`$ for the physical region $`0\kappa ^{\left(n\right)}1`$. This implies that for the Ising model on an $`n`$-lattice, the long-range interactions among different $`n`$-bonds are responsible for the occurrence of the Ising phase transition at the finite temperature, and the critical point is determined by the condition $`\kappa ^{\left(n\right)}=1`$. Thus, we can obtain the bulk critical temperature of the ferromagnetic phase transition from the solution of the condition,

$$A_0^{(n)}2A_1^{(n)}\stackrel{𝑐}{=}0.$$
(24)
Here, for convenience, we use the notation, $`\stackrel{𝑐}{=}`$, to denote the equivalence established only at the critical temperature. With the definitions of $`A_0^{(n)}`$ and $`A_1^{(n)}`$ given by Eqs. (19) and (20), we can rewrite the critical condition as the familiar form, $`\mathrm{sinh}2\eta ^{\left(n\right)}\stackrel{𝑐}{=}1`$, or as

$$\mathrm{exp}\left(\eta ^{\left(n\right)}\right)\stackrel{𝑐}{=}h^{(0)}$$
(25)
with

$$h^{(0)}=\sqrt{1+\sqrt{2}}$$
(26)
for an $`n`$-lattice.

As the consequence of the map of Eq. (10), Eq. (25) yields the critical value of the $`\mathrm{exp}\left(\eta ^{\left(n1\right)}\right)`$ variable as

$$\mathrm{exp}\left(\eta ^{\left(n1\right)}\right)\stackrel{𝑐}{=}h^{(1)},$$
(27)
for $`n1`$ with

$$h^{(1)}=\left(h^{\left(0\right)}+\sqrt{\left(h^{(0)}\right)^21}\right)^{1/2}.$$
(28)
By applying the recursion relation of Eq. (10) continuously, we can establish the result,

$$\mathrm{exp}\left(\eta ^{\left(nk\right)}\right)\stackrel{𝑐}{=}h^{(k)},$$
(29)
for $`1kn`$, where $`h^{\left(k\right)}`$ is determined by the recursion relation

$$h^{(k)}=\left(h^{\left(k1\right)}+\sqrt{\left(h^{(k1)}\right)^21}\right)^{1/2}.$$
(30)
By taking $`k=n`$ in Eq. (29) we obtain the reduced critical temperature of the Ising transition as

$$\left(\frac{k_BT_c}{J}\right)_n=\left[\mathrm{ln}\left(h^{\left(n\right)}\right)\right]^1$$
(31)
for an arbitrary $`n`$-lattice. For $`n=0`$, this gives the well-known result, $`k_BT_c/J=2.269185`$. As $`n`$ approaches $`\mathrm{}`$, the $`h^{\left(n\right)}`$ value can be determined by the fixed points of the map of Eq. (30). In fact, the $`h^{\left(\mathrm{}\right)}`$ value is the repellor of the map of Eq. (30), and this repellor locates at $`h^{\left(\mathrm{}\right)}=1.839287`$. The corresponding critical temperature is $`\left(k_BT_c/J\right)_n\mathrm{}=1.641018`$. However, the nature of the phase transition at the repellor remains to be answered, and this will be discussed in Section 5 where the case of infinite decoration-level is discussed.

The numerical values of $`\left(k_BT_c/J\right)_n`$ for some different $`n`$ values are given in Table 1, and the result as a function of $`n`$ is shown in Fig. 3 with the solid curve given by the fitting function $`\left(k_BT_c/J\right)_n=1.6410+\left(0.6281\right)\mathrm{exp}\left[\left(0.5857\right)n\right]`$. This result indicates that the critical temperature decreases exponentially toward the limiting value $`1.6410`$ as $`n`$ increases. On the other hand, it has been shown$`\left[14\right]`$ that for the cell-decorated triangular Ising model the critical temperature as a function of $`n`$ is given by $`\left(k_BT_c/J\right)_n=4/\mathrm{ln}\left\{4n4\mathrm{ln}\left[\left(n+1\right)/2\right]+1\right\}`$ for $`n1`$, and this leads to $`\left(k_BT_c/J\right)_n=4/\mathrm{ln}\left(4n\right)`$ as $`n\mathrm{}`$. Thus, the decreasing rate of the critical temperature toward the limit of infinite decoration level is more faster for the bond-decorated Ising model in comparing with the cell-decorated Ising model.

## 4 specific heat

In this section, we follow the standard procedure, for example as given by Ref. $`\left[22\right]`$, to calculate the specific heat for an $`n`$-lattice. The analytic result is expressed in terms of the complete elliptic integrals of the first and second kinds. The contributions to the specific heat contain those from the local and the long-range interactions, and the comparison between these two is given quantitatively. The critical behaviors of the specific heats of different $`n`$-lattices are also discussed. In particular, we compare the results with the cell-decorated case.

Using the free energy density $`f^{\left(n\right)}`$ of an $`n`$-lattice given by Eq. (17), we define the dimensionless internal energy per bond as

$$ϵ^{\left(n\right)}=\frac{f^{\left(n\right)}}{\eta }.$$
(32)
Due to the separable contributions to the free energy density, we can write the internal energy density as the sum of two parts,

$$ϵ^{\left(n\right)}=ϵ_b^{\left(n\right)}+ϵ_{bb}^{\left(n\right)},$$
(33)
with $`ϵ_b^{\left(n\right)}=f_b^{\left(n\right)}/\eta `$ and $`ϵ_{bb}^{\left(n\right)}=f_{bb}^{\left(n\right)}/\eta `$. By using $`f_b^{\left(n\right)}`$ and $`f_{bb}^{\left(n\right)}`$ of Eqs. (14) and (23), we can perform the direct differentiation to obtain the corresponding internal energy density as

$$ϵ_b^{\left(n\right)}=\left(\frac{\mathrm{tanh}\left(\eta ^{\left(n\right)}\right)}{4^n}\right)\left(\frac{\eta ^{\left(n\right)}}{\eta }\right)\underset{k=1}{\overset{n}{}}\left[\left(\frac{1}{4^k}\right)\left(\frac{\eta ^{\left(k\right)}}{\eta }\right)\right],$$
(34)
and

$$ϵ_{bb}^{\left(n\right)}=\frac{1}{2N_b^{\left(n\right)}}\left[\left(\frac{}{\eta }\mathrm{ln}A_1^{(n)}\right)\left(\frac{2}{\pi }\right)\left(\frac{}{\eta }\mathrm{ln}\kappa ^{\left(n\right)}\right)K\left(\kappa ^{\left(n\right)}\right)\right],$$
(35)
where the explicit forms of the functions, $`\mathrm{ln}A_1^{(n)}/\eta `$ and $`\mathrm{ln}\kappa ^{\left(n\right)}/\eta `$, can be obtained by direct differentiations of Eqs. (20) and (22),

$$\frac{}{\eta }\mathrm{ln}A_1^{(n)}=\left(\frac{\mathrm{cosh}\left(\eta ^{\left(n\right)}\right)}{\mathrm{sinh}\left(\eta ^{\left(n\right)}\right)}\frac{3\mathrm{sinh}\left(\eta ^{\left(n\right)}\right)}{\mathrm{cosh}\left(\eta ^{\left(n\right)}\right)}\right)\left(\frac{\eta ^{\left(n\right)}}{\eta }\right),$$
(36)
and

$$\frac{}{\eta }\mathrm{ln}\kappa ^{\left(n\right)}=\left(\frac{2\mathrm{cosh}\left(2\eta ^{\left(n\right)}\right)}{\mathrm{sinh}\left(2\eta ^{\left(n\right)}\right)}\frac{4\mathrm{sinh}\left(2\eta ^{\left(n\right)}\right)}{\mathrm{cosh}\left(2\eta ^{\left(n\right)}\right)}\right)\left(\frac{\eta ^{\left(n\right)}}{\eta }\right),$$
(37)
with

$$\frac{\eta ^{\left(k\right)}}{\eta }=2^k\left\{\underset{j=1}{\overset{k}{}}\left[\mathrm{tanh}\left(2\eta ^{\left(j1\right)}\right)\right]\right\}$$
(38)
for $`1kn`$ obtained from the map of Eq. (10), and $`K\left(\kappa ^{\left(n\right)}\right)`$ is the complete elliptic integral of the first kind defined as

$$K\left(\kappa ^{\left(n\right)}\right)=_0^{\pi /2}\frac{d\omega }{\left(1\kappa ^{\left(n\right)2}\mathrm{sin}^2\omega \right)^{1/2}},$$
(39)
We notice that $`K\left(\kappa ^{\left(n\right)}\right)`$ is singular at the critical point $`\kappa ^{\left(n\right)}=1`$. As $`\kappa ^{\left(n\right)}1`$, it can be shown$`\left[22\right]`$ that

$$K\left(\kappa ^{\left(n\right)}\right)\underset{TT_c}{=}\mathrm{ln}\left|\frac{TT_c}{T_c}\right|.$$
(40)
This logarithmic divergence appears in the specific heat, and it becomes the signature of the Ising phase transition in two-dimensions. But, by substituting the critical condition, $`\mathrm{sinh}2\eta ^{\left(n\right)}\stackrel{𝑐}{=}1`$, into Eq. (37), we obtain

$$\frac{}{\eta }\mathrm{ln}\kappa ^{\left(n\right)}\stackrel{𝑐}{=}0.$$
(41)
As a consequence of Eq. (41), we see that $`ϵ_{bb}^{\left(n\right)}`$ of Eq. (35) is a continuous function of temperature $`T`$ even at the critical point, $`T=T_c`$.

To have a qualitative understanding about the varitions of $`ϵ_b^{\left(n\right)}`$ and $`ϵ_{bb}^{\left(n\right)}`$ with respect to the decoration-level $`n`$, we compute the critical values of $`ϵ_b^{\left(n\right)}`$ and $`ϵ_{bb}^{\left(n\right)}`$ for different $`n`$ values. By using the critical condition of Eq. (29) for an $`n`$-lattice, we obtain the critical values of Eqs. (36) and (38) as

$$\frac{}{\eta }\mathrm{ln}A_1^{(n)}\stackrel{𝑐}{=}\left(42\sqrt{2}\right)\left(\frac{\eta ^{\left(n\right)}}{\eta }\right)_c,$$
(42)
and

$$\frac{\eta ^{\left(k\right)}}{\eta }\stackrel{𝑐}{=}2^k\left\{\underset{j=nk+1}{\overset{n}{}}\left[\frac{\left(h^{\left(j\right)}\right)^41}{\left(h^{\left(j\right)}\right)^4+1}\right]\right\}$$
(43)
for $`1kn`$. Here and hereafter, we use the subscript $`c`$ to denote that the quantity in the parenthesis is evaluated at the critical point of an $`n`$-lattice. Then, using Eqs. (41) and (42) we can express the critical value of $`ϵ_b^{\left(n\right)}`$ and $`ϵ_{bb}^{\left(n\right)}`$ given by Eqs. (34) and (35) as

$$ϵ_b^{\left(n\right)}\stackrel{𝑐}{=}\frac{1}{4^n}\left(\frac{\sqrt{2}}{2+\sqrt{2}}\right)\left(\frac{\eta ^{\left(n\right)}}{\eta }\right)_c\underset{k=1}{\overset{n}{}}\left[\left(\frac{1}{4^k}\right)\left(\frac{\eta ^{\left(k\right)}}{\eta }\right)_c\right],$$
(44)
and

$$ϵ_{bb}^{\left(n\right)}\stackrel{𝑐}{=}\frac{1}{4^{n+1}}\left(42\sqrt{2}\right)\left(\frac{\eta ^{\left(n\right)}}{\eta }\right)_c,$$
(45)
with $`\left(\eta ^{\left(k\right)}/\eta \right)_c`$ given by Eq. (43).

The numerical results of Eqs. (44) and (45) for some decoration levels $`n`$ are given in Table 1, and the corresponding curves are shown in Fig. 4. These results indicate that the local correlations among the $`b^{\left(n\right)}`$ bonds of a given $`n`$-bond are enhenced and the long-range correlations among different $`n`$-bonds are reduced as $`n`$ increases. For $`n6`$, the long-range correlations are almost absent, and the system is very closed to that composing of independent $`n`$-bonds.

In calculating the specific heat $`C_V^{\left(n\right)}`$ of an $`n`$-lattice, we use the dimensionless quantity, $`c^{\left(n\right)}=C_V^{\left(n\right)}/k_B`$, defined as

$$c^{(n)}=\eta ^2\frac{ϵ^{\left(n\right)}}{\eta }.$$
(46)
Similar to the internal energy density, we can express $`c^{(n)}`$ as the sum of two parts,

$$c^{\left(n\right)}=c_b^{\left(n\right)}+c_{bb}^{\left(n\right)},$$
(47)
with $`c_b^{\left(n\right)}=\eta ^2ϵ_b^{\left(n\right)}/\eta `$ and $`c_b^{\left(n\right)}=\eta ^2ϵ_b^{\left(n\right)}/\eta `$.

For the direct differentiation of $`ϵ_b^{\left(n\right)}`$ of Eq. (34) with respect to $`\eta `$, we obtain

$$\frac{c_b^{\left(n\right)}}{\eta ^2}=\frac{1}{4^n}\left[\left(\frac{\eta ^{\left(n\right)}/\eta }{\mathrm{cosh}\left(\eta ^{\left(n\right)}\right)}\right)^2+\left(\mathrm{tanh}\left(\eta ^{\left(n\right)}\right)\right)\left(\frac{^2\eta ^{\left(n\right)}}{\eta ^2}\right)\right]+\left\{\underset{k=1}{\overset{n}{}}\left[\left(\frac{1}{4^k}\right)\left(\frac{^2\eta ^{\left(k\right)}}{\eta ^2}\right)\right]\right\},$$
(48)
where the function, $`\eta ^{\left(k\right)}/\eta `$, is given by Eq. (38), and its further derivative yields

$$\frac{^2\eta ^{\left(k\right)}}{\eta ^2}=\left(\frac{\eta ^{\left(k\right)}}{\eta }\right)\left\{\underset{i=1}{\overset{k}{}}\left[\left(\frac{2^{i+1}}{\mathrm{sinh}\left(4\eta ^{\left(i1\right)}\right)}\right)\left(\underset{j=1}{\overset{i1}{}}\mathrm{tanh}\left(2\eta ^{\left(j1\right)}\right)\right)\right]\right\}$$
(49)
for $`k1`$. Note that in Eq. (49) the product for the index $`j`$ in the last parenthesis exists only for $`i2`$.

On the other hand, the direct differentiation of $`ϵ_{bb}^{\left(n\right)}`$ of Eq. (35) with respect to $`\eta `$ yields

$$\frac{c_{bb}^{\left(n\right)}}{\eta ^2}=\frac{1}{2N_b^{\left(n\right)}}\left\{\left(\frac{^2}{\eta ^2}\mathrm{ln}A_1^{\left(n\right)}\right)\frac{2}{\pi }\left[P_0(\kappa ^{\left(n\right)})K\left(\kappa ^{\left(n\right)}\right)+Q_0\left(\kappa ^{\left(n\right)}\right)E\left(\kappa ^{\left(n\right)}\right)\right]\right\},$$
(50)
where the functions, $`P_0\left(\kappa ^{\left(n\right)}\right)`$ and $`Q_0\left(\kappa ^{\left(n\right)}\right)`$, are defined as

$$P_0\left(\kappa ^{\left(n\right)}\right)=\left(\frac{^2}{\eta ^2}\mathrm{ln}\kappa ^{\left(n\right)}\right)\left(\frac{}{\eta }\mathrm{ln}\kappa ^{\left(n\right)}\right)^2$$
(51)
and

$$Q_0\left(\kappa ^{\left(n\right)}\right)=\left(\frac{1}{1\left(\kappa ^{\left(n\right)}\right)^2}\right)\left(\frac{}{\eta }\mathrm{ln}\kappa ^{\left(n\right)}\right)^2,$$
(52)
and $`E\left(\kappa ^{\left(n\right)}\right)`$ is referred as the complete elliptic integral of the second kind and defined as

$$E\left(\kappa ^{\left(n\right)}\right)=_0^{\pi /2}𝑑\omega \left(1\kappa ^{\left(n\right)2}\mathrm{sin}^2\omega \right)^{1/2}.$$
(53)
Note that in obtaining Eq. (50) we have used the well-known formula$`\left[22\right]`$,

$$\frac{K\left(\kappa ^{\left(n\right)}\right)}{\kappa ^{\left(n\right)}}=\left(\frac{1}{\kappa ^{\left(n\right)}}\right)\left[\left(\frac{1}{1\left(\kappa ^{\left(n\right)}\right)^2}\right)E\left(\kappa ^{\left(n\right)}\right)K\left(\kappa ^{\left(n\right)}\right)\right].$$
(54)
We also notice that by further differentiating the results of Eqs. (36) and (37) with respect to $`\eta `$ we obtain

$$\frac{^2}{\eta ^2}\mathrm{ln}A_1^{\left(n\right)}=P_1\left(\eta ^{\left(n\right)}\right)\left(\frac{\eta ^{\left(n\right)}}{\eta }\right)^2+Q_1\left(\eta ^{\left(n\right)}\right)\left(\frac{^2\eta ^{\left(n\right)}}{\eta ^2}\right),$$
(55)
and

$$\frac{^2}{\eta ^2}\mathrm{ln}\kappa ^{\left(n\right)}=P_2\left(\eta ^{\left(n\right)}\right)\left(\frac{\eta ^{\left(n\right)}}{\eta }\right)^2+Q_2\left(\eta ^{\left(n\right)}\right)\left(\frac{^2\eta ^{\left(n\right)}}{\eta ^2}\right),$$
(56)
with

$$P_1\left(\eta ^{\left(n\right)}\right)=\frac{1}{\mathrm{sinh}^2\left(\eta ^{\left(n\right)}\right)}+\frac{3}{\mathrm{cosh}^2\left(\eta ^{\left(n\right)}\right)},$$
(57)
$$Q_1\left(\eta ^{\left(n\right)}\right)=\frac{\mathrm{cosh}\left(\eta ^{\left(n\right)}\right)}{\mathrm{sinh}\left(\eta ^{\left(n\right)}\right)}\frac{3\mathrm{sinh}\left(\eta ^{\left(n\right)}\right)}{\mathrm{cosh}\left(\eta ^{\left(n\right)}\right)},$$
(58)
$$P_2\left(\eta ^{\left(n\right)}\right)=\frac{4}{\mathrm{sinh}^2\left(2\eta ^{\left(n\right)}\right)}+\frac{8}{\mathrm{cosh}^2\left(2\eta ^{\left(n\right)}\right)},$$
(59)
and

$$Q_2\left(\eta ^{\left(n\right)}\right)=\frac{2\mathrm{cosh}\left(2\eta ^{\left(n\right)}\right)}{\mathrm{sinh}\left(2\eta ^{\left(n\right)}\right)}\frac{4\mathrm{sinh}\left(2\eta ^{\left(n\right)}\right)}{\mathrm{cosh}\left(2\eta ^{\left(n\right)}\right)}.$$
(60)
The valuess of $`c_b^{\left(n\right)}`$ and $`c_{bb}^{\left(n\right)}`$ of Eqs. (48) and (50) can be calculated by using the explicit forms of Eqs. (36-38), (49), (55), and (56). The numerical results as a function of reduced temperature $`k_BT/J`$ for some $`n`$ values are shown in Fig. 5. The peak temperature of $`c_b^{\left(n\right)}`$, denoted by $`\left(k_BT_m/J\right)_n`$, marks the occurrence of the local ordering in an $`n`$-bond, and the numerical values of $`\left(k_BT_m/J\right)_n`$ for some decoration levels $`n`$ are given in Table 1. The plot of $`\left(k_BT_m/J\right)_n`$ versus $`n`$ is also shown in the lower part of Fig. 3 with the solid curve given by the fitting function $`\left(k_BT_m/J\right)_n=1.6410\left(0.8063\right)\mathrm{exp}\left[\left(0.7144\right)n\right]`$. Thus, because of the enhancement of the local correlation for the increase of $`n`$, $`\left(k_BT_m/J\right)_n`$ increases as $`n`$ increases. Moreover, the peak temperature of $`c_{bb}^{\left(n\right)}`$ is the critical temperature $`\left(k_BT_c/J\right)_n`$ which signifies the occurrence of the long-range ordering among the Ising spins defined on the primary sites, and the behavior of $`\left(k_BT_c/J\right)_n`$ versus $`n`$ is already given in the upper part of Fig. 3. The results of Fig. 3 indicate that the local ordering occurs at the temperature below the critical temperature, and this is exactly opposite to the case of cell-decorated triangular Ising model$`\left[14\right]`$.

The singular behavior of $`c^{\left(n\right)}`$ arises from $`c_{bb}^{\left(n\right)}`$ of Eq. (50) in the form of the logarithmic divergence given by Eq. (40). Thus, we can express the singular behavior of $`c^{\left(n\right)}`$ as

$$c^{\left(n\right)}\underset{TT_c}{=}A_{\mathrm{sin}g}^{\left(n\right)}\mathrm{ln}\left|\frac{TT_c}{T_c}\right|.$$
(61)
where $`A_{\mathrm{sin}g}^{\left(n\right)}`$, referred as the critical amplitude, is given as

$$A_{\mathrm{sin}g}^{\left(n\right)}=\frac{\eta _c^2}{N_b^{\left(n\right)}\pi }\left[P_0\left(\kappa ^{\left(n\right)}\right)\right]_c,$$
(62)
with $`P_0\left(\kappa ^{\left(n\right)}\right)`$ given by Eq. (51). By using Eq. (41), we can rewrite the expression of $`A_{\mathrm{sin}g}^{\left(n\right)}`$ as

$$A_{\mathrm{sin}g}^{\left(n\right)}=\frac{\eta _c^2}{N_b^{\left(n\right)}\pi }\left(\frac{^2}{\eta ^2}\mathrm{ln}\kappa ^{\left(n\right)}\right)_c.$$
(63)
By taking the critical value of Eq. (56) with $`\left[Q_2\left(\eta ^{\left(n\right)}\right)\right]_c=0`$ and $`\left[P_2\left(\eta ^{\left(n\right)}\right)\right]_c=8`$, we can use the result of Eq. (43) to obtain

$$\frac{^2}{\eta ^2}\mathrm{ln}\kappa ^{\left(n\right)}\stackrel{𝑐}{=}\frac{4^{n+2}}{2}\left\{\underset{k=1}{\overset{n}{}}\left[\frac{\left(h^{\left(k\right)}\right)^41}{\left(h^{\left(k\right)}\right)^4+1}\right]\right\}^2.$$
(64)
Thus, using the reduced critical temperature given by Eq. (31) we can express Eq. (63) as

$$A_{\mathrm{sin}g}^{\left(n\right)}=\frac{4}{\pi }\left[\mathrm{ln}h^{\left(n\right)}\right]^2\left\{\underset{k=1}{\overset{n}{}}\left[\frac{\left(h^{\left(k\right)}\right)^41}{\left(h^{\left(k\right)}\right)^4+1}\right]\right\}^2.$$
(65)
The numerical values of Eq. (65) for some $`n`$ values are given in Table 1. The plot of $`A_{\mathrm{sin}g}^{\left(n\right)}`$ versus $`n`$ is also given in Fig. 6 where the solid curve is given by the fitting function $`A_{\mathrm{sin}g}^{\left(n\right)}=\left(0.2473\right)\mathrm{exp}\left[\left(0.3018\right)n\right]`$. Since the critical amplitude associated with the logrithmic singularity characterizes the width of the critical region, our results indicate that as increasing $`n`$ the critical region is narrowing down to the vanish limit. This characteristic feature in the specific heat can also be seen in the model on anisotropic regular lattices$`[13,14,20,22]`$. We also notice that for the cell-decorated triangular Ising model the critical amplitude behaves as $`A_{\mathrm{sin}g}^{\left(n\right)}3^n\left(n\mathrm{ln}n\right)^2`$ for very large $`n\left[14\right]`$. Thus, as the decoration level $`n`$ increases the critical amplitude approaches zero less quickly for the bond-decorated model in comparing with the cell-decorated model.

## 5 infinite decoration-level

In this section, we discuss the critical properties of the Ising model on an $`n`$-lattice in the limit of infinite $`n`$. First we demonstrate the absence of Ising phase transition. This feature is caused mainly by the fragmentation of $`n`$-bonds. Then, we calculate the critical values of the internal energy and specific heat. The specific heat curve near the critical point becomes a cusp, and the exact critical value is given.

The map of Eq. (10) has three fixed points, $`\eta _f`$, given by the solutions of

$$\eta _f=\mathrm{ln}\left[\mathrm{cosh}\left(2\eta _f\right)\right].$$
(66)
Among the three fixed points, two are attractors locating at $`\eta _f^{a_1}=0`$ and $`\eta _f^{a_2}=\mathrm{}`$ respectively, and one is a repellor locating at $`\eta _f^r=0.609378`$. For the limit of infinite decoration-level, the repellor is the critical point of the bulk system. Since the critical point is a repellor, the critical behavior of the system may be different from that obtained by extrapolating the results of finite-decorated systems to the limit of infinite $`n`$. To show this, first we give the plots of $`\eta ^{\left(n\right)}`$ versus $`\eta `$ for $`n=0`$, $`2`$, $`4`$, and $`6`$ in Fig. 7. The results of Fig. 7 indicate that for the limit of infinite $`n`$ the possible $`\eta ^{\left(n\right)}`$ values are restricted to the fixed points of Eq. (66): $`\eta ^{\left(n\right)}=\eta _f^{a_1}`$ for $`\eta <\eta _f^r`$, $`\eta ^{\left(n\right)}=\eta _f^r`$ for $`\eta =\eta _f^r`$, and $`\eta ^{\left(n\right)}=\eta _f^{a_2}`$ for $`\eta >\eta _f^r`$. This leads to the result that the value of $`\kappa ^{\left(n\right)}`$ defined by Eq. (22) has only two possible values,

$$\kappa ^{\left(\mathrm{}\right)}=0$$
(67)
for $`\eta \eta _f^r`$, and

$$\kappa ^{\left(\mathrm{}\right)}=2\sqrt{\left[\mathrm{exp}(2\eta _f^r)\right]\left[\mathrm{exp}\left(4\eta _f^r\right)\right]}$$
(68)
for $`\eta =\eta _f^r`$. The latter gives the numerical value $`\kappa ^{\left(\mathrm{}\right)}=0.912622`$. Thence, the logarithmic divergence of the Ising transition, caused by the divergence of the elliptic integral $`K\left(\kappa ^{\left(n\right)}\right)`$ of Eq. (39) at the point $`\kappa ^{\left(n\right)}=1`$, is absent for the infinite $`n`$ limit. This indicates that the cross over from a finite-decorated system to an infinite-decorated system may not be a smooth continuation.

To exhibit the critical properties of the infinite-decorated system more explicitly, we then calculate the critical values of the internal energy and specific heat. Since the critical point of the infinite-decorated system is the repellor of the map of Eq. (10), we have the following identities for the derivatives evaluated at the critical point $`\eta =\eta _f^r`$:

$$\frac{\eta ^{\left(i+k\right)}}{\eta ^{\left(i\right)}}\stackrel{}{=}\frac{\eta ^{\left(k\right)}}{\eta },$$
(69)
and

$$\frac{^2\eta ^{\left(i+k\right)}}{\eta ^{\left(i\right)2}}\stackrel{}{=}\frac{^2\eta ^{\left(k\right)}}{\eta ^2}.$$
(70)
Here, we use the notation, $`\stackrel{}{=}`$ , to denote the equivalence established at the repellor. On the other hand, by successively applying the chain rule we have

$$\frac{\eta ^{\left(k\right)}}{\eta }=\underset{i=1}{\overset{k}{}}\frac{\eta ^{\left(i\right)}}{\eta ^{\left(i1\right)}}.$$
(71)
Then, as a result of Eq. (69) we obtain

$$\frac{\eta ^{\left(k\right)}}{\eta }\stackrel{}{=}\left(d^{\left(1\right)}\right)^k,$$
(72)
with $`d^{\left(1\right)}=\left(\eta ^{\left(1\right)}/\eta \right)_{\eta =\eta _f^r}=1.678574`$.

For the second derivative, $`^2\eta ^{\left(k\right)}/\eta ^2`$, we first perform the direct differentiation of Eq. (71) and then use the chain rule to obtain

$$\frac{^2\eta ^{\left(k\right)}}{\eta ^2}=\left(\frac{^2\eta ^{\left(1\right)}}{\eta ^2}\right)\left[D_2^k\left(\eta \right)\right]+\left(\frac{\eta ^{\left(1\right)}}{\eta }\right)^2\left(\frac{^2\eta ^{\left(2\right)}}{\eta ^{\left(1\right)2}}\right)\left[D_3^k\left(\eta \right)\right]+\mathrm{}+\left[D_1^{k1}\left(\eta \right)\right]^2\left(\frac{^2\eta ^{\left(k\right)}}{\eta ^{\left(k1\right)2}}\right),$$
(73)
with

$$D_a^k\left(\eta \right)=\underset{i=a}{\overset{k}{}}\frac{\eta ^{\left(i\right)}}{\eta ^{\left(i1\right)}}$$
(74)
for $`ak`$. Then, by using the identities of Eqs. (69) and (70) we can express Eq. (73) as

$$\frac{^2\eta ^{\left(k\right)}}{\eta ^2}\stackrel{}{=}d^{\left(2\right)}\left[\left(d^{\left(1\right)}\right)^{k1}+\left(d^{\left(1\right)}\right)^k+\mathrm{}+\left(d^{\left(1\right)}\right)^{2\left(k1\right)}\right]$$
(75)
which yields

$$\frac{^2\eta ^{\left(k\right)}}{\eta ^2}\stackrel{}{=}d^{\left(2\right)}\left(d^{\left(1\right)}\right)^{k1}\left(\frac{\left(d^{\left(1\right)}\right)^k1}{d^{\left(1\right)}1}\right),$$
(76)
with $`d^{\left(2\right)}=\left(^2\eta ^{\left(1\right)}/\eta ^2\right)_{\eta =\eta _f^r}=1.182391`$.

To calculate the internal energy density and specific heat at the repellor, we take $`\eta =\eta _f^r`$, Eq. (72) for $`\left(\eta ^{\left(n\right)}/\eta \right)_{\eta =\eta _f^r}`$, and Eq. (76) for $`\left(^2\eta ^{\left(n\right)}/\eta ^2\right)_{\eta =\eta _f^r}`$. Then, for the internal energy density of Eqs. (34) and (35) we obtain the numerical results as

$$ϵ_b^{\left(n\right)}\stackrel{}{=}\left(0.543689\right)\left(a^{\left(1\right)}\right)^na^{\left(1\right)}\left(\frac{1\left(a^{\left(1\right)}\right)^n}{1a^{\left(1\right)}}\right),$$
(77)
and

$$ϵ_{bb}^{\left(n\right)}\stackrel{}{=}\left(0.414742\right)\left(a^{\left(1\right)}\right)^n,$$
(78)
with $`a^{\left(1\right)}=d^{\left(1\right)}/4<1`$. Similarly, for the specific heat of Eqs. (48) and (50) at the repellor we obtain the numerical results as

$$c_b^{\left(n\right)}\stackrel{}{=}\left(0.261574\right)\left(a^{\left(2\right)}\right)^n+\left(0.201894\right)g^{\left(1\right)}\left(1\frac{1}{\left(d^{\left(1\right)}\right)^n}\right)\left(a^{\left(2\right)}\right)^n+\left(0.371342\right)g^{\left(1\right)}\left[g^{\left(2\right)}g^{\left(3\right)}\right],$$
(79)
and

$$c_{bb}^{\left(n\right)}\stackrel{}{=}\left(0.115082\right)\left(a^{\left(2\right)}\right)^n+\left(0.154011\right)g^{\left(1\right)}\left(1\frac{1}{\left(d^{\left(1\right)}\right)^n}\right)\left(a^{\left(2\right)}\right)^n,$$
(80)
with $`a^{\left(2\right)}=\left(d^{\left(1\right)}\right)^2/4<1`$, $`g^{\left(1\right)}=d^{\left(2\right)}/\left(\left(d^{\left(1\right)}\right)^2d^{\left(1\right)}\right)`$, $`g^{\left(2\right)}=a^{\left(2\right)}\left(1\left(a^{\left(2\right)}\right)^n\right)/\left(1a^{\left(2\right)}\right)`$, and $`g^{\left(3\right)}=a^{\left(1\right)}\left(1\left(a^{\left(1\right)}\right)^n\right)/\left(1a^{\left(1\right)}\right)`$.

For the limit of infinite decoration-level, we have $`\left(a^{\left(1\right)}\right)^n0`$, $`\left(a^{\left(2\right)}\right)^n0`$, and $`1/\left(d^{\left(1\right)}\right)^n0`$ as $`n\mathrm{}`$. Then, in this limit Eqs. (77) and (79) reduce to

$$ϵ_b^{\left(\mathrm{}\right)}\stackrel{}{=}\left(\frac{a^{\left(1\right)}}{1a^{\left(1\right)}}\right)=0.723079,$$
(81)
and

$$c_b^{\left(\mathrm{}\right)}\stackrel{}{=}\left(0.371342\right)g^{\left(1\right)}\left[\left(\frac{a^{\left(2\right)}}{1a^{\left(2\right)}}\right)\left(\frac{a^{\left(1\right)}}{1a^{\left(1\right)}}\right)\right]=0.639852,$$
(82)
while the long-range parts vanish, $`ϵ_{bb}^{\left(\mathrm{}\right)}\stackrel{}{=}0`$ and $`c_{bb}^{\left(\mathrm{}\right)}\stackrel{}{=}0`$

From these results we may conclude that the system of an $`n`$-lattice reduces to that composing of independent $`n`$-bonds as $`n\mathrm{}`$. In this limit, our results also show that the critical specific heat becomes a cusp with the height $`c^{\left(\mathrm{}\right)}=0.639852`$ at the critical point. Thus, the corresponding critical exponent of the specific heat, $`\alpha `$, is the same as that of the diamond-hierarchical Ising model, $`\alpha =`$ $`0.67652[9,23]`$.

## 6 summary and discussion

We use the Ising model defined on square lattices with diamond-type bond-decorations to study the ferromagnetic phase transition of an inhomogeneous system. In this model, the decoration level $`n`$ associated with a decorated lattice can be used as a parameter to characterize the degree of the inhomogeneity of a system. Based on the exact solutions of the two-dimensional Ising mobel, we employ the bond-renormalization scheme to obtain the analytic form of the free energy. The free energy and its derived quantities, the internal energy and specific heat, are expressed as the sum of those contributed by the local interactions among the sub-bonds of an $`n`$-bond and the long-range interactions among different $`n`$-bonds. Then, the change in the nature of the phase transition as $`n`$ increases can be viewed as the reflection of the competetion between the local and the long-range interactions.

In general, as $`n`$ increases, the long-range interaction is reduced and the local interaction is enhenced. The logarithmic singularity of the specific heat at the critical point is caused by the long-range interaction. Along with the reduction of the long-range interaction as $`n`$ gets larger, we have the critical temperature lowering down and the width of the critical region becoming narrower. By fitting the data numerically, we obtain the critical temperature changes with the decoration level $`n`$ as $`\left(k_BT_c/J\right)_n=1.6410+\left(0.6281\right)\mathrm{exp}\left[\left(0.5857\right)n\right]`$ and the critical amplitude of the specific heat changes as $`A_{\mathrm{sin}g}^{\left(n\right)}=\left(0.2465\right)4^n\mathrm{exp}\left[\left(1.6530\right)n\right]`$. Moreover, because of the enhancement of the local interactions, the occurring temperature of the local ordering increases with the form, $`\left(k_BT_m/J\right)_n=1.6410\left(0.1081\right)\mathrm{exp}\left[\left(0.4098\right)n\right]`$, as $`n`$ increases. Both $`\left(k_BT_c/J\right)_n`$ and $`\left(k_BT_m/J\right)_n`$ approach to the same limit $`1.6410`$ as $`n`$ approaches the infinite limit. Thus, the distance between $`\left(k_BT_c/J\right)_n`$ and $`\left(k_BT_m/J\right)_n`$ decreases with increasing $`n`$.

The cross over from a finite-decorated system to an infinite-decorated system is not a smooth continuation. For the infinite decoration-level, the critical temperature $`k_BT_c/J`$ is $`1.64101(8)`$ which is exactly the location of the repellor of the renormalization map of the variable $`\eta ^{\left(n\right)}`$. In this limit, the system of an $`n`$-lattice reduces to that of independent $`n`$-bonds, and the specific heat becomes a cusp with the value $`0.63985(0)`$ at the critical point instead of the logarithmic divergence.

It is interesting to compare the results we obtain with those from the cell-decorated triangular Ising model given in Ref. $`\left[14\right]`$. The critical region decreases with increasing $`n`$ for both cases. But, there exists some different features as the followings: $`\left(i\right)`$ The decreasing speed of the critical amplitude of the specific heat with increasing $`n`$ is slower for the bond-decoration case in comparing with that of the cell-decoration. $`\left(ii\right)`$ For the case of bond-decorations, the occurring temperature of the local ordering $`\left(k_BT_m/J\right)_n`$ of an $`n`$-lattice is smaller than the critical temperature $`\left(k_BT_c/J\right)_n`$ and the distance between $`\left(k_BT_m/J\right)_n`$ and $`\left(k_BT_c/J\right)_n`$ decreases with increasing $`n`$. However, for the case of cell-decorations, $`\left(k_BT_m/J\right)_n`$ is larger than $`\left(k_BT_c/J\right)_n`$ and $`\left(k_BT_m/J\right)_n`$ is freezed for $`n5`$. $`\left(iii\right)`$ There is a finite peak arising to the right of the critical point for the specific heat due to the local ordering in the cell-decoration case, but such a peak is not observed in the bond-decoration case.

We may use the idea of the order of ramification to explain the differences described in the above. To characterize the geometric difference between bond-decorated and cell-decorated lattices, we introduce the quantities, $`R_1`$and $`R_2`$, with $`R_1`$ for the cutting-bond number of isolating an $`n`$-bond or $`n`$-cell and $`R_2`$ for the average of the total bond number contained in an $`n`$-bond or $`n`$-cell. We notice that the $`R_1`$ value signifies the ability of establishing the long-range correlations among different $`n`$-bonds or $`n`$-cells and the interactions among $`R_2`$ bonds are referred as the local interactions.

For the cell-decorated triangular Ising model with decoration-level $`n`$, we have $`R_1=12`$ and $`R_2=3^{n+1}`$. Because of $`R_1R_2`$ for $`n3`$, the long-range correlation is very hard to establish, and this leads to the higher temperature for the occurrence of the local ordering in comparing with that of the long-range ordering. For sufficiently large $`n`$ with $`R_1/R_20`$, an $`n`$-cell almost become a closed subsystem, and thence a round peak, contributed by the subsystem, appears above the critical point in the curve of the specific heat.

For the diamond-type bond-decorated Ising model discussed in this work, we have $`R_1=62^n`$ and $`R_2=4^n`$. The values of $`R_1`$ and $`R_2`$ are comparable for $`n5`$, and the long-range effect is almost absent for $`n6`$. Thence, the round peak in the specific heat, which appears in the case of cell-decorations, is absent. Moreover, the coordination numbers of the decorated sites in an $`n`$-bond are highly inhomogeneous, and this is mainly responsible for the lower occurring temperature of the local ordering in comparing with that of the long-range ordering.

Thus, the non-universal features of phase transition for the diamond-type bond-decorated Ising model is quite different from that for the cell-decorated triangular Ising model.

## 7 acknowledgement

This work was partially suppoted by the National Science Council of ROC ( Taiwan ) under the Grand No. NSC 90-2112-M-033-002.

Table 1. The numerical results of the crical temperature $`\left(k_BT_c/J\right)_n`$, the occurring temperature of the local ordering $`\left(k_BT_m/J\right)_n`$, the critical values of the internal energy density $`\left(ϵ^{\left(n\right)}\right)_c`$, the internal energy density due to the local interactions $`\left(ϵ_b^{\left(n\right)}\right)_c`$, and the internal energy density due to the long-range interactions $`\left(ϵ_{bb}^{\left(n\right)}\right)_c`$, and the critical amplitude at the logarithmic divergence in the specific heat $`A_{\mathrm{sin}g}^{\left(n\right)}.`$

| n | $`(kT_c/J)_n`$ | $`(kT_m/J)_n`$ | $`\left(\epsilon ^{(n)}\right)_c`$ | $`\left(\epsilon _b^{(n)}\right)_c`$ | $`\left(\epsilon _{bb}^{(n)}\right)_c`$ | $`A_{\mathrm{sin}g}`$ |
| --- | --- | --- | --- | --- | --- | --- |
| 0 | 2.269185 | 0.833520 | -0.707107 | -0.414214 | -0.292893 | 0.247269 |
| 1 | 1.982072 | 1.255390 | -0.653281 | -0.541190 | -0.112085 | 0.189850 |
| 2 | 1.833652 | 1.449949 | -0.658949 | -0.614280 | -0.044674 | 0.140957 |
| 3 | 1.752268 | 1.540593 | -0.675920 | -0.657720 | -0.018202 | 0.102495 |
| 4 | 1.706095 | 1.586440 | -0.691315 | -0.683810 | -0.007508 | 0.073587 |
| 5 | 1.679371 | 1.610705 | -0.702627 | -0.699510 | -0.003119 | 0.052430 |
| 6 | 1.663721 | 1.623980 | -0.710252 | -0.708950 | -0.001301 | 0.037184 |
| 7 | 1.654492 | 1.631345 | -0.715166 | -0.714626 | -0.000544 | 0.026299 |
| 8 | 1.649027 | 1.635501 | -0.718250 | -0.718022 | -0.000228 | 0.018570 |
| 9 | 1.645783 | 1.637855 | -0.720154 | -0.720055 | -9.54965(-5) | 0.013100 |
| 10 | 1.643854 | 1.639191 | -0.721316 | -0.721280 | -4.00447(-5) | 0.009235 |
| 11 | 1.642707 | 1.639962 | -0.722020 | -0.722003 | -1.67971(-5) | 0.006509 |
| 12 | 1.642024 | 1.640400 | -0.722444 | -0.722433 | -7.04693(-6) | 0.004586 |
| 13 | 1.641617 | 1.640664 | -0.722699 | -0.722697 | -2.95673(-6) | 0.003231 |
| 14 | 1.641375 | 1.640809 | -0.722852 | -0.722849 | -1.24066(-6) | 0.002276 |
| 15 | 1.641231 | 1.640893 | -0.722943 | -0.722939 | -5.20605(-7) | 0.001604 |
| 16 | 1.641145 | 1.640946 | -0.722998 | -0.723000 | -2.18461(-7) | 0.001130 |
| 17 | 1.641093 | 1.640975 | -0.723030 | -0.723030 | -9.1674(-8) | 7.95678(-4) |
| 18 | 1.641063 | 1.640993 | -0.723050 | -0.723050 | -3.84699(-8) | 5.60485(-4) |
| 19 | 1.641045 | 1.641003 | -0.723061 | -0.723060 | -1.61435(-8) | 3.9481(-4) |
| 20 | 1.641034 | 1.641009 | -0.723068 | -0.723070 | -6.7745(-9) | 2.78106(-4) |


format: text

text:
# Simulations of Buoyant Bubbles in Galaxy Clusters

## 1 Introduction

The X-ray surface brightness of many clusters of galaxies shows a strong central peak which is generally interpreted as the signature of a cooling flow (Cowie & Binney 1977, Fabian & Nulsen 1977, Sarazin 1988, Fabian 1994). However, the simple cooling flow model conflicts with a growing number of observations that show that while the temperature is declining in the central region, gas with a temperature below $``$1 keV is significantly less abundant than predicted (Peterson et al. 2001, Tamura et al. 2001). It has been shown that star formation rates as inferred from observations only amount to a few percent of the predicted mass deposition rates of 1 - 1000 $`M_{}`$/yr (e.g. McNamara 1997, Crawford et al. 1999). Moreover, spectral evidence gathered by XMM-Newton has shown inconsistencies of the observed spectra with those predicted by cooling flow models (Böhringer et al. 2002).

Recently, two main candidates for the heating of the gas in the central regions of clusters have emerged: (i) heating by outflows from active galactic nuclei (AGN) (Binney & Tabor 1995, Tabor & Binney 1993, Churazov et al. 2001, Brüggen et al. 2002, Brüggen & Kaiser 2002, Basson & Alexander 2002, Reynolds, Heinz, & Begelman 2001), and (ii) transport of heat from the outer regions of the cluster by thermal conduction (Zakamska & Narayan 2002, Narayan & Medvedev 2001, Fabian, Voigt & Morris 2002, Ruszkowski & Begelman 2002, Voigt et al. 2002, Gruzinov 2002, Bertschinger & Meiksin 1986). A third possibility that has been discussed is heating by magnetic fields through reconnection (e.g. Soker & Sarazin 1990).

In this paper we will consider the first scenario. Radio-loud active galactic nuclei (AGN) drive strong outflows in the form of jets that inflate bubbles or lobes. The lobes are filled with hot plasma, and can heat the cluster gas in various ways. In the case of very energetic jets, the expansion of the lobes will be supersonic, and the resulting strong shock will heat and compress the gas. However, there is mounting evidence that weaker jets, which are presumably much more common, do not lead to efficient shock heating. Nevertheless, they still produce pockets of low-density gas that displace the hot intra-cluster medium (ICM) and thus appear as depressions in the X-ray surface brightness. High-resolution X-ray observations of cooling flow clusters with Chandra have revealed a multitude of X-ray holes on scales $`<`$ 50 kpc, often coincident with patches of radio emission, e.g. in Abell 2052 (Blanton et al. 2001), Perseus (Böhringer et al. 1993, Churazov et al. 2000, Fabian et al. 2000), Cygnus A (Wilson, Young, & Shopbell 2000 and Carilli, Perley, & Harris 1994), Abell 2597 (McNamara et al. 2001), MKW3s (Mazzotta et al. 2002), RBS797 (Schindler et al. 2001), HCG62 (Vrtilek et al. 2000) and Centaurus A (Saxton, Sutherland, & Bicknell 2001). In addition to radio-bright cavities, radio-faint “ghost” cavities have been found, e.g. in Perseus, which presumably date from previous radio episodes of the AGN. The pressure support in the cavities is not well understood. In Abell 2052 the external gas pressure is about 10 times higher than the internal equipartition pressure. The cavities found in MKW3s, however, are filled with hot thermal gas that provides pressure support against the ICM and the radio observations of Hydra A also suggest that the bubbles are in local pressure equilibrium with the surroundings. The ubiquity of cavities in clusters means that they must be relatively long-lived and that the majority must be in rough pressure balance with their surroundings. However, the source of this pressure, the structure of the magnetic fields inside the bubbles and the equation of state of the hot plasma are still uncertain.

Chandra images of Hydra A, Perseus and Abell 2052 show soft X-ray emission form the rims of the bubbles which imply that the gas around the rims of these bubbles is colder and not hotter than the ambient gas (e.g. Fabian et al. 2000, 2001, Nulsen et al. 2002). The existence of these cold rims is taken as evidence that the gas has not been shocked during the inflation of the bubbles but that the cold rims consist of low-entropy gas that has been displaced from the centre of the cluster (Brighenti & Mathews 2002). Therefore, the radio lobes must have expanded gently into the ICM. The buoyant rise and the subsequent mixing of high-entropy gas with the cluster material can lead to a significant re-structuring of the inner regions of a cooling flow (e.g., Brüggen & Kaiser 2002, Basson & Alexander 2002).

Assuming that a typical central galaxy produces 10-100 bubbles in its life time, each carrying an energy of up to $`10^{59}`$ ergs, the total energy in bubbles can be comparable to the thermal energy of the X-ray emitting plasma in the inner cluster regions (McNamara 2002). Moreover, there is evidence that the total power in radio jets is substantially larger than the radio power of the lobes that they feed (Bicknell & Begelman 1996). Thus the question arises of where the energy ultimately goes and what fraction of the energy is deposited into the ICM. This is the question that is addressed in this paper.

Both, the rise and the mixing are complex processes that are difficult to treat on the basis of analytical estimates alone. We therefore performed highly resolved two-dimensional hydrodynamic simulations of buoyant gas in a typical cluster environment similar to previous simulations of lower resolution. Only two-dimensional simulations were performed because our primary aim was to achieve a very high spatial resolution in order to minimize the effects of numerical diffusivity that have hampered previous simulations. This was also aided by the use of an adaptive grid that places resolution elements only where they are needed, such that an effective resolution of 2048 times 2048 elements could be achieved. Experience has shown that large-scale transport properties are not very different in two dimensional and three-dimensional simulations, so we expect our conclusions to hold even for the three-dimensional reality.

In this paper, we will extend the work presented in Brüggen & Kaiser (2002). Apart from a much more detailed analysis, we have extended previous work by simulating an ambipolar jet in the entire cluster compared to merely a section of the cluster. Thus, any artificial surfaces of symmetry are avoided (as also in Reynolds, Heinz, & Begelman 2001; Basson & Alexander 2002) and the amount of global, distributed heating can be studied. Moreover, we ran the simulations for a longer physical time, namely for much longer than the central cooling time of the cluster. We also included radiative cooling and treated the radio plasma and the ambient plasma as separate fluids with different equations of state.

## 2 Model

The initial density was assumed to follow the dark matter, which is given by a modified NFW profile (Navarro et al. 1997)

$$\rho (r)=\frac{\rho (0)r_s}{(r+r_s)(1+r/r_c)^2},$$
(1)
where $`r`$ is the distance from the centre, $`r_c`$ the core radius and $`r_s`$ a softening radius, which we assumed to be proportional to the core radius, i.e. we chose $`r_s=r_c/30`$. The parameters of our model are $`n_\mathrm{e}(0)=0.3`$ cm<sup>-3</sup>, $`r_c=300`$ kpc, $`r_s=r_c/30`$ and the baryon fraction is 10%. The temperature and pressure are calculated assuming initial hydrostatic equilibrium. The initial model of our simulations is shown in Fig. 1. Thus our model has got a central cooling time of around 20 Myrs.

Inspired by a growing number of observations of X-ray cavities in the ICM, we introduce heating at off-centre sites by ad-hoc bubbles. To mimic the energy input by a jet, hot, buoyant gas was injected continuously into two small spherical regions which lie on opposite sides of the centre of the cluster. The injection region had a radius of 0.5 kpc, and the rate of energy input was $`L=6\times 10^{44}`$ which is a typical kinetic luminosity of a radio jet (Owen, Eilek, & Kassim 2000). This rate is also near the maximum rate that allows the bubbles to detach themselves from the injection region in the course of their rise. For higher luminosities a large bubble is inflated that evacuates most of the cluster centre before it can rise and become detached from the injection region. Which mechanism is responsible for tuning the energy of the radio jets remains unclear. The energy input was switched off after 208 Myr which is estimated to be a typical time for a recurring FRI source (Saxton, Sutherland, & Bicknell 2001). We assumed that the temperature of the injected fluid was 100 times greater than the temperature of the ambient fluid at the same location and that the gas had zero initial velocity.

We treated the hot gas in the bubbles and the ambient cluster gas as two separate fluids which obey two separate equations of state. The ICM was assumed to obey a $`\gamma =5/3`$ polytropic equation of state. However, the equation of state within the bubbles is unknown. We assume that the magnetic field within the rarefied bubbles is tangled on small scales and can be approximated together with the relativistic particles by a $`\gamma =4/3`$ equation of state. Thus, the adiabatic index of the injected fluid in our simulation was $`4/3`$. The ambient gas and the gas in the bubbles is allowed to cool via thermal bremsstrahlung and line radiation. The cooling function of the gas was taken from Sarazin (1986) and Raymond (1976) (see Fig. 2). We have neglected the synchrotron losses of the radio plasma assuming it to be very hot so that it hardly cools at all. Finally, we have produced purely hydrodynamic simulations and, therefore, magnetic fields have been ignored. The only effect of magnetic fields that has implicitly been taken into account is the change of the adiabatic index of the injected fluid as described above.

The simulations were performed with the FLASH code using an explicit solver based on the piecewise-parabolic method (PPM). This method is accurate to second order in both time and space and uses a monotonicity constraint and contact steepening. For a detailed description of the FLASH code, the reader is referred to Fryxell et al. (2000). Radiative cooling was treated explicitly with a separate timestep limiter. The code uses an adaptive grid which refines and derefines the computational mesh as needed. The computational domain span 140 kpc $`\times `$ 140 kpc and was centred on the centre of the cluster. The effective resolution in our simulations was 2048 $`\times `$ 2048 zones. On all boundaries outflow conditions were enforced. All simulations were run on 16 processors on an SGI Origin 3000 and on an IBM Regatta system.

## 3 Results and Discussion

Figure 3 shows logarithmic contour plots of the density at 84, 167, 252 and 336 Myr after the start of the energy injection. In the last two pictures the injection of energy has ceased. One can see how the buoyant plume rises upwards while Rayleigh-Taylor instabilities form mushroom-shaped plumes on the bubble surfaces. Kelvin-Helmholtz instabilities at the boundary between the radio plasma and the ICM leads to the formation of tori and the shedding of vortex rings. The tori form when the rarefied gas in the wake of the bubble penetrates the bubble from below and makes contact with the thermal gas on the upper side of the bubble (Saxton, Sutherland, & Bicknell 2001). The rising tori cause the thermal gas to circulate around them. The velocity of the bubbles increase with increasing distance from the centre as expected (Churazov et al. 2001). A stream of denser gas trails the bubbles and marks their path. The largest bubbles are found at the leading edge of the buoyant gas. One can see that as the bubble rises upwards the space between the top bubble and the centre gets scattered with smaller radio bubbles and vortex rings. Thus the bubbles are shredded producing a broad distribution of sizes and velocities where the smaller bubbles rise more slowly. The final fate of the bubbles, i.e. whether they will rise out of the cluster before they are mixed microscopically into the ICM, is difficult to tell from these simulations and will depend crucially on the strength and structure of the magnetic fields inside the bubbles. The simulations show that denser gas is pushed outwards despite the effect of radiative cooling, which would cause the gas to flow inward. It is remarkable that bubbles are found at fairly large distances from the jet axis leading to a broad distribution of bubbles. From an injection region of width of 1 kpc the bubbles spread laterally to distances of up to 30 kpc - which corresponds to 1/10th of the core radius - from the axis.

Thus the positions of X-ray cavities do not always allow simple inferences on the orientation of the jet axis. In some cases bubbles have been observed that are located off-centre from the current jet axis and this has led to suggestions that the jet axis has turned. However, this does not always have to be the case. Some bubbles may merely have migrated away from the axis due to instabilities. This may explain why in some clusters, such as Perseus, the X-ray holes appear to be randomly distributed about the centre.

In Fig. 4 we show a map of the X-ray surface brightness at a time of 100 Myrs after the start of the simulation. For this purpose we have rotated the 2D output around the jet axis and integrated the bremsstrahlung emissivity along lines of sight. Since the bremsstrahlung emissivity goes as $`n_\mathrm{e}^2T^{1/2}`$ the X-ray maps show a very luminous core. Therefore, the cavities are more difficult to observe at the centre but become more prominent once they have reached a height of roughly a tenth of the core radius. The surface brightness image also shows a bias towards larger bubbles while many of the smaller structures remain obscured.

Figure 5 shows a contour plot of the temperature at 120 Myr after the start of the energy injection. A trunk of colder gas is uplifted in the wake of the bubbles. One can note how colder gas has been entrained and is lifted upwards by the bubbles. As mentioned in Sec. I, observations have shown that gas immediately around the bubbles is colder than the ambient gas (e.g. Fabian et al. 2000, 2001, Nulsen et al. 2002). In Fig. 6 we show an enlarged contour plot of the temperature in the vicinity of the bubbles. One can clearly see how entrained, colder gas from the central regions is being pushed out by the buoyant bubbles. The bubbles are surrounded by colder, low-entropy gas which appears bright in X-ray images. These rims have got widths of $``$ 1 - 5 kpc. Their widths decrease with increasing distance from the centre which can be attributed to three factors: First, as the bubbles accelerate, they compress the gas in front of them; second, the material in the rims cools faster than the ambient ICM and third, some of the rim material is being mixed into the ambient ICM. The rims are 20% to 40% cooler than the adjacent thermal gas which is consistent with observations (Nulsen et al. 2002). They also become cooler with their distance from the centre because they cool faster. Our simulations confirm the suggestion by Brighenti & Mathews (2002) who identified the source of the bright rims in low-entropy gas from the cluster centre.

We have found that without cooling the average temperature in the core goes up because hotter material from outer regions is mixed into the core and the expanding bubbles are performing work on the ambient medium (Brüggen & Kaiser 2002). For the simulations presented here we have plotted the temperature evolution of the cluster gas in Fig. 7. The temperature has been averaged over spherical shells of thickness 10 kpc about the cluster centre. The code treats the ambient and the injected material as two separate fluids throughout, which enables us to study the temperature evolution of the original cluster gas alone. Only the ambient cluster gas has been included in the average shown in Fig. 7. In a separate run, we have switched off the energy input by the bubbles and the corresponding temperature evolution is shown in the right panel of Fig. 7. The comparison shows clearly how, over the time scales considered here, the bubbles manage to keep the temperature in the central regions of the cluster constant (or even raise them), whereas in the run without energy input we see the usual radiatively-regulated cooling inflow. Fig. 8 shows the evolution of the specific entropy in the cluster. Again, the entropy has been averaged over spherical shells of thickness 5 kpc. One can notice how in all but the innermost shells the entropy is kept at a fairly constant level or even rises despite the cooling. This effect caused by the rising bubbles persists for a time of more than 400 Myrs which is twice as long as the time for which the jet was active and approximately the sound crossing time of the core radius $`r_c`$. In the cooling-flow model, mass deposition is greatest in the core, where the cooling time drops sharply. The mass that condenses out in the cores ceases to provide pressure support, and thus causes a slow inflow. Therefore, the cooling flow can be delayed significantly if the cooling of gas in the central regions can be halted. In order to gauge what fraction of the energy of the injected material is transferred to the ICM, we have repeated our simulation with the cooling switched off. We then calculate the total internal energy of the ambient gas and find that after 200 Myrs $`10`$ % of the injected energy has been deposited in the cluster gas. This calculation ignores that the injection of hot material might cause enhanced cooling in some of the cluster gas but we estimate that this effect is not very large.

One can conceive of various physical mechanisms whereby the energy of the hot and buoyant fluid is transferred to the ICM: (i) $`pdV`$ work on the ICM by the expanding bubbles (ii) entrainment of ambient fluid, (iii) microscopic mixing between bubble and ambient fluid, (iv) turbulent excitation of sound and gravity waves and (v) thermal conduction.

As pointed out in Churazov et al. (2002) the work done by the rising bubble can be estimated by considering the enthalpy of the bubble during adiabatic expansion

$$H=\frac{\gamma }{\gamma 1}pV=H_0\left(\frac{p}{p_0}\right),$$
(2)
where $`H_0`$ and $`p_0`$ are the enthalpy and pressure at the initial position. Taking a typical pressure profile, approximately half the energy is lost at a distance of $``$ 20 kpc from the origin (Churazov et al. 2002). This work goes mainly into $`pdV`$ work on the ambient medium and kinetic and potential energy of the entrained gas. Given the steep entropy profile in the inner regions of the cluster, large masses of entrained gas cannot travel a long distance with the rising bubbles, and therefore most of the energy will stay in the cluster.

Microscopic mixing is not included in our code and any small-scale mixing that occurs in the simulation must therefore be entirely numerical. Convergence tests with our code have shown that any numerical mixing is negligible.

Sound waves will only have a minor share in the energy lost since the ratio between the energy in sound waves over the dissipated energy for subsonic turbulence goes approximately as (Landau & Lifshitz 1987)

$$\frac{ϵ_s}{ϵ_d}\left(\frac{v}{c_s}\right)^51.$$
(3)
Since most of the turbulent motions induced by the buoyant bubbles are subsonic, this fraction is going to be small. Estimates made in Churazov et al. (2002) show that the energy emitted in sound waves constitutes not more than a small fraction of the internal energy of the bubble.

Finally, thermal conduction is not included in our simulations. The role of thermal conduction in clusters is a matter of contention, with most observations pointing towards a heavily suppressed thermal conductivity in the ICM (e.g. Fabian 1994). Hence, we conclude that $`pdV`$ work by the rising and expanding bubbles as well as turbulent entrainment are the most likely mechanisms whereby the radio plasma deposits energy in the ICM.

Clearly, the lifetime of the activity of the central AGN is short compared to the evolutionary time scale of the cluster gas. Therefore, once the AGN stops supplying energy to the buoyant bubbles, the cluster gas will settle down once again and a full cooling flow may be re-established. The cooling gas flowing to the centre of the cluster may then trigger a further active phase of the AGN. Thus a self-regulating process with cooling periods alternating with brief bursts of AGN activity may be established. The rising bubbles uplift colder material from the vicinity of the AGN and thus disrupt the supply of fresh fuel for the radio jets. This feedback mechanism might automatically regulate the power of the radio jets. The details of this mechanism remain unclear but it is potentially important as it might provide an upper mass limit for galaxies. For more on feedback between gas accretion and radio power see, e.g. Quilis, Bower, & Balogh (2001), Voit & Bryan (2001), Ruszkowski & Begelman (2002), Churazov et al. (2002b).

It has been found that 71 % of all cD galaxies at the centres of cooling-flow clusters show evidence of radio activity (Burns 1990). This fraction is higher than in non-cooling flow galaxies which again points towards the existence of some form of feedback mechanism. If all AGN at the centres of cooling flows go through activity cycles of constant length, $`t_{\mathrm{on}}`$, separated by periods of quiescence, again of constant length, $`t_{\mathrm{off}}`$, then these observations imply that $`t_{\mathrm{off}}<0.41t_{\mathrm{on}}`$. Figure 8 indicates that the simulated cluster starts to relax back to its original state after a time of $`400`$ Myr and returns to a cooling dominated inflow. Thus it gets close to a configuration similar to that at which the AGN activity started. Despite having chosen one particular realisation of the gas distribution inside the cluster, these findings are insensitive to the details of the cluster model. Experiments with other cluster models have yielded very similar results. Here we have only presented results of 2D simulations because we wanted to achive the highest effective resolution that we could afford. However, previous 3D simulations (Brüggen et al. 2002) of lower resolution have shown that the gross features and integrated quantities do not differ greatly from 2D simulations.

Cosmological simulations as well as observations suggest that the gas in many clusters may not be in an equilibrium state, in contrast to what we assumed in our simulation. However, gas flows of low Mach numbers within a cluster are unable to impede the rise of the bubbles, and are more likely to increase than to decrease the mixing. Moreover, clusters that are far from hydrostatic equilibrium are generally not expected to host cooling flows.

In summary, we have found that bubbles can heat the cluster core by input of mechanical energy. The effect of the bubbles on the ICM consists of heating via $`PdV`$ work and redistribution of mass via buoyancy-driven mixing. For realistic parameters the heating is shown to be sufficient to balance radiative cooling and to disrupt, at least temporarily, the cooling flow. Thus a significant fraction ($`10`$ %) of the energy residing in radio lobes can be dissipated in the cooling flow. Clearly, heating by AGN will not account for the entirety of the missing gas in cooling flows but it is shown to be an efficient and common mechanism which introduces heat into the ICM. During their rise in the cluster potential, the bubbles expand and produce depressions in the X-ray surface brightness. Moreover, hydrodynamical instabilities produce a broad distribution of smaller bubbles and vortices. Bright rims are a common feature around buoyant bubbles and consist of colder gas that has been uplifted from the centre. The bright rims seen in our simulations are consistent with those found in observations.

The software used in this work was in part developed by the DOE supported ASCI/Alliances Center for Thermonuclear Flashes at the University of Chicago. I thank Christian Kaiser, Bill Forman and Eugene Churazov for helpful discussions.