text:
# Hollow sphere, a flexible multimode Gravitational Wave antenna

## 1 Introduction

A Gravitational Wave (GW) is a time–dependent perturbation of the background space-time geometry which propagates away from the source, and which is generated by changes in the mass–energy distribution of the latter. Due to the weakness of gravitation, the mathematical description of GWs, so far as detection by a remote observatory is concerned, can be very accurately accomplished by means of a linear approximation in the field equations. More precisely, if $`\eta _{\mu \nu }`$ is the (Lorentzian) background metric, the actual metric $`g_{\mu \nu }(𝐱,t)`$ can be expediently split up as

$$g_{\mu \nu }(𝐱,t)=\eta _{\mu \nu }+h_{\mu \nu }(𝐱,t),|h_{\mu \nu }|1$$
(1)
where x and $`t`$ are quasi–Lorentzian coordinates, and the $`h_{\mu \nu }(𝐱,t)`$ are the GW perturbations.

Since GWs will of course be observed far from their source they can be safely assumed to be plane across the detector’s span. In this case only the six so called electric components, $`R_{0i0j}(𝐱,t)`$, of the Riemann tensor associated with (1) are not identically zero. Although General Relativity theory predicts that at least four of these are zero, as explained in textbooks , alternative theories of gravity (e.g. Brans–Dicke) are generally less restrictive —see  for a complete discussion of all the possibilities. It thus appears that detection of GWs is a strong experimental test on which is the theory which best describes gravity.

Currently operating Weber bars, as well as interferometric detectors, such as LIGO and VIRGO, are however single channel devices, i.e., they can only sense one linear combination of the 5 quadrupole GW amplitudes —the so called antenna pattern . A network of such detectors in joint operation is thus necessary for a complete deconvolution of the GW signal, but this poses non-trivial practical problems of analysis, because there are no two identical devices, noise characteristics change from one to the other, etc.

A simple way out of such limitation is to use a spherical detector, because this is a natural multimode, or tensor, GW antenna . The reason for this is that, having spherical symmetry, an elastic solid has oscillation eigenmodes which come in degenerate $`l`$-pole harmonics, the monopole ($`l=0`$) and quadrupole ones efficiently coupling to the incoming GW . By use of suitable sensor systems  it is theoretically possible to completely deconvolve the GW signal, i.e., to measure all six amplitudes which characterise a general metric GW .

A hollow sphere is an interesting variation of the solid one  which has a number of interesting capabilities as a GW detector, and flexibility to implement them by adjusting the thickness parameter. In this paper I intend to summarise the most relevant features of an ideal hollow spherical GW antenna (i.e., I will assume that it is technically possible to build a working system), and see what interesting physics could be done with it.

## 2 Response and sensitivity

As usual (see  for a thorough discussion), the solid’s response to a GW excitation is given by the solution, with appropriate boundary conditions, to the partial differential equation

$$\varrho \frac{^2𝐮}{t^2}\mu ^2𝐮(\lambda +\mu )(𝐮)=𝐟(𝐱,t)$$
(2)
where $`\varrho `$, $`\lambda `$ and $`\mu `$ are the density and elastic Lamé coefficients, respectively, and

$$𝐟(𝐱,t)=𝐟^{(00)}(𝐱)g^{(00)}(t)+\underset{m=2}{\overset{2}{}}𝐟^{(2m)}(𝐱)g^{(2m)}(t)$$
(3)
is the tidal GW driving density of force. The hollow sphere’s response is then given by

$$𝐮(𝐱,t)=\underset{n=1}{\overset{\mathrm{}}{}}\frac{a_{n0}}{\omega _{n0}}𝐮_{n00}(𝐱)g_{n0}^{(00)}(t)+\underset{n=1}{\overset{\mathrm{}}{}}\frac{a_{n2}}{\omega _{n2}}\left[\underset{m=2}{\overset{2}{}}𝐮_{n2m}(𝐱)g_{n2}^{(2m)}(t)\right]$$
(4)
where $`\omega _{nl}`$ are the oscillation eigenfrequencies of the $`𝐮_{nlm}(𝐱)`$ eigenmode, $`a_{nl}`$ are overlap coefficients and $`g_{nl}^{(lm)}(t)`$ is the convolution between the signal $`g^{(lm)}(t)`$ and the detector mode. Eigenmode wave functions $`𝐮_{nlm}`$ are somewhat complicated functions which are calculated in detail in reference , while the frequencies of the first two quadrupole and first monopole harmonics are plotted in figure 1 for a conceivable materials and dimensions of the system.

The most relevant parameter characterising the solid’s sensitivity to GWs is the cross section for the absorption of incoming GW energy. As shown in reference , the following formula holds for quadrupole modes:

$$\sigma _n\sigma _{\mathrm{abs}}(\omega _{n2})=F_n(\sigma _P;a/R)\frac{GMv_s^2}{c^3}$$
(5)
where $`v_s`$ is the speed of sound in the material, $`\sigma _P`$ its Poisson ratio, and $`F_n(\sigma _P;a/R)`$ is a dimensionless parameter associated with the $`n`$-th quadrupole mode of the hollow sphere. A plot of the first two $`F_n`$’s is given in figure 2. The crossing between the curves happens at $`a/R=0.375`$, which means a hollow sphere is a more efficient GW energy absorber in its second mode than in the first, which favours in turn higher frequency signal reception.

## 3 Physics with a hollow sphere

Let us now briefly assess which sources of GWs can be visible with a hollow spherical detector, under the ideal assumption that technological problems are solved to the effect that the detector’s quantum limit has been reached. If such is the case then peak sensitivity is given by

$$S_h(\omega _n)=\frac{G\mathrm{}}{c^3}\frac{4\pi \beta _n}{\sigma _n}$$
(6)
where $`\omega _n`$ is a quadrupole resonance frequency, and $`\beta _n`$ is the efficiency of the coupling between mechanics and electronics in that mode; the detector’s bandwidth is $`\mathrm{\Delta }\omega _n=\beta _n\omega _n`$.

Take for example a hollow sphere of CuAl, 3 metres of outer diameter and 22 centimetres thick ($`a/R=0.15`$), weighing 40 tons; the first two quadrupole frequencies are then 400 and 1500 Hz, respectively, and the peak sensitivity is essentially identical in both modes, reaching $`10^{46}`$ Hz-1 for a $`\beta `$ factor of 0.1. A system like this would be able to see a millisecond burst with an amplitude of $`10^{22}`$ metres/metre, and a monochromatic signal as weak as $`10^{27}`$ metres/metre at the first resonant frequency if integrated for one year. Such sensitivities would most likely enable the observations of supernova explosions in galaxies out to the Virgo cluster (at a rate of a few per year), and a significant sample of the galactic pulsar population.

A larger detector, 6 metres in outer diameter, 25 centimetres thick, 200 tons of CuAl alloy, has a minimum spectral density of noise of $`3\times 10^{47}`$ Hz-1 at its first resonance (190 Hz), and $`2\times 10^{47}`$ Hz-1 at the second (750 Hz). Such a system can be used for chirp signal detection by a robust method first described for spherical detectors by Coccia and Fafone . The idea is conceptually simple: the GW signal generated e.g. by a compact binary system nearing coalescence sweeps an ascending frequency range which eventually overlaps with the sensitivity bands of the detector. If the binary remains undisrupted  while in the frequency intervals between the first two just quoted resonances then the antenna’s response will peak at those frequencies, with a signal to noise ratio given by

$$SNR=\frac{5\pi 2^{1/3}}{24}\frac{(GM_c)^{5/3}}{c^3}\frac{\mathrm{\Delta }\omega _n}{S_h(\omega _n)}\frac{1}{r^2\omega _n^{7/3}}$$
(7)
where $`r`$ is the distance to the source, $`\omega _n`$ is a detector quadrupole resonance frequency, and $`M_c`$ is the mass parameter

$$M_c(m_1m_2)^{3/5}(m_1+m_2)^{1/5}$$
(8)
of the chirp signal. Let now $`\tau _1`$ and $`\tau _2`$ be the times at which the the signal reaches the detector with frequencies $`\omega _1`$ and $`\omega _2`$, respectively. The mass parameter can then be estimated by the formula

$$M_c=\frac{5^{3/5}}{2^{16/5}}\frac{c^3}{G}\left(\frac{\omega _2^{8/3}\omega _1^{8/3}}{\tau _2\tau _1}\right)^{3/5}$$
(9)
in which the lowest post-Newtonian approximation for the waveform has been assumed . The strong point of this double passage method is its robustness relative to the details of the signal’s phase , which considerably complicates the performance of the matched filter in extracting the signal from the noise in an interferometric detector . It thus appears that both detectors, a sphere and an interferometer, can constructively cooperate in the identification of compact binary coalescences, as the sphere can provide a good estimate of the mass parameter which will simplify and enhance the efficiency of the bank of filters needed in an interferometer .

The hollow sphere is particularly suited for this, as it can be made to have low resonance frequencies, and this favours SNR because the signal stays longer in the lower frequencies —gravitational bremsstrahlung tends to accelerate the inspiral as the binary’s period of revolution diminishes. This fact is clearly reflected by equation (7), where SNR is seen to strongly depend on the 7/3 inverse power of frequency.

Cosmic backgrounds of GWs can also be observed with hollow spheres. Such signals, due to their stochastic nature, are impossible to tell from the detector noise. A minimum of two detectors, whose outputs are cross-correlated, must be used in order to see a common signal which can be identified by a distinctive spectral density, $`S_{\mathrm{GW}}(\omega )`$. With unit SNR, a pair of identical hollow spheres could set a good limit on the density parameter of GWs in the Universe :

$$\mathrm{\Omega }_{\mathrm{GW}}10^9\times \left(\frac{f_n}{200\mathrm{Hz}}\right)^3\left(\frac{S_h^{(1)}(f_n)S_h^{(2)}(f_n)}{10^{96}\mathrm{Hz}^2}\right)^{1/2}\left(\frac{20\mathrm{Hz}}{\mathrm{\Delta }f_n}\right)^{1/2}\left(\frac{1\mathrm{year}}{\tau }\right)^{1/2}$$
(10)
where $`\tau `$ is the integration time and $`\mathrm{\Delta }f_n`$ the detector’s bandwidth around $`f_n`$.

It has recently been conjectured by Amelino-Camelia , on the basis of dimensional arguments, that a cosmic background of GWs could have, in first approximation, a flat spectrum with spectral density

$$S_{\mathrm{GW}}(\omega )\frac{\mathrm{}_{\mathrm{Planck}}}{c}=\left(\frac{G\mathrm{}}{c^5}\right)^{1/2}5\times 10^{44}\mathrm{Hz}^1$$
(11)
It is unclear whether this is a realistic assumption ; if it were however it would mean that Planck scale scale physics would become accessible to the next generation of GW detectors, and hollow spheres would enable the observation of such phenomena with signal to noise ratios of the order of about 100.

## 4 Closing remarks

As already mentioned, spherical elastic bodies are natural multimode GW antennas. The reason for this is the degeneracy of their quadrupole oscillation eigenmodes, which ideally match the structure of GWs. Both a solid and a hollow sphere share this property, of course, and in this paper I have briefly considered some of the more interesting properties of the latter, related to its sensitivity and subsequent capabilities as a powerful GW telescope. I have deliberately left aside the issues of how it could be cast and suspended in practice, as well as questions of how a quantum limited transducer system can be coupled efficiently to the body of the antenna. The paper is thus intended to give a consistent flavour of the best performance one can possibly aim at with such a detector.

There being the possibility to adjust the spherical shell’s thickness, the hollow sphere is a more flexible device than a solid one to meet given specifications. In fact, recent research work  shows that it is advantageous to use a hollow and a solid sphere placed in a joint concentric layout to produce a GW detector which has enhanced sensitivity over a frequency band between the resonances of the inner solid body and the outer hollow one. By an optimised choice of the latter’s thickness, and by means of optical Fabry-Perot cavities acting as non-resonant transducers, such a nested sphere system appears to be a promising device for GW signal detection in the kHz region in bandwidths of several hundred Hz, where it beats the prospects of the very sensitive earth based interferometric GW detectors. And, let me stress this once more, all of it with the ability to potentially deconvolve the five quadrupole amplitudes with a single system.

## Acknowledgements

I acknowledge with thanks support from the Spanish Ministry of Education, contract code BFM2000-0604.

## References


format: text

text:
# A Fast, Accurate and Robust Algorithm For Transferring Radiation in Three-Dimensional Space

## 1 Introduction

Radiative transfer in three dimensional space is a seven-dimensional problem: three spatial dimensions, two angular dimensions, one frequency dimension plus the time dimension. As a result, although the basic physics involved is well understood, a direct computation of radiative transfer in three dimensional space is prohibitively costly.

However, radiation field is known to play a very important role in determining the ionizational and thermodynamic state of the cosmic gas. Rudimentary treatments of radiative transfer from assuming a uniform radiation field (e.g., Cen & Ostriker 1993) to using the somewhat improved method with local optical depth approximation (Gnedin & Ostriker 1997; Cen & Ostriker 1999) have afforded the breakthrough computations of the optically thin regions of the Ly$`\alpha `$ forest with notable successes (Cen et al. 1994; Zhang et al. 1995; Miralda-Escudé et al. 1996; Hernquist et al. 1996). Nevertheless, fluctuations in the ionizing radiation field are expected to exist even in the optically thin regions, since the radiation sources (galaxies and quasars) as well as density fields are known to be significantly clustered at all times. Turning away from optical thin regions, one is faced with transitional regions of optical depth of order unity and optical thick regions largely shielded from external radiation. Moreover, radiation sources themselves most likely reside in dense, at least partially shielded regions often with complicated structures and geometries. For all these regions a proper treatment of radiative transfer is clearly demanded. The predictive power of any theory of galaxy and structure formation, especially with respect to the cosmological reionization, the formation and evolution of Ly$`\alpha `$ forest, damped Ly$`\alpha `$ systems and galaxies/stars, is fundamentally limited until one can accurately treat this essential process of radiative transfer.

Significant progress has been made in recent years in developing practical algorithms for radiative transfer in cosmological applications. Most of the effort has largely been concentrated on two types of approaches: the ray tracing approach (Abel et al. 1999; Rozoumov & Scott 1999; Kessel-Deynet & Burkert 2000) and moment equations approach (Stone, Mihalas & Norman 1992 for two dimensional space; Norman, Paschos, & Abel 1998; Gnedin & Abel 2001). Several Monte Carlo methods have also been explored (Sokasian et al. 2001; Ciardi et al. 2001). The primary limitation of the ray tracing algorithm is its intrinsic high cost to follow a large number of radiation sources. The moment equations approach condenses the intrinsic difficulty of treating long range radiation propagation out as the Eddington tensor. The problem then boils down to evaluating the Eddington tensor accurately and efficiently. The recent suggestion by Gnedin & Abel (2001) of computing the Eddington tensor by ignoring optical depth is novel and clearly worth being explored further.

In this paper we develop an entirely new algorithm for transferring radiation in three dimensional space. In essence, this algorithm computes the optical depth between any pair of points and the source function for each point using convolution techniques. The angular discretization is performed at the receiving site, rather than at the source site as in the normal ray tracing scheme, which guarantees coverage of all space regardless of the fineness of the angular discretization. The overall scaling of the cost with this algorithm is $`N(\mathrm{log}N)^2`$. Unlike conventional ray casting schemes, the computational cost with the present algorithm is independent of the number of sources, suitable for problems encountered in cosmological simulations where a large number of direct, ionizing sources (e.g., small galaxies) as well as secondary, processed ionizing sources (e.g., scattered photons, recombination photons) are present. This paper is organized as follows. In §2 we describe the basic algorithm. In §3 we present a battery of tests to show that the present algorithm is accurate. Discussion and conclusions are given in §4.

## 2 A New Algorithm for Radiative Transfer in Three-Dimensional Space

In the rest frame of a bundle of photons the equation of radiation transfer for those photons (Spitzer 1980) is:

$$\frac{dI_\nu }{ds}=\kappa _\nu I_\nu +j_\nu ,$$
(1)
where $`s`$ is the path length and $`I_\nu (\stackrel{}{x},\stackrel{}{n},t)`$ is the specific intensity with $`I_\nu d\nu d\omega dAdt`$ being the energy during a time interval $`dt`$ passing through an area $`dA`$ about a spatial point $`\stackrel{}{x}`$, within a frequency interval $`d\nu `$, within the solid angle $`d\omega `$ about $`\stackrel{}{n}`$; $`\kappa _\nu (\stackrel{}{x},t)`$ and $`j_\nu (\stackrel{}{x},t)`$ are the opacity and emissivity at $`\stackrel{}{x}`$ at time $`t`$. We can integrate equation (1) and write it in an integral form for the specific flux $`F_\nu `$ defined as

$$F_\nu (\stackrel{}{x},t)I_\nu (\stackrel{}{x},\stackrel{}{n},t)𝑑\omega $$
(2)
at position $`\stackrel{}{x}`$ and time $`t`$:

$$F_\nu (\stackrel{}{x},t)=\frac{1}{4\pi }j_\nu (\stackrel{}{x}+s\stackrel{}{n},t)e^{\tau _\nu (\stackrel{}{x},s\stackrel{}{n},t)}𝑑s𝑑\omega ,$$
(3)
where $`j(\stackrel{}{x}+s\stackrel{}{n},t)`$ is the emissivity at distance $`s`$ from position $`\stackrel{}{x}`$ in the direction $`\stackrel{}{n}`$ (an unit vector) at time $`t`$; the two integrals are over path length and solid angle, respectively; $`\tau _\nu (\stackrel{}{x},s\stackrel{}{n},t)`$ is the optical depth from position $`\stackrel{}{x}`$ to position $`\stackrel{}{x}+s\stackrel{}{n}`$ at time $`t`$:

$$\tau _\nu (\stackrel{}{x},s\stackrel{}{n},t)=_0^r\kappa _\nu (\stackrel{}{x},s^{}\stackrel{}{n},t)𝑑s^{},$$
(4)
where $`\kappa (\stackrel{}{x},s\stackrel{}{n},t)`$ is the opacity at distance $`s`$ from position $`\stackrel{}{x}`$ in the direction $`\stackrel{}{n}`$ at time $`t`$.

The implicit principal assumption that has been made to derive Equation (3) is that the speed of light is infinity, analogous to the conventional treatment of gravitational interactions. Although this assumption may not be necessary, it makes implementation of this method especially simple and we adopt it. This assumption is generally excellent for cosmological simulations, where light crossing time over a simulation box is indeed much smaller than the typical hydrodynamic time step. It should, however, be noted that this assumption does not impose any practical restrictions on the propagation speed of ionization fronts.

Equation (3) may be discretized in the following form:

$$F_\nu (\stackrel{}{x},t)=\underset{s}{}\underset{\stackrel{}{n}}{}S_\nu (\stackrel{}{x},s\stackrel{}{n},t)e^{\tau _\nu (\stackrel{}{x},s\stackrel{}{n},t)}\left(\frac{1}{4\pi s^2}\right),$$
(5)
where the two sums are over the ray path (assuming to be a straight line) and the solid angle, respectively, and the source term $`S_\nu `$ is defined as

$$S_\nu (\stackrel{}{x},s\stackrel{}{n},t)\overline{j}_\nu (\stackrel{}{x},s\stackrel{}{n},t)\mathrm{\Delta }V(\stackrel{}{x},s\stackrel{}{n})$$
(6)
with $`\mathrm{\Delta }V(\stackrel{}{x},s\stackrel{}{n})`$ being the discretization volume element about position $`\stackrel{}{x}+s\stackrel{}{n}`$ and $`\overline{j}_\nu `$ is the mean emissivity over $`\mathrm{\Delta }V(\stackrel{}{x},s\stackrel{}{n})`$. The problem now translates to the evaluation of the two terms on the right hand side in Equation (5), $`S_\nu (\stackrel{}{x},s\stackrel{}{n},t)`$ and $`e^{\tau _\nu (\stackrel{}{x},s\stackrel{}{n},t)}`$. This is the core of the overall computation and how it is computed determines the effectiveness and accuracy of the algorithm. We propose to evaluate these two terms in a novel way.

To proceed further we will now choose the spherical coordinate system for the radiation field about $`\stackrel{}{x}`$, under which we will discretize the angular and radial dimensions in a spherically symmetric fashion. This is to say, one can always find a pair of equal volume elements (of arbitrary domain shape) one about point $`\stackrel{}{x}`$ and the other point $`\stackrel{}{x}`$ related by

$$\mathrm{\Delta }V(\stackrel{}{x},r\stackrel{}{n})=\mathrm{\Delta }V(\stackrel{}{x},r\stackrel{}{n}),$$
(7)
where $`r`$ is the radial distance from $`\stackrel{}{x}`$ along a radial direction $`\stackrel{}{n}`$ or $`\stackrel{}{n}`$. For every vector $`y=r\stackrel{}{n}`$ that lies in $`\mathrm{\Delta }V(\stackrel{}{x},r\stackrel{}{n})`$, one can find an opposite vector (about $`\stackrel{}{x}`$) $`\stackrel{}{y}=r(\stackrel{}{n})`$ that lies in $`\mathrm{\Delta }V(\stackrel{}{x},r\stackrel{}{n})`$. This property will be essential in our implementation of the computation of the two terms in Equation (5). Let us rewrite the source term $`S_\nu `$ in Equation (6) as

$$S_\nu (\stackrel{}{x},r\stackrel{}{n},t)=_{\mathrm{\Delta }V}j_\nu (\stackrel{}{x}+r\stackrel{}{n},t)𝑑V,$$
(8)
where the integration domain is the volume element $`\mathrm{\Delta }V(\stackrel{}{x},r\stackrel{}{n})r^2drd\omega (r,\stackrel{}{n})`$ about position $`\stackrel{}{x}+r\stackrel{}{n}`$. It is now time to make a simple but algorithmically essential modification of Equation (8) by inserting a window function in the integral:

$$S_\nu (\stackrel{}{x},r\stackrel{}{n},t)=j_\nu (\stackrel{}{x}+r\stackrel{}{n},t)g(\stackrel{}{r})d^3\stackrel{}{r},$$
(9)
where $`g(\stackrel{}{r})`$ is a window function about the origin $`\stackrel{}{x}`$ with the following property:

$`g(\stackrel{}{r})`$ $`=`$ $`1\mathrm{for}\mathrm{vectors}\mathrm{within}\mathrm{\Delta }V`$ (10)
$`=`$ $`0\mathrm{otherwise}.`$
Note that the integral in Equation (9) is now over the entire space about $`\stackrel{}{x}`$. The reader is then asked to make the following critical observation: for each window function $`g(\stackrel{}{r})`$ there is another window function $`f(\stackrel{}{r})`$, which is spherically symmetrically paired with it, i.e.,

$$f(\stackrel{}{r})=g(\stackrel{}{r}),$$
(11)
as indicated by Equation (7). Replacing $`g(\stackrel{}{r})`$ with $`f(\stackrel{}{r})`$ in Equation (9) gives,

$$S_\nu (\stackrel{}{x},r\stackrel{}{n},t)=j_\nu (\stackrel{}{x}+r\stackrel{}{n},t)f(\stackrel{}{r})d^3\stackrel{}{r}.$$
(12)
Changing integration variable from $`\stackrel{}{r}`$ to $`\stackrel{}{y}\stackrel{}{x}+\stackrel{}{r}`$ Equation (12) becomes

$$S_\nu (\stackrel{}{x},r\stackrel{}{n},t)=j_\nu (\stackrel{}{y},t)f(\stackrel{}{x}\stackrel{}{y})d^3\stackrel{}{y}.$$
(13)
The integration domain is still over the entire space but the window function $`f(\stackrel{}{x}\stackrel{}{y})`$ now has a moving origin $`\stackrel{}{x}`$. Equation (13) is of the standard convolution form and can be efficiently evaluated for all points $`\stackrel{}{x}`$ simultaneously using the Fast Fourier Transform (FFT) techniques, as the Fourier transforms are related by

$$\underset{¯}{S}_\stackrel{}{k}=\underset{¯}{j}_\stackrel{}{k}\underset{¯}{f}_\stackrel{}{k},$$
(14)
where $`\underset{¯}{S}_\stackrel{}{k}`$, $`\underset{¯}{j}_\stackrel{}{k}`$ and $`\underset{¯}{f}_\stackrel{}{k}`$ are Fourier transforms of $`S_\nu `$, $`j_\nu `$ and $`f`$, respectively.

The same formulation can be constructed for the optical depth $`\tau _\nu (\stackrel{}{x},r\stackrel{}{n},t)`$. Here we instead compute the mean opacity of each volume element

$`\overline{\kappa }_\nu (\stackrel{}{x},r\stackrel{}{n},t)`$ $`=`$ $`{\displaystyle \frac{_{\mathrm{\Delta }V}\kappa _\nu (\stackrel{}{x}+r\stackrel{}{n},t)𝑑V}{\mathrm{\Delta }V}}`$ (15)
$`=`$ $`{\displaystyle \frac{\kappa _\nu (\stackrel{}{y},t)f(\stackrel{}{x}\stackrel{}{y})d^3\stackrel{}{y}}{\mathrm{\Delta }V}}`$
using the same technique. Then we integrate $`\overline{\kappa }_\nu `$ radially along $`\stackrel{}{n}`$ in real space to obtain $`\tau _\nu (\stackrel{}{x},r\stackrel{}{n},t)`$.

$$\tau _\nu (\stackrel{}{x},r\stackrel{}{n},t)=_0^r\overline{\kappa }_\nu (\stackrel{}{x},r^{}\stackrel{}{n},t)𝑑r^{}.$$
(16)
Using the FFT method, $`S_\nu (\stackrel{}{x},r\stackrel{}{n},t)`$ (Equation 13) and $`\overline{\kappa }_\nu (\stackrel{}{x},r\stackrel{}{n},t)`$ (Equation 15) for all grid points of a simulation box with respect to all other points in the simulation box can be evaluated simultaneously with the total cost scaling as $`N(\mathrm{log}N)m_rm_a`$, where $`N`$ is the total number of grid points, $`m_r`$ is the number of radial bins along each cone and $`m_a`$ is the number of angular discretization elements.

The choice of $`m_r`$ and how to discretize the radial direction may be problem dependent. It is, however, illuminating to note that both radiation flux and gravitational force follow the same inverse square law. It is also noted that the effect due to optical depth attenuation tends to relatively suppress the contribution from distant sources. Thus, one may conservatively apply some of the techniques and ideas used in gravitational tree codes (Barnes & Hut 1986). Specifically, one can more coarsely sample more distant regions to discretize the radial direction, the simplest of which would result in $`m_r\mathrm{log}N`$. The choice of $`m_a`$ requires some experiments. We note that, in the present algorithm, any coarse level of angular discretization guarantees that all sources inside the simulation box are covered. This property is in contrast with the conventional ray casting scheme, where a large number of rays from a source is necessary in order to intersect every cell in the simulation box at least once. The primary operational difference is that in the present algorithm the operation is “gathering” flux by the receiving region under consideration, whereas in the ray casting method the operation is “broadcasting” flux to each receiving region. Testing shows that $`m_a=16\times 16=256`$ to cover the full solid angle $`4\pi `$ appears adequate. Combining the three scaling factors gives the total number of operations for evaluating $`S_\nu `$ or $`\tau _\nu `$ for all grid points at each time step:

$$\mathrm{Total}\mathrm{number}\mathrm{of}\mathrm{operations}N_{op}=A_{\mathrm{FFT}}N(\mathrm{log}N)^2m_a,$$
(17)
where $`A_{\mathrm{FFT}}`$ is the normal prefactor of the FFT operation. Note that the algorithm presented here has the desired property that its computational cost is independent of the number of radiation sources. As a matter of fact, each grid point is treated as both a source and sink. This property allows us to compute problems where a large number of radiation sources are present. Examples include reionization of the universe by sub-galactic galaxies and radiation field where diffuse sources such as recombination photons and scattered photons contribute significantly.

We will now implement the algorithm in a three-dimensional simulation box with a uniform mesh for hydrodynamic variables and periodic boundary conditions for both hydro quantities and radiation field. Two separate coordinate systems and two correspondingly separate grids are in use now: the Cartesian coordinate system for the hydrodynamic quantities with a uniform Cartesian grid and the spherical coordinate system for the radiation field about each hydro grid point with a spherical grid. The angular discretization about each hydro grid point is made uniformly across the solid angle $`4\pi `$ with $`m_a`$ elements each having a full solid angle of $`4\pi /m_a`$. The radial discretization is performed using $`m_r`$ bins from $`r=0`$ to $`r_{max}=\sqrt{3}L`$ spaced uniformly in logarithm maximizing the shape of each volume element as a cube, where $`L`$ is the simulation box size. The $`m_rm_a`$ window functions (i.e., $`f`$ in Equations 13,15) for each discretized element of the radiation grid in the spherical coordinates are computed once initially and their Fourier transforms are stored and used for all subsequent time steps.

So far we have operated in a static coordinate system. In cosmological applications one needs to take into account the cosmological effects. We will separate the overall radiation field into two parts: 1) radiation flux due to sources within the simulation box, and 2) external radiation flux. As long as the simulation box is much smaller than the Hubble radius, this division allows us to split these two kinds of sources cleanly in the simple fashion: radiation flux from (1) is treated as local, which is not affected by cosmic expansion, and radiation flux from (2) is diffuse and subject to the cosmological effects. The two kinds of radiation fields are related and computed in a self-consistent way. We keep track of radiation that leaves the simulation box and add it to the diffuse background (see below) and every point in the simulation box is subject to the diffuse background with self-consistently computed attenuation. Thus, the total flux at any point inside the simulation box is

$`F_{\nu ,tot}(\stackrel{}{x},t)=F_{\nu ,in}(\stackrel{}{x},t)+F_{\nu ,ext}(\stackrel{}{x},t),`$ (18)
where $`F_{\nu ,in}(\stackrel{}{x},t)`$ is the contribution from sources inside the simulation box in Equation (3) and $`F_{\nu ,ext}(\stackrel{}{x},t)`$ is the contribution from the diffuse background related to the specific intensity of the diffuse background at time $`t`$, $`I_{\nu ,diff}(t)`$:

$`F_{\nu ,ext}(\stackrel{}{x},t)`$ $`=`$ $`I_{\nu ,diff}(t){\displaystyle e^{\tau _\nu (\stackrel{}{x},r_{upp}\stackrel{}{n},t)}𝑑\omega (\stackrel{}{n})}`$ (19)
$`=`$ $`I_{\nu ,diff}(t){\displaystyle \frac{4\pi }{m_a}}{\displaystyle \underset{l=1}{\overset{m_a}{}}}e^{\tau _{\nu ,l}(\stackrel{}{x},r_{upp}\stackrel{}{n},t)},`$
where the integral or the sum is over the solid angle about position $`\stackrel{}{x}`$; $`\tau _\nu (\stackrel{}{x},r_{upp}\stackrel{}{n},t)`$ is total optical depth along the cone about $`\stackrel{}{n}`$ from position $`\stackrel{}{x}`$ to a distance of $`r_{upp}`$ at time $`t`$ ($`r_{upp}`$ is the upper radius of that cone given the constraint that the simulation box is periodic; $`r_{upp}`$ is orientation-dependent because the hydro simulation box is not spherical about each grid point for the radiation field; the range for $`r_{upp}`$ is $`Lr_{upp}r_{max}=\sqrt{3}L`$). The evolution of the specific intensity of the diffuse background, $`I_{\nu ,diff}(t)`$, is governed (Cen & Ostriker 1993) by

$$\frac{I_{\nu ,diff}(t)}{t}=H(t)\left[\nu \frac{I_{\nu ,diff}(t)}{\nu }3I_{\nu ,diff}(t)\right]+cs_{\nu ,in}(t)+cs_{\nu ,ext}(t)c\overline{\kappa }_{\nu ,diff}(t)I_{\nu ,diff}(t),$$
(20)
where $`H(t)`$ is the Hubble constant, $`c`$ is the speed of light and the mean opacity $`\overline{\kappa }_{\nu ,diff}(t)`$ at time $`t`$ is determined self-consistently by averaging the opacities in the simulation box properly taking into account optical shielding effect:

$`\overline{\kappa }_{\nu ,diff}(t)`$ $`=`$ $`{\displaystyle \frac{1}{4\pi N}}{\displaystyle \underset{\stackrel{}{x}}{}}\kappa _\nu (\stackrel{}{x}){\displaystyle e^{\tau _\nu (\stackrel{}{x},r_{upp}\stackrel{}{n},t)}𝑑\omega (\stackrel{}{n})}`$ (21)
$`=`$ $`{\displaystyle \frac{1}{Nm_a}}{\displaystyle \underset{\stackrel{}{x}}{}}\kappa _\nu (\stackrel{}{x}){\displaystyle \underset{l=1}{\overset{m_a}{}}}e^{\tau _{\nu ,l}(\stackrel{}{x},r_{upp}\stackrel{}{n},t)},`$
where the outer sum is over all $`N`$ grid points of the simulation box each having an opacity $`\kappa _\nu (\stackrel{}{x})`$ and inner sum is over the $`m_a`$ solid angles about each point. The term on the right hand side of Equation (20) due to sources inside the simulation box is $`s_{\nu ,in}`$, which may be computed from emissivity $`j_\nu (\stackrel{}{x})`$ of each hydro grid point of the simulation box:

$`s_{\nu ,in}(t)`$ $`=`$ $`{\displaystyle \frac{1}{4\pi N}}{\displaystyle \underset{\stackrel{}{x}}{}}j_\nu (\stackrel{}{x}){\displaystyle e^{\tau _\nu (\stackrel{}{x},r_{upp}\stackrel{}{n},t)}𝑑\omega (\stackrel{}{n})}`$ (22)
$`=`$ $`{\displaystyle \frac{1}{Nm_a}}{\displaystyle \underset{\stackrel{}{x}}{}}j_\nu (\stackrel{}{x}){\displaystyle \underset{l=1}{\overset{m_a}{}}}e^{\tau _{\nu ,l}(\stackrel{}{x},r_{upp}\stackrel{}{n},t)};`$
this is the radiation due to internal sources that leaves the simulation box. There may be additional sources to the diffuse background. For example, one may wish to take into account rare sources such as quasars that are outside the simulation box. This can be achieved by simply adding the needed terms as $`s_{\nu ,ext}(t)`$.

The formulation of Equation (15) deserves more discussion. The integration over the domain $`\mathrm{\Delta }V`$ is effectively collapsing the angular dimensions to lump together absorbing material into a radial “line” of length $`\mathrm{\Delta }r`$ for that volume element. This could cause an overestimate of opacity for radiation sources that are located inside the same volume element or behind the considered volume element but are spatially displaced perpendicular to the line of sight with respect to the absorbing material. In a more extreme example, a handful of cells of very high optical depth may totally block out radiation downstream, if the mean optical depth is evaluated by linearly averaging optical depth spatially. We propose to modify the original formulation (Equation 15) as follows.

$`<e^{\kappa _\nu (\stackrel{}{x}+s\stackrel{}{n},t)\mathrm{\Delta }L}>`$ $`=`$ $`{\displaystyle \frac{_{\mathrm{\Delta }V}e^{\kappa _\nu (\stackrel{}{x}+s\stackrel{}{n},t)\mathrm{\Delta }L}𝑑V}{\mathrm{\Delta }V}}`$ (23)
$`=`$ $`{\displaystyle \frac{e^{\kappa _\nu (\stackrel{}{y},t)\mathrm{\Delta }L}f(\stackrel{}{x}\stackrel{}{y})d^3\stackrel{}{y}}{\mathrm{\Delta }V}}`$
and

$$\overline{\kappa }_\nu (\stackrel{}{x}+s\stackrel{}{n},t)=\frac{1}{\mathrm{\Delta }L}\mathrm{ln}<e^{\kappa _\nu (\stackrel{}{x}+s\stackrel{}{n},t)\mathrm{\Delta }L}>,$$
(24)
where $`\mathrm{\Delta }L`$ is the size of a hydro cell. This modified formulation is in some sense like averaging the flux instead of the optical depth. This physically motivated modification quite successfully circumvents the forementioned possible problems in realistic situations, as subsequent tests will show (see §3). We adopt this particular scheme to present quantitative results in the next section (§3). We note that in the case of a uniform opacity Equation (15) and Equations (23,24) would yield identical results.

Let us briefly summarize the salient features of the proposed algorithm. What we have shown is that the algorithm is fast, scaling as $`N(\mathrm{log}N)^2`$, where $`N`$ is the number of grid points in a simulation box, independent of the number of radiation sources. Given that computation does not involve any differentiation and all computed integrals are guaranteed to converge in all situations, the algorithm is completely stable and robust. We emphasize that the proposed algorithm, by design, guarantees that flux propagates in the right direction in any situation. The tests presented in the next section will show that the amplitude of flux is also computed with very small errors, indicating that the algorithm is also accurate.

## 3 Tests of the Algorithm

We will present a suite of tests, which are intended to have significant bearings on cosmological applications. We will show that the method performs very well in all tested situations. Clearly, fine tuning of the scheme is worth exploring and is deferred for future work.

For the tests shown we use a $`32^3`$ Cartesian grid with periodic boundary conditions. We use $`m_a=256`$ uniformly spaced angular bin and $`m_r=18`$ logarithmically spaced radial bins for all cases. For the tests shown we do not include cosmological affects; i.e., $`H(t)`$ is set to zero. But the accuracy of the algorithm should not be affected by this simplification.

### 3.1 Static Spherical Optical Depth Distribution

It is guaranteed that the proposed formalism will give the correct flux distribution in the case of no optical attenuation. The simplest non-trivial test then is a source embedded in a static, uniform, non-zero opacity distribution. The flux distribution about a single source is then

$$F(r)=\frac{L}{4\pi r^2}e^{\overline{\kappa }r},$$
(25)
where $`L`$ is the luminosity of the source; $`r`$ is source-centric radius; $`\overline{\kappa }`$ is the opacity. Figure 1 shows the flux distribution in the x-y plane with $`z=16`$, where the luminous source is located at cell $`(16,16,16)`$ and the opacity per cell is $`\overline{\kappa }=0.5`$. The solid contours nearly overlap with the analytic expectations (dashed contours) with some slight deformations from sphericity simply due to the nature of the Cartesian grid. Clearly, the computed flux is accurate.

Since the equation for the flux is linear with respect to the source luminosity, the method would obviously yield results of the same accuracy for any number of distributed sources.

### 3.2 Static Elliptical Optical Depth Distribution

The next step is to test a non-spherical static opacity distribution. An axisymmetric elliptical distribution of opacity of the following form is tested:

$$\kappa (\theta )=\frac{cb}{\sqrt{1(1b^2)\mathrm{cos}^2\theta }},$$
(26)
where $`\theta `$ is the angle with respect to the positive $`z`$ direction, and $`c`$ and $`b`$ are the normalization and the ellipticity parameter of the distribution, respectively; $`\kappa =c`$ at $`\mathrm{cos}\theta =1`$ and $`\kappa =cb`$ at $`\mathrm{cos}\theta =0`$.

Figure 2 shows the computed flux distribution (solid contours) in the x-z plane with $`y=16`$ compared to the analytic result (dashed contours). We use $`c=0.2`$ and $`b=0.2`$ (see Equation 26) for this illustration. We see that the results are quite satisfactory given the Cartesian grid. The agreement with the analytic result is very good except near $`\mathrm{cos}\theta =1`$ where the opacity is the highest and a maximum error of 2 cells is made there. This is caused by a slight underestimate of the optical depth there, due to the averaging scheme of optical depth as indicated by Equations (23,24). We have also tested similar cases by varying either the opacity $`c`$ or ellipticity parameter $`b`$ (see Equation 26) and find that the results have similar accuracies.

### 3.3 Ionization of a Uniform Neutral Medium

We now switch to a higher gear to test the algorithm in a situation where the distribution of opacity changes self-consistently with time due to ionization. For this and subsequent tests we monitor photon number conservation by computing

$$\eta (t)=\frac{N_e(t)+N_{diff}(t)}{N_{ph,tot}(t)},$$
(27)
where $`N_e(t)`$ is the number of free electrons created by time $`t`$ (assuming 100% hydrogen) in the simulation box and $`N_{diff}`$ is the number of photons in the diffuse background in the simulation box volume at time $`t`$ and $`N_{ph,tot}`$ is the total number of photons emitted by time $`t`$ by sources in the simulation box. If photon number conservation is strictly observed, $`\eta `$ would be unity. In all the following tests recombination time is set to infinity to make the problem simple to track. But the accuracy of the method is unchanged by this simplification.

We start with the case of a uniform neutral medium, consisting entirely of hydrogen atoms. An ionizing source of ionizing luminosity $`\dot{N}_{ph}`$ photon/sec is placed at cell $`(16,16,16)`$. The ionizing photons from the source ionize the neutral medium outward and the radius of the ionization front, the surface that separates the interior ionized medium from the exterior neutral medium, evolves with time approximately as:

$`r(t)=13.7({\displaystyle \frac{\dot{N}_{ph}}{10^{51}}})^{1/3}({\displaystyle \frac{n}{10^3\mathrm{cm}^3}})^{1/3}({\displaystyle \frac{t}{10^7\mathrm{yrs}}})^{1/3}\mathrm{kpc},`$ (28)
where $`n`$ is the number density of the neutral hydrogen and $`t`$ is the elapsed time. The reason that this formula (26) is only approximate is that some photons may penetrate further out ahead of the ionization front, if the optical depth per cell is not sufficiently high, and then the “front” is no longer well defined; the formula would be highly accurate for a high opacity distribution. The optical depth per cell is $`\tau _\nu =\sigma _Hn\mathrm{\Delta }x=15.4(n/10^3\mathrm{cm}^3)(\mathrm{\Delta }x/5\mathrm{k}\mathrm{p}\mathrm{c})`$, where the hydrogen photo-ionization cross section $`\sigma _H=10^{18}`$cm<sup>2</sup> is used throughout and $`\mathrm{\Delta }x`$ is the cell size. Thus, for the present test, the optical depth per cell is very large and formula (26) should be very accurate.

Figure 3 shows the contours of the neutral hydrogen fraction in the x-y plane with $`z=16`$ at four epochs. We use $`\dot{N}_{ph}=10^{51}`$photon/sec, $`n=10^3`$cm<sup>-3</sup> and cell size of $`\mathrm{\Delta }x=5`$kpc. The adopted gas density is about 15 times the mean gas density at redshift $`z=6`$. We see that the agreement between the computed results and analytic expectations is excellent at all times. The discrepancy on the radius of the ionization front is evidently no larger than one cell. Given the approximations that have been made in the implementation, especially with respect to the computation of mean opacity (Equations 21,22), it is unclear how well the algorithm conserves photons. Figure 4 shows the degree of conservation of photons as a function of time. It is seen that, except for the first 8 time steps, the number of photons is conserved at better than $`4\%`$ with the average at about $`23\%`$. Timesteps of $`10^7`$yrs are used but the results are insensitive to the timestep; using a timestep 3 times larger yield very similar results except for the upper left corner panel at the beginning of the simulation. Note that the age of the universe at $`z=6`$ is about $`5\times 10^8`$yrs. The indicated convergence of results at this relatively large size of timesteps is clearly advantageous.

### 3.4 Ionization of an $`r^1`$ Neutral Medium

Next, we examine a case of a source sitting in a spherical density distribution, $`n=n_c(r/r_c)^1`$cm<sup>-3</sup>; the central cell at $`r=0`$ is set to a density equal to $`n_c(\mathrm{\Delta }x/r_c)^2`$, where $`\mathrm{\Delta }x`$ is the cell size. This density profile is close to that in the central regions of dark matter halos (Navarro, Frenk, & White 1997) and is of significant cosmological relevance. The expected evolution of the radius of the ionization front is approximately:

$`r(t)=\sqrt{35.2({\displaystyle \frac{\dot{N}_{ph}}{10^{51}\mathrm{sec}^1}})({\displaystyle \frac{n_c}{10^2\mathrm{cm}^3}})^1({\displaystyle \frac{t}{10^7\mathrm{yrs}}})({\displaystyle \frac{r_c}{5\mathrm{k}\mathrm{p}\mathrm{c}}})^1+25\left(({\displaystyle \frac{3}{4\pi }})^{2/3}{\displaystyle \frac{1}{2\pi }}\right)\left({\displaystyle \frac{\mathrm{\Delta }x}{5\mathrm{k}\mathrm{p}\mathrm{c}}}\right)^2}\mathrm{kpc},`$ (29)
where the second (correction) term inside the square root is due to the non-singular core of the modified $`r^1`$ density distribution. The optical depth per cell as a function of source-centric radius is $`\tau (r)=154(n/10^2\mathrm{cm}^3)(\mathrm{\Delta }x/5\mathrm{k}\mathrm{p}\mathrm{c})(r/r_c)^1`$. Therefore, for the adopted values of $`r_c`$, $`n_c`$ and $`\mathrm{\Delta }x`$ (see below), the formula will be very accurate over the range of radii examined and the thickness of the ionization front should be thin and close to one cell.

Figure 5 shows the contours of the neutral hydrogen fraction in the x-y plane with $`z=16`$ at four epochs. We use $`\dot{N}_{ph}=10^{51}`$photon/sec, $`n_c=1.5\times 10^2`$cm<sup>-3</sup>, $`r_c=5`$kpc and cell size of $`\mathrm{\Delta }x=5`$kpc. Note that the adopted $`n_c`$ is approximately $`200`$ times the mean gas density at $`z=6`$. Timestep used is $`10^7`$yrs but the results are insensitive to the timestep. From Figure 5 it is seen that the agreement between the computed results and analytic expectations is excellent at all times and the discrepancy on the radius of the ionization front is no larger than one cell. Figure 6 presents the degree of conservation of photons as a function of time, showing that the number of photons is conserved at better than $`10\%`$ with the average at about $`23\%`$ in this case.

### 3.5 Ionization of an $`r^2`$ Neutral Medium

We turn to a steeper density distribution, $`n=n_c(r/r_c)^2`$cm<sup>-3</sup>; the central cell at $`r=0`$ is set to a density equal to $`n_c(\mathrm{\Delta }x/r_c)^2`$, where $`\mathrm{\Delta }x`$ is the cell size. The density profile is close to that of halos just interior of the virial radius (e.g., Tyson, Kochanski, & dell’Antonio 1998). The expected evolution of the radius of the ionization front is approximately:

$$r(t)=3.4(\frac{\dot{N}_{ph}}{10^{51}})(\frac{n_c}{10^2\mathrm{cm}^3})^1(\frac{t}{10^7\mathrm{yrs}})(\frac{r_c}{5\mathrm{k}\mathrm{p}\mathrm{c}})^1+5\left((\frac{3}{4\pi })^{1/3}\frac{1}{4\pi }\right)(\frac{\mathrm{\Delta }x}{5\mathrm{k}\mathrm{p}\mathrm{c}})\mathrm{kpc},$$
(30)
where the second term on the right hand side is due to the non-singular core imposed on the density profile.

Figure 7 shows the contours of the neutral hydrogen fraction in the x-y plane with $`z=16`$. We use $`\dot{N}_{ph}=10^{51}`$photon/sec, $`n_c=1.5\times 10^2`$cm<sup>-3</sup>, $`r_c=5`$kpc and cell size of $`\mathrm{\Delta }x=5`$kpc at four epochs. Timestep used is $`10^7`$yrs but like in other cases the results are insensitive to the timestep. The computed results and analytic expectations agree very well at all times. For the adopted values of $`n_c`$, $`r_c`$ and $`\mathrm{\Delta }x`$, the optical depth per cell as a function of source-centric radius is $`\tau (r)=154(n/10^2\mathrm{cm}^3)(\mathrm{\Delta }x/5\mathrm{k}\mathrm{p}\mathrm{c})(r/r_c)^2`$ shown in Figure 8; we have $`\tau (r)=3.6`$ for $`r=8`$cells. Thus, the computed shell of the ionization front will be somewhat broadened at late times due to a deeper penetration of ionizing photons outward of the ionization front, clearly evident in Figure 7. Figure 9 shows that the number of photons is conserved at better than $`6\%`$ with the average at about $`14\%`$ for this case.

### 3.6 Ionization of an Elliptical $`r^2`$ Neutral Medium

In realistic situations a source may sit inside a gas cloud that may not be spherical. We examine an elliptical isothermal distribution with the following density distribution

$$n(r,\theta )=n_c(r/r_c)^2\frac{b}{\sqrt{1(1b^2)\mathrm{cos}^2\theta }},$$
(31)
where $`\theta `$ is the angle with respect to the positive $`z`$ direction, and $`b`$ indicates the ellipticity of the distribution. The central cell at $`r=0`$ is designed not to be singular but has a density equal to $`n_c(\mathrm{\Delta }x/r_c)^2`$, where $`\mathrm{\Delta }x`$ is the cell size. The expected, orientation-dependent evolution of the radius of the ionization front follows Equation (30).

Figure 10 shows the contours of the neutral hydrogen fraction in the x-z plane with $`y=16`$ at four epochs. We use $`\dot{N}_{ph}=10^{51}`$photon/sec, $`n_c=1.5\times 10^2`$cm<sup>-3</sup>, $`r_c=5`$kpc and cell size of $`\mathrm{\Delta }x=5`$kpc. Timestep used is $`1.3\times 10^6`$yrs but the results are insensitive to the timestep. It is seen that the algorithm handles the anisotropic density distribution quite well and the agreement between the computed results and analytic expectations is very good at all times, with the computed ionization “fronts” tracing out the analytic results quite nicely with an error no larger than about one cell. The thickening of the computed “fronts” is real due to penetration of photons. Figure 11 shows the degree of conservation of photons. It is seen that the number of photons is conserved at better than $`6\%`$ with the average at about $`14\%`$.

### 3.7 Ionization of a Sharply Divided Anisotropic Neutral Medium

So far we have examined several different cases with smooth density distributions. It is prudent to test the algorithm in a more challenging situation. We pick a case where there are two different density fields separated by a sharp boundary. An axisymmetric density field is prescribed as

$`n(\theta )`$ $`=`$ $`10^3\mathrm{cm}^3\mathrm{for}{\displaystyle \frac{\pi }{4}}\theta {\displaystyle \frac{3\pi }{4}}`$ (32)
$`=`$ $`10^4\mathrm{cm}^3\mathrm{otherwise}.`$
The source sits at cell $`(16,16,16)`$. The expected evolution of the radius of the ionization front follows Equation (28) for each of the two domains with different speeds of the ionization fronts.

Figure 12 shows the contours of the neutral hydrogen fraction in the x-z plane with $`y=16`$ at four epochs. We use $`\dot{N}_{ph}=10^{51}`$photon/sec, $`r_c=5`$kpc and cell size of $`\mathrm{\Delta }x=5`$kpc. Timestep used is $`6.7\times 10^5`$yrs but the results are insensitive to the timestep. The agreement between the computed results and the analytic results are surprisingly good with the radius of the ionization front accurate to about one cell. At the boundaries separating the two regions there is a slight over-shielding on the low density side affected by the high optical depth, high density side. The width of this “affected” region is limited by the size of the solid angle of each discretized angular element, which in the present case is about $`12^2`$ square degrees. Figure 13 shows the degree of conservation of photons, indicating that the number of photons is conserved at better than $`10\%`$ with the average at about $`15\%`$.

### 3.8 Double Sources in a Uniform Density

Up to now we have tested cases with a single radiation source. It is important that the algorithm is capable of handling interactions of multiple radiation sources through the process of percolation, a process that is thought to occur during the cosmological reionization. We will demonstrate this by studying the situation where there are two sources of unequal luminosities. The two sources with luminosities $`\dot{N}_{ph}=10^{51}`$photon/sec, and $`2\times 10^{50}`$photon/sec, respectively, sit in a uniform neutral medium of density $`n=10^3`$cm<sup>-3</sup> initially. A cell size of $`5`$kpc is used for the simulation box. The two sources are located at cell (23,16,16) and cell (13,16,16) with a separation of 10 cells (i.e., $`50`$kpc). Timestep used is $`10^7`$yrs but the results are insensitive to the size of the timestep.

Prior to the overlap of the two HII regions produced by the two sources separately, the evolution of each HII region is separate and the evolution of the radii of the ionization fronts follows Equation (28). Subsequent to the overlap their combined HII region will continue to expand but an analytic solution is not readily available. Figure 14 shows the contours of the neutral hydrogen fraction in the x-y plane with $`z=16`$ at four epochs. It is clear that the algorithm nicely treats the two separate regions without any interference before the overlap, with the two separate ionization fronts traveling at the correct speeds. For the post-overlap era we only show the analytic solution as if overlap has not occurred for the sole purpose of guiding the reader’s eye. The evolution of the ionization front in the post-overlap era is more complicated but the computed results appear to be quite reasonable, with the overlapped region continuing to expand and becoming rounder with time. In this case, one is more keen to check if the total number of photons is conserved. Figure 15 shows the degree of conservation of photons. It is seen the number of photons is conserved at better than $`9\%`$ with the average at about $`12\%`$.

### 3.9 Quadruple Sources in a Uniform Density

Let us make the situation a bit more interesting. We will add two more sources to the case tested in §3.8. The two additional sources have luminosities $`\dot{N}_{ph}=3\times 10^{50}`$photon/sec, and $`\dot{N}_{ph}=5\times 10^{50}`$photon/sec, sitting at cells (20,8,16) and (8,24,16) but turned on with lags of $`1\times 10^8`$yrs and $`2\times 10^8`$yrs (relative to the turn-on time of the initial two sources), respectively.

Figure 16 shows the contours of the neutral hydrogen fraction in the x-y plane with $`z=16`$ at four epochs. For the post-overlap era we only show the analytic solution as if overlap has not occurred for the sake of illustration. We see that the algorithm follows the separate HII regions, as expected. Continuous formation of galaxies in time can evidently be properly followed. Subsequent mergers of HII regions occur as expected. Figure 17 shows the degree of conservation of photons. It is seen the number of photons is conserved at better than $`9\%`$ with the average at about $`13\%`$. The evolution of the ionization front in the post-overlap era is very complicated but the computed results look reasonable, with the overlapped region continuing to expand in a fashion that is expected. The fact that photon number is well conserved and flux is designed to travel in the right direction ensures that the computed results should be accurate.

### 3.10 Shadowing by a Dense Gas Cloud

We will now step back to a case with a single source. But we will include in the uniform density field a dense gas cloud, which is designed to be optically thick to cast a shadow in the direction of radiation propagation. It is expected that shadowing should be common in the real universe since galaxies and other dense clouds (such as damped Lyman alpha systems) are (at least partially) opaque to radiation. A spherical gas cloud of radius $`20`$kpc with a density distribution of $`1\times 10^1(r/20\mathrm{k}\mathrm{p}\mathrm{c})^1`$cm<sup>-3</sup> is centered at cell $`(24,16,16)`$. The rest of the simulation box has a uniform density of $`1\times 10^3`$cm<sup>-3</sup>. The ionizing source is located at cell (16,16,16) with a luminosity of $`10^{51}`$sec<sup>-1</sup>. A cell size of $`\mathrm{\Delta }x=5`$kpc is used for the simulation. Figure 18 shows the neutral hydrogen fraction contours at four epochs. The analytic results are only valid before the ionization front reaches the halo; at subsequent times regions behind the halo within the region delimited by the two dotted lines should remain neutral. We indeed see the shadow cast by the dense gas cloud. The slightly overshadowing near the edges of the shadowed region is due to the limited resolution of angular discretization, also seen in Figure 12.

The relatively complicated situation here does not guarantee that photon number conservation will be obeyed. Figure 19 shows the degree of conservation of photons. Quite reassuringly, we see that the number of photons is conserved at better than $`10\%`$ with the average at about $`12\%`$.

### 3.11 Outside-In Ionization of a Spherical Isothermal Gas Cloud

Let us next examine another class of ionization processes, where external diffuse radiation ionizes isolated overdense regions. This type of ionization process is common in cosmological situations where low density regions become ionized first. In this case we check the photon number conservation in a slightly different way as:

$$\eta (t)=\frac{N_e(t)}{N_{incoming}(t)},$$
(33)
where $`N_e(t)`$ is the number of free electrons created by time $`t`$ in the simulation box and $`N_{incoming}`$ is the number of incoming photons from the diffuse background that are absorbed in the simulation box volume by time $`t`$. If photon number conservation is strictly observed, $`\eta `$ would be unity. Recombination time is set to infinity.

The first simple example is the ionization of a spherical isothermal gas cloud of density $`n=n_c(r/r_c)^2`$cm<sup>-3</sup> with an outer cutoff radius $`r_i`$; the central cell at $`r=0`$ has a density equal to $`n_c(\mathrm{\Delta }x/r_c)^2`$ to avoid singularity, where $`\mathrm{\Delta }x`$ is the cell size. The expected evolution of the radius of the ionization front, which now travels inward towards the center of the gas cloud, is approximately:

$$r(t)=\frac{r_in_cr_c^2}{0.5Fr_it+n_cr_c^2},$$
(34)
where $`F`$ is the diffuse flux assumed to be constant with time. Note that the factor $`0.5`$ in front of flux $`F`$ in Equation (34) is due to the fact that any element on the surface is subject to radiation only at half the total solid angle.

Figure 20 shows the evolution of the ionization front in the x-y plane with $`z=16`$. We use $`r_c=5`$kpc, $`n_c=1.5\times 10^2`$, $`\mathrm{\Delta }x=1`$kpc, $`b=0.5`$ and $`F=1.9\times 10^5`$cm<sup>-2</sup>sec<sup>-1</sup>. We see that the analytic evolution of the ionization front is well tracked by the computed results with error on the radius no larger than one cell. The spiky surface at the earliest time shown reflects the initial condition laid out on a Cartesian grid. Figure 21 shows the photon number conservation as a function of time and again indicates that the method conserves total number of photons very well at $`<1\%`$. The high degree of photon conservation in case is attributable to the fact that there is no obscuration for the cells on the ionization surface that receive ionizing photons.

### 3.12 Outside-In Ionization of an Elliptical Isothermal Gas Cloud

We next consider ionization of an isolated elliptical isothermal gas cloud of the following density profile by a diffuse ionizing background

$$n(r,\theta )=n_c(r/r_c)^2\frac{b}{\sqrt{1(1b^2)\mathrm{cos}^2\theta }};$$
(35)
the central cell at $`r=0`$ has a density equal to $`n_c(\mathrm{\Delta }x/r_c)^2`$ to avoid singularity, where $`\mathrm{\Delta }x`$ is the cell size. The expected evolution of the radius of the ionization front follows Equation (34) with an angle-dependent $`n_c`$ indicated by Equation (35).

Figure 22 shows the evolution of the ionization front in the x-z plane with $`y=16`$. We use $`r_c=5`$kpc, $`n_c=1.5\times 10^2`$, $`\mathrm{\Delta }x=1`$kpc and $`F=1.9\times 10^5`$cm<sup>-2</sup>sec<sup>-1</sup>. We see that the analytic evolution of the ionization front is well traced by the computed results with error on the radius no larger than one cell. Figure 23 shows the photon number conservation as a function of time and again indicates that the method conserves total number of photons very well at $`<1\%`$.

### 3.13 Radiation Propagation In a Realistic Cosmological Density Field

Finally, we turn to the radiation propagation in a realistic density field produced by cosmological simulations. In this case, the situation is much more complicated and no analytic solution is possible. But many of the cases tested above have significant bearings on this realistic case. We take a subbox from one of our latest simulations (Cen et al. 2001) at redshift $`z=6`$. The subbox has a $`32^3`$ grid and a cell size of $`4.65`$kpc proper.

We turn on by hand a galaxy of ionizing luminosity of $`\dot{N}_{ph}=10^{51}`$sec<sup>-1</sup> at cell $`(11,22,21)`$ corresponding to the densest cell in the box, as indicated by the solid dot in Figures (24,25,26). We assume that the gas is composed entirely of hydrogen and at the starting time hydrogen is entirely neutral. For simplicity we have set the recombination time to infinity. It is noted that the galaxy sits in a generic filament, which is embedded in a complex structure containing voids and other interesting structures. Since the subbox is a small region of a large, periodic box, it is not periodic itself. But we treat it as if it were periodic.

Figure 24 and Figure 25 show the distribution of neutral hydrogen at time $`t=(8\times 10^7,4\times 10^8)`$yrs (middle columns), respectively, compared to that at the starting time ($`t=0`$; left column). The most notable feature is that the radiation propagation is highly anisotropic. The ionizing photons quickly clear out widening tunnels towards low density (void) regions in the direction perpendicular to the filament, while radiation front travels at much slower speeds in other directions. Although this is not at all surprising, the figures do remind us of the complex situations encountered in cosmological structures and re-iterate the need for a proper treatment of radiative transfer. We also use a high resolution ray tracing code to obtain numerically the approximate “true” distribution, shown as the right columns. Note that in the ray tracing code one just counts photons and compares to the number of hydrogen atoms along each ray cone; photon penetration ahead of the ionization front is not taken into account and ionization front is taken to be abruptly sharp. Nevertheless, the results obtained with the ray tracing method should be fairly close to the ‘truth”. We see that the agreement between the result obtained with the present method and that with the ray tracing method is excellent.

Figure 26 shows the projected flux distributions (volume-weighted) at the two epochs respectively. We see that flux distributions are quite complex but consistent with the neutral hydrogen distributions. It is seen at the earlier epoch, when substantial blocking or shadowing exists, the contours of the flux distribution tend to orient perpendicular to the contours of the neutral hydrogen density. Note that in the present calculations we have set the recombination time to infinity; in a more involved calculation including recombination it is expected that fluctuations in the flux distributions would be enhanced.

Finally, Figure 27 shows how the total number of photons is conserved. We are once again seeing satisfactory conservation of photons in this realistic simulation at about $`110\%`$ level with a time averaged accuracy of $`35\%`$. This level of photon number conservation is quite adequate for cosmological simulations, since other uncertainties, including galaxy luminosity, radiation escape fraction from the galaxies, star formation efficiency, etc., are usually much larger. Besides, we expect that we may be able to improve the accuracy of the algorithm by more detailed experimenting and careful fine tuning, if necessary.

The fact that photon number is well conserved and flux is designed under the algorithm to propagate in the right directions guarantee that the results should be reliable, as verified by the comparison with the results from a high resolution ray tracing method.

### 3.14 Summary on the Tests

We have performed a variety of tests on the proposed algorithm for transferring radiation in three dimensional space. The tests are designed to be relevant to cosmological applications. It is worthwhile to re-emphasize that the proposed algorithm, by design, guarantees that flux propagates in the right direction in any situation. Thus, the tests should mainly demonstrate how accurately the amplitude of the flux is computed.

To summarize the results we find, the flux is computed very accurately, resulting in ionization fronts that travel at the correct speed with error no larger than one cell in all cases. Photon number is conserved with a maximum error of about $`10\%`$ and an average error of $`15\%`$ over hundreds of time steps for all the cases tested, with $`m_a=256`$ angular elements and $`m_r=18`$ logarithmically spaced radial bins on a $`32^3`$ grid. We note that in the limit of infinite values for $`m_r`$ and $`m_a`$, the present method solves the equation of the radiative transfer exactly and the approximation would represent the truth. In other words, the tests have basically shown that with reasonable number of angular and radial discretization elements the proposed method gives quite accurate results.

Finally, we find that the accuracy of the algorithm does not depend sensitively on the time step used. A time step of size comparable to the typical time step in cosmological hydrodynamic simulations appears sufficient. This feature is desirable with regard to the overall computational cost.

## 4 Discussion and Conclusions

In light of the recent observations of high redshift quasars (e.g., Fan et al. 2001), indicating that the cosmological reionization may be just ending at redshift $`z6`$ (Becker et al. 2001; Barkana 2001; Cen & McDonald 2001), a new and exciting frontier is forming. Inhomogeneous reionization of the universe is now within the observational reach and could potentially provide a great tool to test both cosmological models and astrophysics of galaxy formation (e.g., Barkana & Loeb 2001). It has thus become mandatory to perform cosmological hydrodynamic simulations that tackle this problem with an adequate treatment of radiative transfer. Furthermore, even in the post-reionization era, an initial inhomogeneous reionization process coupled with the significant clustering of radiation sources (galaxies or quasars) as well as the fluctuating density field likely produces fluctuations both in the thermodynamic properties of the cosmic gas and in the ionizing radiation field. Such fluctuations may substantially alter quantitatively (or even qualitatively) the picture of the optically thin regions of the Ly$`\alpha `$ forest on a variety of scales, which would limit our ability of accurately extracting important cosmological information from observations of the Ly$`\alpha `$ forest. Lastly (but not the least), since almost all neutral gas resides and star formation occurs only in dense, optically thick regions, many questions pertaining neutral gas and star formation in the universe can be answered with confidence only when radiative transfer is properly included.

We have developed a fast, accurate and robust algorithm for radiative transfer in three-dimensional space. The core of the algorithm is based on a formulation that the summation of any quantity (such as emissivity or opacity) over any volume can be written in the standard convolution form, which can be computed efficiently using the Fast Fourier Transform techniques. We will name this method “FRT” (Fast Radiative Transfer) method. The overall computational time with this algorithm scales as $`N(\mathrm{log}N)^2`$, where $`N`$ is the number of grid points in a simulation box. Therefore, this is a fast algorithm. We stress that the computational cost with this method is independent of the number of radiation sources, providing a natural match to cosmological simulations where a large number of sources are present. Local sources and diffuse background are naturally split and the evolution of each of the two components as well as interactions between them are handled self-consistently. The devised integral form of the algorithm guarantees that the method is completely stable and robust.

The algorithm is tested on a wide range of analytically tractable problems that have significant bearings on cosmological applications. We find that the algorithm performs very well in all cases, accurately tracking ionization fronts with the error no larger than one cell. Conservation of photons is observed to be at a level of a few percent. The accuracy of the results depends weakly on the size of the time step in all cases; a time step comparable to the typical size of a time step for a cosmological hydrodynamic simulation is sufficient. We also apply the algorithm to a realistic density field produced by a cosmological hydrodynamic simulation and find that conservation of photons is observed at a level of a few percent. These extensive tests indicate that the algorithm is also very accurate.

The tests performed is based on a first implementation of the basic algorithm outlined in the paper. It is clear that there is significant room for further improvement over the specific implementation. The current implementation is on a uniform mesh but we think that variant approaches (such as with adaptive mesh refinement, for example) based on this method are possible and will be explored in the future.

The work is supported in part by grants NAG5-2759, AST91-08103 and ASC93-18185. This implementation of the algorithm will be fine tuned, documented and made available in due course.


text:
# ESTIMATE ON DECELERATION PARAMETER IN A UNIVERSE WITH VARIABLE FINE STRUCTURE “CONSTANT ”.

## Abstract

We determine a cosmological model that includes acceleration of the Universe and time-decreasing fine structure “constant ”, as found recently. We found that the present -day deceleration parameter should be $`q1.1\times 10^5`$. The variation in $`\alpha `$ is attributed to a slowly increasing speed of light.

Pacs 98.80 Hw

ESTIMATE ON DECELERATION PARAMETER IN A UNIVERSE WITH VARIABLE FINE STRUCTURE “CONSTANT ”.

Marcelo S. Berman and Luis A. Trevisan

The latest developments of physics, when we enter this century, are the discovery that the Universe is accelerating, and that the fine structure “constant ” is varying with the age of the Universe.

The first conclusion, based on Supernovae  observations may mean that the value of the cosmological constant is non-null; however, what we know for sure is that that the deceleration parameter,

$$q=\frac{\ddot{R}R}{\dot{R}^2}$$
(1)
is negative. In (1), $`R=R(t)`$ stands for the scale-factor, and overdots are time-derivatives, while the metric is Robertson-Walker’s.

The second finding is that the “constant ”

$$\alpha \frac{e^2}{\frac{h}{2\pi }c},$$
(2)
is varying with time.

Let us call $`F=\mathrm{\Delta }\alpha /(\alpha \mathrm{\Delta }t)`$ (experimental). We can make a rough estimate for $`F`$ by employing the values found by Webb et al for the deviation from the average:

$$\frac{\mathrm{\Delta }\alpha }{\alpha }0.72\times 10^5$$
(3)
spanning $`23\%`$ to $`87\%`$ of the age of the Universe . A rough estimate for $`\mathrm{\Delta }t`$ is given by $`\mathrm{\Delta }t=(0.870.23)t`$ where $`t`$ stands for the age of the Universe. Any competent reader will be able to modify our results when a better estimate shall be available in the literature for $`F.`$

As in cosmology, the charge of the electron, or Planck’s constant, do not play a large-scale rôle, we fix our attention on a time-varying speed of light, $`c=c(t).`$

Webb et al. pointed out that a common property of unified theories, applied to Cosmology, is that they predict space and time-dependence of the coupling constants. Let us suppose that Newton’s gravitational “constant” $`G`$ is also time varying; the natural framework we need is JBD ( Jordan-Brans-Dicke) theory with varying speed, as Barrow has been studying .

It is well accepted that, for the present Universe a power-law of time can represent the scale-factor. We are thus led to consider constant deceleration parameter models, as studied initially by Berman , and Berman and Gomide. On defining $`m=1+q=`$constant, we find,

$$R=R(t)=(mDt)^{\frac{1}{m}}$$
(4)
where $`D=`$constant.

Einstein-de Sitter model, with $`k=0`$ ($`k`$ is the tricurvature), is of this kind, with $`m=3/2.`$ As our time scale for the present Universe duration so far is about 1010 years, and there is not evidence for a significant time variation of $`q,`$ for the present Universe, we find that a constant deceleration parameter model is reasonable for the present Universe. In reference  there is a comment on one single evidence for a changing $`q`$ for the Universe in its present phase, but obviously this is not a conclusive evidence up to now. We shall show that the solution with variable $`\alpha `$ differs from Einstein-de Sitter’s model.

Barrow’s equation read:

$$H^2=\frac{8\pi \rho }{\varphi }\frac{\dot{\varphi }}{\varphi }H+\frac{\omega }{6}\frac{\dot{\varphi }^2}{\varphi ^2}\frac{kc^2(t)}{R^2}$$
(5)
$$\frac{3\ddot{R}}{R}=\frac{8\pi }{(3+2\omega )\varphi }\left[(2+\omega )\rho +\frac{3(1+\omega )p}{c^2(t)}\right]\frac{\ddot{\varphi }}{\varphi }\frac{\omega \dot{\varphi }^2}{\varphi ^2}$$
(6)
$$\ddot{\varphi }+3H\dot{\varphi }=\frac{8\pi }{3+2\omega }\left(\rho \frac{3p}{c^2(t)}\right)$$
(7)
Unlike Barrow, we shall not impose a definite equation of state; the values for energy density $`\rho `$ and cosmic-pressure $`p`$ will arise naturally in our solution.

From the theory of constant deceleration parameter we have:

$$H=(mt)^1$$
(8)
We find the following solution:

$$R=(mDt)^{\frac{1}{m}}$$
(9)
$$\varphi =At^n$$
(10)
$$\rho =ABt^{n2}$$
(11)
$$c=Ct^{\frac{1}{m}1}$$
(12)
$$p=Ft^{\gamma 4+2/m}$$
(13)
where $`A,B,\gamma ,F,C`$ are constants, related by means of equations (4), (5) and (6).

From (11) and (2), we find,

$$\frac{\dot{\alpha }}{\alpha }=Hq$$
(14)
Let us now compare with the experimental values:

$$\frac{\mathrm{\Delta }\alpha }{\alpha \mathrm{\Delta }t}=\frac{\mathrm{\Delta }\alpha }{\alpha \left[0.64\frac{H^1}{m}\right]}\frac{0.72}{0.64}10^5Hm$$
(15)
or

$$\frac{\mathrm{\Delta }\alpha }{\alpha \mathrm{\Delta }t}1.1\times 10^5Hm.$$
(16)
For the age of the Universe, we make use of relation (7). We see that $`q1.1\times 10^{5\text{ }}`$ and $`m0.99999.`$

The speed of light is, then increasing slowly with time. The Universe has a minimal, though finite, acceleration, as can be checked from (8), with the values of $`q`$ and $`m.`$

It has been a happy coincidence, that when the $`\alpha `$ varying experiments were published, we were able to apply the recent Barrow-Jordan-Brans-Dicke equations, in order to derive a value for $`q,`$ which was lacking in the literature. Barrow and Magueijo  explained how a changing $`\alpha `$ could yield the Supernova results with a different argumentation than ours. For up to date references in relation to variable $`G`$ cosmologies, see ref  . We refer to these papers for further information.

Acknowledgments

Both authors thank support by Prof Ramiro Wahrhftifig, Secretary of Science and Technology of the State of Parana, and from our Institutions, especially to Jorge L.Valgas, Roberto Merhy, Mauro K. Nagashima, Carlos Fior, C.R. Kloss, J.L.Buso, and Roberto Almeida. M.S.B. thanks also his two intelectual mentors, Prof. M.M.Som and Prof. F.M. Gomide.


text:
# Cardy’s Formula for some Dependent Percolation Models

## 1 Introduction

It was understood by physicists since the early seventies that critical statistical mechanics models should possess continuum scaling limits with a global conformal invariance that goes beyond pure scale invariance. The phenomenon is particularly interesting in two dimensions, where every analytic function gives rise to a conformal transformation and the local conformal transformations form an infinite dimensional group; in that context, it was first studied by Belavin, Polyakov and Zamolodchikov . For an introduction to the methods of conformal field theory as applied to two-dimensional critical percolation, see .

Until recently, however, there was no rigorous mathematical proof of this phenomenon, with the exception of the Simple Symmetric Random Walk, whose continuum scaling limit is Brownian Motion. Then, S. Smirnov managed to prove  existence, uniqueness and conformal invariance of the continuum scaling limit of critical site percolation on the triangular lattice, obtaining in particular conformal invariance of crossing probabilities and Cardy’s formula for rectangular crossings .

In this paper we show that there are some natural *dependent* percolation models for which conformal invariance of the crossing probabilities and Cardy’s formula can be proved. Our proof relies on Smirnov’s result and on properties of the dependent percolation models which make them, in a sense to be specified later, “small perturbations” of the independent model treated by Smirnov.

The dependent percolation models we consider are the distributions at time $`n1`$ (including the final state as $`n\mathrm{}`$) of a discrete time deterministic dynamical process $`\sigma ^n`$ with state space $`\{1,+1\}^𝕃`$ consisting of assignments of $`1`$ or $`+1`$ to a regular lattice $`𝕃`$. The initial $`\sigma ^0`$ is “uniformly random”, i.e., the distribution of $`\sigma ^0`$ is a Bernoulli(1/2) product measure. The dynamics are those of Domany’s stochastic Ising ferromagnet  at zero temperature. There are two essentially equivalent versions — one where $`𝕃`$ is the hexagonal lattice $``$ and one where it is the triangular lattice $`𝕋`$. We take $``$ and $`𝕋`$ to be regular lattices embedded in $`^2`$ so that the elementary cells of $``$ (resp., $`𝕋`$) are regular hexagons (resp., equilateral triangles). In the first version, $``$, as a bipartite graph, is partitioned into two subsets $`𝒜`$ and $``$ which are alternately updated so that each $`\sigma _x`$ is forced to agree with a majority of its three neighbors (which are in the other subset). In the second version, all sites are updated simultaneously according to a rule based on a deterministic pairing of the six neighbors of every site into three pairs (see the end of Section 2 for a complete explanation). The rule is that $`\sigma _x`$ flips if and only if it disagrees (after the previous update) with both sites in two or more of its three neighbor pairs; thus there is (resp., is not) a flip if the number $`D_x`$ of disagreeing neighbors is $`5`$ (resp., $`3`$) and there is also a flip for some cases of $`D_x=4`$. We note that Cardy’s formula can also be verified for a modified rule in which there is a flip if and only if $`D_x5`$; the case of a modified rule where there is a flip if and only if $`D_x4`$ is an interesting open problem.

## 2 Definition of the model(s) and results

In this section we give a more detailed description of the dependent percolation models and results.

Consider the homogeneous ferromagnet on the hexagonal lattice $``$ with states denoted by $`\sigma =\{\sigma _x\}_x,\sigma _x=\pm 1`$, and with (formal) Hamiltonian

$$=\underset{x,y}{}\sigma _x\sigma _y,$$
(1)
where $`_{x,y}`$ denotes the sum over all pairs of neighbor sites, each pair counted once. The variables $`\sigma _x,\sigma _y`$ are called spins. We write $`𝒩^{}(x)`$ for the set of three neighbors of $`x`$, and indicate with

$$\mathrm{\Delta }_x(\sigma )=2\underset{y𝒩^{}(x)}{}\sigma _x\sigma _y$$
(2)
the change in the Hamiltonian when the spin $`\sigma _x`$ at site $`x`$ is flipped (i.e., changes sign).

Notice that the hexagonal lattice can be partitioned into two subsets $`𝒜`$ and $``$ in such a way that all three neighbors of any site in $`𝒜`$ (resp., $``$) are in $``$ (resp., $`𝒜`$). By placing an edge between any two sites of $`𝒜`$ (resp., $``$) that are next-nearest neighbors in $``$, the subset $`𝒜`$ (resp., $``$) becomes a triangular lattice. (This relation between an hexagonal lattice and its triangular “sublattice,” sometimes expressed in terms of a “star-triangle transformation,” will be used again in Remark 2.1 below.) We now consider the discrete time Markov process $`\sigma ^n,n`$, with state space $`𝒮=\{1,+1\}^{}`$, which is the zero temperature limit of a model of Domany , constructed as follows:

* The initial state $`\sigma ^0`$ is chosen from a symmetric Bernoulli product measure.
* At odd times $`n=1,3,\mathrm{}`$, the spins in the sublattice $`𝒜`$ are updated according to the following rule: $`\sigma _x,x𝒜`$, is flipped if and only if $`\mathrm{\Delta }_x(\sigma )<0`$.
* At even times $`n=2,4,\mathrm{}`$, the spins in the sublattice $``$ are updated according to the same rule as for those of the sublattice $`𝒜`$.

In order to present the main result of this paper, let us denote by $`\sigma ^{\mathrm{}}`$ the final state of the process $`\sigma ^n`$ defined above. $`\sigma ^{\mathrm{}}=lim_n\mathrm{}\sigma ^n`$ exists with probability one, as was proved in , and, like $`\sigma ^n`$ for $`1n<\mathrm{}`$, defines a dependent percolation model on $``$. These are the the main objects of our investigation.

We will call $`\delta `$ the “mesh” of the lattice and consider the continuum scaling limit of the dependent percolation model $`\sigma ^n`$ on $`\delta `$ as $`\delta 0`$. For simplicity of exposition, we will prove Cardy’s formula in the special case of a rectangle, aligned with the coordinate axes and of given cross-ratio $`\eta `$ (a similar approach would work for any domain with a “regular” boundary, but it would involve dealing with more complex deformations of the boundary). Consider a finite rectangle $`=(a,b)(a/2,a/2)\times (b/2,b/2)^2`$ with sides of lengths $`a`$ and $`b`$, such that the cross-ratio $`a/b`$ is $`\eta `$. We say that there is (in $`\sigma ^n`$) a vertical plus-crossing if $`\delta `$ contains a path of $`+1`$ spins from $`\sigma ^n`$ joining the top and bottom sides of the rectangle $``$, and call $`P_\delta (\eta ;n)`$ the probability of such a plus-crossing at time $`n`$. More precisely, there is a vertical plus crossing if there is a path $`x_0,x_1,\mathrm{},x_m,x_{m+1}`$ in $``$ with $`\sigma _{x_j}^n=+1`$ for all $`j`$, with $`\delta x_1,\mathrm{},\delta x_m`$ all in $``$, and with the line segments $`\overline{\delta x_0,\delta x_1}`$ and $`\overline{\delta x_m,\delta x_{m+1}}`$ touching respectively the top side $`[a/2,a/2]\times \{b/2\}`$ and the bottom side $`[a/2,a/2]\times \{b/2\}`$. In the next section we will prove the following result:

###### Theorem 1.

For all $`n1`$ (including $`n=\mathrm{}`$), the limit $`P(\eta ;n)=lim_{\delta 0}P_\delta (\eta ;n)`$ exists and is given by Cardy’s formula:

$$P(\eta ;n)=F_C(\eta )\frac{\mathrm{\Gamma }(\frac{2}{3})}{\mathrm{\Gamma }(\frac{4}{3})\mathrm{\Gamma }(\frac{1}{3})}\eta ^{\frac{1}{3}}{}_{2}{}^{}F_{1}^{}(\frac{1}{3},\frac{2}{3};\frac{4}{3};\eta ).$$
(3)
A stronger result than Theorem 1 can be obtained, i.e., it is possible to prove existence, uniqueness and conformal invariance of the continuum scaling limit, as proven by Smirnov  for independent site percolation on the triangular lattice. Such a result, though, requires more work and will be pursued in a future paper. Here we just note that the proof is based on showing that the limit for our dependent percolation models (on the hexagonal lattice) coincides with that of Smirnov for independent percolation on the triangular lattice, i.e., that the models belong to the same universality class.

The following observations are useful in understanding the behavior of the model and will help in the proof of Theorem 1.

* The values of the spins in the sublattice $`𝒜`$ at time zero are irrelevant, since at time $`1`$, after the first update, those values are uniquely determined by the values of the spins in the sublattice $``$.
* Once the initial spin configuration in the sublattice $``$ is chosen, the dynamics is completely deterministic.
* A site can no longer flip once it belongs to either a loop or “barbell” of constant sign in $``$, where a loop means a simple loop (with no subloops) and a barbell consists of two disjoint loops connected by a path (we regard a loop as a degenerate barbell).

We also note that, by studying the percolation properties of the final state $`\sigma ^{\mathrm{}}`$ on the infinite lattice $``$, it can be shown that every site is in some barbell of constant $`\sigma ^{\mathrm{}}`$-sign .

The discrete time Markov process defined above can be considered a simplified version of a continuous time process where an independent (rate $`1`$) Poisson clock is assigned to each site $`x`$, and the spin at site $`x`$ is updated (with the same rule as in our discrete time process) when the corresponding clock rings. The percolation properties of the final state $`\sigma ^{\mathrm{}}`$ of that process were studied, both rigorously and numerically, in ; the results there (about critical exponents rather than critical crossing probabilities) strongly suggest that that dependent percolation model is also in the same universality class as independent percolation. Similar stochastic processes on different types of lattices have been studied in various papers. See, for example,  for models on $`^d`$ and  for a model on the homogeneous tree of degree three. Such models are also discussed extensively in the physics literature, usually on $`^d`$ (see, for example,  and ). On the hexagonal lattice, the discrete time dynamics is the zero-temperature case of Domany’s dynamics . Numerical simulations have been done by Nienhius  and rigorous results for both the continuous and discrete dynamics have been obtained in , including a detailed analysis of the discrete time (synchronous) case. The analysis of  is at the heart of this paper, and we will refer to and heavily rely on it for the proof of Theorem 1, which is given in the next section.

There is an alternative, but equivalent, way of describing the discrete time dynamics as a deterministic cellular automaton on the triangular lattice $`𝕋`$ (with random initial state). The initial state is again chosen by assigning value $`+1`$ or $`1`$ independently, with equal probability, to each site of the triangular lattice. Given some site $`\overline{x}𝕋`$, group its six $`𝕋`$-neighbors $`y`$ in three disjoint pairs $`\{y_1^{\overline{x}},y_2^{\overline{x}}\},\{y_3^{\overline{x}},y_4^{\overline{x}}\},\{y_5^{\overline{x}},y_6^{\overline{x}}\}`$, so that $`y_1^{\overline{x}}`$ and $`y_2^{\overline{x}}`$ are $`𝕋`$-neighbors, and so on for the other two pairs. Translate this construction to all sites $`x𝕋`$, thus producing three pairs of sites $`\{y_1^x,y_2^x\},\{y_3^x,y_4^x\},\{y_5^x,y_6^x\}`$ associated to each site $`x𝕋`$. (Note that this construction does not need to specify how $`𝕋`$ is embedded in $`^2`$.) Site $`x`$ is updated at times $`m=1,2,\mathrm{}`$ according to the following rule: the spin at site $`x`$ is changed from $`\sigma _x`$ to $`\sigma _x`$ if and only if at least two of its pairs of neighbors have the same sign and this sign is $`\sigma _x`$.

###### Remark 2.1.

This dynamics on the triangular lattice $`𝕋`$ is equivalent to the alternating sublattice dynamics on the hexagonal lattice $``$ when restricted to the sublattice $``$ for even times $`n=2m`$. To see this, start with $`𝕋`$ and construct an hexagonal lattice $`^{}`$ by means of a star-triangle transformation (see, for example, p. 335 of ) such that a site is added at the center of each of the triangles $`(x,y_1^x,y_2^x),(x,y_3^x,y_4^x)`$, and $`(x,y_5^x,y_6^x)`$. $`^{}`$ may be partitioned into two triangular sublattices $`𝒜^{}`$ and $`^{}`$ with $`^{}=𝕋`$. It is now easy to see that the dynamics on $`𝕋`$ for $`m=1,2,\mathrm{}`$ and the alternating sublattice dynamics on $`^{}`$ restricted to $`^{}`$ for even times $`n=2m`$ are the same.

Theorem 1 (and its generalizations) in this context means that, at all times $`m0`$, the crossing probabilities for the states $`\sigma ^m`$ of this cellular automaton on $`𝕋`$ have the same conformally invariant continuum scaling limit as that for critical independent percolation on $`𝕋`$, despite the dependence induced by the cellular automaton dynamics.

## 3 Proof of Theorem 1

In this final section of the paper we prove Theorem 1. We follow the notation of  and start by giving some definitions. Let us consider a loop $`\gamma `$ in the triangular sublattice $``$, written as an ordered sequence of sites $`(y_0,y_1,\mathrm{},y_n)`$ with $`n3`$, which are distinct except that $`y_n=y_0`$. For $`i=1,\mathrm{},n`$, let $`\zeta _i`$ be the unique site in $`𝒜`$ that is an $``$-neighbor of both $`y_{i1}`$ and $`y_i`$. We call $`\gamma `$ an *s-loop* if $`\zeta _1,\mathrm{},\zeta _n`$ are all distinct. Similarly, a (site-self avoiding) path $`(y_0,y_1,\mathrm{},y_n)`$ in $``$, between $`y_0`$ and $`y_n`$, is called an *s-path* if $`\zeta _1,\mathrm{},\zeta _n`$ are all distinct. Notice that any path in $``$ between $`y`$ and $`y^{}`$ (seen as a collection of sites) contains an s-path between $`y`$ and $`y^{}`$. An s-loop of constant sign is stable for the dynamics since at the next update of $`𝒜`$ the presence of the constant sign s-loop in $``$ will produce a stable loop of that sign in the hexagonal lattice. Similarly an s-path of constant sign between $`y`$ and $`y^{}`$ will be stable if $`y`$ and $`y^{}`$ are stable — e.g., if they each belong to an s-loop. A triangular loop $`x_1,x_2,x_3`$ with a common $``$-neighbor $`\zeta 𝒜`$ is called a *star*; it is not an s-loop. A triangular loop in $``$ that is not a star is an s-loop and will be called an *antistar*, while any loop in $``$ that contains more than three sites contains an s-loop.

Before stating a lemma, that will be a main ingredient in the proof of Theorem 1, we need one more definition. For $`(x,x^{})`$ an ordered pair of neighbors in $``$, we define the “partial cluster” $`C_{(x,x^{})}^{}`$ to be the set of sites $`y`$ such that there is a (site-self avoiding) path $`x_0=x^{},x_1,\mathrm{},x_n=y`$ in $``$ of constant sign in $`\sigma ^0`$, with $`x_1x`$ and $`(x_0=x^{},x_1,x)`$ not forming a star. Combining the stability properties of s-loops and s-paths just discussed, we have the following lemma.

###### Lemma 3.1.

An s-path $`(y_0,\mathrm{},y_m)`$ in $``$ of constant sign in $`\sigma ^0`$ is stable (i.e., retains that same sign in $`\sigma ^n`$ for all $`0n\mathrm{}`$) if $`C_{(y_1,y_0)}^{}`$ and $`C_{(y_{m1},y_m)}^{}`$ both contain s-loops.

Proof of Lemma 3.1. The original s-path $`(y_0,\mathrm{},y_m)`$ is stable because either $`y_0`$ and $`y_m`$ both belong to s-loops of constant sign in $`\sigma ^0`$ or else there is a longer s-path of constant sign in $`\sigma ^0`$, between some $`y`$ and $`y^{}`$ (with the original $`(y_0,\mathrm{},y_m)`$ as a subpath), such that both $`y`$ and $`y^{}`$ belong to s-loops of constant sign in $`\sigma ^0`$.

With this preparation, we are now ready to start the proof of Theorem 1. What we will prove, roughly speaking, is that, in the limit $`\delta 0`$, there exists a vertical plus-crossing of $``$ from $`\sigma ^n`$ with $`n1`$, in $`\delta `$, if and only if there exists a vertical plus-crossing of $``$ from $`\sigma ^0`$ in $`\delta `$. Since $``$ is a triangular lattice and the initial state $`\sigma ^0`$ is chosen from a symmetric Bernoulli product measure, this implies that the limit $`P(\eta ;n)=lim_{\delta 0}P_\delta (\eta ;n)`$ exists for $`n1`$ and is the same as in the case of the crossing probability for independent site percolation on the triangular lattice, thus proving the theorem.

Consider two rectangles, $`^{}=(a^{},b^{})`$ with $`b^{}`$ slightly larger than $`b`$ and $`a^{}`$ slightly smaller than $`a`$, and $`^{\prime \prime }=(a^{\prime \prime },b^{\prime \prime })`$ with $`b^{\prime \prime }`$ slightly smaller than $`b`$ and $`a^{\prime \prime }`$ slightly larger than $`a`$. Call $`P_\delta ^{}(a^{},b^{})`$ the probability of a vertical plus-crossing from $`\sigma ^0`$ in $`^{}\delta `$ joining the top and bottom sides of $`^{}`$ and $`P_\delta ^{\prime \prime }(a^{\prime \prime },b^{\prime \prime })`$ the probability of a *horizontal minus*-crossing from $`\sigma ^0`$ in $`^{\prime \prime }\delta `$ joining the left and right sides of $`^{\prime \prime }`$. Note that a vertical plus crossing (on the triangular lattice $`\delta `$) occurs if and only if a horizontal minus-crossing does not occur. Clearly, from  we have

$$P^{}(a^{},b^{})\underset{\delta 0}{lim}P_\delta ^{}(a^{},b^{})=F_C(a^{}/b^{}),$$
(4)
$$\underset{a^{}a,b^{}b}{lim}P^{}(a^{},b^{})=P^{}(a,b)=F_C(\eta ),$$
(5)
and

$$\underset{a^{\prime \prime }a,b^{\prime \prime }b}{lim}\underset{\delta 0}{lim}P_\delta ^{\prime \prime }(a^{\prime \prime },b^{\prime \prime })=1F_C(\eta ).$$
(6)
Any vertical plus-crossing of $`^{}\delta `$ at time $`0`$ yields a vertical plus-crossing by some s-path $`(y_0,\mathrm{},y_m)`$, which then yields at time $`1`$ a vertical plus-crossing of $`\delta `$ by a path $`(y_{k_1},\zeta _{k_1+1},\mathrm{},\zeta _{k_2},y_{k_2})`$, providing $`a^{}<a`$, $`b^{}b`$ and $`\delta `$ is sufficiently small. (The reason we first take $`b^{}>b`$ and then let $`b^{}b`$ is to handle the case of time $`n>1`$, as we shall see.) Therefore, for small $`\delta `$,

$$P_\delta (\eta ;n=1)P_\delta ^{}(a^{},b^{}).$$
(7)
On the other hand, if there is a horizontal minus-crossing of $`^{\prime \prime }\delta `$ at time $`0`$, it produces a horizontal minus-crossing in $`\delta `$ at time $`1`$ (for small $`\delta `$) which blocks any possible vertical plus-crossing in $`\delta `$ at that time; therefore, for small $`\delta `$,

$$P_\delta (\eta ;n=1)1P_\delta ^{\prime \prime }(a^{\prime \prime },b^{\prime \prime }).$$
(8)
Letting $`\delta 0`$ and then $`a^{},a^{\prime \prime }a`$ and $`b^{},b^{\prime \prime }b`$ and using (5)-(8), we conclude that $`P_\delta (\eta ;n=1)`$ converges to Cardy’s formula, $`F_C(\eta )`$, as $`\delta 0`$.

It remains to prove that the same is true for all times $`n2`$. In order to do that, we first have to show that our vertical plus-crossing of $`^{}\delta `$ by $`(y_0,\zeta _1,\mathrm{},\zeta _m,y_m)`$ created at time $`1`$ doesn’t “shrink” too much due to the effect of the dynamics, so that at all later times, including $`n=\mathrm{}`$, there is a vertical plus-crossing of $`\delta `$ by $`(y_{k_1},\zeta _{k_1+1},\mathrm{},\zeta _{k_2},y_{k_2})`$.

To do this by extending the bound (7) to all $`n1`$, at the cost of a correction to the right hand side that tends to zero with $`\delta `$, we apply Lemma 3.1. Noting that each of the partial paths $`(y_0,\mathrm{},y_{k_1})`$ and $`(y_{k_2},\mathrm{},y_m)`$ contains of the order of $`(b^{}b)/\delta `$ sites, we see that the lemma implies that it suffices to show that there is some $`\beta >0`$ and $`K<\mathrm{}`$ such that for any deterministic $`(x,x^{})`$,

$$P(|C_{(x,x^{})}^{}|\mathrm{}\text{ and }C_{(x,x^{})}^{}\text{ contains no antistar})Ke^\beta \mathrm{}.$$
(9)
To prove (9), we partition $``$ into disjoint antistars and denote by $`\tau `$ the collection of these antistars. We do an algorithmic construction of $`C_{(x,x^{})}^{}`$ (as in, e.g., ), where the order of checking the sign of sites is such that when the first site in an antistar from $`\tau `$ is checked (and found to have the same sign as $`x^{}`$), then the other two sites in that antistar are checked next. Without loss of generality, we assume that $`\sigma _x^{}^0=+1`$. Then standard arguments show that the probability in (9) is bounded by $`K(1(\frac{1}{2})^3)^{(\mathrm{}/3)}`$.

To similarly extend the bound (8), one proceeds in the same way, but considering horizontal minus-crossings of $`^{\prime \prime }\delta `$ at time zero which produce horizontal minus-crossings of $`\delta `$ at time $`n1`$. Taking the limits $`\delta 0,a^{}a,b^{}b`$ concludes the proof.