text:
# Non conservative Abelian sandpile model with BTW toppling rule

## I Introduction

Recently Tsuchiya an Katori  have introduced a non conservative Abelian sandpile model with a toppling rule similar to that of the well known Bak, Tang and Wiesenfeld (BTW) sandpile model . The model is defined in a square lattice of size $`L`$ and an integer energy profile $`z_{ij}`$ is considered. Sites with energy below the threshold $`z_c=4\alpha \zeta `$ are stable. If the energy at any given site $`(i,j)`$ exceeds this threshold the site transfer energy to its four nearest neighbors following the toppling rule: $`z_{ij}z_{ij}z_c`$, $`z_{i\pm 1j}z_{i\pm 1j}+\zeta `$ and $`z_{ij\pm 1}z_{ij\pm 1}+\zeta `$. $`\zeta `$ is an integer number and $`\alpha `$ is such that $`z_c`$ is also integer. The boundaries are assumed to be open and the system is perturbed by adding a unit of energy at a site selected at random and letting it evolve until an equilibrium configuration is reached.

On each toppling event an amount of energy $`ϵ=4\zeta (\alpha 1)>0`$ is dissipated. For $`\alpha =1`$ the model is conservative, it is just the BTW model but with a different scale of toppling and energy addition. In the BTW model  the energy added to perturb the system is 1 and on toppling an active site transfers an energy 1 to each neighbor, they are of the same order. In the model defined above the energy added is still 1 but the energy transferred on toppling is $`\zeta `$. In the limit $`\alpha =1`$ and $`\zeta =1`$ the BTW model  is recovered while for $`\alpha =1`$ and $`\zeta 1`$ it is similar to the BTW model but with a uniform driving.

In the BTW model ($`\alpha =1`$ and $`\zeta =1`$) the avalanches can be decomposed in a sequence of sub avalanches called waves  with well defined finite size scaling properties. On the contrary, the distribution of the overall avalanche size $`s`$ is better described using a multi-fractal analysis . The break down of the finite size scaling has been recently shown to be a consequence of the existence of correlations in the sequence of waves .

For $`\alpha >1`$ the model is non conservative but still Abelian . In the thermodynamic limit $`L\mathrm{}`$, exact calculations by Tsuchiya and Katori yield the mean avalanche size (including avalanches with size zero) $`T=ϵ^1`$ . Since $`ϵ=4\zeta (\alpha 1)`$ they concluded that in the limit $`\alpha 1`$ $`T`$ diverges. However, as it is shown in sec. II this conclusion is wrong because $`\alpha `$ cannot goes to zero in an arbitrary way, in order to satisfy the constraint that $`z_c=4\alpha \zeta `$ remains integer. Here it is demonstrated that $`ϵ1`$ and, therefore, $`T1`$ for all possible values of $`\zeta `$ and $`\alpha `$.

The main goal of this work is to investigate the scaling properties of this non conservative BTW like model in the limit $`\alpha 1`$. The main questions are related with the existence or not of criticality in the conservative limit $`\alpha 1`$ and if in this limit one recovers the conventional BTW model $`\alpha =1`$. From the analysis of the energy scales involved in the model ( sec. II) and numerical simulations (sec. III) it is concluded that the model is critical when $`\alpha 1`$ but it does not belong to the universality class of the BTW model. Its relation with other non conservative models with BTW like toppling rule introduced in the literature  is also discussed (sec. III).

## II Scaling laws

In this section some scaling laws are derived based on the energy scales involved in the model. The main idea of this approach is that the balance between input and dissipation of energy determine the scaling of some magnitudes with the dissipation per toppling event, following the general guidelines introduced by Vespiganani et al . For this purpose the avalanches are assumed to be instantaneous and the analysis is focused in the time scale of the driving field. On each step one adds 1 unit of energy and measure the toppling activity and the energy dissipated. On each toppling event an amount of energy $`ϵ=4\zeta (\alpha 1)`$ is locally dissipated while an amount $`4\zeta `$ is transferred to nearest neighbors. For boundary sites part of the energy is also dissipated through the boundary.

Let $`G(r)`$ be the Green function , the average number of toppling events at a distance $`r`$ from the site where the energy was added. Close to $`r=0`$ the effect of local dissipation gives an small contribution and the main energy scale is given by the transport of the energy from active sites to their nearest neighbors. On the contrary, far from $`r=0`$ the effects of the local dissipation becomes more important. How far will depend on certain correlation length $`\xi `$, such that for $`r\xi `$ transport is more important than local dissipation while for $`r\xi `$ the opposite occurs.

Thus, there are two characteristic lengths in this model: the system size $`L`$ and the correlation length $`\xi `$. The analysis developed above is valid in the thermodynamic limit $`L\xi `$. In this case the dissipation through the boundary of the system is negligible in comparison with the energy dissipated on each toppling event. In such a situation the only way to reach an stationary state is to balance the input of energy from the driving field with the energy locally dissipated. Moreover, since $`\xi `$ is the only characteristic length $`G(r)`$ is expected to satisfy the scaling law

$$G(r)=r^{\eta d}(r/\xi ),$$
(1)
where $`\eta `$ is an scaling exponent and $`d`$ is the spatial dimension.

The amount of energy $`\delta E_d(r)`$ locally dissipated inside an hyper-circle of radius $`r`$ is

$$\delta E_d(r)ϵ_0^r𝑑\rho \rho ^{d1}G(\rho )ϵ\xi ^\eta f(r/\xi ),$$
(2)
where $`f(x)=_0^x𝑑yy^{\eta 1}(y)`$ and the second proportionality is obtained using eq. (1). On the other hand, the average energy transported through its boundary $`\delta E_t(r)`$ is given by

$$\delta E_t(r)\zeta r^{d1}\frac{dG}{dr}(r)\zeta \xi ^{\eta 2}g(r/\xi ),$$
(3)
where $`g(x)=d[x^{\eta 2}(x)]/dx`$ and the second proportionality is obtained using eq. (1). The correlation length $`\xi `$ can be defined as the radius $`r`$ at which this two contribution becomes of the same order. With this definition and equating eqs. (2) and (3) with $`r=\xi `$ it results that

$$\xi \left(\frac{\zeta }{ϵ}\right)^\nu ,\nu =\frac{1}{2}.$$
(4)
On the other hand, on each step 1 unit of energy is added and in average the amount $`ϵT`$ is dissipated, where $`T_0^{\mathrm{}}𝑑rr^{d1}G(r)`$ is the mean avalanche size counting even the case when it has size 0. Equating this two contributions it results that

$$T=\frac{1}{ϵ}=\frac{1}{4\zeta (\alpha 1)}.$$
(5)
Moreover, using eq. one obtains (1)

$$\eta =0.$$
(6)
Eq. (5) reproduces the exact result by Tsuchiya and Katori. The present approach is however based on more general arguments an can be easily adapted to any sandpile model with local dissipation. The same argument (energy balance) has been previously used by Vespiganani et al  to understand the scaling properties of other sandpile models with local dissipation. Here, a slightly different approach has been considered where the new parameter $`\zeta `$, the ratio between the energy received from nearest active neighbors and from the driving field, has been taken into account.

From eq. (5) Tsuchiya and Katori concluded that when $`\alpha 1`$ $`T`$ diverges. However, this conclusion is not valid if $`z_c=4\zeta \alpha `$ is restricted to be an integer number. To show this let us write $`\alpha =1+ϵ/4\zeta `$ which follows from eq. (5). But $`4\zeta \alpha =4\zeta +ϵ`$ is restricted to take integer values. With $`\zeta `$ being an integer number the only way to satisfy this requirement is that $`ϵ`$ is also integer, i.e. $`ϵ=1,2,3,\mathrm{}`$. Then, since the smaller non negative integer is 1 it is concluded that $`ϵ1`$ and, therefore, from eq. (5) $`T1`$, i.e. it is bounded.

Nevertheless, the correlation length $`\xi `$ in eq. (4) does not only depend on $`ϵ`$ but also on $`\zeta `$. For fixed $`ϵ`$ it diverges in the limit $`\zeta \mathrm{}`$ and the model is critical. The real control parameter is then $`ϵ_{eff}=ϵ/\zeta `$, i.e. the energy dissipated per toppling event relative to the characteristic energy scale of transport $`\zeta `$. Although this result is in complete agreement with the field theory approach of Vespignani et al  the fact that $`T`$ does not diverges when $`ϵ_{eff}0`$ ($`\zeta \mathrm{}`$) excludes this model from their analysis.

Moreover, in previous sandpile models conservation implies the scaling law $`s\xi ^2`$, where $`s`$ is the mean avalanche size excluding those with size 0 . To investigate the validity of such scaling relation for the present model let us take into account that $`s`$ is related to $`T`$ through the expression

$$s=T/P_a,$$
(7)
where $`P_a`$ is the probability to obtain an avalanche with non zero size. In the models considered by Vespignani et al  $`\zeta =1`$ and, therefore, from eqs. (4), (5) and (7) it results that $`s\xi ^2/P_a`$. Moreover, in this models $`P_a`$ has a finite value and, therefore, one obtains the mentioned scaling law $`s\xi ^2`$.

On the contrary, in the model considered here $`s`$ can not be related to $`\xi `$ using these arguments. For fixed $`ϵ`$ from eqs. (5) and (7) one obtains that $`s1/P_a`$. Thus, from the energy balance invoked above we cannot say anything about the scaling of $`s`$ with $`\xi `$ (an exponent 2 will be an accidental coincidence) and, therefore, this model belongs to a new universality class.

## III Numerical simulations and discussion

In this sections results obtained from numerical simulations of the model studied above are presented. The simulations were performed using $`ϵ=1`$, $`L=4096`$ and $`\zeta `$ ranging from $`2^0`$ to $`2^{10}`$ ($`ϵ_{eff}=1/\zeta `$ ranging from 1 to $`2^{10}`$). For these values the condition $`L\xi `$ was observed to be satisfied. Statistics was taken over $`10^8`$ avalanches after the system reached the stationary state.

Before entering in the analysis of the statistics of the avalanches let us check the validity of eq. (5). The log-log plot of $`s`$ vs. $`\zeta `$ is shown in fig. 1. A clear linear behavior is observed for $`\mathrm{log}_{10}\zeta 5`$ suggesting that above this value simple scaling applies. On top of this points the numerically computed values of $`1/P_a`$ are plotted obtaining an overlap in agreement with eq. (5). If the scaling relation $`s\xi ^2`$ where valid, using eq. (4), $`s\zeta `$. However, a linear fit to this log-log plot gives a slope $`0.9`$.

The fact that this scaling relation thus not hold is clearly shown in fig. 2, where the stationary energy distribution is shown. As it can be seen $`P_a=P_{z_c1}`$ does not decreases as $`1/\zeta `$ but slower, which explains why $`s`$ increases with $`\zeta `$ with a slope smaller than 1. The rest of the distribution scales like $`1/\zeta `$ which is just a consequence of the increase of the density of possible values of $`z`$.

The avalanche statistics will be characterized by: the number of toppling events $`s`$ and steps $`t`$ required to reach an stable configuration, the number of sites $`a`$ ”touched” by the avalanche, and the characteristic radius of the cluster formed by these sites $`r`$. The main goal of the simulations is to determine the probability densities $`p_x(x,\zeta )`$ ($`x=s,t,a,r`$) in the stationary state.

One can easily see that $`s=a`$, in other words sites topples only once within an avalanche. In this model, as a difference with the original BTW model, only one wave of topplings takes place. The first wave if generated from an initial site with height $`z=z_c=4\zeta +ϵ`$. When this site topples it transfers an amount equal to $`z_c`$ to its nearest neighbors and, therefore, ends with energy $`z=0`$. The best we can have to obtain a second toppling at this site is that its four nearest neighbors also become active. In such a case the initial side will receive $`4\zeta <z_c`$ units of energy, which is not enough to make it active again. Hence, no second wave will be obtained yielding $`s=a`$.

Since the waves are known to satisfy well defined finite-size scaling properties and in the present model an avalanche is made by one wave it is expected that the distributions $`p_x(x,\zeta )`$ also satisfy a finite-size scaling. However, the scaling exponents will not necessarily be those obtained for the scaling of waves because, in the present model, conservation does not introduce any scaling relation among the different scaling exponents.

If finite size scaling applies then these densities will satisfy

$$p_x(x,\zeta )=x^{\tau _x}𝒢[x/x_c(\zeta )],$$
(8)
where $`\tau _x`$ is the power law exponent characterizing the self-similar regime and $`x_c`$ is a cut off above which the distribution deviates from a power law and has a fast decay given by $`𝒢`$. The validity of this scaling form is supported by the numerical results. The cut off $`x_c`$ is determined by the existence of the characteristic length $`\xi \zeta ^\nu `$ and is expected to scale as $`x_c\xi ^{D_x}\zeta ^{d_x}`$, where $`d_x=D_x\nu `$ is an effective fractal dimension.

To compute the exponents $`\tau _x`$ and $`d_x`$ the moment analysis technique introduced by De Menech et al  is used. The moments of the probability density in eq. (8) are given by

$$x^q=_0^{\mathrm{}}𝑑xp(x)x^q\zeta ^{\sigma _x(q)},$$
(9)
where

$$\sigma _x(q)=(1\tau _x)d_x+d_xq.$$
(10)
The last equivalence in eq. (9) is valid for values of $`q`$ not to small, for which the precise form of $`p_x(x,\zeta )`$ at small $`x`$ is not important.

$`\sigma _x(q)`$ can be determined from a linear fit to the log-log plot of $`x^q`$ vs. $`\zeta `$. The resulting values using $`\zeta =2^5,2^6,\mathrm{},2^{10}`$ are shown in fig. 3. In all cases ($`x=s,t,r`$) for $`q`$ larger than 1 a well defined linear dependence is observed. From the linear fit (see eq. 10) to these plots the exponents $`\tau _x`$ and $`d_x`$ are computed. The results are shown in table I.

The exponent $`\nu `$ is very close to $`1/2`$ in very good agreement with the scaling arguments of previous section. On the other hand, $`d_s`$ is quite close to 1 which implies that the avalanche size (or area) scale as $`sr^2`$, i.e. avalanches are compact $`D_s=2`$. With this value, the scaling relation $`(2\tau _s)D_s=2`$ yields the power law exponent $`\tau _s=1`$ which is clearly in disagreement with the value computed numerically. The reason for this result is that conservation does not introduce any scaling relation as it generally occurs in sandpile models .

The exponents computed using the moment analysis technique can be checked using rescaled plots of the integrated distribution $`P_x(x,\zeta )=_x^{\mathrm{}}𝑑xp_x(x,\zeta )`$. The resulting plots are shown in figs. 4, 5, and 6. The scaling works quite good supporting the validity of the reported exponents.

In the literature we can find other sandpile models with local dissipation and BTW like toppling rule . In the models considered in  and  the energy profile is continuous and the dissipation rate per toppling event $`ϵ`$ is a control parameter which can take any real value and, therefore, can be tuned to zero. Another feature of these models is that only one wave of toppling can take place and, therefore, for any finite $`ϵ`$ the model is in a different universality class from that of the BTW model.

On the other hand in  the energy profile is discrete as in the original BTW model at the prize of introducing stochasticity in the model. In this case with a probability $`p`$ energy is fully dissipated yielding an average dissipation per toppling event $`ϵ=2dp`$. Clearly $`p`$ may take any real variable between 0 and 1 and, therefore, also in this case the dissipation per toppling event can be fine tuned to zero. As a difference with the models described in the previous paragraph, in this case multiple toppling of a site within an avalanche is possible, which make it closer to the original BTW model. Moreover, the use of finite size scaling techniques can be also questioned and a multi-fractal analysis may be more appropriate , which is another characteristic feature of the BTW model . All this elements together with the numerical results reported in  suggest that in the limit $`p1`$ ($`ϵ0`$) this model belong to the same universality class of the BTW model.

A common feature of all this models  is that $`sϵ^1`$ as predicted by the field theory approach of Vespignani et al , leading to the scaling relation $`(2\tau _s)D_s=2`$. On the contrary, in the present model the scaling of $`s`$ with $`ϵ_{eff}`$ is not known and conservation does not introduce the above scaling relation. Hence, the model introduced by Tsuchiya and Katori belongs to different class among sandpile models.

## IV Summary and conclusions

A non conservative Abelian sandpile model with a BTW like toppling rule has been studied. The model can be though as the only possible generalization of the BTW model to include local dissipation without introducing stochasticity in the toppling rule and keeping a discrete energy profile. However, the scaling approach and numerical the simulations reported here show that it does not belong to the universality class of the BTW model, not even to the universality class of any sandpile model previously considered in the literature.

## Acknowledgements

I thanks R. Pastor Satorras and A. Vespignani for useful comments and discussion. The numerical simulations where performed using the computing facilities at the ICTP.


format: text

text:
# On Non– Efficiency of Quantum Computer

## I Introduction

The standard description of a quantum computer is the following. We have a quantum system with $`N=2^n`$ orthogonal states (computational basis) which can store $`n=\mathrm{log}_2N`$ bits of information typically realized as a collection of $`n`$ 2-level subsystems ($`q`$-bits).

The process of computation is divided into three stages, first an initial (input) state is prepared, then a quantum algorithm is performed which is realized as a sequence of $`𝒢`$ unitary transformation called quantum gates and finally, the output state is measured. The number $`𝒢`$ of involved quantum gates is called sometimes a quantum complexity and one assumes, in analogy to classical digital computers, that the physical time needed to achieve a given task is proportional to $`𝒢`$. If $`𝒢`$ is polynomial in $`n`$ the algorithm is called efficient . The remarkable result of Shor  had an enormous impact on the whole field of quantum information and quantum computing . He constructed a quantum algorithm as a sequence of polynomial in $`\mathrm{log}_2N`$ number of unitary transformations which factorizes numbers smaller than $`N`$ into primes. It is known that the logical complexity $`𝒞`$ of the problem measured by a number of logical steps needed to solve it grows exponentially with $`n=\mathrm{log}_2N`$ for the case of factorization. The possibility of efficient factorization of large numbers by quantum computers must frighten the experts responsible for the safety of information transmission.

The reasoning outlined above possesses several drawbacks. First of all the idea that the physical time of computation $`t_c`$ is proportional to the complexity

$$t_c𝒞\mathrm{or}t_c𝒢$$
$`(1)`$
is only true for the existing digital computers which are ensembles of controlled bistable elements which correspond to Boolean logical values 0,1 and can literary mimic logical operations. Other theoretically conceivable classical computers like, for example, ballistic computer of Fredkin and Toffoli (hard spheres colliding with each other and with fixed reflective barriers) need not satisfy (1). By rescaling masses, distances and initial velocities we can have in principle for a fixed $`N`$ arbitrarily short $`t_c`$ and arbitrarily low energy store needed. Obviously, atomic structure puts limits on the rescaling and classical chaos together with friction make the idea of ballistic computers impractical.

The number of quantum gates $`𝒢`$ which form a quantum algorithm is not unique and does not determine the real physical time of computation. The very idea of quantum gate is not natural in quantum mechanical context. Quantum evolution is continuous in time and governed by a (possibly time-dependent) Hamiltonian which describes the energy of a system. In particular by increasing the energy level spacing of the computer’s Hamiltonian we can speed up its time evolution. Hence the

physical efficiency of computation should be given in terms of the product of a characteristic energy and time (”action parameter”) rather than the computation time alone. A dynamical character and finite duration of the input state preparation procedure and the output state measurement process are also important. As we shall see studying concrete examples there is no clear separation between the preparation-measurement processes and the action of quantum algorithm.

As many authors noticed the decoherence effects due to the interaction with an environment are the main practical obstacles for the operation of quantum computer. As the decoherence typically grows in a nonlinear way with the

energy level spacing this can produce an optimal splitting of the action

parameter into time and energy. This problem is studied in details in  and will be not discussed here.

## II Heisenberg energy-time relations

One of the interpretations of the Heisenberg energy-time uncertainty relation

$$\mathrm{\Delta }E\mathrm{\Delta }t\mathrm{}$$
$`(2)`$
is that a quantum state with spread in energy $`\mathrm{\Delta }E`$ takes time at least $`\mathrm{\Delta }t=\pi \mathrm{}/2\mathrm{\Delta }E`$ to evolve to an orthogonal and hence distinguishable state . Hence, it is quite natural to investigate the

quantum mechanical limitations on the action parameter $`E_ct_c`$ where $`E_c`$ is an energy store needed for the quantum computation and $`t_c`$ is a total computation time including the preparation of an input state and

the measurement of an output. Following  we can propose the inequality of the form

$$E_ct_c>>\mathrm{}𝒞_q$$
$`(3)`$
where $`𝒞_q`$ is a not yet defined ”quantum complexity” of the problem. The ”much larger” symbol in (3) takes into account the probabilistic character of quantum algorithms which should run a certain number of times to achieve a given level of confidence. Although is it not formulated explicitly in

the form of the inequality (3) the general believe among the experts in quantum computation is that the quantum complexity is equal to the minimal number of quantum

gates

$$𝒞_q=𝒢.$$
$`(4)`$
In  the following hypothesis supported by a single example has been proposed

$$𝒞_q𝒞$$
$`(5)`$
what means essentially that there is no quantum  complexity. In particular, if the problem needs an exponential in the input bit size number of logical steps the quantum computer needs also an exponential time or an exponential energy store. It is still better than for the existing digital computers which need both but nevertheless the quantum computation is not practically efficient.

## III Examples

In order to advocate the hypothesis (5) we discuss several examples of quantum algorithms or their essential parts.

### A  State preparation

Consider $`n`$ two level system ($`q`$-bits) with a standard computational basis which consists of products of two $`q`$-bit states $`|0>`$ and $`|1>`$. To prepare an arbitrary input state from this basis starting from a certain fixed initial state we need to rotate some of the 2-level system in a fixed time $`t_c`$ to the orthogonal states. Therefore according to the presented interpretation of (1) we need on the average the energy $`E_c=(n/2)\pi \mathrm{}/2t_c`$ what gives in this case

$$𝒞_qn.$$
$`(6)`$
The similar result can be obtained for the quantum measurement of an output state. As writing or reading of $`n`$-bit messages involves $`n`$ logical steps we have here $`𝒞_q=𝒞`$

### B Grover’s search algorithms

We shall analyze two different Hamiltonian realizations of the Grover’s search algorithm . In the original formulation every element of the database is represented by a state of the standard computational basis of the $`n`$ $`q`$-bit system. One of this states – a searched one– is denoted by $`|x>`$ and we have also a certain standard initial state $`|in>`$ which is usually a uniform superposition of computational basis $`\{|\varphi _j>;j=1,2,\mathrm{},N=2^n\}`$

$$|in>=(N)^{1/2}\underset{j=1}{\overset{N}{}}|\varphi _j>,<s|in>=N^{1/2}.$$
$`(7)`$
These two states are hidden in the dynamics of the system described either by the repeated unitary transformations or the continuous time Hamiltonian evolution. The first Hamiltonian proposed in  has form

$$𝐇_1=E(|x><x|+|in><in|)$$
$`(8)`$
where $`E`$ is an energy scale. The Hamiltonian acts essentially on the two dimensional subspace of the Hilbert space and the corresponding energy difference is equal to

$$\mathrm{}\omega _N=EN^{1/2}.$$
$`(9)`$
One can easily compute that we need a time

$$t_N=\pi N^{1/2}/2E$$
$`(10)`$
to reach from the initial state $`|in>`$ the searched one $`|x>`$. So the quantum complexity of this stage of searching is of the order $`\mathrm{}\omega _Nt_N1`$. Adding preparation and measurement processes we obtain

$$𝒞_qn.$$
$`(11)`$
The second Hamiltonian used with the different rescaling in  reads

$$𝐇_2=iE(|x><in||in><x|).$$
$`(12)`$
Again the problem is essentially two dimensional with the energy difference

$$\mathrm{}\omega _N=2E+o(1/N).$$
$`(12)`$
The time needed to reach the state $`|x>`$ can be estimated by

$$t_N=\pi /4E+o(1/N).$$
$`(13)`$
Here again adding state preparation and measurement we obtain (11).

In the literature on the Grover’s algorithm it is claimed that its quantum complexity is $`\sqrt{N}`$ which is compared with the classical complexity

of the problem claimed to be equal to $`N/2`$. Both statements are incorrect. We have just computed the quantum complexity of the problem equal to $`n=\mathrm{log}_2N`$. The classical analogon of Grover’s search algorithm is not finding an item in a randomly ordered phone book but rather a search for a one heavier ball among $`N`$ otherwise indistinguishable ones. The later problem can be solved in $`\mathrm{log}_2N`$ steps.

### C Quantum telephone directory

We discuss now a true quantum analog of a random telephone directory. It is again a $`n`$ $`q`$-bit system with a computational basis $`\{\varphi _j;j=1,2,\mathrm{},N=2^n\}`$. We fix an initial state to be $`\varphi _1`$ and propose the following time dependent Hamiltonian

$$𝐇(t)=𝐇_0+𝐕\mathrm{cos}(\mathrm{\Omega }t)$$
$`(14)`$
where

$$𝐇_0=\underset{j=1}{\overset{N}{}}E_j|\varphi _j><\varphi _j|$$
$`(15)`$
$`𝐕`$ is a (randomly chosen) weak perturbation and $`\mathrm{\Omega }`$ is a tunable frequency. The energies $`E_j0`$ are not degenerated and provide labels for the states $`\varphi _j`$ (we put $`E_1=0`$). To find a state labeled by $`E_j`$ we tune the frequency to the value $`\mathrm{\Omega }_j=E_j/\mathrm{}`$. The time dependent first-order perturbation calculus  gives us the probability of excitation of the state $`\varphi _k`$

$$p_k(t)=2|<\varphi _1|𝐕|\varphi _k>|^2\frac{\mathrm{sin}^2\left\{\frac{1}{2}(E_kE_j)t/\mathrm{}\right\}}{(E_kE_j)^2}.$$
$`(16)`$
It follows from the formula (16) that we have to wait for a time at least of the order

$$t_j\frac{\mathrm{}}{|E_jE_k|}$$
$`(17)`$
to be sure that the searched state $`\varphi _j`$ has been prepared with a much larger probability than the other neighboring state $`\varphi _k`$. Therefore, on the average we obtain the computation time

$$t_c\mathrm{}N/E_{max}$$
$`(18)`$
where $`E_{max}=\mathrm{max}\{E_j\}`$. Then as $`E_cE_{max}`$ we see again that the quantum complexity coincides with the classical one.

### D  Quantum device factorizing numbers

This model has been introduced in  but we briefly discuss it again for the sake of completeness. In fact this model is very similar to the previous one.

A resonant cavity supports radiation modes with the frequencies being the logarithms of prime numbers times a fixed frequency unit $`\omega `$

$$\omega _q=\omega \mathrm{log}q,q=2,3,5,7,11,13,\mathrm{}$$
$`(19)`$
The second quantization Hamiltonian of the electromagnetic field

$$𝐇=\mathrm{}\omega \underset{q}{}(\mathrm{log}q)a_q^+a_q$$
$`(20)`$
has nondegenerated eigenvalues being proportional to the logarithms of all natural numbers

$$𝐇\psi _N=E_N\psi _N,E_N=\mathrm{}\omega \mathrm{log}N,N=1,2,3,\mathrm{}$$
$`(21)`$
The structure of $`\psi _N`$ reveals the factorization of $`N`$ into prime numbers

$$\psi _N(a_{q_1}^+)^{m_1}(a_{q_2}^+)^{m_2}\mathrm{}(a_{q_r}^+)^{m_r}\psi _1$$
$`(22)`$
where

$$N=(q_1)^{m_1}(q_2)^{m_2}\mathrm{}(q_r)^{m_r}$$
$`(23)`$
and $`\psi _1`$ is a vacuum state. The eq.(22) means that we have $`m_1`$ photons of the frequency $`\omega \mathrm{log}q_1`$, $`m_2`$ photons of the frequency $`\omega \mathrm{log}q_2`$, …, and $`m_r`$ photons of the frequency $`\omega \mathrm{log}q_r`$. Therefore, transferring a given energy portion $`\mathrm{}\omega \mathrm{log}N`$ to the empty cavity and then opening the cavity and counting photons in different modes we obtain the factorization of $`N`$. It can be done similarly to the previous example perturbing the system periodically in time with a tunable frequency $`\mathrm{\Omega }`$ and selecting $`\mathrm{\Omega }=\omega \mathrm{log}N`$.

As the energy level spacing around $`E_N=\mathrm{}\omega \mathrm{log}N`$ is $`\delta E_N\mathrm{}\omega /N`$ it follows from the analogue of the formula (16) that we have to wait for a time at least of the order

$$t_cN\omega ^1$$
$`(24)`$
to select a proper state $`\psi _N`$. The energy used equals $`E_c=\mathrm{}\omega \mathrm{log}N`$ and once again the quantum complexity essentially coincides with the classical one.

### E Shor’s algorithm

The Shor’s factorization algorithm is rather complicated but for our purposes we need only a part of it – the so-called phase shift computation. We follow here the simple presentation in . The phase shift computation consists in applying to the initial state of $`n`$ $`q`$-bit system

$$|in>=2^{n/2}(|0>+|1>)\mathrm{}(|0>+|1>)$$
$`(25)`$
a $`n`$-gates unitary operation

$$𝐔=U^{2^0}U^{2^1}\mathrm{}U^{2^n}$$
$`(26)`$
where

$$U|0>=|0>,U|1>=e^{i\alpha }|1>,\alpha [0,2\pi )$$
$`(27)`$
The unitary $`𝐔`$ can be realized as $`\mathrm{exp}\{i𝐇t_n/\mathrm{}\}`$ with the following Hamiltonian

$$𝐇=\mathrm{}\omega \underset{k=0}{\overset{n1}{}}2^k(|1><1|)_k$$
$`(28)`$
acting for a time $`t_n=\alpha /\omega `$. The averaged energy in the state $`|in>`$ grows exponentially with $`n`$ and is given by

$$\overline{E}_n=<in|𝐇|in>=(\mathrm{}\omega /2)\underset{k=0}{\overset{n1}{}}2^k=\mathrm{}\omega (2^{n1}1/2).$$
$`(29)`$
For the total factorization procedure $`t_c>>t_n`$ and $`E_c>>\overline{E}_n`$. Taking for $`\alpha `$ its average value $`\pi `$ we have $`E_ct_c>>\mathrm{}\mathrm{\hspace{0.17em}2}^n=\mathrm{}N`$ in agreement with the hypothesis (5).

## IV Conclusion

Although the general proof would be very desirable the analysis of the presented examples, in particular the powerful Shor’s algorithm, provides a strong evidence for the Heisenberg-like bound (3)(5) on the efficiency of quantum computations.

###### Acknowledgements.

The author thanks Michał, Paweł and Ryszard Horodecki for discussions. The work is supported by the Grant KBN PB/273/PO3/99/16.


text:
# New type of stable Q balls in the gauge-mediated supersymmetry breaking

\[

## Abstract

We obtain a new type of a stable Q ball in the context of gauge-mediated supersymmetry breaking in minimal supersymmetric standard model. It is so-calledgravity-mediation type of Q ball, but stable against the decay into nucleons, since the energy per unit charge is equal to gravitino mass $`m_{3/2}`$, which can be smaller than nucleon mass in the gauge-mediation mechanism. We consider the cosmological consequences in this new Q-ball scenario, and find that this new type of the Q ball can be considered as the dark matter and the source for the baryon number of the universe simultaneously.

\]

The standard cosmology provides a whole description from a few minutes after “big bang” to now . One of the evidence which supports it is the nucleosynthesis, which successfully predicts cosmological abundances of all light elements. It requires that there is a small asymmetry of the baryons in the universe: $`\eta _B=n_B/n_\gamma 10^{10}`$, where $`n_B`$ and $`n_\gamma `$ denote the number density of the baryon and photon, respectively. This and other observations show that our universe is made almost entirely of matters and devoid of antimatters. Such matter-antimatter asymmetry is produced by baryogenesis, which takes place nonthermally through baryon and CP violating interactions in the very early universe.

On the other hand, inflation solves many problems which cannot be explained in the standard cosmology, such as the homogeneity, flatness, and monopole problems. Inflation thus predicts that $`\mathrm{\Omega }_{\mathrm{tot}}=1`$, where $`\mathrm{\Omega }`$ is the density parameter, the ratio of the density to the critical density $`\rho _c1.9\times 10^{29}h^2`$ g/cm3, and $`h`$ is the Hubble parameter normalized by 100 km/sec/Mpc. However, $`\eta _B10^{10}`$ corresponds to $`\mathrm{\Omega }_Bh^20.02`$, far less than the prediction from the inflation. Even if one does not consider the inflation, $`\mathrm{\Omega }_{\mathrm{matter}}0.2`$ is expected from observations for dynamical properties of galaxies and clusters of galaxies. Therefore, most of the density of the universe has to be in the form of dark matter.

Several mechanisms for creating baryons have been proposed, but none of them explain directly why the universe has almost the same amount of baryons and the dark matter. Their answer is that it is a numerical coincidence. However, the Q ball provides a natural scenario for explaining both of them simultaneously .

In the supersymmetric standard models, Affleck-Dine (AD) mechanism  is the most promising procedure for the baryogenesis. In the minimal supersymmetric standard model (MSSM), there are many flat directions consist of squarks and sleptons , which can be identified as the AD field. Its potential is almost flat but slightly lifted up by effects of supersymmetry (SUSY) breaking. For the mechanism of SUSY breaking, there are two famous scenarios: the gravity- and gauge-mediated SUSY breakings. It was believed that the AD field stays at large field value at the inflationary stage, and, when the Hubble parameter becomes as small as the AD scalar mass after inflation, rolls down homogeneously its potential with rotation, making the baryon number of the universe.

Recently, however, it was revealed that the AD field does not evolve homogeneously, but feels spatial instabilities, which grow nonlinear and form into Q balls . A Q ball is a kind of the nontopological soliton, whose stability is guaranteed by the nonzero charge $`Q`$ . In the context of the AD baryogenesis, the charge $`Q`$ is the baryon number $`B`$. In the gauge-mediated SUSY breaking, a Q ball is stable against the decay into nucleons, provided that its charge is large enough so that its energy per unit charge is less than nucleon mass: $`E_Q/Qm_\varphi Q^{1/4}<1`$ GeV , where $`m_\varphi 1`$ TeV, is the mass of AD field. Therefore, the Q ball itself can be a candidate for the dark matter. On the other hand, in the gravity-mediated SUSY breaking, the energy of a Q ball per unit charge is essentially constant: $`E_Q/Qm_\varphi >`$ 1 GeV . Thus, it should decay into nucleons, and the dark matter will be lightest supersymmetric particles (neutralinos) produced in Q-ball decays. In either case, the dark matter and the baryon number of the universe can be explained simultaneously by the Q-ball formation through the AD mechanism.

In all the previous studies of Q balls in the context of SUSY breaking, the effects of gauge- and gravity-mediations are considered separately. However, it is natural to have both effects in the gauge-mediated SUSY breaking scenario, since the gravity-mediation effects will dominate over the gauge ones at the large field value. Cosmology including AD baryogenesis in such more realistic SUSY breaking scenario was considered in Ref. . There, AD field is regarded as a homogeneously rotating condensate, but we notice that it will form Q balls due to the instabilities of the field. Particular interest is the smallness of the gravitino mass comparing with that in the gravity-mediation scenario. It usually ranges between 100 keV and 1 GeV. Therefore, we can imagine a new type of a stable Q ball: the profile is the same as that in the gravity-mediation, but its energy per unit charge is less that 1 GeV because of the small gravitino mass. In this Letter, we study the cosmological consequences of Q balls (baryogenesis and the dark matter) in the gauge-mediated SUSY breaking, taking into account the gravity-mediation effects at large field value.

To be concrete, let us assume the following potential for the AD field,

$`V(\mathrm{\Phi })`$ $`=`$ $`m_\varphi ^4\mathrm{log}\left(1+{\displaystyle \frac{|\mathrm{\Phi }|^2}{m_\varphi ^2}}\right)`$ (2)
$`+m_{3/2}^2|\mathrm{\Phi }|^2\left[1+K\mathrm{log}\left({\displaystyle \frac{|\mathrm{\Phi }|^2}{M_{}^2}}\right)\right],`$
where $`m_{3/2}`$ is the gravitino mass, $`K(<0)`$ term a one-loop correction, Usually, gaugino contributions to $`K`$ is dominated and $`K`$ becomes negative. However, if $`K>0`$ by some large Yukawa couplings, the AD field can be stabilized and Q balls can be created only for $`\varphi \varphi _{\mathrm{eq}}`$, where $`\varphi _{\mathrm{eq}}`$ is defined in Eq.(3). $`M_{}`$ the renormalization scale, and we assume that the second term should be neglected for small field value. This is nothing but the sum of the potentials for the gauge- and gravity-mediation mechanisms studied previously . However, as we mentioned earlier, the gravitino mass is considerably smaller. The second term will dominate the potential when

$$\varphi \varphi _{\mathrm{eq}}\sqrt{2}\frac{m_\varphi ^2}{m_{3/2}},$$
(3)
where $`\mathrm{\Phi }=\varphi \mathrm{exp}(i\omega t)/\sqrt{2}`$. In this case, a new type of a stable Q ball is produced. Its property is very similar to that in the gravity-mediation, such as

$`R_Q|K|^{1/2}m_{3/2}^1,`$ $`\omega m_{3/2},`$ (4)
$`\varphi |K|^{3/4}m_{3/2}Q^{1/2},`$ $`E_Qm_{3/2}Q,`$ (5)
but, as can be seen from the last equation, it is stable against the decay into nucleons. In the opposite case, the Q-ball properties are the same as in the gauge-mediation only :

$`R_Qm_\varphi ^1Q^{1/4},`$ $`\omega m_\varphi Q^{1/4},`$ (6)
$`\varphi m_\varphi Q^{1/4},`$ $`E_Qm_\varphi Q^{3/4}.`$ (7)
The energy per unit charge can be treated from unified viewpoint. The largest charge of the Q ball formed depends linearly on the initial charge density of the AD field as

$$Q\beta \frac{q(0)}{m_{3/2}^3}\beta \left(\frac{\varphi (0)}{m_{3/2}}\right)^2,$$
(8)
where $`\beta 1`$ , and we use $`q=\omega \varphi ^2m_{3/2}\varphi ^2`$. It can be understood by estimating Q-ball charge in the standard way. The charge is given by

$$Q=d^3x\omega \varphi ^2=\left(\frac{\pi }{2}\right)^{3/2}\omega \varphi _0^2R^3=\beta ^{}\left(\frac{\varphi _0}{m_{3/2}}\right)^2,$$
(9)
where we assume the Gaussian profile ansatz , $`\varphi =\varphi _0\mathrm{exp}(r^2/R^2)`$, which is a very good approximation, and use relations (5), and $`\beta ^{}2\times 10^3(|K|/0.01)^{3/2}`$. The discrepancy of coefficients $`\beta `$ and $`\beta ^{}`$ comes from the fact that $`\varphi (0)\varphi _0`$ and there are more than one Q balls with charges of the same order of magnitude as the largest.

Inserting it into Eq.(3), we obtain

$$m_{3/2}(2\beta )^{1/4}m_\varphi Q^{1/4}.$$
(10)
The right hand side of this equation is identical to the expression for the energy per unit charge of the gauge-mediation besides the factor of order unity. The energy per unit charge of the Q ball is written as

$$\frac{E_Q}{Q}=\{\begin{array}{cc}m_\varphi Q^{1/4}\hfill & \varphi \varphi _{\mathrm{eq}}\hfill \\ m_{3/2}\hfill & \varphi \varphi _{\mathrm{eq}}\hfill \end{array},$$
(11)
which is shown in Fig. 1. The gap on the boundary should disappear and both sides of the curves will be smoothly connected because Q balls formed in this region are not the exact type of either (5) or (7), but will show properties between them.

Since Q balls are stable even for $`\varphi \varphi _{\mathrm{eq}}`$, where the gravity-mediation effect is crucial, the baryon number in the universe should be explained by the baryons evaporated from Q balls, as is the same as for the gauge-mediation type . The evaporation rate of the Q ball is

$$\mathrm{\Gamma }_{\mathrm{evap}}\frac{dQ}{dt}=\alpha \mu T^24\pi R_Q^2,$$
(12)
where $`\mu `$ is a chemical potential of the Q ball, which is estimated as $`\mu \omega `$ because $`\omega `$ is energy of $`\varphi `$-field inside the Q ball. Although the mass of the (free) AD particle $`m_\varphi `$ is affected by thermal corrections, which should be changed as $`m_\varphi m_\varphi (T)T`$, at $`Tm_\varphi `$, the gravitino mass is not affected, since particles coupled to the AD field are decoupled from thermal bath when the AD field has a large vacuum expectation value. At $`Tm_\varphi `$, large numbers of $`\varphi `$-particles are in thermal bath outside Q balls, so $`\alpha 1`$. On the other hand, since only light quarks are in thermal bath at $`Tm_\varphi `$, the corresponding cross section is highly suppressed by the heavy gluino exchanges, and $`\alpha (T/m_\varphi )^2`$.

However, if the rate of the charge diffusion from the “atmosphere” of the Q ball is smaller than the evaporation rate, chemical equilibrium will established there, which results in the suppression of the evaporation . The diffusion rate is

$$\mathrm{\Gamma }_{\mathrm{diff}}\frac{dQ}{dt}=4\pi \zeta R_QD\mu _QT^2,$$
(13)
where $`D=a/T`$ with $`a46`$, and $`\zeta 1`$.

The time scale of charge transportation is determined by the diffusion when $`\mathrm{\Gamma }_{\mathrm{diff}}\mathrm{\Gamma }_{\mathrm{evap}}`$. It holds for $`TT_{}a^{1/3}|K|^{1/6}(m_{3/2}m_\varphi ^2)^{1/3}`$. In this case, using Eqs.(5) and (13), and assuming the radiation-dominated universe, $`t=AM/T^2`$, where $`A0.2`$ and $`M2.4\times 10^{18}`$ GeV, we obtain

$$\frac{dQ}{dT}=\frac{8\pi aAM}{|K|^{1/2}T^2}.$$
(14)
On the other hand, when $`TT_{}`$, the diffusion effect is negligible, and Eq.(14) should be replaced by

$$\frac{dQ}{dT}=\frac{8\pi AMT}{|K|m_{3/2}m_\varphi ^2}.$$
(15)
Therefore, total amount of the charge evaporated from the Q ball is

$$\mathrm{\Delta }Q12\pi AM\left(\frac{a}{|K|}\right)^{2/3}(m_{3/2}m_\varphi ^2)^{1/3}.$$
(16)
Provided that the initial charge of the Q ball is larger than the evaporated charge, we regard that the Q ball survives from evaporation, and contributes to the dark matter of the universe:

$$Q_{\mathrm{init}}\mathrm{\Delta }Q9.8\times 10^{18}\left(\frac{m_{3/2}}{\mathrm{GeV}}\right)^{1/3}\left(\frac{m_\varphi }{\mathrm{TeV}}\right)^{2/3},$$
(17)
where we set $`a=4`$ and $`|K|=0.01`$.

Now we can relate the baryon number and the amount of the dark matter in the universe. As mentioned above, the baryon number of the universe should be explained by the amount of the charge evaporated from Q balls, $`\mathrm{\Delta }Q`$, and the survived Q balls become the dark matter. If we assume that Q balls do not exceed the critical density of the universe, i.e., $`\mathrm{\Omega }_Q1`$, and the baryon-to-photon ratio as $`\eta _B10^{10}`$,

$$Q3.2h^2\times 10^{21}\left(\frac{m_{3/2}}{\mathrm{GeV}}\right)^{4/3}\left(\frac{m_\varphi }{\mathrm{TeV}}\right)^{2/3}.$$
(18)
Rewriting Eq.(10), we have

$$Q2\beta \times 10^{12}\left(\frac{m_\varphi }{\mathrm{TeV}}\right)^4\left(\frac{m_{3/2}}{\mathrm{GeV}}\right)^4.$$
(19)
Combining this constraint with Eqs.(17) and (18), together with $`m_{3/2}1`$ GeV, which implies that the gravity-mediation type of the Q ball is stable against the decay into nucleons, we obtain the allowed region for the new type of the stable Q ball explaining the baryon number of the universe. Figure 2 shows the allowed region on $`(Q,m_{3/2})`$ plane for $`m_\varphi =1`$ TeV. The shaded regions represent that the new type of stable Q balls are created, and the baryon number of the universe can be explained by the mechanism mentioned above. Furthermore, the new type of stable Q balls contribute crucially to the dark matter of the universe at present, if the Q balls have the charge given by the thick line in the figure. Notice that Q balls with very large charge are not allowed because they will overclose the density of the universe.  If the baryons are produced by other mechanism, larger Q balls can be allowed. In this case, however, a nice relation between the densities of the baryon and the dark matter does not hold. Therefore, the initial conditions for the AD field is restricted severely .

One may wonder if these new type of stable Q balls can be detected. When a Q ball collides with nucleons, they enter the surface layer of the Q ball, and dissociate into quarks, which are converted into squarks. In this process, Q balls release $`1`$ GeV energy per collision by emitting soft pions. This process is the basis for the Q-ball detections , which is called (Kusenko-Kuzmin-Shaposhnikov-Tinyakov) KKST process in the literature. It occurs for electrically neutral Q balls (ENQB). For electrically positively charged Q balls (EPCQB), the KKST process is strongly suppressed by Coulomb repulsion, and only electromagnetic processes will take place. For electrically negatively charged Q balls (ENCQB), the both KKST and electromagnetic processes occur, but the former is dominant, which is essentially the same as for ENQBs.

In either case, the detection is more difficult than for the gauge-mediation type of Q balls, since the geometrical cross section is smaller for large $`Q`$, and the Q-ball mass is larger for the same $`Q`$, which results in small flux. With the discussions similar to Ref. , we can restrict the parameter space $`(Q,m_{3/2})`$ by the several experiments. Fig. 3 shows the results for ENQBs. Lower left regions are excluded by the various experiments. The present available data prohibit the stable gravity-mediation type Q balls with large gravitino mass to be both the dark matter and the source of the baryons, and future experiments such as the Telescope Array Project or the OWL-AIRWATCH detector may detect this type of Q balls with an interesting gravitino mass $`100`$ keV.

For EPCQBs with $`Z=1`$, similar constraints as for ENQBs are put by the MACRO experiment .

In summary, we have obtained a new type of a stable Q ball in the context of gauge-mediated SUSY breaking in MSSM. Many properties are the same as the gravity-mediation type of Q ball, but it is stable against the decay into nucleons, since the energy per unit charge is equal to the gravitino mass $`m_{3/2}`$, which can be smaller than nucleon mass of 1 GeV in the gauge-mediation mechanism. We have considered the cosmological consequences in this new Q-ball scenario. Because of its stability, it can be a nice candidate for the dark matter of the universe. In the present case, the baryons are produced only by evaporation from Q balls, since (almost) all the baryons are trapped in Q balls during their formation. We have found that the Q ball with $`Q10^{25}10^{22}`$ can account for both the dark matter and the baryon number of the present universe for $`m_{3/2}10^410^1`$ GeV and $`m_\varphi =1`$ TeV, and such Q balls may be detected by the future experiments.

The authors are grateful to R. Banerjee for useful comments and T. Yoshida for helpful discussions. M.K. is supported in part by the Grant-in-Aid, Priority Area “Supersymmetry and Unified Theory of Elementary Particles”($`\mathrm{\#}707`$).


text:
# Elastic properties of model porous ceramics 11footnote 1Submitted to the Journal of the American Ceramic Society

## 1 Introduction

The elastic properties of two-phase (solid-pore) porous materials depend on the geometrical nature of the pore space and solid phase, as well as the value of porosity . Relevant aspects of porous materials may include pore shape and size as well as the size and type of the interconnections between solid regions. These features, which generally lack precise definition, comprise the microstructure of the material. In order to predict properties, or properly interpret experimental property-porosity relationships, it is necessary to have an accurate method of relating elastic properties to porosity and microstructure. In this paper we use the finite element method to derive simple formulae that relate Young’s modulus and Poisson’s ratio to porosity and microstructure, for three different models of microstructure.

There have been several different approaches to deriving property-porosity relations for porous materials. Formulae derived using the micro-mechanics method  are essentially various methods of approximately extending exact results for small fractions of spherical or ellipsoidal pores to higher porosities. This includes the differential  and self consistent methods  as well as the commonly used semi-empirical correction to the dilute result made by Coble and Kingery  to explain the properties of porous alumina. A drawback of this approach is that the microstructure corresponding to a particular formula is not precisely known; hence agreement or disagreement with data can neither confirm nor reject a particular model. A second problem is that these types of models provide no predictions for the case where the microstructure is comprised of incompletely sintered grains, which is a common morphology in porous ceramics. A second class of results  have been termed minimum solid area (MSA) models. In this approach purely geometrical reasoning is used to predict the elastic moduli based on the weakest points within the structure. Again, the microstructure that corresponds to the MSA predictions is not exactly known. A number of semi-empirical relations have also been proposed , which generally provide a reasonable means of describing data, extrapolating results and comparing data among materials. However, lacking a rigorous connection with microstructure, these results do not offer either predictive or interpretive power. Theoretical bounds  exist for the elastic properties, but the vanishing of the lower bound for porous materials lessens their predictive power when the upper bound does not provide a good estimate. There are numerous other approaches, including the generalized method of cells .

Another approach is to computationally solve the equations of elasticity for digital models of microstructure . In principle this can be done exactly. However, large statistical variations and insufficient resolution, have limited the accuracy of results obtained to date. Only recently have computers been able to handle the large three-dimensional models and number of computations needed to obtain reasonable results. As input to the method, we employ three different microstructural models that broadly cover the types of morphology observed in porous ceramics. The models are based on randomly placed spherical pores, solid spheres, and ellipsoidal pores . The centers of the pores or solid particles are un-correlated which leads to realistic microstructures in which both the pore and solid phase are interconnected. The results, which can be expressed simply by two (or sometimes three) parameter relations, correspond to a particular microstructure and explicitly show how the properties depend on the nature of the porosity. Therefore, the results can be used as a predictive tool for cases where the microstructure of the ceramic is similar to one of the models, and as an interpretive tool if the microstructure is unknown. The numerically exact FEM results are compared with various well-known micro-mechanics and MSA results to determine how close an approximation a particular formula provides for each model. In the FEM, we can freely vary the properties of the solid phase, allowing us to determine the dependence of Young’s modulus and Poisson’s ratio on the solid Poisson’s ratio as well as on the porosity. This question has attracted recent interest in the ceramics literature .

II. Computational Results

A microstructure made up of a digital image is already naturally discretized and so lends itself to numerical computation of many quantities. The finite element method uses a variational formulation of the linear elastic equations, and finds the solution by minimizing the elastic energy via a fast conjugate gradient method. The digital image is assumed to have periodic boundary conditions. Details of the theory and programs used are reported in the papers of Garboczi & Day  and Garboczi .

In order to obtain accurate results using the FEM on models of random porous materials, it is absolutely necessary to estimate and minimize three sources of error: finite size effects, discretization errors, and statistical fluctuations. This has generally not been done in the past, owing to limitations in computer memory and speed. FEM results for random microstructures do not have much meaning without such an error analysis.

The various sources of error are defined in the following way. First, the length scale of the microstructure is fixed, usually by fixing the size of a typical pore (e.g. the spherical pore radius). The size of the system is then controlled by the side length of the cubic sample, denoted $`T`$. The size of $`T`$ compared to the pore size controls how many pores will appear in the computational cell. A real material has many thousands or more such pores. Errors can occur in using a smaller number in a periodic cell to simulate a much larger number. We vary $`T`$ in order to map out this effect.

Once a value of $`T`$ is chosen that minimizes finite size errors but is still computationally possible, we next must consider the discretization error, which comes about because we are using discrete pixels to represent continuum objects. The number of pixels along each edge of the cubical unit cell is $`M`$, giving a resolution of $`dx`$=$`T/M`$ (in units of $`\mu `$m per pixel, if $`T`$ is in $`\mu `$m). For the chosen value of $`T`$, a value of $`M`$ is chosen that also gives acceptable discretization errors, usually on the order of a few percent.

Finally, when computing the properties of random materials, either computationally or experimentally, one must carefully choose the number of samples ($`N_s`$) over which the results need to be averaged to produce acceptable uncertainties. This value is again chosen, within computational constraints, to keep statistical fluctuations within a few percent.

Overlapping Solid Spheres. Realizations of the overlapping solid sphere model  are generated by placing solid spheres at random points in the unit cell. This produces a set of overlapping grains that mimic the microstructure of sintered ceramic composites (see Fig. 1a). The space outside the solid grains is the pore space, with porosity $`\varphi `$. The pore phase is macroscopically connected above porosities of $`\varphi 0.03`$ and the solid phase remains connected for values of $`\varphi `$ below $`0.70`$ . Above $`\varphi `$=0.7, the solid phase is composed of isolated solid particles. So between $`\varphi =0.03`$ and $`\varphi =0.70`$, the overlapping solid sphere model is bi-continuous. In ceramics the porosity is generally less then $`0.40`$, in this bi-continuous regime. We therefore consider the elastic properties for $`0.1\varphi 0.50`$, where the solid Poisson’s ratio, $`\nu _s`$, varied over the range $`0.1\nu _s0.4`$.

To generate the microstructure we chose solid spheres of radius $`r`$=1 $`\mu `$m. Note that the elastic properties are length scale invariant so the results apply to spheres of any radius for which the continuum assumption holds. A preliminary study showed that finite size errors were acceptably small for cubic samples with edge length $`T`$=12 $`\mu `$m. To study the discretization errors we generated one realization of the model with porosity $`\varphi `$=0.5 at $`M`$= 48$`\mathrm{}`$128. The elastic properties depend quite strongly on resolution. We found that the variation of Young’s modulus with $`M`$ could be described by the relation

$$E_{\mathrm{FEM}}(M)E_0+aM^1$$
(1)
where $`E_0`$ can be identified as the continuum value (corresponding to infinitely large $`M`$). The same is true for Poisson’s ratio. Even at $`M=128`$ the finite element code overestimates the ‘exact’ result for the Young’s modulus by 30 %. Therefore, for the overlapping sphere model it is necessary to measure the elastic moduli at three different values of $`M`$ and extrapolate the results to $`M\mathrm{}`$. We chose $`N_s`$=5 samples at each resolution and porosity, except at $`\varphi `$=0.5 where large statistical variations implied a larger number of samples was necessary ($`N_s`$=10). Thus 30 different realizations of the models were considered, each at 3 different discretizations, for a total of $`90`$ models.

The statistical variation in Young’s modulus and Poisson’s ratio for the case $`\nu _s`$=0.2 are shown in Table 1. The error bars shown in the table are equal to twice the standard error (S.E.=$`\sigma /\sqrt{N_s}`$ with $`\sigma `$ the standard deviation). Therefore there is a 95% chance that the “true” result lies between the indicated error bars. The results are accurate to within 20% at $`\varphi `$=0.5; the error decreasing with porosity to less than 10% for $`\varphi 0.30`$. The expected Gaussian distribution of the measured averages implies that the results are actually more accurate than this. For example, the anticipated relative errors are halved if a 68% likelihood threshold is used (i.e., $`\pm `$ one standard error).

In addition to the above results we also computed the elastic moduli of the 90 model microstructures at solid Poisson’s ratios $`\nu _s`$=0.1, 0.3 and 0.4. The statistical variation was not significantly different from the case $`\nu _s`$=0.2. Combined with the data for $`\nu _s`$=0.2 this covers most commonly occuring solids. The scaled Young’s modulus for each value of $`\nu _s`$ is plotted against porosity in Fig. 2. Remarkably, the scaled Young’s modulus of the porous material appears to be practically independent of $`\nu _s`$. This result has been proven to be exact in 2-D  and appears to hold to a very good approximation in 3-D. We found that the Young’s modulus data are well described by an equation of the form

$$\frac{E}{E_s}=\left(1\frac{\varphi }{\varphi }_0\right)^n$$
(2)
with $`n`$=2.23 and $`\varphi _0`$=0.652 and $`0\varphi 0.5`$. Note that $`n`$ and $`\varphi _0`$ are empirical correlation parameters and should not be interpreted as the percolation exponent and threshold, respectively. Percolation concepts are generally valid closer to the threshold $`\varphi _c0.7`$ (for this model) and a higher value of $`n`$ is expected. The computational cost of accurately measuring the elastic properties increases greatly as the percolation threshold is approached.

The Poisson’s ratio of the porous material is shown in Fig. 3 as a function of $`\varphi `$ and $`\nu _s`$, and appears to be a flow diagram , where the Poisson’s ratio asymptotically approaches a fixed point, independently of the value of the solid Poisson’s ratio. This flow diagram has been analytically proven to hold in 2-D, when a percolation threshold exists at which the Young’s modulus goes to zero . This flow diagram also appears to be valid in 3-D as well, within numerical uncertainty. The Poisson’s ratio data shown in Fig. 3 can be roughly described by the simple linear relation,

$$\nu =\nu _s+\frac{\varphi }{\varphi }_0(\nu _0\nu _s)=\nu _0+\left(1\frac{\varphi }{\varphi }_0\right)(\nu _s\nu _0)$$
(3)
with two fitting parameters $`\nu _0`$=0.140 and $`\varphi _0`$=0.472. A more accurate fit is obtained with the three parameter relation,

$$\nu =\nu _0+\left(1\frac{\varphi }{\varphi }_0\right)^m(\nu _s\nu _0).$$
(4)
with $`\nu _0`$=0.140, $`\varphi _0`$=0.500 and $`m`$=1.22.

Overlapping Spherical Pores. The overlapping spherical pore (or swiss cheese) model  is generated by interchanging the roles of the solid and pore phase of the overlapping solid sphere model (see Fig. 1b). The morphology corresponds to isolated spherical pores at low porosity, with the pores becoming macroscopically interconnected at $`\varphi 0.3`$. The solid phase remains connected up to $`\varphi 0.97`$. This type of morphology may arise in ceramics generated with a particulate filler , or where bubbles form in a molten state . We consider solid Poisson’s ratios in the range $`0.1\nu _s0.4`$.

We determined that statistical errors were acceptable for a computational cube of size $`T`$=12 $`\mu `$m with pores of radii $`r`$=1 $`\mu `$m. Using $`M`$=80 pixels, the discretization errors were less than 3 % for $`\varphi =0.5`$ and 2 % for $`\varphi =0.3`$. Therefore, it was not considered necessary to generate samples at different discretizations ($`M`$) and extrapolate the results. As for solid spheres, the Young’s modulus was independent of the solid Poisson’s ratio to a very good approximation. The Young’s modulus can be described by Eq. (2) with $`n=1.65`$ and $`\varphi _0`$=0.818 (Fig. 2). Poisson’s ratio of the porous material is shown in Fig. 4 and is simply described by the linear relation given in Eq. (3) with $`\nu _0`$=0.221 and $`\varphi _0`$=0.840. Again, a flow diagram is observed.

Overlapping Ellipsoidal Pores. A common method of analyzing the effect of pore shape on elastic properties is to study ellipsoidal pores. In analytic formulae, it is possible to treat the limiting cases of needles and platelets, although the difficulty of resolving these fine structures prohibits these limits from being treated with the finite element method. However, the percolation properties of these limiting cases can be computationally studied . To gauge the effect of deviations from spherical shaped pores we considered isotropically oriented overlapping oblate ellipsoidal pores bounded by the surface $`(x/a)^2+(y/b)^2+(z/c)^2=1`$ with $`a`$=$`b`$=1 $`\mu `$m and $`c`$=0.25 $`\mu `$m (see Fig. 1c). For this case, the pore phase becomes connected at porosity $`\varphi =0.2`$ . Statistical errors were found to be acceptable for a computational cube of size $`T`$=10 $`\mu `$m. Using $`M`$=96 pixels the discretization errors were 3 % for $`\varphi =0.5`$ and 2 % for $`\varphi =0.3`$. As for the case of spherical pores, these errors were considered sufficiently small, so that the added computational burden of the extrapolation technique could be again avoided.

Again the Young’s modulus was found to be independent of the solid Poisson’s ratio to a very good approximation. The results, shown in Fig. 2, can be accurately described by Eq. (2) with $`n`$=2.25 and $`\varphi _0`$=0.798. The Poisson’s ratio (Fig. 5) can be roughly fit using the Eq. (3) with $`\nu _0`$=0.166 and $`\varphi _0`$=0.604. A better fit is obtained using Eq. (4) with $`m`$=1.91, $`\nu _0`$=0.161 and $`\varphi _0`$=0.959. A flow diagram similar to that seen before is obtained.

The CPU time and memory required for these computations are an important “experimental” detail. The memory requirement for a given model was $`230\times M^3`$ bytes, where $`M`$ was the edge length in pixels of a cubic unit cell. So for the largest computations carried out, $`M=128`$, the memory requirement was about $`500`$ Mbytes. The amount of CPU time consumed was approximately $`3000`$ hours, divided among different modern workstations.

III. Comparison with micro-mechanical and MSA formulae

In this section we compare a selection of well known theoretical results with the ‘numerically exact’ data computed in the previous section. These results include both analytically exact results (bounds, expansions, dilute limits, composite sphere assemblage), as well as approximate results (effective medium theories, minimum solid area models).

There are several kinds of exact bounds that have been derived for elastic materials . These are equations involving the different phase moduli, the volume fractions of the different phases, and various correlation functions that define the geometry of the composite. The upper bound gives the maximum possible composite elastic moduli, and the lower bound gives the lowest possible composite elastic moduli. The bounds used in this paper are three- point bounds, which have been written out explicitly for overlapping solid spheres and overlapping spherical pores . In the case where one phase has zero elastic moduli, as is true in this paper, the lower bound becomes zero as well, and so only the upper bound is meaningful.

An exact perturbation expansion also exists, where the elastic moduli of a two-phase material is expanded in terms of parameters involving the individual elastic moduli of each phase and geometrical quantities . This expansion has been carried out to three terms explicitly, and it is this truncated form to which we will compare our numerical data. The result is expected to be accurate when the void phase is not interconnected.

Another exact result, which is used later in this section to build the various effective medium theories, is the case of dilute spherical pores for which the exact effective moduli are given by

$`K`$ $`=`$ $`K_m+c_iP^{mi}(K_iK_m)`$ (5)
$`G`$ $`=`$ $`G_m+c_iQ^{mi}(G_iG_m)`$ (6)
with

$`P^{mi}`$ $`=`$ $`{\displaystyle \frac{3K_m+4G_m}{3K_i+4G_m}},Q^{mi}={\displaystyle \frac{G_m+F_m}{G_i+F_m}},`$
$`F_m`$ $`=`$ $`{\displaystyle \frac{G_m}{6}}{\displaystyle \frac{9K_m+8G_m}{K_m+2G_m}}.`$ (7)
Here $`c_i`$ denotes the concentration (volume fraction) of inclusions and the subscripts $`i`$ and $`m`$ on the bulk $`K`$ and shear modulus $`G`$ denote the properties of the inclusion and matrix, respectively. The result is attributed to numerous authors . For a porous matrix $`K_i=G_i=0`$ and the porosity is $`\varphi =c_i`$. The result is strictly valid for small concentrations of inclusions $`\varphi 1`$ (in practice $`\varphi <0.1`$.) Cast in terms of the engineering constants for porous inclusions this result becomes

$`E`$ $`=`$ $`E_m{\displaystyle \frac{3}{2}}\varphi E_m{\displaystyle \frac{94\nu _m5\nu _m^2}{75\nu _m}}+O(\varphi ^2)`$ (8)
$`\nu `$ $`=`$ $`\nu _m{\displaystyle \frac{3}{2}}\varphi {\displaystyle \frac{(5\nu _m1)(1\nu _m^2)}{75\nu _m}}+O(\varphi ^2)`$ (9)
Our prior statement of Eq. (8) inadvertently omitted the factor of $`3/2`$, although the correct result was used in the paper. A non zero quadratic term can be added (as an empirical correction) to ensure that $`E`$=0 at $`\varphi `$=1. This was suggested by Coble and Kingery  for MacKenzie’s  result for spherical pores, which is equivalent to Eqs. (5)-(7) with $`K_i=G_i=0`$.

To adapt the dilute formulas to the case of finite porosity a number of proposals have been made. The approximate equations that result are usually called effective medium theories. The most common approximation is the so-called self consistent method (SCM) of Hill  and Budiansky . In this model the equations of elasticity are solved for a spherical inclusion embedded in a medium of unknown effective moduli. The effective moduli $`K`$ and $`G`$ are then derived. In the dilute case the embedding medium is just the matrix. The Hill-Budiansky result can be stated as

$`c_iP^i(K_iK_{})+c_mP^m(K_mK_{})`$ $`=`$ $`0`$ (10)
$`c_iQ^i(G_iG_{})+c_mQ^m(G_mG_{})`$ $`=`$ $`0`$ (11)
where $`K_{}`$ and $`G_{}`$ denote the effective moduli and $`P^m`$ and $`Q^m`$ are given in Eq. (7). The equations cannot be explicitly solved and numerical methods are necessary (see Hill  and Berryman  for details). In the case of porous inclusions, the moduli vanish at $`\varphi =\frac{1}{2}`$, which is a property not shared with most composites (e.g., the overlapping sphere model). To derive a more realistic result, Christensen and Lo  generalized the SCM (GSCM) to the case of a spherical shell embedded in a matrix of unknown moduli. The result is complicated and not reproduced here.

The differential method (reviewed by McLaughlin ) provides an alternative model using a similar philosophy. Suppose that the effective moduli of a composite medium are known to be $`K_{}`$ and $`G_{}`$. Now if a small additional concentration of inclusions are added, the change in $`K_{}`$ and $`G_{}`$ is approximated to be that which would arise if a dilute concentration of inclusions were added to a uniform, homogeneous matrix with moduli $`K_{}`$ and $`G_{}`$. This leads to a pair of coupled differential equations,

$`{\displaystyle \frac{dK_{}}{dc_i}}`$ $`=`$ $`P^i{\displaystyle \frac{K_iK_{}}{1c_i}};K_{}(c_i=0)=K_m`$ (12)
$`{\displaystyle \frac{dG_{}}{dc_i}}`$ $`=`$ $`Q^i{\displaystyle \frac{G_iG_{}}{1c_i}};G_{}(c_i=0)=G_m.`$ (13)
The dilute result, the self consistent result , and the differential method  can all be extended to the case of spheroidal inclusions. The general results  for $`P^{mi}`$ and $`Q^{mi}`$ have been given by Berryman . In addition to these results Wu  derived a variant of the self consistent method, where $`K_{}`$ and $`G_{}`$, the effective moduli, are found by implicitly solving the equations

$`K_{}`$ $`=`$ $`K_m+c_iP^i(K_iK_m)`$ (14)
$`G_{}`$ $`=`$ $`G_m+c_iQ^i(G_iG_m).`$ (15)
A different type of microstructure is provided by Hashin’s  model of space-filling poly-disperse hollow spheres (the “composite-sphere assemblage”). Although a simple formula exists for the bulk modulus over the full porosity range , exact results for the Young’s moduli are not available. Ramakrishnan and Arunachalam  recently derived the approximation

$`{\displaystyle \frac{E}{E_s}}`$ $`=`$ $`{\displaystyle \frac{(1\varphi )^2}{(1+2\varphi 3\nu _s\varphi )}}`$ (16)
$`\nu `$ $`=`$ $`{\displaystyle \frac{(4\nu _s+3\varphi 7\nu _s\varphi )}{4(1+2\varphi 3\nu _s\varphi )}}.`$ (17)
However, the derivation is not rigorous. In particular, the exact result for the bulk modulus of the model  is around twice that predicted by Eqs. (1617) at $`\varphi =0.5`$. Since Eq. (16) was found to provide reasonable agreement with experimental data for porous ceramics , we compare its predictions to our FEM data below.

The final class of results we consider is provided by the ‘minimum solid area’ (MSA) models  (which have been recently reviewed by Rice ). This approach is based on the assumption that the ratio of the effective moduli to the solid moduli is directly proportional to the minimum ratio of solid contact area to the total cross-sectional area of periodic structures. The approximation derived depends on the particular model considered. We consider two basic models most closely aligned with our FEM data: simple cubic arrays of solid and porous spheres. The latter case provides a particularly simple example of the type of result which can be derived. Suppose the repeat distance of the lattice is $`2h`$ and the sphere radius is $`r`$. The Young’s modulus is assumed to be proportional to the area fraction, giving

$`{\displaystyle \frac{E}{E_s}}={\displaystyle \frac{(2h)^2\pi r^2}{(2h)^2}}=1{\displaystyle \frac{\pi }{4}}\left({\displaystyle \frac{6}{\pi }}\right)^{\frac{2}{3}}\varphi ^{\frac{2}{3}}`$ (18)
since $`\varphi =\frac{1}{6}\pi (r/h)^3`$. The form of the result changes for $`r>h`$ (or $`\varphi >\pi /6`$=0.52) as the spheres begin to coalesce. Rice  has noted the results of many different periodic structures can be approximated by the form $`E/E_s=e^{b\varphi }`$ over a range of porosities. For example, $`b5`$ for the solid sphere model and $`b3`$ for the porous sphere model. It is argued that for a given set of data, $`b`$ can be compared with known values to assess the type of porosity. Often fractions of different types of porosity are assumed to match experimental data making the method an interpretive rather than a predictive tool. Since we have measured $`E`$ for microstructures based on solid sphere contacts and porous spheres we should be able to ascertain the accuracy of the MSA formulae for these cases.

Fig. 6 shows the comparison between the exact three-point bounds  for the overlapping solid sphere and spherical pore case, the truncated expansion  for the overlapping spherical pore case, and the numerical results. Clearly the expansion does better than the three-point bound for the overlapping spherical pore case, though both are fairly close to the numerical results. The bound lies far away from the overlapping solid sphere numerical results, however. For this case, the truncated expansion does not exist. Only the $`\nu _s=0.2`$ data is shown. Using the truncated expansion, one can show that in 3-D, the Young’s modulus is not exactly independent of the solid Poisson’s ratio, but is rather a very good approximation, as was shown earlier in this paper.

In Fig. 7, we compare the FEM data ($`\nu _s=0.2`$) for overlapping spherical pores with dilute and effective medium theory analytic results. At this Poisson’s ratio the SCM and dilute results reduce to $`E/E_s=12\varphi `$ while the differential and dilute results with the Coble–Kingery correction reduce to $`E/E_s=(1\varphi )^2`$. Since the analytic results are based on the case of dilute spherical pores they all match the FEM data at $`\varphi =0.1`$. The dilute and SCM results under-estimate the FEM data at higher porosities because of the aphysical percolation threshold at $`\varphi =\frac{1}{2}`$. The generalized SCM over-estimates the data, while the differential method performs reasonably well over the entire porosity range. The latter observation might have been anticipated given the close association between the definition of the spherical pore model and the assumption of the differential method. At increasing porosities we are simply adding additional spherical pores to a porous matrix. The data for overlapping solid spheres is also shown in the figure, and seen to be quite different from any of the available results. This demonstrates that microstructure (the geometrical nature of the porosity) is an important factor besides the actual value of the porosity.

In Fig. 8, the minimum solid area models and the Ramakrishnan and Arunachalam results  are compared with the data. The MSA model for spherical pores performs reasonably well, although underestimating the FEM data for overlapping random spherical pores at low porosities $`\varphi <0.3`$. The MSA model for solid spheres considerably underestimates these data for $`\varphi <0.3`$. The Ramakrishnan and Arunachalam  approximation falls between the FEM data for $`\varphi >0.1`$ indicating that it corresponds to neither of the microstructures. For purposes of comparison we also report numerical results obtained using the computational generalized method of cells . For a periodic spherical pore the results significantly underestimate the FEM data for overlapping spherical pores at low porosities (and hence the exact dilute result). It is not clear if this is due to the assumptions, or the particular implementation, of the method.

The FEM data for overlapping oblate ellipsoidal pores is compared with the available theories in Fig. 9. The SCM results of Wu  and Berryman  underestimate the porosity as a result of underestimating the physical percolation threshold. The Berryman result performs significantly better than does the Wu result. As for the case of spheres, the differential method matches the data quite closely because of the similarity between the assumptions of the theory and the definition of the model.

We have also compared the Poisson’s ratio predicted by the various self-consistent and differential methods to the FEM data for overlapping spherical and ellipsoidal pores. The theoretical results converge to different fixed points (e.g. Fig. 3) in qualitative agreement with the data. But only the differential method provides reasonable agreement with the FEM data (with absolute error less than 0.02 for $`\varphi 0.4`$ and $`0.1\nu _s0.4`$).

IV. Comparison with Experiment

We now use the FEM results to analyze experimental measurements of the elastic properties of porous ceramic materials. The dependence of the elastic moduli on porosity has been the subject of many studies . Data for porous alumina from numerous studies  are shown in Fig. 10. The Coble-Kingery  material is markedly stiffer than other materials, and is in very good agreement with the FEM results for the overlapping spherical pore model. The pores in the alumina matrix were actually created by the incorporation of a particulate filler , which corresponds well with the definition of the model microstructure. The remaining data closely follow the overlapping solid sphere FEM result for $`\varphi <0.25`$, indicating that the solid alumina phase has the sintered granular morphology exhibited by the model microstructure (Fig.1a). However, Knudsen notes that several of the samples summarized were also created using particulate fillers. At higher porosities the solid sphere result underestimates the data. One reason for this might be that the model contains isolated solid spheres which artificially reduce the actual porosity. This was checked and found not to be the case for the porosities studied. Therefore, the solid connections in these samples of porous alumina are likely stiffer than those found in the solid sphere model at porosities $`\varphi >0.25`$. Overlapping spheres can create very sharp “valleys” between a pair of overlapping solid spheres (see Fig. 1a), which would be rounded off in the sintering process, presumably strengthening the solid-solid connection.

Hunter et al.  have studied the Young’s modulus of several different oxides. In all cases, the porous material was created by sintering a powder of the pure oxide. The results for the Young’s modulus are reproduced in Fig. 11. For low porosities ($`\varphi <0.1`$) all of the data followed the FEM results for overlapping spherical pores. For Gd<sub>2</sub>O<sub>3</sub> the FEM result continues to provide excellent agreement up to the maximum porosity measured ($`\varphi =0.4`$) indicating that the microstructure is similar to that of the model (overlapping pores). In contrast, the data for the other three oxides decreases towards the result for overlapping solid spheres indicating a more granular character.

The data of Walsh et al.  for porous glass is compared with the FEM results for overlapping spherical pores in Fig. 12. The agreement is good for small to moderate porosities ($`\varphi <0.3`$), but the FEM results underestimate the data at higher porosities. Walsh et al. point out that the pores in the glass are actually not interconnected (unlike the overlapping pores of the model). This would account for the increased stiffness. It is interesting that the FEM results begin to deviate from the experimental data at the threshold where the pores become macroscopically connected ($`\varphi `$=0.3). Data for sintered MgAl<sub>2</sub>O<sub>2</sub>  powder is shown in Fig. 13, and is well modeled by the FEM results for overlapping solid spheres. Micrographs of the ceramic indicate a granular structure similar to that of the model microstructure (although the grains appear more like polyhedra, not spheres).

V. Discussion and Conclusions

We have derived empirical theories for the dependence of the Young’s modulus on porosity for three distinct models of porous ceramics, based on careful finite element computations. An advantage of these results over many conventional theories is that they correspond to a priori known microstructures. The dilute result (extended by Coble-Kingery to all porosities), the differential method, and the self consistent methods all have a “built-in” microstructure, but apart from the dilute case, it is not clear what that microstructure is. Therefore, agreement or disagreement with a particular analytic theory neither confirms nor rejects a particular physical model. For the minimum solid area models, the microstructure is exactly known, but the approximation involved in making the Young’s modulus directly proportional to the contact area leads to a similar conclusion. Indeed, we found that the MSA models did not provide quantitative agreement with the moduli of the random microstructures studied. We found that the differential method (Eqs. 12 and 13) gave results in reasonable agreement with computed data for the cases of overlapping spherical and ellipsoidal pores, probably due to the similarities between the assumptions of the model and the definitions of the microstructure. Results for the granular model of overlapping solid spheres were not well modeled by any of the analytic theories, demonstrating the importance of finite element techniques in this case of great physical interest.

We have also generated data that shows the dependence of Poisson’s ratio on porosity and the solid Poisson’s ratio. It is difficult to study this question experimentally because of the inability to vary the Poisson’s ratio of the solid independently, and the well known difficulties of accurately measuring the Poisson’s ratio at moderate to high porosities . At sufficiently high porosities we find that the Poisson’s ratio converges to a fixed non-zero value $`\nu _0`$ irrespective of the solid Poisson’s ratio. For overlapping solid spheres $`\nu _0=0.14`$, spherical pores $`\nu _0=0.22`$ and oblate ellipsoidal pores $`\nu _0=0.16`$. This behavior is exact in two dimensions  and is exhibited by many of the analytic theories in three dimensions. At present the available experimental data cannot confirm this qualitative behavior . We have shown that the Poisson’s ratio does not vanish at high porosities as has been recently argued .

It is not simple to attribute our results to features of the solid-pore morphology such as the size, shape, distribution and connectivity of pores or solid grains, since these features have no obvious definition for complex bi-continuous random microstructures. A few general observations can be made, and interpreted in terms of interrelated geometrical and mechanical features of the models. For a given porosity, the sintered grain structure of the overlapping solid sphere model is relatively weak. The small solid contacts between spheres and the highly interconnected porosity (which becomes macroscopically connected at $`\varphi `$=0.03) lead to a weak structure. We also assume that the valleys which occur between grains will provide sites of large stress concentrations, and consequently, large deformations. In contrast, spherical pores provide high (near optimal) stiffness at a given porosity. The dispersed nature of the porosity (which is macroscopically disconnected for $`\varphi <0.3`$) corresponds to a well connected solid matrix. Ellipsoidal pores tend to weaken a structure more than spherical pores due to a combination of a less well connected solid phase (the pores become macroscopically connected at $`\varphi =0.2`$), and greater stresses and deformations near the high curvature regions of the ellipsoid.

We have compared our FEM results with several sets of previously published experimental data. In cases where the microstructure of the porous ceramics roughly matched that of the models, the agreement was very good. Since the FEM results correspond to a known microstructure, it was possible to explain deviations in terms of specific microstructural features. Thus, comparison of experimental data with the three computational results provides a useful interpretive tool. Note that a given elastic modulus does not correspond to a particular microstructure. Therefore, it is important to corroborate microstructural interpretations obtained from the elastic moduli with information about the particular material (such as a micrograph). In the future it would be useful to extend this work to higher porosities and to other relevant models (such as non-overlapping porous spheres). It is also possible to use statistical microstructural information obtained from two-dimensional micrographs to generate models  that actually mimic physical microstructures.

Acknowledgments A.R. thanks the Fulbright Foundation and Australian Research Council for financial support. We also thank the Partnership for High-Performance Concrete Technology program of the National Institute of Standards and Technology for partial support of this work.