text:
# Exact stationary correlations for a far-from-equilibrium spin chain

## I Introduction

Exact solutions of simple model systems play a key role in statistical mechanics. They provide reliable quantitative information about regular as well as singular behavior and serve as proving grounds for approximation schemes. Moreover, they can set important pointers in areas where a general theoretical framework is still lacking. A prime example is the study of non-equilibrium stationary states (NESS). Similar to systems in thermal equilibrium, NESS are characterized by time-independent macroscopic observables; however, at present these can only be computed with explicit reference to the imposed dynamics: there is, as yet, no equivalent of Gibbs ensemble theory for far-from-equilibrium steady states. As a consequence, most progress to date is made by studying specific models.

Since the Ising model  is one of the most intimately known interacting many-particle systems, many non-equilibrium models depart from it by imposing an external force which drives the system out of equilibrium. Examples include couplings to multiple heat baths , competing Glauber and Kawasaki dynamics , or a current-inducing global bias . In all of these models, reviewed in , the non-equilibrium perturbation violates the detailed balance symmetry of the equilibrium dynamics. On large time and length scales, one can probe how this affects universal behavior, or, adopting a more microscopic but no less fundamental perspective, one can ask how non-universal properties, such as the exact configurational probabilities, are modified. One finds, generically, that non-equilibrium forces have especially profound effects if (i) they couple to the bulk rather than the boundaries  of the system , and (ii) the dynamics satisfies a conservation law . For example, Ising lattice gases with conserved particle number can be driven into two distinct non-equilibrium universality classes, depending on the symmetries of the external forces . In contrast, Ising-like systems with non-conserved dynamics remain in the universality class of the equilibrium Ising model, even if the usual $`_2`$-symmetry of the Ising model is broken . Even then, however, the violation of detailed balance leads to fundamental changes in the configurational probabilities which must be identified and interpreted before we may hope to formalize our understanding of NESS. It is here that exact solutions are expected to be most helpful.

Motivated by these considerations, we recently  investigated a very simple non-equilibrium Ising-like model, namely, a one-dimensional interacting spin chain with spin-flip dynamics coupled to two temperature baths. The rates are a simple but nontrivial generalization of the familiar Glauber  rates: spins at odd (even) sites are coupled to a temperature $`T_o`$ ($`T_e`$). In two dimensions, this model  exhibits an order-disorder phase transition which belongs to the Ising class, according to renormalization group arguments  and Monte Carlo simulations . In one dimension, the two-spin correlation function can be calculated exactly . Its only singularity lies at $`T_o=T_e=0`$, so that the lower critical dimension is $`d=1`$. Seeking an expression for the steady state, we solved the master equation perturbatively, in an expansion in the temperature difference of the two heat baths, up to and including second order terms . To our surprise, the full stationary distribution turned out to be rather complex: At each order, additional spin operators appear, characterized by longer spatial range interactions and higher-order spin products, and lower-order coupling constants acquire corrections. So, at first order, one encounters next-nearest neighbor pair interactions, while at second order, two new terms appear: a next-next-nearest neighbor pair interaction, and a four-spin interaction spanning four nearest-neighbor sites. All of these are allowed by symmetry, and none allowed by symmetry are absent. Given this structure, one can, at least in principle, extrapolate to higher orders in perturbation theory.

Of course, these findings immediately raise an obvious question: how does this relatively complicated stationary distribution generate long-wavelength behavior in the Ising universality class? In this article, we pursue an alternate route towards the answer: Instead of aiming for the stationary state directly, we seek its representation in terms of correlation functions \- since the knowledge of all stationary correlation functions is equivalent to knowing the steady state. Starting from the master equation, we derive a hierarchy of equations of motion for correlation functions of arbitrary numbers of spins. In $`d=1`$, this hierarchy is closed and soluble. Remarkably, we find that - in contrast to the apparent complexity of the stationary distribution itself - the correlation functions are very simple. Apart from two nontrivial amplitudes which reflect the temperature difference between the two sublattices, the structure of correlations is completely analogous to the equilibrium Ising model: specifically, for $`N\mathrm{}`$, arbitrary (even) $`m`$-point correlation functions factorize into a product of $`m/2`$ two-point correlations. Of course, all correlation functions involving an odd number of spins vanish by symmetry.

While some exact analytic results for steady-state distributions are available, they are confined to three classes of systems: first, one-dimensional lattice gas models, restricted to excluded volume interactions, such as the asymmetric exclusion process and its relatives , second, very special one-dimensional spin systems whose master equations are solved by the Ising Boltzmann factor , and third, interacting systems in one or two dimensions on very small lattices so that the number of degrees of freedom remains manageable . To the best of our knowledge, the work presented here is amongst the first complete solutions for non-equilibrium stationary states with nontrivial (nearest-neighbor) interactions and arbitrary number of degrees of freedom.

The article is organized as follows. We first introduce our model and its master equation. Next, we derive the equations of motion for arbitrary correlation functions. Following a brief review of the solution  for the two-point correlation function, we show how four-point correlations can be factorized into two-point correlations in the long chain ( $`N\mathrm{}`$ ) limit. We then postulate that all correlation functions factorize in this manner, and show that this factorization solves the equations of motion. In an appendix, we generalize this solution to finite systems with periodic boundary conditions. We conclude with some comments and open questions.

## II The model

Our model is defined on a one-dimensional ring, with an even number $`N`$ of sites, and periodic boundary conditions. A spin variable, $`\sigma _i=\pm 1`$, denotes the value of the spin at site $`i`$, and nearest-neighbor spins interact according to the usual Ising Hamiltonian,

$$=J\underset{i}{}\sigma _i\sigma _{i+1}$$
(1)
with an exchange coupling $`J`$. The dynamics is a non-equilibrium generalization of the usual Glauber model : spins on even and odd lattice sites experience different temperatures, $`T_e`$ and $`T_o`$. Specifically, a configuration $`\{\sigma \}=\{\sigma _1,\sigma _2,\mathrm{},\sigma _N\}`$ evolves into a new configuration by flipping a randomly selected spin $`\sigma _i`$ with a rate

$$w_i(\sigma _i\sigma _i)=1\frac{\gamma _i}{2}\sigma _i(\sigma _{i1}+\sigma _{i+1}),$$
(2)
where

$$\gamma _i=\{\begin{array}{c}\gamma _e=\mathrm{tanh}(2J/k_BT_e),ieven\hfill \\ \gamma _o=\mathrm{tanh}(2J/k_BT_o),iodd.\hfill \end{array}$$
(3)
Thus, the full time-dependent configurational probability $`p(\{\sigma \};t)`$ evolves according to a master equation,

$$_tp(\{\sigma \};t)=\underset{i=1}{\overset{N}{}}\left[w_i(\sigma _i\sigma _i)p(\{\sigma \};t)+w_i(\sigma _i\sigma _i)p(\{\sigma ^{[i]}\};t)\right]$$
(4)
where $`\{\sigma ^{[i]}\}`$ differs from $`\{\sigma \}`$ by a flip of the $`i`$-th spin. A trivial time scale has been set to unity, and we use dimensionless units for inverse temperature, i.e., $`\beta _eJ/(k_BT_e)`$, etc.

Our goal in the following is to find a representation for the stationary solution of Eqn (4), $`q(\{\sigma \})lim_t\mathrm{}p(\{\sigma \};t)`$. This limit is unique, since Eqn (4) is ergodic: every configuration $`\{\sigma \}`$ can be reached in finite time from every other configuration $`\{\sigma ^{}\}`$ (unless $`T_e=T_o=0`$). For equal temperatures $`TT_e=T_o`$, the steady state is just the (canonical) distribution for the Ising chain,

$$q_o\left(\{\sigma \}\right)=\frac{1}{Z}\mathrm{exp}\left(H/k_BT\right).$$
(5)
It is of course straightforward to compute arbitrary correlation functions for the equilibrium Ising chain. Since the Ising model is invariant under a global spin flip ($`_2`$-symmetry), only correlations of even numbers of spins are nonzero, and are easily expressed in terms of the parameter $`\overline{\omega }=\mathrm{tanh}(J/k_BT)`$. Of course, we may – and always will – order the arguments of an $`m`$-point correlation function ($`m`$ even) without loss of generality, such that $`1k_1<k_2<\mathrm{}<k_mN`$. Then, in the $`N\mathrm{}`$ limit one finds easily:

$`\sigma _{k_1}\mathrm{}\sigma _{k_m}^{eq}`$ $`=`$ $`\overline{\omega }^{\left(k_2k_1\right)+\left(k_4k_3\right)+\mathrm{}+\left(k_mk_{m1}\right)}`$ (6)
$`=`$ $`\sigma _{k_1}\sigma _{k_2}^{eq}\sigma _{k_3}\sigma _{k_4}^{eq}\mathrm{}\sigma _{k_{m1}}\sigma _{k_m}^{eq}\text{ for }m\text{ even and }N\mathrm{},`$ (7)
i.e., a general correlation function factorizes into a product of two-point functions.

We now turn to the non-equilibrium model, characterized by different temperatures, $`T_eT_o`$. The associated stationary state violates detailed balance  and differs from the Boltzmann distribution, Eqn (5). The degree to which detailed balance is violated can be measured by the parameter $`d(\gamma _o\gamma _e)/2`$. Similar to the Ising model, the stationary state is invariant under a global spin flip ($`_2`$). As a consequence, all stationary correlations of an odd number of spins vanish identically, and again, only even correlations need to be discussed. We also note the symmetry under translations by $`n`$ lattice sites, $`𝕋_n`$, combined with $`dd`$ if $`n`$ is odd:

$`q(\{\sigma \};d)`$ $`=`$ $`q(\{\sigma \};d)`$ (8)
$`q(\{\sigma \};d)`$ $`=`$ $`q(𝕋_n\{\sigma \};(1)^nd)`$ (9)
While a direct exact solution of the master equation has proven elusive, a perturbative calculation, in powers of $`d`$, shows  that the stationary distribution for this non-equilibrium model is rather complicated, with longer-range and higher-order spin operators appearing. Specifically, writing $`q\left(\{\sigma \}\right)\stackrel{~}{Z}^1\mathrm{exp}\left[V\left(\{\sigma \}\right)\right]`$, we find that the potential function $`V\left(\{\sigma \}\right)`$, to second order in $`d`$, has the form

$$V\left(\{\sigma \}\right)=\overline{\beta }\underset{i}{}\sigma _i\sigma _{i+1}+d\lambda \underset{i}{}(1)^i\sigma _i\sigma _{i+2}+2\left(d\lambda \right)^2\underset{i}{}\left[\sigma _i\sigma _{i+2}+\sigma _i\sigma _{i+1}\sigma _{i+2}\sigma _{i+3}\mathrm{coth}(2\overline{\beta })\left(\sigma _i\sigma _{i+1}+\sigma _i\sigma _{i+3}\right)\right],$$
(10)
where $`\mathrm{tanh}(2\overline{\beta })\left(\gamma _e+\gamma _o\right)/2`$ and $`\lambda =\frac{1}{8}\mathrm{sinh}(4\overline{\beta })`$. An analysis of the structure of the perturbation series indicates that, at each order, additional spin operators appear, consisting of spatially longer-range interactions and higher-order spin products. It is therefore quite remarkable that the correlation functions turn out to be very similar to those of the Ising model, as we now proceed to show.

## III Equations of motion for arbitrary correlation functions

We begin by deriving the equations of motion for arbitrary time-dependent $`m`$-point spin correlation functions, $`\sigma _{k_1}\sigma _{k_2}\mathrm{}\sigma _{k_m}_t`$, starting from the master equation:

$`_t\sigma _{k_1}\sigma _{k_2}\mathrm{}\sigma _{k_m}_t={\displaystyle \underset{i=1}{\overset{N}{}}}\left\{{\displaystyle \underset{\{\sigma \}}{}}\sigma _{k_1}\sigma _{k_2}\mathrm{}\sigma _{k_m}\left[w_i(\sigma _i\sigma _i)p(\{\sigma \};t)+w_i(\sigma _i\sigma _i)p(\{\sigma ^{[i]}\};t)\right]\right\}`$

Here, the subscript $`_t`$ will be used to distinguish time-dependent averages from their stationary limits, $`lim_t\mathrm{}_t`$. Due to the sum over all configurations, the {} bracket on the right hand side obviously vanishes for all sites $`i`$ which do not belong to the set $`\{k_1,k_2,\mathrm{},k_m\}`$, and one finds easily that

$`_t\sigma _{k_1}\sigma _{k_2}\mathrm{}\sigma _{k_m}_t=2{\displaystyle \underset{i=1}{\overset{m}{}}}\sigma _{k_1}\sigma _{k_2}\mathrm{}\sigma _{k_m}w_{k_i}(\sigma _{k_i}\sigma _{k_i})_t.`$

Inserting Eqn (2) for $`w_{k_i}(\sigma _{k_i}\sigma _{k_i})`$ and taking the infinite time limit, we obtain the equations satisfied by the stationary correlation functions, $`lim_t\mathrm{}\sigma _{k_1}\sigma _{k_2}\mathrm{}\sigma _{k_m}_t\sigma _{k_1}\sigma _{k_2}\mathrm{}\sigma _{k_m}`$:

$$0=2m\sigma _{k_1}\sigma _{k_2}\mathrm{}\sigma _{k_m}+\underset{i=1}{\overset{m}{}}\gamma _{k_i}\sigma _{k_1}..\sigma _{k_{i1}}(\sigma _{k_i+1}+\sigma _{k_i1})\sigma _{k_{i+1}}..\sigma _{k_m}$$
(11)
As mentioned before, all stationary correlations of odd numbers of spins are trivially zero, and we need to focus only on the case of even $`m`$. Moreover, all stationary correlations are translationally invariant, modulo the sublattice structure. As in the Glauber model , the hierarchy of correlation functions is closed: the equation for $`m`$-point correlations does not involve any higher correlations, and it is homogeneous provided there are no nearest-neighbor pairs among the arguments. Otherwise, the $`m`$-point correlations couple to lower correlation functions, which appear as inhomogeneities. In the following section, we motivate an ansatz for $`\sigma _{k_1}\sigma _{k_2}\mathrm{}\sigma _{k_m}`$, and then show that it satisfies Eqn (11) for any choice of $`\{k_1,k_2,\mathrm{},k_m\}`$.

## IV Exact solutions for stationary correlations.

To establish several key relations, we briefly review the exact solution for the two-point functions $`\sigma _i\sigma _j`$ . For this case, Eqn (11) reads

$$0=4\sigma _i\sigma _j+\gamma _i(\sigma _{i+1}+\sigma _{i1})\sigma _j+\gamma _j\sigma _i(\sigma _{j+1}+\sigma _{j1}).$$
(12)
For nearest-neighbor sites, e.g., $`j=i+1`$, this equation becomes inhomogeneous:

$$0=4\sigma _i\sigma _{i+1}+\gamma _i\sigma _{i1}\sigma _{i+1}+\gamma _{i+1}\sigma _i\sigma _{i+2}+\gamma _i+\gamma _{i+1}$$
(13)
For simplicity, we restrict the discussion in this section to correlation functions in the thermodynamic limit $`N\mathrm{}`$. Thus we may label the lattice sites by the integers, $`i`$. Correlations on finite periodic chains and some details for the $`N\mathrm{}`$ limit will be addressed in Appendix A. Following , it is easy to show that the solution is unique and takes the form (for $`i<j`$, without loss of generality)

$$\sigma _i\sigma _j=\sqrt{A_iA_j}\omega ^{ji}$$
(14)
with the spatial decay length controlled by the parameter

$$\omega =\frac{1}{\sqrt{\gamma _e\gamma _o}}\left(1\sqrt{1\gamma _e\gamma _o}\right)$$
(15)
The key difference to the equilibrium Ising model is the emergence of two amplitudes, to match the sublattice symmetries, namely,

$$A_i\{\begin{array}{c}A_e=(\gamma _e+\gamma _o)/(2\gamma _o),\text{if}i\text{even}\hfill \\ A_o=(\gamma _e+\gamma _o)/(2\gamma _e),\text{if}i\text{odd}\hfill \end{array}$$
(16)
For later reference, we note that these relations imply a further identity for products of pair correlations, namely, for any $`i<j<j+1<l`$,

$$0=\left(\gamma _j+\gamma _{j+1}\right)\sigma _i\sigma _l\gamma _j\sigma _i\sigma _{j+1}\sigma _{j+1}\sigma _l\gamma _{j+1}\sigma _i\sigma _j\sigma _j\sigma _l.$$
(17)
This relation will be needed in the discussion of four-point functions.

It is quite remarkable that the three quantities $`\omega `$, $`A_e`$, and $`A_o`$ also determine all higher correlation functions, through a structure which is almost perfectly analogous to the equilibrium case. Deferring the case of finite $`N`$ to Appendix A, we consider only an infinite chain here. Focusing on even correlations, we assert that arbitrary $`m`$-point correlation functions, with $`k_1<k_2<\mathrm{}<k_m`$ and $`m`$ even, are given by

$$\sigma _{k_1}\mathrm{}\sigma _{k_m}=\sqrt{A_{k_1}A_{k_2}\mathrm{}A_{k_m}}\omega ^{\left(k_2k_1\right)+\left(k_4k_3\right)+\mathrm{}+\left(k_mk_{m1}\right)}=\sigma _{k_1}\sigma _{k_2}\sigma _{k_3}\sigma _{k_4}\mathrm{}\sigma _{k_{m1}}\sigma _{k_m}.$$
(18)
i.e., higher correlation functions factorize into two-spin correlations as in the equilibrium case, Eqn (7). However, two features distinguish these correlations from their equilibrium counterparts. First, the spatial dependence is controlled by a different parameter, $`\omega `$. In analogy to $`\overline{\omega }`$ for the equilibrium system, $`\omega `$ defines an effective temperature for the non-equilibrium system, via $`\omega \mathrm{tanh}(J/k_BT_{\mathrm{eff}})`$. $`T_{\mathrm{eff}}`$ diverges with $`T_o`$ or $`T_e`$ but vanishes only if both $`T_o`$ and $`T_e`$ go to zero. As a result, the correlation length $`\xi `$, defined via $`\mathrm{ln}\omega \xi ^1`$, diverges only if both temperatures vanish. Second, and more importantly, we note the appearance of the even/odd amplitudes $`A_e,A_o`$, reflecting the sublattice identities of the two spins. These amplitudes carry the primary information about the non-equilibrium nature of our dynamics.

Before turning to a general proof of Eqn (18), it is instructive to confirm it explicitly for the four-point functions. In this case, factorization implies

$$\sigma _i\sigma _j\sigma _k\sigma _l=\sigma _i\sigma _j\sigma _k\sigma _l,$$
(19)
for $`i<j<k<l`$. When inserting this ansatz into the right hand side of Eqn (11), we need to distinguish whether $`j`$ and $`k`$ are nearest-neighbor sites or not. If $`k>j+1`$, we obtain

$`8\sigma _i\sigma _j\sigma _k\sigma _l+\gamma _i(\sigma _{i+1}+\sigma _{i1})\sigma _j\sigma _k\sigma _l+\gamma _j\sigma _i(\sigma _{j+1}+\sigma _{j1})\sigma _k\sigma _l`$ (21)
$`+\gamma _k\sigma _i\sigma _j(\sigma _{k+1}+\sigma _{k1})\sigma _l+\gamma _l\sigma _i\sigma _j\sigma _k(\sigma _{l+1}+\sigma _{l1})`$
$`=`$ $`\left[4\sigma _i\sigma _j+\gamma _i(\sigma _{i+1}+\sigma _{i1})\sigma _j+\gamma _j\sigma _i(\sigma _{j+1}+\sigma _{j1})\right]\sigma _k\sigma _l`$ (22)
$`+\left[4\sigma _k\sigma _l+\gamma _k(\sigma _{k+1}+\sigma _{k1})\sigma _l+\gamma _l\sigma _j\sigma _k(\sigma _{l+1}+\sigma _{l1})\right]\sigma _i\sigma _j`$
$`=`$ $`0,`$ (23)
since the expressions in the square brackets vanish for the stationary two-point functions, by virtue of Eqns (12) or (13), depending on whether ($`i,j`$) or ($`k,l`$) are nearest neighbors or not. If $`k=j+1`$, we need to consider a slightly different equation, namely

$`8\sigma _i\sigma _j\sigma _{j+1}\sigma _l+\gamma _i(\sigma _{i+1}+\sigma _{i1})\sigma _j\sigma _{j+1}\sigma _l+\gamma _j\sigma _i(\sigma _{j+1}+\sigma _{j1})\sigma _{j+1}\sigma _l`$ (25)
$`+\gamma _{j+1}\sigma _i\sigma _j(\sigma _{j+2}+\sigma _j)\sigma _l+\gamma _l\sigma _i\sigma _j\sigma _{j+1}(\sigma _{l+1}+\sigma _{l1})`$
$`=`$ $`\left[4\sigma _i\sigma _j+\gamma _i(\sigma _{i+1}+\sigma _{i1})\sigma _j+\gamma _j\sigma _i\sigma _{j1}\right]\sigma _{j+1}\sigma _l+\gamma _j\sigma _i\sigma _l`$ ()
$`+\sigma _i\sigma _j\left[4\sigma _{j+1}\sigma _l+\gamma _{j+1}\sigma _{j+2}\sigma _l+\gamma _l\sigma _{j+1}(\sigma _{l+1}+\sigma _{l1})\right]+\gamma _{j+1}\sigma _i\sigma _l`$
$`=`$ $`\gamma _j\sigma _i\sigma _l+\gamma _{j+1}\sigma _i\sigma _l\gamma _j\sigma _i\sigma _{j+1}\sigma _{j+1}\sigma _l\gamma _{j+1}\sigma _i\sigma _j\sigma _j\sigma _l`$ (27)
$`=`$ $`0.`$ (28)
The last equality follows from Eqn (17).

The proof is completed by induction. Assuming that the factorization has been proven for $`(m2)`$-point functions, it is sufficient to show that the ansatz

$$\sigma _{k_1}\mathrm{}\sigma _{k_m}=\sigma _{k_1}\mathrm{}\sigma _{k_{m2}}\sigma _{k_{m1}}\sigma _{k_m}$$
(29)
solves the equations of motion for the $`m`$-point functions. Again, we need to distinguish whether $`k_{m2}`$ and $`k_{m1}`$ are nearest neighbors or not. If they are separated by more than one lattice spacing, i.e., if $`k_{m1}>k_{m2}+1`$, we have

$`2m\sigma _{k_1}\sigma _{k_2}\mathrm{}\sigma _{k_m}+{\displaystyle \underset{i=1}{\overset{m}{}}}\gamma _{k_i}\sigma _{k_1}..\sigma _{k_{i1}}(\sigma _{k_i+1}+\sigma _{k_i1})\sigma _{k_{i+1}}..\sigma _{k_m}`$ (30)
$`=`$ $`\sigma _{k_{m1}}\sigma _{k_m}[2(m2)\sigma _{k_1}\mathrm{}\sigma _{k_{m2}}+{\displaystyle \underset{i=1}{\overset{m2}{}}}\gamma _{k_i}\sigma _{k_1}..\sigma _{k_{i1}}(\sigma _{k_i+1}+\sigma _{k_i1})\sigma _{k_{i+1}}..\sigma _{k_{m2}}]`$ ()
$`+\sigma _{k_1}\mathrm{}\sigma _{k_{m2}}\left[2\sigma _{k_{m1}}\sigma _{k_m}+\gamma _{k_{m1}}(\sigma _{k_{m1}+1}+\sigma _{k_{m1}1})\sigma _{k_m}+\gamma _{k_m}\sigma _{k_{m1}}(\sigma _{k_m+1}+\sigma _{k_m1})\right]`$
$`=`$ $`0`$ (32)
Here, we have used the fact that both square brackets vanish, since they enclose the equations of motion for the two- and $`(m2)`$-point functions, respectively. Next, we consider the case $`k_{m1}=k_{m2}+1`$. Again, we add and subtract the terms missing from the full equations of motion of the lower correlations functions:

$`2m\sigma _{k_1}\sigma _{k_2}\mathrm{}\sigma _{k_m}+{\displaystyle \underset{i=1}{\overset{m}{}}}\gamma _{k_i}\sigma _{k_1}..\sigma _{k_{i1}}(\sigma _{k_i+1}+\sigma _{k_i1})\sigma _{k_{i+1}}..\sigma _{k_m}`$ (33)
$`=`$ $`\gamma _{k_{m2}}\sigma _{k_1}..\sigma _{k_{m3}}\sigma _{k_{m2}+1}\sigma _{k_{m1}}\sigma _{k_m}+\gamma _{k_{m2}}\sigma _{k_1}..\sigma _{k_{m3}}\sigma _{k_m}`$ (35)
$`\gamma _{k_{m1}}\sigma _{k_1}..\sigma _{k_{m2}}\sigma _{k_{m1}1}\sigma _{k_m}+\gamma _{k_{m1}}\sigma _{k_1}..\sigma _{k_{m3}}\sigma _{k_m}`$
$`=`$ $`\sigma _{k_1}..\sigma _{k_{m4}}\{(\gamma _{k_{m2}}+\gamma _{k_{m1}})\sigma _{k_{m3}}\sigma _{k_m}`$ (36)
$`\gamma _{k_{m2}}\sigma _{k_{m3}}\sigma _{k_{m2}+1}\sigma _{k_{m1}}\sigma _{k_m}\gamma _{k_{m1}}\sigma _{k_{m3}}\sigma _{k_{m2}}\sigma _{k_{m1}1}\sigma _{k_m}\}`$
$`=`$ $`0.`$ (37)
To obtain the last two identities, we have factored out an $`(m4)`$-point correlation, and used Eqn (17). This concludes the proof. We emphasize again that the condition $`N1`$ has been imposed for reasons of simplicity alone: expressions for arbitrary correlation functions on finite rings, while somewhat more cumbersome, are easily derived and also follow the Ising pattern, as we will show in Appendix A.

Finally, we discuss one special case of the two temperature Glauber dynamics where the two–point correlations do not have the form of Eqn (14). This is the case, if one of the two parameters $`\gamma _e,\gamma _o`$ is zero, the other finite, i. e. the temperature of one sublattice is infinite. Without loss of generality we assume in the following $`\gamma _o=0`$. One can then easily show directly from the equation of motion (12) that the two-point function reads

$$\sigma _i\sigma _j\{\begin{array}{c}\gamma _e/4,\text{if}ji=1\hfill \\ \gamma _e^2/8,\text{if}ji=2\text{and}i\text{even}\hfill \\ 0,\text{otherwise}\hfill \end{array}$$
(38)
Hence, spins are only correlated over a distance of at most two lattice sites. Turning to the higher correlations, we observe that Eqn (17) also holds for these parameter values. Therefore, $`m`$-point correlations ($`m`$ even) again factorize into two-point correlations.

## V Conclusions

To summarize, we have found an exact solution for all stationary correlation functions of a one-dimensional non-equilibrium Ising spin chain with an even number $`N`$ of sites. The system is globally coupled to two temperature baths: spins on odd (even) lattice sites experience a temperature $`T_o`$ ($`T_e`$) and flip according to a generalization of the familiar Glauber rates. The presence of two different temperatures violates detailed balance and maintains a nontrivial non-equilibrium stationary state. The complete set of correlation functions provides us with the full, exact solution for this steady state. This allows us to reconcile two potentially contradictory earlier findings: while simulations and renormalization group arguments indicate that the long-wavelength behavior of non-conserved two-temperature models should be Ising-like, a perturbative analysis of the stationary distribution $`q\left(\{\sigma \}\right)`$ showed that a zoo of long-range and multi-spin interactions are present in $`\mathrm{ln}q\left(\{\sigma \}\right)`$. Here, we have shown that the exact correlation functions for this model are remarkably Ising-like: they decay exponentially with a characteristic correlation length which diverges only if both both $`T_o`$ and $`T_e`$ vanish. Also, in the $`N\mathrm{}`$ limit, arbitrary (even) $`m`$-point correlations factorize into products of two-point correlations, following exactly the same scheme in both the Ising and the two-temperature chain. The only key difference is the appearance of two nontrivial amplitudes, $`A_o`$ and $`A_e`$, which reflect the sublattice symmetry: for each spin on an odd (even) site, the correlation function carries a factor of $`\sqrt{A_o}`$ ($`\sqrt{A_e}`$). Of course, stationary correlations of odd numbers of spins vanish identically, due to symmetry. In short, the correlations of the model are entirely consistent with its Ising-like long-distance properties.

It is remarkable - and not at all immediately obvious - that these two amplitudes should be the only remnants of the zoo of interactions in $`q\left(\{\sigma \}\right)`$. There is, however, a very simple and elegant representation  of the stationary state in terms of an extended Ising model, consisting of $`2N`$ spins on a comb-like $`d=1`$ lattice with Hamiltonian $`_{aux}=J_is_is_{i+1}J_o_{\text{odd }i}s_i\sigma _iJ_e_{\text{even }i}s_i\sigma _i`$. Here, the spins $`\{s_i\}`$ form an auxiliary set which must be traced out in order to obtain stationary observables associated with the original variables $`\{\sigma _i\}`$. If the interactions $`J`$, $`J_o`$ and $`J_e`$ are tuned appropriately in the complex plane, the $`\sigma `$-correlations of $`_{aux}`$ are identical to those of the two-temperature model. Further, the exact stationary state of the two-temperature model follows as $`q\left(\{\sigma \}\right)\text{Tr}_{\{s_i\}}\mathrm{exp}[_{aux}]`$. Details and generalizations will be presented elsewhere .

It would of course be interesting to investigate other non-equilibrium versions of Glauber dynamics. Will the correlations still have a rather simple structure if a spin chain is coupled, in a translation-invariant manner, to $`n>2`$ different temperatures? Will they still factorize, provided a sufficiently large number of amplitudes is introduced? Clearly, we still need to cover lots of ground before even such simple non-equilibrium systems are fully understood.

Acknowledgements We thank H.J. Hilhorst, R.K.P. Zia, J. Slawny, U.C. Täuber, W. Just, E. Ben-Naim and M. Hastings for fruitful discussions. Financial support from the NSF through the Division of Materials Research and from the DFG through contract no. SCHM 1537 is gratefully acknowledged.

## A Stationary correlations in finite chains

In this Appendix, we address the correlation functions for our model on a finite periodic chain (a ring) of $`N`$ sites, with $`i=1,2,\mathrm{},N`$. To set the scene, we first review the equilibrium case, with uniform temperature $`T`$. The correlation functions of the Ising model on a ring are well known. Using an ordered set of arguments, $`1k_1<k_2<\mathrm{}<k_mN`$ without loss of generality, one finds:

$`\sigma _{k_1}\mathrm{}\sigma _{k_m}^{eq}`$ $`=`$ $`0\text{ for odd }m\text{ },`$ (A1)
$`\sigma _{k_1}\mathrm{}\sigma _{k_m}^{eq}`$ $`=`$ $`{\displaystyle \frac{1}{1+\overline{\omega }^N}}\left\{\overline{\omega }^{\left(k_2k_1\right)+\left(k_4k_3\right)+\mathrm{}+\left(k_mk_{m1}\right)}+\overline{\omega }^{N\left(k_2k_1\right)\left(k_4k_3\right)\mathrm{}\left(k_mk_{m1}\right)}\right\}\text{ for even }m.`$ (A2)
with $`\overline{\omega }=\mathrm{tanh}(J/k_BT)<1`$. The periodicity of the ring implies that $`\sigma _{k_1}\mathrm{}\sigma _{k_m}^{eq}`$ is invariant under $`\overline{\omega }^n\overline{\omega }^{Nn}`$ for any integer $`0nN`$. For $`N\mathrm{}`$, Eqn (A2) obviously reduces to Eqn (7).

Turning to the non-equilibrium case with two temperatures, we first assert that the stationary two-point correlation function, with $`1i<jN`$, is given by

$$\sigma _i\sigma _j=\frac{1}{Z(\omega )}\sqrt{A_iA_j}\left[\omega ^{ji}+\omega ^{N(ji)}\right].$$
(A3)
where the parameter $`\omega `$ and the amplitudes $`A_i\{A_e,A_o\}`$ are simply those of Eqns (15) and (16), and the normalization factor $`Z(\omega )`$ is given by $`Z(\omega )=1+\omega ^N`$. Clearly, $`\sigma _i\sigma _j`$ is periodic with period $`N`$. For easy reference, we fill in some details here which are also relevant to the $`N\mathrm{}`$ case.

Eqn (12) is essentially a second order difference equation in the variable $`ji`$. To find a unique solution, we need two boundary conditions: one of these is Eqn (13) for $`ji=1`$ and the other is the requirement that $`\sigma _i\sigma _j`$ be periodic with period $`N`$ (or, for $`N\mathrm{}`$, that $`\sigma _i\sigma _j`$ vanish for large separations). We therefore expect a superposition of two linearly independent solutions whose coefficients can then be determined. Away from the boundaries, the ansatz $`\sigma _i\sigma _j=\sqrt{A_iA_j}\omega ^{ji}`$ reduces Eqn (12) to three relations, for the combinations $`i,j`$ odd/odd, even/even, and odd/even respectively:

$`0`$ $`=`$ $`4A_o+2\gamma _o\sqrt{A_eA_o}\left(\omega +\omega ^1\right)`$
$`0`$ $`=`$ $`4A_e+2\gamma _e\sqrt{A_eA_o}\left(\omega +\omega ^1\right)`$
$`0`$ $`=`$ $`4\sqrt{A_eA_o}+\left(\gamma _oA_e+\gamma _eA_o\right)\left(\omega +\omega ^1\right)`$
We note that the third equation is simply a linear combination of the first and second, expressing the fact that there are only two amplitudes, and not three, as one might have assumed (cf. ) based on the three types (odd/odd, even/even, odd/even) of pair correlations. Moreover, the symmetry under $`\omega \omega ^1`$ implies that, for each solution $`\omega `$, $`\omega ^1`$ is also a solution. Proceeding to solve the system, we find $`\gamma _eA_o=\gamma _oA_e`$ and

$`\omega +\omega ^1={\displaystyle \frac{2}{\sqrt{\gamma _o\gamma _e}}}\omega _\pm ={\displaystyle \frac{1}{\sqrt{\gamma _e\gamma _o}}}\left(1\pm \sqrt{1\gamma _e\gamma _o}\right)`$

with $`\omega _+=1/\omega _{}`$. The two roots $`\omega _\pm `$ of this quadratic equation provide us with the two anticipated linearly independent solutions. However, to satisfy the inhomogeneous Eqn (13) for both finite and infinite $`N`$, we have to take a convex combination of the two solutions

$$\sigma _i\sigma _j=\alpha \sqrt{A_iA_j}\omega _{}^{ji}+\beta \sqrt{A_iA_j}\omega _+^{ji}\text{with}\alpha +\beta =1,$$
(A4)
which then leads to $`\gamma _eA_o+\gamma _oA_e=\gamma _e+\gamma _o`$. This relation, together with the previous identity $`\gamma _eA_o=\gamma _oA_e`$, determines the values of the two amplitudes $`A_e`$ and $`A_o`$. Demanding periodicity for the ansatz (A4) results in a second identity for the two integration constants $`\alpha ,\beta `$, namely $`\beta =\alpha \omega _{}^N`$. Comparing the result with our assertion, Eqn (A3), we identify $`\omega \omega _{}`$ and $`Z(\omega )1+\omega _{}^N`$. Since $`\omega _{}<1`$, the $`N\mathrm{}`$ limit is also obvious.

To treat the general case, we set up some notation. For any even $`m2`$, we define two auxiliary functions which depend on an ordered set of arguments $`1k_1<k_2<\mathrm{}<k_mN`$,

$`S_m^{(f)}(k_1,k_2,\mathrm{}k_m)`$ $``$ $`\sqrt{A_{k_1}A_{k_2}\mathrm{}A_{k_m}}\omega ^{\left(k_2k_1\right)+\left(k_4k_3\right)+\mathrm{}+\left(k_mk_{m1}\right)}=S_2^{(f)}(k_1,k_2)S_2^{(f)}(k_3,k_4)\mathrm{}S_2^{(f)}(k_{m1},k_m)`$
$`S_m^{(b)}(k_1,k_2,\mathrm{}k_m)`$ $``$ $`\sqrt{A_{k_1}A_{k_2}\mathrm{}A_{k_m}}\omega ^{\left(k_2k_1\right)\left(k_4k_3\right)\mathrm{}\left(k_mk_{m1}\right)}=S_2^{(b)}(k_1,k_2)S_2^{(b)}(k_3,k_4)\mathrm{}S_2^{(b)}(k_{m1},k_m)`$
Thus, we rewrite Eqn (A3) in the form

$`\sigma _i\sigma _j={\displaystyle \frac{1}{Z(\omega )}}\left[S_2^{(f)}(i,j)+\omega ^NS_2^{(b)}(i,j)\right].`$

and recall that both $`S_2^{(f)}(i,j)`$ and $`S_2^{(b)}(i,j)`$ satisfy the equations of motion for the $`\mathrm{i}\mathrm{n}\mathrm{f}\mathrm{i}\mathrm{n}\mathrm{i}\mathrm{t}\mathrm{e}`$ chain, i.e., (12), (13), and (17).

We now proceed as follows: First, we argue that both $`S_m^{(f)}(k_1,k_2,\mathrm{}k_m)`$ and $`S_m^{(b)}(k_1,k_2,\mathrm{}k_m)`$ solve the homogeneous and inhomogeneous equations of motion (11) for $`m`$-point functions on the infinite chain. For $`S_m^{(f)}(k_1,k_2,\mathrm{}k_m)=lim_N\mathrm{}\sigma _{k_1}\mathrm{}\sigma _{k_m}`$, this is immediately obvious. For $`S_m^{(b)}(k_1,k_2,\mathrm{}k_m)`$, we simply need to retrace the inductive proof from Section 4: Since the basic relations (12), (13), and (17) all hold for $`S_2^{(b)}`$, we may replace each $`\sigma _{k_1}\mathrm{}\sigma _{k_m}`$ in Eqns (17-18) by the appropriate $`S_m^{(b)}`$ without violating any of the equalities. Since some of the equations of motion (11) for the $`m`$-point correlations are inhomogeneous, only a convex combination $`\alpha S_m^{(f)}(k_1,k_2,\mathrm{}k_m)+\beta S_m^{(b)}(k_1,k_2,\mathrm{}k_m)`$ with $`\alpha +\beta =1`$ is a solution.

Finally, we have to to satisfy the last remaining property of the full correlation function, namely, periodicity. Given an ordered set of arguments $`1k_1<k_2<\mathrm{}<k_mN`$, we assert that a general correlation function of $`m2`$ spins takes the form

$`\sigma _{k_1}\mathrm{}\sigma _{k_m}`$ $`=`$ $`0\text{ for odd }m\text{ },`$ (A5)
$`\sigma _{k_1}\mathrm{}\sigma _{k_m}`$ $`=`$ $`{\displaystyle \frac{\sqrt{A_{k_1}A_{k_2}\mathrm{}A_{k_m}}}{Z(\omega )}}\left\{\omega ^{\left(k_2k_1\right)+\left(k_4k_3\right)+\mathrm{}+\left(k_mk_{m1}\right)}+\omega ^{N\left[\left(k_2k_1\right)+\left(k_4k_3\right)+\mathrm{}+\left(k_mk_{m1}\right)\right]}\right\}`$ (A6)
$`=`$ $`{\displaystyle \frac{1}{Z(\omega )}}\left\{S_m^{(f)}(k_1,k_2,\mathrm{}k_m)+\omega ^NS_m^{(b)}(k_1,k_2,\mathrm{}k_m)\right\}\text{for even }m`$ (A7)
which is obviously periodic on the ring. Comparing Eqns (A2) and (A6), it is manifest how closely the correlations of the non-equilibrium model mirror those of the Ising chain. Again, the only significant difference is the appearance of the two amplitudes, $`A_e`$ and $`A_o`$.


format: text

text:
# The Hydrodynamics and Chemical Evolution of Starburst-driven Outflows

## 1. Introduction

Massive starbursts, well localized short episodes of violent star formation, are among the intrinsic characteristics of the evolution of galaxies. They are found both at high and intermediate redshifts as well as in galaxies in the Local universe.

The large kinetic luminosity produced by violent bursts of star formation is well known to drastically affect the surrounding interstellar medium (ISM), generating giant superbubbles, supershells and in extreme cases supergalactic winds. Powerful gas outflows from star forming regions are important redistributors of the ISM mass and momentum. On the other hand, the gas ejected by supernova explosions contains the products of massive stars. Therefore gaseous outflows initiated by starburst affect also the chemical evolution of the galactic ISM and the intergalactic medium.

Supergalactic winds are expected to be detectable thanks to their extended X-ray emission. Indeed, diffuse X-ray emission associated with starburst galaxies has been detected by ROSAT, ASCA, BeppoSAX, Chandra and XMM missions in several starburst sources. However, the origin of the soft X-ray emission, its connection with Hα filaments and the chemical composition of the X-ray emitting gas, remain ambiguous.

The paper deals with three central aspects of the resultant outflows. The chemical evolution expected in the interior of superbubbles. The conditions required to establish a supergalactic wind and the intrinsic properties and observational manifestations of fully developed supergalactic winds.

## 2. Hot gas chemical composition

The metallicity of the interstellar and intergalactic media reflects the history of star formation in the Universe. Galactic supershells and superwinds have been advocated as major mechanisms which re-distribute heavy elements throughout galactic discs and the intergalactic space. However, not until recently theoretical models did not care about the hot bubble chemical evolution and associated the hot gas metallicity with the metallicity of the surrounding ISM.

The enrichment caused by supernova products has been incorporated in the models by Silich et al. (2001). The calculations assume an instantaneous burst of star formation with a Salpeter initial mass function with an exponent $`\alpha `$ and incorporate the oxygen and iron yields YO,Fe derived from stellar evolution models. The mass of oxygen (and iron) released by the stellar cluster of total mass MSB with upper and lower cutoff masses Mup and Mlow, during the evolutionary time t is

$$M_{ej}(O,Fe)=\frac{(\alpha 2)M_{SB}}{M_{low}^{2\alpha }M_{up}^{2\alpha }}_{M_{}(t)}^{M_{up}}Y_{O,Fe}(m)m^\alpha dm.$$
(1)
The ejected metals were assumed to be effectively mixed with the stellar hydrogen envelopes and matter evaporated from the shock driven outer shell. The mean hot gas metallicity measured in Solar units is then

$$Z_{O,Fe}=\frac{M_{ej}(O,Fe)/Z_{}(O,Fe)+Z_{ISM}M_{ev}}{M_{ev}+M_{ej}},$$
(2)
where Z(O,Fe) are the oxygen and iron solar content, ZISM is the interstellar gas metallicity, Mej and Mev are the total ejected mass and the mass loaded by thermal evaporation. The results of the calculations are presented in Figure 1.

Figure 1a displays the extreme case without any mass loading and therefore represents the upper limit for the metallicity of a superbubble. Figure 1b represents a more realistic case, with dilution caused by the low metallicity ISM. In this case the oxygen abundance reaches its maximum Z$`{}_{O}{}^{}2`$Z after 6Myr of evolution and then falls slowly to the solar value after 10Myr to finally approach 0.1Z at the end of the starburst activity.

These results are in contradiction with the previous estimates of the X-ray emitting gas metallicities based on the analysis of the ROSAT data. However, as mentioned by Strickland & Stevens (1998) the low, subsolar metallicities most probably come from the single temperature model used to fit the ROSAT spectra. Another possible solution to the low abundance problem has been indicated by Breitschwerdt (2003), who provided a detailed analysis of the non-equilibrium ionization model.

The recent estimates of the hot gas metal abundances in the starburst galaxy NGC 1569 by Martin, Kobulnicky, & Heckman (2002), are based on two temperature model and show an excellent agreement with our theoretical predictions. The $`\alpha `$ element metallicities for a starburst with the age between 10Myr to 20Myr have been best fitted by Z$`{}_{\alpha }{}^{}=1.0`$Z, and abundances higher than solar have not been excluded.

The ultimate fate of metals released within the starburst regions has been intensively debated in recent years (see Kunth; this volume and references therein). It is highly dependent on the kinematic properties of the host galaxy (Silich & Tenorio-Tagle 2001). Supergalactic winds rapidly develop in fast rotating systems with a disc-like ISM distribution and are less likely in slow rotating galaxies with a thick-disc or a spherical gas distribution. The metallicity-flattening relation revealed by Barazza & Binggeli (2002) that shows that rounder dwarf ellipticals tend to be more metal-rich, may support these conclusions.

## 3. Supergalactic wind

Energetic starbursts, evolving in a flat disc-like ISM, are able to drive their associated shock waves to the outskirts of their host galaxies, leading to the development of supergalactic winds. It is usually assumed that within the region of star formation the matter ejected by strong stellar winds and supernova explosions is fully thermalized by means of random collisions. This generates the high central temperature and large overpressure that initiates the metal-rich gas outflow. There, the mean total energy LSB and mass ṀSB deposition rates control, together with the actual size of the star forming region RSB, the properties of the resultant outflow. After crossing r = RSB the gas is immediately accelerated by steep pressure gradients and rapidly reaches its terminal velocity Vt. This is due to a fast conversion of thermal energy, into kinetic energy of the resultant wind. The free wind analytic solution of Chevalier & Clegg (1985) assumes that the thermalized gas freely expands out of the star forming region and then adiabatically cools down approaching an r-4/3 temperature distribution.

### 3.1. The steady state solution

Here we study the true physical properties of such well developed free wind outflows, taking into consideration strong radiative cooling. Following Chevalier & Clegg (1985) we assume a spherically symmetric wind, unaffected by the gravitational pull caused by the central star cluster. The equations that govern the steady gas outflow away from the star forming region are:

$`{\displaystyle \frac{1}{r^2}}{\displaystyle \frac{\mathrm{d}}{\mathrm{d}r}}\left(\rho ur^2\right)=0,`$ (3)
$`\rho u{\displaystyle \frac{\mathrm{d}u}{\mathrm{d}r}}={\displaystyle \frac{\mathrm{d}P}{\mathrm{d}r}},`$ (4)
$`{\displaystyle \frac{1}{r^2}}{\displaystyle \frac{\mathrm{d}}{\mathrm{d}r}}\left[\rho ur^2\left({\displaystyle \frac{u^2}{2}}+{\displaystyle \frac{\gamma }{\gamma 1}}{\displaystyle \frac{P}{\rho }}\right)\right]=Q,`$ (5)
where r is the spherical radius, and $`u(r),\rho (r)`$ and $`P(r)`$ are the wind velocity, density and thermal pressure, respectively. Q is the cooling rate ($`Q=n^2\mathrm{\Lambda }`$) where $`n`$ is the wind number density and $`\mathrm{\Lambda }`$ is the cooling function.

### 3.2. The impact of radiative cooling

A first order of magnitude estimate of whether or not radiative cooling could affect the thermodynamics of superwinds, results from a comparison of the radiative cooling time with the characteristic dynamical time scale :

$$\tau _{cool}(r)=\frac{3kT}{n\mathrm{\Lambda }},\tau _{dyn}(r)=_{R_{SB}}^r\frac{\mathrm{d}r}{u(r)}.$$
(6)
Our calculations show that both the velocity and the density distributions remain practically unaffected by radiative cooling. However, if radiative cooling becomes efficient, the temperature distribution strongly deviates away from the adiabatic solution forcing the gas to soon reach temperature values of the order of 104 K (see Figure 2a).

Figure 2b implies that radiative cooling may become efficient for compact massive star clusters. If the initial wind parameters (LSB and RSB) intersect below the corresponding terminal velocity curve, cooling will be inefficient and deviations from the adiabatic solution would be negligible. However, if the initial wind parameters intersect above the corresponding terminal velocity curve, radiative cooling is expected to become important, causing a major impact on the appearance of superwinds.

### 3.3. The radiative mode observational implications

Cooling provides significant modifications of a superwind internal structure (see Figure 3), bringing the outer boundary of the X-ray emitting zone and inner boundary of the cold envelope closer to the star cluster. This in turn modifies the predicted X-ray and Hα luminosities. The X-ray emission vanishes at large distances and arises only from zones close to the star forming region, thus leading to a drop in the total X-ray luminosity.

The trend of the Hα emission is different. When cooling sets in, much denser superwind layers acquire temperatures below $`10^4`$K. These are to be ionized by the central star cluster and by the soft X-ray photons to become visible in the optical line regime. This causes an increase in the Hα luminosity. For compact star clusters the difference between the adiabatic and radiative model predictions may reach almost two orders of magnitude.

## 4. Conclusions

The results from the calculations presented here imply that:

The metallicity of superbubbles vary with time and can easily exceed the solar value even if the host galaxy ISM has a low metal abundance. The effects of metal contamination are most noticeable during the first 10 Myr of the bubble evolution.

Galactic superwinds driven by compact and powerful starbursts undergo catastrophic cooling close to the star cluster surface and establish a temperature distribution radically different to that predicted by adiabatic calculations.

The fall of the superwind temperature leads to a smaller zone radiating in X-rays and decreases the superwind X-ray luminosity.

At the same time cooling brings the inner radius of the warm ionized gas envelope closer to the star cluster. This increases the estimated Hα luminosity and predicts a low-intensity broad ($`1000`$ km s-1) line emission component.

#### Acknowledgments.

The authors feel honoured to have attended this special meeting and acknowledge financial support from CONACYT (México) grant 36132-E.

## References

Barazza, F.D. & Binggeli, B. 2002, Astro-ph/0209112

Breitschwerdt, D. 2003, in RevMexAA Conf. Ser. Winds, Bubbles and Explosions, (in press)

Chevalier, R.A. & Clegg, A.W. 1985, Nature, 317, 44

Martin, C.L., Kobulnicky, H.A. & Heckman, T.M., 2002, ApJ, 574, 663

Kunth, D. 2003, This volume

Silich, S. & Tenorio-Tagle, G., 2001, ApJ, 552, 91

Silich, S.A., Tenorio-Tagle, G., Terlevich R., Terlevich E. & Netzer H., 2001, MNRAS, 324, 191

Strickland, D.K. & Stevens, I.R. 1998, MNRAS, 297, 747

## Discussion

Casiana Muñoz-Tuñón: Do your results imply that cooling may affect the evolution of superbubbles ?

Silich: No, we have considered the late post-blowout stage, when the reverse shock reaches the outskirts of the galaxy and the outflow approaches a free-wind steady state regime.


text:
# Matching of the continuous gravitational wave in an all sky search

## 1 Introduction

The first generation of kilometer-scale gravitational wave (GW) laser interferometric detectors with sensitivity in the frequency band 10 Hz to few kHz and ultra cryogenic bar detectors sensitive at frequencies around 1 kHz will start collecting data soon. The TAMA 300 (Tsubona 1995) has already done the first large scale data acquisition (Tagoshi et al. 2001), while LIGO (Abramovici et al. 1992) and GEO600 (Danzmann 1995) have recently carried out their first science observations. VIRGO (Bradaschia et al. 1991) may become operational in couple of years. Also, an eighty meter research interferometer ACIGA (McCleland et al. 2000) near Perth, Australia is under construction, hoping that it may be possible to extend it to multi-kilometer scale in the future.

At present, majority of searches are focussed in the detection of chirp and burst signals. However interest for the search of continuous gravitational wave (CGW) signals from the output of detectors is growing (Jaranwoski, et. al 1998, Jaranowski & Królak 1999, 2000; Astone et al. 2002; Brady et. al. 1998, Brady & Creighton 2000) due to the possibility of enhancing the signal-to-noise ratio (S/N) by square root of observation time $`\sqrt{T_{obs}}`$. An optimistic estimate suggests that the Earth based laser interferometric detectors may detect such signals with an observation time of 1-yr.

The strength of the CGW largely depend on the degree of long-lived asymmetry in the source. There are several mechanism for producing such an asymmetry (Pandharipande et al. 1976; Bonazzola & Gourgoulhon 1996; Zimmermann & Szedenits 1979; Zimmermann 1980). The estimates of the asymmetry in neutron stars shows that the amplitude of CGW may be $`10^{25}`$. Hence, long integration will be required to get the signature of signal. But this in turn induces several other problems viz. Doppler modulation and non-stationarity of the detector noise. Consequently data analysis becomes more harder. However, Doppler modulation will provide the information of position of the source in the sky.

It has been realized that with present computing power, coherent all-sky full frequency search in the bandwidth (BW) of the detector of a month long data is computationally prohibitive. Some computationally efficient alternative approaches have been suggested viz. tracking and stacking (Brady and Creighton 2000; Schutz 1998). Tracking  involves the tracking of lines in the time-frequency plane built from the Fourier transform (FT) of one day long stretches of data while stacking involves dividing the data into day long stretches, searching each stretch for signals, and enhancing the detectability by incoherently summing the FT of data stretches. Also, accurate modeling of GW form, optimal data processing and efficient programming are an integral part of all sky-search.

The basic method to analyze the detector output to get the signature of GW signals rely on how efficiently one can Fourier analyze the data. Fourier analysis of the data has the advantage of incorporating the interferometers noise spectral density. The problem for the search of CGW largely depend on how accurately one can take into account the translatory motion of the detector acquired from the motions of the Earth in Solar system barycentre (SSB) frame. It has been shown (Srivastava and Sahay 2002a,b) that amplitude modulation will only redistribute the power of the frequency modulated (FM) signal in five frequency bands $`f\pm 2f_{rot}`$, $`f`$, $`f\pm f_{rot}`$, where $`f`$ and $`f_{rot}`$ are the frequencies of the FM signal and the rotational frequency of the Earth respectively. Hence it is sufficient to consider only FM signal for the analysis of the matching of signals from different locations in the sky .

The study of the matching of signals from locations in the sky that differ in their co-latitude and longitude by $`\pi `$ radians has been made for 1-d data set (Srivastava and Sahay 2002c). They observed symmetries in the sky locations. In this paper we report that the observed symmetries are not generic. We study this problem for longer observation time and the other parameters which affect the symmetries. In the next Section we briefly review the FT of FM signal. In Section 3, using the concept of fitting factor (FF), we investigate the matching of signals in an all sky-search with reference to the Earth based laser interferometric detector by considering the source location as the parameter of signal manifold and templates corresponding to different source locations for longer observation time. We present our conclusions in Section 4.

## 2 Fourier transform of the frequency modulated continuous gravitational wave

The time dependence of the phase of a monochromatic CGW signal of frequency $`f_o`$ observed at detector location is given as (Srivastava and Sahay 2002a)

$`\mathrm{\Phi }(t)`$ $`=`$ $`2\pi f_o[t+{\displaystyle \frac{R_{se}}{c}}\mathrm{sin}\theta \mathrm{cos}\varphi ^{}+{\displaystyle \frac{R_e}{c}}\mathrm{sin}\alpha \{\mathrm{sin}\theta (\mathrm{sin}\beta \mathrm{cos}ϵ\mathrm{sin}\varphi +\mathrm{cos}\varphi \mathrm{cos}\beta )+`$ (1)
$`\mathrm{sin}\beta \mathrm{sin}ϵ\mathrm{cos}\theta \}{\displaystyle \frac{R_{se}}{c}}\mathrm{sin}\theta \mathrm{cos}\varphi {\displaystyle \frac{R_e}{c}}\mathrm{sin}\alpha \{\mathrm{sin}\theta (\mathrm{sin}\beta _o\mathrm{cos}ϵ\mathrm{sin}\varphi +`$
$`\mathrm{cos}\varphi \mathrm{cos}\beta _o)+\mathrm{sin}\beta _o\mathrm{sin}ϵ\mathrm{cos}\theta \}]`$
$`=`$ $`2\pi f_ot+𝒵\mathrm{cos}(a\xi _{rot}\varphi )+𝒩\mathrm{cos}(\xi _{rot}\delta )𝒬`$
where

$$\begin{array}{ccc}𝒫\hfill & =& 2\pi f_o\frac{R_e}{c}\mathrm{sin}\alpha (\mathrm{cos}\beta _o(\mathrm{sin}\theta \mathrm{cos}ϵ\mathrm{sin}\varphi +\mathrm{cos}\theta \mathrm{sin}ϵ)\mathrm{sin}\beta _o\mathrm{sin}\theta \mathrm{cos}\varphi ),\hfill \\ 𝒬\hfill & =& 2\pi f_o\frac{R_e}{c}\mathrm{sin}\alpha (\mathrm{sin}\beta _o(\mathrm{sin}\theta \mathrm{cos}ϵ\mathrm{sin}\varphi +\mathrm{cos}\theta \mathrm{sin}ϵ)+\mathrm{cos}\beta _o\mathrm{sin}\theta \mathrm{cos}\varphi ),\hfill \\ 𝒩\hfill & =& \sqrt{𝒫^2+𝒬^2},\hfill \\ 𝒵\hfill & =& 2\pi f_o\frac{R_{se}}{c}\mathrm{sin}\theta ,\hfill \\ \hfill & =& 𝒵\mathrm{cos}\varphi ,\hfill \end{array}\}$$
(2)
$$\begin{array}{ccc}\delta \hfill & =& \mathrm{tan}^1\frac{𝒫}{𝒬},\hfill \\ \varphi ^{}\hfill & =& w_{orb}t\varphi ,\hfill \\ \beta \hfill & =& \beta _o+w_{rot}t,\hfill \\ \xi _{orb}\hfill & =& w_{orb}t=a\xi _{rot};a=w_{orb}/w_{rot}\mathrm{\hspace{0.33em}1}/365.26,\hfill \\ \xi _{rot}\hfill & =& w_{rot}t\hfill \end{array}\}$$
(3)
where $`R_e`$, $`R_{se}`$, $`w_{rot}`$ and $`w_{orb}`$ represent respectively the Earth’s radius, the average distance between the centre of Earth from the origin of SSB frame, the rotational and the orbital angular velocity of the Earth. $`ϵ`$ and $`c`$ represent the obliquity of the ecliptic and the velocity of light. $`\alpha `$ is the colatitude of the detector. Here $`t`$ represents the time in s elapsed from the instant the Sun is at the Vernal Equinox and $`\beta _o`$ is the local sidereal time at that instant, expressed in radians. $`\theta `$ and $`\varphi `$ denote the celestial colatitude and celestial longitude of the source. These coordinates are related to the right ascension, $`\overline{\alpha }`$ and the declination, $`\overline{\delta }`$ of the source via

$$\begin{array}{ccc}\hfill \mathrm{cos}\theta & =& \mathrm{sin}\overline{\delta }\mathrm{cos}ϵ\mathrm{cos}\overline{\delta }\mathrm{sin}ϵ\mathrm{sin}\overline{\alpha }\hfill \\ \hfill \mathrm{sin}\theta \mathrm{cos}\varphi & =& \mathrm{cos}\overline{\delta }\mathrm{cos}\overline{\alpha }\hfill \\ \hfill \mathrm{sin}\theta \mathrm{sin}\varphi & =& \mathrm{sin}\overline{\delta }\mathrm{sin}ϵ+\mathrm{cos}\overline{\delta }\mathrm{cos}ϵ\mathrm{sin}\overline{\alpha }\hfill \end{array}\}$$
(4)
The two polarisation states of the signal can be taken as

$$h_+(t)=h_{o_+}\mathrm{cos}[\mathrm{\Phi }(t)]$$
(5)
$$h_\times (t)=h_{o_\times }\mathrm{sin}[\mathrm{\Phi }(t)]$$
(6)
$`h_{o_+}`$, $`h_{o_\times }`$ are the time independent amplitude of $`h_+(t)`$, and $`h_\times (t)`$ respectively.

For the analysis of the matching of signals from different locations in the sky, it is sufficient to consider either of the polarisation states of FM signal. We consider the $`+`$ polarisation of the signal of unit amplitude. Therefore, the FT for a data of $`T_{obs}`$ observation time may be expressed as (Srivastava and Sahay 2002b)

$`\stackrel{~}{h}(f)`$ $`=`$ $`{\displaystyle _0^{T_{obs}}}\mathrm{cos}[\mathrm{\Phi }(t)]e^{i2\pi ft}𝑑t`$ (7)
$`{\displaystyle \frac{\nu }{2w_{rot}}}{\displaystyle \underset{k=\mathrm{}}{\overset{k=\mathrm{}}{}}}{\displaystyle \underset{m=\mathrm{}}{\overset{m=\mathrm{}}{}}}e^{i𝒜}[\stackrel{~}{𝒞}i\stackrel{~}{𝒟}];`$
where

$$\begin{array}{ccc}\nu \hfill & =& \frac{f_of}{f_{rot}}\hfill \\ 𝒜\hfill & =& \frac{(k+m)\pi }{2}𝒬\hfill \\ \hfill & =& \frac{J_k(𝒵)J_m(𝒩)}{\nu ^2(ak+m)^2}\hfill \\ \stackrel{~}{𝒞}\hfill & =& \mathrm{sin}\nu \xi _o\mathrm{cos}(ak\xi _o+m\xi _ok\varphi m\delta )\hfill \\ & & \frac{ak+m}{\nu }\{\mathrm{cos}\nu \xi _o\mathrm{sin}(ak\xi _o+m\xi _ok\varphi m\delta )+\mathrm{sin}(k\varphi +m\delta )\}\hfill \\ \stackrel{~}{𝒟}\hfill & =& \mathrm{cos}\nu \xi _o\mathrm{cos}(ak\xi _o+m\xi _ok\varphi m\delta )\hfill \\ & & +\frac{ka+m}{\nu }\mathrm{sin}\nu \xi _o\mathrm{sin}(ak\xi _o+m\xi _ok\varphi m\delta )\mathrm{cos}(k\varphi +m\delta )\hfill \\ \xi _o\hfill & =& w_{rot}T_{obs}\hfill \end{array}\}$$
(8)
J stands for the Bessel function of first kind. The computational strain to compute $`\stackrel{~}{h}(f)`$ from equation (7) can be reduced by $`50\%`$ by using the symmetrical property of the Bessel functions, given as

$`\stackrel{~}{h}(f)`$ $``$ $`{\displaystyle \frac{\nu }{w_{rot}}}[{\displaystyle \frac{J_o(𝒵)J_o(𝒩)}{2\nu ^2}}[\{\mathrm{sin}(+𝒬)\mathrm{sin}(+𝒬\nu \xi _o)\}+`$ (9)
$`i\{\mathrm{cos}(+𝒬)\mathrm{cos}(+𝒬\nu \xi _o)\}]+`$
$`J_o(𝒵){\displaystyle \underset{m=1}{\overset{m=\mathrm{}}{}}}{\displaystyle \frac{J_m(𝒩)}{\nu ^2m^2}}\left[(𝒴𝒰𝒳𝒱)i(𝒳𝒰+𝒴𝒱)\right]+`$
$`{\displaystyle \underset{k=1}{\overset{k=\mathrm{}}{}}}{\displaystyle \underset{m=\mathrm{}}{\overset{m=\mathrm{}}{}}}e^{i𝒜}(\stackrel{~}{𝒞}i\stackrel{~}{𝒟})];`$
$$\begin{array}{ccc}𝒳& =& \mathrm{sin}(+𝒬m\pi /2)\hfill \\ 𝒴& =& \mathrm{cos}(+𝒬m\pi /2)\hfill \\ 𝒰& =& \mathrm{sin}\nu \xi _o\mathrm{cos}m(\xi _o\delta )\frac{m}{\nu }\left\{\mathrm{cos}\nu \xi _o\mathrm{sin}m(\xi _o\delta )\mathrm{sin}m\delta \right\}\hfill \\ 𝒱& =& \mathrm{cos}\nu \xi _o\mathrm{cos}m(\xi _o\delta )+\frac{m}{\nu }\mathrm{sin}\nu \xi _o\mathrm{sin}m(\xi _o\delta )\mathrm{cos}m\delta \hfill \end{array}\}$$
(10)
Equation (9) contains double infinite series of Bessel function. However, we know that the value of Bessel function decreases rapidly as its order exceeds the argument. From equations (2), it can be shown that the order of Bessel function required to compute $`\stackrel{~}{h}(f)`$ in the infinite series are

$$k3133.22\times 10^3\mathrm{sin}\theta \left(\frac{f_o}{1\text{kHz}}\right),$$
(11)
$$m134\left(\frac{f_o}{1\text{kHz}}\right).$$
(12)
## 3 Matching of the continuous gravitational wave

Matched filtering is the most suitable technique for the detection of signals whose wave form is known. The wave forms are used to construct a bank of templates, which represent the expected signal wave form with all possible ranges of its parameters (source location, ellipticity, etc.). In an all sky search by match filtering, Srivastava and Sahay 2002c observed symmetries under the following transformations

$$\theta _T\pi \theta _T0\theta _T\pi ,$$
(13)
$$\varphi _T\pi \varphi _T0\varphi _T\pi ,$$
(14)
$$\varphi _T3\pi \varphi _T\pi \varphi _T2\pi .$$
(15)
The above observations are made by studying few cases for 1-d data set. Hence, it will be important to understand and check the generic nature of the matching of signals under the above transformations. To study this problem we use the formula for FF (Apostolatos 1995) which quantitatively describes the closeness of two signals, given as

FF $`=`$ $`{\displaystyle \frac{h(f)|h_T(f;\theta _T,\varphi _T)}{\sqrt{h_T(f;\theta _T,\varphi _T)|h_T(f;\theta _T,\varphi _T)h(f)|h(f)}}}`$ (16)
where $`h(f)`$ and $`h_T(f;\theta _T,\varphi _T)`$ represent respectively the FTs of the actual signal wave form and the templates. The inner product of two waveform $`h_1`$ and $`h_2`$ is defined as

$`h_1|h_2`$ $`=`$ $`2{\displaystyle _0^{\mathrm{}}}{\displaystyle \frac{\stackrel{~}{h}_1^{}(f)\stackrel{~}{h}_2(f)+\stackrel{~}{h}_1(f)\stackrel{~}{h}_2^{}(f)}{S_n(f)}}𝑑f`$ (17)
$`=`$ $`4{\displaystyle _0^{\mathrm{}}}{\displaystyle \frac{\stackrel{~}{h}_1^{}(f)\stackrel{~}{h}_2(f)}{S_n(f)}}𝑑f`$
where  denotes complex conjugation, $`\stackrel{~}{}`$ denotes the FT of the quantity underneath $`(\stackrel{~}{a}(f)=_{\mathrm{}}^{\mathrm{}}a(t)\text{exp}(2\pi ift)𝑑t)`$ and $`S_n(f)`$ is the spectral density of the detector noise. In our analysis we assume the noise to be stationary and Gaussian. It is remarked that to compute the inner product of the signals one would require the BW of the Doppler modulated signal. From the Doppler shift (Srivastava and Sahay 2002a), the BW of Doppler modulated signal may be given as

$$\text{BW}(1.99115\times 10^4\mathrm{sin}\theta +3.09672\times 10^6)f_o.$$
(18)
### 3.1 Celestial colatitude

Let us consider that GEO600 detector (the position and orientation of the detectors can be found in Jaranowski et al. 1998) receive a CGW signal of frequency $`f_o=0.1`$ Hz (unreasonably low frequency has been chosen for illustrative purposes limited by available computational power) from a source located at $`(\theta ,\varphi )=(25^o,20^o)`$. In order to evaluate the matching of the signals in colatitude, we first maximize FF over $`\varphi `$ by choosing $`\varphi =\varphi _T=20^o`$. Now, we wish to check the symmetries in colatitude represented by equation (13) for the data set $`T_{obs}=120`$ d. For the purpose we maximize the FF over $`\theta `$ by varying $`\theta _T`$ in discrete steps over entire range i.e. $`0^o`$ to $`180^o`$. For the present case it is sufficient to take the ranges of $`k`$ and $`m`$ as $`1`$ to $`345`$ and $`3`$ to $`3`$ respectively and $`\text{BW}=20.1954\times 10^6\text{Hz}`$. The results so obtained are shown in Figure (1). To establish the observed symmetries, we similarly compute the FF for $`T_{obs}=1,2,3\mathrm{}.365`$ d and observe that the matching of signals remains almost same. However, due to obliquity of the ecliptic, the variation in the matching of signals will depend on source frequency, colatitude and detector position and orientation. We check the dependence of the FF on this parameters. The result so obtained for different Earth based laser interferometric detectors are shown in the Tables (1), & (2) and Figures (2) & (3).

The analysis of the matching of signals in $`\theta `$ space shows that

1. for fixed $`f_o`$, the FF

1. is independent of $`T_{obs}`$ and $`\varphi `$.
    2. does not vary significantly with the variation of source location, detector position and orientation \[Tables (1), and (2)\].
2. the FF falls with the source frequency \[Figures (2) and (3)\]. From the figure we find that it may not be relevant to take care of the symmetries in the sky locations for the search of CGW from the output of LIGO-I, GEO600 and TAMA detectors whose lower cut off frequency is 40/75 Hz (Owen and Sathyaprakash 1999).
3. the approximate fall of FF based on the figures (2)and (3) may be given as

$$\text{FF}=A_o+A_1\left(\frac{f_o}{\text{Hz}}\right)A_2\left(\frac{f_o}{\text{Hz}}\right)^2+A_3\left(\frac{f_o}{\text{Hz}}\right)^3$$
(19)
where $`A_o`$, $`A_1`$, $`A_2`$, $`A_3`$ are constants given in Table (2).

### 3.2 Celestial longitude

The Doppler shift due to the motions of the Earth is mainly depend on the colatitude and source frequency and have very less dependence on the longitude. Consequently, grid spacing of the templates for matched filtering for an all sky search will insignificantly depend on longitude (Brady & Creighton 2000). Keeping this in view we similarly check the matching of signals under the transformation given by the equation (14). We chosen the LIGO detector located at Livingston, selected a data set $`T_{obs}=120`$ d, $`(\theta ,\varphi )=(0.5^o,40^o)`$, $`f_o=5`$ Hz. In order to compute the FF, we first maximize equation (16) over $`\theta `$ by selecting $`\theta =\theta _T=0.5^o`$, followed by maximization over $`\varphi `$ in discrete steps over its entire range, $`0^o\varphi 360^o`$. The result so obtained is shown in Figure (4). We also check the mismatch of the signals for different $`\theta ,\varphi `$ and $`f_o`$ by computing the FF for the data set of $`T_{obs}=1,2,\mathrm{}.25/100`$ d. The results so obtained are shown in Figures (5), (6), and (7) respectively. Almost same behavior has been observe for the transformation represented by the equation (15).

From the figures we note that the matching of the signals in longitude decreases with $`T_{obs},f_o,\theta `$ and $`\varphi `$. However, behavior of the matching of signals are similar in nature and may be represented by

$$\text{FF}=B_oB_1T_{obs}+B_2T_{obs}^2B_3T_{obs}^3+B_4T_{obs}^4;T_{obs}=1,2,\mathrm{}.25/100\text{d}$$
(20)
Where, $`B_o`$, $`B_1`$, $`B_2`$, $`B_3`$, $`B_4`$ are the constants given in Table (3). The equations (19) and (20) does not represents the oscillatory part of the figures. This oscillatory behavior is more or less typical in waveform that match well or bad depending on their parameters and is low enough in comparison to the threshold of the detection of GW that one may choose.

## 4 Conclusion

In view of blind all sky search for CGW, we have studied the matching of signals from different source locations assuming the noise to be stationary and Gaussian. For fixed $`f_o`$, we observed that the matching of signals from locations in the sky that differ in their colatitude by $`\pi `$ radians is independent of $`T_{obs}`$, and $`\varphi `$. However, it falls with $`f_o`$. But in longitude the matching of signals falls with $`T_{obs}`$, $`\theta `$, $`\varphi `$, $`f_o`$. We believe that the matching of signals will increase for the real data. This is due to the resolution provided by the Fast Fourier transform (FFT). However, it may not be relevant to account this symmetries for the search of CGW from the output of LIGO-I, GEO600 and TAMA detectors.

This analysis will be more relevant if one performs hierarchical search (Mohanty & Dhurandhar 1996; Mohanty 1998). This search is basically a two step search, in first step the detection threshold is kept low and in the second step a higher threshold is used. The higher threshold is used for those templates which exceed the first step threshold. The study of the matching of signals in different source locations will be also relevant for Laser interferometer space antenna (LISA) (Hough J., 1995). The work has been initiated and may be useful for the data analysis of CGW.

## Acknowledgment

I am thankful to IUCAA for providing hospitality where major part of the work was carried out. I am thankful to Prof. S.V. Dhurandhar, IUCAA, Pune for his useful discussions and Mr. A. Sengupta, IUCAA, Pune for reading the paper.


text:
# Beam Tests of Ionization Chambers for the NuMI Neutrino Beam

## I Introduction

The Neutrinos at the Main Injector (NuMI) beamline at the Fermi National Accelerator Laboratory  will generate an intense $`\nu _\mu `$ beam from the decays of mesons produced in the collisions of 120 GeV protons in a graphite target. The mesons are focused by magnetic “horns” into a 675 m evacuated volume to allow decays to neutrinos. A downstream Aluminum/Steel absorber and bedrock absorb the remnant hadrons and muons in the beam, leaving only neutrinos. The facility is expected deliver beam to neutrino experiments, beginning with the MINOS neutrino oscillation experiment, starting in early 2005.

The meson decays $`\pi /K\mu \nu _\mu `$ produce an energetic muon for every neutrino, allowing monitoring and validation of the neutrino beam focusing to be accomplished by monitoring of the muon flux. As in several previous experiments, the muon flux and remnant hadron flux at the end of the decay volume will be measured by arrays of ionization chambers .

The beam monitoring system will measure the intensity and spatial distribution of the hadron beam at the end of the decay tunnel, upstream of the absorber, and of the muon beam after the absorber and at several stations in the bedrock. The monitoring system will consist of 2 m $`\times `$ 2 m arrays of ionization chambers with 25 cm inter-chamber spacing, with one chamber array in each of the above stations. The hadron and muon fluxes are measures of any targeting or focusing failures. The peak charged particle fluxes in one 9 $`\mu `$s accelerator burst will be 2000, 25, 3, and $`1.5\times 10^6`$/cm2 in the four monitoring stations.

Each ionization chamber will measure the flux of charged particles by using an applied electric field to collect the ionization created in a helium gas volume. The charge measured from each chamber will be proportional to the charged particle flux at that location. By using an array of chambers the spatial distribution of beam intensity can be inferred. While operated without gas amplification the signal from the intense NuMI beam in one 8 cm $`\times `$ 8 cm ionization chamber will be 33000, 1400, 170, and 83 pC at the four stations.

The individual ionization chambers within each array are parallel plate chambers made up of two 4” square ceramic plates with Ag-Pt electrodes. One plate has a single electrode that applies HV bias. The second plate has two electrodes: a central square sense pad measuring 3” $`\times `$ 3”, surrounded by a 1 cm guard ring. The sense pad is connected into the electronics which provides a virtual ground. The guard ring is grounded. The chamber gas is pure helium at atmospheric pressure.

The major limitation of ionization chambers used as beam monitors has been space charge build up inside the chamber. Intense particle fluxes release sufficient ionized charge in the chamber gas so as to create a reverse electric field inside the chamber. With the reduced net field in the chamber, ions require a longer time to reach the collection electrodes, and hence suffer more recombination loss in the gas. Recombination losses increase at larger particle fluxes, resulting in a non-linear performance of the ionization chamber beam monitor at large intensities.

## II Booster Beam Test

A beam test of prototype ionization chambers was undertaken at the Fermilab Booster accelerator, which delivers up to 1012 protons/cm2 of 8 GeV in a 1.56 $`\mu `$sec spill. Tested were two chambers, one with a 1 mm electrode spacing, the other with 2 mm, with continuos gas flow. We studied the shape of the ionization vs. voltage plateau curve at several intensities and the linearity of the chamber response vs. beam intensity at several applied voltages.

For our beam test we placed two ionization chambers in the beamline. The chambers were housed in a stainless steel vessel with .005” Ti beam entrance and exit windows. Electrical feedthroughs were made with stainless steel compression fittings and PEEK plastic.

Two gas mixtures were used in the beam test. The primary gas was pure helium supplied from a cylinder with 99.998% purity. The other gas was a mixture consisting of 98% Helium and 2% Hydrogen, with $``$ 20 p.p.m. of impurities. The gas was first passed through a getter and gas analyzer. Online measurement indicated impurites of $`<`$ 1.5 p.p.m. All gas seals were metallic consisting of compression fittings or copper gaskets compressed by vacuum fittings.

Upstream of our ion chamber vessel Fermilab provided a secondary emission monitor (SEM), to locate the beam. The SEM provided targeting information upstream of the chamber, but was not usable in the analysis because the beam diverged after passing through the SEM, and it was only capable of measuring either the horizontal or vertical profile at any given time. Furthermore, the SEM was insensitive below $`2\times 10^{11}`$ protons per pulse. During the portion of datataking above $`2\times 10^{11}`$ the SEM indicated constant spot size.

A beam toroid was the primary method of measuring the beam intensity delivered to the apparatus. The toroid was supplied with an amplifier and ADC whose least significant bit was $`5\times 10^9`$ protons.

Rigidly attached on the outside of the ion chamber vessel were two beam profile chambers fashioned out of G-10 circuit board and epoxy. Each chamber had a segmented signal electrode composed of 1$`\times `$10 cm2 strips. One chamber provided the vertical profile, the other the horizontal profile. The profile chambers were the primary method for determining beam size. Gas flow was the same as that passed through the vessel. The profile chambers indicated constant beam spot size $``$5 cm2.

The signal from each of the ion chambers and beam profile chambers were read out into a charge integrating amplifier and then into an ADC . The electronics were triggered on an accelerator clock signal shortly before the beam pulse. The charge integration time could be altered and taken between beam spills, allowing measurement of pedestals and backgrounds.

## III Experimental Results

Tests of the chambers consisted of two complementary measurements. The first held the beam intensity constant, while varying the voltages applied to the chambers. The second held the applied voltages constant while adjusting the beam intensity.

The results of the first test are displayed in Figures 2 and 3. Each chamber is exposed to several beam intensities and the voltage varied 0-350 V. The ratio of collected charge to the measured beam intensity is plotted as a function of applied voltage.

The ideal voltage plateau curve would consist of a quick rise to a constant charge collected per proton, independent of voltage and intensity. This constant charge collected would be equal to the amount of charge liberated in the gas per proton. At higher voltages gas amplification is expected and charge collected per proton would increase above the plateau. From the height of the plateau on the 1 mm chamber one may infer a charge ionized per proton of 2.6 pC/107 protons, or 1.6 electron-ion pairs per proton incident on a 1 mm gas gap.

The drifting ions and electrons inside the chamber establish their own electric fields and the relatively slow drift velocities of the ions create a net space charge in the gas. This space charge induced electric field screens the electrodes, slowing the transit of ions and electrons across the electrode gap. As discussed extensively in , operating ion chambers at very large particle fluences modifies the voltage plateau curve discussed above. If the speed of the charges is sufficiently slowed, recombination may take place. This recombination loss is especially evident $`<`$150 V, where the lower voltages result in slower ion drift velocities and longer ion transit times. This loss also increases at higher beam intensity or in the larger 2 mm gap chamber, where space charge buildup should be worse.

Space charge effects were evident at all of the intensities delivered by the Booster. In Figure 2 all of the curves, except possibly one, show a slope in the plateau region, suggesting that there is little or no region where the charge is collected without loss or gain. The only useful voltage looks to be 130-190 V for the 1 mm chamber, where the curves all intersect. In Section IV we analyze this in terms of competing recombination and multiplication effects in the gas.

Figure 4 displays the results of the second test performed in our beam test for the case of the 1 mm chamber in Helium. The charge collected from the ion chamber is plotted as a function of beam intensity for several applied voltages. Here, ideal curves would all join at low intensity. Recombination loss, especially visible for lower voltages, cause the curves to fall off at higher intensity. The linear region is the operating range for the chamber where charge collected is proportional to incident flux.

An analysis of the 200 V Helium curve up to $`20\times 10^{10}`$protons/spill gives a good linear fit ($`\chi ^2`$/NDOF = 456/625). However, the intercept is less than zero: (-0.05 $`\pm `$ 0.012) $`\times 10^3`$ pC. The other data could not be satisfactorily fit to a linear dependence for beam intensities $`>10^{10}`$ protons/spill.

The nonzero intercept in the 200 V data is more apparent when the curve is plotted as a ratio, as in Figure 5 where the ratio of charge collected to beam intensity is plotted versus intensity. Here a straight line with zero intercept in the previous plot corresponds to a horizontal line. In Figure 5 the 200 V curve is obviously curved upward then downward as the intensity is increased from 1 to $`20\times 10^{10}`$ protons/spill. This rise is responsible for the negative intercept in the linear fit. In Section IV we discuss how the interplay of multiplication and recombination can cause such an effect. Similar to the voltage plateau data, the data in Figure 5 indicate a collected charge of $``$3pC/107protons in the 1 mm chamber, or 18 electron-ion pairs per cm in the He.

The similarity between the Helium and Helium-Hydrogen data is interesting. Gas additives with lower excitation potentials are commonly used to increase the drift velocity of electrons , which might naively be expected to reduce recombination losses. However, the recombination loss is not substantially changed by the addition of H2, and in Section IV we discuss how the ion drift more strongly affects space charge build-up, hence recombination loss, in the chamber.

Another parameter that might be expected to vary as a result of the additive is the ionization per proton. The incident protons excite metastable states in the Helium, which are followed by collision of Helium atoms with the H2 additive, yielding strong ionization of the H2 . That the observed ionization/proton is not significantly increased by the H2 additive is perhaps indicative that even the minute impurities ($`<`$1.5 p.p.m.) in the nominally pure Helium gas are sufficient to increase the ionization. The observed 1.6-1.8 ionizations/mm per proton in our data is notably higher than the 0.8-1.0 ionizations/mm inferref from the dE/dx of fast charged particles and the $`w`$=42 eV/ionization  for the purest Helium. Thus, the nominally pure Helium used in typical chamber applications is likely effectively doped.

## IV Simulation

In this section we model the data collected at the Booster beam test using a computer simulation of the pulse development in an ion chamber. The simulation incorporates the known drift of electrons and of Helium ions in electric fields, the effects of volume recombination of charges, and of gas amplification. The simulation follows the charges during and after the beam spill, recording the net charge collected at the chamber electrodes. The simulations indicate a complex interplay between gas amplification and charge recombination as a result of space charge build up. This calculation will then be used to extrapolate to the particle fluxes and 8.6$`\mu `$s spill duration expected in NuMI. Only Helium gas was simulated as literature data for the gas properties of Helium-Hydrogen are minimal, and the beamtest data were not significantly different for the two gases

### IV-A Simulation Method

The differential equations governing the charge flux and electric field evolution inside the ion chamber are nonlinear. Some earlier work analyzed special cases.  demonstrated that space charge accumulation in a continuously-ionized chamber causes the formation of a “dead region” of no electric field, and makes the approximation that all charges in the dead region are lost due to recombination. In , an ion chamber ionized over a short duration is studied, but again the assumption is made that all charges in the dead zone are lost to recombination. That assumption is not valid for our experiment where the beam pulse is of short duration allowing the dead zone to disappear before the charge completely recombines.

The present work considers the time dependant case where ionization is delivered in short pulse of 1.6 $`\mu `$s and allowed to drift out of the chamber. Our simulation is a finite element calculation of one spatial dimension, such that there are series of infinite planes of charge between the electrodes. The electrons and ions are drifted with the velocites discussed in Section IV-B. Space charge is calculated at each step and an image charge is induced such that the potential difference betweeen the electrodes is maintained at the applied voltage.

Simulated electric field distributions as a function of time are shown in Figures 6, 7, and 8. The figures show the field development using only the charge transport, and ignore charge recombination or multiplication.

The electric field simulated in Figure 6 is calculated at an intensity of 1010 ionizations/cm3/$`\mu `$s, where the excess of ions slightly warps the field. In Figure 7, the ion excess is much greater, to the point where it entirely screens the anode from the applied field, creating a dead zone. As ionization continues in the dead zone the ion and electron densities there increase because the charges are effectively trapped. When the beam spill ends the dead zone slowly fades away. Charge is able to escape from the edge of the dead zone because the dead zone is not created by the charge inside, but the excess of ions outside of it.

Figure 8 is calculated at the same intensity as Figure 7, but for the 2 mm gap chamber. The half of the chamber close to the anode has almost exactly the behavior of the 1 mm case. This similarity is a result of the complete screening of the electric fied by space charge which makes the rest of the chamber in the dead zone irrelevant, until the charge has a chance to excape after the ionization period (the beam spill) ends.

### IV-B Gas Properties

The parameters of charge drift, loss, and amplification used in this simulation have been taken from several literature sources. The electron drift velocities as a function of electric field are taken from . For the drift of Helium, we use the mobility for the He$`{}_{}{}^{+}{}_{2}{}^{}`$ ion, $`\mu _{He_2^+}=20`$ cm2/(V$``$s) , following the discussion of  which notes that above pressures of a few Torr, He+ ions tend to collide in the gas and form molecular ions. Previous work  have cited 10 cm2/(V$``$s) as the relevant mobility for He.

We simulated charge loss through two- and three-body volume recombination of electrons ions.  and  parameterize recombination of the form: $`dn_{}/dt=dn_+/dt=rn_{}n_+`$, with $`n_+`$ ($`n_{}`$) the positive (negative) ion densities in the gas, and the measured recombination coefficient $`r=2.4\times 10^8`$ cm3/(ion$``$s).

We did not simulate charge loss through the process of electron attachment to electronegative impurities and subsequent recombination. Using the data from  and  we find that the attachment of electrons to Oxygen in the gas is insignificant for impurity levels less than 30 p.p.m. (to be compared with 1.5 p.p.m. observed in our chambers). Attachment is negligle over the drift region and occurs primarily in the dead zone where the electrons persist for a greater time. There, the electrons will attach with a time constant $`\tau =(67/x)`$ ms, where $`x`$ is the impurity level in p.p.m. The electrons persist in the dead zone time of order 10 $`\mu `$s and the impurities in our beam test were of order 1 p.p.m giving about 0.1% attachment.

Charge amplification in the chamber gas is modeled via $`dN/dx=N\alpha `$, with the Townsend coefficient $`\alpha `$ modelled via $`\alpha /P=A\mathrm{exp}[B/(X/P)]`$, where $`P`$ is the chamber gas pressure, $`X`$ is the applied electric field, and the parameters $`A=3`$ ion pairs/(cm-Torr) and $`B=25`$ V/(cm-Torr) . This parameterization is taken from a source with similar gas purity, which is of importance because the purest Helium gas has significantly lower Townsend coefficient but is sensitive to even a few p.p.b. impurity level .

The ionization rate in the simulation is quoted in units of ionizations/cm3/$`\mu `$s in order to factor out the question of how many ionizations are created by an incident proton. Given $`dE/dx`$ and $`W`$, the average energy loss necessary to create one ionization in the gas, for Helium, values of 8-16/cm have been derived , but this too is sensitive to impurity level. For comparison to our data from the Booster beamtest, 1 ion./cc/$`\mu `$s corresponds to 0.6 protons/spill if the beam spot size of 5 cm2 and $``$16 ionizations/cm in Helium are used.

### IV-C Simulation Results

We found that the simulated space charge accumulation causes not only recombination losses, but the high fields also cause significant multiplication. This result can be seen in Figures 9 and 10 where the charge collected per ionization can actually increase with ionization, until recombination takes over. This should be compared to Figure 5. Both the data and the simulation indicate such an effect. This feature of the curve makes the fit line have a negative intercept, as stated earlier. The lower voltage required in the simulation to generate this effect suggests that multiplication may be overstated in simulation.

The simulation was used to generate voltage plateau curves as well. The results, shown in Figures 11 and 12, can be compared to Figures 2 and 3, respectively. The behavior is marked by greater losses at high intensity and low voltage, and also by a crossing point where the curves all approach each other. The crossing point occurs when the charge collected equals the charge liberated in the gas. In the simulation this occurs at $``$140 V and 210 V for the 1 and 2 mm chambers. In the data, these points occur at 175 V and and 190 V. Another feature is that higher intensities can show gain at lower voltages, resulting in the curves diverging after the crossing point.

The behavior observed in the data is nominally reproduced by the simulation. The crossing point occurs at a higher voltage than in the simulation of the 1 mm chamber, again suggesting that multiplication is overstated there. The crossing point can be moved left and right by simply changing the value of the Townsend coefficient $`\alpha `$, and is particularly sensitive to the parameter $`B`$. Hence the data can constrain the value of $`\alpha `$, and possibly the recombination coefficient $`r`$.

## V Conclusion

We have performed a beam test of Helium- and Helium-Hydrogen-filled ionization chambers at the Fermilab Booster accelerator. We have compared the experimental results to our own calculation of the expected charge collection in such chambers. While we have yet to extrapolate our calculations to the anticipated NuMI beam environment using the constraints of our beam test data, several effects of interest are observed. First, calculation of the expected ionization per charged particle in Helium based $`dE/dx`$ loss and the value $`w=42`$ eV from  may ignore $``$p.p.m. impurities in most chamber gas which result in additional ionization. Second, similar to previous work, our results indicate that space charge effects induce recombination losses in the collected charge at high particle fluxes; however, our results indicate additionally that gas amplification occurs in the ion chambers at high beam intensities due to the large space charge build-up of the electric field. The gas amplification gains compete with the losses due to recombination, effectively extending the range of linear response of the ion chamber with respect to beam intensity.

## VI Acknowledgements

We gratefully acknowledge the excellent facilities and technical help provided by the Fermilab Booster Radiation Damage Facility group, J. Lackey, M. Ferguson, T. Sullivan as well as B. Webber, Fermilab Proton Source Department Head.