text:
# 1 Introduction

## 1 Introduction

In continuum QCD, the conventional explanation for the smallness of the mass of pseudoscalar $`\pi `$ mesons is the following: QCD with $`N_f`$ massless quark flavours has a global chiral $`U(N_f)\times U(N_f)`$ symmetry, which, spontaneously broken, reduces to $`U(N_f)`$ and yields $`N_f^2`$ Goldstone bosons. Explicit breaking of the original chiral symmetry by a small quark mass renders these $`N_f^2`$ Goldstone bosons massive with correspondingly small mass. To accord with nature, one of these Goldstone bosons (the $`\eta `$-particle in the case $`N_f=2`$) has to acquire an an additional mass. To explain this was known as the $`U(1)`$ problem in the continuum. Its resolution there comes from the axial anomaly, whereby the axial symmetry corresponding to the $`U(1)`$ subgroup of $`U(N_f)`$ is explicitly broken by a quantum effect, reducing the number of Goldstone bosons to $`N_f^21`$.

Naive lattice regularization of such a fermionic theory is hindered by the doubling problem, namely that a return to the continuum manifests too many fermionic degrees of freedom. This doubling problem is resolved by the usage of Wilson fermions. However the extra Wilson term that removes the fermion doublers breaks chiral symmetry explicitly. This effect can be traced back to the existence of the axial anomaly in the continuum. For this reason the staggered format has often been the favoured one for the study of models with chiral symmetry breaking .

The Wilson action for free fermions in terms of dimensionless fermionic fields $`\psi (n)`$ defined at the sites $`n`$ of a $`d`$-dimensional lattice is

$$S_F^{(0)}[\overline{\psi },\psi ]=\frac{1}{2\kappa }\underset{n}{}\overline{\psi }(n)\psi (n)\frac{1}{2}\underset{n,\mu }{}\left[\overline{\psi }(n)(r\gamma _\mu )\psi (n+\mu )+\overline{\psi }(n+\mu )(r+\gamma _\mu )\psi (n)\right],$$
(1.1)
where

$$1/2\kappa =\widehat{M}_0+dr.$$
(1.2)
Here, $`\kappa `$ is the hopping parameter, $`r`$ is the Wilson parameter, $`a`$ is the lattice spacing and $`\widehat{M}_0`$ is a dimensionless fermion bare mass parameter. We use $`d=2`$ and $`r=1`$ throughout. This free fermion model is the weak coupling limit of an interactive theory in which a bare parameter $`g`$ measures the coupling of the free theory to some interaction. Even if $`\widehat{M}_0=0`$, now, the Wilson term contributes to the hopping parameter and there is no obvious chiral symmetry. The question arises \- what is the status of the chiral phase transition and the $`U(1)`$ problem on the lattice?

Despite the lack of an obvious chiral symmetry, there exists a host of numerical and analytical evidence for the existence of massless pions in the lattice formulation of QCD. These are believed to exist on a critical line $`\kappa _c(g)`$. In the literature, there are two explanations for the existence of the critical line and the masslessness of the lattice pions.

The first of these was sometimes referred to as the conventional explanation . Although there is no obvious chiral symmetry at non-zero $`r`$, the conventional explanation suggests that tuning $`\kappa `$ effects its recovery in some unknown way. Now, with chiral symmetry recovered at $`\kappa _c(g)`$, the same arguments as in the continuum may be applied.

The second explanation was first forwarded in 1984 by Aoki . Here it is accepted that since there is no chiral symmetry in the lattice formulation of QCD, its spontaneous breaking cannot be responsible for the masslessness of pions. Instead there is an Ising-like second order parity breaking phase transition. In the single flavour case the order parameter for parity symmetry is $`\overline{\psi }i\gamma _5\psi `$, the operator corresponding to the single $`\pi `$ meson. The parity symmetric phase is where $`\pi =\overline{\psi }i\gamma _5\psi =0`$. there is also a phase with long range order where $`\pi 0`$. At the transition between these phases, a correlation length $`\xi `$ diverges. This correlation length is identified as the inverse of the pion mass, which, hence, becomes zero on the phase boundary. Thus the pion is not a Goldstone boson in the Wilson lattice formulation. Aoki also recovered the current algebra relation between pion and quark mass ($`m_\pi ^2m_q(\kappa \kappa _c)`$) by considering the effective meson theory as a scalar field theory in four dimensions with mean field like critical behaviour. In the multiflavour case the parity symmetry breaking is accompanied by a flavour symmetry breaking and, with it, Goldstone bosons in the form of the charged pions. The $`\eta `$ remains massive according to Aoki’s analysis, and the $`U(1)`$ problem in the lattice successfully resolved .

Two main features distinguish Aoki’s QCD phase diagram from the conventional one. Firstly, the existence of the phase transition in Aoki’s scenario is due to parity symmetry breaking as opposed to chiral symmetry breaking in the conventional picture. The order parameter is $`\overline{\psi }i\gamma _5\psi `$ rather than $`\overline{\psi }\psi `$ . Secondly, instead of a single critical line extending from the strongly coupled limit $`g=\mathrm{}`$ to $`\kappa =1/2d`$ in the weakly coupled limit $`g=0`$, Aoki’s picture involves the existence of two such lines and (in QCD) five critical points linked by four cusps in the weakly coupled zone.

Aoki’s QCD phase diagram is based on infinite-volume analyses in the limits of strong and weak coupling and on an analogy to the Gross-Neveu model, which, except for confinement, has features similar to QCD. One of these features is asymptotic freedom, so that in the Gross-Neveu model, as in QCD, the continuum limit is taken in the weakly coupled regime. Aoki’s scenario in the Gross-Neveu model again involves two critical lines spanning the full coupling range, with three critical points at zero coupling, linked by two cusps. This picture is based on saddle point methods . There exists substantial evidence in support of this scenario in the strongly coupled regime . In the weakly coupled regime, however, the evidence has been clear cut  and this is the region where our attention is focused. Recently, also, Creutz  has posed a question as to the size of the Aoki phase. This question is whether the Aoki phase is “squeezed out” between the arms of the cusps at non-zero coupling or whether it only vanishes in the weak coupling limit .

This sets the twofold motivation for this paper. Firstly, a new type of weak coupling expansion is developed . From it, the partition function zeroes of Wilson fermionic models can be extracted in a natural way. This weak coupling technique is then applied to the Gross-Neveu model, where the existence of an Aoki phase was first suggested . We confine our attention to the single flavour Gross-Neveu model and variants thereof. We also address the question of the “squeezing out” of the Aoki phase at weak coupling. This multiplicative approach to the single flavour Gross-Neveu model, shows that the width of the central Aoki cusp is $`𝒪(g^2)`$ while the Aoki phase has not yet emerged at this order from the left and right extremes. Furthermore, that parity symmetry breaking is the phase transition mechanism emerges in a very transparent way.

## 2 The Gross-Neveu Model

The original motivation for the introduction of the Gross-Neveu model in the continuum  was to study a renormalisable quantum field theory involving dynamical spontaneous symmetry breaking. Such models evolved from four dimensional four-fermi models studied by Nambu and Jona-Lasinio  and are essentially their two dimensional equivalents. The Gross-Neveu model is, however, renormalisable and asymptotically free. It is a model of fermions only, which interact through a short range quartic interaction. We start with a generalized Gross-Neveu model, whose action, in euclidean continuum space, is given by

$$S_{\mathrm{GN}}^{(\mathrm{cnm})}=d^2x\{\overline{\psi }(x)(/+M)\psi (x)\frac{g_\varphi ^2}{2}\left(\overline{\psi }(x)\psi (x)\right)^2\frac{g_\pi ^2}{2}\left(\overline{\psi }(x)i\gamma _S\psi (x)\right)^2\},$$
(2.1)
where $`\gamma _S=i^1\gamma _1\gamma _2`$ and $`\psi (x)`$ is a $`2`$ component fermion field. Note that we have allowed for two different four-fermion couplings. This allows for some flexibility to tune in or out the continuous chiral symmetry present in the continuum action .

We use the following representation for the Dirac $`\gamma `$-matrices in two dimensions,

$$\gamma _1^{(d=2)}=\left(\begin{array}{cc}0& 1\\ 1& 0\end{array}\right)\gamma _2^{(d=2)}=\left(\begin{array}{cc}0& i\\ i& 0\end{array}\right),$$
(2.2)
so that the chirality operator is

$$\gamma _S=i^1\gamma _1\gamma _2=\left(\begin{array}{cc}1& 0\\ 0& 1\end{array}\right).$$
(2.3)
Each term in the action (2.1) is invariant under the continuous global $`U(1)`$ symmetry

$$\psi (x)\mathrm{exp}(i\alpha )\psi (x),\overline{\psi }(x)\mathrm{exp}(i\alpha )\overline{\psi }(x).$$
(2.4)
If, further, the fermion mass $`M`$ vanishes, the action (2.1) is also invariant under a discrete global chiral $`Z_2`$ transformation

$$\psi (x)\gamma _S\psi (x),\overline{\psi }(x)\overline{\psi }(x)\gamma _S.$$
(2.5)
This is the symmetry of the original (standard) version of the model, in which the last term of (2.1) is absent (i.e., $`g_\pi =0`$). Finally, if the four fermi couplings are tuned such that $`g_\varphi =g_\pi `$, the discrete chiral symmetry is promoted to a continuous one, namely,

$$\psi (x)\mathrm{exp}(i\theta \gamma _S)\psi (x),\overline{\psi }(x)\overline{\psi }(x)\mathrm{exp}(i\theta \gamma _S).$$
(2.6)
This cannot be broken since there are no Goldstone bosons in two dimensions due to the Mermin-Wagner theorem . Nonetheless, a topological long range order of the Kosterlitz-Thouless type could exist in the model . The Mermin-Wagner theorem refers only to continuous symmetries and does not preclude the spontaneous breaking of a discrete symmetry in two dimensions. In the continuum Gross-Neveu model, the spontaneous breaking of the discrete $`\gamma _S`$ symmetry leads to dynamical fermion mass generation. The mass term explicitly breaks chiral symmetry and is analogous to an external field in the Ising model, say.

Bosonizing the action gives for the partition function,

$$Z_{GN}^{(\mathrm{cnm})}=𝒟\varphi 𝒟\pi 𝒟\overline{\psi }𝒟\psi e^S,$$
(2.7)
where

$$S=d^dx\{\overline{\psi }(x)(/+M)\psi (x)+\frac{1}{2g_\varphi ^2}\varphi ^2(x)+\frac{1}{2g_\pi ^2}\pi ^2(x)+\varphi (x)\overline{\psi }(x)\psi (x)+\pi (x)\overline{\psi }(x)i\gamma _S\psi (x)\},$$
(2.8)
where $`\varphi (x)`$ and $`\pi (x)`$ are auxiliary boson fields. The chiral transformations now represent rotations between these auxiliary fields.

## 3 Lattice Regularization with Wilson Fermions

Lattice regularization of the bosonized Gross-Neveu model, with Wilson fermions, leads to the action

$$S_F^{(\mathrm{W})}[\varphi ,\pi ,\overline{\psi },\psi ]=S_F^{(0)}[\overline{\psi },\psi ]+S_{(\mathrm{int})}[\varphi ,\pi ,\overline{\psi },\psi ]+S_{(\mathrm{bosons})}[\varphi ,\pi ],$$
(3.1)
where

$$S_F^{(0)}=\frac{1}{2\kappa }\underset{n}{}\overline{\psi }(n)\psi (n)\frac{1}{2}\underset{n,\mu }{}\left\{\overline{\psi }(n)(1\gamma _\mu )\psi (n+\widehat{\mu })+\overline{\psi }(n+\widehat{\mu })(1+\gamma _\mu )\psi (n)\right\},$$
(3.2)
$$S_{(\mathrm{int})}=\underset{n}{}\varphi (n)\overline{\psi }(n)\psi (n)+\underset{n}{}\pi (n)\overline{\psi }(n)i\gamma _S\psi (n),$$
(3.3)
and

$$S_{(\mathrm{bosons})}=\frac{1}{2g_\varphi ^2}\underset{n}{}\varphi ^2(n)+\frac{1}{2g_\pi ^2}\underset{n}{}\pi ^2(n).$$
(3.4)
The lattice sites are labeled $`n_\mu =N/2,\mathrm{},N/21`$, and $`N`$ is the number of sites in each of the two directions, which we assume to be even. Appropriate tuning of the two couplings $`g_\varphi ^2`$ and $`g_\pi ^2`$ may allow recovery of chiral symmetry in the continuum limit (see  for discussions).

Lattice Fourier transforms are defined as

$$f(n)=\left(\frac{1}{Na}\right)^2\underset{p}{}\stackrel{~}{f}(p)e^{ip.na},\stackrel{~}{f}(p)=a^2\underset{n}{}f(n)e^{ip.na},$$
(3.5)
where

$$p_\mu =\frac{2\pi }{Na}\widehat{p}_\mu ,$$
(3.6)
and where $`\widehat{p}_\mu `$ are integers or half integers depending on the field type and the boundary conditions. We henceforth drop the tilde on Fourier transformed field variables.

The fermionic part of the action can be expressed in terms of momentum space variables as

$$S_F^{(0)}[\overline{\psi },\psi ]+S_{(\mathrm{int})}[\varphi ,\pi ,\overline{\psi },\psi ]=\frac{1}{a^4}\frac{1}{N^2}\underset{q,p}{}\overline{\psi }(q)M^{(W)}(q,p)\psi (p).$$
(3.7)
Here $`M^{(W)}(q,p)`$ are $`2\times 2`$ matrices and

$$M^{(W)}(q,p)=M^{(0)}(q,p)+M^{(\mathrm{int})}(q,p),$$
(3.8)
with

$`M^{(0)}(q,p)`$ $`=`$ $`\delta _{p,q}M^{(0)}(p),`$ (3.9)
$`=`$ $`\delta _{p,q}\left[{\displaystyle \frac{1}{2\kappa }}{\displaystyle \underset{\mu }{}}\mathrm{cos}p_\mu a+i{\displaystyle \underset{\mu }{}}\gamma _\mu \mathrm{sin}p_\mu a\right],`$ (3.10)
and

$$M^{(\mathrm{int})}(q,p)=\frac{1}{N^2}\underset{n}{}e^{i(pq)na}\left[\varphi (n)+\pi (n)i\gamma _S\right].$$
(3.11)
Integration over the Grassmann variables gives the full partition function

$$Z=𝒟\varphi 𝒟\pi 𝒟\overline{\psi }𝒟\psi \mathrm{exp}\left(S_F^{(\mathrm{W})}\right)detM^{(\mathrm{W})}\underset{\alpha ,p}{}\lambda _\alpha (p),$$
(3.12)
with $`\lambda _\alpha (p)`$ the eigenvalues of the fermion matrix and the expectation values being taken over the bosonic fields. Note that there is no hopping parameter dependence in $`M^{(\mathrm{int})}`$.

In the free fermion case the partition function is simply proportional to

$$detM^{(0)}=\underset{\alpha ,p}{}\lambda _\alpha ^{(0)}(p),$$
(3.13)
where $`\lambda _\alpha ^{(0)}(p)`$ are the eigenvalue solutions of

$$M^{(0)}(p)|\lambda ^{(0)}=\lambda ^{(0)}|\lambda ^{(0)}.$$
(3.14)
Using the representation (2.2) for the Dirac $`\gamma `$-matrices, the solution to this problem is easily found to be

$`|\lambda _\alpha ^{(0)}(p)`$ $`=`$ $`{\displaystyle \frac{1}{\sqrt{2}}}\left(\begin{array}{c}1\\ (1)^\alpha \frac{\mathrm{sin}p_1a+i\mathrm{sin}p_2a}{\sqrt{_{\mu =1}^2\mathrm{sin}^2p_\mu a}}\end{array}\right),`$ (3.17)
$`\lambda _\alpha ^{(0)}(p)`$ $`=`$ $`{\displaystyle \frac{1}{2\kappa }}{\displaystyle \underset{\mu =1}{\overset{2}{}}}\mathrm{cos}p_\mu a+i(1)^\alpha \sqrt{{\displaystyle \underset{\mu =1}{\overset{2}{}}}\mathrm{sin}^2p_\mu a},`$ (3.18)
where $`\alpha =1,2`$. These eigenfunctions form a complete orthonormal set. As is usual for Grassmann variables, we impose antiperiodic boundary conditions in the temporal ($`1`$-) direction and periodic boundary conditions in the spatial ($`2`$-) direction. With these mixed boundary conditions the momenta $`\widehat{p}_\mu `$ for the Fourier transformed fermion fields take the integer or half-integer values, $`\widehat{p}_1=N/2+1/2,\mathrm{},N/21/2`$ and $`\widehat{p}_2=N/2,\mathrm{},N/21`$. Then, the eigenvalues (3.18) in the free fermion case are either two-fold or four-fold degenerate, the former being the case if $`\widehat{p}_2=0`$ or $`N/2`$.

In the free fermion case the Lee-Yang Zeroes  are given by $`\lambda _\alpha ^{(0)}(p)=0`$. From (3.18), this is the case at

$$\frac{1}{2\kappa }=\eta _\alpha (p)=\underset{\mu =1}{\overset{2}{}}\mathrm{cos}p_\mu ai(1)^\alpha \sqrt{\underset{\mu =1}{\overset{2}{}}\mathrm{sin}p_\mu a}.$$
(3.19)
The lowest zeroes (with the smallest imaginary parts) correspond to

$$\widehat{p}=(\pm (N/21/2),N/2),(\pm 1/2,N/2),(\pm (N/21/2),0)\mathrm{and}(\pm 1/2,0),$$
(3.20)
impacting onto the real axis at $`1/2\kappa =2,0,0`$ and $`2`$ respectively. These are precisely the three nadirs of the Aoki cusps in the Gross-Neveu model (see Fig. 1).

Note that the zeroes in the upper half plane are given by $`\alpha =1`$, while their complex conjugates correspond to $`\alpha =2`$. Note, further, that the zeroes in (3.19) are two or four fold degenerate in the momenta. I.e., these zeroes are invariant under $`p_\mu p_\mu `$. This transformation is just a rotation through an angle $`\pi `$ in the space-time plane. The lowest zeroes (3.20), which are responsible for the critical behaviour of the free model, are actually two-fold degenerate. There is also a symmetry under $`p_1p_2`$ which is manifest in the infinite volume limit. This is equivalent to a trivial rotation by $`\pi /2`$ in the $`p_1`$-$`p_2`$ plane, followed by reflection through the $`p_2`$ axis. Since this reflection is through the spatial axis, this transformation is, in fact, parity. I.e., apart from rotations in space and time, the critical points are left unchanged under the parity transformation.

## 4 A New Weak Coupling Expansion

The usual weak coupling expansion of the full determinant for a general fermionic field theory is the Taylor expansion of

$`detM^{(\mathrm{W})}`$ $`=`$ $`detM^{(0)}\times det\left(M_{}^{(0)}{}_{}{}^{1}M^{(\mathrm{W})}\right)`$
$`=`$ $`detM^{(0)}\mathrm{exp}\mathrm{tr}\mathrm{ln}\left(1+M_{}^{(0)}{}_{}{}^{1}M^{(\mathrm{int})}\right).`$ (4.1)
This expansion is additive in nature and, from it, the ratio of full to free fermion determinants may be written,

$$\frac{detM^{(\mathrm{W})}}{detM^{(0)}}=1+\underset{i=1}{\overset{2N^2}{}}\frac{M_{ii}^{(\mathrm{int})}}{\lambda _i^{(0)}}\frac{1}{2}\underset{i,j=1}{\overset{2N^2}{}}\frac{M_{ij}^{(\mathrm{int})}M_{ji}^{(\mathrm{int})}}{\lambda _i^{(0)}\lambda _j^{(0)}}+\frac{1}{2}\underset{i,j=1}{\overset{2N^2}{}}\frac{M_{ii}^{(\mathrm{int})}M_{jj}^{(\mathrm{int})}}{\lambda _i^{(0)}\lambda _j^{(0)}}+\mathrm{}.$$
(4.2)
Here the indices $`i`$ and $`j`$ stand for the combination of Dirac index and momenta $`(\alpha ,p)`$ which label fermionic matrix elements, so that $`M_{ij}^{(\mathrm{int})}M_{(\alpha p)(\beta q)}^{(\mathrm{int})}`$ represents $`\lambda _\alpha ^{(0)}(p)|M^{(\mathrm{int})}(p,q)|\lambda _\beta ^{(0)}(q)`$. The traces in (4.2) are, in fact, the diagrams which contribute to the vacuum polarization tensor.

Setting

$`t_i`$ $`=`$ $`M_{ii}^{(\mathrm{int})},`$ (4.3)
$`s_{ij}`$ $`=`$ $`s_{ji}=M_{ii}^{(\mathrm{int})}M_{jj}^{(\mathrm{int})},`$ (4.4)
$`t_{ij}`$ $`=`$ $`t_{ji}=M_{ij}^{(\mathrm{int})}M_{ji}^{(\mathrm{int})}s_{ij},`$ (4.5)
the ratio of the interactive and free partition functions may be written

$$\frac{detM^{(\mathrm{W})}}{detM^{(0)}}=1+\underset{i=1}{\overset{2N^2}{}}\frac{t_i}{\lambda _i^{(0)}}\frac{1}{2}\underset{i,j=1}{\overset{2N^2}{}}\frac{t_{ij}}{\lambda _i^{(0)}\lambda _j^{(0)}}+\mathrm{}.$$
(4.6)
This Taylor expansion is analytic in $`1/2\kappa `$ with poles at $`\lambda _i^{(0)}=0`$ or $`1/2\kappa =\eta _i^{(0)}`$.

The Wilson fermion matrix $`M^{(\mathrm{W})}`$ is a $`2N^2`$ dimensional square matrix given by (3.8)-(3.11). Its determinant, and the bosonic expectation value thereof, are therefore polynomials of degree $`2N^2`$ with corresponding number of zeroes. As such, the latter may be written (up to an irrelevant constant)

$$detM^{(\mathrm{W})}=\underset{i=1}{\overset{2N^2}{}}\left(1/2\kappa \eta _i\right),$$
(4.7)
where $`\eta _i`$ represents $`\eta _\alpha (p)`$ and are the Lee-Yang zeroes of the interactive model. These are the quantities to be determined at weak coupling.

Writing

$$\mathrm{\Delta }_i=\eta _i\eta _{i}^{}{}_{}{}^{(0)},$$
(4.8)
gives, now, a new type of weak coupling expansion for the ratio of partitions functions, which is ‘multiplicative’ rather than additive in form,

$$\frac{detM^{(\mathrm{W})}}{detM^{(0)}}=\underset{i=1}{\overset{2N^2}{}}\left(\frac{1/2\kappa \eta _i}{\lambda _i^{(0)}}\right)=\underset{i=1}{\overset{2N^2}{}}\left(1\frac{\mathrm{\Delta }_i}{1/2\kappa \eta _i^{(0)}}\right).$$
(4.9)
Note that the expression (4.9), like its additive counterpart, (4.6), analytic in $`1/2\kappa `$ with poles at $`\eta _i^{(0)}`$.

Expanding (4.9) gives

$$\frac{detM^{(\mathrm{W})}}{detM^{(0)}}=1\underset{i=1}{\overset{2N^2}{}}\frac{\mathrm{\Delta }_i}{\lambda _i^{(0)}}+\frac{1}{2}\underset{i=1}{\overset{2N^2}{}}\underset{ji}{\overset{2N^2}{}}\frac{\mathrm{\Delta }_i\mathrm{\Delta }_j}{\lambda _i^{(0)}\lambda _j^{(0)}}+\mathrm{}.$$
(4.10)
Let $`\{n\}`$ denote the $`n^{\mathrm{th}}`$ degeneracy class in the free fermion case, so that the $`D_n`$ eigenvalues $`\lambda _{n_1}^{(0)}=\mathrm{}=\lambda _{n_{D_n}}^{(0)}`$ are identical to $`\lambda _n^{(0)}`$, say, with $`D_n=2`$ or $`4`$. Take the hopping parameter to be complex and arbitrarily close to a free fermion zero,

$$1/2\kappa =\eta _n^{(0)}+ϵ.$$
(4.11)
The additive and multiplicative expressions (4.6) and (4.10) for the ratio of partition functions may now be expanded in $`ϵ^1`$. Indeed, (4.6) gives

$$\frac{detM^{(\mathrm{W})}}{detM^{(0)}}=ϵ^2\frac{1}{2}\underset{n_i,n_j\{n\}}{}t_{n_in_j}+ϵ^1\left\{\underset{n_i\{n\}}{}t_{n_i}\underset{n_i\{n\},j\{n\}}{}\frac{t_{n_ij}}{\eta _n^{(0)}\eta _j^{(0)}+ϵ}\right\}+𝒪\left(ϵ^0\right),$$
(4.12)
while (4.10) yields

$`{\displaystyle \frac{detM^{(\mathrm{W})}}{detM^{(0)}}}`$ $`=`$ $`ϵ^2{\displaystyle \underset{n_i,n_j\{n\},n_in_j}{}}{\displaystyle \frac{\mathrm{\Delta }_{n_i}\mathrm{\Delta }_{n_j}}{2}}`$ (4.13)
$`+ϵ^1\left\{{\displaystyle \underset{n_i\{n\}}{}}\mathrm{\Delta }_{n_i}+{\displaystyle \underset{n_i\{n\},j\{n\}}{}}{\displaystyle \frac{\mathrm{\Delta }_{n_i}\mathrm{\Delta }_j}{\eta _n^{(0)}\eta _j^{(0)}+ϵ}}\right\}+𝒪\left(ϵ^0\right).`$
Equating these two expansions to $`𝒪(ϵ^2)`$ gives

$$\underset{n_i,n_j\{n\},n_in_j}{}\mathrm{\Delta }_{n_i}\mathrm{\Delta }_{n_j}=\underset{n_i,n_j\{n\}}{}t_{n_in_j},$$
(4.14)
while to $`𝒪(ϵ^1)`$ it gives

$$\underset{n_i\{n\}}{}\mathrm{\Delta }_{n_i}\left\{1\underset{j\{n\}}{}\frac{\mathrm{\Delta }_j}{\eta _n^{(0)}\eta _j^{(0)}}\right\}=\underset{n_i\{n\},j\{n\}}{}\frac{t_{n_ij}}{\eta _n^{(0)}\eta _j^{(0)}}\underset{n_i\{n\}}{}t_{n_i},$$
(4.15)
having taken the $`ϵ0`$ limit.

Let

$`t_i`$ $`=`$ $`t_i^{(1)}+t_i^{(2)}+𝒪(g^3),`$ (4.16)
$`t_{ij}`$ $`=`$ $`t_{ij}^{(2)}+𝒪(g^3),`$ (4.17)
$`\mathrm{\Delta }_i`$ $`=`$ $`\eta _i^{(1)}+\eta _i^{(2)}+𝒪(g^3),`$ (4.18)
where $`t_i^{(1)}`$ and $`\eta _i^{(1)}`$ are the order $`g`$ contributions to the expectation values of the matrix elements and zero shifts and where $`t_i^{(2)}`$, $`t_{ij}^{(2)}`$ and $`\eta _i^{(2)}`$ are their order $`g^2`$ equivalents. The $`𝒪(ϵ^1)`$ equation to order $`g`$ is

$$\underset{n_i\{n\}}{}\eta _{n_i}^{(1)}=\underset{n_i\{n\}}{}t_{n_i}^{(1)},$$
(4.19)
and its order $`g^2`$ counterpart is

$$\underset{n_i\{n\}}{}\eta _{n_i}^{(2)}=\underset{n_i\{n\}}{}t_{n_i}^{(2)}+\underset{n_i\{n\},j\{n\}}{}\frac{t_{n_ij}^{(2)}+t_{n_i}^{(1)}t_j^{(1)}}{\eta _n^{(0)}\eta _j^{(0)}}.$$
(4.20)
Also, the $`𝒪(ϵ^2)`$ equation, (4.14), is, now,

$$\underset{n_i\{n\}}{}\left(\eta _{n_i}^{}{}_{}{}^{(1)}\right)^2=\underset{n_i,n_j\{n\}}{}t_{n_in_j}^{(2)}+\left(\underset{n_i\{n\}}{}t_{n_i}^{(1)}\right)^2.$$
(4.21)
With relations (4.19)-(4.21), the multiplicative expression (4.10) recovers (4.6) to $`𝒪(g^2)`$. Thus, equating the $`𝒪\left(ϵ^0\right)`$ contributions to (4.12) and (4.13) yields no extra information.

The partition function zeroes are ‘protocritical points’ in the sense that they have the potential to become true critical points . In the limit of infinite volume, the lowest zeroes impact on to the real hopping parameter axis precipitating the phase transition. The real parts of the lowest zeroes are therefore pseudocritical points in the statistical mechanics sense.

In the free case, the lowest zeroes, and those responsible for criticality, are two fold degenerate. One expects critical behaviour in the weakly coupled case to be governed by their equivalents there. The two equations, (4.19) and (4.21), allow full determination of the first order shifts to two-fold degenerate zeroes. Indeed,

$$\eta _{n_i}^{(1)}=\frac{1}{2}\left\{t_{n_1}^{(1)}t_{n_2}^{(1)}\pm \sqrt{(t_{n_1}^{(1)}+t_{n_2}^{(1)})^2+4t_{n_1n_2}^{(2)}}\right\},$$
(4.22)
where $`n_i\{n\}`$ for $`i=1`$ or $`2`$. The second order equation, (4.20), in the two-fold degenerate case is

$$\eta _{n_1}^{(2)}+\eta _{n_2}^{(2)}=t_{n_1}^{(2)}t_{n_2}^{(2)}+\underset{j\{n\}}{}\frac{(t_{n_1}^{(1)}+t_{n_2}^{(1)})t_j^{(1)}+t_{jn_1}^{(2)}+t_{jn_2}^{(2)}}{\eta _n^{(0)}\eta _j^{(0)}}.$$
(4.23)
To find the individual shifts, let

$`\eta _{n_1}^{(2)}`$ $`=`$ $`\eta _n^{(2)}+\delta ^{(2)},`$ (4.24)
$`\eta _{n_2}^{(2)}`$ $`=`$ $`\eta _n^{(2)}\delta ^{(2)}.`$ (4.25)
Their average, $`\eta _n^{(2)}`$, is determined directly from (4.23). Removing the expectation values over the bosonic fields converts the zeroes to the shifts in the eigenvalues of the fermion matrix in the presence of a small perturbation, $`M^{(\mathrm{int})}`$. The problem of determining such shifts is simply (two-fold degenerate) time independent perturbation theory. Indeed, one finds, for example,

$`\lambda _{n_1}`$ $`=`$ $`\lambda _n^{(0)}{\displaystyle \frac{1}{2}}\left\{M_{n_1n_1}^{(\mathrm{int})}+M_{n_2n_2}^{(\mathrm{int})}\pm \sqrt{(M_{n_1n_1}^{(\mathrm{int})}M_{n_1n_1}^{(\mathrm{int})})^2+4M_{n_1n_2}^{(\mathrm{int})}M_{n_2n_1}^{(\mathrm{int})}}\right\}`$ (4.26)
$`+{\displaystyle \frac{1}{2}}{\displaystyle \underset{j\{n\}}{}}{\displaystyle \frac{M_{jn_1}^{(\mathrm{int})}M_{n_1j}^{(\mathrm{int})}+M_{jn_2}^{(\mathrm{int})}M_{n_2j}^{(\mathrm{int})}}{\lambda _n^{(0)}\lambda _j^{(0)}}}\delta _0^{(2)},`$
in which $`\delta _0^{(2)}=\delta ^{(2)}`$. This recovers time independent perturbation theory if

$$\delta _0^{(2)}=\frac{1}{2}\underset{j\{n\}}{}\frac{M_{jn_1}^{(\mathrm{int})}M_{n_1j}^{(\mathrm{int})}M_{jn_2}^{(\mathrm{int})}M_{n_2j}^{(\mathrm{int})}}{\lambda _n^{(0)}\lambda _j^{(0)}},$$
(4.27)
whence

$$\delta ^{(2)}=\frac{1}{2}\underset{j\{n\}}{}\frac{t_{jn_1}^{(2)}+s_{jn_1}^{(2)}t_{jn_2}^{(2)}s_{jn_2}^{(2)}}{\eta _n^{(0)}\eta _j^{(0)}}.$$
(4.28)
Finally, the full expression for the second order shift in an erstwhile two-fold degenerate zero is

$$\eta _{n_i}^{(2)}=\frac{1}{2}\left(t_{n_1}^{(2)}+t_{n_2}^{(2)}\right)+\frac{1}{2}\underset{j\{n\}}{}\frac{2(t_{jn_i}^{(2)}+s_{jn_i}^{(2)})+t_j^{(1)}t_{n_1}^{(1)}+t_j^{(1)}t_{n_2}^{(1)}s_{jn_1}^{(2)}s_{jn_2}^{(2)}}{\eta _n^{(0)}\eta _j^{(0)}}.$$
(4.29)
## 5 The Zeroes and Phase Diagram of the Gross-Neveu Model

The interactive part of the fermion matrix (3.11) may be split into

$$M^{(\mathrm{int})}(q,p)=M_\varphi ^{(\mathrm{int})}(q,p)+M_\pi ^{(\mathrm{int})}(q,p),$$
(5.1)
where

$`M_\varphi ^{(\mathrm{int})}(q,p)`$ $`=`$ $`{\displaystyle \frac{1}{N^2}}{\displaystyle \underset{n}{}}e^{i(pq)na}\varphi (n)=\left({\displaystyle \frac{1}{Na}}\right)^2\varphi (qp),`$ (5.2)
$`M_\pi ^{(\mathrm{int})}(q,p)`$ $`=`$ $`{\displaystyle \frac{1}{N^2}}{\displaystyle \underset{n}{}}e^{i(pq)na}\pi (n)i\gamma _S=\left({\displaystyle \frac{1}{Na}}\right)^2\pi (qp)i\gamma _5.`$ (5.3)
One notes that the momentum dependency of the bosonic field variables involves even integers, so the bosons have periodic boundary conditions. The generic matrix elements required for the calculation of (4.3), (4.4) and (4.5) are

$`M_{\varphi }^{(\mathrm{int})}{}_{(\alpha p)(\beta q)}{}^{}`$ $`=`$ $`\left({\displaystyle \frac{1}{Na}}\right)^2\varphi (pq){\displaystyle \frac{1}{2}}\times `$ (5.4)
$`\left(1+(1)^{\alpha +\beta }{\displaystyle \frac{_\mu \mathrm{sin}p_\mu a\mathrm{sin}q_\mu a+i(\mathrm{sin}p_1a\mathrm{sin}q_2a\mathrm{sin}p_2a\mathrm{sin}q_1a)}{\sqrt{_\mu \mathrm{sin}^2p_\mu a}\sqrt{_\mu \mathrm{sin}^2p_\mu a}}}\right),`$
$`M_{\pi }^{(\mathrm{int})}{}_{(\alpha p)(\beta q)}{}^{}`$ $`=`$ $`\left({\displaystyle \frac{1}{Na}}\right)^2\pi (pq){\displaystyle \frac{i}{2}}\times `$ (5.5)
$`\left(1(1)^{\alpha +\beta }{\displaystyle \frac{_\mu \mathrm{sin}p_\mu a\mathrm{sin}q_\mu a+i(\mathrm{sin}p_1a\mathrm{sin}q_2a\mathrm{sin}p_2a\mathrm{sin}q_1a)}{\sqrt{_\mu \mathrm{sin}^2p_\mu a}\sqrt{_\mu \mathrm{sin}^2p_\mu a}}}\right).`$
In the (generalized) Gross-Neveu case, the pure bosonic action is given by (3.4). The pure bosonic expectation values in momentum space are thus

$`\varphi (k)`$ $`=`$ $`\pi (k)=0,`$ (5.6)
$`\varphi (k)\varphi (k)`$ $`=`$ $`N^da^{2d}2g_\varphi ^2,`$ (5.7)
$`\pi (k)\pi (k)`$ $`=`$ $`N^da^{2d}2g_\pi ^2.`$ (5.8)
The bosonic expectation values of the matrix elements required in the calculation of the shifts (4.22) and (4.29) are then

$`t_i`$ $``$ $`t_{\alpha ,p}=0,`$ (5.9)
$`s_{ij}`$ $``$ $`s_{(\alpha ,p)(\beta ,q)}={\displaystyle \frac{2g_{\varphi }^{}{}_{}{}^{2}}{N^2}},`$ (5.10)
$`t_{ij}`$ $``$ $`t_{(\alpha ,p)(\beta ,q)}={\displaystyle \frac{g_{\varphi }^{}{}_{}{}^{2}+g_{\pi }^{}{}_{}{}^{2}}{N^2}}\left\{(1)^{\alpha +\beta }{\displaystyle \frac{_\rho \mathrm{sin}p_\rho \mathrm{sin}q_\rho }{\sqrt{_\mu \mathrm{sin}^2p_\mu _\nu \mathrm{sin}^2q_\nu }}}1\right\}.`$ (5.11)
From these equations, together with (4.22) and (4.29), the $`𝒪(g)`$ and $`𝒪(g^2)`$ shifts for the erstwhile two-fold degenerate zeroes, $`\eta _\alpha (\pm |p_1|,p_2)`$ (for $`\widehat{p}_2=0`$ or $`N/2`$), are, respectively,

$`\eta _\alpha ^{(1)}(\pm |p_1|,p_2)`$ $`=`$ $`\pm i{\displaystyle \frac{\sqrt{2}\sqrt{g_\varphi ^2+g_\pi ^2}}{N}},`$ (5.12)
$`\eta _\alpha ^{(2)}(\pm |p_1|,p_2)`$ $`=`$ $`{\displaystyle \frac{g_\varphi ^2+g_\pi ^2}{N^2}}{\displaystyle \underset{(\beta ,q)\{(\alpha ,p)\}}{}}{\displaystyle \frac{1}{\eta _\alpha ^{(0)}(p)\eta _\beta ^{(0)}(q)}}.`$ (5.13)
So the two four-fermi interactions in fact contribute the same amounts to the shifts in the zeroes.

In the thermodynamic limit these lowest zeroes become the true critical points of the theory and their determination amounts to determination of the phase diagram. I.e., the phase diagram is given to order $`g^2`$ by the limit

$$\frac{1}{2\kappa (g)}=\underset{N\mathrm{}}{lim}\left\{\eta _\alpha ^{(0)}(p)+\eta _\alpha ^{(1)}(p)+\eta _\alpha ^{(2)}(p)\right\},$$
(5.14)
where $`p`$ is the momentum corresponding to the lowest zeroes. The first order shift in (5.12) gives the relative separation of the erstwhile two-fold degenerate zeroes and vanishes in the infinite volume limit. The average shift is represented by (5.13) and is second order. The shift in the corresponding critical point is

$$\underset{N\mathrm{}}{lim}\eta _\alpha ^{(2)}(\pm |p_1|,p_2)=(g_\varphi ^2+g_\pi ^2)\underset{N\mathrm{}}{lim}c(N),$$
(5.15)
where

$$c(N)=\frac{1}{N^2}\underset{(\beta ,q)\{(\alpha ,p)\}}{}\frac{1}{\eta _\alpha ^{(0)}(p)\eta _\beta ^{(0)}(q)}.$$
(5.16)
One finds, numerically, that the imaginary contribution to this factor vanishes in the thermodynamic limit, meaning that these zeroes indeed impact on to the real hopping parameter axis. The real part of (5.16) becomes an $`N`$-independent constant whose actual value depends on the free zero from which it evolved. Indeed, (5.16) approaches approximately $`0.77`$ and $`0.77`$ for $`(\widehat{p}_1,\widehat{p}_2)=(\pm 1/2,0)`$ and $`(\pm (N1)/2,N/2)`$ respectively. These correspond to the rightmost and leftmost critical lines (see Fig. 1). Also, (5.16) is approximately $`0.2`$ and $`0.2`$ for $`(\widehat{p}_1,\widehat{p}_2)=(\pm (N1)/2,0)`$ and $`(\pm 1/2,N/2)`$ respectively. These give the two lines that generate the inner cusp. The situation is summarized in Table 1 where the critical hopping parameters in the free and interacting cases and the momentum indices of the corresponding zeroes are listed.

The actual phase diagram for weakly coupling is pictured in Fig. 1 for $`g_\varphi =g_\pi =g`$ (dark lines). The lighter curves are a schematic representation of the expected full phase diagram. One sees the degeneracy of the central free fermion critical point is lifted and two critical lines emerge in the presence of weak bosonic coupling. These are the lines corresponding to the central cusp in Aoki’s phase diagram. From (5.15), the central cusp can only be made vanish at this order, in the unphysical situation of imaginary couplings. The Aoki phase does not yet emerge to $`𝒪(g^2)`$ from the left- and rightmost free critical points.

In the free case, the zeroes and hence the critical points are invariant under momentum inversion $`p_\mu p_\mu `$, corresponding to a rotation in space-time. While this degeneracy is lifted at finite size in the interacting case, it is recovered in the limit of infinite volume. In that limit, the free zeroes and critical points are also invariant under the parity transformation $`p_1p_2`$. This is no longer the case in the presence of interactions. Indeed, the inner pair of critical lines are interchanged under parity. The overall phase structure, however, remains the same.

The situation is similar to the two dimensional Potts model. There, the partition function is invariant under a duality transformation which exchanges the high and low temperature phases. The critical point is that which is invariant under that transformation. Here, the zeroes, and hence the partition function, are invariant under parity. The phase structure is also unchanged by parity. However, parity even and parity odd regions of the phase diagram are interchanged.

## 6 Conclusions

A new type of weak coupling expansion appropriate for Wilson fermionic lattice field theories has been developed. This expansion is multiplicative, but recovers the standard additive expansion. Its multiplicative form allows the Lee-Yang zeroes of the weakly coupled theory to be extracted in a natural way. These zeroes, are protocritical points, which, if they impact on to the real hopping parameter axis, precipitate a phase transition there.

The expansion is applied to the single flavour lattice Gross-Neveu model to track the movement of zeroes and thereby the critical points in the presence of bosonic field variables. This model shares features with QCD, one of which is expected to be the existence of an Aoki phase.

Using the new weak coupling expansion, a phase diagram is obtained in the weakly coupled region which is consistent with that of Aoki. The widths of the Aoki cusps are analytically determined to second order in the couplings. The central cusp cannot be tuned away for real physical couplings. The lateral cusps do not yet emerge at this order. This is the answer to the question posed by Creutz in  for the single flavour Gross-Neveu model.

Finally, while the full phase structure is unaltered by a parity transformation, such an operation has the effect of exchanging the critical lines forming the inner Aoki cusp.


format: text

text:
# Mechanical Feedback: From Stellar Wind Bubbles to Starbursts

## 1 Introduction

Mechanical feedback from massive stars is one of the principal drivers of evolutionary processes in galaxies. Our understanding of the feedback process is grounded in the conventional model for the evolution of stellar wind- and supernova-driven superbubbles (e.g., Pikel’ner 1968; Weaver et al. 1977): hot (log $`T/\mathrm{K}67`$), low-density ($`n0.01\mathrm{cm}^3`$) gas is generated within a double-shock structure, and the pressure of this hot gas, chemically enriched by the stellar products, drives the growth of the thin, radiatively-cooled shell. In the standard model, energy loss from the hot gas is considered to be negligible, yielding simple analytic expressions for the shell radius $`R`$ and expansion velocity $`v`$ as a function of time $`t`$:

$`R`$ $``$ $`(L/n)^{1/5}t^{3/5},`$
$`v=\text{ }\mathrm{}R`$ $``$ $`(L/n)^{1/5}t^{2/5}.`$ (1)
For a stellar wind-driven bubble, the input mechanical luminosity $`L=1/2\text{ }\mathrm{}Mv_{\mathrm{}}^2`$, where $`\text{ }\mathrm{}M`$ and $`v_{\mathrm{}}`$ are the stellar wind mass-loss rate and terminal velocity, respectively. For OB associations, supernovae quickly dominate over winds, in which case $`L=N_{}E_{51}/t_e`$ (e.g., Mac Low & McCray 1988), where $`N_{}`$ and $`E_{51}`$ are the total number of supernovae and supernova energy, respectively, and $`t_e`$ is the time over which the supernovae occur.

The fate of this hot gas is crucial in understanding the phase balance and enrichment of the interstellar and intergalactic media. Is stellar feedback indeed the source of the diffuse, hot gas in the ISM? Do galactic superwinds from starbursts eject gas from galaxies? How does mechanical feedback affect the structure of the ISM and porosity, for example, to ionizing radiation? Do superbubbles trigger renewed star formation? To answer these basic questions, we must first evaluate the relevance of the standard, adiabatic model for shell evolution, which is applied over scales ranging from single stellar wind bubbles to galactic superwinds.

## 2 Single star bubbles

Single star bubbles of isolated OB stars are woefully understudied. Oey & Massey (1994) examined two nebular examples in M33, and spectroscopically classified the late-type O stars. The inferred stellar masses and ages implied wind parameters that were consistent with the observed sizes and shell ages predicted by the adiabatic model. However, the parameters were loosely constrained and the shell kinematics were not observed. H i shells with radii of several tens of pc have been identified as wind-blown bubbles around a number of Galactic O and Of stars (Cappa & Benaglia 1998; Benaglia & Cappa 1999). These are largely consistent with the standard model, and probe a specific subset of fairly evolved stars with old shells that have essentially stopped expanding.

Studies of Wolf-Rayet ring nebulae suggest shells that are too small, equivalent to an overestimate in $`L/n`$ by an order of magnitude (e.g., Treffers & Chu 1982; García-Segura & Mac Low 1995; Drissen et al. 1995). However, the progenitor star produces several wind phases, including both fast and slow winds, with strongly changing $`\text{ }\mathrm{}M`$. Their cumulative effect on shell morphology is even more complex and poorly-understood than the wind evolution itself, and it is therefore unsurprising to find significant discrepancies between the predictions and observations of shell parameters. Hence it is desireable to study more of the simpler, single OB star bubbles.

## 3 Superbubbles

Superbubbles created by OB associations are more prominent than single-star bubbles, and thus have been studied more extensively. Soft X-ray emission has been detected from the interiors of many objects in projection, which is qualitatively consistent with the adiabatic superbubble model. Quantitatively, two classes of X-ray emission have been identified: objects with X-ray luminosity $`L_\mathrm{x}`$ in excess of the standard model’s prediction (Chu & Mac Low 1990; Wang & Helfand 1991), and objects that remain undetected in X-rays (Chu et al. 1995). The X-ray–bright objects are suggested to be overluminous because of SNR impacts on the shell walls. Upper limits on the X-ray–dim objects remain consistent with $`L_\mathrm{x}`$ predicted by the adiabatic model. It will thus be of great interest to determine $`L_\mathrm{x}`$ for these objects with XMM-Newton or Chandra. The existence of hot gas within superbubbles also implies an interface region with the cooler shells where intermediate temperatures and ions should be present. Chu et al. (1994) searched a number sightlines through LMC superbubbles and confirmed the existence of C iv and Si iv absorption in all cases.

A stringent test of the adiabatic model is to compare the predicted and observed shell kinematics where the input mechanical power and other parameters are well-constrained. We obtained spectroscopic classifications of the OB associations in eight LMC superbubbles (Oey & Massey 1995; Oey 1996; Oey & Smedley 1998), thereby providing estimates of $`L`$ in these young, wind-dominated objects. The predicted growth rate for the shells was higher than implied by their observed $`R`$ and $`v`$, equivalent to an overestimate in $`L/n`$ by an order of magnitude. However, even after adjusting $`L/n`$ in the models, over half the objects still showed observed expansion velocities that were typically a factor of two higher than predicted for the given $`R`$. Similar discrepancies were reported for Galactic objects by Saken et al. (1992) and Brown et al. (1995). Figure 1 shows representative examples: the two objects in the left column show predicted $`R`$ (solid curves) and $`v`$ (dotted curves) in agreement with the observed values (horizontal segments); while the objects on the right show anomalously high observed $`v`$. All the models have $`L/n`$ reduced by a factor of 10. The objects with anomalously high $`v`$ also each exhibit anomalously high $`L_\mathrm{x}`$ and \[S ii\]/H$`\alpha `$, consistent with the suggestion by Chu & Mac Low (1990) and Wang & Helfand (1991) that SNR impacts have heated the shell. These would then also accelerate the shell walls, explaining the high observed $`v`$. The observed stellar mass function in these objects also support this scenario, implying the former existence of 1–2 higher-mass stars corresponding to SN progenitors.

While the SNR impact hypothesis appears the most likely explanation for the high-velocity shells, a sudden drop in the ambient density can induce a “mini-blowout” with shell kinematics that can easily reproduce the anomalous observations (Oey & Smedley 1998; Mac Low et al. 1998; Franco & Silich 1999). Indeed, were it not for the high X-ray emission and additional evidence above, it would be impossible to distinguish between the two shell acceleration mechanisms from the kinematics alone.

We therefore see that the ambient density distribution is critical in determining the shell evolution. An underestimate in $`n`$ could, for example, contribute to the growth rate discrepancy described above, that is seen in all the objects. To clarify the effect of the ambient gas distribution, Oey et al. (2001) mapped the H i distribution within a $`40^{}`$ radius of three nebular LMC superbubbles at 30<sup>′′</sup> resolution, using the Australia Telescope Compact Array. The results show neutral environments that vary to an extreme: DEM L25 shows no H i shell component, but appears to be nestled in an H i hole; DEM L50 shows an H i shell with $`10`$ times more mass than the H$`\alpha `$ component, but is otherwise in a large (1.4 kpc) H i void; and DEM L301 is in a non-descript neutral environment with no apparent correspondence to the optical shell. Thus, despite the similar appearance of the H$`\alpha `$ nebulae, the H i distributions could not be more different. It is therefore extremely difficult to infer properties of ambient density distribution without direct, empirical determinations.

Another vital parameter for shell evolution is the ambient pressure, $`P_0`$, which determines whether and when the superbubble growth becomes pressure-confined. While $`P_0`$ is usually relatively unimportant in young, high-pressure superbubbles such as the nebular objects modeled above, it is of vital importance in the mid- to late-stage evolution. It may also be relevant in high-pressure, ionized environments like dense star-forming regions (e.g., García-Segura & Franco 1996). The value of $`P_0`$ ultimately determines the final size of the objects, and conditions relative to blowout. The uniformity and distribution of $`P_0`$ in the multiphase ISM is therefore especially relevant to the global effect of superbubbles on the ISM (see below).

Finally, if the hot gas within superbubbles does not somehow blow out and merge into the hot, ionized medium (HIM), it is likely that the objects will cool and depart from energy conservation. Indeed, cooling of the hot interior has long been one of the principal questions for superbubble evolution and the fate of the hot gas. Thermal conduction at the interface between the cool shell wall and hot gas should cause a high rate of mass-loading into the interior, by evaporation. The evaporated shell material strongly dominates the mass in the hot region, which could be further supplemented by evaporation and ablation from small clouds overrun by the expanding shocks (e.g., Cowie & McKee 1977; McKee et al. 1984; Arthur & Henney 1996). If the interior density is sufficiently increased, radiative cooling will dominate, and the shells will no longer grow adiabatically. In addition, Silich et al. (2001) point out the importance of enhanced metallicity in the superbubble interiors, caused by the stellar and SN yields. Preliminary investigation for individual objects by Oey & Silich (2001, in preparation) shows enhancement in $`L_\mathrm{x}`$ by almost an order of magnitude for low-metallicity ($`Z=0.05Z_{}`$) objects. The elemental yields therefore may also significantly enhance the cooling rate, facilitating the transition from adiabatic to momentum-conserving evolution.

## 4 Global mechanical feedback and superwinds

Quantifying the statistical properties of superbubbles in galaxies allows us to test the evolutionary model and also gain insight on the properties of the ISM. In particular, the so-called interstellar porosity (see below) is a conventional parameter for evaluating the relative importance of the HIM in the multiphase ISM.

### 4.1 Superbubble populations

Oey & Clarke (1997) derived expressions for the differential size distribution $`N(R)dR`$ of superbubbles in a uniform ISM, using the analytic expressions for adiabatic evolution (equation 1). We considered a power-law mechanical luminosity function for the parent OB associations,

$$\varphi (L)dL=AL^\beta dL,$$
(2)
with $`\beta =2`$, based on the typical H ii region luminosity function (H ii LF; e.g., Kennicutt et al. 1989; Oey & Clarke 1998a). The superbubble growth was also taken to be pressure-confined when the interior pressure $`P_i=P_0`$. Star formation is assumed to be coeval within each OB association, with SNe therefore exploding over a period $`t_e=40`$ Myr, the lifetime of the lowest-mass SN progenitors. These simple assumptions yield two evolutionary paths for the superbubbles: objects whose growth stalls at radius $`R_\mathrm{f}`$ due to pressure confinement at time $`t_\mathrm{f}(n/P_0)^{1/2}R_\mathrm{f}`$; and objects whose growth continues until $`t_e`$. The characteristic timescale $`t_e`$ therefore determines the characteristic parameters $`R_e`$ and $`L_e`$, which correspond to conditions for which an object stalls at $`t_e`$.

The resultant size distribution for a constant superbubble creation rate $`\psi `$, and single-valued $`\varphi (L)`$ is:

$$N_{\mathrm{grow}}(R)R^{2/3}.$$
(3)
This is determined exclusively by growing objects, since the stall radius $`R_\mathrm{f}`$ in this case is the same for all the shells, and represents the maximum $`R`$ in the distribution.

If all objects are created in a single burst, but now with a power-law $`\varphi (L)`$ (equation 2), then both growing and stalled objects contribute to the size distribution:

$`N_{\mathrm{grow}}(R)`$ $``$ $`R^{45\beta }R^3`$
$`N_{\mathrm{stall}}(R)`$ $``$ $`R^{12\beta }R^6,`$ (4)
with $`\beta =2`$. The stalled objects dominate the total $`N(R)`$ for $`RR_\mathrm{f}(t_e)`$, beyond which the growing objects dominate, falling off with a much steeper $`R`$-dependence.

For constant $`\psi `$ and power-law $`\varphi (L)`$, the stalled objects again dominate, with a total $`N(R)`$:

$$N(R)=A\psi \frac{L_e^{1\beta }}{R_e}\left(\frac{R}{R_e}\right)^{12\beta }\left[2(t_e+t_s)+\frac{96\beta }{2+3\beta }t_e\left(\frac{R}{R_e}\right)\right],$$
(5)
where $`t_s`$ is an arbitrary period after $`t_e`$ during the which the superbubbles survive before disappearing. Here again, the stalled objects dominate, effectively yielding $`N(R)R^{12\beta }R^3`$. The maximum in $`N(R)`$ is at $`R_\mathrm{f}(L_{\mathrm{min}})`$, corresponding to single-progenitor SNRs. At $`RR_\mathrm{f}(L_{\mathrm{min}})`$, growing, high-$`L`$ objects dominate with $`N(R)R^{2/3}`$ (equation 3). Oey & Clarke (1998b) show schematic representations of $`N(R)`$ for the above three combinations of $`\psi `$ and $`\varphi (L)`$.

We compared this result to the H i shell catalog for the Small Magellanic Cloud (SMC) compiled by Staveley-Smith et al. (1997). This is by far the most complete sample of H i shells obtained for any galaxy, as evidenced by the fact that the relative number counts of H ii regions and H i shells are in excellent agreement with their relative life expectancies. For shells having $`R100`$ pc, the fitted power-law slope $`\alpha =12\beta `$ is $`2.7\pm 0.6`$, in excellent agreement with the general prediction of $`\alpha =3`$. Using $`\beta =1.9`$ from the observed H ii LF for the SMC (Kennicutt et al. 1989) yields a specific SMC prediction of $`\alpha =2.8\pm 0.4`$, still in remarkable agreement with the observed value (Oey & Clarke 1997).

We note that different models for ISM structure yield different predictions for $`N(R)`$. For example, Stanimirović et al. (1999) suggest a possible fractal structure for the neutral ISM. From the same H i dataset of the SMC, they find a fractal dimension implying a size distribution for H i holes of $`\alpha =3.5`$. Within the uncertainties, this prediction cannot be empirically differentiated from our model with $`\alpha =3`$; but it is worth noting that the predictions are in fact intrinsically different.

However, the superbubble size distribution presently is not a sensitive test in determining whether or not the objects evolve adiabatically. If all the internal energy is radiated away, the objects are predicted to follow the momentum-conserving law given by Steigman et al. (1975):

$$R(L/nv_{\mathrm{}})^{1/4}t^{1/2}.$$
(6)
The stall radius $`R_\mathrm{f}`$ in this case is only 1.3 times larger than for the adiabatic model. Furthermore, the size distribution follows the same law, $`N(R)R^{12\beta }`$ (Oey & Clarke 1997). The observations of hot gas are therefore vital confirmation that the adiabatic model applies to a significant fraction of superbubbles.

It is also possible to derive the distribution of expansion velocities $`N(v)dv`$. This is naturally determined only by the growing objects (Oey & Clarke 1998b):

$$N_{\mathrm{grow}}(v)v^{7/2},\beta >1.5.$$
(7)
We again compared this prediction to the SMC H i shell catalog, and found a fitted power-law slope of $`2.9\pm 1.4`$, which is also in agreement with the prediction. Thus, in spite of the crude assumptions in deriving the shell size and velocity distributions, the data suggest that the neutral ISM in the SMC is fully consistent with superbubble activity dominating the structure. Although most other available H i shell catalogs are highly incomplete, preliminary results for a few other galaxies also show agreement with our model for the size distribution (Kim et al. 1999; Mashchenko et al. 1999; Oey & Clarke 1997).

### 4.2 ISM porosity and galactic superwinds

It is straightforward to use the analytic expression for $`N(R)`$ (equation 5) to derive the interstellar porosity parameter $`Q`$, which is given by Oey & Clarke (1997) for 2- and 3-dimensional cases. $`Q`$ is the ratio of: total area or volume occupied by superbubbles, to the total area or volume of the galaxy. Thus it is essentially the filling factor of hot gas, assuming hot gas is contained within all of the superbubbles. Determining $`Q`$ has been the conventional way to evaluate the relative importance of the HIM (e.g., Cox & Smith 1974; McKee & Ostriker 1977; Heiles 1990). Values near unity indicate the HIM dominates the multiphase ISM, and values significantly $`>1`$ imply an outflow, with the galaxy generating more hot gas than it can contain.

We can write $`Q`$ in terms of a galaxy’s star-formation rate (SFR), $`\mathrm{\Psi }`$:

$$Q16\frac{\mathrm{\Psi }(\mathrm{M}_{}\mathrm{yr}^1)}{hR_g^2(\mathrm{kpc}^3)}\frac{1}{P_0},$$
(8)
for $`\beta =2`$, a Salpeter (1955) IMF for stellar masses $`0.1m100\mathrm{M}_{}`$, and $`P_0/k=9500`$. $`R_g`$ and $`h`$ are the galaxy’s gaseous radius and scale height, respectively. We caution that $`Q`$ depends on ambient interstellar parameters, for example, $`P_0^1`$ as indicated.

Table 1 shows the interstellar porosities for all the galaxies in the Local Group. Dwarf galaxy parameters are taken from Mateo (1998), and SFR’s for M31, M33, and the Magellanic Clouds from their H ii LF’s (Kennicutt et al. 1989). In general, we see that almost all the galaxies show low values with $`Q<0.5`$, implying that the HIM is not the dominant ISM component in these cases. The LMC is an exception, showing $`Q1`$, so the HIM appears likely to dominate in this galaxy. The Milky Way results historically have been ambiguous and controversial (e.g., McKee & Ostriker 1977; Slavin & Cox 1993), and we obtain different answers from different methods. Taking the Galaxy’s SFR from the H ii LF of McKee & Williams (1997) and Smith & Kennicutt (1989) yields a low $`Q0.2`$; but using the SN rate of van den Bergh & Tammann (1991) yields $`Q\mathrm{}>1`$.

However, by far the highest porosity is found for IC 10, the one starburst galaxy in the Local Group. IC 10 shows a $`Q`$ value that is over an order of magnitude larger than for any other galaxy examined, therefore fulfilling the breakout criterion for a galactic superwind outflow. It is possible to estimate the mass-loss rate in this outflow $`\text{ }\mathrm{}M_{\mathrm{out}}`$, assuming that the material is largely evaporated from shell walls by thermal conduction. Mac Low & McCray (1988) have estimated the evaporation rate for an individual superbubble:

$$\text{ }\mathrm{}M(R)16\pi \mu m_\mathrm{H}/25kCT^{5/2}R,$$
(9)
where $`\mu m_\mathrm{H}`$ is the mean mass per nucleon and $`C`$ is a conductivity coefficient. Integrating this over the superbubble size distribution yields $`\text{ }\mathrm{}M_{\mathrm{out}}`$:

$$\text{ }\mathrm{}M_{\mathrm{out}}_{R_{\mathrm{min}}}^{R_{\mathrm{max}}}\text{ }\mathrm{}M(R)N(R)𝑑R.$$
(10)
The estimate of $`\text{ }\mathrm{}M_{\mathrm{out}}`$ is necessarily crude, for example, we rather arbitrarily take $`R_{\mathrm{max}}=R_e`$, with the caveat that it may be limited by $`R_g`$.

The maximum $`\text{ }\mathrm{}M_{\mathrm{out}}`$ is obtained for the largest evaporating surface area of superbubbles, which occurs for $`Q1`$. Porosities of $`Q>1`$ imply strong interaction and destruction of the shells, or blowouts, probably reducing the available surface area. For IC 10, the derived normalization for $`N(R)`$ yielding $`Q=1`$ corresponds to $`\text{ }\mathrm{}M_{\mathrm{out}}9\times 10^3\mathrm{M}_{}\mathrm{yr}^1`$. For comparison, the SFR required for $`Q=1`$ (equation 8) is $`\mathrm{\Psi }3\times 10^2\mathrm{M}_{}\mathrm{yr}^1`$, which is only about 3 times more than $`\text{ }\mathrm{}M_{\mathrm{out}}`$! Although the actual SFR for IC 10 is much larger ($`7\times 10^1\mathrm{M}_{}\mathrm{yr}^1`$), yielding the higher original $`Q`$ value, this rough calculation suggests that $`\text{ }\mathrm{}M_{\mathrm{out}}`$ from galactic winds can be of similar magnitude as the parent starburst SFR. Indeed, an extensive absorption-line study of local starburst galaxies by Heckman et al. (2000) also concluded that empirically, the outflow and star-formation rates have the same order of magnitude for their sample.

Finally, we note the extensive body of numerical work on superbubbles and blowout conditions. These are too numerous to discuss here, but we refer the reader to Mac Low (1999) and Strickland & Stevens (2000), who provide overviews of this field. The details of the numerical predictions are presently difficult to confirm empirically, but observations with XMM-Newton and Chandra will be especially helpful to constrain the dominant processes.

## 5 Summary

We have seen that observations of mechanical feedback ranging from individual stellar wind bubbles to galactic superwinds are all roughly consistent with the conventional adiabatic model for bubble evolution. The empirical confirmation of hot, coronal gas supports the relevance of the adiabatic evolution for at least an important proportion of shell phenomena. Although kinematic discrepancies are often found for individual objects, these anomalies have plausible explanations, and presently there are no cases that suggest a need for any major revision of the conventional understanding. It is a remarkable strength that the adiabatic model apparently applies to phenomena over a range of size scales covering at least three orders of magnitude, from pc to kpc.

However, characterizing the dominant parameters and their effects on the shell evolution is still highly problematic. For example, critical ambient ISM conditions like the density, pressure, and ionization distributions remain elusive. The mechanisms and conditions for cooling of the interior energy need to be identified, and the energy budgets reliably determined. Perhaps the most fundamental question is the fate of the hot gas generated within the superbubbles: Does it escape and thereby constitute the HIM? Does it indeed escape from starburst galaxies, and from their gravitational potentials? These issues have vital consequences for almost all galactic evolutionary processes, and our understanding depends on further clarifying the bubble and superbubble process.

MSO gratefully acknowledges travel support from the Graduiertenkolleg and STScI Institute Fellowship. Some material in §4.2 was derived in 1998 at the Aspen Center for Physics.

References

Arthur, S. J. & Henney, W. J. 1996, ApJ, 457, 752

Benaglia, P. & Cappa, C. E. 1999, A&A, 346, 979

Brown, A. G. A., Hartmann, D., & Burton, W. B. 1995, A&A, 300, 903

Cappa, C. E. & Benaglia, P. 1998, AJ, 116, 1906

Chu, Y.-H., Chang, H.-W., Su, Y.-L., & Mac Low, M.-M., 1995, ApJ, 450, 157

Chu, Y.-H. & Mac Low, M-M., 1990, ApJ, 365, 510

Chu, Y.-H., Wakker, B., Mac Low, M.-M., & García-Segura, G. 1994, AJ, 108, 1696

Cowie, L. L. & McKee, C. F. 1977, ApJ, 211, 135

Cox, D. P. & Smith, B. W. 1974, ApJ, 189, L105

Drissen, L., Moffat, A. F. J., Walborn, N. R., & Shara, M. R. 1995, AJ, 110, 2235

García-Segura, G. & Franco, J. 1996, ApJ, 469, 171

García-Segura, G. & Mac Low, M.-M. 1995, ApJ, 455, 145

Heckman, T. M., Lehnert, M. D., Strickland, D. K., & Armus, L. 2000, ApJS, 129, 493

Heiles, C. 1990, ApJ, 354, 483

Kennicutt, R. C., Edgar, B. K., & Hodge, P. W. 1989, ApJ, 337, 761

Kim, S., Dopita, M. A., Staveley-Smith, L., & Bessell, M. S. 1999, AJ, 118, 2797

Mac Low, M.-M., in New Perspectives on the Interstellar Medium, eds. A. R. Taylor, T. L. Landecker, & G. Joncas, San Francisco: Astron. Soc. Pacific, 303

Mac Low, M.-M., Chang, T. H., Chu, Y.-H., Points, S. D., Smith, R. C., & Wakker, B. P. 1998, ApJ, 493, 260

Mac Low, M.-M. & McCray, R. 1988, ApJ, 324, 776

Magnier, E. A., Chu, Y.-H., Points, S. D., Hwang, U., & Smith, R. C. 1996, ApJ, 464, 829

Mashchenko, S. Y., Thilker, D. A., & Braun, R. 1999, A&A, 343, 352

Mateo, M. L. 1998, ARAA, 36, 435

McKee, C. F. & Ostriker, J. P. 1977, ApJ, 218, 148

McKee, C. F., Van Buren, D., & Lazareff, B. 1984, ApJ, 278, L115

McKee, C. F. & Williams, J. P. 1997, ApJ, 476, 144

Oey, M. S., 1996, ApJ, 467, 666

Oey, M. S. & Clarke, C. J. 1997, MNRAS, 289, 570

Oey, M. S. & Clarke, C. J. 1998a, AJ, 115, 1543

Oey, M. S. & Clarke, C. J. 1998b, in Interstellar Turbulence, eds. J. Franco & A. Carramiñana, Cambridge: Cambridge Univ. Press, 112

Oey, M. S., Groves, B., Staveley-Smith, L., & Smith, R. C. 2001, in preparation

Oey, M. S. & Massey, P., 1994, ApJ, 425, 635

Oey, M. S. & Massey, P., 1995, ApJ, 452, 210

Oey, M. S. & Smedley, S. A. 1998, AJ, 116, 1263

Pikel’ner, S. B. 1968, Astrophys. Lett., 2, 97

Saken, J. M., Shull, J. M., Garmany, C. D., Nichols-Bohlin, J., & Fesen, R. A. 1992, ApJ, 397, 537

Salpeter, E. E. 1955, ApJ, 121, 161

Silich, S. & Franco, J. 1999, ApJ, 522, 863

Silich, S. A., Tenorio-Tagle, G., Terlevich, R., Terlevich, E., & Netzer, H. 2001, MNRAS in press, astro-ph/0009092

Slavin, J. D. & Cox, D. P. 1993, ApJ, 417, 187

Smith, T. R. & Kennicutt, R. C. 1989, PASP, 101, 649

Stanimirović, S., Staveley-Smith, L., Dickey, J. M., Sault, R. J., & Snowden, S. L. 1999, MNRAS, 302, 417

Staveley-Smith, L., Sault, R. J., Hatzidimitriou, D., Kesteven, M. J., & McConnell, D. 1997, MNRAS 289, 225

Steigman, G., Strittmatter, P. A., & Williams, R. E. 1975, ApJ, 198, 575

Strickland, D. K. & Stevens, I. R. 2000, MNRAS, 314, 511

Treffers, R. R. & Chu, Y.-H. 1982, ApJ, 254, 569

van den Bergh, S. & Tammann, G. A. 1991, ARAA, 29, 363

Wang, Q. & Helfand, D. J. 1991, ApJ, 373, 497

Weaver, R., McCray, R., Castor, J., Shapiro, P., & Moore, R. 1977, ApJ, 218, 377